From bf099d86bdd34ca91b128dbe08deed17c36a19c3 Mon Sep 17 00:00:00 2001 From: Benjamin Linnik Date: Thu, 23 Nov 2017 15:37:48 +0100 Subject: [PATCH] Analyzer: Bugfix, calibration was broken (by me) --- MABS_run_analyzer/ChargeSpektrum.c | 25 +++++++++++--- MABS_run_analyzer/HistogramType.c | 55 +++++++++++++++++++++++------- MABS_run_analyzer/HistogramType.h | 5 +++ MABS_run_analyzer/Run.c | 29 ++++++++++++---- MABS_run_analyzer/Run.h | 1 + 5 files changed, 92 insertions(+), 23 deletions(-) diff --git a/MABS_run_analyzer/ChargeSpektrum.c b/MABS_run_analyzer/ChargeSpektrum.c index a30464c..b24821b 100644 --- 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 compareHistogramClassVectorClassic; + // variables for cluster threshold analysis + vector compareHistogramClassVectorClusterThreshold; + // variables for calibrated analysis vector compareHistogramClassVectorCalibrated; @@ -80,6 +83,9 @@ void ChargeSpektrum(TString runnumber = "") vector compareHistogramVectorClusterMultiplicity; vector ntupleVectorAvgChargeColl; + // f0 analysis variables + vector compareHistogramAllRunsVectorF0; + for(Int_t runi=0; runi0) @@ -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 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 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 7517852..38df915 100644 --- 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 97f0e4c..4434506 100644 --- 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 074af13..93faba2 100644 --- 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 f0ac8cf..ad48f65 100644 --- 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 * -- 2.43.0