various ways to handle bkg for single histo fits

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);
+        if (!isnan(fRuns[i].fBkgFix[j])) {
-        f.width(12);
+          f.precision(prec);
-        f << left << fRuns[i].fBkgFix[j];
+          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; i<tokens->GetEntries(); i++) {
           ostr = dynamic_cast<TObjString*>(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<fRuns.size(); i++) {
+    if ((fRuns[i].fFitType == MSR_FITTYPE_SINGLE_HISTO) ||
+        (fRuns[i].fFitType == MSR_FITTYPE_ASYM) ||
+        (fRuns[i].fFitType == MSR_FITTYPE_ASYM_RRF)) {
+      bool found;
+      if (fRuns[i].fBkgFitParamNo >= 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;
 }
 

src/classes/PRunSingleHisto.cpp

@@ -117,7 +117,16 @@ double PRunSingleHisto::CalcChiSquare(const std::vector<double>& 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; i<fMsrInfo->GetNoOfFuncs(); i++) {
@@ -137,6 +146,8 @@ double PRunSingleHisto::CalcChiSquare(const std::vector<double>& par)
     }
   }
 
+//cout << endl << "chisq=" << chisq*fRunInfo->fPacking;
+
   return chisq;
 }
 
@@ -172,7 +183,16 @@ double PRunSingleHisto::CalcMaxLikelihood(const std::vector<double>& 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; i<fMsrInfo->GetNoOfFuncs(); 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; i<fMsrInfo->GetNoOfFuncs(); 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; i<fMsrInfo->GetNoOfFuncs(); i++) {
@@ -720,8 +777,18 @@ bool PRunSingleHisto::PrepareViewData()
 //cout << endl << ">> tau = " << tau;
 
   // get background
-  double bkg = par[fRunInfo->fBkgFitParamNo-1];
+  double bkg;
-//cout << endl << ">> bkg = " << 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
+//--------------------------------------------------------------------------
+/**
+ * <p>
+ */
+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; i<end; i++)
+    bkg += runData->fDataBin[histoNo][i];
+  bkg /= static_cast<double>(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;
+}

src/include/PRunAsymmetry.h

@@ -59,10 +59,10 @@ class PRunAsymmetry : public PRunBase
     double fFitStopTime;
     unsigned int fNoOfFitBins;
 
-    vector<double> fForward;     ///< forward histo data
+    PDoubleVector fForward;     ///< forward histo data
-    vector<double> fForwardErr;  ///< forward histo errors
+    PDoubleVector fForwardErr;  ///< forward histo errors
-    vector<double> fBackward;    ///< backward histo data
+    PDoubleVector fBackward;    ///< backward histo data
-    vector<double> fBackwardErr; ///< backward histo errors
+    PDoubleVector fBackwardErr; ///< backward histo errors
 
     bool SubtractFixBkg();
     bool SubtractEstimatedBkg();

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_