From c16916b5634a368bc1940e01eb58c331340ec0aa Mon Sep 17 00:00:00 2001
From: "Bastian M. Wojek" <basti@kth.se>
Date: Sun, 25 May 2008 14:40:13 +0000
Subject: [PATCH] Added a TTrimSPDataHandler method and changed TFitPofB to
 TUserLondon since it has become quite specific.

---
 .../TFitPofB-lib/classes/Makefile.TFitPofB    |  2 +-
 .../TFitPofB-lib/classes/TPofBCalc.cpp        | 26 +++++-
 .../TFitPofB-lib/classes/TPofTCalc.cpp        | 16 +++-
 .../classes/TTrimSPDataHandler.cpp            | 79 ++++++++++++++++++-
 .../classes/{TFitPofB.cpp => TUserLondon.cpp} | 54 ++++++++++---
 src/external/TFitPofB-lib/include/TPofBCalc.h |  5 +-
 src/external/TFitPofB-lib/include/TPofTCalc.h |  6 +-
 .../TFitPofB-lib/include/TTrimSPDataHandler.h |  7 +-
 .../include/{TFitPofB.h => TUserLondon.h}     | 20 ++---
 src/external/TFitPofB-lib/test/test.cpp       | 13 +--
 10 files changed, 185 insertions(+), 43 deletions(-)
 rename src/external/TFitPofB-lib/classes/{TFitPofB.cpp => TUserLondon.cpp} (70%)
 rename src/external/TFitPofB-lib/include/{TFitPofB.h => TUserLondon.h} (60%)

diff --git a/src/external/TFitPofB-lib/classes/Makefile.TFitPofB b/src/external/TFitPofB-lib/classes/Makefile.TFitPofB
index d1865a52..b927d1c2 100644
--- a/src/external/TFitPofB-lib/classes/Makefile.TFitPofB
+++ b/src/external/TFitPofB-lib/classes/Makefile.TFitPofB
@@ -12,7 +12,7 @@ OBJS += TTrimSPDataHandler.o
 OBJS += TBofZCalc.o
 OBJS += TPofBCalc.o
 OBJS += TPofTCalc.o
-OBJS += TFitPofB.o
+OBJS += TUserLondon.o
 
 SHLIB = libTFitPofB.so
 
diff --git a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp
index 4999ce7b..de2aac13 100644
--- a/src/external/TFitPofB-lib/classes/TPofBCalc.cpp
+++ b/src/external/TFitPofB-lib/classes/TPofBCalc.cpp
@@ -5,7 +5,7 @@
   Author: Bastian M. Wojek
   e-mail: bastian.wojek@psi.ch
 
-  2008/05/24
+  2008/05/25
 
 ***************************************************************************/
 
@@ -13,7 +13,11 @@
 #include <cmath>
 #include <iostream>
 #include <fstream>
-//#include <cstdio>
+
+/* USED FOR DEBUGGING-----------------------------------
+#include <cstdio>
+#include <ctime>
+-------------------------------------------------------*/
 
 //-----------
 // Constructor that does the P(B) calculation for given B(z) and n(z)
@@ -46,8 +50,11 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons
 /* USED FOR DEBUGGING-----------------------------------
   cout << "Bmin = " << fBmin << ", Bmax = " << fBmax << endl;
 
+  time_t seconds;
+  seconds = time (NULL);
+
   char debugfile[50];
-  int n = sprintf (debugfile, "test_Bz_%f.dat", fBmin);
+  int n = sprintf (debugfile, "test_Bz_%f_%ld.dat", fBmin, seconds);
 
   if (n > 0) {
     ofstream of(debugfile);
@@ -57,6 +64,19 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons
     }
     of.close();
   }
+
+  char debugfile1[50];
+  int n1 = sprintf (debugfile1, "test_NZ_%f_%ld.dat", para[2], seconds);
+
+  if (n1 > 0) {
+    ofstream of1(debugfile1);
+
+    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();
+  }
+
 ---------------------------------------------------------*/
 
   double nn;
diff --git a/src/external/TFitPofB-lib/classes/TPofTCalc.cpp b/src/external/TFitPofB-lib/classes/TPofTCalc.cpp
index 6077710c..d2613313 100644
--- a/src/external/TFitPofB-lib/classes/TPofTCalc.cpp
+++ b/src/external/TFitPofB-lib/classes/TPofTCalc.cpp
@@ -5,16 +5,20 @@
   Author: Bastian M. Wojek
   e-mail: bastian.wojek@psi.ch
 
-  2008/05/24
+  2008/05/25
 
 ***************************************************************************/
 
 #include "TPofTCalc.h"
 #include "fftw3.h"
 #include <cmath>
-#include <cstdio>
-#include <fstream>
 #include <iostream>
+#include <cstdio>
+
+/* USED FOR DEBUGGING -----------------------
+#include <ctime>
+#include <fstream>
+--------------------------------------------*/
 
 //------------------
 // Constructor of the TPofTCalc class - it creates the FFT plan
@@ -61,11 +65,15 @@ void TPofTCalc::DoFFT(const TPofBCalc &PofB) {
   vector<double> pB(PofB.DataPB());
 
 /* USED FOR DEBUGGING -----------------------
+
+  time_t seconds;
+  seconds = time(NULL);
+
   vector<double> B(PofB.DataB());
   double Bmin(PofB.GetBmin());
 
   char debugfile[50];
-  int n = sprintf (debugfile, "test_PB_%f.dat", Bmin);
+  int n = sprintf (debugfile, "test_PB_%f_%ld.dat", Bmin, seconds);
 
   if (n > 0) {
     ofstream of(debugfile);
diff --git a/src/external/TFitPofB-lib/classes/TTrimSPDataHandler.cpp b/src/external/TFitPofB-lib/classes/TTrimSPDataHandler.cpp
index ec4f6de2..59695b57 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/24
+  2008/05/25
 
 ***************************************************************************/
 
@@ -65,6 +65,8 @@ TTrimSPData::TTrimSPData(const string &path, vector<string> &energyVec) {
 
     }
   }
+
+  fOrigDataNZ = fDataNZ;
 }
 
 //---------------------
@@ -86,7 +88,7 @@ vector<double> TTrimSPData::DataZ(double e) const {
 }
 
 //---------------------
-// Method returning n(z)-vector calculated by trim.SP for given energy[keV]
+// Method returning actual n(z)-vector calculated by trim.SP and potentially altered by the WeightLayers-method for given energy[keV]
 //---------------------
 
 vector<double> TTrimSPData::DataNZ(double e) const {
@@ -102,6 +104,79 @@ vector<double> TTrimSPData::DataNZ(double e) const {
 
 }
 
+//---------------------
+// Method returning original n(z)-vector calculated by trim.SP for given energy[keV]
+//---------------------
+
+vector<double> TTrimSPData::OrigDataNZ(double e) const {
+
+  for(unsigned int i(0); i<fEnergy.size(); i++) {
+    if(!(fEnergy[i] - e)) {
+      return fOrigDataNZ[i];
+    }
+  }
+  // default
+  cout << "No implantation profile available for the specified energy... You get back the first one." << endl;
+  return fOrigDataNZ[0];
+
+}
+
+//---------------------
+// Method putting different weight to different layers of your thin film
+// Parameters: Implantation Energy[keV], Interfaces[nm], Weights [0.0 <= w[i] <= 1.0]
+// Example: 25.0, (50, 100), (1.0, 0.33, 1.0)
+//          at 25keV consider 3 layers, where the first ends after 50nm, the second after 100nm (these are NOT the layer thicknesses!!)
+//          the first and last layers get the full n(z), where only one third of the muons in the second layer will be taken into account
+//---------------------
+
+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;
+    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;
+      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;
+        return;
+      }
+    }
+  }
+
+  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;
+      return;
+    }
+  }
+
+  for(unsigned int i(0); i<fEnergy.size(); i++) {
+    if(!(fEnergy[i] - e)) {
+      unsigned int k(0);
+      for(unsigned int j(0); j<fDataZ[i].size(); j++) {
+        if(k<interface.size()) {
+          if(fDataZ[i][j] < interface[k]*10.0)
+            fDataNZ[i][j] = fOrigDataNZ[i][j]*weight[k];
+          else {
+            k++;
+            fDataNZ[i][j] = fOrigDataNZ[i][j]*weight[k];
+          }
+        }
+        else
+          fDataNZ[i][j] = fOrigDataNZ[i][j]*weight[k];
+      }
+    }
+  }
+
+}
+
 //---------------------
 // Method returning n(z) for given z[nm] and energy[keV]
 //---------------------
diff --git a/src/external/TFitPofB-lib/classes/TFitPofB.cpp b/src/external/TFitPofB-lib/classes/TUserLondon.cpp
similarity index 70%
rename from src/external/TFitPofB-lib/classes/TFitPofB.cpp
rename to src/external/TFitPofB-lib/classes/TUserLondon.cpp
index 653bcac5..eb325573 100644
--- a/src/external/TFitPofB-lib/classes/TFitPofB.cpp
+++ b/src/external/TFitPofB-lib/classes/TUserLondon.cpp
@@ -1,24 +1,25 @@
 /***************************************************************************
 
-  TFitPofB.cpp
+  TUserLondon.cpp
 
   Author: Bastian M. Wojek
   e-mail: bastian.wojek@psi.ch
 
-  2008/05/24
+  2008/05/25
 
 ***************************************************************************/
 
-#include "TFitPofB.h"
+#include "TUserLondon.h"
 #include <iostream>
 using namespace std;
 
 //------------------
-// Constructor of the TFitPofB class -- reading available implantation profiles and
+// Constructor of the TUserLondon class -- reading available implantation profiles and
 // creates (a pointer to) the TPofTCalc object (with the FFT plan)
 //------------------
 
-TFitPofB::TFitPofB(const vector<unsigned int> &parNo, const vector<double> &par) : fCalcNeeded(true) {
+TUserLondon::TUserLondon(const vector<unsigned int> &parNo, const vector<double> &par) 
+  : fCalcNeeded(true), fLastTwoChanged(true), fLastThreeChanged(true) {
 
   for(unsigned int i(0); i<parNo.size(); i++) {
     fPar.push_back(par[parNo[i]-1]);
@@ -47,10 +48,10 @@ TFitPofB::TFitPofB(const vector<unsigned int> &parNo, const vector<double> &par)
 }
 
 //------------------
-// Destructor of the TFitPofB class -- cleaning up
+// Destructor of the TUserLondon class -- cleaning up
 //------------------
 
-TFitPofB::~TFitPofB() {
+TUserLondon::~TUserLondon() {
     fPar.clear();
     delete fImpProfile;
     fImpProfile = 0;
@@ -61,10 +62,10 @@ TFitPofB::~TFitPofB() {
 //------------------
 // 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 TFitPofB
+// Parameters: all the parameters for the function to be fitted through TUserLondon
 //------------------
 
-double TFitPofB::Eval(double t, const vector<double> &par) const {
+double TUserLondon::Eval(double t, const vector<double> &par) const {
 
   // check if any parameter has changed
 
@@ -77,12 +78,18 @@ double TFitPofB::Eval(double t, const vector<double> &par) const {
       par_changed = true;
       if(i == 0 || i == 2 || i == 3 || i == 4) {
         cout << "You are varying the model-flag, dt, dB or E! These parameters have to be fixed! Quitting..." << endl; 
-        exit(0);
+        exit(-1);
       } else if (i == 1) {
         only_phase_changed = true;
       } else {
         only_phase_changed = false;
       }
+      if(i == fPar.size()-3) {
+        fLastThreeChanged = true; 
+      } else if (i == fPar.size()-2 || i == fPar.size()-1) {
+        fLastTwoChanged = true;
+        fLastThreeChanged = true;
+      }
     }
   }
 
@@ -140,6 +147,17 @@ double TFitPofB::Eval(double t, const vector<double> &par) const {
         case 2:
         {
 //          cout << "Found the 1D-London model.2L" << endl;
+          if(fLastTwoChanged) {
+            vector<double> interfaces;
+            interfaces.push_back(par[6]+par[7]);
+
+            vector<double> weights;
+            for(unsigned int i(11); i<par.size(); i++)
+              weights.push_back(par[i]);
+
+            cout << "Weighting has changed, re-calculating n(z) now..." << endl;
+            fImpProfile->WeightLayers(par[4], interfaces, weights);
+          }
           TLondon1D_2L BofZ2(par_for_BofZ);
           TPofBCalc PofB2(BofZ2, *fImpProfile, par_for_PofB);
           fPofT->DoFFT(PofB2);
@@ -148,6 +166,18 @@ double TFitPofB::Eval(double t, const vector<double> &par) const {
         case 3:
         {
 //          cout << "Found the 1D-London model.3L" << endl;
+          if(fLastThreeChanged) {
+            vector<double> interfaces;
+            interfaces.push_back(par[6]+par[7]);
+            interfaces.push_back(par[6]+par[7]+par[8]);
+
+            vector<double> weights;
+            for(unsigned int i(12); i<par.size(); i++)
+              weights.push_back(par[i]);
+
+            cout << "Weighting has changed, re-calculating n(z) now..." << endl;
+            fImpProfile->WeightLayers(par[4], interfaces, weights);
+          }
           TLondon1D_3L BofZ3(par_for_BofZ);
           TPofBCalc PofB3(BofZ3, *fImpProfile, par_for_PofB);
           fPofT->DoFFT(PofB3);
@@ -155,7 +185,7 @@ double TFitPofB::Eval(double t, const vector<double> &par) const {
           break;
         default:
           cout << "The user function you specified with the first parameter, does not exist. Quitting..." << endl;
-          exit(0);
+          exit(-1);
       }
     } else {
       cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl;
@@ -164,6 +194,8 @@ double TFitPofB::Eval(double t, const vector<double> &par) const {
     fPofT->CalcPol(par_for_PofT);
 
     fCalcNeeded = false;
+    fLastTwoChanged = false;
+    fLastThreeChanged = false;
 
   }
 
diff --git a/src/external/TFitPofB-lib/include/TPofBCalc.h b/src/external/TFitPofB-lib/include/TPofBCalc.h
index 8730706e..e59e41c8 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/24
+  2008/05/25
 
 ***************************************************************************/
 
@@ -15,6 +15,8 @@
 #include "TBofZCalc.h"
 #include "TTrimSPDataHandler.h"
 
+#define gBar 0.0135538817
+
 class TPofBCalc {
 
 public:
@@ -35,7 +37,6 @@ private:
   vector<double> fPB;
   double fBmin;
   double fBmax;
-  static const double gBar = 0.0135538817;
 
 };
 
diff --git a/src/external/TFitPofB-lib/include/TPofTCalc.h b/src/external/TFitPofB-lib/include/TPofTCalc.h
index d48e35ed..9420bb8b 100644
--- a/src/external/TFitPofB-lib/include/TPofTCalc.h
+++ b/src/external/TFitPofB-lib/include/TPofTCalc.h
@@ -5,7 +5,7 @@
   Author: Bastian M. Wojek
   e-mail: bastian.wojek@psi.ch
 
-  2008/05/24
+  2008/05/25
 
 ***************************************************************************/
 
@@ -15,6 +15,8 @@
 #include "TPofBCalc.h"
 #include "fftw3.h"
 
+#define PI 3.14159265358979323846
+
 class TPofTCalc {
 
 public:
@@ -35,8 +37,6 @@ private:
   vector<double> fPT;
   double fTBin;
   unsigned int fNFFT;
-  static const double PI = 3.14159265358979323846;
-  static const double gBar = 0.0135538817;
 
 };
 
diff --git a/src/external/TFitPofB-lib/include/TTrimSPDataHandler.h b/src/external/TFitPofB-lib/include/TTrimSPDataHandler.h
index 6510e2c9..ea43f0c2 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/24
+  2008/05/25
 
 ***************************************************************************/
 
@@ -30,12 +30,15 @@ public:
   vector<double> Energy() const {return fEnergy;}
   vector<double> DataZ(double) const;
   vector<double> DataNZ(double) const;
+  vector<double> OrigDataNZ(double) const;
+  void WeightLayers(double, const vector<double>&, const vector<double>&) const;
   double GetNofZ(double, double) const;
 
 private:
   vector<double> fEnergy;
   vector< vector<double> > fDataZ;
-  vector< vector<double> > fDataNZ;
+  mutable vector< vector<double> > fDataNZ;
+  vector< vector<double> > fOrigDataNZ;
 };
 
 #endif // _TTrimSPDataHandler_H_
diff --git a/src/external/TFitPofB-lib/include/TFitPofB.h b/src/external/TFitPofB-lib/include/TUserLondon.h
similarity index 60%
rename from src/external/TFitPofB-lib/include/TFitPofB.h
rename to src/external/TFitPofB-lib/include/TUserLondon.h
index a968ec65..1a73ca03 100644
--- a/src/external/TFitPofB-lib/include/TFitPofB.h
+++ b/src/external/TFitPofB-lib/include/TUserLondon.h
@@ -1,24 +1,24 @@
 /***************************************************************************
 
-  TFitPofB.h
+  TUserLondon.h
 
   Author: Bastian M. Wojek
   e-mail: bastian.wojek@psi.ch
 
-  2008/05/24
+  2008/05/25
 
 ***************************************************************************/
 
-#ifndef _TFitPofB_H_
-#define _TFitPofB_H_
+#ifndef _TUserLondon_H_
+#define _TUserLondon_H_
 
 #include "TPofTCalc.h"
 
-class TFitPofB {
+class TUserLondon {
 
 public:
-  TFitPofB(const vector<unsigned int>& , const vector<double>&);
-  ~TFitPofB();
+  TUserLondon(const vector<unsigned int>& , const vector<double>&);
+  ~TUserLondon();
 
   double Eval(double, const vector<double>&) const;
 
@@ -27,7 +27,9 @@ private:
   TTrimSPData *fImpProfile;
   TPofTCalc *fPofT;
   mutable bool fCalcNeeded;
-  
+  mutable bool fLastTwoChanged;
+  mutable bool fLastThreeChanged;
+
 };
 
-#endif //_TFitPofB_H_
+#endif //_TUserLondon_H_
diff --git a/src/external/TFitPofB-lib/test/test.cpp b/src/external/TFitPofB-lib/test/test.cpp
index 3d39b88f..c191843e 100644
--- a/src/external/TFitPofB-lib/test/test.cpp
+++ b/src/external/TFitPofB-lib/test/test.cpp
@@ -1,4 +1,4 @@
-#include "TFitPofB.h"
+#include "TUserLondon.h"
 #include <iostream>
 #include <fstream>
 
@@ -101,8 +101,8 @@ int main(){
   
   
 */
-  unsigned int parNo_arr[] = {1, 3, 5, 7, 9, 11, 12, 13, 14, 15, 16, 17};
-  double par_arr[] = {3.0, 999.0, 0.0, 999.0, 0.01, 999.0, 0.05, 999.0, 20.0, 999.0, 100.0, 10.0, 80.0, 50.0, 80.0, 180.0, 500.0};
+  unsigned int parNo_arr[] = {1, 3, 5, 7, 9, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20};
+  double par_arr[] = {3.0, 999.0, 0.0, 999.0, 0.01, 999.0, 0.05, 999.0, 25.0, 999.0, 100.0, 10.0, 65.0, 50.0, 75.0, 180.0, 500.0, 1.0, 0.3, 1.0};
   
   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 +113,7 @@ int main(){
     par_vec_sub.push_back(par_vec[parNo_vec[i]-1]);
   }
   
-  TFitPofB fitter(parNo_vec, par_vec);
+  TUserLondon fitter(parNo_vec, par_vec);
 
   ofstream of01("test_fitter01.dat");
   ofstream of02("test_fitter02.dat");
@@ -179,7 +179,6 @@ int main(){
   of07.close();
   
     par_vec_sub[1] += 10.0;
-    par_vec_sub[7] = 75.0;
   
   for (double i(0.); i<12.0; i+=0.003) {
     of08 << i << " " << fitter.Eval(i, par_vec_sub) << endl;
@@ -194,7 +193,9 @@ int main(){
   }
   of09.close();
   
-    par_vec_sub[1] += 10.0;
+    par_vec_sub[1] = 0.0;
+    par_vec_sub[11] = 500.0;
+    par_vec_sub[13] = 0.8;
   
   for (double i(0.); i<12.0; i+=0.003) {
     of10 << i << " " << fitter.Eval(i, par_vec_sub) << endl;