LMU/musrfit
LMU/musrfit
5
0
Fork 0
Code Issues 6 Pull Requests Actions Packages Projects Releases 3 Activity

first full implementation of the sector command.

Browse Source 
Conflicts resolved:
	src/classes/PFitter.cpp
	src/classes/PFitterFcn.cpp
	src/classes/PRunListCollection.cpp
This commit is contained in:
suter_a 2020-02-03 20:54:23 +01:00
parent 1ea877621f
commit 3797c73aeb
6 changed files with 341 additions and 124 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

377
src/classes/PFitter.cpp
View File

@ -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; i<fNoOfRuns; i++) {
fFirst[i] = 0.0;
fChisqRun[i] = 0.0;
fExpectedChisqRun[i] = 0.0;
fNDFRun[i] = 0;
}
}
@ -138,6 +141,27 @@ void PSectorChisq::SetChisq(Double_t chisq, UInt_t idx)
fChisqRun[idx] = chisq;
}
//--------------------------------------------------------------------------
// SetExpectedChisq
//--------------------------------------------------------------------------
/**
* <p>Set 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
//--------------------------------------------------------------------------
/**
* <p>chisq/maxLH of the run with index idx
* <p>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
//--------------------------------------------------------------------------
/**
* <p>NDF of the run with index idx
* <p>Get expected chisq/maxLH of the run with index idx
*
* @param idx index of the run
*
* <b>return:</b> 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
//--------------------------------------------------------------------------
/**
* <p>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<Double_t> param = fMnUserParams.Params();
std::vector<Double_t> 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<Int_t>(fFitterFcnDKS->GetTotalNoOfFittedBins()) - usedParams;
ndf = fFitterFcnDKS->GetTotalNoOfFittedBins() - static_cast<UInt_t>(usedParams);
val = (*fFitterFcnDKS)(param);
} else {
ndf = static_cast<Int_t>(fFitterFcn->GetTotalNoOfFittedBins()) - usedParams;
ndf = fFitterFcn->GetTotalNoOfFittedBins() - static_cast<UInt_t>(usedParams);
val = (*fFitterFcn)(param);
}
if (fUseChi2) {
// calculate expected chisq
Double_t totalExpectedChisq = 0.0;
std::vector<Double_t> 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; i<expectedChisqPerRun.size(); i++) {
if (fDKSReady)
ndf_histo = fFitterFcnDKS->GetNoOfFittedBins(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; i<chisqPerRun.size(); i++) {
if (fDKSReady)
ndf_histo = fFitterFcnDKS->GetNoOfFittedBins(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<Int_t>(fFitterFcnDKS->GetTotalNoOfFittedBins()) - usedParams;
ndf = fFitterFcnDKS->GetTotalNoOfFittedBins() - static_cast<UInt_t>(usedParams);
val = (*fFitterFcnDKS)(param);
} else {
ndf = static_cast<Int_t>(fFitterFcn->GetTotalNoOfFittedBins()) - usedParams;
ndf = fFitterFcn->GetTotalNoOfFittedBins() - static_cast<UInt_t>(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; i<expectedChisqPerRun.size(); i++) {
if (fDKSReady)
ndf_histo = fFitterFcnDKS->GetNoOfFittedBins(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; i<chisqPerRun.size(); i++) {
if (fDKSReady)
ndf_histo = fFitterFcnDKS->GetNoOfFittedBins(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; i<expectedMaxLHPerRun.size(); i++) {
if (fDKSReady)
ndf_histo = fFitterFcnDKS->GetNoOfFittedBins(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; i<expectedMaxLHPerRun.size(); i++) {
if (fDKSReady)
ndf_histo = fFitterFcnDKS->GetNoOfFittedBins(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<Double_t> param;
Double_t totalExpectedChisq = 0.0;
std::vector<Double_t> expectedChisqPerHisto;
std::vector<Double_t> expectedchisqPerRun;
std::vector<UInt_t> ndfPerHisto;
for (UInt_t i=0; i<fParams.size(); i++)
@ -2075,54 +2112,63 @@ Bool_t PFitter::ExecuteSave(Bool_t firstSave)
// CalcExpectedChiSquare handles both, chisq and mlh
if (fDKSReady)
fFitterFcnDKS->CalcExpectedChiSquare(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<Double_t> chisqPerHisto;
std::vector<Double_t> chisqPerRun;
for (UInt_t i=0; i<fRunInfo->GetMsrRunList()->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; i<expectedChisqPerHisto.size(); i++) {
UInt_t ndf_run;
for (UInt_t i=0; i<expectedchisqPerRun.size(); i++) {
if (fDKSReady)
ndf_histo = fFitterFcnDKS->GetNoOfFittedBins(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<chisqPerHisto.size(); i++) {
} else if (chisqPerRun.size() > 1) { // in case expected chisq is not applicable like for asymmetry fits
UInt_t ndf_run = 0;
for (UInt_t i=0; i<chisqPerRun.size(); i++) {
if (fDKSReady)
ndf_histo = fFitterFcnDKS->GetNoOfFittedBins(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<fParams.size(); i++)
error.push_back(fParams[i].fStep);
PrepareSector(param, error);
}
// clean up
param.clear();
expectedChisqPerHisto.clear();
expectedchisqPerRun.clear();
ndfPerHisto.clear();
chisqPerHisto.clear();
chisqPerRun.clear();
std::cout << ">> 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()
//--------------------------------------------------------------------------
/**
* <p>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 &param, PDoubleVector &error)
{
// NOT YET IMPLEMENTED //as35
Double_t val;
UInt_t ndf;
/*
for (UInt_t i=0; i<fSector.size(); i++) {
std::cout << "debug> +++++" << std::endl;
std::cout << "debug> sector " << i << std::endl;
std::cout << "debug> noOfRuns: " << fSector[i].GetNoRuns() << std::endl;
for (UInt_t j=0; j<fSector[i].GetNoRuns(); j++) {
std::cout << "debug> ----" << 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; i<error.size(); i++) {
if (error[i] != 0.0)
usedParams++;
}
PDoublePairVector secFitRange;
secFitRange.resize(1);
if (fUseChi2) {
Double_t totalExpectedChisq = 0.0;
PDoubleVector expectedChisqPerRun;
PDoubleVector chisqPerRun;
for (UInt_t k=0; k<fSector.size(); k++) {
// set sector fit range
secFitRange[0].first = fSector[k].GetTimeRangeFirst(0);
secFitRange[0].second = fSector[k].GetTimeRangeLast();
fRunListCollection->SetFitRange(secFitRange);
// calculate chisq
val = (*fFitterFcn)(param);
fSector[k].SetChisq(val);
// calculate NDF
ndf = static_cast<UInt_t>(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; i<fRunInfo->GetMsrRunList()->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; i<expectedChisqPerRun.size(); i++) {
ndf_run = fFitterFcn->GetNoOfFittedBins(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; i<chisqPerRun.size(); i++) {
ndf_run = fFitterFcn->GetNoOfFittedBins(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; k<fSector.size(); k++) {
// set sector fit range
secFitRange[0].first = fSector[k].GetTimeRangeFirst(0);
secFitRange[0].second = fSector[k].GetTimeRangeLast();
fRunListCollection->SetFitRange(secFitRange);
// calculate maxLH
val = (*fFitterFcn)(param);
fSector[k].SetChisq(val);
// calculate NDF
ndf = static_cast<UInt_t>(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; i<fRunInfo->GetMsrRunList()->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; i<expectedMaxLHPerRun.size(); i++) {
ndf_run = fFitterFcn->GetNoOfFittedBins(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()
*
* <b>return:</b> 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<fSector.size(); i++) {
fout << std::endl;
fout << " Sector " << i+1 << ": FitRange: " << fSector[i].GetTimeRangeFirst(0) << ", " << fSector[i].GetTimeRangeLast() << std::endl;
fout << " chisq = " << fSector[i].GetChisq() << ", NDF = " << fSector[i].GetNDF() << ", chisq/NDF = " << fSector[i].GetChisq()/fSector[i].GetNDF();
if (fSector[i].GetExpectedChisq() > 0) {
fout << std::endl << " expected chisq = " << fSector[i].GetExpectedChisq() << ", NDF = " << fSector[i].GetNDF() << ", expected chisq/NDF = " << fSector[i].GetExpectedChisq()/fSector[i].GetNDF();
for (UInt_t j=0; j<fSector[i].GetNoRuns(); j++) {
fout << std::endl << " run block " << j+1 << " (NDF/red.chisq/red.chisq_e) = (" << fSector[i].GetNDF(j) << "/" << fSector[i].GetChisq(j)/fSector[i].GetNDF(j) << "/" << fSector[i].GetExpectedChisq(j)/fSector[i].GetNDF(j) << ")";
}
} else if (fSector[i].GetNoRuns() > 0) { // in case expected chisq is not applicable like for asymmetry fits
for (UInt_t j=0; j<fSector[i].GetNoRuns(); j++) {
fout << std::endl << " run block " << j+1 << " (NDF/red.chisq) = (" << fSector[i].GetNDF(j) << "/" << fSector[i].GetChisq(j)/fSector[i].GetNDF(j) << ")";
}
}
fout << std::endl << "++++";
}
} else { // max log likelihood
for (UInt_t i=0; i<fSector.size(); i++) {
fout << std::endl;
fout << " Sector " << i+1 << ": FitRange: " << fSector[i].GetTimeRangeFirst(0) << ", " << fSector[i].GetTimeRangeLast() << std::endl;
fout << " maxLH = " << fSector[i].GetChisq() << ", NDF = " << fSector[i].GetNDF() << ", maxLH/NDF = " << fSector[i].GetChisq()/fSector[i].GetNDF();
if (fSector[i].GetExpectedChisq() > 0) {
fout << std::endl << " expected maxLH = " << fSector[i].GetExpectedChisq() << ", NDF = " << fSector[i].GetNDF() << ", expected maxLH/NDF = " << fSector[i].GetExpectedChisq()/fSector[i].GetNDF();
for (UInt_t j=0; j<fSector[i].GetNoRuns(); j++) {
fout << std::endl << " run block " << j+1 << " (NDF/red.maxLH/red.maxLH_e) = (" << fSector[i].GetNDF(j) << "/" << fSector[i].GetChisq(j)/fSector[i].GetNDF(j) << "/" << fSector[i].GetExpectedChisq(j)/fSector[i].GetNDF(j) << ")";
}
} else if (fSector[i].GetNoRuns() > 0) { // in case expected chisq is not applicable like for asymmetry fits
for (UInt_t j=0; j<fSector[i].GetNoRuns(); j++) {
fout << std::endl << " run block " << j+1 << " (NDF/red.maxLH) = (" << fSector[i].GetNDF(j) << "/" << fSector[i].GetChisq(j)/fSector[i].GetNDF(j) << ")";
}
}
fout << std::endl << "++++";
}
}
return true;
}

4
src/classes/PFitterFcn.cpp
View File

@ -124,13 +124,13 @@ void PFitterFcn::CalcExpectedChiSquare(const std::vector<Double_t> &par, Double_
} else { // log max. likelihood
// single histo
for (UInt_t i=0; i<fRunListCollection->GetNoOfSingleHisto(); 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; i<fRunListCollection->GetNoOfMuMinus(); 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;
}

42
src/classes/PRunListCollection.cpp
View File

@ -399,20 +399,21 @@ Double_t PRunListCollection::GetSingleRunChisqExpected(const std::vector<Double_
Double_t expectedChisq = 0.0;
if (idx > 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; i<idx; i++) {
if (fMsrInfo->GetMsrRunList()->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; i<idx; i++) {
if (fMsrInfo->GetMsrRunList()->at(i).GetFitType() == type)
subIdx++;
}
}
// return the chisq of the single run
@ -658,36 +659,37 @@ Double_t PRunListCollection::GetNonMusrMaximumLikelihood(const std::vector<Doubl
}
//--------------------------------------------------------------------------
// GetSingleHistoMaximumLikelihoodExpected (public)
// GetSingleRunMaximumLikelihoodExpected (public)
//--------------------------------------------------------------------------
/**
* <p>Calculates expected mlh of the single histogram with run block index idx of a msr-file.
* <p>Calculates expected mlh of the run block index idx of a msr-file.
*
* <b>return:</b>
* - 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<Double_t>& par, const UInt_t idx) const
Double_t PRunListCollection::GetSingleRunMaximumLikelihoodExpected(const std::vector<Double_t>& 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; i<idx; i++) {
if (fMsrInfo->GetMsrRunList()->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; i<idx; i++) {
if (fMsrInfo->GetMsrRunList()->at(i).GetFitType() == type)
subIdx++;
}
}
// return the mlh of the single run

14
src/classes/PRunSingleHisto.cpp
View File

@ -276,9 +276,9 @@ Double_t PRunSingleHisto::CalcChiSquareExpected(const std::vector<Double_t>& 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<fEndTimeBin; ++i) {
time = fData.GetDataTimeStart() + static_cast<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;
@ -373,9 +373,8 @@ Double_t PRunSingleHisto::CalcMaxLikelihood(const std::vector<Double_t>& par)
time = fData.GetDataTimeStart() + static_cast<Double_t>(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<Double_t>& par)
}
}
return 2.0*mllh;
return normalizer*2.0*mllh;
}
//--------------------------------------------------------------------------
@ -472,9 +471,8 @@ Double_t PRunSingleHisto::CalcMaxLikelihoodExpected(const std::vector<Double_t>&
time = fData.GetDataTimeStart() + static_cast<Double_t>(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<Double_t>&
}
}
return 2.0*mllh;
return normalizer*2.0*mllh;
}
//--------------------------------------------------------------------------

26
src/include/PFitter.h
View File

@ -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<Double_t> fFirst; ///< time stamp for fgb for a given run
std::vector<Double_t> fChisqRun; ///< chisq or maxLH for the sector and run
std::vector<UInt_t> 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<PSectorChisq> 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 &param, PDoubleVector &error);
Bool_t ExecuteSector(std::ofstream &fout);
Double_t MilliTime();

2
src/include/PRunListCollection.h
View File

@ -97,7 +97,7 @@ class PRunListCollection
virtual Double_t GetMuMinusMaximumLikelihood(const std::vector<Double_t>& par) const;
virtual Double_t GetNonMusrMaximumLikelihood(const std::vector<Double_t>& par) const;
virtual Double_t GetSingleHistoMaximumLikelihoodExpected(const std::vector<Double_t>& par, const UInt_t idx) const;
virtual Double_t GetSingleRunMaximumLikelihoodExpected(const std::vector<Double_t>& par, const UInt_t idx) const;
virtual Double_t GetSingleRunMaximumLikelihood(const std::vector<Double_t>& par, const UInt_t idx) const;
virtual UInt_t GetNoOfBinsFitted(const UInt_t idx) const;