diff --git a/src/external/TFitPofB-lib/classes/TBofZCalc.cpp b/src/external/TFitPofB-lib/classes/TBofZCalc.cpp
index 99332251..8c9c16bb 100644
--- a/src/external/TFitPofB-lib/classes/TBofZCalc.cpp
+++ b/src/external/TFitPofB-lib/classes/TBofZCalc.cpp
@@ -52,8 +52,8 @@ double TBofZCalc::GetBofZ(double zz) const {
   }
 
   if (!found || i == 0) {
-    cout << "B(z) cannot be calculated for z = " << zz << " !" << endl;
-    cout << "Check your theory function!" << endl;
+    cout << "TBofZCalc::GetBofZ: B(z) cannot be calculated for z = " << zz << " !" << endl;
+    cout << "TBofZCalc::GetBofZ: Check your theory function!" << endl;
     return -1.0;
   }
 
diff --git a/src/external/TFitPofB-lib/classes/TFitPofBStartupHandler.cpp b/src/external/TFitPofB-lib/classes/TFitPofBStartupHandler.cpp
index 28bb1ae3..85ae238a 100644
--- a/src/external/TFitPofB-lib/classes/TFitPofBStartupHandler.cpp
+++ b/src/external/TFitPofB-lib/classes/TFitPofBStartupHandler.cpp
@@ -187,7 +187,7 @@ void TFitPofBStartupHandler::OnComment(const char *str)
  */
 void TFitPofBStartupHandler::OnWarning(const char *str)
 {
-  cout << endl << "TFitPofBStartupHandler **WARNING** " << str;
+  cout << endl << "TFitPofBStartupHandler::OnWarning: TFitPofBStartupHandler **WARNING** " << str;
   cout << endl;
 }
 
@@ -201,7 +201,7 @@ void TFitPofBStartupHandler::OnWarning(const char *str)
  */
 void TFitPofBStartupHandler::OnError(const char *str)
 {
-  cout << endl << "TFitPofBStartupHandler **ERROR** " << str;
+  cout << endl << "TFitPofBStartupHandler::OnError: TFitPofBStartupHandler **ERROR** " << str;
   cout << endl;
 }
 
@@ -215,7 +215,7 @@ void TFitPofBStartupHandler::OnError(const char *str)
  */
 void TFitPofBStartupHandler::OnFatalError(const char *str)
 {
-  cout << endl << "TFitPofBStartupHandler **FATAL ERROR** " << str;
+  cout << endl << "TFitPofBStartupHandler::OnFatalError: TFitPofBStartupHandler **FATAL ERROR** " << str;
   cout << endl;
 }
 
@@ -244,16 +244,16 @@ void TFitPofBStartupHandler::CheckLists()
   // check if anything was set, and if not set some default stuff
 
   // check if any data path is given
-  cout << endl << ">> check data path ...";
+  cout << endl << "TFitPofBStartupHandler::CheckLists: check data path ...";
   if (!fDataPath.size()) {
-    cout << endl << ">> This is not going to work, you have to set a valid data path where to find the rge-files in the xml-file!" << endl;
+    cout << endl << "TFitPofBStartupHandler::CheckLists: This is not going to work, you have to set a valid data path where to find the rge-files in the xml-file!" << endl;
     exit(-1);
   }
 
   // check if any energies are given
-  cout << endl << ">> check energy list ..." << endl;
+  cout << endl << "TFitPofBStartupHandler::CheckLists: check energy list ..." << endl;
   if (!fEnergyList.size()) {
-    cout << endl << ">> Energy list empty! Setting the default list." << endl;
+    cout << endl << "TFitPofBStartupHandler::CheckLists: Energy list empty! Setting the default list." << endl;
     char eChar[5];
     for(unsigned int i(0); i<33; i++) {
       for(unsigned int j(0); j<10; j++) {
@@ -264,30 +264,30 @@ void TFitPofBStartupHandler::CheckLists()
   }
 
   // check if delta_t is given, if not set default
-  cout << endl << ">> check specified time resolution ..." << endl;
+  cout << endl << "TFitPofBStartupHandler::CheckLists: check specified time resolution ..." << endl;
   if(!fDeltat) {
-    cout << endl << ">> You did not specify the time resolution. Setting the default." << endl;
+    cout << endl << "TFitPofBStartupHandler::CheckLists: You did not specify the time resolution. Setting the default." << endl;
     fDeltat = 0.01;
   }
 
   // check if delta_B is given, if not set default
-  cout << endl << ">> check specified field resolution ..." << endl;
+  cout << endl << "TFitPofBStartupHandler::CheckLists: check specified field resolution ..." << endl;
   if(!fDeltaB) {
-    cout << endl << ">> You did not specify the field resolution. Setting the default." << endl;
+    cout << endl << "TFitPofBStartupHandler::CheckLists: You did not specify the field resolution. Setting the default." << endl;
     fDeltaB = 0.05;
   }
 
   // check if any wisdom-file is specified
-  cout << endl << ">> check wisdom-file ..." << endl;
+  cout << endl << "TFitPofBStartupHandler::CheckLists: check wisdom-file ..." << endl;
   if (!fWisdomFile.size()) {
-    cout << endl << ">> You did not specify a wisdom file. Setting the default." << endl;
+    cout << endl << "TFitPofBStartupHandler::CheckLists: You did not specify a wisdom file. Setting the default." << endl;
     fWisdomFile = "WordsOfWisdom.dat";
   }
 
   // check if any number of steps for the theory function is specified
-  cout << endl << ">> check number of steps for theory ..." << endl;
+  cout << endl << "TFitPofBStartupHandler::CheckLists: check number of steps for theory ..." << endl;
   if (!fNSteps) {
-    cout << endl << ">> You did not specify the number of steps for the theory. Setting the default." << endl;
+    cout << endl << "TFitPofBStartupHandler::CheckLists: You did not specify the number of steps for the theory. Setting the default." << endl;
     fNSteps = 3000;
   }
 
diff --git a/src/external/TFitPofB-lib/classes/TLondon1D.cpp b/src/external/TFitPofB-lib/classes/TLondon1D.cpp
index 06971a5e..42e261d0 100644
--- a/src/external/TFitPofB-lib/classes/TLondon1D.cpp
+++ b/src/external/TFitPofB-lib/classes/TLondon1D.cpp
@@ -1113,7 +1113,7 @@ double TLondon1D3LSub::operator()(double t, const vector<double> &par) const {
       } else {
         only_phase_changed = false;
       }
-      if (i == fPar.size()-5 || i == fPar.size()-5 || i == fPar.size()-3 || i == fPar.size()-2)
+      if (i == fPar.size()-5 || i == fPar.size()-4 || i == fPar.size()-3 || i == fPar.size()-2)
         fWeightsChanged = true;
     }
   }
diff --git a/src/external/TFitPofB-lib/classes/TPofTCalc.cpp b/src/external/TFitPofB-lib/classes/TPofTCalc.cpp
index 273baa82..978c1cf8 100644
--- a/src/external/TFitPofB-lib/classes/TPofTCalc.cpp
+++ b/src/external/TFitPofB-lib/classes/TPofTCalc.cpp
@@ -68,7 +68,7 @@ TPofTCalc::TPofTCalc (const string &wisdom, const vector<double> &par) : fWisdom
   fFFTin = (double *)malloc(sizeof(double) * fNFFT);
   fFFTout = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (fNFFT/2+1));
 
-  cout << "Check for the FFT plan..." << endl;
+  cout << "TPofTCalc::TPofTCalc: Check for the FFT plan..." << endl;
 
   // Load wisdom from file
 
@@ -77,14 +77,14 @@ TPofTCalc::TPofTCalc (const string &wisdom, const vector<double> &par) : fWisdom
   FILE *wordsOfWisdomR;
   wordsOfWisdomR = fopen(fWisdom.c_str(), "r");
   if (wordsOfWisdomR == NULL) {
-    cout << "Couldn't open wisdom file ..." << endl;
+    cout << "TPofTCalc::TPofTCalc: Couldn't open wisdom file ..." << endl;
   } else {
     wisdomLoaded = fftw_import_wisdom_from_file(wordsOfWisdomR);
     fclose(wordsOfWisdomR);
   }
 
   if (!wisdomLoaded) {
-    cout << "No wisdom is imported..." << endl;
+    cout << "TPofTCalc::TPofTCalc: No wisdom is imported..." << endl;
   }
 
   fFFTplan = fftw_plan_dft_r2c_1d(fNFFT, fFFTin, fFFTout, FFTW_EXHAUSTIVE);
@@ -157,10 +157,10 @@ void TPofTCalc::CalcPol(const vector<double> &par) {
 
 //---------------------
 // Method for generating fake LEM decay histograms from p(B)
-// Parameters: par(dt, dB, timeres, channels, asyms, phases, t0s, N0s, bgs), output filename
+// Parameters: output filename, par(dt, dB, timeres, channels, asyms, phases, t0s, N0s, bgs)
 //---------------------
 
-void TPofTCalc::FakeData(const vector<double> &par, const string &rootOutputFileName) {
+void TPofTCalc::FakeData(const string &rootOutputFileName, const vector<double> &par) {
 
   //determine the number of histograms to be built
   unsigned int numHist(0);
@@ -168,16 +168,30 @@ void TPofTCalc::FakeData(const vector<double> &par, const string &rootOutputFile
     numHist=(par.size()-4)/5;
 
   if(!numHist){
-    cout << "The number of parameters for the histogram creation is not correct. Do nothing." << endl;
+    cout << "TPofTCalc::FakeData: The number of parameters for the histogram creation is not correct. Do nothing." << endl;
     return;
   }
 
-  cout << numHist << " histograms to be built" << endl;
+  cout << "TPofTCalc::FakeData: " << numHist << " histograms to be built" << endl;
 
-  vector<double> param;
-  param.push_back(0.0);
-  param.push_back(par[0]);
-  param.push_back(par[1]);
+  vector<unsigned int> t0;
+  vector<double> asy0;
+  vector<double> phase0;
+  vector<double> N0;
+  vector<double> bg;
+
+  for(unsigned int i(0); i<numHist; i++) {
+    t0.push_back(int(par[i+4+numHist*2]));
+    asy0.push_back(par[i+4]);
+    phase0.push_back(par[i+4+numHist]);
+    N0.push_back(par[i+4+numHist*3]);
+    bg.push_back(par[i+4+numHist*4]);
+  }
+
+  vector<double> param; // Parameters for TPofTCalc::CalcPol
+  param.push_back(0.0); // phase
+  param.push_back(par[0]); // dt
+  param.push_back(par[1]); // dB
 
   vector< vector<double> > asy;
   vector<double> asydata;
@@ -185,7 +199,7 @@ void TPofTCalc::FakeData(const vector<double> &par, const string &rootOutputFile
   double pol(0.0);
 
   for(unsigned int i(0); i<numHist; i++) {
-    param[0]=par[i+numHist*1+4];
+    param[0]=phase0[i];
     // calculate asymmetry
     CalcPol(param);
     for(unsigned int j(0); j<par[3]; j++){
@@ -193,14 +207,14 @@ void TPofTCalc::FakeData(const vector<double> &par, const string &rootOutputFile
       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(par[i+numHist*0+4]*pol);
+          asydata.push_back(asy0[i]*pol);
           break;
         }
       }
     }
     asy.push_back(asydata);
     asydata.clear();
-    cout << "Asymmetry " << i+1 << "/" << numHist << " calculated!" << endl;
+    cout << "TPofTCalc::FakeData: " << i+1 << "/" << numHist << " calculated!" << endl;
   }
 
   // calculate the histograms
@@ -210,20 +224,20 @@ void TPofTCalc::FakeData(const vector<double> &par, const string &rootOutputFile
 
   for (unsigned int i(0); i<numHist; i++) {    // loop over all histos
     for (unsigned int j(0); j<par[3]; j++) { // loop over time
-      if (j < par[i+numHist*2+4]) // j<t0
-        value = par[i+numHist*4+4]; // background
+      if (j < t0[i]) // j<t0
+        value = bg[i]; // background
       else
-        value = par[i+numHist*3+4]*exp(-par[2]*(j-par[i+numHist*2+4])/tauMu)*(1.0+asy[i][j- int(par[i+numHist*2+4])])+par[i+numHist*4+4];
+        value = N0[i]*exp(-par[2]*double(j-t0[i])/tauMu)*(1.0+asy[i][j-t0[i]])+bg[i];
       data.push_back(value);
     }
     histo.push_back(data);
     data.clear();
-    cout << endl << ">> histo " << i+1 << "/" << numHist << " done ...";
+    cout << "TPofTCalc::FakeData: " << i+1 << "/" << numHist << " done ...";
   }
 
   // add Poisson noise to the histograms
 
-  cout << endl << ">> add Poisson noise ..." << endl;
+  cout << endl << "TPofTCalc::FakeData: Adding Poisson noise ..." << endl;
 
   TH1F* theoHisto;
   TH1F* fakeHisto;
@@ -269,7 +283,7 @@ void TPofTCalc::FakeData(const vector<double> &par, const string &rootOutputFile
   runHeader->SetNHist(histoData.size());
   double *t0array = new double[histoData.size()];
   for (unsigned int i(0); i<histoData.size(); i++)
-    t0array[i] = par[i+numHist*2+4];
+    t0array[i] = t0[i];
   runHeader->SetTimeZero(t0array);
   if (t0array)
     delete t0array;
@@ -316,7 +330,13 @@ void TPofTCalc::FakeData(const vector<double> &par, const string &rootOutputFile
   histoData.clear();
   histoDataPPC.clear();
 
-  cout << endl << ">> DONE." << endl;
+  t0.clear();
+  asy0.clear();
+  phase0.clear();
+  N0.clear();
+  bg.clear();
+
+  cout << endl << "TPofTCalc::FakeData: DONE." << endl;
 
   return;
 }
@@ -331,7 +351,7 @@ double TPofTCalc::Eval(double t) const {
       return fPT[i]+(fPT[i+1]-fPT[i])/(fT[i+1]-fT[i])*(t-fT[i]);
   }
 
-  cout << "No data for the time " << t << " us available! Returning -999.0 ..." << endl;
+  cout << "TPofTCalc::Eval: No data for the time " << t << " us available! Returning -999.0 ..." << endl;
   return -999.0;
 }
 
@@ -345,7 +365,7 @@ TPofTCalc::~TPofTCalc() {
   FILE *wordsOfWisdomW;
   wordsOfWisdomW = fopen(fWisdom.c_str(), "w");
   if (wordsOfWisdomW == NULL) {
-    cout << "couldn't open file ... No wisdom is exported..." << endl;
+    cout << "TPofTCalc::~TPofTCalc(): Could not open file ... No wisdom is exported..." << endl;
   }
 
   fftw_export_wisdom_to_file(wordsOfWisdomW);
diff --git a/src/external/TFitPofB-lib/classes/TTrimSPDataHandler.cpp b/src/external/TFitPofB-lib/classes/TTrimSPDataHandler.cpp
index 86e4db89..10e308cc 100644
--- a/src/external/TFitPofB-lib/classes/TTrimSPDataHandler.cpp
+++ b/src/external/TFitPofB-lib/classes/TTrimSPDataHandler.cpp
@@ -5,7 +5,7 @@
   Author: Bastian M. Wojek
   e-mail: bastian.wojek@psi.ch
 
-  2008/09/02
+  2008/11/21
 
 ***************************************************************************/
 
@@ -14,6 +14,7 @@
 #include <fstream>
 #include <string>
 #include <cmath>
+#include <cassert>
 
 using namespace std;
 
@@ -25,6 +26,13 @@ using namespace std;
 //          vector<string> energyVec(energyArr, energyArr+(sizeof(energyArr)/sizeof(energyArr[0])));
 //
 // This will read the files "/home/user/TrimSP/SomeSample-02_1.rge", "/home/user/TrimSP/SomeSample-02_5.rge" and so on.
+//
+// Alternative supported energy formats in the energyVec are, e.g. for E=2.1keV:
+// <some_path>02-1.rge
+// <some_path>02.1.rge
+// <some_path>021.rge
+//
+// <some_path>21.rge is explicitly not possible since it is not clear, if this denotes 2.1keV or 21.0keV!
 //--------------------
 
 TTrimSPData::TTrimSPData(const string &path, vector<string> &energyVec) {
@@ -39,10 +47,24 @@ TTrimSPData::TTrimSPData(const string &path, vector<string> &energyVec) {
 
     ifstream *rgeFile = new ifstream(energyStr.c_str());
     if(! *rgeFile) {
-      cout << "rge-file not found! Try next energy..." << endl;
+      cout << "TTrimSPData::TTrimSPData: rge-file not found! Try next energy..." << endl;
       delete rgeFile;
+      rgeFile = 0;
     } else {
-      fEnergy.push_back(atof(energyVec[i].replace(2,1,".").c_str()));
+      if (energyVec[i].length() == 4)
+        fEnergy.push_back(atof(energyVec[i].replace(2,1,".").c_str()));
+      else if (energyVec[i].length() == 3) {
+        energyVec[i].insert(energyVec[i].end()-1, 1, '.');
+        fEnergy.push_back(atof(energyVec[i].c_str()));
+      } else {
+        cout << "TTrimSPData::TTrimSPData: The energy cannot be correctly extracted from the rge-file name!" << endl;
+        cout << "TTrimSPData::TTrimSPData: Please use file names in one of the following formats, e.g. for E=2.1keV use:" << endl;
+        cout << "TTrimSPData::TTrimSPData: <some_path>02_1.rge" << endl;
+        cout << "TTrimSPData::TTrimSPData: <some_path>02-1.rge" << endl;
+        cout << "TTrimSPData::TTrimSPData: <some_path>02.1.rge" << endl;
+        cout << "TTrimSPData::TTrimSPData: <some_path>021.rge" << endl;
+        assert(false);
+      }
 
       while(*rgeFile >> word)
         if(word == "PARTICLES") break;
@@ -53,19 +75,29 @@ TTrimSPData::TTrimSPData(const string &path, vector<string> &energyVec) {
         vnzz.push_back(nzz);
       }
 
+      fDZ.push_back(vzz[1]-vzz[0]);
+
+      while(zz < 2100.0){
+        zz += *(fDZ.end()-1);
+        vzz.push_back(zz);
+        vnzz.push_back(0.0);
+      }
+
       fDataZ.push_back(vzz);
       fDataNZ.push_back(vnzz);
 
+
       rgeFile->close();
       delete rgeFile;
       rgeFile = 0;
+
       vzz.clear();
       vnzz.clear();
 
     }
   }
 
-  cout << "Read in " << fDataNZ.size() << " implantation profiles in total." << endl;
+  cout << "TTrimSPData::TTrimSPData: Read in " << fDataNZ.size() << " implantation profiles in total." << endl;
 
   fOrigDataNZ = fDataNZ;
 
@@ -87,7 +119,7 @@ vector<double> TTrimSPData::DataZ(double e) const {
     }
   }
   // default
-  cout << "No implantation profile available for the specified energy... You get back the first one." << endl;
+  cout << "TTrimSPData::DataZ: No implantation profile available for the specified energy... You get back the first one." << endl;
   return fDataZ[0];
 
 }
@@ -105,7 +137,7 @@ vector<double> TTrimSPData::DataNZ(double e) const {
     }
   }
   // default
-  cout << "No implantation profile available for the specified energy... You get back the first one." << endl;
+  cout << "TTrimSPData::DataNZ: No implantation profile available for the specified energy... You get back the first one." << endl;
   return fDataNZ[0];
 
 }
@@ -122,7 +154,7 @@ vector<double> TTrimSPData::OrigDataNZ(double e) const {
     }
   }
   // default
-  cout << "No implantation profile available for the specified energy... You get back the first one." << endl;
+  cout << "TTrimSPData::OrigDataNZ: No implantation profile available for the specified energy... You get back the first one." << endl;
   return fOrigDataNZ[0];
 
 }
@@ -135,7 +167,7 @@ vector<double> TTrimSPData::OrigDataNZ(double e) const {
 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;
+    cout << "TTrimSPData::LayerFraction: No such layer available according to your specified interfaces... Returning 0.0!" << endl;
     return 0.0;
   }
 
@@ -162,13 +194,13 @@ double TTrimSPData::LayerFraction(double e, unsigned int layno, const vector<dou
             layerNumber += fDataNZ[i][j];
       }
       // fraction of muons in layer layno
-      cout << "Fraction of muons in layer " << layno << ": " << layerNumber/totalNumber << endl;
+      // 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;
+  cout << "TTrimSPData::LayerFraction: No implantation profile available for the specified energy... Returning 0.0" << endl;
   return 0.0;
 
 }
@@ -184,19 +216,19 @@ double TTrimSPData::LayerFraction(double e, unsigned int layno, const vector<dou
 void TTrimSPData::WeightLayers(double e, const vector<double>& interface, const vector<double>& weight) const {
 
   if(weight.size()-interface.size()-1) {
-    cout << "For the weighting the number of interfaces has to be one less than the number of weights!" << endl;
-    cout << "No weighting of the implantation profile will be done unless you take care of that!" << endl;
+    cout << "TTrimSPData::WeightLayers: For the weighting the number of interfaces has to be one less than the number of weights!" << endl;
+    cout << "TTrimSPData::WeightLayers: No weighting of the implantation profile will be done unless you take care of that!" << endl;
     return;
   }
 
   for(unsigned int i(0); i<interface.size(); i++) {
     if (interface[i]<0.0) {
-      cout << "One of your layer interfaces has a negative coordinate! - No weighting will be done!" << endl;
+      cout << "TTrimSPData::WeightLayers: One of your layer interfaces has a negative coordinate! - No weighting will be done!" << endl;
       return;
     }
     else if (i>1) {
       if (interface[i]<interface[i-1]) {
-        cout << "The specified interfaces appear to be not in ascending order! - No weighting will be done!" << endl;
+        cout << "TTrimSPData::WeightLayers: The specified interfaces appear to be not in ascending order! - No weighting will be done!" << endl;
         return;
       }
     }
@@ -204,7 +236,7 @@ void TTrimSPData::WeightLayers(double e, const vector<double>& interface, const
 
   for(unsigned int i(0); i<weight.size(); i++) {
     if (weight[i]>1.0 || weight[i]<0.0) {
-      cout << "At least one of the specified weights is out of range - no weighting will be done!" << endl;
+      cout << "TTrimSPData::WeightLayers: At least one of the specified weights is out of range - no weighting will be done!" << endl;
       return;
     }
   }
@@ -243,7 +275,7 @@ void TTrimSPData::WeightLayers(double e, const vector<double>& interface, const
     }
   }
 
-  cout << "No implantation profile available for the specified energy... No weighting done." << endl;
+  cout << "TTrimSPData::WeightLayers: No implantation profile available for the specified energy... No weighting done." << endl;
   return;
 }
 
@@ -262,7 +294,7 @@ double TTrimSPData::GetNofZ(double zz, double e) const {
       break;
     }
     if(i == fEnergy.size() - 1) {
-      cout << "No implantation profile available for the specified energy... Quitting!" << endl;
+      cout << "TTrimSPData::GetNofZ: No implantation profile available for the specified energy... Quitting!" << endl;
       exit(-1);
     }
   }
@@ -299,7 +331,7 @@ void TTrimSPData::Normalize(double e) const {
       double nZsum = 0.0;
       for (unsigned int j(0); j<fDataZ[i].size(); j++)
         nZsum += fDataNZ[i][j];
-      nZsum *= (fDataZ[i][1]-fDataZ[i][0]);
+      nZsum *= fDZ[i];
       for (unsigned int j(0); j<fDataZ[i].size(); j++)
         fDataNZ[i][j] /= nZsum;
 
@@ -308,7 +340,7 @@ void TTrimSPData::Normalize(double e) const {
     }
   }
   // default
-  cout << "No implantation profile available for the specified energy... No normalization done." << endl;
+  cout << "TTrimSPData::Normalize: No implantation profile available for the specified energy... No normalization done." << endl;
   return;
 
 }
@@ -324,10 +356,32 @@ bool TTrimSPData::IsNormalized(double e) const {
     }
   }
 
-  cout << "No implantation profile available for the specified energy... Returning false! Check your code!" << endl;
+  cout << "TTrimSPData::IsNormalized: No implantation profile available for the specified energy... Returning false! Check your code!" << endl;
   return false;
 }
 
+//---------------------
+// Calculate the mean range in (nm) for a given energy e
+//---------------------
+
+double TTrimSPData::MeanRange(double e) const {
+  for(unsigned int i(0); i<fEnergy.size(); i++) {
+    if(!(fEnergy[i] - e)) {
+      if (!fIsNormalized[i])
+        Normalize(e);
+      double mean(0.0);
+      for(unsigned int j(0); j<fDataNZ[i].size(); j++){
+        mean += fDataNZ[i][j]*fDataZ[i][j];
+      }
+      mean *= fDZ[i]/10.0;
+      return mean;
+    }
+  }
+
+  cout << "TTrimSPData::MeanRange: No implantation profile available for the specified energy... Returning -1! Check your code!" << endl;
+  return -1.;
+}
+
 //---------------------
 // Convolve the n(z)-vector calculated by trim.SP for a given energy e [keV] with a gaussian exp(-z^2/(2*w^2))
 // No normalization is done!
@@ -357,10 +411,12 @@ void TTrimSPData::ConvolveGss(double w, double e) const {
         fDataNZ[i][k] = nn;
       }
 
+      fIsNormalized[i] = false;
+
       return;
     }
   }
 
-  cout << "No implantation profile available for the specified energy... No convolution done!" << endl;
+  cout << "TTrimSPData::ConvolveGss: No implantation profile available for the specified energy... No convolution done!" << endl;
   return;
 }
diff --git a/src/external/TFitPofB-lib/include/TPofTCalc.h b/src/external/TFitPofB-lib/include/TPofTCalc.h
index ad1f0dac..83cf8bb8 100644
--- a/src/external/TFitPofB-lib/include/TPofTCalc.h
+++ b/src/external/TFitPofB-lib/include/TPofTCalc.h
@@ -30,7 +30,7 @@ public:
   vector<double> DataPT() const {return fPT;}
   void DoFFT(const TPofBCalc&);
   void CalcPol(const vector<double>&);
-  void FakeData(const vector<double>&, const string&); 
+  void FakeData(const string&, const vector<double>&); 
   double Eval(double) const;
 
 private:
diff --git a/src/external/TFitPofB-lib/include/TTrimSPDataHandler.h b/src/external/TFitPofB-lib/include/TTrimSPDataHandler.h
index cd46b510..14ac49b9 100644
--- a/src/external/TFitPofB-lib/include/TTrimSPDataHandler.h
+++ b/src/external/TFitPofB-lib/include/TTrimSPDataHandler.h
@@ -5,7 +5,7 @@
   Author: Bastian M. Wojek
   e-mail: bastian.wojek@psi.ch
 
-  2008/09/02
+  2008/11/21
 
 ***************************************************************************/
 
@@ -25,6 +25,10 @@ public:
   ~TTrimSPData() {
     fDataZ.clear();
     fDataNZ.clear();
+    fOrigDataNZ.clear();
+    fEnergy.clear();
+    fDZ.clear();
+    fIsNormalized.clear();
   }
 
   vector<double> Energy() const {return fEnergy;}
@@ -37,9 +41,11 @@ public:
   void Normalize(double) const;
   bool IsNormalized(double) const;
   void ConvolveGss(double, double) const;
+  double MeanRange(double) const;
 
 private:
   vector<double> fEnergy;
+  vector<double> fDZ;
   vector< vector<double> > fDataZ;
   mutable vector< vector<double> > fDataNZ;
   vector< vector<double> > fOrigDataNZ;