From ed0bfd516bd745647b2cfce5b937b870dedef4c3 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Fri, 31 Jan 2020 16:39:52 +0100 Subject: [PATCH 1/5] started to implement a chisq/maxLH sector command. Currently only the skeleton is in place. --- src/classes/PFitter.cpp | 330 ++++++++++++++++++++++++++++++++-------- src/include/PFitter.h | 40 +++++ 2 files changed, 308 insertions(+), 62 deletions(-) diff --git a/src/classes/PFitter.cpp b/src/classes/PFitter.cpp index c9f5db4b..2b301c78 100644 --- a/src/classes/PFitter.cpp +++ b/src/classes/PFitter.cpp @@ -67,6 +67,110 @@ #include "PFitter.h" + +//+++ PSectorChisq class +++++++++++++++++++++++++++++++++++++++++++++++++++ + +//-------------------------------------------------------------------------- +// Constructor +//-------------------------------------------------------------------------- +/** + *

Constructor. + */ +PSectorChisq::PSectorChisq(UInt_t noOfRuns) : fNoOfRuns(noOfRuns) +{ + // init + fFirst = 0.0; + fLast = 0.0; + fChisq = 0.0; + fNDF = 0; + fChisqRun.resize(fNoOfRuns); + fNDFRun.resize(fNoOfRuns); + for (UInt_t i=0; iSet the time range for one sector + * + * @param first time stamp of the fgb + * @param last time stamp of the requested sector end + */ +void PSectorChisq::SetTimeRange(Double_t first, Double_t last) +{ + // NOT YET IMPLEMENTED //as35 +} + +//-------------------------------------------------------------------------- +// SetChisq +//-------------------------------------------------------------------------- +/** + *

Set the chisq/maxLH for a given detector with index idx. + * + * @param chisq chisq/maxLH to be set + * @param idx index of the run to be set + */ +void PSectorChisq::SetChisq(Double_t chisq, UInt_t idx) +{ + // NOT YET IMPLEMENTED //as35 +} + +//-------------------------------------------------------------------------- +// SetNDF +//-------------------------------------------------------------------------- +/** + *

Set the NDF for a given detector with index idx. + * + * @param ndf to be set + * @param idx index of the run to be set + */ +void PSectorChisq::SetNDF(Double_t ndf, UInt_t idx) +{ + // NOT YET IMPLEMENTED //as35 +} + +//-------------------------------------------------------------------------- +// GetChisq +//-------------------------------------------------------------------------- +/** + *

chisq/maxLH of the run with index idx + * + * @param idx index of the run + * + * return: chisq/maxLH of the requested run or -1.0 if the idx is out-of-range. + */ +Double_t PSectorChisq::GetChisq(UInt_t idx) +{ + if (idx >= fNoOfRuns) + return -1.0; + + return fChisqRun[idx]; +} + +//-------------------------------------------------------------------------- +// GetChisq +//-------------------------------------------------------------------------- +/** + *

NDF of the run with index idx + * + * @param idx index of the run + * + * return: NDF of the requested run or 0.0 if the idx is out-of-range. + */ +UInt_t PSectorChisq::GetNDF(UInt_t idx) +{ + if (idx >= fNoOfRuns) + return 0; + + return fNDFRun[idx]; +} + +//+++ PFitter class ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + //-------------------------------------------------------------------------- // Constructor //-------------------------------------------------------------------------- @@ -103,12 +207,20 @@ PFitter::PFitter(PMsrHandler *runInfo, PRunListCollection *runListCollection, Bo fPrintLevel = 1.0; // keep all the fit ranges in case RANGE command is present + PDoublePair rangeGlob; + PMsrGlobalBlock *global = fRunInfo->GetMsrGlobal(); + rangeGlob.first = global->GetFitRange(0); + rangeGlob.second = global->GetFitRange(1); + 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); + if (range.first == PMUSR_UNDEFINED) + fOriginalFitRange.push_back(rangeGlob); + else + fOriginalFitRange.push_back(range); } // check msr minuit commands @@ -338,28 +450,37 @@ Bool_t PFitter::CheckCommands() PIntPair cmd; PMsrLines::iterator it; UInt_t cmdLineNo = 0; + TString line; + Ssiz_t pos; for (it = fCmdLines.begin(); it != fCmdLines.end(); ++it) { if (it == fCmdLines.begin()) cmdLineNo = 0; else cmdLineNo++; - if (it->fLine.Contains("COMMANDS", TString::kIgnoreCase)) { + + // strip potential comments + line = it->fLine; + pos = line.First('#'); + if (pos > 0) // comment present + line.Remove(pos,line.Length()-pos); + + if (line.Contains("COMMANDS", TString::kIgnoreCase)) { continue; - } else if (it->fLine.Contains("SET BATCH", TString::kIgnoreCase)) { // needed for backward compatibility + } else if (line.Contains("SET BATCH", TString::kIgnoreCase)) { // needed for backward compatibility continue; - } else if (it->fLine.Contains("END RETURN", TString::kIgnoreCase)) { // needed for backward compatibility + } else if (line.Contains("END RETURN", TString::kIgnoreCase)) { // needed for backward compatibility continue; - } else if (it->fLine.Contains("CHI_SQUARE", TString::kIgnoreCase)) { + } else if (line.Contains("CHI_SQUARE", TString::kIgnoreCase)) { continue; - } else if (it->fLine.Contains("MAX_LIKELIHOOD", TString::kIgnoreCase)) { + } else if (line.Contains("MAX_LIKELIHOOD", TString::kIgnoreCase)) { continue; - } else if (it->fLine.Contains("SCALE_N0_BKG", TString::kIgnoreCase)) { + } else if (line.Contains("SCALE_N0_BKG", TString::kIgnoreCase)) { continue; - } else if (it->fLine.Contains("INTERACTIVE", TString::kIgnoreCase)) { + } else if (line.Contains("INTERACTIVE", TString::kIgnoreCase)) { cmd.first = PMN_INTERACTIVE; cmd.second = cmdLineNo; fCmdList.push_back(cmd); - } else if (it->fLine.Contains("CONTOURS", TString::kIgnoreCase)) { + } else if (line.Contains("CONTOURS", TString::kIgnoreCase)) { cmd.first = PMN_CONTOURS; cmd.second = cmdLineNo; fCmdList.push_back(cmd); @@ -369,7 +490,7 @@ Bool_t PFitter::CheckCommands() TString str; UInt_t ival; - tokens = it->fLine.Tokenize(", \t"); + tokens = line.Tokenize(", \t"); for (Int_t i=0; iGetEntries(); i++) { ostr = dynamic_cast(tokens->At(i)); @@ -379,7 +500,7 @@ Bool_t PFitter::CheckCommands() // check that token is a UInt_t if (!str.IsDigit()) { std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo; - std::cerr << std::endl << ">> " << it->fLine.Data(); + std::cerr << std::endl << ">> " << line.Data(); std::cerr << std::endl << ">> parameter number is not number!"; std::cerr << std::endl << ">> command syntax for CONTOURS is: CONTOURS parameter-X parameter-Y [# of points]"; std::cerr << std::endl; @@ -394,7 +515,7 @@ Bool_t PFitter::CheckCommands() // check that parameter is within range if ((ival < 1) || (ival > fParams.size())) { std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo; - std::cerr << std::endl << ">> " << it->fLine.Data(); + std::cerr << std::endl << ">> " << line.Data(); std::cerr << std::endl << ">> parameter number is out of range [1," << fParams.size() << "]!"; std::cerr << std::endl << ">> command syntax for CONTOURS is: CONTOURS parameter-X parameter-Y [# of points]"; std::cerr << std::endl; @@ -414,7 +535,7 @@ Bool_t PFitter::CheckCommands() // check that token is a UInt_t if (!str.IsDigit()) { std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo; - std::cerr << std::endl << ">> " << it->fLine.Data(); + std::cerr << std::endl << ">> " << line.Data(); std::cerr << std::endl << ">> number of points is not number!"; std::cerr << std::endl << ">> command syntax for CONTOURS is: CONTOURS parameter-X parameter-Y [# of points]"; std::cerr << std::endl; @@ -428,7 +549,7 @@ Bool_t PFitter::CheckCommands() ival = str.Atoi(); if ((ival < 1) || (ival > 100)) { std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo; - std::cerr << std::endl << ">> " << it->fLine.Data(); + std::cerr << std::endl << ">> " << line.Data(); std::cerr << std::endl << ">> number of scan points is out of range [1,100]!"; std::cerr << std::endl << ">> command syntax for CONTOURS is: CONTOURS parameter-X parameter-Y [# of points]"; std::cerr << std::endl; @@ -447,11 +568,11 @@ Bool_t PFitter::CheckCommands() delete tokens; tokens = nullptr; } - } else if (it->fLine.Contains("EIGEN", TString::kIgnoreCase)) { + } else if (line.Contains("EIGEN", TString::kIgnoreCase)) { cmd.first = PMN_EIGEN; cmd.second = cmdLineNo; fCmdList.push_back(cmd); - } else if (it->fLine.Contains("FIT_RANGE", TString::kIgnoreCase)) { + } else if (line.Contains("FIT_RANGE", TString::kIgnoreCase)) { // check the 5 options: // (i) FIT_RANGE RESET, // (ii) FIT_RANGE start end, @@ -462,7 +583,7 @@ Bool_t PFitter::CheckCommands() TObjString *ostr; TString str; - tokens = it->fLine.Tokenize(", \t"); + tokens = line.Tokenize(", \t"); if (tokens->GetEntries() == 2) { // should only be RESET ostr = dynamic_cast(tokens->At(1)); @@ -473,7 +594,7 @@ Bool_t PFitter::CheckCommands() fCmdList.push_back(cmd); } else { std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo; - std::cerr << std::endl << ">> " << it->fLine.Data(); + std::cerr << std::endl << ">> " << line.Data(); std::cerr << std::endl << ">> Syntax: FIT_RANGE RESET | FIT_RANGE start end | FIT_RANGE s1 e1 s2 e2 .. sN eN,"; std::cerr << std::endl << ">> with N the number of runs in the msr-file." << std::endl; std::cerr << std::endl << ">> Found " << str.Data() << ", instead of RESET" << std::endl; @@ -487,7 +608,7 @@ Bool_t PFitter::CheckCommands() } else if ((tokens->GetEntries() > 1) && (static_cast(tokens->GetEntries()) % 2) == 1) { if ((tokens->GetEntries() > 3) && ((static_cast(tokens->GetEntries())-1)) != 2*fRunInfo->GetMsrRunList()->size()) { std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo; - std::cerr << std::endl << ">> " << it->fLine.Data(); + std::cerr << std::endl << ">> " << line.Data(); std::cerr << std::endl << ">> Syntax: FIT_RANGE RESET | FIT_RANGE | FIT_RANGE .. |"; std::cerr << std::endl << ">> FIT_RANGE fgb+ lgb- | FIT_RANGE fgb+ lgb- fgb+ lgb- ... fgb+ lgb-,"; std::cerr << std::endl << ">> with N the number of runs in the msr-file."; @@ -520,7 +641,7 @@ Bool_t PFitter::CheckCommands() fCmdList.push_back(cmd); } else { std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo; - std::cerr << std::endl << ">> " << it->fLine.Data(); + std::cerr << std::endl << ">> " << line.Data(); std::cerr << std::endl << ">> Syntax: FIT_RANGE RESET | FIT_RANGE | FIT_RANGE .. |"; std::cerr << std::endl << ">> FIT_RANGE fgb+ lgb- | FIT_RANGE fgb+ lgb- fgb+ lgb- ... fgb+ lgb-,"; std::cerr << std::endl << ">> with N the number of runs in the msr-file."; @@ -535,7 +656,7 @@ Bool_t PFitter::CheckCommands() } } else { std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo; - std::cerr << std::endl << ">> " << it->fLine.Data(); + std::cerr << std::endl << ">> " << line.Data(); std::cerr << std::endl << ">> Syntax: FIT_RANGE RESET | FIT_RANGE | FIT_RANGE .. |"; std::cerr << std::endl << ">> FIT_RANGE fgb+ lgb- | FIT_RANGE fgb+ lgb- fgb+ lgb- ... fgb+ lgb-,"; std::cerr << std::endl << ">> with N the number of runs in the msr-file."; @@ -551,14 +672,14 @@ Bool_t PFitter::CheckCommands() delete tokens; tokens = nullptr; } - } else if (it->fLine.Contains("FIX", TString::kIgnoreCase)) { + } else if (line.Contains("FIX", TString::kIgnoreCase)) { // check if the given set of parameters (number or names) is present TObjArray *tokens = nullptr; TObjString *ostr; TString str; UInt_t ival; - tokens = it->fLine.Tokenize(", \t"); + tokens = line.Tokenize(", \t"); for (Int_t i=1; iGetEntries(); i++) { ostr = dynamic_cast(tokens->At(i)); @@ -569,7 +690,7 @@ Bool_t PFitter::CheckCommands() // check that ival is in the parameter list if (ival > fParams.size()) { std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo; - std::cerr << std::endl << ">> " << it->fLine.Data(); + std::cerr << std::endl << ">> " << line.Data(); std::cerr << std::endl << ">> Parameter " << ival << " is out of the Parameter Range [1," << fParams.size() << "]"; std::cerr << std::endl; fIsValid = false; @@ -590,7 +711,7 @@ Bool_t PFitter::CheckCommands() } if (!found) { std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo; - std::cerr << std::endl << ">> " << it->fLine.Data(); + std::cerr << std::endl << ">> " << line.Data(); std::cerr << std::endl << ">> Parameter '" << str.Data() << "' is NOT present as a parameter name"; std::cerr << std::endl; fIsValid = false; @@ -612,46 +733,46 @@ Bool_t PFitter::CheckCommands() cmd.first = PMN_FIX; cmd.second = cmdLineNo; fCmdList.push_back(cmd); - } else if (it->fLine.Contains("HESSE", TString::kIgnoreCase)) { + } else if (line.Contains("HESSE", TString::kIgnoreCase)) { fIsScanOnly = false; cmd.first = PMN_HESSE; cmd.second = cmdLineNo; fCmdList.push_back(cmd); - } else if (it->fLine.Contains("MACHINE_PRECISION", TString::kIgnoreCase)) { + } else if (line.Contains("MACHINE_PRECISION", TString::kIgnoreCase)) { cmd.first = PMN_MACHINE_PRECISION; cmd.second = cmdLineNo; fCmdList.push_back(cmd); - } else if (it->fLine.Contains("MIGRAD", TString::kIgnoreCase)) { + } else if (line.Contains("MIGRAD", TString::kIgnoreCase)) { fIsScanOnly = false; cmd.first = PMN_MIGRAD; cmd.second = cmdLineNo; fCmdList.push_back(cmd); - } else if (it->fLine.Contains("MINIMIZE", TString::kIgnoreCase)) { + } else if (line.Contains("MINIMIZE", TString::kIgnoreCase)) { fIsScanOnly = false; cmd.first = PMN_MINIMIZE; cmd.second = cmdLineNo; fCmdList.push_back(cmd); - } else if (it->fLine.Contains("MINOS", TString::kIgnoreCase)) { + } else if (line.Contains("MINOS", TString::kIgnoreCase)) { fIsScanOnly = false; cmd.first = PMN_MINOS; cmd.second = cmdLineNo; fCmdList.push_back(cmd); - } else if (it->fLine.Contains("MNPLOT", TString::kIgnoreCase)) { + } else if (line.Contains("MNPLOT", TString::kIgnoreCase)) { cmd.first = PMN_PLOT; cmd.second = cmdLineNo; fCmdList.push_back(cmd); - } else if (it->fLine.Contains("PRINT_LEVEL", TString::kIgnoreCase)) { + } else if (line.Contains("PRINT_LEVEL", TString::kIgnoreCase)) { cmd.first = PMN_PRINT; cmd.second = cmdLineNo; fCmdList.push_back(cmd); - } else if (it->fLine.Contains("RELEASE", TString::kIgnoreCase)) { + } else if (line.Contains("RELEASE", TString::kIgnoreCase)) { // check if the given set of parameters (number or names) is present TObjArray *tokens = nullptr; TObjString *ostr; TString str; UInt_t ival; - tokens = it->fLine.Tokenize(", \t"); + tokens = line.Tokenize(", \t"); for (Int_t i=1; iGetEntries(); i++) { ostr = dynamic_cast(tokens->At(i)); @@ -662,7 +783,7 @@ Bool_t PFitter::CheckCommands() // check that ival is in the parameter list if (ival > fParams.size()) { std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo; - std::cerr << std::endl << ">> " << it->fLine.Data(); + std::cerr << std::endl << ">> " << line.Data(); std::cerr << std::endl << ">> Parameter " << ival << " is out of the Parameter Range [1," << fParams.size() << "]"; std::cerr << std::endl; fIsValid = false; @@ -683,7 +804,7 @@ Bool_t PFitter::CheckCommands() } if (!found) { std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo; - std::cerr << std::endl << ">> " << it->fLine.Data(); + std::cerr << std::endl << ">> " << line.Data(); std::cerr << std::endl << ">> Parameter '" << str.Data() << "' is NOT present as a parameter name"; std::cerr << std::endl; fIsValid = false; @@ -703,15 +824,15 @@ Bool_t PFitter::CheckCommands() cmd.first = PMN_RELEASE; cmd.second = cmdLineNo; fCmdList.push_back(cmd); - } else if (it->fLine.Contains("RESTORE", TString::kIgnoreCase)) { + } else if (line.Contains("RESTORE", TString::kIgnoreCase)) { cmd.first = PMN_RESTORE; cmd.second = cmdLineNo; fCmdList.push_back(cmd); - } else if (it->fLine.Contains("SAVE", TString::kIgnoreCase)) { + } else if (line.Contains("SAVE", TString::kIgnoreCase)) { cmd.first = PMN_SAVE; cmd.second = cmdLineNo; fCmdList.push_back(cmd); - } else if (it->fLine.Contains("SCAN", TString::kIgnoreCase)) { + } else if (line.Contains("SCAN", TString::kIgnoreCase)) { cmd.first = PMN_SCAN; cmd.second = cmdLineNo; fCmdList.push_back(cmd); @@ -721,7 +842,7 @@ Bool_t PFitter::CheckCommands() TString str; UInt_t ival; - tokens = it->fLine.Tokenize(", \t"); + tokens = line.Tokenize(", \t"); for (Int_t i=0; iGetEntries(); i++) { ostr = dynamic_cast(tokens->At(i)); @@ -730,7 +851,7 @@ Bool_t PFitter::CheckCommands() // check that token is a UInt_t if (!str.IsDigit()) { std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo; - std::cerr << std::endl << ">> " << it->fLine.Data(); + std::cerr << std::endl << ">> " << line.Data(); std::cerr << std::endl << ">> parameter number is not number!"; std::cerr << std::endl << ">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]"; std::cerr << std::endl; @@ -745,7 +866,7 @@ Bool_t PFitter::CheckCommands() // check that parameter is within range if ((ival < 1) || (ival > fParams.size())) { std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo; - std::cerr << std::endl << ">> " << it->fLine.Data(); + std::cerr << std::endl << ">> " << line.Data(); std::cerr << std::endl << ">> parameter number is out of range [1," << fParams.size() << "]!"; std::cerr << std::endl << ">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]"; std::cerr << std::endl; @@ -765,7 +886,7 @@ Bool_t PFitter::CheckCommands() // check that token is a UInt_t if (!str.IsDigit()) { std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo; - std::cerr << std::endl << ">> " << it->fLine.Data(); + std::cerr << std::endl << ">> " << line.Data(); std::cerr << std::endl << ">> number of points is not number!"; std::cerr << std::endl << ">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]"; std::cerr << std::endl; @@ -779,7 +900,7 @@ Bool_t PFitter::CheckCommands() ival = str.Atoi(); if ((ival < 1) || (ival > 100)) { std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo; - std::cerr << std::endl << ">> " << it->fLine.Data(); + std::cerr << std::endl << ">> " << line.Data(); std::cerr << std::endl << ">> number of scan points is out of range [1,100]!"; std::cerr << std::endl << ">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]"; std::cerr << std::endl; @@ -797,7 +918,7 @@ Bool_t PFitter::CheckCommands() // check that token is a Double_t if (!str.IsFloat()) { std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo; - std::cerr << std::endl << ">> " << it->fLine.Data(); + std::cerr << std::endl << ">> " << line.Data(); std::cerr << std::endl << ">> low is not a floating point number!"; std::cerr << std::endl << ">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]"; std::cerr << std::endl; @@ -815,7 +936,7 @@ Bool_t PFitter::CheckCommands() // check that token is a Double_t if (!str.IsFloat()) { std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo; - std::cerr << std::endl << ">> " << it->fLine.Data(); + std::cerr << std::endl << ">> " << line.Data(); std::cerr << std::endl << ">> high is not a floating point number!"; std::cerr << std::endl << ">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]"; std::cerr << std::endl; @@ -834,16 +955,16 @@ Bool_t PFitter::CheckCommands() delete tokens; tokens = nullptr; } - } else if (it->fLine.Contains("SIMPLEX", TString::kIgnoreCase)) { + } else if (line.Contains("SIMPLEX", TString::kIgnoreCase)) { cmd.first = PMN_SIMPLEX; cmd.second = cmdLineNo; fCmdList.push_back(cmd); - } else if (it->fLine.Contains("STRATEGY", TString::kIgnoreCase)) { + } else if (line.Contains("STRATEGY", TString::kIgnoreCase)) { TObjArray *tokens = nullptr; TObjString *ostr; TString str; - tokens = it->fLine.Tokenize(" \t"); + tokens = line.Tokenize(" \t"); if (tokens->GetEntries() == 2) { ostr = dynamic_cast(tokens->At(1)); str = ostr->GetString(); @@ -866,19 +987,89 @@ Bool_t PFitter::CheckCommands() delete tokens; tokens = nullptr; } - } else if (it->fLine.Contains("USER_COVARIANCE", TString::kIgnoreCase)) { + } else if (line.Contains("USER_COVARIANCE", TString::kIgnoreCase)) { cmd.first = PMN_USER_COVARIANCE; cmd.second = cmdLineNo; fCmdList.push_back(cmd); - } else if (it->fLine.Contains("USER_PARAM_STATE", TString::kIgnoreCase)) { + } else if (line.Contains("USER_PARAM_STATE", TString::kIgnoreCase)) { cmd.first = PMN_USER_PARAM_STATE; cmd.second = cmdLineNo; fCmdList.push_back(cmd); + } else if (line.Contains("SECTOR", TString::kIgnoreCase)) { + cmd.first = PMN_SECTOR; + cmd.second = cmdLineNo; + fCmdList.push_back(cmd); + + // check if the given sector arguments are valid time stamps, i.e. doubles and value < lgb time stamp + TObjArray *tokens = nullptr; + TObjString *ostr; + TString str; + + tokens = line.Tokenize(" ,\t"); + + if (tokens->GetEntries() == 1) { // no sector time stamps given -> issue an error + std::cerr << std::endl << ">> PFitter::CheckCommands(): **FATAL ERROR** in line " << it->fLineNo; + std::cerr << std::endl << ">> " << line.Data(); + std::cerr << std::endl << ">> At least one sector time stamp is expected."; + std::cerr << std::endl << ">> Will stop ..."; + std::cerr << std::endl; + // cleanup + if (tokens) { + delete tokens; + tokens = nullptr; + } + fIsValid = false; + break; + } + + Double_t dval; + for (Int_t i=1; iGetEntries(); i++) { + ostr = dynamic_cast(tokens->At(i)); + str = ostr->GetString(); + if (str.IsFloat()) { + dval = str.Atof(); + // check that the sector time stamp is smaller than all lgb time stamps + for (UInt_t j=0; j fOriginalFitRange[j].second) { + std::cerr << std::endl << ">> PFitter::CheckCommands(): **FATAL ERROR** in line " << it->fLineNo; + std::cerr << std::endl << ">> " << line.Data(); + std::cerr << std::endl << ">> The sector time stamp " << dval << " is > as the lgb time stamp (" << fOriginalFitRange[j].second << ") of run " << j << "."; + std::cerr << std::endl << ">> Will stop ..."; + std::cerr << std::endl; + // cleanup + if (tokens) { + delete tokens; + tokens = nullptr; + } + fIsValid = false; + return fIsValid; + } + } + } else { // sector element is NOT a float + std::cerr << std::endl << ">> PFitter::CheckCommands(): **FATAL ERROR** in line " << it->fLineNo; + std::cerr << std::endl << ">> " << line.Data(); + std::cerr << std::endl << ">> The sector time stamp '" << str << "' is not a number."; + std::cerr << std::endl << ">> Will stop ..."; + std::cerr << std::endl; + // cleanup + if (tokens) { + delete tokens; + tokens = nullptr; + } + fIsValid = false; + break; + } + } + + if (tokens) { + delete tokens; + tokens = nullptr; + } } else { // unkown command - std::cerr << std::endl << ">> PFitter::CheckCommands(): **FATAL ERROR** in line " << it->fLineNo << " an unkown command is found:"; - std::cerr << std::endl << ">> " << it->fLine.Data(); - std::cerr << std::endl << ">> Will stop ..."; - std::cerr << std::endl; + std::cerr << std::endl << ">> PFitter::CheckCommands(): **FATAL ERROR** in line " << it->fLineNo << " an unkown command is found:"; + std::cerr << std::endl << ">> " << line.Data(); + std::cerr << std::endl << ">> Will stop ..."; + std::cerr << std::endl; fIsValid = false; break; } @@ -890,17 +1081,17 @@ Bool_t PFitter::CheckCommands() Bool_t releaseFlag = false; Bool_t minimizerFlag = false; for (it = fCmdLines.begin(); it != fCmdLines.end(); ++it) { - if (it->fLine.Contains("FIX", TString::kIgnoreCase)) + if (line.Contains("FIX", TString::kIgnoreCase)) fixFlag = true; - else if (it->fLine.Contains("RELEASE", TString::kIgnoreCase) || - it->fLine.Contains("RESTORE", TString::kIgnoreCase)) + else if (line.Contains("RELEASE", TString::kIgnoreCase) || + line.Contains("RESTORE", TString::kIgnoreCase)) releaseFlag = true; - else if (it->fLine.Contains("MINIMIZE", TString::kIgnoreCase) || - it->fLine.Contains("MIGRAD", TString::kIgnoreCase) || - it->fLine.Contains("SIMPLEX", TString::kIgnoreCase)) { + else if (line.Contains("MINIMIZE", TString::kIgnoreCase) || + line.Contains("MIGRAD", TString::kIgnoreCase) || + line.Contains("SIMPLEX", TString::kIgnoreCase)) { if (releaseFlag) minimizerFlag = true; - } else if (it->fLine.Contains("MINOS", TString::kIgnoreCase)) { + } else if (line.Contains("MINOS", TString::kIgnoreCase)) { if (fixFlag && releaseFlag && !minimizerFlag) { std::cerr << std::endl << ">> PFitter::CheckCommands(): **WARNING** RELEASE/RESTORE command present"; std::cerr << std::endl << ">> without minimizer command (MINIMIZE/MIGRAD/SIMPLEX) between"; @@ -1969,6 +2160,21 @@ Bool_t PFitter::ExecuteSimplex() return true; } +//-------------------------------------------------------------------------- +// ExecuteSector +//-------------------------------------------------------------------------- +/** + *

Collects all the necessary chisq/maxLH for the given sectors. + * + * return: if the sector command was successful, otherwise return flase. + */ +Bool_t PFitter::ExecuteSector() +{ + // NOT YET IMPLEMENTED //as35 + + return true; +} + //-------------------------------------------------------------------------- // MilliTime //-------------------------------------------------------------------------- diff --git a/src/include/PFitter.h b/src/include/PFitter.h index d6e12f7b..ea015686 100644 --- a/src/include/PFitter.h +++ b/src/include/PFitter.h @@ -58,7 +58,44 @@ #define PMN_USER_COVARIANCE 17 #define PMN_USER_PARAM_STATE 18 #define PMN_PRINT 19 +#define PMN_SECTOR 20 +//----------------------------------------------------------------------------- +/** + *

The PSectorChisq class is needed to store the chisq/maxLH of a sector. + * A sector is a time window from fgb to fLast (fLast < lgb). It also stores + * these information for each run of the msr-file. + */ +class PSectorChisq +{ + public: + PSectorChisq(UInt_t noOfRuns); + + void SetTimeRange(Double_t first, Double_t last); + void SetChisq(Double_t chisq) { fChisq = chisq; } + void SetChisq(Double_t chisq, UInt_t idx); + void SetNDF(Double_t ndf) { fNDF = ndf; } + void SetNDF(Double_t ndf, UInt_t idx); + + Double_t GetTimeRangeFirst() { return fFirst; } + Double_t GetTimeRangeLast() { return fLast; } + Double_t GetChisq() { return fChisq; } + UInt_t GetNDF() { return fNDF; } + UInt_t GetNoRuns() { return fNoOfRuns; } + Double_t GetChisq(UInt_t idx); + UInt_t GetNDF(UInt_t idx); + + private: + UInt_t fNoOfRuns; ///< number of runs presesent + Double_t fFirst; ///< time stamp for fgb + Double_t fLast; ///< requested time stamp + Double_t fChisq; ///< chisq or maxLH for the sector + UInt_t fNDF; ///< NDF for the sector + std::vector fChisqRun; ///< chisq or maxLH for the sector and run + std::vector fNDFRun; ///< NDF for the sector and run +}; + +//----------------------------------------------------------------------------- /** *

Interface class to minuit2. */ @@ -108,6 +145,8 @@ class PFitter PStringVector fElapsedTime; + std::vector fSector; ///< stores all chisq/maxLH sector information + // commands Bool_t CheckCommands(); Bool_t SetParameters(); @@ -126,6 +165,7 @@ class PFitter Bool_t ExecuteScan(); Bool_t ExecuteSave(Bool_t first); Bool_t ExecuteSimplex(); + Bool_t ExecuteSector(); Double_t MilliTime(); }; From 1fb1753d8238bbdfd9907d1b2146d6c349fa0cce Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Sun, 2 Feb 2020 17:12:08 +0100 Subject: [PATCH 2/5] more work on the SECTOR command. CalcNoOfFitBins() needed to be changed from protected to public. --- src/classes/PFitter.cpp | 240 ++++++++++++++++++++++++++--- src/classes/PRunAsymmetry.cpp | 2 +- src/classes/PRunAsymmetryBNMR.cpp | 2 +- src/classes/PRunAsymmetryRRF.cpp | 2 +- src/classes/PRunListCollection.cpp | 24 ++- src/classes/PRunMuMinus.cpp | 2 +- src/classes/PRunSingleHisto.cpp | 2 +- src/classes/PRunSingleHistoRRF.cpp | 2 +- src/include/PFitter.h | 9 +- src/include/PMsrHandler.h | 2 + src/include/PRunAsymmetry.h | 3 +- src/include/PRunAsymmetryBNMR.h | 3 +- src/include/PRunAsymmetryRRF.h | 3 +- src/include/PRunMuMinus.h | 3 +- src/include/PRunSingleHisto.h | 3 +- src/include/PRunSingleHistoRRF.h | 3 +- 16 files changed, 263 insertions(+), 42 deletions(-) diff --git a/src/classes/PFitter.cpp b/src/classes/PFitter.cpp index 2b301c78..ccba9226 100644 --- a/src/classes/PFitter.cpp +++ b/src/classes/PFitter.cpp @@ -79,30 +79,38 @@ PSectorChisq::PSectorChisq(UInt_t noOfRuns) : fNoOfRuns(noOfRuns) { // init - fFirst = 0.0; fLast = 0.0; fChisq = 0.0; fNDF = 0; + fFirst.resize(fNoOfRuns); fChisqRun.resize(fNoOfRuns); fNDFRun.resize(fNoOfRuns); for (UInt_t i=0; iSet the time range for one sector + *

Set the time of the fgb of a given RUN * * @param first time stamp of the fgb - * @param last time stamp of the requested sector end + * @param idx index of the RUN */ -void PSectorChisq::SetTimeRange(Double_t first, Double_t last) +void PSectorChisq::SetRunFirstTime(Double_t first, UInt_t idx) { - // NOT YET IMPLEMENTED //as35 + if (idx > fNoOfRuns) { + std::cerr << "**WARNING** from PSectorChisq::SetRunFirstTime. It tries to set" << std::endl; + std::cerr << " a fgb time stamp with idx=" << idx << " which is larger than #RUNS=" << fNoOfRuns << "." << std::endl; + std::cerr << " Will ignore it, but you better check what is going on!" << std::endl; + return; + } + + fFirst[idx] = first; } //-------------------------------------------------------------------------- @@ -116,7 +124,14 @@ void PSectorChisq::SetTimeRange(Double_t first, Double_t last) */ void PSectorChisq::SetChisq(Double_t chisq, UInt_t idx) { - // NOT YET IMPLEMENTED //as35 + if (idx > fNoOfRuns) { + std::cerr << "**WARNING** from PSectorChisq::SetChisq. It tries to set" << std::endl; + std::cerr << " a chisq with idx=" << idx << " which is larger than #RUNS=" << fNoOfRuns << "." << std::endl; + std::cerr << " Will ignore it, but you better check what is going on!" << std::endl; + return; + } + + fChisqRun[idx] = chisq; } //-------------------------------------------------------------------------- @@ -130,7 +145,33 @@ void PSectorChisq::SetChisq(Double_t chisq, UInt_t idx) */ void PSectorChisq::SetNDF(Double_t ndf, UInt_t idx) { - // NOT YET IMPLEMENTED //as35 + if (idx > fNoOfRuns) { + std::cerr << "**WARNING** from PSectorChisq::SetNDF. It tries to set" << std::endl; + std::cerr << " a NDF with idx=" << idx << " which is larger than #RUNS=" << fNoOfRuns << "." << std::endl; + std::cerr << " Will ignore it, but you better check what is going on!" << std::endl; + return; + } + + fNDFRun[idx] = ndf; +} + +//-------------------------------------------------------------------------- +// GetTimeRangeFirst +//-------------------------------------------------------------------------- +/** + *

Get the fgb time of RUN with index idx. If idx is out-of-range + * PMUSR_UNDEFINED is returned. + * + * @param idx index of the RUN + * + * return: return the fgb time of RUN with index idx. + */ +Double_t PSectorChisq::GetTimeRangeFirst(UInt_t idx) +{ + if (idx > fNoOfRuns) + return PMUSR_UNDEFINED; + + return fFirst[idx]; } //-------------------------------------------------------------------------- @@ -191,6 +232,8 @@ PFitter::PFitter(PMsrHandler *runInfo, PRunListCollection *runListCollection, Bo fStrategy = 1; // 0=low, 1=default, 2=high + fSectorFlag = false; + fParams = *(runInfo->GetMsrParamList()); fCmdLines = *runInfo->GetMsrCommands(); @@ -276,6 +319,11 @@ Bool_t PFitter::DoFit() // check if only chisq/maxLH shall be calculated once if (fChisqOnly) { + + // check for sector command + if (fSectorFlag) + PrepareSector(); + std::vector param = fMnUserParams.Params(); std::vector error = fMnUserParams.Errors(); Int_t usedParams = 0; @@ -288,12 +336,12 @@ Bool_t PFitter::DoFit() if (fUseChi2) { // calculate expected chisq Double_t totalExpectedChisq = 0.0; - std::vector expectedChisqPerHisto; - fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto); + std::vector expectedChisqPerRun; + fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerRun); // calculate chisq per run - std::vector chisqPerHisto; + std::vector chisqPerRun; for (UInt_t i=0; iGetMsrRunList()->size(); i++) { - chisqPerHisto.push_back(fRunListCollection->GetSingleRunChisq(param, i)); + chisqPerRun.push_back(fRunListCollection->GetSingleRunChisq(param, i)); } std::cout << std::endl << std::endl << ">> chisq = " << val << ", NDF = " << ndf << ", chisq/NDF = " << val/ndf; @@ -301,25 +349,137 @@ Bool_t PFitter::DoFit() if (totalExpectedChisq != 0.0) { std::cout << std::endl << ">> expected chisq = " << totalExpectedChisq << ", NDF = " << ndf << ", expected chisq/NDF = " << totalExpectedChisq/ndf; UInt_t ndf_histo = 0; - for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); if (ndf_histo > 0) - std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq/red.chisq_e) = (" << ndf_histo << "/" << chisqPerHisto[i]/ndf_histo << "/" << expectedChisqPerHisto[i]/ndf_histo << ")"; + std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq/red.chisq_e) = (" << ndf_histo << "/" << chisqPerRun[i]/ndf_histo << "/" << expectedChisqPerRun[i]/ndf_histo << ")"; } - } else if (chisqPerHisto.size() > 0) { // in case expected chisq is not applicable like for asymmetry fits + } else if (chisqPerRun.size() > 0) { // in case expected chisq is not applicable like for asymmetry fits UInt_t ndf_histo = 0; - for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); if (ndf_histo > 0) - std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq) = (" << ndf_histo << "/" << chisqPerHisto[i]/ndf_histo << ")"; + std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq) = (" << ndf_histo << "/" << chisqPerRun[i]/ndf_histo << ")"; } } // clean up - chisqPerHisto.clear(); - expectedChisqPerHisto.clear(); + chisqPerRun.clear(); + expectedChisqPerRun.clear(); + + if (fSectorFlag) { + PDoublePairVector secFitRange; + secFitRange.resize(1); + for (UInt_t k=0; kSetFitRange(secFitRange); + // calculate chisq + val = (*fFitterFcn)(param); + // calculate NDF + ndf = static_cast(fFitterFcn->GetTotalNoOfFittedBins()) - usedParams; + // calculate expected chisq + totalExpectedChisq = 0.0; + fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerRun); + // calculate chisq per run + for (UInt_t i=0; iGetMsrRunList()->size(); i++) { + chisqPerRun.push_back(fRunListCollection->GetSingleRunChisq(param, i)); + } + + std::cout << std::endl; + std::cout << "++++" << std::endl; + std::cout << ">> Sector " << k << ": FitRange: " << secFitRange[0].first << ", " << secFitRange[0].second << std::endl; + std::cout << ">> chisq = " << val << ", NDF = " << ndf << ", chisq/NDF = " << val/ndf; + + if (totalExpectedChisq != 0.0) { + std::cout << std::endl << ">> expected chisq = " << totalExpectedChisq << ", NDF = " << ndf << ", expected chisq/NDF = " << totalExpectedChisq/ndf; + UInt_t ndf_histo = 0; + for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_histo > 0) + std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq/red.chisq_e) = (" << ndf_histo << "/" << chisqPerRun[i]/ndf_histo << "/" << expectedChisqPerRun[i]/ndf_histo << ")"; + } + } else if (chisqPerRun.size() > 0) { // in case expected chisq is not applicable like for asymmetry fits + UInt_t ndf_histo = 0; + for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_histo > 0) + std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq) = (" << ndf_histo << "/" << chisqPerRun[i]/ndf_histo << ")"; + } + } + // clean up + chisqPerRun.clear(); + expectedChisqPerRun.clear(); + } + } } else { // max. log likelihood + // calculate expected maxLH + Double_t totalExpectedMaxLH = 0.0; + std::vector expectedMaxLHPerRun; + fFitterFcn->CalcExpectedChiSquare(param, totalExpectedMaxLH, expectedMaxLHPerRun); + // calculate maxLH per run + std::vector maxLHPerRun; + for (UInt_t i=0; iGetMsrRunList()->size(); i++) { + maxLHPerRun.push_back(fRunListCollection->GetSingleRunMaximumLikelihood(param, i)); + } + std::cout << std::endl << std::endl << ">> maxLH = " << val << ", NDF = " << ndf << ", maxLH/NDF = " << val/ndf; + + if (totalExpectedMaxLH != 0.0) { + std::cout << std::endl << ">> expected maxLH = " << totalExpectedMaxLH << ", NDF = " << ndf << ", expected maxLH/NDF = " << totalExpectedMaxLH/ndf; + UInt_t ndf_histo = 0; + for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_histo > 0) + std::cout << std::endl << ">> run block " << i+1 << ": (NDF/maxLH.chisq/maxLH.chisq_e) = (" << ndf_histo << "/" << maxLHPerRun[i]/ndf_histo << "/" << expectedMaxLHPerRun[i]/ndf_histo << ")"; + } + } + + // clean up + maxLHPerRun.clear(); + expectedMaxLHPerRun.clear(); + + if (fSectorFlag) { + PDoublePairVector secFitRange; + secFitRange.resize(1); + for (UInt_t k=0; kSetFitRange(secFitRange); + // calculate maxLH + val = (*fFitterFcn)(param); + // calculate NDF + ndf = static_cast(fFitterFcn->GetTotalNoOfFittedBins()) - usedParams; + // calculate expected maxLH + totalExpectedMaxLH = 0.0; + fFitterFcn->CalcExpectedChiSquare(param, totalExpectedMaxLH, expectedMaxLHPerRun); + // calculate maxLH per run + for (UInt_t i=0; iGetMsrRunList()->size(); i++) { + maxLHPerRun.push_back(fRunListCollection->GetSingleRunMaximumLikelihood(param, i)); + } + + std::cout << std::endl; + std::cout << "++++" << std::endl; + std::cout << ">> Sector " << k << ": FitRange: " << secFitRange[0].first << ", " << secFitRange[0].second << std::endl; + std::cout << ">> maxLH = " << val << ", NDF = " << ndf << ", maxLH/NDF = " << val/ndf; + + if (totalExpectedMaxLH != 0.0) { + std::cout << std::endl << ">> expected maxLH = " << totalExpectedMaxLH << ", NDF = " << ndf << ", expected maxLH/NDF = " << totalExpectedMaxLH/ndf; + UInt_t ndf_histo = 0; + for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_histo > 0) + std::cout << std::endl << ">> run block " << i+1 << ": (NDF/maxLH.chisq/maxLH.chisq_e) = (" << ndf_histo << "/" << maxLHPerRun[i]/ndf_histo << "/" << expectedMaxLHPerRun[i]/ndf_histo << ")"; + } + } + + // clean up + maxLHPerRun.clear(); + expectedMaxLHPerRun.clear(); + } + } } std::cout << std::endl << std::endl; return true; @@ -396,6 +556,9 @@ Bool_t PFitter::DoFit() case PMN_SCAN: status = ExecuteScan(); break; + case PMN_SECTOR: + // nothing to be done here + break; case PMN_SIMPLEX: status = ExecuteSimplex(); break; @@ -996,6 +1159,7 @@ Bool_t PFitter::CheckCommands() cmd.second = cmdLineNo; fCmdList.push_back(cmd); } else if (line.Contains("SECTOR", TString::kIgnoreCase)) { + fSectorFlag = true; cmd.first = PMN_SECTOR; cmd.second = cmdLineNo; fCmdList.push_back(cmd); @@ -1019,11 +1183,15 @@ Bool_t PFitter::CheckCommands() tokens = nullptr; } fIsValid = false; + fSectorFlag = false; break; } Double_t dval; for (Int_t i=1; iGetEntries(); i++) { + // keep time range of sector + PSectorChisq sec(fRunInfo->GetNoOfRuns()); + // get parse tokens ostr = dynamic_cast(tokens->At(i)); str = ostr->GetString(); if (str.IsFloat()) { @@ -1042,9 +1210,13 @@ Bool_t PFitter::CheckCommands() tokens = nullptr; } fIsValid = false; + fSectorFlag = false; return fIsValid; } + sec.SetRunFirstTime(fOriginalFitRange[j].first, j); // keep fgb time stamp for sector } + sec.SetSectorTime(dval); + fSector.push_back(sec); } else { // sector element is NOT a float std::cerr << std::endl << ">> PFitter::CheckCommands(): **FATAL ERROR** in line " << it->fLineNo; std::cerr << std::endl << ">> " << line.Data(); @@ -1057,6 +1229,7 @@ Bool_t PFitter::CheckCommands() tokens = nullptr; } fIsValid = false; + fSectorFlag = false; break; } } @@ -1216,7 +1389,7 @@ Bool_t PFitter::ExecuteFitRange(UInt_t lineNo) PMsrRunList *runList = fRunInfo->GetMsrRunList(); - // execute command, no error checking needed since this has been already carried out in CheckCommands() + // 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 @@ -2160,11 +2333,36 @@ Bool_t PFitter::ExecuteSimplex() return true; } +//-------------------------------------------------------------------------- +// PrepareSector +//-------------------------------------------------------------------------- +/** + *

Collect all the necessary chisq/maxLH sector information. + */ +void PFitter::PrepareSector() +{ + // NOT YET IMPLEMENTED //as35 + +/* + for (UInt_t i=0; i +++++" << std::endl; + std::cout << "debug> sector " << i << std::endl; + std::cout << "debug> noOfRuns: " << fSector[i].GetNoRuns() << std::endl; + for (UInt_t j=0; j ----" << std::endl; + std::cout << "debug> RUN BLOCK " << j << std::endl; + std::cout << "debug> Sector Time Range: " << fSector[i].GetTimeRangeFirst(j) << ", " << fSector[i].GetTimeRangeLast() << std::endl; + } + } +*/ +} + //-------------------------------------------------------------------------- // ExecuteSector //-------------------------------------------------------------------------- /** - *

Collects all the necessary chisq/maxLH for the given sectors. + *

Write all chisq/maxLH sector information to MINUIT.OUTPUT and dump it + * to stdout. * * return: if the sector command was successful, otherwise return flase. */ diff --git a/src/classes/PRunAsymmetry.cpp b/src/classes/PRunAsymmetry.cpp index 2616ea1b..eb6edba7 100644 --- a/src/classes/PRunAsymmetry.cpp +++ b/src/classes/PRunAsymmetry.cpp @@ -378,7 +378,7 @@ void PRunAsymmetry::SetFitRangeBin(const TString fitRange) } //-------------------------------------------------------------------------- -// CalcNoOfFitBins (protected) +// CalcNoOfFitBins (public) //-------------------------------------------------------------------------- /** *

Calculate the number of fitted bins for the current fit range. diff --git a/src/classes/PRunAsymmetryBNMR.cpp b/src/classes/PRunAsymmetryBNMR.cpp index 0235c989..9b95c074 100644 --- a/src/classes/PRunAsymmetryBNMR.cpp +++ b/src/classes/PRunAsymmetryBNMR.cpp @@ -400,7 +400,7 @@ void PRunAsymmetryBNMR::SetFitRangeBin(const TString fitRange) } //-------------------------------------------------------------------------- -// CalcNoOfFitBins (protected) +// CalcNoOfFitBins (public) //-------------------------------------------------------------------------- /** *

Calculate the number of fitted bins for the current fit range. diff --git a/src/classes/PRunAsymmetryRRF.cpp b/src/classes/PRunAsymmetryRRF.cpp index 64d2fdc2..0c39f0e7 100644 --- a/src/classes/PRunAsymmetryRRF.cpp +++ b/src/classes/PRunAsymmetryRRF.cpp @@ -370,7 +370,7 @@ void PRunAsymmetryRRF::SetFitRangeBin(const TString fitRange) } //-------------------------------------------------------------------------- -// CalcNoOfFitBins (protected) +// CalcNoOfFitBins (public) //-------------------------------------------------------------------------- /** *

Calculate the number of fitted bins for the current fit range. diff --git a/src/classes/PRunListCollection.cpp b/src/classes/PRunListCollection.cpp index 103912ae..e456ccd2 100644 --- a/src/classes/PRunListCollection.cpp +++ b/src/classes/PRunListCollection.cpp @@ -206,18 +206,30 @@ void PRunListCollection::SetFitRange(const TString fitRange) */ void PRunListCollection::SetFitRange(const PDoublePairVector fitRange) { - for (UInt_t i=0; iSetFitRange(fitRange); - for (UInt_t i=0; iCalcNoOfFitBins(); // needed to update fStartTimeBin, fEndTimeBin + } + for (UInt_t i=0; iSetFitRange(fitRange); - for (UInt_t i=0; iCalcNoOfFitBins(); // needed to update fStartTimeBin, fEndTimeBin + } + for (UInt_t i=0; iSetFitRange(fitRange); - for (UInt_t i=0; iCalcNoOfFitBins(); // needed to update fStartTimeBin, fEndTimeBin + } + for (UInt_t i=0; iSetFitRange(fitRange); - for (UInt_t i=0; iCalcNoOfFitBins(); // needed to update fStartTimeBin, fEndTimeBin + } + for (UInt_t i=0; iSetFitRange(fitRange); - for (UInt_t i=0; iCalcNoOfFitBins(); // needed to update fStartTimeBin, fEndTimeBin + } + for (UInt_t i=0; iSetFitRange(fitRange); + fRunMuMinusList[i]->CalcNoOfFitBins(); // needed to update fStartTimeBin, fEndTimeBin + } for (UInt_t i=0; iSetFitRange(fitRange); } diff --git a/src/classes/PRunMuMinus.cpp b/src/classes/PRunMuMinus.cpp index 83101a8e..38240625 100644 --- a/src/classes/PRunMuMinus.cpp +++ b/src/classes/PRunMuMinus.cpp @@ -382,7 +382,7 @@ void PRunMuMinus::SetFitRangeBin(const TString fitRange) } //-------------------------------------------------------------------------- -// CalcNoOfFitBins (private) +// CalcNoOfFitBins (public) //-------------------------------------------------------------------------- /** *

Calculate the number of fitted bins for the current fit range. diff --git a/src/classes/PRunSingleHisto.cpp b/src/classes/PRunSingleHisto.cpp index 4931703f..b3d15f09 100644 --- a/src/classes/PRunSingleHisto.cpp +++ b/src/classes/PRunSingleHisto.cpp @@ -659,7 +659,7 @@ void PRunSingleHisto::SetFitRangeBin(const TString fitRange) } //-------------------------------------------------------------------------- -// CalcNoOfFitBins (protected) +// CalcNoOfFitBins (public) //-------------------------------------------------------------------------- /** *

Calculate the number of fitted bins for the current fit range. diff --git a/src/classes/PRunSingleHistoRRF.cpp b/src/classes/PRunSingleHistoRRF.cpp index ec9ab895..e13c1c8d 100644 --- a/src/classes/PRunSingleHistoRRF.cpp +++ b/src/classes/PRunSingleHistoRRF.cpp @@ -388,7 +388,7 @@ void PRunSingleHistoRRF::SetFitRangeBin(const TString fitRange) } //-------------------------------------------------------------------------- -// CalcNoOfFitBins (protected) +// CalcNoOfFitBins (public) //-------------------------------------------------------------------------- /** *

Calculate the number of fitted bins for the current fit range. diff --git a/src/include/PFitter.h b/src/include/PFitter.h index ea015686..56cd6176 100644 --- a/src/include/PFitter.h +++ b/src/include/PFitter.h @@ -71,13 +71,14 @@ class PSectorChisq public: PSectorChisq(UInt_t noOfRuns); - void SetTimeRange(Double_t first, Double_t last); + void SetRunFirstTime(Double_t first, UInt_t idx); + void SetSectorTime(Double_t last) { fLast = last; } void SetChisq(Double_t chisq) { fChisq = chisq; } void SetChisq(Double_t chisq, UInt_t idx); void SetNDF(Double_t ndf) { fNDF = ndf; } void SetNDF(Double_t ndf, UInt_t idx); - Double_t GetTimeRangeFirst() { return fFirst; } + Double_t GetTimeRangeFirst(UInt_t idx); Double_t GetTimeRangeLast() { return fLast; } Double_t GetChisq() { return fChisq; } UInt_t GetNDF() { return fNDF; } @@ -87,10 +88,10 @@ class PSectorChisq private: UInt_t fNoOfRuns; ///< number of runs presesent - Double_t fFirst; ///< time stamp for fgb Double_t fLast; ///< requested time stamp Double_t fChisq; ///< chisq or maxLH for the sector UInt_t fNDF; ///< NDF for the sector + std::vector fFirst; ///< time stamp for fgb for a given run std::vector fChisqRun; ///< chisq or maxLH for the sector and run std::vector fNDFRun; ///< NDF for the sector and run }; @@ -145,6 +146,7 @@ class PFitter PStringVector fElapsedTime; + Bool_t fSectorFlag; ///< sector command present flag std::vector fSector; ///< stores all chisq/maxLH sector information // commands @@ -165,6 +167,7 @@ class PFitter Bool_t ExecuteScan(); Bool_t ExecuteSave(Bool_t first); Bool_t ExecuteSimplex(); + void PrepareSector(); Bool_t ExecuteSector(); Double_t MilliTime(); diff --git a/src/include/PMsrHandler.h b/src/include/PMsrHandler.h index a283b3a1..e7824fb6 100644 --- a/src/include/PMsrHandler.h +++ b/src/include/PMsrHandler.h @@ -66,6 +66,8 @@ class PMsrHandler virtual TString* GetMsrFileDirectoryPath() { return &fMsrFileDirectoryPath; } + virtual UInt_t GetNoOfRuns() { return fRuns.size(); } + virtual UInt_t GetNoOfParams() { return fParam.size(); } virtual const TString& GetFileName() const { return fFileName; } diff --git a/src/include/PRunAsymmetry.h b/src/include/PRunAsymmetry.h index 2d028cbb..798e26da 100644 --- a/src/include/PRunAsymmetry.h +++ b/src/include/PRunAsymmetry.h @@ -56,8 +56,9 @@ class PRunAsymmetry : public PRunBase virtual Int_t GetEndTimeBin() { return fEndTimeBin; } virtual Int_t GetPacking() { return fPacking; } - protected: virtual void CalcNoOfFitBins(); + + protected: virtual Bool_t PrepareData(); virtual Bool_t PrepareFitData(); virtual Bool_t PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]); diff --git a/src/include/PRunAsymmetryBNMR.h b/src/include/PRunAsymmetryBNMR.h index 688c84c8..bc827273 100644 --- a/src/include/PRunAsymmetryBNMR.h +++ b/src/include/PRunAsymmetryBNMR.h @@ -57,8 +57,9 @@ class PRunAsymmetryBNMR : public PRunBase virtual Int_t GetEndTimeBin() { return fEndTimeBin; } virtual Int_t GetPacking() { return fPacking; } - protected: virtual void CalcNoOfFitBins(); + + protected: virtual Bool_t PrepareData(); virtual Bool_t PrepareFitData(); virtual Bool_t PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]); diff --git a/src/include/PRunAsymmetryRRF.h b/src/include/PRunAsymmetryRRF.h index 78971e76..4aef2598 100644 --- a/src/include/PRunAsymmetryRRF.h +++ b/src/include/PRunAsymmetryRRF.h @@ -55,8 +55,9 @@ class PRunAsymmetryRRF : public PRunBase virtual Int_t GetStartTimeBin() { return fStartTimeBin; } virtual Int_t GetEndTimeBin() { return fEndTimeBin; } - protected: virtual void CalcNoOfFitBins(); + + protected: virtual Bool_t PrepareData(); virtual Bool_t PrepareFitData(); virtual Bool_t PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]); diff --git a/src/include/PRunMuMinus.h b/src/include/PRunMuMinus.h index 8c9a0e20..f95056d1 100644 --- a/src/include/PRunMuMinus.h +++ b/src/include/PRunMuMinus.h @@ -55,8 +55,9 @@ class PRunMuMinus : public PRunBase virtual Int_t GetEndTimeBin() { return fEndTimeBin; } virtual Int_t GetPacking() { return fPacking; } - protected: virtual void CalcNoOfFitBins(); + + protected: virtual Bool_t PrepareData(); virtual Bool_t PrepareFitData(PRawRunData* runData, const UInt_t histoNo); virtual Bool_t PrepareRawViewData(PRawRunData* runData, const UInt_t histoNo); diff --git a/src/include/PRunSingleHisto.h b/src/include/PRunSingleHisto.h index 6b3a02d1..2152cae6 100644 --- a/src/include/PRunSingleHisto.h +++ b/src/include/PRunSingleHisto.h @@ -59,8 +59,9 @@ class PRunSingleHisto : public PRunBase virtual Int_t GetPacking() { return fPacking; } virtual Bool_t GetScaleN0AndBkg() { return fScaleN0AndBkg; } - protected: virtual void CalcNoOfFitBins(); + + protected: virtual Bool_t PrepareData(); virtual Bool_t PrepareFitData(PRawRunData* runData, const UInt_t histoNo); virtual Bool_t PrepareRawViewData(PRawRunData* runData, const UInt_t histoNo); diff --git a/src/include/PRunSingleHistoRRF.h b/src/include/PRunSingleHistoRRF.h index 5631bfb4..8f92dcfe 100644 --- a/src/include/PRunSingleHistoRRF.h +++ b/src/include/PRunSingleHistoRRF.h @@ -54,8 +54,9 @@ class PRunSingleHistoRRF : public PRunBase virtual Int_t GetStartTimeBin() { return fStartTimeBin; } virtual Int_t GetEndTimeBin() { return fEndTimeBin; } - protected: virtual void CalcNoOfFitBins(); + + protected: virtual Bool_t PrepareData(); virtual Bool_t PrepareFitData(PRawRunData* runData, const UInt_t histoNo); virtual Bool_t PrepareViewData(PRawRunData* runData, const UInt_t histoNo); From fe7a1b79204f0661aaddfd3155325eb603e96c69 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Mon, 3 Feb 2020 20:54:23 +0100 Subject: [PATCH 3/5] first full implementation of the sector command. --- src/classes/PFitter.cpp | 355 +++++++++++++++++++++++------ src/classes/PFitterFcn.cpp | 4 +- src/classes/PRunListCollection.cpp | 69 +++--- src/classes/PRunSingleHisto.cpp | 14 +- src/include/PFitter.h | 26 ++- src/include/PRunListCollection.h | 4 +- 6 files changed, 345 insertions(+), 127 deletions(-) diff --git a/src/classes/PFitter.cpp b/src/classes/PFitter.cpp index ccba9226..9841a36f 100644 --- a/src/classes/PFitter.cpp +++ b/src/classes/PFitter.cpp @@ -81,13 +81,16 @@ PSectorChisq::PSectorChisq(UInt_t noOfRuns) : fNoOfRuns(noOfRuns) // init fLast = 0.0; fChisq = 0.0; + fExpectedChisq = 0.0; fNDF = 0; fFirst.resize(fNoOfRuns); fChisqRun.resize(fNoOfRuns); + fExpectedChisqRun.resize(fNoOfRuns); fNDFRun.resize(fNoOfRuns); for (UInt_t i=0; iSet the expected chisq/maxLH for a given detector with index idx. + * + * @param chisq expected chisq/maxLH to be set + * @param idx index of the run to be set + */ +void PSectorChisq::SetExpectedChisq(Double_t chisq, UInt_t idx) +{ + if (idx > fNoOfRuns) { + std::cerr << "**WARNING** from PSectorChisq::SetExpectedChisq. It tries to set" << std::endl; + std::cerr << " a chisq with idx=" << idx << " which is larger than #RUNS=" << fNoOfRuns << "." << std::endl; + std::cerr << " Will ignore it, but you better check what is going on!" << std::endl; + return; + } + + fExpectedChisqRun[idx] = chisq; +} + //-------------------------------------------------------------------------- // SetNDF //-------------------------------------------------------------------------- @@ -143,7 +167,7 @@ void PSectorChisq::SetChisq(Double_t chisq, UInt_t idx) * @param ndf to be set * @param idx index of the run to be set */ -void PSectorChisq::SetNDF(Double_t ndf, UInt_t idx) +void PSectorChisq::SetNDF(UInt_t ndf, UInt_t idx) { if (idx > fNoOfRuns) { std::cerr << "**WARNING** from PSectorChisq::SetNDF. It tries to set" << std::endl; @@ -178,7 +202,7 @@ Double_t PSectorChisq::GetTimeRangeFirst(UInt_t idx) // GetChisq //-------------------------------------------------------------------------- /** - *

chisq/maxLH of the run with index idx + *

Get chisq/maxLH of the run with index idx * * @param idx index of the run * @@ -193,10 +217,28 @@ Double_t PSectorChisq::GetChisq(UInt_t idx) } //-------------------------------------------------------------------------- -// GetChisq +// GetExpectedChisq //-------------------------------------------------------------------------- /** - *

NDF of the run with index idx + *

Get expected chisq/maxLH of the run with index idx + * + * @param idx index of the run + * + * return: chisq/maxLH of the requested run or -1.0 if the idx is out-of-range. + */ +Double_t PSectorChisq::GetExpectedChisq(UInt_t idx) +{ + if (idx >= fNoOfRuns) + return -1.0; + + return fExpectedChisqRun[idx]; +} + +//-------------------------------------------------------------------------- +// GetNDF +//-------------------------------------------------------------------------- +/** + *

Get NDF of the run with index idx * * @param idx index of the run * @@ -319,11 +361,6 @@ Bool_t PFitter::DoFit() // check if only chisq/maxLH shall be calculated once if (fChisqOnly) { - - // check for sector command - if (fSectorFlag) - PrepareSector(); - std::vector param = fMnUserParams.Params(); std::vector error = fMnUserParams.Errors(); Int_t usedParams = 0; @@ -331,12 +368,12 @@ Bool_t PFitter::DoFit() if (error[i] != 0.0) usedParams++; } - Int_t ndf = static_cast(fFitterFcn->GetTotalNoOfFittedBins()) - usedParams; + UInt_t ndf = static_cast(fFitterFcn->GetTotalNoOfFittedBins()) - usedParams; Double_t val = (*fFitterFcn)(param); if (fUseChi2) { // calculate expected chisq Double_t totalExpectedChisq = 0.0; - std::vector expectedChisqPerRun; + PDoubleVector expectedChisqPerRun; fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerRun); // calculate chisq per run std::vector chisqPerRun; @@ -348,18 +385,18 @@ Bool_t PFitter::DoFit() if (totalExpectedChisq != 0.0) { std::cout << std::endl << ">> expected chisq = " << totalExpectedChisq << ", NDF = " << ndf << ", expected chisq/NDF = " << totalExpectedChisq/ndf; - UInt_t ndf_histo = 0; + UInt_t ndf_run = 0; for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); - if (ndf_histo > 0) - std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq/red.chisq_e) = (" << ndf_histo << "/" << chisqPerRun[i]/ndf_histo << "/" << expectedChisqPerRun[i]/ndf_histo << ")"; + ndf_run = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_run > 0) + std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq/red.chisq_e) = (" << ndf_run << "/" << chisqPerRun[i]/ndf_run << "/" << expectedChisqPerRun[i]/ndf_run << ")"; } } else if (chisqPerRun.size() > 0) { // in case expected chisq is not applicable like for asymmetry fits - UInt_t ndf_histo = 0; + UInt_t ndf_run = 0; for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); - if (ndf_histo > 0) - std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq) = (" << ndf_histo << "/" << chisqPerRun[i]/ndf_histo << ")"; + ndf_run = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_run > 0) + std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq) = (" << ndf_run << "/" << chisqPerRun[i]/ndf_run << ")"; } } @@ -378,7 +415,7 @@ Bool_t PFitter::DoFit() // calculate chisq val = (*fFitterFcn)(param); // calculate NDF - ndf = static_cast(fFitterFcn->GetTotalNoOfFittedBins()) - usedParams; + ndf = static_cast(fFitterFcn->GetTotalNoOfFittedBins()) - usedParams; // calculate expected chisq totalExpectedChisq = 0.0; fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerRun); @@ -389,23 +426,23 @@ Bool_t PFitter::DoFit() std::cout << std::endl; std::cout << "++++" << std::endl; - std::cout << ">> Sector " << k << ": FitRange: " << secFitRange[0].first << ", " << secFitRange[0].second << std::endl; + std::cout << ">> Sector " << k+1 << ": FitRange: " << secFitRange[0].first << ", " << secFitRange[0].second << std::endl; std::cout << ">> chisq = " << val << ", NDF = " << ndf << ", chisq/NDF = " << val/ndf; if (totalExpectedChisq != 0.0) { std::cout << std::endl << ">> expected chisq = " << totalExpectedChisq << ", NDF = " << ndf << ", expected chisq/NDF = " << totalExpectedChisq/ndf; - UInt_t ndf_histo = 0; + UInt_t ndf_run = 0; for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); - if (ndf_histo > 0) - std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq/red.chisq_e) = (" << ndf_histo << "/" << chisqPerRun[i]/ndf_histo << "/" << expectedChisqPerRun[i]/ndf_histo << ")"; + ndf_run = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_run > 0) + std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq/red.chisq_e) = (" << ndf_run << "/" << chisqPerRun[i]/ndf_run << "/" << expectedChisqPerRun[i]/ndf_run << ")"; } } else if (chisqPerRun.size() > 0) { // in case expected chisq is not applicable like for asymmetry fits - UInt_t ndf_histo = 0; + UInt_t ndf_run = 0; for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); - if (ndf_histo > 0) - std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq) = (" << ndf_histo << "/" << chisqPerRun[i]/ndf_histo << ")"; + ndf_run = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_run > 0) + std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq) = (" << ndf_run << "/" << chisqPerRun[i]/ndf_run << ")"; } } // clean up @@ -428,11 +465,11 @@ Bool_t PFitter::DoFit() if (totalExpectedMaxLH != 0.0) { std::cout << std::endl << ">> expected maxLH = " << totalExpectedMaxLH << ", NDF = " << ndf << ", expected maxLH/NDF = " << totalExpectedMaxLH/ndf; - UInt_t ndf_histo = 0; + UInt_t ndf_run = 0; for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); - if (ndf_histo > 0) - std::cout << std::endl << ">> run block " << i+1 << ": (NDF/maxLH.chisq/maxLH.chisq_e) = (" << ndf_histo << "/" << maxLHPerRun[i]/ndf_histo << "/" << expectedMaxLHPerRun[i]/ndf_histo << ")"; + ndf_run = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_run > 0) + std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.maxLH/red.maxLH_e) = (" << ndf_run << "/" << maxLHPerRun[i]/ndf_run << "/" << expectedMaxLHPerRun[i]/ndf_run << ")"; } } @@ -462,16 +499,16 @@ Bool_t PFitter::DoFit() std::cout << std::endl; std::cout << "++++" << std::endl; - std::cout << ">> Sector " << k << ": FitRange: " << secFitRange[0].first << ", " << secFitRange[0].second << std::endl; + std::cout << ">> Sector " << k+1 << ": FitRange: " << secFitRange[0].first << ", " << secFitRange[0].second << std::endl; std::cout << ">> maxLH = " << val << ", NDF = " << ndf << ", maxLH/NDF = " << val/ndf; if (totalExpectedMaxLH != 0.0) { std::cout << std::endl << ">> expected maxLH = " << totalExpectedMaxLH << ", NDF = " << ndf << ", expected maxLH/NDF = " << totalExpectedMaxLH/ndf; - UInt_t ndf_histo = 0; + UInt_t ndf_run = 0; for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); - if (ndf_histo > 0) - std::cout << std::endl << ">> run block " << i+1 << ": (NDF/maxLH.chisq/maxLH.chisq_e) = (" << ndf_histo << "/" << maxLHPerRun[i]/ndf_histo << "/" << expectedMaxLHPerRun[i]/ndf_histo << ")"; + ndf_run = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_run > 0) + std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.maxLH/red.maxLH_e) = (" << ndf_run << "/" << maxLHPerRun[i]/ndf_run << "/" << expectedMaxLHPerRun[i]/ndf_run << ")"; } } @@ -1938,60 +1975,69 @@ Bool_t PFitter::ExecuteSave(Bool_t firstSave) // calculate expected chisq std::vector param; Double_t totalExpectedChisq = 0.0; - std::vector expectedChisqPerHisto; + std::vector expectedchisqPerRun; std::vector ndfPerHisto; for (UInt_t i=0; iCalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto); + fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedchisqPerRun); // calculate chisq per run - std::vector chisqPerHisto; + std::vector chisqPerRun; for (UInt_t i=0; iGetMsrRunList()->size(); i++) { if (fUseChi2) - chisqPerHisto.push_back(fRunListCollection->GetSingleRunChisq(param, i)); + chisqPerRun.push_back(fRunListCollection->GetSingleRunChisq(param, i)); else - chisqPerHisto.push_back(fRunListCollection->GetSingleRunMaximumLikelihood(param, i)); + chisqPerRun.push_back(fRunListCollection->GetSingleRunMaximumLikelihood(param, i)); } if (totalExpectedChisq != 0.0) { // i.e. applicable for single histogram fits only // get the ndf's of the histos - UInt_t ndf_histo; - for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); - ndfPerHisto.push_back(ndf_histo); + UInt_t ndf_run; + for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + ndfPerHisto.push_back(ndf_run); } // feed the msr-file handler PMsrStatisticStructure *statistics = fRunInfo->GetMsrStatistic(); if (statistics) { - statistics->fMinPerHisto = chisqPerHisto; + statistics->fMinPerHisto = chisqPerRun; statistics->fMinExpected = totalExpectedChisq; - statistics->fMinExpectedPerHisto = expectedChisqPerHisto; + statistics->fMinExpectedPerHisto = expectedchisqPerRun; statistics->fNdfPerHisto = ndfPerHisto; } - } else if (chisqPerHisto.size() > 1) { // in case expected chisq is not applicable like for asymmetry fits - UInt_t ndf_histo = 0; - for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); - ndfPerHisto.push_back(ndf_histo); + } else if (chisqPerRun.size() > 1) { // in case expected chisq is not applicable like for asymmetry fits + UInt_t ndf_run = 0; + for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + ndfPerHisto.push_back(ndf_run); } // feed the msr-file handler PMsrStatisticStructure *statistics = fRunInfo->GetMsrStatistic(); if (statistics) { - statistics->fMinPerHisto = chisqPerHisto; + statistics->fMinPerHisto = chisqPerRun; statistics->fNdfPerHisto = ndfPerHisto; } } + // check if sector command has been requested + if (fSectorFlag) { + PDoubleVector error; + for (UInt_t i=0; i> PFitter::ExecuteSave(): will write minuit2 output file ..." << std::endl; @@ -2026,7 +2072,7 @@ Bool_t PFitter::ExecuteSave(Bool_t firstSave) fout << std::endl; fElapsedTime.clear(); - // write global informations + // write global information fout << std::endl << " Fval() = " << mnState.Fval() << ", Edm() = " << mnState.Edm() << ", NFcn() = " << mnState.NFcn(); fout << std::endl; fout << std::endl << "*************************************************************************"; @@ -2246,6 +2292,41 @@ Bool_t PFitter::ExecuteSave(Bool_t firstSave) } else { fout << std::endl << " no correlation coefficients available"; } + + fout << std::endl; + fout << std::endl << "*************************************************************************"; + fout << std::endl << " chisq/maxLH RESULT "; + fout << std::endl << "*************************************************************************"; + PMsrStatisticStructure *statistics = fRunInfo->GetMsrStatistic(); + + // get time range and write it + Double_t fitStartTime = PMUSR_UNDEFINED, fitEndTime = PMUSR_UNDEFINED; + // first check if there is a global block with a valid time range + PMsrGlobalBlock *global = fRunInfo->GetMsrGlobal(); + fitStartTime = global->GetFitRange(0); + fitEndTime = global->GetFitRange(1); + if (fitStartTime == PMUSR_UNDEFINED) { // no global time range, hence take the one from the first run block + PMsrRunList *run = fRunInfo->GetMsrRunList(); + fitStartTime = run->at(0).GetFitRange(0); + fitEndTime = run->at(0).GetFitRange(1); + } + fout.setf(std::ios::floatfield); + fout << std::endl << " Time Range: " << fitStartTime << ", " << fitEndTime << std::endl; + if (fUseChi2) { + fout.setf(std::ios::fixed, std::ios::floatfield); + fout << std::endl << " chisq = " << std::setprecision(4) << statistics->fMin << ", NDF = " << statistics->fNdf << ", chisq/NDF = " << std::setprecision(6) << statistics->fMin/statistics->fNdf; + if (statistics->fMinExpected > 0) + fout << std::endl << " chisq_e = " << std::setprecision(4) << statistics->fMinExpected << ", NDF = " << statistics->fNdf << ", chisq_e/NDF = " << std::setprecision(6) << statistics->fMinExpected/statistics->fNdf; + } else { // maxLH + fout.setf(std::ios::fixed, std::ios::floatfield); + fout << std::endl << " maxLH = " << std::setprecision(4) << statistics->fMin << ", NDF = " << statistics->fNdf << ", maxLH/NDF = " << std::setprecision(6) << statistics->fMin/statistics->fNdf; + if (statistics->fMinExpected > 0) + fout << std::endl << " maxLH_e = " << std::setprecision(4) << statistics->fMinExpected << ", NDF = " << statistics->fNdf << ", maxLH_e/NDF = " << std::setprecision(6) << statistics->fMinExpected/statistics->fNdf; + } + + if (fSectorFlag) + ExecuteSector(fout); + fout << std::endl; fout << std::endl << "*************************************************************************"; fout << std::endl << " DONE "; @@ -2338,23 +2419,112 @@ Bool_t PFitter::ExecuteSimplex() //-------------------------------------------------------------------------- /** *

Collect all the necessary chisq/maxLH sector information. + * + * @param param parameter value vector of the converged fit. + * @param error step value vector of the converged fit. */ -void PFitter::PrepareSector() +void PFitter::PrepareSector(PDoubleVector ¶m, PDoubleVector &error) { - // NOT YET IMPLEMENTED //as35 + Double_t val; + UInt_t ndf; -/* - for (UInt_t i=0; i +++++" << std::endl; - std::cout << "debug> sector " << i << std::endl; - std::cout << "debug> noOfRuns: " << fSector[i].GetNoRuns() << std::endl; - for (UInt_t j=0; j ----" << std::endl; - std::cout << "debug> RUN BLOCK " << j << std::endl; - std::cout << "debug> Sector Time Range: " << fSector[i].GetTimeRangeFirst(j) << ", " << fSector[i].GetTimeRangeLast() << std::endl; + Int_t usedParams = 0; + for (UInt_t i=0; iSetFitRange(secFitRange); + // calculate chisq + val = (*fFitterFcn)(param); + fSector[k].SetChisq(val); + // calculate NDF + ndf = static_cast(fFitterFcn->GetTotalNoOfFittedBins()) - usedParams; + fSector[k].SetNDF(ndf); + // calculate expected chisq + totalExpectedChisq = 0.0; + fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerRun); + fSector[k].SetExpectedChisq(totalExpectedChisq); + // calculate chisq per run + for (UInt_t i=0; iGetMsrRunList()->size(); i++) { + chisqPerRun.push_back(fRunListCollection->GetSingleRunChisq(param, i)); + fSector[k].SetChisq(chisqPerRun[i], i); + fSector[k].SetExpectedChisq(expectedChisqPerRun[i], i); + } + + if (totalExpectedChisq != 0.0) { + UInt_t ndf_run = 0; + for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_run > 0) { + fSector[k].SetNDF(ndf_run, i); + } + } + } else if (chisqPerRun.size() > 0) { // in case expected chisq is not applicable like for asymmetry fits + UInt_t ndf_run = 0; + for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_run > 0) { + fSector[k].SetNDF(ndf_run, i); + } + } + } + // clean up + chisqPerRun.clear(); + expectedChisqPerRun.clear(); + } + } else { + Double_t totalExpectedMaxLH = 0.0; + PDoubleVector expectedMaxLHPerRun; + PDoubleVector maxLHPerRun; + for (UInt_t k=0; kSetFitRange(secFitRange); + // calculate maxLH + val = (*fFitterFcn)(param); + fSector[k].SetChisq(val); + // calculate NDF + ndf = static_cast(fFitterFcn->GetTotalNoOfFittedBins()) - usedParams; + fSector[k].SetNDF(ndf); + // calculate expected maxLH + totalExpectedMaxLH = 0.0; + fFitterFcn->CalcExpectedChiSquare(param, totalExpectedMaxLH, expectedMaxLHPerRun); + fSector[k].SetExpectedChisq(totalExpectedMaxLH); + // calculate maxLH per run + for (UInt_t i=0; iGetMsrRunList()->size(); i++) { + maxLHPerRun.push_back(fRunListCollection->GetSingleRunMaximumLikelihood(param, i)); + fSector[k].SetChisq(maxLHPerRun[i], i); + fSector[k].SetExpectedChisq(expectedMaxLHPerRun[i], i); + } + + if (totalExpectedMaxLH != 0.0) { + UInt_t ndf_run = 0; + for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_run > 0) { + fSector[k].SetNDF(ndf_run, i); + } + } + } + + // clean up + maxLHPerRun.clear(); + expectedMaxLHPerRun.clear(); } } -*/ } //-------------------------------------------------------------------------- @@ -2366,9 +2536,50 @@ void PFitter::PrepareSector() * * return: if the sector command was successful, otherwise return flase. */ -Bool_t PFitter::ExecuteSector() +Bool_t PFitter::ExecuteSector(std::ofstream &fout) { - // NOT YET IMPLEMENTED //as35 + fout << std::endl; + fout << std::endl << "*************************************************************************"; + fout << std::endl << " SECTOR RESULTS"; + fout << std::endl << "*************************************************************************"; + + if (fUseChi2) { + for (UInt_t i=0; i 0) { // in case expected chisq is not applicable like for asymmetry fits + for (UInt_t j=0; j 0) { // in case expected chisq is not applicable like for asymmetry fits + for (UInt_t j=0; j &par, Double_ if (fUseChi2) { // single histo for (UInt_t i=0; iGetNoOfSingleHisto(); i++) { - value = fRunListCollection->GetSingleHistoChisqExpected(par, i); // calculate the expected chisq for single histo run block 'i' + value = fRunListCollection->GetSingleRunChisqExpected(par, i); // calculate the expected chisq for single histo run block 'i' expectedChisqPerRun.push_back(value); totalExpectedChisq += value; } } else { // log max. likelihood // single histo for (UInt_t i=0; iGetNoOfSingleHisto(); i++) { - value = fRunListCollection->GetSingleHistoMaximumLikelihoodExpected(par, i); // calculate the expected mlh for single histo run block 'i' + value = fRunListCollection->GetSingleRunMaximumLikelihoodExpected(par, i); // calculate the expected mlh for single histo run block 'i' expectedChisqPerRun.push_back(value); totalExpectedChisq += value; } diff --git a/src/classes/PRunListCollection.cpp b/src/classes/PRunListCollection.cpp index e456ccd2..c201976d 100644 --- a/src/classes/PRunListCollection.cpp +++ b/src/classes/PRunListCollection.cpp @@ -382,38 +382,39 @@ Double_t PRunListCollection::GetNonMusrChisq(const std::vector& par) c } //-------------------------------------------------------------------------- -// GetSingleHistoChisqExpected (public) +// GetSingleRunChisqExpected (public) //-------------------------------------------------------------------------- /** - *

Calculates expected chi-square of the single histogram with run block index idx of a msr-file. + *

Calculates expected chi-square of the run block index idx of a msr-file. * * return: - * - expected chi-square of for a single histogram + * - expected chi-square of for a single run block * * \param par fit parameter vector * \param idx run block index */ -Double_t PRunListCollection::GetSingleHistoChisqExpected(const std::vector& par, const UInt_t idx) const +Double_t PRunListCollection::GetSingleRunChisqExpected(const std::vector& par, const UInt_t idx) const { Double_t expectedChisq = 0.0; if (idx > fMsrInfo->GetMsrRunList()->size()) { - std::cerr << ">> PRunListCollection::GetSingleHistoChisqExpected() **ERROR** idx=" << idx << " is out of range [0.." << fMsrInfo->GetMsrRunList()->size() << "[" << std::endl << std::endl; + std::cerr << ">> PRunListCollection::GetSingleRunChisqExpected() **ERROR** idx=" << idx << " is out of range [0.." << fMsrInfo->GetMsrRunList()->size() << "[" << std::endl << std::endl; return expectedChisq; } + UInt_t subIdx = 0; Int_t type = fMsrInfo->GetMsrRunList()->at(idx).GetFitType(); if (type == -1) { // i.e. not found in the RUN block, try the GLOBAL block type = fMsrInfo->GetMsrGlobal()->GetFitType(); + subIdx = idx; + } else { // found in the RUN block + // count how many entries of this fit-type are present up to idx + for (UInt_t i=0; iGetMsrRunList()->at(i).GetFitType() == type) + subIdx++; + } } - // count how many entries of this fit-type are present up to idx - UInt_t subIdx = 0; - for (UInt_t i=0; iGetMsrRunList()->at(i).GetFitType() == type) || - (fMsrInfo->GetMsrRunList()->at(i).GetFitType() == -1)) // the -1 is needed if there is a global section - subIdx++; - } // return the chisq of the single run switch (type) { @@ -661,36 +662,37 @@ Double_t PRunListCollection::GetNonMusrMaximumLikelihood(const std::vectorCalculates expected mlh of the single histogram with run block index idx of a msr-file. + *

Calculates expected mlh of the run block index idx of a msr-file. * * return: - * - expected mlh of for a single histogram + * - expected mlh of for a single run block * * \param par fit parameter vector * \param idx run block index */ -Double_t PRunListCollection::GetSingleHistoMaximumLikelihoodExpected(const std::vector& par, const UInt_t idx) const +Double_t PRunListCollection::GetSingleRunMaximumLikelihoodExpected(const std::vector& par, const UInt_t idx) const { Double_t expected_mlh = 0.0; if (idx > fMsrInfo->GetMsrRunList()->size()) { - std::cerr << ">> PRunListCollection::GetSingleHistoMaximumLikelihoodExpected() **ERROR** idx=" << idx << " is out of range [0.." << fMsrInfo->GetMsrRunList()->size() << "[" << std::endl << std::endl; + std::cerr << ">> PRunListCollection::GetSingleRunMaximumLikelihoodExpected() **ERROR** idx=" << idx << " is out of range [0.." << fMsrInfo->GetMsrRunList()->size() << "[" << std::endl << std::endl; return expected_mlh; } + UInt_t subIdx = 0; Int_t type = fMsrInfo->GetMsrRunList()->at(idx).GetFitType(); if (type == -1) { // i.e. not found in the RUN block, try the GLOBAL block type = fMsrInfo->GetMsrGlobal()->GetFitType(); - } - - // count how many entries of this fit-type are present up to idx - UInt_t subIdx = 0; - for (UInt_t i=0; iGetMsrRunList()->at(i).GetFitType() == type) - subIdx++; + subIdx = idx; + } else { // found in the RUN block + // count how many entries of this fit-type are present up to idx + for (UInt_t i=0; iGetMsrRunList()->at(i).GetFitType() == type) + subIdx++; + } } // return the mlh of the single run @@ -772,16 +774,17 @@ UInt_t PRunListCollection::GetNoOfBinsFitted(const UInt_t idx) const return result; } - Int_t type = fMsrInfo->GetMsrRunList()->at(idx).GetFitType(); - if (type == -1) { // i.e. not forun in the RUN block, try the GLOBAL block - type = fMsrInfo->GetMsrGlobal()->GetFitType(); - } - - // count how many entries of this fit-type are present up to idx UInt_t subIdx = 0; - for (UInt_t i=0; iGetMsrRunList()->at(i).GetFitType() == type) - subIdx++; + Int_t type = fMsrInfo->GetMsrRunList()->at(idx).GetFitType(); + if (type == -1) { // i.e. not found in the RUN block, try the GLOBAL block + type = fMsrInfo->GetMsrGlobal()->GetFitType(); + subIdx = idx; + } else { // found in the RUN block + // count how many entries of this fit-type are present up to idx + for (UInt_t i=0; iGetMsrRunList()->at(i).GetFitType() == type) + subIdx++; + } } // return the chisq of the single run diff --git a/src/classes/PRunSingleHisto.cpp b/src/classes/PRunSingleHisto.cpp index b3d15f09..5d114688 100644 --- a/src/classes/PRunSingleHisto.cpp +++ b/src/classes/PRunSingleHisto.cpp @@ -276,9 +276,9 @@ Double_t PRunSingleHisto::CalcChiSquareExpected(const std::vector& par 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) + #pragma omp parallel for default(shared) private(i,time,theo,diff) schedule(dynamic,chunk) reduction(+:chisq) #endif - for (i=fStartTimeBin; i < fEndTimeBin; ++i) { + for (i=fStartTimeBin; i(i)*fData.GetDataTimeStep(); theo = N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg; diff = fData.GetValue()->at(i) - theo; @@ -373,9 +373,8 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector& par) time = fData.GetDataTimeStart() + static_cast(i)*fData.GetDataTimeStep(); // calculate theory for the given parameter set theo = N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg; - theo *= normalizer; - data = normalizer*fData.GetValue()->at(i); + data = fData.GetValue()->at(i); if (theo <= 0.0) { std::cerr << ">> PRunSingleHisto::CalcMaxLikelihood: **WARNING** NEGATIVE theory!!" << std::endl; @@ -389,7 +388,7 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector& par) } } - return 2.0*mllh; + return normalizer*2.0*mllh; } //-------------------------------------------------------------------------- @@ -472,9 +471,8 @@ Double_t PRunSingleHisto::CalcMaxLikelihoodExpected(const std::vector& time = fData.GetDataTimeStart() + static_cast(i)*fData.GetDataTimeStep(); // calculate theory for the given parameter set theo = N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg; - theo *= normalizer; - data = normalizer*fData.GetValue()->at(i); + data = fData.GetValue()->at(i); if (theo <= 0.0) { std::cerr << ">> PRunSingleHisto::CalcMaxLikelihood: **WARNING** NEGATIVE theory!!" << std::endl; @@ -486,7 +484,7 @@ Double_t PRunSingleHisto::CalcMaxLikelihoodExpected(const std::vector& } } - return 2.0*mllh; + return normalizer*2.0*mllh; } //-------------------------------------------------------------------------- diff --git a/src/include/PFitter.h b/src/include/PFitter.h index 56cd6176..b0a55959 100644 --- a/src/include/PFitter.h +++ b/src/include/PFitter.h @@ -75,25 +75,31 @@ class PSectorChisq void SetSectorTime(Double_t last) { fLast = last; } void SetChisq(Double_t chisq) { fChisq = chisq; } void SetChisq(Double_t chisq, UInt_t idx); - void SetNDF(Double_t ndf) { fNDF = ndf; } - void SetNDF(Double_t ndf, UInt_t idx); + void SetExpectedChisq(Double_t expChisq) { fExpectedChisq = expChisq; } + void SetExpectedChisq(Double_t chisq, UInt_t idx); + void SetNDF(UInt_t ndf) { fNDF = ndf; } + void SetNDF(UInt_t ndf, UInt_t idx); Double_t GetTimeRangeFirst(UInt_t idx); Double_t GetTimeRangeLast() { return fLast; } Double_t GetChisq() { return fChisq; } - UInt_t GetNDF() { return fNDF; } - UInt_t GetNoRuns() { return fNoOfRuns; } Double_t GetChisq(UInt_t idx); + Double_t GetExpectedChisq() { return fExpectedChisq; } + Double_t GetExpectedChisq(UInt_t idx); + UInt_t GetNDF() { return fNDF; } UInt_t GetNDF(UInt_t idx); + UInt_t GetNoRuns() { return fNoOfRuns; } private: UInt_t fNoOfRuns; ///< number of runs presesent Double_t fLast; ///< requested time stamp Double_t fChisq; ///< chisq or maxLH for the sector + Double_t fExpectedChisq; ///< keep the expected chisq or maxLH for the sector UInt_t fNDF; ///< NDF for the sector - std::vector fFirst; ///< time stamp for fgb for a given run - std::vector fChisqRun; ///< chisq or maxLH for the sector and run - std::vector fNDFRun; ///< NDF for the sector and run + PDoubleVector fFirst; ///< time stamp for fgb for a given run + PDoubleVector fChisqRun; ///< chisq or maxLH for the sector and run + PDoubleVector fExpectedChisqRun; ///< expected chisq or maxLH for the sector and run + PUIntVector fNDFRun; ///< NDF for the sector and run }; //----------------------------------------------------------------------------- @@ -146,7 +152,7 @@ class PFitter PStringVector fElapsedTime; - Bool_t fSectorFlag; ///< sector command present flag + Bool_t fSectorFlag; ///< sector command present flag std::vector fSector; ///< stores all chisq/maxLH sector information // commands @@ -167,8 +173,8 @@ class PFitter Bool_t ExecuteScan(); Bool_t ExecuteSave(Bool_t first); Bool_t ExecuteSimplex(); - void PrepareSector(); - Bool_t ExecuteSector(); + void PrepareSector(PDoubleVector ¶m, PDoubleVector &error); + Bool_t ExecuteSector(std::ofstream &fout); Double_t MilliTime(); }; diff --git a/src/include/PRunListCollection.h b/src/include/PRunListCollection.h index 7f41c1c1..19aef805 100644 --- a/src/include/PRunListCollection.h +++ b/src/include/PRunListCollection.h @@ -67,7 +67,7 @@ class PRunListCollection virtual Double_t GetMuMinusChisq(const std::vector& par) const; virtual Double_t GetNonMusrChisq(const std::vector& par) const; - virtual Double_t GetSingleHistoChisqExpected(const std::vector& par, const UInt_t idx) const; + virtual Double_t GetSingleRunChisqExpected(const std::vector& par, const UInt_t idx) const; virtual Double_t GetSingleRunChisq(const std::vector& par, const UInt_t idx) const; virtual Double_t GetSingleHistoMaximumLikelihood(const std::vector& par) const; @@ -78,7 +78,7 @@ class PRunListCollection virtual Double_t GetMuMinusMaximumLikelihood(const std::vector& par) const; virtual Double_t GetNonMusrMaximumLikelihood(const std::vector& par) const; - virtual Double_t GetSingleHistoMaximumLikelihoodExpected(const std::vector& par, const UInt_t idx) const; + virtual Double_t GetSingleRunMaximumLikelihoodExpected(const std::vector& par, const UInt_t idx) const; virtual Double_t GetSingleRunMaximumLikelihood(const std::vector& par, const UInt_t idx) const; virtual UInt_t GetNoOfBinsFitted(const UInt_t idx) const; From e9247ed3d5ea8de683f13e501230a41c47529ed6 Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Mon, 3 Feb 2020 21:03:45 +0100 Subject: [PATCH 4/5] increment version. --- CMakeLists.txt | 2 +- ChangeLog | 6 ++++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index d0555e3f..01c1e94e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -5,7 +5,7 @@ if (CMAKE_VERSION GREATER_EQUAL 3.12) cmake_policy(SET CMP0075 NEW) endif (CMAKE_VERSION GREATER_EQUAL 3.12) -project(musrfit VERSION 1.5.1 LANGUAGES C CXX) +project(musrfit VERSION 1.5.2 LANGUAGES C CXX) #--- musrfit specific options ------------------------------------------------- option(nexus "build optional NeXus support. Needed for ISIS" OFF) diff --git a/ChangeLog b/ChangeLog index 6d003e47..7a2092b7 100644 --- a/ChangeLog +++ b/ChangeLog @@ -12,6 +12,12 @@ or https://bitbucket.org/muonspin/musrfit/commits/all +Release of V1.5.2, 2020/02/03 +============================= + +Implemented a SECTOR command. It allows to get chisq/maxLH information from +different time sector slots. + Release of V1.5.0, 2019/05/15 ============================= From c3e3af03113ce07a9c0f7824369df1c1d9ba138a Mon Sep 17 00:00:00 2001 From: Andreas Suter Date: Tue, 4 Feb 2020 19:37:31 +0100 Subject: [PATCH 5/5] update of the docu including the SECTOR cmd. --- doc/html/.buildinfo | 2 +- .../Makefile.TMyLibrary | 60 + doc/html/_sources/acknowledgement.rst.txt | 22 + doc/html/_sources/any2many.rst.txt | 8 + doc/html/_sources/bugtracking.rst.txt | 9 + doc/html/_sources/cite.rst.txt | 18 + doc/html/_sources/index.rst.txt | 31 + doc/html/_sources/msr2data.rst.txt | 374 + doc/html/_sources/mupp.rst.txt | 166 + doc/html/_sources/musr-root.rst.txt | 681 + doc/html/_sources/musredit.rst.txt | 515 + doc/html/_sources/setup-dks.rst.txt | 275 + doc/html/_sources/setup-standard.rst.txt | 1392 +++ doc/html/_sources/tutorial.rst.txt | 417 + doc/html/_sources/user-libs.rst.txt | 747 ++ doc/html/_sources/user-manual.rst.txt | 2394 ++++ doc/html/_static/basic.css | 282 +- doc/html/_static/doctools.js | 108 +- doc/html/_static/documentation_options.js | 10 + doc/html/_static/file.png | Bin 392 -> 286 bytes doc/html/_static/jquery-3.2.1.js | 10253 ++++++++++++++++ doc/html/_static/jquery.js | 6 +- doc/html/_static/language_data.js | 297 + doc/html/_static/minus.png | Bin 199 -> 90 bytes doc/html/_static/nature.css | 17 +- doc/html/_static/plus.png | Bin 199 -> 90 bytes doc/html/_static/searchtools.js | 360 +- doc/html/_static/underscore-1.3.1.js | 999 ++ doc/html/acknowledgement.html | 96 +- doc/html/any2many.html | 80 +- doc/html/bugtracking.html | 74 +- doc/html/cite.html | 90 +- doc/html/genindex.html | 1283 +- doc/html/index.html | 102 +- doc/html/msr2data.html | 444 +- doc/html/mupp.html | 236 +- doc/html/musr-root.html | 881 +- doc/html/musredit.html | 489 +- doc/html/objects.inv | Bin 1957 -> 2154 bytes doc/html/search.html | 54 +- doc/html/searchindex.js | 2 +- doc/html/setup-dks.html | 230 +- doc/html/setup-standard.html | 934 +- doc/html/tutorial.html | 384 +- doc/html/user-libs.html | 533 +- doc/html/user-manual.html | 2746 +++-- 46 files changed, 23269 insertions(+), 4832 deletions(-) create mode 100644 doc/html/_downloads/105384a3f798b4479ee2ecadd2e3112e/Makefile.TMyLibrary create mode 100644 doc/html/_sources/acknowledgement.rst.txt create mode 100644 doc/html/_sources/any2many.rst.txt create mode 100644 doc/html/_sources/bugtracking.rst.txt create mode 100644 doc/html/_sources/cite.rst.txt create mode 100644 doc/html/_sources/index.rst.txt create mode 100644 doc/html/_sources/msr2data.rst.txt create mode 100644 doc/html/_sources/mupp.rst.txt create mode 100644 doc/html/_sources/musr-root.rst.txt create mode 100644 doc/html/_sources/musredit.rst.txt create mode 100644 doc/html/_sources/setup-dks.rst.txt create mode 100644 doc/html/_sources/setup-standard.rst.txt create mode 100644 doc/html/_sources/tutorial.rst.txt create mode 100644 doc/html/_sources/user-libs.rst.txt create mode 100644 doc/html/_sources/user-manual.rst.txt create mode 100644 doc/html/_static/documentation_options.js create mode 100644 doc/html/_static/jquery-3.2.1.js create mode 100644 doc/html/_static/language_data.js create mode 100644 doc/html/_static/underscore-1.3.1.js diff --git a/doc/html/.buildinfo b/doc/html/.buildinfo index 7c5e4ed5..6bd6cce3 100644 --- a/doc/html/.buildinfo +++ b/doc/html/.buildinfo @@ -1,4 +1,4 @@ # Sphinx build info version 1 # This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done. -config: b0f54e143f8f1d163af95476812112f7 +config: 08712e994533f427d7efeeb2574ca61a tags: 645f666f9bcd5a90fca523b33c5a78b7 diff --git a/doc/html/_downloads/105384a3f798b4479ee2ecadd2e3112e/Makefile.TMyLibrary b/doc/html/_downloads/105384a3f798b4479ee2ecadd2e3112e/Makefile.TMyLibrary new file mode 100644 index 00000000..73a85c8e --- /dev/null +++ b/doc/html/_downloads/105384a3f798b4479ee2ecadd2e3112e/Makefile.TMyLibrary @@ -0,0 +1,60 @@ +#--------------------------------------------------- +# get compilation flags from root-config + +ROOTCFLAGS = $(shell $(ROOTSYS)/bin/root-config --cflags) + +#--------------------------------------------------- + +OS = LINUX +CXX = g++ +CXXFLAGS = -O3 -Wall -Wno-trigraphs -fPIC +LOCALINCLUDE = . +ROOTINCLUDE = $(ROOTSYS)/include +INCLUDES = -I$(LOCALINCLUDE) -I$(ROOTINCLUDE) +LD = g++ +LDFLAGS = +SOFLAGS = -O -shared + +# the output from the root-config script: +CXXFLAGS += $(ROOTCFLAGS) +LDFLAGS += + +# some definitions: headers (used to generate *Dict* stuff), sources, objects,... +OBJS = +OBJS += TMyFunction.o TMyLibraryDict.o + +SHLIB = libTMyLibrary.so + +# make the shared lib: +# +all: $(SHLIB) + +$(SHLIB): $(OBJS) + @echo "---> Building shared library $(SHLIB) ..." + /bin/rm -f $(SHLIB) + $(LD) $(OBJS) $(SOFLAGS) -o $(SHLIB) + @echo "done" + +# clean up: remove all object file (and core files) +# semicolon needed to tell make there is no source +# for this target! +# +clean:; @rm -f $(OBJS) *Dict* core* + @echo "---> removing $(OBJS)" + +# +$(OBJS): %.o: %.cpp + $(CXX) $(INCLUDES) $(CXXFLAGS) -c $< + +# Generate the ROOT CINT dictionary + +TMyLibraryDict.cpp: TMyFunction.h TMyLibraryLinkDef.h + @echo "Generating dictionary $@..." + rootcint -f $@ -c -p -I$(ROOTINCLUDE) $^ + +install: all + @echo "Installing shared lib: libTApproximation.so" +ifeq ($(OS),LINUX) + cp -pv $(SHLIB) $(ROOTSYS)/lib + cp -pv $(LOCALINCLUDE)/*.h $(ROOTSYS)/include +endif diff --git a/doc/html/_sources/acknowledgement.rst.txt b/doc/html/_sources/acknowledgement.rst.txt new file mode 100644 index 00000000..7e8c47fa --- /dev/null +++ b/doc/html/_sources/acknowledgement.rst.txt @@ -0,0 +1,22 @@ +.. include:: +.. index:: acknowledgment + +.. _acknowledgment: + +Acknowledgements +================ + + +**Bastian M. Wojek** + I am very much indebted to BMW for his rigorous testing of ``musrfit``, his many useful suggestions, contributions, and for the + largest part of the user manual of ``musrfit`` which makes it accessible to a broader audience! Many thanks Bastian! + +**Uldis Locans** + I am very much indebted to Uldis work on :ref:`DKS ` enabling the GPU support for ``musrfit``. His kind, calm, and + extremely competent way to deal with his projects as well as to deal with the chaos of physicists way to think is admirable. Many thanks Uldis! + +**Zaher Salman** + Thanks for his beta-NMR and web-interface contributions to ``musrfit``! + +**Robert Scheuermann** + Thanks for his constant contructive input on ``musrfit``! diff --git a/doc/html/_sources/any2many.rst.txt b/doc/html/_sources/any2many.rst.txt new file mode 100644 index 00000000..6b5fd201 --- /dev/null +++ b/doc/html/_sources/any2many.rst.txt @@ -0,0 +1,8 @@ +.. include:: +.. index:: any2many + +any2many - a Universal |mgr|\SR-file-format converter +===================================================== + +``any2many`` allows to convert most |mgr|\SR-file-formats from one to the other. +For a detailed description see :ref:`here `. \ No newline at end of file diff --git a/doc/html/_sources/bugtracking.rst.txt b/doc/html/_sources/bugtracking.rst.txt new file mode 100644 index 00000000..55a5ff54 --- /dev/null +++ b/doc/html/_sources/bugtracking.rst.txt @@ -0,0 +1,9 @@ +.. index:: bugtracking +.. _bugtracking: + +Bugtracking +=========== + +For reporting bugs or requesting new features and improvements please use +the `bitbucket-repo `_ (preferred) +or send an e-mail to A. Suter at PSI. diff --git a/doc/html/_sources/cite.rst.txt b/doc/html/_sources/cite.rst.txt new file mode 100644 index 00000000..5a71e6f1 --- /dev/null +++ b/doc/html/_sources/cite.rst.txt @@ -0,0 +1,18 @@ +.. include:: +.. index:: cite +.. _cite: + +How to Cite ``musrfit``? +======================== + +Since quite some effort is going into the development and maintenance of the ``musrfit`` package, you should at least acknowledge it in your publication if you have used it to analyze your data. Even better of course is to cite it properly by the reference given beneath + + * A.\ Suter, B.M. Wojek, "Musrfit: A Free Platform-Independent Framework for |mgr|\SR Data Analysis", Physics Procedia **30**, 69 (2012). ``_ + +The GPU high speed ``musrfit`` version is utilizing ``DKS``. In case you are using this version, please also add the following citations + + * A.\ Adelmann, U. Locans, A. Suter, "The Dynamic Kernel Scheduler—Part 1", Computer Physics Communications **207**, 83 (2016). ``_ + * U.\ Locans, *et al.*, "Real-time computation of parameter fitting and image reconstruction using graphical processing units", Computer Physics Communications **215**, 71 (2017). ``_ + * U.\ Locans and A.\ Suter, "Musrfit – Real Time Parameter Fitting Using GPUs", JPS Conf. Proc. *21*, 011051 (2018). ``_ + + diff --git a/doc/html/_sources/index.rst.txt b/doc/html/_sources/index.rst.txt new file mode 100644 index 00000000..b4ed7d2f --- /dev/null +++ b/doc/html/_sources/index.rst.txt @@ -0,0 +1,31 @@ +.. musrfit docu documentation master file, created by + sphinx-quickstart on Sun Jun 17 11:00:32 2018. + You can adapt this file completely to your liking, but it should at least + contain the root `toctree` directive. + +Welcome to the musrfit documentation! +===================================== + +.. toctree:: + :maxdepth: 2 + + cite + tutorial + user-manual + user-libs + setup-standard + setup-dks + musredit + mupp + msr2data + any2many + musr-root + acknowledgement + bugtracking + + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`search` diff --git a/doc/html/_sources/msr2data.rst.txt b/doc/html/_sources/msr2data.rst.txt new file mode 100644 index 00000000..8e77bbca --- /dev/null +++ b/doc/html/_sources/msr2data.rst.txt @@ -0,0 +1,374 @@ +.. include:: +.. index:: msr2data + +.. _msr2data: + +msr2data - A Program for Automatically Processing Multiple ``musrfit`` msr Files +================================================================================ + +``msr2data`` (originally written by B. M. Wojek) is a program implemented in ``C++``. Its purpose is +to process multiple msr files (input files for ``musrfit``) with the same parameters and summarize the fitting +results either in a *TRIUMF DB* [#f1]_ or a *column ASCII* file. This allows essentially to + +#. Collect the fit parameters. +#. Generate *new* input msr files based on old ones. + +.. [#f1] For an abridged description of this format see `here `_. The DB files + produced by ``msr2data`` can be viewed for instance with :ref:`mupp ` or |mgr|\View `see here `_, however, + they are not completely backward-compatible to the original ``db language`` since the parameter names can be longer than five or + six characters! In order to establish this backward compatibility (if needed) the user has to ensure the correct length of the + parameter names in the msr files! + +.. _msr2data-basic-usage: + +Basic Types of Usage +-------------------- + +Apart from numerous :ref:`optional parameters ` that might be set, in principle there are four different ways of calling ``msr2data``. +These differ in how the list of runs which should be processed is supplied: + +**msr2data [optional parameters]** + A single run number. +**msr2data [optional parameters]** + An interval of run numbers is specified through the first and the last run number. The condition ```` < ```` is not necessary. +**msr2data \[ \] [optional parameters]** + Where ```` is one or a combination of the following: + + #. ``, , , ... `` : run numbers, *e.g.* 123 124, + #. ``-`` : a range, *e.g.* 123-125 -> 123 124 125, + #. ``::`` : a sequence, *e.g.* 123:127:2 -> 123 125 127. ```` has to be a positive integer. + #. A ```` can also combine (1)-(3), *e.g.* 123 128-130 133, etc. + +**msr2data [optional parameters]** + An ASCII file containing a list of run numbers and optional external parameters is passed to ``msr2data``. For the structure of the ASCII file + see :ref:`below `. + +All four basic types of calling ``msr2data`` contain the *mandatory* file-name ```` passed right after the list of runs. The meaning of +this ```` should become clear after giving examples for all four cases: + +.. code-block:: bash + + $ msr2data 8472 _tf_h13 + +generates the DB file ``out.db`` (can be changed by using the -o option) from ``8472_tf_h13.msr``. + +.. code-block:: bash + + $ msr2data 8472 8474 _tf_h13 + +generates the DB file ``out.db`` (can be changed by using the -o option) from ``8472_tf_h13.msr``, ``8473_tf_h13.msr``, and ``8474_tf_h13.msr``. + +.. code-block:: bash + + $ msr2data [8472 8470] _tf_h13 + +generates the DB file ``out.db`` (can be changed by using the -o option) from ``8472_tf_h13.msr`` and ``8470_tf_h13.msr``. + +.. code-block:: bash + + $ msr2data [8470:8474:2] _tf_h13 + +generates the DB file ``out.db`` (can be changed by using the -o option) from ``8470_tf_h13.msr``, ``8472_tf_h13.msr``, and ``8474_tf_h13.msr``. + +.. _run-list-file_structure: + +Run List File Structure ++++++++++++++++++++++++ + +.. code-block:: bash + + $ msr2data run.list _tf_h13 + +generates the DB file ``out.db`` (can be changed by using the -o option) from all runs listed in the ASCII file ``run.list`` in the working directory. +In this file it is also possible to include *external* parameters which should be put in the resulting DB file. The structure of the ``run.list`` is the following: + +:: + + RUN VAR1 VAR2 VAR3 ... + 8460 200 27.1 46.2 ... + 8472 205 27.1 46.3 ... + 8453 210 27.2 45.9 ... + · · · · + · · · · + · · · · + +*The first not commented and not empty line determines the parameter names and labels and has to be present!* + +It is allowed to add comments (with a preceding '#') or empty lines to the run-list file. + +The following should be mentioned together with the above examples: + +* The output files in the examples above are only newly created if they did *not* exist before invoking ``msr2data``. + If the files were already present the msr file data would be appended! +* If the files have been newly created, also the DB file header is written. If the files were present before, only + the data blocks are appended. The output of the header can either be forced or completely suppressed with the ``header`` + and ``noheader`` options as shall be seen later. +* If the ``musrfit`` output files do not have an ```` as specified above like ``8472.msr`` one has to call ``msr2data`` like in the following example: + + .. code-block:: bash + + $ msr2data 8472 8460 "" + +.. _msr2data-opt-param: + +Optional Parameters +------------------- + +As mentioned already above there are some optional parameters which change the behavior of ``msr2data`` and can be passed in any order. Here is a complete list: + +**data** + The output file format is changed to a simple column ASCII file (default output file name: out.dat). +**new** + An existing output file is deleted before new information is written to it. +**header** + Force the output of the file header even if the output file was present before. +**noheader** + The output of the file header is suppressed—also if the output file is newly created. + If either both or none of the header options are given, ``msr2data`` writes the file header only to new files + and it solely appends the data blocks to an existing output file assuming that the header is present already. +**nosummary** + There will be no attempt to read additional information like the temperature or the applied magnetic field from + the data files even if these information were present there. +**paramList ** + option used to select the parameters which shall be exported. ```` is a list of parameter numbers to be exported. + Allowed lists are: ``-``, *e.g.* ``1-16`` will export parameters 1 to 16. Space separated numbers, *e.g.:* ``1 3 5``. + A combination of both is possible, *e.g.* ``1-16 19 31 62``, and so on. +**-o, -o ** + The processed data will be written to the file ```` instead of the default ``out.db`` or ``out.dat``. + If ```` is equal to none (case-insensitive) the parameter data are not appended to any output file. +**fit** + Additionally to the final data collection ``msr2data`` will invoke ``musrfit`` to fit the specified runs. + All msr files are assumed to be present, none is newly generated! +**fit-