diff --git a/src/classes/PFitter.cpp b/src/classes/PFitter.cpp index 9040b60f..b0ad8f73 100644 --- a/src/classes/PFitter.cpp +++ b/src/classes/PFitter.cpp @@ -85,13 +85,16 @@ PSectorChisq::PSectorChisq(UInt_t noOfRuns) : fNoOfRuns(noOfRuns) // init fLast = 0.0; fChisq = 0.0; + fExpectedChisq = 0.0; fNDF = 0; fFirst.resize(fNoOfRuns); fChisqRun.resize(fNoOfRuns); + fExpectedChisqRun.resize(fNoOfRuns); fNDFRun.resize(fNoOfRuns); for (UInt_t i=0; iSet the expected chisq/maxLH for a given detector with index idx. + * + * @param chisq expected chisq/maxLH to be set + * @param idx index of the run to be set + */ +void PSectorChisq::SetExpectedChisq(Double_t chisq, UInt_t idx) +{ + if (idx > fNoOfRuns) { + std::cerr << "**WARNING** from PSectorChisq::SetExpectedChisq. It tries to set" << std::endl; + std::cerr << " a chisq with idx=" << idx << " which is larger than #RUNS=" << fNoOfRuns << "." << std::endl; + std::cerr << " Will ignore it, but you better check what is going on!" << std::endl; + return; + } + + fExpectedChisqRun[idx] = chisq; +} + //-------------------------------------------------------------------------- // SetNDF //-------------------------------------------------------------------------- @@ -147,7 +171,7 @@ void PSectorChisq::SetChisq(Double_t chisq, UInt_t idx) * @param ndf to be set * @param idx index of the run to be set */ -void PSectorChisq::SetNDF(Double_t ndf, UInt_t idx) +void PSectorChisq::SetNDF(UInt_t ndf, UInt_t idx) { if (idx > fNoOfRuns) { std::cerr << "**WARNING** from PSectorChisq::SetNDF. It tries to set" << std::endl; @@ -182,7 +206,7 @@ Double_t PSectorChisq::GetTimeRangeFirst(UInt_t idx) // GetChisq //-------------------------------------------------------------------------- /** - *

chisq/maxLH of the run with index idx + *

Get chisq/maxLH of the run with index idx * * @param idx index of the run * @@ -197,10 +221,28 @@ Double_t PSectorChisq::GetChisq(UInt_t idx) } //-------------------------------------------------------------------------- -// GetChisq +// GetExpectedChisq //-------------------------------------------------------------------------- /** - *

NDF of the run with index idx + *

Get expected chisq/maxLH of the run with index idx + * + * @param idx index of the run + * + * return: chisq/maxLH of the requested run or -1.0 if the idx is out-of-range. + */ +Double_t PSectorChisq::GetExpectedChisq(UInt_t idx) +{ + if (idx >= fNoOfRuns) + return -1.0; + + return fExpectedChisqRun[idx]; +} + +//-------------------------------------------------------------------------- +// GetNDF +//-------------------------------------------------------------------------- +/** + *

Get NDF of the run with index idx * * @param idx index of the run * @@ -354,11 +396,6 @@ Bool_t PFitter::DoFit() // check if only chisq/maxLH shall be calculated once if (fChisqOnly) { - - // check for sector command - if (fSectorFlag) - PrepareSector(); - std::vector param = fMnUserParams.Params(); std::vector error = fMnUserParams.Errors(); Int_t usedParams = 0; @@ -366,20 +403,20 @@ Bool_t PFitter::DoFit() if (error[i] != 0.0) usedParams++; } - Int_t ndf = 0; + UInt_t ndf = 0; Double_t val = 0.0; if (fDKSReady) { - ndf = static_cast(fFitterFcnDKS->GetTotalNoOfFittedBins()) - usedParams; + ndf = fFitterFcnDKS->GetTotalNoOfFittedBins() - static_cast(usedParams); val = (*fFitterFcnDKS)(param); } else { - ndf = static_cast(fFitterFcn->GetTotalNoOfFittedBins()) - usedParams; + ndf = fFitterFcn->GetTotalNoOfFittedBins() - static_cast(usedParams); val = (*fFitterFcn)(param); } if (fUseChi2) { // calculate expected chisq Double_t totalExpectedChisq = 0.0; - std::vector expectedChisqPerRun; + PDoubleVector expectedChisqPerRun; if (fDKSReady) fFitterFcnDKS->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerRun); else @@ -394,25 +431,25 @@ Bool_t PFitter::DoFit() if (totalExpectedChisq != 0.0) { std::cout << std::endl << ">> expected chisq = " << totalExpectedChisq << ", NDF = " << ndf << ", expected chisq/NDF = " << totalExpectedChisq/ndf; - UInt_t ndf_histo = 0; + UInt_t ndf_run = 0; for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + ndf_run = fFitterFcnDKS->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); else - ndf_histo = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); - if (ndf_histo > 0) - std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq/red.chisq_e) = (" << ndf_histo << "/" << chisqPerRun[i]/ndf_histo << "/" << expectedChisqPerRun[i]/ndf_histo << ")"; + ndf_run = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_run > 0) + std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq/red.chisq_e) = (" << ndf_run << "/" << chisqPerRun[i]/ndf_run << "/" << expectedChisqPerRun[i]/ndf_run << ")"; } } else if (chisqPerRun.size() > 0) { // in case expected chisq is not applicable like for asymmetry fits - UInt_t ndf_histo = 0; + UInt_t ndf_run = 0; for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + ndf_run = fFitterFcnDKS->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); else - ndf_histo = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + ndf_run = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); - if (ndf_histo > 0) - std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq) = (" << ndf_histo << "/" << chisqPerRun[i]/ndf_histo << ")"; + if (ndf_run > 0) + std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq) = (" << ndf_run << "/" << chisqPerRun[i]/ndf_run << ")"; } } @@ -430,10 +467,10 @@ Bool_t PFitter::DoFit() fRunListCollection->SetFitRange(secFitRange); // calculate chisq and NDF if (fDKSReady) { - ndf = static_cast(fFitterFcnDKS->GetTotalNoOfFittedBins()) - usedParams; + ndf = fFitterFcnDKS->GetTotalNoOfFittedBins() - static_cast(usedParams); val = (*fFitterFcnDKS)(param); } else { - ndf = static_cast(fFitterFcn->GetTotalNoOfFittedBins()) - usedParams; + ndf = fFitterFcn->GetTotalNoOfFittedBins() - static_cast(usedParams); val = (*fFitterFcn)(param); } // calculate expected chisq @@ -449,29 +486,29 @@ Bool_t PFitter::DoFit() std::cout << std::endl; std::cout << "++++" << std::endl; - std::cout << ">> Sector " << k << ": FitRange: " << secFitRange[0].first << ", " << secFitRange[0].second << std::endl; + std::cout << ">> Sector " << k+1 << ": FitRange: " << secFitRange[0].first << ", " << secFitRange[0].second << std::endl; std::cout << ">> chisq = " << val << ", NDF = " << ndf << ", chisq/NDF = " << val/ndf; if (totalExpectedChisq != 0.0) { std::cout << std::endl << ">> expected chisq = " << totalExpectedChisq << ", NDF = " << ndf << ", expected chisq/NDF = " << totalExpectedChisq/ndf; - UInt_t ndf_histo = 0; + UInt_t ndf_run = 0; for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + ndf_run = fFitterFcnDKS->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); else - ndf_histo = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); - if (ndf_histo > 0) - std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq/red.chisq_e) = (" << ndf_histo << "/" << chisqPerRun[i]/ndf_histo << "/" << expectedChisqPerRun[i]/ndf_histo << ")"; + ndf_run = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_run > 0) + std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq/red.chisq_e) = (" << ndf_run << "/" << chisqPerRun[i]/ndf_run << "/" << expectedChisqPerRun[i]/ndf_run << ")"; } } else if (chisqPerRun.size() > 0) { // in case expected chisq is not applicable like for asymmetry fits - UInt_t ndf_histo = 0; + UInt_t ndf_run = 0; for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + ndf_run = fFitterFcnDKS->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); else - ndf_histo = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); - if (ndf_histo > 0) - std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq) = (" << ndf_histo << "/" << chisqPerRun[i]/ndf_histo << ")"; + ndf_run = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_run > 0) + std::cout << std::endl << ">> run block " << i+1 << ": (NDF/red.chisq) = (" << ndf_run << "/" << chisqPerRun[i]/ndf_run << ")"; } } // clean up @@ -497,14 +534,14 @@ Bool_t PFitter::DoFit() if (totalExpectedMaxLH != 0.0) { std::cout << std::endl << ">> expected maxLH = " << totalExpectedMaxLH << ", NDF = " << ndf << ", expected maxLH/NDF = " << totalExpectedMaxLH/ndf; - UInt_t ndf_histo = 0; + UInt_t ndf_run = 0; for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + ndf_run = fFitterFcnDKS->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); else - ndf_histo = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); - if (ndf_histo > 0) - std::cout << std::endl << ">> run block " << i+1 << ": (NDF/maxLH.chisq/maxLH.chisq_e) = (" << ndf_histo << "/" << maxLHPerRun[i]/ndf_histo << "/" << expectedMaxLHPerRun[i]/ndf_histo << ")"; + ndf_run = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_run > 0) + std::cout << std::endl << ">> run block " << i+1 << ": (NDF/maxLH.chisq/maxLH.chisq_e) = (" << ndf_run << "/" << maxLHPerRun[i]/ndf_run << "/" << expectedMaxLHPerRun[i]/ndf_run << ")"; } } @@ -541,19 +578,19 @@ Bool_t PFitter::DoFit() std::cout << std::endl; std::cout << "++++" << std::endl; - std::cout << ">> Sector " << k << ": FitRange: " << secFitRange[0].first << ", " << secFitRange[0].second << std::endl; + std::cout << ">> Sector " << k+1 << ": FitRange: " << secFitRange[0].first << ", " << secFitRange[0].second << std::endl; std::cout << ">> maxLH = " << val << ", NDF = " << ndf << ", maxLH/NDF = " << val/ndf; if (totalExpectedMaxLH != 0.0) { std::cout << std::endl << ">> expected maxLH = " << totalExpectedMaxLH << ", NDF = " << ndf << ", expected maxLH/NDF = " << totalExpectedMaxLH/ndf; - UInt_t ndf_histo = 0; + UInt_t ndf_run = 0; for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + ndf_run = fFitterFcnDKS->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); else - ndf_histo = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); - if (ndf_histo > 0) - std::cout << std::endl << ">> run block " << i+1 << ": (NDF/maxLH.chisq/maxLH.chisq_e) = (" << ndf_histo << "/" << maxLHPerRun[i]/ndf_histo << "/" << expectedMaxLHPerRun[i]/ndf_histo << ")"; + ndf_run = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_run > 0) + std::cout << std::endl << ">> run block " << i+1 << ": (NDF/maxLH.chisq/maxLH.chisq_e) = (" << ndf_run << "/" << maxLHPerRun[i]/ndf_run << "/" << expectedMaxLHPerRun[i]/ndf_run << ")"; } } @@ -2067,7 +2104,7 @@ Bool_t PFitter::ExecuteSave(Bool_t firstSave) // calculate expected chisq std::vector param; Double_t totalExpectedChisq = 0.0; - std::vector expectedChisqPerHisto; + std::vector expectedchisqPerRun; std::vector ndfPerHisto; for (UInt_t i=0; iCalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto); + fFitterFcnDKS->CalcExpectedChiSquare(param, totalExpectedChisq, expectedchisqPerRun); else - fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerHisto); + fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedchisqPerRun); // calculate chisq per run - std::vector chisqPerHisto; + std::vector chisqPerRun; for (UInt_t i=0; iGetMsrRunList()->size(); i++) { if (fUseChi2) - chisqPerHisto.push_back(fRunListCollection->GetSingleRunChisq(param, i)); + chisqPerRun.push_back(fRunListCollection->GetSingleRunChisq(param, i)); else - chisqPerHisto.push_back(fRunListCollection->GetSingleRunMaximumLikelihood(param, i)); + chisqPerRun.push_back(fRunListCollection->GetSingleRunMaximumLikelihood(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; - for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + ndf_run = fFitterFcnDKS->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); else - ndf_histo = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); - ndfPerHisto.push_back(ndf_histo); + ndf_run = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + ndfPerHisto.push_back(ndf_run); } // feed the msr-file handler PMsrStatisticStructure *statistics = fRunInfo->GetMsrStatistic(); if (statistics) { - statistics->fMinPerHisto = chisqPerHisto; + statistics->fMinPerHisto = chisqPerRun; statistics->fMinExpected = totalExpectedChisq; - statistics->fMinExpectedPerHisto = expectedChisqPerHisto; + statistics->fMinExpectedPerHisto = expectedchisqPerRun; 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; i 1) { // in case expected chisq is not applicable like for asymmetry fits + UInt_t ndf_run = 0; + for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + ndf_run = fFitterFcnDKS->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); else - ndf_histo = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); - ndfPerHisto.push_back(ndf_histo); + ndf_run = fFitterFcn->GetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + ndfPerHisto.push_back(ndf_run); } } + // check if sector command has been requested + if (fSectorFlag) { + PDoubleVector error; + for (UInt_t i=0; i> PFitter::ExecuteSave(): will write minuit2 output file ..." << std::endl; @@ -2194,7 +2240,7 @@ Bool_t PFitter::ExecuteSave(Bool_t firstSave) fout << std::endl; fElapsedTime.clear(); - // write global informations + // write global information fout << std::endl << " Fval() = " << mnState.Fval() << ", Edm() = " << mnState.Edm() << ", NFcn() = " << mnState.NFcn(); fout << std::endl; fout << std::endl << "*************************************************************************"; @@ -2414,6 +2460,41 @@ Bool_t PFitter::ExecuteSave(Bool_t firstSave) } else { fout << std::endl << " no correlation coefficients available"; } + + fout << std::endl; + fout << std::endl << "*************************************************************************"; + fout << std::endl << " chisq/maxLH RESULT "; + fout << std::endl << "*************************************************************************"; + PMsrStatisticStructure *statistics = fRunInfo->GetMsrStatistic(); + + // get time range and write it + Double_t fitStartTime = PMUSR_UNDEFINED, fitEndTime = PMUSR_UNDEFINED; + // first check if there is a global block with a valid time range + PMsrGlobalBlock *global = fRunInfo->GetMsrGlobal(); + fitStartTime = global->GetFitRange(0); + fitEndTime = global->GetFitRange(1); + if (fitStartTime == PMUSR_UNDEFINED) { // no global time range, hence take the one from the first run block + PMsrRunList *run = fRunInfo->GetMsrRunList(); + fitStartTime = run->at(0).GetFitRange(0); + fitEndTime = run->at(0).GetFitRange(1); + } + fout.setf(std::ios::floatfield); + fout << std::endl << " Time Range: " << fitStartTime << ", " << fitEndTime << std::endl; + if (fUseChi2) { + fout.setf(std::ios::fixed, std::ios::floatfield); + fout << std::endl << " chisq = " << std::setprecision(4) << statistics->fMin << ", NDF = " << statistics->fNdf << ", chisq/NDF = " << std::setprecision(6) << statistics->fMin/statistics->fNdf; + if (statistics->fMinExpected > 0) + fout << std::endl << " chisq_e = " << std::setprecision(4) << statistics->fMinExpected << ", NDF = " << statistics->fNdf << ", chisq_e/NDF = " << std::setprecision(6) << statistics->fMinExpected/statistics->fNdf; + } else { // maxLH + fout.setf(std::ios::fixed, std::ios::floatfield); + fout << std::endl << " maxLH = " << std::setprecision(4) << statistics->fMin << ", NDF = " << statistics->fNdf << ", maxLH/NDF = " << std::setprecision(6) << statistics->fMin/statistics->fNdf; + if (statistics->fMinExpected > 0) + fout << std::endl << " maxLH_e = " << std::setprecision(4) << statistics->fMinExpected << ", NDF = " << statistics->fNdf << ", maxLH_e/NDF = " << std::setprecision(6) << statistics->fMinExpected/statistics->fNdf; + } + + if (fSectorFlag) + ExecuteSector(fout); + fout << std::endl; fout << std::endl << "*************************************************************************"; fout << std::endl << " DONE "; @@ -2515,23 +2596,112 @@ Bool_t PFitter::ExecuteSimplex() //-------------------------------------------------------------------------- /** *

Collect all the necessary chisq/maxLH sector information. + * + * @param param parameter value vector of the converged fit. + * @param error step value vector of the converged fit. */ -void PFitter::PrepareSector() +void PFitter::PrepareSector(PDoubleVector ¶m, PDoubleVector &error) { - // NOT YET IMPLEMENTED //as35 + Double_t val; + UInt_t ndf; -/* - for (UInt_t i=0; i +++++" << std::endl; - std::cout << "debug> sector " << i << std::endl; - std::cout << "debug> noOfRuns: " << fSector[i].GetNoRuns() << std::endl; - for (UInt_t j=0; j ----" << std::endl; - std::cout << "debug> RUN BLOCK " << j << std::endl; - std::cout << "debug> Sector Time Range: " << fSector[i].GetTimeRangeFirst(j) << ", " << fSector[i].GetTimeRangeLast() << std::endl; + Int_t usedParams = 0; + for (UInt_t i=0; iSetFitRange(secFitRange); + // calculate chisq + val = (*fFitterFcn)(param); + fSector[k].SetChisq(val); + // calculate NDF + ndf = static_cast(fFitterFcn->GetTotalNoOfFittedBins()) - usedParams; + fSector[k].SetNDF(ndf); + // calculate expected chisq + totalExpectedChisq = 0.0; + fFitterFcn->CalcExpectedChiSquare(param, totalExpectedChisq, expectedChisqPerRun); + fSector[k].SetExpectedChisq(totalExpectedChisq); + // calculate chisq per run + for (UInt_t i=0; iGetMsrRunList()->size(); i++) { + chisqPerRun.push_back(fRunListCollection->GetSingleRunChisq(param, i)); + fSector[k].SetChisq(chisqPerRun[i], i); + fSector[k].SetExpectedChisq(expectedChisqPerRun[i], i); + } + + if (totalExpectedChisq != 0.0) { + UInt_t ndf_run = 0; + for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_run > 0) { + fSector[k].SetNDF(ndf_run, i); + } + } + } else if (chisqPerRun.size() > 0) { // in case expected chisq is not applicable like for asymmetry fits + UInt_t ndf_run = 0; + for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_run > 0) { + fSector[k].SetNDF(ndf_run, i); + } + } + } + // clean up + chisqPerRun.clear(); + expectedChisqPerRun.clear(); + } + } else { + Double_t totalExpectedMaxLH = 0.0; + PDoubleVector expectedMaxLHPerRun; + PDoubleVector maxLHPerRun; + for (UInt_t k=0; kSetFitRange(secFitRange); + // calculate maxLH + val = (*fFitterFcn)(param); + fSector[k].SetChisq(val); + // calculate NDF + ndf = static_cast(fFitterFcn->GetTotalNoOfFittedBins()) - usedParams; + fSector[k].SetNDF(ndf); + // calculate expected maxLH + totalExpectedMaxLH = 0.0; + fFitterFcn->CalcExpectedChiSquare(param, totalExpectedMaxLH, expectedMaxLHPerRun); + fSector[k].SetExpectedChisq(totalExpectedMaxLH); + // calculate maxLH per run + for (UInt_t i=0; iGetMsrRunList()->size(); i++) { + maxLHPerRun.push_back(fRunListCollection->GetSingleRunMaximumLikelihood(param, i)); + fSector[k].SetChisq(maxLHPerRun[i], i); + fSector[k].SetExpectedChisq(expectedMaxLHPerRun[i], i); + } + + if (totalExpectedMaxLH != 0.0) { + UInt_t ndf_run = 0; + for (UInt_t i=0; iGetNoOfFittedBins(i) - fRunInfo->GetNoOfFitParameters(i); + if (ndf_run > 0) { + fSector[k].SetNDF(ndf_run, i); + } + } + } + + // clean up + maxLHPerRun.clear(); + expectedMaxLHPerRun.clear(); } } -*/ } //-------------------------------------------------------------------------- @@ -2543,9 +2713,50 @@ void PFitter::PrepareSector() * * return: if the sector command was successful, otherwise return flase. */ -Bool_t PFitter::ExecuteSector() +Bool_t PFitter::ExecuteSector(std::ofstream &fout) { - // NOT YET IMPLEMENTED //as35 + fout << std::endl; + fout << std::endl << "*************************************************************************"; + fout << std::endl << " SECTOR RESULTS"; + fout << std::endl << "*************************************************************************"; + + if (fUseChi2) { + for (UInt_t i=0; i 0) { // in case expected chisq is not applicable like for asymmetry fits + for (UInt_t j=0; j 0) { // in case expected chisq is not applicable like for asymmetry fits + for (UInt_t j=0; j &par, Double_ } else { // log max. likelihood // single histo for (UInt_t i=0; iGetNoOfSingleHisto(); i++) { - value = fRunListCollection->GetSingleHistoMaximumLikelihoodExpected(par, i); // calculate the expected mlh for single histo run block 'i' + value = fRunListCollection->GetSingleRunMaximumLikelihoodExpected(par, i); // calculate the expected mlh for single histo run block 'i' expectedChisqPerRun.push_back(value); totalExpectedChisq += value; } // mu minus for (UInt_t i=0; iGetNoOfMuMinus(); i++) { - value = fRunListCollection->GetSingleHistoMaximumLikelihoodExpected(par, i); // calculate the expected maxLH for mu minus run block 'i' + value = fRunListCollection->GetSingleRunMaximumLikelihoodExpected(par, i); // calculate the expected maxLH for mu minus run block 'i' expectedChisqPerRun.push_back(value); totalExpectedChisq += value; } diff --git a/src/classes/PRunListCollection.cpp b/src/classes/PRunListCollection.cpp index 11cfcb03..534d9bcf 100644 --- a/src/classes/PRunListCollection.cpp +++ b/src/classes/PRunListCollection.cpp @@ -399,20 +399,21 @@ Double_t PRunListCollection::GetSingleRunChisqExpected(const std::vector fMsrInfo->GetMsrRunList()->size()) { - std::cerr << ">> PRunListCollection::GetSingleHistoChisqExpected() **ERROR** idx=" << idx << " is out of range [0.." << fMsrInfo->GetMsrRunList()->size() << "[" << std::endl << std::endl; + std::cerr << ">> PRunListCollection::GetSingleRunChisqExpected() **ERROR** idx=" << idx << " is out of range [0.." << fMsrInfo->GetMsrRunList()->size() << "[" << std::endl << std::endl; return expectedChisq; } + UInt_t subIdx = 0; Int_t type = fMsrInfo->GetMsrRunList()->at(idx).GetFitType(); if (type == -1) { // i.e. not found in the RUN block, try the GLOBAL block type = fMsrInfo->GetMsrGlobal()->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++; + subIdx = idx; + } else { // found in the RUN block + // count how many entries of this fit-type are present up to idx + for (UInt_t i=0; iGetMsrRunList()->at(i).GetFitType() == type) + subIdx++; + } } // return the chisq of the single run @@ -658,36 +659,37 @@ Double_t PRunListCollection::GetNonMusrMaximumLikelihood(const std::vectorCalculates expected mlh of the single histogram with run block index idx of a msr-file. + *

Calculates expected mlh of the run block index idx of a msr-file. * * return: - * - expected mlh of for a single histogram + * - expected mlh of for a single run block * * \param par fit parameter vector * \param idx run block index */ -Double_t PRunListCollection::GetSingleHistoMaximumLikelihoodExpected(const std::vector& par, const UInt_t idx) const +Double_t PRunListCollection::GetSingleRunMaximumLikelihoodExpected(const std::vector& par, const UInt_t idx) const { Double_t expected_mlh = 0.0; if (idx > fMsrInfo->GetMsrRunList()->size()) { - std::cerr << ">> PRunListCollection::GetSingleHistoMaximumLikelihoodExpected() **ERROR** idx=" << idx << " is out of range [0.." << fMsrInfo->GetMsrRunList()->size() << "[" << std::endl << std::endl; + std::cerr << ">> PRunListCollection::GetSingleRunMaximumLikelihoodExpected() **ERROR** idx=" << idx << " is out of range [0.." << fMsrInfo->GetMsrRunList()->size() << "[" << std::endl << std::endl; return expected_mlh; } + UInt_t subIdx = 0; Int_t type = fMsrInfo->GetMsrRunList()->at(idx).GetFitType(); if (type == -1) { // i.e. not found in the RUN block, try the GLOBAL block type = fMsrInfo->GetMsrGlobal()->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++; + subIdx = idx; + } else { // found in the RUN block + // count how many entries of this fit-type are present up to idx + for (UInt_t i=0; iGetMsrRunList()->at(i).GetFitType() == type) + subIdx++; + } } // return the mlh of the single run diff --git a/src/classes/PRunSingleHisto.cpp b/src/classes/PRunSingleHisto.cpp index b3d15f09..5d114688 100644 --- a/src/classes/PRunSingleHisto.cpp +++ b/src/classes/PRunSingleHisto.cpp @@ -276,9 +276,9 @@ Double_t PRunSingleHisto::CalcChiSquareExpected(const std::vector& par Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs(); if (chunk < 10) chunk = 10; - #pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq) + #pragma omp parallel for default(shared) private(i,time,theo,diff) schedule(dynamic,chunk) reduction(+:chisq) #endif - for (i=fStartTimeBin; i < fEndTimeBin; ++i) { + for (i=fStartTimeBin; i(i)*fData.GetDataTimeStep(); theo = N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg; diff = fData.GetValue()->at(i) - theo; @@ -373,9 +373,8 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector& par) time = fData.GetDataTimeStart() + static_cast(i)*fData.GetDataTimeStep(); // calculate theory for the given parameter set theo = N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg; - theo *= normalizer; - data = normalizer*fData.GetValue()->at(i); + data = fData.GetValue()->at(i); if (theo <= 0.0) { std::cerr << ">> PRunSingleHisto::CalcMaxLikelihood: **WARNING** NEGATIVE theory!!" << std::endl; @@ -389,7 +388,7 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector& par) } } - return 2.0*mllh; + return normalizer*2.0*mllh; } //-------------------------------------------------------------------------- @@ -472,9 +471,8 @@ Double_t PRunSingleHisto::CalcMaxLikelihoodExpected(const std::vector& time = fData.GetDataTimeStart() + static_cast(i)*fData.GetDataTimeStep(); // calculate theory for the given parameter set theo = N0*TMath::Exp(-time/tau)*(1.0+fTheory->Func(time, par, fFuncValues))+bkg; - theo *= normalizer; - data = normalizer*fData.GetValue()->at(i); + data = fData.GetValue()->at(i); if (theo <= 0.0) { std::cerr << ">> PRunSingleHisto::CalcMaxLikelihood: **WARNING** NEGATIVE theory!!" << std::endl; @@ -486,7 +484,7 @@ Double_t PRunSingleHisto::CalcMaxLikelihoodExpected(const std::vector& } } - return 2.0*mllh; + return normalizer*2.0*mllh; } //-------------------------------------------------------------------------- diff --git a/src/include/PFitter.h b/src/include/PFitter.h index 724c1f60..797e02d5 100644 --- a/src/include/PFitter.h +++ b/src/include/PFitter.h @@ -78,25 +78,31 @@ class PSectorChisq void SetSectorTime(Double_t last) { fLast = last; } void SetChisq(Double_t chisq) { fChisq = chisq; } void SetChisq(Double_t chisq, UInt_t idx); - void SetNDF(Double_t ndf) { fNDF = ndf; } - void SetNDF(Double_t ndf, UInt_t idx); + void SetExpectedChisq(Double_t expChisq) { fExpectedChisq = expChisq; } + void SetExpectedChisq(Double_t chisq, UInt_t idx); + void SetNDF(UInt_t ndf) { fNDF = ndf; } + void SetNDF(UInt_t ndf, UInt_t idx); Double_t GetTimeRangeFirst(UInt_t idx); Double_t GetTimeRangeLast() { return fLast; } Double_t GetChisq() { return fChisq; } - UInt_t GetNDF() { return fNDF; } - UInt_t GetNoRuns() { return fNoOfRuns; } Double_t GetChisq(UInt_t idx); + Double_t GetExpectedChisq() { return fExpectedChisq; } + Double_t GetExpectedChisq(UInt_t idx); + UInt_t GetNDF() { return fNDF; } UInt_t GetNDF(UInt_t idx); + UInt_t GetNoRuns() { return fNoOfRuns; } private: UInt_t fNoOfRuns; ///< number of runs presesent Double_t fLast; ///< requested time stamp Double_t fChisq; ///< chisq or maxLH for the sector + Double_t fExpectedChisq; ///< keep the expected chisq or maxLH for the sector UInt_t fNDF; ///< NDF for the sector - std::vector fFirst; ///< time stamp for fgb for a given run - std::vector fChisqRun; ///< chisq or maxLH for the sector and run - std::vector fNDFRun; ///< NDF for the sector and run + PDoubleVector fFirst; ///< time stamp for fgb for a given run + PDoubleVector fChisqRun; ///< chisq or maxLH for the sector and run + PDoubleVector fExpectedChisqRun; ///< expected chisq or maxLH for the sector and run + PUIntVector fNDFRun; ///< NDF for the sector and run }; //----------------------------------------------------------------------------- @@ -152,7 +158,7 @@ class PFitter PStringVector fElapsedTime; - Bool_t fSectorFlag; ///< sector command present flag + Bool_t fSectorFlag; ///< sector command present flag std::vector fSector; ///< stores all chisq/maxLH sector information // commands @@ -173,8 +179,8 @@ class PFitter Bool_t ExecuteScan(); Bool_t ExecuteSave(Bool_t first); Bool_t ExecuteSimplex(); - void PrepareSector(); - Bool_t ExecuteSector(); + void PrepareSector(PDoubleVector ¶m, PDoubleVector &error); + Bool_t ExecuteSector(std::ofstream &fout); Double_t MilliTime(); diff --git a/src/include/PRunListCollection.h b/src/include/PRunListCollection.h index 1ea14848..5bd49ccd 100644 --- a/src/include/PRunListCollection.h +++ b/src/include/PRunListCollection.h @@ -97,7 +97,7 @@ class PRunListCollection virtual Double_t GetMuMinusMaximumLikelihood(const std::vector& par) const; virtual Double_t GetNonMusrMaximumLikelihood(const std::vector& par) const; - virtual Double_t GetSingleHistoMaximumLikelihoodExpected(const std::vector& par, const UInt_t idx) const; + virtual Double_t GetSingleRunMaximumLikelihoodExpected(const std::vector& par, const UInt_t idx) const; virtual Double_t GetSingleRunMaximumLikelihood(const std::vector& par, const UInt_t idx) const; virtual UInt_t GetNoOfBinsFitted(const UInt_t idx) const;