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

Fixed bugs and added a few comments after first testing.

Browse Source
This commit is contained in:
Bastian M. Wojek
2008-05-24 15:05:18 +00:00
parent 728ec00816
commit a8046f80d0
12 changed files with 642 additions and 175 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@ -5,47 +5,70 @@
Author: Bastian M. Wojek
e-mail: bastian.wojek@psi.ch
2008/05/23
2008/05/24
***************************************************************************/
#include "TPofBCalc.h"
#include <cmath>
#include <iostream>
#include <fstream>
//#include <cstdio>
//-----------
// Constructor that does the P(B) calculation for given B(z) and n(z)
// Parameters: dt[us], dB[G], Energy[keV]
//-----------
TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, const vector<double> &para ) {
// Parameters: dt[us], dB[G], Energy[keV]
fBmin = BofZ.GetBmin();
fBmax = BofZ.GetBmax();
double BB, BBnext;
double zm, zp, zNextm, zNextp, dz;
// fill not used Bs before Bmin with 0.0
for ( BB = 0.0 ; BB < BofZ.GetBmin() ; BB += para[1] ) {
for ( BB = 0.0 ; BB < fBmin ; BB += para[1] ) {
fB.push_back(BB);
fPB.push_back(0.0);
}
unsigned int firstZerosEnd(fB.size());
// calculate p(B) from B(z)
vector<double> bofzZ(BofZ.DataZ());
vector<double> bofzBZ(BofZ.DataBZ());
double ddZ(BofZ.GetdZ());
double Bmax(BofZ.GetBmax());
/* USED FOR DEBUGGING-----------------------------------
cout << "Bmin = " << fBmin << ", Bmax = " << fBmax << endl;
char debugfile[50];
int n = sprintf (debugfile, "test_Bz_%f.dat", fBmin);
if (n > 0) {
ofstream of(debugfile);
for (unsigned int i(0); i<bofzZ.size(); i++) {
of << bofzZ[i] << " " << bofzBZ[i] << endl;
}
of.close();
}
---------------------------------------------------------*/
double nn;
for ( ; BB <= Bmax ; BB += para[1]) {
for ( ; BB <= fBmax ; BB += para[1]) {
BBnext = BB + para[1];
fB.push_back(BB);
fPB.push_back(0.0);
for ( unsigned int j(0); j < bofzZ.size() - 1; j++ ) {
if ( bofzBZ[j] >= BB && bofzBZ[j+1] <= BB ) {
zm = (BB-bofzBZ[j])*ddZ/(bofzBZ[j+1]-bofzBZ[j]) + bofzZ[j];
for (unsigned int k(0); k < j; k++) {
if ( ( bofzBZ[j-k] <= BBnext && bofzBZ[j-k-1] >= BBnext ) ) {
// cout << "1 " << j << " " << k << endl;
@ -53,17 +76,17 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons
break;
}
}
dz = zNextm-zm;
nn = dataTrimSP.GetNofZ(zm, para[2]);
if (nn != -1.0) {
// cout << "zNext = " << zNextm << ", zm = " << zm << ", dz = " << dz << endl;
*(fPB.end()-1) += nn*fabs(dz/para[1]);
}
} else if (bofzBZ[j] <= BB && bofzBZ[j+1] >= BB ) {
zp = (BB-bofzBZ[j])*ddZ/(bofzBZ[j+1]-bofzBZ[j]) + bofzZ[j];
for (unsigned int k(0); k < bofzZ.size() - j - 1; k++) {
if ( ( bofzBZ[j+k] <= BBnext && bofzBZ[j+k+1] >= BBnext ) ) {
// cout << "2 " << j << " " << k << endl;
@ -71,7 +94,7 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons
break;
}
}
dz = zNextp-zp;
nn = dataTrimSP.GetNofZ(zp, para[2]);
if (nn != -1.0) {
@ -80,33 +103,33 @@ TPofBCalc::TPofBCalc( const TBofZCalc &BofZ, const TTrimSPData &dataTrimSP, cons
}
}
}
}
unsigned int lastZerosStart(fB.size());
// fill not used Bs after Bext with 0.0
double BmaxFFT(1.0/gBar/para[0]);
//// cout << "N = " << int(BmaxFFT/para[1]+1.0) << endl;
// cout << "N = " << int(BmaxFFT/para[1]+1.0) << endl;
for ( ; BB <= BmaxFFT ; BB += para[1] ) {
fB.push_back(BB);
fPB.push_back(0.0);
}
// make sure that we have an even number of elements in p(B) for FFT
// make sure that we have an even number of elements in p(B) for FFT, so we do not have to care later
if (fB.size() % 2) {
fB.push_back(BB);
fPB.push_back(0.0);
}
cout << "size of B = " << fB.size() << ", size of p(B) = " << fPB.size() << endl;
// normalize pB
// cout << "size of B = " << fB.size() << ", size of p(B) = " << fPB.size() << endl;
// normalize p(B)
double pBsum = 0.0;
for (unsigned int i(firstZerosEnd); i<lastZerosStart; i++)
pBsum += fPB[i];