LMU/musrfit
LMU/musrfit
5
0
Fork 0
Code Issues 6 Pull Requests Actions Packages Projects Releases 3 Activity

start to work on a more efficient version of MM PDepthProfile user function.

Browse Source
This commit is contained in:
suter_a 2022-12-21 07:55:33 +01:00
parent 7bf2cfd8c1
commit 348f02b217
2 changed files with 141 additions and 117 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
src/external/DepthProfile/inc/PDepthProfile.h vendored
View File

@ -53,6 +53,7 @@ class PDepthProfileGlobal
mutable std::vector<Double_t> fPreviousParam; mutable std::vector<Double_t> fPreviousParam;
PRgeHandler *fRgeHandler{nullptr}; PRgeHandler *fRgeHandler{nullptr};
PRgeDataList fCfd;
ClassDef(PDepthProfileGlobal, 1) ClassDef(PDepthProfileGlobal, 1)
}; };

257
src/external/DepthProfile/src/PDepthProfile.cpp vendored
View File

@ -2,8 +2,8 @@
PDepthProfile.cpp PDepthProfile.cpp
Author: Andreas Suter Authors: Maria Martins, Andreas Suter
e-mail: andreas.suter@psi.ch e-mail: maria.martins@psi.ch, andreas.suter@psi.ch
***************************************************************************/ ***************************************************************************/
@ -46,14 +46,28 @@ ClassImp(PDepthProfileGlobal)
* <p>Constructor. Reads the necessary rge-files based on the depth_profile_startup.xml * <p>Constructor. Reads the necessary rge-files based on the depth_profile_startup.xml
*/ */
PDepthProfileGlobal::PDepthProfileGlobal() { PDepthProfileGlobal::PDepthProfileGlobal() {
// load all the TRIM.SP rge-files // load all the TRIM.SP rge-files
fRgeHandler = new PRgeHandler("./depth_profile_startup.xml"); fRgeHandler = new PRgeHandler("./depth_profile_startup.xml");
if (!fRgeHandler->IsValid()) { if (!fRgeHandler->IsValid()) {
std::cout << std::endl << ">> PDepthProfileGlobal::PDepthProfileGlobal **PANIC ERROR**"; 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 << ">> rge data handler too unhappy. Will terminate unfriendly, sorry.";
std::cout << std::endl; std::cout << std::endl;
fValid = false; fValid = false;
}
// calculate cumulative frequency distribution of all the rge-files
PRgeDataList rgeData = fRgeHandler->GetRgeData();
fCfd.resize(fRgeHandler->GetNoOfRgeDataSets());
for (unsigned int i=0; i<fCfd.size(); i++) {
fCfd[i].energy = rgeData[i].energy;
fCfd[i].depth = rgeData[i].depth;
fCfd[i].nn.resize(rgeData[i].nn.size());
double dval=0.0;
for (unsigned int j=0; j<fCfd[i].nn.size(); j++) {
dval += rgeData[i].nn[j];
fCfd[i].nn[j] = dval;
} }
}
} }
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
@ -63,10 +77,10 @@ PDepthProfileGlobal::PDepthProfileGlobal() {
* <p>Clean up the rge-handler. * <p>Clean up the rge-handler.
*/ */
PDepthProfileGlobal::~PDepthProfileGlobal() { PDepthProfileGlobal::~PDepthProfileGlobal() {
if (fRgeHandler) { if (fRgeHandler) {
delete fRgeHandler; delete fRgeHandler;
fRgeHandler = nullptr; fRgeHandler = nullptr;
} }
} }
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@ -80,10 +94,10 @@ ClassImp(PDepthProfile)
* <p>Clean up the global part. * <p>Clean up the global part.
*/ */
PDepthProfile::~PDepthProfile() { PDepthProfile::~PDepthProfile() {
if ((fDepthProfileGlobal != 0) && fInvokedGlobal) { if ((fDepthProfileGlobal != nullptr) && fInvokedGlobal) {
delete fDepthProfileGlobal; delete fDepthProfileGlobal;
fDepthProfileGlobal = nullptr; fDepthProfileGlobal = nullptr;
} }
} }
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
@ -98,31 +112,31 @@ PDepthProfile::~PDepthProfile() {
* \param idx * \param idx
*/ */
void PDepthProfile::SetGlobalPart(std::vector<void *> &globalPart, UInt_t idx) { void PDepthProfile::SetGlobalPart(std::vector<void *> &globalPart, UInt_t idx) {
fIdxGlobal = static_cast<Int_t>(idx); fIdxGlobal = static_cast<Int_t>(idx);
if ((Int_t) globalPart.size() <= fIdxGlobal) { if ((Int_t) globalPart.size() <= fIdxGlobal) {
fDepthProfileGlobal = new PDepthProfileGlobal(); fDepthProfileGlobal = new PDepthProfileGlobal();
if (fDepthProfileGlobal == nullptr) { if (fDepthProfileGlobal == nullptr) {
fValid = false; fValid = false;
std::cerr << std::endl std::cerr << std::endl
<< ">> PDepthProfile::SetGlobalPart(): **ERROR** Couldn't invoke global user function object, sorry ..." << ">> PDepthProfile::SetGlobalPart(): **ERROR** Couldn't invoke global user function object, sorry ..."
<< std::endl; << std::endl;
} else if (!fDepthProfileGlobal->IsValid()) { } else if (!fDepthProfileGlobal->IsValid()) {
fValid = false; fValid = false;
std::cerr << std::endl std::cerr << std::endl
<< ">> PDepthProfile::SetGlobalPart(): **ERROR** initialization of global user function object failed, sorry ..." << ">> PDepthProfile::SetGlobalPart(): **ERROR** initialization of global user function object failed, sorry ..."
<< std::endl; << std::endl;
} else {
fValid = true;
fInvokedGlobal = true;
globalPart.resize(fIdxGlobal + 1);
globalPart[fIdxGlobal] = dynamic_cast<PDepthProfileGlobal *>(fDepthProfileGlobal);
}
} else { } else {
fValid = true; fValid = true;
fDepthProfileGlobal = (PDepthProfileGlobal * ) fInvokedGlobal = true;
globalPart[fIdxGlobal]; globalPart.resize(fIdxGlobal + 1);
globalPart[fIdxGlobal] = dynamic_cast<PDepthProfileGlobal *>(fDepthProfileGlobal);
} }
} else {
fValid = true;
fDepthProfileGlobal = (PDepthProfileGlobal * )
globalPart[fIdxGlobal];
}
} }
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
@ -134,13 +148,19 @@ void PDepthProfile::SetGlobalPart(std::vector<void *> &globalPart, UInt_t idx) {
* <b>return:</b> * <b>return:</b>
*/ */
Bool_t PDepthProfile::GlobalPartIsValid() const { Bool_t PDepthProfile::GlobalPartIsValid() const {
return (fValid && fDepthProfileGlobal->IsValid()); return (fValid && fDepthProfileGlobal->IsValid());
} }
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
// GetStoppingProbability() // GetStoppingProbability()
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
/**
* @brief PDepthProfileGlobal::GetStoppingProbability
* @param a
* @param b
* @param energy
* @return
*/
double PDepthProfileGlobal::GetStoppingProbability(double a, double b, Double_t energy) const { double PDepthProfileGlobal::GetStoppingProbability(double a, double b, Double_t energy) const {
// calculation of stopping probability for a given z interval and experimental energy // calculation of stopping probability for a given z interval and experimental energy
@ -175,94 +195,97 @@ double PDepthProfileGlobal::GetStoppingProbability(double a, double b, Double_t
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
// operator() // operator()
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
/**
* @brief PDepthProfile::operator ()
* @param t
* @param param
* @return
*/
Double_t PDepthProfile::operator()(Double_t t, const std::vector <Double_t> &param) const { Double_t PDepthProfile::operator()(Double_t t, const std::vector <Double_t> &param) const {
//verify number of parameters: 2n+1 // verify number of parameters: 2n+1, i.e. it has to be an odd number of parameters
// parameters: {E,f1, f2, ..., f_n, x1, ..., x_(n-1)} // parameters: {f1, f2, ..., f_n, x1, ..., x_(n-1)}
assert(param.size() > 2); assert(param.size() > 2);
assert(((param.size() - 1) % 2) == 0); assert(((param.size() - 1) % 2) == 0);
//number of steps: n+1
int n = (param.size() - 1) / 2;
std::vector<double> parameters;
if (fPreviousParam.size() == 0) { //number of steps: n+1
for (int i = 0; i < param.size(); i++) { int n = (param.size() - 1) / 2;
fPreviousParam.push_back(param[i]); 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]); 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 { } else {
for (int i = n + 1; i < param.size() ; i++) { if (param[i] < param[i-1]) {
//std::cout << "param[i]: " << param[i] << "param[i+1]: " << param[i + 1] << std::endl; parameters.push_back(fPreviousParam[i-1]);
if (i!=(param.size() -1)){ } else {
if (param[i] > param[i + 1]) { //std::cout<<"bem" << std::endl;
parameters.push_back(fPreviousParam[i+1]); parameters.push_back(param[i]);
} 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;
} }
}
} }
for (int i = 0; i < param.size(); i++) {
Double_t energy = t; fPreviousParam.push_back(parameters[i]);
std::vector<double> probability; std::cout << "param[i]: " << fPreviousParam[i] << std::endl;
// 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; Double_t energy = t;
std::vector<double> probability;
for (int j = 0; j < n + 1; j++) { // get stopping probability
fit = fit + parameters[j] * probability[j]; //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);
}
/* std::cout << "Energy " << energy << std::endl; double fit=0;
std::cout << "FRACTION " << fit << std::endl;*/
return fit;
for (int j=0; j<n+1; j++) {
fit = fit + parameters[j] * probability[j];
}
return fit;
} }