Added the NGL calculation to libTFitPofB - some tests need to be done and some more checks have to be implemented

src/external/TFitPofB-lib/classes/TVortex.cpp

@@ -39,6 +39,7 @@ using namespace std;
 #include "TFitPofBStartupHandler.h"
 
 ClassImp(TBulkTriVortexLondon)
+ClassImp(TBulkTriVortexNGL)
 
 //------------------
 // Destructor of the TBulkTriVortexLondon class -- cleaning up
@@ -57,6 +58,23 @@ TBulkTriVortexLondon::~TBulkTriVortexLondon() {
     fParForPofT.clear();
 }
 
+//------------------
+// Destructor of the TBulkTriVortexNGL class -- cleaning up
+//------------------
+
+TBulkTriVortexNGL::~TBulkTriVortexNGL() {
+    delete fPofT;
+    fPofT = 0;
+    delete fPofB;
+    fPofT = 0;
+    delete fVortex;
+    fVortex = 0;
+    fPar.clear();
+    fParForVortex.clear();
+    fParForPofB.clear();
+    fParForPofT.clear();
+}
+
 //------------------
 // Constructor of the TBulkTriVortexLondon class
 // creates (a pointer to) the TPofTCalc object (with the FFT plan)
@@ -194,3 +212,141 @@ double TBulkTriVortexLondon::operator()(double t, const vector<double> &par) con
   return fPofT->Eval(t);
 
 }
+
+//------------------
+// Constructor of the TBulkTriVortexNGL class
+// creates (a pointer to) the TPofTCalc object (with the FFT plan)
+//------------------
+
+TBulkTriVortexNGL::TBulkTriVortexNGL() : 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();
+
+    fParForVortex.resize(3); // field, lambda, xi
+
+    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
+
+    TBulkTriVortexNGLFieldCalc *x = new TBulkTriVortexNGLFieldCalc(fWisdom, fGridSteps);
+    fVortex = x;
+    x = 0;
+
+    TPofBCalc *y = new TPofBCalc(fParForPofB);
+    fPofB = y;
+    y = 0;
+
+    TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT);
+    fPofT = z;
+    z = 0;
+
+    // clean up
+    if (saxParser) {
+      delete saxParser;
+      saxParser = 0;
+    }
+    if (startupHandler) {
+      delete startupHandler;
+      startupHandler = 0;
+    }
+}
+
+//------------------
+// TBulkTriVortexNGL-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 TBulkTriVortexNGL (phase, appl.field, lambda, xi, [not implemented: bkg weight])
+//------------------
+
+double TBulkTriVortexNGL::operator()(double t, const vector<double> &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(0); i < 3; 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 < 3; 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);
+
+}

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

@@ -62,8 +62,8 @@ public:
 protected:
   vector<double> fParam;
   unsigned int fSteps;  // number of steps, the "unit cell" of the vortex lattice is devided in (in x- and y- direction)
-  fftw_complex *fFFTin; // Fourier components of the field
+  mutable fftw_complex *fFFTin; // Fourier components of the field
-  double *fFFTout; // fields in a "unit cell" of the vortex lattice
+  mutable double *fFFTout; // fields in a "unit cell" of the vortex lattice
   fftw_plan fFFTplan;
   bool fUseWisdom;
   string fWisdom;
@@ -97,40 +97,29 @@ public:
   ~TBulkTriVortexNGLFieldCalc();
 
   void CalculateGrid() const;
-  fftw_complex* GetTempMatrix() const {return fTempMatrix;}
+
   fftw_complex* GetAkMatrix() const {return fFFTin;}
   fftw_complex* GetBkMatrix() const {return fBkMatrix;}
   double* GetOmegaMatrix() const {return fOmegaMatrix;}
   double* GetBMatrix() const {return fFFTout;}
   fftw_complex* GetOmegaDiffMatrix() const {return fOmegaDiffMatrix;}
-//  double* GetOmegaDiffYMatrix() const {return fOmegaDiffYMatrix;}
   fftw_complex* GetQMatrix() const {return fQMatrix;}
-//  double* GetQyMatrix() const {return fQyMatrix;}
-//  double* GetAbrikosovCheck() const {return fAbrikosovCheck;}
 
 private:
 
   void CalculateGradient() const;
+  void CalculateSumAk() const;
   void FillAbrikosovCoefficients() const;
-  void ManipulateFourierCoefficients(fftw_complex*, const double&, const double&) const;
+  void ManipulateFourierCoefficientsA() const;
+  void ManipulateFourierCoefficientsB() const;
   void ManipulateFourierCoefficientsForQx() const;
   void ManipulateFourierCoefficientsForQy() const;
 
-  fftw_complex *fTempMatrix;
-  fftw_complex *fBkMatrix;
-//  mutable fftw_complex *fRealSpaceMatrix;
-//  mutable double *fAbrikosovCheck;
   mutable double *fOmegaMatrix;
-//  mutable double *fOmegaSqMatrix;
   mutable fftw_complex *fOmegaDiffMatrix;
-//  mutable double *fOmegaDiffYMatrix;
+  mutable fftw_complex *fBkMatrix;
-//  mutable double *fOmegaDDiffXMatrix;
-//  mutable double *fOmegaDDiffYMatrix;
   mutable fftw_complex *fQMatrix;
-//  mutable double *fQyMatrix;
   mutable fftw_complex *fQMatrixA;
-//  mutable double *fQyMatrixA;
-//  mutable double *fGMatrix;
 
   mutable double *fCheckAkConvergence;
   mutable double *fCheckBkConvergence;

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

@@ -59,4 +59,28 @@ private:
   ClassDef(TBulkTriVortexLondon,1)
 };
 
+class TBulkTriVortexNGL : public PUserFcnBase {
+
+public:
+  TBulkTriVortexNGL();
+  ~TBulkTriVortexNGL();
+
+  double operator()(double, const vector<double>&) const;
+
+private:
+  mutable vector<double> fPar;
+  TBulkTriVortexNGLFieldCalc *fVortex;
+  TPofBCalc *fPofB;
+  TPofTCalc *fPofT;
+  mutable bool fCalcNeeded;
+  mutable bool fFirstCall;
+  mutable vector<double> fParForVortex;
+  mutable vector<double> fParForPofB;
+  mutable vector<double> fParForPofT;
+  string fWisdom;
+  unsigned int fGridSteps;
+
+  ClassDef(TBulkTriVortexNGL,1)
+};
+
 #endif //_TLondon1D_H_

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

@@ -37,6 +37,7 @@
 #pragma link off all functions;
 
 #pragma link C++ class TBulkTriVortexLondon+;
+#pragma link C++ class TBulkTriVortexNGL+;
 
 #endif //__CINT__
 // root dictionary stuff --------------------------------------------------