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

added FIT_RANGE RESET | start end | s1 e1 s2 e2 .. sN eN command to the COMMAND block

Browse Source
This commit is contained in:
nemu
2010-10-20 09:24:12 +00:00
parent 616740f299
commit 5cd4283d30
16 changed files with 416 additions and 66 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
ChangeLog
View File

@ -4,8 +4,9 @@
# $Id$ # $Id$
#--------------------------------------------------------------------- #---------------------------------------------------------------------
changes since 0.6.0 musrfit 0.7.0 - changes since 0.6.0
=================== ===================================
NEW added FIT_RANGE RESET | start end | s1 e1 s2 e2 .. sN eN command to the COMMAND block
NEW added FIX/RELEASE/RESTORE minuit2 command to the COMMAND block NEW added FIX/RELEASE/RESTORE minuit2 command to the COMMAND block
NEW implemented more checks on the integrity of the msr-file NEW implemented more checks on the integrity of the msr-file
NEW [MUSRGUI/MUSREDIT] msr2data default option flags in musredit_startup.xml NEW [MUSRGUI/MUSREDIT] msr2data default option flags in musredit_startup.xml

4
configure.ac
View File

@ -1,5 +1,5 @@
AC_PREREQ(2.59) AC_PREREQ(2.59)
AC_INIT(musrfit, 0.6.0, andreas.suter@psi.ch) AC_INIT(musrfit, 0.7.0, andreas.suter@psi.ch)
AC_CONFIG_AUX_DIR(admin) AC_CONFIG_AUX_DIR(admin)
AC_CANONICAL_HOST AC_CANONICAL_HOST
#AC_MSG_RESULT([${host} ${host_cpu} ${host_vendor} ${host_os}]) #AC_MSG_RESULT([${host} ${host_cpu} ${host_vendor} ${host_os}])
@ -28,7 +28,7 @@ dnl -----------------------------------------------
#release versioning #release versioning
MUSR_MAJOR_VERSION=0 MUSR_MAJOR_VERSION=0
MUSR_MINOR_VERSION=6 MUSR_MINOR_VERSION=7
MUSR_MICRO_VERSION=0 MUSR_MICRO_VERSION=0
#release versioning #release versioning

225
src/classes/PFitter.cpp
View File

@ -69,7 +69,7 @@ using namespace std;
* \param chisq_only flag: true=calculate chisq only (no fitting) * \param chisq_only flag: true=calculate chisq only (no fitting)
*/ */
PFitter::PFitter(PMsrHandler *runInfo, PRunListCollection *runListCollection, Bool_t chisq_only) : PFitter::PFitter(PMsrHandler *runInfo, PRunListCollection *runListCollection, Bool_t chisq_only) :
fChisqOnly(chisq_only), fRunInfo(runInfo) fChisqOnly(chisq_only), fRunInfo(runInfo), fRunListCollection(runListCollection)
{ {
// initialize variables // initialize variables
fIsScanOnly = true; fIsScanOnly = true;
@ -92,6 +92,15 @@ PFitter::PFitter(PMsrHandler *runInfo, PRunListCollection *runListCollection, Bo
fScanLow = 0.0; // minuit2 default, i.e. 2 std deviations fScanLow = 0.0; // minuit2 default, i.e. 2 std deviations
fScanHigh = 0.0; // minuit2 default, i.e. 2 std deviations fScanHigh = 0.0; // minuit2 default, i.e. 2 std deviations
// keep all the fit ranges in case RANGE command is present
PMsrRunList *runs = fRunInfo->GetMsrRunList();
PDoublePair range;
for (UInt_t i=0; i<runs->size(); i++) {
range.first = (*runs)[i].GetFitRange(0);
range.second = (*runs)[i].GetFitRange(1);
fOriginalFitRange.push_back(range);
}
// check msr minuit commands // check msr minuit commands
if (!CheckCommands()) { if (!CheckCommands()) {
return; return;
@ -183,6 +192,9 @@ Bool_t PFitter::DoFit()
case PMN_CONTOURS: case PMN_CONTOURS:
status = ExecuteContours(); status = ExecuteContours();
break; break;
case PMN_FIT_RANGE:
status = ExecuteFitRange(fCmdList[i].second);
break;
case PMN_FIX: case PMN_FIX:
status = ExecuteFix(fCmdList[i].second); status = ExecuteFix(fCmdList[i].second);
break; break;
@ -318,6 +330,10 @@ Bool_t PFitter::CheckCommands()
cerr << endl << ">> command syntax for CONTOURS is: CONTOURS parameter-X parameter-Y [# of points]"; cerr << endl << ">> command syntax for CONTOURS is: CONTOURS parameter-X parameter-Y [# of points]";
cerr << endl; cerr << endl;
fIsValid = false; fIsValid = false;
if (tokens) {
delete tokens;
tokens = 0;
}
break; break;
} }
ival = str.Atoi(); ival = str.Atoi();
@ -329,6 +345,10 @@ Bool_t PFitter::CheckCommands()
cerr << endl << ">> command syntax for CONTOURS is: CONTOURS parameter-X parameter-Y [# of points]"; cerr << endl << ">> command syntax for CONTOURS is: CONTOURS parameter-X parameter-Y [# of points]";
cerr << endl; cerr << endl;
fIsValid = false; fIsValid = false;
if (tokens) {
delete tokens;
tokens = 0;
}
break; break;
} }
// keep parameter // keep parameter
@ -345,6 +365,10 @@ Bool_t PFitter::CheckCommands()
cerr << endl << ">> command syntax for CONTOURS is: CONTOURS parameter-X parameter-Y [# of points]"; cerr << endl << ">> command syntax for CONTOURS is: CONTOURS parameter-X parameter-Y [# of points]";
cerr << endl; cerr << endl;
fIsValid = false; fIsValid = false;
if (tokens) {
delete tokens;
tokens = 0;
}
break; break;
} }
ival = str.Atoi(); ival = str.Atoi();
@ -355,6 +379,10 @@ Bool_t PFitter::CheckCommands()
cerr << endl << ">> command syntax for CONTOURS is: CONTOURS parameter-X parameter-Y [# of points]"; cerr << endl << ">> command syntax for CONTOURS is: CONTOURS parameter-X parameter-Y [# of points]";
cerr << endl; cerr << endl;
fIsValid = false; fIsValid = false;
if (tokens) {
delete tokens;
tokens = 0;
}
break; break;
} }
fScanNoPoints = ival; fScanNoPoints = ival;
@ -369,6 +397,90 @@ Bool_t PFitter::CheckCommands()
cmd.first = PMN_EIGEN; cmd.first = PMN_EIGEN;
cmd.second = cmdLineNo; cmd.second = cmdLineNo;
fCmdList.push_back(cmd); fCmdList.push_back(cmd);
} else if (it->fLine.Contains("FIT_RANGE", TString::kIgnoreCase)) {
// check the 3 options: FIT_RANGE RESET, FIT_RANGE start end, FIT_RANGE start1 end1 start2 end2 ... startN endN
TObjArray *tokens = 0;
TObjString *ostr;
TString str;
tokens = it->fLine.Tokenize(", \t");
if (tokens->GetEntries() == 2) {
ostr = dynamic_cast<TObjString*>(tokens->At(1));
str = ostr->GetString();
if (str.Contains("RESET"), TString::kIgnoreCase) {
cmd.first = PMN_FIT_RANGE;
cmd.second = cmdLineNo;
fCmdList.push_back(cmd);
} else {
cerr << endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
cerr << endl << ">> " << it->fLine.Data();
cerr << endl << ">> Syntax: FIT_RANGE RESET | FIT_RANGE start end | FIT_RANGE s1 e1 s2 e2 .. sN eN,";
cerr << endl << ">> with N the number of runs in the msr-file." << endl;
cerr << endl << ">> Found " << str.Data() << ", instead of RESET" << endl;
fIsValid = false;
if (tokens) {
delete tokens;
tokens = 0;
}
break;
}
} else if ((tokens->GetEntries() > 1) && (static_cast<UInt_t>(tokens->GetEntries()) % 2) == 1) {
if ((tokens->GetEntries() > 3) && ((static_cast<UInt_t>(tokens->GetEntries())-1)) != 2*fRunInfo->GetMsrRunList()->size()) {
cerr << endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
cerr << endl << ">> " << it->fLine.Data();
cerr << endl << ">> Syntax: FIT_RANGE RESET | FIT_RANGE start end | FIT_RANGE s1 e1 s2 e2 .. sN eN,";
cerr << endl << ">> with N the number of runs in the msr-file.";
cerr << endl << ">> Found N=" << (tokens->GetEntries()-1)/2 << ", # runs in msr-file=" << fRunInfo->GetMsrRunList()->size() << endl;
fIsValid = false;
if (tokens) {
delete tokens;
tokens = 0;
}
break;
} else {
// check that all range entries are numbers
Int_t n=1;
do {
ostr = dynamic_cast<TObjString*>(tokens->At(n));
str = ostr->GetString();
} while ((++n < tokens->GetEntries()) && str.IsFloat());
if (str.IsFloat()) { // everything is fine, last string was a floating point number
cmd.first = PMN_FIT_RANGE;
cmd.second = cmdLineNo;
fCmdList.push_back(cmd);
} else {
cerr << endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
cerr << endl << ">> " << it->fLine.Data();
cerr << endl << ">> Syntax: FIT_RANGE RESET | FIT_RANGE <start> <end> | FIT_RANGE <s1> <e1> <s2> <e2> .. <sN> <eN>,";
cerr << endl << ">> with N the number of runs in the msr-file.";
cerr << endl << ">> Found token '" << str.Data() << "', which is not a floating point number." << endl;
fIsValid = false;
if (tokens) {
delete tokens;
tokens = 0;
}
break;
}
}
} else {
cerr << endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
cerr << endl << ">> " << it->fLine.Data();
cerr << endl << ">> Syntax: FIT_RANGE RESET | FIT_RANGE start end | FIT_RANGE s1 e1 s2 e2 .. sN eN,";
cerr << endl << ">> with N the number of runs in the msr-file." << endl;
fIsValid = false;
if (tokens) {
delete tokens;
tokens = 0;
}
break;
}
if (tokens) {
delete tokens;
tokens = 0;
}
} else if (it->fLine.Contains("FIX", TString::kIgnoreCase)) { } else if (it->fLine.Contains("FIX", TString::kIgnoreCase)) {
// check if the given set of parameters (number or names) is present // check if the given set of parameters (number or names) is present
TObjArray *tokens = 0; TObjArray *tokens = 0;
@ -391,6 +503,10 @@ Bool_t PFitter::CheckCommands()
cerr << endl << ">> Parameter " << ival << " is out of the Parameter Range [1," << fParams.size() << "]"; cerr << endl << ">> Parameter " << ival << " is out of the Parameter Range [1," << fParams.size() << "]";
cerr << endl; cerr << endl;
fIsValid = false; fIsValid = false;
if (tokens) {
delete tokens;
tokens = 0;
}
break; break;
} }
} else { // token might be a parameter name } else { // token might be a parameter name
@ -408,6 +524,10 @@ Bool_t PFitter::CheckCommands()
cerr << endl << ">> Parameter '" << str.Data() << "' is NOT present as a parameter name"; cerr << endl << ">> Parameter '" << str.Data() << "' is NOT present as a parameter name";
cerr << endl; cerr << endl;
fIsValid = false; fIsValid = false;
if (tokens) {
delete tokens;
tokens = 0;
}
break; break;
} }
} }
@ -472,6 +592,10 @@ Bool_t PFitter::CheckCommands()
cerr << endl << ">> Parameter " << ival << " is out of the Parameter Range [1," << fParams.size() << "]"; cerr << endl << ">> Parameter " << ival << " is out of the Parameter Range [1," << fParams.size() << "]";
cerr << endl; cerr << endl;
fIsValid = false; fIsValid = false;
if (tokens) {
delete tokens;
tokens = 0;
}
break; break;
} }
} else { // token might be a parameter name } else { // token might be a parameter name
@ -489,6 +613,10 @@ Bool_t PFitter::CheckCommands()
cerr << endl << ">> Parameter '" << str.Data() << "' is NOT present as a parameter name"; cerr << endl << ">> Parameter '" << str.Data() << "' is NOT present as a parameter name";
cerr << endl; cerr << endl;
fIsValid = false; fIsValid = false;
if (tokens) {
delete tokens;
tokens = 0;
}
break; break;
} }
} }
@ -533,6 +661,10 @@ Bool_t PFitter::CheckCommands()
cerr << endl << ">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]"; cerr << endl << ">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]";
cerr << endl; cerr << endl;
fIsValid = false; fIsValid = false;
if (tokens) {
delete tokens;
tokens = 0;
}
break; break;
} }
ival = str.Atoi(); ival = str.Atoi();
@ -544,6 +676,10 @@ Bool_t PFitter::CheckCommands()
cerr << endl << ">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]"; cerr << endl << ">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]";
cerr << endl; cerr << endl;
fIsValid = false; fIsValid = false;
if (tokens) {
delete tokens;
tokens = 0;
}
break; break;
} }
// keep parameter // keep parameter
@ -560,6 +696,10 @@ Bool_t PFitter::CheckCommands()
cerr << endl << ">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]"; cerr << endl << ">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]";
cerr << endl; cerr << endl;
fIsValid = false; fIsValid = false;
if (tokens) {
delete tokens;
tokens = 0;
}
break; break;
} }
ival = str.Atoi(); ival = str.Atoi();
@ -570,6 +710,10 @@ Bool_t PFitter::CheckCommands()
cerr << endl << ">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]"; cerr << endl << ">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]";
cerr << endl; cerr << endl;
fIsValid = false; fIsValid = false;
if (tokens) {
delete tokens;
tokens = 0;
}
break; break;
} }
fScanNoPoints = ival; fScanNoPoints = ival;
@ -584,6 +728,10 @@ Bool_t PFitter::CheckCommands()
cerr << endl << ">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]"; cerr << endl << ">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]";
cerr << endl; cerr << endl;
fIsValid = false; fIsValid = false;
if (tokens) {
delete tokens;
tokens = 0;
}
break; break;
} }
fScanLow = str.Atof(); fScanLow = str.Atof();
@ -598,6 +746,10 @@ Bool_t PFitter::CheckCommands()
cerr << endl << ">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]"; cerr << endl << ">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]";
cerr << endl; cerr << endl;
fIsValid = false; fIsValid = false;
if (tokens) {
delete tokens;
tokens = 0;
}
break; break;
} }
fScanHigh = str.Atof(); fScanHigh = str.Atof();
@ -659,6 +811,7 @@ Bool_t PFitter::CheckCommands()
cerr << endl << ">> Will stop ..."; cerr << endl << ">> Will stop ...";
cerr << endl; cerr << endl;
fIsValid = false; fIsValid = false;
break;
} }
} }
@ -776,6 +929,72 @@ Bool_t PFitter::ExecuteContours()
return true; return true;
} }
//--------------------------------------------------------------------------
// ExecuteFitRange
//--------------------------------------------------------------------------
/**
* <p>Change the fit range via command block.
*
* \param lineNo the line number of the command block
*
* <b>return:</b> true if done, otherwise returns false.
*/
Bool_t PFitter::ExecuteFitRange(UInt_t lineNo)
{
cout << ">> PFitter::ExecuteFitRange(): " << fCmdLines[lineNo].fLine.Data() << endl;
TObjArray *tokens = 0;
TObjString *ostr;
TString str;
tokens = fCmdLines[lineNo].fLine.Tokenize(", \t");
PMsrRunList *runList = fRunInfo->GetMsrRunList();
// execute command, no error checking needed since this has been already carried out in CheckCommands()
if (tokens->GetEntries() == 2) { // reset command
fRunListCollection->SetFitRange(fOriginalFitRange);
} else if (tokens->GetEntries() == 3) { // single fit range for all runs
Double_t start = 0.0, end = 0.0;
PDoublePair fitRange;
PDoublePairVector fitRangeVector;
ostr = dynamic_cast<TObjString*>(tokens->At(1));
str = ostr->GetString();
start = str.Atof();
ostr = dynamic_cast<TObjString*>(tokens->At(2));
str = ostr->GetString();
end = str.Atof();
fitRange.first = start;
fitRange.second = end;
fitRangeVector.push_back(fitRange);
fRunListCollection->SetFitRange(fitRangeVector);
} else { // individual fit ranges for each run
Double_t start = 0.0, end = 0.0;
PDoublePair fitRange;
PDoublePairVector fitRangeVector;
for (UInt_t i=0; i<runList->size(); i++) {
ostr = dynamic_cast<TObjString*>(tokens->At(2*i+1));
str = ostr->GetString();
start = str.Atof();
ostr = dynamic_cast<TObjString*>(tokens->At(2*i+2));
str = ostr->GetString();
end = str.Atof();
fitRange.first = start;
fitRange.second = end;
fitRangeVector.push_back(fitRange);
}
fRunListCollection->SetFitRange(fitRangeVector);
}
return true;
}
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
// ExecuteFix // ExecuteFix
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
@ -788,14 +1007,14 @@ Bool_t PFitter::ExecuteContours()
*/ */
Bool_t PFitter::ExecuteFix(UInt_t lineNo) Bool_t PFitter::ExecuteFix(UInt_t lineNo)
{ {
cout << ">> PFitter::ExecuteFix(): " << fCmdLines[lineNo].fLine.Data() << endl;
TObjArray *tokens = 0; TObjArray *tokens = 0;
TObjString *ostr; TObjString *ostr;
TString str; TString str;
tokens = fCmdLines[lineNo].fLine.Tokenize(", \t"); tokens = fCmdLines[lineNo].fLine.Tokenize(", \t");
cout << ">> PFitter::ExecuteFix(): " << fCmdLines[lineNo].fLine.Data() << endl;
// Check if there is already a function minimum, i.e. migrad, minimization, or simplex has been called previously. // Check if there is already a function minimum, i.e. migrad, minimization, or simplex has been called previously.
// If so, update minuit2 user parameters // If so, update minuit2 user parameters
if (fFcnMin != 0) { if (fFcnMin != 0) {

33
src/classes/PRunAsymmetry.cpp
View File

@ -44,8 +44,6 @@
*/ */
PRunAsymmetry::PRunAsymmetry() : PRunBase() PRunAsymmetry::PRunAsymmetry() : PRunBase()
{ {
fFitStartTime = 0.0;
fFitStopTime = 0.0;
fNoOfFitBins = 0; fNoOfFitBins = 0;
} }
@ -159,7 +157,7 @@ Double_t PRunAsymmetry::CalcChiSquare(const std::vector<Double_t>& par)
Double_t time; Double_t time;
for (UInt_t i=0; i<fData.GetValue()->size(); i++) { for (UInt_t i=0; i<fData.GetValue()->size(); i++) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
if ((time>=fFitStartTime) && (time<=fFitStopTime)) { if ((time>=fFitStartTime) && (time<=fFitEndTime)) {
switch (fAlphaBetaTag) { switch (fAlphaBetaTag) {
case 1: // alpha == 1, beta == 1 case 1: // alpha == 1, beta == 1
asymFcnValue = fTheory->Func(time, par, fFuncValues); asymFcnValue = fTheory->Func(time, par, fFuncValues);
@ -206,6 +204,27 @@ Double_t PRunAsymmetry::CalcMaxLikelihood(const std::vector<Double_t>& par)
return 1.0; return 1.0;
} }
//--------------------------------------------------------------------------
// GetNoOfFitBins (public)
//--------------------------------------------------------------------------
/**
* <p>Calculate the number of fitted bins for the current fit range.
*
* <b>return:</b> number of fitted bins.
*/
UInt_t PRunAsymmetry::GetNoOfFitBins()
{
Double_t time;
fNoOfFitBins=0;
for (UInt_t i=0; i<fData.GetValue()->size(); i++) {
time = fData.GetDataTimeStart() + (Double_t)i * fData.GetDataTimeStep();
if ((time >= fFitStartTime) && (time <= fFitEndTime))
fNoOfFitBins++;
}
return fNoOfFitBins;
}
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
// CalcTheory // CalcTheory
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
@ -298,10 +317,6 @@ Bool_t PRunAsymmetry::PrepareData()
// keep the time resolution in (us) // keep the time resolution in (us)
fTimeResolution = runData->GetTimeResolution()/1.0e3; fTimeResolution = runData->GetTimeResolution()/1.0e3;
// keep start/stop time for fit
fFitStartTime = fRunInfo->GetFitRange(0);
fFitStopTime = fRunInfo->GetFitRange(1);
// collect histogram numbers // collect histogram numbers
PUIntVector forwardHistoNo; PUIntVector forwardHistoNo;
PUIntVector backwardHistoNo; PUIntVector backwardHistoNo;
@ -918,7 +933,7 @@ Bool_t PRunAsymmetry::PrepareFitData(PRawRunData* runData, UInt_t histoNo[2])
fNoOfFitBins=0; fNoOfFitBins=0;
for (UInt_t i=0; i<fData.GetValue()->size(); i++) { for (UInt_t i=0; i<fData.GetValue()->size(); i++) {
time = fData.GetDataTimeStart() + (Double_t)i * fData.GetDataTimeStep(); time = fData.GetDataTimeStart() + (Double_t)i * fData.GetDataTimeStep();
if ((time >= fFitStartTime) && (time <= fFitStopTime)) if ((time >= fFitStartTime) && (time <= fFitEndTime))
fNoOfFitBins++; fNoOfFitBins++;
} }
@ -1147,7 +1162,7 @@ Bool_t PRunAsymmetry::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2])
fNoOfFitBins=0; fNoOfFitBins=0;
for (UInt_t i=0; i<fData.GetValue()->size(); i++) { for (UInt_t i=0; i<fData.GetValue()->size(); i++) {
time = fData.GetDataTimeStart() + (Double_t)i * fData.GetDataTimeStep(); time = fData.GetDataTimeStart() + (Double_t)i * fData.GetDataTimeStep();
if ((time >= fFitStartTime) && (time <= fFitStopTime)) if ((time >= fFitStartTime) && (time <= fFitEndTime))
fNoOfFitBins++; fNoOfFitBins++;
} }

51
src/classes/PRunBase.cpp
View File

@ -56,6 +56,9 @@ PRunBase::PRunBase()
fRawData = 0; fRawData = 0;
fTimeResolution = -1.0; fTimeResolution = -1.0;
fFitStartTime = 0.0;
fFitEndTime = 0.0;
fValid = true; fValid = true;
fHandleTag = kEmpty; fHandleTag = kEmpty;
} }
@ -106,6 +109,10 @@ PRunBase::PRunBase(PMsrHandler *msrInfo, PRunDataHandler *rawData, UInt_t runNo,
cerr << endl << "**SEVERE ERROR** PRunBase::PRunBase: Theory is not valid :-(, will quit" << endl; cerr << endl << "**SEVERE ERROR** PRunBase::PRunBase: Theory is not valid :-(, will quit" << endl;
exit(0); exit(0);
} }
// set fit time ranges
fFitStartTime = fRunInfo->GetFitRange(0);
fFitEndTime = fRunInfo->GetFitRange(1);
} }
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
@ -123,8 +130,50 @@ PRunBase::~PRunBase()
fFuncValues.clear(); fFuncValues.clear();
} }
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
// CleanUp // SetFitRange (public)
//--------------------------------------------------------------------------
/**
* <p> Sets the current fit range, and recalculated the number of fitted bins
*
* \param fitRange vector with fit ranges
*/
void PRunBase::SetFitRange(PDoublePairVector fitRange)
{
Double_t start=0.0, end=0.0;
assert(fitRange.size()); // make sure fitRange is not empty
if (fitRange.size()==1) { // one fit range for all
start = fitRange[0].first;
end = fitRange[0].second;
} else {
// check that fRunNo is found within fitRange
if (fRunNo < static_cast<Int_t>(fitRange.size())) { // fRunNo present
start = fitRange[fRunNo].first;
end = fitRange[fRunNo].second;
} else { // fRunNo NOT present
cerr << endl << ">> PRunBase::SetFitRange(): **ERROR** msr-file run entry " << fRunNo << " not present in fit range vector.";
cerr << endl << ">> Will not do anything! Please check, this shouldn't happen." << endl;
return;
}
}
// check that start is before end
if (start > end) {
cerr << endl << ">> PRunBase::SetFitRange(): **WARNING** start=" << start << " is > as end=" << end;
cerr << endl << ">> Will swap them, since otherwise chisq/logLH == 0.0" << endl;
fFitStartTime = end;
fFitEndTime = start;
} else {
fFitStartTime = start;
fFitEndTime = end;
}
}
//--------------------------------------------------------------------------
// CleanUp (public)
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
/** /**
* <p> Clean up all locally allocate memory * <p> Clean up all locally allocate memory

22
src/classes/PRunListCollection.cpp
View File

@ -127,6 +127,28 @@ Bool_t PRunListCollection::Add(Int_t runNo, EPMusrHandleTag tag)
return success; return success;
} }
//--------------------------------------------------------------------------
// SetFitRange (public)
//--------------------------------------------------------------------------
/**
* <p>Set the current fit range. If fitRange.size()==1 the given fit range will be used for all the runs,
* otherwise fitRange.size()==the number of runs in the msr-file, and for each run there will be an induvidual
* fit range.
*
* \param fitRange vector holding the fit range(s).
*/
void PRunListCollection::SetFitRange(const PDoublePairVector fitRange)
{
for (UInt_t i=0; i<fRunSingleHistoList.size(); i++)
fRunSingleHistoList[i]->SetFitRange(fitRange);
for (UInt_t i=0; i<fRunAsymmetryList.size(); i++)
fRunAsymmetryList[i]->SetFitRange(fitRange);
for (UInt_t i=0; i<fRunMuMinusList.size(); i++)
fRunMuMinusList[i]->SetFitRange(fitRange);
for (UInt_t i=0; i<fRunNonMusrList.size(); i++)
fRunNonMusrList[i]->SetFitRange(fitRange);
}
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
// GetSingleHistoChisq // GetSingleHistoChisq
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------

15
src/classes/PRunMuMinus.cpp
View File

@ -41,8 +41,6 @@
*/ */
PRunMuMinus::PRunMuMinus() : PRunBase() PRunMuMinus::PRunMuMinus() : PRunBase()
{ {
fFitStartTime = 0.0;
fFitStopTime = 0.0;
fNoOfFitBins = 0; fNoOfFitBins = 0;
fHandleTag = kEmpty; fHandleTag = kEmpty;
@ -116,6 +114,19 @@ Double_t PRunMuMinus::CalcMaxLikelihood(const std::vector<Double_t>& par)
return 1.0; return 1.0;
} }
//--------------------------------------------------------------------------
// GetNoOfFitBins (public)
//--------------------------------------------------------------------------
/**
* <p>Calculate the number of fitted bins for the current fit range.
*
* <b>return:</b> number of fitted bins.
*/
UInt_t PRunMuMinus::GetNoOfFitBins()
{
return fNoOfFitBins;
}
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
// CalcTheory // CalcTheory
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------

31
src/classes/PRunNonMusr.cpp
View File

@ -41,8 +41,6 @@
*/ */
PRunNonMusr::PRunNonMusr() : PRunBase() PRunNonMusr::PRunNonMusr() : PRunBase()
{ {
fFitStartTime = 0.0;
fFitStopTime = 0.0;
fNoOfFitBins = 0; fNoOfFitBins = 0;
fHandleTag = kEmpty; fHandleTag = kEmpty;
@ -111,7 +109,7 @@ Double_t PRunNonMusr::CalcChiSquare(const std::vector<Double_t>& par)
Double_t x; Double_t x;
for (UInt_t i=0; i<fData.GetValue()->size(); i++) { for (UInt_t i=0; i<fData.GetValue()->size(); i++) {
x = fData.GetX()->at(i); x = fData.GetX()->at(i);
if ((x>=fFitStartTime) && (x<=fFitStopTime)) { if ((x>=fFitStartTime) && (x<=fFitEndTime)) {
diff = fData.GetValue()->at(i) - fTheory->Func(x, par, fFuncValues); diff = fData.GetValue()->at(i) - fTheory->Func(x, par, fFuncValues);
chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i)); chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i));
} }
@ -145,6 +143,27 @@ void PRunNonMusr::CalcTheory()
{ {
} }
//--------------------------------------------------------------------------
// GetNoOfFitBins (public)
//--------------------------------------------------------------------------
/**
* <p>Calculate the number of fitted bins for the current fit range.
*
* <b>return:</b> number of fitted bins.
*/
UInt_t PRunNonMusr::GetNoOfFitBins()
{
fNoOfFitBins=0;
Double_t x;
for (UInt_t i=0; i<fData.GetValue()->size(); i++) {
x = fData.GetX()->at(i);
if ((x >= fFitStartTime) && (x <= fFitEndTime))
fNoOfFitBins++;
}
return fNoOfFitBins;
}
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
// PrepareData // PrepareData
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
@ -187,10 +206,6 @@ Bool_t PRunNonMusr::PrepareFitData()
{ {
Bool_t success = true; Bool_t success = true;
// keep start/stop time for fit: here the meaning is of course start x, stop x
fFitStartTime = fRunInfo->GetFitRange(0);
fFitStopTime = fRunInfo->GetFitRange(1);
// get x-, y-index // get x-, y-index
UInt_t xIndex = GetXIndex(); UInt_t xIndex = GetXIndex();
UInt_t yIndex = GetYIndex(); UInt_t yIndex = GetYIndex();
@ -222,7 +237,7 @@ Bool_t PRunNonMusr::PrepareFitData()
Double_t x; Double_t x;
for (UInt_t i=0; i<fData.GetValue()->size(); i++) { for (UInt_t i=0; i<fData.GetValue()->size(); i++) {
x = fData.GetX()->at(i); x = fData.GetX()->at(i);
if ((x >= fFitStartTime) && (x <= fFitStopTime)) if ((x >= fFitStartTime) && (x <= fFitEndTime))
fNoOfFitBins++; fNoOfFitBins++;
} }

37
src/classes/PRunSingleHisto.cpp
View File

@ -43,8 +43,6 @@
*/ */
PRunSingleHisto::PRunSingleHisto() : PRunBase() PRunSingleHisto::PRunSingleHisto() : PRunBase()
{ {
fFitStartTime = 0.0;
fFitStopTime = 0.0;
fNoOfFitBins = 0; fNoOfFitBins = 0;
} }
@ -136,7 +134,7 @@ Double_t PRunSingleHisto::CalcChiSquare(const std::vector<Double_t>& par)
Double_t time; Double_t time;
for (UInt_t i=0; i<fData.GetValue()->size(); i++) { for (UInt_t i=0; i<fData.GetValue()->size(); i++) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
if ((time>=fFitStartTime) && (time<=fFitStopTime)) { if ((time>=fFitStartTime) && (time<=fFitEndTime)) {
diff = fData.GetValue()->at(i) - diff = fData.GetValue()->at(i) -
(N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg); (N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg);
chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i)); chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i));
@ -204,7 +202,7 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector<Double_t>& par)
Double_t time; Double_t time;
for (UInt_t i=0; i<fData.GetValue()->size(); i++) { for (UInt_t i=0; i<fData.GetValue()->size(); i++) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
if ((time>=fFitStartTime) && (time<=fFitStopTime)) { if ((time>=fFitStartTime) && (time<=fFitEndTime)) {
// calculate theory for the given parameter set // calculate theory for the given parameter set
theo = N0*TMath::Exp(-time/tau)*(1+fTheory->Func(time, par, fFuncValues))+bkg; theo = N0*TMath::Exp(-time/tau)*(1+fTheory->Func(time, par, fFuncValues))+bkg;
// check if data value is not too small // check if data value is not too small
@ -284,6 +282,27 @@ void PRunSingleHisto::CalcTheory()
par.clear(); par.clear();
} }
//--------------------------------------------------------------------------
// GetNoOfFitBins (public)
//--------------------------------------------------------------------------
/**
* <p>Calculate the number of fitted bins for the current fit range.
*
* <b>return:</b> number of fitted bins.
*/
UInt_t PRunSingleHisto::GetNoOfFitBins()
{
fNoOfFitBins=0;
Double_t time;
for (UInt_t i=0; i<fData.GetValue()->size(); i++) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
if ((time >= fFitStartTime) && (time <= fFitEndTime))
fNoOfFitBins++;
}
return fNoOfFitBins;
}
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
// PrepareData // PrepareData
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
@ -511,10 +530,6 @@ Bool_t PRunSingleHisto::PrepareData()
*/ */
Bool_t PRunSingleHisto::PrepareFitData(PRawRunData* runData, const UInt_t histoNo) Bool_t PRunSingleHisto::PrepareFitData(PRawRunData* runData, const UInt_t histoNo)
{ {
// keep start/stop time for fit
fFitStartTime = fRunInfo->GetFitRange(0);
fFitStopTime = fRunInfo->GetFitRange(1);
// transform raw histo data. This is done the following way (for details see the manual): // transform raw histo data. This is done the following way (for details see the manual):
// for the single histo fit, just the rebinned raw data are copied // for the single histo fit, just the rebinned raw data are copied
@ -622,7 +637,7 @@ Bool_t PRunSingleHisto::PrepareFitData(PRawRunData* runData, const UInt_t histoN
Double_t time; Double_t time;
for (UInt_t i=0; i<fData.GetValue()->size(); i++) { for (UInt_t i=0; i<fData.GetValue()->size(); i++) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
if ((time >= fFitStartTime) && (time <= fFitStopTime)) if ((time >= fFitStartTime) && (time <= fFitEndTime))
fNoOfFitBins++; fNoOfFitBins++;
} }
@ -721,7 +736,7 @@ Bool_t PRunSingleHisto::PrepareRawViewData(PRawRunData* runData, const UInt_t hi
Double_t time; Double_t time;
for (UInt_t i=0; i<fData.GetValue()->size(); i++) { for (UInt_t i=0; i<fData.GetValue()->size(); i++) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
if ((time >= fFitStartTime) && (time <= fFitStopTime)) if ((time >= fFitStartTime) && (time <= fFitEndTime))
fNoOfFitBins++; fNoOfFitBins++;
} }
@ -1013,7 +1028,7 @@ Bool_t PRunSingleHisto::PrepareViewData(PRawRunData* runData, const UInt_t histo
fNoOfFitBins=0; fNoOfFitBins=0;
for (UInt_t i=0; i<fData.GetValue()->size(); i++) { for (UInt_t i=0; i<fData.GetValue()->size(); i++) {
time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep();
if ((time >= fFitStartTime) && (time <= fFitStopTime)) if ((time >= fFitStartTime) && (time <= fFitEndTime))
fNoOfFitBins++; fNoOfFitBins++;
} }

37
src/include/PFitter.h
View File

@ -43,22 +43,23 @@
#define PMN_INTERACTIVE 0 #define PMN_INTERACTIVE 0
#define PMN_CONTOURS 1 #define PMN_CONTOURS 1
#define PMN_EIGEN 2 #define PMN_EIGEN 2
#define PMN_FIX 3 #define PMN_FIT_RANGE 3
#define PMN_HESSE 4 #define PMN_FIX 4
#define PMN_MACHINE_PRECISION 5 #define PMN_HESSE 5
#define PMN_MIGRAD 6 #define PMN_MACHINE_PRECISION 6
#define PMN_MINIMIZE 7 #define PMN_MIGRAD 7
#define PMN_MINOS 8 #define PMN_MINIMIZE 8
#define PMN_PLOT 9 #define PMN_MINOS 9
#define PMN_RELEASE 10 #define PMN_PLOT 10
#define PMN_RESTORE 11 #define PMN_RELEASE 11
#define PMN_SAVE 12 #define PMN_RESTORE 12
#define PMN_SCAN 13 #define PMN_SAVE 13
#define PMN_SIMPLEX 14 #define PMN_SCAN 14
#define PMN_STRATEGY 15 #define PMN_SIMPLEX 15
#define PMN_USER_COVARIANCE 16 #define PMN_STRATEGY 16
#define PMN_USER_PARAM_STATE 17 #define PMN_USER_COVARIANCE 17
#define PMN_PRINT 18 #define PMN_USER_PARAM_STATE 18
#define PMN_PRINT 19
/** /**
* <p>Interface class to minuit2. * <p>Interface class to minuit2.
@ -84,6 +85,7 @@ class PFitter
UInt_t fStrategy; ///< fitting strategy (see minuit2 manual). UInt_t fStrategy; ///< fitting strategy (see minuit2 manual).
PMsrHandler *fRunInfo; ///< pointer to the msr-file handler PMsrHandler *fRunInfo; ///< pointer to the msr-file handler
PRunListCollection *fRunListCollection; ///< pointer to the run list collection
PMsrParamList fParams; ///< msr-file parameters PMsrParamList fParams; ///< msr-file parameters
@ -103,11 +105,14 @@ class PFitter
Double_t fScanHigh; ///< scan range high. default=0.0 which means 2 std dev. (see MnScan/MnContours in the minuit2 user manual) Double_t fScanHigh; ///< scan range high. default=0.0 which means 2 std dev. (see MnScan/MnContours in the minuit2 user manual)
PDoublePairVector fScanData; ///< keeps the scan/contour data PDoublePairVector fScanData; ///< keeps the scan/contour data
PDoublePairVector fOriginalFitRange; ///< keeps the original fit range in case there is a range command in the COMMAND block
// commands // commands
Bool_t CheckCommands(); Bool_t CheckCommands();
Bool_t SetParameters(); Bool_t SetParameters();
Bool_t ExecuteContours(); Bool_t ExecuteContours();
Bool_t ExecuteFitRange(UInt_t lineNo);
Bool_t ExecuteFix(UInt_t lineNo); Bool_t ExecuteFix(UInt_t lineNo);
Bool_t ExecuteHesse(); Bool_t ExecuteHesse();
Bool_t ExecuteMigrad(); Bool_t ExecuteMigrad();

4
src/include/PRunAsymmetry.h
View File

@ -49,7 +49,7 @@ class PRunAsymmetry : public PRunBase
virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par); virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par);
virtual void CalcTheory(); virtual void CalcTheory();
virtual UInt_t GetNoOfFitBins() { return fNoOfFitBins; } ///< returns the number of bins to be fitted. virtual UInt_t GetNoOfFitBins();
protected: protected:
virtual Bool_t PrepareData(); virtual Bool_t PrepareData();
@ -60,8 +60,6 @@ class PRunAsymmetry : public PRunBase
private: private:
UInt_t fAlphaBetaTag; ///< \f$ 1 \to \alpha = \beta = 1\f$; \f$ 2 \to \alpha \neq 1, \beta = 1\f$; \f$ 3 \to \alpha = 1, \beta \neq 1\f$; \f$ 4 \to \alpha \neq 1, \beta \neq 1\f$. UInt_t fAlphaBetaTag; ///< \f$ 1 \to \alpha = \beta = 1\f$; \f$ 2 \to \alpha \neq 1, \beta = 1\f$; \f$ 3 \to \alpha = 1, \beta \neq 1\f$; \f$ 4 \to \alpha \neq 1, \beta \neq 1\f$.
Double_t fFitStartTime; ///< fit start time
Double_t fFitStopTime; ///< fit stop time
UInt_t fNoOfFitBins; ///< number of bins to be be fitted UInt_t fNoOfFitBins; ///< number of bins to be be fitted
PDoubleVector fForward; ///< forward histo data PDoubleVector fForward; ///< forward histo data

4
src/include/PRunBase.h
View File

@ -55,6 +55,7 @@ class PRunBase
virtual Double_t CalcChiSquare(const vector<Double_t>& par) = 0; ///< pure virtual, i.e. needs to be implemented by the deriving class!! virtual Double_t CalcChiSquare(const vector<Double_t>& par) = 0; ///< pure virtual, i.e. needs to be implemented by the deriving class!!
virtual Double_t CalcMaxLikelihood(const vector<Double_t>& par) = 0; ///< pure virtual, i.e. needs to be implemented by the deriving class!! virtual Double_t CalcMaxLikelihood(const vector<Double_t>& par) = 0; ///< pure virtual, i.e. needs to be implemented by the deriving class!!
virtual void SetFitRange(PDoublePairVector fitRange);
virtual void CalcTheory() = 0; ///< pure virtual, i.e. needs to be implemented by the deriving class!! virtual void CalcTheory() = 0; ///< pure virtual, i.e. needs to be implemented by the deriving class!!
@ -79,6 +80,9 @@ class PRunBase
Double_t fTimeResolution; ///< time resolution in (us) Double_t fTimeResolution; ///< time resolution in (us)
PIntVector fT0s; ///< all t0's of a run! The derived classes will handle it PIntVector fT0s; ///< all t0's of a run! The derived classes will handle it
Double_t fFitStartTime; ///< fit start time
Double_t fFitEndTime; ///< fit end time
PDoubleVector fFuncValues; ///< is keeping the values of the functions from the FUNCTIONS block PDoubleVector fFuncValues; ///< is keeping the values of the functions from the FUNCTIONS block
PTheory *fTheory; ///< theory needed to calculate chi-square PTheory *fTheory; ///< theory needed to calculate chi-square

2
src/include/PRunListCollection.h
View File

@ -56,6 +56,8 @@ class PRunListCollection
virtual Bool_t Add(Int_t runNo, EPMusrHandleTag tag); virtual Bool_t Add(Int_t runNo, EPMusrHandleTag tag);
virtual void SetFitRange(const PDoublePairVector fitRange);
virtual Double_t GetSingleHistoChisq(const std::vector<Double_t>& par) const; virtual Double_t GetSingleHistoChisq(const std::vector<Double_t>& par) const;
virtual Double_t GetAsymmetryChisq(const std::vector<Double_t>& par) const; virtual Double_t GetAsymmetryChisq(const std::vector<Double_t>& par) const;
virtual Double_t GetMuMinusChisq(const std::vector<Double_t>& par) const; virtual Double_t GetMuMinusChisq(const std::vector<Double_t>& par) const;

4
src/include/PRunMuMinus.h
View File

@ -48,14 +48,12 @@ class PRunMuMinus : public PRunBase
virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par); virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par);
virtual void CalcTheory(); virtual void CalcTheory();
virtual UInt_t GetNoOfFitBins() { return fNoOfFitBins; } virtual UInt_t GetNoOfFitBins();
protected: protected:
virtual Bool_t PrepareData(); virtual Bool_t PrepareData();
private: private:
Double_t fFitStartTime; ///< fit start time
Double_t fFitStopTime; ///< fit end time
UInt_t fNoOfFitBins; ///< number of bins to be fitted UInt_t fNoOfFitBins; ///< number of bins to be fitted
}; };

4
src/include/PRunNonMusr.h
View File

@ -49,7 +49,7 @@ class PRunNonMusr : public PRunBase
virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par); virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par);
virtual void CalcTheory(); virtual void CalcTheory();
virtual UInt_t GetNoOfFitBins() { return fNoOfFitBins; } ///< returns the number of fitted bins virtual UInt_t GetNoOfFitBins();
virtual UInt_t GetXIndex(); virtual UInt_t GetXIndex();
virtual UInt_t GetYIndex(); virtual UInt_t GetYIndex();
@ -62,8 +62,6 @@ class PRunNonMusr : public PRunBase
private: private:
PRawRunData *fRawRunData; ///< raw run data handler PRawRunData *fRawRunData; ///< raw run data handler
Double_t fFitStartTime; ///< fit start time
Double_t fFitStopTime; ///< fit stop time
UInt_t fNoOfFitBins; ///< number of bins to be be fitted UInt_t fNoOfFitBins; ///< number of bins to be be fitted
}; };

4
src/include/PRunSingleHisto.h
View File

@ -48,7 +48,7 @@ class PRunSingleHisto : public PRunBase
virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par); virtual Double_t CalcMaxLikelihood(const std::vector<Double_t>& par);
virtual void CalcTheory(); virtual void CalcTheory();
virtual UInt_t GetNoOfFitBins() { return fNoOfFitBins; } ///< returns the number of bins fitted virtual UInt_t GetNoOfFitBins();
protected: protected:
virtual Bool_t PrepareData(); virtual Bool_t PrepareData();
@ -57,8 +57,6 @@ class PRunSingleHisto : public PRunBase
virtual Bool_t PrepareViewData(PRawRunData* runData, const UInt_t histoNo); virtual Bool_t PrepareViewData(PRawRunData* runData, const UInt_t histoNo);
private: private:
Double_t fFitStartTime; ///< fit start time
Double_t fFitStopTime; ///< fit stop time
UInt_t fNoOfFitBins; ///< number of bins to be fitted UInt_t fNoOfFitBins; ///< number of bins to be fitted
Double_t fBackground; ///< needed if background range is given (units: 1/bin) Double_t fBackground; ///< needed if background range is given (units: 1/bin)