Submission of a few minor changes of the last time

src/classes/PFitter.cpp

@@ -1128,6 +1128,7 @@ Bool_t PFitter::ExecuteSave()
       else
         hcorr->Draw("COLZ");
       ccorr->Write("ccorr", TObject::kOverwrite, sizeof(ccorr));
+      hcorr->Write("hcorr", TObject::kOverwrite, sizeof(hcorr));
       ff.Close();
       // clean up
       if (ccorr) {

src/external/BMWIntegrator/BMWIntegrator.h

@@ -61,10 +61,7 @@ class TIntegrator {
 };
 
 inline TIntegrator::TIntegrator() : fFunc(0) {
-  ROOT::Math::GSLIntegrator *integrator = new ROOT::Math::GSLIntegrator(ROOT::Math::Integration::kADAPTIVE,ROOT::Math::Integration::kGAUSS31);
+  fIntegrator = new ROOT::Math::GSLIntegrator(ROOT::Math::Integration::kADAPTIVE,ROOT::Math::Integration::kGAUSS31);
-  fIntegrator = integrator;
-  integrator = 0;
-  delete integrator;
 }
 
 inline TIntegrator::~TIntegrator(){
@@ -104,10 +101,7 @@ class TMCIntegrator {
 };
 
 inline TMCIntegrator::TMCIntegrator() : fFunc(0) {
-  ROOT::Math::GSLMCIntegrator *integrator = new ROOT::Math::GSLMCIntegrator(ROOT::Math::MCIntegration::kMISER, 1.E-6, 1.E-4, 500000);
+  fMCIntegrator = new ROOT::Math::GSLMCIntegrator(ROOT::Math::MCIntegration::kMISER, 1.E-6, 1.E-4, 500000);
-  fMCIntegrator = integrator;
-  integrator = 0;
-  delete integrator;
 }
 
 inline TMCIntegrator::~TMCIntegrator(){

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

@@ -286,6 +286,130 @@ void TBulkTriVortexLondonFieldCalc::CalculateGrid() const {
 
 }
 
+TBulkSqVortexLondonFieldCalc::TBulkSqVortexLondonFieldCalc(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)
+    fftw_plan_with_nthreads(2);
+#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 TBulkSqVortexLondonFieldCalc::CalculateGrid() const {
+  // SetParameters - method has to be called from the user before the calculation!!
+  if (fParam.size() < 3) {
+    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]) {
+    cout << endl << "The field, penetration depth and coherence length have to have finite values in order to calculate B(x,y)!" << endl;
+    return;
+  }
+
+  double field(fabs(fParam[0])), lambda(fabs(fParam[1])), xi(fabs(fParam[2]));
+  double Hc2(getHc2(xi));
+
+  double latConstSq(sqrt(fluxQuantum/field));
+  double xisq_2_scaled(2.0*pow(xi*PI/latConstSq,2.0)), lambdasq_scaled(4.0*pow(lambda*PI/latConstSq,2.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) || (lambda < xi/sqrt(2.0))) {
+    int m;
+    #pragma omp parallel for default(shared) private(m) schedule(dynamic)
+    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 Gsq, ll;
+  int k, l, lNFFT_2;
+
+  for (l = 0; l < NFFT_2; ++l) {
+    lNFFT_2 = l*(NFFT_2 + 1);
+    ll = static_cast<double>(l*l);
+    for (k = 0; k <= NFFT_2; ++k) {
+      Gsq = static_cast<double>(k*k) + ll;
+      fFFTin[lNFFT_2 + k][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq);
+      fFFTin[lNFFT_2 + k][1] = 0.0;
+    }
+  }
+
+  for (l = NFFT_2; l < NFFT; ++l) {
+    lNFFT_2 = l*(NFFT_2 + 1);
+    ll = static_cast<double>((NFFT-l)*(NFFT-l));
+    for (k = 0; k <= NFFT_2; ++k) {
+      Gsq = static_cast<double>(k*k) + ll;
+      fFFTin[lNFFT_2 + k][0] = exp(-xisq_2_scaled*Gsq)/(1.0+lambdasq_scaled*Gsq);
+      fFFTin[lNFFT_2 + k][1] = 0.0;
+    }
+  }
+
+  // Do the Fourier transform to get B(x,y)
+
+  fftw_execute(fFFTplan);
+
+  // Multiply by the applied field
+  #pragma omp parallel for default(shared) private(l) schedule(dynamic)
+  for (l = 0; l < NFFTsq; l++) {
+    fFFTout[l] *= field;
+  }
+
+  // Set the flag which shows that the calculation has been done
+
+  fGridExists = true;
+  return;
+
+}
+
+
 
 TBulkTriVortexMLFieldCalc::TBulkTriVortexMLFieldCalc(const string& wisdom, const unsigned int steps) {
   fWisdom = wisdom;

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

@@ -49,7 +49,7 @@ ClassImpQ(TFitPofBStartupHandler)
 /**
  * <p>
  */
-TFitPofBStartupHandler::TFitPofBStartupHandler() : fDeltat(0.), fDeltaB(0.), fNSteps(0)
+TFitPofBStartupHandler::TFitPofBStartupHandler() : fDeltat(0.), fDeltaB(0.), fNSteps(0), fGridSteps(0)
 {
 }
 

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

@@ -209,18 +209,11 @@ TLondon1DHS::TLondon1DHS() : fCalcNeeded(true), fFirstCall(true) {
     fParForPofB.push_back(0.005); // Bkg-width
     fParForPofB.push_back(0.0); // Bkg-weight
 
-    TTrimSPData *x = new TTrimSPData(rge_path, energy_vec);
+    fImpProfile = new TTrimSPData(rge_path, energy_vec);
-    fImpProfile = x;
-    x = 0;
 
-    TPofBCalc *y = new TPofBCalc(fParForPofB);
+    fPofB = new TPofBCalc(fParForPofB);
-    fPofB = y;
-    y = 0;
-
-    TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT);
-    fPofT = z;
-    z = 0;
 
+    fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT);
 
     // clean up
     if (saxParser) {
@@ -370,17 +363,11 @@ TLondon1D1L::TLondon1D1L() : fCalcNeeded(true), fFirstCall(true), fCallCounter(0
     fParForPofB.push_back(startupHandler->GetDeltaB());
     fParForPofB.push_back(0.0);
 
-    TTrimSPData *x = new TTrimSPData(rge_path, energy_vec);
+    fImpProfile = new TTrimSPData(rge_path, energy_vec);
-    fImpProfile = x;
-    x = 0;
 
-    TPofBCalc *y = new TPofBCalc(fParForPofB);
+    fPofB = new TPofBCalc(fParForPofB);
-    fPofB = y;
-    y = 0;
 
-    TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT);
+    fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT);
-    fPofT = z;
-    z = 0;
 
     // clean up
     if (saxParser) {
@@ -524,18 +511,11 @@ TLondon1D2L::TLondon1D2L() : fCalcNeeded(true), fFirstCall(true), fLastTwoChange
     fParForPofB.push_back(startupHandler->GetDeltaB());
     fParForPofB.push_back(0.0);
 
-    TTrimSPData *x = new TTrimSPData(rge_path, energy_vec);
+    fImpProfile = new TTrimSPData(rge_path, energy_vec);
-    fImpProfile = x;
-    x = 0;
 
-    TPofBCalc *y = new TPofBCalc(fParForPofB);
+    fPofB = new TPofBCalc(fParForPofB);
-    fPofB = y;
-    y = 0;
-
-    TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT);
-    fPofT = z;
-    z = 0;
 
+    fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT);
 
     // clean up
     if (saxParser) {
@@ -684,17 +664,11 @@ TProximity1D1LHS::TProximity1D1LHS() : fCalcNeeded(true), fFirstCall(true) {
     fParForPofB.push_back(0.01); // Bkg-width
     fParForPofB.push_back(0.0); // Bkg-weight
 
-    TTrimSPData *x = new TTrimSPData(rge_path, energy_vec);
+    fImpProfile = new TTrimSPData(rge_path, energy_vec);
-    fImpProfile = x;
-    x = 0;
 
-    TPofBCalc *y = new TPofBCalc(fParForPofB);
+    fPofB = new TPofBCalc(fParForPofB);
-    fPofB = y;
-    y = 0;
 
-    TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT);
+    fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT);
-    fPofT = z;
-    z = 0;
 
     // clean up
     if (saxParser) {
@@ -851,17 +825,11 @@ TProximity1D1LHSGss::TProximity1D1LHSGss() : fCalcNeeded(true), fFirstCall(true)
     fParForPofB.push_back(0.0);
 //    fParForPofB.push_back(0.0);
 
-    TTrimSPData *x = new TTrimSPData(rge_path, energy_vec);
+    fImpProfile = new TTrimSPData(rge_path, energy_vec);
-    fImpProfile = x;
-    x = 0;
 
-    TPofBCalc *y = new TPofBCalc(fParForPofB);
+    fPofB = new TPofBCalc(fParForPofB);
-    fPofB = y;
-    y = 0;
 
-    TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT);
+    fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT);
-    fPofT = z;
-    z = 0;
 
     // clean up
     if (saxParser) {
@@ -991,18 +959,11 @@ TLondon1D3L::TLondon1D3L() : fCalcNeeded(true), fFirstCall(true), fLastThreeChan
     fParForPofB.push_back(startupHandler->GetDeltaB());
     fParForPofB.push_back(0.0);
 
-    TTrimSPData *x = new TTrimSPData(rge_path, energy_vec);
+    fImpProfile = new TTrimSPData(rge_path, energy_vec);
-    fImpProfile = x;
-    x = 0;
 
-    TPofBCalc *y = new TPofBCalc(fParForPofB);
+    fPofB = new TPofBCalc(fParForPofB);
-    fPofB = y;
-    y = 0;
-
-    TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT);
-    fPofT = z;
-    z = 0;
 
+    fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT);
 
     // clean up
     if (saxParser) {
@@ -1162,17 +1123,11 @@ TLondon1D3LS::TLondon1D3LS() : fCalcNeeded(true), fFirstCall(true), fLastThreeCh
     fParForPofB.push_back(startupHandler->GetDeltaB());
     fParForPofB.push_back(0.0);
 
-    TTrimSPData *x = new TTrimSPData(rge_path, energy_vec);
+    fImpProfile = new TTrimSPData(rge_path, energy_vec);
-    fImpProfile = x;
-    x = 0;
 
-    TPofBCalc *y = new TPofBCalc(fParForPofB);
+    fPofB = new TPofBCalc(fParForPofB);
-    fPofB = y;
-    y = 0;
 
-    TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT);
+    fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT);
-    fPofT = z;
-    z = 0;
 
     // clean up
     if (saxParser) {
@@ -1314,15 +1269,9 @@ double TLondon1D3LS::operator()(double t, const vector<double> &par) const {
 //     fParForPofB.push_back(startupHandler->GetDeltaB());
 //     fParForPofB.push_back(0.0);
 // 
-//     TTrimSPData *x = new TTrimSPData(rge_path, energy_vec);
+//     fImpProfile = new TTrimSPData(rge_path, energy_vec);
-//     fImpProfile = x;
-//     x = 0;
-//     delete x;
 // 
-//     TPofTCalc *y = new TPofTCalc(fWisdom, fParForPofT);
+//     fPofT = new TPofTCalc(fWisdom, fParForPofT);
-//     fPofT = y;
-//     y = 0;
-//     delete y;
 // 
 //     // clean up
 //     if (saxParser) {
@@ -1485,17 +1434,11 @@ TLondon1D3LSub::TLondon1D3LSub() : fCalcNeeded(true), fFirstCall(true), fWeights
     fParForPofB.push_back(startupHandler->GetDeltaB());
     fParForPofB.push_back(0.0);
 
-    TTrimSPData *x = new TTrimSPData(rge_path, energy_vec);
+    fImpProfile = new TTrimSPData(rge_path, energy_vec);
-    fImpProfile = x;
-    x = 0;
 
-    TPofBCalc *y = new TPofBCalc(fParForPofB);
+    fPofB = new TPofBCalc(fParForPofB);
-    fPofB = y;
-    y = 0;
 
-    TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT);
+    fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT);
-    fPofT = z;
-    z = 0;
 
     // clean up
     if (saxParser) {

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

@@ -388,12 +388,12 @@ void TPofBCalc::Calculate(const TBulkVortexFieldCalc *vortexLattice, const vecto
   int a3(static_cast<int>(floor(fBmax/fDB)));
   int a4(static_cast<int>(ceil(fBmax/fDB)));
 
-  unsigned int firstZerosEnd ((a1 < a2) ? a1 : ((a1 > 0) ? (a1 - 1) : 0));
+  //unsigned int firstZerosEnd ((a1 < a2) ? a1 : ((a1 > 0) ? (a1 - 1) : 0));
   unsigned int lastZerosStart ((a3 < a4) ? a4 : (a4 + 1));
   unsigned int numberOfSteps(vortexLattice->GetNumberOfSteps());
   unsigned int numberOfStepsSq(numberOfSteps*numberOfSteps);
   unsigned int numberOfSteps_2(numberOfSteps/2);
-  unsigned int numberOfStepsSq_2(numberOfStepsSq/2);
+  //unsigned int numberOfStepsSq_2(numberOfStepsSq/2);
 
   if (lastZerosStart >= fPBSize)
     lastZerosStart = fPBSize - 1;
@@ -405,30 +405,76 @@ void TPofBCalc::Calculate(const TBulkVortexFieldCalc *vortexLattice, const vecto
   }
 
   double *vortexFields = vortexLattice->DataB();
-  unsigned int fill_index, counter(0);
+  unsigned int fill_index;
 
-  for (unsigned int j(0); j < numberOfStepsSq_2; j++) {
+  if (para.size() == 7 && para[6] == 1.0 && para[5] != 0.0 && vortexLattice->IsTriangular()) {
-    fill_index = static_cast<unsigned int>(ceil(fabs((vortexFields[j]/fDB))));
+    // weight distribution with Gaussian around vortex-cores
-    if (fill_index >= fPBSize)
+    double Rsq1, Rsq2, Rsq3, Rsq4, Rsq5, Rsq6, sigmaSq(-0.5*para[5]*para[5]);
-      fPB[fPBSize - 1] += 1.0;
+    for (unsigned int j(0); j < numberOfSteps_2; ++j) {
-    else
+      for (unsigned int i(0); i < numberOfSteps_2; ++i) {
-      fPB[fill_index] += 1.0;
+        fill_index = static_cast<unsigned int>(ceil(fabs((vortexFields[i + numberOfSteps*j]/fDB))));
-    counter++;
+        if (fill_index < fPBSize) {
-    if (counter == numberOfSteps_2) { // sum only over the first quadrant of B(x,y)
+          Rsq1 = static_cast<double>(3*i*i + j*j)/static_cast<double>(numberOfStepsSq);
-      counter = 0;
+          Rsq2 = static_cast<double>(3*(numberOfSteps_2 - i)*(numberOfSteps_2 - i) \
-      j += numberOfSteps_2;
+               + (numberOfSteps_2 - j)*(numberOfSteps_2 - j))/static_cast<double>(numberOfStepsSq);
+          Rsq3 = static_cast<double>(3*(numberOfSteps - i)*(numberOfSteps - i) \
+               + j*j)/static_cast<double>(numberOfStepsSq);
+          Rsq4 = static_cast<double>(3*(numberOfSteps_2 - i)*(numberOfSteps_2 - i) \
+               + (numberOfSteps_2 + j)*(numberOfSteps_2 + j))/static_cast<double>(numberOfStepsSq);
+          Rsq5 = static_cast<double>(3*i*i \
+               + (numberOfSteps - j)*(numberOfSteps - j))/static_cast<double>(numberOfStepsSq);
+          Rsq6 = static_cast<double>(3*(numberOfSteps_2 + i)*(numberOfSteps_2 + i) \
+               + (numberOfSteps_2 - j)*(numberOfSteps_2 - j))/static_cast<double>(numberOfStepsSq);
+          fPB[fill_index] += exp(sigmaSq*Rsq1) + exp(sigmaSq*Rsq2) + exp(sigmaSq*Rsq3) \
+                           + exp(sigmaSq*Rsq4) + exp(sigmaSq*Rsq5) + exp(sigmaSq*Rsq6);
+        }
+      }
+    }
+  } else if (para.size() == 7 && para[6] == 2.0 && para[5] != 0.0 && vortexLattice->IsTriangular()) {
+    // weight distribution with Lorentzian around vortex-cores
+    double Rsq1, Rsq2, Rsq3, Rsq4, Rsq5, Rsq6, sigmaSq(para[5]*para[5]);
+    for (unsigned int j(0); j < numberOfSteps_2; ++j) {
+      for (unsigned int i(0); i < numberOfSteps_2; ++i) {
+        fill_index = static_cast<unsigned int>(ceil(fabs((vortexFields[i + numberOfSteps*j]/fDB))));
+        if (fill_index < fPBSize) {
+          Rsq1 = static_cast<double>(3*i*i + j*j)/static_cast<double>(numberOfStepsSq);
+          Rsq2 = static_cast<double>(3*(numberOfSteps_2 - i)*(numberOfSteps_2 - i) \
+               + (numberOfSteps_2 - j)*(numberOfSteps_2 - j))/static_cast<double>(numberOfStepsSq);
+          Rsq3 = static_cast<double>(3*(numberOfSteps - i)*(numberOfSteps - i) \
+               + j*j)/static_cast<double>(numberOfStepsSq);
+          Rsq4 = static_cast<double>(3*(numberOfSteps_2 - i)*(numberOfSteps_2 - i) \
+               + (numberOfSteps_2 + j)*(numberOfSteps_2 + j))/static_cast<double>(numberOfStepsSq);
+          Rsq5 = static_cast<double>(3*i*i \
+               + (numberOfSteps - j)*(numberOfSteps - j))/static_cast<double>(numberOfStepsSq);
+          Rsq6 = static_cast<double>(3*(numberOfSteps_2 + i)*(numberOfSteps_2 + i) \
+               + (numberOfSteps_2 - j)*(numberOfSteps_2 - j))/static_cast<double>(numberOfStepsSq);
+          fPB[fill_index] += 1.0/(1.0+sigmaSq*Rsq1) + 1.0/(1.0+sigmaSq*Rsq2) + 1.0/(1.0+sigmaSq*Rsq3) \
+                           + 1.0/(1.0+sigmaSq*Rsq4) + 1.0/(1.0+sigmaSq*Rsq5) + 1.0/(1.0+sigmaSq*Rsq6);
+        }
+      }
+    }
+  } else {
+    for (unsigned int j(0); j < numberOfSteps_2; ++j) {
+      for (unsigned int i(0); i < numberOfSteps_2; ++i) {
+        fill_index = static_cast<unsigned int>(ceil(fabs((vortexFields[i + numberOfSteps*j]/fDB))));
+        if (fill_index < fPBSize) {
+          fPB[fill_index] += 1.0;
+        }
+      }
     }
   }
 
   vortexFields = 0;
 
-  double normalizer(static_cast<double>(numberOfStepsSq_2/2)*fDB);
+  // normalize P(B)
-
+  double sum(0.0);
+  for (unsigned int i(0); i < fPBSize; ++i)
+    sum += fPB[i];
+  sum *= fDB;
   int i;
 #pragma omp parallel for default(shared) private(i) schedule(dynamic)
-  for (i = firstZerosEnd; i <= static_cast<int>(lastZerosStart); i++) {
+  for (i = 0; i < static_cast<int>(fPBSize); ++i)
-    fPB[i] /= normalizer;
+    fPB[i] /= sum;
-  }
 // end pragma omp parallel
 
   if(para.size() == 5)

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

@@ -39,6 +39,7 @@ using namespace std;
 #include "TFitPofBStartupHandler.h"
 
 ClassImp(TBulkTriVortexLondon)
+ClassImp(TBulkSqVortexLondon)
 ClassImp(TBulkTriVortexML)
 ClassImp(TBulkTriVortexAGL)
 ClassImp(TBulkTriVortexNGL)
@@ -60,6 +61,23 @@ TBulkTriVortexLondon::~TBulkTriVortexLondon() {
     fParForPofT.clear();
 }
 
+//------------------
+// Destructor of the TBulkSqVortexLondon class -- cleaning up
+//------------------
+
+TBulkSqVortexLondon::~TBulkSqVortexLondon() {
+    delete fPofT;
+    fPofT = 0;
+    delete fPofB;
+    fPofT = 0;
+    delete fVortex;
+    fVortex = 0;
+    fPar.clear();
+    fParForVortex.clear();
+    fParForPofB.clear();
+    fParForPofT.clear();
+}
+
 //------------------
 // Destructor of the TBulkTriVortexML class -- cleaning up
 //------------------
@@ -148,18 +166,69 @@ TBulkTriVortexLondon::TBulkTriVortexLondon() : fCalcNeeded(true), fFirstCall(tru
     fParForPofB.push_back(0.0); // Bkg-Field
     fParForPofB.push_back(0.005); // Bkg-width
     fParForPofB.push_back(0.0); // Bkg-weight
+    fParForPofB.push_back(0.0); // vortex-weighting
+    fParForPofB.push_back(0.0); // vortex-weighting: 0.0 homogeneous, 1.0 Gaussian, 2.0 Lorentzian
 
-    TBulkTriVortexLondonFieldCalc *x = new TBulkTriVortexLondonFieldCalc(fWisdom, fGridSteps);
+    fVortex = new TBulkTriVortexLondonFieldCalc(fWisdom, fGridSteps);
-    fVortex = x;
-    x = 0;
 
-    TPofBCalc *y = new TPofBCalc(fParForPofB);
+    fPofB = new TPofBCalc(fParForPofB);
-    fPofB = y;
-    y = 0;
 
-    TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT);
+    fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT);
-    fPofT = z;
+
-    z = 0;
+    // clean up
+    if (saxParser) {
+      delete saxParser;
+      saxParser = 0;
+    }
+    if (startupHandler) {
+      delete startupHandler;
+      startupHandler = 0;
+    }
+}
+
+//------------------
+// Constructor of the TBulkSqVortexLondon class
+// creates (a pointer to) the TPofTCalc object (with the FFT plan)
+//------------------
+
+TBulkSqVortexLondon::TBulkSqVortexLondon() : 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(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
+
+    fVortex = new TBulkSqVortexLondonFieldCalc(fWisdom, fGridSteps);
+
+    fPofB = new TPofBCalc(fParForPofB);
+
+    fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT);
 
     // clean up
     if (saxParser) {
@@ -180,6 +249,91 @@ TBulkTriVortexLondon::TBulkTriVortexLondon() : fCalcNeeded(true), fFirstCall(tru
 
 double TBulkTriVortexLondon::operator()(double t, const vector<double> &par) const {
 
+  assert(par.size() == 4 || par.size() == 5 || par.size() == 7); // normal, +BkgWeight, +VortexWeighting
+
+  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)
+      if (par.size() == 7) {
+        fParForPofB[5] = par[5];
+        assert((par[6] == 0.0) || (par[6] == 1.0) || (par[6] == 2.0));
+        fParForPofB[6] = par[6];
+      }
+
+      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);
+
+}
+
+//------------------
+// TBulkSqVortexLondon-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 TBulkSqVortexLondon (phase, av.field, lambda, xi, [not implemented: bkg weight])
+//------------------
+
+double TBulkSqVortexLondon::operator()(double t, const vector<double> &par) const {
+
   assert(par.size() == 4 || par.size() == 5);
 
   if(t<0.0)
@@ -252,6 +406,7 @@ double TBulkTriVortexLondon::operator()(double t, const vector<double> &par) con
 
 }
 
+
 //------------------
 // Constructor of the TBulkTriVortexML class
 // creates (a pointer to) the TPofTCalc object (with the FFT plan)
@@ -290,17 +445,11 @@ TBulkTriVortexML::TBulkTriVortexML() : fCalcNeeded(true), fFirstCall(true) {
     fParForPofB.push_back(0.005); // Bkg-width
     fParForPofB.push_back(0.0); // Bkg-weight
 
-    TBulkTriVortexMLFieldCalc *x = new TBulkTriVortexMLFieldCalc(fWisdom, fGridSteps);
+    fVortex = new TBulkTriVortexMLFieldCalc(fWisdom, fGridSteps);
-    fVortex = x;
-    x = 0;
 
-    TPofBCalc *y = new TPofBCalc(fParForPofB);
+    fPofB = new TPofBCalc(fParForPofB);
-    fPofB = y;
-    y = 0;
 
-    TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT);
+    fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT);
-    fPofT = z;
-    z = 0;
 
     // clean up
     if (saxParser) {
@@ -432,17 +581,11 @@ TBulkTriVortexAGL::TBulkTriVortexAGL() : fCalcNeeded(true), fFirstCall(true) {
     fParForPofB.push_back(0.005); // Bkg-width
     fParForPofB.push_back(0.0); // Bkg-weight
 
-    TBulkTriVortexAGLFieldCalc *x = new TBulkTriVortexAGLFieldCalc(fWisdom, fGridSteps);
+    fVortex = new TBulkTriVortexAGLFieldCalc(fWisdom, fGridSteps);
-    fVortex = x;
-    x = 0;
 
-    TPofBCalc *y = new TPofBCalc(fParForPofB);
+    fPofB = new TPofBCalc(fParForPofB);
-    fPofB = y;
-    y = 0;
 
-    TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT);
+    fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT);
-    fPofT = z;
-    z = 0;
 
     // clean up
     if (saxParser) {
@@ -574,17 +717,11 @@ TBulkTriVortexNGL::TBulkTriVortexNGL() : fCalcNeeded(true), fFirstCall(true) {
     fParForPofB.push_back(0.005); // Bkg-width
     fParForPofB.push_back(0.0); // Bkg-weight
 
-    TBulkTriVortexNGLFieldCalc *x = new TBulkTriVortexNGLFieldCalc(fWisdom, fGridSteps);
+    fVortex = new TBulkTriVortexNGLFieldCalc(fWisdom, fGridSteps);
-    fVortex = x;
-    x = 0;
 
-    TPofBCalc *y = new TPofBCalc(fParForPofB);
+    fPofB = new TPofBCalc(fParForPofB);
-    fPofB = y;
-    y = 0;
 
-    TPofTCalc *z = new TPofTCalc(fPofB, fWisdom, fParForPofT);
+    fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT);
-    fPofT = z;
-    z = 0;
 
     // clean up
     if (saxParser) {

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

@@ -58,6 +58,7 @@ public:
   virtual bool GridExists() const {return fGridExists;}
   virtual void UnsetGridExists() const {fGridExists = false;}
   virtual unsigned int GetNumberOfSteps() const {return fSteps;}
+  virtual bool IsTriangular() const = 0;
 
 protected:
   vector<double> fParam;
@@ -82,6 +83,23 @@ public:
   ~TBulkTriVortexLondonFieldCalc() {}
 
   void CalculateGrid() const;
+  bool IsTriangular() const {return true;}
+
+};
+
+//--------------------
+// Class for square vortex lattice, London model with Gaussian Cutoff
+//--------------------
+
+class TBulkSqVortexLondonFieldCalc : public TBulkVortexFieldCalc {
+
+public:
+
+  TBulkSqVortexLondonFieldCalc(const string&, const unsigned int steps = 256);
+  ~TBulkSqVortexLondonFieldCalc() {}
+
+  void CalculateGrid() const;
+   bool IsTriangular() const {return false;}
 
 };
 
@@ -97,6 +115,7 @@ public:
   ~TBulkTriVortexMLFieldCalc() {}
 
   void CalculateGrid() const;
+  bool IsTriangular() const {return true;}
 
 };
 
@@ -112,6 +131,7 @@ public:
   ~TBulkTriVortexAGLFieldCalc() {}
 
   void CalculateGrid() const;
+  bool IsTriangular() const {return true;}
 
 };
 
@@ -134,6 +154,7 @@ public:
   double* GetBMatrix() const {return fFFTout;}
   fftw_complex* GetOmegaDiffMatrix() const {return fOmegaDiffMatrix;}
   fftw_complex* GetQMatrix() const {return fQMatrix;}
+  bool IsTriangular() const {return true;}
 
 private:
 

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

@@ -59,6 +59,30 @@ private:
   ClassDef(TBulkTriVortexLondon,1)
 };
 
+class TBulkSqVortexLondon : public PUserFcnBase {
+
+public:
+  TBulkSqVortexLondon();
+  ~TBulkSqVortexLondon();
+
+  double operator()(double, const vector<double>&) const;
+
+private:
+  mutable vector<double> fPar;
+  TBulkSqVortexLondonFieldCalc *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(TBulkSqVortexLondon,1)
+};
+
 class TBulkTriVortexML : public PUserFcnBase {
 
 public:

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 TBulkSqVortexLondon+;
 #pragma link C++ class TBulkTriVortexML+;
 #pragma link C++ class TBulkTriVortexAGL+;
 #pragma link C++ class TBulkTriVortexNGL+;

src/external/libLFRelaxation/TLFRelaxation.cpp

@@ -53,10 +53,7 @@ ClassImp(TLFSGInterpolation)
 // LF Static Gaussian KT
 
 TLFStatGssKT::TLFStatGssKT() {
-  TIntSinGss *intSinGss = new TIntSinGss();
+  fIntSinGss = new TIntSinGss();
-  fIntSinGss = intSinGss;
-  intSinGss = 0;
-  delete intSinGss;
 }
 
 TLFStatGssKT::~TLFStatGssKT() {
@@ -89,10 +86,7 @@ double TLFStatGssKT::operator()(double t, const vector<double> &par) const {
 // LF Static Lorentzian KT
 
 TLFStatLorKT::TLFStatLorKT() {
-  TIntBesselJ0Exp *intBesselJ0Exp = new TIntBesselJ0Exp();
+  fIntBesselJ0Exp = new TIntBesselJ0Exp();
-  fIntBesselJ0Exp = intBesselJ0Exp;
-  intBesselJ0Exp = 0;
-  delete intBesselJ0Exp;
 }
 
 TLFStatLorKT::~TLFStatLorKT() {