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

Added a TTrimSPDataHandler method and changed TFitPofB to TUserLondon since it has become quite specific.

Browse Source
This commit is contained in:
Bastian M. Wojek 2008-05-25 14:40:13 +00:00
parent a8046f80d0
commit c16916b563
10 changed files with 185 additions and 43 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/external/TFitPofB-lib/classes/Makefile.TFitPofB vendored
View File

@ -12,7 +12,7 @@ OBJS += TTrimSPDataHandler.o
OBJS += TBofZCalc.o
OBJS += TPofBCalc.o
OBJS += TPofTCalc.o
OBJS += TFitPofB.o
OBJS += TUserLondon.o
SHLIB = libTFitPofB.so

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

@ -5,7 +5,7 @@
Author: Bastian M. Wojek
e-mail: bastian.wojek@psi.ch
2008/05/24
2008/05/25
***************************************************************************/
@ -13,7 +13,11 @@
#include <cmath>
#include <iostream>
#include <fstream>
//#include <cstdio>
/* USED FOR DEBUGGING-----------------------------------
#include <cstdio>
#include <ctime>
-------------------------------------------------------*/
//-----------
// Constructor that does the P(B) calculation for given B(z) and n(z)
@ -46,8 +50,11 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons
/* USED FOR DEBUGGING-----------------------------------
cout << "Bmin = " << fBmin << ", Bmax = " << fBmax << endl;
time_t seconds;
seconds = time (NULL);
char debugfile[50];
int n = sprintf (debugfile, "test_Bz_%f.dat", fBmin);
int n = sprintf (debugfile, "test_Bz_%f_%ld.dat", fBmin, seconds);
if (n > 0) {
ofstream of(debugfile);
@ -57,6 +64,19 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons
}
of.close();
}
char debugfile1[50];
int n1 = sprintf (debugfile1, "test_NZ_%f_%ld.dat", para[2], seconds);
if (n1 > 0) {
ofstream of1(debugfile1);
for (unsigned int i(0); i<dataTrimSP.DataZ(para[2]).size(); i++) {
of1 << dataTrimSP.DataZ(para[2])[i] << " " << dataTrimSP.DataNZ(para[2])[i] << " " << dataTrimSP.OrigDataNZ(para[2])[i] << endl;
}
of1.close();
}
---------------------------------------------------------*/
double nn;

16
src/external/TFitPofB-lib/classes/TPofTCalc.cpp vendored
View File

@ -5,16 +5,20 @@
Author: Bastian M. Wojek
e-mail: bastian.wojek@psi.ch
2008/05/24
2008/05/25
***************************************************************************/
#include "TPofTCalc.h"
#include "fftw3.h"
#include <cmath>
#include <cstdio>
#include <fstream>
#include <iostream>
#include <cstdio>
/* USED FOR DEBUGGING -----------------------
#include <ctime>
#include <fstream>
--------------------------------------------*/
//------------------
// Constructor of the TPofTCalc class - it creates the FFT plan
@ -61,11 +65,15 @@ void TPofTCalc::DoFFT(const TPofBCalc &PofB) {
vector<double> pB(PofB.DataPB());
/* USED FOR DEBUGGING -----------------------
time_t seconds;
seconds = time(NULL);
vector<double> B(PofB.DataB());
double Bmin(PofB.GetBmin());
char debugfile[50];
int n = sprintf (debugfile, "test_PB_%f.dat", Bmin);
int n = sprintf (debugfile, "test_PB_%f_%ld.dat", Bmin, seconds);
if (n > 0) {
ofstream of(debugfile);

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

@ -5,7 +5,7 @@
Author: Bastian M. Wojek
e-mail: bastian.wojek@psi.ch
2008/05/24
2008/05/25
***************************************************************************/
@ -65,6 +65,8 @@ TTrimSPData::TTrimSPData(const string &path, vector<string> &energyVec) {
}
}
fOrigDataNZ = fDataNZ;
}
//---------------------
@ -86,7 +88,7 @@ vector<double> TTrimSPData::DataZ(double e) const {
}
//---------------------
// Method returning n(z)-vector calculated by trim.SP for given energy[keV]
// Method returning actual n(z)-vector calculated by trim.SP and potentially altered by the WeightLayers-method for given energy[keV]
//---------------------
vector<double> TTrimSPData::DataNZ(double e) const {
@ -102,6 +104,79 @@ vector<double> TTrimSPData::DataNZ(double e) const {
}
//---------------------
// Method returning original n(z)-vector calculated by trim.SP for given energy[keV]
//---------------------
vector<double> TTrimSPData::OrigDataNZ(double e) const {
for(unsigned int i(0); i<fEnergy.size(); i++) {
if(!(fEnergy[i] - e)) {
return fOrigDataNZ[i];
}
}
// default
cout << "No implantation profile available for the specified energy... You get back the first one." << endl;
return fOrigDataNZ[0];
}
//---------------------
// Method putting different weight to different layers of your thin film
// Parameters: Implantation Energy[keV], Interfaces[nm], Weights [0.0 <= w[i] <= 1.0]
// Example: 25.0, (50, 100), (1.0, 0.33, 1.0)
// at 25keV consider 3 layers, where the first ends after 50nm, the second after 100nm (these are NOT the layer thicknesses!!)
// the first and last layers get the full n(z), where only one third of the muons in the second layer will be taken into account
//---------------------
void TTrimSPData::WeightLayers(double e, const vector<double>& interface, const vector<double>& weight) const {
if(weight.size()-interface.size()-1) {
cout << "For the weighting the number of interfaces has to be one less than the number of weights!" << endl;
cout << "No weighting of the implantation profile will be done unless you take care of that!" << endl;
return;
}
for(unsigned int i(0); i<interface.size(); i++) {
if (interface[i]<0.0) {
cout << "One of your layer interfaces has a negative coordinate! - No weighting will be done!" << endl;
return;
}
else if (i>1) {
if (interface[i]<interface[i-1]) {
cout << "The specified interfaces appear to be not in ascending order! - No weighting will be done!" << endl;
return;
}
}
}
for(unsigned int i(0); i<weight.size(); i++) {
if (weight[i]>1.0 || weight[i]<0.0) {
cout << "At least one of the specified weights is out of range - no weighting will be done!" << endl;
return;
}
}
for(unsigned int i(0); i<fEnergy.size(); i++) {
if(!(fEnergy[i] - e)) {
unsigned int k(0);
for(unsigned int j(0); j<fDataZ[i].size(); j++) {
if(k<interface.size()) {
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];
}
}
}
}
//---------------------
// Method returning n(z) for given z[nm] and energy[keV]
//---------------------

54
src/external/TFitPofB-lib/classes/TFitPofB.cpp → src/external/TFitPofB-lib/classes/TUserLondon.cpp vendored
View File

@ -1,24 +1,25 @@
/***************************************************************************
TFitPofB.cpp
TUserLondon.cpp
Author: Bastian M. Wojek
e-mail: bastian.wojek@psi.ch
2008/05/24
2008/05/25
***************************************************************************/
#include "TFitPofB.h"
#include "TUserLondon.h"
#include <iostream>
using namespace std;
//------------------
// Constructor of the TFitPofB class -- reading available implantation profiles and
// Constructor of the TUserLondon class -- reading available implantation profiles and
// creates (a pointer to) the TPofTCalc object (with the FFT plan)
//------------------
TFitPofB::TFitPofB(const vector<unsigned int> &parNo, const vector<double> &par) : fCalcNeeded(true) {
TUserLondon::TUserLondon(const vector<unsigned int> &parNo, const vector<double> &par)
: fCalcNeeded(true), fLastTwoChanged(true), fLastThreeChanged(true) {
for(unsigned int i(0); i<parNo.size(); i++) {
fPar.push_back(par[parNo[i]-1]);
@ -47,10 +48,10 @@ TFitPofB::TFitPofB(const vector<unsigned int> &parNo, const vector<double> &par)
}
//------------------
// Destructor of the TFitPofB class -- cleaning up
// Destructor of the TUserLondon class -- cleaning up
//------------------
TFitPofB::~TFitPofB() {
TUserLondon::~TUserLondon() {
fPar.clear();
delete fImpProfile;
fImpProfile = 0;
@ -61,10 +62,10 @@ TFitPofB::~TFitPofB() {
//------------------
// 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 TFitPofB
// Parameters: all the parameters for the function to be fitted through TUserLondon
//------------------
double TFitPofB::Eval(double t, const vector<double> &par) const {
double TUserLondon::Eval(double t, const vector<double> &par) const {
// check if any parameter has changed
@ -77,12 +78,18 @@ double TFitPofB::Eval(double t, const vector<double> &par) const {
par_changed = true;
if(i == 0 || i == 2 || i == 3 || i == 4) {
cout << "You are varying the model-flag, dt, dB or E! These parameters have to be fixed! Quitting..." << endl;
exit(0);
exit(-1);
} else if (i == 1) {
only_phase_changed = true;
} else {
only_phase_changed = false;
}
if(i == fPar.size()-3) {
fLastThreeChanged = true;
} else if (i == fPar.size()-2 || i == fPar.size()-1) {
fLastTwoChanged = true;
fLastThreeChanged = true;
}
}
}
@ -140,6 +147,17 @@ double TFitPofB::Eval(double t, const vector<double> &par) const {
case 2:
{
// cout << "Found the 1D-London model.2L" << endl;
if(fLastTwoChanged) {
vector<double> interfaces;
interfaces.push_back(par[6]+par[7]);
vector<double> weights;
for(unsigned int i(11); i<par.size(); i++)
weights.push_back(par[i]);
cout << "Weighting has changed, re-calculating n(z) now..." << endl;
fImpProfile->WeightLayers(par[4], interfaces, weights);
}
TLondon1D_2L BofZ2(par_for_BofZ);
TPofBCalc PofB2(BofZ2, *fImpProfile, par_for_PofB);
fPofT->DoFFT(PofB2);
@ -148,6 +166,18 @@ double TFitPofB::Eval(double t, const vector<double> &par) const {
case 3:
{
// cout << "Found the 1D-London model.3L" << endl;
if(fLastThreeChanged) {
vector<double> interfaces;
interfaces.push_back(par[6]+par[7]);
interfaces.push_back(par[6]+par[7]+par[8]);
vector<double> weights;
for(unsigned int i(12); i<par.size(); i++)
weights.push_back(par[i]);
cout << "Weighting has changed, re-calculating n(z) now..." << endl;
fImpProfile->WeightLayers(par[4], interfaces, weights);
}
TLondon1D_3L BofZ3(par_for_BofZ);
TPofBCalc PofB3(BofZ3, *fImpProfile, par_for_PofB);
fPofT->DoFFT(PofB3);
@ -155,7 +185,7 @@ double TFitPofB::Eval(double t, const vector<double> &par) const {
break;
default:
cout << "The user function you specified with the first parameter, does not exist. Quitting..." << endl;
exit(0);
exit(-1);
}
} else {
cout << "Only the phase parameter has changed, (re-)calculating P(t) now..." << endl;
@ -164,6 +194,8 @@ double TFitPofB::Eval(double t, const vector<double> &par) const {
fPofT->CalcPol(par_for_PofT);
fCalcNeeded = false;
fLastTwoChanged = false;
fLastThreeChanged = false;
}

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

@ -5,7 +5,7 @@
Author: Bastian M. Wojek
e-mail: bastian.wojek@psi.ch
2008/05/24
2008/05/25
***************************************************************************/
@ -15,6 +15,8 @@
#include "TBofZCalc.h"
#include "TTrimSPDataHandler.h"
#define gBar 0.0135538817
class TPofBCalc {
public:
@ -35,7 +37,6 @@ private:
vector<double> fPB;
double fBmin;
double fBmax;
static const double gBar = 0.0135538817;
};

6
src/external/TFitPofB-lib/include/TPofTCalc.h vendored
View File

@ -5,7 +5,7 @@
Author: Bastian M. Wojek
e-mail: bastian.wojek@psi.ch
2008/05/24
2008/05/25
***************************************************************************/
@ -15,6 +15,8 @@
#include "TPofBCalc.h"
#include "fftw3.h"
#define PI 3.14159265358979323846
class TPofTCalc {
public:
@ -35,8 +37,6 @@ private:
vector<double> fPT;
double fTBin;
unsigned int fNFFT;
static const double PI = 3.14159265358979323846;
static const double gBar = 0.0135538817;
};

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

@ -5,7 +5,7 @@
Author: Bastian M. Wojek
e-mail: bastian.wojek@psi.ch
2008/05/24
2008/05/25
***************************************************************************/
@ -30,12 +30,15 @@ public:
vector<double> Energy() const {return fEnergy;}
vector<double> DataZ(double) const;
vector<double> DataNZ(double) const;
vector<double> OrigDataNZ(double) const;
void WeightLayers(double, const vector<double>&, const vector<double>&) const;
double GetNofZ(double, double) const;
private:
vector<double> fEnergy;
vector< vector<double> > fDataZ;
vector< vector<double> > fDataNZ;
mutable vector< vector<double> > fDataNZ;
vector< vector<double> > fOrigDataNZ;
};
#endif // _TTrimSPDataHandler_H_

18
src/external/TFitPofB-lib/include/TFitPofB.h → src/external/TFitPofB-lib/include/TUserLondon.h vendored
View File

@ -1,24 +1,24 @@
/***************************************************************************
TFitPofB.h
TUserLondon.h
Author: Bastian M. Wojek
e-mail: bastian.wojek@psi.ch
2008/05/24
2008/05/25
***************************************************************************/
#ifndef _TFitPofB_H_
#define _TFitPofB_H_
#ifndef _TUserLondon_H_
#define _TUserLondon_H_
#include "TPofTCalc.h"
class TFitPofB {
class TUserLondon {
public:
TFitPofB(const vector<unsigned int>& , const vector<double>&);
~TFitPofB();
TUserLondon(const vector<unsigned int>& , const vector<double>&);
~TUserLondon();
double Eval(double, const vector<double>&) const;
@ -27,7 +27,9 @@ private:
TTrimSPData *fImpProfile;
TPofTCalc *fPofT;
mutable bool fCalcNeeded;
mutable bool fLastTwoChanged;
mutable bool fLastThreeChanged;
};
#endif //_TFitPofB_H_
#endif //_TUserLondon_H_

13
src/external/TFitPofB-lib/test/test.cpp vendored
View File

@ -1,4 +1,4 @@
#include "TFitPofB.h"
#include "TUserLondon.h"
#include <iostream>
#include <fstream>
@ -101,8 +101,8 @@ int main(){
*/
unsigned int parNo_arr[] = {1, 3, 5, 7, 9, 11, 12, 13, 14, 15, 16, 17};
double par_arr[] = {3.0, 999.0, 0.0, 999.0, 0.01, 999.0, 0.05, 999.0, 20.0, 999.0, 100.0, 10.0, 80.0, 50.0, 80.0, 180.0, 500.0};
unsigned int parNo_arr[] = {1, 3, 5, 7, 9, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20};
double par_arr[] = {3.0, 999.0, 0.0, 999.0, 0.01, 999.0, 0.05, 999.0, 25.0, 999.0, 100.0, 10.0, 65.0, 50.0, 75.0, 180.0, 500.0, 1.0, 0.3, 1.0};
vector<unsigned int> parNo_vec(parNo_arr, parNo_arr+(sizeof(parNo_arr)/sizeof(parNo_arr[0])));
vector<double> par_vec(par_arr, par_arr+(sizeof(par_arr)/sizeof(par_arr[0])));
@ -113,7 +113,7 @@ int main(){
par_vec_sub.push_back(par_vec[parNo_vec[i]-1]);
}
TFitPofB fitter(parNo_vec, par_vec);
TUserLondon fitter(parNo_vec, par_vec);
ofstream of01("test_fitter01.dat");
ofstream of02("test_fitter02.dat");
@ -179,7 +179,6 @@ int main(){
of07.close();
par_vec_sub[1] += 10.0;
par_vec_sub[7] = 75.0;
for (double i(0.); i<12.0; i+=0.003) {
of08 << i << " " << fitter.Eval(i, par_vec_sub) << endl;
@ -194,7 +193,9 @@ int main(){
}
of09.close();
par_vec_sub[1] += 10.0;
par_vec_sub[1] = 0.0;
par_vec_sub[11] = 500.0;
par_vec_sub[13] = 0.8;
for (double i(0.); i<12.0; i+=0.003) {
of10 << i << " " << fitter.Eval(i, par_vec_sub) << endl;