From a8046f80d0b11d8fe55e71c03a422533c1e9e8a5 Mon Sep 17 00:00:00 2001 From: "Bastian M. Wojek" Date: Sat, 24 May 2008 15:05:18 +0000 Subject: [PATCH] Fixed bugs and added a few comments after first testing. --- .../TFitPofB-lib/classes/TBofZCalc.cpp | 82 ++++-- .../TFitPofB-lib/classes/TFitPofB.cpp | 144 ++++++--- .../TFitPofB-lib/classes/TPofBCalc.cpp | 85 ++++-- .../TFitPofB-lib/classes/TPofTCalc.cpp | 106 ++++--- .../classes/TTrimSPDataHandler.cpp | 60 ++-- src/external/TFitPofB-lib/include/TBofZCalc.h | 12 +- src/external/TFitPofB-lib/include/TFitPofB.h | 6 +- src/external/TFitPofB-lib/include/TPofBCalc.h | 6 +- src/external/TFitPofB-lib/include/TPofTCalc.h | 5 +- .../TFitPofB-lib/include/TTrimSPDataHandler.h | 2 +- src/external/TFitPofB-lib/test/Makefile.test | 33 +++ src/external/TFitPofB-lib/test/test.cpp | 276 ++++++++++++++++++ 12 files changed, 642 insertions(+), 175 deletions(-) create mode 100644 src/external/TFitPofB-lib/test/Makefile.test create mode 100644 src/external/TFitPofB-lib/test/test.cpp diff --git a/src/external/TFitPofB-lib/classes/TBofZCalc.cpp b/src/external/TFitPofB-lib/classes/TBofZCalc.cpp index 958c9be6..21f976b5 100644 --- a/src/external/TFitPofB-lib/classes/TBofZCalc.cpp +++ b/src/external/TFitPofB-lib/classes/TBofZCalc.cpp @@ -5,7 +5,7 @@ Author: Bastian M. Wojek e-mail: bastian.wojek@psi.ch - 2008/05/23 + 2008/05/24 ***************************************************************************/ @@ -14,23 +14,34 @@ #include #include +//------------------ +// Method returning the field minimum of the model +//------------------ double TBofZCalc::GetBmin() const { vector::const_iterator iter; iter = min_element(fBZ.begin(),fBZ.end()); - + return *iter; } +//------------------ +// Method returning the field maximum of the model +//------------------ + double TBofZCalc::GetBmax() const { vector::const_iterator iter; iter = max_element(fBZ.begin(),fBZ.end()); - + return *iter; } +//------------------ +// Method returning the field B(z) at a given z according to the model +//------------------ + double TBofZCalc::GetBofZ(double zz) const { - + bool found = false; unsigned int i; for (i=0; i ¶m) { - - // Parameters for 1D-London screening in a thin superconducting film - // Bext[G], deadlayer[nm], thickness[nm], lambda[nm] - + unsigned int n(2000); // number of steps for the calculation - + double N(cosh(param[2]/2.0/param[3])); - + fDZ = param[2]/double(n); double ZZ, BBz; - + for (unsigned int j(0); j ¶m) { - - // Parameters for 1D-London screening in a thin superconducting film, two layers, two lambda - // Bext[G], deadlayer[nm], thickness1[nm], thickness2[nm], lambda1[nm], lambda2[nm] - - unsigned int n(2000); // number of steps for the calculation - + + unsigned int n(4000); // number of steps for the calculation + double N1(param[5]*cosh(param[3]/param[5])*sinh(param[2]/param[4]) + param[4]*cosh(param[2]/param[4])*sinh(param[3]/param[5])); double N2(4.0*N1); - + fDZ = (param[2]+param[3])/double(n); double ZZ, BBz; - + for (unsigned int j(0); j ¶m) { } fBZ.push_back(BBz); } - + } +//------------------ +// Constructor of the TLondon1D_3L class +// 1D-London screening in a thin superconducting film, three layers, two lambdas (top and bottom layer: lambda1) +// Parameters: Bext[G], deadlayer[nm], thickness1[nm], thickness2[nm], thickness3[nm], lambda1[nm], lambda2[nm] +//------------------ + TLondon1D_3L::TLondon1D_3L(const vector ¶m) { - - // Parameters for 1D-London screening in a thin superconducting film, three layers, two lambdas (top and bottom layer: lambda1) - // Bext[G], deadlayer[nm], thickness1[nm], thickness2[nm], thickness3[nm], lambda1[nm], lambda2[nm] - - unsigned int n(2000); // number of steps for the calculation - + + unsigned int n(4000); // number of steps for the calculation + double N1(8.0*(param[5]*param[6]*cosh(param[3]/param[6])*sinh((param[2]+param[4])/param[5]) + ((param[5]*param[5]*cosh(param[2]/param[5])*cosh(param[4]/param[5])) + (param[6]*param[6]*sinh(param[2]/param[5])*sinh(param[4]/param[5])))*sinh(param[3]/param[6]))); - + double N2(2.0*param[5]*param[6]*cosh(param[3]/param[6])*sinh((param[2]+param[4])/param[5]) + 2.0*(param[5]*param[5]*cosh(param[2]/param[5])*cosh(param[4]/param[5]) + param[6]*param[6]*sinh(param[2]/param[5])*sinh(param[4]/param[5]))*sinh(param[3]/param[6])); - + double N3(8.0*(param[5]*param[6]*cosh(param[3]/param[6])*sinh((param[2]+param[4])/param[5]) + (param[5]*param[5]*cosh(param[2]/param[5])*cosh(param[4]/param[5]) + param[6]*param[6]*sinh(param[2]/param[5])*sinh(param[4]/param[5]))*sinh(param[3]/param[6]))); - + fDZ = (param[2]+param[3]+param[4])/double(n); double ZZ, BBz; - + for (unsigned int j(0); j using namespace std; +//------------------ +// Constructor of the TFitPofB class -- reading available implantation profiles and +// creates (a pointer to) the TPofTCalc object (with the FFT plan) +//------------------ + TFitPofB::TFitPofB(const vector &parNo, const vector &par) : fCalcNeeded(true) { - + for(unsigned int i(0); i energy_vec(energy_arr, energy_arr+(sizeof(energy_arr)/sizeof(energy_arr[0]))); vector par_for_PofT; - + for (unsigned int i(1); i<4; i++) par_for_PofT.push_back(fPar[i]); @@ -33,7 +38,7 @@ TFitPofB::TFitPofB(const vector &parNo, const vector &par) fImpProfile = x; x = 0; delete x; - + TPofTCalc *y = new TPofTCalc(par_for_PofT); fPofT = y; y = 0; @@ -41,81 +46,128 @@ TFitPofB::TFitPofB(const vector &parNo, const vector &par) } +//------------------ +// Destructor of the TFitPofB class -- cleaning up +//------------------ + TFitPofB::~TFitPofB() { fPar.clear(); - fImpProfile = 0; delete fImpProfile; - fPofT = 0; + fImpProfile = 0; delete fPofT; + fPofT = 0; } +//------------------ +// Method that calls the procedures to create B(z), p(B) and P(t) +// It finally returns P(t) for a given t. +// Parameters: all the parameters for the function to be fitted through TFitPofB +//------------------ + double TFitPofB::Eval(double t, const vector &par) const { // check if any parameter has changed bool par_changed(false); + bool only_phase_changed(false); for (unsigned int i(0); i par_for_BofZ; - vector par_for_PofB; vector par_for_PofT; - for (unsigned int i(5); iDoFFT(PofB1, par_for_PofT); + if(!only_phase_changed) { + + cout << " Parameters have changed, (re-)calculating p(B) and P(t) now..." << endl; + + vector par_for_BofZ; + vector par_for_PofB; + +// cout << "par_for_BofZ: "; + + for (unsigned int i(5); iDoFFT(PofB2, par_for_PofT); +// cout << endl; + +// cout << "par_for_PofB: "; + + for (unsigned int i(2); i<5; i++) { + par_for_PofB.push_back(par[i]); +// cout << par[i] << " "; } - break; - case 3: - { - TLondon1D_3L BofZ3(par_for_BofZ); - TPofBCalc PofB3(BofZ3, *fImpProfile, par_for_PofB); - fPofT->DoFFT(PofB3, par_for_PofT); - } - break; - default: - cout << "The user function you specified, does not exist. Quitting..." << endl; - exit(0); +// cout << endl; + + switch (int(par[0])) { + case 1: + { +// cout << "Found the 1D-London model." << endl; + TLondon1D_1L BofZ1(par_for_BofZ); + TPofBCalc PofB1(BofZ1, *fImpProfile, par_for_PofB); + fPofT->DoFFT(PofB1); + } + break; + case 2: + { +// cout << "Found the 1D-London model.2L" << endl; + TLondon1D_2L BofZ2(par_for_BofZ); + TPofBCalc PofB2(BofZ2, *fImpProfile, par_for_PofB); + fPofT->DoFFT(PofB2); + } + break; + case 3: + { +// cout << "Found the 1D-London model.3L" << endl; + TLondon1D_3L BofZ3(par_for_BofZ); + TPofBCalc PofB3(BofZ3, *fImpProfile, par_for_PofB); + fPofT->DoFFT(PofB3); + } + break; + default: + cout << "The user function you specified with the first parameter, does not exist. Quitting..." << endl; + exit(0); } + } else { + cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl; + } + + fPofT->CalcPol(par_for_PofT); + + fCalcNeeded = false; } - + return fPofT->Eval(t); - + } diff --git a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp index bed045dd..4999ce7b 100644 --- a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp +++ b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp @@ -5,47 +5,70 @@ Author: Bastian M. Wojek e-mail: bastian.wojek@psi.ch - 2008/05/23 + 2008/05/24 ***************************************************************************/ #include "TPofBCalc.h" #include #include +#include +//#include + +//----------- +// Constructor that does the P(B) calculation for given B(z) and n(z) +// Parameters: dt[us], dB[G], Energy[keV] +//----------- TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, const vector ¶ ) { - // Parameters: dt[us], dB[G], Energy[keV] - + fBmin = BofZ.GetBmin(); + fBmax = BofZ.GetBmax(); + double BB, BBnext; double zm, zp, zNextm, zNextp, dz; // fill not used Bs before Bmin with 0.0 - - for ( BB = 0.0 ; BB < BofZ.GetBmin() ; BB += para[1] ) { + + for ( BB = 0.0 ; BB < fBmin ; BB += para[1] ) { fB.push_back(BB); fPB.push_back(0.0); } - + unsigned int firstZerosEnd(fB.size()); - + // calculate p(B) from B(z) - + vector bofzZ(BofZ.DataZ()); vector bofzBZ(BofZ.DataBZ()); double ddZ(BofZ.GetdZ()); - double Bmax(BofZ.GetBmax()); - + +/* USED FOR DEBUGGING----------------------------------- + cout << "Bmin = " << fBmin << ", Bmax = " << fBmax << endl; + + char debugfile[50]; + int n = sprintf (debugfile, "test_Bz_%f.dat", fBmin); + + if (n > 0) { + ofstream of(debugfile); + + for (unsigned int i(0); i= BB && bofzBZ[j+1] <= BB ) { zm = (BB-bofzBZ[j])*ddZ/(bofzBZ[j+1]-bofzBZ[j]) + bofzZ[j]; - + for (unsigned int k(0); k < j; k++) { if ( ( bofzBZ[j-k] <= BBnext && bofzBZ[j-k-1] >= BBnext ) ) { // cout << "1 " << j << " " << k << endl; @@ -53,17 +76,17 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons break; } } - + dz = zNextm-zm; nn = dataTrimSP.GetNofZ(zm, para[2]); if (nn != -1.0) { // cout << "zNext = " << zNextm << ", zm = " << zm << ", dz = " << dz << endl; *(fPB.end()-1) += nn*fabs(dz/para[1]); } - + } else if (bofzBZ[j] <= BB && bofzBZ[j+1] >= BB ) { zp = (BB-bofzBZ[j])*ddZ/(bofzBZ[j+1]-bofzBZ[j]) + bofzZ[j]; - + for (unsigned int k(0); k < bofzZ.size() - j - 1; k++) { if ( ( bofzBZ[j+k] <= BBnext && bofzBZ[j+k+1] >= BBnext ) ) { // cout << "2 " << j << " " << k << endl; @@ -71,7 +94,7 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons break; } } - + dz = zNextp-zp; nn = dataTrimSP.GetNofZ(zp, para[2]); if (nn != -1.0) { @@ -80,33 +103,33 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons } } } - + } - + unsigned int lastZerosStart(fB.size()); - + // fill not used Bs after Bext with 0.0 - + double BmaxFFT(1.0/gBar/para[0]); - -//// cout << "N = " << int(BmaxFFT/para[1]+1.0) << endl; - + +// cout << "N = " << int(BmaxFFT/para[1]+1.0) << endl; + for ( ; BB <= BmaxFFT ; BB += para[1] ) { fB.push_back(BB); fPB.push_back(0.0); } - - // make sure that we have an even number of elements in p(B) for FFT - + + // make sure that we have an even number of elements in p(B) for FFT, so we do not have to care later + if (fB.size() % 2) { fB.push_back(BB); fPB.push_back(0.0); } - - cout << "size of B = " << fB.size() << ", size of p(B) = " << fPB.size() << endl; - // normalize pB - +// cout << "size of B = " << fB.size() << ", size of p(B) = " << fPB.size() << endl; + + // normalize p(B) + double pBsum = 0.0; for (unsigned int i(firstZerosEnd); i #include +#include #include +//------------------ +// Constructor of the TPofTCalc class - it creates the FFT plan +// Parameters: phase, dt, dB +//------------------ + TPofTCalc::TPofTCalc (const vector &par) { fNFFT = ( int(1.0/gBar/par[1]/par[2]+1.0) % 2 ) ? int(1.0/gBar/par[1]/par[2]+2.0) : int(1.0/gBar/par[1]/par[2]+1.0); fTBin = 1.0/gBar/double(fNFFT-1)/par[2]; -// tBin = 1.0/gBar/(*(PofB.DataB().end()-1)); -// nFFT = PofB.DataB().size(); - + fFFTin = (double *)malloc(sizeof(double) * fNFFT); fFFTout = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (fNFFT/2+1)); - + cout << "Check for the FFT plan..." << endl; - + // Load wisdom from file - + int wisdomLoaded(0); - + FILE *wordsOfWisdomR; wordsOfWisdomR = fopen("WordsOfWisdom.dat", "r"); if (wordsOfWisdomR == NULL) { @@ -39,19 +43,40 @@ TPofTCalc::TPofTCalc (const vector &par) { wisdomLoaded = fftw_import_wisdom_from_file(wordsOfWisdomR); fclose(wordsOfWisdomR); } - + if (!wisdomLoaded) { cout << "No wisdom is imported..." << endl; } - + fFFTplan = fftw_plan_dft_r2c_1d(fNFFT, fFFTin, fFFTout, FFTW_EXHAUSTIVE); - + } -void TPofTCalc::DoFFT(const TPofBCalc &PofB, const vector &par) { - +//-------------- +// Method that does the FFT of a given p(B) +//-------------- + +void TPofTCalc::DoFFT(const TPofBCalc &PofB) { + vector pB(PofB.DataPB()); - + +/* USED FOR DEBUGGING ----------------------- + vector B(PofB.DataB()); + double Bmin(PofB.GetBmin()); + + char debugfile[50]; + int n = sprintf (debugfile, "test_PB_%f.dat", Bmin); + + if (n > 0) { + ofstream of(debugfile); + + for (unsigned int i(0); i &par) { fFFTout[i][0] = 0.0; fFFTout[i][1] = 0.0; } - + cout << "perform the Fourier transform..." << endl; - + fftw_execute(fFFTplan); - // Calculate polarisation - +} + +//--------------------- +// Method for calculating the muon spin polarization P(t) from the Fourier transformed p(B) +// Parameters: phase, dt, dB +//--------------------- + +void TPofTCalc::CalcPol(const vector &par) { + double sinph(sin(par[0]*PI/180.0)), cosph(cos(par[0]*PI/180.0)); - - double polTemp(0.0); - + + double polTemp(0.0); + fT.clear(); + fPT.clear(); + for (unsigned int i(0); i energyVec(energyArr, energyArr+(sizeof(energyArr)/sizeof(energyArr[0]))); +// +// This will read the files "/home/user/TrimSP/SomeSample-02_1.rge", "/home/user/TrimSP/SomeSample-02_5.rge" and so on. +//-------------------- TTrimSPData::TTrimSPData(const string &path, vector &energyVec) { - //string energyArr[] = {"02_1", "02_5", "03_5", "05_0", "07_5", "10_0", "12_5", "15_0", "17_5", "19_0", "20_0", "22_5", "25_0"}; - double zz(0.0), nzz(0.0); vector vzz, vnzz; string word, energyStr; - + for(unsigned int i(0); i> word) if(word == "PARTICLES") break; @@ -49,57 +53,63 @@ TTrimSPData::TTrimSPData(const string &path, vector &energyVec) { vzz.push_back(zz); vnzz.push_back(nzz); } - + fDataZ.push_back(vzz); fDataNZ.push_back(vnzz); - + rgeFile->close(); delete rgeFile; rgeFile = 0; vzz.clear(); vnzz.clear(); - + } } } -// vector DataZ(double) -- returning z-vector calculated by trim.SP for energy given energy +//--------------------- +// Method returning z-vector calculated by trim.SP for given energy[keV] +//--------------------- vector TTrimSPData::DataZ(double e) const { - + for(unsigned int i(0); i DataNZ(double) -- returning n(z)-vector calculated by trim.SP for given energy +//--------------------- +// Method returning n(z)-vector calculated by trim.SP for given energy[keV] +//--------------------- vector TTrimSPData::DataNZ(double e) const { - + for(unsigned int i(0); i z, nz; - + for(unsigned int i(0); i using namespace std; +//-------------------- // Base class for any kind of theory function B(z) +// In principle only constructors for the different models have to be implemented +//-------------------- class TBofZCalc { @@ -41,7 +44,9 @@ protected: double fDZ; }; +//-------------------- // Class "for Meissner screening" in a thin superconducting film +//-------------------- class TLondon1D_1L : public TBofZCalc { @@ -51,8 +56,9 @@ public: }; - +//-------------------- // Class "for Meissner screening" in a thin superconducting film - bilayer with two different lambdas +//-------------------- class TLondon1D_2L : public TBofZCalc { @@ -62,7 +68,9 @@ public: }; +//-------------------- // Class "for Meissner screening" in a thin superconducting film - tri-layer with two different lambdas +//-------------------- class TLondon1D_3L : public TBofZCalc { diff --git a/src/external/TFitPofB-lib/include/TFitPofB.h b/src/external/TFitPofB-lib/include/TFitPofB.h index b8a275d0..a968ec65 100644 --- a/src/external/TFitPofB-lib/include/TFitPofB.h +++ b/src/external/TFitPofB-lib/include/TFitPofB.h @@ -5,7 +5,7 @@ Author: Bastian M. Wojek e-mail: bastian.wojek@psi.ch - 2008/05/23 + 2008/05/24 ***************************************************************************/ @@ -19,9 +19,9 @@ class TFitPofB { public: TFitPofB(const vector& , const vector&); ~TFitPofB(); - + double Eval(double, const vector&) const; - + private: mutable vector fPar; TTrimSPData *fImpProfile; diff --git a/src/external/TFitPofB-lib/include/TPofBCalc.h b/src/external/TFitPofB-lib/include/TPofBCalc.h index 8268c3bf..8730706e 100644 --- a/src/external/TFitPofB-lib/include/TPofBCalc.h +++ b/src/external/TFitPofB-lib/include/TPofBCalc.h @@ -5,7 +5,7 @@ Author: Bastian M. Wojek e-mail: bastian.wojek@psi.ch - 2008/05/23 + 2008/05/24 ***************************************************************************/ @@ -27,10 +27,14 @@ public: vector DataB() const {return fB;} vector DataPB() const {return fPB;} + double GetBmin() const {return fBmin;} + double GetBmax() const {return fBmax;} private: vector fB; vector fPB; + double fBmin; + double fBmax; static const double gBar = 0.0135538817; }; diff --git a/src/external/TFitPofB-lib/include/TPofTCalc.h b/src/external/TFitPofB-lib/include/TPofTCalc.h index 23f2fc48..d48e35ed 100644 --- a/src/external/TFitPofB-lib/include/TPofTCalc.h +++ b/src/external/TFitPofB-lib/include/TPofTCalc.h @@ -5,7 +5,7 @@ Author: Bastian M. Wojek e-mail: bastian.wojek@psi.ch - 2008/05/23 + 2008/05/24 ***************************************************************************/ @@ -23,7 +23,8 @@ public: vector DataT() const {return fT;} vector DataPT() const {return fPT;} - void DoFFT(const TPofBCalc&, const vector&); + void DoFFT(const TPofBCalc&); + void CalcPol(const vector&); double Eval(double) const; private: diff --git a/src/external/TFitPofB-lib/include/TTrimSPDataHandler.h b/src/external/TFitPofB-lib/include/TTrimSPDataHandler.h index c163e8f9..6510e2c9 100644 --- a/src/external/TFitPofB-lib/include/TTrimSPDataHandler.h +++ b/src/external/TFitPofB-lib/include/TTrimSPDataHandler.h @@ -5,7 +5,7 @@ Author: Bastian M. Wojek e-mail: bastian.wojek@psi.ch - 2008/05/23 + 2008/05/24 ***************************************************************************/ diff --git a/src/external/TFitPofB-lib/test/Makefile.test b/src/external/TFitPofB-lib/test/Makefile.test new file mode 100644 index 00000000..3b4e826c --- /dev/null +++ b/src/external/TFitPofB-lib/test/Makefile.test @@ -0,0 +1,33 @@ +CXX = g++ +CXXFLAGS = -g -Wall +LOCALINCLUDE = ../include/ +INCLUDES = -I$(LOCALINCLUDE) +LD = g++ +LDFLAGS = -g -L../classes/ -lTFitPofB -lfftw3 -lm + +EXEC = test + +# some definitions: headers, sources, objects,... +OBJS = +OBJS += $(EXEC).o + +# make the executable: +# +all: $(EXEC) + +$(EXEC): $(OBJS) + @echo "---> Building $(EXEC) ..." + $(LD) $(LDFLAGS) $(OBJS) -o $(EXEC) + @echo "done" + +# clean up: remove all object file (and core files) +# semicolon needed to tell make there is no source +# for this target! +# +clean:; @rm -f $(OBJS) + @echo "---> removing $(OBJS)" + +# +$(OBJS): %.o: %.cpp + $(CXX) $(INCLUDES) $(CXXFLAGS) -c $< + diff --git a/src/external/TFitPofB-lib/test/test.cpp b/src/external/TFitPofB-lib/test/test.cpp new file mode 100644 index 00000000..3d39b88f --- /dev/null +++ b/src/external/TFitPofB-lib/test/test.cpp @@ -0,0 +1,276 @@ +#include "TFitPofB.h" +#include +#include + +using namespace std; + +int main(){ +/* string rge_path("/home/l_wojek/nt/wojek/g/Bastian/ImplantationDepth/YBCO_PBCO-"); + string energy_arr[] = {"02_1", "02_5", "03_5", "05_0", "07_5", "10_0", "12_5", "15_0", "17_5", "19_0", "20_0", "22_5", "25_0"}; + + vector energy_vec(energy_arr, energy_arr+(sizeof(energy_arr)/sizeof(energy_arr[0]))); + + TTrimSPData calcData(rge_path, energy_vec); + + vector energies(calcData.Energy()); + for (unsigned int i(0); i z(calcData.DataZ(2.5)); + vector nz(calcData.DataNZ(2.5)); + + vector z2(calcData.DataZ(25.0)); + vector nz2(calcData.DataNZ(25.0)); + + ofstream of("test_out1.dat"); + for (unsigned int i(0); i parameter(param,param+8); + vector param_for_BofZ; + vector param_for_PofB; + vector param_for_PofT; + + for (unsigned int i(0); i<4; i++) + param_for_BofZ.push_back(parameter[i]); + + for (unsigned int i(5); i<8; i++) + param_for_PofB.push_back(parameter[i]); + + for (unsigned int i(4); i<7; i++) + param_for_PofT.push_back(parameter[i]); + + TLondon1D_1L bofz(param_for_BofZ); + + cout << "Bmin = " << bofz.GetBmin() << endl; + cout << "Bmax = " << bofz.GetBmax() << endl; + cout << "dZ = " << bofz.GetdZ() << endl; + + ofstream of5("test_Bz.dat"); + for (double i(0); i<50.; i+=0.1) { + of5 << i << " " << bofz.GetBofZ(i) << endl; + } + of5.close(); + + ofstream of6("test_zBz.dat"); + for (unsigned int i(0); i hurgaB(pofb.DataB()); + vector hurgaPB(pofb.DataPB()); + + ofstream of7("test_BpB.dat"); + for (unsigned int i(0); i parNo_vec(parNo_arr, parNo_arr+(sizeof(parNo_arr)/sizeof(parNo_arr[0]))); + vector par_vec(par_arr, par_arr+(sizeof(par_arr)/sizeof(par_arr[0]))); + + vector par_vec_sub; + + for(unsigned int i(0); i data01, data02, data03, data04, data05, data06, data07, data08, data09, data10; + + for (double i(0.); i<12.0; i+=0.003) + data01.push_back(fitter.Eval(i, par_vec_sub)); + + + par_vec_sub[1] += 10.0; + par_vec_sub[8] -= 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data02.push_back(fitter.Eval(i, par_vec_sub)); + + + par_vec_sub[1] += 10.0; + par_vec_sub[8] -= 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data03.push_back(fitter.Eval(i, par_vec_sub)); + + par_vec_sub[1] += 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data04.push_back(fitter.Eval(i, par_vec_sub)); + + + + par_vec_sub[1] += 10.0; + par_vec_sub[8] -= 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data05.push_back(fitter.Eval(i, par_vec_sub)); + + par_vec_sub[1] += 10.0; + par_vec_sub[8] -= 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data06.push_back(fitter.Eval(i, par_vec_sub)); + + par_vec_sub[1] += 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data07.push_back(fitter.Eval(i, par_vec_sub)); + + par_vec_sub[1] += 10.0; + par_vec_sub[7] = 190.0; + + for (double i(0.); i<12.0; i+=0.003) + data08.push_back(fitter.Eval(i, par_vec_sub)); + + par_vec_sub[1] += 10.0; + par_vec_sub[8] -= 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data09.push_back(fitter.Eval(i, par_vec_sub)); + + par_vec_sub[1] += 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data10.push_back(fitter.Eval(i, par_vec_sub)); + + */ + + parNo_vec.clear(); + par_vec.clear(); + par_vec_sub.clear(); + + + + + + return EXIT_SUCCESS; +}