Added a first try of globalization to libTFitPofB - very chaotic and mainly copy/paste of examples

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<double>(l*l);
+    for (k = 0; k < NFFT_2; k += 2) {
+      kk = static_cast<double>(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<double>(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<double>((NFFT-l)*(NFFT-l));
+    for (k = 0; k < NFFT_2; k += 2) {
+      kk = static_cast<double>(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<double>(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<double>(l*l);
+    for (k = 0; k < NFFT_2; k += 2) {
+      kk = static_cast<double>((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<double>((NFFT-l)*(NFFT-l));
+    for (k = 0; k < NFFT_2; k += 2) {
+      kk = static_cast<double>((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;
+
+}
+

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<double> &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<double> &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); i<fPar.size(); i++) {
+    if( fPar[i]-par[i] ) {
+      fPar[i] = par[i];
+      par_changed = true;
+      if (i == 0) {
+        only_phase_changed = true;
+      } else {
+        only_phase_changed = false;
+      }
+    }
+  }
+
+  if (par_changed)
+    fCalcNeeded = true;
+
+  // if model parameters have changed, recalculate B(x,y), P(B) and P(t)
+
+  if (fCalcNeeded) {
+
+    fParForPofT[0] = par[0]; // phase
+
+    if(!only_phase_changed) {
+
+//      cout << " Parameters have changed, (re-)calculating p(B) and P(t) now..." << endl;
+
+      for (unsigned int i(0); i < 5; i++) {
+        fParForVortex[i] = par[i+1];
+      }
+
+      fParForPofB[2] = par[1]; // Bkg-Field
+      //fParForPofB[3] = 0.005; // Bkg-width (in principle zero)
+
+      fVortex->SetParameters(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;
+}
+
+//------------------------------------------------------------------------------------
+/**
+ * <p> Constructor.
+ */
+TBulkAnisotropicTriVortexLondon::TBulkAnisotropicTriVortexLondon()
+{
+  fValid = false;
+  fInvokedGlobal = false;
+  fIdxGlobal = -1;
+  fGlobalUserFcn = 0;
+}
+
+//------------------------------------------------------------------------------------
+/**
+ * <p> Destructor.
+ */
+TBulkAnisotropicTriVortexLondon::~TBulkAnisotropicTriVortexLondon()
+{
+  if ((fGlobalUserFcn != 0) && fInvokedGlobal) {
+    delete fGlobalUserFcn;
+    fGlobalUserFcn = 0;
+  }
+}
+
+//------------------------------------------------------------------------------------
+/**
+ * <p> 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<void *> &globalPart, UInt_t idx)
+{
+  fIdxGlobal = static_cast<Int_t>(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<TBulkAnisotropicTriVortexLondonGlobal*>(fGlobalUserFcn);
+      fValid = true;
+      fInvokedGlobal = true;
+    }
+  } else { // global user function already present hence just retrieve a pointer to it
+    fValid = true;
+    fGlobalUserFcn = (TBulkAnisotropicTriVortexLondonGlobal*)globalPart[fIdxGlobal];
+  }
+}
+
+//------------------------------------------------------------------------------------
+/**
+ * <p> Used to check if the user function is OK.
+ *
+ * <b>return:</b> true if both, the user function and the global user function object are valid
+ */
+bool TBulkAnisotropicTriVortexLondon::GlobalPartIsValid() const
+{
+  return (fValid && fGlobalUserFcn->IsValid());
+}
+
+//------------------------------------------------------------------------------------
+/**
+ * <p> function operator
+ *
+ * <b>return:</b> 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<double> &param) 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);
+}

src/external/TFitPofB-lib/include/TBulkTriVortexFieldCalc.h

@@ -86,6 +86,22 @@ public:
 
 };
 
+/**
+ * <p>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;}
+
+};
+
 /**
  * <p>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

src/external/TFitPofB-lib/include/TVortex.h

@@ -200,4 +200,72 @@ private:
   ClassDef(TBulkTriVortexNGL,1)
 };
 
-#endif //_TLondon1D_H_
+/**
+ * <p>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<double>&) 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<double> 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<double> fParForVortex; ///< parameters for the calculation of B(x,y)
+  mutable vector<double> fParForPofB; ///< parameters for the calculation of P(B)
+  mutable vector<double> 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)
+};
+
+/**
+ * <p>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<void *> &globalPart, unsigned int idx);
+  virtual bool GlobalPartIsValid() const;
+
+  double operator()(double, const vector<double>&) const;
+
+private:
+  bool fValid;
+  bool fInvokedGlobal;
+  int  fIdxGlobal;
+  TBulkAnisotropicTriVortexLondonGlobal *fGlobalUserFcn;
+
+  ClassDef(TBulkAnisotropicTriVortexLondon,1)
+};
+
+#endif //_TVortex_H_

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 --------------------------------------------------