diff --git a/src/external/TFitPofB-lib/classes/Makefile.TFitPofB b/src/external/TFitPofB-lib/classes/Makefile.TFitPofB index f2469b5c..be2e475e 100644 --- a/src/external/TFitPofB-lib/classes/Makefile.TFitPofB +++ b/src/external/TFitPofB-lib/classes/Makefile.TFitPofB @@ -24,10 +24,12 @@ LDFLAGS += # some definitions: headers (used to generate *Dict* stuff), sources, objects,... OBJS = OBJS += TTrimSPDataHandler.o +OBJS += TBulkVortexFieldCalc.o OBJS += TBofZCalc.o OBJS += TPofBCalc.o OBJS += TPofTCalc.o OBJS += TFitPofBStartupHandler.o TFitPofBStartupHandlerDict.o +OBJS += TVortex.o TVortexDict.o OBJS += TLondon1D.o TLondon1DDict.o OBJS += TSkewedGss.o TSkewedGssDict.o @@ -56,6 +58,10 @@ $(OBJS): %.o: %.cpp # Generate the ROOT CINT dictionary +TVortexDict.cpp: ../include/TVortex.h ../include/TVortexLinkDef.h + @echo "Generating dictionary $@..." + rootcint -f $@ -c -p -I$(MUSRFITINCLUDE) $^ + TLondon1DDict.cpp: ../include/TLondon1D.h ../include/TLondon1DLinkDef.h @echo "Generating dictionary $@..." rootcint -f $@ -c -p -I$(MUSRFITINCLUDE) $^ diff --git a/src/external/TFitPofB-lib/classes/TBulkVortexFieldCalc.cpp b/src/external/TFitPofB-lib/classes/TBulkVortexFieldCalc.cpp new file mode 100644 index 00000000..08613310 --- /dev/null +++ b/src/external/TFitPofB-lib/classes/TBulkVortexFieldCalc.cpp @@ -0,0 +1,101 @@ +/*************************************************************************** + + TBulkVortexFieldCalc.cpp + + Author: Bastian M. Wojek, Alexander Maisuradze + e-mail: bastian.wojek@psi.ch + + 2009/10/17 + +***************************************************************************/ + +#include "TBulkVortexFieldCalc.h" + +#include +#include + +#define PI 3.14159265358979323846 +#define TWOPI 6.28318530717958647692 + +const double fluxQuantum(2.067833667e7); // 10e14 times CGS units %% in CGS units should be 10^-7 + // in this case this is Gauss per square nm + +double getXi(const double hc2) { // get xi given Hc2 in Gauss + if (hc2) + return sqrt(fluxQuantum/(TWOPI*hc2)); + else + return 0.0; +} + +double getHc2(const double xi) { // get Hc2 given xi in nm + if (xi) + return fluxQuantum/(TWOPI*xi*xi); + else + return 0.0; +} + +TBulkTriVortexLondonFieldCalc::TBulkTriVortexLondonFieldCalc(const vector& param, const unsigned int steps, const unsigned int Ncomp) + : fNFourierComp(Ncomp) { + fSteps = steps; + fParam = param; + fB.clear(); +} + +double TBulkTriVortexLondonFieldCalc::GetBatPoint(double relX, double relY) const { + + double latConstTr(sqrt(fluxQuantum/fParam[2])*sqrt(sqrt(4.0/3.0))); + double Hc2(getHc2(fParam[1])); + + double xCoord(latConstTr*relX); // in nanometers + double yCoord(latConstTr*relY); + + if ((fParam[2] < Hc2) && (fParam[0] > fParam[1]/sqrt(2.0))) { + double KLatTr(4.0*PI/(latConstTr*sqrt(3.0))); + double fourierSum(0.0); + for(int mm = -fNFourierComp; mm <= static_cast(fNFourierComp); mm++) { + for(int nn = -fNFourierComp; nn <= static_cast(fNFourierComp); nn++) { + fourierSum += cos(KLatTr*(xCoord*mm*(0.5*sqrt(3.0)) + yCoord*mm*0.5 + yCoord*nn))*exp(-(0.5*fParam[1]*fParam[1]*KLatTr*KLatTr)* \ + (0.75*mm*mm + (nn + 0.5*mm)*(nn + 0.5*mm)))/(1.0+(fParam[0]*KLatTr*fParam[0]*KLatTr)*(0.75*mm*mm + (nn + 0.5*mm)*(nn + 0.5*mm))); + } + } + return fParam[2]*fourierSum; + } + else + return fParam[2]; + +} + +void TBulkTriVortexLondonFieldCalc::CalculateGrid() const { + + double bb(cos(PI/6.0)*0.5); + double yStep(bb/static_cast(fSteps)); + double xStep(0.5/static_cast(fSteps)); + + fB.resize(fSteps); + for (unsigned int i(0); i < fSteps; i++) + fB[i].resize(fSteps); + + for(unsigned int i(0); i < fSteps; i++) { + int j; +#pragma omp parallel for default(shared) private(j) schedule(dynamic) + for(j=0; j < static_cast(fSteps); j++) { + fB[i][j] = GetBatPoint(static_cast(i)*xStep, static_cast(j)*yStep); + } +// end pragma omp parallel + } + + return; + +} + +double TBulkTriVortexLondonFieldCalc::GetBmin() const { + return GetBatPoint(0.5, 0.5*tan(PI/6.0)); +} + +double TBulkTriVortexLondonFieldCalc::GetBmax() const { + if (!fB.empty()) + return fB[0][0]; + else + return GetBatPoint(0.0, 0.0); +} + diff --git a/src/external/TFitPofB-lib/classes/TFitPofBStartupHandler.cpp b/src/external/TFitPofB-lib/classes/TFitPofBStartupHandler.cpp index 6f77f18a..454d66ae 100644 --- a/src/external/TFitPofB-lib/classes/TFitPofBStartupHandler.cpp +++ b/src/external/TFitPofB-lib/classes/TFitPofBStartupHandler.cpp @@ -107,8 +107,12 @@ void TFitPofBStartupHandler::OnStartElement(const char *str, const TList *attrib fKey = eDeltaB; } else if (!strcmp(str, "wisdom")) { fKey = eWisdomFile; - } if (!strcmp(str, "N_theory")) { + } else if (!strcmp(str, "N_theory")) { fKey = eNSteps; + } else if (!strcmp(str, "N_VortexGrid")) { + fKey = eGridSteps; + } else if (!strcmp(str, "N_VortexFourier")) { + fKey = eVortexFourierComp; } } @@ -157,9 +161,17 @@ void TFitPofBStartupHandler::OnCharacters(const char *str) fWisdomFile = str; break; case eNSteps: - // convert str to int and assign it to the deltat-member + // convert str to int and assign it to the NSteps-member fNSteps = atoi(str); break; + case eGridSteps: + // convert str to int and assign it to the GridSteps-member + fGridSteps = atoi(str); + break; + case eVortexFourierComp: + // convert str to int and assign it to the VortexFourierComp-member + fVortexFourierComp = atoi(str); + break; default: break; } @@ -292,6 +304,20 @@ void TFitPofBStartupHandler::CheckLists() fNSteps = 3000; } + // check if any number of steps for the theory function is specified + cout << endl << "TFitPofBStartupHandler::CheckLists: check number of steps for Vortex grid ..." << endl; + if (!fGridSteps) { + cout << endl << "TFitPofBStartupHandler::CheckLists: You did not specify the number of steps for the grid. Setting the default." << endl; + fGridSteps = 200; + } + + // check if any number of steps for the theory function is specified + cout << endl << "TFitPofBStartupHandler::CheckLists: check number of steps for Vortex lattice Fourier components ..." << endl; + if (!fVortexFourierComp) { + cout << endl << "TFitPofBStartupHandler::CheckLists: You did not specify the number of Vortex lattice Fourier components. Setting the default." << endl; + fVortexFourierComp = 17; + } + } // end --------------------------------------------------------------------- diff --git a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp index d4a089ef..6b37f88a 100644 --- a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp +++ b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp @@ -16,6 +16,8 @@ #include #include +#include + /* USED FOR DEBUGGING----------------------------------- #include #include @@ -368,6 +370,81 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons } +//----------- +// Constructor that does the P(B) calculation for a bulk vortex lattice +// Parameters: dt[us], dB[G] [, Bbg[G], width[us^{-1}], weight[1] ] +//----------- + +TPofBCalc::TPofBCalc( const TBulkVortexFieldCalc &vortexLattice, const vector ¶ ) : fDT(para[0]), fDB(para[1]) +{ + + fBmin = vortexLattice.GetBmin(); + fBmax = vortexLattice.GetBmax(); + + unsigned int a1(static_cast(floor(fBmin/fDB))); + unsigned int a2(static_cast(ceil(fBmin/fDB))); + unsigned int a3(static_cast(floor(fBmax/fDB))); + unsigned int a4(static_cast(ceil(fBmax/fDB))); + + unsigned int firstZerosEnd ((a1 != a2) ? a1 : (a1 - 1)); + unsigned int lastZerosStart ((a3 != a4) ? a4 : (a4 + 1)); + + // fill the B vector + unsigned int indexBmaxFFT(static_cast(ceil(1.0/(gBar*fDT*fDB)))); + fB.resize(indexBmaxFFT, 0.0); + fPB.resize(indexBmaxFFT, 0.0); + + int i; +#pragma omp parallel for default(shared) private(i) schedule(dynamic) + for (i = 0; i < static_cast(indexBmaxFFT); i++) + fB[i] = fDB*static_cast(i); +// end pragma omp parallel + + // 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(*(fB.end()-1)+fDB); + fPB.push_back(0.0); + } + + vector< vector > vortexFields(vortexLattice.DataB()); + + if (vortexFields.empty()) { + vortexLattice.CalculateGrid(); + vortexFields = vortexLattice.DataB(); + } + + unsigned int numberOfPoints(vortexFields.size()); + + for (unsigned int j(0); j < numberOfPoints; j++){ + for (unsigned int k(0); k < numberOfPoints; k++){ + fPB[static_cast(ceil(vortexFields[j][k]/fDB))] += 1.0; + } + } + + double normalizer(static_cast(numberOfPoints*numberOfPoints)*fDB); + +#pragma omp parallel for default(shared) private(i) schedule(dynamic) + for (i = firstZerosEnd; i <= static_cast(lastZerosStart); i++) { + fPB[i] /= normalizer; + } +// end pragma omp parallel + +// // normalize p(B) +// +// double pBsum = 0.0; +// for (unsigned int i(firstZerosEnd); i +#include +#include +using namespace std; + +#include +#include "TFitPofBStartupHandler.h" + +ClassImp(TBulkTriVortexLondon) + +//------------------ +// Destructor of the TBulkTriVortexLondon class -- cleaning up +//------------------ + +TBulkTriVortexLondon::~TBulkTriVortexLondon() { + fPar.clear(); + fParForVortex.clear(); + fParForPofB.clear(); + fParForPofT.clear(); + delete fPofT; + fPofT = 0; +} + +//------------------ +// Constructor of the TBulkTriVortexLondon class +// creates (a pointer to) the TPofTCalc object (with the FFT plan) +//------------------ + +TBulkTriVortexLondon::TBulkTriVortexLondon() : fCalcNeeded(true), fFirstCall(true) { + + // read startup file + string startup_path_name("TFitPofB_startup.xml"); + + TSAXParser *saxParser = new TSAXParser(); + TFitPofBStartupHandler *startupHandler = new TFitPofBStartupHandler(); + saxParser->ConnectToHandler("TFitPofBStartupHandler", startupHandler); + int status (saxParser->ParseFile(startup_path_name.c_str())); + // check for parse errors + if (status) { // error + cout << endl << "**WARNING** reading/parsing TFitPofB_startup.xml failed." << endl; + } + + fGridSteps = startupHandler->GetGridSteps(); + fWisdom = startupHandler->GetWisdomFile(); + fVortexFourierComp = startupHandler->GetVortexFourierComp(); + + fParForVortex.clear(); + + fParForPofT.push_back(0.0); + fParForPofT.push_back(startupHandler->GetDeltat()); + fParForPofT.push_back(startupHandler->GetDeltaB()); + + fParForPofB.push_back(startupHandler->GetDeltat()); + fParForPofB.push_back(startupHandler->GetDeltaB()); + + fParForPofB.push_back(0.0); // Bkg-Field + fParForPofB.push_back(0.005); // Bkg-width + fParForPofB.push_back(0.0); // Bkg-weight + + TPofTCalc *y = new TPofTCalc(fWisdom, fParForPofT); + fPofT = y; + y = 0; + + // clean up + if (saxParser) { + delete saxParser; + saxParser = 0; + } + if (startupHandler) { + delete startupHandler; + startupHandler = 0; + } +} + +//------------------ +// TBulkTriVortexLondon-Method that calls the procedures to create B(x,y), 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 TBulkTriVortexLondon (phase, appl.field, lambda, xi) +//------------------ + +double TBulkTriVortexLondon::operator()(double t, const vector &par) const { + + assert(par.size() == 4 || par.size() == 5); + + if(t<0.0) + return cos(par[0]*0.017453293); + + // check if the function is called the first time and if yes, read in parameters + + if(fFirstCall){ + fPar = par; + + for (unsigned int i(1); i < 4; i++){ + fParForVortex.push_back(fPar[i]); + } + fFirstCall = false; + } + + // check if any parameter has changed + + bool par_changed(false); + bool only_phase_changed(false); + + for (unsigned int i(0); iDoFFT(PofB); + + }/* else { + cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl; + }*/ + + fPofT->CalcPol(fParForPofT); + + fCalcNeeded = false; + } + + return fPofT->Eval(t); + +} diff --git a/src/external/TFitPofB-lib/include/TBulkVortexFieldCalc.h b/src/external/TFitPofB-lib/include/TBulkVortexFieldCalc.h new file mode 100644 index 00000000..552e982b --- /dev/null +++ b/src/external/TFitPofB-lib/include/TBulkVortexFieldCalc.h @@ -0,0 +1,67 @@ +/*************************************************************************** + + TBulkVortexFieldCalc.h + + Author: Bastian M. Wojek, Alexander Maisuradze + e-mail: bastian.wojek@psi.ch + + 2009/10/17 + +***************************************************************************/ + +#ifndef _TBulkVortexFieldCalc_H_ +#define _TBulkVortexFieldCalc_H_ + +#include +using namespace std; + +//-------------------- +// Base class for any kind of vortex symmetry +//-------------------- + +class TBulkVortexFieldCalc { + +public: + + TBulkVortexFieldCalc() {} + + virtual ~TBulkVortexFieldCalc() { + for(unsigned int i(0); i < fB.size(); i++) + fB[i].clear(); + fB.clear(); + fParam.clear(); + } + + virtual vector< vector > DataB() const {return fB;} + virtual void CalculateGrid() const = 0; + virtual double GetBatPoint(double, double) const = 0; + virtual double GetBmin() const = 0; + virtual double GetBmax() const = 0; + +protected: + unsigned int fSteps; // number of steps, the "unit cell" of the vortex lattice is devided in (in x- and y- direction) + vector fParam; // lambda, xi, appl.field + mutable vector< vector > fB; // fields in a "unit cell" of the vortex lattice +}; + +//-------------------- +// Class for triangular vortex lattice, London model +//-------------------- + +class TBulkTriVortexLondonFieldCalc : public TBulkVortexFieldCalc { + +public: + + TBulkTriVortexLondonFieldCalc(const vector&, const unsigned int steps = 200, const unsigned int Ncomp = 17); + ~TBulkTriVortexLondonFieldCalc() {} + + void CalculateGrid() const; + double GetBatPoint(double, double) const; + double GetBmin() const; + double GetBmax() const; + +private: + unsigned int fNFourierComp; +}; + +#endif // _TBulkVortexFieldCalc_H_ diff --git a/src/external/TFitPofB-lib/include/TFitPofBStartupHandler.h b/src/external/TFitPofB-lib/include/TFitPofBStartupHandler.h index bdd6e373..e3a8ce6e 100644 --- a/src/external/TFitPofB-lib/include/TFitPofBStartupHandler.h +++ b/src/external/TFitPofB-lib/include/TFitPofBStartupHandler.h @@ -65,9 +65,11 @@ class TFitPofBStartupHandler : public TQObject { virtual const double GetDeltaB() const { return fDeltaB; } virtual const string GetWisdomFile() const { return fWisdomFile; } virtual const unsigned int GetNSteps() const { return fNSteps; } + virtual const unsigned int GetGridSteps() const { return fGridSteps; } + virtual const unsigned int GetVortexFourierComp() const { return fVortexFourierComp; } private: - enum EKeyWords {eEmpty, eComment, eDataPath, eEnergy, eEnergyList, eDeltat, eDeltaB, eWisdomFile, eNSteps}; + enum EKeyWords {eEmpty, eComment, eDataPath, eEnergy, eEnergyList, eDeltat, eDeltaB, eWisdomFile, eNSteps, eGridSteps, eVortexFourierComp}; EKeyWords fKey; @@ -77,6 +79,8 @@ class TFitPofBStartupHandler : public TQObject { double fDeltaB; string fWisdomFile; unsigned int fNSteps; + unsigned int fGridSteps; + unsigned int fVortexFourierComp; ClassDef(TFitPofBStartupHandler, 1) }; diff --git a/src/external/TFitPofB-lib/include/TPofBCalc.h b/src/external/TFitPofB-lib/include/TPofBCalc.h index 101b5ccf..ecd0ee0a 100644 --- a/src/external/TFitPofB-lib/include/TPofBCalc.h +++ b/src/external/TFitPofB-lib/include/TPofBCalc.h @@ -14,6 +14,7 @@ #include "TBofZCalc.h" #include "TTrimSPDataHandler.h" +#include "TBulkVortexFieldCalc.h" #define gBar 0.0135538817 #define pi 3.14159265358979323846 @@ -25,6 +26,7 @@ public: TPofBCalc( const string&, const vector& ); TPofBCalc( const TBofZCalc&, const TTrimSPData&, const vector&, unsigned int ); TPofBCalc( const TBofZCalcInverse&, const TTrimSPData&, const vector& ); + TPofBCalc( const TBulkVortexFieldCalc&, const vector& ); TPofBCalc( const vector&, const vector& , double dt = 0.01 ); ~TPofBCalc() { fB.clear(); diff --git a/src/external/TFitPofB-lib/include/TVortex.h b/src/external/TFitPofB-lib/include/TVortex.h new file mode 100644 index 00000000..46d7312e --- /dev/null +++ b/src/external/TFitPofB-lib/include/TVortex.h @@ -0,0 +1,41 @@ +/*************************************************************************** + + TVortex.h + + Author: Bastian M. Wojek + e-mail: bastian.wojek@psi.ch + + 2009/10/17 + +***************************************************************************/ + +#ifndef _TVortex_H_ +#define _TVortex_H_ + +#include "PUserFcnBase.h" +#include "TPofTCalc.h" + +class TBulkTriVortexLondon : public PUserFcnBase { + +public: + TBulkTriVortexLondon(); + ~TBulkTriVortexLondon(); + + double operator()(double, const vector&) const; + +private: + mutable vector fPar; + TPofTCalc *fPofT; + mutable bool fCalcNeeded; + mutable bool fFirstCall; + mutable vector fParForPofT; + mutable vector fParForVortex; + mutable vector fParForPofB; + string fWisdom; + unsigned int fGridSteps; + unsigned int fVortexFourierComp; + + ClassDef(TBulkTriVortexLondon,1) +}; + +#endif //_TLondon1D_H_ diff --git a/src/external/TFitPofB-lib/include/TVortexLinkDef.h b/src/external/TFitPofB-lib/include/TVortexLinkDef.h new file mode 100644 index 00000000..8e9fdfc4 --- /dev/null +++ b/src/external/TFitPofB-lib/include/TVortexLinkDef.h @@ -0,0 +1,23 @@ +/*************************************************************************** + + TVortexLinkDef.h + + Author: Bastian M. Wojek + e-mail: bastian.wojek@psi.ch + + 2009/10/17 + +***************************************************************************/ + +// root dictionary stuff -------------------------------------------------- +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class TBulkTriVortexLondon+; + +#endif //__CINT__ +// root dictionary stuff -------------------------------------------------- + diff --git a/src/external/TFitPofB-lib/test/Makefile.testVortex b/src/external/TFitPofB-lib/test/Makefile.testVortex new file mode 100644 index 00000000..e7195be1 --- /dev/null +++ b/src/external/TFitPofB-lib/test/Makefile.testVortex @@ -0,0 +1,49 @@ +#--------------------------------------------------- +# get compilation and library flags from root-config + +ROOTCFLAGS = $(shell $(ROOTSYS)/bin/root-config --cflags) +ROOTLIBS = $(shell $(ROOTSYS)/bin/root-config --libs) + +#--------------------------------------------------- + +CXX = g++-4.4.0 +CXXFLAGS = -O3 -Wall +LOCALINCLUDE = ../include +ROOTINCLUDE = $(ROOTSYS)/include +INCLUDES = -I$(LOCALINCLUDE) -I$(ROOTINCLUDE) +LD = g++-4.4.0 +LDFLAGS = -O3 -L../classes -lTFitPofB -lfftw3 -lm -fopenmp -lPMusr + +# the output from the root-config script: +CXXFLAGS += $(ROOTCFLAGS) +LDFLAGS += + +# the ROOT libraries +LIBS = $(ROOTLIBS) -lXMLParser -lMathMore + +EXEC = testVortex + +# some definitions: headers, sources, objects,... +OBJS = +OBJS += $(EXEC).o + +# make the executable: +# +all: $(EXEC) + +$(EXEC): $(OBJS) + @echo "---> Building $(EXEC) ..." + $(LD) $(LDFLAGS) $(OBJS) -o $(EXEC) $(LIBS) + @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/TFitPofB_startup.xml b/src/external/TFitPofB-lib/test/TFitPofB_startup.xml index dfcf3a2d..cdc2e709 100644 --- a/src/external/TFitPofB-lib/test/TFitPofB_startup.xml +++ b/src/external/TFitPofB-lib/test/TFitPofB_startup.xml @@ -5,12 +5,16 @@ Defines path/prefix and energies (keV, format: %02u_%1u) of TrimSP-rge-files, path/name to the FFTW-wisdom-file and time/field binning (us/G) N_theory determines the number of points in "real space" where the theory function will be calculated + N_VortexGrid determines the number of points in which each axis of the "real space" will be divided when calculating the vortex B(x,y) + N_VortexFourier defines the number of Fourier components which should be summed up in order to get B(x,y) /home/l_wojek/TrimSP/YBCOxtal/YBCOxtal-500000- /home/l_wojek/analysis/WordsOfWisdom.dat 0.01 0.1 5000 + 200 + 17 03_0 03_6 diff --git a/src/external/TFitPofB-lib/test/testVortex.cpp b/src/external/TFitPofB-lib/test/testVortex.cpp new file mode 100644 index 00000000..d109e1e6 --- /dev/null +++ b/src/external/TFitPofB-lib/test/testVortex.cpp @@ -0,0 +1,565 @@ +#include "TPofTCalc.h" +#include +#include + +using namespace std; + +#include +#include + +int main(){ + + vector parForVortex; + parForVortex.push_back(150.0); //lambda + parForVortex.push_back(3.0); //xi + parForVortex.push_back(500.0); //app.field + + vector parForPofB; + parForPofB.push_back(0.01); //dt + parForPofB.push_back(1.0); //dB + + TBulkTriVortexLondonFieldCalc *vortexLattice = new TBulkTriVortexLondonFieldCalc(parForVortex, 250, 17); + TPofBCalc *PofB = new TPofBCalc(*vortexLattice, parForPofB); + + vector b(PofB->DataB()); + vector pb(PofB->DataPB()); + + double test(0.0); + + ofstream ofx("testVector.dat"); + for (unsigned int i(0); i par_vec(par_arr, par_arr+(sizeof(par_arr)/sizeof(par_arr[0]))); + + + vector parForBofZ; + for (unsigned int i(2); i parForPofB; + parForPofB.push_back(0.01); //dt + parForPofB.push_back(0.1); //dB + parForPofB.push_back(par_vec[1]); + parForPofB.push_back(par_vec[2]); // Bkg-Field + parForPofB.push_back(0.005); // Bkg-width (in principle zero) + vector interfaces; + interfaces.push_back(par_vec[3]); // dead layer + parForPofB.push_back(calcData.LayerFraction(par_vec[1], 1, interfaces)); // Fraction of muons in the deadlayer + + + vector zonk; + vector donk; + vector hurgaB; + vector hurgaPB; + + for (unsigned int u(0); u<10 ; u++) { + + parForBofZ[par_vec.size()-3] = par_vec[4] + double(u)*5.; + + TLondon1D_HS *BofZ = new TLondon1D_HS(parForBofZ); + TPofBCalc *PofB = new TPofBCalc(*BofZ, calcData, parForPofB); + + if(u == 0){ + hurgaB = PofB->DataB(); + hurgaPB = PofB->DataPB(); + PofB->ConvolveGss(7.0); + donk = PofB->DataPB(); + } + + delete BofZ; + delete PofB; + } + +// TPofBCalc sumPB(hurgaB, zonk); + +// sumPB.AddBackground(par_vec[2], 0.15, 0.1); +// sumPB.ConvolveGss(7.0); +// donk = sumPB.DataPB(); + + char debugfile[50]; + int n = sprintf (debugfile, "testInverseExp-BpB_B-%.4f_l-%.3f_E-%.1f_normal.dat", par_vec[2], par_vec[4], par_vec[1]); + + if (n > 0) { + ofstream of7(debugfile); + for (unsigned int i(0); i parForPofT; +// parForPofT.push_back(par_vec[0]); //phase +// parForPofT.push_back(0.01); //dt +// parForPofT.push_back(0.02); //dB + +// TLondon1D_HS BofZ(parForBofZ); + +// vector zz; +// vector Bz; +// for(double i(0.5); i<2001; i+=1.0) { +// zz.push_back(i/10.0); +// if(!(i/10.0 < parForBofZ[1])) { +// Bz.push_back(BofZ.GetBofZ(i/10.0)); +// } +// else { +// Bz.push_back(parForBofZ[0]); +// } +// } +/* + char debugfile1[50]; + int nn = sprintf (debugfile1, "testInverseExp-Bz_z-%.4f_l-%.3f_E-%.1f_normal.dat", par_vec[2], par_vec[9], par_vec[1]); + if (nn > 0) { + ofstream of01(debugfile1); + for (unsigned int i(0); i<2000; i++) { + of01 << zz[i] << " " << Bz[i] << endl; + } + of01.close(); + } + */ +// TPofBCalc PofB(BofZ, calcData, parForPofB); + + +// double t1(0.0); +// get start time +// struct timeval tv_start, tv_stop; +// gettimeofday(&tv_start, 0); +// +// TPofBCalc PofB(BofZ, calcData, parForPofB); +// +// gettimeofday(&tv_stop, 0); +// t1 = (tv_stop.tv_sec - tv_start.tv_sec)*1000.0 + (tv_stop.tv_usec - tv_start.tv_usec)/1000.0; + +// cout << "p(B) calculation time with derivatives of the inverse function: " << t1 << " ms" << endl; + +// gettimeofday(&tv_start, 0); +// +// BofZ.Calculate(); +// TPofBCalc PofB1(BofZ, calcData, parForPofB, 1); +// +// gettimeofday(&tv_stop, 0); +// t1 = (tv_stop.tv_sec - tv_start.tv_sec)*1000.0 + (tv_stop.tv_usec - tv_start.tv_usec)/1000.0; +// +// cout << "p(B) calculation time without derivatives of the inverse function: " << t1 << " ms" << endl; + +// PofB.ConvolveGss(1.17); + +// PofB.AddBackground(par_vec[2], 0.2, calcData.LayerFraction(E, 4, interfaces)); +// PofB.ConvolveGss(1.17); + +// vector hurgaB(PofB.DataB()); +// vector hurgaPB(PofB.DataPB()); +// vector hurgaPB1(PofB1.DataPB()); +// +// char debugfile[50]; +// int n = sprintf (debugfile, "testInverseExp-BpB_B-%.4f_l-%.3f_E-%.1f_normal.dat", par_vec[2], par_vec[4], par_vec[1]); +// +// if (n > 0) { +// ofstream of7(debugfile); +// for (unsigned int i(0); i 0) { +// ofstream of8(debugfilex); +// for (double i(0.); i<12.0; i+=0.003) { +// of8 << i << " " << poft.Eval(i) << endl; +// } +// of8.close(); +// } +/* + PofB.ConvolveGss(8.8); + + vector hurgaB1(PofB.DataB()); + vector hurgaPB1(PofB.DataPB()); + + n = sprintf (debugfile, "BpB_B-%.4f_l-%.3f_E-%.1f_broadend8.8G.dat", par_vec[2], par_vec[4], par_vec[1]); + + if (n > 0) { + ofstream of1(debugfile); + for (unsigned int i(0); i z2(calcData.DataZ(par_vec[1])); + vector nz2(calcData.DataNZ(par_vec[1])); + + ofstream of2("Implantation-profile-broad.dat"); + for (unsigned int i(0); i hurgaB2(PofB2.DataB()); + vector hurgaPB2(PofB2.DataPB()); + + n = sprintf (debugfile, "BpB_B-%.4f_l-%.3f_E-%.1f_broadend10nm.dat", par_vec[2], par_vec[4], par_vec[1]); + + if (n>0) { + ofstream of8(debugfile); + for (unsigned int i(0); i hurgaB3(PofB2.DataB()); + vector hurgaPB3(PofB2.DataPB()); + + n = sprintf (debugfile, "BpB_B-%.4f_l-%.3f_E-%.1f_broadend10nm+7.5G.dat", par_vec[2], par_vec[4], par_vec[1]); + + if (n > 0) { + ofstream of9(debugfile); + 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 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 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 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 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[0] += 10.0; + par_vec_sub[10] -= 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data02.push_back(fitter.Eval(i, par_vec_sub)); + + + par_vec_sub[0] += 10.0; + par_vec_sub[10] -= 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data03.push_back(fitter.Eval(i, par_vec_sub)); + + par_vec_sub[0] += 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data04.push_back(fitter.Eval(i, par_vec_sub)); + + + + par_vec_sub[0] += 10.0; + par_vec_sub[10] -= 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data05.push_back(fitter.Eval(i, par_vec_sub)); + + par_vec_sub[0] += 10.0; + par_vec_sub[10] -= 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data06.push_back(fitter.Eval(i, par_vec_sub)); + + par_vec_sub[0] += 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data07.push_back(fitter.Eval(i, par_vec_sub)); + + par_vec_sub[0] += 10.0; + par_vec_sub[10] = 190.0; + + for (double i(0.); i<12.0; i+=0.003) + data08.push_back(fitter.Eval(i, par_vec_sub)); + + par_vec_sub[0] += 10.0; + par_vec_sub[10] -= 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data09.push_back(fitter.Eval(i, par_vec_sub)); + + par_vec_sub[0] += 10.0; + + for (double i(0.); i<12.0; i+=0.003) + data10.push_back(fitter.Eval(i, par_vec_sub)); + + */ + + +// par_vec.clear(); + + + + + + return 0; +} diff --git a/src/external/libLFRelaxation/Makefile.libLFRelaxation b/src/external/libLFRelaxation/Makefile.libLFRelaxation index 21563e4d..2c74da52 100644 --- a/src/external/libLFRelaxation/Makefile.libLFRelaxation +++ b/src/external/libLFRelaxation/Makefile.libLFRelaxation @@ -1,63 +1,169 @@ #--------------------------------------------------- -# get compilation flags from root-config +# Makefile.libGapIntegrals +# +# Author: Bastian M. Wojek +# e-mail: bastian.wojek@psi.ch +# +# $Id$ +# +#--------------------------------------------------- + +#--------------------------------------------------- +# get compilation and library flags from root-config ROOTCFLAGS = $(shell $(ROOTSYS)/bin/root-config --cflags) -# ROOTLIBS = $(shell $(ROOTSYS)/bin/root-config --libs) +ROOTLIBS = $(shell $(ROOTSYS)/bin/root-config --libs) +ROOTGLIBS = $(shell $(ROOTSYS)/bin/root-config --glibs) #--------------------------------------------------- +# depending on the architecture, choose the compiler, +# linker, and the flags to use +# -OS = LINUX +ARCH = $(shell $(ROOTSYS)/bin/root-config --arch) + +ifeq ($(ARCH),linux) +OS = LINUX +endif +ifeq ($(ARCH),linuxx8664gcc) +OS = LINUX +endif +ifeq ($(ARCH),win32gcc) +OS = WIN32GCC +endif +ifeq ($(ARCH),macosx) +OS = DARWIN +endif + +# -- Linux +ifeq ($(OS),LINUX) CXX = g++-4.4.0 CXXFLAGS = -O3 -Wall -Wno-trigraphs -fPIC -MUSRFITINCLUDE = ../../include -#MUSRFITINCLUDE = /home/l_wojek/rep/analysis/musrfit/src/include -LOCALINCLUDE = . -ROOTINCLUDE = $(ROOTSYS)/include -INCLUDES = -I$(LOCALINCLUDE) -I$(MUSRFITINCLUDE) -I$(ROOTINCLUDE) +PMUSRPATH = ../../include +MNPATH = $(ROOTSYS)/include +GSLPATH = /usr/include/gsl +CUBAPATH = /usr/local/include +INCLUDES = -I$(PMUSRPATH) -I$(MNPATH) -I$(GSLPATH) -I$(CUBAPATH) LD = g++-4.4.0 -LDFLAGS = -O3 +LDFLAGS = -O SOFLAGS = -shared +SHLIB = libLFRelaxation.so +endif + +# -- Windows/Cygwin +ifeq ($(OS),WIN32GCC) +CXX = g++ +CXXFLAGS = -O3 -Wall -Wno-trigraphs -D_DLL +PMUSRPATH = ../../include +MNPATH = $(ROOTSYS)/include +GSLPATH = /usr/include/gsl +CUBAPATH = /usr/local/include +INCLUDES = -I$(PMUSRPATH) -I$(MNPATH) -I$(GSLPATH) -I$(CUBAPATH) +LD = g++ +LDFLAGS = -O -Wl,--enable-auto-import -Wl,--enable-runtime-pseudo-reloc +SOFLAGS = -shared -Wl,--export-all-symbols +SHLIB = libLFRelaxation.dll +endif + +# -- MacOSX/Darwin +ifeq ($(OS),DARWIN) +CXX = g++ +CXXFLAGS = -O3 -Wall -Wno-trigraphs -fPIC +PMUSRPATH = ../../include +MNPATH = $(ROOTSYS)/include +FINKPATH = /sw/include +GSLPATH = $(FINKPATH)/gsl +CUBAPATH = $(FINKPATH) +INCLUDES = -I$(PMUSRPATH) -I$(MNPATH) -I$(GSLPATH) -I$(CUBAPATH) +LD = g++ +LDFLAGS = -O -Xlinker -bind_at_load +SOFLAGS = -dynamiclib -flat_namespace -undefined suppress -Wl,-x +SHLIB = libLFRelaxation.dylib +endif # the output from the root-config script: CXXFLAGS += $(ROOTCFLAGS) LDFLAGS += +# the ROOT libraries (G = graphic) +LIBS = $(ROOTLIBS) +GLIBS = $(ROOTGLIBS) + +# GSL lib +GSLLIB = -lgslcblas -lgsl +# Cuba lib +CUBALIB = -L/usr/local/lib -lcuba -lm +# PMusr lib +PMUSRLIB = -L$(ROOTSYS)/lib -lPMusr + +ifeq ($(OS),WIN32GCC) +# GSL lib +GSLLIB = -L/usr/lib -lgslcblas -lgsl +# Cuba lib +CUBALIB = -L/usr/local/lib -lcuba -lm +# PMusr lib +PMUSRLIB = -L$(ROOTSYS)/lib -lPMusr -lMathMore +endif + +ifeq ($(OS),DARWIN) +# GSL lib +GSLLIB = -L/sw/lib -lgslcblas -lgsl +# Cuba lib +CUBALIB = -L/sw/lib -lcuba -lm +# PMusr lib +PMUSRLIB = -L$(ROOTSYS)/lib -lPMusr -lMathMore +endif + # some definitions: headers (used to generate *Dict* stuff), sources, objects,... OBJS = OBJS += TLFRelaxation.o TLFRelaxationDict.o -SHLIB = libLFRelaxation.so +EXTOBJS = BMWIntegrator.o + +# make the shared libs: -# make the shared lib: -# all: $(SHLIB) -$(SHLIB): $(OBJS) +$(SHLIB): $(EXTOBJS) $(OBJS) @echo "---> Building shared library $(SHLIB) ..." /bin/rm -f $(SHLIB) - $(LD) $(LDFLAGS) $(OBJS) $(SOFLAGS) -o $(SHLIB) + $(LD) $(SOFLAGS) $(LDFLAGS) -o $(SHLIB) $(EXTOBJS) $(OBJS) $(LIBS) $(PMUSRLIB) $(GSLLIB) $(CUBALIB) @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) $(SHLIB) *Dict* core* *~ - @echo "---> removing $(OBJS) $(SHLIB)" +clean:; @rm -f $(OBJS) $(EXTOBJS) $(SHLIB) *Dict* core* *~ + @echo "---> removing $(OBJS) $(EXTOBJS) $(SHLIB)" # +BMWIntegrator.o: ../BMWIntegrator/BMWIntegrator.cpp + $(CXX) $(INCLUDES) $(CXXFLAGS) -c $< + $(OBJS): %.o: %.cpp $(CXX) $(INCLUDES) $(CXXFLAGS) -c $< # Generate the ROOT CINT dictionary -TLFRelaxationDict.cpp: ./TLFRelaxation.h ./TLFRelaxationLinkDef.h +TLFRelaxationDict.cpp: ../BMWIntegrator/BMWIntegrator.h ./TLFRelaxation.h ./TLFRelaxationLinkDef.h @echo "Generating dictionary $@..." - rootcint -f $@ -c -p -I$(MUSRFITINCLUDE) $^ + rootcint -f $@ -c -p $^ install: all - @echo "Installing shared lib: libLFRelaxation.so" + @echo "Installing shared lib: $(SHLIB)" ifeq ($(OS),LINUX) cp -pv $(SHLIB) $(ROOTSYS)/lib - cp -pv $(LOCALINCLUDE)/*.h $(ROOTSYS)/include + cp -pv ../BMWIntegrator/BMWIntegrator.h TLFRelaxation.h $(ROOTSYS)/include endif +ifeq ($(OS),WIN32GCC) + cp -pv $(SHLIB) $(ROOTSYS)/bin + ln -sf $(ROOTSYS)/bin/$(SHLIB) $(ROOTSYS)/lib/$(SHLIB) + cp -pv ../BMWIntegrator/BMWIntegrator.h TLFRelaxation.h $(ROOTSYS)/include +endif +ifeq ($(OS),DARWIN) + cp -pv $(SHLIB) $(ROOTSYS)/lib + cp -pv ../BMWIntegrator/BMWIntegrator.h TLFRelaxation.h $(ROOTSYS)/include +endif + +cleaninstall: clean install diff --git a/src/external/libLFRelaxation/TLFRelaxation.cpp b/src/external/libLFRelaxation/TLFRelaxation.cpp index 9fde97ce..abd1a3f4 100644 --- a/src/external/libLFRelaxation/TLFRelaxation.cpp +++ b/src/external/libLFRelaxation/TLFRelaxation.cpp @@ -9,11 +9,31 @@ ***************************************************************************/ +/*************************************************************************** + * Copyright (C) 2009 by Bastian M. Wojek * + * bastian.wojek@psi.ch * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + #include #include #include -#include "TIntegrator.h" +#include "../BMWIntegrator/BMWIntegrator.h" #include "TLFRelaxation.h" using namespace std; @@ -26,6 +46,7 @@ ClassImp(TLFStatGssKT) ClassImp(TLFStatLorKT) ClassImp(TLFDynGssKT) ClassImp(TLFDynLorKT) +ClassImp(TLFSGInterpolation) // LF Static Gaussian KT @@ -96,7 +117,17 @@ double TLFStatLorKT::operator()(double t, const vector &par) const { double coeff2(coeff1*coeff1); double coeff3((1.0+coeff2)*par[1]); - return (1.0-(coeff1*TMath::Exp(-par[1]*t)*TMath::BesselJ1(omegaL*t))-(coeff2*(TMath::BesselJ0(omegaL*t)*TMath::Exp(-par[1]*t)-1.0))-coeff3*fIntBesselJ0Exp->IntegrateFunc(0.,t)); + double w0t(omegaL*t); + double j1, j0; + if (fabs(w0t) < 0.001) { // check zero time limits of the spherical bessel functions j0(x) and j1(x) + j0 = 1.0; + j1 = 0.0; + } else { + j0 = TMath::Sin(w0t)/w0t; + j1 = (TMath::Sin(w0t)-w0t*TMath::Cos(w0t))/(w0t*w0t); + } + + return (1.0-(coeff1*TMath::Exp(-par[1]*t)*j1)-(coeff2*(j0*TMath::Exp(-par[1]*t)-1.0))-coeff3*fIntBesselJ0Exp->IntegrateFunc(0.,t)); } // LF Dynamic Gaussian KT @@ -188,7 +219,7 @@ double TLFDynGssKT::operator()(double t, const vector &par) const { double sigsqtsq(sigsq*t*t); return 0.33333333333333333333+0.66666666666666666667*(1.0-sigsqtsq)*exp(-0.5*sigsqtsq); // ZF static Gauss KT } - return exp(-2.0*sigsq/nusq*(enut-1.0+nut)); // ZF Abragam Delta^2->2*Delta^2 +// return exp(-2.0*sigsq/nusq*(enut-1.0+nut)); // ZF Abragam Delta^2->2*Delta^2 (only here for TESTING - DO NOT return this value in an production environment!!!!!!) } if((par[2] >= 5.0*par[1]) || (omegaL >= 10.0*par[1])) { @@ -272,7 +303,7 @@ double TLFDynGssKT::operator()(double t, const vector &par) const { // LF Dynamic Lorentz KT -TLFDynLorKT::TLFDynLorKT() : fCalcNeeded(true), fFirstCall(true), fWisdom("/home/l_wojek/analysis/WordsOfWisdom.dat"), fNSteps(524288), fDt(0.000040), fCounter(0), fA(1.0), fL1(0.0), fL2(0.0) { +TLFDynLorKT::TLFDynLorKT() : fCalcNeeded(true), fFirstCall(true), fWisdom("/home/l_wojek/analysis/WordsOfWisdom.dat"), fNSteps(524288), fDt(0.000040), fCounter(0), fL1(0.0), fL2(0.0) { // Calculate d_omega and C for given NFFT and dt fDw = TMath::Pi()/fNSteps/fDt; fC = 2.0*TMath::Log(double(fNSteps))/(double(fNSteps-1)*fDt); @@ -283,14 +314,14 @@ TLFDynLorKT::TLFDynLorKT() : fCalcNeeded(true), fFirstCall(true), fWisdom("/home FILE *wordsOfWisdomR; wordsOfWisdomR = fopen(fWisdom.c_str(), "r"); if (wordsOfWisdomR == NULL) { - cout << "TLFDynLorKT::TLFDynGssKT: Couldn't open wisdom file ..." << endl; + cout << "TLFDynLorKT::TLFDynLorKT: Couldn't open wisdom file ..." << endl; } else { wisdomLoaded = fftw_import_wisdom_from_file(wordsOfWisdomR); fclose(wordsOfWisdomR); } if (!wisdomLoaded) { - cout << "TLFDynLorKT::TLFDynGssKT: No wisdom is imported..." << endl; + cout << "TLFDynLorKT::TLFDynLorKT: No wisdom is imported..." << endl; } // END of WisdomLoading @@ -307,7 +338,7 @@ TLFDynLorKT::~TLFDynLorKT() { FILE *wordsOfWisdomW; wordsOfWisdomW = fopen(fWisdom.c_str(), "w"); if (wordsOfWisdomW == NULL) { - cout << "TLFDynGssKT::~TLFDynGssKT: Could not open file ... No wisdom is exported..." << endl; + cout << "TLFDynLorKT::~TLFDynLorKT: Could not open file ... No wisdom is exported..." << endl; } else { fftw_export_wisdom_to_file(wordsOfWisdomW); fclose(wordsOfWisdomW); @@ -322,10 +353,6 @@ TLFDynLorKT::~TLFDynLorKT() { cout << "TLFDynLorKT::~TLFDynLorKT(): " << fCounter << " full FFT cyles needed..." << endl; } -const double TLFDynLorKT::fX[16]={1.61849, 1.27059, 0.890315, 6.72608, 0.327258, 0.981975, 1.5964, 2.31023, 2.37689, 5.63945, 0.235734, 1.10617, 1.1664, 0.218402, 0.266562, 0.471331}; - -const double TLFDynLorKT::fY[20]={1.54699, 0.990321, 0.324923, 0.928399, 2.11155, 0.615106, 1.49907, 0.679423, 6.17621, 1.38907, 0.979608, 0.0593631, 0.806427, 3.33144, 1.05059, 1.29922, 0.254661, 0.284348, 0.991419, 0.375213}; - double TLFDynLorKT::operator()(double t, const vector &par) const { assert(par.size() == 3); // three parameters nuL=gbar*B,sigma,fluct.rate nu @@ -349,6 +376,8 @@ double TLFDynLorKT::operator()(double t, const vector &par) const { } double omegaL(TWOPI*par[0]); // Larmor frequency + double a(par[1]); // static width + double nu(par[2]); // hopping rate if(fabs(par[0])<0.00135538817){ // Zero Field if(!par[2]){ // nu=0 @@ -357,37 +386,43 @@ double TLFDynLorKT::operator()(double t, const vector &par) const { } } - // semi-analytical approximations for {nu >= 4a} and {omegaL >= 10a} + // semi-analytical approximation for {nu >= 5a} and {omegaL >= 30a} - if(par[2]){ // nu!=0 - if(par[2] >= 4.0*par[1]){ // nu >= 4a - if(fCalcNeeded){ - double wLnu(omegaL/par[2]); // omegaL/nu - double wLnusq(wLnu*wLnu); // (omegaL/nu)^2 - double anu(par[1]/par[2]); // a/nu + // check if hopping > 5 * damping, of Larmor angular frequency is > 30 * damping (BMW limit) + if((par[2] > 5.0*par[1]) || (omegaL > 30.0*par[1])){ + if(fCalcNeeded){ - fA=1.0+wLnu*anu*(-fX[0]/(wLnusq+fX[1])+fX[2]/(wLnusq*wLnu+fX[3])+fX[4]/(wLnusq*wLnusq+fX[5])); - fL1=par[1]*(fX[6]/(wLnu+fX[7])+fX[8]/(wLnusq+fX[9])+fX[10]/(wLnusq*wLnusq+fX[11])); - fL2=par[2]*(fX[12]*tanh(fX[13]*wLnu)+fX[14]*wLnu+fX[15]); - - fCalcNeeded=false; + // 'c' and 'd' are parameters BMW obtained by fitting large parameter space LF-curves to the model below + const double c[7] = {1.15331, 1.64826, -0.71763, 3.0, 0.386683, -5.01876, 2.41854}; + const double d[4] = {2.44056, 2.92063, 1.69581, 0.667277}; + double w0N[4]; + double nuN[4]; + w0N[0] = omegaL; + nuN[0] = nu; + for (unsigned int i=1; i<4; i++) { + w0N[i] = omegaL * w0N[i-1]; + nuN[i] = nu * nuN[i-1]; } - return fA*exp(-fL1*t)+(1.0-fA)*exp(-fL2*t)*cos(omegaL*t); - } - if(omegaL >= 10.0*par[1]){ // omegaL >= 10a - if(fCalcNeeded){ - double wLnu(omegaL/par[2]); // omegaL/nu - double wLnusq(wLnu*wLnu); // (omegaL/nu)^2 - double anu(par[1]/par[2]); // a/nu + double denom(w0N[3]+d[0]*w0N[2]*nuN[0]+d[1]*w0N[1]*nuN[1]+d[2]*w0N[0]*nuN[2]+d[3]*nuN[3]); + fL1 = (c[0]*w0N[2]+c[1]*w0N[1]*nuN[0]+c[2]*w0N[0]*nuN[1])/denom; + fL2 = (c[3]*w0N[2]+c[4]*w0N[1]*nuN[0]+c[5]*w0N[0]*nuN[1]+c[6]*nuN[2])/denom; - fA=1.0-fY[0]*pow(anu,fY[1])/(wLnu+fY[2])+fY[3]*pow(anu,fY[4])/(wLnusq-fY[5])+fY[6]*pow(anu,fY[7])/(wLnusq*wLnu+fY[8]); - fL1=par[1]*(fY[9]/(pow(wLnu,fY[10])+fY[11])-fY[12]/(pow(wLnu,fY[13])+fY[14])); - fL2=par[2]*(fY[15]*tanh(fY[16]*wLnu)+fY[17]*pow(wLnu,fY[18])+fY[19]); - - fCalcNeeded=false; - } - return fA*exp(-fL1*t)+(1.0-fA)*exp(-fL2*t)*cos(omegaL*t); + fCalcNeeded=false; } + + double w0t(omegaL*t); + double j1, j0; + if (fabs(w0t) < 0.001) { // check zero time limits of the spherical bessel functions j0(x) and j1(x) + j0 = 1.0; + j1 = 0.0; + } else { + j0 = TMath::Sin(w0t)/w0t; + j1 = (TMath::Sin(w0t)-w0t*TMath::Cos(w0t))/(w0t*w0t); + } + + double Gamma_t = -4.0/3.0*a*(fL1*(1.0-j0*TMath::Exp(-nu*t))+fL2*j1*TMath::Exp(-nu*t)+(1.0-fL2*omegaL/3.0-fL1*nu)*t); + + return TMath::Exp(Gamma_t); } // if no approximation can be used -> Laplace transform @@ -411,15 +446,17 @@ double TLFDynLorKT::operator()(double t, const vector &par) const { for(unsigned int i(1); i &par) const { // return fFFTtime[int(t/fDt)]; return fDw*TMath::Exp(fC*t)/TMath::Pi()*fFFTtime[int(t/fDt)]; } + + +// De Renzi Spin Glass Interpolation + +TLFSGInterpolation::TLFSGInterpolation() { + TIntSGInterpolation *integral = new TIntSGInterpolation(); + fIntegral = integral; + integral = 0; +} + +TLFSGInterpolation::~TLFSGInterpolation() { + delete fIntegral; + fIntegral = 0; +} + +double TLFSGInterpolation::operator()(double t, const vector &par) const { + + assert(par.size() == 4); // two parameters nu=gbar*B [MHz], a [1/us], lambda [1/us], beta [1] + + if((t <= 0.0) || (par[3] <= 0.0)) + return 1.0; + + double expo(0.5*par[1]*par[1]*t*t/par[3]+par[2]*t); + + if(TMath::Abs(par[0])<0.00135538817){ + return (0.33333333333333333333*TMath::Exp(-TMath::Power(par[2]*t,par[3])) + \ + 0.66666666666666666667*(1.0-par[1]*par[1]*t*t/TMath::Power(expo,(1.0-par[3])))*TMath::Exp(-TMath::Power(expo,par[3]))); + } + + vector intpar(par); + intpar.push_back(t); + fIntegral->SetParameters(intpar); + + double omegaL(TMath::TwoPi()*par[0]); + + return (TMath::Exp(-TMath::Power(par[2]*t,par[3])) + 2.0*par[1]*par[1]/(omegaL*omegaL) * \ + ((TMath::Cos(omegaL*t)-TMath::Sin(omegaL*t)/(omegaL*t))*TMath::Exp(-TMath::Power(expo,par[3]))/TMath::Power(expo,(1.0-par[3])) + \ + fIntegral->IntegrateFunc(0.,t))); +} diff --git a/src/external/libLFRelaxation/TLFRelaxation.h b/src/external/libLFRelaxation/TLFRelaxation.h index f541ee10..c5324b74 100644 --- a/src/external/libLFRelaxation/TLFRelaxation.h +++ b/src/external/libLFRelaxation/TLFRelaxation.h @@ -9,6 +9,26 @@ ***************************************************************************/ +/*************************************************************************** + * Copyright (C) 2009 by Bastian M. Wojek * + * bastian.wojek@psi.ch * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + #ifndef _TLFRelaxation_H_ #define _TLFRelaxation_H_ @@ -34,7 +54,7 @@ using namespace std; //#include "TMath.h" #include "PUserFcnBase.h" #include "fftw3.h" -#include "TIntegrator.h" +#include "../BMWIntegrator/BMWIntegrator.h" class TLFStatGssKT : public PUserFcnBase { @@ -112,14 +132,25 @@ private: double *fFFTtime; fftw_complex *fFFTfreq; mutable unsigned int fCounter; - static const double fX[16]; - static const double fY[20]; - mutable double fA; mutable double fL1; mutable double fL2; ClassDef(TLFDynLorKT,1) }; +class TLFSGInterpolation : public PUserFcnBase { + +public: + TLFSGInterpolation(); + ~TLFSGInterpolation(); + + double operator()(double, const vector&) const; + +private: + TIntSGInterpolation *fIntegral; + + ClassDef(TLFSGInterpolation,1) +}; + #endif //_LFRelaxation_H_ diff --git a/src/external/libLFRelaxation/TLFRelaxationLinkDef.h b/src/external/libLFRelaxation/TLFRelaxationLinkDef.h index 0c21ff24..5085946e 100644 --- a/src/external/libLFRelaxation/TLFRelaxationLinkDef.h +++ b/src/external/libLFRelaxation/TLFRelaxationLinkDef.h @@ -9,6 +9,26 @@ ***************************************************************************/ +/*************************************************************************** + * Copyright (C) 2009 by Bastian M. Wojek * + * bastian.wojek@psi.ch * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program; if not, write to the * + * Free Software Foundation, Inc., * + * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * + ***************************************************************************/ + // root dictionary stuff -------------------------------------------------- #ifdef __CINT__ @@ -20,6 +40,7 @@ #pragma link C++ class TLFStatLorKT+; #pragma link C++ class TLFDynGssKT+; #pragma link C++ class TLFDynLorKT+; +#pragma link C++ class TLFSGInterpolation+; #endif //__CINT__ // root dictionary stuff --------------------------------------------------