From 81f59e3538e68b3ff7909a9dbd8ea4b4c85551a8 Mon Sep 17 00:00:00 2001 From: "Bastian M. Wojek" Date: Thu, 11 Nov 2010 23:53:52 +0000 Subject: [PATCH] Added a first try of globalization to libTFitPofB - very chaotic and mainly copy/paste of examples --- .../classes/TBulkTriVortexFieldCalc.cpp | 186 +++++++++++++- src/external/TFitPofB-lib/classes/TVortex.cpp | 238 ++++++++++++++++++ .../include/TBulkTriVortexFieldCalc.h | 16 ++ src/external/TFitPofB-lib/include/TVortex.h | 70 +++++- .../TFitPofB-lib/include/TVortexLinkDef.h | 2 + 5 files changed, 509 insertions(+), 3 deletions(-) diff --git a/src/external/TFitPofB-lib/classes/TBulkTriVortexFieldCalc.cpp b/src/external/TFitPofB-lib/classes/TBulkTriVortexFieldCalc.cpp index 9268e4c8..49d62a50 100644 --- a/src/external/TFitPofB-lib/classes/TBulkTriVortexFieldCalc.cpp +++ b/src/external/TFitPofB-lib/classes/TBulkTriVortexFieldCalc.cpp @@ -2,7 +2,7 @@ TBulkTriVortexFieldCalc.cpp - Author: Bastian M. Wojek, Alexander Maisuradze + Author: Bastian M. Wojek e-mail: bastian.wojek@psi.ch 2009/10/17 @@ -10,7 +10,7 @@ ***************************************************************************/ /*************************************************************************** - * Copyright (C) 2009 by Bastian M. Wojek, Alexander Maisuradze * + * Copyright (C) 2009 by Bastian M. Wojek * * * * * * This program is free software; you can redistribute it and/or modify * @@ -1917,3 +1917,185 @@ void TBulkTriVortexNGLFieldCalc::CalculateGrid() const { return; } + +TBulkAnisotropicTriVortexLondonFieldCalc::TBulkAnisotropicTriVortexLondonFieldCalc(const string& wisdom, const unsigned int steps) { + fWisdom = wisdom; + if (steps % 2) { + fSteps = steps + 1; + } else { + fSteps = steps; + } + fParam.resize(3); + fGridExists = false; + +#ifdef HAVE_LIBFFTW3_THREADS + int init_threads(fftw_init_threads()); + if (init_threads) { +#ifdef HAVE_GOMP + fftw_plan_with_nthreads(omp_get_num_procs()); +#else + fftw_plan_with_nthreads(2); +#endif /* HAVE_GOMP */ + } +#endif /* HAVE_LIBFFTW3_THREADS */ + + fFFTin = new fftw_complex[(fSteps/2 + 1) * fSteps]; + fFFTout = new double[fSteps*fSteps]; + +// cout << "Check for the FFT plan..." << endl; + +// Load wisdom from file if it exists and should be used + + fUseWisdom = true; + int wisdomLoaded(0); + + FILE *wordsOfWisdomR; + wordsOfWisdomR = fopen(fWisdom.c_str(), "r"); + if (wordsOfWisdomR == NULL) { + fUseWisdom = false; + } else { + wisdomLoaded = fftw_import_wisdom_from_file(wordsOfWisdomR); + fclose(wordsOfWisdomR); + } + + if (!wisdomLoaded) { + fUseWisdom = false; + } + +// create the FFT plan + + if (fUseWisdom) + fFFTplan = fftw_plan_dft_c2r_2d(fSteps, fSteps, fFFTin, fFFTout, FFTW_EXHAUSTIVE); + else + fFFTplan = fftw_plan_dft_c2r_2d(fSteps, fSteps, fFFTin, fFFTout, FFTW_ESTIMATE); +} + +void TBulkAnisotropicTriVortexLondonFieldCalc::CalculateGrid() const { + // SetParameters - method has to be called from the user before the calculation!! + if (fParam.size() < 5) { + cout << endl << "The SetParameters-method has to be called before B(x,y) can be calculated!" << endl; + return; + } + if (!fParam[0] || !fParam[1] || !fParam[2] || !fParam[3] || !fParam[4]) { + cout << endl << "The field, penetration depths and coherence lengths have to have finite values in order to calculate B(x,y)!" \ + << endl; + return; + } + + double field(fabs(fParam[0])), lambdaX(fabs(fParam[1])), lambdaY(fabs(fParam[2])), xiX(fabs(fParam[3])), xiY(fabs(fParam[4])); + double Hc2(fluxQuantum/(TWOPI*xiX*xiY)); + + double latConstTr(sqrt(2.0*fluxQuantum/(field*sqrt3))); + double xiXsq_2_scaled(field/fluxQuantum*pow(xiX*PI,2.0)*sqrt3*lambdaY/lambdaX); + double xiYsq_2_scaled(field/fluxQuantum*pow(xiY*PI,2.0)*lambdaX/(lambdaY*sqrt3)); + double lambdaXsq_scaled(2.0*field/fluxQuantum*PI*PI/(sqrt3*lambdaY)*pow(lambdaX,3.0)); + double lambdaYsq_scaled(2.0*field/fluxQuantum*PI*PI*sqrt3/lambdaX*pow(lambdaY,3.0)); + + const int NFFT(fSteps); + const int NFFT_2(fSteps/2); + const int NFFTsq(fSteps*fSteps); + + // fill the field Fourier components in the matrix + + // ... but first check that the field is not larger than Hc2 and that we are dealing with a type II SC + if ((field >= Hc2) || (sqrt(lambdaX*lambdaY) < sqrt(xiX*xiY)/sqrt(2.0))) { + int m; + #ifdef HAVE_GOMP + #pragma omp parallel for default(shared) private(m) schedule(dynamic) + #endif + for (m = 0; m < NFFTsq; ++m) { + fFFTout[m] = field; + } + // Set the flag which shows that the calculation has been done + fGridExists = true; + return; + } + + // ... now fill in the Fourier components if everything was okay above + double kk, ll; + int k, l, lNFFT_2; + + for (l = 0; l < NFFT_2; l += 2) { + lNFFT_2 = l*(NFFT_2 + 1); + ll = static_cast(l*l); + for (k = 0; k < NFFT_2; k += 2) { + kk = static_cast(k*k); + fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk); + fFFTin[lNFFT_2 + k][1] = 0.0; + fFFTin[lNFFT_2 + k + 1][0] = 0.0; + fFFTin[lNFFT_2 + k + 1][1] = 0.0; + } + k = NFFT_2; + kk = static_cast(k*k); + fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk); + fFFTin[lNFFT_2 + k][1] = 0.0; + } + + + for (l = NFFT_2; l < NFFT; l += 2) { + lNFFT_2 = l*(NFFT_2 + 1); + ll = static_cast((NFFT-l)*(NFFT-l)); + for (k = 0; k < NFFT_2; k += 2) { + kk = static_cast(k*k); + fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk); + fFFTin[lNFFT_2 + k][1] = 0.0; + fFFTin[lNFFT_2 + k + 1][0] = 0.0; + fFFTin[lNFFT_2 + k + 1][1] = 0.0; + } + k = NFFT_2; + kk = static_cast(k*k); + fFFTin[lNFFT_2 + k][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk); + fFFTin[lNFFT_2 + k][1] = 0.0; + } + + // intermediate rows + + for (l = 1; l < NFFT_2; l += 2) { + lNFFT_2 = l*(NFFT_2 + 1); + ll = static_cast(l*l); + for (k = 0; k < NFFT_2; k += 2) { + kk = static_cast((k + 1)*(k + 1)); + fFFTin[lNFFT_2 + k][0] = 0.0; + fFFTin[lNFFT_2 + k][1] = 0.0; + fFFTin[lNFFT_2 + k + 1][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk); + fFFTin[lNFFT_2 + k + 1][1] = 0.0; + } + k = NFFT_2; + fFFTin[lNFFT_2 + k][0] = 0.0; + fFFTin[lNFFT_2 + k][1] = 0.0; + } + + for (l = NFFT_2 + 1; l < NFFT; l += 2) { + lNFFT_2 = l*(NFFT_2 + 1); + ll = static_cast((NFFT-l)*(NFFT-l)); + for (k = 0; k < NFFT_2; k += 2) { + kk = static_cast((k+1)*(k+1)); + fFFTin[lNFFT_2 + k][0] = 0.0; + fFFTin[lNFFT_2 + k][1] = 0.0; + fFFTin[lNFFT_2 + k + 1][0] = exp(-(xiXsq_2_scaled*kk + xiYsq_2_scaled*ll))/(1.0+lambdaXsq_scaled*ll+lambdaYsq_scaled*kk); + fFFTin[lNFFT_2 + k + 1][1] = 0.0; + } + k = NFFT_2; + fFFTin[lNFFT_2 + k][0] = 0.0; + fFFTin[lNFFT_2 + k][1] = 0.0; + } + + // Do the Fourier transform to get B(x,y) + + fftw_execute(fFFTplan); + + // Multiply by the applied field + #ifdef HAVE_GOMP + #pragma omp parallel for default(shared) private(l) schedule(dynamic) + #endif + for (l = 0; l < NFFTsq; ++l) { + fFFTout[l] *= field; + } + + // Set the flag which shows that the calculation has been done + + fGridExists = true; + return; + +} + diff --git a/src/external/TFitPofB-lib/classes/TVortex.cpp b/src/external/TFitPofB-lib/classes/TVortex.cpp index 3fdc3ce4..fb02d5d5 100644 --- a/src/external/TFitPofB-lib/classes/TVortex.cpp +++ b/src/external/TFitPofB-lib/classes/TVortex.cpp @@ -43,6 +43,8 @@ ClassImp(TBulkSqVortexLondon) ClassImp(TBulkTriVortexML) ClassImp(TBulkTriVortexAGL) ClassImp(TBulkTriVortexNGL) +ClassImp(TBulkAnisotropicTriVortexLondonGlobal) +ClassImp(TBulkAnisotropicTriVortexLondon) //------------------ // Destructor of the TBulkTriVortexLondon class -- cleaning up @@ -819,3 +821,239 @@ double TBulkTriVortexNGL::operator()(double t, const vector &par) const return fPofT->Eval(t); } + +//------------------ +// Destructor of the TBulkAnisotropicTriVortexLondonGlobal class -- cleaning up +//------------------ + +TBulkAnisotropicTriVortexLondonGlobal::~TBulkAnisotropicTriVortexLondonGlobal() { + delete fPofT; + fPofT = 0; + delete fPofB; + fPofT = 0; + delete fVortex; + fVortex = 0; + fPar.clear(); + fParForVortex.clear(); + fParForPofB.clear(); + fParForPofT.clear(); +} + +//------------------ +// Constructor of the TBulkAnisotropicTriVortexLondonGlobal class +// creates (a pointer to) the TPofTCalc object (with the FFT plan) +//------------------ + +TBulkAnisotropicTriVortexLondonGlobal::TBulkAnisotropicTriVortexLondonGlobal() : 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 + cerr << endl << "**ERROR** reading/parsing TFitPofB_startup.xml failed." \ + << endl << "**ERROR** Please make sure that the file exists in the local directory and it is set up correctly!" \ + << endl; + assert(false); + } + + fGridSteps = startupHandler->GetGridSteps(); + fWisdom = startupHandler->GetWisdomFile(); + + fParForVortex.resize(5); // field, lambdaX, lambdaY, xiX, xiY + + 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 + + fVortex = new TBulkAnisotropicTriVortexLondonFieldCalc(fWisdom, fGridSteps); + + fPofB = new TPofBCalc(fParForPofB); + + fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT); + + // clean up + if (saxParser) { + delete saxParser; + saxParser = 0; + } + if (startupHandler) { + delete startupHandler; + startupHandler = 0; + } +} + +void TBulkAnisotropicTriVortexLondonGlobal::CalcPofB(const vector &par) const { + + assert(par.size() == 6); +/* + 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(0); i < 5; i++) { + fParForVortex[i] = fPar[i+1]; + } + fFirstCall = false; + } + + // check if any parameter has changed + + bool par_changed(false); + bool only_phase_changed(false); + + for (unsigned int i(0); iSetParameters(fParForVortex); + fVortex->CalculateGrid(); + fPofB->UnsetPBExists(); + fPofB->Calculate(fVortex, fParForPofB); + fPofT->DoFFT(); + + }/* else { + cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl; + }*/ + + fPofT->CalcPol(fParForPofT); + + fCalcNeeded = false; + } +/* + return fPofT->Eval(t); +*/ + fValid = true; + return; +} + +//------------------------------------------------------------------------------------ +/** + *

Constructor. + */ +TBulkAnisotropicTriVortexLondon::TBulkAnisotropicTriVortexLondon() +{ + fValid = false; + fInvokedGlobal = false; + fIdxGlobal = -1; + fGlobalUserFcn = 0; +} + +//------------------------------------------------------------------------------------ +/** + *

Destructor. + */ +TBulkAnisotropicTriVortexLondon::~TBulkAnisotropicTriVortexLondon() +{ + if ((fGlobalUserFcn != 0) && fInvokedGlobal) { + delete fGlobalUserFcn; + fGlobalUserFcn = 0; + } +} + +//------------------------------------------------------------------------------------ +/** + *

Used to invoke/retrieve the proper global user function + * + * \param globalPart reference to the global user function object vector + * \param idx global user function index within the theory tree + */ +void TBulkAnisotropicTriVortexLondon::SetGlobalPart(vector &globalPart, UInt_t idx) +{ + fIdxGlobal = static_cast(idx); + + if ((Int_t)globalPart.size() <= fIdxGlobal) { // global user function not present, invoke it + fGlobalUserFcn = new TBulkAnisotropicTriVortexLondonGlobal(); + if (fGlobalUserFcn == 0) { // global user function object couldn't be invoked -> error + fValid = false; + cerr << endl << ">> TBulkAnisotropicTriVortexLondon::SetGlobalPart(): **ERROR** Couldn't invoke global user function object, sorry ..." << endl; + } else { // global user function object could be invoked -> resize to global user function vector and keep the pointer to the corresponding object + globalPart.resize(fIdxGlobal+1); + globalPart[fIdxGlobal] = dynamic_cast(fGlobalUserFcn); + fValid = true; + fInvokedGlobal = true; + } + } else { // global user function already present hence just retrieve a pointer to it + fValid = true; + fGlobalUserFcn = (TBulkAnisotropicTriVortexLondonGlobal*)globalPart[fIdxGlobal]; + } +} + +//------------------------------------------------------------------------------------ +/** + *

Used to check if the user function is OK. + * + * return: true if both, the user function and the global user function object are valid + */ +bool TBulkAnisotropicTriVortexLondon::GlobalPartIsValid() const +{ + return (fValid && fGlobalUserFcn->IsValid()); +} + +//------------------------------------------------------------------------------------ +/** + *

function operator + * + * return: function value + * + * \param t time in \f$(\mu\mathrm{s})\f$, or x-axis value for non-muSR fit + * \param param parameter vector + */ +double TBulkAnisotropicTriVortexLondon::operator()(double t, const std::vector ¶m) const +{ + // expected parameters: phase, field, lambdaX, lambdaY, xiX, xiY + + assert(param.size() == 6); + assert(fGlobalUserFcn); + + if(t<0.0) + return cos(param[0]*0.017453293); + + // call the global user function object + fGlobalUserFcn->CalcPofB(param); + + return fGlobalUserFcn->GetPolarizationPointer()->Eval(t); +} diff --git a/src/external/TFitPofB-lib/include/TBulkTriVortexFieldCalc.h b/src/external/TFitPofB-lib/include/TBulkTriVortexFieldCalc.h index c623f096..534a2675 100644 --- a/src/external/TFitPofB-lib/include/TBulkTriVortexFieldCalc.h +++ b/src/external/TFitPofB-lib/include/TBulkTriVortexFieldCalc.h @@ -86,6 +86,22 @@ public: }; +/** + *

Class for the calculation of the spatial field distribution B(x,y) within a 2D "triangular vortex lattice" + * using the London model with a Gaussian cutoff for an anisotropic superconductor with the field applied along one of the principal axes + */ +class TBulkAnisotropicTriVortexLondonFieldCalc : public TBulkVortexFieldCalc { + +public: + + TBulkAnisotropicTriVortexLondonFieldCalc(const string&, const unsigned int steps = 256); + ~TBulkAnisotropicTriVortexLondonFieldCalc() {} + + void CalculateGrid() const; + bool IsTriangular() const {return true;} + +}; + /** *

Class for the calculation of the spatial field distribution B(x,y) within a 2D square vortex lattice * using the London model with a Gaussian cutoff diff --git a/src/external/TFitPofB-lib/include/TVortex.h b/src/external/TFitPofB-lib/include/TVortex.h index 368d6936..10b2b39d 100644 --- a/src/external/TFitPofB-lib/include/TVortex.h +++ b/src/external/TFitPofB-lib/include/TVortex.h @@ -200,4 +200,72 @@ private: ClassDef(TBulkTriVortexNGL,1) }; -#endif //_TLondon1D_H_ +/** + *

Class implementing the global part of the musrfit user function interface for calculating muon spin depolarization functions + * originating from the spatial field distribution B(x,y) within a 2D "triangular vortex lattice" + * calculated using the London model with a Gaussian cutoff for an anisotropic superconductor with the field applied along + * one of the principal axes + */ +class TBulkAnisotropicTriVortexLondonGlobal { + +public: + // default constructor and destructor + TBulkAnisotropicTriVortexLondonGlobal(); + ~TBulkAnisotropicTriVortexLondonGlobal(); + + bool IsValid() { return fValid; } + + // the following function will check if something needs to be calculated, which + // is the case if param != fPrevParam + void CalcPofB(const vector&) const; + + // this routine will return the calculated values, e.g. B(z,E) for TMyFunction::operator()() + // (...) here sketches only that some parameters are likley to be fed + const TPofTCalc* GetPolarizationPointer() const { return fPofT; } + +private: + mutable bool fValid; + mutable vector fPar; + TBulkAnisotropicTriVortexLondonFieldCalc *fVortex; ///< spatial field distribution B(x,y) within the vortex lattice + TPofBCalc *fPofB; ///< static field distribution P(B) + TPofTCalc *fPofT; ///< muon spin polarization p(t) + mutable bool fCalcNeeded; ///< tag needed to avoid unnecessary calculations if the core parameters were unchanged + mutable bool fFirstCall; ///< tag for checking if the function operator is called the first time + mutable vector fParForVortex; ///< parameters for the calculation of B(x,y) + mutable vector fParForPofB; ///< parameters for the calculation of P(B) + mutable vector fParForPofT; ///< parameters for the calculation of p(t) + string fWisdom; ///< file name of the FFTW wisdom file + unsigned int fGridSteps; ///< number of points in x- and y-direction for which B(x,y) is calculated + + // definition of the class for the ROOT-dictionary + ClassDef(TBulkAnisotropicTriVortexLondonGlobal,1) +}; + +/** + *

Class implementing the musrfit user function interface for calculating muon spin depolarization functions + * originating from the spatial field distribution B(x,y) within a 2D "triangular vortex lattice" + * calculated using the London model with a Gaussian cutoff for an anisotropic superconductor with the field applied along + * one of the principal axes + */ +class TBulkAnisotropicTriVortexLondon : public PUserFcnBase { + +public: + TBulkAnisotropicTriVortexLondon(); + ~TBulkAnisotropicTriVortexLondon(); + + virtual bool NeedGlobalPart() const { return true; } + virtual void SetGlobalPart(vector &globalPart, unsigned int idx); + virtual bool GlobalPartIsValid() const; + + double operator()(double, const vector&) const; + +private: + bool fValid; + bool fInvokedGlobal; + int fIdxGlobal; + TBulkAnisotropicTriVortexLondonGlobal *fGlobalUserFcn; + + ClassDef(TBulkAnisotropicTriVortexLondon,1) +}; + +#endif //_TVortex_H_ diff --git a/src/external/TFitPofB-lib/include/TVortexLinkDef.h b/src/external/TFitPofB-lib/include/TVortexLinkDef.h index b2efbe95..7ec587ea 100644 --- a/src/external/TFitPofB-lib/include/TVortexLinkDef.h +++ b/src/external/TFitPofB-lib/include/TVortexLinkDef.h @@ -41,6 +41,8 @@ #pragma link C++ class TBulkTriVortexML+; #pragma link C++ class TBulkTriVortexAGL+; #pragma link C++ class TBulkTriVortexNGL+; +#pragma link C++ class TBulkAnisotropicTriVortexLondon+; +#pragma link C++ class TBulkAnisotropicTriVortexLondonGlobal+; #endif //__CINT__ // root dictionary stuff --------------------------------------------------