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

implementation of depth step fits

Browse Source
This commit is contained in:
martins_m
2022-08-19 17:50:25 +02:00
parent 26c1f49554
commit 8fed497adc
137 changed files with 19332 additions and 57 deletions
Download Patch File Download Diff File Expand all files Collapse all files

13
src/classes/PRgeHandler.cpp
View File

@@ -489,15 +489,16 @@ Double_t PRgeHandler::GetZmax(const Double_t energy)
{
int idx=-1;
for (int i=0; i<fData.size(); i++) {
if (fabs(fData[i].energy-energy) < 1.0) {
if (fabs(fData[i].energy-energy)*0.001 < 0.9){
idx = i;
break;
}
}
}
if (idx != -1)
return GetZmax(idx);
return GetZmax(idx);
return -1.0;
return -1;
}
//--------------------------------------------------------------------------
@@ -531,7 +532,7 @@ Double_t PRgeHandler::Get_n(const Double_t energy, const Double_t z)
{
int idx=-1;
for (int i=0; i<fData.size(); i++) {
if (fabs(fData[i].energy-energy) < 1.0) {
if (fabs(fData[i].energy-energy)*0.001 < 0.90) {
idx = i;
break;
}
@@ -593,7 +594,7 @@ Int_t PRgeHandler::GetEnergyIndex(const Double_t energy)
{
int idx=-1;
for (int i=0; i<fData.size(); i++) {
if (fabs(fData[i].energy-energy) < 1.0) {
if (fabs(fData[i].energy-energy)*0.001 < 0.90) {
idx = i;
break;
}

37
src/external/DepthProfile/depth_profile_startup.xml vendored Normal file
View File

@@ -0,0 +1,37 @@
<?xml version="1.0" encoding="UTF-8"?>
<nonlocal xmlns="http://nemu.web.psi.ch/musrfit/nonlocal">
<comment>
TrimSp information
</comment>
<trim_sp>
<data_path>./Trimsp/Si10/</data_path>
<rge_fln_pre>Si10_100nm_2.2_E</rge_fln_pre>
<energy_list>
<energy>1000</energy>
<energy>2000</energy>
<energy>3000</energy>
<energy>4000</energy>
<energy>5000</energy>
<energy>6000</energy>
<energy>7000</energy>
<energy>8000</energy>
<energy>9000</energy>
<energy>10000</energy>
<energy>11000</energy>
<energy>12000</energy>
<energy>13000</energy>
<energy>14000</energy>
<energy>15000</energy>
<energy>16000</energy>
<energy>17000</energy>
<energy>17000</energy>
<energy>18000</energy>
<energy>19000</energy>
<energy>20000</energy>
<energy>21000</energy>
<energy>22000</energy>
<energy>23000</energy>
<energy>24000</energy>
</energy_list>
</trim_sp>
</nonlocal>

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

@@ -44,9 +44,11 @@ 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;
private:
Bool_t fValid{true};
mutable std::vector<Double_t> fPreviousParam;
@@ -72,6 +74,7 @@ class PDepthProfile : public PUserFcnBase
Bool_t fInvokedGlobal{false};
Int_t fIdxGlobal;
mutable std::vector<Double_t> fPreviousParam;
PDepthProfileGlobal *fDepthProfileGlobal{nullptr};
// definition of the class for the ROOT dictionary

8
src/external/DepthProfile/src/.idea/.gitignore generated vendored Normal file
View File

@@ -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

1
src/external/DepthProfile/src/.idea/.name generated vendored Normal file
View File

@@ -0,0 +1 @@
PDepthProfile.cpp

8
src/external/DepthProfile/src/.idea/modules.xml generated vendored Normal file
View File

@@ -0,0 +1,8 @@
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ProjectModuleManager">
<modules>
<module fileurl="file://$PROJECT_DIR$/.idea/src.iml" filepath="$PROJECT_DIR$/.idea/src.iml" />
</modules>
</component>
</project>

8
src/external/DepthProfile/src/.idea/src.iml generated vendored Normal file
View File

@@ -0,0 +1,8 @@
<?xml version="1.0" encoding="UTF-8"?>
<module type="CPP_MODULE" version="4">
<component name="NewModuleRootManager">
<content url="file://$MODULE_DIR$" />
<orderEntry type="inheritedJdk" />
<orderEntry type="sourceFolder" forTests="false" />
</component>
</module>

2
src/external/DepthProfile/src/CMakeLists.txt vendored
View File

@@ -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)

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

@@ -45,16 +45,15 @@ ClassImp(PDepthProfileGlobal)
/**
* <p>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()
/**
* <p>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)
/**
* <p>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<void*> &globalPart, UInt_t idx)
{
fIdxGlobal = static_cast<Int_t>(idx);
void PDepthProfile::SetGlobalPart(std::vector<void *> &globalPart, UInt_t idx) {
fIdxGlobal = static_cast<Int_t>(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<PDepthProfileGlobal *>(fDepthProfileGlobal);
}
} else {
fValid = true;
fInvokedGlobal = true;
globalPart.resize(fIdxGlobal+1);
globalPart[fIdxGlobal] = dynamic_cast<PDepthProfileGlobal*>(fDepthProfileGlobal);
fValid = true;
fDepthProfileGlobal = (PDepthProfileGlobal * )
globalPart[fIdxGlobal];
}
} else {
fValid = true;
fDepthProfileGlobal = (PDepthProfileGlobal*)globalPart[fIdxGlobal];
}
}
//--------------------------------------------------------------------------
@@ -132,19 +133,136 @@ void PDepthProfile::SetGlobalPart(std::vector<void*> &globalPart, UInt_t idx)
*
* <b>return:</b>
*/
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<double> 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()
//--------------------------------------------------------------------------
/**
* <p>Yet to be implemented.
*/
Double_t PDepthProfile::operator()(Double_t t, const std::vector<Double_t> &param) const
{
return 0.0;
Double_t PDepthProfile::operator()(Double_t t, const std::vector <Double_t> &param) 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<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;
}
}
Double_t energy = t;
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];
}
/* std::cout << "Energy " << energy << std::endl;
std::cout << "FRACTION " << fit << std::endl;*/
return fit;
}

247
src/external/DepthProfile/src/PDepthProfile.cpp.backup vendored Normal file
View File

@@ -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 <cassert>
#include <cmath>
#include <iostream>
#include <vector>
#include <TSAXParser.h>
#include <TMath.h>
#include "PDepthProfile.h"
ClassImp(PDepthProfileGlobal)
//--------------------------------------------------------------------------
// Constructor (PDepthProfileGlobal)
//--------------------------------------------------------------------------
/**
* <p>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)
//--------------------------------------------------------------------------
/**
* <p>Clean up the rge-handler.
*/
PDepthProfileGlobal::~PDepthProfileGlobal()
{
if (fRgeHandler) {
delete fRgeHandler;
fRgeHandler = nullptr;
}
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
ClassImp(PDepthProfile)
//--------------------------------------------------------------------------
// Destructor (PDepthProfile)
//--------------------------------------------------------------------------
/**
* <p>Clean up the global part.
*/
PDepthProfile::~PDepthProfile()
{
if ((fDepthProfileGlobal != 0) && fInvokedGlobal) {
delete fDepthProfileGlobal;
fDepthProfileGlobal = nullptr;
}
}
//--------------------------------------------------------------------------
// SetGlobalPart (public)
//--------------------------------------------------------------------------
/**
* <p>
*
* <b>return:</b>
*
* \param globalPart
* \param idx
*/
void PDepthProfile::SetGlobalPart(std::vector<void*> &globalPart, UInt_t idx)
{
fIdxGlobal = static_cast<Int_t>(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<PDepthProfileGlobal*>(fDepthProfileGlobal);
}
} else {
fValid = true;
fDepthProfileGlobal = (PDepthProfileGlobal*)globalPart[fIdxGlobal];
}
}
//--------------------------------------------------------------------------
// GlobalPartIsValid (public)
//--------------------------------------------------------------------------
/**
* <p>
*
* <b>return:</b>
*/
Bool_t PDepthProfile::GlobalPartIsValid() const
{
return (fValid && fDepthProfileGlobal->IsValid());
}
//--------------------------------------------------------------------------
// GetStoppingProbability()
//--------------------------------------------------------------------------
/**
* <p>Yet to be implemented.
*/
std::array<double,24> 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<double> z;
double x=a;
double xMax=b;
if (x==0) {
x=0.5;
}
if (b>zMax){
xMax=zMax;
}
int j=0;
while (x<xMax){
z[j]=a;
x++;
j++;
}
vector<double> n;
double prob[24];
double sum=0;
for (int i=0; i<24;i++){
for (int j=0; j<z.size(); j++){
//n[j]=fRgeHandler.Get_n(energy[i], z[j]);
if ((j==0) || (j==z.size()-1)){
prob=0.5*fRgeHandler->Get_n(energy[i], z[j];
} else {
prob=fRgeHandler->Get_n(energy[i], z[j];
}
sum=sum+prob;
}
prob[i]=sum;
}
return prob;
}
//--------------------------------------------------------------------------
// operator()
//--------------------------------------------------------------------------
/**
* <p>Yet to be implemented.
*/
Double_t PDepthProfile::operator()(Double_t t, const std::vector<Double_t> &param) 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;i<sizeof(energy); i++){
energy[i]=(i+1)*1000;
}
double probability [sizeof(energy)][n+1];
// get stopping probability
int l=0;
for (int i=n+1; i<2*n+1; i++){
double prob[24];
if (i==n+1){
prob=GetStoppingProbability(0.0,param(i));
} else if (i==2*n) {
prob=GetStoppingProbability(i,fRgeHandler->GetZmax(energy[23]));
} else {
prob=GetStoppingProbability(i,i+1);
}
for (int k=0;k<sizeof(energy); k++){
probability[k][l]=prob[k];
}
l++;
}
double fit[24];
double intermediate;
for (int i=0; i<sizeof(energy); i++){
intermediate=0;
for (int j=0; j<n+1; j++){
intermediate=intermediate+param[j]*probability[i][j];
}
fit[i]=intermediate;
}
return 0.0;
}