diff --git a/ChangeLog b/ChangeLog index b35eaeaa..36b85967 100644 --- a/ChangeLog +++ b/ChangeLog @@ -6,6 +6,7 @@ changes since 0.9.0 =================================== +NEW added chisq per run block information (all fit types). The same output rules as for the expected chisq are in place. NEW calculate expected chisq (Pearson's chisq) for single histogram fits. It always will send this information to the stdout but only writes it to the msr-file if a corresponding flag in the musrfit_startup.xml is enabling it. The following expected chisq's are calculated: total expected chisq, and the per histo expected chisq. (MUSR-194) diff --git a/src/classes/PFitter.cpp b/src/classes/PFitter.cpp index eb011658..008fd1cb 100644 --- a/src/classes/PFitter.cpp +++ b/src/classes/PFitter.cpp @@ -174,6 +174,11 @@ Bool_t PFitter::DoFit() Double_t totalExpectedChisq = 0.0; std::vector expectedChisqPerHisto; fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto); + // calculate chisq per run + std::vector chisqPerHisto; + for (UInt_t i=0; iGetMsrRunList()->size(); i++) { + chisqPerHisto.push_back(fRunListCollection->GetSingleRunChisq(param, i)); + } cout << endl << endl << ">> chisq = " << val << ", NDF = " << ndf << ", chisq/NDF = " << val/ndf; @@ -183,14 +188,20 @@ Bool_t PFitter::DoFit() for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); if (ndf_histo > 0) - cout << endl << ">> run block " << i+1 << ": expected chisq = " << expectedChisqPerHisto[i] << ", ndf = " << ndf_histo << ", expected chisq/NDF = " << expectedChisqPerHisto[i]/ndf_histo; - else - cout << endl << ">> run block " << i+1 << ": expected chisq = " << expectedChisqPerHisto[i]; + cout << endl << ">> run block " << i+1 << ": (NDF/red.chisq/red.chisq_e) = (" << ndf_histo << "/" << chisqPerHisto[i]/ndf_histo << "/" << expectedChisqPerHisto[i]/ndf_histo << ")"; + } + } else if (chisqPerHisto.size() > 0) { // in case expected chisq is not applicable like for asymmetry fits + UInt_t ndf_histo = 0; + for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_histo > 0) + cout << endl << ">> run block " << i+1 << ": (NDF/red.chisq) = (" << ndf_histo << "/" << chisqPerHisto[i]/ndf_histo << ")"; } - - // clean up - expectedChisqPerHisto.clear(); } + + // clean up + chisqPerHisto.clear(); + expectedChisqPerHisto.clear(); } else { // max. log likelihood cout << endl << endl << ">> maxLH = " << val << ", NDF = " << ndf << ", maxLH/NDF = " << val/ndf; } @@ -1456,6 +1467,12 @@ Bool_t PFitter::ExecuteSave() fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto); + // calculate chisq per run + std::vector chisqPerHisto; + for (UInt_t i=0; iGetMsrRunList()->size(); i++) { + chisqPerHisto.push_back(fRunListCollection->GetSingleRunChisq(param, i)); + } + if (totalExpectedChisq != 0.0) { // i.e. applicable for single histogram fits only // get the ndf's of the histos UInt_t ndf_histo; @@ -1467,11 +1484,31 @@ Bool_t PFitter::ExecuteSave() // feed the msr-file handler PMsrStatisticStructure *statistics = fRunInfo->GetMsrStatistic(); if (statistics) { + statistics->fMinPerHisto = chisqPerHisto; statistics->fMinExpected = totalExpectedChisq; statistics->fMinExpectedPerHisto = expectedChisqPerHisto; statistics->fNdfPerHisto = ndfPerHisto; } + } else if (chisqPerHisto.size() > 1) { // in case expected chisq is not applicable like for asymmetry fits + UInt_t ndf_histo = 0; + for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + ndfPerHisto.push_back(ndf_histo); + } + + // feed the msr-file handler + PMsrStatisticStructure *statistics = fRunInfo->GetMsrStatistic(); + if (statistics) { + statistics->fMinPerHisto = chisqPerHisto; + statistics->fNdfPerHisto = ndfPerHisto; + } } + + // clean up + param.clear(); + expectedChisqPerHisto.clear(); + ndfPerHisto.clear(); + chisqPerHisto.clear(); } cout << ">> PFitter::ExecuteSave(): will write minuit2 output file ..." << endl; diff --git a/src/classes/PFitterFcn.cpp b/src/classes/PFitterFcn.cpp index 7c1203f6..187a63d4 100644 --- a/src/classes/PFitterFcn.cpp +++ b/src/classes/PFitterFcn.cpp @@ -92,30 +92,11 @@ Double_t PFitterFcn::operator()(const std::vector& par) const return value; } - -//-------------------------------------------------------------------------- -// GetNoOfFittedBins() -//-------------------------------------------------------------------------- -/** - *

Get the number of fitted bins of the run block idx. - * - * \param idx index of the run block - */ -UInt_t PFitterFcn::GetNoOfFittedBins(const UInt_t idx) -{ - UInt_t result = 0; - - if (idx < fRunListCollection->GetNoOfSingleHisto()) - result = fRunListCollection->GetNoOfBinsFitted(idx); - - return result; -} - //-------------------------------------------------------------------------- // CalcExpectedChiSquare() //-------------------------------------------------------------------------- /** - *

Calculates the expected chisq, if applicable. + *

Calculates the expected chisq, expected chisq per run, and chisq per run, if applicable. * * \param par * \param totalExpectedChisq expected chisq for all run blocks @@ -130,6 +111,8 @@ void PFitterFcn::CalcExpectedChiSquare(const std::vector &par, Double_ // only do something for chisq if (fUseChi2) { Double_t value = 0.0; + + // single histo for (UInt_t i=0; iGetNoOfSingleHisto(); i++) { value = fRunListCollection->GetSingleHistoChisqExpected(par, i); // calculate the expected chisq for single histo run block 'i' expectedChisqPerRun.push_back(value); diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index 7a013f2d..2ee74488 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -1019,15 +1019,8 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) for (UInt_t i=0; i 0) { - str.Form(" run block %d: expected chisq=%.1lf, NDF=%d, expected chisq/NDF=%lf", - i+1, fStatistic.fMinExpectedPerHisto[i], fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i] /fStatistic.fNdfPerHisto[i]); - if (fWriteExpectedChisq) - fout << str.Data() << endl; - - if (messages) - cout << str.Data() << endl; - } else { - str.Form(" run block %d: expected chisq=%.1lf", i+1, fStatistic.fMinExpectedPerHisto[i]); + str.Form(" run block %d: (NDF/red.chisq/red.chisq_e) = (%d/%lf/%lf)", + i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i]/fStatistic.fNdfPerHisto[i]); if (fWriteExpectedChisq) fout << str.Data() << endl; @@ -1035,6 +1028,16 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) cout << str.Data() << endl; } } + } else if (fStatistic.fNdfPerHisto.size() > 1) { // check if per run chisq needs to be written + for (UInt_t i=0; i 0) { - str.Form(" run block %d: expected chisq=%.1lf, NDF=%d, expected chisq/NDF=%lf", - i+1, fStatistic.fMinExpectedPerHisto[i], fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i] /fStatistic.fNdfPerHisto[i]); - if (fWriteExpectedChisq) - fout << str.Data() << endl; - - if (messages) - cout << str.Data() << endl; - } else { - str.Form(" run block %d: expected chisq=%.1lf", i+1, fStatistic.fMinExpectedPerHisto[i]); + str.Form(" run block %d: (NDF/red.chisq/red.chisq_e) = (%d/%lf/%lf)", + i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i]/fStatistic.fNdfPerHisto[i]); if (fWriteExpectedChisq) fout << str.Data() << endl; @@ -1083,6 +1079,16 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) cout << str.Data() << endl; } } + } else if (fStatistic.fNdfPerHisto.size() > 1) { // check if per run chisq needs to be written + for (UInt_t i=0; i 0) { - str.Form(" run block %d: expected chisq=%.1lf, NDF=%d, expected chisq/NDF=%lf", - i+1, fStatistic.fMinExpectedPerHisto[i], fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i] /fStatistic.fNdfPerHisto[i]); - if (fWriteExpectedChisq) - fout << str.Data() << endl; - - if (messages) - cout << str.Data() << endl; - } else { - str.Form(" run block %d: expected chisq=%.1lf", i+1, fStatistic.fMinExpectedPerHisto[i]); + str.Form(" run block %d: (NDF/red.chisq/red.chisq_e) = (%d/%lf/%lf)", + i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinPerHisto[i]/fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i]/fStatistic.fNdfPerHisto[i]); if (fWriteExpectedChisq) fout << str.Data() << endl; @@ -1153,6 +1152,16 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) cout << str.Data() << endl; } } + } else if (fStatistic.fNdfPerHisto.size() > 1) { // check if per run chisq needs to be written + for (UInt_t i=0; i 0) { - str.Form(" run block %d: expected chisq=%.1lf, NDF=%d, expected chisq/NDF=%lf", - i+1, fStatistic.fMinExpectedPerHisto[i], fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i] /fStatistic.fNdfPerHisto[i]); - if (fWriteExpectedChisq) - fout << str.Data() << endl; - - if (messages) - cout << str.Data() << endl; - } else { - str.Form(" run block %d: expected chisq=%.1lf", i+1, fStatistic.fMinExpectedPerHisto[i]); + str.Form(" run block %d: (NDF/red.chisq/red.chisq_e) =(%d/%lf/%lf)", + i+1, fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i]/fStatistic.fNdfPerHisto[i], fStatistic.fMinExpectedPerHisto[i]/fStatistic.fNdfPerHisto[i]); if (fWriteExpectedChisq) fout << str.Data() << endl; @@ -1207,6 +1209,16 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) cout << str.Data() << endl; } } + } else if (fStatistic.fNdfPerHisto.size() > 1) { // check if per run chisq needs to be written + for (UInt_t i=0; iCalculate the number of fit parameters for single histo fits. + *

Calculate the number of fit parameters. * * \param idx run block index */ @@ -3993,7 +4005,7 @@ UInt_t PMsrHandler::GetNoOfFitParameters(UInt_t idx) return 0; } - // get N0 parameter, possible parameter number or function + // get N0 parameter, possible parameter number or function (single histo fit) if (fRuns[idx].GetNormParamNo() != -1) { if (fRuns[idx].GetNormParamNo() < MSR_PARAM_FUN_OFFSET) // parameter paramVector.push_back(fRuns[idx].GetNormParamNo()); @@ -4001,10 +4013,18 @@ UInt_t PMsrHandler::GetNoOfFitParameters(UInt_t idx) funVector.push_back(fRuns[idx].GetNormParamNo() - MSR_PARAM_FUN_OFFSET); } - // get background parameter, for the case the background is fitted. + // get background parameter, for the case the background is fitted (single histo fit) if (fRuns[idx].GetBkgFitParamNo() != -1) paramVector.push_back(fRuns[idx].GetBkgFitParamNo()); + // get alpha parameter if present (asymmetry fit) + if (fRuns[idx].GetAlphaParamNo() != -1) + paramVector.push_back(fRuns[idx].GetAlphaParamNo()); + + // get beta parameter if present (asymmetry fit) + if (fRuns[idx].GetBetaParamNo() != -1) + paramVector.push_back(fRuns[idx].GetBetaParamNo()); + // go through the theory block and collect parameters // possible entries: number -> parameter, fun -> function, map -> maps for (UInt_t i=0; i& par) return chisq; } +//-------------------------------------------------------------------------- +// CalcChiSquareExpected (public) +//-------------------------------------------------------------------------- +/** + *

Calculate expected chi-square. Currently not implemented since not clear what to be done. + * + * return: + * - chisq value == 0.0 + * + * \param par parameter vector iterated by minuit2 + */ +Double_t PRunAsymmetry::CalcChiSquareExpected(const std::vector& par) +{ + return 0.0; +} + //-------------------------------------------------------------------------- // CalcMaxLikelihood //-------------------------------------------------------------------------- @@ -244,13 +260,6 @@ Double_t PRunAsymmetry::CalcMaxLikelihood(const std::vector& par) */ UInt_t PRunAsymmetry::GetNoOfFitBins() { -// Double_t time; -// fNoOfFitBins=0; -// for (UInt_t i=0; isize(); i++) { -// time = fData.GetDataTimeStart() + (Double_t)i * fData.GetDataTimeStep(); -// if ((time >= fFitStartTime) && (time <= fFitEndTime)) -// fNoOfFitBins++; -// } CalcNoOfFitBins(); return fNoOfFitBins; diff --git a/src/classes/PRunListCollection.cpp b/src/classes/PRunListCollection.cpp index b2467780..29200912 100644 --- a/src/classes/PRunListCollection.cpp +++ b/src/classes/PRunListCollection.cpp @@ -80,7 +80,7 @@ PRunListCollection::~PRunListCollection() } //-------------------------------------------------------------------------- -// Add +// Add (public) //-------------------------------------------------------------------------- /** *

Adds a processed set of data to the handler. @@ -150,7 +150,7 @@ void PRunListCollection::SetFitRange(const PDoublePairVector fitRange) } //-------------------------------------------------------------------------- -// GetSingleHistoChisq +// GetSingleHistoChisq (public) //-------------------------------------------------------------------------- /** *

Calculates chi-square of all single histogram runs of a msr-file. @@ -171,29 +171,7 @@ Double_t PRunListCollection::GetSingleHistoChisq(const std::vector& pa } //-------------------------------------------------------------------------- -// GetSingleHistoChisqExpected -//-------------------------------------------------------------------------- -/** - *

Calculates expected chi-square of the single histogram with run block index idx of a msr-file. - * - * return: - * - expected chi-square of for a single histogram - * - * \param par fit parameter vector - * \param idx run block index - */ -Double_t PRunListCollection::GetSingleHistoChisqExpected(const std::vector& par, const UInt_t idx) const -{ - Double_t expectedChisq = 0.0; - - if (idx < fRunSingleHistoList.size()) - expectedChisq = fRunSingleHistoList[idx]->CalcChiSquareExpected(par); - - return expectedChisq; -} - -//-------------------------------------------------------------------------- -// GetAsymmetryChisq +// GetAsymmetryChisq (public) //-------------------------------------------------------------------------- /** *

Calculates chi-square of all asymmetry runs of a msr-file. @@ -214,7 +192,7 @@ Double_t PRunListCollection::GetAsymmetryChisq(const std::vector& par) } //-------------------------------------------------------------------------- -// GetMuMinusChisq +// GetMuMinusChisq (public) //-------------------------------------------------------------------------- /** *

Calculates chi-square of all mu minus runs of a msr-file. @@ -235,7 +213,7 @@ Double_t PRunListCollection::GetMuMinusChisq(const std::vector& par) c } //-------------------------------------------------------------------------- -// GetNonMusrChisq +// GetNonMusrChisq (public) //-------------------------------------------------------------------------- /** *

Calculates chi-square of all non-muSR runs of a msr-file. @@ -256,7 +234,109 @@ Double_t PRunListCollection::GetNonMusrChisq(const std::vector& par) c } //-------------------------------------------------------------------------- -// GetSingleHistoMaximumLikelihood +// GetSingleHistoChisqExpected (public) +//-------------------------------------------------------------------------- +/** + *

Calculates expected chi-square of the single histogram with run block index idx of a msr-file. + * + * return: + * - expected chi-square of for a single histogram + * + * \param par fit parameter vector + * \param idx run block index + */ +Double_t PRunListCollection::GetSingleHistoChisqExpected(const std::vector& par, const UInt_t idx) const +{ + Double_t expectedChisq = 0.0; + + if (idx > fMsrInfo->GetMsrRunList()->size()) { + cerr << ">> PRunListCollection::GetSingleHistoChisqExpected() **ERROR** idx=" << idx << " is out of range [0.." << fMsrInfo->GetMsrRunList()->size() << "[" << endl << endl; + return expectedChisq; + } + + Int_t type = fMsrInfo->GetMsrRunList()->at(idx).GetFitType(); + + // count how many entries of this fit-type are present up to idx + UInt_t subIdx = 0; + for (UInt_t i=0; iGetMsrRunList()->at(i).GetFitType() == type) + subIdx++; + } + + // return the chisq of the single run + switch (type) { + case PRUN_SINGLE_HISTO: + expectedChisq = fRunSingleHistoList[subIdx]->CalcChiSquareExpected(par); + break; + case PRUN_ASYMMETRY: + expectedChisq = fRunAsymmetryList[subIdx]->CalcChiSquareExpected(par); + break; + case PRUN_MU_MINUS: + expectedChisq = fRunMuMinusList[subIdx]->CalcChiSquareExpected(par); + break; + case PRUN_NON_MUSR: + expectedChisq = fRunNonMusrList[subIdx]->CalcChiSquareExpected(par); + break; + default: + break; + } + + return expectedChisq; +} + +//-------------------------------------------------------------------------- +// GetSingleRunChisq (public) +//-------------------------------------------------------------------------- +/** + *

Calculates chi-square of a single run-block entry of the msr-file. + * + * return: + * - chi-square of single run-block entry with index idx + * + * \param par fit parameter vector + * \param idx run block index + */ +Double_t PRunListCollection::GetSingleRunChisq(const std::vector& par, const UInt_t idx) const +{ + Double_t chisq = 0.0; + + if (idx > fMsrInfo->GetMsrRunList()->size()) { + cerr << ">> PRunListCollection::GetSingleRunChisq() **ERROR** idx=" << idx << " is out of range [0.." << fMsrInfo->GetMsrRunList()->size() << "[" << endl << endl; + return chisq; + } + + Int_t type = fMsrInfo->GetMsrRunList()->at(idx).GetFitType(); + + // count how many entries of this fit-type are present up to idx + UInt_t subIdx = 0; + for (UInt_t i=0; iGetMsrRunList()->at(i).GetFitType() == type) + subIdx++; + } + + // return the chisq of the single run + switch (type) { + case PRUN_SINGLE_HISTO: + chisq = fRunSingleHistoList[subIdx]->CalcChiSquare(par); + break; + case PRUN_ASYMMETRY: + chisq = fRunAsymmetryList[subIdx]->CalcChiSquare(par); + break; + case PRUN_MU_MINUS: + chisq = fRunMuMinusList[subIdx]->CalcChiSquare(par); + break; + case PRUN_NON_MUSR: + chisq = fRunNonMusrList[subIdx]->CalcChiSquare(par); + break; + default: + break; + } + + return chisq; +} + +//-------------------------------------------------------------------------- +// GetSingleHistoMaximumLikelihood (public) //-------------------------------------------------------------------------- /** *

Calculates log max-likelihood of all single histogram runs of a msr-file. @@ -277,7 +357,7 @@ Double_t PRunListCollection::GetSingleHistoMaximumLikelihood(const std::vector Since it is not clear yet how to handle asymmetry fits with max likelihood @@ -299,7 +379,7 @@ Double_t PRunListCollection::GetAsymmetryMaximumLikelihood(const std::vectorCalculates log max-likelihood of all mu minus runs of a msr-file. @@ -320,7 +400,7 @@ Double_t PRunListCollection::GetMuMinusMaximumLikelihood(const std::vector Since it is not clear yet how to handle non musr fits with max likelihood @@ -342,7 +422,7 @@ Double_t PRunListCollection::GetNonMusrMaximumLikelihood(const std::vectorNumber of bins in run block idx to be fitted. Only used for single histogram @@ -357,15 +437,44 @@ UInt_t PRunListCollection::GetNoOfBinsFitted(const UInt_t idx) const { UInt_t result = 0; - if (idx < fRunSingleHistoList.size()) - result = fRunSingleHistoList[idx]->GetNoOfFitBins(); + if (idx > fMsrInfo->GetMsrRunList()->size()) { + cerr << ">> PRunListCollection::GetNoOfBinsFitted() **ERROR** idx=" << idx << " is out of range [0.." << fMsrInfo->GetMsrRunList()->size() << "[" << endl << endl; + return result; + } + + Int_t type = fMsrInfo->GetMsrRunList()->at(idx).GetFitType(); + + // count how many entries of this fit-type are present up to idx + UInt_t subIdx = 0; + for (UInt_t i=0; iGetMsrRunList()->at(i).GetFitType() == type) + subIdx++; + } + + // return the chisq of the single run + switch (type) { + case PRUN_SINGLE_HISTO: + result = fRunSingleHistoList[subIdx]->GetNoOfFitBins(); + break; + case PRUN_ASYMMETRY: + result = fRunAsymmetryList[subIdx]->GetNoOfFitBins(); + break; + case PRUN_MU_MINUS: + result = fRunMuMinusList[subIdx]->GetNoOfFitBins(); + break; + case PRUN_NON_MUSR: + result = fRunNonMusrList[subIdx]->GetNoOfFitBins(); + break; + default: + break; + } return result; } //-------------------------------------------------------------------------- -// GetTotalNoOfBinsFitted +// GetTotalNoOfBinsFitted (public) //-------------------------------------------------------------------------- /** *

Counts the total number of bins to be fitted. @@ -393,7 +502,7 @@ UInt_t PRunListCollection::GetTotalNoOfBinsFitted() const } //-------------------------------------------------------------------------- -// GetSingleHisto +// GetSingleHisto (public) //-------------------------------------------------------------------------- /** *

Get a processed single histogram data set. @@ -436,7 +545,7 @@ PRunData* PRunListCollection::GetSingleHisto(UInt_t index, EDataSwitch tag) } //-------------------------------------------------------------------------- -// GetAsymmetry +// GetAsymmetry (public) //-------------------------------------------------------------------------- /** *

Get a processed asymmetry data set. @@ -479,7 +588,7 @@ PRunData* PRunListCollection::GetAsymmetry(UInt_t index, EDataSwitch tag) } //-------------------------------------------------------------------------- -// GetMuMinus +// GetMuMinus (public) //-------------------------------------------------------------------------- /** *

Get a processed mu minus data set. @@ -519,7 +628,7 @@ PRunData* PRunListCollection::GetMuMinus(UInt_t index, EDataSwitch tag) } //-------------------------------------------------------------------------- -// GetNonMusr +// GetNonMusr (public) //-------------------------------------------------------------------------- /** *

Get a processed non-muSR data set. @@ -559,7 +668,7 @@ PRunData* PRunListCollection::GetNonMusr(UInt_t index, EDataSwitch tag) } //-------------------------------------------------------------------------- -// GetTemp +// GetTemp (public) //-------------------------------------------------------------------------- /** *

Get the temperature from the data set. @@ -575,7 +684,7 @@ const PDoublePairVector* PRunListCollection::GetTemp(const TString &runName) con } //-------------------------------------------------------------------------- -// GetField +// GetField (public) //-------------------------------------------------------------------------- /** *

Get the magnetic field from the data set. @@ -591,7 +700,7 @@ Double_t PRunListCollection::GetField(const TString &runName) const } //-------------------------------------------------------------------------- -// GetEnergy +// GetEnergy (public) //-------------------------------------------------------------------------- /** *

Get the muon implantation energy from the data set. @@ -607,7 +716,7 @@ Double_t PRunListCollection::GetEnergy(const TString &runName) const } //-------------------------------------------------------------------------- -// GetSetup +// GetSetup (public) //-------------------------------------------------------------------------- /** *

Get the setup information from the data set. @@ -623,7 +732,7 @@ const Char_t* PRunListCollection::GetSetup(const TString &runName) const } //-------------------------------------------------------------------------- -// GetXAxisTitle +// GetXAxisTitle (public) //-------------------------------------------------------------------------- /** *

Get the x-axis title (used with non-muSR fit). @@ -656,7 +765,7 @@ const Char_t* PRunListCollection::GetXAxisTitle(const TString &runName, const UI } //-------------------------------------------------------------------------- -// GetYAxisTitle +// GetYAxisTitle (public) //-------------------------------------------------------------------------- /** *

Get the y-axis title (used with non-muSR fit). diff --git a/src/classes/PRunMuMinus.cpp b/src/classes/PRunMuMinus.cpp index 15df0dfe..74a4003a 100644 --- a/src/classes/PRunMuMinus.cpp +++ b/src/classes/PRunMuMinus.cpp @@ -96,6 +96,22 @@ Double_t PRunMuMinus::CalcChiSquare(const std::vector& par) return chisq; } +//-------------------------------------------------------------------------- +// CalcChiSquareExpected (public) +//-------------------------------------------------------------------------- +/** + *

Calculate expected chi-square. Currently not implemented. + * + * return: + * - chisq value == 0.0 + * + * \param par parameter vector iterated by minuit2 + */ +Double_t PRunMuMinus::CalcChiSquareExpected(const std::vector& par) +{ + return 0.0; +} + //-------------------------------------------------------------------------- // CalcMaxLikelihood //-------------------------------------------------------------------------- diff --git a/src/classes/PRunNonMusr.cpp b/src/classes/PRunNonMusr.cpp index 6f5dd875..fe0e1155 100644 --- a/src/classes/PRunNonMusr.cpp +++ b/src/classes/PRunNonMusr.cpp @@ -118,6 +118,22 @@ Double_t PRunNonMusr::CalcChiSquare(const std::vector& par) return chisq; } +//-------------------------------------------------------------------------- +// CalcChiSquareExpected (public) +//-------------------------------------------------------------------------- +/** + *

Calculate expected chi-square. Currently not implemented since not clear what to be done. + * + * return: + * - chisq value == 0.0 + * + * \param par parameter vector iterated by minuit2 + */ +Double_t PRunNonMusr::CalcChiSquareExpected(const std::vector& par) +{ + return 0.0; +} + //-------------------------------------------------------------------------- // CalcMaxLikelihood //-------------------------------------------------------------------------- diff --git a/src/classes/PStartupHandler.cpp b/src/classes/PStartupHandler.cpp index 6416205f..0b3052d9 100644 --- a/src/classes/PStartupHandler.cpp +++ b/src/classes/PStartupHandler.cpp @@ -166,7 +166,7 @@ void PStartupHandler::OnStartDocument() fFourierDefaults.fPlotRange[1] = -1.0; fFourierDefaults.fPhaseIncrement = 1.0; - fWriteExpectedChisq = false; + fWritePerRunBlockChisq = false; } //-------------------------------------------------------------------------- @@ -195,8 +195,8 @@ void PStartupHandler::OnStartElement(const Char_t *str, const TList *attributes) { if (!strcmp(str, "data_path")) { fKey = eDataPath; - } else if (!strcmp(str, "write_expected_chisq")) { - fKey = eWriteExpectedChisq; + } else if (!strcmp(str, "write_per_run_block_chisq")) { + fKey = eWritePerRunBlockChisq; } else if (!strcmp(str, "marker")) { fKey = eMarker; } else if (!strcmp(str, "color")) { @@ -251,10 +251,10 @@ void PStartupHandler::OnCharacters(const Char_t *str) // add str to the path list fDataPathList.push_back(str); break; - case eWriteExpectedChisq: + case eWritePerRunBlockChisq: tstr = TString(str); if (tstr.BeginsWith("y") || tstr.BeginsWith("Y")) - fWriteExpectedChisq = true; + fWritePerRunBlockChisq = true; break; case eMarker: // check that str is a number diff --git a/src/include/PFitterFcn.h b/src/include/PFitterFcn.h index 0e9215dc..59135144 100644 --- a/src/include/PFitterFcn.h +++ b/src/include/PFitterFcn.h @@ -51,7 +51,7 @@ class PFitterFcn : public ROOT::Minuit2::FCNBase Double_t operator()(const std::vector &par) const; UInt_t GetTotalNoOfFittedBins() { return fRunListCollection->GetTotalNoOfBinsFitted(); } - UInt_t GetNoOfFittedBins(const UInt_t idx); + UInt_t GetNoOfFittedBins(const UInt_t idx) { return fRunListCollection->GetNoOfBinsFitted(idx); } void CalcExpectedChiSquare(const std::vector &par, Double_t &totalExpectedChisq, std::vector &expectedChisqPerRun); private: diff --git a/src/include/PMusr.h b/src/include/PMusr.h index 7fd56a6d..9dfe21db 100644 --- a/src/include/PMusr.h +++ b/src/include/PMusr.h @@ -583,7 +583,8 @@ typedef struct { PMsrLines fStatLines; ///< statistics block in msr-file clear text TString fDate; ///< string holding fitting date and time Bool_t fChisq; ///< flag telling if min = chi2 or min = max.likelihood - Double_t fMin; ///< chi2 or max. likelihood + Double_t fMin; ///< chisq or max. likelihood + PDoubleVector fMinPerHisto; ///< chisq or max. likelihood per histo UInt_t fNdf; ///< number of degrees of freedom Double_t fMinExpected; ///< expected total chi2 or max. likelihood PDoubleVector fMinExpectedPerHisto; ///< expected pre histo chi2 or max. likelihood diff --git a/src/include/PRunAsymmetry.h b/src/include/PRunAsymmetry.h index bb3ce942..e6c2741f 100644 --- a/src/include/PRunAsymmetry.h +++ b/src/include/PRunAsymmetry.h @@ -46,6 +46,7 @@ class PRunAsymmetry : public PRunBase virtual ~PRunAsymmetry(); virtual Double_t CalcChiSquare(const std::vector& par); + virtual Double_t CalcChiSquareExpected(const std::vector& par); virtual Double_t CalcMaxLikelihood(const std::vector& par); virtual void CalcTheory(); diff --git a/src/include/PRunListCollection.h b/src/include/PRunListCollection.h index ce93454f..05ee6cf5 100644 --- a/src/include/PRunListCollection.h +++ b/src/include/PRunListCollection.h @@ -59,11 +59,13 @@ class PRunListCollection virtual void SetFitRange(const PDoublePairVector fitRange); virtual Double_t GetSingleHistoChisq(const std::vector& par) const; - virtual Double_t GetSingleHistoChisqExpected(const std::vector& par, const UInt_t idx) const; virtual Double_t GetAsymmetryChisq(const std::vector& par) const; virtual Double_t GetMuMinusChisq(const std::vector& par) const; virtual Double_t GetNonMusrChisq(const std::vector& par) const; + virtual Double_t GetSingleHistoChisqExpected(const std::vector& par, const UInt_t idx) const; + virtual Double_t GetSingleRunChisq(const std::vector& par, const UInt_t idx) const; + virtual Double_t GetSingleHistoMaximumLikelihood(const std::vector& par) const; virtual Double_t GetAsymmetryMaximumLikelihood(const std::vector& par) const; virtual Double_t GetMuMinusMaximumLikelihood(const std::vector& par) const; diff --git a/src/include/PRunMuMinus.h b/src/include/PRunMuMinus.h index 4adabfcc..6e316b33 100644 --- a/src/include/PRunMuMinus.h +++ b/src/include/PRunMuMinus.h @@ -45,6 +45,7 @@ class PRunMuMinus : public PRunBase virtual ~PRunMuMinus(); virtual Double_t CalcChiSquare(const std::vector& par); + virtual Double_t CalcChiSquareExpected(const std::vector& par); virtual Double_t CalcMaxLikelihood(const std::vector& par); virtual void CalcTheory(); diff --git a/src/include/PRunNonMusr.h b/src/include/PRunNonMusr.h index 35a6cd76..4f33fff6 100644 --- a/src/include/PRunNonMusr.h +++ b/src/include/PRunNonMusr.h @@ -46,6 +46,7 @@ class PRunNonMusr : public PRunBase virtual ~PRunNonMusr(); virtual Double_t CalcChiSquare(const std::vector& par); + virtual Double_t CalcChiSquareExpected(const std::vector& par); virtual Double_t CalcMaxLikelihood(const std::vector& par); virtual void CalcTheory(); diff --git a/src/include/PStartupHandler.h b/src/include/PStartupHandler.h index 17abbce9..dedde505 100644 --- a/src/include/PStartupHandler.h +++ b/src/include/PStartupHandler.h @@ -75,14 +75,14 @@ class PStartupHandler : public TObject, public TQObject virtual void CheckLists(); - virtual const Bool_t GetWriteExpectedChisq() { return fWriteExpectedChisq; } ///< returns the write_expected_chisq flag + virtual const Bool_t GetWritePerRunBlockChisq() { return fWritePerRunBlockChisq; } ///< returns the write_per_run_block_chisq flag virtual PMsrFourierStructure GetFourierDefaults() { return fFourierDefaults; } ///< returns the Fourier defaults virtual const PStringVector GetDataPathList() const { return fDataPathList; } ///< returns the search data path list virtual const PIntVector GetMarkerList() const { return fMarkerList; } ///< returns the marker list virtual const PIntVector GetColorList() const { return fColorList; } ///< returns the color list private: - enum EKeyWords {eEmpty, eComment, eDataPath, eWriteExpectedChisq, + enum EKeyWords {eEmpty, eComment, eDataPath, eWritePerRunBlockChisq, eFourierSettings, eUnits, eFourierPower, eApodization, ePlot, ePhase, ePhaseIncrement, eRootSettings, eMarkerList, eMarker, eColorList, eColor}; @@ -90,7 +90,7 @@ class PStartupHandler : public TObject, public TQObject Bool_t fStartupFileFound; ///< startup file found flag TString fStartupFilePath; ///< full musrfit_startup.xml startup file paths - Bool_t fWriteExpectedChisq; ///< flag showing if the expected chisq shall be written to the msr-file + Bool_t fWritePerRunBlockChisq; ///< flag showing if per run block chisq and the expected chisq shall be written to the msr-file PMsrFourierStructure fFourierDefaults; ///< Fourier defaults PStringVector fDataPathList; ///< search data path list PIntVector fMarkerList; ///< marker list diff --git a/src/musrfit.cpp b/src/musrfit.cpp index f3810eb5..c17b584a 100644 --- a/src/musrfit.cpp +++ b/src/musrfit.cpp @@ -463,7 +463,7 @@ int main(int argc, char *argv[]) } // read msr-file - PMsrHandler *msrHandler = new PMsrHandler(filename, startupHandler->GetWriteExpectedChisq()); + PMsrHandler *msrHandler = new PMsrHandler(filename, startupHandler->GetWritePerRunBlockChisq()); status = msrHandler->ReadMsrFile(); if (status != PMUSR_SUCCESS) { switch (status) { diff --git a/src/musrfit_startup.xml b/src/musrfit_startup.xml index 8e7e1492..cd27b09a 100644 --- a/src/musrfit_startup.xml +++ b/src/musrfit_startup.xml @@ -13,7 +13,7 @@ /afs/psi.ch/project/bulkmusr/data/gpd /afs/psi.ch/project/bulkmusr/data/ltf /afs/psi.ch/project/bulkmusr/data/alc - n + n Gauss 0