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

Removed the strange energy-loops in the TrimSPDataHandler and use generic algorithms instead

Browse Source
This commit is contained in:
Bastian M. Wojek 2009-05-15 13:36:15 +00:00
parent 908cf0a478
commit 3dd8b3e73f
2 changed files with 178 additions and 112 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

285
src/external/TFitPofB-lib/classes/TTrimSPDataHandler.cpp vendored
View File

@ -5,7 +5,7 @@
Author: Bastian M. Wojek Author: Bastian M. Wojek
e-mail: bastian.wojek@psi.ch e-mail: bastian.wojek@psi.ch
2008/11/21 2009/05/15
***************************************************************************/ ***************************************************************************/
@ -15,6 +15,7 @@
#include <string> #include <string>
#include <cmath> #include <cmath>
#include <cassert> #include <cassert>
#include <algorithm>
using namespace std; using namespace std;
@ -33,10 +34,15 @@ using namespace std;
// <some_path>021.rge // <some_path>021.rge
// //
// <some_path>21.rge is explicitly not possible since it is not clear, if this denotes 2.1keV or 21.0keV! // <some_path>21.rge is explicitly not possible since it is not clear, if this denotes 2.1keV or 21.0keV!
//
// Also always use the same format within one energyVec - otherwise sorting of the vector will not work properly!
//-------------------- //--------------------
TTrimSPData::TTrimSPData(const string &path, vector<string> &energyVec) { TTrimSPData::TTrimSPData(const string &path, vector<string> &energyVec) {
// sort the energies in ascending orders - this might be useful for later applications (energy-interpolations etc.)
sort(energyVec.begin(), energyVec.end());
double zz(0.0), nzz(0.0); double zz(0.0), nzz(0.0);
vector<double> vzz, vnzz; vector<double> vzz, vnzz;
string word, energyStr; string word, energyStr;
@ -104,6 +110,31 @@ TTrimSPData::TTrimSPData(const string &path, vector<string> &energyVec) {
for(unsigned int i(0); i<fEnergy.size();i++) for(unsigned int i(0); i<fEnergy.size();i++)
fIsNormalized.push_back(false); fIsNormalized.push_back(false);
fEnergyIter = fEnergy.end();
}
void TTrimSPData::UseHighResolution(double e) {
fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e);
if(fEnergyIter != fEnergy.end()) {
unsigned int i(fEnergyIter - fEnergy.begin());
vector<double> vecZ;
vector<double> vecNZ;
for(double zz(1.); zz<2100.; zz+=1.) {
vecZ.push_back(zz);
vecNZ.push_back(GetNofZ(zz/10.0, e));
}
fDataZ[i] = vecZ;
fDataNZ[i] = vecNZ;
fOrigDataNZ[i] = vecNZ;
fDZ[i] = 1.;
fIsNormalized[i] = false;
return;
}
cout << "TTrimSPData::DataZ: No implantation profile available for the specified energy... Nothing happens." << endl;
return;
} }
//--------------------- //---------------------
@ -112,16 +143,15 @@ TTrimSPData::TTrimSPData(const string &path, vector<string> &energyVec) {
vector<double> TTrimSPData::DataZ(double e) const { vector<double> TTrimSPData::DataZ(double e) const {
for(unsigned int i(0); i<fEnergy.size(); i++) { fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e);
// cout << tEnergy[i] << " " << e << " " << tEnergy[i] - e << endl;
if(!(fEnergy[i] - e)) { if(fEnergyIter != fEnergy.end()) {
return fDataZ[i]; unsigned int i(fEnergyIter - fEnergy.begin());
} return fDataZ[i];
} }
// default // default
cout << "TTrimSPData::DataZ: No implantation profile available for the specified energy... You get back the first one." << endl; cout << "TTrimSPData::DataZ: No implantation profile available for the specified energy... You get back the first one." << endl;
return fDataZ[0]; return fDataZ[0];
} }
//--------------------- //---------------------
@ -131,10 +161,11 @@ vector<double> TTrimSPData::DataZ(double e) const {
vector<double> TTrimSPData::DataNZ(double e) const { vector<double> TTrimSPData::DataNZ(double e) const {
for(unsigned int i(0); i<fEnergy.size(); i++) { fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e);
if(!(fEnergy[i] - e)) {
return fDataNZ[i]; if(fEnergyIter != fEnergy.end()) {
} unsigned int i(fEnergyIter - fEnergy.begin());
return fDataNZ[i];
} }
// default // default
cout << "TTrimSPData::DataNZ: No implantation profile available for the specified energy... You get back the first one." << endl; cout << "TTrimSPData::DataNZ: No implantation profile available for the specified energy... You get back the first one." << endl;
@ -148,10 +179,11 @@ vector<double> TTrimSPData::DataNZ(double e) const {
vector<double> TTrimSPData::OrigDataNZ(double e) const { vector<double> TTrimSPData::OrigDataNZ(double e) const {
for(unsigned int i(0); i<fEnergy.size(); i++) { fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e);
if(!(fEnergy[i] - e)) {
return fOrigDataNZ[i]; if(fEnergyIter != fEnergy.end()) {
} unsigned int i(fEnergyIter - fEnergy.begin());
return fOrigDataNZ[i];
} }
// default // default
cout << "TTrimSPData::OrigDataNZ: No implantation profile available for the specified energy... You get back the first one." << endl; cout << "TTrimSPData::OrigDataNZ: No implantation profile available for the specified energy... You get back the first one." << endl;
@ -171,32 +203,33 @@ double TTrimSPData::LayerFraction(double e, unsigned int layno, const vector<dou
return 0.0; return 0.0;
} }
for(unsigned int i(0); i<fEnergy.size(); i++) { fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e);
if(!(fEnergy[i] - e)) {
// Because we do not know if the implantation profile is normalized or not, do not care about this and calculate the fraction from the beginning if(fEnergyIter != fEnergy.end()) {
// Total "number of muons" unsigned int i(fEnergyIter - fEnergy.begin());
double totalNumber(0.0); // Because we do not know if the implantation profile is normalized or not, do not care about this and calculate the fraction from the beginning
// Total "number of muons"
double totalNumber(0.0);
for(unsigned int j(0); j<fDataZ[i].size(); j++)
totalNumber += fDataNZ[i][j];
// "number of muons" in layer layno
double layerNumber(0.0);
if(!(layno-1)){
for(unsigned int j(0); j<fDataZ[i].size(); j++) for(unsigned int j(0); j<fDataZ[i].size(); j++)
totalNumber += fDataNZ[i][j]; if(fDataZ[i][j] < interface[0]*10.0)
// "number of muons" in layer layno layerNumber += fDataNZ[i][j];
double layerNumber(0.0); } else if(!(layno-interface.size()-1)){
if(!(layno-1)){ for(unsigned int j(0); j<fDataZ[i].size(); j++)
for(unsigned int j(0); j<fDataZ[i].size(); j++) if(fDataZ[i][j] >= *(interface.end()-1)*10.0)
if(fDataZ[i][j] < interface[0]*10.0) layerNumber += fDataNZ[i][j];
layerNumber += fDataNZ[i][j]; } else {
} else if(!(layno-interface.size()-1)){ for(unsigned int j(0); j<fDataZ[i].size(); j++)
for(unsigned int j(0); j<fDataZ[i].size(); j++) if(fDataZ[i][j] >= interface[layno-2]*10.0 && fDataZ[i][j] < interface[layno-1]*10.0)
if(fDataZ[i][j] >= *(interface.end()-1)*10.0) layerNumber += fDataNZ[i][j];
layerNumber += fDataNZ[i][j];
} else {
for(unsigned int j(0); j<fDataZ[i].size(); j++)
if(fDataZ[i][j] >= interface[layno-2]*10.0 && fDataZ[i][j] < interface[layno-1]*10.0)
layerNumber += fDataNZ[i][j];
}
// fraction of muons in layer layno
// cout << "Fraction of muons in layer " << layno << ": " << layerNumber/totalNumber << endl;
return layerNumber/totalNumber;
} }
// fraction of muons in layer layno
// cout << "Fraction of muons in layer " << layno << ": " << layerNumber/totalNumber << endl;
return layerNumber/totalNumber;
} }
// default // default
@ -241,38 +274,39 @@ void TTrimSPData::WeightLayers(double e, const vector<double>& interface, const
} }
} }
fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e);
// If all weights are equal to one, use the original n(z) vector
for(unsigned int i(0); i<weight.size(); i++) { for(unsigned int i(0); i<weight.size(); i++) {
if(weight[i]-1.0) if(weight[i]-1.0)
break; break;
if(i == weight.size() - 1) { if(i == weight.size() - 1) {
for(unsigned int j(0); j<fEnergy.size(); j++) { if(fEnergyIter != fEnergy.end()) {
if(!(fEnergy[j] - e)) { unsigned int j(fEnergyIter - fEnergy.begin());
fDataNZ[j] = fOrigDataNZ[j]; fDataNZ[j] = fOrigDataNZ[j];
fIsNormalized[j] = false; fIsNormalized[j] = false;
return; return;
}
} }
} }
} }
for(unsigned int i(0); i<fEnergy.size(); i++) { if(fEnergyIter != fEnergy.end()) {
if(!(fEnergy[i] - e)) { unsigned int i(fEnergyIter - fEnergy.begin());
unsigned int k(0); unsigned int k(0);
for(unsigned int j(0); j<fDataZ[i].size(); j++) { for(unsigned int j(0); j<fDataZ[i].size(); j++) {
if(k<interface.size()) { if(k<interface.size()) {
if(fDataZ[i][j] < interface[k]*10.0) 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]; fDataNZ[i][j] = fOrigDataNZ[i][j]*weight[k];
else {
k++;
fDataNZ[i][j] = fOrigDataNZ[i][j]*weight[k];
}
} }
fIsNormalized[i] = false; else
return; fDataNZ[i][j] = fOrigDataNZ[i][j]*weight[k];
} }
fIsNormalized[i] = false;
return;
} }
cout << "TTrimSPData::WeightLayers: No implantation profile available for the specified energy... No weighting done." << endl; cout << "TTrimSPData::WeightLayers: No implantation profile available for the specified energy... No weighting done." << endl;
@ -287,16 +321,15 @@ double TTrimSPData::GetNofZ(double zz, double e) const {
vector<double> z, nz; vector<double> z, nz;
for(unsigned int i(0); i<fEnergy.size(); i++) { fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e);
if(!(fEnergy[i] - e)) {
z = fDataZ[i]; if(fEnergyIter != fEnergy.end()) {
nz = fDataNZ[i]; unsigned int i(fEnergyIter - fEnergy.begin());
break; z = fDataZ[i];
} nz = fDataNZ[i];
if(i == fEnergy.size() - 1) { } else {
cout << "TTrimSPData::GetNofZ: No implantation profile available for the specified energy... Quitting!" << endl; cout << "TTrimSPData::GetNofZ: No implantation profile available for the specified energy... Quitting!" << endl;
exit(-1); exit(-1);
}
} }
if(zz < 0) if(zz < 0)
@ -326,18 +359,19 @@ double TTrimSPData::GetNofZ(double zz, double e) const {
void TTrimSPData::Normalize(double e) const { void TTrimSPData::Normalize(double e) const {
for(unsigned int i(0); i<fEnergy.size(); i++) { fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e);
if(!(fEnergy[i] - e)) {
double nZsum = 0.0;
for (unsigned int j(0); j<fDataZ[i].size(); j++)
nZsum += fDataNZ[i][j];
nZsum *= fDZ[i];
for (unsigned int j(0); j<fDataZ[i].size(); j++)
fDataNZ[i][j] /= nZsum;
fIsNormalized[i] = true; if(fEnergyIter != fEnergy.end()) {
return; unsigned int i(fEnergyIter - fEnergy.begin());
} double nZsum = 0.0;
for (unsigned int j(0); j<fDataZ[i].size(); j++)
nZsum += fDataNZ[i][j];
nZsum *= fDZ[i];
for (unsigned int j(0); j<fDataZ[i].size(); j++)
fDataNZ[i][j] /= nZsum;
fIsNormalized[i] = true;
return;
} }
// default // default
cout << "TTrimSPData::Normalize: No implantation profile available for the specified energy... No normalization done." << endl; cout << "TTrimSPData::Normalize: No implantation profile available for the specified energy... No normalization done." << endl;
@ -350,10 +384,11 @@ void TTrimSPData::Normalize(double e) const {
//--------------------- //---------------------
bool TTrimSPData::IsNormalized(double e) const { bool TTrimSPData::IsNormalized(double e) const {
for(unsigned int i(0); i<fEnergy.size(); i++) { fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e);
if(!(fEnergy[i] - e)) {
return fIsNormalized[i]; if(fEnergyIter != fEnergy.end()) {
} unsigned int i(fEnergyIter - fEnergy.begin());
return fIsNormalized[i];
} }
cout << "TTrimSPData::IsNormalized: No implantation profile available for the specified energy... Returning false! Check your code!" << endl; cout << "TTrimSPData::IsNormalized: No implantation profile available for the specified energy... Returning false! Check your code!" << endl;
@ -365,23 +400,50 @@ bool TTrimSPData::IsNormalized(double e) const {
//--------------------- //---------------------
double TTrimSPData::MeanRange(double e) const { double TTrimSPData::MeanRange(double e) const {
for(unsigned int i(0); i<fEnergy.size(); i++) { fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e);
if(!(fEnergy[i] - e)) {
if (!fIsNormalized[i]) if(fEnergyIter != fEnergy.end()) {
Normalize(e); unsigned int i(fEnergyIter - fEnergy.begin());
double mean(0.0); if (!fIsNormalized[i])
for(unsigned int j(0); j<fDataNZ[i].size(); j++){ Normalize(e);
mean += fDataNZ[i][j]*fDataZ[i][j]; double mean(0.0);
} for(unsigned int j(0); j<fDataNZ[i].size(); j++){
mean *= fDZ[i]/10.0; mean += fDataNZ[i][j]*fDataZ[i][j];
return mean;
} }
mean *= fDZ[i]/10.0;
return mean;
} }
cout << "TTrimSPData::MeanRange: No implantation profile available for the specified energy... Returning -1! Check your code!" << endl; cout << "TTrimSPData::MeanRange: No implantation profile available for the specified energy... Returning -1! Check your code!" << endl;
return -1.; return -1.;
} }
//---------------------
// Find the peak range in (nm) for a given energy e
//---------------------
double TTrimSPData::PeakRange(double e) const {
fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e);
if(fEnergyIter != fEnergy.end()) {
unsigned int i(fEnergyIter - fEnergy.begin());
vector<double>::const_iterator nziter;
nziter = max_element(fDataNZ[i].begin(),fDataNZ[i].end());
if(nziter != fDataNZ[i].end()){
unsigned int j(nziter - fDataNZ[i].begin());
return fDataZ[i][j];
}
cout << "TTrimSPData::PeakRange: No maximum found in the implantation profile... Returning -1! Please check the profile!" << endl;
return -1.;
}
cout << "TTrimSPData::PeakRange: No implantation profile available for the specified energy... Returning -1! Check your code!" << endl;
return -1.;
}
//--------------------- //---------------------
// Convolve the n(z)-vector calculated by trim.SP for a given energy e [keV] with a gaussian exp(-z^2/(2*w^2)) // Convolve the n(z)-vector calculated by trim.SP for a given energy e [keV] with a gaussian exp(-z^2/(2*w^2))
// No normalization is done! // No normalization is done!
@ -394,27 +456,28 @@ void TTrimSPData::ConvolveGss(double w, double e) const {
vector<double> z, nz, gss; vector<double> z, nz, gss;
double nn; double nn;
for(unsigned int i(0); i<fEnergy.size(); i++) { fEnergyIter = find(fEnergy.begin(), fEnergy.end(), e);
if(!(fEnergy[i] - e)) {
z = fDataZ[i];
nz = fOrigDataNZ[i];
for(unsigned int k(0); k<z.size(); k++) { if(fEnergyIter != fEnergy.end()) {
gss.push_back(exp(-z[k]*z[k]/200.0/w/w)); unsigned int i(fEnergyIter - fEnergy.begin());
} z = fDataZ[i];
nz = fOrigDataNZ[i];
for(unsigned int k(0); k<nz.size(); k++) { for(unsigned int k(0); k<z.size(); k++) {
nn = 0.0; gss.push_back(exp(-z[k]*z[k]/200.0/w/w));
for(unsigned int j(0); j<nz.size(); j++) {
nn += nz[j]*gss[abs(int(k)-int(j))];
}
fDataNZ[i][k] = nn;
}
fIsNormalized[i] = false;
return;
} }
for(unsigned int k(0); k<nz.size(); k++) {
nn = 0.0;
for(unsigned int j(0); j<nz.size(); j++) {
nn += nz[j]*gss[abs(int(k)-int(j))];
}
fDataNZ[i][k] = nn;
}
fIsNormalized[i] = false;
return;
} }
cout << "TTrimSPData::ConvolveGss: No implantation profile available for the specified energy... No convolution done!" << endl; cout << "TTrimSPData::ConvolveGss: No implantation profile available for the specified energy... No convolution done!" << endl;

5
src/external/TFitPofB-lib/include/TTrimSPDataHandler.h vendored
View File

@ -5,7 +5,7 @@
Author: Bastian M. Wojek Author: Bastian M. Wojek
e-mail: bastian.wojek@psi.ch e-mail: bastian.wojek@psi.ch
2008/11/21 2009/05/15
***************************************************************************/ ***************************************************************************/
@ -35,6 +35,7 @@ public:
vector<double> DataZ(double) const; vector<double> DataZ(double) const;
vector<double> DataNZ(double) const; vector<double> DataNZ(double) const;
vector<double> OrigDataNZ(double) const; vector<double> OrigDataNZ(double) const;
void UseHighResolution(double);
void WeightLayers(double, const vector<double>&, const vector<double>&) const; void WeightLayers(double, const vector<double>&, const vector<double>&) const;
double LayerFraction(double, unsigned int, const vector<double>&) const; double LayerFraction(double, unsigned int, const vector<double>&) const;
double GetNofZ(double, double) const; double GetNofZ(double, double) const;
@ -42,6 +43,7 @@ public:
bool IsNormalized(double) const; bool IsNormalized(double) const;
void ConvolveGss(double, double) const; void ConvolveGss(double, double) const;
double MeanRange(double) const; double MeanRange(double) const;
double PeakRange(double) const;
private: private:
vector<double> fEnergy; vector<double> fEnergy;
@ -50,6 +52,7 @@ private:
mutable vector< vector<double> > fDataNZ; mutable vector< vector<double> > fDataNZ;
vector< vector<double> > fOrigDataNZ; vector< vector<double> > fOrigDataNZ;
mutable vector<bool> fIsNormalized; mutable vector<bool> fIsNormalized;
mutable vector<double>::const_iterator fEnergyIter;
}; };
#endif // _TTrimSPDataHandler_H_ #endif // _TTrimSPDataHandler_H_