diff --git a/ChangeLog b/ChangeLog index fba002e7..bec8ba22 100644 --- a/ChangeLog +++ b/ChangeLog @@ -6,6 +6,9 @@ changes since 0.9.0 =================================== +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) NEW the chi^2 calculation in single-histogram and asymmetry fits is parallelized if musrfit is built using a compiler supporting OpenMP (e.g. GCC >= 4.2) Using --disable-omp this feature can be disabled on the configure level. diff --git a/src/classes/PFitter.cpp b/src/classes/PFitter.cpp index 27cea37e..eb011658 100644 --- a/src/classes/PFitter.cpp +++ b/src/classes/PFitter.cpp @@ -170,7 +170,27 @@ Bool_t PFitter::DoFit() Int_t ndf = fFitterFcn->GetTotalNoOfFittedBins() - usedParams; Double_t val = (*fFitterFcn)(param); if (fUseChi2) { + // calculate expected chisq + Double_t totalExpectedChisq = 0.0; + std::vector expectedChisqPerHisto; + fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto); + cout << endl << endl << ">> chisq = " << val << ", NDF = " << ndf << ", chisq/NDF = " << val/ndf; + + if (totalExpectedChisq != 0.0) { + cout << endl << ">> expected chisq = " << totalExpectedChisq << ", NDF = " << ndf << ", expected chisq/NDF = " << totalExpectedChisq/ndf; + UInt_t ndf_histo = 0; + for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_histo > 0) + 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]; + } + + // clean up + expectedChisqPerHisto.clear(); + } } else { // max. log likelihood cout << endl << endl << ">> maxLH = " << val << ", NDF = " << ndf << ", maxLH/NDF = " << val/ndf; } @@ -1421,6 +1441,39 @@ Bool_t PFitter::ExecuteSave() return false; } + // handle expected chisq if applicable + if (fUseChi2) { + fParams = *(fRunInfo->GetMsrParamList()); // get the update parameters back + + // calculate expected chisq + std::vector param; + Double_t totalExpectedChisq = 0.0; + std::vector expectedChisqPerHisto; + std::vector ndfPerHisto; + + for (UInt_t i=0; iCalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto); + + if (totalExpectedChisq != 0.0) { // i.e. applicable for single histogram fits only + // get the ndf's of the histos + UInt_t ndf_histo; + for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + ndfPerHisto.push_back(ndf_histo); + } + + // feed the msr-file handler + PMsrStatisticStructure *statistics = fRunInfo->GetMsrStatistic(); + if (statistics) { + statistics->fMinExpected = totalExpectedChisq; + statistics->fMinExpectedPerHisto = expectedChisqPerHisto; + statistics->fNdfPerHisto = ndfPerHisto; + } + } + } + cout << ">> PFitter::ExecuteSave(): will write minuit2 output file ..." << endl; ofstream fout; diff --git a/src/classes/PFitterFcn.cpp b/src/classes/PFitterFcn.cpp index aef48dab..7c1203f6 100644 --- a/src/classes/PFitterFcn.cpp +++ b/src/classes/PFitterFcn.cpp @@ -65,6 +65,8 @@ PFitterFcn::~PFitterFcn() { } +//-------------------------------------------------------------------------- +// operator() //-------------------------------------------------------------------------- /** *

Minuit2 interface function call routine. This is the function which should be minimized. @@ -90,3 +92,48 @@ 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. + * + * \param par + * \param totalExpectedChisq expected chisq for all run blocks + * \param expectedChisqPerRun expected chisq vector for all the run blocks + */ +void PFitterFcn::CalcExpectedChiSquare(const std::vector &par, Double_t &totalExpectedChisq, std::vector &expectedChisqPerRun) +{ + // init expected chisq related variables + totalExpectedChisq = 0.0; + expectedChisqPerRun.clear(); + + // only do something for chisq + if (fUseChi2) { + Double_t value = 0.0; + 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); + totalExpectedChisq += value; + } + } +} diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index 8b291f67..7a013f2d 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -52,7 +52,8 @@ using namespace std; * * \param fileName name of a msr-file. */ -PMsrHandler::PMsrHandler(const Char_t *fileName) : fFileName(fileName) +PMsrHandler::PMsrHandler(const Char_t *fileName, const Bool_t writeExpectedChisq) : + fWriteExpectedChisq(writeExpectedChisq), fFileName(fileName) { // init variables fMsrBlockCounter = 0; @@ -64,6 +65,9 @@ PMsrHandler::PMsrHandler(const Char_t *fileName) : fFileName(fileName) fStatistic.fChisq = true; fStatistic.fMin = -1.0; fStatistic.fNdf = 0; + fStatistic.fMinExpected = 0.0; + fStatistic.fMinExpectedPerHisto.clear(); + fStatistic.fNdfPerHisto.clear(); fFuncHandler = 0; @@ -95,6 +99,8 @@ PMsrHandler::~PMsrHandler() fCommands.clear(); fPlots.clear(); fStatistic.fStatLines.clear(); + fStatistic.fMinExpectedPerHisto.clear(); + fStatistic.fNdfPerHisto.clear(); fParamInUse.clear(); if (fFuncHandler) { @@ -996,21 +1002,46 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) } else if (sstr.BeginsWith("chisq") || sstr.BeginsWith("maxLH")) { partialStatisticBlockFound = false; if (fStatistic.fValid) { // valid fit result - if (fStatistic.fChisq) - str = " chisq = "; - else - str = " maxLH = "; - str += fStatistic.fMin; - str += ", NDF = "; - str += fStatistic.fNdf; - if (fStatistic.fChisq) - str += ", chisq/NDF = "; - else - str += ", maxLH/NDF = "; - str += fStatistic.fMin / fStatistic.fNdf; - fout << str.Data() << endl; - if (messages) - cout << endl << str.Data() << endl; + if (fStatistic.fChisq) { + str.Form(" chisq = %.1lf, NDF = %d, chisq/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf); + fout << str.Data() << endl; + if (messages) + cout << endl << str.Data() << endl; + + // check if expected chisq needs to be written + if (fStatistic.fMinExpected != 0.0) { + str.Form(" expected chisq = %.1lf, NDF = %d, expected chisq/NDF = %lf", + fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf); + if (fWriteExpectedChisq) + fout << str.Data() << endl; + if (messages) + cout << endl << str.Data() << endl; + + 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]); + if (fWriteExpectedChisq) + fout << str.Data() << endl; + + if (messages) + cout << str.Data() << endl; + } + } + } + } else { // maxLH + str.Form(" maxLH = %.1lf, NDF = %d, maxLH/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf); + fout << str.Data() << endl; + if (messages) + cout << endl << str.Data() << endl; + } } else { fout << "*** FIT DID NOT CONVERGE ***" << endl; if (messages) @@ -1020,22 +1051,41 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) partialStatisticBlockFound = false; if (fStatistic.fValid) { // valid fit result if (fStatistic.fChisq) { // chisq - str = " chisq = "; - str += fStatistic.fMin; - str += ", NDF = "; - str += fStatistic.fNdf; - str += ", chisq/NDF = "; - str += fStatistic.fMin / fStatistic.fNdf; + str.Form(" chisq = %.1lf, NDF = %d, chisq/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf); fout << str.Data() << endl; if (messages) cout << endl << str.Data() << endl; + + // check if expected chisq needs to be written + if (fStatistic.fMinExpected != 0.0) { + str.Form(" expected chisq = %.1lf, NDF = %d, expected chisq/NDF = %lf", + fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf); + if (fWriteExpectedChisq) + fout << str.Data() << endl; + if (messages) + cout << str.Data() << endl; + + 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]); + if (fWriteExpectedChisq) + fout << str.Data() << endl; + + if (messages) + cout << str.Data() << endl; + } + } + } } else { // max. log. liklihood - str = " maxLH = "; - str += fStatistic.fMin; - str += ", NDF = "; - str += fStatistic.fNdf; - str += ", maxLH/NDF = "; - str += fStatistic.fMin / fStatistic.fNdf; + str.Form(" maxLH = %.1lf, NDF = %d, maxLH/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf); fout << str.Data() << endl; if (messages) cout << endl << str.Data() << endl; @@ -1046,11 +1096,15 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) cout << endl << "*** FIT DID NOT CONVERGE ***" << endl; } } else { - if (str.Length() > 0) - fout << str.Data() << endl; - else // only write endl if not eof is reached. This is preventing growing msr-files, i.e. more and more empty lines at the end of the file + if (str.Length() > 0) { + sstr = str; + sstr.Remove(TString::kLeading, ' '); + if (!sstr.BeginsWith("expected chisq") && !sstr.BeginsWith("run block")) + fout << str.Data() << endl; + } else { // only write endl if not eof is reached. This is preventing growing msr-files, i.e. more and more empty lines at the end of the file if (!fin.eof()) fout << endl; + } } break; default: @@ -1066,28 +1120,47 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) TDatime dt; fout << "STATISTIC --- " << dt.AsSQLString() << endl; if (fStatistic.fValid) { // valid fit result - if (fStatistic.fChisq) { // chisq - str = " chisq = "; - str += fStatistic.fMin; - str += ", NDF = "; - str += fStatistic.fNdf; - str += ", chisq/NDF = "; - str += fStatistic.fMin / fStatistic.fNdf; + if (fStatistic.fChisq) { + str.Form(" chisq = %.1lf, NDF = %d, chisq/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf); fout << str.Data() << endl; if (messages) cout << endl << str.Data() << endl; - } else { // max. log. liklihood - str = " maxLH = "; - str += fStatistic.fMin; - str += ", NDF = "; - str += fStatistic.fNdf; - str += ", maxLH/NDF = "; - str += fStatistic.fMin / fStatistic.fNdf; + + // check if expected chisq needs to be written + if (fStatistic.fMinExpected != 0.0) { + str.Form(" expected chisq = %.1lf, NDF = %d, expected chisq/NDF = %lf", + fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf); + if (fWriteExpectedChisq) + fout << str.Data() << endl; + if (messages) + cout << str.Data() << endl; + + 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]); + if (fWriteExpectedChisq) + fout << str.Data() << endl; + + if (messages) + cout << str.Data() << endl; + } + } + } + } else { // maxLH + str.Form(" maxLH = %.1lf, NDF = %d, maxLH/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf); fout << str.Data() << endl; if (messages) cout << endl << str.Data() << endl; } - } else { + } else { fout << "*** FIT DID NOT CONVERGE ***" << endl; if (messages) cout << endl << "*** FIT DID NOT CONVERGE ***" << endl; @@ -1102,28 +1175,47 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) fout << "STATISTIC --- " << dt.AsSQLString() << endl; if (fStatistic.fValid) { // valid fit result if (fStatistic.fChisq) { // chisq - str = " chisq = "; - str += fStatistic.fMin; - str += ", NDF = "; - str += fStatistic.fNdf; - str += ", chisq/NDF = "; - str += fStatistic.fMin / fStatistic.fNdf; + str.Form(" chisq = %.1lf, NDF = %d, chisq/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf); fout << str.Data() << endl; if (messages) cout << endl << str.Data() << endl; + + // check if expected chisq needs to be written + if (fStatistic.fMinExpected != 0.0) { + str.Form(" expected chisq = %.1lf, NDF = %d, expected chisq/NDF = %lf", + fStatistic.fMinExpected, fStatistic.fNdf, fStatistic.fMinExpected/fStatistic.fNdf); + if (fWriteExpectedChisq) + fout << str.Data() << endl; + if (messages) + cout << str.Data() << endl; + + 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]); + if (fWriteExpectedChisq) + fout << str.Data() << endl; + + if (messages) + cout << str.Data() << endl; + } + } + } } else { // max. log. liklihood - str = " maxLH = "; - str += fStatistic.fMin; - str += ", NDF = "; - str += fStatistic.fNdf; - str += ", maxLH/NDF = "; - str += fStatistic.fMin / fStatistic.fNdf; + str.Form(" maxLH = %.1lf, NDF = %d, maxLH/NDF = %lf", fStatistic.fMin, fStatistic.fNdf, fStatistic.fMin / fStatistic.fNdf); fout << str.Data() << endl; if (messages) cout << endl << str.Data() << endl; } } else { - fout << "*** FIT DID NOT CONVERGE ***" << endl; + fout << "*** FIT DID NOT CONVERGE (4) ***" << endl; if (messages) cout << endl << "*** FIT DID NOT CONVERGE ***" << endl; } @@ -3813,7 +3905,8 @@ Bool_t PMsrHandler::HandleStatisticEntry(PMsrLines &lines) tstr.Remove(TString::kLeading, ' '); if (tstr.Length() > 0) { if (!tstr.BeginsWith("#") && !tstr.BeginsWith("STATISTIC") && !tstr.BeginsWith("chisq") && - !tstr.BeginsWith("maxLH") && !tstr.BeginsWith("*** FIT DID NOT CONVERGE ***")) { + !tstr.BeginsWith("maxLH") && !tstr.BeginsWith("*** FIT DID NOT CONVERGE ***") && + !tstr.BeginsWith("expected chisq") && !tstr.BeginsWith("run block")) { cerr << endl << ">> PMsrHandler::HandleStatisticEntry: **SYNTAX ERROR** in line " << lines[i].fLineNo; cerr << endl << ">> '" << lines[i].fLine.Data() << "'"; cerr << endl << ">> not a valid STATISTIC block line"; @@ -3872,6 +3965,185 @@ Bool_t PMsrHandler::HandleStatisticEntry(PMsrLines &lines) return true; } + +//-------------------------------------------------------------------------- +// GetNoOfFitParameters (public) +//-------------------------------------------------------------------------- +/** + *

Calculate the number of fit parameters for single histo fits. + * + * \param idx run block index + */ +UInt_t PMsrHandler::GetNoOfFitParameters(UInt_t idx) +{ + UInt_t noOfFitParameters = 0; + PUIntVector paramVector; + PUIntVector funVector; + PUIntVector mapVector; + TObjArray *tokens = 0; + TObjString *ostr = 0; + TString str; + UInt_t k, dval; + Int_t status, pos; + + // check that idx is valid + if (idx >= fRuns.size()) { + cerr << endl << ">> PMsrHandler::GetNoOfFitParameters() **ERROR** idx=" << idx << ", out of range fRuns.size()=" << fRuns.size(); + cerr << endl; + return 0; + } + + // get N0 parameter, possible parameter number or function + if (fRuns[idx].GetNormParamNo() != -1) { + if (fRuns[idx].GetNormParamNo() < MSR_PARAM_FUN_OFFSET) // parameter + paramVector.push_back(fRuns[idx].GetNormParamNo()); + else // function + funVector.push_back(fRuns[idx].GetNormParamNo() - MSR_PARAM_FUN_OFFSET); + } + + // get background parameter, for the case the background is fitted. + if (fRuns[idx].GetBkgFitParamNo() != -1) + paramVector.push_back(fRuns[idx].GetBkgFitParamNo()); + + // go through the theory block and collect parameters + // possible entries: number -> parameter, fun -> function, map -> maps + for (UInt_t i=0; i= 0) + str.Resize(pos); + // tokenize + tokens = str.Tokenize(" \t"); + if (!tokens) { + mapVector.clear(); + funVector.clear(); + paramVector.clear(); + return 0; + } + + for (Int_t j=0; jGetEntries(); j++) { + ostr = dynamic_cast(tokens->At(j)); + str = ostr->GetString(); + // check for parameter number + if (str.IsDigit()) { + dval = str.Atoi(); + paramVector.push_back(dval); + } + + // check for map + if (str.Contains("map")) { + status = sscanf(str.Data(), "map%d", &dval); + if (status == 1) { + mapVector.push_back(dval); + } + } + + // check for function + if (str.Contains("fun")) { + status = sscanf(str.Data(), "fun%d", &dval); + if (status == 1) { + funVector.push_back(dval); + } + } + } + + delete tokens; + tokens = 0; + } + + // go through the function block and collect parameters + for (UInt_t i=0; i> PMsrHandler::GetNoOfFitParameters() **ERROR** couldn't find fun" << funVector[i]; + cerr << endl << endl; + + // clean up + mapVector.clear(); + funVector.clear(); + paramVector.clear(); + + return 0; + } + + // remove potential comments + str = fFunctions[k].fLine; + pos = str.Index('#'); + if (pos >= 0) + str.Resize(pos); + + // tokenize + tokens = str.Tokenize(" \t"); + if (!tokens) { + mapVector.clear(); + funVector.clear(); + paramVector.clear(); + return 0; + } + + // filter out parameters and maps + for (Int_t j=0; jGetEntries(); j++) { + ostr = dynamic_cast(tokens->At(j)); + str = ostr->GetString(); + + // check for parameter + if (str.BeginsWith("par")) { + status = sscanf(str.Data(), "par%d", &dval); + if (status == 1) + paramVector.push_back(dval); + } + + // check for map + if (str.BeginsWith("map")) { + status = sscanf(str.Data(), "map%d", &dval); + if (status == 1) + mapVector.push_back(dval); + } + } + } + + // go through the map and collect parameters + for (UInt_t i=0; i& pa return chisq; } +//-------------------------------------------------------------------------- +// 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 //-------------------------------------------------------------------------- @@ -319,6 +341,29 @@ Double_t PRunListCollection::GetNonMusrMaximumLikelihood(const std::vectorNumber of bins in run block idx to be fitted. Only used for single histogram + * fitting together with the expected chisq. + * + * return: + * - number of bins fitted. + * + * \param idx run block index + */ +UInt_t PRunListCollection::GetNoOfBinsFitted(const UInt_t idx) const +{ + UInt_t result = 0; + + if (idx < fRunSingleHistoList.size()) + result = fRunSingleHistoList[idx]->GetNoOfFitBins(); + + return result; +} + + //-------------------------------------------------------------------------- // GetTotalNoOfBinsFitted //-------------------------------------------------------------------------- diff --git a/src/classes/PRunSingleHisto.cpp b/src/classes/PRunSingleHisto.cpp index 8c787c64..b41487a4 100644 --- a/src/classes/PRunSingleHisto.cpp +++ b/src/classes/PRunSingleHisto.cpp @@ -185,6 +185,99 @@ Double_t PRunSingleHisto::CalcChiSquare(const std::vector& par) return chisq; } +//-------------------------------------------------------------------------- +// CalcChiSquareExpected (public) +//-------------------------------------------------------------------------- +/** + *

Calculate expected chi-square. + * + * return: + * - chisq value + * + * \param par parameter vector iterated by minuit2 + */ +Double_t PRunSingleHisto::CalcChiSquareExpected(const std::vector& par) +{ + Double_t chisq = 0.0; + Double_t diff = 0.0; + Double_t theo = 0.0; + + Double_t N0 = 0.0; + + // check if norm is a parameter or a function + if (fRunInfo->GetNormParamNo() < MSR_PARAM_FUN_OFFSET) { // norm is a parameter + N0 = par[fRunInfo->GetNormParamNo()-1]; + } else { // norm is a function + // get function number + UInt_t funNo = fRunInfo->GetNormParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + N0 = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + + // get tau + Double_t tau; + if (fRunInfo->GetLifetimeParamNo() != -1) + tau = par[fRunInfo->GetLifetimeParamNo()-1]; + else + tau = PMUON_LIFETIME; + + // get background + Double_t bkg; + if (fRunInfo->GetBkgFitParamNo() == -1) { // bkg not fitted + if (fRunInfo->GetBkgFix(0) == PMUSR_UNDEFINED) { // no fixed background given (background interval) + bkg = fBackground; + } else { // fixed bkg given + bkg = fRunInfo->GetBkgFix(0); + } + } else { // bkg fitted + bkg = par[fRunInfo->GetBkgFitParamNo()-1]; + } + + // calculate functions + for (Int_t i=0; iGetNoOfFuncs(); i++) { + Int_t funcNo = fMsrInfo->GetFuncNo(i); + fFuncValues[i] = fMsrInfo->EvalFunc(funcNo, *fRunInfo->GetMap(), par); + } + + // calculate chi square + Double_t time(1.0); + Int_t i, N(static_cast(fData.GetValue()->size())); + + // In order not to have an IF in the next loop, determine the start and end bins for the fit range now + Int_t startTimeBin = static_cast(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())); + if (startTimeBin < 0) + startTimeBin = 0; + Int_t endTimeBin = static_cast(floor((fFitEndTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())) + 1; + if (endTimeBin > N) + endTimeBin = N; + + // Calculate the theory function once to ensure one function evaluation for the current set of parameters. + // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once + // for a given set of parameters---which should be done outside of the parallelized loop. + // For all other functions it means a tiny and acceptable overhead. + time = fTheory->Func(time, par, fFuncValues); + + #ifdef HAVE_GOMP + Int_t chunk = (endTimeBin - startTimeBin)/omp_get_num_procs(); + if (chunk < 10) + chunk = 10; + #pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq) + #endif + for (i=startTimeBin; i < endTimeBin; ++i) { + time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); + theo = N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg; + diff = fData.GetValue()->at(i) - theo; + chisq += diff*diff / theo; + } + + // the correction factor is need since the data scales like pack*t_res, + // whereas the error scales like sqrt(pack*t_res) + if (fScaleN0AndBkg) + chisq *= fRunInfo->GetPacking() * (fTimeResolution * 1.0e3); + + return chisq; +} + //-------------------------------------------------------------------------- // CalcMaxLikelihood (public) //-------------------------------------------------------------------------- @@ -360,14 +453,6 @@ void PRunSingleHisto::CalcTheory() */ UInt_t PRunSingleHisto::GetNoOfFitBins() { -// fNoOfFitBins=0; -// -// Double_t time; -// 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/PStartupHandler.cpp b/src/classes/PStartupHandler.cpp index 160facd7..6416205f 100644 --- a/src/classes/PStartupHandler.cpp +++ b/src/classes/PStartupHandler.cpp @@ -165,6 +165,8 @@ void PStartupHandler::OnStartDocument() fFourierDefaults.fPlotRange[0] = -1.0; fFourierDefaults.fPlotRange[1] = -1.0; fFourierDefaults.fPhaseIncrement = 1.0; + + fWriteExpectedChisq = false; } //-------------------------------------------------------------------------- @@ -193,6 +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, "marker")) { fKey = eMarker; } else if (!strcmp(str, "color")) { @@ -247,6 +251,11 @@ void PStartupHandler::OnCharacters(const Char_t *str) // add str to the path list fDataPathList.push_back(str); break; + case eWriteExpectedChisq: + tstr = TString(str); + if (tstr.BeginsWith("y") || tstr.BeginsWith("Y")) + fWriteExpectedChisq = true; + break; case eMarker: // check that str is a number tstr = TString(str); diff --git a/src/classes/PTheory.cpp b/src/classes/PTheory.cpp index 25441666..4b1ae7ce 100644 --- a/src/classes/PTheory.cpp +++ b/src/classes/PTheory.cpp @@ -307,8 +307,6 @@ PTheory::PTheory(PMsrHandler *msrInfo, UInt_t runNo, const Bool_t hasParent) : f fUserParam.resize(fParamNo.size()); } -//cout << endl << "debug> fUserFcn=" << fUserFcn << ", fUserFcn->NeedGlobalPart()=" << fUserFcn->NeedGlobalPart() << ", gGlobalUserFcn=" << gGlobalUserFcn << endl; - // check if the global part of the user function is needed if (fUserFcn->NeedGlobalPart()) { fUserFcn->SetGlobalPart(gGlobalUserFcn, fUserFcnIdx); // invoke or retrieve global user function object diff --git a/src/include/PFitterFcn.h b/src/include/PFitterFcn.h index b6d9e876..0e9215dc 100644 --- a/src/include/PFitterFcn.h +++ b/src/include/PFitterFcn.h @@ -48,9 +48,11 @@ class PFitterFcn : public ROOT::Minuit2::FCNBase ~PFitterFcn(); Double_t Up() const { return fUp; } - Double_t operator()(const std::vector&) const; + Double_t operator()(const std::vector &par) const; UInt_t GetTotalNoOfFittedBins() { return fRunListCollection->GetTotalNoOfBinsFitted(); } + UInt_t GetNoOfFittedBins(const UInt_t idx); + void CalcExpectedChiSquare(const std::vector &par, Double_t &totalExpectedChisq, std::vector &expectedChisqPerRun); private: Double_t fUp; ///< for chisq == 1.0, i.e. errors are 1 std. deviation errors. for log max-likelihood == 0.5, i.e. errors are 1 std. deviation errors (for details see the minuit2 user manual). diff --git a/src/include/PMsrHandler.h b/src/include/PMsrHandler.h index d12c9f8a..003d2bc1 100644 --- a/src/include/PMsrHandler.h +++ b/src/include/PMsrHandler.h @@ -47,7 +47,7 @@ class PMsrHandler { public: - PMsrHandler(const Char_t *fileName); + PMsrHandler(const Char_t *fileName, const Bool_t writeExpectedChisq=false); virtual ~PMsrHandler(); virtual Int_t ReadMsrFile(); @@ -95,6 +95,7 @@ class PMsrHandler virtual Double_t EvalFunc(UInt_t i, vector map, vector param) { return fFuncHandler->Eval(i,map,param); } + virtual UInt_t GetNoOfFitParameters(UInt_t idx); virtual Int_t ParameterInUse(UInt_t paramNo); virtual Bool_t CheckRunBlockIntegrity(); virtual Bool_t CheckUniquenessOfParamNames(UInt_t &parX, UInt_t &parY); @@ -105,6 +106,8 @@ class PMsrHandler virtual void CheckMaxLikelihood(); private: + Bool_t fWriteExpectedChisq; ///< flag shows if expected chisq shall be written to the msr-file + TString fFileName; ///< file name of the msr-file TString fMsrFileDirectoryPath; ///< msr-file directory path TString fTitle; ///< holds the title string of the msr-file diff --git a/src/include/PMusr.h b/src/include/PMusr.h index ba615df1..7fd56a6d 100644 --- a/src/include/PMusr.h +++ b/src/include/PMusr.h @@ -582,9 +582,12 @@ typedef struct { Bool_t fValid; ///< flag showing if the statistics block is valid, i.e. a fit took place which converged 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.likelyhood - Double_t fMin; ///< chi2 or max. likelyhood + Bool_t fChisq; ///< flag telling if min = chi2 or min = max.likelihood + Double_t fMin; ///< chi2 or max. likelihood 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 + PUIntVector fNdfPerHisto; ///< number of degrees of freedom per histo } PMsrStatisticStructure; //------------------------------------------------------------- diff --git a/src/include/PRunBase.h b/src/include/PRunBase.h index 48d6dd1b..650adfcc 100644 --- a/src/include/PRunBase.h +++ b/src/include/PRunBase.h @@ -74,8 +74,6 @@ class PRunBase PMsrRunBlock *fRunInfo; ///< run info used to filter out needed infos of a run PRunDataHandler *fRawData; ///< holds the raw run data - PIntVector fParamNo; ///< vector of parameter numbers for the specifc run - PRunData fData; ///< data to be fitted, viewed, i.e. binned data Double_t fTimeResolution; ///< time resolution in (us) PIntVector fT0s; ///< all t0's of a run! The derived classes will handle it diff --git a/src/include/PRunListCollection.h b/src/include/PRunListCollection.h index d424cc31..ce93454f 100644 --- a/src/include/PRunListCollection.h +++ b/src/include/PRunListCollection.h @@ -59,6 +59,7 @@ 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; @@ -68,6 +69,7 @@ class PRunListCollection virtual Double_t GetMuMinusMaximumLikelihood(const std::vector& par) const; virtual Double_t GetNonMusrMaximumLikelihood(const std::vector& par) const; + virtual UInt_t GetNoOfBinsFitted(const UInt_t idx) const; virtual UInt_t GetTotalNoOfBinsFitted() const; virtual UInt_t GetNoOfSingleHisto() const { return fRunSingleHistoList.size(); } ///< returns the number of single histogram data sets present in the msr-file @@ -91,10 +93,10 @@ class PRunListCollection PMsrHandler *fMsrInfo; ///< pointer to the msr-file handler PRunDataHandler *fData; ///< pointer to the run-data handler - vector fRunSingleHistoList; ///< stores all precessed single histogram data - vector fRunAsymmetryList; ///< stores all precessed asymmetry data - vector fRunMuMinusList; ///< stores all precessed mu-minus data - vector fRunNonMusrList; ///< stores all precessed non-muSR data + vector fRunSingleHistoList; ///< stores all processed single histogram data + vector fRunAsymmetryList; ///< stores all processed asymmetry data + vector fRunMuMinusList; ///< stores all processed mu-minus data + vector fRunNonMusrList; ///< stores all processed non-muSR data }; #endif // _PRUNLISTCOLLECTION_H_ diff --git a/src/include/PRunSingleHisto.h b/src/include/PRunSingleHisto.h index db2c0f96..3b147d10 100644 --- a/src/include/PRunSingleHisto.h +++ b/src/include/PRunSingleHisto.h @@ -45,6 +45,7 @@ class PRunSingleHisto : public PRunBase virtual ~PRunSingleHisto(); 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 3854c303..17abbce9 100644 --- a/src/include/PStartupHandler.h +++ b/src/include/PStartupHandler.h @@ -75,13 +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 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, + enum EKeyWords {eEmpty, eComment, eDataPath, eWriteExpectedChisq, eFourierSettings, eUnits, eFourierPower, eApodization, ePlot, ePhase, ePhaseIncrement, eRootSettings, eMarkerList, eMarker, eColorList, eColor}; @@ -89,6 +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 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 9e9a58ef..f3810eb5 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); + PMsrHandler *msrHandler = new PMsrHandler(filename, startupHandler->GetWriteExpectedChisq()); status = msrHandler->ReadMsrFile(); if (status != PMUSR_SUCCESS) { switch (status) { diff --git a/src/musrfit_startup.xml b/src/musrfit_startup.xml index af0f8010..8e7e1492 100644 --- a/src/musrfit_startup.xml +++ b/src/musrfit_startup.xml @@ -13,6 +13,7 @@ /afs/psi.ch/project/bulkmusr/data/gpd /afs/psi.ch/project/bulkmusr/data/ltf /afs/psi.ch/project/bulkmusr/data/alc + n Gauss 0