diff --git a/src/ToDo.txt b/src/ToDo.txt index 24446324..32d7f386 100644 --- a/src/ToDo.txt +++ b/src/ToDo.txt @@ -89,6 +89,7 @@ short term: * PFitter.cpp: the value of the parabolic error in the MINUIT2.OUTPUT is wrong (only MIGRAD is called and neither HESSE nor MINOS). Suspect the problem in the ExecuteSave() routine. Needs to be checked. + **DONE** 08-09-01 * PFitter.cpp: implement HESSE diff --git a/src/classes/PFitter.cpp b/src/classes/PFitter.cpp index 7a2f7a15..a931ed5a 100644 --- a/src/classes/PFitter.cpp +++ b/src/classes/PFitter.cpp @@ -134,8 +134,7 @@ bool PFitter::DoFit() cout << endl; break; case PMN_HESSE: - cout << endl << "**WARNING** from PFitter::DoFit() : the command HESSE is not yet implemented."; - cout << endl; + status = ExecuteHesse(); break; case PMN_MACHINE_PRECISION: cout << endl << "**WARNING** from PFitter::DoFit() : the command MACHINE_PRECISION is not yet implemented."; @@ -270,18 +269,17 @@ bool PFitter::CheckCommands() */ bool PFitter::SetParameters() { - PMsrParamList::iterator it; - - for (it = fParams.begin(); it != fParams.end(); ++it) { + for (unsigned int i=0; ifStep == 0.0) { // add fixed parameter - fMnUserParams.Add(it->fName.Data(), it->fValue); + 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 (it->fNoOfParams > 5) { // boundaries given - fMnUserParams.Add(it->fName.Data(), it->fValue, it->fStep, it->fLowerBoundary, it->fUpperBoundary); + if (fParams[i].fNoOfParams > 5) { // boundaries given + fMnUserParams.Add(fParams[i].fName.Data(), fParams[i].fValue, fParams[i].fStep, + fParams[i].fLowerBoundary, fParams[i].fUpperBoundary); } else { // no boundaries given - fMnUserParams.Add(it->fName.Data(), it->fValue, it->fStep); + fMnUserParams.Add(fParams[i].fName.Data(), fParams[i].fValue, fParams[i].fStep); } } } @@ -298,6 +296,24 @@ bool PFitter::SetParameters() return true; } +//-------------------------------------------------------------------------- +// ExecuteHesse +//-------------------------------------------------------------------------- +/** + *

+ * + */ +bool PFitter::ExecuteHesse() +{ + cout << "PFitter::ExecuteHesse(): will call hesse ..." << endl; + + // if already some minimization is done use the minuit2 output as input + if (fFcnMin) + fMnUserParams = fFcnMin->UserParameters(); + + return true; +} + //-------------------------------------------------------------------------- // ExecuteMigrad //-------------------------------------------------------------------------- @@ -534,7 +550,7 @@ bool PFitter::ExecuteSave() fout.setf(ios::left, ios::adjustfield); fout.precision(6); fout.width(11); - fout << fMnUserParams.Error(i); + fout << mnState.Error(i); // write minos errors if (fParams[i].fPosErrorPresent) { fout.setf(ios::left, ios::adjustfield); @@ -554,7 +570,7 @@ bool PFitter::ExecuteSave() fout << "---"; } // write limits - if (!isnan(fParams[i].fLowerBoundary)) { + if (fParams[i].fNoOfParams > 5) { fout.setf(ios::left, ios::adjustfield); fout.width(7); fout << fParams[i].fLowerBoundary; diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index cca350d9..67d5137c 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -506,15 +506,16 @@ int PMsrHandler::WriteMsrLogFile() CheckAndWriteComment(f, ++lineNo); } // backgr.fix - if (!isnan(fRuns[i].fBkgFix[0])) { + if (fRuns[i].fBkgFixPresent[0]) { f.width(15); f << endl << left << "backgr.fix"; - for (unsigned int j=0; j<2; j++) { - if (!isnan(fRuns[i].fBkgFix[j])) { - f.precision(prec); - f.width(12); - f << left << fRuns[i].fBkgFix[j]; - } + f.precision(prec); + f.width(12); + f << left << fRuns[i].fBkgFix[0]; + if (fRuns[i].fBkgFixPresent[1]) { + f.precision(prec); + f.width(12); + f << left << fRuns[i].fBkgFix[1]; } CheckAndWriteComment(f, ++lineNo); } @@ -830,21 +831,21 @@ bool PMsrHandler::HandleFitParameterEntry(PMsrLines &lines) TObjString *ostr; TString str; - // init param structure - param.fNoOfParams = -1; - param.fNo = -1; - param.fName = TString(""); - param.fValue = nan("NAN"); - param.fStep = nan("NAN"); - param.fPosErrorPresent = nan("NAN"); - param.fPosError = nan("NAN"); - param.fLowerBoundary = nan("NAN"); - param.fUpperBoundary = nan("NAN"); - // fill param structure iter = lines.begin(); while ((iter != lines.end()) && !error) { + // init param structure + param.fNoOfParams = -1; + param.fNo = -1; + param.fName = TString(""); + param.fValue = 0.0; + param.fStep = 0.0; + param.fPosErrorPresent = false; + param.fPosError = 0.0; + param.fLowerBoundary = 0.0; + param.fUpperBoundary = 0.0; + tokens = iter->fLine.Tokenize(" \t"); if (!tokens) { cout << endl << "SEVERE ERROR: Couldn't tokenize Parameters in line " << iter->fLineNo; @@ -991,6 +992,7 @@ bool PMsrHandler::HandleFitParameterEntry(PMsrLines &lines) cout << endl; } else { // everything is OK, therefore add the parameter to the parameter list fParam.push_back(param); +cout << endl << ">> PMsrHandler::HandleFitParameterEntry: i=" << fParam.size() << ", param.fLowerBoundary=" << param.fLowerBoundary << ", param.fUpperBoundary=" << param.fUpperBoundary; } // clean up @@ -1509,7 +1511,7 @@ bool PMsrHandler::HandleRunEntry(PMsrLines &lines) bool found; if (fRuns[i].fBkgFitParamNo >= 0) { // check if backgr.fit is given found = true; - } else if (!isnan(fRuns[i].fBkgFix[0])) { // check if backgr.fix is given + } else if (fRuns[i].fBkgFixPresent[0]) { // check if backgr.fix is given found = true; } else if (fRuns[i].fBkgRange[0] >= 0) { // check if background window is given found = true; @@ -1555,8 +1557,10 @@ void PMsrHandler::InitRunParameterStructure(PMsrRunStructure ¶m) param.fMap.clear(); // empty list param.fForwardHistoNo = -1; param.fBackwardHistoNo = -1; - for (int i=0; i<2; i++) - param.fBkgFix[i] = nan("NAN"); + for (int i=0; i<2; i++) { + param.fBkgFixPresent[i] = false; + param.fBkgFix[i] = 0.0; + } for (int i=0; i<4; i++) param.fBkgRange[i] = -1; for (int i=0; i<4; i++) diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index 220a6e99..4b839c5c 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -594,7 +594,7 @@ void PMusrCanvas::UpdateInfoPad() // temperature if present tstr += TString("T="); dval = fRunList->GetTemp(runs[runNo].fRunName); - if (isnan(dval)) { + if (dval == -9.9e99) { tstr += TString("??,"); } else { sprintf(sval, "%0.2lf", dval); @@ -603,7 +603,7 @@ void PMusrCanvas::UpdateInfoPad() // field if present tstr += TString("B="); dval = fRunList->GetField(runs[runNo].fRunName); - if (isnan(dval)) { + if (dval == -9.9e99) { tstr += TString("??,"); } else { sprintf(sval, "%0.2lf", dval); @@ -612,7 +612,7 @@ void PMusrCanvas::UpdateInfoPad() // energy if present tstr += TString("E="); dval = fRunList->GetEnergy(runs[runNo].fRunName); - if (isnan(dval)) { + if (dval == -9.9e99) { tstr += TString("??,"); } else { sprintf(sval, "%0.2lf", dval); diff --git a/src/classes/PRunAsymmetry.cpp b/src/classes/PRunAsymmetry.cpp index 28715c44..1d50f7be 100644 --- a/src/classes/PRunAsymmetry.cpp +++ b/src/classes/PRunAsymmetry.cpp @@ -356,7 +356,7 @@ bool PRunAsymmetry::PrepareData() } // subtract background from histogramms ------------------------------------------ - if (isnan(fRunInfo->fBkgFix[0])) { // no fixed background given + if (!fRunInfo->fBkgFixPresent[0]) { // no fixed background given if (fRunInfo->fBkgRange[0] != 0) { if (!SubtractEstimatedBkg()) return false; diff --git a/src/classes/PRunDataHandler.cpp b/src/classes/PRunDataHandler.cpp index d59b78bb..13319df6 100644 --- a/src/classes/PRunDataHandler.cpp +++ b/src/classes/PRunDataHandler.cpp @@ -488,10 +488,10 @@ bool PRunDataHandler::ReadNemuFile() runData.fRunName = TString(""); runData.fRunTitle = TString(""); runData.fSetup = TString(""); - runData.fField = nan("NAN"); - runData.fTemp = nan("NAN"); - runData.fEnergy = nan("NAN"); - runData.fTimeResolution = nan("NAN"); + runData.fField = -9.9e99; + runData.fTemp = -9.9e99; + runData.fEnergy = -9.9e99; + runData.fTimeResolution = 0.0; // open file ifstream f; @@ -578,7 +578,7 @@ bool PRunDataHandler::ReadNemuFile() f.getline(instr, sizeof(instr)); } while (headerInfo && !f.eof()); - if ((groups == 0) || (channels == 0) || isnan(runData.fTimeResolution)) { + if ((groups == 0) || (channels == 0) || runData.fTimeResolution == 0.0) { cout << endl << "PRunDataHandler::ReadNemuFile(): essential header informations are missing!"; f.close(); return false; diff --git a/src/classes/PRunSingleHisto.cpp b/src/classes/PRunSingleHisto.cpp index 5c3ff886..14ebb96f 100644 --- a/src/classes/PRunSingleHisto.cpp +++ b/src/classes/PRunSingleHisto.cpp @@ -119,7 +119,7 @@ double PRunSingleHisto::CalcChiSquare(const std::vector& par) // get background double bkg; if (fRunInfo->fBkgFitParamNo == -1) { // bkg not fitted - if (isnan(fRunInfo->fBkgFix[0])) { // no fixed background given (background interval) + if (!fRunInfo->fBkgFixPresent[0]) { // no fixed background given (background interval) bkg = fBackground; } else { // fixed bkg given bkg = fRunInfo->fBkgFix[0]; @@ -185,7 +185,7 @@ double PRunSingleHisto::CalcMaxLikelihood(const std::vector& par) // get background double bkg; if (fRunInfo->fBkgFitParamNo == -1) { // bkg not fitted - if (isnan(fRunInfo->fBkgFix[0])) { // no fixed background given (background interval) + if (!fRunInfo->fBkgFixPresent[0]) { // no fixed background given (background interval) bkg = fBackground; } else { // fixed bkg given bkg = fRunInfo->fBkgFix[0]; @@ -259,7 +259,7 @@ void PRunSingleHisto::CalcTheory() // get background double bkg; if (fRunInfo->fBkgFitParamNo == -1) { // bkg not fitted - if (isnan(fRunInfo->fBkgFix[0])) { // no fixed background given (background interval) + if (!fRunInfo->fBkgFixPresent[0]) { // no fixed background given (background interval) bkg = fBackground; } else { // fixed bkg given bkg = fRunInfo->fBkgFix[0]; @@ -412,7 +412,7 @@ bool PRunSingleHisto::PrepareFitData() // check how the background shall be handled if (fRunInfo->fBkgFitParamNo == -1) { // bkg shall **NOT** be fitted // subtract background from histogramms ------------------------------------------ - if (isnan(fRunInfo->fBkgFix[0])) { // no fixed background given + if (!fRunInfo->fBkgFixPresent[0]) { // no fixed background given if (fRunInfo->fBkgRange[0] != 0) { if (!EstimateBkg(histoNo)) return false; @@ -628,7 +628,7 @@ cout << endl << ">> data start time = " << fData.fDataTimeStart; // get background double bkg; if (fRunInfo->fBkgFitParamNo == -1) { // bkg not fitted - if (isnan(fRunInfo->fBkgFix[0])) { // no fixed background given (background interval) + if (!fRunInfo->fBkgFixPresent[0]) { // no fixed background given (background interval) if (!EstimateBkg(histoNo)) return false; bkg = fBackground; @@ -800,7 +800,7 @@ bool PRunSingleHisto::PrepareViewData() // get background double bkg; if (fRunInfo->fBkgFitParamNo == -1) { // bkg not fitted - if (isnan(fRunInfo->fBkgFix[0])) { // no fixed background given (background interval) + if (!fRunInfo->fBkgFixPresent[0]) { // no fixed background given (background interval) if (!EstimateBkg(histoNo)) return false; bkg = fBackground; diff --git a/src/include/PFitter.h b/src/include/PFitter.h index 921eaea7..bcbac3b5 100644 --- a/src/include/PFitter.h +++ b/src/include/PFitter.h @@ -87,6 +87,7 @@ class PFitter bool CheckCommands(); bool SetParameters(); + bool ExecuteHesse(); bool ExecuteMigrad(); bool ExecuteMinimize(); bool ExecuteMinos(); diff --git a/src/include/PMusr.h b/src/include/PMusr.h index eed04aac..61edf595 100644 --- a/src/include/PMusr.h +++ b/src/include/PMusr.h @@ -234,6 +234,7 @@ typedef struct { PIntVector fMap; ///< int fForwardHistoNo; ///< int fBackwardHistoNo; ///< + bool fBkgFixPresent[2]; ///< flag showing if a fixed background is present double fBkgFix[2]; ///< int fBkgRange[4]; ///< int fDataRange[4]; ///<