47 #include "Minuit2/FunctionMinimum.h" 48 #include "Minuit2/MnContours.h" 49 #include "Minuit2/MnHesse.h" 50 #include "Minuit2/MnMinimize.h" 51 #include "Minuit2/MnMigrad.h" 52 #include "Minuit2/MnMinos.h" 53 #include "Minuit2/MnPlot.h" 54 #include "Minuit2/MnPrint.h" 55 #include "Minuit2/MnScan.h" 56 #include "Minuit2/MnSimplex.h" 57 #include "Minuit2/MnUserParameterState.h" 58 #include "Minuit2/MinosError.h" 65 #include <TObjArray.h> 66 #include <TObjString.h> 110 std::cerr <<
"**WARNING** from PSectorChisq::SetRunFirstTime. It tries to set" << std::endl;
111 std::cerr <<
" a fgb time stamp with idx=" << idx <<
" which is larger than #RUNS=" <<
fNoOfRuns <<
"." << std::endl;
112 std::cerr <<
" Will ignore it, but you better check what is going on!" << std::endl;
131 std::cerr <<
"**WARNING** from PSectorChisq::SetChisq. It tries to set" << std::endl;
132 std::cerr <<
" a chisq with idx=" << idx <<
" which is larger than #RUNS=" <<
fNoOfRuns <<
"." << std::endl;
133 std::cerr <<
" Will ignore it, but you better check what is going on!" << std::endl;
152 std::cerr <<
"**WARNING** from PSectorChisq::SetExpectedChisq. It tries to set" << std::endl;
153 std::cerr <<
" a chisq with idx=" << idx <<
" which is larger than #RUNS=" <<
fNoOfRuns <<
"." << std::endl;
154 std::cerr <<
" Will ignore it, but you better check what is going on!" << std::endl;
173 std::cerr <<
"**WARNING** from PSectorChisq::SetNDF. It tries to set" << std::endl;
174 std::cerr <<
" a NDF with idx=" << idx <<
" which is larger than #RUNS=" <<
fNoOfRuns <<
"." << std::endl;
175 std::cerr <<
" Will ignore it, but you better check what is going on!" << std::endl;
268 fChisqOnly(chisq_only), fRunInfo(runInfo), fRunListCollection(runListCollection)
299 for (UInt_t i=0; i<runs->size(); i++) {
300 range.first = (*runs)[i].GetFitRange(0);
301 range.second = (*runs)[i].GetFitRange(1);
345 for (
unsigned int i=0; i<
fPhase.size(); i++)
352 TObjArray *tok =
nullptr;
353 TObjString *ostr =
nullptr;
356 for (
unsigned int i=0; i<theo->size(); i++) {
358 TString line = theo->at(i).fLine;
359 if (line.Contains(
"TFieldCos") || line.Contains(
"tf ") ||
360 line.Contains(
"bessel") || line.Contains(
"b ") ||
361 line.Contains(
"skewedGss") || line.Contains(
"skg ") ||
362 line.Contains(
"staticNKTF") || line.Contains(
"snktf ") ||
363 line.Contains(
"dynamicNKTF") || line.Contains(
"dnktf ")) {
366 if (line.Contains(
"internFld") || line.Contains(
"if ") ||
367 line.Contains(
"internBsl") || line.Contains(
"ib ")) {
370 if (line.Contains(
"muMinusExpTF") || line.Contains(
"mmsetf ")) {
378 tok = line.Tokenize(
" \t");
379 if (tok ==
nullptr) {
380 std::cerr <<
"PFitter::GetPhaseParams(): **ERROR** couldn't tokenize theory line string." << std::endl;
383 if (tok->GetEntries() > pos) {
384 ostr =
dynamic_cast<TObjString*
>(tok->At(pos));
385 str = ostr->GetString();
392 if (str.Contains(
"fun")) {
394 for (
int i=0; i<parVec.size(); i++) {
396 fPhase[parVec[i]-1] =
true;
398 }
else if (str.Contains(
"map")) {
400 for (
int i=0; i<parVec.size(); i++) {
402 fPhase[parVec[i]-1] =
true;
405 int idx = str.Atoi();
407 std::cerr <<
"PFitter::GetPhaseParams(): **ERROR** str=" << str.View() <<
" is not an integer!" << std::endl;
412 std::cerr <<
"PFitter::GetPhaseParams(): **ERROR** idx=" << idx <<
" is > #param = " <<
fRunInfo->
GetNoOfParams() <<
"!" << std::endl;
435 TObjArray *tok =
nullptr;
436 TObjString *ostr =
nullptr;
439 for (
int i=0; i<funList->size(); i++) {
440 if (funList->at(i).fLine.Contains(funStr)) {
442 tok = funList->at(i).fLine.Tokenize(
" =+-*/");
443 if (tok ==
nullptr) {
444 std::cerr <<
"PFitter::GetParFromFun(): **ERROR** couldn't tokenize function string." << std::endl;
448 for (
int j=1; j<tok->GetEntries(); j++) {
449 ostr =
dynamic_cast<TObjString*
>(tok->At(j));
450 str = ostr->GetString();
452 if (str.Contains(
"par")) {
454 Ssiz_t idx = str.Index(
"par");
459 }
while (isdigit(str[idx++]));
460 parVec.push_back(parStr.Atoi());
463 if (str.Contains(
"map")) {
465 Ssiz_t idx = str.Index(
"map");
467 TString mapStr(
"map");
470 }
while (isdigit(str[idx++]));
472 for (
int k=0; k<mapParVec.size(); k++) {
473 parVec.push_back(mapParVec[k]);
501 TString str = mapStr;
506 std::cerr <<
"PFitter::GetParFromMap(): **ERROR** couldn't get propper index from mapX!" << std::endl;
514 if (runList ==
nullptr) {
515 std::cerr <<
"PFitter::GetParFromMap(): **ERROR** couldn't get required run list information!" << std::endl;
520 for (
int i=0; i<runList->size(); i++) {
521 map = runList->at(i).GetMap();
522 if (map ==
nullptr) {
523 std::cerr <<
"PFitter::GetParFromMap(): **ERROR** couldn't get required map information (idx=" << i <<
")!" << std::endl;
527 if (idx >= map->size()) {
528 std::cerr <<
"PFitter::GetParFromMap(): **ERROR** requested map index (idx=" << idx <<
") out-of-range (" << map->size() <<
")!" << std::endl;
532 parVec.push_back(map->at(idx));
555 Int_t usedParams = 0;
556 for (UInt_t i=0; i<error.size(); i++) {
560 UInt_t ndf =
static_cast<int>(
fFitterFcn->GetTotalNoOfFittedBins()) - usedParams;
561 Double_t val = (*fFitterFcn)(param);
564 Double_t totalExpectedChisq = 0.0;
566 fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerRun);
568 std::vector<Double_t> chisqPerRun;
573 std::cout << std::endl << std::endl <<
">> chisq = " << val <<
", NDF = " << ndf <<
", chisq/NDF = " << val/ndf;
575 if (totalExpectedChisq != 0.0) {
576 std::cout << std::endl <<
">> expected chisq = " << totalExpectedChisq <<
", NDF = " << ndf <<
", expected chisq/NDF = " << totalExpectedChisq/ndf;
578 for (UInt_t i=0; i<expectedChisqPerRun.size(); i++) {
581 std::cout << std::endl <<
">> run block " << i+1 <<
": (NDF/red.chisq/red.chisq_e) = (" << ndf_run <<
"/" << chisqPerRun[i]/ndf_run <<
"/" << expectedChisqPerRun[i]/ndf_run <<
")";
583 }
else if (chisqPerRun.size() > 0) {
585 for (UInt_t i=0; i<chisqPerRun.size(); i++) {
588 std::cout << std::endl <<
">> run block " << i+1 <<
": (NDF/red.chisq) = (" << ndf_run <<
"/" << chisqPerRun[i]/ndf_run <<
")";
594 expectedChisqPerRun.clear();
598 secFitRange.resize(1);
599 for (UInt_t k=0; k<
fSector.size(); k++) {
601 secFitRange[0].first =
fSector[k].GetTimeRangeFirst(0);
602 secFitRange[0].second =
fSector[k].GetTimeRangeLast();
605 val = (*fFitterFcn)(param);
607 ndf =
static_cast<UInt_t
>(
fFitterFcn->GetTotalNoOfFittedBins()) - usedParams;
609 totalExpectedChisq = 0.0;
610 fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerRun);
616 std::cout << std::endl;
617 std::cout <<
"++++" << std::endl;
618 std::cout <<
">> Sector " << k+1 <<
": FitRange: " << secFitRange[0].first <<
", " << secFitRange[0].second << std::endl;
619 std::cout <<
">> chisq = " << val <<
", NDF = " << ndf <<
", chisq/NDF = " << val/ndf;
621 if (totalExpectedChisq != 0.0) {
622 std::cout << std::endl <<
">> expected chisq = " << totalExpectedChisq <<
", NDF = " << ndf <<
", expected chisq/NDF = " << totalExpectedChisq/ndf;
624 for (UInt_t i=0; i<expectedChisqPerRun.size(); i++) {
627 std::cout << std::endl <<
">> run block " << i+1 <<
": (NDF/red.chisq/red.chisq_e) = (" << ndf_run <<
"/" << chisqPerRun[i]/ndf_run <<
"/" << expectedChisqPerRun[i]/ndf_run <<
")";
629 }
else if (chisqPerRun.size() > 0) {
631 for (UInt_t i=0; i<chisqPerRun.size(); i++) {
634 std::cout << std::endl <<
">> run block " << i+1 <<
": (NDF/red.chisq) = (" << ndf_run <<
"/" << chisqPerRun[i]/ndf_run <<
")";
639 expectedChisqPerRun.clear();
644 Double_t totalExpectedMaxLH = 0.0;
645 std::vector<Double_t> expectedMaxLHPerRun;
646 fFitterFcn->CalcExpectedChiSquare(param, totalExpectedMaxLH, expectedMaxLHPerRun);
648 std::vector<Double_t> maxLHPerRun;
653 std::cout << std::endl << std::endl <<
">> maxLH = " << val <<
", NDF = " << ndf <<
", maxLH/NDF = " << val/ndf;
655 if (totalExpectedMaxLH != 0.0) {
656 std::cout << std::endl <<
">> expected maxLH = " << totalExpectedMaxLH <<
", NDF = " << ndf <<
", expected maxLH/NDF = " << totalExpectedMaxLH/ndf;
658 for (UInt_t i=0; i<expectedMaxLHPerRun.size(); i++) {
661 std::cout << std::endl <<
">> run block " << i+1 <<
": (NDF/red.maxLH/red.maxLH_e) = (" << ndf_run <<
"/" << maxLHPerRun[i]/ndf_run <<
"/" << expectedMaxLHPerRun[i]/ndf_run <<
")";
667 expectedMaxLHPerRun.clear();
671 secFitRange.resize(1);
672 for (UInt_t k=0; k<
fSector.size(); k++) {
674 secFitRange[0].first =
fSector[k].GetTimeRangeFirst(0);
675 secFitRange[0].second =
fSector[k].GetTimeRangeLast();
678 val = (*fFitterFcn)(param);
680 ndf =
static_cast<int>(
fFitterFcn->GetTotalNoOfFittedBins()) - usedParams;
682 totalExpectedMaxLH = 0.0;
683 fFitterFcn->CalcExpectedChiSquare(param, totalExpectedMaxLH, expectedMaxLHPerRun);
689 std::cout << std::endl;
690 std::cout <<
"++++" << std::endl;
691 std::cout <<
">> Sector " << k+1 <<
": FitRange: " << secFitRange[0].first <<
", " << secFitRange[0].second << std::endl;
692 std::cout <<
">> maxLH = " << val <<
", NDF = " << ndf <<
", maxLH/NDF = " << val/ndf;
694 if (totalExpectedMaxLH != 0.0) {
695 std::cout << std::endl <<
">> expected maxLH = " << totalExpectedMaxLH <<
", NDF = " << ndf <<
", expected maxLH/NDF = " << totalExpectedMaxLH/ndf;
697 for (UInt_t i=0; i<expectedMaxLHPerRun.size(); i++) {
700 std::cout << std::endl <<
">> run block " << i+1 <<
": (NDF/red.maxLH/red.maxLH_e) = (" << ndf_run <<
"/" << maxLHPerRun[i]/ndf_run <<
"/" << expectedMaxLHPerRun[i]/ndf_run <<
")";
706 expectedMaxLHPerRun.clear();
710 std::cout << std::endl << std::endl;
716 std::cout << std::endl <<
">> Number of available threads for the function optimization: " << omp_get_max_threads() << std::endl;
721 std::cout << std::endl <<
">> Chi Square fit will be executed" << std::endl;
723 std::cout << std::endl <<
">> Maximum Likelihood fit will be executed" << std::endl;
727 for (UInt_t i=0; i<
fParams.size(); i++) {
732 Bool_t firstSave =
true;
733 for (UInt_t i=0; i<
fCmdList.size(); i++) {
736 std::cerr << std::endl <<
"**WARNING** from PFitter::DoFit() : the command INTERACTIVE is not yet implemented.";
737 std::cerr << std::endl;
749 std::cerr << std::endl <<
"**WARNING** from PFitter::DoFit() : the command EIGEN is not yet implemented.";
750 std::cerr << std::endl;
756 std::cerr << std::endl <<
"**WARNING** from PFitter::DoFit() : the command MACHINE_PRECISION is not yet implemented.";
757 std::cerr << std::endl;
792 std::cerr << std::endl <<
"**WARNING** from PFitter::DoFit() : the command USER_COVARIANCE is not yet implemented.";
793 std::cerr << std::endl;
796 std::cerr << std::endl <<
"**WARNING** from PFitter::DoFit() : the command USER_PARAM_STATE is not yet implemented.";
797 std::cerr << std::endl;
803 std::cerr << std::endl <<
"**PANIC ERROR**: PFitter::DoFit(): You should never have reached this point";
804 std::cerr << std::endl;
840 PMsrLines::iterator it;
841 UInt_t cmdLineNo = 0;
852 pos = line.First(
'#');
854 line.Remove(pos,line.Length()-pos);
856 if (line.Contains(
"COMMANDS", TString::kIgnoreCase)) {
858 }
else if (line.Contains(
"SET BATCH", TString::kIgnoreCase)) {
860 }
else if (line.Contains(
"END RETURN", TString::kIgnoreCase)) {
862 }
else if (line.Contains(
"CHI_SQUARE", TString::kIgnoreCase)) {
864 }
else if (line.Contains(
"MAX_LIKELIHOOD", TString::kIgnoreCase)) {
866 }
else if (line.Contains(
"SCALE_N0_BKG", TString::kIgnoreCase)) {
868 }
else if (line.Contains(
"INTERACTIVE", TString::kIgnoreCase)) {
870 cmd.second = cmdLineNo;
872 }
else if (line.Contains(
"CONTOURS", TString::kIgnoreCase)) {
874 cmd.second = cmdLineNo;
877 TObjArray *tokens =
nullptr;
882 tokens = line.Tokenize(
", \t");
884 for (Int_t i=0; i<tokens->GetEntries(); i++) {
885 ostr =
dynamic_cast<TObjString*
>(tokens->At(i));
886 str = ostr->GetString();
888 if ((i==1) || (i==2)) {
890 if (!str.IsDigit()) {
891 std::cerr << std::endl <<
">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
892 std::cerr << std::endl <<
">> " << line.Data();
893 std::cerr << std::endl <<
">> parameter number is not number!";
894 std::cerr << std::endl <<
">> command syntax for CONTOURS is: CONTOURS parameter-X parameter-Y [# of points]";
895 std::cerr << std::endl;
905 if ((ival < 1) || (ival >
fParams.size())) {
906 std::cerr << std::endl <<
">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
907 std::cerr << std::endl <<
">> " << line.Data();
908 std::cerr << std::endl <<
">> parameter number is out of range [1," <<
fParams.size() <<
"]!";
909 std::cerr << std::endl <<
">> command syntax for CONTOURS is: CONTOURS parameter-X parameter-Y [# of points]";
910 std::cerr << std::endl;
925 if (!str.IsDigit()) {
926 std::cerr << std::endl <<
">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
927 std::cerr << std::endl <<
">> " << line.Data();
928 std::cerr << std::endl <<
">> number of points is not number!";
929 std::cerr << std::endl <<
">> command syntax for CONTOURS is: CONTOURS parameter-X parameter-Y [# of points]";
930 std::cerr << std::endl;
939 if ((ival < 1) || (ival > 100)) {
940 std::cerr << std::endl <<
">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
941 std::cerr << std::endl <<
">> " << line.Data();
942 std::cerr << std::endl <<
">> number of scan points is out of range [1,100]!";
943 std::cerr << std::endl <<
">> command syntax for CONTOURS is: CONTOURS parameter-X parameter-Y [# of points]";
944 std::cerr << std::endl;
960 }
else if (line.Contains(
"EIGEN", TString::kIgnoreCase)) {
962 cmd.second = cmdLineNo;
964 }
else if (line.Contains(
"FIT_RANGE", TString::kIgnoreCase)) {
971 TObjArray *tokens =
nullptr;
975 tokens = line.Tokenize(
", \t");
977 if (tokens->GetEntries() == 2) {
978 ostr =
dynamic_cast<TObjString*
>(tokens->At(1));
979 str = ostr->GetString();
980 if (str.Contains(
"RESET", TString::kIgnoreCase)) {
982 cmd.second = cmdLineNo;
985 std::cerr << std::endl <<
">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
986 std::cerr << std::endl <<
">> " << line.Data();
987 std::cerr << std::endl <<
">> Syntax: FIT_RANGE RESET | FIT_RANGE start end | FIT_RANGE s1 e1 s2 e2 .. sN eN,";
988 std::cerr << std::endl <<
">> with N the number of runs in the msr-file." << std::endl;
989 std::cerr << std::endl <<
">> Found " << str.Data() <<
", instead of RESET" << std::endl;
997 }
else if ((tokens->GetEntries() > 1) && (static_cast<UInt_t>(tokens->GetEntries()) % 2) == 1) {
998 if ((tokens->GetEntries() > 3) && ((static_cast<UInt_t>(tokens->GetEntries())-1)) != 2*
fRunInfo->
GetMsrRunList()->size()) {
999 std::cerr << std::endl <<
">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
1000 std::cerr << std::endl <<
">> " << line.Data();
1001 std::cerr << std::endl <<
">> Syntax: FIT_RANGE RESET | FIT_RANGE <start> <end> | FIT_RANGE <s1> <e1> <s2> <e2> .. <sN> <eN> |";
1002 std::cerr << std::endl <<
">> FIT_RANGE fgb+<n0> lgb-<n1> | FIT_RANGE fgb+<n00> lgb-<n01> fgb+<n10> lgb-<n11> ... fgb+<nN0> lgb-<nN1>,";
1003 std::cerr << std::endl <<
">> with N the number of runs in the msr-file.";
1004 std::cerr << std::endl <<
">> Found N=" << (tokens->GetEntries()-1)/2 <<
", # runs in msr-file=" <<
fRunInfo->
GetMsrRunList()->size() << std::endl;
1014 for (Int_t n=1; n<tokens->GetEntries(); n++) {
1015 ostr =
dynamic_cast<TObjString*
>(tokens->At(n));
1016 str = ostr->GetString();
1017 if (!str.IsFloat()) {
1018 if ((n%2 == 1) && (!str.Contains(
"fgb", TString::kIgnoreCase)))
1020 if ((n%2 == 0) && (!str.Contains(
"lgb", TString::kIgnoreCase)))
1029 cmd.second = cmdLineNo;
1032 std::cerr << std::endl <<
">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
1033 std::cerr << std::endl <<
">> " << line.Data();
1034 std::cerr << std::endl <<
">> Syntax: FIT_RANGE RESET | FIT_RANGE <start> <end> | FIT_RANGE <s1> <e1> <s2> <e2> .. <sN> <eN> |";
1035 std::cerr << std::endl <<
">> FIT_RANGE fgb+<n0> lgb-<n1> | FIT_RANGE fgb+<n00> lgb-<n01> fgb+<n10> lgb-<n11> ... fgb+<nN0> lgb-<nN1>,";
1036 std::cerr << std::endl <<
">> with N the number of runs in the msr-file.";
1037 std::cerr << std::endl <<
">> Found token '" << str.Data() <<
"', which is not a floating point number." << std::endl;
1047 std::cerr << std::endl <<
">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
1048 std::cerr << std::endl <<
">> " << line.Data();
1049 std::cerr << std::endl <<
">> Syntax: FIT_RANGE RESET | FIT_RANGE <start> <end> | FIT_RANGE <s1> <e1> <s2> <e2> .. <sN> <eN> |";
1050 std::cerr << std::endl <<
">> FIT_RANGE fgb+<n0> lgb-<n1> | FIT_RANGE fgb+<n00> lgb-<n01> fgb+<n10> lgb-<n11> ... fgb+<nN0> lgb-<nN1>,";
1051 std::cerr << std::endl <<
">> with N the number of runs in the msr-file.";
1064 }
else if (line.Contains(
"FIX", TString::kIgnoreCase)) {
1066 TObjArray *tokens =
nullptr;
1071 tokens = line.Tokenize(
", \t");
1073 for (Int_t i=1; i<tokens->GetEntries(); i++) {
1074 ostr =
dynamic_cast<TObjString*
>(tokens->At(i));
1075 str = ostr->GetString();
1077 if (str.IsDigit()) {
1081 std::cerr << std::endl <<
">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
1082 std::cerr << std::endl <<
">> " << line.Data();
1083 std::cerr << std::endl <<
">> Parameter " << ival <<
" is out of the Parameter Range [1," <<
fParams.size() <<
"]";
1084 std::cerr << std::endl;
1094 Bool_t found =
false;
1095 for (UInt_t j=0; j<
fParams.size(); j++) {
1096 if (
fParams[j].fName.CompareTo(str, TString::kIgnoreCase) == 0) {
1102 std::cerr << std::endl <<
">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
1103 std::cerr << std::endl <<
">> " << line.Data();
1104 std::cerr << std::endl <<
">> Parameter '" << str.Data() <<
"' is NOT present as a parameter name";
1105 std::cerr << std::endl;
1123 cmd.second = cmdLineNo;
1125 }
else if (line.Contains(
"HESSE", TString::kIgnoreCase)) {
1128 cmd.second = cmdLineNo;
1130 }
else if (line.Contains(
"MACHINE_PRECISION", TString::kIgnoreCase)) {
1132 cmd.second = cmdLineNo;
1134 }
else if (line.Contains(
"MIGRAD", TString::kIgnoreCase)) {
1137 cmd.second = cmdLineNo;
1139 }
else if (line.Contains(
"MINIMIZE", TString::kIgnoreCase)) {
1142 cmd.second = cmdLineNo;
1144 }
else if (line.Contains(
"MINOS", TString::kIgnoreCase)) {
1147 cmd.second = cmdLineNo;
1149 }
else if (line.Contains(
"MNPLOT", TString::kIgnoreCase)) {
1151 cmd.second = cmdLineNo;
1153 }
else if (line.Contains(
"PRINT_LEVEL", TString::kIgnoreCase)) {
1155 cmd.second = cmdLineNo;
1157 }
else if (line.Contains(
"RELEASE", TString::kIgnoreCase)) {
1159 TObjArray *tokens =
nullptr;
1164 tokens = line.Tokenize(
", \t");
1166 for (Int_t i=1; i<tokens->GetEntries(); i++) {
1167 ostr =
dynamic_cast<TObjString*
>(tokens->At(i));
1168 str = ostr->GetString();
1170 if (str.IsDigit()) {
1174 std::cerr << std::endl <<
">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
1175 std::cerr << std::endl <<
">> " << line.Data();
1176 std::cerr << std::endl <<
">> Parameter " << ival <<
" is out of the Parameter Range [1," <<
fParams.size() <<
"]";
1177 std::cerr << std::endl;
1187 Bool_t found =
false;
1188 for (UInt_t j=0; j<
fParams.size(); j++) {
1189 if (
fParams[j].fName.CompareTo(str, TString::kIgnoreCase) == 0) {
1195 std::cerr << std::endl <<
">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
1196 std::cerr << std::endl <<
">> " << line.Data();
1197 std::cerr << std::endl <<
">> Parameter '" << str.Data() <<
"' is NOT present as a parameter name";
1198 std::cerr << std::endl;
1214 cmd.second = cmdLineNo;
1216 }
else if (line.Contains(
"RESTORE", TString::kIgnoreCase)) {
1218 cmd.second = cmdLineNo;
1220 }
else if (line.Contains(
"SAVE", TString::kIgnoreCase)) {
1222 cmd.second = cmdLineNo;
1224 }
else if (line.Contains(
"SCAN", TString::kIgnoreCase)) {
1226 cmd.second = cmdLineNo;
1229 TObjArray *tokens =
nullptr;
1234 tokens = line.Tokenize(
", \t");
1236 for (Int_t i=0; i<tokens->GetEntries(); i++) {
1237 ostr =
dynamic_cast<TObjString*
>(tokens->At(i));
1238 str = ostr->GetString();
1241 if (!str.IsDigit()) {
1242 std::cerr << std::endl <<
">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
1243 std::cerr << std::endl <<
">> " << line.Data();
1244 std::cerr << std::endl <<
">> parameter number is not number!";
1245 std::cerr << std::endl <<
">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]";
1246 std::cerr << std::endl;
1256 if ((ival < 1) || (ival >
fParams.size())) {
1257 std::cerr << std::endl <<
">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
1258 std::cerr << std::endl <<
">> " << line.Data();
1259 std::cerr << std::endl <<
">> parameter number is out of range [1," <<
fParams.size() <<
"]!";
1260 std::cerr << std::endl <<
">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]";
1261 std::cerr << std::endl;
1276 if (!str.IsDigit()) {
1277 std::cerr << std::endl <<
">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
1278 std::cerr << std::endl <<
">> " << line.Data();
1279 std::cerr << std::endl <<
">> number of points is not number!";
1280 std::cerr << std::endl <<
">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]";
1281 std::cerr << std::endl;
1290 if ((ival < 1) || (ival > 100)) {
1291 std::cerr << std::endl <<
">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
1292 std::cerr << std::endl <<
">> " << line.Data();
1293 std::cerr << std::endl <<
">> number of scan points is out of range [1,100]!";
1294 std::cerr << std::endl <<
">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]";
1295 std::cerr << std::endl;
1308 if (!str.IsFloat()) {
1309 std::cerr << std::endl <<
">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
1310 std::cerr << std::endl <<
">> " << line.Data();
1311 std::cerr << std::endl <<
">> low is not a floating point number!";
1312 std::cerr << std::endl <<
">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]";
1313 std::cerr << std::endl;
1326 if (!str.IsFloat()) {
1327 std::cerr << std::endl <<
">> PFitter::CheckCommands: **ERROR** in line " << it->fLineNo;
1328 std::cerr << std::endl <<
">> " << line.Data();
1329 std::cerr << std::endl <<
">> high is not a floating point number!";
1330 std::cerr << std::endl <<
">> command syntax for SCAN is: SCAN [parameter no [# of points [low high]]]";
1331 std::cerr << std::endl;
1347 }
else if (line.Contains(
"SIMPLEX", TString::kIgnoreCase)) {
1349 cmd.second = cmdLineNo;
1351 }
else if (line.Contains(
"STRATEGY", TString::kIgnoreCase)) {
1352 TObjArray *tokens =
nullptr;
1356 tokens = line.Tokenize(
" \t");
1357 if (tokens->GetEntries() == 2) {
1358 ostr =
dynamic_cast<TObjString*
>(tokens->At(1));
1359 str = ostr->GetString();
1360 if (str.CompareTo(
"0") == 0) {
1362 }
else if (str.CompareTo(
"1") == 0) {
1364 }
else if (str.CompareTo(
"2") == 0) {
1366 }
else if (str.CompareTo(
"LOW") == 0) {
1368 }
else if (str.CompareTo(
"DEFAULT") == 0) {
1370 }
else if (str.CompareTo(
"HIGH") == 0) {
1379 }
else if (line.Contains(
"USER_COVARIANCE", TString::kIgnoreCase)) {
1381 cmd.second = cmdLineNo;
1383 }
else if (line.Contains(
"USER_PARAM_STATE", TString::kIgnoreCase)) {
1385 cmd.second = cmdLineNo;
1387 }
else if (line.Contains(
"SECTOR", TString::kIgnoreCase)) {
1390 cmd.second = cmdLineNo;
1394 TObjArray *tokens =
nullptr;
1398 tokens = line.Tokenize(
" ,\t");
1400 if (tokens->GetEntries() == 1) {
1401 std::cerr << std::endl <<
">> PFitter::CheckCommands(): **FATAL ERROR** in line " << it->fLineNo;
1402 std::cerr << std::endl <<
">> " << line.Data();
1403 std::cerr << std::endl <<
">> At least one sector time stamp is expected.";
1404 std::cerr << std::endl <<
">> Will stop ...";
1405 std::cerr << std::endl;
1417 for (Int_t i=1; i<tokens->GetEntries(); i++) {
1421 ostr =
dynamic_cast<TObjString*
>(tokens->At(i));
1422 str = ostr->GetString();
1423 if (str.IsFloat()) {
1428 std::cerr << std::endl <<
">> PFitter::CheckCommands(): **FATAL ERROR** in line " << it->fLineNo;
1429 std::cerr << std::endl <<
">> " << line.Data();
1430 std::cerr << std::endl <<
">> The sector time stamp " << dval <<
" is > as the lgb time stamp (" <<
fOriginalFitRange[j].second <<
") of run " << j <<
".";
1431 std::cerr << std::endl <<
">> Will stop ...";
1432 std::cerr << std::endl;
1444 sec.SetSectorTime(dval);
1447 std::cerr << std::endl <<
">> PFitter::CheckCommands(): **FATAL ERROR** in line " << it->fLineNo;
1448 std::cerr << std::endl <<
">> " << line.Data();
1449 std::cerr << std::endl <<
">> The sector time stamp '" << str <<
"' is not a number.";
1450 std::cerr << std::endl <<
">> Will stop ...";
1451 std::cerr << std::endl;
1468 std::cerr << std::endl <<
">> PFitter::CheckCommands(): **FATAL ERROR** in line " << it->fLineNo <<
" an unkown command is found:";
1469 std::cerr << std::endl <<
">> " << line.Data();
1470 std::cerr << std::endl <<
">> Will stop ...";
1471 std::cerr << std::endl;
1479 Bool_t fixFlag =
false;
1480 Bool_t releaseFlag =
false;
1481 Bool_t minimizerFlag =
false;
1483 if (line.Contains(
"FIX", TString::kIgnoreCase))
1485 else if (line.Contains(
"RELEASE", TString::kIgnoreCase) ||
1486 line.Contains(
"RESTORE", TString::kIgnoreCase))
1488 else if (line.Contains(
"MINIMIZE", TString::kIgnoreCase) ||
1489 line.Contains(
"MIGRAD", TString::kIgnoreCase) ||
1490 line.Contains(
"SIMPLEX", TString::kIgnoreCase)) {
1492 minimizerFlag =
true;
1493 }
else if (line.Contains(
"MINOS", TString::kIgnoreCase)) {
1494 if (fixFlag && releaseFlag && !minimizerFlag) {
1495 std::cerr << std::endl <<
">> PFitter::CheckCommands(): **WARNING** RELEASE/RESTORE command present";
1496 std::cerr << std::endl <<
">> without minimizer command (MINIMIZE/MIGRAD/SIMPLEX) between";
1497 std::cerr << std::endl <<
">> RELEASE/RESTORE and MINOS. Behaviour might be different to the";
1498 std::cerr << std::endl <<
">> expectation of the user ?!?" << std::endl;
1501 releaseFlag =
false;
1502 minimizerFlag =
false;
1520 for (UInt_t i=0; i<
fParams.size(); i++) {
1522 if (
fParams[i].fStep == 0.0) {
1526 if (
fParams[i].fNoOfParams > 5) {
1527 if (
fParams[i].fLowerBoundaryPresent &&
1528 fParams[i].fUpperBoundaryPresent) {
1531 }
else if (
fParams[i].fLowerBoundaryPresent &&
1532 !
fParams[i].fUpperBoundaryPresent) {
1546 for (UInt_t i=0; i<
fParams.size(); i++) {
1550 std::cerr << std::endl <<
"**WARNING** : Parameter No " << i+1 <<
" is not used at all, will fix it";
1551 std::cerr << std::endl;
1568 std::cout <<
">> PFitter::ExecuteContours() ..." << std::endl;
1572 std::cerr << std::endl <<
"**WARNING**: CONTOURS musn't be called before any minimization (MINIMIZE/MIGRAD/SIMPLEX) is done!!";
1573 std::cerr << std::endl;
1579 std::cerr << std::endl <<
"**ERROR**: CONTOURS cannot started since the previous minimization failed :-(";
1580 std::cerr << std::endl;
1603 std::cout <<
">> PFitter::ExecuteFitRange(): " <<
fCmdLines[lineNo].fLine.Data() << std::endl;
1605 if (
fCmdLines[lineNo].fLine.Contains(
"fgb", TString::kIgnoreCase)) {
1610 TObjArray *tokens =
nullptr;
1614 tokens =
fCmdLines[lineNo].fLine.Tokenize(
", \t");
1619 if (tokens->GetEntries() == 2) {
1621 }
else if (tokens->GetEntries() == 3) {
1622 Double_t start = 0.0, end = 0.0;
1626 ostr =
dynamic_cast<TObjString*
>(tokens->At(1));
1627 str = ostr->GetString();
1629 ostr =
dynamic_cast<TObjString*
>(tokens->At(2));
1630 str = ostr->GetString();
1633 fitRange.first = start;
1634 fitRange.second = end;
1635 fitRangeVector.push_back(fitRange);
1639 Double_t start = 0.0, end = 0.0;
1643 for (UInt_t i=0; i<runList->size(); i++) {
1644 ostr =
dynamic_cast<TObjString*
>(tokens->At(2*i+1));
1645 str = ostr->GetString();
1647 ostr =
dynamic_cast<TObjString*
>(tokens->At(2*i+2));
1648 str = ostr->GetString();
1651 fitRange.first = start;
1652 fitRange.second = end;
1653 fitRangeVector.push_back(fitRange);
1674 std::cout <<
">> PFitter::ExecuteFix(): " <<
fCmdLines[lineNo].fLine.Data() << std::endl;
1676 TObjArray *tokens =
nullptr;
1680 tokens =
fCmdLines[lineNo].fLine.Tokenize(
", \t");
1682 for (Int_t i=1; i<tokens->GetEntries(); i++) {
1683 ostr =
dynamic_cast<TObjString*
>(tokens->At(i));
1684 str = ostr->GetString();
1686 if (str.IsDigit()) {
1712 std::cout <<
">> PFitter::ExecuteHesse(): will call hesse ..." << std::endl;
1715 ROOT::Minuit2::MnHesse hesse;
1718 UInt_t maxfcn = std::numeric_limits<UInt_t>::max();
1721 Double_t start=0.0, end=0.0;
1725 std::cout <<
">> PFitter::ExecuteMinimize(): execution time for Hesse = " << std::setprecision(3) << (end-start)/1.0e3 <<
" sec." << std::endl;
1726 TString str = TString::Format(
"Hesse: %.3f sec", (end-start)/1.0e3);
1728 if (!mnState.IsValid()) {
1729 std::cerr << std::endl <<
">> PFitter::ExecuteHesse(): **WARNING** Hesse encountered a problem! The state found is invalid.";
1730 std::cerr << std::endl;
1733 if (!mnState.HasCovariance()) {
1734 std::cerr << std::endl <<
">> PFitter::ExecuteHesse(): **WARNING** Hesse encountered a problem! No covariance matrix available.";
1735 std::cerr << std::endl;
1740 for (UInt_t i=0; i<
fParams.size(); i++) {
1746 std::cout << mnState << std::endl;
1761 std::cout <<
">> PFitter::ExecuteMigrad(): will call migrad ..." << std::endl;
1769 UInt_t maxfcn = std::numeric_limits<UInt_t>::max();
1771 Double_t tolerance = 0.1;
1773 Double_t start=0.0, end=0.0;
1775 ROOT::Minuit2::FunctionMinimum min = migrad(maxfcn, tolerance);
1777 std::cout <<
">> PFitter::ExecuteMinimize(): execution time for Migrad = " << std::setprecision(3) << (end-start)/1.0e3 <<
" sec." << std::endl;
1778 TString str = TString::Format(
"Migrad: %.3f sec", (end-start)/1.0e3);
1780 if (!min.IsValid()) {
1781 std::cerr << std::endl <<
">> PFitter::ExecuteMigrad(): **WARNING**: Fit did not converge, sorry ...";
1782 std::cerr << std::endl;
1788 fFcnMin.reset(
new ROOT::Minuit2::FunctionMinimum(min));
1795 for (UInt_t i=0; i<
fParams.size(); i++) {
1796 Double_t dval = min.UserState().Value(i);
1798 Int_t m = (Int_t)(dval/360.0);
1799 dval = dval - m*360.0;
1807 Double_t minVal = min.Fval();
1808 UInt_t ndf =
fFitterFcn->GetTotalNoOfFittedBins();
1810 for (UInt_t i=0; i<
fParams.size(); i++) {
1811 if ((min.UserState().Error(i) != 0.0) && !
fMnUserParams.Parameters().at(i).IsFixed())
1822 std::cout << *
fFcnMin << std::endl;
1837 std::cout <<
">> PFitter::ExecuteMinimize(): will call minimize ..." << std::endl;
1845 UInt_t maxfcn = std::numeric_limits<UInt_t>::max();
1848 Double_t tolerance = 0.1;
1850 Double_t start=0.0, end=0.0;
1852 ROOT::Minuit2::FunctionMinimum min = minimize(maxfcn, tolerance);
1854 std::cout <<
">> PFitter::ExecuteMinimize(): execution time for Minimize = " << std::setprecision(3) << (end-start)/1.0e3 <<
" sec." << std::endl;
1855 TString str = TString::Format(
"Minimize: %.3f sec", (end-start)/1.0e3);
1857 if (!min.IsValid()) {
1858 std::cerr << std::endl <<
">> PFitter::ExecuteMinimize(): **WARNING**: Fit did not converge, sorry ...";
1859 std::cerr << std::endl;
1865 fFcnMin.reset(
new ROOT::Minuit2::FunctionMinimum(min));
1872 for (UInt_t i=0; i<
fParams.size(); i++) {
1873 Double_t dval = min.UserState().Value(i);
1875 Int_t m = (Int_t)(dval/360.0);
1876 dval = dval - m*360.0;
1884 Double_t minVal = min.Fval();
1885 UInt_t ndf =
fFitterFcn->GetTotalNoOfFittedBins();
1887 for (UInt_t i=0; i<
fParams.size(); i++) {
1888 if ((min.UserState().Error(i) != 0.0) && !
fMnUserParams.Parameters().at(i).IsFixed())
1899 std::cout << *
fFcnMin << std::endl;
1914 std::cout <<
">> PFitter::ExecuteMinos(): will call minos ..." << std::endl;
1918 std::cerr << std::endl <<
"**ERROR**: MINOS musn't be called before any minimization (MINIMIZE/MIGRAD/SIMPLEX) is done!!";
1919 std::cerr << std::endl;
1925 std::cerr << std::endl <<
"**ERROR**: MINOS cannot started since the previous minimization failed :-(";
1926 std::cerr << std::endl;
1931 Double_t start=0.0, end=0.0;
1935 for (UInt_t i=0; i<
fParams.size(); i++) {
1941 std::cout <<
">> PFitter::ExecuteMinos(): calculate errors for " <<
fParams[i].fName << std::endl;
1944 ROOT::Minuit2::MinosError err = minos.Minos(i);
1946 if (err.IsValid()) {
1957 std::cerr << std::endl <<
">> PFitter::ExecuteMinos(): **WARNING** Parameter " <<
fMnUserParams.Name(i) <<
" (ParamNo " << i+1 <<
") is fixed!";
1958 std::cerr << std::endl <<
">> Will set STEP to zero, i.e. making it a constant parameter";
1959 std::cerr << std::endl;
1966 std::cout <<
">> PFitter::ExecuteMinimize(): execution time for Minos = " << std::setprecision(3) << (end-start)/1.0e3 <<
" sec." << std::endl;
1967 TString str = TString::Format(
"Minos: %.3f sec", (end-start)/1.0e3);
1983 std::cout <<
">> PFitter::ExecutePlot() ..." << std::endl;
1985 ROOT::Minuit2::MnPlot plot;
2003 std::cout <<
">> PFitter::ExecutePrintLevel(): " <<
fCmdLines[lineNo].fLine.Data() << std::endl;
2005 TObjArray *tokens =
nullptr;
2009 tokens =
fCmdLines[lineNo].fLine.Tokenize(
", \t");
2011 if (tokens->GetEntries() < 2) {
2012 std::cerr << std::endl <<
"**ERROR** from PFitter::ExecutePrintLevel(): SYNTAX: PRINT_LEVEL <N>, where <N>=0-3" << std::endl << std::endl;
2016 ostr = (TObjString*)tokens->At(1);
2017 str = ostr->GetString();
2020 if (str.IsDigit()) {
2022 if ((ival >=0) && (ival <= 3)) {
2025 std::cerr << std::endl <<
"**ERROR** from PFitter::ExecutePrintLevel(): SYNTAX: PRINT_LEVEL <N>, where <N>=0-3";
2026 std::cerr << std::endl <<
" found <N>=" << ival << std::endl << std::endl;
2030 std::cerr << std::endl <<
"**ERROR** from PFitter::ExecutePrintLevel(): SYNTAX: PRINT_LEVEL <N>, where <N>=0-3" << std::endl << std::endl;
2034 #ifdef ROOT_GRTEQ_24 2035 ROOT::Minuit2::MnPrint::SetGlobalLevel(
fPrintLevel);
2061 TObjArray *tokens =
nullptr;
2065 tokens =
fCmdLines[lineNo].fLine.Tokenize(
", \t");
2067 std::cout <<
">> PFitter::ExecuteRelease(): " <<
fCmdLines[lineNo].fLine.Data() << std::endl;
2069 for (Int_t i=1; i<tokens->GetEntries(); i++) {
2070 ostr =
dynamic_cast<TObjString*
>(tokens->At(i));
2071 str = ostr->GetString();
2073 if (str.IsDigit()) {
2103 std::cout <<
"PFitter::ExecuteRestore(): release all fixed parameters (RESTORE) ..." << std::endl;
2105 for (UInt_t i=0; i<
fMnUserParams.Parameters().size(); i++) {
2125 std::cout <<
">> PFitter::ExecuteScan(): will call scan ..." << std::endl;
2154 std::cout << std::endl <<
">> PFitter::ExecuteSave(): nothing to be saved ...";
2158 ROOT::Minuit2::MnUserParameterState mnState =
fFcnMin->UserState();
2161 if (!mnState.IsValid()) {
2162 std::cerr << std::endl <<
">> PFitter::ExecuteSave: **WARNING** Minuit2 User Parameter State is not valid, i.e. nothing to be saved";
2163 std::cerr << std::endl;
2171 std::vector<Double_t> param;
2172 Double_t totalExpectedChisq = 0.0;
2173 std::vector<Double_t> expectedchisqPerRun;
2174 std::vector<UInt_t> ndfPerHisto;
2176 for (UInt_t i=0; i<
fParams.size(); i++)
2177 param.push_back(
fParams[i].fValue);
2180 fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedchisqPerRun);
2183 std::vector<Double_t> chisqPerRun;
2191 if (totalExpectedChisq != 0.0) {
2194 for (UInt_t i=0; i<expectedchisqPerRun.size(); i++) {
2196 ndfPerHisto.push_back(ndf_run);
2207 }
else if (chisqPerRun.size() > 1) {
2209 for (UInt_t i=0; i<chisqPerRun.size(); i++) {
2211 ndfPerHisto.push_back(ndf_run);
2225 for (UInt_t i=0; i<
fParams.size(); i++)
2226 error.push_back(
fParams[i].fStep);
2233 expectedchisqPerRun.clear();
2234 ndfPerHisto.clear();
2235 chisqPerRun.clear();
2237 std::cout <<
">> PFitter::ExecuteSave(): will write minuit2 output file ..." << std::endl;
2243 fout.open(
"MINUIT2.OUTPUT", std::iostream::out);
2245 fout.open(
"MINUIT2.OUTPUT", std::iostream::out | std::iostream::app);
2247 if (!fout.is_open()) {
2248 std::cerr << std::endl <<
"**ERROR** PFitter::ExecuteSave() couldn't open MINUIT2.OUTPUT file";
2249 std::cerr << std::endl;
2255 fout << std::endl <<
"*************************************************************************";
2256 fout << std::endl <<
" musrfit MINUIT2 output file from " <<
fRunInfo->
GetFileName().Data() <<
" - " << dt.AsSQLString();
2257 fout << std::endl <<
"*************************************************************************";
2261 fout << std::endl <<
" elapsed times:";
2266 fout << std::endl <<
"*************************************************************************";
2271 fout << std::endl <<
" Fval() = " << mnState.Fval() <<
", Edm() = " << mnState.Edm() <<
", NFcn() = " << mnState.NFcn();
2273 fout << std::endl <<
"*************************************************************************";
2277 Int_t maxLength = 10;
2278 for (UInt_t i=0; i<
fParams.size(); i++) {
2279 if (
fParams[i].fName.Length() > maxLength)
2280 maxLength =
fParams[i].fName.Length() + 1;
2285 fout << std::endl <<
" PARAMETERS";
2286 fout << std::endl <<
"-------------------------------------------------------------------------";
2287 fout << std::endl <<
" ";
2288 for (Int_t j=0; j<=maxLength-4; j++)
2290 fout <<
"Parabolic Minos";
2291 fout << std::endl <<
" No Name";
2292 for (Int_t j=0; j<=maxLength-4; j++)
2294 fout <<
"Value Error Negative Positive Limits";
2295 for (UInt_t i=0; i<
fParams.size(); i++) {
2297 fout.setf(std::ios::right, std::ios::adjustfield);
2299 fout << std::endl << i+1 <<
" ";
2301 fout <<
fParams[i].fName.Data();
2302 for (Int_t j=0; j<=maxLength-
fParams[i].fName.Length(); j++)
2305 fout.setf(std::ios::left, std::ios::adjustfield);
2308 fout <<
fParams[i].fValue <<
" ";
2310 fout.setf(std::ios::left, std::ios::adjustfield);
2318 if (
fParams[i].fPosErrorPresent) {
2319 fout.setf(std::ios::left, std::ios::adjustfield);
2323 fout.setf(std::ios::left, std::ios::adjustfield);
2326 fout <<
fParams[i].fPosError <<
" ";
2328 fout.setf(std::ios::left, std::ios::adjustfield);
2331 fout.setf(std::ios::left, std::ios::adjustfield);
2336 if (
fParams[i].fNoOfParams > 5) {
2337 fout.setf(std::ios::left, std::ios::adjustfield);
2339 if (
fParams[i].fLowerBoundaryPresent)
2340 fout <<
fParams[i].fLowerBoundary;
2343 fout.setf(std::ios::left, std::ios::adjustfield);
2345 if (
fParams[i].fUpperBoundaryPresent)
2346 fout <<
fParams[i].fUpperBoundary;
2350 fout.setf(std::ios::left, std::ios::adjustfield);
2353 fout.setf(std::ios::left, std::ios::adjustfield);
2359 fout << std::endl <<
"*************************************************************************";
2362 fout << std::endl <<
" COVARIANCE MATRIX";
2363 fout << std::endl <<
"-------------------------------------------------------------------------";
2364 if (mnState.HasCovariance()) {
2365 ROOT::Minuit2::MnUserCovariance cov = mnState.Covariance();
2366 fout << std::endl <<
"from " << cov.Nrow() <<
" free parameters";
2367 for (UInt_t i=0; i<cov.Nrow(); i++) {
2369 for (UInt_t j=0; j<i; j++) {
2370 fout.setf(std::ios::left, std::ios::adjustfield);
2372 if (cov(i,j) > 0.0) {
2382 fout << std::endl <<
" no covariance matrix available";
2385 fout << std::endl <<
"*************************************************************************";
2388 fout << std::endl <<
" CORRELATION COEFFICIENTS";
2389 fout << std::endl <<
"-------------------------------------------------------------------------";
2390 if (mnState.HasGlobalCC() && mnState.HasCovariance()) {
2391 ROOT::Minuit2::MnGlobalCorrelationCoeff corr = mnState.GlobalCC();
2392 ROOT::Minuit2::MnUserCovariance cov = mnState.Covariance();
2394 fout << std::endl <<
" No Global ";
2395 for (UInt_t i=0; i<
fParams.size(); i++) {
2398 fout.setf(std::ios::left, std::ios::adjustfield);
2405 if (parNo.size() != cov.Nrow()) {
2406 std::cerr << std::endl <<
"**SEVERE ERROR** in PFitter::ExecuteSave(): minuit2 and musrfit book keeping to not correspond! Unable to write correlation matrix.";
2407 std::cerr << std::endl;
2409 TString title(
"Minuit2 Output Correlation Matrix for ");
2412 title += dt.AsSQLString();
2413 std::unique_ptr<TCanvas> ccorr = std::make_unique<TCanvas>(
"ccorr",
"title", 500, 500);
2414 std::unique_ptr<TH2D> hcorr = std::make_unique<TH2D>(
"hcorr", title, cov.Nrow(), 0.0, cov.Nrow(), cov.Nrow(), 0.0, cov.Nrow());
2416 for (UInt_t i=0; i<cov.Nrow(); i++) {
2418 fout << std::endl <<
" ";
2419 fout.setf(std::ios::left, std::ios::adjustfield);
2423 fout.setf(std::ios::left, std::ios::adjustfield);
2426 fout << corr.GlobalCC()[i];
2428 for (UInt_t j=0; j<cov.Nrow(); j++) {
2429 fout.setf(std::ios::left, std::ios::adjustfield);
2434 hcorr->Fill((Double_t)i,(Double_t)i,1.0);
2438 std::cerr << std::endl <<
"**SEVERE ERROR** in PFitter::ExecuteSave(): parameter no " << parNo[i]+1 <<
" has an error == 0. Cannot correctly handle the correlation matrix.";
2439 std::cerr << std::endl;
2442 std::cerr << std::endl <<
"**SEVERE ERROR** in PFitter::ExecuteSave(): parameter no " << parNo[j]+1 <<
" has an error == 0. Cannot correctly handle the correlation matrix.";
2443 std::cerr << std::endl;
2448 hcorr->Fill((Double_t)i,(Double_t)j,dval);
2450 if (dval < 1.0e-2) {
2462 fout << dval <<
" ";
2467 TFile ff(
"MINUIT2.root",
"recreate");
2469 if (cov.Nrow() <= 6)
2470 hcorr->Draw(
"COLZTEXT");
2472 hcorr->Draw(
"COLZ");
2473 ccorr->Write(
"ccorr", TObject::kOverwrite,
sizeof(ccorr));
2474 hcorr->Write(
"hcorr", TObject::kOverwrite,
sizeof(hcorr));
2479 fout << std::endl <<
" no correlation coefficients available";
2483 fout << std::endl <<
"*************************************************************************";
2484 fout << std::endl <<
" chisq/maxLH RESULT ";
2485 fout << std::endl <<
"*************************************************************************";
2496 fitStartTime = run->at(0).GetFitRange(0);
2497 fitEndTime = run->at(0).GetFitRange(1);
2499 fout.setf(std::ios::fixed, std::ios::floatfield);
2500 fout << std::endl <<
" Time Range: " << fitStartTime <<
", " << fitEndTime << std::endl;
2502 fout.setf(std::ios::fixed, std::ios::floatfield);
2503 fout << std::endl <<
" chisq = " << std::setprecision(4) << statistics->
fMin <<
", NDF = " << statistics->
fNdf <<
", chisq/NDF = " << std::setprecision(6) << statistics->
fMin/statistics->
fNdf;
2505 fout << std::endl <<
" chisq_e = " << std::setprecision(4) << statistics->
fMinExpected <<
", NDF = " << statistics->
fNdf <<
", chisq_e/NDF = " << std::setprecision(6) << statistics->
fMinExpected/statistics->
fNdf;
2507 fout.setf(std::ios::fixed, std::ios::floatfield);
2508 fout << std::endl <<
" maxLH = " << std::setprecision(4) << statistics->
fMin <<
", NDF = " << statistics->
fNdf <<
", maxLH/NDF = " << std::setprecision(6) << statistics->
fMin/statistics->
fNdf;
2510 fout << std::endl <<
" maxLH_e = " << std::setprecision(4) << statistics->
fMinExpected <<
", NDF = " << statistics->
fNdf <<
", maxLH_e/NDF = " << std::setprecision(6) << statistics->
fMinExpected/statistics->
fNdf;
2517 fout << std::endl <<
"*************************************************************************";
2518 fout << std::endl <<
" DONE ";
2519 fout << std::endl <<
"*************************************************************************";
2520 fout << std::endl << std::endl;
2538 std::cout <<
">> PFitter::ExecuteSimplex(): will call simplex ..." << std::endl;
2546 UInt_t maxfcn = std::numeric_limits<UInt_t>::max();
2548 Double_t tolerance = 0.1;
2550 Double_t start=0.0, end=0.0;
2552 ROOT::Minuit2::FunctionMinimum min = simplex(maxfcn, tolerance);
2554 std::cout <<
">> PFitter::ExecuteMinimize(): execution time for Simplex = " << std::setprecision(3) << (end-start)/1.0e3 <<
" sec." << std::endl;
2555 TString str = TString::Format(
"Simplex: %.3f sec", (end-start)/1.0e3);
2557 if (!min.IsValid()) {
2558 std::cerr << std::endl <<
">> PFitter::ExecuteSimplex(): **WARNING**: Fit did not converge, sorry ...";
2559 std::cerr << std::endl;
2565 fFcnMin.reset(
new ROOT::Minuit2::FunctionMinimum(min));
2572 for (UInt_t i=0; i<
fParams.size(); i++) {
2573 Double_t dval = min.UserState().Value(i);
2575 Int_t m = (Int_t)(dval/360.0);
2576 dval = dval - m*360.0;
2584 Double_t minVal = min.Fval();
2585 UInt_t ndf =
fFitterFcn->GetTotalNoOfFittedBins();
2587 for (UInt_t i=0; i<
fParams.size(); i++) {
2588 if ((min.UserState().Error(i) != 0.0) && !
fMnUserParams.Parameters().at(i).IsFixed())
2599 std::cout << *
fFcnMin << std::endl;
2618 Int_t usedParams = 0;
2619 for (UInt_t i=0; i<error.size(); i++) {
2620 if (error[i] != 0.0)
2625 secFitRange.resize(1);
2628 Double_t totalExpectedChisq = 0.0;
2631 for (UInt_t k=0; k<
fSector.size(); k++) {
2633 secFitRange[0].first =
fSector[k].GetTimeRangeFirst(0);
2634 secFitRange[0].second =
fSector[k].GetTimeRangeLast();
2637 val = (*fFitterFcn)(param);
2640 ndf =
static_cast<UInt_t
>(
fFitterFcn->GetTotalNoOfFittedBins()) - usedParams;
2643 totalExpectedChisq = 0.0;
2644 fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerRun);
2645 fSector[k].SetExpectedChisq(totalExpectedChisq);
2649 fSector[k].SetChisq(chisqPerRun[i], i);
2650 fSector[k].SetExpectedChisq(expectedChisqPerRun[i], i);
2653 if (totalExpectedChisq != 0.0) {
2655 for (UInt_t i=0; i<expectedChisqPerRun.size(); i++) {
2658 fSector[k].SetNDF(ndf_run, i);
2661 }
else if (chisqPerRun.size() > 0) {
2663 for (UInt_t i=0; i<chisqPerRun.size(); i++) {
2666 fSector[k].SetNDF(ndf_run, i);
2671 chisqPerRun.clear();
2672 expectedChisqPerRun.clear();
2675 Double_t totalExpectedMaxLH = 0.0;
2678 for (UInt_t k=0; k<
fSector.size(); k++) {
2680 secFitRange[0].first =
fSector[k].GetTimeRangeFirst(0);
2681 secFitRange[0].second =
fSector[k].GetTimeRangeLast();
2684 val = (*fFitterFcn)(param);
2687 ndf =
static_cast<UInt_t
>(
fFitterFcn->GetTotalNoOfFittedBins()) - usedParams;
2690 totalExpectedMaxLH = 0.0;
2691 fFitterFcn->CalcExpectedChiSquare(param, totalExpectedMaxLH, expectedMaxLHPerRun);
2692 fSector[k].SetExpectedChisq(totalExpectedMaxLH);
2696 fSector[k].SetChisq(maxLHPerRun[i], i);
2697 fSector[k].SetExpectedChisq(expectedMaxLHPerRun[i], i);
2700 if (totalExpectedMaxLH != 0.0) {
2702 for (UInt_t i=0; i<expectedMaxLHPerRun.size(); i++) {
2705 fSector[k].SetNDF(ndf_run, i);
2711 maxLHPerRun.clear();
2712 expectedMaxLHPerRun.clear();
2729 fout << std::endl <<
"*************************************************************************";
2730 fout << std::endl <<
" SECTOR RESULTS";
2731 fout << std::endl <<
"*************************************************************************";
2734 for (UInt_t i=0; i<
fSector.size(); i++) {
2736 fout <<
" Sector " << i+1 <<
": FitRange: " <<
fSector[i].GetTimeRangeFirst(0) <<
", " <<
fSector[i].GetTimeRangeLast() << std::endl;
2737 fout <<
" chisq = " <<
fSector[i].GetChisq() <<
", NDF = " <<
fSector[i].GetNDF() <<
", chisq/NDF = " <<
fSector[i].GetChisq()/
fSector[i].GetNDF();
2739 if (
fSector[i].GetExpectedChisq() > 0) {
2740 fout << std::endl <<
" expected chisq = " <<
fSector[i].GetExpectedChisq() <<
", NDF = " <<
fSector[i].GetNDF() <<
", expected chisq/NDF = " <<
fSector[i].GetExpectedChisq()/
fSector[i].GetNDF();
2741 for (UInt_t j=0; j<
fSector[i].GetNoRuns(); j++) {
2742 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) <<
")";
2744 }
else if (
fSector[i].GetNoRuns() > 0) {
2745 for (UInt_t j=0; j<
fSector[i].GetNoRuns(); j++) {
2746 fout << std::endl <<
" run block " << j+1 <<
" (NDF/red.chisq) = (" <<
fSector[i].GetNDF(j) <<
"/" <<
fSector[i].GetChisq(j)/
fSector[i].GetNDF(j) <<
")";
2749 fout << std::endl <<
"++++";
2752 for (UInt_t i=0; i<
fSector.size(); i++) {
2754 fout <<
" Sector " << i+1 <<
": FitRange: " <<
fSector[i].GetTimeRangeFirst(0) <<
", " <<
fSector[i].GetTimeRangeLast() << std::endl;
2755 fout <<
" maxLH = " <<
fSector[i].GetChisq() <<
", NDF = " <<
fSector[i].GetNDF() <<
", maxLH/NDF = " <<
fSector[i].GetChisq()/
fSector[i].GetNDF();
2757 if (
fSector[i].GetExpectedChisq() > 0) {
2758 fout << std::endl <<
" expected maxLH = " <<
fSector[i].GetExpectedChisq() <<
", NDF = " <<
fSector[i].GetNDF() <<
", expected maxLH/NDF = " <<
fSector[i].GetExpectedChisq()/
fSector[i].GetNDF();
2759 for (UInt_t j=0; j<
fSector[i].GetNoRuns(); j++) {
2760 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) <<
")";
2762 }
else if (
fSector[i].GetNoRuns() > 0) {
2763 for (UInt_t j=0; j<
fSector[i].GetNoRuns(); j++) {
2764 fout << std::endl <<
" run block " << j+1 <<
" (NDF/red.maxLH) = (" <<
fSector[i].GetNDF(j) <<
"/" <<
fSector[i].GetChisq(j)/
fSector[i].GetNDF(j) <<
")";
2767 fout << std::endl <<
"++++";
2785 gettimeofday(&now,
nullptr);
2787 return ((Double_t)now.tv_sec * 1.0e6 + (Double_t)now.tv_usec)/1.0e3;
Bool_t fIsValid
flag. true: the fit is valid.
PFitter(PMsrHandler *runInfo, PRunListCollection *runListCollection, Bool_t chisq_only=false)
virtual void SetFitRange(const PDoublePairVector fitRange)
std::pair< Double_t, Double_t > PDoublePair
Double_t GetTimeRangeFirst(UInt_t idx)
Bool_t ExecuteSave(Bool_t first)
PDoubleVector fFirst
time stamp for fgb for a given run
std::unique_ptr< PFitterFcn > fFitterFcn
pointer to the fitter function object
PIntPairVector fCmdList
command list, first=cmd, second=cmd line index
std::vector< bool > fPhase
flag array in which an entry is true if the related parameter value is a phase
Bool_t fIsScanOnly
flag. true: scan along some parameters (no fitting).
Bool_t fChisqOnly
flag. true: calculate chi^2 only (no fitting).
std::vector< PMsrRunBlock > PMsrRunList
void PrepareSector(PDoubleVector ¶m, PDoubleVector &error)
PDoubleVector fMinPerHisto
chisq or max. likelihood per histo
PMsrLines fCmdLines
all the Minuit commands from the msr-file
#define PMN_USER_COVARIANCE
PIntVector GetParFromMap(const TString mapStr)
PStringVector fElapsedTime
PDoubleVector fChisqRun
chisq or maxLH for the sector and run
Double_t fMin
chisq or max. likelihood
void SetChisq(Double_t chisq)
Bool_t ExecuteRelease(UInt_t lineNo)
PMsrHandler * fRunInfo
pointer to the msr-file handler
PDoublePairVector fScanData
keeps the scan/contour data
Bool_t fValid
flag showing if the statistics block is valid, i.e. a fit took place which converged ...
Bool_t ExecuteFix(UInt_t lineNo)
std::vector< Int_t > PIntVector
void SetExpectedChisq(Double_t expChisq)
UInt_t fStrategy
fitting strategy (see minuit2 manual).
virtual void SetMsrStatisticMin(Double_t min)
PDoubleVector fExpectedChisqRun
expected chisq or maxLH for the sector and run
std::vector< Double_t > PDoubleVector
virtual void SetMsrStatisticNdf(UInt_t ndf)
Bool_t ExecuteSector(std::ofstream &fout)
virtual PMsrRunList * GetMsrRunList()
UInt_t fNDF
NDF for the sector.
virtual PMsrStatisticStructure * GetMsrStatistic()
Double_t fExpectedChisq
keep the expected chisq or maxLH for the sector
virtual UInt_t GetNoOfRuns()
Bool_t fUseChi2
flag. true: chi^2 fit. false: log-max-likelihood
virtual UInt_t GetNoOfParams()
ROOT::Minuit2::MnUserParameters fMnUserParams
minuit2 input parameter list
virtual Int_t ParameterInUse(UInt_t paramNo)
PMsrParamList fParams
msr-file parameters
virtual PMsrGlobalBlock * GetMsrGlobal()
std::vector< PMsrLineStructure > PMsrLines
Double_t fScanHigh
scan range high. default=0.0 which means 2 std dev. (see MnScan/MnContours in the minuit2 user manual...
virtual Double_t GetSingleRunChisq(const std::vector< Double_t > &par, const UInt_t idx) const
PDoubleVector fMinExpectedPerHisto
expected pre histo chi2 or max. likelihood
virtual Bool_t SetMsrParamStep(UInt_t i, Double_t value)
UInt_t fScanNoPoints
number of points in a scan/contour (see MnScan/MnContours in the minuit2 user manual) ...
Double_t GetExpectedChisq()
Bool_t fSectorFlag
sector command present flag
virtual Bool_t SetMsrParamValue(UInt_t i, Double_t value)
Bool_t fChisq
flag telling if min = chi2 or min = max.likelihood
#define PMN_MACHINE_PRECISION
#define PMN_USER_PARAM_STATE
virtual PMsrLines * GetMsrTheory()
virtual UInt_t GetNoOfFitParameters(UInt_t idx)
Double_t fMinExpected
expected total chi2 or max. likelihood
std::vector< PDoublePair > PDoublePairVector
UInt_t fNoOfRuns
number of runs presesent
PSectorChisq(UInt_t noOfRuns)
PRunListCollection * fRunListCollection
pointer to the run list collection
Bool_t ExecuteFitRange(UInt_t lineNo)
PUIntVector fNDFRun
NDF for the sector and run.
void SetRunFirstTime(Double_t first, UInt_t idx)
UInt_t fScanParameter[2]
scan parameter. idx=0: used for scan and contour, idx=1: used for contour (see MnScan/MnContours in t...
virtual PMsrParamList * GetMsrParamList()
Double_t fScanLow
scan range low. default=0.0 which means 2 std dev. (see MnScan/MnContours in the minuit2 user manual)...
Double_t fChisq
chisq or maxLH for the sector
Bool_t fScanAll
flag. false: single parameter scan, true: not implemented yet (see MnScan/MnContours in the minuit2 u...
UInt_t fNdf
number of degrees of freedom
PIntVector GetParFromFun(const TString funStr)
virtual PMsrLines * GetMsrFunctions()
UInt_t fPrintLevel
tag, showing the level of messages whished. 0=minimum, 1=standard, 2=maximum
virtual PMsrLines * GetMsrCommands()
virtual Bool_t SetMsrParamPosErrorPresent(UInt_t i, Bool_t value)
Bool_t fConverged
flag. true: the fit has converged.
virtual Double_t GetFitRange(UInt_t idx)
Double_t fLast
requested time stamp
PUIntVector fNdfPerHisto
number of degrees of freedom per histo
Bool_t ExecutePrintLevel(UInt_t lineNo)
virtual Bool_t SetMsrParamPosError(UInt_t i, Double_t value)
PDoublePairVector fOriginalFitRange
keeps the original fit range in case there is a range command in the COMMAND block ...
virtual Double_t GetSingleRunMaximumLikelihood(const std::vector< Double_t > &par, const UInt_t idx) const
std::pair< Int_t, Int_t > PIntPair
std::vector< PSectorChisq > fSector
stores all chisq/maxLH sector information
virtual const TString & GetFileName() const
std::unique_ptr< ROOT::Minuit2::FunctionMinimum > fFcnMin
function minimum object