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

Harmonized a few function calls in libFitPofB

Browse Source
This commit is contained in:
Bastian M. Wojek 2010-08-04 13:12:14 +00:00
parent 7536ca6fc6
commit 599c83a17c
4 changed files with 177 additions and 324 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

463
src/external/TFitPofB-lib/classes/TLondon1D.cpp vendored
View File

@ -46,7 +46,6 @@ ClassImp(TProximity1D1LHSGss)
ClassImp(TLondon1D3L) ClassImp(TLondon1D3L)
ClassImp(TLondon1D3LS) ClassImp(TLondon1D3LS)
// ClassImp(TLondon1D4L) // ClassImp(TLondon1D4L)
ClassImp(TLondon1D3LSub)
ClassImp(TLondon1D3Lestimate) ClassImp(TLondon1D3Lestimate)
@ -158,19 +157,6 @@ TLondon1D3LS::~TLondon1D3LS() {
// fPofT = 0; // fPofT = 0;
// } // }
TLondon1D3LSub::~TLondon1D3LSub() {
fPar.clear();
fParForBofZ.clear();
fParForPofB.clear();
fParForPofT.clear();
delete fImpProfile;
fImpProfile = 0;
delete fPofB;
fPofB = 0;
delete fPofT;
fPofT = 0;
}
//------------------ //------------------
// Constructor of the TLondon1DHS class -- reading available implantation profiles and // Constructor of the TLondon1DHS class -- reading available implantation profiles and
// creates (a pointer to) the TPofTCalc object (with the FFT plan) // creates (a pointer to) the TPofTCalc object (with the FFT plan)
@ -234,13 +220,12 @@ TLondon1DHS::TLondon1DHS() : fCalcNeeded(true), fFirstCall(true) {
double TLondon1DHS::operator()(double t, const vector<double> &par) const { double TLondon1DHS::operator()(double t, const vector<double> &par) const {
assert(par.size() == 5 || par.size() == 6); assert(par.size() == 5);
if(t<0.0) if(t<0.0)
return cos(par[0]*0.017453293); return cos(par[0]*0.017453293);
bool dead_layer_changed(false); bool dead_layer_changed(false);
bool width_changed(false);
// check if the function is called the first time and if yes, read in parameters // check if the function is called the first time and if yes, read in parameters
@ -252,8 +237,6 @@ double TLondon1DHS::operator()(double t, const vector<double> &par) const {
} }
fFirstCall = false; fFirstCall = false;
dead_layer_changed = true; dead_layer_changed = true;
if(par.size() == 6)
width_changed = true;
} }
// check if any parameter has changed // check if any parameter has changed
@ -271,8 +254,6 @@ double TLondon1DHS::operator()(double t, const vector<double> &par) const {
only_phase_changed = false; only_phase_changed = false;
if (i == 3) { if (i == 3) {
dead_layer_changed = true; dead_layer_changed = true;
} else if (i == 5) {
width_changed = true;
} }
} }
} }
@ -294,10 +275,6 @@ double TLondon1DHS::operator()(double t, const vector<double> &par) const {
for (unsigned int i(2); i<fPar.size(); i++) for (unsigned int i(2); i<fPar.size(); i++)
fParForBofZ[i-2] = par[i]; fParForBofZ[i-2] = par[i];
if(width_changed) { // Convolution of the implantation profile with Gaussian
fImpProfile->ConvolveGss(par[5], par[1]);
}
fParForPofB[2] = par[1]; // energy fParForPofB[2] = par[1]; // energy
fParForPofB[3] = par[2]; // Bkg-Field fParForPofB[3] = par[2]; // Bkg-Field
//fParForPofB[4] = 0.005; // Bkg-width (in principle zero) //fParForPofB[4] = 0.005; // Bkg-width (in principle zero)
@ -361,7 +338,10 @@ TLondon1D1L::TLondon1D1L() : fCalcNeeded(true), fFirstCall(true), fCallCounter(0
fParForPofB.push_back(startupHandler->GetDeltat()); fParForPofB.push_back(startupHandler->GetDeltat());
fParForPofB.push_back(startupHandler->GetDeltaB()); fParForPofB.push_back(startupHandler->GetDeltaB());
fParForPofB.push_back(0.0); fParForPofB.push_back(0.0); // Energy
fParForPofB.push_back(0.0); // Bkg-Field
fParForPofB.push_back(0.005); // Bkg-width
fParForPofB.push_back(0.0); // Bkg-weight
fImpProfile = new TTrimSPData(rge_path, energy_vec); fImpProfile = new TTrimSPData(rge_path, energy_vec);
@ -388,30 +368,31 @@ TLondon1D1L::TLondon1D1L() : fCalcNeeded(true), fFirstCall(true), fCallCounter(0
double TLondon1D1L::operator()(double t, const vector<double> &par) const { double TLondon1D1L::operator()(double t, const vector<double> &par) const {
assert(par.size() == 6); assert(par.size() == 6 || par.size() == 8);
// Debugging // Debugging
// Count the number of function calls // Count the number of function calls
// fCallCounter++; // fCallCounter++;
if(t<0.0) if(t<0.0)
return cos(par[0]*0.017453293); return cos(par[0]*0.017453293);
bool bkg_fraction_changed(false);
bool weights_changed(false);
// check if the function is called the first time and if yes, read in parameters // check if the function is called the first time and if yes, read in parameters
if(fFirstCall){ if(fFirstCall){
fPar = par; fPar = par;
/* for (unsigned int i(0); i<fPar.size(); i++){
cout << "fPar[" << i << "] = " << fPar[i] << endl;
}
*/
for (unsigned int i(2); i<fPar.size(); i++){ for (unsigned int i(2); i<fPar.size(); i++){
fParForBofZ.push_back(fPar[i]); fParForBofZ.push_back(fPar[i]);
// cout << "fParForBofZ[" << i-2 << "] = " << fParForBofZ[i-2] << endl;
} }
fFirstCall=false;
// cout << this << endl; fFirstCall = false;
bkg_fraction_changed = true;
if (par.size() == 8)
weights_changed = true;
} }
// check if any parameter has changed // check if any parameter has changed
@ -427,6 +408,10 @@ double TLondon1D1L::operator()(double t, const vector<double> &par) const {
only_phase_changed = true; only_phase_changed = true;
} else { } else {
only_phase_changed = false; only_phase_changed = false;
if (i == 3 || i == 4)
bkg_fraction_changed = true;
if (i == 6 || i == 7)
weights_changed = true;
} }
} }
} }
@ -448,6 +433,33 @@ double TLondon1D1L::operator()(double t, const vector<double> &par) const {
fParForBofZ[i-2] = par[i]; fParForBofZ[i-2] = par[i];
fParForPofB[2] = par[1]; // energy fParForPofB[2] = par[1]; // energy
fParForPofB[3] = par[2]; // Bkg-Field
//fParForPofB[4] = 0.005; // Bkg-width (in principle zero)
if(weights_changed) {
vector<double> interfaces;
interfaces.push_back(par[3]+par[4]);
vector<double> weights;
for(unsigned int i(6); i<8; i++)
weights.push_back(par[i]);
// cout << "Weighting has changed, re-calculating n(z) now..." << endl;
fImpProfile->WeightLayers(par[1], interfaces, weights);
interfaces.clear();
weights.clear();
}
if(bkg_fraction_changed || weights_changed) {
vector<double> interfaces;
interfaces.push_back(par[3]);// dead layer
interfaces.push_back(par[3] + par[4]);// dead layer + first layer
fParForPofB[5] = fImpProfile->LayerFraction(par[1], 1, interfaces) + // Fraction of muons in the deadlayer
fImpProfile->LayerFraction(par[1], 3, interfaces); // Fraction of muons in the substrate
interfaces.clear();
}
TLondon1D_1L BofZ1(fParForBofZ); TLondon1D_1L BofZ1(fParForBofZ);
fPofB->UnsetPBExists(); fPofB->UnsetPBExists();
@ -482,7 +494,7 @@ double TLondon1D1L::operator()(double t, const vector<double> &par) const {
// creates (a pointer to) the TPofTCalc object (with the FFT plan) // creates (a pointer to) the TPofTCalc object (with the FFT plan)
//------------------ //------------------
TLondon1D2L::TLondon1D2L() : fCalcNeeded(true), fFirstCall(true), fLastTwoChanged(true) { TLondon1D2L::TLondon1D2L() : fCalcNeeded(true), fFirstCall(true) {
// read startup file // read startup file
string startup_path_name("TFitPofB_startup.xml"); string startup_path_name("TFitPofB_startup.xml");
@ -509,7 +521,10 @@ TLondon1D2L::TLondon1D2L() : fCalcNeeded(true), fFirstCall(true), fLastTwoChange
fParForPofB.push_back(startupHandler->GetDeltat()); fParForPofB.push_back(startupHandler->GetDeltat());
fParForPofB.push_back(startupHandler->GetDeltaB()); fParForPofB.push_back(startupHandler->GetDeltaB());
fParForPofB.push_back(0.0); fParForPofB.push_back(0.0); // Energy
fParForPofB.push_back(0.0); // Bkg-Field
fParForPofB.push_back(0.005); // Bkg-width
fParForPofB.push_back(0.0); // Bkg-weight
fImpProfile = new TTrimSPData(rge_path, energy_vec); fImpProfile = new TTrimSPData(rge_path, energy_vec);
@ -536,25 +551,26 @@ TLondon1D2L::TLondon1D2L() : fCalcNeeded(true), fFirstCall(true), fLastTwoChange
double TLondon1D2L::operator()(double t, const vector<double> &par) const { double TLondon1D2L::operator()(double t, const vector<double> &par) const {
assert(par.size() == 10); assert(par.size() == 8 || par.size() == 11);
if(t<0.0) if(t<0.0)
return cos(par[0]*0.017453293); return cos(par[0]*0.017453293);
bool bkg_fraction_changed(false);
bool weights_changed(false);
// check if the function is called the first time and if yes, read in parameters // check if the function is called the first time and if yes, read in parameters
if(fFirstCall){ if(fFirstCall){
fPar = par; fPar = par;
/* for (unsigned int i(0); i<fPar.size(); i++){
cout << "fPar[" << i << "] = " << fPar[i] << endl;
}
*/
for (unsigned int i(2); i<fPar.size(); i++){ for (unsigned int i(2); i<fPar.size(); i++){
fParForBofZ.push_back(fPar[i]); fParForBofZ.push_back(fPar[i]);
// cout << "fParForBofZ[" << i-2 << "] = " << fParForBofZ[i-2] << endl;
} }
fFirstCall=false; fFirstCall=false;
bkg_fraction_changed = true;
if (par.size() == 11)
weights_changed = true;
} }
// check if any parameter has changed // check if any parameter has changed
@ -570,9 +586,11 @@ double TLondon1D2L::operator()(double t, const vector<double> &par) const {
only_phase_changed = true; only_phase_changed = true;
} else { } else {
only_phase_changed = false; only_phase_changed = false;
if (i == 3 || i == 4 || i == 5)
bkg_fraction_changed = true;
if (i == 8 || i == 9 || i == 10)
weights_changed = true;
} }
if (i == fPar.size()-2 || i == fPar.size()-1)
fLastTwoChanged = true;
} }
} }
@ -593,17 +611,32 @@ double TLondon1D2L::operator()(double t, const vector<double> &par) const {
fParForBofZ[i-2] = par[i]; fParForBofZ[i-2] = par[i];
fParForPofB[2] = par[1]; // energy fParForPofB[2] = par[1]; // energy
fParForPofB[3] = par[2]; // Bkg-Field
//fParForPofB[4] = 0.005; // Bkg-width (in principle zero)
if(fLastTwoChanged) { if(weights_changed) {
vector<double> interfaces; vector<double> interfaces;
interfaces.push_back(par[3]+par[4]); interfaces.push_back(par[3]+par[4]);
interfaces.push_back(par[3]+par[4]+par[5]);
vector<double> weights; vector<double> weights;
for(unsigned int i(par.size()-2); i<par.size(); i++) for(unsigned int i(8); i<11; i++)
weights.push_back(par[i]); weights.push_back(par[i]);
// cout << "Weighting has changed, re-calculating n(z) now..." << endl; // cout << "Weighting has changed, re-calculating n(z) now..." << endl;
fImpProfile->WeightLayers(par[1], interfaces, weights); fImpProfile->WeightLayers(par[1], interfaces, weights);
interfaces.clear();
weights.clear();
}
if(bkg_fraction_changed || weights_changed) {
vector<double> interfaces;
interfaces.push_back(par[3]);// dead layer
interfaces.push_back(par[3] + par[4] + par[5]);// dead layer + first layer + second layer
fParForPofB[5] = fImpProfile->LayerFraction(par[1], 1, interfaces) + // Fraction of muons in the deadlayer
fImpProfile->LayerFraction(par[1], 3, interfaces); // Fraction of muons in the substrate
interfaces.clear();
} }
TLondon1D_2L BofZ2(fParForBofZ); TLondon1D_2L BofZ2(fParForBofZ);
@ -619,7 +652,6 @@ double TLondon1D2L::operator()(double t, const vector<double> &par) const {
fPofT->CalcPol(fParForPofT); fPofT->CalcPol(fParForPofT);
fCalcNeeded = false; fCalcNeeded = false;
fLastTwoChanged = false;
} }
return fPofT->Eval(t); return fPofT->Eval(t);
@ -659,7 +691,7 @@ TProximity1D1LHS::TProximity1D1LHS() : fCalcNeeded(true), fFirstCall(true) {
fParForPofB.push_back(startupHandler->GetDeltat()); fParForPofB.push_back(startupHandler->GetDeltat());
fParForPofB.push_back(startupHandler->GetDeltaB()); fParForPofB.push_back(startupHandler->GetDeltaB());
fParForPofB.push_back(0.0); fParForPofB.push_back(0.0); // Energy
fParForPofB.push_back(0.0); // Bkg-Field fParForPofB.push_back(0.0); // Bkg-Field
fParForPofB.push_back(0.01); // Bkg-width fParForPofB.push_back(0.01); // Bkg-width
fParForPofB.push_back(0.0); // Bkg-weight fParForPofB.push_back(0.0); // Bkg-weight
@ -689,31 +721,23 @@ TProximity1D1LHS::TProximity1D1LHS() : fCalcNeeded(true), fFirstCall(true) {
double TProximity1D1LHS::operator()(double t, const vector<double> &par) const { double TProximity1D1LHS::operator()(double t, const vector<double> &par) const {
assert(par.size() == 8); assert(par.size() == 7);
if(t<0.0) if(t<0.0)
return cos(par[0]*0.017453293); return cos(par[0]*0.017453293);
// check if the function is called the first time and if yes, read in parameters // check if the function is called the first time and if yes, read in parameters
bool width_changed(false);
bool dead_layer_changed(false); bool dead_layer_changed(false);
if(fFirstCall){ if(fFirstCall){
fPar = par; fPar = par;
// for (unsigned int i(0); i<fPar.size(); i++){
// cout << "fPar[" << i << "] = " << fPar[i] << endl;
// }
for (unsigned int i(2); i<fPar.size(); i++){ for (unsigned int i(2); i<fPar.size(); i++){
fParForBofZ.push_back(fPar[i]); fParForBofZ.push_back(fPar[i]);
// cout << "fParForBofZ[" << i-2 << "] = " << fParForBofZ[i-2] << endl;
} }
fFirstCall = false; fFirstCall = false;
width_changed = true;
dead_layer_changed = true; dead_layer_changed = true;
// cout << this << endl;
} }
// check if any parameter has changed // check if any parameter has changed
@ -729,9 +753,6 @@ double TProximity1D1LHS::operator()(double t, const vector<double> &par) const {
only_phase_changed = true; only_phase_changed = true;
} else { } else {
only_phase_changed = false; only_phase_changed = false;
if (i == 7){
width_changed = true;
}
if (i == 4){ if (i == 4){
dead_layer_changed = true; dead_layer_changed = true;
} }
@ -755,10 +776,6 @@ double TProximity1D1LHS::operator()(double t, const vector<double> &par) const {
for (unsigned int i(2); i<fPar.size(); i++) for (unsigned int i(2); i<fPar.size(); i++)
fParForBofZ[i-2] = par[i]; fParForBofZ[i-2] = par[i];
if(width_changed) { // Convolution of the implantation profile with Gaussian
fImpProfile->ConvolveGss(par[7], par[1]);
}
fParForPofB[2] = par[1]; // energy fParForPofB[2] = par[1]; // energy
fParForPofB[3] = par[2]; // Bkg-Field fParForPofB[3] = par[2]; // Bkg-Field
//fParForPofB[4] = 0.005; // Bkg-width (in principle zero) //fParForPofB[4] = 0.005; // Bkg-width (in principle zero)
@ -823,7 +840,9 @@ TProximity1D1LHSGss::TProximity1D1LHSGss() : fCalcNeeded(true), fFirstCall(true)
fParForPofB.push_back(startupHandler->GetDeltat()); fParForPofB.push_back(startupHandler->GetDeltat());
fParForPofB.push_back(startupHandler->GetDeltaB()); fParForPofB.push_back(startupHandler->GetDeltaB());
fParForPofB.push_back(0.0); fParForPofB.push_back(0.0);
// fParForPofB.push_back(0.0); fParForPofB.push_back(0.0); // Bkg-Field
fParForPofB.push_back(0.01); // Bkg-width
fParForPofB.push_back(0.0); // Bkg-weight
fImpProfile = new TTrimSPData(rge_path, energy_vec); fImpProfile = new TTrimSPData(rge_path, energy_vec);
@ -845,7 +864,7 @@ TProximity1D1LHSGss::TProximity1D1LHSGss() : fCalcNeeded(true), fFirstCall(true)
//------------------ //------------------
// TProximity1D1LHS-Method that calls the procedures to create B(z), p(B) and P(t) // TProximity1D1LHS-Method that calls the procedures to create B(z), p(B) and P(t)
// It finally returns P(t) for a given t. // It finally returns P(t) for a given t.
// Parameters: all the parameters for the function to be fitted through TProximity1D1LHS // Parameters: all the parameters for the function to be fitted through TProximity1D1LHSGss
//------------------ //------------------
double TProximity1D1LHSGss::operator()(double t, const vector<double> &par) const { double TProximity1D1LHSGss::operator()(double t, const vector<double> &par) const {
@ -857,19 +876,16 @@ double TProximity1D1LHSGss::operator()(double t, const vector<double> &par) cons
// check if the function is called the first time and if yes, read in parameters // check if the function is called the first time and if yes, read in parameters
bool dead_layer_changed(false);
if(fFirstCall){ if(fFirstCall){
fPar = par; fPar = par;
// for (unsigned int i(0); i<fPar.size(); i++){
// cout << "fPar[" << i << "] = " << fPar[i] << endl;
// }
for (unsigned int i(2); i<fPar.size(); i++){ for (unsigned int i(2); i<fPar.size(); i++){
fParForBofZ.push_back(fPar[i]); fParForBofZ.push_back(fPar[i]);
// cout << "fParForBofZ[" << i-2 << "] = " << fParForBofZ[i-2] << endl;
} }
fFirstCall=false; fFirstCall = false;
// cout << this << endl; dead_layer_changed = true;
} }
// check if any parameter has changed // check if any parameter has changed
@ -885,6 +901,9 @@ double TProximity1D1LHSGss::operator()(double t, const vector<double> &par) cons
only_phase_changed = true; only_phase_changed = true;
} else { } else {
only_phase_changed = false; only_phase_changed = false;
if (i == 4){
dead_layer_changed = true;
}
} }
} }
} }
@ -906,12 +925,22 @@ double TProximity1D1LHSGss::operator()(double t, const vector<double> &par) cons
fParForBofZ[i-2] = par[i]; fParForBofZ[i-2] = par[i];
fParForPofB[2] = par[1]; // energy fParForPofB[2] = par[1]; // energy
fParForPofB[3] = par[2]; // Bkg-Field
//fParForPofB[4] = 0.005; // Bkg-width (in principle zero)
if(dead_layer_changed){
vector<double> interfaces;
interfaces.push_back(par[4]);// dead layer
fParForPofB[5] = fImpProfile->LayerFraction(par[1], 1, interfaces); // Fraction of muons in the deadlayer
interfaces.clear();
}
TProximity1D_1LHSGss BofZ(fParForBofZ); TProximity1D_1LHSGss BofZ(fParForBofZ);
fPofB->UnsetPBExists(); fPofB->UnsetPBExists();
fPofB->Calculate(&BofZ, fImpProfile, fParForPofB); fPofB->Calculate(&BofZ, fImpProfile, fParForPofB);
fPofT->DoFFT(); fPofT->DoFFT();
}/* else { }/* else {
cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl; cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl;
}*/ }*/
@ -930,7 +959,7 @@ double TProximity1D1LHSGss::operator()(double t, const vector<double> &par) cons
// creates (a pointer to) the TPofTCalc object (with the FFT plan) // creates (a pointer to) the TPofTCalc object (with the FFT plan)
//------------------ //------------------
TLondon1D3L::TLondon1D3L() : fCalcNeeded(true), fFirstCall(true), fLastThreeChanged(true) { TLondon1D3L::TLondon1D3L() : fCalcNeeded(true), fFirstCall(true) {
// read startup file // read startup file
string startup_path_name("TFitPofB_startup.xml"); string startup_path_name("TFitPofB_startup.xml");
@ -957,7 +986,10 @@ TLondon1D3L::TLondon1D3L() : fCalcNeeded(true), fFirstCall(true), fLastThreeChan
fParForPofB.push_back(startupHandler->GetDeltat()); fParForPofB.push_back(startupHandler->GetDeltat());
fParForPofB.push_back(startupHandler->GetDeltaB()); fParForPofB.push_back(startupHandler->GetDeltaB());
fParForPofB.push_back(0.0); fParForPofB.push_back(0.0); // Energy
fParForPofB.push_back(0.0); // Bkg-Field
fParForPofB.push_back(0.005); // Bkg-width
fParForPofB.push_back(0.0); // Bkg-weight
fImpProfile = new TTrimSPData(rge_path, energy_vec); fImpProfile = new TTrimSPData(rge_path, energy_vec);
@ -984,25 +1016,26 @@ TLondon1D3L::TLondon1D3L() : fCalcNeeded(true), fFirstCall(true), fLastThreeChan
double TLondon1D3L::operator()(double t, const vector<double> &par) const { double TLondon1D3L::operator()(double t, const vector<double> &par) const {
assert(par.size() == 13); assert(par.size() == 10 || par.size() == 14);
if(t<0.0) if(t<0.0)
return cos(par[0]*0.017453293); return cos(par[0]*0.017453293);
bool bkg_fraction_changed(false);
bool weights_changed(false);
// check if the function is called the first time and if yes, read in parameters // check if the function is called the first time and if yes, read in parameters
if(fFirstCall){ if(fFirstCall){
fPar = par; fPar = par;
/* for (unsigned int i(0); i<fPar.size(); i++){
cout << "fPar[" << i << "] = " << fPar[i] << endl;
}
*/
for (unsigned int i(2); i<fPar.size(); i++){ for (unsigned int i(2); i<fPar.size(); i++){
fParForBofZ.push_back(fPar[i]); fParForBofZ.push_back(fPar[i]);
// cout << "fParForBofZ[" << i-2 << "] = " << fParForBofZ[i-2] << endl;
} }
fFirstCall=false; fFirstCall=false;
bkg_fraction_changed = true;
if (par.size() == 14)
weights_changed = true;
} }
// check if any parameter has changed // check if any parameter has changed
@ -1018,9 +1051,11 @@ double TLondon1D3L::operator()(double t, const vector<double> &par) const {
only_phase_changed = true; only_phase_changed = true;
} else { } else {
only_phase_changed = false; only_phase_changed = false;
if (i == 3 || i == 4 || i == 5 || i == 6)
bkg_fraction_changed = true;
if (i == 10 || i == 11 || i == 12 || i == 13)
weights_changed = true;
} }
if (i == fPar.size()-3 || i == fPar.size()-2 || i == fPar.size()-1)
fLastThreeChanged = true;
} }
} }
@ -1041,32 +1076,33 @@ double TLondon1D3L::operator()(double t, const vector<double> &par) const {
fParForBofZ[i-2] = par[i]; fParForBofZ[i-2] = par[i];
fParForPofB[2] = par[1]; // energy fParForPofB[2] = par[1]; // energy
fParForPofB[3] = par[2]; // Bkg-Field
//fParForPofB[4] = 0.005; // Bkg-width (in principle zero)
/* DEBUG --------------------------- if(weights_changed) {
for(unsigned int i(0); i<fParForBofZ.size(); i++) {
cout << "ParForBofZ[" << i << "] = " << fParForBofZ[i] << endl;
}
for(unsigned int i(0); i<fParForPofB.size(); i++) {
cout << "ParForPofB[" << i << "] = " << fParForPofB[i] << endl;
}
for(unsigned int i(0); i<fParForPofT.size(); i++) {
cout << "ParForPofT[" << i << "] = " << fParForPofT[i] << endl;
}
------------------------------------*/
if(fLastThreeChanged) {
vector<double> interfaces; vector<double> interfaces;
interfaces.push_back(par[3]+par[4]); interfaces.push_back(par[3]+par[4]);
interfaces.push_back(par[3]+par[4]+par[5]); interfaces.push_back(par[3]+par[4]+par[5]);
interfaces.push_back(par[3]+par[4]+par[5]+par[6]);
vector<double> weights; vector<double> weights;
for(unsigned int i(par.size()-3); i<par.size(); i++) for(unsigned int i(10); i<14; i++)
weights.push_back(par[i]); weights.push_back(par[i]);
// cout << "Weighting has changed, re-calculating n(z) now..." << endl; // cout << "Weighting has changed, re-calculating n(z) now..." << endl;
fImpProfile->WeightLayers(par[1], interfaces, weights); fImpProfile->WeightLayers(par[1], interfaces, weights);
interfaces.clear();
weights.clear();
}
if(bkg_fraction_changed || weights_changed) {
vector<double> interfaces;
interfaces.push_back(par[3]);// dead layer
interfaces.push_back(par[3] + par[4] + par[5] + par[6]);// dead layer + first layer + second layer + third layer
fParForPofB[5] = fImpProfile->LayerFraction(par[1], 1, interfaces) + // Fraction of muons in the deadlayer
fImpProfile->LayerFraction(par[1], 3, interfaces); // Fraction of muons in the substrate
interfaces.clear();
} }
TLondon1D_3L BofZ3(fParForBofZ); TLondon1D_3L BofZ3(fParForBofZ);
@ -1074,7 +1110,6 @@ double TLondon1D3L::operator()(double t, const vector<double> &par) const {
fPofB->Calculate(&BofZ3, fImpProfile, fParForPofB); fPofB->Calculate(&BofZ3, fImpProfile, fParForPofB);
fPofT->DoFFT(); fPofT->DoFFT();
}/* else { }/* else {
cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl; cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl;
}*/ }*/
@ -1082,7 +1117,6 @@ double TLondon1D3L::operator()(double t, const vector<double> &par) const {
fPofT->CalcPol(fParForPofT); fPofT->CalcPol(fParForPofT);
fCalcNeeded = false; fCalcNeeded = false;
fLastThreeChanged = false;
} }
return fPofT->Eval(t); return fPofT->Eval(t);
@ -1094,7 +1128,7 @@ double TLondon1D3L::operator()(double t, const vector<double> &par) const {
// creates (a pointer to) the TPofTCalc object (with the FFT plan) // creates (a pointer to) the TPofTCalc object (with the FFT plan)
//------------------ //------------------
TLondon1D3LS::TLondon1D3LS() : fCalcNeeded(true), fFirstCall(true), fLastThreeChanged(true) { TLondon1D3LS::TLondon1D3LS() : fCalcNeeded(true), fFirstCall(true) {
// read startup file // read startup file
string startup_path_name("TFitPofB_startup.xml"); string startup_path_name("TFitPofB_startup.xml");
@ -1121,7 +1155,10 @@ TLondon1D3LS::TLondon1D3LS() : fCalcNeeded(true), fFirstCall(true), fLastThreeCh
fParForPofB.push_back(startupHandler->GetDeltat()); fParForPofB.push_back(startupHandler->GetDeltat());
fParForPofB.push_back(startupHandler->GetDeltaB()); fParForPofB.push_back(startupHandler->GetDeltaB());
fParForPofB.push_back(0.0); fParForPofB.push_back(0.0); // Energy
fParForPofB.push_back(0.0); // Bkg-Field
fParForPofB.push_back(0.005); // Bkg-width
fParForPofB.push_back(0.0); // Bkg-weight
fImpProfile = new TTrimSPData(rge_path, energy_vec); fImpProfile = new TTrimSPData(rge_path, energy_vec);
@ -1148,25 +1185,26 @@ TLondon1D3LS::TLondon1D3LS() : fCalcNeeded(true), fFirstCall(true), fLastThreeCh
double TLondon1D3LS::operator()(double t, const vector<double> &par) const { double TLondon1D3LS::operator()(double t, const vector<double> &par) const {
assert(par.size() == 12); assert(par.size() == 9 || par.size() == 13);
if(t<0.0) if(t<0.0)
return cos(par[0]*0.017453293); return cos(par[0]*0.017453293);
bool bkg_fraction_changed(false);
bool weights_changed(false);
// check if the function is called the first time and if yes, read in parameters // check if the function is called the first time and if yes, read in parameters
if(fFirstCall){ if(fFirstCall){
fPar = par; fPar = par;
/* for (unsigned int i(0); i<fPar.size(); i++){
cout << "fPar[" << i << "] = " << fPar[i] << endl;
}
*/
for (unsigned int i(2); i<fPar.size(); i++){ for (unsigned int i(2); i<fPar.size(); i++){
fParForBofZ.push_back(fPar[i]); fParForBofZ.push_back(fPar[i]);
// cout << "fParForBofZ[" << i-2 << "] = " << fParForBofZ[i-2] << endl;
} }
fFirstCall=false; fFirstCall=false;
bkg_fraction_changed = true;
if (par.size() == 13)
weights_changed = true;
} }
// check if any parameter has changed // check if any parameter has changed
@ -1182,9 +1220,11 @@ double TLondon1D3LS::operator()(double t, const vector<double> &par) const {
only_phase_changed = true; only_phase_changed = true;
} else { } else {
only_phase_changed = false; only_phase_changed = false;
if (i == 3 || i == 4 || i == 5 || i == 6)
bkg_fraction_changed = true;
if (i == 9 || i == 10 || i == 11 || i == 12)
weights_changed = true;
} }
if (i == fPar.size()-3 || i == fPar.size()-2 || i == fPar.size()-1)
fLastThreeChanged = true;
} }
} }
@ -1205,23 +1245,38 @@ double TLondon1D3LS::operator()(double t, const vector<double> &par) const {
fParForBofZ[i-2] = par[i]; fParForBofZ[i-2] = par[i];
fParForPofB[2] = par[1]; // energy fParForPofB[2] = par[1]; // energy
fParForPofB[3] = par[2]; // Bkg-Field
//fParForPofB[4] = 0.005; // Bkg-width (in principle zero)
if(fLastThreeChanged) { if(weights_changed) {
vector<double> interfaces; vector<double> interfaces;
interfaces.push_back(par[3]+par[4]); interfaces.push_back(par[3]+par[4]);
interfaces.push_back(par[3]+par[4]+par[5]); interfaces.push_back(par[3]+par[4]+par[5]);
interfaces.push_back(par[3]+par[4]+par[5]+par[6]);
vector<double> weights; vector<double> weights;
for(unsigned int i(par.size()-3); i<par.size(); i++) for(unsigned int i(9); i<13; i++)
weights.push_back(par[i]); weights.push_back(par[i]);
// cout << "Weighting has changed, re-calculating n(z) now..." << endl; // cout << "Weighting has changed, re-calculating n(z) now..." << endl;
fImpProfile->WeightLayers(par[1], interfaces, weights); fImpProfile->WeightLayers(par[1], interfaces, weights);
interfaces.clear();
weights.clear();
} }
TLondon1D_3LS BofZ3S(fParForBofZ); if(bkg_fraction_changed || weights_changed) {
vector<double> interfaces;
interfaces.push_back(par[3]);// dead layer
interfaces.push_back(par[3] + par[4] + par[5] + par[6]);// dead layer + first layer + second layer + third layer
fParForPofB[5] = fImpProfile->LayerFraction(par[1], 1, interfaces) + // Fraction of muons in the deadlayer
fImpProfile->LayerFraction(par[1], 3, interfaces); // Fraction of muons in the substrate
interfaces.clear();
}
TLondon1D_3LS BofZ3(fParForBofZ);
fPofB->UnsetPBExists(); fPofB->UnsetPBExists();
fPofB->Calculate(&BofZ3S, fImpProfile, fParForPofB); fPofB->Calculate(&BofZ3, fImpProfile, fParForPofB);
fPofT->DoFFT(); fPofT->DoFFT();
}/* else { }/* else {
@ -1231,7 +1286,6 @@ double TLondon1D3LS::operator()(double t, const vector<double> &par) const {
fPofT->CalcPol(fParForPofT); fPofT->CalcPol(fParForPofT);
fCalcNeeded = false; fCalcNeeded = false;
fLastThreeChanged = false;
} }
return fPofT->Eval(t); return fPofT->Eval(t);
@ -1396,183 +1450,6 @@ double TLondon1D3LS::operator()(double t, const vector<double> &par) const {
// //
// } // }
//------------------
// Constructor of the TLondon1D3LSub class -- reading available implantation profiles and
// creates (a pointer to) the TPofTCalc object (with the FFT plan)
//------------------
TLondon1D3LSub::TLondon1D3LSub() : fCalcNeeded(true), fFirstCall(true), fWeightsChanged(true) {
// omp_set_nested(1);
// omp_set_dynamic(1);
// omp_set_num_threads(4);
// read startup file
string startup_path_name("TFitPofB_startup.xml");
TSAXParser *saxParser = new TSAXParser();
TFitPofBStartupHandler *startupHandler = new TFitPofBStartupHandler();
saxParser->ConnectToHandler("TFitPofBStartupHandler", startupHandler);
int status (saxParser->ParseFile(startup_path_name.c_str()));
// check for parse errors
if (status) { // error
cerr << endl << "**ERROR** reading/parsing TFitPofB_startup.xml failed." \
<< endl << "**ERROR** Please make sure that the file exists in the local directory and it is set up correctly!" \
<< endl;
assert(false);
}
fNSteps = startupHandler->GetNSteps();
fWisdom = startupHandler->GetWisdomFile();
string rge_path(startupHandler->GetDataPath());
map<double, string> energy_vec(startupHandler->GetEnergies());
fParForPofT.push_back(0.0);
fParForPofT.push_back(startupHandler->GetDeltat());
fParForPofT.push_back(startupHandler->GetDeltaB());
fParForPofB.push_back(startupHandler->GetDeltat());
fParForPofB.push_back(startupHandler->GetDeltaB());
fParForPofB.push_back(0.0);
fImpProfile = new TTrimSPData(rge_path, energy_vec);
fPofB = new TPofBCalc(fParForPofB);
fPofT = new TPofTCalc(fPofB, fWisdom, fParForPofT);
// clean up
if (saxParser) {
delete saxParser;
saxParser = 0;
}
if (startupHandler) {
delete startupHandler;
startupHandler = 0;
}
}
//------------------
// TLondon1D3LSub-Method that calls the procedures to create B(z), p(B) and P(t)
// It finally returns P(t) for a given t.
// Parameters: all the parameters for the function to be fitted through TLondon1D3LSub
//------------------
double TLondon1D3LSub::operator()(double t, const vector<double> &par) const {
assert(par.size() == 15);
if(t<0.0)
return cos(par[0]*0.017453293);
// check if the function is called the first time and if yes, read in parameters
//#pragma omp critical
//{
if(fFirstCall){
fPar = par;
/* for (unsigned int i(0); i<fPar.size(); i++){
cout << "fPar[" << i << "] = " << fPar[i] << endl;
}
*/
for (unsigned int i(2); i<fPar.size(); i++){
fParForBofZ.push_back(fPar[i]);
// cout << "fParForBofZ[" << i-2 << "] = " << fParForBofZ[i-2] << endl;
}
fFirstCall=false;
}
// check if any parameter has changed
bool par_changed(false);
bool only_phase_changed(false);
for (unsigned int i(0); i<fPar.size(); i++) {
if( fPar[i]-par[i] ) {
fPar[i] = par[i];
par_changed = true;
if (i == 0) {
only_phase_changed = true;
} else {
only_phase_changed = false;
}
if (i == fPar.size()-5 || i == fPar.size()-4 || i == fPar.size()-3 || i == fPar.size()-2)
fWeightsChanged = true;
}
}
if (par_changed)
fCalcNeeded = true;
// if model parameters have changed, recalculate B(z), P(B) and P(t)
if (fCalcNeeded) {
fParForPofT[0] = par[0]; // phase
if(!only_phase_changed) {
// cout << " Parameters have changed, (re-)calculating p(B) and P(t) now..." << endl;
for (unsigned int i(2); i<par.size(); i++)
fParForBofZ[i-2] = par[i];
fParForPofB[2] = par[1]; // energy
/* DEBUG ---------------------------
for(unsigned int i(0); i<fParForBofZ.size(); i++) {
cout << "ParForBofZ[" << i << "] = " << fParForBofZ[i] << endl;
}
for(unsigned int i(0); i<fParForPofB.size(); i++) {
cout << "ParForPofB[" << i << "] = " << fParForPofB[i] << endl;
}
for(unsigned int i(0); i<fParForPofT.size(); i++) {
cout << "ParForPofT[" << i << "] = " << fParForPofT[i] << endl;
}
------------------------------------*/
vector<double> interfaces;
interfaces.push_back(par[3]+par[4]);
interfaces.push_back(par[3]+par[4]+par[5]);
interfaces.push_back(par[3]+par[4]+par[5]+par[6]);
if(fWeightsChanged) {
vector<double> weights;
for(unsigned int i(par.size()-5); i<(par.size()-1); i++)
weights.push_back(par[i]);
// cout << "Weighting has changed, re-calculating n(z) now..." << endl;
fImpProfile->WeightLayers(par[1], interfaces, weights);
}
TLondon1D_3L BofZ3(fParForBofZ);
fPofB->UnsetPBExists();
fPofB->Calculate(&BofZ3, fImpProfile, fParForPofB);
// Add background contribution from the substrate
fPofB->AddBackground(par[2], par[14], fImpProfile->LayerFraction(par[1], 4, interfaces));
// FourierTransform of P(B)
fPofT->DoFFT();
}/* else {
cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl;
}*/
fPofT->CalcPol(fParForPofT);
fCalcNeeded = false;
fWeightsChanged = false;
}
//}
return fPofT->Eval(t);
}
double TLondon1D3Lestimate::operator()(double z, const vector<double>& par) const { double TLondon1D3Lestimate::operator()(double z, const vector<double>& par) const {

8
src/external/TFitPofB-lib/classes/TPofBCalc.cpp vendored
View File

@ -205,7 +205,7 @@ void TPofBCalc::Calculate(const TBofZCalcInverse *BofZ, const TTrimSPData *dataT
for (i = firstZerosEnd; i<=lastZerosStart; i++) for (i = firstZerosEnd; i<=lastZerosStart; i++)
fPB[i] /= pBsum; fPB[i] /= pBsum;
if(para.size() == 6) if(para.size() == 6 && para[5] != 0.0)
AddBackground(para[3], para[4], para[5]); AddBackground(para[3], para[4], para[5]);
} }
@ -433,6 +433,7 @@ void TPofBCalc::Calculate(const TBulkVortexFieldCalc *vortexLattice, const vecto
} else if (para.size() == 7 && para[6] == 2.0 && para[5] != 0.0 && vortexLattice->IsTriangular()) { } else if (para.size() == 7 && para[6] == 2.0 && para[5] != 0.0 && vortexLattice->IsTriangular()) {
// weight distribution with Lorentzian around vortex-cores // weight distribution with Lorentzian around vortex-cores
double Rsq1, Rsq2, Rsq3, Rsq4, Rsq5, Rsq6, sigmaSq(para[5]*para[5]); double Rsq1, Rsq2, Rsq3, Rsq4, Rsq5, Rsq6, sigmaSq(para[5]*para[5]);
// ofstream of("LorentzWeight.dat");
for (unsigned int j(0); j < numberOfSteps_2; ++j) { for (unsigned int j(0); j < numberOfSteps_2; ++j) {
for (unsigned int i(0); i < numberOfSteps_2; ++i) { for (unsigned int i(0); i < numberOfSteps_2; ++i) {
fill_index = static_cast<unsigned int>(ceil(fabs((vortexFields[i + numberOfSteps*j]/fDB)))); fill_index = static_cast<unsigned int>(ceil(fabs((vortexFields[i + numberOfSteps*j]/fDB))));
@ -450,9 +451,14 @@ void TPofBCalc::Calculate(const TBulkVortexFieldCalc *vortexLattice, const vecto
+ (numberOfSteps_2 - j)*(numberOfSteps_2 - j))/static_cast<double>(numberOfStepsSq); + (numberOfSteps_2 - j)*(numberOfSteps_2 - j))/static_cast<double>(numberOfStepsSq);
fPB[fill_index] += 1.0/(1.0+sigmaSq*Rsq1) + 1.0/(1.0+sigmaSq*Rsq2) + 1.0/(1.0+sigmaSq*Rsq3) \ fPB[fill_index] += 1.0/(1.0+sigmaSq*Rsq1) + 1.0/(1.0+sigmaSq*Rsq2) + 1.0/(1.0+sigmaSq*Rsq3) \
+ 1.0/(1.0+sigmaSq*Rsq4) + 1.0/(1.0+sigmaSq*Rsq5) + 1.0/(1.0+sigmaSq*Rsq6); + 1.0/(1.0+sigmaSq*Rsq4) + 1.0/(1.0+sigmaSq*Rsq5) + 1.0/(1.0+sigmaSq*Rsq6);
// of << 1.0/(1.0+sigmaSq*Rsq1) + 1.0/(1.0+sigmaSq*Rsq2) + 1.0/(1.0+sigmaSq*Rsq3) \
// + 1.0/(1.0+sigmaSq*Rsq4) + 1.0/(1.0+sigmaSq*Rsq5) + 1.0/(1.0+sigmaSq*Rsq6) << " ";
} }
} }
// of << endl;
} }
// of.close();
} else { } else {
for (unsigned int j(0); j < numberOfSteps_2; ++j) { for (unsigned int j(0); j < numberOfSteps_2; ++j) {
for (unsigned int i(0); i < numberOfSteps_2; ++i) { for (unsigned int i(0); i < numberOfSteps_2; ++i) {

29
src/external/TFitPofB-lib/include/TLondon1D.h vendored
View File

@ -107,7 +107,6 @@ private:
mutable vector<double> fParForPofB; mutable vector<double> fParForPofB;
string fWisdom; string fWisdom;
unsigned int fNSteps; unsigned int fNSteps;
mutable bool fLastTwoChanged;
ClassDef(TLondon1D2L,1) ClassDef(TLondon1D2L,1)
}; };
@ -184,7 +183,6 @@ private:
mutable vector<double> fParForPofB; mutable vector<double> fParForPofB;
string fWisdom; string fWisdom;
unsigned int fNSteps; unsigned int fNSteps;
mutable bool fLastThreeChanged;
ClassDef(TLondon1D3L,1) ClassDef(TLondon1D3L,1)
}; };
@ -211,7 +209,6 @@ private:
mutable vector<double> fParForPofB; mutable vector<double> fParForPofB;
string fWisdom; string fWisdom;
unsigned int fNSteps; unsigned int fNSteps;
mutable bool fLastThreeChanged;
ClassDef(TLondon1D3LS,1) ClassDef(TLondon1D3LS,1)
}; };
@ -241,32 +238,6 @@ private:
// ClassDef(TLondon1D4L,1) // ClassDef(TLondon1D4L,1)
// }; // };
class TLondon1D3LSub : public PUserFcnBase {
public:
// default constructor
TLondon1D3LSub();
~TLondon1D3LSub();
double operator()(double, const vector<double>&) const;
private:
mutable vector<double> fPar;
TTrimSPData *fImpProfile;
TPofBCalc *fPofB;
TPofTCalc *fPofT;
mutable bool fCalcNeeded;
mutable bool fFirstCall;
mutable vector<double> fParForPofT;
mutable vector<double> fParForBofZ;
mutable vector<double> fParForPofB;
string fWisdom;
unsigned int fNSteps;
mutable bool fWeightsChanged;
ClassDef(TLondon1D3LSub,1)
};
// Class for fitting directly B(z) without any P(B)-calculation // Class for fitting directly B(z) without any P(B)-calculation
class TLondon1D3Lestimate : public PUserFcnBase { class TLondon1D3Lestimate : public PUserFcnBase {

1
src/external/TFitPofB-lib/include/TLondon1DLinkDef.h vendored
View File

@ -44,7 +44,6 @@
#pragma link C++ class TLondon1D3L+; #pragma link C++ class TLondon1D3L+;
#pragma link C++ class TLondon1D3LS+; #pragma link C++ class TLondon1D3LS+;
//#pragma link C++ class TLondon1D4L+; //#pragma link C++ class TLondon1D4L+;
#pragma link C++ class TLondon1D3LSub+;
#pragma link C++ class TLondon1D3Lestimate; #pragma link C++ class TLondon1D3Lestimate;