From: Benjamin Linnik <blinnik@jspc61.x-matter.uni-frankfurt.de>
Date: Wed, 16 Aug 2017 13:21:39 +0000 (+0200)
Subject: Analyzer: New gauss fitting routines, ajustments for depletion volume studies
X-Git-Url: https://jspc29.x-matter.uni-frankfurt.de/git/?a=commitdiff_plain;h=00994fb0bf36f5557381fcc6253416c9c465ce9c;p=radhard.git

Analyzer: New gauss fitting routines, ajustments for depletion volume studies
---

diff --git a/MABS_run_analyzer/ChargeSpektrum.c b/MABS_run_analyzer/ChargeSpektrum.c
index 25d58f6..fbb844f 100644
--- a/MABS_run_analyzer/ChargeSpektrum.c
+++ b/MABS_run_analyzer/ChargeSpektrum.c
@@ -118,11 +118,30 @@ void ChargeSpektrum(TString runnumber = "")
 //              runs[runi]->plot1DHistogram(runs[runi]->histogram, runs[runi]->histogramwoRTS->pixeltimefiredsorted, "", 0);
 
 // Integralberechnung vom Veto Peak
-                    runs[runi]->plot1DHistogram(runs[runi]->histogram, runs[runi]->histogram->Seed, "GaussTail", true);
-                    runs[runi]->plot1DHistogram(runs[runi]->histogram->normalized, runs[runi]->histogram->normalized->Seed, "GaussTail", true);
+//                     runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS, runs[runi]->histogramwoRTS->Seed, "GaussTail", true);
+//                     runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS, runs[runi]->histogramwoRTS->Sum, "GaussTail", true);
 					
-					compareHistogramVector.push_back(runs[runi]->histogram->normalized->Seed);
-					compareHistogramVector.push_back(runs[runi]->histogram->Seed);
+                    // 					compareHistogramVector.push_back(runs[runi]->histogram->normalized->Seed);
+//                     compareHistogramVector.push_back(runs[runi]->histogram->normalized->Seed);
+//                     compareHistogramVector.push_back(runs[runi]->histogramthreshold->normalized->Seed);
+//                     compareHistogramVector.push_back(runs[runi]->histogramwoRTS->normalized->Seed);
+//                     compareHistogramVector.push_back(runs[runi]->histogramwoRTS->normalized->Sum);
+//                     compareHistogramVector.push_back(runs[runi]->histogramwoRTSAggresive->normalized->Seed);
+                    
+                    runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS, runs[runi]->histogramwoRTS->Seed, "GaussTail", true, true, false, 500);
+                    runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS, runs[runi]->histogramwoRTS->Sum, "GaussTail", true, true, false, 500);
+                    
+//                     compareHistogramVector2.push_back(runs[runi]->histogramwoRTS->normalized->Seed);
+//                     compareHistogramVector3.push_back(runs[runi]->histogramwoRTS->normalized->Sum);
+                    
+//                     compareHistogramVector3.push_back(runs[runi]->histogram->pixeltimefired);
+//                     compareHistogramVector3.push_back(runs[runi]->histogramwoRTS->pixeltimefired);
+//                     compareHistogramVector3.push_back(runs[runi]->histogramwoRTSAggresive->pixeltimefired);
+                    //                     
+                    //                     compareHistogramVector5.push_back(runs[runi]->histogram->pixeltimefiredsorted);
+                    
+//                     compareHistogramVector2.push_back(runs[runi]->histogramwoRTS->Seed);
+//                     compareHistogramVector3.push_back(runs[runi]->histogramwoRTS->Sum);
 
 //                     runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS->normalized->calibrated, runs[runi]->histogramwoRTS->normalized->calibrated->Sum, "gaus", true);
                     
@@ -178,8 +197,6 @@ void ChargeSpektrum(TString runnumber = "")
 //                     
 // //                         compareHistogramClassVector.push_back(runs[runi]->histogram);
 //                     //                         compareHistogramVector2.push_back(runs[runi]->histogram->LeakageCurrentInPixel);
-//                     compareHistogramVector3.push_back(runs[runi]->histogram->LeakageCurrentInPixelSorted);
-//                     compareHistogramVector4.push_back(runs[runi]->histogramwoRTS->LeakageCurrentInPixelSorted);
 //                     
 //                     compareHistogramVector5.push_back(runs[runi]->histogram->pixeltimefiredsorted);
 //                     compareHistogramVector6.push_back(runs[runi]->histogramwoRTS->pixeltimefiredsorted);
diff --git a/MABS_run_analyzer/ChargeSpektrumFunctions.c b/MABS_run_analyzer/ChargeSpektrumFunctions.c
index 4c0bd2c..909ad0e 100644
--- a/MABS_run_analyzer/ChargeSpektrumFunctions.c
+++ b/MABS_run_analyzer/ChargeSpektrumFunctions.c
@@ -400,7 +400,7 @@ Bool_t CompareHistograms(vector<TH1F*>* ptCompareHistogramVector, TString titles
             
             curhistogramclone->SetAxisRange(0,lastbin1*1.1,"X");
             //             curhistogramclone->GetYaxis()->SetRangeUser(5,heighestval1*4);
-            //             gPad->SetLogy(1);
+                         gPad->SetLogy(1);
             owntitle->Clear();
             owntitle->AddText(trimRunnumberAtBegin(curhistogramclone->GetName()));
             owntitle->Draw("SAME");
diff --git a/MABS_run_analyzer/HistogramType.c b/MABS_run_analyzer/HistogramType.c
index 662affa..93f7ae5 100644
--- a/MABS_run_analyzer/HistogramType.c
+++ b/MABS_run_analyzer/HistogramType.c
@@ -126,8 +126,10 @@ Bool_t HistogramType::calibrateHistograms( Float_t gotgain ) {
     calibrated->posSum = posSum * gain;
     calibrated->posVeto = posVeto * gain;
     calibrated->integralSeed = integralSeed;
+    calibrated->integralSeedErr = integralSeedErr;
     calibrated->integralVeto = integralVeto;
     calibrated->integralSum = integralVeto;
+    calibrated->integralSumErr = integralSumErr;
     calibrated->posSeedPerc = posSeedPerc;
     calibrated->sigmaSeedPerc = sigmaSeedPerc;
     calibrated->avgNoise = avgNoise * gain;
@@ -143,9 +145,9 @@ Bool_t HistogramType::calibrateHistograms( Float_t gotgain ) {
     calibrated->RTSpixel = RTSpixel;
     calibrated->percentageofRTSpixel = percentageofRTSpixel;
     calibrated->avgLeakageCurrentInChip = avgLeakageCurrentInChip * gain;
-    calibrated->medianLeakageCurrent = medianLeakageCurrent * gain / 3.7 / pow10(-3) * (1.6*pow10(-18)) * pow10(15);
-    calibrated->medianLeakageCurrentMinus = medianLeakageCurrentMinus * gain / 3.7 / pow10(-3) * (1.6*pow10(-18)) * pow10(15);
-    calibrated->medianLeakageCurrentPlus = medianLeakageCurrentPlus * gain / 3.7 / pow10(-3) * (1.6*pow10(-18)) * pow10(15);
+    calibrated->medianLeakageCurrent = medianLeakageCurrent * gain / 3.7 / pow10(-3) * (1.6*pow10(-19)) * pow10(15);
+    calibrated->medianLeakageCurrentMinus = medianLeakageCurrentMinus * gain / 3.7 / pow10(-3) * (1.6*pow10(-19)) * pow10(15);
+    calibrated->medianLeakageCurrentPlus = medianLeakageCurrentPlus * gain / 3.7 / pow10(-3) * (1.6*pow10(-19)) * pow10(15);
     
     calibrated->iscalibrated = true;
     
@@ -225,9 +227,11 @@ Bool_t HistogramType::normalizeHistograms( Int_t got_frames_found ) {
     normalized->posSum = posSum;
     normalized->posVeto = posVeto;
     normalized->integralSeed = integralSeed/frames_found*pow10(6);
+    normalized->integralSeedErr = integralSeedErr/frames_found*pow10(6);
     normalized->integralVeto= integralVeto/frames_found*pow10(6);
     normalized->integralVeto= integralVeto/frames_found*pow10(6);
     normalized->integralSum = integralSum/frames_found*pow10(6);
+    normalized->integralSumErr = integralSumErr/frames_found*pow10(6);
     normalized->posSeedPerc = posSeedPerc;
     normalized->sigmaSeedPerc = sigmaSeedPerc;
     normalized->noisethresholdborder = noisethresholdborder;
@@ -278,10 +282,13 @@ void HistogramType::normalizeHistogramXAxis(TH1F* &histogrampointernew, TH1F* &h
 }
 
 
-Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType, Bool_t verbose, Float_t fitstart) {
+//Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType, TFitResult* fit_result, Bool_t verbose, Float_t fitstart) {
+Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType, TFitResult* fit_result_giveback, Bool_t verbose, Float_t fitstart) {
+    TFitResultPtr fit_result_ptr;
     Double_t* parameters = 0;
     Float_t posMax = 0;
     Float_t integralPeak = 0;
+    Double_t integralPeakError = 0;
     Float_t posMax2 = 0;
     Float_t posMaxValHist = histogrampointer->GetXaxis()->GetXmax();
     // empirical value for USB system is 50, should be overwritten in 99% of cases
@@ -297,6 +304,8 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType,
     
     if (fitFuncType.EqualTo("gaus")) 
     {
+        parameters = (Double_t *)calloc(10, sizeof(Double_t));
+        
         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,noiseborder,posMaxValHist);
@@ -347,11 +356,11 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType,
             int fittries = 0;
             do {
 //                cout << fittries <<  ": New range: " <<  histogrampointer->GetXaxis()->GetBinCenter(min+fittries*min/20) << " to " << posMaxValHist << endl;
-                histogrampointer->Fit(fitFunc, "N,Q,W", "", histogrampointer->GetXaxis()->GetBinCenter(min+fittries*min/20), posMaxValHist);    
+                fit_result_ptr = histogrampointer->Fit(fitFunc, "NQWS", "", histogrampointer->GetXaxis()->GetBinCenter(min+fittries*min/20), posMaxValHist);    
 //                 cout << colorred << gMinuit->fCstatu.Data() << endlr;            
             } while (gMinuit->fCstatu.Contains("FAILED") && fittries++ < 8);
             posMax = histogrampointer->GetXaxis()->GetBinCenter(histogrampointer->GetMaximumBin()); // Methode 1
-            integralPeak = histogrampointer->Integral(min+fittries*min/20, histogrampointer->GetXaxis()->FindBin(posMaxValHist), "width");
+            integralPeak = histogrampointer->IntegralAndError(min+fittries*min/20, histogrampointer->GetXaxis()->FindBin(posMaxValHist), integralPeakError, "width");
 //             if (verbose)
 //             {
 //                 cout << coloryellow << "min " << histogrampointer->GetXaxis()->GetBinCenter( min ) << endlr;
@@ -364,21 +373,97 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType,
             if (verbose)
               fitFunc->DrawCopy("same");
         } else {
-            histogrampointer->Fit(fitFunc, "N,Q,W", "", noiseborder, posMaxValHist);
-            noiseborder = FindBorderToPeak(histogrampointer, noiseborder,fitFunc->GetParameter(1), verbose); // starting point of histogram integration
-            integralPeak = histogrampointer->Integral(histogrampointer->GetXaxis()->FindBin(noiseborder), histogrampointer->GetXaxis()->FindBin(posMaxValHist), "width");
+            // modified the gaussian approuch to be more powerful, escpecially the
+            // integral is calculated in a +- 2 sigma region.
+            
+            const Double_t def_amplitude=15;
+            const Double_t def_peakcenter=histogrampointer->GetBinCenter(histogrampointer->GetMaximumBin());
+            const Double_t def_gausssig=21;
+            // set start values for some parameters
+            fitFunc->SetParName(0,"amplitude of peak");
+            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->SetParName(2,"Gaussian sigma");
+            fitFunc->SetParameter(2,def_gausssig);
+            fitFunc->SetParLimits(2,0,100);
+                            
+            fitFunc->SetLineWidth(4);
+            fitFunc->SetLineColor(kGreen);
+            
+            int fittries = 0;
+            Bool_t failed = false;
+            do {
+                if (failed)
+                {
+                    fitFunc->SetParameter(0,def_amplitude*(1.0-0.1*++fittries));
+                    fitFunc->SetParameter(1,def_peakcenter);
+                    fitFunc->SetParameter(2,def_gausssig);
+                } else fittries = 100;         
+                failed = false;
+                fit_result_ptr = histogrampointer->Fit(fitFunc, "NMWQS", "", noiseborder, posMaxValHist);  
+                if (gMinuit == nullptr) { 
+                    failed = true;
+                } else 
+                {
+                    if (gMinuit->fCstatu == "NOT POSDEF" || gMinuit->fCstatu.Contains("FAILED")) {
+                        failed = true;
+                    }
+                }                
+            } while (fittries < 10);   
+            // get parameters
+            for (Int_t pari=0; pari<3; pari++) {
+                //cout << colorcyan << fitFunc->GetParameter(pari) << endlr;
+                parameters[pari]=fitFunc->GetParameter(pari);
+            }
+            if (failed || ( parameters[2]/parameters[0] > 100 ) )
+            {
+                fitFunc->SetParameter(0,def_amplitude);
+                fitFunc->SetParameter(1,def_peakcenter);
+                fitFunc->SetParameter(2,def_gausssig);
+                fit_result_ptr = histogrampointer->Fit(fitFunc, "NMWQS", "", def_peakcenter*0.7, def_peakcenter*1.3);  
+            }
+            if (failed)
+            {
+                cout << colorred << " Could not find " << histogrampointer->GetName() << " peak" << endlr;
+                parameters = (Double_t *)calloc(8, sizeof(Double_t));
+                return parameters;
+            }
+            
+            // get parameters
+            for (Int_t pari=0; pari<3; pari++) {
+                //cout << colorcyan << fitFunc->GetParameter(pari) << endlr;
+                parameters[pari]=fitFunc->GetParameter(pari);
+            }
+            // https://suchideas.com/articles/maths/applied/histogram-errors/#The_Sensible_Way
+            // make an error estimate, in case of rare events one can use the poisson distribution
+            // error bars become +- sqrt(n) for each bin, where n is the bin content.
+                       
+            parameters[7] = parameters[1] - 2*parameters[2] ; // starting point of histogram integration
+            parameters[8] = parameters[1] + 2*parameters[2] ; // end point of histogram integration
+            integralPeak  = histogrampointer->IntegralAndError(histogrampointer->GetXaxis()->FindBin( parameters[7]  ), histogrampointer->GetXaxis()->FindBin( parameters[8]), integralPeakError, "width"); // integral value of histogram (NOT fit), normaliezed with bin size!
             posMax = fitFunc->GetMaximumX(); // Methode 2
             fitFunc->SetLineStyle(1); // normal for the following fits
             if (verbose)
                 fitFunc->Draw("SAME");
         }
-        parameters = (Double_t *)calloc(5, sizeof(Double_t));
+
+        
+        // parameters are 
+        // parameters[0]: amplitude of peak
+        // parameters[1]: x position of peak
+        // parameters[2]: sigma of gaussian
+        // parameters[3]: integral of peak
+        // parameters[6]: integral of peak
+        // parameters[7]: start of integration
+        // parameters[8]: end of integratio
         for (Int_t pari=0; pari<3; pari++) {
             //cout << colorcyan << fitFunc->GetParameter(pari) << endlr;
             parameters[pari]=fitFunc->GetParameter(pari);
         }
-        parameters[3]=integralPeak;
-        parameters[4] = noiseborder;
+        parameters[6]=integralPeak;
+        parameters[9]=integralPeakError;
         Float_t sigma = fitFunc->GetParameter(2);
         posMax2 = fitFunc->GetMaximumX(); // Methode 2
 //         if (verbose)
@@ -417,7 +502,20 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType,
                     cout << coloryellow << "  Please check " << histogrampointer->GetName() << " peak position: " << colorreset << colorwhite << posMax  << endlr;
                 }
 //             }
-        }       
+        }  
+        if (verbose) {
+			cout << colorcyan << "Depletion voltage: " << labbook->depletionV << endlr; 			
+            cout << colorcyan << "Integral from bin : " << histogrampointer->FindBin(parameters[7])  << " to " << histogrampointer->GetXaxis()->FindBin(parameters[8]) << endlr;  
+            cout << colorcyan << "Integral from val : " << parameters[7] << " to " << parameters[8] << endlr;      
+            cout << colorcyan << "Integral value without bg: " << parameters[6] << endlr;                
+            cout << colorcyan << "Integral error: " << parameters[9] << endlr;           
+			cout << colorcyan << "Gauss Sigma: " << parameters[2] << endlr; 			
+            cout << colorcyan << "Number of events : " << frames_found << endlr; 
+            
+            (fit_result_ptr.Get())->Print("V"); 
+            
+//             cout <<  fitFunc->GetParError(0) << endlr;             // retrieve the error for the parameter 0            
+        }
         
 //         TString  legendEntry = TString(Form("%s","Gaus fit"));
 //         TLegend* leg = new TLegend(0.5,0.5,0.89,0.89);//(0.6,0.7,0.89,0.89);
@@ -441,7 +539,7 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType,
         Float_t fitMax3 = 1000;
         Float_t minFitMax = 1000;
         Float_t maxFitMax = 1000;
-        histogrampointer->Fit(fitFunc, "N,Q,W", "", noiseborder, posMaxValHist);
+        fit_result_ptr = histogrampointer->Fit(fitFunc, "N,Q,W,S", "", noiseborder, posMaxValHist);
         for (Int_t pari=0; pari<3; pari++)
             parameters[pari]=fitFunc->GetParameter(pari);
 //         posMax = fitFunc->GetParameter(1); // new method, just use MPV of fit, as noise border is well known. Compare with old method anyway and warn if something suspecious.
@@ -454,12 +552,12 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType,
 //         cout << colorred << "noiseborder: " << noiseborder << endlr;
         fitMax1 = fitFunc->GetMaximumX();
 //                       fitFunc->DrawCopy("same");         
-        histogrampointer->Fit(fitFunc, "N,Q,W", "", noiseborder, fitMax1*1.1);
+        fit_result_ptr = histogrampointer->Fit(fitFunc, "N,Q,W,S", "", noiseborder, fitMax1*1.1);
         fitMax2 = fitFunc->GetMaximumX();
         fitFunc->SetLineColor(kBlue);
         fitFunc->SetLineStyle(2); // dashed
 //                       fitFunc->DrawCopy("same");
-        histogrampointer->Fit(fitFunc, "N,Q,W", "", fitMax1*0.9, posMaxValHist);
+        fit_result_ptr = histogrampointer->Fit(fitFunc, "N,Q,W,S", "", fitMax1*0.9, posMaxValHist);
         //             histogrampointer->Fit(fitFunc, "N,Q,W", "", fitMax1, histogrampointer->GetBinCenter(bini));
         fitMax3 = fitFunc->GetMaximumX();
         fitFunc->SetLineColor(kGreen);
@@ -507,7 +605,7 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType,
         //             fitFcn->FixParameter(3,-22.43016284);
         
         //histogrampointer->Fit("fitFcn","V+","ep");
-        histogrampointer->Fit("fitFcn","Q,W","ep");        
+        fit_result_ptr = histogrampointer->Fit("fitFcn","Q,W,S","ep");        
         parameters = (Double_t *)calloc(4, sizeof(Double_t));
         for (Int_t pari=0; pari<4; pari++)
             parameters[pari]=fitFcn->GetParameter(pari);
@@ -532,7 +630,7 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType,
         
 //         histogrampointer->GetXaxis()->SetRangeUser(noiseborder,posMaxValHist);
         TF1* fitFunc = new TF1("fitFunc",GaussTail,noiseborder,posMaxValHist,6);
-        parameters = (Double_t *)calloc(9, sizeof(Double_t));  
+        parameters = (Double_t *)calloc(10, sizeof(Double_t));  
         const Double_t def_amplitude=459.951;
         const Double_t def_peakcenter=histogrampointer->GetBinCenter(histogrampointer->GetMaximumBin());
 //                  cout << colorcyan << "def_peakcenter: " << def_peakcenter << endlr;
@@ -553,6 +651,7 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType,
         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);
@@ -564,13 +663,16 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType,
         fitFunc->SetParameter(5,def_bgoffs);
         fitFunc->SetParLimits(5,0,def_bgoffs*4);
         
-        parameters = (Double_t *)calloc(9, sizeof(Double_t));
+        // TODO: This fix disables the background
+        fitFunc->FixParameter(4,0);
+       // fitFunc->FixParameter(5,0);
+        
         int fittries = 0;
         Bool_t failed = false;
         do {
             failed = false;
 //                            cout << fittries <<  ": New range: " <<  histogrampointer->GetXaxis()->GetBinCenter(min+fittries*min/20) << " to " << posMaxValHist << endl;
-            histogrampointer->Fit(fitFunc, "N,M,W,Q", "", noiseborder, posMaxValHist);
+            fit_result_ptr = histogrampointer->Fit(fitFunc, "N,M,W,Q,S", "", noiseborder, posMaxValHist);
 //             cout << colorcyan << " AFTER fit " << endlr;   
             if (gMinuit == nullptr) { 
                 failed = true;
@@ -588,6 +690,10 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType,
                 fitFunc->SetParameter(2,def_gausssig);
                 fitFunc->SetParameter(4,def_bgslope);
                 fitFunc->SetParameter(5,def_bgoffs);
+                
+                // TODO: This fix disables the background
+                fitFunc->FixParameter(4,0);
+//                 fitFunc->FixParameter(5,0);
             }
             else
                 fittries = 100;
@@ -603,7 +709,12 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType,
             fitFunc->SetParameter(3,def_distgauss);
             fitFunc->SetParameter(4,def_bgslope);
             fitFunc->SetParameter(5,def_bgoffs);
-            histogrampointer->Fit(fitFunc, "N,M,W,Q", "", noiseborder, posMaxValHist);
+            
+            // TODO: This fix disables the background
+            fitFunc->FixParameter(4,0);
+//             fitFunc->FixParameter(5,0);
+            
+            fit_result_ptr = histogrampointer->Fit(fitFunc, "N,M,W,Q,S", "", noiseborder, posMaxValHist);
             if (gMinuit == nullptr) { 
                 failed = true;
             } else 
@@ -617,19 +728,31 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType,
         if (failed)
         {
             cout << colorred << " Could not find " << histogrampointer->GetName() << " peak" << endlr;
-            parameters = (Double_t *)calloc(8, sizeof(Double_t));
+            parameters = (Double_t *)calloc(10, sizeof(Double_t));
             return parameters;
         }
 //         fitFunc->FixParameter(1,parameters[1]+histogrampointer->GetBinWidth(0));            
 //         histogrampointer->Fit(fitFunc, "N,M,W,Q", "", noiseborder, posMaxValHist);           
+
+        // parameters are 
+        // parameters[0]: amplitude of peak
+        // parameters[1]: x position of peak
+        // parameters[2]: sigma of gaussian
+        // parameters[6]: integral of peak
+        // parameters[7]: start of integral
+        // parameters[8]: end of integratio
+        // parameters[9]: error of integral
         for (Int_t pari=0; pari<6; pari++)
             parameters[pari]=fitFunc->GetParameter(pari);     
         //parameters[1] = parameters[1] - histogrampointer->GetBinWidth(0)*2;
         
         //parameters[7] = FindBorderToPeak(histogrampointer, noiseborder,def_peakcenter, verbose); // starting point of histogram integration
-        parameters[7] = parameters[1] - 2*parameters[2] ; // starting point of histogram integration
-        parameters[8] = parameters[1] + 2*parameters[2] ; // end point of histogram integration
-        parameters[6] = histogrampointer->Integral(histogrampointer->GetXaxis()->FindBin( parameters[1] - 2*parameters[2]  ), histogrampointer->GetXaxis()->FindBin( parameters[1] + 2*parameters[2]), "width"); // integral value of histogram (NOT fit)
+        parameters[7] = parameters[1] - 3*parameters[2] ; // starting point of histogram integration
+        parameters[8] = parameters[1] + 3*parameters[2] ; // end point of histogram integration
+        integralPeak  = histogrampointer->IntegralAndError(histogrampointer->GetXaxis()->FindBin( parameters[7]  ), histogrampointer->GetXaxis()->FindBin( parameters[8]), integralPeakError, "width"); // integral value of histogram (NOT fit),
+                
+        parameters[6]=integralPeak;
+        parameters[9]=integralPeakError;
         
         TF1 *bgfct = new TF1("f1","[0] +[1]*x",0,posMaxValHist);
         bgfct->SetParameters(parameters[5],parameters[4]);     
@@ -644,7 +767,10 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType,
             cout << colorcyan << "Integral value without bg: " << parameters[6] << endlr;              
 			cout << colorcyan << "Integral value with bg: " << parameters[6] + integralbg << endlr;          
 			cout << colorcyan << "Gauss Sigma: " << parameters[2] << endlr; 			
-            cout << colorcyan << "Number of events : " << frames_found << endlr;                             
+            cout << colorcyan << "Number of events : " << frames_found << endlr;     
+            
+             (fit_result_ptr.Get())->Print("V"); 
+             
         }
 
 //             DEBUG
@@ -657,7 +783,21 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType,
             fitFunc->Draw("SAME");
         
     }
+
+
+    fit_result_giveback = fit_result_ptr.Get();
+//     cout << "Correct address:" << endl;
+//     cout << fit_result_ptr.Get() << endl;
+//     
+//     fit_result_ptr_giveback = new TFitResultPtr();
+//     
+//     fit_result_ptr_giveback = &fit_result_ptr;
+//     cout << "Gave pointer address:" << endl;
+//     cout << fit_result_ptr_giveback->Get() << endl;
     
+//     TFitResult resultnew = new TFitResult(fit_result_ptr.Get());
+//     cout << "saved address:" << endl;
+//     cout << &resultnew << endl;
     return parameters;
 }
 
diff --git a/MABS_run_analyzer/HistogramType.h b/MABS_run_analyzer/HistogramType.h
index ba758fb..8a512cd 100644
--- a/MABS_run_analyzer/HistogramType.h
+++ b/MABS_run_analyzer/HistogramType.h
@@ -12,6 +12,7 @@
 #include <TMath.h>
 #include <TStyle.h>
 #include <TMinuit.h>
+#include <TFitResult.h>
 
 #include "help.h"
 
@@ -152,12 +153,16 @@ public:
     
     /// fitted position of the most probable value of the seed spectrum
     Float_t posSeed=0;
-    /// fintegral over the Seed peak, 4 sigma area
+    /// fintegral over the Seed peak, 2 sigma area
     Double_t integralSeed=0;
+    /// fintegral error over the Seed peak, poisson dist. assumed -> sqrt(N) for each bin
+    Double_t integralSeedErr=0;
     /// fitted position of the most probable value in the over cluster summed spectrum
     Float_t posSum=0;
     /// fintegral over the sum peak, normalized to number of events
     Double_t integralSum=0;
+    /// fintegral error over the Seed peak, poisson dist. assumed -> sqrt(N) for each bin
+    Double_t integralSumErr=0;
     /// fitted position of the calibration peak of Fe55-beta-photons in the seed spectrum, from a run best suited to the current
     Float_t posVeto=0;
     /// fintegral over the veto peak, normalized to number of events
@@ -258,6 +263,9 @@ public:
     /**
      * @brief rescales one specific histogram bin content from ADU to electrons */     
     void calibrate2DHistogramCounts(TH2F* &histogrampointernew, TH2F* &histogrampointerold );
+
+
+    
     /**
      * @brief intern function to calculate and plot fit curve to given histogram
      * 
@@ -265,7 +273,8 @@ public:
      * @return peak position of the fit
      * 
      */
-    Double_t* FitPerform(TH1F* histogrampointer, TString fitFuncType, Bool_t verbose = 0, Float_t fitstart = -1);
+    Double_t* FitPerform(TH1F* histogrampointer, TString fitFuncType, TFitResult* r = 0, Bool_t verbose = 0, Float_t fitstart = -1);
+    //Double_t* FitPerform(TH1F* histogrampointer, TString fitFuncType, TFitResultPtr* r = 0, Bool_t verbose = 0, Float_t fitstart = -1);
     /** @brief calculates Charge Collection Efficiency if Fe55 run */
     Bool_t calculteCCE(Bool_t verbose=false);
     /**
@@ -322,4 +331,4 @@ public:
     
 };
 
-#endif
\ No newline at end of file
+#endif
diff --git a/MABS_run_analyzer/MAPS.c b/MABS_run_analyzer/MAPS.c
index 6f96291..99a22bc 100644
--- a/MABS_run_analyzer/MAPS.c
+++ b/MABS_run_analyzer/MAPS.c
@@ -1584,7 +1584,8 @@ void MAPS::hitana() {
 //                     cout << endl;
 
                     // if cluster charge > clusternoise * const
-                    if(fOrderCode=="FSBB" || fOrderCode=="Mi19")
+                    //if(fOrderCode=="FSBB" || fOrderCode=="Mi19") // if you need CCE_4, uncomment this and comment the line below for Mi19
+                    if(fOrderCode=="FSBB")
                     {
                         if (1.0*chargesumincluster4 > noisesumincluster4*6.0)
                             fFrameInfo.pixelthreshold[fHits] = Hitlist[hit];
diff --git a/MABS_run_analyzer/MAPS.h b/MABS_run_analyzer/MAPS.h
index 53c8c83..10939a0 100644
--- a/MABS_run_analyzer/MAPS.h
+++ b/MABS_run_analyzer/MAPS.h
@@ -498,4 +498,4 @@ public:
     /// Size of cluster. Don't change without code modification
     const Int_t clustersize = 5;
 };
-#endif
\ No newline at end of file
+#endif
diff --git a/MABS_run_analyzer/Run.c b/MABS_run_analyzer/Run.c
index a29e976..d9399a1 100644
--- a/MABS_run_analyzer/Run.c
+++ b/MABS_run_analyzer/Run.c
@@ -91,6 +91,7 @@ Run::Run(Int_t runnumber, Int_t loopi, TString customtitle)
             } else {
                 storepathRAWLinux=Form("/d/garlic/%s/%d",labbook.chipGen.Data(),runnumber); // default empty path
             }
+            
             storepathRootLinux = Form("/d/garlic/%s/rootFiles",labbook.chipGen.Data()); // default empty path
             labbook.posVetoDB = (rowsql->GetField(10) != NULL)?atof(rowsql->GetField(10)):-1;
             labbook.posSeedDB = (rowsql->GetField(11) != NULL)?atof(rowsql->GetField(11)):-1;
@@ -176,10 +177,10 @@ Run::Run(Int_t runnumber, Int_t loopi, TString customtitle)
             HistogramClassVector.push_back(histogramwoRTSthreshold);
             
             
-            //             histogramwoRTSAggresive = new HistogramType(" more RTS cleaned", &cursystemparam, &cursensorinfo, humanreadablestr, &labbook, this, rootcolors[plotStyle], rootlinestyle[plotStyle]);
-//             histogramwoRTSAggresive->maskRTSpixel = true;
-//             histogramwoRTSAggresive->RTSthreshold = 1.5;
-//             HistogramClassVector.push_back(histogramwoRTSAggresive);
+             histogramwoRTSAggresive = new HistogramType(" more RTS cleaned", &cursystemparam, &cursensorinfo, humanreadablestr, &labbook, this, rootcolors[plotStyle], rootlinestyle[plotStyle]);
+             histogramwoRTSAggresive->maskRTSpixel = true;
+             histogramwoRTSAggresive->RTSthreshold = 1.5;
+             HistogramClassVector.push_back(histogramwoRTSAggresive);
             //             histogramwoRTSAggresivethreshold = new HistogramType(" Threshold, more RTS cleaned", &cursystemparam, &cursensorinfo, humanreadablestr, &labbook, this, rootcolors[plotStyle], rootlinestyle[plotStyle] );
 //             histogramwoRTSAggresivethreshold->maskRTSpixel = true;
 //             histogramwoRTSAggresivethreshold->RTSthreshold = 1.0;
@@ -254,7 +255,7 @@ void Run::setSystemSpecificParameters()
     systemparam systemparamPXI (800*16,800,75,150,150);
     systemparam systemparamFSBB (2800,2800/4,25,10,100);
     systemparam systemparamUSBMi19 (400/*maxbin*/,400/1/*nbins*/, 25/*vetothreshold*/, 10/*maxbinnoise*/, 100/*nbinsnoise*/);
-    systemparam systemparamPipper2 (1000,500,15,10,100);
+    systemparam systemparamPipper2 (1000,200,15,10,100);
     if (labbook.system.EqualTo("USB")  && labbook.chipGen.EqualTo("Mi34") )
         cursystemparam = systemparamUSB;
     else if (labbook.system.EqualTo("USB")  && labbook.chipGen.EqualTo("FSBB") )
@@ -527,6 +528,7 @@ Bool_t Run::FindRTSPixelToMask(HistogramType* HistogramTypepointer)
         }
     }
     
+    
     // calculate the mean firing times and bin fireing times in a propability histogram, first approach
     for (Int_t pixeli=0; pixeli < HistogramTypepointer->pixeltimefired->GetNbinsX(); pixeli++) {
         if (HistogramTypepointer->pixeltimefired->GetBinContent(pixeli) > 0) {
@@ -536,10 +538,58 @@ Bool_t Run::FindRTSPixelToMask(HistogramType* HistogramTypepointer)
 //             HistogramTypepointer->pixeltimefiredsorted->Fill(2000);
         }
     }   
-//     TH1F* pixeltimefiredWithOverflow = ShowOverflow(HistogramTypepointer->pixeltimefiredsorted);
+    meanpixeltimesfired /= numberofactivepixel; // Very rough estimate of a mean firing time
+    
+    // vrey rough estimate on standard deviation
+    for (Int_t pixeli=0; pixeli < HistogramTypepointer->pixeltimefired->GetNbinsX(); pixeli++) {
+        if (HistogramTypepointer->pixeltimefired->GetBinContent(pixeli) > 0)
+            stdeviation += pow(HistogramTypepointer->pixeltimefired->GetBinContent(pixeli)-meanpixeltimesfired,2);
+    }
+    stdeviation /= numberofactivepixel;
+    stdeviation = sqrt(stdeviation);
+    cout << colorcyan << "Mean firing rate first approximation: " << meanpixeltimesfired << endlr;
+    cout << colorcyan << "Deviation value first approximation: " << stdeviation << endlr;
+    cout << colorcyan << "number of considered pixel: " << numberofactivepixel << endlr;
+    cout << colorcyan << endlr;        
+    
+    // better estimate on average firing time
+    Double_t meanpixeltimesfired2, stdeviation2;
+    uint numberofactivepixel2;
+    Bool_t RTSpixel;
+    u_int numberofconsideredpixel;
+    u_int poscounter;
+    
+    
+    for (UInt_t i=0; i < 3; ++i) // 3 times remove most fired pixel
+    {
+        meanpixeltimesfired2 = 0;
+        numberofactivepixel2=0;
+        HistogramTypepointer->RTSpixel.clear();
+        cout << "i " << i << endl;
+        cout << "HistogramTypepointer->pixeltimefired->GetBinContent(0) " << HistogramTypepointer->pixeltimefired->GetBinContent(0) << endl;
+        cout << "meanpixeltimesfired " << meanpixeltimesfired << endl;
+        cout << "HistogramTypepointer->RTSthreshold " << HistogramTypepointer->RTSthreshold << endl;
+        cout << "stdeviation " << stdeviation << endl;
+        for (Int_t pixeli=0; pixeli < HistogramTypepointer->pixeltimefired->GetNbinsX(); pixeli++) {
+            if (abs(HistogramTypepointer->pixeltimefired->GetBinContent(pixeli)-meanpixeltimesfired) > HistogramTypepointer->RTSthreshold * stdeviation) {
+                HistogramTypepointer->RTSpixel.push_back(pixeli);
+                RTSpixelHits+=HistogramTypepointer->pixeltimefired->GetBinContent(pixeli);
+            }
+            else
+            {
+                meanpixeltimesfired2 += HistogramTypepointer->pixeltimefired->GetBinContent(pixeli);
+                numberofactivepixel2++;
+            }
+            
+        }
+        if (numberofactivepixel2 > 0) {
+            meanpixeltimesfired2 /= numberofactivepixel2;
+            meanpixeltimesfired = meanpixeltimesfired2;
+        }
     
 //     Double_t const probabilities[] = {1-erf(HistogramTypepointer->RTSthreshold/sqrt(2)), 0.5, erf(HistogramTypepointer->RTSthreshold/sqrt(2))}; // sigma/2 from gaus to the left and to the right //{0.17, 0.5, 1-0.17};     
 //     HistogramTypepointer->pixeltimefiredsorted->GetQuantiles( 3, leakagequantiles, probabilities);
+    
 //     HistogramTypepointer->medianLeakageCurrent = leakagequantiles[1];
 //     HistogramTypepointer->medianLeakageCurrentPlus = leakagequantiles[2] - leakagequantiles[1];
 //     HistogramTypepointer->medianLeakageCurrentMinus = leakagequantiles[1] - leakagequantiles[0];
@@ -553,92 +603,84 @@ Bool_t Run::FindRTSPixelToMask(HistogramType* HistogramTypepointer)
 //     cout << "last value above 1 in overflow TH1F: " << pixeltimefiredWithOverflow->GetBinCenter(pixeltimefiredWithOverflow->FindLastBinAbove(4)) << endl;
     //     cout << colorcyan << "Deviation value after one more loop: " << stdeviation << endlr;
     
-    HistogramTypepointer->RTSpixel.clear();
-    for (Int_t pixeli=0; pixeli < HistogramTypepointer->pixeltimefired->GetNbinsX(); pixeli++) {
-        if (HistogramTypepointer->pixeltimefired->GetBinContent(pixeli) > leakagequantiles[2]) {
-//             cout << coloryellow << numberToString(5) << " mal drüber, Böser pixel #" << pixeli << " Wert: " << HistogramTypepointer->pixeltimefired->GetBinContent(pixeli) << endlr;
-            HistogramTypepointer->RTSpixel.push_back(pixeli);
-        }
-    }
-//     cout << colorred << " number of RTS pixel in quantile run: " << HistogramTypepointer->RTSpixel.size() << endlr;
+//     HistogramTypepointer->RTSpixel.clear();
+//     for (Int_t pixeli=0; pixeli < HistogramTypepointer->pixeltimefired->GetNbinsX(); pixeli++) {
+//         if (HistogramTypepointer->pixeltimefired->GetBinContent(pixeli) > leakagequantiles[2]) {
+// //              cout << coloryellow << numberToString(5) << " mal drüber, Böser pixel #" << pixeli << " Wert: " << HistogramTypepointer->pixeltimefired->GetBinContent(pixeli) << endlr;
+//             HistogramTypepointer->RTSpixel.push_back(pixeli);
+//         }
+//     }
+//      cout << colorred << " number of RTS pixel in quantile run: " << HistogramTypepointer->RTSpixel.size() << endlr;
     
     
-    meanpixeltimesfired /= numberofactivepixel; // Very rough estimate of a mean firing time
+//     meanpixeltimesfired /= HistogramTypepointer->pixeltimefired->GetNbinsX(); // Very rough estimate of a mean firing time
 //     meanpixeltimesfired = leakagequantiles[1];
 //     cout << colorcyan << "Mean value own: " << meanpixeltimesfired << endlr;
         
-    // vrey rough estimate on standard deviation
-    for (Int_t pixeli=0; pixeli < HistogramTypepointer->pixeltimefired->GetNbinsX(); pixeli++) {
-        if (HistogramTypepointer->pixeltimefired->GetBinContent(pixeli) > 0)
-            stdeviation += pow(HistogramTypepointer->pixeltimefired->GetBinContent(pixeli)-meanpixeltimesfired,2);
-    }
-    stdeviation /= numberofactivepixel;
-    stdeviation = sqrt(stdeviation);
-//     cout << colorcyan << "Deviation value own 1: " << stdeviation << endlr;
-    
-    // better estimate on average firing time
-    Double_t meanpixeltimesfired2 = 0;
-    uint numberofactivepixel2=0;
-    HistogramTypepointer->RTSpixel.clear();
-    for (Int_t pixeli=0; pixeli < HistogramTypepointer->pixeltimefired->GetNbinsX(); pixeli++) {
-        if (abs(HistogramTypepointer->pixeltimefired->GetBinContent(pixeli)-meanpixeltimesfired) > 2 * stdeviation) {
-            HistogramTypepointer->RTSpixel.push_back(pixeli);
-        }
-        else
-        {
-            meanpixeltimesfired2 += HistogramTypepointer->pixeltimefired->GetBinContent(pixeli);
-            numberofactivepixel2++;
-        }
-    }
-    meanpixeltimesfired2 /= numberofactivepixel2;
-//     cout << colorred << " number of RTS pixel in first run: " << HistogramTypepointer->RTSpixel.size() << endlr;
-//     cout << colorcyan << " mean value new: " << meanpixeltimesfired2 << endl;
-    
-//  Estimate new standard deviation    
-    Bool_t RTSpixel =false;
-    u_int numberofconsideredpixel=0;
-    u_int poscounter = 0;
-    stdeviation = 0;
-    for (u_int pixeli=0; (int)pixeli<cursensorinfo.columns*cursensorinfo.rows ; pixeli++) {
-        if (HistogramTypepointer->pixeltimefired->GetBinContent(pixeli) > 0) {
-            RTSpixel = false;
-            for (u_int RTSpixeli=poscounter; RTSpixeli < HistogramTypepointer->RTSpixel.size(); RTSpixeli++) {
-                if (pixeli == HistogramTypepointer->RTSpixel[RTSpixeli]) {
-                    poscounter = RTSpixeli;
-                    RTSpixel = true;
+    //  Estimate new standard deviation    
+        RTSpixel =false;
+        numberofconsideredpixel=0;
+        poscounter = 0;
+        stdeviation2 = 0;
+        for (UInt_t pixeli=0; pixeli < (UInt_t)HistogramTypepointer->pixeltimefired->GetNbinsX(); pixeli++) {
+            if (HistogramTypepointer->pixeltimefired->GetBinContent(pixeli) > 0) {
+                RTSpixel = false;
+                for (u_int RTSpixeli=0; RTSpixeli < (u_int) HistogramTypepointer->RTSpixel.size(); RTSpixeli++) {
+                    if (pixeli == HistogramTypepointer->RTSpixel[RTSpixeli]) {
+                        poscounter = RTSpixeli;
+                        RTSpixel = true;
+                    }
+                }
+                if (!RTSpixel) {
+                    stdeviation2 += pow(HistogramTypepointer->pixeltimefired->GetBinContent(pixeli)-meanpixeltimesfired2,2);
+                    numberofconsideredpixel++;
                 }
             }
-            if (!RTSpixel) {
-                stdeviation += pow(HistogramTypepointer->pixeltimefired->GetBinContent(pixeli)-meanpixeltimesfired2,2);
-                numberofconsideredpixel++;
-            }
-        }
-    }
-    stdeviation /= numberofconsideredpixel;
-    stdeviation = sqrt(stdeviation); 
-//     cout << colorcyan << "Deviation value own 2: " << stdeviation << endlr;
-    
-    
-    // better estimate on RTS pixel
-    meanpixeltimesfired = 0;
-    numberofconsideredpixel = 0;
-    HistogramTypepointer->RTSpixel.clear();
-    for (Int_t pixeli=0; pixeli < HistogramTypepointer->pixeltimefired->GetNbinsX(); pixeli++) {
-        totalHits+=HistogramTypepointer->pixeltimefired->GetBinContent(pixeli);
-        if (HistogramTypepointer->pixeltimefired->GetBinContent(pixeli)-meanpixeltimesfired2 > HistogramTypepointer->RTSthreshold * stdeviation) {
-            HistogramTypepointer->RTSpixel.push_back(pixeli);
-            RTSpixelHits+=HistogramTypepointer->pixeltimefired->GetBinContent(pixeli);
-//             cout << coloryellow << numberToString(HistogramTypepointer->RTSthreshold) << " mal drüber, Böser pixel #" << pixeli << " Wert: " << HistogramTypepointer->pixeltimefired->GetBinContent(pixeli) << endlr;
         }
-        else
-        {
-            meanpixeltimesfired += HistogramTypepointer->pixeltimefired->GetBinContent(pixeli);
-            numberofconsideredpixel++;
+        if (numberofconsideredpixel > 0) {
+            stdeviation2 /= numberofconsideredpixel;
+            stdeviation2 = sqrt(stdeviation); 
+            stdeviation = stdeviation2;
         }
+        cout << colorred << " number of RTS pixel in " << i << "  run: " << HistogramTypepointer->RTSpixel.size() << endlr;
+        cout << colorred << " number of considered pixel: " << numberofconsideredpixel << endlr;
+        cout << colorcyan << " mean value: " << meanpixeltimesfired << endl;
+        cout << colorcyan << " Deviation value: " << stdeviation << endlr;
+        
+        
+        
+        //     Double_t const probabilities[] = {1-erf(HistogramTypepointer->RTSthreshold/sqrt(2)), 0.5, erf(HistogramTypepointer->RTSthreshold/sqrt(2))}; // sigma/2 from gaus to the left and to the right //{0.17, 0.5, 1-0.17};     
+        //     HistogramTypepointer->pixeltimefiredsorted->GetQuantiles( 3, leakagequantiles, probabilities);
+        
+        //     HistogramTypepointer->medianLeakageCurrent = leakagequantiles[1];
+        //     HistogramTypepointer->medianLeakageCurrentPlus = leakagequantiles[2] - leakagequantiles[1];
+        //     HistogramTypepointer->medianLeakageCurrentMinus = leakagequantiles[1] - leakagequantiles[0];
+        //     cout << "Median is at: " << HistogramTypepointer->medianLeakageCurrent << endl;
+        //     cout << "Upper limit is at: " << leakagequantiles[2] << endl;
+        //     cout << "Lower limit is at: " << leakagequantiles[0] << endl;
+        //     cout << "Plus: " << HistogramTypepointer->medianLeakageCurrentPlus << endl;
+        //     cout << "minus: " << HistogramTypepointer->medianLeakageCurrentMinus << endl;
+        
+        //     cout << "last value above 1: " << HistogramTypepointer->pixeltimefiredsorted->GetBinCenter(HistogramTypepointer->pixeltimefiredsorted->FindLastBinAbove(4)) << endl; 
+        //     cout << "last value above 1 in overflow TH1F: " << pixeltimefiredWithOverflow->GetBinCenter(pixeltimefiredWithOverflow->FindLastBinAbove(4)) << endl;
+        //     cout << colorcyan << "Deviation value after one more loop: " << stdeviation << endlr;
+        
+        //     HistogramTypepointer->RTSpixel.clear();
+        //     for (Int_t pixeli=0; pixeli < HistogramTypepointer->pixeltimefired->GetNbinsX(); pixeli++) {
+        //         if (HistogramTypepointer->pixeltimefired->GetBinContent(pixeli) > leakagequantiles[2]) {
+        // //              cout << coloryellow << numberToString(5) << " mal drüber, Böser pixel #" << pixeli << " Wert: " << HistogramTypepointer->pixeltimefired->GetBinContent(pixeli) << endlr;
+        //             HistogramTypepointer->RTSpixel.push_back(pixeli);
+        //         }
+        //     }
+        //      cout << colorred << " number of RTS pixel in quantile run: " << HistogramTypepointer->RTSpixel.size() << endlr;
+        
+        
+        //     meanpixeltimesfired /= HistogramTypepointer->pixeltimefired->GetNbinsX(); // Very rough estimate of a mean firing time
+        //     meanpixeltimesfired = leakagequantiles[1];
+        //     cout << colorcyan << "Mean value own: " << meanpixeltimesfired << endlr;
     }
-    meanpixeltimesfired /= numberofconsideredpixel;
-//     cout << colorred << " number of RTS pixel in second run: " << HistogramTypepointer->RTSpixel.size() << endlr;
-//     cout << colorcyan << " mean value new new: " << meanpixeltimesfired << endl;
+    
+ 
 
     HistogramTypepointer->percentageofRTSpixel = (HistogramTypepointer->RTSpixel.size()*1.0)/(numberofactivepixel*1.0)*100.0;
     cout << "  cutted away evertyhing with more then " << std::right<< Form("%.1f",HistogramTypepointer->RTSthreshold * stdeviation+meanpixeltimesfired2) << " hits" << endlr;
@@ -963,9 +1005,7 @@ void Run::updateDatabase() {
 	} else {
 		histogramclassToUseForDB = histogramwoRTSthreshold;
 	}
-		
-	
-	
+						
     string sqlupdatequery = "";	
     constructUpdateString(&sqlupdatequery, "Gain",       histogramclassToUseForDB->normalized->gain);
     constructUpdateString(&sqlupdatequery, "SumPeak",    histogramclassToUseForDB->normalized->posSum, 4);
@@ -973,9 +1013,11 @@ void Run::updateDatabase() {
     constructUpdateString(&sqlupdatequery, "AvgF0",   labbook.averageF0, 6);
     constructUpdateString(&sqlupdatequery, "SigmaF0",   labbook.sigmaF0, 6);
     constructUpdateString(&sqlupdatequery, "SeedIntegralMinusBG2Sigma",   histogramclassToUseForDB->normalized->integralSeed, 10);
+    constructUpdateString(&sqlupdatequery, "SeedIntegralMinusBG2SigmaErr",   histogramclassToUseForDB->normalized->integralSeedErr, 5);
     constructUpdateString(&sqlupdatequery, "VetoPeak",   histogramclassToUseForDB->normalized->posVeto, 4);
     constructUpdateString(&sqlupdatequery, "VetoIntegral",   histogramclassToUseForDB->normalized->integralVeto, 10);
     constructUpdateString(&sqlupdatequery, "SumIntegral",   histogramclassToUseForDB->normalized->integralSum, 10);
+    constructUpdateString(&sqlupdatequery, "SumIntegralErr",   histogramclassToUseForDB->normalized->integralSumErr, 5);
     if (histogramclassToUseForDB->normalized->calibrated != 0)
 //         if (histogramthreshold->normalized->calibrated->avgNoise < 100)
             constructUpdateString(&sqlupdatequery, "Avg.Noise",  histogramclassToUseForDB->normalized->calibrated->avgNoise);
@@ -1039,7 +1081,7 @@ string Run::to_str_w_prec(const Float_t a_value, int precision = 3)
 Bool_t Run::binLeakageCurrent(HistogramType* oneHistogramClass)
 {
     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 const probabilities[] = {0.25, 0.5, 0.75}; // sigma/2 from gaus to the left and to the right //{0.17, 0.5, 1-0.17}; 
+//     Double_t const probabilities[] = {0.25, 0.5, 0.75}; // sigma/2 from gaus to the left and to the right, an area of 68% must be within the interesting range //{0.17, 0.5, 1-0.17}; 
     Double_t pedestals [cursensorinfo.columns*cursensorinfo.rows];
     for (Int_t pixeli=0; pixeli<cursensorinfo.columns*cursensorinfo.rows ; pixeli++) { // loop over all pixel
         pedestals[pixeli]=0;
@@ -1243,7 +1285,8 @@ Bool_t Run::binSeedSumVeto()
                     // sum histogram
                     pixelSum = 0;
                     notSeedSum = 0;
-                    if(labbook.chipGen=="FSBB" || labbook.chipGen=="Mi19")
+                    //if(labbook.chipGen=="FSBB" || labbook.chipGen=="Mi19")
+                    if(labbook.chipGen=="FSBB")
                     {	
                         Float_t clusterArray[processed->clustersize*processed->clustersize];// temp variable clusterArray necessary, because Sort only accepts 1-dim arrays
                         Int_t index[processed->clustersize*processed->clustersize];
@@ -1347,7 +1390,7 @@ Bool_t Run::binSeedSumVeto()
                             if (processed->fFrameInfo.pixel[hiti] == histogramwoRTS->RTSpixel[RTSpixeli])
                             {
                                 RTSpixel = true;
-//                                 cout << "not binned! RTS pixel #" << processed->fFrameInfo.pixel[hiti] << " with charge: " << processed->fFrameInfo.p[12][hiti] << endl;
+//                                  cout << "not binned! RTS pixel #" << processed->fFrameInfo.pixel[hiti] << " with charge: " << processed->fFrameInfo.p[12][hiti] << endl;
                             }
                         }
                         if (!RTSpixel) {
@@ -1365,10 +1408,12 @@ Bool_t Run::binSeedSumVeto()
                             if (processed->fFrameInfo.pixel[hiti] == histogramwoRTSAggresive->RTSpixel[RTSpixeli])
                             {
                                 RTSpixel = true;
+//                                 cout << "not binned! RTS pixel #" << processed->fFrameInfo.pixel[hiti] << " with charge: " << processed->fFrameInfo.p[12][hiti] << endl;
                                 //                                 cout << "not binned! RTS pixel #" << processed->fFrameInfo.pixel[hiti] << " with charge: " << processed->fFrameInfo.p[12][hiti] << endl;
                             }
                         }
                         if (!RTSpixel) {
+//                             cout << "binned! pixel #" << processed->fFrameInfo.pixel[hiti] << " with charge: " << processed->fFrameInfo.p[12][hiti] << endl;
                             histogramwoRTSAggresive->Seed->Fill(processed->fFrameInfo.p[12][hiti]);
                             histogramwoRTSAggresive->Sum->Fill(pixelSum);
                             histogramwoRTSAggresive->SeedPerc->Fill(processed->fFrameInfo.p[12][hiti]/pixelSum*100);  
@@ -1380,9 +1425,10 @@ Bool_t Run::binSeedSumVeto()
             } // end loop over hits in frame
         }
     } // end loop over all frames
+    cout << colorred << " histogramwoRTSAggresive->Seed size: " << histogramwoRTSAggresive->Seed->GetEntries() << endlr;
     //     gROOT->SetBatch(kTRUE);
     for (vector<HistogramType*>::iterator curHistogramClass = HistogramClassVector.begin(); curHistogramClass != HistogramClassVector.end(); curHistogramClass++)  {
-        Double_t* parameters = (Double_t *)calloc(10, sizeof(Double_t)); // allocate 10 parameters for safety, maximum 9 used at the moment     
+        Double_t* parameters = (Double_t *)calloc(11, sizeof(Double_t)); // allocate 11 parameters for safety, maximum 10 used at the moment     
         (*curHistogramClass)->FindNoisethresholdborder((*curHistogramClass)->Seed, false, false);
         if (labbook.runnumber == 343056)
             (*curHistogramClass)->noisethresholdborder = 34;
@@ -1403,14 +1449,16 @@ Bool_t Run::binSeedSumVeto()
             (*curHistogramClass)->integralVeto = parameters[6];
         }
         if (labbook.chipGen.EqualTo("Pipper2")) 
-            parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->Seed, "GaussTail");
+            parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->Seed, "GaussTail", 0, false, 500);
         else
             parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->Seed, "landau");
         (*curHistogramClass)->integralSeed = parameters[6];
         (*curHistogramClass)->posSeed = parameters[1];
-        parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->Sum, "gaus");
+        (*curHistogramClass)->integralSeedErr = parameters[9];
+        parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->Sum, "GaussTail", 0, false, 500); //TODO change back to gauss
         (*curHistogramClass)->posSum = parameters[1];
-        (*curHistogramClass)->integralSum = parameters[3];
+        (*curHistogramClass)->integralSum = parameters[6];
+        (*curHistogramClass)->integralSumErr = parameters[9];
         
         parameters =(*curHistogramClass)->FitPerform((*curHistogramClass)->SeedPerc, "landau");
         (*curHistogramClass)->posSeedPerc = parameters[1];
@@ -2057,13 +2105,13 @@ Bool_t Run::plotClusterDistribution(HistogramType* HistogramTypepointer)
     
     // create lorentz fit of a slice
     TH1F* xslicetroughcluster = new TH1F(Form("%s_xslice", HistogramTypepointer->histAvgCluster->GetTitle()),"X slice",HistogramTypepointer->histAvgCluster->GetXaxis()->GetNbins(),HistogramTypepointer->histAvgCluster->GetXaxis()->GetBinLowEdge(1),HistogramTypepointer->histAvgCluster->GetXaxis()->GetBinUpEdge(HistogramTypepointer->histAvgCluster->GetXaxis()->GetNbins()));
-    Double_t* xslicetroughclusterpar = (Double_t *)calloc(4, sizeof(Double_t)); // allocate 6 parameters for safety, maximum 4 used at the moment 
+    Double_t* xslicetroughclusterpar = (Double_t *)calloc(11, sizeof(Double_t)); // allocate 11 parameters for safety, maximum 10 used at the moment 
     Int_t middlebin = (int)(HistogramTypepointer->histAvgCluster->GetXaxis()->GetNbins()/2.0 + 0.5);
     for (Int_t bini=1; bini <= HistogramTypepointer->histAvgCluster->GetXaxis()->GetNbins(); bini++)
         xslicetroughcluster->SetBinContent(bini,HistogramTypepointer->histAvgCluster->GetBinContent(bini,middlebin));
     
     TH1F* yslicetroughcluster = new TH1F(Form("%s_yslice", HistogramTypepointer->histAvgCluster->GetTitle()),"Y slice",HistogramTypepointer->histAvgCluster->GetYaxis()->GetNbins(),HistogramTypepointer->histAvgCluster->GetYaxis()->GetBinLowEdge(1),HistogramTypepointer->histAvgCluster->GetYaxis()->GetBinUpEdge(HistogramTypepointer->histAvgCluster->GetYaxis()->GetNbins()));
-    Double_t* yslicetroughclusterpar = (Double_t *)calloc(4, sizeof(Double_t)); // allocate 6 parameters for safety, maximum 4 used at the moment 
+    Double_t* yslicetroughclusterpar = (Double_t *)calloc(11, sizeof(Double_t)); // allocate 11 parameters for safety, maximum 10 used at the moment 
     middlebin = (int)(HistogramTypepointer->histAvgCluster->GetYaxis()->GetNbins()/2.0 + 0.5);
     for (Int_t bini=1; bini <= HistogramTypepointer->histAvgCluster->GetYaxis()->GetNbins(); bini++)
         yslicetroughcluster->SetBinContent(bini,HistogramTypepointer->histAvgCluster->GetBinContent(middlebin,bini));
@@ -2111,6 +2159,10 @@ TCanvas* Run::plot1DHistogram(HistogramType* HistogramTypepointer, TH1F* onehist
     if (onehistogram != nullptr)
     {
         Int_t random = random1->Rndm()*100000;
+        TFitResult* fit_result;
+//         TFitResultPtr* fit_result_ptr;
+        fit_result = new TFitResult();
+//         fit_result_ptr = new TFitResultPtr();
         TString canvastitle = Form("%s %s", onehistogram->GetName(), runcode.Data());
         TString canvasname =  Form("%s %s %d", onehistogram->GetName(), runcode.Data(), random);
         TCanvas* canvas = new TCanvas(canvasname, canvastitle, 900, 700);
@@ -2123,11 +2175,11 @@ TCanvas* Run::plot1DHistogram(HistogramType* HistogramTypepointer, TH1F* onehist
             onehistogram=ShowOverflow(onehistogram);
         onehistogram->Draw();
         
-        Double_t* parameters = (Double_t *)calloc(10, sizeof(Double_t)); // allocate 10 parameters for safety, maximum 9 used at the moment     
+        Double_t* parameters = (Double_t *)calloc(11, sizeof(Double_t)); // allocate 10 parameters for safety, maximum 9 used at the moment     
         Float_t maxValue = 0;
         if (fitFuncType!="")
         {
-            parameters = HistogramTypepointer->FitPerform(onehistogram, fitFuncType, verbose, fitstart);
+            parameters = HistogramTypepointer->FitPerform(onehistogram, fitFuncType, fit_result, verbose, fitstart);
             maxValue = parameters[1];
             plotVerticalLine(onehistogram, maxValue);
         }
@@ -2152,7 +2204,7 @@ TCanvas* Run::plot1DHistogram(HistogramType* HistogramTypepointer, TH1F* onehist
 //         canvas->Draw();
 //                 canvas->Update();
                 
-        if (fitFuncType=="GaussTail")
+        if (fitFuncType=="GaussTail" || fitFuncType=="gaus")
         {
             Float_t integralstart = parameters[7];
             Float_t integralend = parameters[8];
@@ -2171,22 +2223,27 @@ TCanvas* Run::plot1DHistogram(HistogramType* HistogramTypepointer, TH1F* onehist
                 cout << "Distance under   : " <<  parameters[3] << endl;
                 cout << "background slope : " <<  parameters[4] << endl;
                 cout << "background offset: " <<  parameters[5] << endl;
+                cout << "FWHM             : " <<  parameters[2]*sqrt(8*log(2)) << endl;
+//                 cout << "got address2:" << endl;
+//                 cout <<  fit_result << endl;
+//                  fit_result->Print("V"); 
+                
                 //TF1 *f1 = new TF1("f1","[0] +[1]*x +gaus(2)",0,5);
                 TF1 *f1 = new TF1("f1","[0] +[1]*x",0,1000);
                 f1->SetParameters(parameters[5],parameters[4]);
                 f1->Draw("SAME");          
             }
         }      
-        if (fitFuncType=="gaus")
-        {
-            Float_t integralstart = parameters[4];
-            Float_t integralval = parameters[3];
-            TString integrallbl = Form("Integral: %.0f", integralval);
-//             cout << colorcyan << "  " << HistogramTypepointer->histogramdescription << endlr;
-//             cout << colorcyan << ": " << parameters[4] << endlr;
-            plotVerticalLineHeigher(onehistogram, integralstart);
-            plotTextAtVal(onehistogram, maxValue, integrallbl);
-        }
+//         if (fitFuncType=="gaus")
+//         {
+//             Float_t integralstart = parameters[4];
+//             Float_t integralval = parameters[3];
+//             TString integrallbl = Form("Integral: %.0f", integralval);
+// //             cout << colorcyan << "  " << HistogramTypepointer->histogramdescription << endlr;
+// //             cout << colorcyan << ": " << parameters[4] << endlr;
+//             plotVerticalLineHeigher(onehistogram, integralstart);
+//             plotTextAtVal(onehistogram, maxValue, integrallbl);
+//         }
                 
 //         canvas -> SaveAs( savepathresults + "/" + runcode + " " + onehistogram->GetName() + ".eps");
             
diff --git a/MABS_run_analyzer/Run.h b/MABS_run_analyzer/Run.h
index f127117..1703891 100644
--- a/MABS_run_analyzer/Run.h
+++ b/MABS_run_analyzer/Run.h
@@ -22,6 +22,7 @@
 #include <TH2F.h>
 #include <TF1.h>
 #include <TImageDump.h>
+#include <TFitResult.h>
 
 #include <TMath.h>
 #include <TF1.h>
diff --git a/MABS_run_analyzer/help.h b/MABS_run_analyzer/help.h
index 12a44e8..a396ae7 100644
--- a/MABS_run_analyzer/help.h
+++ b/MABS_run_analyzer/help.h
@@ -603,4 +603,4 @@ void *sortThread(void *ptr)
 }
 
 //####################################################################
-#endif
\ No newline at end of file
+#endif