diff --git a/src/external/TFitPofB-lib/classes/TBofZCalc.cpp b/src/external/TFitPofB-lib/classes/TBofZCalc.cpp
index 5b70b3a2..55821d29 100644
--- a/src/external/TFitPofB-lib/classes/TBofZCalc.cpp
+++ b/src/external/TFitPofB-lib/classes/TBofZCalc.cpp
@@ -403,6 +403,182 @@ vector< pair<double, double> > TLondon1D_2L::GetInverseAndDerivative(double BB)
   return inv;
 }
 
+//------------------
+// Constructor of the TProximity1D_1LHS class
+// 1D-Proximity screening in a thin superconducting film, 1layer+bulk, linear + exponential decay
+// Parameters: Bext[G], B1[G], deadlayer[nm], thickness1[nm], lambda[nm]
+//------------------
+
+TProximity1D_1LHS::TProximity1D_1LHS(const vector<double> &param, unsigned int steps)
+{
+  fSteps = steps;
+  fDZ = 200./double(steps);
+  fMinTag = -1;
+  fMinZ = -1.0;
+  fMinB = -1.0;
+
+  // replaced the former assertions of positive lengths by fixed parameter limits
+
+  fParam.resize(param.size());
+  fParam[0] = param[0]; // Bext
+  (param[0] != param[1]) ? fParam[1] = param[1] : fParam[1] = param[1]-0.01; // B1
+  (param[2] >= 0.) ? fParam[2] = param[2] : fParam[2] = 0.; // deadlayer
+  (param[3] >= 0.) ? fParam[3] = param[3] : fParam[3] = 0.; // thickness1
+  (param[3] >= 1.) ? fParam[4] = param[4] : fParam[4] = 1.; // lambda
+
+  fInterfaces[0]=fParam[2];
+  fInterfaces[1]=fParam[2]+fParam[3];
+
+  SetBmin();
+}
+
+double TProximity1D_1LHS::GetBofZ(double ZZ) const
+{
+  if(ZZ < fInterfaces[0])
+    return fParam[0];
+
+  if (ZZ < fInterfaces[1]) {
+    return fParam[0]-(fParam[0]-fParam[1])/fParam[3]*(ZZ-fParam[2]);
+  } else {
+    return fParam[1]*exp((fParam[2]+fParam[3]-ZZ)/fParam[4]);
+  }
+}
+
+double TProximity1D_1LHS::GetBmax() const
+{
+  // return applied field
+  return fParam[0];
+}
+
+double TProximity1D_1LHS::GetBmin() const
+{
+  // return field minimum
+  return fMinB;
+}
+
+void TProximity1D_1LHS::SetBmin()
+{
+  fMinTag = 2;
+  fMinZ = 200.;
+  fMinB = GetBofZ(fMinZ);
+  return;
+}
+
+vector< pair<double, double> > TProximity1D_1LHS::GetInverseAndDerivative(double BB) const
+{
+  vector< pair<double, double> > inv;
+
+  if(BB <= fMinB || BB > fParam[0])
+    return inv;
+
+  double inverse[2];
+  pair<double, double> invAndDerivative;
+
+  inverse[0]=(fParam[0]*(fParam[2]+fParam[3])-fParam[1]*fParam[2]-BB*fParam[3])/(fParam[0]-fParam[1]);
+  inverse[1]=fParam[2]+fParam[3]-fParam[4]*log(BB/fParam[1]);
+
+  if(inverse[0] > fInterfaces[0] && inverse[0] <= fInterfaces[1]) {
+    invAndDerivative.first = inverse[0];
+    invAndDerivative.second = fParam[3]/(fParam[1]-fParam[0]);
+    inv.push_back(invAndDerivative);
+  }
+  if(inverse[1] > fInterfaces[1]) {
+    invAndDerivative.first = inverse[1];
+    invAndDerivative.second = -fParam[4]/BB;
+    inv.push_back(invAndDerivative);
+  }
+
+  return inv;
+}
+
+//------------------
+// Constructor of the TProximity1D_1LHSGss class
+// 1D-Proximity screening in a thin superconducting film, 1layer+bulk, linear + exponential decay
+// Parameters: Bext[G], B1[G], deadlayer[nm], thickness1[nm], sigma[nm]
+//------------------
+
+TProximity1D_1LHSGss::TProximity1D_1LHSGss(const vector<double> &param, unsigned int steps)
+{
+  fSteps = steps;
+  fDZ = 200./double(steps);
+  fMinTag = -1;
+  fMinZ = -1.0;
+  fMinB = -1.0;
+
+  // replaced the former assertions of positive lengths by fixed parameter limits
+
+  fParam.resize(param.size());
+  fParam[0] = param[0]; // Bext
+  (param[0] != param[1]) ? fParam[1] = param[1] : fParam[1] = param[1]-0.01; // B1
+  (param[2] >= 0.) ? fParam[2] = param[2] : fParam[2] = 0.; // deadlayer
+  (param[3] >= 0.) ? fParam[3] = param[3] : fParam[3] = 0.; // thickness1
+  (param[3] >= 1.) ? fParam[4] = param[4] : fParam[4] = 1.; // sigma
+
+  fInterfaces[0]=fParam[2];
+  fInterfaces[1]=fParam[2]+fParam[3];
+
+  SetBmin();
+}
+
+double TProximity1D_1LHSGss::GetBofZ(double ZZ) const
+{
+  if(ZZ < fInterfaces[0])
+    return fParam[0];
+
+  if (ZZ < fInterfaces[1]) {
+    return fParam[0]-(fParam[0]-fParam[1])/fParam[3]*(ZZ-fParam[2]);
+  } else {
+    return fParam[1]*exp(-(ZZ-fParam[2]-fParam[3])*(ZZ-fParam[2]-fParam[3])/(4.0*fParam[4]*fParam[4]));
+  }
+}
+
+double TProximity1D_1LHSGss::GetBmax() const
+{
+  // return applied field
+  return fParam[0];
+}
+
+double TProximity1D_1LHSGss::GetBmin() const
+{
+  // return field minimum
+  return fMinB;
+}
+
+void TProximity1D_1LHSGss::SetBmin()
+{
+  fMinTag = 2;
+  fMinZ = 200.;
+  fMinB = GetBofZ(fMinZ);
+  return;
+}
+
+vector< pair<double, double> > TProximity1D_1LHSGss::GetInverseAndDerivative(double BB) const
+{
+  vector< pair<double, double> > inv;
+
+  if(BB <= fMinB || BB > fParam[0])
+    return inv;
+
+  double inverse[2];
+  pair<double, double> invAndDerivative;
+
+  inverse[0]=(fParam[0]*(fParam[2]+fParam[3])-fParam[1]*fParam[2]-BB*fParam[3])/(fParam[0]-fParam[1]);
+  inverse[1]=fParam[2]+fParam[3]+fParam[4]*sqrt(2.0*log(fParam[1]/BB));
+
+  if(inverse[0] > fInterfaces[0] && inverse[0] <= fInterfaces[1]) {
+    invAndDerivative.first = inverse[0];
+    invAndDerivative.second = fParam[3]/(fParam[1]-fParam[0]);
+    inv.push_back(invAndDerivative);
+  }
+  if(inverse[1] > fInterfaces[1]) {
+    invAndDerivative.first = inverse[1];
+    invAndDerivative.second = -fParam[4]/(BB*sqrt(2.0*log(fParam[1]/BB)));
+    inv.push_back(invAndDerivative);
+  }
+
+  return inv;
+}
+
 //------------------
 // Constructor of the TLondon1D_3L class
 // 1D-London screening in a thin superconducting film, three layers, three lambdas
diff --git a/src/external/TFitPofB-lib/classes/TLondon1D.cpp b/src/external/TFitPofB-lib/classes/TLondon1D.cpp
index 988c629d..030f4168 100644
--- a/src/external/TFitPofB-lib/classes/TLondon1D.cpp
+++ b/src/external/TFitPofB-lib/classes/TLondon1D.cpp
@@ -12,6 +12,7 @@
 #include "TLondon1D.h"
 #include <iostream>
 #include <cassert>
+#include <cmath>
 using namespace std;
 
 #include <TSAXParser.h>
@@ -20,6 +21,8 @@ using namespace std;
 ClassImp(TLondon1DHS)
 ClassImp(TLondon1D1L)
 ClassImp(TLondon1D2L)
+ClassImp(TProximity1D1LHS)
+ClassImp(TProximity1D1LHSGss)
 ClassImp(TLondon1D3L)
 ClassImp(TLondon1D3LS)
 // ClassImp(TLondon1D4L)
@@ -66,6 +69,28 @@ TLondon1D2L::~TLondon1D2L() {
     fPofT = 0;
 }
 
+TProximity1D1LHS::~TProximity1D1LHS() {
+    fPar.clear();
+    fParForBofZ.clear();
+    fParForPofB.clear();
+    fParForPofT.clear();
+    delete fImpProfile;
+    fImpProfile = 0;
+    delete fPofT;
+    fPofT = 0;
+}
+
+TProximity1D1LHSGss::~TProximity1D1LHSGss() {
+    fPar.clear();
+    fParForBofZ.clear();
+    fParForPofB.clear();
+    fParForPofT.clear();
+    delete fImpProfile;
+    fImpProfile = 0;
+    delete fPofT;
+    fPofT = 0;
+}
+
 TLondon1D3L::~TLondon1D3L() {
     fPar.clear();
     fParForBofZ.clear();
@@ -175,7 +200,7 @@ double TLondon1DHS::operator()(double t, const vector<double> &par) const {
   assert(par.size() == 5);
 
   if(t<0.0)
-    return 0.0;
+    return cos(par[0]*0.017453293);
 
   // check if the function is called the first time and if yes, read in parameters
 
@@ -316,7 +341,7 @@ double TLondon1D1L::operator()(double t, const vector<double> &par) const {
 //  fCallCounter++;
 
   if(t<0.0)
-    return 0.0;
+    return cos(par[0]*0.017453293);
 
   // check if the function is called the first time and if yes, read in parameters
 
@@ -460,7 +485,7 @@ double TLondon1D2L::operator()(double t, const vector<double> &par) const {
   assert(par.size() == 10);
 
   if(t<0.0)
-    return 0.0;
+    return cos(par[0]*0.017453293);
 
   // check if the function is called the first time and if yes, read in parameters
 
@@ -545,6 +570,285 @@ double TLondon1D2L::operator()(double t, const vector<double> &par) const {
 
 }
 
+//------------------
+// Constructor of the TProximity1D1LHS class -- reading available implantation profiles and
+// creates (a pointer to) the TPofTCalc object (with the FFT plan)
+//------------------
+
+TProximity1D1LHS::TProximity1D1LHS() : 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;
+    }
+
+    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);
+//    fParForPofB.push_back(0.0);
+
+    TTrimSPData *x = new TTrimSPData(rge_path, energy_vec);
+    fImpProfile = x;
+    x = 0;
+
+    TPofTCalc *y = new TPofTCalc(fWisdom, fParForPofT);
+    fPofT = y;
+    y = 0;
+
+    // clean up
+    if (saxParser) {
+      delete saxParser;
+      saxParser = 0;
+    }
+    if (startupHandler) {
+      delete startupHandler;
+      startupHandler = 0;
+    }
+}
+
+//------------------
+// TProximity1D1LHS-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 TProximity1D1LHS
+//------------------
+
+double TProximity1D1LHS::operator()(double t, const vector<double> &par) const {
+
+  assert(par.size() == 8);
+
+  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
+
+  bool width_changed(false);
+
+  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;
+    width_changed = true;
+//  cout << this << endl;
+  }
+
+  // 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 == 7){
+          width_changed = 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<fPar.size(); i++)
+        fParForBofZ[i-2] = par[i];
+
+      fParForPofB[2] = par[1]; // energy
+
+      if(width_changed) { // Convolution of the implantation profile with Gaussian
+        fImpProfile->ConvolveGss(par[7], par[1]);
+        width_changed = false;
+      }
+
+      TProximity1D_1LHS BofZ(fParForBofZ);
+      TPofBCalc PofB(BofZ, *fImpProfile, fParForPofB);
+      fPofT->DoFFT(PofB);
+
+    }/* else {
+      cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl;
+    }*/
+
+    fPofT->CalcPol(fParForPofT);
+
+    fCalcNeeded = false;
+  }
+
+  return fPofT->Eval(t);
+
+}
+
+//------------------
+// Constructor of the TProximity1D1LHSGss class -- reading available implantation profiles and
+// creates (a pointer to) the TPofTCalc object (with the FFT plan)
+//------------------
+
+TProximity1D1LHSGss::TProximity1D1LHSGss() : 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;
+    }
+
+    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);
+//    fParForPofB.push_back(0.0);
+
+    TTrimSPData *x = new TTrimSPData(rge_path, energy_vec);
+    fImpProfile = x;
+    x = 0;
+
+    TPofTCalc *y = new TPofTCalc(fWisdom, fParForPofT);
+    fPofT = y;
+    y = 0;
+
+    // clean up
+    if (saxParser) {
+      delete saxParser;
+      saxParser = 0;
+    }
+    if (startupHandler) {
+      delete startupHandler;
+      startupHandler = 0;
+    }
+}
+
+//------------------
+// TProximity1D1LHS-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 TProximity1D1LHS
+//------------------
+
+double TProximity1D1LHSGss::operator()(double t, const vector<double> &par) const {
+
+  assert(par.size() == 7);
+
+  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<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;
+//  cout << this << endl;
+  }
+
+  // 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(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<fPar.size(); i++)
+        fParForBofZ[i-2] = par[i];
+
+      fParForPofB[2] = par[1]; // energy
+
+      TProximity1D_1LHSGss BofZ(fParForBofZ);
+      TPofBCalc PofB(BofZ, *fImpProfile, fParForPofB);
+      fPofT->DoFFT(PofB);
+
+    }/* else {
+      cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl;
+    }*/
+
+    fPofT->CalcPol(fParForPofT);
+
+    fCalcNeeded = false;
+  }
+
+  return fPofT->Eval(t);
+
+}
+
 //------------------
 // Constructor of the TLondon1D3L class -- reading available implantation profiles and
 // creates (a pointer to) the TPofTCalc object (with the FFT plan)
@@ -608,7 +912,7 @@ double TLondon1D3L::operator()(double t, const vector<double> &par) const {
   assert(par.size() == 13);
 
   if(t<0.0)
-    return 0.0;
+    return cos(par[0]*0.017453293);
 
   // check if the function is called the first time and if yes, read in parameters
 
@@ -771,7 +1075,7 @@ double TLondon1D3LS::operator()(double t, const vector<double> &par) const {
   assert(par.size() == 12);
 
   if(t<0.0)
-    return 0.0;
+    return cos(par[0]*0.017453293);
 
   // check if the function is called the first time and if yes, read in parameters
 
@@ -1088,7 +1392,7 @@ double TLondon1D3LSub::operator()(double t, const vector<double> &par) const {
   assert(par.size() == 15);
 
   if(t<0.0)
-    return 0.0;
+    return cos(par[0]*0.017453293);
 
   // check if the function is called the first time and if yes, read in parameters
 
diff --git a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp
index 34618f83..e9bc0520 100644
--- a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp
+++ b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp
@@ -372,12 +372,12 @@ void TPofBCalc::AddBackground(double B, double s, double w) {
     return;
 
   unsigned int sizePB(fPB.size());
-  double BsSq(s*s/gBar/gBar/4.0/pi/pi);
+  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));
+    bg.push_back(exp(-(fB[i]-B)*(fB[i]-B)/(2.0*BsSq)));
   }
 
   // normalize background
diff --git a/src/external/TFitPofB-lib/classes/TPofTCalc.cpp b/src/external/TFitPofB-lib/classes/TPofTCalc.cpp
index 72ada3ef..8dc5dbf5 100644
--- a/src/external/TFitPofB-lib/classes/TPofTCalc.cpp
+++ b/src/external/TFitPofB-lib/classes/TPofTCalc.cpp
@@ -193,7 +193,8 @@ void TPofTCalc::FakeData(const string &rootOutputFileName, const vector<double>
 
   cout << "TPofTCalc::FakeData: " << numHist << " histograms to be built" << endl;
 
-  vector<unsigned int> t0;
+  int nChannels = int(par[3]);
+  vector<int> t0;
   vector<double> asy0;
   vector<double> phase0;
   vector<double> N0;
@@ -213,44 +214,51 @@ void TPofTCalc::FakeData(const string &rootOutputFileName, const vector<double>
   param.push_back(par[1]); // dB
 
   vector< vector<double> > asy;
-  vector<double> asydata;
+  vector<double> asydata(nChannels);
   double ttime(0.0);
-  double pol(0.0);
+  int j,k;
 
   for(unsigned int i(0); i<numHist; i++) {
     param[0]=phase0[i];
     // calculate asymmetry
     CalcPol(param);
-    for(unsigned int j(0); j<par[3]; j++){
+
+#pragma omp parallel for default(shared) private(j,ttime,k) schedule(dynamic)
+    for(j=0; j<nChannels; j++){
       ttime=j*par[2];
-      for(unsigned int k(0); k<fT.size()-1; k++){
-        if (ttime < fT[k+1]) {
-          pol=fPT[k]+(fPT[k+1]-fPT[k])/(fT[k+1]-fT[k])*(ttime-fT[k]);
-          asydata.push_back(asy0[i]*pol);
-          break;
-        }
-      }
+      k = floor(ttime/fTBin);
+      asydata[j]=asy0[i]*(fPT[k]+(fPT[k+1]-fPT[k])/fTBin*(ttime-fT[k]));
     }
+// end omp
+
+//       for(unsigned int k(0); k<fT.size()-1; k++){
+//         if (ttime < fT[k+1]) {
+//           pol=fPT[k]+(fPT[k+1]-fPT[k])/(fT[k+1]-fT[k])*(ttime-fT[k]);
+//           asydata.push_back(asy0[i]*pol);
+//           break;
+//         }
+//       }
     asy.push_back(asydata);
-    asydata.clear();
+//    asydata.clear();
     cout << "TPofTCalc::FakeData: " << i+1 << "/" << numHist << " calculated!" << endl;
   }
 
   // calculate the histograms
   vector< vector<double> > histo;
-  vector<double> data;
-  double value;
+  vector<double> data(nChannels);
 
   for (unsigned int i(0); i<numHist; i++) {    // loop over all histos
-    for (unsigned int j(0); j<par[3]; j++) { // loop over time
+
+#pragma omp parallel for default(shared) private(j) schedule(dynamic)
+    for (j = 0; j<nChannels; j++) { // loop over time
       if (j < t0[i]) // j<t0
-        value = bg[i]; // background
+        data[j] = bg[i]; // background
       else
-        value = N0[i]*exp(-par[2]*double(j-t0[i])/tauMu)*(1.0+asy[i][j-t0[i]])+bg[i];
-      data.push_back(value);
+        data[j] = N0[i]*exp(-par[2]*double(j-t0[i])/tauMu)*(1.0+asy[i][j-t0[i]])+bg[i];
     }
+// end omp
+
     histo.push_back(data);
-    data.clear();
     cout << "TPofTCalc::FakeData: " << i+1 << "/" << numHist << " done ...";
   }
 
@@ -275,8 +283,11 @@ void TPofTCalc::FakeData(const string &rootOutputFileName, const vector<double>
     name  += i;
     fakeHisto = new TH1F(name.Data(), name.Data(), int(par[3]), -par[2]/2.0, (par[3]+0.5)*par[2]);
     // fill theoHisto
-    for (unsigned int j(0); j<par[3]; j++)
+#pragma omp parallel for default(shared) private(j) schedule(dynamic)
+    for (j = 0; j<nChannels; j++)
       theoHisto->SetBinContent(j, histo[i][j]);
+// end omp
+
     // fill fakeHisto
     fakeHisto->FillRandom(theoHisto, (int)theoHisto->Integral());
 
@@ -298,14 +309,16 @@ void TPofTCalc::FakeData(const string &rootOutputFileName, const vector<double>
   runHeader->SetRunTitle("Fake Data");
   float fval = par[2]*1000.; //us->ns
   runHeader->SetTimeResolution(fval);
-  runHeader->SetNChannels(int(par[3]));
+  runHeader->SetNChannels(nChannels);
   runHeader->SetNHist(histoData.size());
   double *t0array = new double[histoData.size()];
   for (unsigned int i(0); i<histoData.size(); i++)
     t0array[i] = t0[i];
   runHeader->SetTimeZero(t0array);
-  if (t0array)
+  if (t0array) {
     delete t0array;
+    t0array = 0;
+  }
   runInfoFolder->Add(runHeader);
 
   // create decay histo folder and content
@@ -348,6 +361,13 @@ void TPofTCalc::FakeData(const string &rootOutputFileName, const vector<double>
   }
   histoData.clear();
   histoDataPPC.clear();
+  fakeHisto = 0;
+
+  delete histoFolder; histoFolder = 0;
+  delete decayAnaModule; decayAnaModule = 0;
+
+  delete runInfoFolder; runInfoFolder = 0;
+  delete runHeader; runHeader = 0;
 
   t0.clear();
   asy0.clear();
diff --git a/src/external/TFitPofB-lib/classes/TTrimSPDataHandler.cpp b/src/external/TFitPofB-lib/classes/TTrimSPDataHandler.cpp
index ecc90653..a707ecfb 100644
--- a/src/external/TFitPofB-lib/classes/TTrimSPDataHandler.cpp
+++ b/src/external/TFitPofB-lib/classes/TTrimSPDataHandler.cpp
@@ -434,7 +434,7 @@ double TTrimSPData::PeakRange(double e) const {
 
     if(nziter != fDataNZ[i].end()){
       unsigned int j(nziter - fDataNZ[i].begin());
-      return fDataZ[i][j];
+      return fDataZ[i][j]/10.0;
     }
     cout << "TTrimSPData::PeakRange: No maximum found in the implantation profile... Returning -1! Please check the profile!" << endl;
     return -1.;
diff --git a/src/external/TFitPofB-lib/include/TBofZCalc.h b/src/external/TFitPofB-lib/include/TBofZCalc.h
index 60a6f98d..edb2fc3d 100644
--- a/src/external/TFitPofB-lib/include/TBofZCalc.h
+++ b/src/external/TFitPofB-lib/include/TBofZCalc.h
@@ -124,6 +124,53 @@ private:
   double fCoeff[4];
 };
 
+
+//--------------------
+// Class "for proximity screening" in a thin superconducting film - one layer + superconducting half space
+//--------------------
+
+class TProximity1D_1LHS : public TBofZCalcInverse {
+
+public:
+
+  TProximity1D_1LHS(const vector<double>&, unsigned int steps = 3000);
+  double GetBofZ(double) const;
+  double GetBmin() const;
+  double GetBmax() const;
+  vector< pair<double, double> > GetInverseAndDerivative(double) const;
+
+private:
+  void SetBmin();
+
+  int fMinTag;
+  double fMinZ;
+  double fMinB;
+  double fInterfaces[2];
+};
+
+//--------------------
+// Class "for proximity screening" in a thin superconducting film - one layer + superconducting half space
+//--------------------
+
+class TProximity1D_1LHSGss : public TBofZCalcInverse {
+
+public:
+
+  TProximity1D_1LHSGss(const vector<double>&, unsigned int steps = 3000);
+  double GetBofZ(double) const;
+  double GetBmin() const;
+  double GetBmax() const;
+  vector< pair<double, double> > GetInverseAndDerivative(double) const;
+
+private:
+  void SetBmin();
+
+  int fMinTag;
+  double fMinZ;
+  double fMinB;
+  double fInterfaces[2];
+};
+
 //--------------------
 // Class "for Meissner screening" in a thin superconducting film - tri-layer with three different lambdas
 //--------------------
diff --git a/src/external/TFitPofB-lib/include/TLondon1D.h b/src/external/TFitPofB-lib/include/TLondon1D.h
index 9faf330b..014281d3 100644
--- a/src/external/TFitPofB-lib/include/TLondon1D.h
+++ b/src/external/TFitPofB-lib/include/TLondon1D.h
@@ -89,6 +89,55 @@ private:
   ClassDef(TLondon1D2L,1)
 };
 
+class TProximity1D1LHS : public PUserFcnBase {
+
+public:
+  // default constructor
+  TProximity1D1LHS();
+  ~TProximity1D1LHS();
+
+  double operator()(double, const vector<double>&) const;
+
+private:
+  mutable vector<double> fPar;
+  TTrimSPData *fImpProfile;
+  TPofTCalc *fPofT;
+  mutable bool fCalcNeeded;
+  mutable bool fFirstCall;
+  mutable vector<double> fParForPofT;
+  mutable vector<double> fParForBofZ;
+  mutable vector<double> fParForPofB;
+  string fWisdom;
+  unsigned int fNSteps;
+
+  ClassDef(TProximity1D1LHS,1)
+};
+
+class TProximity1D1LHSGss : public PUserFcnBase {
+
+public:
+  // default constructor
+  TProximity1D1LHSGss();
+  ~TProximity1D1LHSGss();
+
+  double operator()(double, const vector<double>&) const;
+
+private:
+  mutable vector<double> fPar;
+  TTrimSPData *fImpProfile;
+  TPofTCalc *fPofT;
+  mutable bool fCalcNeeded;
+  mutable bool fFirstCall;
+  mutable vector<double> fParForPofT;
+  mutable vector<double> fParForBofZ;
+  mutable vector<double> fParForPofB;
+  string fWisdom;
+  unsigned int fNSteps;
+
+  ClassDef(TProximity1D1LHSGss,1)
+};
+
+
 class TLondon1D3L : public PUserFcnBase {
 
 public:
diff --git a/src/external/TFitPofB-lib/include/TLondon1DLinkDef.h b/src/external/TFitPofB-lib/include/TLondon1DLinkDef.h
index 8b723dac..f8ff47f6 100644
--- a/src/external/TFitPofB-lib/include/TLondon1DLinkDef.h
+++ b/src/external/TFitPofB-lib/include/TLondon1DLinkDef.h
@@ -19,6 +19,8 @@
 #pragma link C++ class TLondon1DHS+;
 #pragma link C++ class TLondon1D1L+;
 #pragma link C++ class TLondon1D2L+;
+#pragma link C++ class TProximity1D1LHS+;
+#pragma link C++ class TProximity1D1LHSGss+;
 #pragma link C++ class TLondon1D3L+;
 #pragma link C++ class TLondon1D3LS+;
 //#pragma link C++ class TLondon1D4L+;
diff --git a/src/external/TFitPofB-lib/test/test.cpp b/src/external/TFitPofB-lib/test/test.cpp
index 5f378f2b..7f85b1f4 100644
--- a/src/external/TFitPofB-lib/test/test.cpp
+++ b/src/external/TFitPofB-lib/test/test.cpp
@@ -67,7 +67,7 @@ int main(){
   }
   
   char debugfile1[50];
-  int nn = sprintf (debugfile1, "4Ltest-Bz_z-%.4f_l-%.3f_E-%.1f_normal.dat", par_vec[2], par_vec[11], par_vec[1]);
+  int nn = sprintf (debugfile1, "4Ltest-Bz_z-%.4f_l-%.3f_E-%.1f_normal.dat", par_vec[2], par_vec[9], par_vec[1]);
   if (nn > 0) {
     ofstream of01(debugfile1);
     for (unsigned int i(0); i<2000; i++) {
@@ -77,15 +77,16 @@ int main(){
   }
 
   TPofBCalc PofB(BofZ, calcData, parForPofB);
-
-  PofB.AddBackground(par_vec[2], 0.01, calcData.LayerFraction(E, 4, interfaces));
 //  PofB.ConvolveGss(1.17);
 
+  PofB.AddBackground(par_vec[2], 0.2, calcData.LayerFraction(E, 4, interfaces));
+  PofB.ConvolveGss(1.17);
+
   vector<double> hurgaB(PofB.DataB());
   vector<double> hurgaPB(PofB.DataPB());
   
   char debugfile[50];
-  int n = sprintf (debugfile, "4Ltest-BpB_B-%.4f_l-%.3f_E-%.1f_normal.dat", par_vec[2], par_vec[11], par_vec[1]);
+  int n = sprintf (debugfile, "4Ltest-BpB_B-%.4f_l-%.3f_E-%.1f_normal.dat", par_vec[2], par_vec[9], par_vec[1]);
 
   if (n > 0) {
     ofstream of7(debugfile);
@@ -101,7 +102,7 @@ int main(){
   poft.CalcPol(parForPofT);
   
   char debugfilex[50];
-  int nx = sprintf (debugfilex, "4Ltest-P_t-%.4f_l-%.3f_E-%.1f_normal.dat", par_vec[2], par_vec[11], par_vec[1]);
+  int nx = sprintf (debugfilex, "4Ltest-P_t-%.4f_l-%.3f_E-%.1f_normal.dat", par_vec[2], par_vec[9], par_vec[1]);
   
   if (nx > 0) { 
     ofstream of8(debugfilex);