From 6128464d5bdfae16600dbe0d3f816bb7422d6c4c Mon Sep 17 00:00:00 2001 From: nemu Date: Wed, 30 Apr 2008 11:46:32 +0000 Subject: [PATCH] various ways to handle bkg for single histo fits --- src/classes/PMsrHandler.cpp | 45 +++++++-- src/classes/PRunSingleHisto.cpp | 157 ++++++++++++++++++++++++++++++-- src/include/PRunAsymmetry.h | 8 +- src/include/PRunSingleHisto.h | 3 + 4 files changed, 193 insertions(+), 20 deletions(-) diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index cf7a998c..f06fd76e 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -510,9 +510,11 @@ int PMsrHandler::WriteMsrLogFile() f.width(15); f << endl << left << "backgr.fix"; for (unsigned int j=0; j<2; j++) { - f.precision(prec); - f.width(12); - f << left << fRuns[i].fBkgFix[j]; + if (!isnan(fRuns[i].fBkgFix[j])) { + f.precision(prec); + f.width(12); + f << left << fRuns[i].fBkgFix[j]; + } } CheckAndWriteComment(f, ++lineNo); } @@ -1053,10 +1055,11 @@ bool PMsrHandler::HandleRunEntry(PMsrLines &lines) // RUN line ---------------------------------------------- if (iter->fLine.BeginsWith("run", TString::kIgnoreCase)) { - if (!first) // not the first run in the list + if (!first) { // not the first run in the list fRuns.push_back(param); - else + } else { first = false; + } InitRunParameterStructure(param); @@ -1252,10 +1255,10 @@ bool PMsrHandler::HandleRunEntry(PMsrLines &lines) // backgr.fix ---------------------------------------------- if (iter->fLine.BeginsWith("backgr.fix", TString::kIgnoreCase)) { - if (tokens->GetEntries() != 3) { + if ((tokens->GetEntries() != 2) && (tokens->GetEntries() != 3)) { error = true; } else { - for (int i=1; i<3; i++) { + for (int i=1; iGetEntries(); i++) { ostr = dynamic_cast(tokens->At(i)); str = ostr->GetString(); if (str.IsFloat()) @@ -1439,13 +1442,39 @@ bool PMsrHandler::HandleRunEntry(PMsrLines &lines) if (error) { --iter; - cout << endl << "ERROR in line " << iter->fLineNo << ":"; + cout << endl << "**ERROR** in line " << iter->fLineNo << ":"; cout << endl << iter->fLine.Data(); cout << endl << "RUN block syntax is too complex to print it here. Please check the manual."; } else { // save last run found fRuns.push_back(param); } + // check if for fittypes: single histo, asymmetry, RRF any background info is given + for (unsigned int i=0; i= 0) { // check if backgr.fit is given + found = true; + } else if (!isnan(fRuns[i].fBkgFix[0])) { // check if backgr.fix is given + found = true; + } else if (fRuns[i].fBkgRange[0] >= 0) { // check if background window is given + found = true; + } else { + found = false; + } + if (!found) { + cout << endl << "**ERROR** for run " << fRuns[i].fRunName << ", forward " << fRuns[i].fForwardHistoNo; + cout << endl << " no background information found!"; + cout << endl << " Either of the tags 'backgr.fit', 'backgr.fix', 'background'"; + cout << endl << " with data is needed."; + cout << endl; + return false; + } + } + } + return !error; } diff --git a/src/classes/PRunSingleHisto.cpp b/src/classes/PRunSingleHisto.cpp index 43de7b2d..c1567c70 100644 --- a/src/classes/PRunSingleHisto.cpp +++ b/src/classes/PRunSingleHisto.cpp @@ -117,7 +117,16 @@ double PRunSingleHisto::CalcChiSquare(const std::vector& par) tau = PMUON_LIFETIME; // get background - double bkg = par[fRunInfo->fBkgFitParamNo-1]; + double bkg; + if (fRunInfo->fBkgFitParamNo == -1) { // bkg not fitted + if (isnan(fRunInfo->fBkgFix[0])) { // no fixed background given (background interval) + bkg = fBackground; + } else { // fixed bkg given + bkg = fRunInfo->fBkgFix[0]; + } + } else { // bkg fitted + bkg = par[fRunInfo->fBkgFitParamNo-1]; + } // calculate functions for (int i=0; iGetNoOfFuncs(); i++) { @@ -137,6 +146,8 @@ double PRunSingleHisto::CalcChiSquare(const std::vector& par) } } +//cout << endl << "chisq=" << chisq*fRunInfo->fPacking; + return chisq; } @@ -172,7 +183,16 @@ double PRunSingleHisto::CalcMaxLikelihood(const std::vector& par) tau = PMUON_LIFETIME; // get background - double bkg = par[fRunInfo->fBkgFitParamNo-1]; + double bkg; + if (fRunInfo->fBkgFitParamNo == -1) { // bkg not fitted + if (isnan(fRunInfo->fBkgFix[0])) { // no fixed background given (background interval) + bkg = fBackground; + } else { // fixed bkg given + bkg = fRunInfo->fBkgFix[0]; + } + } else { // bkg fitted + bkg = par[fRunInfo->fBkgFitParamNo-1]; + } // calculate functions for (int i=0; iGetNoOfFuncs(); i++) { @@ -237,7 +257,16 @@ void PRunSingleHisto::CalcTheory() tau = PMUON_LIFETIME; // get background - double bkg = par[fRunInfo->fBkgFitParamNo-1]; + double bkg; + if (fRunInfo->fBkgFitParamNo == -1) { // bkg not fitted + if (isnan(fRunInfo->fBkgFix[0])) { // no fixed background given (background interval) + bkg = fBackground; + } else { // fixed bkg given + bkg = fRunInfo->fBkgFix[0]; + } + } else { // bkg fitted + bkg = par[fRunInfo->fBkgFitParamNo-1]; + } // calculate functions for (int i=0; iGetNoOfFuncs(); i++) { @@ -380,8 +409,23 @@ bool PRunSingleHisto::PrepareFitData() return false; } + // check how the background shall be handled + if (fRunInfo->fBkgFitParamNo == -1) { // bkg shall **NOT** be fitted + // subtract background from histogramms ------------------------------------------ + if (isnan(fRunInfo->fBkgFix[0])) { // no fixed background given + if (fRunInfo->fBkgRange[0] != 0) { + if (!EstimateBkg(histoNo)) + return false; + } else { // no background given to do the job + cout << endl << "PRunSingleHisto::PrepareData(): Neither fix background nor background bins are given!"; + cout << endl << "One of the two is needed! Will quit ..."; + return false; + } + } + } + // everything looks fine, hence fill data set - double value = 0.0; + double value = 0.0; // data start at data_start-t0 // time shifted so that packing is included correctly, i.e. t0 == t0 after packing fData.fDataTimeStart = fTimeResolution*(((double)start-t0)+(double)fRunInfo->fPacking/2.0); @@ -396,6 +440,7 @@ bool PRunSingleHisto::PrepareFitData() fData.fError.push_back(1.0); else fData.fError.push_back(TMath::Sqrt(value/fRunInfo->fPacking)); + // reset values value = 0.0; } value += runData->fDataBin[histoNo][i]; @@ -509,7 +554,7 @@ bool PRunSingleHisto::PrepareRawViewData() } // everything looks fine, hence fill data set - double value = 0.0; + double value = 0.0; // data start at data_start-t0 // time shifted so that packing is included correctly, i.e. t0 == t0 after packing fData.fDataTimeStart = fTimeResolution*(((double)start-t0)+(double)fRunInfo->fPacking/2.0); @@ -524,6 +569,7 @@ bool PRunSingleHisto::PrepareRawViewData() fData.fError.push_back(1.0); else fData.fError.push_back(TMath::Sqrt(value/fRunInfo->fPacking)); + // reset values value = 0.0; } value += runData->fDataBin[histoNo][i]; @@ -565,7 +611,18 @@ bool PRunSingleHisto::PrepareRawViewData() tau = PMUON_LIFETIME; // get background - double bkg = par[fRunInfo->fBkgFitParamNo-1]; + double bkg; + if (fRunInfo->fBkgFitParamNo == -1) { // bkg not fitted + if (isnan(fRunInfo->fBkgFix[0])) { // no fixed background given (background interval) + if (!EstimateBkg(histoNo)) + return false; + bkg = fBackground; + } else { // fixed bkg given + bkg = fRunInfo->fBkgFix[0]; + } + } else { // bkg fitted + bkg = par[fRunInfo->fBkgFitParamNo-1]; + } // calculate functions for (int i=0; iGetNoOfFuncs(); i++) { @@ -720,8 +777,18 @@ bool PRunSingleHisto::PrepareViewData() //cout << endl << ">> tau = " << tau; // get background - double bkg = par[fRunInfo->fBkgFitParamNo-1]; -//cout << endl << ">> bkg = " << bkg; + double bkg; + if (fRunInfo->fBkgFitParamNo == -1) { // bkg not fitted + if (isnan(fRunInfo->fBkgFix[0])) { // no fixed background given (background interval) + if (!EstimateBkg(histoNo)) + return false; + bkg = fBackground; + } else { // fixed bkg given + bkg = fRunInfo->fBkgFix[0]; + } + } else { // bkg fitted + bkg = par[fRunInfo->fBkgFitParamNo-1]; + } double value = 0.0; double expval; @@ -774,3 +841,77 @@ bool PRunSingleHisto::PrepareViewData() return true; } + +//-------------------------------------------------------------------------- +// EstimatBkg +//-------------------------------------------------------------------------- +/** + *

+ */ +bool PRunSingleHisto::EstimateBkg(unsigned int histoNo) +{ + double beamPeriod = 0.0; + + // check if data are from PSI, RAL, or TRIUMF + if (fRunInfo->fInstitute.Contains("psi")) + beamPeriod = ACCEL_PERIOD_PSI; + else if (fRunInfo->fInstitute.Contains("ral")) + beamPeriod = ACCEL_PERIOD_RAL; + else if (fRunInfo->fInstitute.Contains("triumf")) + beamPeriod = ACCEL_PERIOD_TRIUMF; + else + beamPeriod = 0.0; + + // check if start and end are in proper order + unsigned int start = fRunInfo->fBkgRange[0]; + unsigned int end = fRunInfo->fBkgRange[1]; + if (end < start) { + cout << endl << "PRunSingleHisto::EstimatBkg(): end = " << end << " > start = " << start << "! Will swap them!"; + unsigned int keep = end; + end = start; + start = keep; + } + + // calculate proper background range + if (beamPeriod != 0.0) { + double beamPeriodBins = beamPeriod/fRunInfo->fPacking; + unsigned int periods = (unsigned int)((double)(end - start + 1) / beamPeriodBins); + end = start + (unsigned int)round((double)periods*beamPeriodBins); + cout << endl << "PRunSingleHisto::EstimatBkg(): Background " << start << ", " << end; + if (end == start) + end = fRunInfo->fBkgRange[1]; + } + + // get the proper run + PRawRunData* runData = fRawData->GetRunData(fRunInfo->fRunName); + + // check if start is within histogram bounds + if ((start < 0) || (start >= runData->fDataBin[histoNo].size())) { + cout << endl << "PRunSingleHisto::EstimatBkg(): background bin values out of bound!"; + cout << endl << " histo lengths = " << runData->fDataBin[histoNo].size(); + cout << endl << " background start = " << start; + return false; + } + + // check if end is within histogram bounds + if ((end < 0) || (end >= runData->fDataBin[histoNo].size())) { + cout << endl << "PRunSingleHisto::EstimatBkg(): background bin values out of bound!"; + cout << endl << " histo lengths = " << runData->fDataBin[histoNo].size(); + cout << endl << " background end = " << end; + return false; + } + + // calculate background + double bkg = 0.0; + + // forward + for (unsigned int i=start; ifDataBin[histoNo][i]; + bkg /= static_cast(end - start + 1); + + fBackground = bkg; // keep background (per bin) for chisq, max.log.likelihood + + cout << endl << ">> fRunInfo->fRunName=" << fRunInfo->fRunName.Data() << ", histNo=" << histoNo << ", bkg=" << bkg; + + return true; +} diff --git a/src/include/PRunAsymmetry.h b/src/include/PRunAsymmetry.h index 15d81aa5..43df615b 100644 --- a/src/include/PRunAsymmetry.h +++ b/src/include/PRunAsymmetry.h @@ -59,10 +59,10 @@ class PRunAsymmetry : public PRunBase double fFitStopTime; unsigned int fNoOfFitBins; - vector fForward; ///< forward histo data - vector fForwardErr; ///< forward histo errors - vector fBackward; ///< backward histo data - vector fBackwardErr; ///< backward histo errors + PDoubleVector fForward; ///< forward histo data + PDoubleVector fForwardErr; ///< forward histo errors + PDoubleVector fBackward; ///< backward histo data + PDoubleVector fBackwardErr; ///< backward histo errors bool SubtractFixBkg(); bool SubtractEstimatedBkg(); diff --git a/src/include/PRunSingleHisto.h b/src/include/PRunSingleHisto.h index 4752e119..4f5320d5 100644 --- a/src/include/PRunSingleHisto.h +++ b/src/include/PRunSingleHisto.h @@ -57,6 +57,9 @@ class PRunSingleHisto : public PRunBase double fFitStartTime; double fFitStopTime; unsigned int fNoOfFitBins; + double fBackground; ///< needed if background range is given. In units per bin + + bool EstimateBkg(unsigned int histoNo); }; #endif // _PRUNSINGLEHISTO_H_