fPreviousParam;
PDepthProfileGlobal *fDepthProfileGlobal{nullptr};
// definition of the class for the ROOT dictionary
diff --git a/src/external/DepthProfile/src/.idea/.gitignore b/src/external/DepthProfile/src/.idea/.gitignore
new file mode 100644
index 00000000..13566b81
--- /dev/null
+++ b/src/external/DepthProfile/src/.idea/.gitignore
@@ -0,0 +1,8 @@
+# Default ignored files
+/shelf/
+/workspace.xml
+# Editor-based HTTP Client requests
+/httpRequests/
+# Datasource local storage ignored files
+/dataSources/
+/dataSources.local.xml
diff --git a/src/external/DepthProfile/src/.idea/.name b/src/external/DepthProfile/src/.idea/.name
new file mode 100644
index 00000000..1084bbed
--- /dev/null
+++ b/src/external/DepthProfile/src/.idea/.name
@@ -0,0 +1 @@
+PDepthProfile.cpp
\ No newline at end of file
diff --git a/src/external/DepthProfile/src/.idea/modules.xml b/src/external/DepthProfile/src/.idea/modules.xml
new file mode 100644
index 00000000..f669a0e5
--- /dev/null
+++ b/src/external/DepthProfile/src/.idea/modules.xml
@@ -0,0 +1,8 @@
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/src/external/DepthProfile/src/.idea/src.iml b/src/external/DepthProfile/src/.idea/src.iml
new file mode 100644
index 00000000..bc2cd874
--- /dev/null
+++ b/src/external/DepthProfile/src/.idea/src.iml
@@ -0,0 +1,8 @@
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/src/external/DepthProfile/src/CMakeLists.txt b/src/external/DepthProfile/src/CMakeLists.txt
index 6aa386c9..9e3562fc 100644
--- a/src/external/DepthProfile/src/CMakeLists.txt
+++ b/src/external/DepthProfile/src/CMakeLists.txt
@@ -50,7 +50,7 @@ target_include_directories(
)
#--- add library dependencies -------------------------------------------------
-target_link_libraries(PDepthProfile ${FFTW3_LIBRARY} ${ROOT_LIBRARIES} PUserFcnBase)
+target_link_libraries(PDepthProfile ${FFTW3_LIBRARY} ${ROOT_LIBRARIES} PRgeHandler PUserFcnBase)
#--- install PDepthProfile solib ----------------------------------------------
install(TARGETS PDepthProfile DESTINATION lib)
diff --git a/src/external/DepthProfile/src/PDepthProfile.cpp b/src/external/DepthProfile/src/PDepthProfile.cpp
index b4d9795f..d0d24ce7 100644
--- a/src/external/DepthProfile/src/PDepthProfile.cpp
+++ b/src/external/DepthProfile/src/PDepthProfile.cpp
@@ -45,16 +45,15 @@ ClassImp(PDepthProfileGlobal)
/**
* Constructor. Reads the necessary rge-files based on the depth_profile_startup.xml
*/
-PDepthProfileGlobal::PDepthProfileGlobal()
-{
- // load all the TRIM.SP rge-files
- fRgeHandler = new PRgeHandler("./depth_profile_startup.xml");
- if (!fRgeHandler->IsValid()) {
- std::cout << std::endl << ">> PDepthProfileGlobal::PDepthProfileGlobal **PANIC ERROR**";
- std::cout << std::endl << ">> rge data handler too unhappy. Will terminate unfriendly, sorry.";
- std::cout << std::endl;
- fValid = false;
- }
+PDepthProfileGlobal::PDepthProfileGlobal() {
+ // load all the TRIM.SP rge-files
+ fRgeHandler = new PRgeHandler("./depth_profile_startup.xml");
+ if (!fRgeHandler->IsValid()) {
+ std::cout << std::endl << ">> PDepthProfileGlobal::PDepthProfileGlobal **PANIC ERROR**";
+ std::cout << std::endl << ">> rge data handler too unhappy. Will terminate unfriendly, sorry.";
+ std::cout << std::endl;
+ fValid = false;
+ }
}
//--------------------------------------------------------------------------
@@ -63,12 +62,11 @@ PDepthProfileGlobal::PDepthProfileGlobal()
/**
*
Clean up the rge-handler.
*/
-PDepthProfileGlobal::~PDepthProfileGlobal()
-{
- if (fRgeHandler) {
- delete fRgeHandler;
- fRgeHandler = nullptr;
- }
+PDepthProfileGlobal::~PDepthProfileGlobal() {
+ if (fRgeHandler) {
+ delete fRgeHandler;
+ fRgeHandler = nullptr;
+ }
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -81,12 +79,11 @@ ClassImp(PDepthProfile)
/**
*
Clean up the global part.
*/
-PDepthProfile::~PDepthProfile()
-{
- if ((fDepthProfileGlobal != 0) && fInvokedGlobal) {
- delete fDepthProfileGlobal;
- fDepthProfileGlobal = nullptr;
- }
+PDepthProfile::~PDepthProfile() {
+ if ((fDepthProfileGlobal != 0) && fInvokedGlobal) {
+ delete fDepthProfileGlobal;
+ fDepthProfileGlobal = nullptr;
+ }
}
//--------------------------------------------------------------------------
@@ -100,28 +97,32 @@ PDepthProfile::~PDepthProfile()
* \param globalPart
* \param idx
*/
-void PDepthProfile::SetGlobalPart(std::vector &globalPart, UInt_t idx)
-{
- fIdxGlobal = static_cast(idx);
+void PDepthProfile::SetGlobalPart(std::vector &globalPart, UInt_t idx) {
+ fIdxGlobal = static_cast(idx);
- if ((Int_t)globalPart.size() <= fIdxGlobal) {
- fDepthProfileGlobal = new PDepthProfileGlobal();
- if (fDepthProfileGlobal == nullptr) {
- fValid = false;
- std::cerr << std::endl << ">> PDepthProfile::SetGlobalPart(): **ERROR** Couldn't invoke global user function object, sorry ..." << std::endl;
- } else if (!fDepthProfileGlobal->IsValid()) {
- fValid = false;
- std::cerr << std::endl << ">> PDepthProfile::SetGlobalPart(): **ERROR** initialization of global user function object failed, sorry ..." << std::endl;
+ if ((Int_t) globalPart.size() <= fIdxGlobal) {
+ fDepthProfileGlobal = new PDepthProfileGlobal();
+ if (fDepthProfileGlobal == nullptr) {
+ fValid = false;
+ std::cerr << std::endl
+ << ">> PDepthProfile::SetGlobalPart(): **ERROR** Couldn't invoke global user function object, sorry ..."
+ << std::endl;
+ } else if (!fDepthProfileGlobal->IsValid()) {
+ fValid = false;
+ std::cerr << std::endl
+ << ">> PDepthProfile::SetGlobalPart(): **ERROR** initialization of global user function object failed, sorry ..."
+ << std::endl;
+ } else {
+ fValid = true;
+ fInvokedGlobal = true;
+ globalPart.resize(fIdxGlobal + 1);
+ globalPart[fIdxGlobal] = dynamic_cast(fDepthProfileGlobal);
+ }
} else {
- fValid = true;
- fInvokedGlobal = true;
- globalPart.resize(fIdxGlobal+1);
- globalPart[fIdxGlobal] = dynamic_cast(fDepthProfileGlobal);
+ fValid = true;
+ fDepthProfileGlobal = (PDepthProfileGlobal * )
+ globalPart[fIdxGlobal];
}
- } else {
- fValid = true;
- fDepthProfileGlobal = (PDepthProfileGlobal*)globalPart[fIdxGlobal];
- }
}
//--------------------------------------------------------------------------
@@ -132,19 +133,136 @@ void PDepthProfile::SetGlobalPart(std::vector &globalPart, UInt_t idx)
*
* return:
*/
-Bool_t PDepthProfile::GlobalPartIsValid() const
-{
- return (fValid && fDepthProfileGlobal->IsValid());
+Bool_t PDepthProfile::GlobalPartIsValid() const {
+ return (fValid && fDepthProfileGlobal->IsValid());
}
+//--------------------------------------------------------------------------
+// GetStoppingProbability()
+//--------------------------------------------------------------------------
+
+double PDepthProfileGlobal::GetStoppingProbability(double a, double b, Double_t energy) const {
+
+ // calculation of stopping probability for a given z interval and experimental energy
+ // \int n(z, E) dz ~ dz [1/2 n_0 + n_1 + n_2 + ... + n_(N-1) + 1/2 n_N]
+ energy = energy * 1000;
+ double zMax = fRgeHandler->GetZmax(energy);
+
+ std::vector z;
+ 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;
+}
+
+
//--------------------------------------------------------------------------
// operator()
//--------------------------------------------------------------------------
-/**
- * Yet to be implemented.
- */
-Double_t PDepthProfile::operator()(Double_t t, const std::vector ¶m) const
-{
- return 0.0;
+
+Double_t PDepthProfile::operator()(Double_t t, const std::vector ¶m) const {
+ //verify number of parameters: 2n+1
+ // parameters: {E,f1, f2, ..., f_n, x1, ..., x_(n-1)}
+ assert(param.size() > 2);
+ assert(((param.size() - 1) % 2) == 0);
+ //number of steps: n+1
+ int n = (param.size() - 1) / 2;
+ std::vector 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;
+ }
+ }
+
+ Double_t energy = t;
+ std::vector 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];
+ }
+
+ /* std::cout << "Energy " << energy << std::endl;
+ std::cout << "FRACTION " << fit << std::endl;*/
+ return fit;
+
}
+
diff --git a/src/external/DepthProfile/src/PDepthProfile.cpp.backup b/src/external/DepthProfile/src/PDepthProfile.cpp.backup
new file mode 100644
index 00000000..80ef1e1d
--- /dev/null
+++ b/src/external/DepthProfile/src/PDepthProfile.cpp.backup
@@ -0,0 +1,247 @@
+/***************************************************************************
+
+ PDepthProfile.cpp
+
+ Author: Andreas Suter
+ e-mail: andreas.suter@psi.ch
+
+***************************************************************************/
+
+/***************************************************************************
+ * Copyright (C) 2009-2022 by Andreas Suter *
+ * andreas.suter@psi.ch *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ***************************************************************************/
+
+#include
+#include
+
+#include
+#include
+
+#include
+#include
+
+#include "PDepthProfile.h"
+
+ClassImp(PDepthProfileGlobal)
+
+//--------------------------------------------------------------------------
+// Constructor (PDepthProfileGlobal)
+//--------------------------------------------------------------------------
+/**
+ * Constructor. Reads the necessary rge-files based on the depth_profile_startup.xml
+ */
+PDepthProfileGlobal::PDepthProfileGlobal()
+{
+ // load all the TRIM.SP rge-files
+ fRgeHandler = new PRgeHandler("./depth_profile_startup.xml");
+ if (!fRgeHandler->IsValid()) {
+ std::cout << std::endl << ">> PDepthProfileGlobal::PDepthProfileGlobal **PANIC ERROR**";
+ std::cout << std::endl << ">> rge data handler too unhappy. Will terminate unfriendly, sorry.";
+ std::cout << std::endl;
+ fValid = false;
+ }
+}
+
+//--------------------------------------------------------------------------
+// Destructor (PDepthProfileGlobal)
+//--------------------------------------------------------------------------
+/**
+ *
Clean up the rge-handler.
+ */
+PDepthProfileGlobal::~PDepthProfileGlobal()
+{
+ if (fRgeHandler) {
+ delete fRgeHandler;
+ fRgeHandler = nullptr;
+ }
+}
+
+//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ClassImp(PDepthProfile)
+
+//--------------------------------------------------------------------------
+// Destructor (PDepthProfile)
+//--------------------------------------------------------------------------
+/**
+ *
Clean up the global part.
+ */
+PDepthProfile::~PDepthProfile()
+{
+ if ((fDepthProfileGlobal != 0) && fInvokedGlobal) {
+ delete fDepthProfileGlobal;
+ fDepthProfileGlobal = nullptr;
+ }
+}
+
+//--------------------------------------------------------------------------
+// SetGlobalPart (public)
+//--------------------------------------------------------------------------
+/**
+ *
+ *
+ * return:
+ *
+ * \param globalPart
+ * \param idx
+ */
+void PDepthProfile::SetGlobalPart(std::vector &globalPart, UInt_t idx)
+{
+ fIdxGlobal = static_cast(idx);
+
+ if ((Int_t)globalPart.size() <= fIdxGlobal) {
+ fDepthProfileGlobal = new PDepthProfileGlobal();
+ if (fDepthProfileGlobal == nullptr) {
+ fValid = false;
+ std::cerr << std::endl << ">> PDepthProfile::SetGlobalPart(): **ERROR** Couldn't invoke global user function object, sorry ..." << std::endl;
+ } else if (!fDepthProfileGlobal->IsValid()) {
+ fValid = false;
+ std::cerr << std::endl << ">> PDepthProfile::SetGlobalPart(): **ERROR** initialization of global user function object failed, sorry ..." << std::endl;
+ } else {
+ fValid = true;
+ fInvokedGlobal = true;
+ globalPart.resize(fIdxGlobal+1);
+ globalPart[fIdxGlobal] = dynamic_cast(fDepthProfileGlobal);
+ }
+ } else {
+ fValid = true;
+ fDepthProfileGlobal = (PDepthProfileGlobal*)globalPart[fIdxGlobal];
+ }
+}
+
+//--------------------------------------------------------------------------
+// GlobalPartIsValid (public)
+//--------------------------------------------------------------------------
+/**
+ *
+ *
+ * return:
+ */
+Bool_t PDepthProfile::GlobalPartIsValid() const
+{
+ return (fValid && fDepthProfileGlobal->IsValid());
+}
+
+//--------------------------------------------------------------------------
+// GetStoppingProbability()
+//--------------------------------------------------------------------------
+/**
+ *
Yet to be implemented.
+ */
+
+std::array PDepthProfile::GetStoppingProbability(double a, double b, double *energy) const
+{
+
+ // calculation of stopping probability for a given interval and experimental energy range
+ // \int n(z, E) dz ~ dz [1/2 n_0 + n_1 + n_2 + ... + n_(N-1) + 1/2 n_N]
+
+
+ double zMax=fRgeHandler->GetZmax(energy[23]);
+ vector z;
+ double x=a;
+ double xMax=b;
+ if (x==0) {
+ x=0.5;
+ }
+ if (b>zMax){
+ xMax=zMax;
+
+ }
+ int j=0;
+
+ while (x n;
+ double prob[24];
+ double sum=0;
+ for (int i=0; i<24;i++){
+ for (int j=0; jGet_n(energy[i], z[j];
+ } else {
+ prob=fRgeHandler->Get_n(energy[i], z[j];
+ }
+ sum=sum+prob;
+ }
+ prob[i]=sum;
+ }
+
+ return prob;
+}
+
+
+//--------------------------------------------------------------------------
+// operator()
+//--------------------------------------------------------------------------
+/**
+ * Yet to be implemented.
+ */
+Double_t PDepthProfile::operator()(Double_t t, const std::vector ¶m) const
+{
+ //verify number of parameters: 2n+1
+ // parameters: {f1, f2, ..., f_n, x1, ..., x_(n-1)}
+ assert(param.size() >2);
+ assert(((param.size()-1)%2) == 0);
+ int n=(param.size()-1)/2;
+
+ double energy[24];
+ for (int i=0;iGetZmax(energy[23]));
+
+ } else {
+ prob=GetStoppingProbability(i,i+1);
+ }
+ for (int k=0;k