Analyzer: Bugfix, calibration was broken (by me)

diff --git a/MABS_run_analyzer/ChargeSpektrum.c b/MABS_run_analyzer/ChargeSpektrum.c
index a30464cd2ff39d8bef1193bd06ca8181c9e42779..b24821b6cb608557c9702643d371f582b4e81f02 100644 (file)
--- a/MABS_run_analyzer/ChargeSpektrum.c
+++ b/MABS_run_analyzer/ChargeSpektrum.c
@@ -44,6 +44,9 @@ void ChargeSpektrum(TString runnumber = "")
     // variables for classical analysis
     vector<HistogramType*> compareHistogramClassVectorClassic;
     
+    // variables for cluster threshold analysis
+    vector<HistogramType*> compareHistogramClassVectorClusterThreshold;
+    
     // variables for calibrated analysis
     vector<HistogramType*> compareHistogramClassVectorCalibrated;
     
@@ -80,6 +83,9 @@ void ChargeSpektrum(TString runnumber = "")
     vector<TH1FO*> compareHistogramVectorClusterMultiplicity;
     vector<TNtupleO*> ntupleVectorAvgChargeColl;
     
+    // f0 analysis variables  
+    vector<TH1FO*> compareHistogramAllRunsVectorF0;
+    
     for(Int_t runi=0; runi<numberRuns; runi++) // loop over runs read from file or given as aparameters
     {
         if (runList[runi]>0)
@@ -143,11 +149,17 @@ void ChargeSpektrum(TString runnumber = "")
                     
                     // use only lower case when matching anaylsis types
                     
-                    if (analysisType.Contains("clas")) { // classic analysis
+                    if (analysisType.Contains("clas")) { // classic analysis                        
                         compareHistogramClassVectorClassic.push_back(runs[runi]->histogram->normalized);
                         if (runi+1 == numberRuns) 
                             compareHistogramClassVectorVector.push_back(compareHistogramClassVectorClassic);
-                    }           
+                    }   
+                    
+                    if (analysisType.Contains("clusthresh")) { // cluster charge > 2 * noise in cluster
+                        compareHistogramClassVectorClusterThreshold.push_back(runs[runi]->histogramthreshold->normalized);
+                        if (runi+1 == numberRuns) 
+                            compareHistogramClassVectorVector.push_back(compareHistogramClassVectorClusterThreshold);
+                    }    
                     
                     if (analysisType.Contains("dynthresh") || analysisType.Contains("dynamicalthresh")) { // fixed threshold analysis
                         vector<HistogramType*> compareHistogramClassVectorDynamicalThreshold;
@@ -158,7 +170,7 @@ void ChargeSpektrum(TString runnumber = "")
                     }   
                     
                     if (analysisType.Contains("cali")) { // calibrated analysis
-                        compareHistogramClassVectorCalibrated.push_back(runs[runi]->histogram->calibrated->normalized);
+                        compareHistogramClassVectorCalibrated.push_back(runs[runi]->histogram->normalized->calibrated);
                         if (runi+1 == numberRuns)
                             compareHistogramClassVectorVector.push_back(compareHistogramClassVectorCalibrated);
                     }  
@@ -391,9 +403,12 @@ void ChargeSpektrum(TString runnumber = "")
                         // Plot F0 signal
                         vector<TH1FO*> compareHistogramVectorF0;
                         compareHistogramVectorF0.push_back(runs[runi]->averageF0DistrEnd);
-                        compareHistogramVectorVector.push_back(compareHistogramVectorF0);    
+                        compareHistogramVectorF0.push_back(runs[runi]->averageF0DistrStart);
+                        CompareHistograms(&compareHistogramVectorF0);                        
+                        compareHistogramVectorF0.clear();
                         
-                        runs[runi]->plot1DHistogram( runs[runi]->averageF0DistrEnd, "gaus", true);
+                        compareHistogramAllRunsVectorF0.push_back(runs[runi]->averageF0Distr);
+                        runs[runi]->plot1DHistogram( runs[runi]->averageF0Distr, "gaus", true);
                     }
 
                     // write all histograms to a dat file, can be imported by ORIGIN or Excel. Leave this line ;)

diff --git a/MABS_run_analyzer/HistogramType.c b/MABS_run_analyzer/HistogramType.c
index 751785269ad1e1b7f0886b7d1b1cf56b9f168177..38df91541fade7141730fecdb345104a92c1177c 100644 (file)
--- a/MABS_run_analyzer/HistogramType.c
+++ b/MABS_run_analyzer/HistogramType.c
@@ -233,11 +233,11 @@ void HistogramType::calibrateHistogram(TH1FO* &histogrampointernew, TH1FO* &hist
     histogrampointernew->GetXaxis()->SetTitle("Collected charge [e]");    
     int nbins = histogrampointernew->GetXaxis()->GetNbins();
     double new_bins[nbins+1];
-    if (scalefactor != 0)
+    if (scalefactor == 0)
         scalefactor = gain;
     for(int i=0; i <= nbins; i++){
         new_bins[i] = histogrampointernew->GetBinLowEdge(i)*scalefactor;
-//         cout << "histogrampointernew->GetBinCenter(i): "  << histogrampointernew->GetBinLowEdge(i) << " * gain " << gain << " = " << histogrampointernew->GetBinLowEdge(i)*gain << endl;
+//          cout << "histogrampointernew->GetBinCenter(i): "  << histogrampointernew->GetBinLowEdge(i) << " * gain " << gain << " = " << histogrampointernew->GetBinLowEdge(i)*gain << endl;
     }    
     histogrampointernew->SetBins(nbins, new_bins);    
     histogrampointernew->GetYaxis()->SetTitle(Form("Entries [1/%.1f e]",histogrampointernew->GetBinWidth(1)));
@@ -251,7 +251,7 @@ void HistogramType::calibrateHistogramYAxis(TH1FO* &histogrampointernew, TH1FO*
     histogrampointernew->SetTitle(Form("%s Calibr.", histogrampointerold->GetTitle()));
     histogrampointernew->GetYaxis()->SetTitle("Collected charge [e]");    
     int nbins = histogrampointernew->GetXaxis()->GetNbins();
-    if (scalefactor != 0)
+    if (scalefactor == 0)
         scalefactor = gain;
     for(int x=0; x <= nbins; x++){
         histogrampointernew->SetBinContent(x,histogrampointerold->GetBinContent(x)*scalefactor);
@@ -769,7 +769,11 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType
         const Double_t def_bgslope=0;
         const Double_t def_bgoffs=histogrampointer->GetBinContent(histogrampointer->FindBin((noiseborder+def_peakcenter)/2));
 //          cout << colorcyan << "histogrampointer->FindBin((noiseborder+def_peakcenter)/2): " << histogrampointer->FindBin((noiseborder+def_peakcenter)/2) << endlr;
-                         
+        
+        if (verbose) {
+            cout << colorcyan << "def_amplitude: " << def_amplitude << endlr;
+            cout << colorcyan << "def_peakcenter: " << def_peakcenter << endlr;
+        }
         fitFunc->SetLineWidth(4);
         fitFunc->SetLineColor(kGreen);
         // set start values for some parameters
@@ -777,14 +781,14 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType
         fitFunc->SetParameter(0,def_amplitude);
         fitFunc->SetParName(1,"peak centroid");
         fitFunc->SetParameter(1,def_peakcenter);
-        fitFunc->SetParLimits(1,def_peakcenter*0.8,def_peakcenter*1.2);
+//         fitFunc->SetParLimits(1,def_peakcenter*0.8,def_peakcenter*1.2);
         fitFunc->SetParName(2,"Gaussian sigma");
         fitFunc->SetParameter(2,def_gausssig);
         fitFunc->SetParLimits(2,2,50);
         fitFunc->SetParLimits(2,0,10000);
         fitFunc->SetParName(3,"Distance from Gauss");
         fitFunc->SetParameter(3,def_distgauss);
-        fitFunc->SetParLimits(3,0,def_distgauss*2);
+        fitFunc->SetParLimits(3,0,def_distgauss*4);
         fitFunc->SetParName(4,"background slope");
         fitFunc->SetParameter(4,def_bgslope);
 //        fitFunc->SetParLimits(4,def_bgslope-0.1,def_bgslope+0.1);
@@ -793,8 +797,8 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType
         fitFunc->SetParLimits(5,0,def_bgoffs*4);
         
         // TODO: This fix disables the background
-        fitFunc->FixParameter(4,0);
-        fitFunc->FixParameter(5,0);
+//         fitFunc->FixParameter(4,0);
+//         fitFunc->FixParameter(5,0);
         
         int fittries = 0;
         Bool_t failed = false;
@@ -823,7 +827,32 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType
             else
                 fittries = 100;
 //              cout << colorcyan << gMinuit->fCstatu.Data() << " fit tries: " << fittries << endlr;            
-        } while (fittries < 6);
+        } while (fittries < 6);        
+        fittries = 0;
+        do {
+            failed = false;
+            fit_result_ptr = histogrampointer->Fit(fitFunc, "N,M,W,Q,S", "", posMaxValHist-posMaxValHist/6*(fittries+1), posMaxValHist);
+            //             cout << colorcyan << " AFTER fit " << endlr;   
+            if (gMinuit == nullptr) { 
+                failed = true;
+            } else 
+            {
+                if (gMinuit->fCstatu == "NOT POSDEF" || gMinuit->fCstatu.Contains("FAILED")) {
+                    failed = true;
+                }
+            }
+            if (failed)
+            {
+                fitFunc->SetParameter(0,def_amplitude);
+                fitFunc->SetParameter(1,def_peakcenter);
+                fitFunc->SetParameter(2,def_gausssig);
+                fitFunc->SetParameter(3,def_distgauss);
+                fitFunc->SetParameter(4,def_bgslope);
+                fitFunc->SetParameter(5,def_bgoffs);
+            }
+            else
+                fittries = 100;        
+        } while (fittries < 6); 
         if (failed)
         {
             failed = false;
@@ -870,7 +899,7 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType
          
         if (failed)
         {
-            cout << colorred << " Could not find " << histogrampointer->GetName() << " peak" << endlr;
+            cout << colorred << " Could not find " << histogrampointer->GetName() << " calibration peak" << endlr;
             parameters = (Double_t *)calloc(10, sizeof(Double_t));
             return parameters;
         }
@@ -917,13 +946,15 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType
         }
 
 //             DEBUG
-//         TCanvas* canvas = new TCanvas(histogrampointer->GetName(), histogrampointer->GetName(), 900, 700);
-//         histogrampointer->Draw();
 //         fitFunc->Draw("SAME");
 
         fitFunc->SetLineStyle(1); // normal for the following fits
         if (verbose)
+        {
+            TCanvas* canvas = new TCanvas(histogrampointer->GetName(), histogrampointer->GetName(), 900, 700);
+            histogrampointer->Draw();
             fitFunc->Draw("SAME");
+        }
         
     }
 

diff --git a/MABS_run_analyzer/HistogramType.h b/MABS_run_analyzer/HistogramType.h
index 97f0e4c79c6cb46530c0d3b235df87b0d8b0e7c6..4434506f0dc20135f573f5f5376b0900f09f9431 100644 (file)
--- a/MABS_run_analyzer/HistogramType.h
+++ b/MABS_run_analyzer/HistogramType.h
@@ -18,6 +18,11 @@
 
 class HistogramType;
 
+
+/**
+ * @brief A class which stores one Histogram, derived from TH1F, see Cern root documentation for more information
+ * 
+ */
 class TH1FO: public TH1F {
 public:
     HistogramType* itsHistogramType = 0;

diff --git a/MABS_run_analyzer/Run.c b/MABS_run_analyzer/Run.c
index 074af13102fb0d4ba54d27f2363556377c0e7258..93faba2f4ef019f25b42ae3cc4b9312092c1f04d 100644 (file)
--- a/MABS_run_analyzer/Run.c
+++ b/MABS_run_analyzer/Run.c
@@ -723,9 +723,9 @@ Bool_t Run::setNoisethresholdborderE(Float_t noisethresholdborder)
 Bool_t Run::rescaleHistogramClasses()
 {
     float_t vetopeakposition = -1;
-    if ( histogram->posVeto > 0 )
+    if ( histogramdynamicalthreshold->posVeto > 0 )
     {
-        vetopeakposition = histogram->posVeto;
+        vetopeakposition = histogramdynamicalthreshold->posVeto;
         cout << " Use calibration obtained from this run, position veto: " << vetopeakposition << endlr;
     }
     else if ( labbook.posVetoDB > 0 )
@@ -1044,7 +1044,7 @@ void Run::updateDatabase() {
     constructUpdateString(&sqlupdatequery, "SigmaF0",   labbook.sigmaF0, 6, 0, 10000);
     constructUpdateString(&sqlupdatequery, "SeedIntegralMinusBG2Sigma",   histogramclassToUseForDB->normalized->integralSeed, 10, 0, 1e7);
     constructUpdateString(&sqlupdatequery, "SeedIntegralMinusBG2SigmaErr",   histogramclassToUseForDB->normalized->integralSeedErr, 5, 0, 1e7);
-    constructUpdateString(&sqlupdatequery, "VetoPeak",   histogramclassToUseForDB->normalized->posVeto, 4, 0, 1000);
+    constructUpdateString(&sqlupdatequery, "VetoPeak",   histogramdynamicalthreshold->normalized->posVeto, 4, 0, 1000);
     constructUpdateString(&sqlupdatequery, "VetoIntegral",   histogramclassToUseForDB->normalized->integralVeto, 10, 0, 1e7);
     constructUpdateString(&sqlupdatequery, "SumIntegral",   histogramclassToUseForDB->normalized->integralSum, 10, 0, 1e7);
     constructUpdateString(&sqlupdatequery, "SumIntegralErr",   histogramclassToUseForDB->normalized->integralSumErr, 5, 0, 1e7);
@@ -1244,6 +1244,10 @@ Bool_t Run::binNoise(HistogramType* oneHistogramClass)
 
 Bool_t Run::binF0()
 {
+    // for better performace analyze 1000 frames at the befginng of a run (skip the very first 1000 frames, they are "Dirty" and last 1000 frames
+    // afterwards sum them up for an 'average' F0 distribution over the run
+    
+    
     averageF0Hist = (TH1FO*) new TH1F(Form("%d AvgF0 for whole chip",labbook.runnumber), Form("%d AvgF0 for whole chip",labbook.runnumber), processed->fHitTree->GetEntries(), 0, processed->fHitTree->GetEntries());  
     averageF0Hist->SetLineStyle(rootlinestyle[plotStyle]); 
     averageF0Hist->SetLineColor(rootcolors[plotStyle]);
@@ -1299,12 +1303,17 @@ Bool_t Run::binF0()
     averageF0DistrStart->itsHistogramType = histogram;  
     averageF0PixelwiseEnd = (TH1FO*) averageF0PixelwiseStart->Clone();
     averageF0PixelwiseEnd->SetNameTitle(Form("%d AvgF0 per Pixel, last 1000 frames",labbook.runnumber),Form("%d AvgF0 per Pixel, last 1000 frames",labbook.runnumber));
-    averageF0PixelwiseEnd->SetLineColor(12);
+    averageF0PixelwiseEnd->SetLineColor(1+100);
     averageF0PixelwiseEnd->itsHistogramType = histogram;  
     averageF0DistrEnd = (TH1FO*) averageF0DistrStart->Clone();
     averageF0DistrEnd->SetNameTitle(Form("%d AvgF0, last 1000 frames",labbook.runnumber),Form("%d Average F0, last 1000 frames",labbook.runnumber));
-    averageF0DistrEnd->SetLineColor(12);
-    averageF0DistrEnd->itsHistogramType = histogram;  
+    averageF0DistrEnd->SetLineColor(1+100);
+    averageF0DistrEnd->itsHistogramType = histogram;
+    averageF0Distr = (TH1FO*) averageF0DistrStart->Clone();
+    averageF0Distr->SetNameTitle(Form("%d AvgF0",labbook.runnumber),Form("%d Average F0",labbook.runnumber));
+    averageF0Distr->SetLineColor(2);
+    averageF0Distr->itsHistogramType = histogram;  
+    
 
     Double_t const probabilities[] = {0.158655254, 0.5, 1-0.158655254}; // sigma/2 from gaus to the left and to the right //{0.17, 0.5, 1-0.17};
     Double_t avgF0ForPixel [cursensorinfo.columns*cursensorinfo.rows];
@@ -1356,6 +1365,14 @@ Bool_t Run::binF0()
     }
     
     
+    //  sum up end and beginning
+    averageF0Distr->Add(averageF0DistrStart);
+    averageF0Distr->Add(averageF0DistrEnd);
+    for (Int_t bini=1; bini <= averageF0Distr->GetXaxis()->GetNbins(); bini++) { // loop over all pixel
+        averageF0Distr->SetBinContent(bini, averageF0Distr->GetBinContent(bini)/2);
+    }
+    
+    
 //     TCanvas *c1 = new TCanvas("c1","c1",600,400);
 //     averageF0DistrStart->Draw("");
 //     averageF0DistrEnd->Draw("SAME");

diff --git a/MABS_run_analyzer/Run.h b/MABS_run_analyzer/Run.h
index f0ac8cf23a573e3f7c5ba9cb25b0e9b58b2d42a6..ad48f653be9be465d589964072c5acd04de1070b 100644 (file)
--- a/MABS_run_analyzer/Run.h
+++ b/MABS_run_analyzer/Run.h
@@ -413,6 +413,7 @@ public:
         
     TH1FO* averageF0DistrStart = 0; // Over first 1000 frames
     TH1FO* averageF0DistrEnd = 0;   // Over last  1000 frames 
+    TH1FO* averageF0Distr = 0;   // Over all frames 
             
     /** @brief Plots all histograms from #Run::HistogramType into one canvas 
      *