From 5cd4283d30111e1fed86d26fc9bdf1906724304d Mon Sep 17 00:00:00 2001 From: nemu Date: Wed, 20 Oct 2010 09:24:12 +0000 Subject: [PATCH] added FIT_RANGE RESET | start end | s1 e1 s2 e2 .. sN eN command to the COMMAND block --- ChangeLog | 5 +- configure.ac | 4 +- src/classes/PFitter.cpp | 225 ++++++++++++++++++++++++++++- src/classes/PRunAsymmetry.cpp | 33 +++-- src/classes/PRunBase.cpp | 51 ++++++- src/classes/PRunListCollection.cpp | 22 +++ src/classes/PRunMuMinus.cpp | 15 +- src/classes/PRunNonMusr.cpp | 31 +++- src/classes/PRunSingleHisto.cpp | 37 +++-- src/include/PFitter.h | 37 +++-- src/include/PRunAsymmetry.h | 4 +- src/include/PRunBase.h | 4 + src/include/PRunListCollection.h | 2 + src/include/PRunMuMinus.h | 4 +- src/include/PRunNonMusr.h | 4 +- src/include/PRunSingleHisto.h | 4 +- 16 files changed, 416 insertions(+), 66 deletions(-) diff --git a/ChangeLog b/ChangeLog index e9dd95be..19000fd5 100644 --- a/ChangeLog +++ b/ChangeLog @@ -4,8 +4,9 @@ # $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 implemented more checks on the integrity of the msr-file NEW [MUSRGUI/MUSREDIT] msr2data default option flags in musredit_startup.xml diff --git a/configure.ac b/configure.ac index cfc3e7b3..52a71f1b 100644 --- a/configure.ac +++ b/configure.ac @@ -1,5 +1,5 @@ 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_CANONICAL_HOST #AC_MSG_RESULT([${host} ${host_cpu} ${host_vendor} ${host_os}]) @@ -28,7 +28,7 @@ dnl ----------------------------------------------- #release versioning MUSR_MAJOR_VERSION=0 -MUSR_MINOR_VERSION=6 +MUSR_MINOR_VERSION=7 MUSR_MICRO_VERSION=0 #release versioning diff --git a/src/classes/PFitter.cpp b/src/classes/PFitter.cpp index e93f5139..d18e647a 100644 --- a/src/classes/PFitter.cpp +++ b/src/classes/PFitter.cpp @@ -69,7 +69,7 @@ using namespace std; * \param chisq_only flag: true=calculate chisq only (no fitting) */ PFitter::PFitter(PMsrHandler *runInfo, PRunListCollection *runListCollection, Bool_t chisq_only) : - fChisqOnly(chisq_only), fRunInfo(runInfo) + fChisqOnly(chisq_only), fRunInfo(runInfo), fRunListCollection(runListCollection) { // initialize variables fIsScanOnly = true; @@ -92,6 +92,15 @@ PFitter::PFitter(PMsrHandler *runInfo, PRunListCollection *runListCollection, Bo fScanLow = 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; isize(); i++) { + range.first = (*runs)[i].GetFitRange(0); + range.second = (*runs)[i].GetFitRange(1); + fOriginalFitRange.push_back(range); + } + // check msr minuit commands if (!CheckCommands()) { return; @@ -183,6 +192,9 @@ Bool_t PFitter::DoFit() case PMN_CONTOURS: status = ExecuteContours(); break; + case PMN_FIT_RANGE: + status = ExecuteFitRange(fCmdList[i].second); + break; case PMN_FIX: status = ExecuteFix(fCmdList[i].second); 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; fIsValid = false; + if (tokens) { + delete tokens; + tokens = 0; + } break; } 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; fIsValid = false; + if (tokens) { + delete tokens; + tokens = 0; + } break; } // 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; fIsValid = false; + if (tokens) { + delete tokens; + tokens = 0; + } break; } 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; fIsValid = false; + if (tokens) { + delete tokens; + tokens = 0; + } break; } fScanNoPoints = ival; @@ -369,6 +397,90 @@ Bool_t PFitter::CheckCommands() cmd.first = PMN_EIGEN; cmd.second = cmdLineNo; 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(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(tokens->GetEntries()) % 2) == 1) { + if ((tokens->GetEntries() > 3) && ((static_cast(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(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 | FIT_RANGE .. ,"; + 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)) { // check if the given set of parameters (number or names) is present 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; fIsValid = false; + if (tokens) { + delete tokens; + tokens = 0; + } break; } } 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; fIsValid = false; + if (tokens) { + delete tokens; + tokens = 0; + } break; } } @@ -472,6 +592,10 @@ Bool_t PFitter::CheckCommands() cerr << endl << ">> Parameter " << ival << " is out of the Parameter Range [1," << fParams.size() << "]"; cerr << endl; fIsValid = false; + if (tokens) { + delete tokens; + tokens = 0; + } break; } } 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; fIsValid = false; + if (tokens) { + delete tokens; + tokens = 0; + } 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; fIsValid = false; + if (tokens) { + delete tokens; + tokens = 0; + } break; } 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; fIsValid = false; + if (tokens) { + delete tokens; + tokens = 0; + } break; } // 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; fIsValid = false; + if (tokens) { + delete tokens; + tokens = 0; + } break; } 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; fIsValid = false; + if (tokens) { + delete tokens; + tokens = 0; + } break; } 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; fIsValid = false; + if (tokens) { + delete tokens; + tokens = 0; + } break; } 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; fIsValid = false; + if (tokens) { + delete tokens; + tokens = 0; + } break; } fScanHigh = str.Atof(); @@ -659,6 +811,7 @@ Bool_t PFitter::CheckCommands() cerr << endl << ">> Will stop ..."; cerr << endl; fIsValid = false; + break; } } @@ -776,6 +929,72 @@ Bool_t PFitter::ExecuteContours() return true; } +//-------------------------------------------------------------------------- +// ExecuteFitRange +//-------------------------------------------------------------------------- +/** + *

Change the fit range via command block. + * + * \param lineNo the line number of the command block + * + * return: 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(tokens->At(1)); + str = ostr->GetString(); + start = str.Atof(); + ostr = dynamic_cast(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; isize(); i++) { + ostr = dynamic_cast(tokens->At(2*i+1)); + str = ostr->GetString(); + start = str.Atof(); + ostr = dynamic_cast(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 //-------------------------------------------------------------------------- @@ -788,14 +1007,14 @@ Bool_t PFitter::ExecuteContours() */ Bool_t PFitter::ExecuteFix(UInt_t lineNo) { + cout << ">> PFitter::ExecuteFix(): " << fCmdLines[lineNo].fLine.Data() << endl; + TObjArray *tokens = 0; TObjString *ostr; TString str; 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. // If so, update minuit2 user parameters if (fFcnMin != 0) { diff --git a/src/classes/PRunAsymmetry.cpp b/src/classes/PRunAsymmetry.cpp index 3efef5ad..0a3ba780 100644 --- a/src/classes/PRunAsymmetry.cpp +++ b/src/classes/PRunAsymmetry.cpp @@ -44,8 +44,6 @@ */ PRunAsymmetry::PRunAsymmetry() : PRunBase() { - fFitStartTime = 0.0; - fFitStopTime = 0.0; fNoOfFitBins = 0; } @@ -159,7 +157,7 @@ Double_t PRunAsymmetry::CalcChiSquare(const std::vector& par) Double_t time; for (UInt_t i=0; isize(); i++) { time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); - if ((time>=fFitStartTime) && (time<=fFitStopTime)) { + if ((time>=fFitStartTime) && (time<=fFitEndTime)) { switch (fAlphaBetaTag) { case 1: // alpha == 1, beta == 1 asymFcnValue = fTheory->Func(time, par, fFuncValues); @@ -206,6 +204,27 @@ Double_t PRunAsymmetry::CalcMaxLikelihood(const std::vector& par) return 1.0; } +//-------------------------------------------------------------------------- +// GetNoOfFitBins (public) +//-------------------------------------------------------------------------- +/** + *

Calculate the number of fitted bins for the current fit range. + * + * return: number of fitted bins. + */ +UInt_t PRunAsymmetry::GetNoOfFitBins() +{ + Double_t time; + fNoOfFitBins=0; + for (UInt_t i=0; isize(); i++) { + time = fData.GetDataTimeStart() + (Double_t)i * fData.GetDataTimeStep(); + if ((time >= fFitStartTime) && (time <= fFitEndTime)) + fNoOfFitBins++; + } + + return fNoOfFitBins; +} + //-------------------------------------------------------------------------- // CalcTheory //-------------------------------------------------------------------------- @@ -298,10 +317,6 @@ Bool_t PRunAsymmetry::PrepareData() // keep the time resolution in (us) fTimeResolution = runData->GetTimeResolution()/1.0e3; - // keep start/stop time for fit - fFitStartTime = fRunInfo->GetFitRange(0); - fFitStopTime = fRunInfo->GetFitRange(1); - // collect histogram numbers PUIntVector forwardHistoNo; PUIntVector backwardHistoNo; @@ -918,7 +933,7 @@ Bool_t PRunAsymmetry::PrepareFitData(PRawRunData* runData, UInt_t histoNo[2]) fNoOfFitBins=0; for (UInt_t i=0; isize(); i++) { time = fData.GetDataTimeStart() + (Double_t)i * fData.GetDataTimeStep(); - if ((time >= fFitStartTime) && (time <= fFitStopTime)) + if ((time >= fFitStartTime) && (time <= fFitEndTime)) fNoOfFitBins++; } @@ -1147,7 +1162,7 @@ Bool_t PRunAsymmetry::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]) fNoOfFitBins=0; for (UInt_t i=0; isize(); i++) { time = fData.GetDataTimeStart() + (Double_t)i * fData.GetDataTimeStep(); - if ((time >= fFitStartTime) && (time <= fFitStopTime)) + if ((time >= fFitStartTime) && (time <= fFitEndTime)) fNoOfFitBins++; } diff --git a/src/classes/PRunBase.cpp b/src/classes/PRunBase.cpp index aa10a28a..1711b8a1 100644 --- a/src/classes/PRunBase.cpp +++ b/src/classes/PRunBase.cpp @@ -56,6 +56,9 @@ PRunBase::PRunBase() fRawData = 0; fTimeResolution = -1.0; + fFitStartTime = 0.0; + fFitEndTime = 0.0; + fValid = true; 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; exit(0); } + + // set fit time ranges + fFitStartTime = fRunInfo->GetFitRange(0); + fFitEndTime = fRunInfo->GetFitRange(1); } //-------------------------------------------------------------------------- @@ -123,8 +130,50 @@ PRunBase::~PRunBase() fFuncValues.clear(); } + //-------------------------------------------------------------------------- -// CleanUp +// SetFitRange (public) +//-------------------------------------------------------------------------- +/** + *

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(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) //-------------------------------------------------------------------------- /** *

Clean up all locally allocate memory diff --git a/src/classes/PRunListCollection.cpp b/src/classes/PRunListCollection.cpp index d8eaf188..a4a14fbf 100644 --- a/src/classes/PRunListCollection.cpp +++ b/src/classes/PRunListCollection.cpp @@ -127,6 +127,28 @@ Bool_t PRunListCollection::Add(Int_t runNo, EPMusrHandleTag tag) return success; } +//-------------------------------------------------------------------------- +// SetFitRange (public) +//-------------------------------------------------------------------------- +/** + *

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; iSetFitRange(fitRange); + for (UInt_t i=0; iSetFitRange(fitRange); + for (UInt_t i=0; iSetFitRange(fitRange); + for (UInt_t i=0; iSetFitRange(fitRange); +} + //-------------------------------------------------------------------------- // GetSingleHistoChisq //-------------------------------------------------------------------------- diff --git a/src/classes/PRunMuMinus.cpp b/src/classes/PRunMuMinus.cpp index f8cabe00..e8f14726 100644 --- a/src/classes/PRunMuMinus.cpp +++ b/src/classes/PRunMuMinus.cpp @@ -41,8 +41,6 @@ */ PRunMuMinus::PRunMuMinus() : PRunBase() { - fFitStartTime = 0.0; - fFitStopTime = 0.0; fNoOfFitBins = 0; fHandleTag = kEmpty; @@ -116,6 +114,19 @@ Double_t PRunMuMinus::CalcMaxLikelihood(const std::vector& par) return 1.0; } +//-------------------------------------------------------------------------- +// GetNoOfFitBins (public) +//-------------------------------------------------------------------------- +/** + *

Calculate the number of fitted bins for the current fit range. + * + * return: number of fitted bins. + */ +UInt_t PRunMuMinus::GetNoOfFitBins() +{ + return fNoOfFitBins; +} + //-------------------------------------------------------------------------- // CalcTheory //-------------------------------------------------------------------------- diff --git a/src/classes/PRunNonMusr.cpp b/src/classes/PRunNonMusr.cpp index eca151fa..6f5dd875 100644 --- a/src/classes/PRunNonMusr.cpp +++ b/src/classes/PRunNonMusr.cpp @@ -41,8 +41,6 @@ */ PRunNonMusr::PRunNonMusr() : PRunBase() { - fFitStartTime = 0.0; - fFitStopTime = 0.0; fNoOfFitBins = 0; fHandleTag = kEmpty; @@ -111,7 +109,7 @@ Double_t PRunNonMusr::CalcChiSquare(const std::vector& par) Double_t x; for (UInt_t i=0; isize(); 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); chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i)); } @@ -145,6 +143,27 @@ void PRunNonMusr::CalcTheory() { } +//-------------------------------------------------------------------------- +// GetNoOfFitBins (public) +//-------------------------------------------------------------------------- +/** + *

Calculate the number of fitted bins for the current fit range. + * + * return: number of fitted bins. + */ +UInt_t PRunNonMusr::GetNoOfFitBins() +{ + fNoOfFitBins=0; + Double_t x; + for (UInt_t i=0; isize(); i++) { + x = fData.GetX()->at(i); + if ((x >= fFitStartTime) && (x <= fFitEndTime)) + fNoOfFitBins++; + } + + return fNoOfFitBins; +} + //-------------------------------------------------------------------------- // PrepareData //-------------------------------------------------------------------------- @@ -187,10 +206,6 @@ Bool_t PRunNonMusr::PrepareFitData() { 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 UInt_t xIndex = GetXIndex(); UInt_t yIndex = GetYIndex(); @@ -222,7 +237,7 @@ Bool_t PRunNonMusr::PrepareFitData() Double_t x; for (UInt_t i=0; isize(); i++) { x = fData.GetX()->at(i); - if ((x >= fFitStartTime) && (x <= fFitStopTime)) + if ((x >= fFitStartTime) && (x <= fFitEndTime)) fNoOfFitBins++; } diff --git a/src/classes/PRunSingleHisto.cpp b/src/classes/PRunSingleHisto.cpp index ed8b9f6a..9c21a470 100644 --- a/src/classes/PRunSingleHisto.cpp +++ b/src/classes/PRunSingleHisto.cpp @@ -43,8 +43,6 @@ */ PRunSingleHisto::PRunSingleHisto() : PRunBase() { - fFitStartTime = 0.0; - fFitStopTime = 0.0; fNoOfFitBins = 0; } @@ -136,7 +134,7 @@ Double_t PRunSingleHisto::CalcChiSquare(const std::vector& par) Double_t time; for (UInt_t i=0; isize(); i++) { time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); - if ((time>=fFitStartTime) && (time<=fFitStopTime)) { + if ((time>=fFitStartTime) && (time<=fFitEndTime)) { diff = fData.GetValue()->at(i) - (N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg); chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i)); @@ -204,7 +202,7 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector& par) Double_t time; for (UInt_t i=0; isize(); i++) { 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 theo = N0*TMath::Exp(-time/tau)*(1+fTheory->Func(time, par, fFuncValues))+bkg; // check if data value is not too small @@ -284,6 +282,27 @@ void PRunSingleHisto::CalcTheory() par.clear(); } +//-------------------------------------------------------------------------- +// GetNoOfFitBins (public) +//-------------------------------------------------------------------------- +/** + *

Calculate the number of fitted bins for the current fit range. + * + * return: number of fitted bins. + */ +UInt_t PRunSingleHisto::GetNoOfFitBins() +{ + fNoOfFitBins=0; + Double_t time; + for (UInt_t i=0; isize(); i++) { + time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); + if ((time >= fFitStartTime) && (time <= fFitEndTime)) + fNoOfFitBins++; + } + + return fNoOfFitBins; +} + //-------------------------------------------------------------------------- // PrepareData //-------------------------------------------------------------------------- @@ -511,10 +530,6 @@ Bool_t PRunSingleHisto::PrepareData() */ 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): // 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; for (UInt_t i=0; isize(); i++) { time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); - if ((time >= fFitStartTime) && (time <= fFitStopTime)) + if ((time >= fFitStartTime) && (time <= fFitEndTime)) fNoOfFitBins++; } @@ -721,7 +736,7 @@ Bool_t PRunSingleHisto::PrepareRawViewData(PRawRunData* runData, const UInt_t hi Double_t time; for (UInt_t i=0; isize(); i++) { time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); - if ((time >= fFitStartTime) && (time <= fFitStopTime)) + if ((time >= fFitStartTime) && (time <= fFitEndTime)) fNoOfFitBins++; } @@ -1013,7 +1028,7 @@ Bool_t PRunSingleHisto::PrepareViewData(PRawRunData* runData, const UInt_t histo fNoOfFitBins=0; for (UInt_t i=0; isize(); i++) { time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); - if ((time >= fFitStartTime) && (time <= fFitStopTime)) + if ((time >= fFitStartTime) && (time <= fFitEndTime)) fNoOfFitBins++; } diff --git a/src/include/PFitter.h b/src/include/PFitter.h index 43a803b0..f90d4f26 100644 --- a/src/include/PFitter.h +++ b/src/include/PFitter.h @@ -43,22 +43,23 @@ #define PMN_INTERACTIVE 0 #define PMN_CONTOURS 1 #define PMN_EIGEN 2 -#define PMN_FIX 3 -#define PMN_HESSE 4 -#define PMN_MACHINE_PRECISION 5 -#define PMN_MIGRAD 6 -#define PMN_MINIMIZE 7 -#define PMN_MINOS 8 -#define PMN_PLOT 9 -#define PMN_RELEASE 10 -#define PMN_RESTORE 11 -#define PMN_SAVE 12 -#define PMN_SCAN 13 -#define PMN_SIMPLEX 14 -#define PMN_STRATEGY 15 -#define PMN_USER_COVARIANCE 16 -#define PMN_USER_PARAM_STATE 17 -#define PMN_PRINT 18 +#define PMN_FIT_RANGE 3 +#define PMN_FIX 4 +#define PMN_HESSE 5 +#define PMN_MACHINE_PRECISION 6 +#define PMN_MIGRAD 7 +#define PMN_MINIMIZE 8 +#define PMN_MINOS 9 +#define PMN_PLOT 10 +#define PMN_RELEASE 11 +#define PMN_RESTORE 12 +#define PMN_SAVE 13 +#define PMN_SCAN 14 +#define PMN_SIMPLEX 15 +#define PMN_STRATEGY 16 +#define PMN_USER_COVARIANCE 17 +#define PMN_USER_PARAM_STATE 18 +#define PMN_PRINT 19 /** *

Interface class to minuit2. @@ -84,6 +85,7 @@ class PFitter UInt_t fStrategy; ///< fitting strategy (see minuit2 manual). PMsrHandler *fRunInfo; ///< pointer to the msr-file handler + PRunListCollection *fRunListCollection; ///< pointer to the run list collection 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) 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 Bool_t CheckCommands(); Bool_t SetParameters(); Bool_t ExecuteContours(); + Bool_t ExecuteFitRange(UInt_t lineNo); Bool_t ExecuteFix(UInt_t lineNo); Bool_t ExecuteHesse(); Bool_t ExecuteMigrad(); diff --git a/src/include/PRunAsymmetry.h b/src/include/PRunAsymmetry.h index a6e901c8..9f77cf69 100644 --- a/src/include/PRunAsymmetry.h +++ b/src/include/PRunAsymmetry.h @@ -49,7 +49,7 @@ class PRunAsymmetry : public PRunBase virtual Double_t CalcMaxLikelihood(const std::vector& par); virtual void CalcTheory(); - virtual UInt_t GetNoOfFitBins() { return fNoOfFitBins; } ///< returns the number of bins to be fitted. + virtual UInt_t GetNoOfFitBins(); protected: virtual Bool_t PrepareData(); @@ -60,8 +60,6 @@ class PRunAsymmetry : public PRunBase 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$. - Double_t fFitStartTime; ///< fit start time - Double_t fFitStopTime; ///< fit stop time UInt_t fNoOfFitBins; ///< number of bins to be be fitted PDoubleVector fForward; ///< forward histo data diff --git a/src/include/PRunBase.h b/src/include/PRunBase.h index f6fb3252..48d6dd1b 100644 --- a/src/include/PRunBase.h +++ b/src/include/PRunBase.h @@ -55,6 +55,7 @@ class PRunBase virtual Double_t CalcChiSquare(const vector& par) = 0; ///< pure virtual, i.e. needs to be implemented by the deriving class!! virtual Double_t CalcMaxLikelihood(const vector& 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!! @@ -79,6 +80,9 @@ class PRunBase Double_t fTimeResolution; ///< time resolution in (us) 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 PTheory *fTheory; ///< theory needed to calculate chi-square diff --git a/src/include/PRunListCollection.h b/src/include/PRunListCollection.h index c997fdd3..d424cc31 100644 --- a/src/include/PRunListCollection.h +++ b/src/include/PRunListCollection.h @@ -56,6 +56,8 @@ class PRunListCollection virtual Bool_t Add(Int_t runNo, EPMusrHandleTag tag); + virtual void SetFitRange(const PDoublePairVector fitRange); + virtual Double_t GetSingleHistoChisq(const std::vector& par) const; virtual Double_t GetAsymmetryChisq(const std::vector& par) const; virtual Double_t GetMuMinusChisq(const std::vector& par) const; diff --git a/src/include/PRunMuMinus.h b/src/include/PRunMuMinus.h index 49b8d8c7..4adabfcc 100644 --- a/src/include/PRunMuMinus.h +++ b/src/include/PRunMuMinus.h @@ -48,14 +48,12 @@ class PRunMuMinus : public PRunBase virtual Double_t CalcMaxLikelihood(const std::vector& par); virtual void CalcTheory(); - virtual UInt_t GetNoOfFitBins() { return fNoOfFitBins; } + virtual UInt_t GetNoOfFitBins(); protected: virtual Bool_t PrepareData(); private: - Double_t fFitStartTime; ///< fit start time - Double_t fFitStopTime; ///< fit end time UInt_t fNoOfFitBins; ///< number of bins to be fitted }; diff --git a/src/include/PRunNonMusr.h b/src/include/PRunNonMusr.h index 5f41d642..35a6cd76 100644 --- a/src/include/PRunNonMusr.h +++ b/src/include/PRunNonMusr.h @@ -49,7 +49,7 @@ class PRunNonMusr : public PRunBase virtual Double_t CalcMaxLikelihood(const std::vector& par); 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 GetYIndex(); @@ -62,8 +62,6 @@ class PRunNonMusr : public PRunBase private: 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 }; diff --git a/src/include/PRunSingleHisto.h b/src/include/PRunSingleHisto.h index c22d325b..38b01f6b 100644 --- a/src/include/PRunSingleHisto.h +++ b/src/include/PRunSingleHisto.h @@ -48,7 +48,7 @@ class PRunSingleHisto : public PRunBase virtual Double_t CalcMaxLikelihood(const std::vector& par); virtual void CalcTheory(); - virtual UInt_t GetNoOfFitBins() { return fNoOfFitBins; } ///< returns the number of bins fitted + virtual UInt_t GetNoOfFitBins(); protected: virtual Bool_t PrepareData(); @@ -57,8 +57,6 @@ class PRunSingleHisto : public PRunBase virtual Bool_t PrepareViewData(PRawRunData* runData, const UInt_t histoNo); private: - Double_t fFitStartTime; ///< fit start time - Double_t fFitStopTime; ///< fit stop time UInt_t fNoOfFitBins; ///< number of bins to be fitted Double_t fBackground; ///< needed if background range is given (units: 1/bin)