Analysis: Improved Veto peak finding routine, added user warnings

diff --git a/MABS_run_analyzer/Run.c b/MABS_run_analyzer/Run.c
index a4b3e975f069372aa80c91e58d7f10bc210f3f0e..77de290dbaabfbe68f817d0b5ba45bd61ad9020d 100644 (file)
--- a/MABS_run_analyzer/Run.c
+++ b/MABS_run_analyzer/Run.c
@@ -211,7 +211,7 @@ Bool_t Run::analyzeRun(Bool_t force)
         }
         else
         {
-            cout << "\033[1;33mSkipped analysis of run " << labbook.runnumber << ", I think the root file for this run allready exists.\033[0m" << endl;
+            cout << "Skipped analysis of run " << labbook.runnumber << ", I think the root file for this run allready exists." << endl;
             cout << colorwhite << "initOldRootFile():" << endlr;
             error = processed->initOldRootFile();
         }
@@ -762,59 +762,72 @@ Bool_t Run::plotAllHistogramsCalibrated()
     return 1;
 }
 
-Float_t Run::FitPerform(TH1F* histogrampointer, TString fitFuncType)
+Float_t Run::FitPerform(TH1F* histogrampointer, TString fitFuncType, Bool_t verbose)
 {
     Float_t posMax = 0;
     Float_t posMaxValHist = histogrampointer->GetXaxis()->GetXmax();
+    Float_t noiseborder = posMaxValHist/15;
     
     if (doFits) 
     {        
-        histogrampointer->GetXaxis()->SetRange(histogrampointer->GetXaxis()->FindBin(posMaxValHist/15),histogrampointer->GetXaxis()->FindBin(posMaxValHist));   // look only for maxima with x greater than posMaxValHist/15, cut away noise
+        histogrampointer->GetXaxis()->SetRange(histogrampointer->GetXaxis()->FindBin(noiseborder),histogrampointer->GetXaxis()->FindBin(posMaxValHist));   // look only for maxima with x greater than noiseborder, cut away noise
 //         Int_t xValMax = histogrampointer->GetBinCenter(histogrampointer->GetMaximumBin());
-        TF1* fitFunc = new TF1("fitFunc",fitFuncType,posMaxValHist/15,posMaxValHist);
+        TF1* fitFunc = new TF1("fitFunc",fitFuncType,noiseborder,posMaxValHist);
         
         if (fitFuncType.Contains("gaus")) 
         {
             if (TString(histogrampointer->GetName()).Contains("Veto")) 
             { 
-//                 histogrampointer->GetXaxis()->SetRange(200,histogrampointer->GetXaxis()->FindBin(posMaxValHist));   // look only for maxima with x greater than posMaxValHist/15, cut away noise
-                
-                Float_t peak11 = histogrampointer->GetXaxis()->GetBinCenter(histogrampointer->GetMaximumBin());
                 Float_t peak1 = histogrampointer->GetMaximumBin();
-                Float_t peak1val = histogrampointer->GetBinContent(histogrampointer->GetMaximumBin());
-                Float_t avg = 0;
-                for(Int_t bin=histogrampointer->GetXaxis()->FindBin(posMaxValHist/15);bin<histogrampointer->FindLastBinAbove(0);bin++)
+                Float_t peak2;
+                if (histogrampointer->GetBinCenter(peak1)<200)
                 {
-                    avg += histogrampointer->GetBinContent(bin);
+                    Float_t avg = 0;
+                    for(Int_t bin=histogrampointer->GetXaxis()->FindBin(noiseborder);bin<histogrampointer->FindLastBinAbove(0);bin++)
+                        avg += histogrampointer->GetBinContent(bin);
+                    avg /= histogrampointer->FindLastBinAbove(0) - histogrampointer->GetXaxis()->FindBin(noiseborder);
+                    peak2 = histogrampointer->FindLastBinAbove(avg/4);
                 }
-                avg /= histogrampointer->FindLastBinAbove(0) - histogrampointer->GetXaxis()->FindBin(posMaxValHist/15);
-                Float_t peak2 = histogrampointer->FindLastBinAbove(avg/4);
-                Float_t peak21 = histogrampointer->GetXaxis()->GetBinCenter(peak2);
-                Float_t peak2val = histogrampointer->GetBinContent(peak2);
-                histogrampointer->GetXaxis()->SetRange(peak1,peak2);   // look only for maxima with x greater than posMaxValHist/15, cut away noise
-                //      
-                cout << coloryellow << "peak1 " << peak11 << endlr;
-                cout << coloryellow << "peak1val " << peak1val << endlr;
-                cout << coloryellow << "peak2 " << peak21 << endlr;
-                cout << coloryellow << "peak2val " << peak1val << endlr;
-                cout << coloryellow << "avg " << avg << endlr;
-                
+                else
+                {
+                    peak2 = peak1;
+                    peak1 = histogrampointer->GetXaxis()->FindBin(noiseborder);
+                }             
+                histogrampointer->GetXaxis()->SetRange(peak1,peak2-1);   // look only for maxima with x greater than noiseborder, cut away noise   
                 Float_t min = histogrampointer->GetMinimumBin();
+                histogrampointer->GetXaxis()->SetRange(min,histogrampointer->FindLastBinAbove(0));   // look only for maxima with x greater than noiseborder, cut away noise
+                // DEBUG 
+                if (verbose)
+                {
+                cout << coloryellow << "peak1 " << histogrampointer->GetXaxis()->GetBinCenter(peak1) << endlr;
+                cout << coloryellow << "peak1val " << histogrampointer->GetBinContent(peak1) << endlr;
+                cout << coloryellow << "peak2 " << histogrampointer->GetXaxis()->GetBinCenter(peak2) << endlr;
+                cout << coloryellow << "peak2val " << histogrampointer->GetBinContent(peak2) << endlr;
                 cout << coloryellow << "min " << histogrampointer->GetXaxis()->GetBinCenter( min ) << endlr;
-                histogrampointer->GetXaxis()->SetRange(min,histogrampointer->FindLastBinAbove(0));   // look only for maxima with x greater than posMaxValHist/15, cut away noise
-                
-                
-                
+                }
                 histogrampointer->Fit(fitFunc, "N,Q,W", "", histogrampointer->GetXaxis()->GetBinCenter(min), posMaxValHist);
+                histogrampointer->GetXaxis()->SetRange(histogrampointer->GetXaxis()->FindBin(noiseborder),histogrampointer->GetXaxis()->FindBin(posMaxValHist));   // look only for maxima with x greater than noiseborder, cut away noise
             } else {
-                histogrampointer->Fit(fitFunc, "N,Q,W", "", posMaxValHist/15, posMaxValHist);
+                histogrampointer->Fit(fitFunc, "N,Q,W", "", noiseborder, posMaxValHist);
             }
             Float_t sigma = fitFunc->GetParameter(2);
-            if (sigma > 50)
+            if (sigma > 260)
             {
-                cout << coloryellow << "sigma suspiciously height when fitting " << histogrampointer->GetName() << " spectrum: " << sigma  << endlr;
-                return 0;
+                if (verbose)
+                {
+                    cout << "Sigma suspiciously height when fitting " << histogrampointer->GetName() << " spectrum: " << sigma  << endl;
+                    cout << colorred << "  Could not find " << histogrampointer->GetName() << " peak" << endlr;
+                }
+                return posMax;
             }
+            else if (sigma > 40)
+            {
+                if (verbose)
+                {
+                    cout << "Sigma suspiciously height when fitting " << histogrampointer->GetName() << " spectrum: " << sigma  << endl;
+                    cout << coloryellow << "  Please check " << histogrampointer->GetName() << " peak position: " << colorreset << colorwhite << fitFunc->GetMaximumX()  << endlr;
+                }
+            }            
             posMax = fitFunc->GetMaximumX();
             fitFunc->DrawCopy("same");
             TString  legendEntry = TString(Form("%s",runcode.Data()));
@@ -833,13 +846,10 @@ Float_t Run::FitPerform(TH1F* histogrampointer, TString fitFuncType)
             Float_t fitMax3 = 1000;
             Float_t minFitMax = 1000;
             Float_t maxFitMax = 1000;
-                        
-            cout << "posMaxValHist/15: " << posMaxValHist/15 << endl;
-            
-            histogrampointer->Fit(fitFunc, "N,Q,W", "", posMaxValHist/15, posMaxValHist);
+            histogrampointer->Fit(fitFunc, "N,Q,W", "", noiseborder, posMaxValHist);
             fitMax1 = fitFunc->GetMaximumX();
             fitFunc->DrawCopy("same");         
-            histogrampointer->Fit(fitFunc, "N,Q,W", "", posMaxValHist/15, fitMax1*1.1);
+            histogrampointer->Fit(fitFunc, "N,Q,W", "", noiseborder, fitMax1*1.1);
             fitMax2 = fitFunc->GetMaximumX();
             fitFunc->SetLineColor(kBlue);
             fitFunc->SetLineStyle(2); // dashed
@@ -895,7 +905,7 @@ TCanvas* Run::plot1DHistogram(TH1F* onehistogram, TString fitFuncType, TString t
     TCanvas* canvas = new TCanvas(canvasname, canvastitle, 900, 700);
     onehistogram->SetTitle(titlestr);
     onehistogram->Draw();
-    Float_t maxValue = FitPerform(onehistogram, fitFuncType);
+    Float_t maxValue = FitPerform(onehistogram, fitFuncType, false);
     plotVerticalLine(onehistogram, maxValue);
     TLegend* leg = new TLegend(0.8,0.8,0.89,0.89);//(0.6,0.7,0.89,0.89);
     leg->SetFillColor(0);