Run analyzer: Added possibility to define noise threshold cut manually and overwrite...

diff --git a/MABS_run_analyzer/ChargeSpektrum.c b/MABS_run_analyzer/ChargeSpektrum.c
index 5b470d0dca78c28fe9f10621a6b7df8617a3a013..40c1d839bf044d7a078012a084e4b084b75e29d7 100644 (file)
--- a/MABS_run_analyzer/ChargeSpektrum.c
+++ b/MABS_run_analyzer/ChargeSpektrum.c
@@ -132,6 +132,7 @@ void ChargeSpektrum(TString runnumber = "")
                 runs[runi]->useDynamicalNoise(true);
                 runs[runi]->useCommonModeFilter(false);
                 runs[runi]->setFixedThresholdValueElectrons(1440);
+                runs[runi]->setNoisethresholdborderADU(40);
                 
 //                 runs[runi]->analyzeFrame(57826);
 //                 runs[runi]->analyzeFrame(57827);

diff --git a/MABS_run_analyzer/HistogramType.c b/MABS_run_analyzer/HistogramType.c
index 78a8f70075b46c2a0921dda46be0059eace574da..148b1abf50a4997a48d1bf06a599f279e637ffdc 100644 (file)
--- a/MABS_run_analyzer/HistogramType.c
+++ b/MABS_run_analyzer/HistogramType.c
@@ -381,78 +381,82 @@ Bool_t HistogramType::calculteCCE(Bool_t verbose) {
     return 1;
 }
 
-Bool_t HistogramType::FindNoisethresholdborder(TH1F* histogrampointer, Bool_t verbose) {
-    Float_t posMaxValHist = histogrampointer->GetXaxis()->GetXmax();
-    
-    TH1F* smoothedcurce = (TH1F*)histogrampointer->Clone();
-    //     smoothedcurce->Smooth(4); // not working method, rescaling is easier and more reliable
-    Int_t rebinningfactor = 4;
-    smoothedcurce->RebinX(rebinningfactor);    
-    smoothedcurce->SetAxisRange(0,histogrampointer->GetBinCenter(histogrampointer->GetNbinsX()));
-    
-    if (verbose) {
-        cout << colorwhite << histogrampointer->GetName() << endlr;
-        TString canvastitle = Form("%s Noisethreshold border", histogrampointer->GetName());
-        TString canvasname =  Form("%s Noisethreshold border", histogrampointer->GetName());
-        TCanvas* canvas = new TCanvas(canvasname, canvastitle, 900, 700);
-        smoothedcurce->SetLineColor(color+1);
-        smoothedcurce->Draw();
-        histogrampointer->Draw("SAME");
-        canvas->Update();
-    }
-    
-    //     Int_t bini =smoothedcurce->GetMaximumBin();
-    //     thresholdbincurcandidate = bini;
-    //     if (verbose)
-    //         cout << "GetMaximumBin: smoothedcurce->GetBinContent(" << bini << "): " << smoothedcurce->GetBinContent(bini) << endl;
-    Int_t bini = 0;
-    Int_t falling = 0;
-    Int_t noisepeakcandidate = 0;
-    
-    do {
-        Float_t curval=smoothedcurce->GetBinContent(bini++);
-        if (curval*0.95 <= smoothedcurce->GetBinContent(bini))
-        {
-            falling = 0;
-            noisepeakcandidate = bini;
-        }
-        else
-        {
-            falling++;
-            if (verbose)
-                cout << "falling at " <<  smoothedcurce->GetXaxis()->GetBinCenter(bini) << "   as " << curval  << " < " << smoothedcurce->GetBinContent(bini) << endl; // debug
+Bool_t HistogramType::FindNoisethresholdborder(TH1F* histogrampointer, Bool_t overwrite, Bool_t verbose) {
+    if (overwrite || noisethresholdborder == -1)
+    {
+        Float_t posMaxValHist = histogrampointer->GetXaxis()->GetXmax();
+        
+        TH1F* smoothedcurce = (TH1F*)histogrampointer->Clone();
+        //     smoothedcurce->Smooth(4); // not working method, rescaling is easier and more reliable
+        Int_t rebinningfactor = 2;
+        smoothedcurce->RebinX(rebinningfactor);    
+        smoothedcurce->SetAxisRange(0,histogrampointer->GetBinCenter(histogrampointer->GetNbinsX()));
+        
+        if (verbose) {
+            cout << colorwhite << histogrampointer->GetName() << endlr;
+            TString canvastitle = Form("%s Noisethreshold border", histogrampointer->GetName());
+            TString canvasname =  Form("%s Noisethreshold border", histogrampointer->GetName());
+            TCanvas* canvas = new TCanvas(canvasname, canvastitle, 900, 700);
+            smoothedcurce->SetLineColor(color+1);
+            smoothedcurce->Draw();
+            histogrampointer->Draw("SAME");
+            canvas->Update();
         }
-    } while (falling < 2 && bini<smoothedcurce->GetNbinsX());        
-    if (verbose) {
-        cout << "Noise peak at: smoothedcurce->GetBinContent(" << noisepeakcandidate << "): " << smoothedcurce->GetBinContent(noisepeakcandidate) << endl;
-    }
-    Int_t thresholdbincurcandidate = 0;
-    bini = noisepeakcandidate;
-    
-    Int_t rising = 0;
-    do {
-        Float_t curval=smoothedcurce->GetBinContent(bini++);
-        if (curval*0.95 <= smoothedcurce->GetBinContent(bini))
-        {
-            rising++;
-            if (verbose)
-                cout << "rising at bin " << bini << ":" <<  smoothedcurce->GetXaxis()->GetBinCenter(bini) << "   as " << curval  << " < " << smoothedcurce->GetBinContent(bini) << endl; // debug
+        
+        //     Int_t bini =smoothedcurce->GetMaximumBin();
+        //     thresholdbincurcandidate = bini;
+        //     if (verbose)
+        //         cout << "GetMaximumBin: smoothedcurce->GetBinContent(" << bini << "): " << smoothedcurce->GetBinContent(bini) << endl;
+        Int_t bini = 0;
+        Int_t falling = 0;
+        Int_t noisepeakcandidate = 0;
+        
+        do {
+            Float_t curval=smoothedcurce->GetBinContent(bini++);
+            if (curval*0.95 <= smoothedcurce->GetBinContent(bini))
+            {
+                falling = 0;
+                noisepeakcandidate = bini;
+            }
+            else
+            {
+                falling++;
+                if (verbose)
+                    cout << "falling at " <<  smoothedcurce->GetXaxis()->GetBinCenter(bini) << "   as " << curval  << " < " << smoothedcurce->GetBinContent(bini) << endl; // debug
+            }
+        } while (falling < 2 && bini<smoothedcurce->GetNbinsX());        
+        if (verbose) {
+            cout << "Noise peak at: smoothedcurce->GetBinContent(" << noisepeakcandidate << "): " << smoothedcurce->GetBinContent(noisepeakcandidate) << endl;
         }
-        else
-        {
-            rising = 0;
-            thresholdbincurcandidate = bini;
-            if (verbose)
-                cout << "falling at bin " << bini << ":" <<  smoothedcurce->GetXaxis()->GetBinCenter(bini) << "   as " << curval  << " < " << smoothedcurce->GetBinContent(bini) << endl; // debug
+        Int_t thresholdbincurcandidate = 0;
+        bini = noisepeakcandidate;
+        
+        Int_t rising = 0;
+        do {
+            Float_t curval=smoothedcurce->GetBinContent(bini++);
+            if (curval*0.95 <= smoothedcurce->GetBinContent(bini))
+            {
+                rising++;
+                if (verbose)
+                    cout << "rising at bin " << bini << ":" <<  smoothedcurce->GetXaxis()->GetBinCenter(bini) << "   as " << curval  << " < " << smoothedcurce->GetBinContent(bini) << endl; // debug
+            }
+            else
+            {
+                rising = 0;
+                thresholdbincurcandidate = bini;
+                if (verbose)
+                    cout << "falling at bin " << bini << ":" <<  smoothedcurce->GetXaxis()->GetBinCenter(bini) << "   as " << curval  << " < " << smoothedcurce->GetBinContent(bini) << endl; // debug
+            }
+        } while (rising < 3 && bini<smoothedcurce->GetNbinsX());
+        thresholdbincurcandidate *= rebinningfactor;
+        
+        noisethresholdborder = histogrampointer->GetXaxis()->GetBinUpEdge(thresholdbincurcandidate);
+        if (verbose) {
+            cout << "     Noise threshold at " << noisethresholdborder << endl;
         }
-    } while (rising < 3 && bini<smoothedcurce->GetNbinsX());
-    thresholdbincurcandidate *= rebinningfactor;
-    
-    noisethresholdborder = histogrampointer->GetXaxis()->GetBinUpEdge(thresholdbincurcandidate);
-    if (verbose) {
-        cout << "     Noise threshold at " << noisethresholdborder << endl;
+        return 0;
     }
-    return 0;
+    return 1;
 }
 
 Bool_t HistogramType::integrateSr90Spectra(TH1F* histogrampointer, Int_t frames_found, Float_t thresholdborder, Bool_t verbose) {
@@ -461,7 +465,7 @@ Bool_t HistogramType::integrateSr90Spectra(TH1F* histogrampointer, Int_t frames_
     if (thresholdborder < 0)
     {
         if (noisethresholdborder < 0) {
-            FindNoisethresholdborder(histogrampointer, verbose);
+            FindNoisethresholdborder(histogrampointer, true, verbose);
         }
         thresholdbincurcandidate = histogrampointer->GetXaxis()->FindBin(noisethresholdborder);
     }

diff --git a/MABS_run_analyzer/HistogramType.h b/MABS_run_analyzer/HistogramType.h
index 76c7a1704f9d9b2b05e266ce0d628629f1339f75..43e310bd6b1903a2e17cfe9e72d783393efebea5 100644 (file)
--- a/MABS_run_analyzer/HistogramType.h
+++ b/MABS_run_analyzer/HistogramType.h
@@ -216,7 +216,7 @@ public:
      * @return true if succesfull
      * 
      */
-    Bool_t FindNoisethresholdborder(TH1F* histogrampointer, Bool_t verbose=false);
+    Bool_t FindNoisethresholdborder(TH1F* histogrampointer, Bool_t overwrite=false, Bool_t verbose=false);
     
     /**
      * @brief calculates StoN ratio, if beta source was used

diff --git a/MABS_run_analyzer/Run.c b/MABS_run_analyzer/Run.c
index f1d67cb07eefa912583f7ce831022ccacd8784bb..f9e121b7f08b63cab63d98a0a36ddbbb025ea65c 100644 (file)
--- a/MABS_run_analyzer/Run.c
+++ b/MABS_run_analyzer/Run.c
@@ -380,7 +380,7 @@ Bool_t Run::analyzeRun(Bool_t force)
         }
         else
         {
-            cout << "Skipped analysis of run " << labbook.runnumber << ", I think the root file for this run allready exists." << endl;
+            cout << endl <<"Skipped analysis of run " << labbook.runnumber << ", I think the root file for this run allready exists." << endl;
             cout << colorwhite << "initOldRootFile():" << endlr;
         }
         if (processed->initOldRootFile()) return 1;
@@ -394,7 +394,16 @@ Bool_t Run::analyzeRun(Bool_t force)
         cout << "---------- loop over all classes --------" << endl;
         for (vector<HistogramType*>::iterator curHistogramClass = HistogramClassVector.begin(); curHistogramClass != HistogramClassVector.end(); curHistogramClass++)  {
             cout << "Processing histograms in class: <" << colorwhite << (*curHistogramClass)->histogramdescription << colorreset << " >" << endlr;
-            binNoise((*curHistogramClass));     
+            binNoise((*curHistogramClass));                 
+            if (noisethresholdborderADU > -1)
+                (*curHistogramClass)->noisethresholdborder = noisethresholdborderADU;   
+            if (noisethresholdborderE > -1) {
+                if (labbook.gainDB > 0) {
+                    (*curHistogramClass)->noisethresholdborder = noisethresholdborderE / labbook.gainDB;       
+                } else {
+                    cout << colorred << "Could not set noise threshold border in units of electrons to " << noisethresholdborderE << "e, no calibration done yet. Please rerun after first analysis." <<endlr;
+                }
+            }
             if (labbook.source.Contains("Fe")||labbook.source.Contains("Cd")) {
                 cout << colorwhite << " calculateCCE():" << endlr;
                 (*curHistogramClass)->calculteCCE(true);
@@ -407,10 +416,11 @@ Bool_t Run::analyzeRun(Bool_t force)
             }
         }
         cout << "---------------------------------" << endl;
-        cout << colorwhite << "integrateSr90Spectra():" << endlr;
-//         histogramthreshold->noisethresholdborder = 35.0; // TODO remove the overwriting of automatic noise threshold setting
-        histogramthreshold->integrateSr90Spectra(histogramthreshold->Seed, frames_found, -1, true);
-        histogramfixedthreshold->integrateSr90Spectra(histogramfixedthreshold->Seed, frames_found, 0);
+        if (labbook.source.Contains("Sr90")) {
+            cout << colorwhite << "integrateSr90Spectra():" << endlr;
+            histogramthreshold->integrateSr90Spectra(histogramthreshold->Seed, frames_found, -1, true);
+            histogramfixedthreshold->integrateSr90Spectra(histogramfixedthreshold->Seed, frames_found, 0);
+        }
         rescaleHistogramClasses();
         cout << colorwhite << "updateDatabase():" << endlr;
         updateDatabase();
@@ -422,6 +432,20 @@ Bool_t Run::analyzeRun(Bool_t force)
     return 1;
 }
 
+Bool_t Run::setNoisethresholdborderADU(Float_t noisethresholdborder)
+{
+    noisethresholdborderADU = noisethresholdborder;
+    cout << " Noise threshold border: " << colorwhite << noisethresholdborder  << " ADU"<< endlr;
+    return 0;
+}
+
+Bool_t Run::setNoisethresholdborderE(Float_t noisethresholdborder)
+{
+    noisethresholdborderADU = noisethresholdborder;
+    cout << " Noise threshold border: " << colorwhite << noisethresholdborder  << " e"<< colorreset << "  <-- only used if GAIN is allready saved in database, rerun after first analysis to make sure" << endl;
+    return 0;
+}
+
 Bool_t Run::rescaleHistogramClasses()
 {
     float_t vetopeakposition = -1;
@@ -845,7 +869,7 @@ Bool_t Run::binSeedSumVeto()
 
     for (vector<HistogramType*>::iterator curHistogramClass = HistogramClassVector.begin(); curHistogramClass != HistogramClassVector.end(); curHistogramClass++)  {
         Double_t* parameters = (Double_t *)calloc(3, sizeof(Double_t));
-        (*curHistogramClass)->FindNoisethresholdborder((*curHistogramClass)->Seed, false);
+        (*curHistogramClass)->FindNoisethresholdborder((*curHistogramClass)->Seed, false, false);
         if (labbook.source.Contains("Fe")||labbook.source.Contains("Cd"))
             parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->Veto, "gaus", true);
         (*curHistogramClass)->posVeto = parameters[1];

diff --git a/MABS_run_analyzer/Run.h b/MABS_run_analyzer/Run.h
index cba89f500cdf4f4423068dc87aef2bc1106e16e1..f02ccd053b8819407424e493c89ee15bddc8a2ad 100644 (file)
--- a/MABS_run_analyzer/Run.h
+++ b/MABS_run_analyzer/Run.h
@@ -205,6 +205,12 @@ private:
     Bool_t dynamicalNoise = 1;
     /** @brief set this variable to true, to use the common mode noise filter in MAPS:analyzeRun() */
     Bool_t commonModeFilter = 1;
+        
+    /** @brief overwrite value from which the integral should be taken, overwrite noise cut in electrons */
+    Float_t noisethresholdborderE = -1;
+    
+    /** @brief overwrite value from which the integral should be taken, overwrite noise cut in ADU */
+    Float_t noisethresholdborderADU = -1;
     
     Bool_t rescaleHistogramClasses();
        
@@ -305,6 +311,16 @@ public:
      */
     Bool_t setDynamicalNoiseMode(TString mode);
     
+    /**
+     * @brief set noise threshold border in ADU
+     */
+    Bool_t setNoisethresholdborderADU(Float_t noisethresholdborder);
+        
+    /**
+     * @brief set noise threshold border in e
+     */
+    Bool_t setNoisethresholdborderE(Float_t noisethresholdborder);
+    
     /// stores information from the SQL database of a given run
     labbooksctruct labbook;