From 35198e46da780674955524032c5c44368f35138d Mon Sep 17 00:00:00 2001
From: Suter Andreas <andreas.suter@psi.ch>
Date: Fri, 6 May 2016 16:58:43 +0200
Subject: [PATCH 1/2] marginal improved error estimate for single histo RRF

---
 src/classes/PRunSingleHistoRRF.cpp | 14 ++++++++++----
 src/include/PRunSingleHistoRRF.h   |  1 +
 2 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/src/classes/PRunSingleHistoRRF.cpp b/src/classes/PRunSingleHistoRRF.cpp
index 547e6442..707b6643 100644
--- a/src/classes/PRunSingleHistoRRF.cpp
+++ b/src/classes/PRunSingleHistoRRF.cpp
@@ -58,7 +58,8 @@ using namespace std;
 PRunSingleHistoRRF::PRunSingleHistoRRF() : PRunBase()
 {
   fNoOfFitBins  = 0;
-  fBackground = 0;
+  fBackground = 0.0;
+  fBkgErr = 1.0;
   fRRFPacking = -1;
 
   // the 2 following variables are need in case fit range is given in bins, and since
@@ -612,7 +613,7 @@ Bool_t PRunSingleHistoRRF::PrepareFitData(PRawRunData* runData, const UInt_t his
     exp_t_tau = exp(time_tau);
     fForward[i] *= exp_t_tau;
     fM.push_back(fForward[i]);   // i.e. M(t) = [N(t)-Nbkg] exp(+t/tau); needed to estimate N0 later on
-    fMerr.push_back(exp_t_tau*sqrt(rawNt[i]-fBackground));
+    fMerr.push_back(exp_t_tau*sqrt(rawNt[i]+fBkgErr*fBkgErr));
   }
 
   // calculate weights
@@ -620,7 +621,7 @@ Bool_t PRunSingleHistoRRF::PrepareFitData(PRawRunData* runData, const UInt_t his
     if (fMerr[i] > 0.0)
       fW.push_back(1.0/(fMerr[i]*fMerr[i]));
     else
-      fW.push_back(0.0);
+      fW.push_back(1.0);
   }
   // now fForward = exp(+t/tau) [N(t)-Nbkg] = M(t)
 
@@ -1188,7 +1189,12 @@ Bool_t PRunSingleHistoRRF::EstimateBkg(UInt_t histoNo)
 
   fBackground = bkg;  // keep background (per bin)
 
-  cout << endl << "info> fBackground=" << fBackground << endl;
+  bkg = 0.0;
+  for (UInt_t i=start; i<end; i++)
+    bkg += pow(fForward[i]-fBackground, 2.0);
+  fBkgErr = sqrt(bkg/(static_cast<Double_t>(end - start)));
+
+  cout << endl << "info> fBackground=" << fBackground << "(" << fBkgErr << ")" << endl;
 
   fRunInfo->SetBkgEstimated(fBackground, 0);
 
diff --git a/src/include/PRunSingleHistoRRF.h b/src/include/PRunSingleHistoRRF.h
index 670b90db..44e7b9a5 100644
--- a/src/include/PRunSingleHistoRRF.h
+++ b/src/include/PRunSingleHistoRRF.h
@@ -65,6 +65,7 @@ class PRunSingleHistoRRF : public PRunBase
 
     UInt_t fNoOfFitBins;    ///< number of bins to be fitted
     Double_t fBackground;   ///< needed if background range is given (units: 1/bin)
+    Double_t fBkgErr;       ///< estimate error on the estimated background
     Int_t fRRFPacking;      ///< RRF packing for this particular run. Given in the GLOBAL-block.
 
     Int_t fGoodBins[2];     ///< keep first/last good bins. 0=fgb, 1=lgb

From 2f91278871e2e81e9fa177d0e9814896222d46de Mon Sep 17 00:00:00 2001
From: Suter Andreas <andreas.suter@psi.ch>
Date: Fri, 6 May 2016 18:54:36 +0200
Subject: [PATCH 2/2] fixed a minor bug in the background handling. For now
 change over to un-weighted N0 estimate which seems to be slightly better.
 This only could mean, that the error estimate if M is NOT optimal. Eventually
 this needs some deeper investigation.

---
 src/classes/PRunSingleHistoRRF.cpp | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/src/classes/PRunSingleHistoRRF.cpp b/src/classes/PRunSingleHistoRRF.cpp
index 707b6643..d43546ae 100644
--- a/src/classes/PRunSingleHistoRRF.cpp
+++ b/src/classes/PRunSingleHistoRRF.cpp
@@ -593,6 +593,9 @@ Bool_t PRunSingleHistoRRF::PrepareFitData(PRawRunData* runData, const UInt_t his
       if (!EstimateBkg(histoNo))
         return false;
     }
+    // subtract background from fForward
+    for (UInt_t i=0; i<fForward.size(); i++)
+      fForward[i] -= fBackground;
   } else { // fixed background given
     for (UInt_t i=0; i<fForward.size(); i++) {
       fForward[i] -= fRunInfo->GetBkgFix(0);
@@ -1093,15 +1096,17 @@ Double_t PRunSingleHistoRRF::EstimateN0(Double_t &errN0, Double_t freqMax)
 {
   // endBin is estimated such that the number of full cycles (according to the maximum frequency of the data)
   // is approximately the time fN0EstimateEndTime.
-  Int_t endBin = (Int_t)round(fN0EstimateEndTime / fTimeResolution * ceil(freqMax)/freqMax);
+  Int_t endBin = (Int_t)round(ceil(fN0EstimateEndTime*freqMax/TMath::TwoPi()) * (TMath::TwoPi()/freqMax) / fTimeResolution);
 
   Double_t n0 = 0.0;
   Double_t wN = 0.0;
   for (Int_t i=0; i<endBin; i++) {
-    n0 += fW[i]*fM[i];
+//    n0 += fW[i]*fM[i];
+    n0 += fM[i];
     wN += fW[i];
   }
-  n0 /= wN;
+//  n0 /= wN;
+  n0 /= endBin;
 
   errN0 = 0.0;
   for (Int_t i=0; i<endBin; i++) {