LMU/musrfit
5
1
Fork 0
Code Issues 5 Pull Requests Actions Packages Projects Releases 4 Activity
Files
f2ab8dbeb9020c9ab7ed3ff8faedf601efd493d0
musrfit/src/classes/PFitter.cpp
T
suter_a 66aa847b58 PFitter: replace TObjArray/TObjString with PStringUtils 
Replace ROOT's TString::Tokenize() + TObjArray/TObjString token
handling with the dependency-free PStringUtils::Split() across all
command/theory parsing in PFitter (GetPhaseParams, GetParFromFun,
CheckCommands and the Execute* helpers). Split() mirrors Tokenize()
semantics (delimiter-set, skips empty tokens), so token counts and
indices are unchanged. Each token is still copied into a TString,
so the downstream Atoi/Atof/IsFloat/IsDigit/Contains/CompareTo logic
stays as-is.

Using a std::vector<std::string> removes the manual TObjArray
cleanup and incidentally fixes three pre-existing leaks: the
tokens array in ExecuteFitRange was never freed, and the early
return paths in ExecutePrintLevel and the SECTOR check in
CheckCommands skipped the cleanup.

Build of libPMusr and musrfit is clean; full ctest suite passes
(85/85).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
2026-06-06 15:38:15 +02:00

2904 lines
111 KiB
C++
Raw Blame History

/***************************************************************************
PFitter.cpp
Author: Andreas Suter
e-mail: andreas.suter@psi.ch
***************************************************************************/
/***************************************************************************
* Copyright (C) 2007-2026 by Andreas Suter *
* andreas.suter@psi.ch *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program; if not, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#ifdef HAVE_GOMP
#include <omp.h>
#endif
#include <iostream>
#include <iomanip>
#include <fstream>
#include <limits>
#include <cmath>
#include <boost/variant/variant.hpp>
#include <sys/time.h>
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnContours.h"
#include "Minuit2/MnHesse.h"
#include "Minuit2/MnMinimize.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/MnMinos.h"
#include "Minuit2/MnPlot.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/MnScan.h"
#include "Minuit2/MnSimplex.h"
#include "Minuit2/MnStrategy.h"
#include "Minuit2/MnUserParameterState.h"
#include "Minuit2/MinosError.h"
#include <TCanvas.h>
#include <TH2.h>
#include <TFile.h>
#include <TDatime.h>
#include <TString.h>
#include "PFitter.h"
#include "PStringUtils.h"
//+++ PSectorChisq class +++++++++++++++++++++++++++++++++++++++++++++++++++
//--------------------------------------------------------------------------
// Constructor
//--------------------------------------------------------------------------
/**
* <p>Constructor.
*/
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; i<fNoOfRuns; i++) {
fFirst[i] = 0.0;
fChisqRun[i] = 0.0;
fExpectedChisqRun[i] = 0.0;
fNDFRun[i] = 0;
}
}
//--------------------------------------------------------------------------
// SetRunFirstTime
//--------------------------------------------------------------------------
/**
* <p>Set the time of the fgb of a given RUN
*
* @param first time stamp of the fgb
* @param idx index of the RUN
*/
void PSectorChisq::SetRunFirstTime(Double_t first, UInt_t idx)
{
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;
}
//--------------------------------------------------------------------------
// SetChisq
//--------------------------------------------------------------------------
/**
* <p>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)
{
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;
}
//--------------------------------------------------------------------------
// SetExpectedChisq
//--------------------------------------------------------------------------
/**
* <p>Set 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
//--------------------------------------------------------------------------
/**
* <p>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(UInt_t ndf, UInt_t idx)
{
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
//--------------------------------------------------------------------------
/**
* <p>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
*
* <b>return:</b> 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];
}
//--------------------------------------------------------------------------
// GetChisq
//--------------------------------------------------------------------------
/**
* <p>Get chisq/maxLH of the run with index idx
*
* @param idx index of the run
*
* <b>return:</b> 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];
}
//--------------------------------------------------------------------------
// GetExpectedChisq
//--------------------------------------------------------------------------
/**
* <p>Get expected chisq/maxLH of the run with index idx
*
* @param idx index of the run
*
* <b>return:</b> 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
//--------------------------------------------------------------------------
/**
* <p>Get NDF of the run with index idx
*
* @param idx index of the run
*
* <b>return:</b> 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
//--------------------------------------------------------------------------
/**
* @brief Constructor for the fitting engine.
*
* Initializes the fitter with MSR configuration and preprocessed data.
* Sets up the fitting environment including:
* - Parameter lists and command queues
* - Fit mode (χ² vs. likelihood)
* - Strategy and print level defaults
* - Original fit ranges for RANGE commands
* - Phase parameter identification
*
* The constructor validates the COMMANDS block and creates the objective
* function (PFitterFcn) but does not start the fit. Call DoFit() to execute.
*
* @param runInfo Pointer to MSR file handler containing fit configuration
* @param runListCollection Pointer to preprocessed data collection
* @param chisq_only If true, only evaluate χ² without fitting (useful for validation)
* @param yaml_out If true, generate YAML output file with fit results
*
* @note Sets fIsValid = false if command validation fails. Check IsValid() before DoFit().
*
* @par Initialization sequence:
* 1. Copy parameters and commands from MSR handler
* 2. Store original fit ranges (for FIT_RANGE command)
* 3. Parse and validate COMMANDS block
* 4. Identify phase parameters (for angle wrapping)
* 5. Create objective function interface
*/
PFitter::PFitter(PMsrHandler *runInfo, PRunListCollection *runListCollection, Bool_t chisq_only, Bool_t yaml_out) :
fChisqOnly(chisq_only), fYamlOut(yaml_out), fRunInfo(runInfo), fRunListCollection(runListCollection)
{
// initialize variables
fIsScanOnly = true;
fConverged = false;
fUseChi2 = true; // chi^2 is the default
fStrategy = 1; // 0=low, 1=default, 2=high
fSectorFlag = false;
fParams = *(runInfo->GetMsrParamList());
fCmdLines = *runInfo->GetMsrCommands();
// init class variables
fScanAll = true;
fScanParameter[0] = 0;
fScanParameter[1] = 0;
fScanNoPoints = 41; // minuit2 default
fScanLow = 0.0; // minuit2 default, i.e. 2 std deviations
fScanHigh = 0.0; // minuit2 default, i.e. 2 std deviations
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; i<runs->size(); i++) {
range.first = (*runs)[i].GetFitRange(0);
range.second = (*runs)[i].GetFitRange(1);
if (range.first == PMUSR_UNDEFINED)
fOriginalFitRange.push_back(rangeGlob);
else
fOriginalFitRange.push_back(range);
}
// check msr minuit commands
if (!CheckCommands()) {
return;
}
// create phase bool array
GetPhaseParams();
// create fit function object
fFitterFcn = std::make_unique<PFitterFcn>(runListCollection, fUseChi2);
}
//--------------------------------------------------------------------------
// Destructor
//--------------------------------------------------------------------------
/**
* @brief Destructor - Cleans up dynamically allocated resources.
*
* Frees memory used by command lists, scan data, and timing information.
* Smart pointers (fFitterFcn, fFcnMin) are automatically cleaned up.
*/
PFitter::~PFitter()
{
fCmdList.clear();
fScanData.clear();
fElapsedTime.clear();
}
//--------------------------------------------------------------------------
// GetPhaseParams (private)
//--------------------------------------------------------------------------
/**
* @brief Identifies which parameters represent phase angles.
*
* Scans the THEORY block to detect parameters used as phases in standard
* muSR functions. Phase parameters are flagged to enable automatic wrapping
* to the interval [-360°, +360°] during fitting, preventing meaningless
* phase values outside this range.
*
* Recognized phase-containing functions:
* - TFieldCos/tf: phase is 1st parameter
* - bessel/b: phase is 1st parameter
* - skewedGss/skg: phase is 1st parameter
* - staticNKTF/snktf: phase is 1st parameter
* - dynamicNKTF/dnktf: phase is 1st parameter
* - internFld/if: phase is 2nd parameter
* - internBsl/ib: phase is 2nd parameter
* - muMinusExpTF/mmsetf: phase is 5th parameter
*
* Phase references can be:
* - Direct parameter numbers (e.g., "7" → par7)
* - Function references (e.g., "fun3" → all params in function 3)
* - Map references (e.g., "map2" → parameters mapped across runs)
*
* @note Populates fPhase vector where fPhase[i] = true means parameter i+1 is a phase.
* @note User-defined functions cannot be automatically analyzed for phases.
*/
void PFitter::GetPhaseParams()
{
fPhase.resize(fRunInfo->GetNoOfParams());
for (unsigned int i=0; i<fPhase.size(); i++)
fPhase[i] = false;
// analyze theory block for parameters. Phases are present in the following
// default functions:
// user functions cannot be checked!
PMsrLines *theo = fRunInfo->GetMsrTheory();
TString str;
int pos = -1;
for (unsigned int i=0; i<theo->size(); i++) {
pos = -1;
TString line = theo->at(i).fLine;
if (line.Contains("TFieldCos") || line.Contains("tf ") ||
line.Contains("bessel") || line.Contains("b ") ||
line.Contains("skewedGss") || line.Contains("skg ") ||
line.Contains("staticNKTF") || line.Contains("snktf ") ||
line.Contains("dynamicNKTF") || line.Contains("dnktf ")) { // phase is 1st param
pos = 1;
}
if (line.Contains("internFld") || line.Contains("if ") ||
line.Contains("internBsl") || line.Contains("ib ")) { // phase is 2nd param
pos = 2;
}
if (line.Contains("muMinusExpTF") || line.Contains("mmsetf ")) { // phase is 5th param
pos = 5;
}
if (pos == -1)
continue;
// extract phase token
std::vector<std::string> tok = PStringUtils::Split(line.Data(), " \t");
if (static_cast<int>(tok.size()) > pos) {
str = tok[pos].c_str();
}
// decode phase token. It can be funX, mapX, or a number
if (str.Contains("fun")) { // function
PIntVector parVec = GetParFromFun(str);
for (int i=0; i<parVec.size(); i++) {
if (parVec[i] <= fRunInfo->GetNoOfParams())
fPhase[parVec[i]-1] = true;
}
} else if (str.Contains("map")) { // map
PIntVector parVec = GetParFromMap(str);
for (int i=0; i<parVec.size(); i++) {
if (parVec[i] <= fRunInfo->GetNoOfParams())
fPhase[parVec[i]-1] = true;
}
} else { // must be a number
int idx = str.Atoi();
if (idx == 0) { // something went wrong, str is not an integer
std::cerr << "PFitter::GetPhaseParams(): **ERROR** str=" << str.View() << " is not an integer!" << std::endl;
return;
}
idx -= 1; // param start at 1, vector at 0
if (idx >= fRunInfo->GetNoOfParams()) { // idx is out-of-range
std::cerr << "PFitter::GetPhaseParams(): **ERROR** idx=" << idx << " is > #param = " << fRunInfo->GetNoOfParams() << "!" << std::endl;
return;
}
fPhase[idx] = true;
}
}
}
//--------------------------------------------------------------------------
// GetParFromFun (private)
//--------------------------------------------------------------------------
/**
* @brief Extracts parameter numbers from a FUNCTIONS block entry.
*
* Parses a function definition to find all parameters it references.
* Recursively handles nested map references within the function.
*
* @param funStr Function identifier (e.g., "fun1", "fun23")
* @return Vector of parameter numbers (1-indexed) used in the function
*
* @par Example:
* If FUNCTIONS block contains:
* @code
* fun1 = par3 * par7 + map2
* @endcode
* Then GetParFromFun("fun1") returns [3, 7, ...params from map2...]
*
* @note Returns empty vector if function not found or parsing fails.
*/
PIntVector PFitter::GetParFromFun(const TString funStr)
{
PIntVector parVec;
PMsrLines *funList = fRunInfo->GetMsrFunctions();
TString str;
for (int i=0; i<funList->size(); i++) {
if (funList->at(i).fLine.Contains(funStr)) {
// tokenize function string
std::vector<std::string> tok = PStringUtils::Split(funList->at(i).fLine.Data(), " =+-*/");
for (int j=1; j<static_cast<int>(tok.size()); j++) {
str = tok[j].c_str();
// parse tok for parX
if (str.Contains("par")) {
// find start idx of par in token
Ssiz_t idx = str.Index("par");
idx += 3;
TString parStr("");
do {
parStr += str[idx];
} while (isdigit(str[idx++]));
parVec.push_back(parStr.Atoi());
}
// parse tok for mapX
if (str.Contains("map")) {
// find start idx of par in token
Ssiz_t idx = str.Index("map");
idx += 3;
TString mapStr("map");
do {
mapStr += str[idx];
} while (isdigit(str[idx++]));
PIntVector mapParVec = GetParFromMap(mapStr);
for (int k=0; k<mapParVec.size(); k++) {
parVec.push_back(mapParVec[k]);
}
}
}
}
}
return parVec;
}
//--------------------------------------------------------------------------
// GetParFromMap (private)
//--------------------------------------------------------------------------
/**
* @brief Extracts parameter numbers from a map reference.
*
* Parses "mapX" to find which parameters are mapped to the X-th position
* across all RUN blocks. Maps allow different runs to use different parameters
* for the same theoretical component, enabling multi-run fits with run-dependent
* parameters.
*
* @param mapStr Map identifier (e.g., "map1", "map5")
* @return Vector of parameter numbers (1-indexed) mapped at this position
*
* @par Example:
* If RUN blocks have:
* @code
* RUN 1: map 7 8 9
* RUN 2: map 10 11 12
* @endcode
* Then GetParFromMap("map1") returns [7, 10] (first map entry per run)
*
* @note Returns empty vector if map index is out of range or invalid.
*/
PIntVector PFitter::GetParFromMap(const TString mapStr)
{
PIntVector parVec;
TString str = mapStr;
str.Remove(0,3); // remove map from string
int idx=str.Atoi();
if (idx == 0) {
std::cerr << "PFitter::GetParFromMap(): **ERROR** couldn't get propper index from mapX!" << std::endl;
return parVec;
}
idx -= 1; // map starts at 1, map vector at 0
// go through all the runs and collect the parameters from the map vectors
PMsrRunList *runList = fRunInfo->GetMsrRunList();
if (runList == nullptr) {
std::cerr << "PFitter::GetParFromMap(): **ERROR** couldn't get required run list information!" << std::endl;
return parVec;
}
PIntVector *map = nullptr;
for (int i=0; i<runList->size(); i++) {
map = runList->at(i).GetMap();
if (map == nullptr) {
std::cerr << "PFitter::GetParFromMap(): **ERROR** couldn't get required map information (idx=" << i << ")!" << std::endl;
parVec.clear();
return parVec;
}
if (idx >= map->size()) {
std::cerr << "PFitter::GetParFromMap(): **ERROR** requested map index (idx=" << idx << ") out-of-range (" << map->size() << ")!" << std::endl;
parVec.clear();
return parVec;
}
parVec.push_back(map->at(idx));
}
return parVec;
}
//--------------------------------------------------------------------------
// DoFit
//--------------------------------------------------------------------------
/**
* @brief Main entry point for executing the fit.
*
* This is the primary method that orchestrates the entire fitting process:
* 1. Transfers parameters from MSR to Minuit2
* 2. If chisq_only mode: calculates χ²/maxLH and returns
* 3. Otherwise: executes COMMANDS block sequentially
* 4. Updates MSR file with fit results and statistics
*
* @return true if fit completed successfully (or chisq calculated), false on errors
*
* @par Execution modes:
* - **chisq_only = true:** Evaluates objective function once with current parameters,
* reports χ²/maxLH and NDF, useful for validating parameter sets
* - **chisq_only = false:** Runs full fit according to COMMANDS block, which typically
* includes MIGRAD (minimization), HESSE (error matrix), optionally MINOS (asymmetric errors)
*
* @par COMMANDS execution:
* Commands are executed in the order they appear in the MSR file. Common sequence:
* @code
* SET PRINT 1 # Set verbosity
* MIGRAD # Find minimum
* HESSE # Calculate symmetric errors
* MINOS # Calculate asymmetric errors (optional)
* SAVE # Save results to MSR file
* @endcode
*
* @note Check IsValid() before calling. Check HasConverged() after completion.
* @note Updates fRunInfo with final parameters, errors, and fit statistics.
* @note Prints detailed χ² breakdown per run if multiple runs are fitted.
*
* @see ExecuteMigrad(), ExecuteHesse(), ExecuteMinos(), ExecuteSave()
*/
Bool_t PFitter::DoFit()
{
// feed minuit parameters
SetParameters();
// check if only chisq/maxLH shall be calculated once
if (fChisqOnly) {
std::vector<Double_t> param = fMnUserParams.Params();
std::vector<Double_t> error = fMnUserParams.Errors();
Int_t usedParams = 0;
for (UInt_t i=0; i<error.size(); i++) {
if (error[i] != 0.0)
usedParams++;
}
UInt_t ndf = static_cast<int>(fFitterFcn->GetTotalNoOfFittedBins()) - usedParams;
Double_t val = (*fFitterFcn)(param);
if (fUseChi2) {
// calculate expected chisq
Double_t totalExpectedChisq = 0.0;
PDoubleVector expectedChisqPerRun;
fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerRun);
// calculate chisq per run
std::vector<Double_t> chisqPerRun;
for (UInt_t i=0; i<fRunInfo->GetMsrRunList()->size(); i++) {
chisqPerRun.push_back(fRunListCollection->GetSingleRunChisq(param, i));
}
std::cout << std::endl << std::endl << ">> 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_run = 0;
for (UInt_t i=0; i<expectedChisqPerRun.size(); i++) {
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_run = 0;
for (UInt_t i=0; i<chisqPerRun.size(); i++) {
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
chisqPerRun.clear();
expectedChisqPerRun.clear();
if (fSectorFlag) {
PDoublePairVector secFitRange;
secFitRange.resize(1);
for (UInt_t k=0; k<fSector.size(); k++) {
// set sector fit range
secFitRange[0].first = fSector[k].GetTimeRangeFirst(0);
secFitRange[0].second = fSector[k].GetTimeRangeLast();
fRunListCollection->SetFitRange(secFitRange);
// calculate chisq
val = (*fFitterFcn)(param);
// calculate NDF
ndf = static_cast<UInt_t>(fFitterFcn->GetTotalNoOfFittedBins()) - usedParams;
// calculate expected chisq
totalExpectedChisq = 0.0;
fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerRun);
// calculate chisq per run
for (UInt_t i=0; i<fRunInfo->GetMsrRunList()->size(); i++) {
chisqPerRun.push_back(fRunListCollection->GetSingleRunChisq(param, i));
}
std::cout << std::endl;
std::cout << "++++" << 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_run = 0;
for (UInt_t i=0; i<expectedChisqPerRun.size(); i++) {
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_run = 0;
for (UInt_t i=0; i<chisqPerRun.size(); i++) {
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
chisqPerRun.clear();
expectedChisqPerRun.clear();
}
}
} else { // max. log likelihood
// calculate expected maxLH
Double_t totalExpectedMaxLH = 0.0;
std::vector<Double_t> expectedMaxLHPerRun;
fFitterFcn->CalcExpectedChiSquare(param, totalExpectedMaxLH, expectedMaxLHPerRun);
// calculate maxLH per run
std::vector<Double_t> maxLHPerRun;
for (UInt_t i=0; i<fRunInfo->GetMsrRunList()->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_run = 0;
for (UInt_t i=0; i<expectedMaxLHPerRun.size(); i++) {
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 << ")";
}
}
// clean up
maxLHPerRun.clear();
expectedMaxLHPerRun.clear();
if (fSectorFlag) {
PDoublePairVector secFitRange;
secFitRange.resize(1);
for (UInt_t k=0; k<fSector.size(); k++) {
// set sector fit range
secFitRange[0].first = fSector[k].GetTimeRangeFirst(0);
secFitRange[0].second = fSector[k].GetTimeRangeLast();
fRunListCollection->SetFitRange(secFitRange);
// calculate maxLH
val = (*fFitterFcn)(param);
// calculate NDF
ndf = static_cast<int>(fFitterFcn->GetTotalNoOfFittedBins()) - usedParams;
// calculate expected maxLH
totalExpectedMaxLH = 0.0;
fFitterFcn->CalcExpectedChiSquare(param, totalExpectedMaxLH, expectedMaxLHPerRun);
// calculate maxLH per run
for (UInt_t i=0; i<fRunInfo->GetMsrRunList()->size(); i++) {
maxLHPerRun.push_back(fRunListCollection->GetSingleRunMaximumLikelihood(param, i));
}
std::cout << std::endl;
std::cout << "++++" << 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_run = 0;
for (UInt_t i=0; i<expectedMaxLHPerRun.size(); i++) {
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 << ")";
}
}
// clean up
maxLHPerRun.clear();
expectedMaxLHPerRun.clear();
}
}
}
std::cout << std::endl << std::endl;
return true;
}
// debugging information
#ifdef HAVE_GOMP
std::cout << std::endl << ">> Number of available threads for the function optimization: " << omp_get_max_threads() << std::endl;
#endif
// real fit wanted
if (fUseChi2)
std::cout << std::endl << ">> Chi Square fit will be executed" << std::endl;
else
std::cout << std::endl << ">> Maximum Likelihood fit will be executed" << std::endl;
Bool_t status = true;
// init positive errors to default false, if minos is called, it will be set true there
for (UInt_t i=0; i<fParams.size(); i++) {
fRunInfo->SetMsrParamPosErrorPresent(i, false);
}
// walk through the command list and execute them
Bool_t firstSave = true;
for (UInt_t i=0; i<fCmdList.size(); i++) {
switch (fCmdList[i].first) {
case PMN_INTERACTIVE:
std::cerr << std::endl << "**WARNING** from PFitter::DoFit() : the command INTERACTIVE is not yet implemented.";
std::cerr << std::endl;
break;
case PMN_CONTOURS:
status = ExecuteContours();
break;
case PMN_FIT_RANGE:
status = ExecuteFitRange(fCmdList[i].second);
break;
case PMN_FIX:
status = ExecuteFix(fCmdList[i].second);
break;
case PMN_EIGEN:
std::cerr << std::endl << "**WARNING** from PFitter::DoFit() : the command EIGEN is not yet implemented.";
std::cerr << std::endl;
break;
case PMN_HESSE:
status = ExecuteHesse();
break;
case PMN_MACHINE_PRECISION:
std::cerr << std::endl << "**WARNING** from PFitter::DoFit() : the command MACHINE_PRECISION is not yet implemented.";
std::cerr << std::endl;
break;
case PMN_MIGRAD:
status = ExecuteMigrad();
break;
case PMN_MINIMIZE:
status = ExecuteMinimize();
break;
case PMN_MINOS:
status = ExecuteMinos();
break;
case PMN_PLOT:
status = ExecutePlot();
break;
case PMN_RELEASE:
status = ExecuteRelease(fCmdList[i].second);
break;
case PMN_RESTORE:
status = ExecuteRestore();
break;
case PMN_SAVE:
status = ExecuteSave(firstSave);
if (firstSave)
firstSave = false;
break;
case PMN_SCAN:
status = ExecuteScan();
break;
case PMN_SECTOR:
// nothing to be done here
break;
case PMN_SIMPLEX:
status = ExecuteSimplex();
break;
case PMN_USER_COVARIANCE:
std::cerr << std::endl << "**WARNING** from PFitter::DoFit() : the command USER_COVARIANCE is not yet implemented.";
std::cerr << std::endl;
break;
case PMN_USER_PARAM_STATE:
std::cerr << std::endl << "**WARNING** from PFitter::DoFit() : the command USER_PARAM_STATE is not yet implemented.";
std::cerr << std::endl;
break;
case PMN_PRINT:
status = ExecutePrintLevel(fCmdList[i].second);
break;
default:
std::cerr << std::endl << "**PANIC ERROR**: PFitter::DoFit(): You should never have reached this point";
std::cerr << std::endl;
exit(0);
}
// check if command has been successful
if (!status)
break;
}
if (IsValid()) {
fRunInfo->GetMsrStatistic()->fValid = true;
} else {
fRunInfo->GetMsrStatistic()->fValid = false;
}
return true;
}
//--------------------------------------------------------------------------
// CheckCommands
//--------------------------------------------------------------------------
/**
* @brief Validates COMMANDS block and builds execution queue.
*
* Parses all command lines from the MSR file's COMMANDS block, validates
* syntax and parameters, and constructs an ordered execution list. This
* ensures commands are executable before starting the fit.
*
* @return true if all commands are valid and parseable, false on any error
*
* @par Validation checks:
* - Command keywords are recognized (MIGRAD, HESSE, MINOS, etc.)
* - Numeric arguments are valid (parameter indices, ranges, point counts)
* - Parameter references are within bounds
* - SCAN/CONTOURS have required arguments
* - FIX/RELEASE specify valid parameter lists
*
* @par Side effects:
* - Populates fCmdList with (command_ID, line_number) pairs
* - Sets fIsValid = false if any command is invalid
* - Configures scan parameters (fScanParameter, fScanNoPoints, etc.)
* - Detects SECTOR command presence (sets fSectorFlag)
*
* @note Called automatically by constructor. Errors are reported to stderr.
* @note Some legacy commands (SET BATCH, END RETURN) are silently ignored.
*
* @see PMN_* command constants in PFitter.h
*/
Bool_t PFitter::CheckCommands()
{
fIsValid = true;
// check if chisq or log max likelihood fit
fUseChi2 = fRunInfo->GetMsrStatistic()->fChisq;
// walk through the msr-file COMMAND block
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++;
// 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 (line.Contains("SET BATCH", TString::kIgnoreCase)) { // needed for backward compatibility
continue;
} else if (line.Contains("END RETURN", TString::kIgnoreCase)) { // needed for backward compatibility
continue;
} else if (line.Contains("CHI_SQUARE", TString::kIgnoreCase)) {
continue;
} else if (line.Contains("MAX_LIKELIHOOD", TString::kIgnoreCase)) {
continue;
} else if (line.Contains("SCALE_N0_BKG", TString::kIgnoreCase)) {
continue;
} else if (line.Contains("INTERACTIVE", TString::kIgnoreCase)) {
cmd.first = PMN_INTERACTIVE;
cmd.second = cmdLineNo;
fCmdList.push_back(cmd);
} else if (line.Contains("CONTOURS", TString::kIgnoreCase)) {
cmd.first = PMN_CONTOURS;
cmd.second = cmdLineNo;
fCmdList.push_back(cmd);
// filter out possible parameters for scan
std::vector<std::string> tokens;
TString str;
UInt_t ival;
tokens = PStringUtils::Split(line.Data(), ", \t");
for (Int_t i=0; i<static_cast<int>(tokens.size()); i++) {
str = tokens[i].c_str();
if ((i==1) || (i==2)) { // parX / parY
// 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 << ">> " << 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;
fIsValid = false;
break;
}
ival = str.Atoi();
// 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 << ">> " << 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;
fIsValid = false;
break;
}
// keep parameter
fScanParameter[i-1] = ival-1; // internally parameter number starts at 0
fScanAll = false;
}
if (i==3) {
// 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 << ">> " << 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;
fIsValid = false;
break;
}
ival = str.Atoi();
if ((ival < 1) || (ival > 100)) {
std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
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;
fIsValid = false;
break;
}
fScanNoPoints = ival;
}
}
} else if (line.Contains("EIGEN", TString::kIgnoreCase)) {
cmd.first = PMN_EIGEN;
cmd.second = cmdLineNo;
fCmdList.push_back(cmd);
} else if (line.Contains("FIT_RANGE", TString::kIgnoreCase)) {
// check the 5 options:
// (i) FIT_RANGE RESET,
// (ii) FIT_RANGE start end,
// (iii) FIT_RANGE start1 end1 start2 end2 ... startN endN
// (iv) FIT_RANGE fgb+n0 lgb-n1
// (v) FIT_RANGE fgb+n00 lgb-n01 fgb+n10 lgb-n11 ... fgb+nN0 lgb-nN1
std::vector<std::string> tokens;
TString str;
tokens = PStringUtils::Split(line.Data(), ", \t");
if (static_cast<int>(tokens.size()) == 2) { // should only be RESET
str = tokens[1].c_str();
if (str.Contains("RESET", TString::kIgnoreCase)) {
cmd.first = PMN_FIT_RANGE;
cmd.second = cmdLineNo;
fCmdList.push_back(cmd);
} else {
std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
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;
fIsValid = false;
break;
}
} else if ((static_cast<int>(tokens.size()) > 1) && (static_cast<UInt_t>(static_cast<int>(tokens.size())) % 2) == 1) {
if ((static_cast<int>(tokens.size()) > 3) && ((static_cast<UInt_t>(static_cast<int>(tokens.size()))-1)) != 2*fRunInfo->GetMsrRunList()->size()) {
std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
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 << ">> FIT_RANGE fgb+<n0> lgb-<n1> | FIT_RANGE fgb+<n00> lgb-<n01> fgb+<n10> lgb-<n11> ... fgb+<nN0> lgb-<nN1>,";
std::cerr << std::endl << ">> with N the number of runs in the msr-file.";
std::cerr << std::endl << ">> Found N=" << (static_cast<int>(tokens.size())-1)/2 << ", # runs in msr-file=" << fRunInfo->GetMsrRunList()->size() << std::endl;
fIsValid = false;
break;
} else {
// check that all range entries are numbers or fgb+n0 / lgb-n1
Bool_t ok = true;
for (Int_t n=1; n<static_cast<int>(tokens.size()); n++) {
str = tokens[n].c_str();
if (!str.IsFloat()) {
if ((n%2 == 1) && (!str.Contains("fgb", TString::kIgnoreCase)))
ok = false;
if ((n%2 == 0) && (!str.Contains("lgb", TString::kIgnoreCase)))
ok = false;
}
if (!ok)
break;
}
if (ok) { // everything is fine
cmd.first = PMN_FIT_RANGE;
cmd.second = cmdLineNo;
fCmdList.push_back(cmd);
} else {
std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
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 << ">> FIT_RANGE fgb+<n0> lgb-<n1> | FIT_RANGE fgb+<n00> lgb-<n01> fgb+<n10> lgb-<n11> ... fgb+<nN0> lgb-<nN1>,";
std::cerr << std::endl << ">> with N the number of runs in the msr-file.";
std::cerr << std::endl << ">> Found token '" << str.Data() << "', which is not a floating point number." << std::endl;
fIsValid = false;
break;
}
}
} else {
std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
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 << ">> FIT_RANGE fgb+<n0> lgb-<n1> | FIT_RANGE fgb+<n00> lgb-<n01> fgb+<n10> lgb-<n11> ... fgb+<nN0> lgb-<nN1>,";
std::cerr << std::endl << ">> with N the number of runs in the msr-file.";
fIsValid = false;
break;
}
} else if (line.Contains("FIX", TString::kIgnoreCase)) {
// check if the given set of parameters (number or names) is present
std::vector<std::string> tokens;
TString str;
UInt_t ival;
tokens = PStringUtils::Split(line.Data(), ", \t");
for (Int_t i=1; i<static_cast<int>(tokens.size()); i++) {
str = tokens[i].c_str();
if (str.IsDigit()) { // token might be a parameter number
ival = str.Atoi();
// 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 << ">> " << line.Data();
std::cerr << std::endl << ">> Parameter " << ival << " is out of the Parameter Range [1," << fParams.size() << "]";
std::cerr << std::endl;
fIsValid = false;
break;
}
} else { // token might be a parameter name
// check if token is present as parameter name
Bool_t found = false;
for (UInt_t j=0; j<fParams.size(); j++) {
if (fParams[j].fName.CompareTo(str, TString::kIgnoreCase) == 0) { // found
found = true;
break;
}
}
if (!found) {
std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
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;
break;
}
}
}
// everything looks fine, feed the command list
cmd.first = PMN_FIX;
cmd.second = cmdLineNo;
fCmdList.push_back(cmd);
} else if (line.Contains("HESSE", TString::kIgnoreCase)) {
fIsScanOnly = false;
cmd.first = PMN_HESSE;
cmd.second = cmdLineNo;
fCmdList.push_back(cmd);
} else if (line.Contains("MACHINE_PRECISION", TString::kIgnoreCase)) {
cmd.first = PMN_MACHINE_PRECISION;
cmd.second = cmdLineNo;
fCmdList.push_back(cmd);
} else if (line.Contains("MIGRAD", TString::kIgnoreCase)) {
fIsScanOnly = false;
cmd.first = PMN_MIGRAD;
cmd.second = cmdLineNo;
fCmdList.push_back(cmd);
} else if (line.Contains("MINIMIZE", TString::kIgnoreCase)) {
fIsScanOnly = false;
cmd.first = PMN_MINIMIZE;
cmd.second = cmdLineNo;
fCmdList.push_back(cmd);
} else if (line.Contains("MINOS", TString::kIgnoreCase)) {
fIsScanOnly = false;
cmd.first = PMN_MINOS;
cmd.second = cmdLineNo;
fCmdList.push_back(cmd);
} else if (line.Contains("MNPLOT", TString::kIgnoreCase)) {
cmd.first = PMN_PLOT;
cmd.second = cmdLineNo;
fCmdList.push_back(cmd);
} else if (line.Contains("PRINT_LEVEL", TString::kIgnoreCase)) {
cmd.first = PMN_PRINT;
cmd.second = cmdLineNo;
fCmdList.push_back(cmd);
} else if (line.Contains("RELEASE", TString::kIgnoreCase)) {
// check if the given set of parameters (number or names) is present
std::vector<std::string> tokens;
TString str;
UInt_t ival;
tokens = PStringUtils::Split(line.Data(), ", \t");
for (Int_t i=1; i<static_cast<int>(tokens.size()); i++) {
str = tokens[i].c_str();
if (str.IsDigit()) { // token might be a parameter number
ival = str.Atoi();
// 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 << ">> " << line.Data();
std::cerr << std::endl << ">> Parameter " << ival << " is out of the Parameter Range [1," << fParams.size() << "]";
std::cerr << std::endl;
fIsValid = false;
break;
}
} else { // token might be a parameter name
// check if token is present as parameter name
Bool_t found = false;
for (UInt_t j=0; j<fParams.size(); j++) {
if (fParams[j].fName.CompareTo(str, TString::kIgnoreCase) == 0) { // found
found = true;
break;
}
}
if (!found) {
std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
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;
break;
}
}
}
cmd.first = PMN_RELEASE;
cmd.second = cmdLineNo;
fCmdList.push_back(cmd);
} else if (line.Contains("RESTORE", TString::kIgnoreCase)) {
cmd.first = PMN_RESTORE;
cmd.second = cmdLineNo;
fCmdList.push_back(cmd);
} else if (line.Contains("SAVE", TString::kIgnoreCase)) {
cmd.first = PMN_SAVE;
cmd.second = cmdLineNo;
fCmdList.push_back(cmd);
} else if (line.Contains("SCAN", TString::kIgnoreCase)) {
cmd.first = PMN_SCAN;
cmd.second = cmdLineNo;
fCmdList.push_back(cmd);
// filter out possible parameters for scan
std::vector<std::string> tokens;
TString str;
UInt_t ival;
tokens = PStringUtils::Split(line.Data(), ", \t");
for (Int_t i=0; i<static_cast<int>(tokens.size()); i++) {
str = tokens[i].c_str();
if (i==1) { // get parameter number
// 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 << ">> " << 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;
fIsValid = false;
break;
}
ival = str.Atoi();
// 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 << ">> " << 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;
fIsValid = false;
break;
}
// keep parameter
fScanParameter[0] = ival-1; // internally parameter number starts at 0
fScanAll = false;
}
if (i==2) { // get number of points
// 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 << ">> " << 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;
fIsValid = false;
break;
}
ival = str.Atoi();
if ((ival < 1) || (ival > 100)) {
std::cerr << std::endl << ">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
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;
fIsValid = false;
break;
}
fScanNoPoints = ival;
}
if (i==3) { // get low
// 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 << ">> " << 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;
fIsValid = false;
break;
}
fScanLow = str.Atof();
}
if (i==4) { // get high
// 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 << ">> " << 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;
fIsValid = false;
break;
}
fScanHigh = str.Atof();
}
}
} else if (line.Contains("SIMPLEX", TString::kIgnoreCase)) {
cmd.first = PMN_SIMPLEX;
cmd.second = cmdLineNo;
fCmdList.push_back(cmd);
} else if (line.Contains("STRATEGY", TString::kIgnoreCase)) {
std::vector<std::string> tokens;
TString str;
tokens = PStringUtils::Split(line.Data(), " \t");
if (static_cast<int>(tokens.size()) == 2) {
str = tokens[1].c_str();
if (str.CompareTo("0") == 0) { // low
fStrategy = 0;
} else if (str.CompareTo("1") == 0) { // default
fStrategy = 1;
} else if (str.CompareTo("2") == 0) { // high
fStrategy = 2;
} else if (str.CompareTo("LOW") == 0) { // low
fStrategy = 0;
} else if (str.CompareTo("DEFAULT") == 0) { // default
fStrategy = 1;
} else if (str.CompareTo("HIGH") == 0) { // high
fStrategy = 2;
}
}
} else if (line.Contains("USER_COVARIANCE", TString::kIgnoreCase)) {
cmd.first = PMN_USER_COVARIANCE;
cmd.second = cmdLineNo;
fCmdList.push_back(cmd);
} 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)) {
fSectorFlag = true;
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
std::vector<std::string> tokens;
TString str;
tokens = PStringUtils::Split(line.Data(), " ,\t");
if (static_cast<int>(tokens.size()) == 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;
fIsValid = false;
fSectorFlag = false;
break;
}
Double_t dval;
for (Int_t i=1; i<static_cast<int>(tokens.size()); i++) {
// keep time range of sector
PSectorChisq sec(fRunInfo->GetNoOfRuns());
// get parse tokens
str = tokens[i].c_str();
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.size(); j++) {
if (dval > 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;
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();
std::cerr << std::endl << ">> The sector time stamp '" << str << "' is not a number.";
std::cerr << std::endl << ">> Will stop ...";
std::cerr << std::endl;
fIsValid = false;
fSectorFlag = false;
break;
}
}
} else { // unkown command
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;
}
}
// Check that in case release/restore is present, that it is followed by a minimizer before minos is called.
// If this is not the case, place a warning
Bool_t fixFlag = false;
Bool_t releaseFlag = false;
Bool_t minimizerFlag = false;
for (it = fCmdLines.begin(); it != fCmdLines.end(); ++it) {
if (line.Contains("FIX", TString::kIgnoreCase))
fixFlag = true;
else if (line.Contains("RELEASE", TString::kIgnoreCase) ||
line.Contains("RESTORE", TString::kIgnoreCase))
releaseFlag = true;
else if (line.Contains("MINIMIZE", TString::kIgnoreCase) ||
line.Contains("MIGRAD", TString::kIgnoreCase) ||
line.Contains("SIMPLEX", TString::kIgnoreCase)) {
if (releaseFlag)
minimizerFlag = true;
} 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";
std::cerr << std::endl << ">> RELEASE/RESTORE and MINOS. Behaviour might be different to the";
std::cerr << std::endl << ">> expectation of the user ?!?" << std::endl;
}
fixFlag = false;
releaseFlag = false;
minimizerFlag = false;
}
}
return fIsValid;
}
//--------------------------------------------------------------------------
// SetParameters
//--------------------------------------------------------------------------
/**
* <p>Feeds the internal minuit2 fit parameters. It also makes sure that unused parameters
* are fixed.
*
* <b>return:</b> true.
*/
Bool_t PFitter::SetParameters()
{
for (UInt_t i=0; i<fParams.size(); i++) {
// check if parameter is fixed
if (fParams[i].fStep == 0.0) { // add fixed parameter
fMnUserParams.Add(fParams[i].fName.Data(), fParams[i].fValue);
} else { // add free parameter
// check if boundaries are given
if (fParams[i].fNoOfParams > 5) { // boundaries given
if (fParams[i].fLowerBoundaryPresent &&
fParams[i].fUpperBoundaryPresent) { // upper and lower boundaries given
fMnUserParams.Add(fParams[i].fName.Data(), fParams[i].fValue, fParams[i].fStep,
fParams[i].fLowerBoundary, fParams[i].fUpperBoundary);
} else if (fParams[i].fLowerBoundaryPresent &&
!fParams[i].fUpperBoundaryPresent) { // lower boundary limited
fMnUserParams.Add(fParams[i].fName.Data(), fParams[i].fValue, fParams[i].fStep);
fMnUserParams.SetLowerLimit(fParams[i].fName.Data(), fParams[i].fLowerBoundary);
} else { // upper boundary limited
fMnUserParams.Add(fParams[i].fName.Data(), fParams[i].fValue, fParams[i].fStep);
fMnUserParams.SetUpperLimit(fParams[i].fName.Data(), fParams[i].fUpperBoundary);
}
} else { // no boundaries given
fMnUserParams.Add(fParams[i].fName.Data(), fParams[i].fValue, fParams[i].fStep);
}
}
}
// check if there is an unused parameter, if so, fix it
for (UInt_t i=0; i<fParams.size(); i++) {
// parameter not used in the whole theory and not already fixed!!
if ((fRunInfo->ParameterInUse(i) == 0) && (fParams[i].fStep != 0.0)) {
fMnUserParams.Fix(i); // fix the unused parameter so that minuit will not vary it
std::cerr << std::endl << "**WARNING** : Parameter No " << i+1 << " is not used at all, will fix it";
std::cerr << std::endl;
}
}
return true;
}
//--------------------------------------------------------------------------
// ExecuteContours
//--------------------------------------------------------------------------
/**
* <p>Execute the minuit2 contour command. Makes sure that a valid minuit2 minimum is present.
*
* <b>return:</b> true if the contour command could be executed successfully, otherwise returns false.
*/
Bool_t PFitter::ExecuteContours()
{
std::cout << ">> PFitter::ExecuteContours() ..." << std::endl;
// if already some minimization is done use the minuit2 output as input
if (!fFcnMin) {
std::cerr << std::endl << "**WARNING**: CONTOURS musn't be called before any minimization (MINIMIZE/MIGRAD/SIMPLEX) is done!!";
std::cerr << std::endl;
return false;
}
// check if minimum was valid
if (!fFcnMin->IsValid()) {
std::cerr << std::endl << "**ERROR**: CONTOURS cannot started since the previous minimization failed :-(";
std::cerr << std::endl;
return false;
}
ROOT::Minuit2::MnContours contours((*fFitterFcn), *fFcnMin);
fScanData = contours(fScanParameter[0], fScanParameter[1], fScanNoPoints);
return true;
}
//--------------------------------------------------------------------------
// ExecuteFitRange
//--------------------------------------------------------------------------
/**
* <p>Change the fit range via command block.
*
* \param lineNo the line number of the command block
*
* <b>return:</b> true if done, otherwise returns false.
*/
Bool_t PFitter::ExecuteFitRange(UInt_t lineNo)
{
std::cout << ">> PFitter::ExecuteFitRange(): " << fCmdLines[lineNo].fLine.Data() << std::endl;
if (fCmdLines[lineNo].fLine.Contains("fgb", TString::kIgnoreCase)) { // fit range given in bins
fRunListCollection->SetFitRange(fCmdLines[lineNo].fLine);
return true;
}
std::vector<std::string> tokens;
TString str;
tokens = PStringUtils::Split(fCmdLines[lineNo].fLine.Data(), ", \t");
PMsrRunList *runList = fRunInfo->GetMsrRunList();
// execute command, no error checking needed since this has been already carried out in CheckCommands()
if (static_cast<int>(tokens.size()) == 2) { // reset command
fRunListCollection->SetFitRange(fOriginalFitRange);
} else if (static_cast<int>(tokens.size()) == 3) { // single fit range for all runs
Double_t start = 0.0, end = 0.0;
PDoublePair fitRange;
PDoublePairVector fitRangeVector;
str = tokens[1].c_str();
start = str.Atof();
str = tokens[2].c_str();
end = str.Atof();
fitRange.first = start;
fitRange.second = end;
fitRangeVector.push_back(fitRange);
fRunListCollection->SetFitRange(fitRangeVector);
} else { // individual fit ranges for each run
Double_t start = 0.0, end = 0.0;
PDoublePair fitRange;
PDoublePairVector fitRangeVector;
for (UInt_t i=0; i<runList->size(); i++) {
str = tokens[2*i+1].c_str();
start = str.Atof();
str = tokens[2*i+2].c_str();
end = str.Atof();
fitRange.first = start;
fitRange.second = end;
fitRangeVector.push_back(fitRange);
}
fRunListCollection->SetFitRange(fitRangeVector);
}
return true;
}
//--------------------------------------------------------------------------
// ExecuteFix
//--------------------------------------------------------------------------
/**
* <p>Fix parameter list given at lineNo of the command block.
*
* \param lineNo the line number of the command block
*
* <b>return:</b> true if done, otherwise returns false.
*/
Bool_t PFitter::ExecuteFix(UInt_t lineNo)
{
std::cout << ">> PFitter::ExecuteFix(): " << fCmdLines[lineNo].fLine.Data() << std::endl;
std::vector<std::string> tokens;
TString str;
tokens = PStringUtils::Split(fCmdLines[lineNo].fLine.Data(), ", \t");
for (Int_t i=1; i<static_cast<int>(tokens.size()); i++) {
str = tokens[i].c_str();
if (str.IsDigit()) { // token is a parameter number
fMnUserParams.Fix(static_cast<UInt_t>(str.Atoi())-1);
} else { // token is a parameter name
fMnUserParams.Fix(str.Data());
}
}
return true;
}
//--------------------------------------------------------------------------
// ExecuteHesse
//--------------------------------------------------------------------------
/**
* <p>Execute the minuit2 hesse command.
*
* <b>return:</b> true if the hesse command could be executed successfully, otherwise returns false.
*/
Bool_t PFitter::ExecuteHesse()
{
std::cout << ">> PFitter::ExecuteHesse(): will call hesse ..." << std::endl;
// create the hesse object
ROOT::Minuit2::MnHesse hesse;
// specify maximal number of function calls
UInt_t maxfcn = std::numeric_limits<UInt_t>::max();
// call hesse
Double_t start=0.0, end=0.0;
start=MilliTime();
ROOT::Minuit2::MnUserParameterState mnState = hesse((*fFitterFcn), fMnUserParams, maxfcn);
end=MilliTime();
std::cout << ">> PFitter::ExecuteMinimize(): execution time for Hesse = " << std::setprecision(3) << (end-start)/1.0e3 << " sec." << std::endl;
TString str = TString::Format("Hesse: %.3f sec", (end-start)/1.0e3);
fElapsedTime.push_back(str);
if (!mnState.IsValid()) {
std::cerr << std::endl << ">> PFitter::ExecuteHesse(): **WARNING** Hesse encountered a problem! The state found is invalid.";
std::cerr << std::endl;
return false;
}
if (!mnState.HasCovariance()) {
std::cerr << std::endl << ">> PFitter::ExecuteHesse(): **WARNING** Hesse encountered a problem! No covariance matrix available.";
std::cerr << std::endl;
return false;
}
// fill parabolic errors
for (UInt_t i=0; i<fParams.size(); i++) {
fRunInfo->SetMsrParamStep(i, mnState.Error(i));
fRunInfo->SetMsrParamPosErrorPresent(i, false);
}
if (fPrintLevel >= 2)
std::cout << mnState << std::endl;
return true;
}
//--------------------------------------------------------------------------
// ExecuteMigrad
//--------------------------------------------------------------------------
/**
* <p>Execute the minuit2 migrad command.
*
* <b>return:</b> true if the migrad command could be executed successfully, otherwise returns false.
*/
Bool_t PFitter::ExecuteMigrad()
{
std::cout << ">> PFitter::ExecuteMigrad(): will call migrad ..." << std::endl;
// create migrad object
// strategy is by default = 'default'
ROOT::Minuit2::MnMigrad migrad((*fFitterFcn), fMnUserParams, ROOT::Minuit2::MnStrategy{fStrategy});
// minimize
// maxfcn is MINUIT2 Default maxfcn
UInt_t maxfcn = std::numeric_limits<UInt_t>::max();
// tolerance = MINUIT2 Default tolerance
Double_t tolerance = 0.1;
// keep track of elapsed time
Double_t start=0.0, end=0.0;
start=MilliTime();
ROOT::Minuit2::FunctionMinimum min = migrad(maxfcn, tolerance);
end=MilliTime();
std::cout << ">> PFitter::ExecuteMinimize(): execution time for Migrad = " << std::setprecision(3) << (end-start)/1.0e3 << " sec." << std::endl;
TString str = TString::Format("Migrad: %.3f sec", (end-start)/1.0e3);
fElapsedTime.push_back(str);
if (!min.IsValid()) {
std::cerr << std::endl << ">> PFitter::ExecuteMigrad(): **WARNING**: Fit did not converge, sorry ...";
std::cerr << std::endl;
fIsValid = false;
return false;
}
// keep FunctionMinimum object
fFcnMin.reset(new ROOT::Minuit2::FunctionMinimum(min));
// keep user parameters
if (fFcnMin)
fMnUserParams = fFcnMin->UserParameters();
// fill run info
for (UInt_t i=0; i<fParams.size(); i++) {
Double_t dval = min.UserState().Value(i);
if (fPhase[i]) {
Int_t m = (Int_t)(dval/360.0);
dval = dval - m*360.0;
}
fRunInfo->SetMsrParamValue(i, dval);
fRunInfo->SetMsrParamStep(i, min.UserState().Error(i));
fRunInfo->SetMsrParamPosErrorPresent(i, false);
}
// handle statistics
Double_t minVal = min.Fval();
UInt_t ndf = fFitterFcn->GetTotalNoOfFittedBins();
// subtract number of varied parameters from total no of fitted bins -> ndf
for (UInt_t i=0; i<fParams.size(); i++) {
if ((min.UserState().Error(i) != 0.0) && !fMnUserParams.Parameters().at(i).IsFixed())
ndf -= 1;
}
// feed run info with new statistics info
fRunInfo->SetMsrStatisticMin(minVal);
fRunInfo->SetMsrStatisticNdf(ndf);
fConverged = true;
if (fPrintLevel >= 2)
std::cout << *fFcnMin << std::endl;
return true;
}
//--------------------------------------------------------------------------
// ExecuteMinimize
//--------------------------------------------------------------------------
/**
* <p>Execute the minuit2 minimize command.
*
* <b>return:</b> true if the minimize command could be executed successfully, otherwise returns false.
*/
Bool_t PFitter::ExecuteMinimize()
{
std::cout << ">> PFitter::ExecuteMinimize(): will call minimize ..." << std::endl;
// create minimizer object
// strategy is by default = 'default'
ROOT::Minuit2::MnMinimize minimize((*fFitterFcn), fMnUserParams, ROOT::Minuit2::MnStrategy{fStrategy});
// minimize
// maxfcn is MINUIT2 Default maxfcn
UInt_t maxfcn = std::numeric_limits<UInt_t>::max();
// tolerance = MINUIT2 Default tolerance
Double_t tolerance = 0.1;
// keep track of elapsed time
Double_t start=0.0, end=0.0;
start = MilliTime();
ROOT::Minuit2::FunctionMinimum min = minimize(maxfcn, tolerance);
end = MilliTime();
std::cout << ">> PFitter::ExecuteMinimize(): execution time for Minimize = " << std::setprecision(3) << (end-start)/1.0e3 << " sec." << std::endl;
TString str = TString::Format("Minimize: %.3f sec", (end-start)/1.0e3);
fElapsedTime.push_back(str);
if (!min.IsValid()) {
std::cerr << std::endl << ">> PFitter::ExecuteMinimize(): **WARNING**: Fit did not converge, sorry ...";
std::cerr << std::endl;
fIsValid = false;
return false;
}
// keep FunctionMinimum object
fFcnMin.reset(new ROOT::Minuit2::FunctionMinimum(min));
// keep user parameters
if (fFcnMin)
fMnUserParams = fFcnMin->UserParameters();
// fill run info
for (UInt_t i=0; i<fParams.size(); i++) {
Double_t dval = min.UserState().Value(i);
if (fPhase[i]) {
Int_t m = (Int_t)(dval/360.0);
dval = dval - m*360.0;
}
fRunInfo->SetMsrParamValue(i, dval);
fRunInfo->SetMsrParamStep(i, min.UserState().Error(i));
fRunInfo->SetMsrParamPosErrorPresent(i, false);
}
// handle statistics
Double_t minVal = min.Fval();
UInt_t ndf = fFitterFcn->GetTotalNoOfFittedBins();
// subtract number of varied parameters from total no of fitted bins -> ndf
for (UInt_t i=0; i<fParams.size(); i++) {
if ((min.UserState().Error(i) != 0.0) && !fMnUserParams.Parameters().at(i).IsFixed())
ndf -= 1;
}
// feed run info with new statistics info
fRunInfo->SetMsrStatisticMin(minVal);
fRunInfo->SetMsrStatisticNdf(ndf);
fConverged = true;
if (fPrintLevel >= 2)
std::cout << *fFcnMin << std::endl;
return true;
}
//--------------------------------------------------------------------------
// ExecuteMinos
//--------------------------------------------------------------------------
/**
* <p>Execute the minuit2 minos command.
*
* <b>return:</b> true if the minos command could be executed successfully, otherwise returns false.
*/
Bool_t PFitter::ExecuteMinos()
{
std::cout << ">> PFitter::ExecuteMinos(): will call minos ..." << std::endl;
// if already some minimization is done use the minuit2 output as input
if (!fFcnMin) {
std::cerr << std::endl << "**ERROR**: MINOS musn't be called before any minimization (MINIMIZE/MIGRAD/SIMPLEX) is done!!";
std::cerr << std::endl;
return false;
}
// check if minimum was valid
if (!fFcnMin->IsValid()) {
std::cerr << std::endl << "**ERROR**: MINOS cannot started since the previous minimization failed :-(";
std::cerr << std::endl;
return false;
}
// make minos analysis
Double_t start=0.0, end=0.0;
start=MilliTime();
ROOT::Minuit2::MnMinos minos((*fFitterFcn), (*fFcnMin));
for (UInt_t i=0; i<fParams.size(); i++) {
// only try to call minos if the parameter is not fixed!!
// the 1st condition is from an user fixed variable,
// the 2nd condition is from an all together unused variable
// the 3rd condition is a variable fixed via the FIX command
if ((fMnUserParams.Error(i) != 0.0) && (fRunInfo->ParameterInUse(i) != 0) && (!fMnUserParams.Parameters().at(i).IsFixed())) {
std::cout << ">> PFitter::ExecuteMinos(): calculate errors for " << fParams[i].fName << std::endl;
// 1-sigma MINOS errors
ROOT::Minuit2::MinosError err = minos.Minos(i);
if (err.IsValid()) {
// fill msr-file structure
fRunInfo->SetMsrParamStep(i, err.Lower());
fRunInfo->SetMsrParamPosError(i, err.Upper());
fRunInfo->SetMsrParamPosErrorPresent(i, true);
} else {
fRunInfo->SetMsrParamPosErrorPresent(i, false);
}
}
if (fMnUserParams.Parameters().at(i).IsFixed()) {
std::cerr << std::endl << ">> PFitter::ExecuteMinos(): **WARNING** Parameter " << fMnUserParams.Name(i) << " (ParamNo " << i+1 << ") is fixed!";
std::cerr << std::endl << ">> Will set STEP to zero, i.e. making it a constant parameter";
std::cerr << std::endl;
fRunInfo->SetMsrParamStep(i, 0.0);
fRunInfo->SetMsrParamPosErrorPresent(i, false);
}
}
end=MilliTime();
std::cout << ">> PFitter::ExecuteMinimize(): execution time for Minos = " << std::setprecision(3) << (end-start)/1.0e3 << " sec." << std::endl;
TString str = TString::Format("Minos: %.3f sec", (end-start)/1.0e3);
fElapsedTime.push_back(str);
return true;
}
//--------------------------------------------------------------------------
// ExecutePlot
//--------------------------------------------------------------------------
/**
* <p>Execute the minuit2 plot command.
*
* <b>return:</b> true.
*/
Bool_t PFitter::ExecutePlot()
{
std::cout << ">> PFitter::ExecutePlot() ..." << std::endl;
ROOT::Minuit2::MnPlot plot;
plot(fScanData);
return true;
}
//--------------------------------------------------------------------------
// ExecutePrintLevel
//--------------------------------------------------------------------------
/**
* <p>Set the print level.
*
* \param lineNo the line number of the command block
*
* <b>return:</b> true if done, otherwise returns false.
*/
Bool_t PFitter::ExecutePrintLevel(UInt_t lineNo)
{
std::cout << ">> PFitter::ExecutePrintLevel(): " << fCmdLines[lineNo].fLine.Data() << std::endl;
std::vector<std::string> tokens;
TString str;
tokens = PStringUtils::Split(fCmdLines[lineNo].fLine.Data(), ", \t");
if (static_cast<int>(tokens.size()) < 2) {
std::cerr << std::endl << "**ERROR** from PFitter::ExecutePrintLevel(): SYNTAX: PRINT_LEVEL <N>, where <N>=0-3" << std::endl << std::endl;
return false;
}
str = tokens[1].c_str();
Int_t ival;
if (str.IsDigit()) {
ival = str.Atoi();
if ((ival >=0) && (ival <= 3)) {
fPrintLevel = (UInt_t) ival;
} else {
std::cerr << std::endl << "**ERROR** from PFitter::ExecutePrintLevel(): SYNTAX: PRINT_LEVEL <N>, where <N>=0-3";
std::cerr << std::endl << " found <N>=" << ival << std::endl << std::endl;
return false;
}
} else {
std::cerr << std::endl << "**ERROR** from PFitter::ExecutePrintLevel(): SYNTAX: PRINT_LEVEL <N>, where <N>=0-3" << std::endl << std::endl;
return false;
}
#ifdef ROOT_GRTEQ_24
ROOT::Minuit2::MnPrint::SetGlobalLevel(fPrintLevel);
#else
ROOT::Minuit2::MnPrint::SetLevel(fPrintLevel);
#endif
return true;
}
//--------------------------------------------------------------------------
// ExecuteRelease
//--------------------------------------------------------------------------
/**
* <p>Release parameter list given at lineNo of the command block.
*
* \param lineNo the line number of the command block
*
* <b>return:</b> true if done, otherwise returns false.
*/
Bool_t PFitter::ExecuteRelease(UInt_t lineNo)
{
std::vector<std::string> tokens;
TString str;
tokens = PStringUtils::Split(fCmdLines[lineNo].fLine.Data(), ", \t");
std::cout << ">> PFitter::ExecuteRelease(): " << fCmdLines[lineNo].fLine.Data() << std::endl;
for (Int_t i=1; i<static_cast<int>(tokens.size()); i++) {
str = tokens[i].c_str();
if (str.IsDigit()) { // token is a parameter number
fMnUserParams.Release(static_cast<UInt_t>(str.Atoi())-1);
// set the error to 2% of the value when releasing
fMnUserParams.SetError(static_cast<UInt_t>(str.Atoi())-1, 0.02*fMnUserParams.Value(static_cast<UInt_t>(str.Atoi())-1));
} else { // token is a parameter name
fMnUserParams.Release(str.Data());
// set the error to 2% of the value when releasing
fMnUserParams.SetError(str.Data(), 0.02*fMnUserParams.Value(str.Data()));
}
}
return true;
}
//--------------------------------------------------------------------------
// ExecuteRestore
//--------------------------------------------------------------------------
/**
* <p>Release all fixed parameters
*
* <b>return:</b> true.
*/
Bool_t PFitter::ExecuteRestore()
{
std::cout << "PFitter::ExecuteRestore(): release all fixed parameters (RESTORE) ..." << std::endl;
for (UInt_t i=0; i<fMnUserParams.Parameters().size(); i++) {
if (fMnUserParams.Parameters().at(i).IsFixed()) {
fMnUserParams.Release(i);
fMnUserParams.SetError(i, 0.02*fMnUserParams.Value(i));
}
}
return true;
}
//--------------------------------------------------------------------------
// ExecuteScan
//--------------------------------------------------------------------------
/**
* <p>Execute the minuit2 scan command.
*
* <b>return:</b> true.
*/
Bool_t PFitter::ExecuteScan()
{
std::cout << ">> PFitter::ExecuteScan(): will call scan ..." << std::endl;
ROOT::Minuit2::MnScan scan((*fFitterFcn), fMnUserParams);
if (fScanAll) { // not clear at the moment what to be done here
// TO BE IMPLEMENTED
} else { // single parameter scan
fScanData = scan.Scan(fScanParameter[0], fScanNoPoints, fScanLow, fScanHigh);
}
fConverged = true;
return true;
}
//--------------------------------------------------------------------------
// ExecuteSave
//--------------------------------------------------------------------------
/**
* <p>Execute the save command.
*
* \param firstSave flag indication if this is the first save call and hence write a fresh MINUIT2.OUTPUT
*
* <b>return:</b> true if the valid minuit2 state is found, otherwise returns false.
*/
Bool_t PFitter::ExecuteSave(Bool_t firstSave)
{
// if any minimization was done, otherwise get out immediately
if (!fFcnMin) {
std::cout << std::endl << ">> PFitter::ExecuteSave(): nothing to be saved ...";
return false;
}
ROOT::Minuit2::MnUserParameterState mnState = fFcnMin->UserState();
// check if the user parameter state is valid
if (!mnState.IsValid()) {
std::cerr << std::endl << ">> PFitter::ExecuteSave: **WARNING** Minuit2 User Parameter State is not valid, i.e. nothing to be saved";
std::cerr << std::endl;
return false;
}
// handle expected chisq if applicable
fParams = *(fRunInfo->GetMsrParamList()); // get the update parameters back
// calculate expected chisq
std::vector<Double_t> param;
std::vector<Double_t> err;
Double_t totalExpectedChisq = 0.0;
std::vector<Double_t> expectedchisqPerRun;
std::vector<UInt_t> ndfPerHisto;
for (UInt_t i=0; i<fParams.size(); i++) {
param.push_back(fParams[i].fValue);
err.push_back(fParams[i].fStep);
}
// CalcExpectedChiSquare handles both, chisq and mlh
Bool_t ok;
PDoubleVector par_r = ParamRound(param, err, ok);
if (!ok)
par_r = param;
fFitterFcn->CalcExpectedChiSquare(par_r, totalExpectedChisq, expectedchisqPerRun);
// calculate chisq per run
std::vector<Double_t> chisqPerRun;
for (UInt_t i=0; i<fRunInfo->GetMsrRunList()->size(); i++) {
if (fUseChi2)
chisqPerRun.push_back(fRunListCollection->GetSingleRunChisq(par_r, i));
else
chisqPerRun.push_back(fRunListCollection->GetSingleRunMaximumLikelihood(par_r, i));
}
if (totalExpectedChisq != 0.0) { // i.e. applicable for single histogram fits only
// get the ndf's of the histos
UInt_t ndf_run;
for (UInt_t i=0; i<expectedchisqPerRun.size(); i++) {
ndf_run = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i);
ndfPerHisto.push_back(ndf_run);
}
// feed the msr-file handler
PMsrStatisticStructure *statistics = fRunInfo->GetMsrStatistic();
if (statistics) {
statistics->fMinPerHisto = chisqPerRun;
statistics->fMinExpected = totalExpectedChisq;
statistics->fMinExpectedPerHisto = expectedchisqPerRun;
statistics->fNdfPerHisto = ndfPerHisto;
}
} 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; i<chisqPerRun.size(); i++) {
ndf_run = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i);
ndfPerHisto.push_back(ndf_run);
}
// feed the msr-file handler
PMsrStatisticStructure *statistics = fRunInfo->GetMsrStatistic();
if (statistics) {
statistics->fMinPerHisto = chisqPerRun;
statistics->fNdfPerHisto = ndfPerHisto;
}
}
// check if sector command has been requested
if (fSectorFlag) {
PDoubleVector error;
for (UInt_t i=0; i<fParams.size(); i++)
error.push_back(fParams[i].fStep);
PrepareSector(param, error);
}
// clean up
param.clear();
expectedchisqPerRun.clear();
ndfPerHisto.clear();
chisqPerRun.clear();
std::cout << ">> PFitter::ExecuteSave(): will write minuit2 output file ..." << std::endl;
std::ofstream fout;
// open minuit2 output file
if (firstSave)
fout.open("MINUIT2.OUTPUT", std::iostream::out);
else
fout.open("MINUIT2.OUTPUT", std::iostream::out | std::iostream::app);
if (!fout.is_open()) {
std::cerr << std::endl << "**ERROR** PFitter::ExecuteSave() couldn't open MINUIT2.OUTPUT file";
std::cerr << std::endl;
return false;
}
// write header
TDatime dt;
fout << std::endl << "*************************************************************************";
fout << std::endl << " musrfit MINUIT2 output file from " << fRunInfo->GetFileName().Data() << " - " << dt.AsSQLString();
fout << std::endl << "*************************************************************************";
fout << std::endl;
// write elapsed times
fout << std::endl << " elapsed times:";
for (UInt_t i=0; i<fElapsedTime.size(); i++) {
fout << std::endl << " " << fElapsedTime[i];
}
fout << std::endl;
fout << std::endl << "*************************************************************************";
fout << std::endl;
fElapsedTime.clear();
// write global information
fout << std::endl << " Fval() = " << mnState.Fval() << ", Edm() = " << mnState.Edm() << ", NFcn() = " << mnState.NFcn();
fout << std::endl;
fout << std::endl << "*************************************************************************";
fout << std::endl;
// identifiy the longest variable name for proper formating reasons
Int_t maxLength = 10;
for (UInt_t i=0; i<fParams.size(); i++) {
if (fParams[i].fName.Length() > maxLength)
maxLength = fParams[i].fName.Length() + 1;
}
// write parameters
fParams = *(fRunInfo->GetMsrParamList()); // get the update parameters back
fout << std::endl << " PARAMETERS";
fout << std::endl << "-------------------------------------------------------------------------";
fout << std::endl << " ";
for (Int_t j=0; j<=maxLength-4; j++)
fout << " ";
fout << "Parabolic Minos";
fout << std::endl << " No Name";
for (Int_t j=0; j<=maxLength-4; j++)
fout << " ";
fout << "Value Error Negative Positive Limits";
for (UInt_t i=0; i<fParams.size(); i++) {
// write no
fout.setf(std::ios::right, std::ios::adjustfield);
fout.width(3);
fout << std::endl << i+1 << " ";
// write name
fout << fParams[i].fName.Data();
for (Int_t j=0; j<=maxLength-fParams[i].fName.Length(); j++)
fout << " ";
// write value
fout.setf(std::ios::left, std::ios::adjustfield);
fout.precision(6);
fout.width(10);
fout << fParams[i].fValue << " ";
// write parabolic error
fout.setf(std::ios::left, std::ios::adjustfield);
fout.precision(6);
fout.width(10);
if (fParams[i].fStep != 0.0)
fout << fMnUserParams.Error(i) << " ";
else
fout << "---";
// write minos errors
if (fParams[i].fPosErrorPresent) {
fout.setf(std::ios::left, std::ios::adjustfield);
fout.precision(6);
fout.width(12);
fout << fParams[i].fStep;
fout.setf(std::ios::left, std::ios::adjustfield);
fout.precision(6);
fout.width(11);
fout << fParams[i].fPosError << " ";
} else {
fout.setf(std::ios::left, std::ios::adjustfield);
fout.width(12);
fout << "---";
fout.setf(std::ios::left, std::ios::adjustfield);
fout.width(12);
fout << "---";
}
// write limits
if (fParams[i].fNoOfParams > 5) {
fout.setf(std::ios::left, std::ios::adjustfield);
fout.width(7);
if (fParams[i].fLowerBoundaryPresent)
fout << fParams[i].fLowerBoundary;
else
fout << "---";
fout.setf(std::ios::left, std::ios::adjustfield);
fout.width(7);
if (fParams[i].fUpperBoundaryPresent)
fout << fParams[i].fUpperBoundary;
else
fout << "---";
} else {
fout.setf(std::ios::left, std::ios::adjustfield);
fout.width(7);
fout << "---";
fout.setf(std::ios::left, std::ios::adjustfield);
fout.width(7);
fout << "---";
}
}
fout << std::endl;
fout << std::endl << "*************************************************************************";
// write covariance matrix
fout << std::endl << " COVARIANCE MATRIX";
fout << std::endl << "-------------------------------------------------------------------------";
if (mnState.HasCovariance()) {
ROOT::Minuit2::MnUserCovariance cov = mnState.Covariance();
fout << std::endl << "from " << cov.Nrow() << " free parameters";
for (UInt_t i=0; i<cov.Nrow(); i++) {
fout << std::endl;
for (UInt_t j=0; j<i; j++) {
fout.setf(std::ios::left, std::ios::adjustfield);
fout.precision(6);
if (cov(i,j) > 0.0) {
fout << " ";
fout.width(13);
} else {
fout.width(14);
}
fout << cov(i,j);
}
}
} else {
fout << std::endl << " no covariance matrix available";
}
fout << std::endl;
fout << std::endl << "*************************************************************************";
// write correlation matrix
fout << std::endl << " CORRELATION COEFFICIENTS";
fout << std::endl << "-------------------------------------------------------------------------";
if (mnState.IsValid() && mnState.HasCovariance()) {
ROOT::Minuit2::MnGlobalCorrelationCoeff corr = mnState.GlobalCC();
ROOT::Minuit2::MnUserCovariance cov = mnState.Covariance();
PIntVector parNo;
fout << std::endl << " No Global ";
for (UInt_t i=0; i<fParams.size(); i++) {
// only free parameters, i.e. not fixed, and not unsed ones!
if ((fParams[i].fStep != 0) && (fRunInfo->ParameterInUse(i) > 0) && (!fMnUserParams.Parameters().at(i).IsFixed())) {
fout.setf(std::ios::left, std::ios::adjustfield);
fout.width(9);
fout << i+1;
parNo.push_back(i);
}
}
// check that there is a correspondens between minuit2 and musrfit book keeping
if (parNo.size() != cov.Nrow()) {
std::cerr << std::endl << "**SEVERE ERROR** in PFitter::ExecuteSave(): minuit2 and musrfit book keeping to not correspond! Unable to write correlation matrix.";
std::cerr << std::endl;
} else { // book keeping is OK
TString title("Minuit2 Output Correlation Matrix for ");
title += fRunInfo->GetFileName();
title += " - ";
title += dt.AsSQLString();
std::unique_ptr<TCanvas> ccorr = std::make_unique<TCanvas>("ccorr", "title", 500, 500);
std::unique_ptr<TH2D> hcorr = std::make_unique<TH2D>("hcorr", title, cov.Nrow(), 0.0, cov.Nrow(), cov.Nrow(), 0.0, cov.Nrow());
Double_t dval;
for (UInt_t i=0; i<cov.Nrow(); i++) {
// parameter number
fout << std::endl << " ";
fout.setf(std::ios::left, std::ios::adjustfield);
fout.width(5);
fout << parNo[i]+1;
// global correlation coefficient
fout.setf(std::ios::left, std::ios::adjustfield);
fout.precision(6);
fout.width(12);
fout << corr.GlobalCC()[i];
// correlations matrix
for (UInt_t j=0; j<cov.Nrow(); j++) {
fout.setf(std::ios::left, std::ios::adjustfield);
// fout.precision(4);
if (i==j) {
fout.width(9);
fout << " 1.0 ";
hcorr->Fill((Double_t)i,(Double_t)i,1.0);
} else {
// check that errors are none zero
if (fMnUserParams.Error(parNo[i]) == 0.0) {
std::cerr << std::endl << "**SEVERE ERROR** in PFitter::ExecuteSave(): parameter no " << parNo[i]+1 << " has an error == 0. Cannot correctly handle the correlation matrix.";
std::cerr << std::endl;
dval = 0.0;
} else if (fMnUserParams.Error(parNo[j]) == 0.0) {
std::cerr << std::endl << "**SEVERE ERROR** in PFitter::ExecuteSave(): parameter no " << parNo[j]+1 << " has an error == 0. Cannot correctly handle the correlation matrix.";
std::cerr << std::endl;
dval = 0.0;
} else {
dval = cov(i,j)/(fMnUserParams.Error(parNo[i])*fMnUserParams.Error(parNo[j]));
}
hcorr->Fill((Double_t)i,(Double_t)j,dval);
// handle precision, ugly but ...
if (dval < 1.0e-2) {
fout.precision(2);
} else {
fout.precision(4);
}
// handle sign
if (dval > 0.0) {
fout << " ";
fout.width(7);
} else {
fout.width(8);
}
fout << dval << " ";
}
}
}
// write correlation matrix into a root file
TFile ff("MINUIT2.root", "recreate");
ccorr->Draw();
if (cov.Nrow() <= 6)
hcorr->Draw("COLZTEXT");
else
hcorr->Draw("COLZ");
ccorr->Write("ccorr", TObject::kOverwrite, sizeof(ccorr));
hcorr->Write("hcorr", TObject::kOverwrite, sizeof(hcorr));
ff.Close();
if (fYamlOut) {
// write the fit results to an easy-to-read/parse yaml file
// note: the block names follow those used by Python library iminuit
// https://github.com/scikit-hep/iminuit
// dynamically name the yaml output file
// https://stackoverflow.com/a/25389052
std::string yaml_filename(fRunInfo->GetFileName().Data());
const std::string msr_ext(".msr");
yaml_filename.replace(yaml_filename.find(msr_ext), msr_ext.length(),
".yaml");
// define yaml's 2-space indentation
const std::string yaml_indent(" ");
// open the yaml file for writing
std::ofstream yaml_file(yaml_filename);
// number formatting of the output
yaml_file << std::scientific << std::setprecision(16);
// write the parameter values
yaml_file << "values:\n";
for (unsigned int i = 0; i < fParams.size(); ++i) {
yaml_file << yaml_indent << fParams[i].fName.Data() << ": "
<< fParams[i].fValue << "\n";
}
// write the parabolic errors
yaml_file << "errors:\n";
for (unsigned int i = 0; i < fParams.size(); ++i) {
yaml_file << yaml_indent << fParams[i].fName.Data() << ": "
<< fMnUserParams.Error(i) << "\n";
}
// write the minos errors
yaml_file << "mnerrors:\n";
for (unsigned int i = 0; i < fParams.size(); ++i) {
// use boost's implementation of a variant, which can be streamed by
// default - see: https://stackoverflow.com/q/47168477
boost::variant<double, std::string> positive_error, negative_error;
if (fParams[i].fPosErrorPresent) {
positive_error = fParams[i].fPosError;
negative_error = fParams[i].fStep;
} else {
positive_error = "null";
negative_error = "null";
}
yaml_file << yaml_indent << fParams[i].fName.Data() << ":\n";
yaml_file << yaml_indent << yaml_indent
<< "positive: " << positive_error << "\n";
yaml_file << yaml_indent << yaml_indent
<< "negative: " << negative_error << "\n";
}
// write the parameter limits
yaml_file << "limits:\n";
for (unsigned int i = 0; i < fParams.size(); ++i) {
// use boost's implementation of a variant, which can be streamed by
// default - see: https://stackoverflow.com/q/47168477
boost::variant<double, std::string> upper_limit, lower_limit;
if (fParams[i].fLowerBoundaryPresent) {
lower_limit = fParams[i].fLowerBoundary;
} else {
lower_limit = "null";
}
if (fParams[i].fUpperBoundaryPresent) {
upper_limit = fParams[i].fUpperBoundary;
} else {
upper_limit = "null";
}
yaml_file << yaml_indent << fParams[i].fName.Data() << ":\n";
yaml_file << yaml_indent << yaml_indent << "lower: " << lower_limit
<< "\n";
yaml_file << yaml_indent << yaml_indent << "upper: " << upper_limit
<< "\n";
}
// write if the parameter is fixed
yaml_file << "fixed:\n";
for (unsigned int i = 0; i < fParams.size(); ++i) {
std::string is_fixed = fParams[i].fStep == 0.0 ? "true" : "false";
yaml_file << yaml_indent << fParams[i].fName.Data() << ": " << is_fixed
<< "\n";
}
// write the covariance matrix (omitting fixed parameters)
yaml_file << "covariance:\n";
for (unsigned int i = 0; i < cov.Nrow(); ++i) {
yaml_file << yaml_indent << mnState.Name(parNo[i]) << ":\n";
for (unsigned int j = 0; j < cov.Nrow(); ++j) {
yaml_file << yaml_indent << yaml_indent << mnState.Name(parNo[j])
<< ": " << cov(i, j) << "\n";
}
}
// write the correlation matrix (omitting fixed parameters)
yaml_file << "correlation:\n";
for (unsigned int i = 0; i < cov.Nrow(); ++i) {
yaml_file << yaml_indent << mnState.Name(parNo[i]) << ":\n";
for (unsigned int j = 0; j < cov.Nrow(); ++j) {
double correlation =
i == j ? 1.0
: cov(i, j) / (fMnUserParams.Error(parNo[i]) *
fMnUserParams.Error(parNo[j]));
yaml_file << yaml_indent << yaml_indent << mnState.Name(parNo[j])
<< ": " << correlation << "\n";
}
}
// close the yaml file
yaml_file.close();
}
}
parNo.clear(); // clean up
} 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::fixed, 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 ";
fout << std::endl << "*************************************************************************";
fout << std::endl << std::endl;
// close MINUIT2.OUTPUT file
fout.close();
return true;
}
//--------------------------------------------------------------------------
// ExecuteSimplex
//--------------------------------------------------------------------------
/**
* <p>Execute the minuit2 simplex command.
*
* <b>return:</b> true if the simplex command could be executed successfully, otherwise returns false.
*/
Bool_t PFitter::ExecuteSimplex()
{
std::cout << ">> PFitter::ExecuteSimplex(): will call simplex ..." << std::endl;
// create minimizer object
// strategy is by default = 'default'
ROOT::Minuit2::MnSimplex simplex((*fFitterFcn), fMnUserParams, ROOT::Minuit2::MnStrategy{fStrategy});
// minimize
// maxfcn is 10*MINUIT2 Default maxfcn
UInt_t maxfcn = std::numeric_limits<UInt_t>::max();
// tolerance = MINUIT2 Default tolerance
Double_t tolerance = 0.1;
// keep track of elapsed time
Double_t start=0.0, end=0.0;
start=MilliTime();
ROOT::Minuit2::FunctionMinimum min = simplex(maxfcn, tolerance);
end=MilliTime();
std::cout << ">> PFitter::ExecuteMinimize(): execution time for Simplex = " << std::setprecision(3) << (end-start)/1.0e3 << " sec." << std::endl;
TString str = TString::Format("Simplex: %.3f sec", (end-start)/1.0e3);
fElapsedTime.push_back(str);
if (!min.IsValid()) {
std::cerr << std::endl << ">> PFitter::ExecuteSimplex(): **WARNING**: Fit did not converge, sorry ...";
std::cerr << std::endl;
fIsValid = false;
return false;
}
// keep FunctionMinimum object
fFcnMin.reset(new ROOT::Minuit2::FunctionMinimum(min));
// keep user parameters
if (fFcnMin)
fMnUserParams = fFcnMin->UserParameters();
// fill run info
for (UInt_t i=0; i<fParams.size(); i++) {
Double_t dval = min.UserState().Value(i);
if (fPhase[i]) {
Int_t m = (Int_t)(dval/360.0);
dval = dval - m*360.0;
}
fRunInfo->SetMsrParamValue(i, dval);
fRunInfo->SetMsrParamStep(i, min.UserState().Error(i));
fRunInfo->SetMsrParamPosErrorPresent(i, false);
}
// handle statistics
Double_t minVal = min.Fval();
UInt_t ndf = fFitterFcn->GetTotalNoOfFittedBins();
// subtract number of varied parameters from total no of fitted bins -> ndf
for (UInt_t i=0; i<fParams.size(); i++) {
if ((min.UserState().Error(i) != 0.0) && !fMnUserParams.Parameters().at(i).IsFixed())
ndf -= 1;
}
// feed run info with new statistics info
fRunInfo->SetMsrStatisticMin(minVal);
fRunInfo->SetMsrStatisticNdf(ndf);
fConverged = true;
if (fPrintLevel >= 2)
std::cout << *fFcnMin << std::endl;
return true;
}
//--------------------------------------------------------------------------
// PrepareSector
//--------------------------------------------------------------------------
/**
* <p>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(PDoubleVector &param, PDoubleVector &error)
{
Double_t val;
UInt_t ndf;
Int_t usedParams = 0;
for (UInt_t i=0; i<error.size(); i++) {
if (error[i] != 0.0)
usedParams++;
}
PDoublePairVector secFitRange;
secFitRange.resize(1);
if (fUseChi2) {
Double_t totalExpectedChisq = 0.0;
PDoubleVector expectedChisqPerRun;
PDoubleVector chisqPerRun;
for (UInt_t k=0; k<fSector.size(); k++) {
// set sector fit range
secFitRange[0].first = fSector[k].GetTimeRangeFirst(0);
secFitRange[0].second = fSector[k].GetTimeRangeLast();
fRunListCollection->SetFitRange(secFitRange);
// calculate chisq
val = (*fFitterFcn)(param);
fSector[k].SetChisq(val);
// calculate NDF
ndf = static_cast<UInt_t>(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; i<fRunInfo->GetMsrRunList()->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; i<expectedChisqPerRun.size(); i++) {
ndf_run = fFitterFcn->GetNoOfFittedBins(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; i<chisqPerRun.size(); i++) {
ndf_run = fFitterFcn->GetNoOfFittedBins(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; k<fSector.size(); k++) {
// set sector fit range
secFitRange[0].first = fSector[k].GetTimeRangeFirst(0);
secFitRange[0].second = fSector[k].GetTimeRangeLast();
fRunListCollection->SetFitRange(secFitRange);
// calculate maxLH
val = (*fFitterFcn)(param);
fSector[k].SetChisq(val);
// calculate NDF
ndf = static_cast<UInt_t>(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; i<fRunInfo->GetMsrRunList()->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; i<expectedMaxLHPerRun.size(); i++) {
ndf_run = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i);
if (ndf_run > 0) {
fSector[k].SetNDF(ndf_run, i);
}
}
}
// clean up
maxLHPerRun.clear();
expectedMaxLHPerRun.clear();
}
}
}
//--------------------------------------------------------------------------
// ExecuteSector
//--------------------------------------------------------------------------
/**
* <p>Write all chisq/maxLH sector information to MINUIT.OUTPUT and dump it
* to stdout.
*
* <b>return:</b> if the sector command was successful, otherwise return flase.
*/
Bool_t PFitter::ExecuteSector(std::ofstream &fout)
{
fout << std::endl;
fout << std::endl << "*************************************************************************";
fout << std::endl << " SECTOR RESULTS";
fout << std::endl << "*************************************************************************";
if (fUseChi2) {
for (UInt_t i=0; i<fSector.size(); i++) {
fout << std::endl;
fout << " Sector " << i+1 << ": FitRange: " << fSector[i].GetTimeRangeFirst(0) << ", " << fSector[i].GetTimeRangeLast() << std::endl;
fout << " chisq = " << fSector[i].GetChisq() << ", NDF = " << fSector[i].GetNDF() << ", chisq/NDF = " << fSector[i].GetChisq()/fSector[i].GetNDF();
if (fSector[i].GetExpectedChisq() > 0) {
fout << std::endl << " expected chisq = " << fSector[i].GetExpectedChisq() << ", NDF = " << fSector[i].GetNDF() << ", expected chisq/NDF = " << fSector[i].GetExpectedChisq()/fSector[i].GetNDF();
for (UInt_t j=0; j<fSector[i].GetNoRuns(); j++) {
fout << std::endl << " run block " << j+1 << " (NDF/red.chisq/red.chisq_e) = (" << fSector[i].GetNDF(j) << "/" << fSector[i].GetChisq(j)/fSector[i].GetNDF(j) << "/" << fSector[i].GetExpectedChisq(j)/fSector[i].GetNDF(j) << ")";
}
} else if (fSector[i].GetNoRuns() > 0) { // in case expected chisq is not applicable like for asymmetry fits
for (UInt_t j=0; j<fSector[i].GetNoRuns(); j++) {
fout << std::endl << " run block " << j+1 << " (NDF/red.chisq) = (" << fSector[i].GetNDF(j) << "/" << fSector[i].GetChisq(j)/fSector[i].GetNDF(j) << ")";
}
}
fout << std::endl << "++++";
}
} else { // max log likelihood
for (UInt_t i=0; i<fSector.size(); i++) {
fout << std::endl;
fout << " Sector " << i+1 << ": FitRange: " << fSector[i].GetTimeRangeFirst(0) << ", " << fSector[i].GetTimeRangeLast() << std::endl;
fout << " maxLH = " << fSector[i].GetChisq() << ", NDF = " << fSector[i].GetNDF() << ", maxLH/NDF = " << fSector[i].GetChisq()/fSector[i].GetNDF();
if (fSector[i].GetExpectedChisq() > 0) {
fout << std::endl << " expected maxLH = " << fSector[i].GetExpectedChisq() << ", NDF = " << fSector[i].GetNDF() << ", expected maxLH/NDF = " << fSector[i].GetExpectedChisq()/fSector[i].GetNDF();
for (UInt_t j=0; j<fSector[i].GetNoRuns(); j++) {
fout << std::endl << " run block " << j+1 << " (NDF/red.maxLH/red.maxLH_e) = (" << fSector[i].GetNDF(j) << "/" << fSector[i].GetChisq(j)/fSector[i].GetNDF(j) << "/" << fSector[i].GetExpectedChisq(j)/fSector[i].GetNDF(j) << ")";
}
} else if (fSector[i].GetNoRuns() > 0) { // in case expected chisq is not applicable like for asymmetry fits
for (UInt_t j=0; j<fSector[i].GetNoRuns(); j++) {
fout << std::endl << " run block " << j+1 << " (NDF/red.maxLH) = (" << fSector[i].GetNDF(j) << "/" << fSector[i].GetChisq(j)/fSector[i].GetNDF(j) << ")";
}
}
fout << std::endl << "++++";
}
}
return true;
}
//--------------------------------------------------------------------------
// MilliTime
//--------------------------------------------------------------------------
/**
* <p>
*
* <b>return:</b>
*/
Double_t PFitter::MilliTime()
{
struct timeval now;
gettimeofday(&now, nullptr);
return ((Double_t)now.tv_sec * 1.0e6 + (Double_t)now.tv_usec)/1.0e3;
}
//--------------------------------------------------------------------------
// ParamRound (private)
//--------------------------------------------------------------------------
/**
* <p>Rounds the parameter vector value according to the given error estimate,
* so that the msr-file value and the fitter result are consistent with each
* other. This means that musrfit -c, and musrfit -e -t should give essentially
* the same values of expected chisq (up to small rounding values).
*
* @param par parameter value vector
* @param err error value vector
* @param ok true if size of par and err are identically, otherwise false.
*
* @return rounded parameter value vector, compatible with the msr-file output.
*/
PDoubleVector PFitter::ParamRound(const PDoubleVector &par, const PDoubleVector &err, Bool_t &ok)
{
PDoubleVector par_r;
par_r.resize(par.size());
ok = true;
if (par.size() != err.size()) {
// error msg
ok = false;
return par_r;
}
int exp;
double dval;
for (unsigned int i=0; i<par.size(); i++) {
if (err[i] != 0.0) {
exp = floor(log10(fabs(err[i]))-1);
dval = round(par[i]*pow(10.0, -exp))/pow(10.0, -exp);
par_r[i] = dval;
} else {
par_r[i] = par[i];
}
}
return par_r;
}
//-------------------------------------------------------------------------------------------------
// end
//-------------------------------------------------------------------------------------------------
Reference in New Issue View Git Blame Copy Permalink