some first bug-fixing and cleanup of the code.

src/external/DepthProfile/inc/PDepthProfile.h

@@ -44,8 +44,8 @@ class PDepthProfileGlobal
     Bool_t IsValid() { return fValid; }
     virtual Int_t GetEnergyIndex(const Double_t energy) { return fRgeHandler->GetEnergyIndex(energy); }
     virtual Double_t GetMuonStoppingDensity(const Int_t energyIndex, const Double_t z) const { return fRgeHandler->Get_n(energyIndex, z); }
-    virtual double GetStoppingProbability(double a, double b, Double_t energy) const;
+    virtual Double_t GetStoppingProbability(double a, double b, Double_t energy) const;
-    virtual double GetZmax(Double_t energy) const { return fRgeHandler->GetZmax(energy); }
+    virtual Double_t GetZmax(const Double_t energy);
 
   private:
     Bool_t fValid{true};

src/external/DepthProfile/src/PDepthProfile.cpp

@@ -101,15 +101,15 @@ PDepthProfile::~PDepthProfile() {
 }
 
 //--------------------------------------------------------------------------
-// SetGlobalPart (public)
+// SetGlobalPart : public
 //--------------------------------------------------------------------------
 /**
  * <p>
  *
- * <b>return:</b>
+ * @param globalPart
+ * @param idx
  *
- * \param globalPart
+ * @return
- * \param idx
  */
 void PDepthProfile::SetGlobalPart(std::vector<void *> &globalPart, UInt_t idx) {
   fIdxGlobal = static_cast<Int_t>(idx);
@@ -140,22 +140,23 @@ void PDepthProfile::SetGlobalPart(std::vector<void *> &globalPart, UInt_t idx) {
 }
 
 //--------------------------------------------------------------------------
-// GlobalPartIsValid (public)
+// GlobalPartIsValid : public
 //--------------------------------------------------------------------------
 /**
- * <p>
+ * <p>Returns true if the global part is valid. It is also checking if all
+ * muon stopping profiles have been loaded, if present.
  *
- * <b>return:</b>
+ * @return true if the global part is valid, false otherwise.
  */
 Bool_t PDepthProfile::GlobalPartIsValid() const {
   return (fValid && fDepthProfileGlobal->IsValid());
 }
 
 //--------------------------------------------------------------------------
-// GetPosIdx()
+// GetPosIdx() : private
 //--------------------------------------------------------------------------
 /**
- * <p>Get the position index for a given distance 'a' and a given implantation
+ * <p>Get the position index for a given distance 'z' and a given implantation
  * energy index.
  *
  * @param z distance for which the index is requested
@@ -165,7 +166,7 @@ Bool_t PDepthProfile::GlobalPartIsValid() const {
  */
 Int_t PDepthProfileGlobal::GetPosIdx(const double z, const Int_t energyIdx) const
 {
-  Int_t idx=0;
+  Int_t idx=-1;
 
   for (UInt_t i=0; i<fCfd[energyIdx].depth.size(); i++) {
     if (z <= fCfd[energyIdx].depth[i]) {
@@ -173,12 +174,14 @@ Int_t PDepthProfileGlobal::GetPosIdx(const double z, const Int_t energyIdx) cons
       break;
     }
   }
+  if (idx==-1) // z larger than any found 'z' value, hence put it to zMax
+    idx = fCfd[energyIdx].depth[fCfd[energyIdx].depth.size()-1];
 
   return idx;
 }
 
 //--------------------------------------------------------------------------
-// GetStoppingProbability()
+// GetStoppingProbability() : public
 //--------------------------------------------------------------------------
 /**
  * <p>Calculates the stopping probability from a to b for a given implantation
@@ -188,13 +191,12 @@ Int_t PDepthProfileGlobal::GetPosIdx(const double z, const Int_t energyIdx) cons
  *
  * @param a starting point in distance
  * @param b stopping point in distance
- * @param energy implantation energy in keV
+ * @param energy implantation energy in (eV)
  *
  * @return stopping probability p.
  */
-double PDepthProfileGlobal::GetStoppingProbability(double a, double b, Double_t energy) const {
+Double_t PDepthProfileGlobal::GetStoppingProbability(double a, double b, Double_t energy) const {
 
-  energy = energy * 1000;
   Int_t idx = fRgeHandler->GetEnergyIndex(energy);
   if (idx == -1) {
     std::cerr << "**WARNING** couldn't find energy " << energy << " (eV) in the rge-files provided." << std::endl;
@@ -205,37 +207,30 @@ double PDepthProfileGlobal::GetStoppingProbability(double a, double b, Double_t
   const Int_t idx_b = GetPosIdx(b, idx); // get index for distance b from cfd(z, E)
 
   return fCfd[idx].nn[idx_b] - fCfd[idx].nn[idx_a];
+}
 
-/* //as35
+//--------------------------------------------------------------------------
-    double zMax = fRgeHandler->GetZmax(energy);
+// GetZmax() : public
+//--------------------------------------------------------------------------
+/**
+ * <p>Returns the maximal distance present in an rge-dataset for a given energy.
+ *
+ * @param energy in (eV) of the rge-dataset
+ *
+ * @return maximal distance in the rge-dataset, or 1e4 if energy is not found.
+ */
+Double_t PDepthProfileGlobal::GetZmax(const Double_t energy)
+{
+  const Int_t idx=GetEnergyIndex(energy);
+  if (idx == -1) return 10000.0;
+  const Int_t maxDepthIdx = fCfd[idx].depth.size()-1;
 
-    std::vector<double> z;
+  return fCfd[idx].depth[maxDepthIdx];
-    double x = a;
-    double xMax = b;
-    std::cout  << " a: " << x << "xmax: " << xMax<< std::endl;
-
-    while (x < xMax) {
-        //z[j]=x;
-        z.push_back(x);
-        x++;
-    }
-
-    double probability = 0;
-    double prob = 0;
-    for (int j = a; j < b - 1; j++) {
-        prob = fRgeHandler->Get_n(energy, j);
-        if ((j == a) || (j == b - 1)) {
-        }
-        probability = probability + prob;
-    }
-
-    return probability;
-*/ //as35
 }
 
 
 //--------------------------------------------------------------------------
-// operator()
+// operator() : public
 //--------------------------------------------------------------------------
 /**
  * <p> calculate the stopping fraction amplitude f(E) for the given parameter:
@@ -256,6 +251,8 @@ Double_t PDepthProfile::operator()(Double_t energy, const std::vector <Double_t>
   // number of steps: n+1
   int n = (param.size() - 1) / 2;
 
+  energy *= 1000.0; // keV -> eV
+
   double fE = param[0]*fDepthProfileGlobal->GetStoppingProbability(0.0, param[n+1], energy); // z=0..x_1
 
   for (int i=n+1; i<param.size()-1; i++) {
@@ -265,83 +262,6 @@ Double_t PDepthProfile::operator()(Double_t energy, const std::vector <Double_t>
   double zMax = fDepthProfileGlobal->GetZmax(energy);
   fE += param[n]*fDepthProfileGlobal->GetStoppingProbability(param[param.size()-1], zMax, energy); // z=x_(n-1)..infinity
 
-/* //as35
-  std::vector<double> parameters;
-
-  if (fPreviousParam.size() == 0) {
-    for (int i = 0; i < param.size(); i++) {
-      fPreviousParam.push_back(param[i]);
-      parameters.push_back(param[i]);
-    }
-  } else {
-    for (int i = 0;  i <(n+1); i++) {
-      if (param[i]>=0 & param[i] <=1) {
-        parameters.push_back(param[i]);
-      } else{
-        parameters.push_back(fPreviousParam[i]);
-      }
-    }
-    if (n==1) {
-      if (param[n + 1]>0) {
-        parameters.push_back(param[n+1]);
-      } else {
-        parameters.push_back(fPreviousParam[n+1]);
-      }
-    } else {
-      for (int i=n+1; i<param.size(); i++) {
-        //std::cout << "param[i]: " << param[i] << "param[i+1]: " << param[i + 1] << std::endl;
-        if (i != (param.size() -1)){
-          if (param[i] > param[i+1]) {
-            parameters.push_back(fPreviousParam[i+1]);
-          } else {
-            //std::cout<<"bem" << std::endl;
-            parameters.push_back(param[i]);
-          }
-        } else {
-          if (param[i] < param[i-1]) {
-            parameters.push_back(fPreviousParam[i-1]);
-          } else {
-            //std::cout<<"bem" << std::endl;
-            parameters.push_back(param[i]);
-          }
-        }
-      }
-    }
-    for (int i = 0; i < param.size(); i++) {
-      fPreviousParam.push_back(parameters[i]);
-      std::cout << "param[i]: " << fPreviousParam[i] << std::endl;
-    }
-  }
-
-  std::vector<double> probability;
-
-  // get stopping probability
-  //int l=0;
-  for (int i=0; i<n+1; i++) {
-    double probE;
-    double a;
-    double b;
-    if (i == 0) { //prob between 0 and x_1
-      a = 0;
-      b = parameters[n + 1];
-    } else if (i == n) { //prob between x_(n-1) and inf
-      a = parameters[2 * n];
-      b = 179;
-    } else { //prob between x_1-x_2, ..., x_n-x_(n-1)
-      a = parameters[n+i];
-      b = parameters[n+1+i];
-    }
-    probE = fDepthProfileGlobal->GetStoppingProbability(a, b, energy);
-    probability.push_back(probE);
-  }
-
-  double fit=0;
-
-  for (int j=0; j<n+1; j++) {
-    fit = fit + parameters[j] * probability[j];
-  }
-*/ //as35
-
   return fE;
 }