Implemented some more London-multilayer-scenarios for testing purposes

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

@@ -195,3 +195,49 @@ TLondon1D_3LS::TLondon1D_3LS(unsigned int steps, const vector<double> &param) {
     fBZ.push_back(BBz);
   }
 }
+
+//------------------
+// Constructor of the TLondon1D_4L class
+// 1D-London screening in a thin superconducting film, four layers, four lambdas
+// Parameters: Bext[G], deadlayer[nm], thickness1[nm], thickness2[nm], thickness3[nm], thickness4[nm],
+// lambda1[nm], lambda2[nm], lambda3[nm], lambda4[nm]
+//------------------
+
+TLondon1D_4L::TLondon1D_4L(unsigned int steps, const vector<double> &param) {
+
+  double N1((param[6]+exp(2.0*param[2]/param[6])*(param[6]-param[7])+param[7])*(-1.0*param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])-param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])-param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))+exp(2.0*param[3]/param[7])*(param[6]-param[7]+exp(2.0*param[2]/param[6])*(param[6]+param[7]))*(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9])));
+
+  double N11(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])*(param[7]*cosh(param[3]/param[7])-param[6]*sinh(param[3]/param[7]))+param[7]*(-1.0*param[6]*cosh(param[3]/param[7])+param[7]*sinh(param[3]/param[7]))*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])*(-1.0*param[6]*cosh(param[3]/param[7])+param[7]*sinh(param[3]/param[7]))+param[8]*(param[7]*cosh(param[3]/param[7])-param[6]*sinh(param[3]/param[7]))*sinh(param[4]/param[8]))*sinh(param[5]/param[9]));
+
+  double N12(exp(2.0*(param[2]+param[1])/param[6])*(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])*(param[7]*cosh(param[3]/param[7])+param[6]*sinh(param[3]/param[7]))+param[7]*(param[6]*cosh(param[3]/param[7])+param[7]*sinh(param[3]/param[7]))*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])*(param[6]*cosh(param[3]/param[7])+param[7]*sinh(param[3]/param[7]))+param[8]*(param[7]*cosh(param[3]/param[7])+param[6]*sinh(param[3]/param[7]))*sinh(param[4]/param[8]))*sinh(param[5]/param[9])));
+
+  double N2((param[6]+param[7]+(param[6]-param[7])*exp(2.0*param[2]/param[6]))*(-1.0*param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])-param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])-param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))+exp(2.0*param[3]/param[7])*(param[6]-param[7]+exp(2.0*param[2]/param[6])*(param[6]+param[7]))*(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9])));
+
+  double N21(exp(param[3]/param[7])*param[8]*param[9]*(param[6]*cosh(param[2]/param[6])+param[7]*sinh(param[2]/param[6]))+param[6]*param[9]*cosh(param[5]/param[9])*(-1.0*param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+param[6]*param[8]*(param[7]*cosh(param[4]/param[8])-param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]));
+
+  double N22(exp((2.0*param[2]+param[3]+2.0*param[1])/param[7])*(-1.0*param[6]*param[8]*param[9]*cosh(param[2]/param[6])+param[7]*param[8]*param[9]*sinh(param[2]/param[6])+exp(param[3]/param[7])*param[6]*param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+exp(param[3]/param[7])*param[6]*param[8]*(param[7]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9])));
+
+  double N3(2.0*((param[6]+exp(2.0*param[2]/param[6])*(param[6]-param[7])+param[7])*(-1.0*param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])-param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])-param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))+exp(2.0*param[3]/param[7])*(param[6]-param[7]+exp(2.0*param[2]/param[6])*(param[6]+param[7]))*(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))));
+
+  double N4(4.0*((param[6]+exp(2.0*param[2]/param[6])*(param[6]-param[7])+param[7])*(-1.0*param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])-param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])-param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))+exp(2.0*param[3]/param[7])*(param[6]-param[7]+exp(2.0*param[2]/param[6])*(param[6]+param[7]))*(param[9]*cosh(param[5]/param[9])*(param[8]*cosh(param[4]/param[8])+param[7]*sinh(param[4]/param[8]))+param[8]*(param[7]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))*sinh(param[5]/param[9]))));
+
+  double N42(-1.0*param[6]*param[7]*param[8]+exp(param[5]/param[9])*(param[6]*cosh(param[2]/param[6])*(param[8]*sinh(param[3]/param[7])*(param[9]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))+param[7]*cosh(param[3]/param[7])*(param[8]*cosh(param[4]/param[8])+param[9]*sinh(param[4]/param[8])))+param[7]*sinh(param[2]/param[6])*(param[8]*cosh(param[3]/param[7])*(param[9]*cosh(param[4]/param[8])+param[8]*sinh(param[4]/param[8]))+param[7]*sinh(param[3]/param[7])*(param[8]*cosh(param[4]/param[8])+param[9]*sinh(param[4]/param[8])))));
+
+  fDZ = (param[2]+param[3]+param[4]+param[5])/double(steps);
+  double ZZ, BBz;
+
+  for (unsigned int j(0); j<steps; j++) {
+    ZZ = param[1] + (double)j*fDZ;
+    fZ.push_back(ZZ);
+    if (ZZ < param[1]+param[2]) {
+      BBz = (2.0*param[0]*exp(-(param[1]+ZZ)/param[6]+param[3]/param[7])*(-1.0*exp(param[2]/param[6])*(exp(2.0*param[1]/param[6])-exp(2.0*ZZ/param[6]))*param[7]*param[8]*param[9]-exp(2.0*ZZ/param[6])*N11+N12))/N1;
+    } else if (ZZ < param[1]+param[2]+param[3]) {
+      BBz = (2.0*param[0]*exp(param[2]/param[6]-(param[2]+param[1]+ZZ)/param[7])*(exp(2.0*ZZ/param[7])*N21+N22))/N2;
+    } else if (ZZ < param[1]+param[2]+param[3]+param[4]) {
+      BBz = (param[0]*exp(-1.0*(param[2]+param[3]+param[1]+ZZ)/param[8]-param[5]/param[9])*(2.0*exp(param[2]/param[6]+param[3]/param[7]+(2.0*param[2]+2.0*param[3]+param[4]+2.0*param[1])/param[8])*param[6]*param[7]*(param[9]-param[8]+exp(2.0*param[5]/param[9])*(param[8]+param[9]))+4.0*exp(param[2]/param[6]+param[3]/param[7]-param[4]/param[8]+param[5]/param[9])*(-1.0*exp((2.0*param[2]+2.0*param[3]+param[4]+2.0*param[1])/param[8])*param[9]*(param[7]*sinh(param[2]/param[6])*(-1.0*param[8]*cosh(param[3]/param[7])+param[7]*sinh(param[3]/param[7]))+param[6]*cosh(param[2]/param[6])*(param[7]*cosh(param[3]/param[7])-param[8]*sinh(param[3]/param[7])))+exp(2.0*ZZ/param[8])*(exp(param[4]/param[8])*param[9]*(param[7]*sinh(param[2]/param[6])*(param[8]*cosh(param[3]/param[7])+param[7]*sinh(param[3]/param[7]))+param[6]*cosh(param[2]/param[6])*(param[7]*cosh(param[3]/param[7])+param[8]*sinh(param[3]/param[7])))+param[6]*param[7]*(-1.0*param[9]*cosh(param[5]/param[9])+param[8]*sinh(param[5]/param[9]))))))/N3;
+    } else {
+      BBz = (param[0]*exp(-1.0*(param[2]+param[3]+param[4]+param[1]+ZZ)/param[9])*(-1.0*exp(-1.0*param[4]/param[8]+2.0*(param[2]+param[3]+param[4]+param[1])/param[9])*((-1.0+exp(2.0*param[2]/param[6]))*param[7]*(exp(2.0*param[4]/param[8])*(param[8]-param[7]+exp(2.0*param[3]/param[7])*(param[7]+param[8]))*(param[8]-param[9])+exp(2.0*param[3]/param[7])*(param[7]-param[8])*(param[8]+param[9])-(param[7]+param[8])*(param[8]+param[9]))+param[6]*(-8.0*exp(param[2]/param[6]+param[3]/param[7]+param[4]/param[8]+param[5]/param[9])*param[7]*param[8]+(1.0+exp(2.0*param[2]/param[6]))*(-1.0+exp(2.0*param[3]/param[7]))*param[8]*(-1.0*param[8]+exp(2.0*param[4]/param[8])*(param[8]-param[9])-param[9])+(1.0+exp(2.0*param[2]/param[6]))*(1.0+exp(2.0*param[3]/param[7]))*param[7]*(param[8]+exp(2.0*param[4]/param[8])*(param[8]-param[9])+param[9])))+8.0*exp(param[2]/param[6]+param[3]/param[7]-(param[5]-2.0*ZZ)/param[9])*N42))/N4;
+    }
+    fBZ.push_back(BBz);
+  }
+}

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

@@ -5,12 +5,13 @@
   Author: Bastian M. Wojek
   e-mail: bastian.wojek@psi.ch
 
-  2008/08/28
+  2008/11/07
 
 ***************************************************************************/
 
 #include "TLondon1D.h"
 #include <iostream>
+#include <cassert>
 using namespace std;
 
 #include <TSAXParser.h>
@@ -21,6 +22,8 @@ ClassImp(TLondon1D1L)
 ClassImp(TLondon1D2L)
 ClassImp(TLondon1D3L)
 ClassImp(TLondon1D3LS)
+ClassImp(TLondon1D4L)
+ClassImp(TLondon1D3LSub)
 
 
 //------------------
@@ -82,6 +85,28 @@ TLondon1D3LS::~TLondon1D3LS() {
     fPofT = 0;
 }
 
+TLondon1D4L::~TLondon1D4L() {
+    fPar.clear();
+    fParForBofZ.clear();
+    fParForPofB.clear();
+    fParForPofT.clear();
+    delete fImpProfile;
+    fImpProfile = 0;
+    delete fPofT;
+    fPofT = 0;
+}
+
+TLondon1D3LSub::~TLondon1D3LSub() {
+    fPar.clear();
+    fParForBofZ.clear();
+    fParForPofB.clear();
+    fParForPofT.clear();
+    delete fImpProfile;
+    fImpProfile = 0;
+    delete fPofT;
+    fPofT = 0;
+}
+
 //------------------
 // Constructor of the TLondon1DHS class -- reading available implantation profiles and
 // creates (a pointer to) the TPofTCalc object (with the FFT plan)
@@ -143,6 +168,9 @@ TLondon1DHS::TLondon1DHS() : fCalcNeeded(true), fFirstCall(true) {
 //------------------
 
 double TLondon1DHS::operator()(double t, const vector<double> &par) const {
+
+  assert(par.size() == 5);
+
   if(t<0.0)
     return 0.0;
 
@@ -278,9 +306,11 @@ TLondon1D1L::TLondon1D1L() : fCalcNeeded(true), fFirstCall(true), fCallCounter(0
 
 double TLondon1D1L::operator()(double t, const vector<double> &par) const {
 
+  assert(par.size() == 6);
+
   // Debugging
   // Count the number of function calls
-  fCallCounter++;
+//  fCallCounter++;
 
   if(t<0.0)
     return 0.0;
@@ -350,15 +380,15 @@ double TLondon1D1L::operator()(double t, const vector<double> &par) const {
     fCalcNeeded = false;
   }
   
-  // Debugging start
+//   // Debugging start
-  if (!(fCallCounter%10000)){
+//   if (!(fCallCounter%10000)){
-    cout << fCallCounter-1 << "\t";
+//     cout << fCallCounter-1 << "\t";
-    for (unsigned int i(0); i<fPar.size(); i++){
+//     for (unsigned int i(0); i<fPar.size(); i++){
-      cout << fPar[i] << "\t";
+//       cout << fPar[i] << "\t";
-    }
+//     }
-    cout << endl;
+//     cout << endl;
-  }
+//   }
-  // Debugging end
+//   // Debugging end
 
   return fPofT->Eval(t);
 
@@ -423,6 +453,9 @@ TLondon1D2L::TLondon1D2L() : fCalcNeeded(true), fFirstCall(true), fLastTwoChange
 //------------------
 
 double TLondon1D2L::operator()(double t, const vector<double> &par) const {
+
+  assert(par.size() == 10);
+
   if(t<0.0)
     return 0.0;
 
@@ -568,6 +601,9 @@ TLondon1D3L::TLondon1D3L() : fCalcNeeded(true), fFirstCall(true), fLastThreeChan
 //------------------
 
 double TLondon1D3L::operator()(double t, const vector<double> &par) const {
+
+  assert(par.size() == 13);
+
   if(t<0.0)
     return 0.0;
 
@@ -728,6 +764,9 @@ TLondon1D3LS::TLondon1D3LS() : fCalcNeeded(true), fFirstCall(true), fLastThreeCh
 //------------------
 
 double TLondon1D3LS::operator()(double t, const vector<double> &par) const {
+
+  assert(par.size() == 12);
+
   if(t<0.0)
     return 0.0;
 
@@ -814,3 +853,337 @@ double TLondon1D3LS::operator()(double t, const vector<double> &par) const {
   return fPofT->Eval(t);
 
 }
+
+//------------------
+// Constructor of the TLondon1D4L class -- reading available implantation profiles and
+// creates (a pointer to) the TPofTCalc object (with the FFT plan)
+//------------------
+
+TLondon1D4L::TLondon1D4L() : fCalcNeeded(true), fFirstCall(true), fLastFourChanged(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;
+    }
+
+    fNSteps = startupHandler->GetNSteps();
+    fWisdom = startupHandler->GetWisdomFile();
+    string rge_path(startupHandler->GetDataPath());
+    vector<string> energy_vec(startupHandler->GetEnergyList());
+
+    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);
+
+    TTrimSPData *x = new TTrimSPData(rge_path, energy_vec);
+    fImpProfile = x;
+    x = 0;
+    delete x;
+
+    TPofTCalc *y = new TPofTCalc(fWisdom, fParForPofT);
+    fPofT = y;
+    y = 0;
+    delete y;
+
+    // clean up
+    if (saxParser) {
+      delete saxParser;
+      saxParser = 0;
+    }
+    if (startupHandler) {
+      delete startupHandler;
+      startupHandler = 0;
+    }
+}
+
+//------------------
+// TLondon1D4L-Method that calls the procedures to create B(z), 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 TLondon1D4L
+//------------------
+
+double TLondon1D4L::operator()(double t, const vector<double> &par) const {
+
+  assert(par.size() == 16);
+
+  if(t<0.0)
+    return 0.0;
+
+  // 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<fPar.size(); i++){
+      cout << "fPar[" << i << "] = " << fPar[i] << endl;
+    }
+*/
+    for (unsigned int i(2); i<fPar.size(); i++){
+      fParForBofZ.push_back(fPar[i]);
+//      cout << "fParForBofZ[" << i-2 << "] = " << fParForBofZ[i-2] << endl;
+    }
+    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 (i == fPar.size()-4 || i == fPar.size()-3 || i == fPar.size()-2 || i == fPar.size()-1)
+        fLastFourChanged = true;
+    }
+  }
+
+  if (par_changed)
+    fCalcNeeded = true;
+
+  // if model parameters have changed, recalculate B(z), 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(2); i<par.size(); i++)
+        fParForBofZ[i-2] = par[i];
+
+      fParForPofB[2] = par[1]; // energy
+
+/* DEBUG ---------------------------
+      for(unsigned int i(0); i<fParForBofZ.size(); i++) {
+        cout << "ParForBofZ[" << i << "] = " << fParForBofZ[i] << endl;
+      }
+
+      for(unsigned int i(0); i<fParForPofB.size(); i++) {
+        cout << "ParForPofB[" << i << "] = " << fParForPofB[i] << endl;
+      }
+
+      for(unsigned int i(0); i<fParForPofT.size(); i++) {
+        cout << "ParForPofT[" << i << "] = " << fParForPofT[i] << endl;
+      }
+------------------------------------*/
+
+      if(fLastFourChanged) {
+        vector<double> interfaces;
+        interfaces.push_back(par[3]+par[4]);
+        interfaces.push_back(par[3]+par[4]+par[5]);
+        interfaces.push_back(par[3]+par[4]+par[5]+par[6]);
+
+        vector<double> weights;
+        for(unsigned int i(par.size()-4); i<par.size(); i++)
+          weights.push_back(par[i]);
+
+//        cout << "Weighting has changed, re-calculating n(z) now..." << endl;
+        fImpProfile->WeightLayers(par[1], interfaces, weights);
+      }
+
+      TLondon1D_4L BofZ4(fNSteps, fParForBofZ);
+      TPofBCalc PofB4(BofZ4, *fImpProfile, fParForPofB);
+      fPofT->DoFFT(PofB4);
+
+    }/* else {
+      cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl;
+    }*/
+
+    fPofT->CalcPol(fParForPofT);
+
+    fCalcNeeded = false;
+    fLastFourChanged = false;
+  }
+
+  return fPofT->Eval(t);
+
+}
+
+//------------------
+// Constructor of the TLondon1D3LSub class -- reading available implantation profiles and
+// creates (a pointer to) the TPofTCalc object (with the FFT plan)
+//------------------
+
+TLondon1D3LSub::TLondon1D3LSub() : fCalcNeeded(true), fFirstCall(true), fWeightsChanged(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;
+    }
+
+    fNSteps = startupHandler->GetNSteps();
+    fWisdom = startupHandler->GetWisdomFile();
+    string rge_path(startupHandler->GetDataPath());
+    vector<string> energy_vec(startupHandler->GetEnergyList());
+
+    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);
+
+    TTrimSPData *x = new TTrimSPData(rge_path, energy_vec);
+    fImpProfile = x;
+    x = 0;
+    delete x;
+
+    TPofTCalc *y = new TPofTCalc(fWisdom, fParForPofT);
+    fPofT = y;
+    y = 0;
+    delete y;
+
+    // clean up
+    if (saxParser) {
+      delete saxParser;
+      saxParser = 0;
+    }
+    if (startupHandler) {
+      delete startupHandler;
+      startupHandler = 0;
+    }
+}
+
+//------------------
+// TLondon1D3L-Method that calls the procedures to create B(z), 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 TLondon1D3L
+//------------------
+
+double TLondon1D3LSub::operator()(double t, const vector<double> &par) const {
+
+  assert(par.size() == 15);
+
+  if(t<0.0)
+    return 0.0;
+
+  // 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<fPar.size(); i++){
+      cout << "fPar[" << i << "] = " << fPar[i] << endl;
+    }
+*/
+    for (unsigned int i(2); i<fPar.size(); i++){
+      fParForBofZ.push_back(fPar[i]);
+//      cout << "fParForBofZ[" << i-2 << "] = " << fParForBofZ[i-2] << endl;
+    }
+    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 (i == fPar.size()-5 || i == fPar.size()-5 || i == fPar.size()-3 || i == fPar.size()-2)
+        fWeightsChanged = true;
+    }
+  }
+
+  if (par_changed)
+    fCalcNeeded = true;
+
+  // if model parameters have changed, recalculate B(z), 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(2); i<par.size(); i++)
+        fParForBofZ[i-2] = par[i];
+
+      fParForPofB[2] = par[1]; // energy
+
+/* DEBUG ---------------------------
+      for(unsigned int i(0); i<fParForBofZ.size(); i++) {
+        cout << "ParForBofZ[" << i << "] = " << fParForBofZ[i] << endl;
+      }
+
+      for(unsigned int i(0); i<fParForPofB.size(); i++) {
+        cout << "ParForPofB[" << i << "] = " << fParForPofB[i] << endl;
+      }
+
+      for(unsigned int i(0); i<fParForPofT.size(); i++) {
+        cout << "ParForPofT[" << i << "] = " << fParForPofT[i] << endl;
+      }
+------------------------------------*/
+
+      vector<double> interfaces;
+      interfaces.push_back(par[3]+par[4]);
+      interfaces.push_back(par[3]+par[4]+par[5]);
+      interfaces.push_back(par[3]+par[4]+par[5]+par[6]);
+
+      if(fWeightsChanged) {
+        vector<double> weights;
+        for(unsigned int i(par.size()-5); i<(par.size()-1); i++)
+          weights.push_back(par[i]);
+
+//        cout << "Weighting has changed, re-calculating n(z) now..." << endl;
+        fImpProfile->WeightLayers(par[1], interfaces, weights);
+      }
+
+      TLondon1D_3L BofZ3(fNSteps, fParForBofZ);
+      TPofBCalc PofB3(BofZ3, *fImpProfile, fParForPofB);
+
+      // Add background contribution from the substrate
+      PofB3.AddBackground(par[2], par[14], fImpProfile->LayerFraction(par[1], 4, interfaces));
+
+      // FourierTransform of P(B)
+      fPofT->DoFFT(PofB3);
+
+    }/* else {
+      cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl;
+    }*/
+
+    fPofT->CalcPol(fParForPofT);
+
+    fCalcNeeded = false;
+    fWeightsChanged = false;
+  }
+
+  return fPofT->Eval(t);
+
+}
+

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

@@ -193,6 +193,42 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons
 
 }
 
+// TPofBCalc::AddBackground, Parameters: field B[G], width s[us], weight w
+
+void TPofBCalc::AddBackground(double B, double s, double w) {
+
+  if(!s || w<0.0 || w>1.0 || B<0.0)
+    return;
+
+  unsigned int sizePB(fPB.size());
+  double BsSq(s*s/gBar/gBar/4.0/pi/pi);
+
+  // calculate Gaussian background
+  vector<double> bg;
+  for(unsigned int i(0); i < sizePB; i++) {
+    bg.push_back(exp(-(fB[i]-B)*(fB[i]-B)/2.0/BsSq));
+  }
+
+  // normalize background
+  double bgsum(0.0);
+  for (unsigned int i(0); i < sizePB; i++)
+    bgsum += bg[i];
+  bgsum *= fDB;
+  for (unsigned int i(0); i < sizePB; i++)
+    bg[i] /= bgsum;
+
+  // add background to P(B)
+  for (unsigned int i(0); i < sizePB; i++)
+    fPB[i] = (1.0 - w)*fPB[i] + w*bg[i];
+
+//   // check if normalization is still valid
+//   double pBsum(0.0);
+//   for (unsigned int i(0); i < sizePB; i++)
+//     pBsum += fPB[i];
+//
+//   cout << "pBsum = " << pBsum << endl;
+
+}
 
 void TPofBCalc::ConvolveGss(double w) {
 

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

@@ -127,6 +127,52 @@ vector<double> TTrimSPData::OrigDataNZ(double e) const {
 
 }
 
+//---------------------
+// Method returning fraction of muons implanted in the specified layer for a given energy[keV]
+// Parameters: Energy[keV], LayerNumber[1], Interfaces[nm]
+//---------------------
+
+double TTrimSPData::LayerFraction(double e, unsigned int layno, const vector<double>& interface) const {
+
+  if(layno < 1 && layno > (interface.size()+1)) {
+    cout << "No such layer available according to your specified interfaces... Returning 0.0!" << endl;
+    return 0.0;
+  }
+
+  for(unsigned int i(0); i<fEnergy.size(); i++) {
+    if(!(fEnergy[i] - e)) {
+      // Because we do not know if the implantation profile is normalized or not, do not care about this and calculate the fraction from the beginning
+      // Total "number of muons"
+      double totalNumber(0.0);
+      for(unsigned int j(0); j<fDataZ[i].size(); j++)
+        totalNumber += fDataNZ[i][j];
+      // "number of muons" in layer layno
+      double layerNumber(0.0);
+      if(!(layno-1)){
+        for(unsigned int j(0); j<fDataZ[i].size(); j++)
+          if(fDataZ[i][j] < interface[0]*10.0)
+            layerNumber += fDataNZ[i][j];
+      } else if(!(layno-interface.size()-1)){
+        for(unsigned int j(0); j<fDataZ[i].size(); j++)
+          if(fDataZ[i][j] >= *(interface.end()-1)*10.0)
+            layerNumber += fDataNZ[i][j];
+      } else {
+        for(unsigned int j(0); j<fDataZ[i].size(); j++)
+          if(fDataZ[i][j] >= interface[layno-2]*10.0 && fDataZ[i][j] < interface[layno-1]*10.0)
+            layerNumber += fDataNZ[i][j];
+      }
+      // fraction of muons in layer layno
+      cout << "Fraction of muons in layer " << layno << ": " << layerNumber/totalNumber << endl;
+      return layerNumber/totalNumber;
+    }
+  }
+
+  // default
+  cout << "No implantation profile available for the specified energy... Returning 0.0" << endl;
+  return 0.0;
+
+}
+
 //---------------------
 // Method putting different weight to different layers of your thin film
 // Parameters: Implantation Energy[keV], Interfaces[nm], Weights [0.0 <= w[i] <= 1.0]