diff --git a/src/external/TFitPofB-lib/classes/TLondon1D.cpp b/src/external/TFitPofB-lib/classes/TLondon1D.cpp
index 075c690a..60c153f2 100644
--- a/src/external/TFitPofB-lib/classes/TLondon1D.cpp
+++ b/src/external/TFitPofB-lib/classes/TLondon1D.cpp
@@ -113,6 +113,7 @@ TLondon1DHS::TLondon1DHS() : fCalcNeeded(true), fFirstCall(true) {
     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;
@@ -196,6 +197,7 @@ double TLondon1DHS::operator()(double t, const vector<double> &par) const {
         fParForBofZ[i-2] = par[i];
 
       fParForPofB[2] = par[1]; // energy
+//      fParForPofB[3] = par[3]; // deadlayer for convolution of field profile
 
       TLondon1D_HS BofZ(fNSteps, fParForBofZ);
       TPofBCalc PofB(BofZ, *fImpProfile, fParForPofB);
@@ -220,7 +222,7 @@ double TLondon1DHS::operator()(double t, const vector<double> &par) const {
 // creates (a pointer to) the TPofTCalc object (with the FFT plan)
 //------------------
 
-TLondon1D1L::TLondon1D1L() : fCalcNeeded(true), fFirstCall(true) {
+TLondon1D1L::TLondon1D1L() : fCalcNeeded(true), fFirstCall(true), fCallCounter(0) {
 
     // read startup file
     string startup_path_name("TFitPofB_startup.xml");
@@ -275,6 +277,11 @@ TLondon1D1L::TLondon1D1L() : fCalcNeeded(true), fFirstCall(true) {
 //------------------
 
 double TLondon1D1L::operator()(double t, const vector<double> &par) const {
+
+  // Debugging
+  // Count the number of function calls
+  fCallCounter++;
+
   if(t<0.0)
     return 0.0;
 
@@ -283,16 +290,16 @@ double TLondon1D1L::operator()(double t, const vector<double> &par) const {
   if(fFirstCall){
     fPar = par;
 
-    for (unsigned int i(0); i<fPar.size(); i++){
+/*    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;
+//      cout << "fParForBofZ[" << i-2 << "] = " << fParForBofZ[i-2] << endl;
     }
     fFirstCall=false;
-  cout << this << endl;
+//  cout << this << endl;
   }
 
   // check if any parameter has changed
@@ -342,6 +349,16 @@ double TLondon1D1L::operator()(double t, const vector<double> &par) const {
 
     fCalcNeeded = false;
   }
+  
+  // Debugging start
+  if (!(fCallCounter%10000)){
+    cout << fCallCounter-1 << "\t";
+    for (unsigned int i(0); i<fPar.size(); i++){
+      cout << fPar[i] << "\t";
+    }
+    cout << endl;
+  }
+  // Debugging end
 
   return fPofT->Eval(t);
 
@@ -414,13 +431,13 @@ double TLondon1D2L::operator()(double t, const vector<double> &par) const {
   if(fFirstCall){
     fPar = par;
 
-    for (unsigned int i(0); i<fPar.size(); i++){
+/*    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;
+//      cout << "fParForBofZ[" << i-2 << "] = " << fParForBofZ[i-2] << endl;
     }
     fFirstCall=false;
   }
@@ -470,7 +487,7 @@ double TLondon1D2L::operator()(double t, const vector<double> &par) const {
         for(unsigned int i(par.size()-2); i<par.size(); i++)
           weights.push_back(par[i]);
 
-        cout << "Weighting has changed, re-calculating n(z) now..." << endl;
+//        cout << "Weighting has changed, re-calculating n(z) now..." << endl;
         fImpProfile->WeightLayers(par[1], interfaces, weights);
       }
 
@@ -559,13 +576,13 @@ double TLondon1D3L::operator()(double t, const vector<double> &par) const {
   if(fFirstCall){
     fPar = par;
 
-    for (unsigned int i(0); i<fPar.size(); i++){
+/*    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;
+//      cout << "fParForBofZ[" << i-2 << "] = " << fParForBofZ[i-2] << endl;
     }
     fFirstCall=false;
   }
@@ -630,7 +647,7 @@ double TLondon1D3L::operator()(double t, const vector<double> &par) const {
         for(unsigned int i(par.size()-3); i<par.size(); i++)
           weights.push_back(par[i]);
 
-        cout << "Weighting has changed, re-calculating n(z) now..." << endl;
+//        cout << "Weighting has changed, re-calculating n(z) now..." << endl;
         fImpProfile->WeightLayers(par[1], interfaces, weights);
       }
 
@@ -719,13 +736,13 @@ double TLondon1D3LS::operator()(double t, const vector<double> &par) const {
   if(fFirstCall){
     fPar = par;
 
-    for (unsigned int i(0); i<fPar.size(); i++){
+/*    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;
+//      cout << "fParForBofZ[" << i-2 << "] = " << fParForBofZ[i-2] << endl;
     }
     fFirstCall=false;
   }
@@ -776,7 +793,7 @@ double TLondon1D3LS::operator()(double t, const vector<double> &par) const {
         for(unsigned int i(par.size()-3); i<par.size(); i++)
           weights.push_back(par[i]);
 
-        cout << "Weighting has changed, re-calculating n(z) now..." << endl;
+//        cout << "Weighting has changed, re-calculating n(z) now..." << endl;
         fImpProfile->WeightLayers(par[1], interfaces, weights);
       }
 
diff --git a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp
index f672e91d..0f63063c 100644
--- a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp
+++ b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp
@@ -5,11 +5,11 @@
   Author: Bastian M. Wojek
   e-mail: bastian.wojek@psi.ch
 
-  2008/06/03
+  2008/09/04
 
 ***************************************************************************/
 
-#include "TPofBCalc.h"
+#include "TPofTCalc.h"
 #include <cmath>
 #include <iostream>
 #include <fstream>
@@ -17,14 +17,14 @@
 /* USED FOR DEBUGGING-----------------------------------
 #include <cstdio>
 #include <ctime>
--------------------------------------------------------*/
+/-------------------------------------------------------*/
 
 //-----------
 // Constructor that does the P(B) calculation for given B(z) and n(z)
 // Parameters: dt[us], dB[G], Energy[keV]
 //-----------
 
-TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, const vector<double> &para ) {
+TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, const vector<double> &para ) : fDT(para[0]), fDB(para[1]) {
 
   fBmin = BofZ.GetBmin();
   fBmax = BofZ.GetBmax();
@@ -34,7 +34,7 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons
 
   // fill not used Bs before Bmin with 0.0
 
-  for ( BB = 0.0 ; BB < fBmin ; BB += para[1] ) {
+  for ( BB = 0.0 ; BB < fBmin ; BB += fDB ) {
     fB.push_back(BB);
     fPB.push_back(0.0);
   }
@@ -69,23 +69,43 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons
   char debugfile1[50];
   int n1 = sprintf (debugfile1, "test_NZ_%ld_%f.dat", seconds, para[2]);
 
+  char debugfile2[50];
+  int n2 = sprintf (debugfile2, "test_NZgss_%ld_%f.dat", seconds, para[2]);
+
   if (n1 > 0) {
     ofstream of1(debugfile1);
 //    assure(of1, debugfile1);
 
+    dataTrimSP.Normalize(para[2]);
+
     for (unsigned int i(0); i<dataTrimSP.DataZ(para[2]).size(); i++) {
       of1 << dataTrimSP.DataZ(para[2])[i] << " " << dataTrimSP.DataNZ(para[2])[i] << " " << dataTrimSP.OrigDataNZ(para[2])[i] << endl;
     }
     of1.close();
   }
 
----------------------------------------------------------*/
+  if (n2 > 0) {
+    ofstream of2(debugfile2);
+//    assure(of1, debugfile1);
+
+    dataTrimSP.ConvolveGss(10.0,para[2]);
+    dataTrimSP.Normalize(para[2]);
+    
+    for (unsigned int i(0); i<dataTrimSP.DataZ(para[2]).size(); i++) {
+      of2 << dataTrimSP.DataZ(para[2])[i] << " " << dataTrimSP.DataNZ(para[2])[i] << " " << dataTrimSP.OrigDataNZ(para[2])[i] << endl;
+    }
+    of2.close();
+  }
+
+/---------------------------------------------------------*/
+
+//  dataTrimSP.ConvolveGss(10.0,para[2]); // convolve implantation profile by gaussian
 
   double nn;
   bool zNextFound(false);
 
-  for ( ; BB <= fBmax ; BB += para[1]) {
-    BBnext = BB + para[1];
+  for ( ; BB <= fBmax ; BB += fDB) {
+    BBnext = BB + fDB;
     fB.push_back(BB);
     fPB.push_back(0.0);
 
@@ -109,7 +129,7 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons
           nn = dataTrimSP.GetNofZ(zm, para[2]);
           if (nn != -1.0) {
 //          cout << "zNext = " << zNextm << ", zm = " << zm << ", dz = " << dz << endl;
-            *(fPB.end()-1) += nn*fabs(dz/para[1]);
+            *(fPB.end()-1) += nn*fabs(dz/fDB);
           }
         }
 
@@ -132,7 +152,7 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons
           nn = dataTrimSP.GetNofZ(zp, para[2]);
           if (nn != -1.0) {
 //            cout << "zNext = " << zNextp << ", zp = " << zp << ", dz = " << dz << endl;
-            *(fPB.end()-1) += nn*fabs(dz/para[1]);
+            *(fPB.end()-1) += nn*fabs(dz/fDB);
           }
         }
      }
@@ -144,11 +164,11 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons
 
   // fill not used Bs after Bext with 0.0
 
-  double BmaxFFT(1.0/gBar/para[0]);
+  double BmaxFFT(1.0/gBar/fDT);
 
 //  cout << "N = " << int(BmaxFFT/para[1]+1.0) << endl;
 
-  for ( ; BB <= BmaxFFT ; BB += para[1] ) {
+  for ( ; BB <= BmaxFFT ; BB += fDB ) {
     fB.push_back(BB);
     fPB.push_back(0.0);
   }
@@ -167,8 +187,82 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons
   double pBsum = 0.0;
   for (unsigned int i(firstZerosEnd); i<lastZerosStart; i++)
     pBsum += fPB[i];
-  pBsum *= para[1];
+  pBsum *= fDB;
   for (unsigned int i(firstZerosEnd); i<lastZerosStart; i++)
     fPB[i] /= pBsum;
 
 }
+
+
+void TPofBCalc::ConvolveGss(double w) {
+
+  if(!w)
+    return;
+
+  unsigned int NFFT;
+  double TBin;
+  fftw_plan FFTplanToTimeDomain;
+  fftw_plan FFTplanToFieldDomain;
+  fftw_complex *FFTout;
+  double *FFTin;
+
+  NFFT = ( int(1.0/gBar/fDT/fDB+1.0) % 2 ) ? int(1.0/gBar/fDT/fDB+2.0) : int(1.0/gBar/fDT/fDB+1.0);
+  TBin = 1.0/gBar/double(NFFT-1)/fDB;
+
+  FFTin = (double *)malloc(sizeof(double) * NFFT);
+  FFTout = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (NFFT/2+1));
+
+  // do the FFT to time domain
+
+  FFTplanToTimeDomain = fftw_plan_dft_r2c_1d(NFFT, FFTin, FFTout, FFTW_ESTIMATE);
+
+  for (unsigned int i(0); i<NFFT; i++) {
+    FFTin[i] = fPB[i];
+  }
+  for (unsigned int i(0); i<NFFT/2+1; i++) {
+    FFTout[i][0] = 0.0;
+    FFTout[i][1] = 0.0;
+  }
+
+  fftw_execute(FFTplanToTimeDomain);
+
+  // multiply everything by a gaussian
+
+  double GssInTimeDomain;
+  double expo(-2.0*PI*PI*gBar*gBar*w*w*TBin*TBin);
+
+  for (unsigned int i(0); i<NFFT/2+1; i++) {
+    GssInTimeDomain = exp(expo*double(i)*double(i));
+    FFTout[i][0] *= GssInTimeDomain;
+    FFTout[i][1] *= GssInTimeDomain;
+  }
+
+  // FFT back to the field domain
+
+  FFTplanToFieldDomain = fftw_plan_dft_c2r_1d(NFFT, FFTout, FFTin, FFTW_ESTIMATE);
+
+  fftw_execute(FFTplanToFieldDomain);
+
+  for (unsigned int i(0); i<NFFT; i++) {
+    fPB[i] = FFTin[i];
+  }
+
+  // cleanup
+
+  fftw_destroy_plan(FFTplanToTimeDomain);
+  fftw_destroy_plan(FFTplanToFieldDomain);
+  free(FFTin);
+  fftw_free(FFTout);
+//  fftw_cleanup();
+
+  // normalize p(B)
+
+  double pBsum = 0.0;
+  for (unsigned int i(0); i<NFFT; i++)
+    pBsum += fPB[i];
+  pBsum *= fDB;
+  for (unsigned int i(0); i<NFFT; i++)
+    fPB[i] /= pBsum;
+
+  return;
+}
diff --git a/src/external/TFitPofB-lib/classes/TPofTCalc.cpp b/src/external/TFitPofB-lib/classes/TPofTCalc.cpp
index 9e0b312e..1ad25e0b 100644
--- a/src/external/TFitPofB-lib/classes/TPofTCalc.cpp
+++ b/src/external/TFitPofB-lib/classes/TPofTCalc.cpp
@@ -5,7 +5,7 @@
   Author: Bastian M. Wojek
   e-mail: bastian.wojek@psi.ch
 
-  2008/08/28
+  2008/09/02
 
 ***************************************************************************/
 
@@ -18,7 +18,7 @@
 /* USED FOR DEBUGGING -----------------------
 #include <ctime>
 #include <fstream>
---------------------------------------------*/
+/--------------------------------------------*/
 
 //------------------
 // Constructor of the TPofTCalc class - it creates the FFT plan
@@ -84,7 +84,7 @@ void TPofTCalc::DoFFT(const TPofBCalc &PofB) {
     }
     of.close();
   }
---------------------------------------------*/
+/--------------------------------------------*/
 
   for (unsigned int i(0); i<fNFFT; i++) {
     fFFTin[i] = pB[i];
diff --git a/src/external/TFitPofB-lib/classes/TTrimSPDataHandler.cpp b/src/external/TFitPofB-lib/classes/TTrimSPDataHandler.cpp
index d10d8bc2..42122ca8 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/05/30
+  2008/09/02
 
 ***************************************************************************/
 
@@ -221,6 +221,9 @@ double TTrimSPData::GetNofZ(double zz, double e) const {
     }
   }
 
+  if(zz < 0)
+    return 0.0;
+
   bool found = false;
   unsigned int i;
   for (i=0; i<z.size(); i++) {
@@ -231,7 +234,7 @@ double TTrimSPData::GetNofZ(double zz, double e) const {
   }
 
   if (!found)
-    return -1.0;
+    return 0.0;
 
   if (i == 0)
     return nz[0]*10.0*zz/z[0];
@@ -243,7 +246,7 @@ double TTrimSPData::GetNofZ(double zz, double e) const {
 // Method normalizing the n(z)-vector calculated by trim.SP for a given energy[keV]
 //---------------------
 
-void TTrimSPData::Normalize(double e) {
+void TTrimSPData::Normalize(double e) const {
 
   for(unsigned int i(0); i<fEnergy.size(); i++) {
     if(!(fEnergy[i] - e)) {
@@ -278,3 +281,40 @@ bool TTrimSPData::IsNormalized(double e) const {
   cout << "No implantation profile available for the specified energy... Returning false! Check your code!" << endl;
   return false;
 }
+
+//---------------------
+// 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!
+//---------------------
+
+void TTrimSPData::ConvolveGss(double w, double e) const {
+  if(!w)
+    return;
+
+  vector<double> z, nz, gss;
+  double nn;
+
+  for(unsigned int i(0); i<fEnergy.size(); i++) {
+    if(!(fEnergy[i] - e)) {
+      z = fDataZ[i];
+      nz = fOrigDataNZ[i];
+
+      for(unsigned int k(0); k<z.size(); k++) {
+        gss.push_back(exp(-z[k]*z[k]/200.0/w/w));
+      }
+
+      for(unsigned int k(0); k<nz.size(); k++) {
+        nn = 0.0;
+        for(unsigned int j(0); j<nz.size(); j++) {
+          nn += nz[j]*gss[abs(int(k)-int(j))];
+        }
+        fDataNZ[i][k] = nn;
+      }
+
+      return;
+    }
+  }
+
+  cout << "No implantation profile available for the specified energy... No convolution done!" << endl;
+  return;
+}
diff --git a/src/external/TFitPofB-lib/include/TLondon1D.h b/src/external/TFitPofB-lib/include/TLondon1D.h
index 0c9b337e..8bf3cdd7 100644
--- a/src/external/TFitPofB-lib/include/TLondon1D.h
+++ b/src/external/TFitPofB-lib/include/TLondon1D.h
@@ -18,7 +18,7 @@
 class TLondon1DHS : public PUserFcnBase {
 
 public:
-  // default conctructor
+  // default constructor
   TLondon1DHS();
   ~TLondon1DHS();
 
@@ -42,7 +42,7 @@ private:
 class TLondon1D1L : public PUserFcnBase {
 
 public:
-  // default conctructor
+  // default constructor
   TLondon1D1L();
   ~TLondon1D1L();
 
@@ -59,6 +59,7 @@ private:
   mutable vector<double> fParForPofB;
   string fWisdom;
   unsigned int fNSteps;
+  mutable unsigned int fCallCounter;
 
   ClassDef(TLondon1D1L,1)
 };
@@ -66,7 +67,7 @@ private:
 class TLondon1D2L : public PUserFcnBase {
 
 public:
-  // default conctructor
+  // default constructor
   TLondon1D2L();
   ~TLondon1D2L();
 
@@ -91,7 +92,7 @@ private:
 class TLondon1D3L : public PUserFcnBase {
 
 public:
-  // default conctructor
+  // default constructor
   TLondon1D3L();
   ~TLondon1D3L();
 
@@ -117,7 +118,7 @@ private:
 class TLondon1D3LS : public PUserFcnBase {
 
 public:
-  // default conctructor
+  // default constructor
   TLondon1D3LS();
   ~TLondon1D3LS();
 
diff --git a/src/external/TFitPofB-lib/include/TPofBCalc.h b/src/external/TFitPofB-lib/include/TPofBCalc.h
index e59e41c8..9d256605 100644
--- a/src/external/TFitPofB-lib/include/TPofBCalc.h
+++ b/src/external/TFitPofB-lib/include/TPofBCalc.h
@@ -5,7 +5,7 @@
   Author: Bastian M. Wojek
   e-mail: bastian.wojek@psi.ch
 
-  2008/05/25
+  2008/09/04
 
 ***************************************************************************/
 
@@ -31,13 +31,15 @@ public:
   vector<double> DataPB() const {return fPB;}
   double GetBmin() const {return fBmin;}
   double GetBmax() const {return fBmax;}
+  void ConvolveGss(double);
 
 private:
   vector<double> fB;
   vector<double> fPB;
   double fBmin;
   double fBmax;
-
+  double fDT;
+  double fDB;
 };
 
 #endif // _TPofBCalc_H_
diff --git a/src/external/TFitPofB-lib/include/TTrimSPDataHandler.h b/src/external/TFitPofB-lib/include/TTrimSPDataHandler.h
index 8289596d..c87e7f65 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/05/26
+  2008/09/02
 
 ***************************************************************************/
 
@@ -33,8 +33,9 @@ public:
   vector<double> OrigDataNZ(double) const;
   void WeightLayers(double, const vector<double>&, const vector<double>&) const;
   double GetNofZ(double, double) const;
-  void Normalize(double);
+  void Normalize(double) const;
   bool IsNormalized(double) const;
+  void ConvolveGss(double, double) const;
 
 private:
   vector<double> fEnergy;
diff --git a/src/external/TFitPofB-lib/test/test.cpp b/src/external/TFitPofB-lib/test/test.cpp
index 20b17fb9..5b1ff431 100644
--- a/src/external/TFitPofB-lib/test/test.cpp
+++ b/src/external/TFitPofB-lib/test/test.cpp
@@ -5,41 +5,105 @@
 using namespace std;
 
 int main(){
-/*  string rge_path("/home/l_wojek/nt/wojek/g/Bastian/ImplantationDepth/YBCO_PBCO-");
-  string energy_arr[] = {"02_1", "02_5", "03_5", "05_0", "07_5", "10_0", "12_5", "15_0", "17_5", "19_0", "20_0", "22_5", "25_0"};
-  
+  string rge_path("/home/l_wojek/TrimSP/YBCOxtal/YBCOxtal-500000-");
+  string energy_arr[] = {"03_0", "03_6", "05_0", "05_3", "07_0", "07_7", "08_0", "09_0", "10_0", "10_2", "12_0", "14_1", "16_0", "16_4", "18_0", "19_7", "20_0", "22_0", "24_0", "24_6"};
+    
   vector<string> energy_vec(energy_arr, energy_arr+(sizeof(energy_arr)/sizeof(energy_arr[0])));
 
   TTrimSPData calcData(rge_path, energy_vec);
   
-  vector<double> energies(calcData.Energy());
+/*  vector<double> energies(calcData.Energy());
   for (unsigned int i(0); i<energies.size(); i++)
     cout << energies[i] << endl;
-  
-  vector<double> z(calcData.DataZ(2.5));
-  vector<double> nz(calcData.DataNZ(2.5));
-  
-  vector<double> z2(calcData.DataZ(25.0));
-  vector<double> nz2(calcData.DataNZ(25.0));
-  
-  ofstream of("test_out1.dat");
+*/
+  calcData.Normalize(22.0);
+
+  vector<double> z(calcData.DataZ(22.0));
+  vector<double> nz(calcData.DataNZ(22.0));
+
+  ofstream of("Implantation-profile-normal.dat");
   for (unsigned int i(0); i<z.size(); i++) {
     of << z[i] << " " << nz[i] << endl;
   }
   of.close();
+   
+  double par_arr[] = {24.4974, 22.0, 95.8253, 7.62096, 143.215};
+  vector<double> par_vec(par_arr, par_arr+(sizeof(par_arr)/sizeof(par_arr[0])));
   
-  ofstream of2("test_out2.dat");
+  vector<double> parForBofZ;
+  for (unsigned int i(2); i<par_vec.size(); i++)
+    parForBofZ.push_back(par_vec[i]);
+
+  vector<double> parForPofB;
+  parForPofB.push_back(0.01); //dt
+  parForPofB.push_back(0.2); //dB
+  parForPofB.push_back(par_vec[1]);
+
+  TLondon1D_HS BofZ(3000, parForBofZ);
+
+  TPofBCalc PofB(BofZ, calcData, parForPofB);
+
+  vector<double> hurgaB(PofB.DataB());
+  vector<double> hurgaPB(PofB.DataPB());
+  
+  ofstream of7("BpB-normal.dat");
+  for (unsigned int i(0); i<hurgaB.size(); i++) {
+    of7 << hurgaB[i] << " " << hurgaPB[i] << endl;
+  }
+  of7.close();
+  
+  PofB.ConvolveGss(5.0);
+  
+  vector<double> hurgaB1(PofB.DataB());
+  vector<double> hurgaPB1(PofB.DataPB());
+  
+  ofstream of1("BpB-field-broad.dat");
+  for (unsigned int i(0); i<hurgaB1.size(); i++) {
+    of1 << hurgaB1[i] << " " << hurgaPB1[i] << endl;
+  }
+  of1.close();
+
+  calcData.ConvolveGss(15.0, par_vec[1]);
+  calcData.Normalize(par_vec[1]);
+  
+  TPofBCalc PofB2(BofZ, calcData, parForPofB);
+  
+  vector<double> z2(calcData.DataZ(par_vec[1]));
+  vector<double> nz2(calcData.DataNZ(par_vec[1]));
+  
+  ofstream of2("Implantation-profile-broad.dat");
   for (unsigned int i(0); i<z2.size(); i++) {
     of2 << z2[i] << " " << nz2[i] << endl;
   }
   of2.close();
   
-  ofstream of3("test_out4.dat");
-  for (unsigned int i(0); i<z.size(); i++) {
-    of3 << i << " " << calcData.GetNofZ(double(i), 25.0) << endl;
-  }
-  of3.close();
+  vector<double> hurgaB2(PofB2.DataB());
+  vector<double> hurgaPB2(PofB2.DataPB());
   
+  ofstream of8("BpB-profile-broad.dat");
+  for (unsigned int i(0); i<hurgaB2.size(); i++) {
+    of8 << hurgaB2[i] << " " << hurgaPB2[i] << endl;
+  }
+  of8.close();
+  
+  PofB2.ConvolveGss(5.0);
+  
+  vector<double> hurgaB3(PofB2.DataB());
+  vector<double> hurgaPB3(PofB2.DataPB());
+  
+  ofstream of9("BpB-profile+field-broad.dat");
+  for (unsigned int i(0); i<hurgaB3.size(); i++) {
+    of9 << hurgaB3[i] << " " << hurgaPB3[i] << endl;
+  }
+  of9.close();
+  
+  z.clear();
+  nz.clear();
+  z2.clear();
+  nz2.clear();
+  
+/*
+
   double param[8] = {100.0, 5.0, 50.0, 100.0, 40.0, 0.01, 0.1, 15.0};
   vector<double> parameter(param,param+8);
   vector<double> param_for_BofZ;
@@ -97,10 +161,10 @@ int main(){
   of8.close();
 
 */
-/**************** Test TLondon1DHS *********************************/
+/**************** Test TLondon1DHS *********************************
 
 //  unsigned int parNo_arr[] = {1, 2, 5, 7, 9, 10, 11, 12};
-  double par_arr[] = {20.0, 24.6, 100.0, 15.0, 140.0};
+  double par_arr[] = {24.4974, 22.0, 95.8253, 7.62096, 143.215};
   
 //  vector<unsigned int> parNo_vec(parNo_arr, parNo_arr+(sizeof(parNo_arr)/sizeof(parNo_arr[0])));
   vector<double> par_vec(par_arr, par_arr+(sizeof(par_arr)/sizeof(par_arr[0])));
@@ -113,7 +177,7 @@ int main(){
   
   TLondon1DHS fitter;
 
-/************************************************************************/
+************************************************************************/
 
 
 /**************** Test TLondon1D1L *********************************
@@ -189,83 +253,83 @@ int main(){
 
 ************************************************************************/
 
-  ofstream of01("test_fitter01.dat");
-  ofstream of02("test_fitter02.dat");
-  ofstream of03("test_fitter03.dat");
-  ofstream of04("test_fitter04.dat");
-  ofstream of05("test_fitter05.dat");
-  ofstream of06("test_fitter06.dat");
-  ofstream of07("test_fitter07.dat");
-  ofstream of08("test_fitter08.dat");
-//  ofstream of09("test_fitter09.dat");
-//  ofstream of10("test_fitter10.dat");
+//  ofstream of01("test_fitter01.dat");
+//   ofstream of02("test_fitter02.dat");
+//   ofstream of03("test_fitter03.dat");
+//   ofstream of04("test_fitter04.dat");
+//   ofstream of05("test_fitter05.dat");
+//   ofstream of06("test_fitter06.dat");
+//   ofstream of07("test_fitter07.dat");
+//   ofstream of08("test_fitter08.dat");
+//   ofstream of09("test_fitter09.dat");
+//   ofstream of10("test_fitter10.dat");
   
 
-  for (double i(0.); i<12.0; i+=0.003) {
-    of01 << i << " " << fitter(i, par_vec) << endl;
-  }
-  of01.close();
+//   for (double i(0.); i<12.0; i+=0.003) {
+//     of01 << i << " " << fitter(i, par_vec) << endl;
+//   }
+//   of01.close();
   
-  par_vec[1] = 7.7;
-  
-  for (double i(0.); i<12.0; i+=0.003) {
-    of02 << i << " " << fitter(i, par_vec) << endl;
-  }
-  of02.close();
-
-  par_vec[0] = 0.0;
-  
-  for (double i(0.); i<12.0; i+=0.003) {
-    of03 << i << " " << fitter(i, par_vec) << endl;
-  }
-  of03.close();
-  
-    par_vec[2] = 200.0;
-  
-  for (double i(0.); i<12.0; i+=0.003) {
-    of04 << i << " " << fitter(i, par_vec) << endl;
-  }
-  of04.close();
-
-    par_vec[4] = 100.0;
-  
-  for (double i(0.); i<12.0; i+=0.003) {
-    of05 << i << " " << fitter(i, par_vec) << endl;
-  }
-  of05.close();
-
-  par_vec[0] = 20.0;
-  par_vec[1] = 24.6;
-  par_vec[2] = 96.5;
-  par_vec[3] = 15.0;
-  par_vec[4] = 130.0;
-  
-  for (double i(0.); i<12.0; i+=0.003) {
-    of06 << i << " " << fitter(i, par_vec) << endl;
-  }
-  of06.close();
-  
-  par_vec[0] = 20.0;
-  par_vec[1] = 24.6;
-  par_vec[2] = 96.5;
-  par_vec[3] = 15.0;
-  par_vec[4] = 140.0;
-  
-  for (double i(0.); i<12.0; i+=0.003) {
-    of07 << i << " " << fitter(i, par_vec) << endl;
-  }
-  of07.close();
-  
-  par_vec[0] = 20.0;
-  par_vec[1] = 24.6;
-  par_vec[2] = 96.5;
-  par_vec[3] = 20.0;
-  par_vec[4] = 130.0;
-  
-  for (double i(0.); i<12.0; i+=0.003) {
-    of08 << i << " " << fitter(i, par_vec) << endl;
-  }
-  of08.close();
+//   par_vec[1] = 7.7;
+//   
+//   for (double i(0.); i<12.0; i+=0.003) {
+//     of02 << i << " " << fitter(i, par_vec) << endl;
+//   }
+//   of02.close();
+// 
+//   par_vec[0] = 0.0;
+//   
+//   for (double i(0.); i<12.0; i+=0.003) {
+//     of03 << i << " " << fitter(i, par_vec) << endl;
+//   }
+//   of03.close();
+//   
+//     par_vec[2] = 200.0;
+//   
+//   for (double i(0.); i<12.0; i+=0.003) {
+//     of04 << i << " " << fitter(i, par_vec) << endl;
+//   }
+//   of04.close();
+// 
+//     par_vec[4] = 100.0;
+//   
+//   for (double i(0.); i<12.0; i+=0.003) {
+//     of05 << i << " " << fitter(i, par_vec) << endl;
+//   }
+//   of05.close();
+// 
+//   par_vec[0] = 20.0;
+//   par_vec[1] = 24.6;
+//   par_vec[2] = 96.5;
+//   par_vec[3] = 15.0;
+//   par_vec[4] = 130.0;
+//   
+//   for (double i(0.); i<12.0; i+=0.003) {
+//     of06 << i << " " << fitter(i, par_vec) << endl;
+//   }
+//   of06.close();
+//   
+//   par_vec[0] = 20.0;
+//   par_vec[1] = 24.6;
+//   par_vec[2] = 96.5;
+//   par_vec[3] = 15.0;
+//   par_vec[4] = 140.0;
+//   
+//   for (double i(0.); i<12.0; i+=0.003) {
+//     of07 << i << " " << fitter(i, par_vec) << endl;
+//   }
+//   of07.close();
+//   
+//   par_vec[0] = 20.0;
+//   par_vec[1] = 24.6;
+//   par_vec[2] = 96.5;
+//   par_vec[3] = 20.0;
+//   par_vec[4] = 130.0;
+//   
+//   for (double i(0.); i<12.0; i+=0.003) {
+//     of08 << i << " " << fitter(i, par_vec) << endl;
+//   }
+//   of08.close();
 /*
     par_vec_sub[0] = 0.0;
     par_vec_sub[7] = 1000.0;
@@ -347,7 +411,7 @@ int main(){
   */
   
 
-  par_vec.clear();
+//  par_vec.clear();