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

merge master

Browse Source
This commit is contained in:
suter_a
2016-08-09 14:05:40 +02:00
parent a9744e273b 5c69c21f56
commit f93ef114e0
128 changed files with 1421 additions and 517 deletions
Download Patch File Download Diff File Expand all files Collapse all files

320
src/classes/PMsr2Data.cpp
View File

@@ -1643,13 +1643,23 @@ bool PMsr2Data::PrepareNewSortedInputFile(unsigned int tempRun) const
* - -2 if a fit has not converged (and the data is not appended to the output file)
*
* \param outfile name of the DB/ASCII output file
* \param paramList parameter list which shall be written to the output file
* \param db DB or plain ASCII output
* \param withHeader write output file header or not
* \param global global mode or not
* \param counter counter used within the global mode to determine how many runs have been processed already
*/
int PMsr2Data::WriteOutput(const string &outfile, bool db, unsigned int withHeader, bool global, unsigned int counter) const
int PMsr2Data::WriteOutput(const string &outfile, const vector<unsigned int>& paramList, bool db,
unsigned int withHeader, bool global, unsigned int counter) const
{
// make sure that the parameter number of the parameter list stays within proper bounds
for (unsigned int i=0; i<paramList.size(); i++) {
if (paramList[i] > fMsrHandler->GetMsrParamList()->size()) {
cerr << "msr2data: **ERROR** found parameter " << paramList[i] << " which is out of bound (>" << fMsrHandler->GetMsrParamList()->size() << ")." << endl;
return -1;
}
}
if (!to_lower_copy(outfile).compare("none")) {
fRunVectorIter++;
return PMUSR_SUCCESS;
@@ -2008,7 +2018,7 @@ int PMsr2Data::WriteOutput(const string &outfile, bool db, unsigned int withHead
if (fDataHandler) {
for (unsigned int i(0); i < dataParamLabels.size(); ++i) {
outFile << dataParamLabels[i] << endl;
outFile << dataParamLabels[i] << endl;
}
}
@@ -2022,22 +2032,26 @@ int PMsr2Data::WriteOutput(const string &outfile, bool db, unsigned int withHead
if (global) {
string tempName;
for (unsigned int i(0); i < fNumGlobalParam; ++i) {
outFile << (*msrParamList)[i].fName.Data() << endl;
if (InParameterList(i, paramList))
outFile << (*msrParamList)[i].fName.Data() << endl;
}
for (unsigned int i(0); i < fNumSpecParam; ++i) {
tempName = (*msrParamList)[fNumGlobalParam + fNumSpecParam*counter + i].fName.Data();
string::size_type loc = tempName.rfind(curRunNumber.str());
if (loc == tempName.length() - fRunNumberDigits) {
outFile << tempName.substr(0, loc) << endl;
} else {
cerr << endl << ">> msr2data: **ERROR** The run index of some parameter does not match the run number being processed!";
cerr << endl << ">> msr2data: **ERROR** The output will be flawed!";
cerr << endl;
if (InParameterList(fNumGlobalParam + fNumSpecParam*counter + i, paramList)) {
tempName = (*msrParamList)[fNumGlobalParam + fNumSpecParam*counter + i].fName.Data();
string::size_type loc = tempName.rfind(curRunNumber.str());
if (loc == tempName.length() - fRunNumberDigits) {
outFile << tempName.substr(0, loc) << endl;
} else {
cerr << endl << ">> msr2data: **ERROR** The run index of some parameter does not match the run number being processed!";
cerr << endl << ">> msr2data: **ERROR** The output will be flawed!";
cerr << endl;
}
}
}
} else {
for (unsigned int i(0); i < msrNoOfParams; ++i) {
outFile << (*msrParamList)[i].fName.Data() << endl;
if (InParameterList(i, paramList))
outFile << (*msrParamList)[i].fName.Data() << endl;
}
}
@@ -2073,22 +2087,26 @@ int PMsr2Data::WriteOutput(const string &outfile, bool db, unsigned int withHead
if (global) {
string tempName;
for (unsigned int i(0); i < fNumGlobalParam; ++i) {
outFile << " " << (*msrParamList)[i].fName.Data();
if (InParameterList(i, paramList))
outFile << " " << (*msrParamList)[i].fName.Data();
}
for (unsigned int i(0); i < fNumSpecParam; ++i) {
tempName = (*msrParamList)[fNumGlobalParam + fNumSpecParam*counter + i].fName.Data();
string::size_type loc = tempName.rfind(curRunNumber.str());
if (loc == tempName.length() - fRunNumberDigits) {
outFile << " " << tempName.substr(0, loc);
} else {
cerr << endl << ">> msr2data: **ERROR** The run index of some parameter does not match the run number being processed!";
cerr << endl << ">> msr2data: **ERROR** The output will be flawed!";
cerr << endl;
if (InParameterList(i, paramList)) {
tempName = (*msrParamList)[fNumGlobalParam + fNumSpecParam*counter + i].fName.Data();
string::size_type loc = tempName.rfind(curRunNumber.str());
if (loc == tempName.length() - fRunNumberDigits) {
outFile << " " << tempName.substr(0, loc);
} else {
cerr << endl << ">> msr2data: **ERROR** The run index of some parameter does not match the run number being processed!";
cerr << endl << ">> msr2data: **ERROR** The output will be flawed!";
cerr << endl;
}
}
}
} else {
for (unsigned int i(0); i < msrNoOfParams; ++i) {
outFile << " " << (*msrParamList)[i].fName.Data();
if (InParameterList(i, paramList))
outFile << " " << (*msrParamList)[i].fName.Data();
}
}
@@ -2132,61 +2150,67 @@ int PMsr2Data::WriteOutput(const string &outfile, bool db, unsigned int withHead
string tempName;
unsigned int idx;
for (unsigned int i(0); i < fNumGlobalParam; ++i) {
outFile << (*msrParamList)[i].fName.Data() << " = ";
if ((*msrParamList)[i].fPosErrorPresent) {
WriteValue(outFile, (*msrParamList)[i].fValue, (*msrParamList)[i].fPosError, outFile.width(), db);
outFile << ", ";
} else {
outFile << (*msrParamList)[i].fValue << ", ";
}
if (InParameterList(i, paramList)) {
outFile << (*msrParamList)[i].fName.Data() << " = ";
if ((*msrParamList)[i].fPosErrorPresent) {
WriteValue(outFile, (*msrParamList)[i].fValue, (*msrParamList)[i].fPosError, outFile.width(), db);
outFile << ", ";
} else {
outFile << (*msrParamList)[i].fValue << ", ";
}
if ((*msrParamList)[i].fPosErrorPresent)
outFile << (*msrParamList)[i].fPosError << ", ";
else
outFile << fabs((*msrParamList)[i].fStep) << ", ";
outFile << fabs((*msrParamList)[i].fStep) << ",\\" << endl;
if ((*msrParamList)[i].fPosErrorPresent)
outFile << (*msrParamList)[i].fPosError << ", ";
else
outFile << fabs((*msrParamList)[i].fStep) << ", ";
outFile << fabs((*msrParamList)[i].fStep) << ",\\" << endl;
}
}
for (unsigned int i(0); i < fNumSpecParam; ++i) {
idx = fNumGlobalParam + fNumSpecParam*counter + i;
tempName = (*msrParamList)[idx].fName.Data();
string::size_type loc = tempName.rfind(curRunNumber.str());
if (loc == tempName.length() - fRunNumberDigits) {
outFile << tempName.substr(0, loc) << " = ";
if ((*msrParamList)[idx].fPosErrorPresent) {
WriteValue(outFile, (*msrParamList)[idx].fValue, (*msrParamList)[idx].fPosError, outFile.width(), db);
outFile << ", ";
} else {
outFile << (*msrParamList)[idx].fValue << ", ";
}
if ((*msrParamList)[idx].fPosErrorPresent) {
WriteValue(outFile, (*msrParamList)[idx].fPosError, (*msrParamList)[idx].fPosError, outFile.width(), db);
outFile << ", ";
} else {
if (InParameterList(idx, paramList)) {
tempName = (*msrParamList)[idx].fName.Data();
string::size_type loc = tempName.rfind(curRunNumber.str());
if (loc == tempName.length() - fRunNumberDigits) {
outFile << tempName.substr(0, loc) << " = ";
if ((*msrParamList)[idx].fPosErrorPresent) {
WriteValue(outFile, (*msrParamList)[idx].fValue, (*msrParamList)[idx].fPosError, outFile.width(), db);
outFile << ", ";
} else {
outFile << (*msrParamList)[idx].fValue << ", ";
}
if ((*msrParamList)[idx].fPosErrorPresent) {
WriteValue(outFile, (*msrParamList)[idx].fPosError, (*msrParamList)[idx].fPosError, outFile.width(), db);
outFile << ", ";
} else {
WriteValue(outFile, (*msrParamList)[idx].fStep, (*msrParamList)[idx].fStep, outFile.width(), db);
outFile << ", ";
}
WriteValue(outFile, (*msrParamList)[idx].fStep, (*msrParamList)[idx].fStep, outFile.width(), db);
outFile << ", ";
outFile << ",\\" << endl;
}
WriteValue(outFile, (*msrParamList)[idx].fStep, (*msrParamList)[idx].fStep, outFile.width(), db);
outFile << ",\\" << endl;
}
}
} else {
for (unsigned int i(0); i < msrNoOfParams; ++i) {
outFile << (*msrParamList)[i].fName.Data() << " = ";
if ((*msrParamList)[i].fPosErrorPresent) {
WriteValue(outFile, (*msrParamList)[i].fValue, (*msrParamList)[i].fPosError, outFile.width(), db);
outFile << ", ";
} else {
outFile << (*msrParamList)[i].fValue << ", ";
}
if ((*msrParamList)[i].fPosErrorPresent) {
WriteValue(outFile, (*msrParamList)[i].fPosError, (*msrParamList)[i].fPosError, outFile.width(), db);
outFile << ", ";
} else {
if (InParameterList(i, paramList)) {
outFile << (*msrParamList)[i].fName.Data() << " = ";
if ((*msrParamList)[i].fPosErrorPresent) {
WriteValue(outFile, (*msrParamList)[i].fValue, (*msrParamList)[i].fPosError, outFile.width(), db);
outFile << ", ";
} else {
outFile << (*msrParamList)[i].fValue << ", ";
}
if ((*msrParamList)[i].fPosErrorPresent) {
WriteValue(outFile, (*msrParamList)[i].fPosError, (*msrParamList)[i].fPosError, outFile.width(), db);
outFile << ", ";
} else {
WriteValue(outFile, (*msrParamList)[i].fStep, (*msrParamList)[i].fStep, outFile.width(), db);
outFile << ", ";
}
WriteValue(outFile, (*msrParamList)[i].fStep, (*msrParamList)[i].fStep, outFile.width(), db);
outFile << ", ";
outFile << ",\\" << endl;
}
WriteValue(outFile, (*msrParamList)[i].fStep, (*msrParamList)[i].fStep, outFile.width(), db);
outFile << ",\\" << endl;
}
}
@@ -2221,30 +2245,38 @@ int PMsr2Data::WriteOutput(const string &outfile, bool db, unsigned int withHead
string s;
if (global) {
for (unsigned int i(0); i < fNumGlobalParam; ++i) {
s = (*msrParamList)[i].fName.Data();
length = s.length();
if (length > maxlength)
maxlength = length;
}
for (unsigned int i(0); i < fNumSpecParam; ++i) {
s = (*msrParamList)[fNumGlobalParam + fNumSpecParam*counter + i].fName.Data();
string::size_type loc = s.rfind(curRunNumber.str());
if (loc == s.length() - fRunNumberDigits) {
length = s.length() - fRunNumberDigits;
if (InParameterList(i, paramList)) {
s = (*msrParamList)[i].fName.Data();
length = s.length();
if (length > maxlength)
maxlength = length;
} else {
cerr << endl << ">> msr2data: **ERROR** The run index of some parameter does not match the run number being processed!";
cerr << endl << ">> msr2data: **ERROR** The output will be flawed!";
cerr << endl;
}
}
unsigned int idx;
for (unsigned int i(0); i < fNumSpecParam; ++i) {
idx = fNumGlobalParam + fNumSpecParam*counter + i;
if (InParameterList(idx, paramList)) {
s = (*msrParamList)[idx].fName.Data();
string::size_type loc = s.rfind(curRunNumber.str());
if (loc == s.length() - fRunNumberDigits) {
length = s.length() - fRunNumberDigits;
if (length > maxlength)
maxlength = length;
} else {
cerr << endl << ">> msr2data: **ERROR** The run index of some parameter does not match the run number being processed!";
cerr << endl << ">> msr2data: **ERROR** The output will be flawed!";
cerr << endl;
}
}
}
} else {
for (unsigned int i(0); i < msrNoOfParams; ++i) {
s = (*msrParamList)[i].fName.Data();
length = s.length();
if (length > maxlength)
maxlength = length;
if (InParameterList(i, paramList)) {
s = (*msrParamList)[i].fName.Data();
length = s.length();
if (length > maxlength)
maxlength = length;
}
}
}
if (maxlength < 13)
@@ -2274,31 +2306,39 @@ int PMsr2Data::WriteOutput(const string &outfile, bool db, unsigned int withHead
// in the GLOBAL mode write only global parameters and those which belong to the actual run - in the NORMAL mode write all parameters
if (global) {
for (unsigned int i(0); i < fNumGlobalParam; ++i) {
s = (*msrParamList)[i].fName.Data();
outFile << setw(maxlength) << left << s \
<< setw(maxlength + 6) << left << s + "PosErr" \
<< setw(maxlength + 6) << left << s + "NegErr";
}
for (unsigned int i(0); i < fNumSpecParam; ++i) {
s = (*msrParamList)[fNumGlobalParam + fNumSpecParam*counter + i].fName.Data();
string::size_type loc = s.rfind(curRunNumber.str());
if (loc == s.length() - fRunNumberDigits) {
s = s.substr(0, loc);
if (InParameterList(i, paramList)) {
s = (*msrParamList)[i].fName.Data();
outFile << setw(maxlength) << left << s \
<< setw(maxlength + 6) << left << s + "PosErr" \
<< setw(maxlength + 6) << left << s + "NegErr";
} else {
cerr << endl << ">> msr2data: **ERROR** The run index of some parameter does not match the run number being processed!";
cerr << endl << ">> msr2data: **ERROR** The output will be flawed!";
cerr << endl;
}
}
unsigned int idx;
for (unsigned int i(0); i < fNumSpecParam; ++i) {
idx = fNumGlobalParam + fNumSpecParam*counter + i;
if (InParameterList(idx, paramList)) {
s = (*msrParamList)[idx].fName.Data();
string::size_type loc = s.rfind(curRunNumber.str());
if (loc == s.length() - fRunNumberDigits) {
s = s.substr(0, loc);
outFile << setw(maxlength) << left << s \
<< setw(maxlength + 6) << left << s + "PosErr" \
<< setw(maxlength + 6) << left << s + "NegErr";
} else {
cerr << endl << ">> msr2data: **ERROR** The run index of some parameter does not match the run number being processed!";
cerr << endl << ">> msr2data: **ERROR** The output will be flawed!";
cerr << endl;
}
}
}
} else {
for (unsigned int i(0); i < msrNoOfParams; ++i) {
s = (*msrParamList)[i].fName.Data();
outFile << setw(maxlength) << left << s \
<< setw(maxlength + 6) << left << s + "PosErr" \
<< setw(maxlength + 6) << left << s + "NegErr";
if (InParameterList(i, paramList)) {
s = (*msrParamList)[i].fName.Data();
outFile << setw(maxlength) << left << s \
<< setw(maxlength + 6) << left << s + "PosErr" \
<< setw(maxlength + 6) << left << s + "NegErr";
}
}
}
s.clear();
@@ -2340,46 +2380,52 @@ int PMsr2Data::WriteOutput(const string &outfile, bool db, unsigned int withHead
// in the GLOBAL mode write only global parameters and those which belong to the actual run - in the NORMAL mode write all parameters
if (global) {
for (unsigned int i(0); i < fNumGlobalParam; ++i) {
if ((*msrParamList)[i].fPosErrorPresent)
WriteValue(outFile, (*msrParamList)[i].fValue, (*msrParamList)[i].fPosError, maxlength, db);
else
WriteValue(outFile, (*msrParamList)[i].fValue, maxlength);
if (InParameterList(i, paramList)) {
if ((*msrParamList)[i].fPosErrorPresent)
WriteValue(outFile, (*msrParamList)[i].fValue, (*msrParamList)[i].fPosError, maxlength, db);
else
WriteValue(outFile, (*msrParamList)[i].fValue, maxlength);
if ((*msrParamList)[i].fPosErrorPresent)
WriteValue(outFile, (*msrParamList)[i].fPosError, (*msrParamList)[i].fPosError, maxlength, db);
else
WriteValue(outFile, fabs((*msrParamList)[i].fStep), (*msrParamList)[i].fStep, maxlength, db);
if ((*msrParamList)[i].fPosErrorPresent)
WriteValue(outFile, (*msrParamList)[i].fPosError, (*msrParamList)[i].fPosError, maxlength, db);
else
WriteValue(outFile, fabs((*msrParamList)[i].fStep), (*msrParamList)[i].fStep, maxlength, db);
WriteValue(outFile, fabs((*msrParamList)[i].fStep), (*msrParamList)[i].fStep, maxlength, db);
}
}
unsigned int idx;
for (unsigned int i(0); i < fNumSpecParam; ++i) {
idx = fNumGlobalParam + fNumSpecParam*counter + i;
if ((*msrParamList)[idx].fPosErrorPresent)
WriteValue(outFile, (*msrParamList)[idx].fValue, (*msrParamList)[idx].fPosError, maxlength, db);
else
WriteValue(outFile, (*msrParamList)[idx].fValue, maxlength);
if (InParameterList(idx, paramList)) {
if ((*msrParamList)[idx].fPosErrorPresent)
WriteValue(outFile, (*msrParamList)[idx].fValue, (*msrParamList)[idx].fPosError, maxlength, db);
else
WriteValue(outFile, (*msrParamList)[idx].fValue, maxlength);
if ((*msrParamList)[idx].fPosErrorPresent)
WriteValue(outFile, (*msrParamList)[idx].fPosError, (*msrParamList)[idx].fPosError, maxlength, db);
else
WriteValue(outFile, (*msrParamList)[idx].fStep, (*msrParamList)[idx].fStep, maxlength, db);
if ((*msrParamList)[idx].fPosErrorPresent)
WriteValue(outFile, (*msrParamList)[idx].fPosError, (*msrParamList)[idx].fPosError, maxlength, db);
else
WriteValue(outFile, (*msrParamList)[idx].fStep, (*msrParamList)[idx].fStep, maxlength, db);
WriteValue(outFile, (*msrParamList)[idx].fStep, (*msrParamList)[idx].fStep, maxlength, db);
}
}
} else {
for (unsigned int i(0); i < msrNoOfParams; ++i) {
if ((*msrParamList)[i].fPosErrorPresent)
WriteValue(outFile, (*msrParamList)[i].fValue, (*msrParamList)[i].fPosError, maxlength, db);
else
WriteValue(outFile, (*msrParamList)[i].fValue, maxlength);
if (InParameterList(i, paramList)) {
if ((*msrParamList)[i].fPosErrorPresent)
WriteValue(outFile, (*msrParamList)[i].fValue, (*msrParamList)[i].fPosError, maxlength, db);
else
WriteValue(outFile, (*msrParamList)[i].fValue, maxlength);
if ((*msrParamList)[i].fPosErrorPresent)
WriteValue(outFile, (*msrParamList)[i].fPosError, (*msrParamList)[i].fPosError, maxlength, db);
else
WriteValue(outFile, fabs((*msrParamList)[i].fStep), fabs((*msrParamList)[i].fStep), maxlength, db);
if ((*msrParamList)[i].fPosErrorPresent)
WriteValue(outFile, (*msrParamList)[i].fPosError, (*msrParamList)[i].fPosError, maxlength, db);
else
WriteValue(outFile, fabs((*msrParamList)[i].fStep), fabs((*msrParamList)[i].fStep), maxlength, db);
WriteValue(outFile, fabs((*msrParamList)[i].fStep), fabs((*msrParamList)[i].fStep), maxlength, db);
}
}
}
@@ -2505,6 +2551,24 @@ int PMsr2Data::GetFirstSignificantDigit(const double &value) const
return prec+1;
}
//-------------------------------------------------------------
/**
*
*/
bool PMsr2Data::InParameterList(const unsigned int &paramValue, const vector<unsigned int> &paramList) const
{
// if paramList.size() == 0, i.e. use ALL parameters
if (paramList.size() == 0)
return true;
for (unsigned int i=0; i<paramList.size(); i++) {
if (paramValue+1 == paramList[i])
return true;
}
return false;
}
//-------------------------------------------------------------------------------------------------------
// end
//-------------------------------------------------------------------------------------------------------

30
src/classes/PRunAsymmetryRRF.cpp
View File

@@ -187,15 +187,7 @@ Double_t PRunAsymmetryRRF::CalcChiSquare(const std::vector<Double_t>& par)
// calculate chi square
Double_t time(1.0);
Int_t i, N(static_cast<Int_t>(fData.GetValue()->size()));
// In order not to have an IF in the next loop, determine the start and end bins for the fit range now
Int_t startTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (startTimeBin < 0)
startTimeBin = 0;
Int_t endTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (endTimeBin > N)
endTimeBin = N;
Int_t i;
// Calculate the theory function once to ensure one function evaluation for the current set of parameters.
// This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
@@ -204,12 +196,12 @@ Double_t PRunAsymmetryRRF::CalcChiSquare(const std::vector<Double_t>& par)
asymFcnValue = fTheory->Func(time, par, fFuncValues);
#ifdef HAVE_GOMP
Int_t chunk = (endTimeBin - startTimeBin)/omp_get_num_procs();
Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
if (chunk < 10)
chunk = 10;
#pragma omp parallel for default(shared) private(i,time,diff,asymFcnValue,a,b,f) schedule(dynamic,chunk) reduction(+:chisq)
#endif
for (i=startTimeBin; i<endTimeBin; ++i) {
for (i=fStartTimeBin; i<fEndTimeBin; ++i) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
switch (fAlphaBetaTag) {
case 1: // alpha == 1, beta == 1
@@ -387,15 +379,15 @@ void PRunAsymmetryRRF::SetFitRangeBin(const TString fitRange)
void PRunAsymmetryRRF::CalcNoOfFitBins()
{
// In order not having to loop over all bins and to stay consistent with the chisq method, calculate the start and end bins explicitly
Int_t startTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (startTimeBin < 0)
startTimeBin = 0;
Int_t endTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (endTimeBin > static_cast<Int_t>(fData.GetValue()->size()))
endTimeBin = fData.GetValue()->size();
fStartTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (fStartTimeBin < 0)
fStartTimeBin = 0;
fEndTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (fEndTimeBin > static_cast<Int_t>(fData.GetValue()->size()))
fEndTimeBin = fData.GetValue()->size();
if (endTimeBin > startTimeBin)
fNoOfFitBins = endTimeBin - startTimeBin;
if (fEndTimeBin > fStartTimeBin)
fNoOfFitBins = fEndTimeBin - fStartTimeBin;
else
fNoOfFitBins = 0;
}

76
src/classes/PRunSingleHistoRRF.cpp
View File

@@ -58,13 +58,16 @@ using namespace std;
PRunSingleHistoRRF::PRunSingleHistoRRF() : PRunBase()
{
fNoOfFitBins = 0;
fBackground = 0;
fBackground = 0.0;
fBkgErr = 1.0;
fRRFPacking = -1;
// the 2 following variables are need in case fit range is given in bins, and since
// the fit range can be changed in the command block, these variables need to be accessible
fGoodBins[0] = -1;
fGoodBins[1] = -1;
fN0EstimateEndTime = 1.0; // end time in (us) over which N0 is estimated.
}
//--------------------------------------------------------------------------
@@ -112,6 +115,8 @@ PRunSingleHistoRRF::PRunSingleHistoRRF(PMsrHandler *msrInfo, PRunDataHandler *ra
fGoodBins[0] = -1;
fGoodBins[1] = -1;
fN0EstimateEndTime = 1.0; // end time in (us) over which N0 is estimated.
if (!PrepareData()) {
cerr << endl << ">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: Couldn't prepare data for fitting!";
cerr << endl << ">> This is very bad :-(, will quit ...";
@@ -155,15 +160,7 @@ Double_t PRunSingleHistoRRF::CalcChiSquare(const std::vector<Double_t>& par)
// calculate chi square
Double_t time(1.0);
Int_t i, N(static_cast<Int_t>(fData.GetValue()->size()));
// In order not to have an IF in the next loop, determine the start and end bins for the fit range now
Int_t startTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (startTimeBin < 0)
startTimeBin = 0;
Int_t endTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (endTimeBin > N)
endTimeBin = N;
Int_t i;
// Calculate the theory function once to ensure one function evaluation for the current set of parameters.
// This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
@@ -172,12 +169,12 @@ Double_t PRunSingleHistoRRF::CalcChiSquare(const std::vector<Double_t>& par)
time = fTheory->Func(time, par, fFuncValues);
#ifdef HAVE_GOMP
Int_t chunk = (endTimeBin - startTimeBin)/omp_get_num_procs();
Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
if (chunk < 10)
chunk = 10;
#pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq)
#endif
for (i=startTimeBin; i<endTimeBin; ++i) {
for (i=fStartTimeBin; i<fEndTimeBin; ++i) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
diff = fData.GetValue()->at(i) - fTheory->Func(time, par, fFuncValues);
chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i));
@@ -211,15 +208,7 @@ Double_t PRunSingleHistoRRF::CalcChiSquareExpected(const std::vector<Double_t>&
// calculate chi square
Double_t time(1.0);
Int_t i, N(static_cast<Int_t>(fData.GetValue()->size()));
// In order not to have an IF in the next loop, determine the start and end bins for the fit range now
Int_t startTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (startTimeBin < 0)
startTimeBin = 0;
Int_t endTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (endTimeBin > N)
endTimeBin = N;
Int_t i;
// Calculate the theory function once to ensure one function evaluation for the current set of parameters.
// This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once
@@ -228,12 +217,12 @@ Double_t PRunSingleHistoRRF::CalcChiSquareExpected(const std::vector<Double_t>&
time = fTheory->Func(time, par, fFuncValues);
#ifdef HAVE_GOMP
Int_t chunk = (endTimeBin - startTimeBin)/omp_get_num_procs();
Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs();
if (chunk < 10)
chunk = 10;
#pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq)
#endif
for (i=startTimeBin; i < endTimeBin; ++i) {
for (i=fStartTimeBin; i < fEndTimeBin; ++i) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
theo = fTheory->Func(time, par, fFuncValues);
diff = fData.GetValue()->at(i) - theo;
@@ -408,15 +397,15 @@ void PRunSingleHistoRRF::SetFitRangeBin(const TString fitRange)
void PRunSingleHistoRRF::CalcNoOfFitBins()
{
// In order not having to loop over all bins and to stay consistent with the chisq method, calculate the start and end bins explicitly
Int_t startTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (startTimeBin < 0)
startTimeBin = 0;
Int_t endTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (endTimeBin > static_cast<Int_t>(fData.GetValue()->size()))
endTimeBin = fData.GetValue()->size();
fStartTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
if (fStartTimeBin < 0)
fStartTimeBin = 0;
fEndTimeBin = static_cast<Int_t>(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1;
if (fEndTimeBin > static_cast<Int_t>(fData.GetValue()->size()))
fEndTimeBin = fData.GetValue()->size();
if (endTimeBin > startTimeBin)
fNoOfFitBins = endTimeBin - startTimeBin;
if (fEndTimeBin > fStartTimeBin)
fNoOfFitBins = fEndTimeBin - fStartTimeBin;
else
fNoOfFitBins = 0;
}
@@ -604,6 +593,9 @@ Bool_t PRunSingleHistoRRF::PrepareFitData(PRawRunData* runData, const UInt_t his
if (!EstimateBkg(histoNo))
return false;
}
// subtract background from fForward
for (UInt_t i=0; i<fForward.size(); i++)
fForward[i] -= fBackground;
} else { // fixed background given
for (UInt_t i=0; i<fForward.size(); i++) {
fForward[i] -= fRunInfo->GetBkgFix(0);
@@ -624,7 +616,7 @@ Bool_t PRunSingleHistoRRF::PrepareFitData(PRawRunData* runData, const UInt_t his
exp_t_tau = exp(time_tau);
fForward[i] *= exp_t_tau;
fM.push_back(fForward[i]); // i.e. M(t) = [N(t)-Nbkg] exp(+t/tau); needed to estimate N0 later on
fMerr.push_back(exp_t_tau*sqrt(rawNt[i]-fBackground));
fMerr.push_back(exp_t_tau*sqrt(rawNt[i]+fBkgErr*fBkgErr));
}
// calculate weights
@@ -632,7 +624,7 @@ Bool_t PRunSingleHistoRRF::PrepareFitData(PRawRunData* runData, const UInt_t his
if (fMerr[i] > 0.0)
fW.push_back(1.0/(fMerr[i]*fMerr[i]));
else
fW.push_back(0.0);
fW.push_back(1.0);
}
// now fForward = exp(+t/tau) [N(t)-Nbkg] = M(t)
@@ -1071,6 +1063,9 @@ Double_t PRunSingleHistoRRF::GetMainFrequency(PDoubleVector &data)
if (power->GetBinContent(i)>power->GetBinContent(i+1))
continue;
}
// ignore everything below 10 MHz
if (power->GetBinCenter(i) < 10.0)
continue;
// check for maximum
if (power->GetBinContent(i) > maxFreqVal) {
maxFreqVal = power->GetBinContent(i);
@@ -1101,15 +1096,17 @@ Double_t PRunSingleHistoRRF::EstimateN0(Double_t &errN0, Double_t freqMax)
{
// endBin is estimated such that the number of full cycles (according to the maximum frequency of the data)
// is approximately the time fN0EstimateEndTime.
Int_t endBin = (Int_t)round(fN0EstimateEndTime / fTimeResolution * ceil(freqMax)/freqMax);
Int_t endBin = (Int_t)round(ceil(fN0EstimateEndTime*freqMax/TMath::TwoPi()) * (TMath::TwoPi()/freqMax) / fTimeResolution);
Double_t n0 = 0.0;
Double_t wN = 0.0;
for (Int_t i=0; i<endBin; i++) {
n0 += fW[i]*fM[i];
// n0 += fW[i]*fM[i];
n0 += fM[i];
wN += fW[i];
}
n0 /= wN;
// n0 /= wN;
n0 /= endBin;
errN0 = 0.0;
for (Int_t i=0; i<endBin; i++) {
@@ -1197,7 +1194,12 @@ Bool_t PRunSingleHistoRRF::EstimateBkg(UInt_t histoNo)
fBackground = bkg; // keep background (per bin)
cout << endl << "info> fBackground=" << fBackground << endl;
bkg = 0.0;
for (UInt_t i=start; i<end; i++)
bkg += pow(fForward[i]-fBackground, 2.0);
fBkgErr = sqrt(bkg/(static_cast<Double_t>(end - start)));
cout << endl << "info> fBackground=" << fBackground << "(" << fBkgErr << ")" << endl;
fRunInfo->SetBkgEstimated(fBackground, 0);

46
src/classes/PStartupHandler.cpp
View File

@@ -97,7 +97,6 @@ PStartupHandler::PStartupHandler()
Char_t *home=0;
Char_t musrpath[128];
Char_t startup_path_name[128];
Bool_t found = false;
strncpy(musrpath, "", sizeof(musrpath));
@@ -106,32 +105,39 @@ PStartupHandler::PStartupHandler()
if (StartupFileExists(startup_path_name)) {
fStartupFileFound = true;
fStartupFilePath = TString(startup_path_name);
} else { // startup file is not found in the current directory
}
if (!fStartupFileFound) { // startup file not found in the current directory
// check if the startup file is found under $HOME/.musrfit
home = getenv("HOME");
if (home != 0) {
sprintf(musrpath, "%s/.musrfit", home);
found = true;
}
pmusrpath = getenv("MUSRFITPATH");
if (!found) {
// check if the MUSRFITPATH system variable is set
if (pmusrpath != 0) {
if (strcmp(pmusrpath, "")) { // MUSRFITPATH variable set but empty
found = true;
}
sprintf(startup_path_name, "%s/.musrfit/musrfit_startup.xml", home);
if (StartupFileExists(startup_path_name)) {
fStartupFilePath = TString(startup_path_name);
fStartupFileFound = true;
}
}
if (!found) { // MUSRFITPATH not set or empty, will try default one
home = getenv("ROOTSYS");
}
if (!fStartupFileFound) { // startup file not found in $HOME/.musrfit
// check if the MUSRFITPATH system variable is set
pmusrpath = getenv("MUSRFITPATH");
if (pmusrpath != 0) {
sprintf(startup_path_name, "%s/musrfit_startup.xml", pmusrpath);
if (StartupFileExists(startup_path_name)) {
fStartupFilePath = TString(startup_path_name);
fStartupFileFound = true;
}
}
}
if (!fStartupFileFound) { // MUSRFITPATH not set or empty, will try $ROOTSYS/bin
home = getenv("ROOTSYS");
if (home != 0) {
sprintf(musrpath, "%s/bin", home);
cerr << endl << "**WARNING** MUSRFITPATH environment variable not set will try " << musrpath << endl;
}
sprintf(startup_path_name, "%s/musrfit_startup.xml", musrpath);
fStartupFilePath = TString(startup_path_name);
if (StartupFileExists(startup_path_name)) {
fStartupFileFound = true;
sprintf(startup_path_name, "%s/musrfit_startup.xml", musrpath);
if (StartupFileExists(startup_path_name)) {
fStartupFilePath = TString(startup_path_name);
fStartupFileFound = true;
}
}
}
}