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

fixed some annoying little bug in the RGE handler. Added some checks if wrong energies are given. Fixed a typo

Browse Source
This commit is contained in:
suter_a 2012-07-05 15:31:20 +00:00
parent aa2e6cead0
commit de69ee8826
3 changed files with 21 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
src/external/Nonlocal/PNL_PippardFitter.cpp vendored
View File

@ -200,7 +200,6 @@ void PNL_PippardFitterGlobal::CalculateField(const std::vector<Double_t> &param)
Double_t xiP = XiP_T(param[4], param[3], param[1]);
Double_t preFactor = pow(xiP/(LambdaL_T(param[5], param[1])),2.0)*xiP/param[4];
// calculate the fFieldq vector, which is x/(x^2 + alpha k(x)), with alpha = xiP(T)^3/(lambdaL(T)^2 xiP(0)), and
// k(x) = 3/2 [(1+x^2) arctan(x) - x]/x^3, see lab-book p.137
Double_t x;
@ -443,6 +442,10 @@ Double_t PNL_PippardFitter::operator()(Double_t t, const std::vector<Double_t> &
// calculate field if parameter have changed
fPippardFitterGlobal->CalculateField(param);
Int_t energyIndex = fPippardFitterGlobal->GetEnergyIndex(param[0]);
if (energyIndex == -1) { // energy not found
cerr << endl << ">> PNL_PippardFitter::operator() energy " << param[0] << " not found. Will terminate." << endl;
assert(0);
}
// calcualte polarization
Bool_t done = false;
@ -450,11 +453,12 @@ Double_t PNL_PippardFitter::operator()(Double_t t, const std::vector<Double_t> &
Double_t z=0.0;
Int_t terminate = 0;
Double_t dz = 1.0;
do {
if (z < param[8]) { // z < dead-layer
dPol = fPippardFitterGlobal->GetMuoneStoppingDensity(energyIndex, z) * cos(GAMMA_MU * param[6] * t + param[7] * DEGREE2RAD);
dPol = fPippardFitterGlobal->GetMuonStoppingDensity(energyIndex, z) * cos(GAMMA_MU * param[6] * t + param[7] * DEGREE2RAD);
} else {
dPol = fPippardFitterGlobal->GetMuoneStoppingDensity(energyIndex, z) * cos(GAMMA_MU * param[6] * fPippardFitterGlobal->GetMagneticField(z-param[8]) * t + param[7] * DEGREE2RAD);
dPol = fPippardFitterGlobal->GetMuonStoppingDensity(energyIndex, z) * cos(GAMMA_MU * param[6] * fPippardFitterGlobal->GetMagneticField(z-param[8]) * t + param[7] * DEGREE2RAD);
}
z += dz;
pol += dPol;

2
src/external/Nonlocal/PNL_PippardFitter.h vendored
View File

@ -51,7 +51,7 @@ class PNL_PippardFitterGlobal
Bool_t IsValid() { return fValid; }
virtual void CalculateField(const std::vector<Double_t> &param) const;
virtual Int_t GetEnergyIndex(const Double_t energy) { return fRgeHandler->GetRgeEnergyIndex(energy); }
virtual Double_t GetMuoneStoppingDensity(const Int_t energyIndex, const Double_t z) const { return fRgeHandler->GetRgeValue(energyIndex, z); }
virtual Double_t GetMuonStoppingDensity(const Int_t energyIndex, const Double_t z) const { return fRgeHandler->GetRgeValue(energyIndex, z); }
virtual Double_t GetMagneticField(const Double_t z) const;
private:

38
src/external/Nonlocal/PNL_RgeHandler.cpp vendored
View File

@ -91,9 +91,17 @@ Double_t PNL_RgeHandler::GetRgeValue(const Int_t index, const Double_t dist)
{
Double_t rgeVal = 0.0;
UInt_t distIdx = (UInt_t)(dist/(fRgeDataList[index].stoppingDistance[1]-fRgeDataList[index].stoppingDistance[0]));
// find the proper position index
Int_t distIdx = -2;
for (UInt_t i=0; i<fRgeDataList[index].stoppingDistance.size(); i++) {
if (fRgeDataList[index].stoppingDistance[i] > dist) {
distIdx = i-1;
break;
}
}
if ((distIdx >= fRgeDataList[index].stoppingDistance.size()) || (distIdx < 0)) {
// get the proper n(z) value
if (distIdx < 0) {
rgeVal = 0.0;
} else {
rgeVal = fRgeDataList[index].stoppingAmplitude[distIdx] +
@ -101,17 +109,6 @@ Double_t PNL_RgeHandler::GetRgeValue(const Int_t index, const Double_t dist)
(dist-fRgeDataList[index].stoppingDistance[distIdx])/(fRgeDataList[index].stoppingDistance[distIdx+1]-fRgeDataList[index].stoppingDistance[distIdx]);
}
/*
cout << endl << "distIdx = " << distIdx << ", fRgeDataList[index].stoppingDistance.size() = " << fRgeDataList[index].stoppingDistance.size();
cout << endl << "fRgeDataList[index].stoppingAmplitude[distIdx] = " << fRgeDataList[index].stoppingAmplitude[distIdx];
cout << endl << "fRgeDataList[index].stoppingAmplitude[distIdx+1] = " << fRgeDataList[index].stoppingAmplitude[distIdx+1];
cout << endl << "dist = " << dist;
cout << endl << "fRgeDataList[index].stoppingDistance[distIdx] = " << fRgeDataList[index].stoppingDistance[distIdx];
cout << endl << "fRgeDataList[index].stoppingDistance[distIdx+1] = " << fRgeDataList[index].stoppingDistance[distIdx+1];
cout << endl << "rgeVal = " << rgeVal;
cout << endl << "----" << endl;
*/
return rgeVal;
}
@ -150,8 +147,7 @@ Bool_t PNL_RgeHandler::LoadRgeData(const PStringVector &rgeDataPathList, const P
{
ifstream fin;
PNL_RgeData data;
Int_t idx=0;
TString dataName, tstr;
TString tstr;
char line[512];
int result;
double dist, val;
@ -170,28 +166,20 @@ Bool_t PNL_RgeHandler::LoadRgeData(const PStringVector &rgeDataPathList, const P
data.energy = rgeDataEnergyList[i]/1000.0;
// read msr-file
idx = 0;
while (!fin.eof()) {
// read a line
fin.getline(line, sizeof(line));
idx++;
// ignore first line
if (idx == 1)
continue;
// ignore empty lines
tstr = line;
if (tstr.IsWhitespace())
continue;
// get values
result = sscanf(line, "%lf %lf", &dist, &val);
// check if data are valid
// check if data are valid, otherwise ignore it
if (result != 2) {
fin.close();
return false;
continue;
}
// feed fRgeDataList