diff --git a/src/external/TFitPofB-lib/classes/TLondon1D.cpp b/src/external/TFitPofB-lib/classes/TLondon1D.cpp
index a86882e1..8350716a 100644
--- a/src/external/TFitPofB-lib/classes/TLondon1D.cpp
+++ b/src/external/TFitPofB-lib/classes/TLondon1D.cpp
@@ -46,7 +46,6 @@ ClassImp(TProximity1D1LHSGss)
 ClassImp(TLondon1D3L)
 ClassImp(TLondon1D3LS)
 // ClassImp(TLondon1D4L)
-ClassImp(TLondon1D3LSub)
 
 ClassImp(TLondon1D3Lestimate)
 
@@ -158,19 +157,6 @@ TLondon1D3LS::~TLondon1D3LS() {
 //     fPofT = 0;
 // }
 
-TLondon1D3LSub::~TLondon1D3LSub() {
-    fPar.clear();
-    fParForBofZ.clear();
-    fParForPofB.clear();
-    fParForPofT.clear();
-    delete fImpProfile;
-    fImpProfile = 0;
-    delete fPofB;
-    fPofB = 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)
@@ -234,13 +220,12 @@ TLondon1DHS::TLondon1DHS() : fCalcNeeded(true), fFirstCall(true) {
 
 double TLondon1DHS::operator()(double t, const vector<double> &par) const {
 
-  assert(par.size() == 5 || par.size() == 6);
+  assert(par.size() == 5);
 
   if(t<0.0)
     return cos(par[0]*0.017453293);
 
   bool dead_layer_changed(false);
-  bool width_changed(false);
 
   // check if the function is called the first time and if yes, read in parameters
 
@@ -252,8 +237,6 @@ double TLondon1DHS::operator()(double t, const vector<double> &par) const {
     }
     fFirstCall = false;
     dead_layer_changed = true;
-    if(par.size() == 6)
-      width_changed = true;
   }
 
   // check if any parameter has changed
@@ -271,8 +254,6 @@ double TLondon1DHS::operator()(double t, const vector<double> &par) const {
         only_phase_changed = false;
         if (i == 3) {
           dead_layer_changed = true;
-        } else if (i == 5) {
-          width_changed = true;
         }
       }
     }
@@ -294,10 +275,6 @@ double TLondon1DHS::operator()(double t, const vector<double> &par) const {
       for (unsigned int i(2); i<fPar.size(); i++)
         fParForBofZ[i-2] = par[i];
 
-      if(width_changed) { // Convolution of the implantation profile with Gaussian
-        fImpProfile->ConvolveGss(par[5], par[1]);
-      }
-
       fParForPofB[2] = par[1]; // energy
       fParForPofB[3] = par[2]; // Bkg-Field
       //fParForPofB[4] = 0.005; // Bkg-width (in principle zero)
@@ -361,7 +338,10 @@ TLondon1D1L::TLondon1D1L() : fCalcNeeded(true), fFirstCall(true), fCallCounter(0
 
     fParForPofB.push_back(startupHandler->GetDeltat());
     fParForPofB.push_back(startupHandler->GetDeltaB());
-    fParForPofB.push_back(0.0);
+    fParForPofB.push_back(0.0); // Energy
+    fParForPofB.push_back(0.0); // Bkg-Field
+    fParForPofB.push_back(0.005); // Bkg-width
+    fParForPofB.push_back(0.0); // Bkg-weight
 
     fImpProfile = new TTrimSPData(rge_path, energy_vec);
 
@@ -388,30 +368,31 @@ TLondon1D1L::TLondon1D1L() : fCalcNeeded(true), fFirstCall(true), fCallCounter(0
 
 double TLondon1D1L::operator()(double t, const vector<double> &par) const {
 
-  assert(par.size() == 6);
+  assert(par.size() == 6 || par.size() == 8);
 
   // Debugging
   // Count the number of function calls
-//  fCallCounter++;
+  //  fCallCounter++;
 
   if(t<0.0)
     return cos(par[0]*0.017453293);
 
+  bool bkg_fraction_changed(false);
+  bool weights_changed(false);
+
   // 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;
+
+    fFirstCall = false;
+    bkg_fraction_changed = true;
+    if (par.size() == 8)
+      weights_changed = true;
   }
 
   // check if any parameter has changed
@@ -427,6 +408,10 @@ double TLondon1D1L::operator()(double t, const vector<double> &par) const {
         only_phase_changed = true;
       } else {
         only_phase_changed = false;
+        if (i == 3 || i == 4)
+          bkg_fraction_changed = true;
+        if (i == 6 || i == 7)
+          weights_changed = true;
       }
     }
   }
@@ -448,6 +433,33 @@ double TLondon1D1L::operator()(double t, const vector<double> &par) const {
         fParForBofZ[i-2] = par[i];
 
       fParForPofB[2] = par[1]; // energy
+      fParForPofB[3] = par[2]; // Bkg-Field
+      //fParForPofB[4] = 0.005; // Bkg-width (in principle zero)
+
+      if(weights_changed) {
+        vector<double> interfaces;
+        interfaces.push_back(par[3]+par[4]);
+
+        vector<double> weights;
+        for(unsigned int i(6); i<8; i++)
+          weights.push_back(par[i]);
+
+//        cout << "Weighting has changed, re-calculating n(z) now..." << endl;
+        fImpProfile->WeightLayers(par[1], interfaces, weights);
+
+        interfaces.clear();
+        weights.clear();
+      }
+
+      if(bkg_fraction_changed || weights_changed) {
+        vector<double> interfaces;
+        interfaces.push_back(par[3]);// dead layer
+        interfaces.push_back(par[3] + par[4]);// dead layer + first layer
+        fParForPofB[5] = fImpProfile->LayerFraction(par[1], 1, interfaces) + // Fraction of muons in the deadlayer
+                         fImpProfile->LayerFraction(par[1], 3, interfaces); // Fraction of muons in the substrate
+        interfaces.clear();
+      }
+
 
       TLondon1D_1L BofZ1(fParForBofZ);
       fPofB->UnsetPBExists();
@@ -482,7 +494,7 @@ double TLondon1D1L::operator()(double t, const vector<double> &par) const {
 // creates (a pointer to) the TPofTCalc object (with the FFT plan)
 //------------------
 
-TLondon1D2L::TLondon1D2L() : fCalcNeeded(true), fFirstCall(true), fLastTwoChanged(true) {
+TLondon1D2L::TLondon1D2L() : fCalcNeeded(true), fFirstCall(true) {
   // read startup file
     string startup_path_name("TFitPofB_startup.xml");
 
@@ -509,7 +521,10 @@ TLondon1D2L::TLondon1D2L() : fCalcNeeded(true), fFirstCall(true), fLastTwoChange
 
     fParForPofB.push_back(startupHandler->GetDeltat());
     fParForPofB.push_back(startupHandler->GetDeltaB());
-    fParForPofB.push_back(0.0);
+    fParForPofB.push_back(0.0); // Energy
+    fParForPofB.push_back(0.0); // Bkg-Field
+    fParForPofB.push_back(0.005); // Bkg-width
+    fParForPofB.push_back(0.0); // Bkg-weight
 
     fImpProfile = new TTrimSPData(rge_path, energy_vec);
 
@@ -536,25 +551,26 @@ TLondon1D2L::TLondon1D2L() : fCalcNeeded(true), fFirstCall(true), fLastTwoChange
 
 double TLondon1D2L::operator()(double t, const vector<double> &par) const {
 
-  assert(par.size() == 10);
+  assert(par.size() == 8 || par.size() == 11);
 
   if(t<0.0)
     return cos(par[0]*0.017453293);
 
+  bool bkg_fraction_changed(false);
+  bool weights_changed(false);
+
   // 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;
+    bkg_fraction_changed = true;
+    if (par.size() == 11)
+      weights_changed = true;
   }
 
   // check if any parameter has changed
@@ -570,9 +586,11 @@ double TLondon1D2L::operator()(double t, const vector<double> &par) const {
         only_phase_changed = true;
       } else {
         only_phase_changed = false;
+        if (i == 3 || i == 4 || i == 5)
+          bkg_fraction_changed = true;
+        if (i == 8 || i == 9 || i == 10)
+          weights_changed = true;
       }
-      if (i == fPar.size()-2 || i == fPar.size()-1)
-        fLastTwoChanged = true;
     }
   }
 
@@ -593,17 +611,32 @@ double TLondon1D2L::operator()(double t, const vector<double> &par) const {
         fParForBofZ[i-2] = par[i];
 
       fParForPofB[2] = par[1]; // energy
+      fParForPofB[3] = par[2]; // Bkg-Field
+      //fParForPofB[4] = 0.005; // Bkg-width (in principle zero)
 
-      if(fLastTwoChanged) {
+      if(weights_changed) {
         vector<double> interfaces;
         interfaces.push_back(par[3]+par[4]);
+        interfaces.push_back(par[3]+par[4]+par[5]);
 
         vector<double> weights;
-        for(unsigned int i(par.size()-2); i<par.size(); i++)
+        for(unsigned int i(8); i<11; i++)
           weights.push_back(par[i]);
 
 //        cout << "Weighting has changed, re-calculating n(z) now..." << endl;
         fImpProfile->WeightLayers(par[1], interfaces, weights);
+
+        interfaces.clear();
+        weights.clear();
+      }
+
+      if(bkg_fraction_changed || weights_changed) {
+        vector<double> interfaces;
+        interfaces.push_back(par[3]);// dead layer
+        interfaces.push_back(par[3] + par[4] + par[5]);// dead layer + first layer + second layer
+        fParForPofB[5] = fImpProfile->LayerFraction(par[1], 1, interfaces) + // Fraction of muons in the deadlayer
+                         fImpProfile->LayerFraction(par[1], 3, interfaces); // Fraction of muons in the substrate
+        interfaces.clear();
       }
 
       TLondon1D_2L BofZ2(fParForBofZ);
@@ -619,7 +652,6 @@ double TLondon1D2L::operator()(double t, const vector<double> &par) const {
     fPofT->CalcPol(fParForPofT);
 
     fCalcNeeded = false;
-    fLastTwoChanged = false;
   }
 
   return fPofT->Eval(t);
@@ -659,7 +691,7 @@ TProximity1D1LHS::TProximity1D1LHS() : fCalcNeeded(true), fFirstCall(true) {
 
     fParForPofB.push_back(startupHandler->GetDeltat());
     fParForPofB.push_back(startupHandler->GetDeltaB());
-    fParForPofB.push_back(0.0);
+    fParForPofB.push_back(0.0); // Energy
     fParForPofB.push_back(0.0); // Bkg-Field
     fParForPofB.push_back(0.01); // Bkg-width
     fParForPofB.push_back(0.0); // Bkg-weight
@@ -689,31 +721,23 @@ TProximity1D1LHS::TProximity1D1LHS() : fCalcNeeded(true), fFirstCall(true) {
 
 double TProximity1D1LHS::operator()(double t, const vector<double> &par) const {
 
-  assert(par.size() == 8);
+  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
 
-  bool width_changed(false);
   bool dead_layer_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;
     dead_layer_changed = true;
-//  cout << this << endl;
   }
 
   // check if any parameter has changed
@@ -729,9 +753,6 @@ double TProximity1D1LHS::operator()(double t, const vector<double> &par) const {
         only_phase_changed = true;
       } else {
         only_phase_changed = false;
-        if (i == 7){
-          width_changed = true;
-        }
         if (i == 4){
           dead_layer_changed = true;
         }
@@ -755,10 +776,6 @@ double TProximity1D1LHS::operator()(double t, const vector<double> &par) const {
       for (unsigned int i(2); i<fPar.size(); i++)
         fParForBofZ[i-2] = par[i];
 
-      if(width_changed) { // Convolution of the implantation profile with Gaussian
-        fImpProfile->ConvolveGss(par[7], par[1]);
-      }
-
       fParForPofB[2] = par[1]; // energy
       fParForPofB[3] = par[2]; // Bkg-Field
       //fParForPofB[4] = 0.005; // Bkg-width (in principle zero)
@@ -823,7 +840,9 @@ TProximity1D1LHSGss::TProximity1D1LHSGss() : fCalcNeeded(true), fFirstCall(true)
     fParForPofB.push_back(startupHandler->GetDeltat());
     fParForPofB.push_back(startupHandler->GetDeltaB());
     fParForPofB.push_back(0.0);
-//    fParForPofB.push_back(0.0);
+    fParForPofB.push_back(0.0); // Bkg-Field
+    fParForPofB.push_back(0.01); // Bkg-width
+    fParForPofB.push_back(0.0); // Bkg-weight
 
     fImpProfile = new TTrimSPData(rge_path, energy_vec);
 
@@ -845,7 +864,7 @@ TProximity1D1LHSGss::TProximity1D1LHSGss() : fCalcNeeded(true), fFirstCall(true)
 //------------------
 // 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
+// Parameters: all the parameters for the function to be fitted through TProximity1D1LHSGss
 //------------------
 
 double TProximity1D1LHSGss::operator()(double t, const vector<double> &par) const {
@@ -857,19 +876,16 @@ double TProximity1D1LHSGss::operator()(double t, const vector<double> &par) cons
 
   // check if the function is called the first time and if yes, read in parameters
 
+  bool dead_layer_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;
-//  cout << this << endl;
+    fFirstCall = false;
+    dead_layer_changed = true;
   }
 
   // check if any parameter has changed
@@ -885,6 +901,9 @@ double TProximity1D1LHSGss::operator()(double t, const vector<double> &par) cons
         only_phase_changed = true;
       } else {
         only_phase_changed = false;
+        if (i == 4){
+          dead_layer_changed = true;
+        }
       }
     }
   }
@@ -906,12 +925,22 @@ double TProximity1D1LHSGss::operator()(double t, const vector<double> &par) cons
         fParForBofZ[i-2] = par[i];
 
       fParForPofB[2] = par[1]; // energy
+      fParForPofB[3] = par[2]; // Bkg-Field
+      //fParForPofB[4] = 0.005; // Bkg-width (in principle zero)
+
+      if(dead_layer_changed){
+        vector<double> interfaces;
+        interfaces.push_back(par[4]);// dead layer
+        fParForPofB[5] = fImpProfile->LayerFraction(par[1], 1, interfaces); // Fraction of muons in the deadlayer
+        interfaces.clear();
+      }
 
       TProximity1D_1LHSGss BofZ(fParForBofZ);
       fPofB->UnsetPBExists();
       fPofB->Calculate(&BofZ, fImpProfile, fParForPofB);
       fPofT->DoFFT();
 
+
     }/* else {
       cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl;
     }*/
@@ -930,7 +959,7 @@ double TProximity1D1LHSGss::operator()(double t, const vector<double> &par) cons
 // creates (a pointer to) the TPofTCalc object (with the FFT plan)
 //------------------
 
-TLondon1D3L::TLondon1D3L() : fCalcNeeded(true), fFirstCall(true), fLastThreeChanged(true) {
+TLondon1D3L::TLondon1D3L() : fCalcNeeded(true), fFirstCall(true) {
     // read startup file
     string startup_path_name("TFitPofB_startup.xml");
 
@@ -957,7 +986,10 @@ TLondon1D3L::TLondon1D3L() : fCalcNeeded(true), fFirstCall(true), fLastThreeChan
 
     fParForPofB.push_back(startupHandler->GetDeltat());
     fParForPofB.push_back(startupHandler->GetDeltaB());
-    fParForPofB.push_back(0.0);
+    fParForPofB.push_back(0.0); // Energy
+    fParForPofB.push_back(0.0); // Bkg-Field
+    fParForPofB.push_back(0.005); // Bkg-width
+    fParForPofB.push_back(0.0); // Bkg-weight
 
     fImpProfile = new TTrimSPData(rge_path, energy_vec);
 
@@ -984,25 +1016,26 @@ TLondon1D3L::TLondon1D3L() : fCalcNeeded(true), fFirstCall(true), fLastThreeChan
 
 double TLondon1D3L::operator()(double t, const vector<double> &par) const {
 
-  assert(par.size() == 13);
+  assert(par.size() == 10 || par.size() == 14);
 
   if(t<0.0)
     return cos(par[0]*0.017453293);
 
+  bool bkg_fraction_changed(false);
+  bool weights_changed(false);
+
   // 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;
+    bkg_fraction_changed = true;
+    if (par.size() == 14)
+      weights_changed = true;
   }
 
   // check if any parameter has changed
@@ -1018,9 +1051,11 @@ double TLondon1D3L::operator()(double t, const vector<double> &par) const {
         only_phase_changed = true;
       } else {
         only_phase_changed = false;
+        if (i == 3 || i == 4 || i == 5 || i == 6)
+          bkg_fraction_changed = true;
+        if (i == 10 || i == 11 || i == 12 || i == 13)
+          weights_changed = true;
       }
-      if (i == fPar.size()-3 || i == fPar.size()-2 || i == fPar.size()-1)
-        fLastThreeChanged = true;
     }
   }
 
@@ -1041,32 +1076,33 @@ double TLondon1D3L::operator()(double t, const vector<double> &par) const {
         fParForBofZ[i-2] = par[i];
 
       fParForPofB[2] = par[1]; // energy
+      fParForPofB[3] = par[2]; // Bkg-Field
+      //fParForPofB[4] = 0.005; // Bkg-width (in principle zero)
 
-/* 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(fLastThreeChanged) {
+      if(weights_changed) {
         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()-3); i<par.size(); i++)
+        for(unsigned int i(10); i<14; i++)
           weights.push_back(par[i]);
 
 //        cout << "Weighting has changed, re-calculating n(z) now..." << endl;
         fImpProfile->WeightLayers(par[1], interfaces, weights);
+
+        interfaces.clear();
+        weights.clear();
+      }
+
+      if(bkg_fraction_changed || weights_changed) {
+        vector<double> interfaces;
+        interfaces.push_back(par[3]);// dead layer
+        interfaces.push_back(par[3] + par[4] + par[5] + par[6]);// dead layer + first layer + second layer + third layer
+        fParForPofB[5] = fImpProfile->LayerFraction(par[1], 1, interfaces) + // Fraction of muons in the deadlayer
+                         fImpProfile->LayerFraction(par[1], 3, interfaces); // Fraction of muons in the substrate
+        interfaces.clear();
       }
 
       TLondon1D_3L BofZ3(fParForBofZ);
@@ -1074,7 +1110,6 @@ double TLondon1D3L::operator()(double t, const vector<double> &par) const {
       fPofB->Calculate(&BofZ3, fImpProfile, fParForPofB);
       fPofT->DoFFT();
 
-
     }/* else {
       cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl;
     }*/
@@ -1082,7 +1117,6 @@ double TLondon1D3L::operator()(double t, const vector<double> &par) const {
     fPofT->CalcPol(fParForPofT);
 
     fCalcNeeded = false;
-    fLastThreeChanged = false;
   }
 
   return fPofT->Eval(t);
@@ -1094,7 +1128,7 @@ double TLondon1D3L::operator()(double t, const vector<double> &par) const {
 // creates (a pointer to) the TPofTCalc object (with the FFT plan)
 //------------------
 
-TLondon1D3LS::TLondon1D3LS() : fCalcNeeded(true), fFirstCall(true), fLastThreeChanged(true) {
+TLondon1D3LS::TLondon1D3LS() : fCalcNeeded(true), fFirstCall(true) {
     // read startup file
     string startup_path_name("TFitPofB_startup.xml");
 
@@ -1121,7 +1155,10 @@ TLondon1D3LS::TLondon1D3LS() : fCalcNeeded(true), fFirstCall(true), fLastThreeCh
 
     fParForPofB.push_back(startupHandler->GetDeltat());
     fParForPofB.push_back(startupHandler->GetDeltaB());
-    fParForPofB.push_back(0.0);
+    fParForPofB.push_back(0.0); // Energy
+    fParForPofB.push_back(0.0); // Bkg-Field
+    fParForPofB.push_back(0.005); // Bkg-width
+    fParForPofB.push_back(0.0); // Bkg-weight
 
     fImpProfile = new TTrimSPData(rge_path, energy_vec);
 
@@ -1148,25 +1185,26 @@ TLondon1D3LS::TLondon1D3LS() : fCalcNeeded(true), fFirstCall(true), fLastThreeCh
 
 double TLondon1D3LS::operator()(double t, const vector<double> &par) const {
 
-  assert(par.size() == 12);
+  assert(par.size() == 9 || par.size() == 13);
 
   if(t<0.0)
     return cos(par[0]*0.017453293);
 
+  bool bkg_fraction_changed(false);
+  bool weights_changed(false);
+
   // 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;
+    bkg_fraction_changed = true;
+    if (par.size() == 13)
+      weights_changed = true;
   }
 
   // check if any parameter has changed
@@ -1182,9 +1220,11 @@ double TLondon1D3LS::operator()(double t, const vector<double> &par) const {
         only_phase_changed = true;
       } else {
         only_phase_changed = false;
+        if (i == 3 || i == 4 || i == 5 || i == 6)
+          bkg_fraction_changed = true;
+        if (i == 9 || i == 10 || i == 11 || i == 12)
+          weights_changed = true;
       }
-      if (i == fPar.size()-3 || i == fPar.size()-2 || i == fPar.size()-1)
-        fLastThreeChanged = true;
     }
   }
 
@@ -1205,23 +1245,38 @@ double TLondon1D3LS::operator()(double t, const vector<double> &par) const {
         fParForBofZ[i-2] = par[i];
 
       fParForPofB[2] = par[1]; // energy
+      fParForPofB[3] = par[2]; // Bkg-Field
+      //fParForPofB[4] = 0.005; // Bkg-width (in principle zero)
 
-      if(fLastThreeChanged) {
+      if(weights_changed) {
         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()-3); i<par.size(); i++)
+        for(unsigned int i(9); i<13; i++)
           weights.push_back(par[i]);
 
 //        cout << "Weighting has changed, re-calculating n(z) now..." << endl;
         fImpProfile->WeightLayers(par[1], interfaces, weights);
+
+        interfaces.clear();
+        weights.clear();
       }
 
-      TLondon1D_3LS BofZ3S(fParForBofZ);
+      if(bkg_fraction_changed || weights_changed) {
+        vector<double> interfaces;
+        interfaces.push_back(par[3]);// dead layer
+        interfaces.push_back(par[3] + par[4] + par[5] + par[6]);// dead layer + first layer + second layer + third layer
+        fParForPofB[5] = fImpProfile->LayerFraction(par[1], 1, interfaces) + // Fraction of muons in the deadlayer
+                         fImpProfile->LayerFraction(par[1], 3, interfaces); // Fraction of muons in the substrate
+        interfaces.clear();
+      }
+
+      TLondon1D_3LS BofZ3(fParForBofZ);
       fPofB->UnsetPBExists();
-      fPofB->Calculate(&BofZ3S, fImpProfile, fParForPofB);
+      fPofB->Calculate(&BofZ3, fImpProfile, fParForPofB);
       fPofT->DoFFT();
 
     }/* else {
@@ -1231,7 +1286,6 @@ double TLondon1D3LS::operator()(double t, const vector<double> &par) const {
     fPofT->CalcPol(fParForPofT);
 
     fCalcNeeded = false;
-    fLastThreeChanged = false;
   }
 
   return fPofT->Eval(t);
@@ -1396,183 +1450,6 @@ double TLondon1D3LS::operator()(double t, const vector<double> &par) const {
 // 
 // }
 
-//------------------
-// 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) {
-//    omp_set_nested(1);
-//    omp_set_dynamic(1);
-//    omp_set_num_threads(4);
-
-    // 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);
-    }
-
-    fNSteps = startupHandler->GetNSteps();
-    fWisdom = startupHandler->GetWisdomFile();
-    string rge_path(startupHandler->GetDataPath());
-    map<double, string> energy_vec(startupHandler->GetEnergies());
-
-    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);
-
-    fImpProfile = new TTrimSPData(rge_path, energy_vec);
-
-    fPofB = new TPofBCalc(fParForPofB);
-
-    fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT);
-
-    // clean up
-    if (saxParser) {
-      delete saxParser;
-      saxParser = 0;
-    }
-    if (startupHandler) {
-      delete startupHandler;
-      startupHandler = 0;
-    }
-}
-
-//------------------
-// TLondon1D3LSub-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 TLondon1D3LSub
-//------------------
-
-double TLondon1D3LSub::operator()(double t, const vector<double> &par) const {
-
-  assert(par.size() == 15);
-
-  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
-
-//#pragma omp critical
-//{
-  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()-4 || 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(fParForBofZ);
-      fPofB->UnsetPBExists();
-      fPofB->Calculate(&BofZ3, fImpProfile, fParForPofB);
-
-      // Add background contribution from the substrate
-      fPofB->AddBackground(par[2], par[14], fImpProfile->LayerFraction(par[1], 4, interfaces));
-
-      // FourierTransform of P(B)
-      fPofT->DoFFT();
-
-    }/* 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);
-
-}
-
 
 
 double TLondon1D3Lestimate::operator()(double z, const vector<double>& par) const {
diff --git a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp
index 9c3ccaed..d6e5fef6 100644
--- a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp
+++ b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp
@@ -205,7 +205,7 @@ void TPofBCalc::Calculate(const TBofZCalcInverse *BofZ, const TTrimSPData *dataT
   for (i = firstZerosEnd; i<=lastZerosStart; i++)
     fPB[i] /= pBsum;
 
-  if(para.size() == 6)
+  if(para.size() == 6 && para[5] != 0.0)
     AddBackground(para[3], para[4], para[5]);
 
 }
@@ -433,6 +433,7 @@ void TPofBCalc::Calculate(const TBulkVortexFieldCalc *vortexLattice, const vecto
   } 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]);
+//    ofstream of("LorentzWeight.dat");
     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))));
@@ -450,9 +451,14 @@ void TPofBCalc::Calculate(const TBulkVortexFieldCalc *vortexLattice, const vecto
                + (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);
+
+//          of << 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) << " ";
         }
       }
+//      of << endl;
     }
+//    of.close();
   } else {
     for (unsigned int j(0); j < numberOfSteps_2; ++j) {
       for (unsigned int i(0); i < numberOfSteps_2; ++i) {
diff --git a/src/external/TFitPofB-lib/include/TLondon1D.h b/src/external/TFitPofB-lib/include/TLondon1D.h
index 87032ce7..62c7de6a 100644
--- a/src/external/TFitPofB-lib/include/TLondon1D.h
+++ b/src/external/TFitPofB-lib/include/TLondon1D.h
@@ -107,7 +107,6 @@ private:
   mutable vector<double> fParForPofB;
   string fWisdom;
   unsigned int fNSteps;
-  mutable bool fLastTwoChanged;
 
   ClassDef(TLondon1D2L,1)
 };
@@ -184,7 +183,6 @@ private:
   mutable vector<double> fParForPofB;
   string fWisdom;
   unsigned int fNSteps;
-  mutable bool fLastThreeChanged;
 
   ClassDef(TLondon1D3L,1)
 };
@@ -211,7 +209,6 @@ private:
   mutable vector<double> fParForPofB;
   string fWisdom;
   unsigned int fNSteps;
-  mutable bool fLastThreeChanged;
 
   ClassDef(TLondon1D3LS,1)
 };
@@ -241,32 +238,6 @@ private:
 //   ClassDef(TLondon1D4L,1)
 // };
 
-class TLondon1D3LSub : public PUserFcnBase {
-
-public:
-  // default constructor
-  TLondon1D3LSub();
-  ~TLondon1D3LSub();
-
-  double operator()(double, const vector<double>&) const;
-
-private:
-  mutable vector<double> fPar;
-  TTrimSPData *fImpProfile;
-  TPofBCalc *fPofB;
-  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;
-  mutable bool fWeightsChanged;
-
-  ClassDef(TLondon1D3LSub,1)
-};
-
 // Class for fitting directly B(z) without any P(B)-calculation
 
 class TLondon1D3Lestimate : public PUserFcnBase {
diff --git a/src/external/TFitPofB-lib/include/TLondon1DLinkDef.h b/src/external/TFitPofB-lib/include/TLondon1DLinkDef.h
index bd816267..bcd64887 100644
--- a/src/external/TFitPofB-lib/include/TLondon1DLinkDef.h
+++ b/src/external/TFitPofB-lib/include/TLondon1DLinkDef.h
@@ -44,7 +44,6 @@
 #pragma link C++ class TLondon1D3L+;
 #pragma link C++ class TLondon1D3LS+;
 //#pragma link C++ class TLondon1D4L+;
-#pragma link C++ class TLondon1D3LSub+;
 
 #pragma link C++ class TLondon1D3Lestimate;