From ba441a179cfe7f4c20ae44b04b2dc70494e50e47 Mon Sep 17 00:00:00 2001 From: Benjamin Linnik Date: Mon, 21 Sep 2015 11:56:17 +0200 Subject: [PATCH] Run analyzer: Added the new class HistogramType.c, now all histograms and cuts are stored and managed in this class. Every cut is an instance of the class HistogramType. --- MABS_run_analyzer/ChargeSpektrumFunctions.c | 71 +- MABS_run_analyzer/HistogramType.c | 419 +++++++++ MABS_run_analyzer/HistogramType.h | 202 +++++ MABS_run_analyzer/MAPS.c | 52 +- MABS_run_analyzer/MAPS.h | 9 - MABS_run_analyzer/Run.c | 954 +++++--------------- MABS_run_analyzer/Run.h | 260 +----- MABS_run_analyzer/help.h | 43 + 8 files changed, 992 insertions(+), 1018 deletions(-) create mode 100644 MABS_run_analyzer/HistogramType.c create mode 100644 MABS_run_analyzer/HistogramType.h diff --git a/MABS_run_analyzer/ChargeSpektrumFunctions.c b/MABS_run_analyzer/ChargeSpektrumFunctions.c index 01de8da..d529180 100644 --- a/MABS_run_analyzer/ChargeSpektrumFunctions.c +++ b/MABS_run_analyzer/ChargeSpektrumFunctions.c @@ -15,12 +15,7 @@ Bool_t writeObservableToFile(TString); Bool_t writeObservableToFile(); void writeObservableToFileHistSelect(TString); -const TString colorwhite = "\033[1;29m"; -const TString colorred = "\033[1;31m"; -const TString coloryellow = "\033[1;33m"; -const TString colorreset = "\033[0m"; -const TString endlr = "\033[0m\n"; -TString ownpath = ""; +TString ownpath = ""; void setCustomPath(TString setOwnpath) { @@ -77,8 +72,8 @@ Bool_t writeObservableToFile(TString histogramtype) { for(Int_t runi=0;runihistogramthresholdCalibrated.calibrated) - runs[runi]->plothistogramstructpointer = &runs[runi]->histogramthresholdCalibrated; + if (runs[runi]->histogramthreshold->iscalibrated) + runs[runi]->plothistogramclasspointer = runs[runi]->histogramthreshold->calibrated; else return 1; } @@ -86,18 +81,18 @@ Bool_t writeObservableToFile(TString histogramtype) else { for(Int_t runi=0;runiplothistogramstructpointer = &runs[runi]->histogramthreshold; + runs[runi]->plothistogramclasspointer = runs[runi]->histogramthreshold; } } else if (histogramtype.Contains("calibrated") || histogramtype.Contains("electron")) { for(Int_t runi=0;runihistogramCalibrated.calibrated) - runs[runi]->plothistogramstructpointer = &runs[runi]->histogramCalibrated; + if (runs[runi]->histogram->iscalibrated) + runs[runi]->plothistogramclasspointer = runs[runi]->histogram->calibrated; else return 1; } } else { for(Int_t runi=0;runiplothistogramstructpointer = &runs[runi]->histogram; } } + runs[runi]->plothistogramclasspointer = runs[runi]->histogram; } } writeObservableToFileHistSelect(histogramtype); return 0; } @@ -105,16 +100,16 @@ void writeObservableToFileHistSelect(TString histogramtype) { if (histogramtype.Contains("Seed") || histogramtype.Contains("seed")) { for(Int_t runi=0;runiplothistogrampointer = &runs[runi]->plothistogramstructpointer->Seed; } } + runs[runi]->plothistogrampointer = runs[runi]->plothistogramclasspointer->Seed; } } else if (histogramtype.Contains("Sum") || histogramtype.Contains("sum")) { for(Int_t runi=0;runiplothistogrampointer = &runs[runi]->plothistogramstructpointer->Sum; } } + runs[runi]->plothistogrampointer = runs[runi]->plothistogramclasspointer->Sum; } } else if (histogramtype.Contains("Veto") || histogramtype.Contains("veto")) { for(Int_t runi=0;runiplothistogrampointer = &runs[runi]->plothistogramstructpointer->Veto; } } + runs[runi]->plothistogrampointer = runs[runi]->plothistogramclasspointer->Veto; } } else { for(Int_t runi=0;runiplothistogrampointer = &runs[runi]->plothistogramstructpointer->Seed; } } + runs[runi]->plothistogrampointer = runs[runi]->plothistogramclasspointer->Seed; } } writeObservableToFile(); } @@ -126,7 +121,7 @@ Bool_t writeObservableToFile() { headerInfo+=runs[runi]->runcode+"\t\t\t"; } - TH1F* plothistogrampointer = *runs[0]->plothistogrampointer; + TH1F* plothistogrampointer = runs[0]->plothistogrampointer; TString runnumberListe=""; TString header=""; for(Int_t runi=0;runierror) { - plothistogrampointer = *runs[runi]->plothistogrampointer; + plothistogrampointer = runs[runi]->plothistogrampointer; runnumberListe+=Form("%d_",runs[runi]->labbook.runnumber); header+=Form("%d\t%s\t\t\t", runs[runi]->labbook.runnumber, plothistogrampointer->GetName()); } @@ -153,7 +148,7 @@ Bool_t writeObservableToFile() { if (!runs[runi]->error) { - plothistogrampointer = *runs[runi]->plothistogrampointer; + plothistogrampointer = runs[runi]->plothistogrampointer; TString outline; Double_t binContentNorm=plothistogrampointer->GetBinContent(bini)/runs[runi]->labbook.frames_foundDB*10000000; Double_t binSumme=0.0; @@ -181,27 +176,27 @@ Bool_t plotAllRuns(TString histogramtype) if (histogramtype.Contains("calibrated")) { for(Int_t runi=0;runihistogramthresholdCalibrated.calibrated) - runs[runi]->plothistogramstructpointer = &runs[runi]->histogramthresholdCalibrated; + if (runs[runi]->histogramthreshold->iscalibrated) + runs[runi]->plothistogramclasspointer = runs[runi]->histogramthreshold->calibrated; else return 1; } else { for(Int_t runi=0;runiplothistogramstructpointer = &runs[runi]->histogramthreshold; + runs[runi]->plothistogramclasspointer = runs[runi]->histogramthreshold; } } else if (histogramtype.Contains("calibrated") || histogramtype.Contains("electron")) { for(Int_t runi=0;runihistogramCalibrated.calibrated) - runs[runi]->plothistogramstructpointer = &runs[runi]->histogramCalibrated; + if (runs[runi]->histogram->iscalibrated) + runs[runi]->plothistogramclasspointer = runs[runi]->histogram->calibrated; else return 1; } } else { for(Int_t runi=0;runiplothistogramstructpointer = &runs[runi]->histogram; } } + runs[runi]->plothistogramclasspointer = runs[runi]->histogram; } } plotAllRuns(); return 0; } @@ -234,24 +229,24 @@ Bool_t plotAllRuns() if (!runs[runi]->error) { canvas->cd(1); - runs[runi]->plothistogramstructpointer->Seed->Draw("SAME"); - lastbin = runs[runi]->plothistogramstructpointer->Seed->GetBinCenter(runs[runi]->plothistogramstructpointer->Seed->FindLastBinAbove(2,1)); - runs[runi]->plothistogramstructpointer->Seed->SetAxisRange(0,lastbin*1.1,"X"); + runs[runi]->plothistogramclasspointer->Seed->Draw("SAME"); + lastbin = runs[runi]->plothistogramclasspointer->Seed->GetBinCenter(runs[runi]->plothistogramclasspointer->Seed->FindLastBinAbove(2,1)); + runs[runi]->plothistogramclasspointer->Seed->SetAxisRange(0,lastbin*1.1,"X"); gPad->SetLogy(1); - legendEntry = Form("%s", runs[runi]->plothistogramstructpointer->Seed->GetTitle()); - leg1->AddEntry(runs[runi]->plothistogramstructpointer->Veto, legendEntry, "l"); + legendEntry = Form("%s", runs[runi]->plothistogramclasspointer->Seed->GetTitle()); + leg1->AddEntry(runs[runi]->plothistogramclasspointer->Veto, legendEntry, "l"); leg1->Draw("SAME"); canvas->cd(2); - runs[runi]->plothistogramstructpointer->Sum->Draw("SAME"); - lastbin = runs[runi]->plothistogramstructpointer->Sum->GetBinCenter(runs[runi]->plothistogramstructpointer->Sum->FindLastBinAbove(2,1)); - runs[runi]->plothistogramstructpointer->Sum->SetAxisRange(0,lastbin*1.1,"X"); + runs[runi]->plothistogramclasspointer->Sum->Draw("SAME"); + lastbin = runs[runi]->plothistogramclasspointer->Sum->GetBinCenter(runs[runi]->plothistogramclasspointer->Sum->FindLastBinAbove(2,1)); + runs[runi]->plothistogramclasspointer->Sum->SetAxisRange(0,lastbin*1.1,"X"); canvas->cd(3); - runs[runi]->plothistogramstructpointer->Veto->Draw("SAME"); - runs[runi]->plothistogramstructpointer->Veto->SetAxisRange(runs[runi]->plothistogramstructpointer->posVeto*0.7,runs[runi]->plothistogramstructpointer->posVeto*1.4,"X"); + runs[runi]->plothistogramclasspointer->Veto->Draw("SAME"); + runs[runi]->plothistogramclasspointer->Veto->SetAxisRange(runs[runi]->plothistogramclasspointer->posVeto*0.7,runs[runi]->plothistogramclasspointer->posVeto*1.4,"X"); canvas->cd(4); - runs[runi]->plothistogramstructpointer->Noise->Draw("SAME"); - legendEntry = Form("%s, Noise: %.2f", runs[runi]->labbook.matrix.Data(), runs[runi]->plothistogramstructpointer->avgNoise); - leg2->AddEntry(runs[runi]->plothistogramstructpointer->Veto, legendEntry, "l"); + runs[runi]->plothistogramclasspointer->Noise->Draw("SAME"); + legendEntry = Form("%s, Noise: %.2f", runs[runi]->labbook.matrix.Data(), runs[runi]->plothistogramclasspointer->avgNoise); + leg2->AddEntry(runs[runi]->plothistogramclasspointer->Veto, legendEntry, "l"); leg2->Draw("SAME"); } } diff --git a/MABS_run_analyzer/HistogramType.c b/MABS_run_analyzer/HistogramType.c new file mode 100644 index 0000000..a2d2dcf --- /dev/null +++ b/MABS_run_analyzer/HistogramType.c @@ -0,0 +1,419 @@ +/** + * @file HistogramType.c + * @brief A class to store the histograms used in #Run.c + * + */ + +#ifndef __HISTOGRAMTYPE__C +#define __HISTOGRAMTYPE__C + +#include "HistogramType.h" +using namespace std; + + +//***************** +// CONSTRUCTOR AND DECONSTRUCTOR +//***************** + +HistogramType::~HistogramType( void) { + +} + +HistogramType::HistogramType(TString suffix, systemparam* gotsystempar, sensorinfostruct* gotsensorinfo, Bool_t threshold, Int_t gotcolor, Int_t gotstyle ) { + histogramdescription = suffix; + cursystempar = gotsystempar; + cursensorinfo = gotsensorinfo; + setThresholdType(threshold); + initHistograms(gotcolor, gotstyle); +}; + + +//***************** +// GETTER AND SETTER +//***************** + +void HistogramType::setThresholdType(Bool_t threshold) { + if (threshold) { + isthresholdclustertype = 1; } + else { + isthresholdclustertype = 0; } +} + +//***************** +// INIT FUNCTIONS +//***************** + +void HistogramType::initHistograms(Int_t gotcolor, Int_t gotstyle) { + color = gotcolor; + style = gotstyle; + initHistogram(Seed, "Seed" + histogramdescription, color, style); + initHistogram(Sum, "Sum" + histogramdescription, color, style); + initHistogram(Veto, "Veto" + histogramdescription, color, style); + initHistogram(Noise, "Noise" + histogramdescription, color, style); + + Noise->SetBins(cursystempar->nbinsnoise, 0, cursystempar->maxbinnoise); +} + +void HistogramType::initHistogram(TH1F* &histogrampointer, TString prefix, Int_t color, Int_t style) { + histogrampointer=new TH1F(prefix.Data(), prefix.Data(), cursystempar->nbins, 0, cursystempar->maxbin); + histogrampointer->SetLineStyle(style); + histogrampointer->SetLineColor(color); + histogrampointer->SetStats(kTRUE); + histogrampointer->SetStats(111111111); + histogrampointer->SetLineWidth(3); + histogrampointer->GetXaxis()->SetTitle("Q_coll [ADU]"); + histogrampointer->GetYaxis()->SetTitle(Form("Entries [1/%.1f ADU]",histogrampointer->GetBinWidth(1))); + histogrampointer->GetXaxis()->CenterTitle(); + histogrampointer->GetYaxis()->CenterTitle(); +} + +void HistogramType::initCluster(TString suffix, Float_t xcoord_min_step, Float_t xcoord_abs_min, Float_t xcoord_abs_max, Float_t ycoord_min_step, Float_t ycoord_abs_min, Float_t ycoord_abs_max) { + TString prefix = "Cluster" + suffix; + Int_t sizearrrotx= (int)((xcoord_abs_max-xcoord_abs_min)/xcoord_min_step+0.5); + Int_t sizearrroty = (int)((ycoord_abs_max-ycoord_abs_min)/ycoord_min_step+0.5); + // cout << "xcoord_abs_min: " << xcoord_abs_min << ", xcoord_abs_max: " << xcoord_abs_max << ", xcoord_min_step: " << xcoord_min_step << endl; + // cout << "ycoord_abs_min: " << ycoord_abs_min << ", ycoord_abs_max: " << ycoord_abs_max << ", ycoord_min_step: " << ycoord_min_step << endl; + histAvgCluster = new TH2F(prefix.Data(), "Avg. cluster distribution; x; y",sizearrrotx, xcoord_abs_min,xcoord_abs_max,sizearrroty,ycoord_abs_min,ycoord_abs_max); + // you can increase the resolution, especially usefull if you observe an rotated (originally staggered) pixel cluster + // histogramstructpointer->histAvgCluster = new TH2F(prefix.Data(), Form("Avg. cluster distribution; x; y %s", humanreadablestr.Data()),2*sizearrrotx, xcoord_abs_min,xcoord_abs_max,sizearrroty*2,ycoord_abs_min,ycoord_abs_max); +} + +//***************** +// METHODS APPLYABLE TO HISTOGRAMS +//***************** + +Bool_t HistogramType::calibrateHistograms( Float_t gotgain ) { + if ( gotgain <= 0 ) + { + cout << coloryellow << "Cannot rescale run from [ADU] to [e] units, no gain provided." << endlr; + return 1; + } + gain = gotgain; + + calibrated = new HistogramType(histogramdescription+" calibrated", cursystempar, cursensorinfo, isthresholdclustertype, color, style); + if (Seed != 0) calibrateHistogram(calibrated->Seed, Seed); + if (Sum != 0) calibrateHistogram(calibrated->Sum, Sum); + if (Veto != 0) calibrateHistogram(calibrated->Veto, Veto); + if (Noise != 0) calibrateHistogram(calibrated->Noise, Noise); + if (histAvgCluster != 0) calibrate2DHistogramCounts(calibrated->histAvgCluster, histAvgCluster); + + calibrated->posSeed = posSeed * gain; + calibrated->posSum = posSum * gain; + calibrated->posVeto = posVeto * gain; + calibrated->avgNoise = avgNoise * gain; + calibrated->avgNoisePlus = avgNoisePlus * gain; + calibrated->avgNoiseMinus = avgNoiseMinus * gain; + calibrated->sr90IntegralVal = sr90IntegralVal * gain; + + calibrated->iscalibrated = true; + + return 1; +} + +void HistogramType::calibrateHistogram(TH1F* &histogrampointernew, TH1F* &histogrampointerold) { + histogrampointernew = (TH1F*)histogrampointerold->Clone(); + histogrampointernew->SetName(Form("%sC", histogrampointerold->GetName())); + histogrampointernew->SetTitle(Form("%sC", histogrampointerold->GetTitle())); + histogrampointernew->GetXaxis()->SetTitle("Q_coll [e]"); + int nbins = histogrampointernew->GetXaxis()->GetNbins(); + double new_bins[nbins+1]; + for(int i=0; i <= nbins; i++){ + new_bins[i] = histogrampointernew->GetBinCenter(i)*gain; + } + histogrampointernew->SetBins(nbins, new_bins); + histogrampointernew->GetYaxis()->SetTitle(Form("Entries [1/%.1f e]",histogrampointernew->GetBinWidth(1))); +} + +void HistogramType::calibrate2DHistogramCounts(TH2F* &histogrampointernew, TH2F* &histogrampointerold ) { + histogrampointernew = (TH2F*)histogrampointerold->Clone(); + histogrampointernew->SetName(Form("%sC", histogrampointerold->GetName())); + histogrampointernew->SetTitle(Form("%sC", histogrampointerold->GetTitle())); + histogrampointernew->GetZaxis()->SetTitle("Q_coll [e]"); + int nbinsx = histogrampointernew->GetXaxis()->GetNbins(); + int nbinsy = histogrampointernew->GetYaxis()->GetNbins(); + + for(int y=0; y <= nbinsy; y++){ + for(int x=0; x <= nbinsx; x++){ + histogrampointernew->SetBinContent(x,y,histogrampointerold->GetBinContent(x,y)*gain); + } + } +} + +Float_t HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType, Bool_t verbose, Double_t* parameters) { + Float_t posMax = 0; + Float_t posMax2 = 0; + Float_t posMaxValHist = histogrampointer->GetXaxis()->GetXmax(); + // posMaxValHist/10 for USB system, the value is 90 + Float_t noiseborder = 90; + if (noisethresholdborder>0) + noiseborder = noisethresholdborder; + else if (noisethresholdborder > 0) + noiseborder = noisethresholdborder; + + if (fitFuncType.Contains("gaus")) + { + 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); + + if (TString(histogrampointer->GetName()).Contains("Veto")) + { + Float_t peak1 = histogrampointer->GetMaximumBin(); + Float_t peak2; + if (histogrampointer->GetBinCenter(peak1)<3.3*noiseborder) // vermutlich ist veto peak nicht am höchsten + { + Float_t avg = 0; + for(Int_t bin=histogrampointer->GetXaxis()->FindBin(noiseborder);binFindLastBinAbove(0);bin++) + avg += histogrampointer->GetBinContent(bin); + avg /= histogrampointer->FindLastBinAbove(0) - histogrampointer->GetXaxis()->FindBin(noiseborder); + peak2 = histogrampointer->FindLastBinAbove(avg/3); + } + 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->Fit(fitFunc, "N,Q,W", "", histogrampointer->GetXaxis()->GetBinCenter(min), posMaxValHist); + posMax = histogrampointer->GetXaxis()->GetBinCenter(histogrampointer->GetMaximumBin()); // Methode 1 + 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", "", noiseborder, posMaxValHist); + posMax = fitFunc->GetMaximumX(); // Methode 2 + } + Float_t sigma = fitFunc->GetParameter(2); + posMax2 = fitFunc->GetMaximumX(); // Methode 2 + Float_t peakposdifperc = abs(posMax-posMax2)/min(posMax,posMax2); + if (sigma > 260 || peakposdifperc>0.3) + { + if (verbose) + { + if (sigma > 260) + { + cout << "Sigma suspiciously height when fitting " << histogrampointer->GetName() << " spectrum: " << sigma << endl; + cout << colorred << " Could not find " << histogrampointer->GetName() << " peak" << endlr; + } + else if (peakposdifperc>0.3) + { + cout << "Maximum peak position and fit gaussian peak position doesn't match in " << histogrampointer->GetName() << " spectrum. Difference: " << peakposdifperc*100 <<" % " << endl; + cout << colorred << " Could not find " << histogrampointer->GetName() << " peak" << endlr; + } + } + return 0; + } + else if (sigma > 40 || peakposdifperc>0.1) + { + if (verbose) + { + if (sigma > 260) + { + cout << "Sigma or suspiciously height when fitting " << histogrampointer->GetName() << " spectrum: " << sigma << endl; + cout << coloryellow << " Please check " << histogrampointer->GetName() << " peak position: " << colorreset << colorwhite << posMax << endlr; + } + else if (peakposdifperc>0.1) + { + cout << "Maximum peak position and fit gaussian peak position in " << histogrampointer->GetName() << " have difference of " << peakposdifperc*100 <<" % " << endl; + cout << coloryellow << " Please check " << histogrampointer->GetName() << " peak position: " << colorreset << colorwhite << posMax << endlr; + } + } + } + // fitFunc->DrawCopy("same"); + 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); + // leg->SetHeader();//"Legend Title"); + leg->SetFillStyle(0); + leg->AddEntry((TObject*) 0, legendEntry, ""); + leg->SetTextSize(0.05); + // leg->Draw(); + + } + else if (fitFuncType=="landau") + { + 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); + + Float_t fitMax1 = 1000; + Float_t fitMax2 = 1000; + Float_t fitMax3 = 1000; + Float_t minFitMax = 1000; + Float_t maxFitMax = 1000; + histogrampointer->Fit(fitFunc, "N,Q,W", "", noiseborder, posMaxValHist); + fitMax1 = fitFunc->GetMaximumX(); + // fitFunc->DrawCopy("same"); + histogrampointer->Fit(fitFunc, "N,Q,W", "", 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); + // histogrampointer->Fit(fitFunc, "N,Q,W", "", fitMax1, histogrampointer->GetBinCenter(bini)); + fitMax3 = fitFunc->GetMaximumX(); + fitFunc->SetLineColor(kGreen); + // fitFunc->DrawCopy("same"); + fitFunc->SetLineStyle(1); // normal for the following fits + + // Sort the three fits and save error estimation + minFitMax = TMath::Min(TMath::Min(fitMax1,fitMax2),fitMax3); + maxFitMax = TMath::Max(TMath::Max(fitMax1,fitMax2),fitMax3); + posMax = fitMax1 + fitMax2 + fitMax3 - minFitMax - maxFitMax; + + //fitLandauErrorLeft.push_back(posMax - minFitMax); + //fitLandauErrorRight.push_back(maxFitMax - posMax); + } else if (fitFuncType=="lorentz") + { + // create a TF1 with the range from 0 to 3 and 6 parameters + TF1 *fitFcn = new TF1("fitFcn",lorentzianPeak,0,160,4); + fitFcn->SetNpx(1000); + fitFcn->SetLineWidth(4); + fitFcn->SetLineColor(color); + // set start values for some parameters + fitFcn->SetParName(0,"Area"); + fitFcn->SetParameter(0,21655); // Area + fitFcn->SetParName(1,"Width"); + fitFcn->SetParameter(1,34); // width + fitFcn->SetParName(2,"Peak pos."); + fitFcn->SetParameter(2,80); // peak position x + fitFcn->SetParName(3,"Y shift"); + fitFcn->SetParameter(3,-22); // y-shift + gStyle->SetOptFit(0111); + + // fitFcn->FixParameter(0,21655.90045); + // fitFcn->FixParameter(1,34.31885101); + fitFcn->FixParameter(2,histogrampointer->GetXaxis()->GetBinCenter((int)(histogrampointer->GetXaxis()->GetNbins()/2.0+0.5))); + // fitFcn->FixParameter(3,-22.43016284); + + //histogrampointer->Fit("fitFcn","V+","ep"); + histogrampointer->Fit("fitFcn","Q,W","ep"); + for (Int_t pari=0; pari<4; pari++) + parameters[pari]=fitFcn->GetParameter(pari); + if (verbose) + { + Printf("width: %.6f ",fitFcn->GetParameter(1)); + cout << "width: " << fitFcn->GetParameter(1) << endl; + cout << "peak: " << fitFcn->GetParameter(2) << endl; + cout << "y-shift:" << fitFcn->GetParameter(3) << endl; + } + } + + return posMax; +} + +Bool_t HistogramType::calculteCCE() { + if (posSeed > 0 && posVeto > 0) + CCE_in_Perc_1 = posSeed / posVeto * 100.0; + if (posSum > 0 && posVeto > 0) + CCE_in_Perc_25 = posSum / posVeto * 100.0; + if (CCE_in_Perc_1 > 0 || CCE_in_Perc_25 > 0) + return 0; + return 1; +} + +Bool_t HistogramType::integrateSr90Spectra(TH1F* histogrampointer, Int_t frames_found, Float_t thresholdborder, Bool_t verbose) { + Float_t posMaxValHist = histogrampointer->GetXaxis()->GetXmax(); + + TH1F* smoothedcurce = (TH1F*)histogrampointer->Clone(); + smoothedcurce->Smooth(4); + + // Int_t random = random1->Rndm()*100000; + // TString canvastitle = Form("%s %s sdfasdf", histogrampointer->GetName(), runcode.Data()); + // TString canvasname = Form("%s %s %d dfgsdfgsdfg", histogrampointer->GetName(), runcode.Data(), random); + // TCanvas* canvas = new TCanvas(canvasname, canvastitle, 900, 700); + // smoothedcurce->Draw(); + + Int_t thresholdbincurcandidate = 0; + if (thresholdborder < 0) + { + Int_t rising = 0; + Int_t bini =smoothedcurce->GetMaximumBin(); + Float_t curval = smoothedcurce->GetXaxis()->GetBinCenter(bini); + // cout << "smoothedcurce->GetXaxis()->GetBinCenter(" << bini << "): " << curval << endl; + + do { + curval=smoothedcurce->GetBinContent(bini++); + if (curval*0.95 <= smoothedcurce->GetBinContent(bini)) + { + rising++; + // cout << "rising at " << smoothedcurce->GetXaxis()->GetBinCenter(bini) << " as " << curval << " < " << smoothedcurce->GetBinContent(bini) << endl; // debug + } + else + { + rising = 0; + thresholdbincurcandidate = bini; + } + } while (rising < 3 && biniGetNbinsX()); + } + else + { + thresholdbincurcandidate = histogrampointer->GetXaxis()->FindBin(thresholdborder); + } + + + // histogrampointer->GetXaxis()->SetRange(histogrampointer->GetXaxis()->FindBin(0),histogrampointer->GetXaxis()->FindBin(posMaxValHist)); + // Float_t posNoiseMax= histogrampointer->GetMaximum(); + // histogrampointer->GetXaxis()->SetRange(posNoiseMax,histogrampointer->GetXaxis()->FindBin(posMaxValHist)); + // Float_t posChargeMax= histogrampointer->GetMaximum(); + + noisethresholdborder = histogrampointer->GetXaxis()->GetBinUpEdge(thresholdbincurcandidate); + + if (verbose) { + cout << " Integrating " << histogrampointer->GetTitle() << endlr; + cout << " Noise threshold at " << noisethresholdborder; + TString histtitle=histogrampointer->GetXaxis()->GetTitle(); + if (histtitle.Contains("[e]")) + { + cout << " e" << endl; + if (noisethresholdborder > 300) { + cout << coloryellow << " This noise threshold seems too high, please check manually the " << histogrampointer->GetName() << " spectrum." << endlr; + } + } + else + { + cout << " ADU" << endl; + if (noisethresholdborder > 75) { + cout << coloryellow << " This noise threshold seems too high, please check manually the " << histogrampointer->GetName() << " spectrum." << endlr; + } + } + } + sr90IntegralVal = histogrampointer->IntegralAndError(thresholdbincurcandidate,histogrampointer->GetNbinsX(), sr90IntegralErr); + // cout << "Integrate from bin " << thresholdbincurcandidate << " to " << histogrampointer->GetNbinsX() << endl; + sr90IntegralErr /= sr90IntegralVal/100; +// cout << " Unscaled integral is " << Form("%e",sr90IntegralVal) << ", error " << sr90IntegralErr << "%" << endl; +// cout << " "; + if (frames_found>0) + { + sr90IntegralVal/=frames_found; + if (verbose) + cout << " Scaled"; + } + if (verbose) + cout << " Integral is " << Form("%e",sr90IntegralVal) << ", error " << sr90IntegralErr << "%" << endl; + + return 0; +} + + +//***************** +// OTHER HELPFUL FUNCTIONS +//***************** + +// Lorenzian Peak function +Double_t HistogramType::lorentzianPeak(Double_t *x, Double_t *par) { + return par[3]+(0.5*par[0]*par[1]/TMath::Pi()) / + TMath::Max( 1.e-10,(x[0]-par[2])*(x[0]-par[2]) + + .25*par[1]*par[1]); +} + +#endif + diff --git a/MABS_run_analyzer/HistogramType.h b/MABS_run_analyzer/HistogramType.h new file mode 100644 index 0000000..f8deeab --- /dev/null +++ b/MABS_run_analyzer/HistogramType.h @@ -0,0 +1,202 @@ +#ifndef __HISTOGRAMTYPE__H +#define __HISTOGRAMTYPE__H + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "help.h" + +/** + * @file HistogramType.h + * @brief Has all the raw and processed information about a specific run + * + * This class makes SQL requests to the server, gets database results to a given runnumber, + * stores information about a specific run, writes the analyzed data and results back + * to the databases and stores the analyzed data as histograms. + * + */ +class HistogramType { + +private: + + //***************** + // FITTING RESULTS OF TH1F HISTOGRAMS + //***************** + + + //***************** + // GENERAL HISTOGRAM PROPERTIES + //***************** + + Int_t color; + Int_t style; + + + //***************** + // GENERAL SYSTEM PROPERTIES + //***************** + /// system information to use in analysis, is the system read out by USB or PXI? Number of bins and limits differ + systemparam* cursystempar; + sensorinfostruct* cursensorinfo; + + //***************** + // PRIVATE METHODS APPLYABLE TO HISTOGRAMS + //***************** + /** + * @brief rescales one specific histogram from ADU to electrons */ + void calibrateHistogram(TH1F* &histogrampointernew, TH1F* &histogrampointerold); + + //***************** + // OTHER + //***************** + + + //***************** + // PRIVATE METHODS + //***************** + + /** + * @brief Lorenzian Peak function */ + static Double_t lorentzianPeak(Double_t *x, Double_t *par); + + +public: + + //***************** + // CONSTRUCTOR AND DECONSTRUCTOR + //***************** + + ~HistogramType(void); + /** @brief constructor */ + HistogramType(TString suffix, systemparam* gotsystempar, sensorinfostruct* gotsensorinfo, Bool_t threshold, Int_t gotcolor, Int_t gotstyle ); + + //***************** + // TH HISTOGRAMS + //***************** + + /// Seed spectrum, only the seed pixel is considered when this histogram is build + TH1F* Seed = 0; + /// Sum spectrum, the charge of whole cluster is summed up in binned into this TH1F histogram + TH1F* Sum = 0; + /// Veto spectrum, used to find better calibration peak, only entries where Sum over not seed pixels is below @c Run::vetothreshold are binned into this histogram + TH1F* Veto = 0; + /// Noise histogram + TH1F* Noise = 0; + /// histogram with cluster data, rotated if staggered + TH2F* histAvgCluster = 0; + + //***************** + // FITTING RESULTS OF TH1F HISTOGRAMS + //***************** + + /// fitted position of the most probable value of the seed spectrum + Float_t posSeed=0; + /// fitted position of the most probable value in the over cluster summed spectrum + Float_t posSum=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; + + /// Average Noise value + Float_t avgNoise = 0; + /// Positive Noise sigma + Float_t avgNoisePlus = 0; + /// Negative Noise sigma + Float_t avgNoiseMinus = 0; + + /// Integral value, after integrating from #noisethresholdborder to maxbin. + Double_t sr90IntegralVal = -1; + Double_t sr90IntegralErr = -1; + + /// threshold for integrating the Sr90 spectrum, is set dynamically in @c integrateSr90Spectra() + Float_t noisethresholdborder = -1; + + //***************** + // GETTER AND SETTER + //***************** + void setThresholdType(Bool_t threshold); + + //***************** + // INIT FUNCTIONS + //***************** + /// + /** @brief initializes all histograms, set binwidth, bin amount and names + * * + * @param histogramstruct Pointer to a structure of type #Run::histogramstruct + * @param suffix A string which will be appended to the generated title of each histograms + */ + void initHistograms(Int_t gotcolor, Int_t gotstyle); + /// init one specific histogram + void initHistogram(TH1F* &histogrampointer, TString prefix, Int_t color, Int_t style); + /** @brief initializes the TH2F cluster for the binning for a specific structure of type #Run::histogramstruct one points to*/ + void initCluster(TString suffix, Float_t xcoord_min_step, Float_t xcoord_abs_min, Float_t xcoord_abs_max, Float_t ycoord_min_step, Float_t ycoord_abs_min, Float_t ycoord_abs_max); + + + //***************** + // GENERAL HISTOGRAM PROPERTIES + //***************** + /// type in here what the histogram is intended for or how it is calculated, will be added to filenames + TString histogramdescription = ""; + /// The gain used to rescale the histograms + Float_t gain = 0; + /// set to true, if bins are in electrons, otherwise in ADU + Bool_t iscalibrated = false; + /// set to true, if threshold clusters are used + Bool_t isthresholdclustertype = false; + /// number of hits/clusters used to generate all distributions + Double_t numberofhits = 0; + /// Charge collection efficciency of the cluster in percent + Float_t CCE_in_Perc_25=0; + /// Charge collection efficciency of the seed pixel in percent + Float_t CCE_in_Perc_1=0; + + //***************** + // METHODS APPLYABLE TO HISTOGRAMS + //***************** + + /** + * @brief rescales all histograms from ADU to electrons */ + Bool_t calibrateHistograms( Float_t gotgain ); + /** + * @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 + * + * @param histogrampointer histogram pointer to calculate fit to + * @return peak position of the fit + * + */ + Float_t FitPerform(TH1F* histogrampointer, TString fitFuncType, Bool_t verbose = 0, Double_t* parameters = 0); + /** @brief calculates Charge Collection Efficiency if Fe55 run */ + Bool_t calculteCCE(); + /** + * @brief find the border between the noise and the signal in Sr90 runs + * + * writes the found threshold into the private variable @c Run::posNoiseThreshold + * + * @param histogrampointer pointer to the histogram structure, threshold will be searched in seed spectra + * @param thresholdborder value from which the integral will be calculated, if not set, the algorithms tries to find the best value + * @return true if succesfull + * + */ + Bool_t integrateSr90Spectra(TH1F* histogrampointer, Int_t frames_found=-1, Float_t thresholdborder=-1, Bool_t verbose=true); + + + //***************** + // OTHER + //***************** + + /// pointer to calibrated HistogramType instance of this class + HistogramType* calibrated = 0; + +}; + +#endif \ No newline at end of file diff --git a/MABS_run_analyzer/MAPS.c b/MABS_run_analyzer/MAPS.c index ca788a0..f9c1582 100644 --- a/MABS_run_analyzer/MAPS.c +++ b/MABS_run_analyzer/MAPS.c @@ -78,6 +78,11 @@ Bool_t MAPS::initNewRootFile() { Bool_t MAPS::initOldRootFile() { fSave = false; + if ( !checkConf(fPixelsData) ) { // TODO FileEvNbInConfig + defaultConf(); + } + int MaxFiles = TMath::Ceil((Float_t) FileTotalEvNbInConfig/FileEvNbInConfig); + checkDataFiles(MaxFiles); fOutputFile = new TFile(fRootFile,"READ"); if (fOutputFile->IsZombie()) { @@ -202,14 +207,15 @@ Bool_t MAPS::initMapsRun( ) { fOutDir = run->storepathOutLinux; // default ouput directory is input directory fRunNumber = run->labbook.runnumber; - fRows = run->sensorinfocur.rows; - fColumns = run->sensorinfocur.columns; + fRows = run->cursensorinfo.rows; + fColumns = run->cursensorinfo.columns; fPixels = fRows*fColumns; fMatrix = run->labbook.matrix; fOrderCode = run->labbook.chipGen; fSystem = run->labbook.system; fConfigFile = fInDir+Form("/RUN_%i_i.rz",fRunNumber); fRootFile = fOutDir+Form("/RUN_%i_i.root",fRunNumber); + fEventsSum = run->labbook.frames_foundDB; cout<<"================================================================="<SetCanvasColor(0); gStyle->SetFrameFillColor(10); gStyle->SetOptStat(0); - + fEventsSum = 1000000000; // remove me TODO if( Frames > fEventsSum ) { Frames = fEventsSum; printf("Changed 'Number Frames' to: %u\n", Frames ); @@ -2062,21 +2068,30 @@ void MAPS::plotPixSignal(UInt_t Start, UInt_t Frames, Int_t Pixel) { printf("Changed 'First Event' to: %u\n", Start ); } - TH2F *h1 = new TH2F("Frame 0 vs T" , "Frame 0 vs T" , Frames, Start, Start+Frames, 2*16384, -16384, 16384); - TH2F *h2 = new TH2F("Frame 1 vs T" , "Frame 1 vs T" , Frames, Start, Start+Frames, 2*16384, -16384, 16384); - TH2F *h3 = new TH2F("CDS vs T" , "CDS vs T" , Frames, Start, Start+Frames, 2*16384, -16384, 16384); +// TH2F *h1 = new TH2F("Frame 0 vs T" , "Frame 0 vs T" , Frames, Start, Start+Frames, 2*16384, -16384, 16384); +// TH2F *h2 = new TH2F("Frame 1 vs T" , "Frame 1 vs T" , Frames, Start, Start+Frames, 2*16384, -16384, 16384); +// TH2F *h3 = new TH2F("CDS vs T" , "CDS vs T" , Frames, Start, Start+Frames, 2*16384, -16384, 16384); + TH1F *h1 = new TH1F("Frame 0 vs T" , "Frame 0 vs T" , Frames, Start, Start+Frames); + TH1F *h2 = new TH1F("Frame 1 vs T" , "Frame 1 vs T" , Frames, Start, Start+Frames); + TH1F *h3 = new TH1F("CDS vs T" , "CDS vs T" , Frames, Start, Start+Frames); - TH1F *h4 = new TH1F("Frames 0 phisto" , "Frames 0 phisto" , 2*16384, -16384, 16384); - TH1F *h5 = new TH1F("Frames 1 phisto" , "Frames 1 phisto" , 2*16384, -16384, 16384); - TH1F *h6 = new TH1F("CDS phisto" , "CDS phisto" , 2*16384, -16384, 16384); + TH1F *h4 = new TH1F("Frames 0 phisto" , "Frames 0 phisto" , 2*16384+1, -16384, 16384); + TH1F *h5 = new TH1F("Frames 1 phisto" , "Frames 1 phisto" , 2*16384+1, -16384, 16384); + TH1F *h6 = new TH1F("CDS phisto" , "CDS phisto" , 2*16384+1, -16384, 16384); +// h5->SetBit(TH1::kCanRebin); +// h6->SetBit(TH1::kCanRebin); +// h4->Rebin(k) for(UInt_t i=Start; iFill( i,fF0matrix[Pixel]) ; - h2->Fill( i,fF1matrix[Pixel]) ; - h3->Fill( i,fCdsmatrix[Pixel] ); +// h1->Fill( i,fF0matrix[Pixel]) ; +// h2->Fill( i,fF1matrix[Pixel]) ; +// h3->Fill( i,fCdsmatrix[Pixel] ); + h1->SetBinContent( i,fF0matrix[Pixel]) ; + h2->SetBinContent( i,fF1matrix[Pixel]) ; + h3->SetBinContent( i,fCdsmatrix[Pixel] ); h4->Fill( fF0matrix[Pixel] ); h5->Fill( fF1matrix[Pixel] ); @@ -2085,17 +2100,18 @@ void MAPS::plotPixSignal(UInt_t Start, UInt_t Frames, Int_t Pixel) { } cm4->cd(1); - h1->Draw("colz"); + // h1->Draw("colz"); + h1->Draw(""); h1->GetXaxis()->SetTitle("Frame#"); h1->GetYaxis()->SetTitle("Signal"); cm4->cd(2); - h2->Draw("colz"); + h2->Draw(""); h2->GetXaxis()->SetTitle("Frame#"); h2->GetYaxis()->SetTitle("Signal"); cm4->cd(3); - h3->Draw("colz"); + h3->Draw(""); h3->GetXaxis()->SetTitle("Frame#"); h3->GetYaxis()->SetTitle("Signal"); @@ -2119,7 +2135,7 @@ void MAPS::plotPixSignal(UInt_t Start, UInt_t Frames, Int_t Pixel) { // cm4->Close(); cout<<"\rPIXELSIGNALS plotted! "<initNewRootFile()) return 1; - frames_found = processed->GetNumberFrames(); /// progress meter, temporal variable ULong_t progress_tmp=-1; /// progress meter @@ -368,22 +377,26 @@ Bool_t Run::analyzeRun(Bool_t force) cout << colorwhite << "initOldRootFile():" << endlr; } if (processed->initOldRootFile()) return 1; - cout << colorwhite << "binNoise():" << endlr; - binNoise(); + frames_found = processed->GetNumberFrames(); +// cout << colorwhite << "plotPixSignal():"<< flush << endlr; +// processed->plotPixSignal(0,10000000,351); cout << colorwhite << "binSeedSumVeto():" << endlr; - binSeedSumVeto(); + binSeedSumVeto(); cout << colorwhite << "binCluster():" << endlr; binCluster(); - cout << colorwhite << "rescaleHistograms():" << endlr; - histogramCalibrated.calibrated = rescaleHistograms(); - histogramthresholdCalibrated.calibrated = histogramCalibrated.calibrated; - histogramfixedthresholdCalibrated.calibrated = histogramCalibrated.calibrated; - cout << colorwhite << "calculateCCE():" << endlr; - calculteCCE(); - if (labbook.source.Contains("Sr90")) { - cout << colorwhite << "integrateSr90Spectra():" << endlr; - integrateSr90Spectra(&histogramthresholdCalibrated, histogramthresholdCalibrated.Seed); + for (vector::iterator curHistogramClass = HistogramClassVector.begin(); curHistogramClass != HistogramClassVector.end(); curHistogramClass++) { + binNoise((*curHistogramClass)); + if (labbook.source.Contains("Fe")||labbook.source.Contains("Cd")) { + cout << colorwhite << "calculateCCE():" << endlr; + (*curHistogramClass)->calculteCCE(); + } + if (labbook.source.Contains("Sr90")) { + cout << colorwhite << "integrateSr90Spectra():" << endlr; + (*curHistogramClass)->integrateSr90Spectra((*curHistogramClass)->Seed, frames_found, -1, false); + } } + histogramfixedthreshold->integrateSr90Spectra(histogramfixedthreshold->Seed, frames_found, 0); + rescaleHistogramClasses(); cout << colorwhite << "updateDatabase():" << endlr; updateDatabase(); cout << colorwhite << "delete MAPS class:" << endlr; @@ -394,24 +407,12 @@ Bool_t Run::analyzeRun(Bool_t force) return 1; } -Bool_t Run::calculteCCE() -{ - if (labbook.source.Contains("Fe")||labbook.source.Contains("Cd")) - { - if (histogram.posSeed > 0 && histogram.posVeto > 0) - CCE_in_Perc_1 = histogram.posSeed / histogram.posVeto * 100.0; - if (histogram.posSum > 0 && histogram.posVeto > 0) - CCE_in_Perc_25 = histogram.posSum / histogram.posVeto * 100.0; - } - return 0; -} - -Bool_t Run::rescaleHistograms() +Bool_t Run::rescaleHistogramClasses() { float_t vetopeakposition = -1; - if ( histogram.posVeto > 0 ) + if ( histogram->posVeto > 0 ) { - vetopeakposition = histogram.posVeto; + vetopeakposition = histogram->posVeto; cout << colorwhite << "Use calibration obtained from this run, position veto: " << vetopeakposition << endlr; } else if ( labbook.posVetoDB > 0 ) @@ -431,88 +432,11 @@ Bool_t Run::rescaleHistograms() cout << coloryellow << "Cannot rescale run from [ADU] to [e] units, no calibration peak found." << endlr; return 0; } - gain = 1640.0/vetopeakposition; - - rescaleHistogram(histogramCalibrated.Seed, histogram.Seed); - rescaleHistogram(histogramCalibrated.Sum, histogram.Sum); - rescaleHistogram(histogramCalibrated.Veto, histogram.Veto); - rescaleHistogram(histogramCalibrated.Noise, histogram.Noise); - - if (histogram.histAvgCluster != nullptr) - { - rescale2DHistogramCounts(histogramCalibrated.histAvgCluster, histogram.histAvgCluster); + Float_t gain = 1640.0/vetopeakposition; + for (vector::iterator curHistogramClass = HistogramClassVector.begin(); curHistogramClass != HistogramClassVector.end(); curHistogramClass++) { + (*curHistogramClass)->calibrateHistograms(gain);; } - - histogramCalibrated.posSeed = histogram.posSeed * gain; - histogramCalibrated.posSum = histogram.posSum * gain; - histogramCalibrated.posVeto = histogram.posVeto * gain; - histogramCalibrated.avgNoise = histogram.avgNoise * gain; - histogramCalibrated.avgNoisePlus = histogram.avgNoisePlus * gain; - histogramCalibrated.avgNoiseMinus = histogram.avgNoiseMinus * gain; - - rescaleHistogram(histogramthresholdCalibrated.Seed, histogramthreshold.Seed); - rescaleHistogram(histogramthresholdCalibrated.Sum, histogramthreshold.Sum); - rescaleHistogram(histogramthresholdCalibrated.Veto, histogramthreshold.Veto); - rescaleHistogram(histogramthresholdCalibrated.Noise, histogramthreshold.Noise); - - histogramthresholdCalibrated.posSeed = histogramthreshold.posSeed * gain; - histogramthresholdCalibrated.posSum = histogramthreshold.posSum * gain; - histogramthresholdCalibrated.posVeto = histogramthreshold.posVeto * gain; - histogramthresholdCalibrated.avgNoise = histogramthreshold.avgNoise * gain; - histogramthresholdCalibrated.avgNoisePlus = histogramthreshold.avgNoisePlus * gain; - histogramthresholdCalibrated.avgNoiseMinus = histogramthreshold.avgNoiseMinus * gain; - - rescaleHistogram(histogramfixedthresholdCalibrated.Seed, histogramfixedthreshold.Seed); - rescaleHistogram(histogramfixedthresholdCalibrated.Sum, histogramfixedthreshold.Sum); - rescaleHistogram(histogramfixedthresholdCalibrated.Veto, histogramfixedthreshold.Veto); - rescaleHistogram(histogramfixedthresholdCalibrated.Noise, histogramfixedthreshold.Noise); - - histogramfixedthresholdCalibrated.posSeed = histogramfixedthreshold.posSeed * gain; - histogramfixedthresholdCalibrated.posSum = histogramfixedthreshold.posSum * gain; - histogramfixedthresholdCalibrated.posVeto = histogramfixedthreshold.posVeto * gain; - histogramfixedthresholdCalibrated.avgNoise = histogramfixedthreshold.avgNoise * gain; - histogramfixedthresholdCalibrated.avgNoisePlus = histogramfixedthreshold.avgNoisePlus * gain; - histogramfixedthresholdCalibrated.avgNoiseMinus = histogramfixedthreshold.avgNoiseMinus * gain; - - if (histogramthreshold.histAvgCluster != nullptr) - { - rescale2DHistogramCounts(histogramthresholdCalibrated.histAvgCluster, histogramthreshold.histAvgCluster); - } - - return 1; -} - -void Run::rescaleHistogram(TH1F* &histogrampointernew, TH1F* &histogrampointerold ) -{ - histogrampointernew = (TH1F*)histogrampointerold->Clone(); - histogrampointernew->SetName(Form("%sC", histogrampointerold->GetName())); - histogrampointernew->SetTitle(Form("%sC", histogrampointerold->GetTitle())); - histogrampointernew->GetXaxis()->SetTitle("Q_coll [e]"); - int nbins = histogrampointernew->GetXaxis()->GetNbins(); - double new_bins[nbins+1]; - for(int i=0; i <= nbins; i++){ - new_bins[i] = histogrampointernew->GetBinCenter(i)*gain; - } - histogrampointernew->SetBins(nbins, new_bins); - histogrampointernew->GetYaxis()->SetTitle(Form("Entries [1/%.1f e]",histogrampointernew->GetBinWidth(1))); -} - - -void Run::rescale2DHistogramCounts(TH2F* &histogrampointernew, TH2F* &histogrampointerold ) -{ - histogrampointernew = (TH2F*)histogrampointerold->Clone(); - histogrampointernew->SetName(Form("%sC", histogrampointerold->GetName())); - histogrampointernew->SetTitle(Form("%sC", histogrampointerold->GetTitle())); - histogrampointernew->GetZaxis()->SetTitle("Q_coll [e]"); - int nbinsx = histogrampointernew->GetXaxis()->GetNbins(); - int nbinsy = histogrampointernew->GetYaxis()->GetNbins(); - - double new_bins[nbinsx][nbinsy]; - for(int y=0; y < nbinsy; y++){ - for(int x=0; x < nbinsx; x++){ - histogrampointernew->SetBinContent(x,y,histogrampointerold->GetBinContent(x,y)*gain); - } - } + return 1; } Bool_t Run::analyzeFrame(Int_t frame) @@ -647,12 +571,12 @@ void Run::selectSubMatrixFSBB(TString matrixname) { if (matrixname.EqualTo("A0")) { - setSubMatrixBorders(0, sensorinfocur.columns/2-2, 0, sensorinfocur.rows, false); + setSubMatrixBorders(0, cursensorinfo.columns/2-2, 0, cursensorinfo.rows, false); runcodesuffix += "_A0_"; humanreadablesuffix += "A0, "; } else if (matrixname.EqualTo("A1")) { - setSubMatrixBorders(sensorinfocur.columns/2+2, sensorinfocur.columns, 0, sensorinfocur.rows, false); + setSubMatrixBorders(cursensorinfo.columns/2+2, cursensorinfo.columns, 0, cursensorinfo.rows, false); runcodesuffix += "_A1_"; humanreadablesuffix += "A1, "; } @@ -726,7 +650,7 @@ void Run::getVetoPeakPositionFromFe55Run() void Run::constructUpdateString(string *sqlupdatequery, const string databasevaluename, const Double_t value, const int precision=3) { - // cout << colorred << databasevaluename << " : " << value << endlr; +// cout << colorred << databasevaluename << " : " << value << endlr; if (!std::isinf(value)) { if (value>0) @@ -742,18 +666,18 @@ void Run::updateDatabase() { string sqlupdatequery = ""; - constructUpdateString(&sqlupdatequery, "Gain", gain); - constructUpdateString(&sqlupdatequery, "SumPeak", histogram.posSum, 4); - constructUpdateString(&sqlupdatequery, "SeedPeak", histogram.posSeed, 4); - constructUpdateString(&sqlupdatequery, "VetoPeak", histogram.posVeto, 4); - constructUpdateString(&sqlupdatequery, "Avg.Noise", histogramCalibrated.avgNoise); - constructUpdateString(&sqlupdatequery, "Avg.Noise+", histogramCalibrated.avgNoisePlus, 2); - constructUpdateString(&sqlupdatequery, "Avg.Noise-", histogramCalibrated.avgNoiseMinus, 2); - constructUpdateString(&sqlupdatequery, "CCE_1", CCE_in_Perc_1); - constructUpdateString(&sqlupdatequery, "CCE_25", CCE_in_Perc_25); - constructUpdateString(&sqlupdatequery, "Avg.NoiseADC", histogram.avgNoise); + constructUpdateString(&sqlupdatequery, "Gain", histogram->gain); + constructUpdateString(&sqlupdatequery, "SumPeak", histogram->posSum, 4); + constructUpdateString(&sqlupdatequery, "SeedPeak", histogram->posSeed, 4); + constructUpdateString(&sqlupdatequery, "VetoPeak", histogram->posVeto, 4); + constructUpdateString(&sqlupdatequery, "Avg.Noise", histogram->calibrated->avgNoise); + constructUpdateString(&sqlupdatequery, "Avg.Noise+", histogram->calibrated->avgNoisePlus, 2); + constructUpdateString(&sqlupdatequery, "Avg.Noise-", histogram->calibrated->avgNoiseMinus, 2); + constructUpdateString(&sqlupdatequery, "CCE_1", histogram->CCE_in_Perc_1); + constructUpdateString(&sqlupdatequery, "CCE_25", histogram->CCE_in_Perc_25); + constructUpdateString(&sqlupdatequery, "Avg.NoiseADC", histogram->avgNoise); constructUpdateString(&sqlupdatequery, "Frames_found", frames_found,100000000); - constructUpdateString(&sqlupdatequery, "Sr90IntegralVal", histogramthresholdCalibrated.sr90IntegralVal,1000000000); + constructUpdateString(&sqlupdatequery, "Sr90IntegralVal", histogramfixedthreshold->calibrated->sr90IntegralVal,1000000000); if (sqlupdatequery.length()>0) { @@ -781,27 +705,21 @@ string Run::to_str_w_prec(const Float_t a_value, int precision = 3) return out.str(); } -Bool_t Run::binNoise() +Bool_t Run::binNoise(HistogramType* oneHistogramClass) { Float_t noise; TBranch* noiseBranch; Double_t const probabilities[] = {0.3415/2, 0.5, 1-0.3415/2}; // sigma/2 from gaus to the left and to the right //{0.17, 0.5, 1-0.17}; - histogram.Noise->Reset(); + oneHistogramClass->Noise->Reset(); processed->fNoiseTree->SetBranchAddress("noise", &noise, &noiseBranch); for (Int_t cnt=0; cntfNoiseTree->GetEntries(); cnt++) { processed->fNoiseTree->GetEntry(cnt); - histogram.Noise->Fill(noise); + oneHistogramClass->Noise->Fill(noise); } - histogram.Noise->GetQuantiles( 3, noisequantiles, probabilities); - histogram.avgNoise = noisequantiles[1]; - histogram.avgNoisePlus = noisequantiles[2] - noisequantiles[1]; - histogram.avgNoiseMinus = noisequantiles[1] - noisequantiles[0]; - histogramthreshold.avgNoise = histogram.avgNoise; - histogramthreshold.avgNoisePlus = histogram.avgNoisePlus; - histogramthreshold.avgNoiseMinus = histogram.avgNoiseMinus; - histogramfixedthreshold.avgNoise = histogram.avgNoise; - histogramfixedthreshold.avgNoisePlus = histogram.avgNoisePlus; - histogramfixedthreshold.avgNoiseMinus = histogram.avgNoiseMinus; + oneHistogramClass->Noise->GetQuantiles( 3, noisequantiles, probabilities); + oneHistogramClass->avgNoise = noisequantiles[1]; + oneHistogramClass->avgNoisePlus = noisequantiles[2] - noisequantiles[1]; + oneHistogramClass->avgNoiseMinus = noisequantiles[1] - noisequantiles[0]; // if (labbook.system == "PXI") // for (int j=0; j<3; j++) // noisequantiles[j] /= 16.0; // TODO analyze PXI scales @@ -809,12 +727,7 @@ Bool_t Run::binNoise() } Bool_t Run::binSeedSumVeto() -{ -// /// pixel number of seed pixel, position on sensor -// UInt_t seedPixel[10000]; -// /// Array of processed->clustersize * processed->clustersize clusters, seed pixel in the middle -// Float_t pixelcluster[processed->clustersize*processed->clustersize][10000]; - +{ /// collected charge in cluster Float_t pixelSum = 0; Float_t notSeedSum = 0; @@ -827,8 +740,8 @@ Bool_t Run::binSeedSumVeto() { for(Int_t hiti=0; (unsigned int)hitifFrameInfo.hits;hiti++) { - uint pixel_column_x = processed->fFrameInfo.pixel[hiti]%sensorinfocur.columns; // column of seed - uint pixel_row_y = processed->fFrameInfo.pixel[hiti]/sensorinfocur.columns; // row of seed + uint pixel_column_x = processed->fFrameInfo.pixel[hiti]%cursensorinfo.columns; // column of seed + uint pixel_row_y = processed->fFrameInfo.pixel[hiti]/cursensorinfo.columns; // row of seed if (false) { cout << "submatrix_x_start: " << submatrix_x_start << ", submatrix_x_end: " << submatrix_x_end << endl; cout << "submatrix_y_start: " << submatrix_y_start << ", submatrix_y_end: " << submatrix_y_end << endl; @@ -836,10 +749,10 @@ Bool_t Run::binSeedSumVeto() } if (pixel_column_x >= submatrix_x_start && pixel_column_x < submatrix_x_end && pixel_row_y >= submatrix_y_start && pixel_row_y < submatrix_y_end) // Diode sitzt oben im SeedPixel, da nach PitchY angeordnet { - histogram.numberofhits++; + histogram->numberofhits++; // histogram with the single pixel - histogram.Seed->Fill(processed->fFrameInfo.p[12][hiti]); + histogram->Seed->Fill(processed->fFrameInfo.p[12][hiti]); // sum histogram pixelSum = 0; @@ -870,48 +783,43 @@ Bool_t Run::binSeedSumVeto() notSeedSum += processed->fFrameInfo.p[clusteri][hiti]; } } - histogram.Sum->Fill(pixelSum); + histogram->Sum->Fill(pixelSum); // veto spectrum - if (TMath::Abs(notSeedSum) < systemparamcur.vetothreshold && (labbook.source.Contains("Fe")||labbook.source.Contains("Cd"))) - histogram.Veto->Fill(processed->fFrameInfo.p[12][hiti]); // histogram with the single pixel + if (TMath::Abs(notSeedSum) < cursystemparam.vetothreshold && (labbook.source.Contains("Fe")||labbook.source.Contains("Cd"))) + histogram->Veto->Fill(processed->fFrameInfo.p[12][hiti]); // histogram with the single pixel if (processed->fFrameInfo.pixelthreshold[hiti]>0) { - histogramthreshold.numberofhits++; + histogramthreshold->numberofhits++; - histogramthreshold.Seed->Fill(processed->fFrameInfo.p[12][hiti]); - histogramthreshold.Sum->Fill(pixelSum); - if (TMath::Abs(notSeedSum) < systemparamcur.vetothreshold && (labbook.source.Contains("Fe")||labbook.source.Contains("Cd"))) - histogramthreshold.Veto->Fill(processed->fFrameInfo.p[12][hiti]); // histogram with the single pixel + histogramthreshold->Seed->Fill(processed->fFrameInfo.p[12][hiti]); + histogramthreshold->Sum->Fill(pixelSum); + if (TMath::Abs(notSeedSum) < cursystemparam.vetothreshold && (labbook.source.Contains("Fe")||labbook.source.Contains("Cd"))) + histogramthreshold->Veto->Fill(processed->fFrameInfo.p[12][hiti]); // histogram with the single pixel } if (processed->fFrameInfo.pixelfixedthreshold[hiti]>0) { - histogramfixedthreshold.numberofhits++; + histogramfixedthreshold->numberofhits++; - histogramfixedthreshold.Seed->Fill(processed->fFrameInfo.p[12][hiti]); - histogramfixedthreshold.Sum->Fill(pixelSum); - if (TMath::Abs(notSeedSum) < systemparamcur.vetothreshold && (labbook.source.Contains("Fe")||labbook.source.Contains("Cd"))) - histogramfixedthreshold.Veto->Fill(processed->fFrameInfo.p[12][hiti]); // histogram with the single pixel + histogramfixedthreshold->Seed->Fill(processed->fFrameInfo.p[12][hiti]); + histogramfixedthreshold->Sum->Fill(pixelSum); + if (TMath::Abs(notSeedSum) < cursystemparam.vetothreshold && (labbook.source.Contains("Fe")||labbook.source.Contains("Cd"))) + histogramfixedthreshold->Veto->Fill(processed->fFrameInfo.p[12][hiti]); // histogram with the single pixel } } } } } // gROOT->SetBatch(kTRUE); - if (labbook.source.Contains("Fe")||labbook.source.Contains("Cd")) - histogram.posVeto=FitPerform(&histogram, histogram.Veto, "gaus", true); - histogram.posSeed=FitPerform(&histogram, histogram.Seed, "landau", true); - histogram.posSum=FitPerform(&histogram, histogram.Sum, "gaus", true); - if (labbook.source.Contains("Fe")||labbook.source.Contains("Cd")) - histogramthreshold.posVeto=FitPerform(&histogram, histogram.Veto, "gaus", false); - histogramthreshold.posSeed=FitPerform(&histogram, histogram.Seed, "landau", false); - histogramthreshold.posSum=FitPerform(&histogram, histogram.Sum, "gaus", false); - if (labbook.source.Contains("Fe")||labbook.source.Contains("Cd")) - histogramfixedthreshold.posVeto=FitPerform(&histogramfixedthreshold, histogramfixedthreshold.Veto, "gaus", false); - histogramfixedthreshold.posSeed=FitPerform(&histogramfixedthreshold, histogramfixedthreshold.Seed, "landau", false); - histogramfixedthreshold.posSum=FitPerform(&histogramfixedthreshold, histogramfixedthreshold.Sum, "gaus", false); + + for (vector::iterator curHistogramClass = HistogramClassVector.begin(); curHistogramClass != HistogramClassVector.end(); curHistogramClass++) { + if (labbook.source.Contains("Fe")||labbook.source.Contains("Cd")) + (*curHistogramClass)->posVeto=(*curHistogramClass)->FitPerform((*curHistogramClass)->Veto, "gaus", true); + (*curHistogramClass)->posSeed=(*curHistogramClass)->FitPerform((*curHistogramClass)->Seed, "landau", true); + (*curHistogramClass)->posSum=(*curHistogramClass)->FitPerform((*curHistogramClass)->Sum, "gaus", true); + } // gROOT->SetBatch(kFALSE); return 0; } @@ -925,10 +833,10 @@ Bool_t Run::binCluster() findRotatedClusterDimension(rotateangle); } else - { + { initClusters(curpixelinfo.pitchY, 0,curpixelinfo.pitchY*5,curpixelinfo.pitchX,0,curpixelinfo.pitchX*5); // prefix = "Cluster" + suffix; - // histogramstructpointer->histAvgCluster = new TH2F(prefix.Data(), Form("Avg. cluster distribution; x; y %s", humanreadablestr.Data()),5, 0,curpixelinfo.pitchY*5,5,0,curpixelinfo.pitchX*5); // TODO: remove hardcoded 5 + // HistogramTypepointer->histAvgCluster = new TH2F(prefix.Data(), Form("Avg. cluster distribution; x; y %s", humanreadablestr.Data()),5, 0,curpixelinfo.pitchY*5,5,0,curpixelinfo.pitchX*5); // TODO: remove hardcoded 5 } @@ -942,8 +850,8 @@ Bool_t Run::binCluster() { for(Int_t hiti=0; (unsigned int)hitifFrameInfo.hits;hiti++) { - uint pixel_column_x = processed->fFrameInfo.pixel[hiti]%sensorinfocur.columns; // column of seed - uint pixel_row_y = processed->fFrameInfo.pixel[hiti]/sensorinfocur.columns; // row of seed + uint pixel_column_x = processed->fFrameInfo.pixel[hiti]%cursensorinfo.columns; // column of seed + uint pixel_row_y = processed->fFrameInfo.pixel[hiti]/cursensorinfo.columns; // row of seed if (false) { cout << "submatrix_x_start: " << submatrix_x_start << ", submatrix_x_end: " << submatrix_x_end << endl; cout << "submatrix_y_start: " << submatrix_y_start << ", submatrix_y_end: " << submatrix_y_end << endl; @@ -975,11 +883,11 @@ Bool_t Run::binCluster() } // cout << colorwhite << "oben: xcoord: " << xcoord << ", ycoord: " << ycoord << endl; - histogram.histAvgCluster->Fill(xcoord,ycoord,processed->fFrameInfo.p[clusteri][hiti]); + histogram->histAvgCluster->Fill(xcoord,ycoord,processed->fFrameInfo.p[clusteri][hiti]); if (processed->fFrameInfo.pixelthreshold[hiti]>0) - histogramthreshold.histAvgCluster->Fill(xcoord,ycoord,processed->fFrameInfo.p[clusteri][hiti]); + histogramthreshold->histAvgCluster->Fill(xcoord,ycoord,processed->fFrameInfo.p[clusteri][hiti]); if (processed->fFrameInfo.pixelfixedthreshold[hiti]>0) - histogramfixedthreshold.histAvgCluster->Fill(xcoord,ycoord,processed->fFrameInfo.p[clusteri][hiti]); + histogramfixedthreshold->histAvgCluster->Fill(xcoord,ycoord,processed->fFrameInfo.p[clusteri][hiti]); } } else // Diode sitzt unten im SeedPixel, da nach PitchY angeordnet @@ -1006,11 +914,11 @@ Bool_t Run::binCluster() } // cout << colorwhite << "unten: xcoord: " << xcoord << ", ycoord: " << ycoord << endl; - histogram.histAvgCluster->Fill(xcoord,ycoord,processed->fFrameInfo.p[clusteri][hiti]); + histogram->histAvgCluster->Fill(xcoord,ycoord,processed->fFrameInfo.p[clusteri][hiti]); if (processed->fFrameInfo.pixelthreshold[hiti]>0) - histogramthreshold.histAvgCluster->Fill(xcoord,ycoord,processed->fFrameInfo.p[clusteri][hiti]); + histogramthreshold->histAvgCluster->Fill(xcoord,ycoord,processed->fFrameInfo.p[clusteri][hiti]); if (processed->fFrameInfo.pixelfixedthreshold[hiti]>0) - histogramfixedthreshold.histAvgCluster->Fill(xcoord,ycoord,processed->fFrameInfo.p[clusteri][hiti]); + histogramfixedthreshold->histAvgCluster->Fill(xcoord,ycoord,processed->fFrameInfo.p[clusteri][hiti]); } } } @@ -1018,11 +926,11 @@ Bool_t Run::binCluster() } } - for (Int_t clusterx = 1; clusterx <=histogram.histAvgCluster->GetXaxis()->GetNbins();clusterx++) { - for (Int_t clustery = 1; clustery <=histogram.histAvgCluster->GetYaxis()->GetNbins();clustery++) { - histogram.histAvgCluster->SetBinContent(clusterx,clustery,histogram.histAvgCluster->GetBinContent(clusterx,clustery)/histogram.numberofhits); - histogramthreshold.histAvgCluster->SetBinContent(clusterx,clustery,histogramthreshold.histAvgCluster->GetBinContent(clusterx,clustery)/histogramthreshold.numberofhits); - histogramfixedthreshold.histAvgCluster->SetBinContent(clusterx,clustery,histogramfixedthreshold.histAvgCluster->GetBinContent(clusterx,clustery)/histogramfixedthreshold.numberofhits); + for (Int_t clusterx = 1; clusterx <=histogram->histAvgCluster->GetXaxis()->GetNbins();clusterx++) { + for (Int_t clustery = 1; clustery <=histogram->histAvgCluster->GetYaxis()->GetNbins();clustery++) { + histogram->histAvgCluster->SetBinContent(clusterx,clustery,histogram->histAvgCluster->GetBinContent(clusterx,clustery)/histogram->numberofhits); + histogramthreshold->histAvgCluster->SetBinContent(clusterx,clustery,histogramthreshold->histAvgCluster->GetBinContent(clusterx,clustery)/histogramthreshold->numberofhits); + histogramfixedthreshold->histAvgCluster->SetBinContent(clusterx,clustery,histogramfixedthreshold->histAvgCluster->GetBinContent(clusterx,clustery)/histogramfixedthreshold->numberofhits); } } @@ -1141,84 +1049,33 @@ Bool_t Run:: findRotatedClusterDimension(Float_t rotateangle) return 0; } -Bool_t Run::plotNoise() -{ - if (!error) - { - TString legendEntry = Form("Noise: %.2f + %.2f - %.2f", histogram.avgNoise, histogram.avgNoisePlus, histogram.avgNoiseMinus ); - noisecanvas = plot1DHistogram(&histogram, histogram.Noise, "landau", "", legendEntry); - return 0; - } - return 1; -} - -Bool_t Run::plotSeed() -{ - if (!error) - { - plot1DHistogram(&histogram, histogram.Seed, "landau"); - return 0; - } - return 1; -} - -Bool_t Run::plotSeedThreshold() -{ - if (!error) - { - plot1DHistogram(&histogramthreshold, histogramthreshold.Seed, "landau"); - return 0; - } - return 1; -} - -Bool_t Run::plotSeedThresholdCalibrated() -{ - if (!error) - { - plot1DHistogram(&histogramthresholdCalibrated, histogramthresholdCalibrated.Seed, "landau"); - return 0; - } - return 1; -} - -Bool_t Run::plotSum() -{ - if (!error) - { - plot1DHistogram(&histogram, histogram.Sum, "gaus"); - return 0; - } - return 1; -} - -Bool_t Run::plotVeto() -{ - if (!error) - { - plot1DHistogram(&histogram, histogram.Veto, "gaus"); - return 0; - } - return 1; -} - +// Bool_t Run::plotNoise() +// { +// if (!error) +// { +// TString legendEntry = Form("Noise: %.2f + %.2f - %.2f", histogram->avgNoise, histogram->avgNoisePlus, histogram->avgNoiseMinus ); +// noisecanvas = plot1DHistogram(histogram, histogram->Noise, "landau", "", legendEntry); +// return 0; +// } +// return 1; +// } Bool_t Run::plotCompareHistograms() { if (!error) { - if ( !plothistogramstructpointerarray.size() ) // XOR operator != + if ( !compareHistogramClassVector.size() ) // XOR operator != { - cout << colorred << "plotCompareHistograms(): No plots to compare, please push them into the vector: Run::plothistogramstructpointerarray." << endlr; + cout << colorred << "plotCompareHistograms(): No plots to compare, please push them into the vector: Run::compareHistogramClassVector." << endlr; return 1; } Bool_t calibrated = true; Bool_t uncalibrated = true; - for (UInt_t histogrami=0; histogrami < plothistogramstructpointerarray.size(); histogrami++) + for (UInt_t histogrami=0; histogrami < compareHistogramClassVector.size(); histogrami++) { - calibrated = calibrated && plothistogramstructpointerarray.at(histogrami)->calibrated; - uncalibrated = uncalibrated && !plothistogramstructpointerarray.at(histogrami)->calibrated; - // cout << colorcyan << plothistogramstructpointerarray.at(histogrami)->histogramdescription << endlr; + calibrated = calibrated && compareHistogramClassVector.at(histogrami)->calibrated; + uncalibrated = uncalibrated && !compareHistogramClassVector.at(histogrami)->calibrated; + // cout << colorcyan << compareHistogramClassVector.at(histogrami)->histogramdescription << endlr; } // cout << "Calibrated " << (calibrated?"Yes":"No") << endl; // cout << "Uncalibrated " << (uncalibrated?"Yes":"No") << endl; @@ -1229,7 +1086,7 @@ Bool_t Run::plotCompareHistograms() } // legend entries - Float_t height = plothistogramstructpointerarray.size() * 0.04; + Float_t height = compareHistogramClassVector.size() * 0.04; TLegend* leg1 = new TLegend(0.4,0.89-height,0.89,0.89);//(0.6,0.7,0.89,0.89); leg1->SetTextSize(0.02); leg1->SetFillColor(0); leg1->SetBorderSize(0); @@ -1247,12 +1104,12 @@ Bool_t Run::plotCompareHistograms() TString canvasname = Form("%s%d",runcode.Data(),random); TCanvas* canvas = new TCanvas(canvasname, canvastitle, 1200, 800); canvas->Divide(2,2); - for (UInt_t histogrami=0; histogrami < plothistogramstructpointerarray.size(); histogrami++) + for (UInt_t histogrami=0; histogrami < compareHistogramClassVector.size(); histogrami++) { - histogramstruct* curhistogramstructp = plothistogramstructpointerarray.at(histogrami); + HistogramType* curhistogramclassp = compareHistogramClassVector.at(histogrami); TH1F* curhistogramclone; canvas->cd(1); - curhistogramclone = (TH1F*) curhistogramstructp->Seed->Clone(); + curhistogramclone = (TH1F*) curhistogramclassp->Seed->Clone(); curhistogramclone->SetLineColor(rootcolors[histogrami]); curhistogramclone->Draw("SAME"); legendEntry = Form("%s", curhistogramclone->GetTitle()); @@ -1262,7 +1119,7 @@ Bool_t Run::plotCompareHistograms() curhistogramclone->SetAxisRange(0,lastbin1*1.1,"X"); gPad->SetLogy(1); canvas->cd(2); - curhistogramclone = (TH1F*) curhistogramstructp->Sum->Clone(); + curhistogramclone = (TH1F*) curhistogramclassp->Sum->Clone(); curhistogramclone->SetLineColor(rootcolors[histogrami]); curhistogramclone->Draw("SAME"); leg1->Draw("SAME"); @@ -1270,17 +1127,17 @@ Bool_t Run::plotCompareHistograms() curhistogramclone->SetAxisRange(0,lastbin2*1.1,"X"); gPad->SetLogy(1); canvas->cd(3); - curhistogramclone = (TH1F*) curhistogramstructp->Veto->Clone(); + curhistogramclone = (TH1F*) curhistogramclassp->Veto->Clone(); curhistogramclone->SetLineColor(rootcolors[histogrami]); curhistogramclone->Draw("SAME"); leg1->Draw("SAME"); lastbin3 = (curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,1))>lastbin3)?curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,1)):lastbin3; curhistogramclone->SetAxisRange(0,lastbin3*1.1,"X"); canvas->cd(4); - curhistogramclone = (TH1F*) curhistogramstructp->Noise->Clone(); + curhistogramclone = (TH1F*) curhistogramclassp->Noise->Clone(); curhistogramclone->SetLineColor(rootcolors[histogrami]); curhistogramclone->Draw(); - legendEntry = Form("Noise: %.2f + %.2f - %.2f",curhistogramstructp->avgNoise, curhistogramstructp->avgNoisePlus, curhistogramstructp->avgNoiseMinus); + legendEntry = Form("Noise: %.2f + %.2f - %.2f",curhistogramclassp->avgNoise, curhistogramclassp->avgNoisePlus, curhistogramclassp->avgNoiseMinus); leg2->AddEntry(curhistogramclone, legendEntry, "l"); leg2->Draw(); } @@ -1290,39 +1147,39 @@ Bool_t Run::plotCompareHistograms() return 1; } -Bool_t Run::plotAllHistograms(histogramstruct* histogramstructpointer) +Bool_t Run::plotAllHistograms(HistogramType* HistogramTypepointer) { if (!error) { Float_t lastbin; Int_t random = random1->Rndm()*1000000; - TString canvastitle = Form("%s %s", runcode.Data(), histogramstructpointer->histogramdescription.Data()); - if (histogramstructpointer->calibrated) + TString canvastitle = Form("%s %s", runcode.Data(), HistogramTypepointer->histogramdescription.Data()); + if (HistogramTypepointer->calibrated) canvastitle += "_el"; TString canvasname = Form("%s%d",runcode.Data(),random); TCanvas* canvas = new TCanvas(canvasname, canvastitle, 1200, 800); canvas->Divide(2,2); canvas->cd(1); - histogramstructpointer->Seed->Draw(""); - plotVerticalLine(histogramstructpointer->Seed, histogramstructpointer->posSeed); - lastbin = histogramstructpointer->Seed->GetBinCenter(histogramstructpointer->Seed->FindLastBinAbove(2,1)); - histogramstructpointer->Seed->SetAxisRange(0,lastbin*1.1,"X"); + HistogramTypepointer->Seed->Draw(""); + plotVerticalLine(HistogramTypepointer->Seed, HistogramTypepointer->posSeed); + lastbin = HistogramTypepointer->Seed->GetBinCenter(HistogramTypepointer->Seed->FindLastBinAbove(2,1)); + HistogramTypepointer->Seed->SetAxisRange(0,lastbin*1.1,"X"); // gPad->SetLogy(1); canvas->cd(2); - histogramstructpointer->Sum->Draw(); - plotVerticalLine(histogramstructpointer->Sum, histogramstructpointer->posSum); - lastbin = histogramstructpointer->Sum->GetBinCenter(histogramstructpointer->Sum->FindLastBinAbove(2,1)); - histogramstructpointer->Sum->SetAxisRange(0,lastbin*1.1,"X"); + HistogramTypepointer->Sum->Draw(); + plotVerticalLine(HistogramTypepointer->Sum, HistogramTypepointer->posSum); + lastbin = HistogramTypepointer->Sum->GetBinCenter(HistogramTypepointer->Sum->FindLastBinAbove(2,1)); + HistogramTypepointer->Sum->SetAxisRange(0,lastbin*1.1,"X"); canvas->cd(3); - histogramstructpointer->Veto->Draw(); + HistogramTypepointer->Veto->Draw(); if (labbook.source.Contains("Fe")||labbook.source.Contains("Cd")) - plotVerticalLine(histogramstructpointer->Veto, histogramstructpointer->posVeto); - histogramstructpointer->Veto->SetAxisRange(histogramstructpointer->posVeto*0.7,histogramstructpointer->posVeto*1.4,"X"); + plotVerticalLine(HistogramTypepointer->Veto, HistogramTypepointer->posVeto); + HistogramTypepointer->Veto->SetAxisRange(HistogramTypepointer->posVeto*0.7,HistogramTypepointer->posVeto*1.4,"X"); canvas->cd(4); - histogramstructpointer->Noise->Draw(); - TString legendEntry = Form("Noise: %.2f + %.2f - %.2f",histogramstructpointer->avgNoise, histogramstructpointer->avgNoisePlus, histogramstructpointer->avgNoiseMinus); + HistogramTypepointer->Noise->Draw(); + TString legendEntry = Form("Noise: %.2f + %.2f - %.2f",HistogramTypepointer->avgNoise, HistogramTypepointer->avgNoisePlus, HistogramTypepointer->avgNoiseMinus); TLegend* leg = new TLegend(0.5,0.8,0.89,0.89);//(0.6,0.7,0.89,0.89); - leg->AddEntry(histogramstructpointer->Veto, legendEntry, "l"); + leg->AddEntry(HistogramTypepointer->Veto, legendEntry, "l"); leg->SetTextSize(0.03); leg->Draw(); @@ -1346,13 +1203,13 @@ Bool_t Run::plotAllHistograms(histogramstruct* histogramstructpointer) return 1; } -Bool_t Run::plotClusterDistribution(histogramstruct* histogramstructpointer) +Bool_t Run::plotClusterDistribution(HistogramType* HistogramTypepointer) { Int_t random = random1->Rndm()*100000; TString canvastitle = Form("%s_Cluster", runcode.Data()); - if (histogramstructpointer->calibrated) + if (HistogramTypepointer->calibrated) canvastitle += "_el"; - if (histogramstructpointer->thresholdcluster) + if (HistogramTypepointer->isthresholdclustertype) canvastitle += "_trsh"; TString canvasname = Form("%s%d",runcode.Data(),random); TCanvas* Canvas_1 = new TCanvas(canvasname, canvastitle, 1200, 800); @@ -1365,9 +1222,9 @@ Bool_t Run::plotClusterDistribution(histogramstruct* histogramstructpointer) // Canvas_1->SetPhi(32.46973); Canvas_1->TAttPad::SetFrameBorderMode(0); Canvas_1->TAttPad::SetFrameBorderMode(0); - histogramstructpointer->histAvgCluster->SetContour(20); + HistogramTypepointer->histAvgCluster->SetContour(20); - TPaletteAxis *palette = new TPaletteAxis(0.8596098,-0.8025656,0.8949913,0.9836416,histogramstructpointer->histAvgCluster); + TPaletteAxis *palette = new TPaletteAxis(0.8596098,-0.8025656,0.8949913,0.9836416,HistogramTypepointer->histAvgCluster); palette->SetLabelColor(1); palette->SetLabelFont(42); palette->SetLabelOffset(0.005); @@ -1376,7 +1233,7 @@ Bool_t Run::plotClusterDistribution(histogramstruct* histogramstructpointer) palette->SetTitleSize(0.035); palette->SetFillColor(100); palette->SetFillStyle(1001); - histogramstructpointer->histAvgCluster->GetListOfFunctions()->Add(palette,"br"); + HistogramTypepointer->histAvgCluster->GetListOfFunctions()->Add(palette,"br"); TPaveStats *ptstats = new TPaveStats(0.7013721,0.7794677,0.850686,0.9226869,"brNDC"); ptstats->SetTextFont(42); @@ -1384,37 +1241,37 @@ Bool_t Run::plotClusterDistribution(histogramstruct* histogramstructpointer) ptstats->SetOptStat(1011); ptstats->SetOptFit(0); ptstats->Draw(); - histogramstructpointer->histAvgCluster->GetListOfFunctions()->Add(ptstats); - ptstats->SetParent(histogramstructpointer->histAvgCluster); + HistogramTypepointer->histAvgCluster->GetListOfFunctions()->Add(ptstats); + ptstats->SetParent(HistogramTypepointer->histAvgCluster); Int_t ci; // for color index setting ci = TColor::GetColor("#000099"); - histogramstructpointer->histAvgCluster->SetLineColor(ci); - histogramstructpointer->histAvgCluster->GetXaxis()->SetTitle(" x [#mum]"); - histogramstructpointer->histAvgCluster->GetXaxis()->SetLabelFont(42); - histogramstructpointer->histAvgCluster->GetXaxis()->SetLabelOffset(0.004); - histogramstructpointer->histAvgCluster->GetXaxis()->SetLabelSize(0.025); - histogramstructpointer->histAvgCluster->GetXaxis()->SetTitleSize(0.035); - histogramstructpointer->histAvgCluster->GetXaxis()->SetTitleFont(42); - histogramstructpointer->histAvgCluster->GetXaxis()->SetTitleOffset(1.5); - histogramstructpointer->histAvgCluster->GetYaxis()->SetTitle(" y [#mum]"); - histogramstructpointer->histAvgCluster->GetYaxis()->SetLabelFont(42); - histogramstructpointer->histAvgCluster->GetYaxis()->SetLabelOffset(-0.016); - histogramstructpointer->histAvgCluster->GetYaxis()->SetLabelSize(0.025); - histogramstructpointer->histAvgCluster->GetYaxis()->SetTitleSize(0.035); - histogramstructpointer->histAvgCluster->GetYaxis()->SetTitleFont(42); - histogramstructpointer->histAvgCluster->GetYaxis()->SetTitleOffset(1.5); - if (histogramstructpointer->calibrated) - histogramstructpointer->histAvgCluster->GetZaxis()->SetTitle(" Charge [e]"); + HistogramTypepointer->histAvgCluster->SetLineColor(ci); + HistogramTypepointer->histAvgCluster->GetXaxis()->SetTitle(" x [#mum]"); + HistogramTypepointer->histAvgCluster->GetXaxis()->SetLabelFont(42); + HistogramTypepointer->histAvgCluster->GetXaxis()->SetLabelOffset(0.004); + HistogramTypepointer->histAvgCluster->GetXaxis()->SetLabelSize(0.025); + HistogramTypepointer->histAvgCluster->GetXaxis()->SetTitleSize(0.035); + HistogramTypepointer->histAvgCluster->GetXaxis()->SetTitleFont(42); + HistogramTypepointer->histAvgCluster->GetXaxis()->SetTitleOffset(1.5); + HistogramTypepointer->histAvgCluster->GetYaxis()->SetTitle(" y [#mum]"); + HistogramTypepointer->histAvgCluster->GetYaxis()->SetLabelFont(42); + HistogramTypepointer->histAvgCluster->GetYaxis()->SetLabelOffset(-0.016); + HistogramTypepointer->histAvgCluster->GetYaxis()->SetLabelSize(0.025); + HistogramTypepointer->histAvgCluster->GetYaxis()->SetTitleSize(0.035); + HistogramTypepointer->histAvgCluster->GetYaxis()->SetTitleFont(42); + HistogramTypepointer->histAvgCluster->GetYaxis()->SetTitleOffset(1.5); + if (HistogramTypepointer->calibrated) + HistogramTypepointer->histAvgCluster->GetZaxis()->SetTitle(" Charge [e]"); else - histogramstructpointer->histAvgCluster->GetZaxis()->SetTitle(" Charge [ADU]"); - histogramstructpointer->histAvgCluster->GetZaxis()->SetTitleOffset(1.5); - histogramstructpointer->histAvgCluster->GetZaxis()->SetLabelFont(42); - histogramstructpointer->histAvgCluster->GetZaxis()->SetLabelSize(0.025); - histogramstructpointer->histAvgCluster->GetZaxis()->SetTitleSize(0.035); - histogramstructpointer->histAvgCluster->GetZaxis()->SetTitleFont(42); - histogramstructpointer->histAvgCluster->GetZaxis()->SetNdivisions(5,5,0); - histogramstructpointer->histAvgCluster->Draw("LEGO2Z "); + HistogramTypepointer->histAvgCluster->GetZaxis()->SetTitle(" Charge [ADU]"); + HistogramTypepointer->histAvgCluster->GetZaxis()->SetTitleOffset(1.5); + HistogramTypepointer->histAvgCluster->GetZaxis()->SetLabelFont(42); + HistogramTypepointer->histAvgCluster->GetZaxis()->SetLabelSize(0.025); + HistogramTypepointer->histAvgCluster->GetZaxis()->SetTitleSize(0.035); + HistogramTypepointer->histAvgCluster->GetZaxis()->SetTitleFont(42); + HistogramTypepointer->histAvgCluster->GetZaxis()->SetNdivisions(5,5,0); + HistogramTypepointer->histAvgCluster->Draw("LEGO2Z "); TPaveText *pt = new TPaveText(0.3177154,0.9348119,0.6822846,0.995,"blNDC"); pt->SetName("title"); @@ -1437,28 +1294,28 @@ Bool_t Run::plotClusterDistribution(histogramstruct* histogramstructpointer) f->Append(img); f->Write(); - write2DHistogramToFile(histogramstructpointer->histAvgCluster); + write2DHistogramToFile(HistogramTypepointer->histAvgCluster); // gStyle->SetPaperSize(10.,10.); // canvas->Print(savepathresults + "/" + canvastitle + ".tex"); // create lorentz fit of a slice - TH1F* xslicetroughcluster = new TH1F(Form("%s_xslice", histogramstructpointer->histAvgCluster->GetTitle()),"X slice",histogramstructpointer->histAvgCluster->GetXaxis()->GetNbins(),histogramstructpointer->histAvgCluster->GetXaxis()->GetBinLowEdge(1),histogramstructpointer->histAvgCluster->GetXaxis()->GetBinUpEdge(histogramstructpointer->histAvgCluster->GetXaxis()->GetNbins())); + 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[4]; - Int_t middlebin = (int)(histogramstructpointer->histAvgCluster->GetXaxis()->GetNbins()/2.0 + 0.5); - for (Int_t bini=1; bini <= histogramstructpointer->histAvgCluster->GetXaxis()->GetNbins(); bini++) - xslicetroughcluster->SetBinContent(bini,histogramstructpointer->histAvgCluster->GetBinContent(bini,middlebin)); + 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", histogramstructpointer->histAvgCluster->GetTitle()),"Y slice",histogramstructpointer->histAvgCluster->GetYaxis()->GetNbins(),histogramstructpointer->histAvgCluster->GetYaxis()->GetBinLowEdge(1),histogramstructpointer->histAvgCluster->GetYaxis()->GetBinUpEdge(histogramstructpointer->histAvgCluster->GetYaxis()->GetNbins())); + 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[4]; - middlebin = (int)(histogramstructpointer->histAvgCluster->GetYaxis()->GetNbins()/2.0 + 0.5); - for (Int_t bini=1; bini <= histogramstructpointer->histAvgCluster->GetYaxis()->GetNbins(); bini++) - yslicetroughcluster->SetBinContent(bini,histogramstructpointer->histAvgCluster->GetBinContent(middlebin,bini)); + 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)); random = random1->Rndm()*100000; canvastitle = Form("%s Cluster slices", runcode.Data()); - if (histogramstructpointer->calibrated) + if (HistogramTypepointer->calibrated) canvastitle += "_el"; - if (histogramstructpointer->thresholdcluster) + if (HistogramTypepointer->isthresholdclustertype) canvastitle += "_trsh"; canvasname = Form("%s%d",runcode.Data(),random); TCanvas* canvas = new TCanvas(canvasname, canvastitle, 1200, 700); @@ -1469,13 +1326,13 @@ Bool_t Run::plotClusterDistribution(histogramstruct* histogramstructpointer) canvas->cd(1); TPad* firstpart = (TPad*)canvas->GetPad(1); firstpart->SetGrid(); - FitPerform(&histogram, xslicetroughcluster,"lorentz", false, &xslicetroughclusterpar[0]); + HistogramTypepointer->FitPerform(xslicetroughcluster,"lorentz", false, &xslicetroughclusterpar[0]); canvas->cd(2); TPad* secondpart = (TPad*)canvas->GetPad(2); secondpart->SetGrid(); - FitPerform(&histogram, yslicetroughcluster,"lorentz", false, &yslicetroughclusterpar[0]); + HistogramTypepointer->FitPerform(yslicetroughcluster,"lorentz", false, &yslicetroughclusterpar[0]); canvas->Draw(); @@ -1491,322 +1348,6 @@ Bool_t Run::plotClusterDistribution(histogramstruct* histogramstructpointer) return 0; } -Bool_t Run::plotClusterDistribution() -{ - if (!error) - { - plotClusterDistribution(&histogram); - return 0; - } - return 1; -} - -Bool_t Run::plotAllHistograms() -{ - if (!error) - { - plotAllHistograms(&histogram); - return 0; - } - return 1; -} - -Bool_t Run::plotAllHistogramsCalibrated() -{ - if (!error) - { - if (histogramCalibrated.calibrated) { - plotAllHistograms(&histogramCalibrated); - return 0; - } - } - return 1; -} - -Bool_t Run::plotAllHistogramsThresholdCluster() -{ - if (!error) - { - plotAllHistograms(&histogramthreshold); - return 0; - } - return 1; -} - -Bool_t Run::plotAllHistogramsThresholdClusterCalibrated() -{ - if (!error) - { - plotAllHistograms(&histogramthresholdCalibrated); - return 0; - } - return 1; -} - -Bool_t Run::integrateSr90Spectra(histogramstruct* histogramstructpointer, TH1F* histogrampointer, Float_t thresholdborder, Bool_t verbose) -{ - Float_t posMaxValHist = histogrampointer->GetXaxis()->GetXmax(); - - TH1F* smoothedcurce = (TH1F*)histogrampointer->Clone(); - smoothedcurce->Smooth(4); - -// Int_t random = random1->Rndm()*100000; -// TString canvastitle = Form("%s %s sdfasdf", histogrampointer->GetName(), runcode.Data()); -// TString canvasname = Form("%s %s %d dfgsdfgsdfg", histogrampointer->GetName(), runcode.Data(), random); -// TCanvas* canvas = new TCanvas(canvasname, canvastitle, 900, 700); -// smoothedcurce->Draw(); - - Int_t thresholdbincurcandidate = 0; - if (thresholdborder < 0) - { - Int_t rising = 0; - Int_t bini =smoothedcurce->GetMaximumBin(); - Float_t curval = smoothedcurce->GetXaxis()->GetBinCenter(bini); -// cout << "smoothedcurce->GetXaxis()->GetBinCenter(" << bini << "): " << curval << endl; - - do { - curval=smoothedcurce->GetBinContent(bini++); - if (curval*0.95 <= smoothedcurce->GetBinContent(bini)) - { - rising++; -// cout << "rising at " << smoothedcurce->GetXaxis()->GetBinCenter(bini) << " as " << curval << " < " << smoothedcurce->GetBinContent(bini) << endl; // debug - } - else - { - rising = 0; - thresholdbincurcandidate = bini; - } - } while (rising < 3 && biniGetNbinsX()); - } - else - { - thresholdbincurcandidate = histogrampointer->GetXaxis()->FindBin(thresholdborder); - } - - -// histogrampointer->GetXaxis()->SetRange(histogrampointer->GetXaxis()->FindBin(0),histogrampointer->GetXaxis()->FindBin(posMaxValHist)); -// Float_t posNoiseMax= histogrampointer->GetMaximum(); -// histogrampointer->GetXaxis()->SetRange(posNoiseMax,histogrampointer->GetXaxis()->FindBin(posMaxValHist)); -// Float_t posChargeMax= histogrampointer->GetMaximum(); - - histogramstructpointer->noisethresholdborder = histogrampointer->GetXaxis()->GetBinUpEdge(thresholdbincurcandidate); - - if (verbose) { - cout << " Integrating " << histogrampointer->GetTitle() << endlr; - cout << " Noise threshold at " << noisethresholdborder; - TString histtitle=histogrampointer->GetXaxis()->GetTitle(); - if (histtitle.Contains("[e]")) - { - cout << " e" << endl; - if (noisethresholdborder > 300) { - cout << coloryellow << " This noise threshold seems too high, please check manually the " << histogrampointer->GetName() << " spectrum." << endlr; - } - } - else - { - cout << " ADU" << endl; - if (noisethresholdborder > 75) { - cout << coloryellow << " This noise threshold seems too high, please check manually the " << histogrampointer->GetName() << " spectrum." << endlr; - } - } - } - histogramstructpointer->sr90IntegralVal = histogrampointer->IntegralAndError(thresholdbincurcandidate,histogrampointer->GetNbinsX(), histogramstructpointer->sr90IntegralErr); - // cout << "Integrate from bin " << thresholdbincurcandidate << " to " << histogrampointer->GetNbinsX() << endl; - histogramstructpointer->sr90IntegralErr /= histogramstructpointer->sr90IntegralVal/100; - cout << " Unscaled integral is " << Form("%e",histogramstructpointer->sr90IntegralVal) << ", error " << histogramstructpointer->sr90IntegralErr << "%" << endl; - cout << " "; - if (labbook.frames_foundDB>0 || frames_found>0) - { - histogramstructpointer->sr90IntegralVal/=(labbook.frames_foundDB>0)?labbook.frames_foundDB:frames_found; - cout << "Scaled "; - } - cout << "Integral is " << Form("%e",histogramstructpointer->sr90IntegralVal) << ", error " << histogramstructpointer->sr90IntegralErr << "%" << endl; - - return 0; -} - -Float_t Run::FitPerform(histogramstruct* histogramstructpointer, TH1F* histogrampointer, TString fitFuncType, Bool_t verbose, Double_t* parameters) -{ - Float_t posMax = 0; - Float_t posMax2 = 0; - Float_t posMaxValHist = histogrampointer->GetXaxis()->GetXmax(); - // posMaxValHist/10 for USB system, the value is 90 - Float_t noiseborder = 90; - if (histogramstructpointer->noisethresholdborder>0) - noiseborder = histogramstructpointer->noisethresholdborder; - else if (noisethresholdborder > 0) - noiseborder = noisethresholdborder; - - if (doFits) - { - if (fitFuncType.Contains("gaus")) - { - 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); - - if (TString(histogrampointer->GetName()).Contains("Veto")) - { - Float_t peak1 = histogrampointer->GetMaximumBin(); - Float_t peak2; - if (histogrampointer->GetBinCenter(peak1)<3.3*noiseborder) // vermutlich ist veto peak nicht am höchsten - { - Float_t avg = 0; - for(Int_t bin=histogrampointer->GetXaxis()->FindBin(noiseborder);binFindLastBinAbove(0);bin++) - avg += histogrampointer->GetBinContent(bin); - avg /= histogrampointer->FindLastBinAbove(0) - histogrampointer->GetXaxis()->FindBin(noiseborder); - peak2 = histogrampointer->FindLastBinAbove(avg/3); - } - 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->Fit(fitFunc, "N,Q,W", "", histogrampointer->GetXaxis()->GetBinCenter(min), posMaxValHist); - posMax = histogrampointer->GetXaxis()->GetBinCenter(histogrampointer->GetMaximumBin()); // Methode 1 - 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", "", noiseborder, posMaxValHist); - posMax = fitFunc->GetMaximumX(); // Methode 2 - } - Float_t sigma = fitFunc->GetParameter(2); - posMax2 = fitFunc->GetMaximumX(); // Methode 2 - Float_t peakposdifperc = abs(posMax-posMax2)/min(posMax,posMax2); - if (sigma > 260 || peakposdifperc>0.3) - { - if (verbose) - { - if (sigma > 260) - { - cout << "Sigma suspiciously height when fitting " << histogrampointer->GetName() << " spectrum: " << sigma << endl; - cout << colorred << " Could not find " << histogrampointer->GetName() << " peak" << endlr; - } - else if (peakposdifperc>0.3) - { - cout << "Maximum peak position and fit gaussian peak position doesn't match in " << histogrampointer->GetName() << " spectrum. Difference: " << peakposdifperc*100 <<" % " << endl; - cout << colorred << " Could not find " << histogrampointer->GetName() << " peak" << endlr; - } - } - return 0; - } - else if (sigma > 40 || peakposdifperc>0.1) - { - if (verbose) - { - if (sigma > 260) - { - cout << "Sigma or suspiciously height when fitting " << histogrampointer->GetName() << " spectrum: " << sigma << endl; - cout << coloryellow << " Please check " << histogrampointer->GetName() << " peak position: " << colorreset << colorwhite << posMax << endlr; - } - else if (peakposdifperc>0.1) - { - cout << "Maximum peak position and fit gaussian peak position in " << histogrampointer->GetName() << " have difference of " << peakposdifperc*100 <<" % " << endl; - cout << coloryellow << " Please check " << histogrampointer->GetName() << " peak position: " << colorreset << colorwhite << posMax << endlr; - } - } - } -// fitFunc->DrawCopy("same"); - TString legendEntry = TString(Form("%s",runcode.Data())); - TLegend* leg = new TLegend(0.5,0.5,0.89,0.89);//(0.6,0.7,0.89,0.89); - // leg->SetHeader();//"Legend Title"); - leg->SetFillStyle(0); - leg->AddEntry((TObject*) 0, legendEntry, ""); - leg->SetTextSize(0.05); -// leg->Draw(); - - } - else if (fitFuncType=="landau") - { - 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); - - Float_t fitMax1 = 1000; - Float_t fitMax2 = 1000; - Float_t fitMax3 = 1000; - Float_t minFitMax = 1000; - Float_t maxFitMax = 1000; - histogrampointer->Fit(fitFunc, "N,Q,W", "", noiseborder, posMaxValHist); - fitMax1 = fitFunc->GetMaximumX(); -// fitFunc->DrawCopy("same"); - histogrampointer->Fit(fitFunc, "N,Q,W", "", 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); -// histogrampointer->Fit(fitFunc, "N,Q,W", "", fitMax1, histogrampointer->GetBinCenter(bini)); - fitMax3 = fitFunc->GetMaximumX(); - fitFunc->SetLineColor(kGreen); -// fitFunc->DrawCopy("same"); - fitFunc->SetLineStyle(1); // normal for the following fits - - // Sort the three fits and save error estimation - minFitMax = TMath::Min(TMath::Min(fitMax1,fitMax2),fitMax3); - maxFitMax = TMath::Max(TMath::Max(fitMax1,fitMax2),fitMax3); - posMax = fitMax1 + fitMax2 + fitMax3 - minFitMax - maxFitMax; - - //fitLandauErrorLeft.push_back(posMax - minFitMax); - //fitLandauErrorRight.push_back(maxFitMax - posMax); - } else if (fitFuncType=="lorentz") - { - // create a TF1 with the range from 0 to 3 and 6 parameters - TF1 *fitFcn = new TF1("fitFcn",lorentzianPeak,0,160,4); - fitFcn->SetNpx(1000); - fitFcn->SetLineWidth(4); - fitFcn->SetLineColor(rootcolors[plotStyle]); - // set start values for some parameters - fitFcn->SetParName(0,"Area"); - fitFcn->SetParameter(0,21655); // Area - fitFcn->SetParName(1,"Width"); - fitFcn->SetParameter(1,34); // width - fitFcn->SetParName(2,"Peak pos."); - fitFcn->SetParameter(2,80); // peak position x - fitFcn->SetParName(3,"Y shift"); - fitFcn->SetParameter(3,-22); // y-shift - gStyle->SetOptFit(0111); - -// fitFcn->FixParameter(0,21655.90045); -// fitFcn->FixParameter(1,34.31885101); - fitFcn->FixParameter(2,histogrampointer->GetXaxis()->GetBinCenter((int)(histogrampointer->GetXaxis()->GetNbins()/2.0+0.5))); -// fitFcn->FixParameter(3,-22.43016284); - - //histogrampointer->Fit("fitFcn","V+","ep"); - histogrampointer->Fit("fitFcn","Q,W","ep"); - for (Int_t pari=0; pari<4; pari++) - parameters[pari]=fitFcn->GetParameter(pari); - if (verbose) - { - Printf("width: %.6f ",fitFcn->GetParameter(1)); - cout << "width: " << fitFcn->GetParameter(1) << endl; - cout << "peak: " << fitFcn->GetParameter(2) << endl; - cout << "y-shift:" << fitFcn->GetParameter(3) << endl; - } - } - - } - - return posMax; -} - -// Lorenzian Peak function -Double_t Run::lorentzianPeak(Double_t *x, Double_t *par) { - return par[3]+(0.5*par[0]*par[1]/TMath::Pi()) / - TMath::Max( 1.e-10,(x[0]-par[2])*(x[0]-par[2]) - + .25*par[1]*par[1]); -} void Run::plotVerticalLine(TH1F* histogrampointer, Float_t xVal) { if (xVal > 0) @@ -1829,7 +1370,7 @@ void Run::plotVerticalLine(TH1F* histogrampointer, Float_t xVal) { } } -TCanvas* Run::plot1DHistogram(histogramstruct* histogramstructpointer, TH1F* onehistogram, TString fitFuncType, TString titlestr, TString legendstr) +TCanvas* Run::plot1DHistogram(HistogramType* HistogramTypepointer, TH1F* onehistogram, TString fitFuncType, TString titlestr, TString legendstr) { if (onehistogram != nullptr) { @@ -1841,7 +1382,7 @@ TCanvas* Run::plot1DHistogram(histogramstruct* histogramstructpointer, TH1F* one TCanvas* canvas = new TCanvas(canvasname, canvastitle, 900, 700); onehistogram->SetTitle(titlestr); onehistogram->Draw(); - Float_t maxValue = FitPerform(histogramstructpointer, onehistogram, fitFuncType); + Float_t maxValue = HistogramTypepointer->FitPerform(onehistogram, fitFuncType); 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); @@ -1920,22 +1461,39 @@ Bool_t Run::writeAllHistogramsToFile() TString filename= savepathresults + "/" + runcode + " histograms.dat"; fstream* fout = new fstream(filename,ios::out); - TString header = Form("#%s %lu frames\n", runcode.Data(), (labbook.frames_foundDB>0)?labbook.frames_foundDB:frames_found); - header += Form("#bin [ADU]\tbin [e]\tSeed\tSum\tVeto\tSeed thrsh\tSum thrsh\tbin noise [ADU]\tbin noise [e]\tnoise\n"); - header += Form("#posVeto, run: %.1f, DB: %.1f, Fe55 DB (%d, %.1f): %.1f\n", histogram.posVeto, labbook.posVetoDB, Fe55run.posVetorunnumber, Fe55run.temperature, Fe55run.posVeto); - header += Form("#posSeed, run: %.1f, DB: %.1f\n", histogram.posSeed, labbook.posSeedDB); - header += Form("#posSum, run: %.1f, DB: %.1f\n", histogram.posSum, labbook.posSumDB); - header += Form("#Noise [ADU]: %.2f +%.2f -%.2f\n", histogram.avgNoise, histogram.avgNoisePlus, histogram.avgNoiseMinus ); - if (histogramCalibrated.calibrated) - header += Form("#Noise [e]: %.2f +%.2f -%.2f",histogramCalibrated.avgNoise, histogramCalibrated.avgNoisePlus, histogramCalibrated.avgNoiseMinus ); + TString header = Form("#%s %lu frames\n", runcode.Data(), frames_found); + header += Form("#posVeto, run: %.1f, DB: %.1f, Fe55 DB (%d, %.1f): %.1f\n", histogram->posVeto, labbook.posVetoDB, Fe55run.posVetorunnumber, Fe55run.temperature, Fe55run.posVeto); + header += Form("#posSeed, run: %.1f, DB: %.1f\n", histogram->posSeed, labbook.posSeedDB); + header += Form("#posSum, run: %.1f, DB: %.1f\n", histogram->posSum, labbook.posSumDB); + if (labbook.source.Contains("Sr90")) { + header += Form("#Sr90 integral: "); + for (vector::iterator curHistogramClass = HistogramClassVector.begin(); curHistogramClass != HistogramClassVector.end(); curHistogramClass++) + header += Form("%s: %.6f\t", (*curHistogramClass)->histogramdescription.Data(), histogram->sr90IntegralVal); + header += Form("\n"); + } + header += Form("#Noise [ADU]: %.2f +%.2f -%.2f\n", histogram->avgNoise, histogram->avgNoisePlus, histogram->avgNoiseMinus ); + if (histogram->calibrated != 0) + header += Form("#Noise [e]: %.2f +%.2f -%.2f\n",histogram->calibrated->avgNoise, histogram->calibrated->avgNoisePlus, histogram->calibrated->avgNoiseMinus ); + header += Form("#\t\t"); + for (vector::iterator curHistogramClass = HistogramClassVector.begin(); curHistogramClass != HistogramClassVector.end(); curHistogramClass++) { + header += Form("%s\t\t\t", (*curHistogramClass)->histogramdescription.Data()); + } + header += Form("\n#bin [ADU]\tbin [e]\t"); + for (vector::iterator curHistogramClass = HistogramClassVector.begin(); curHistogramClass != HistogramClassVector.end(); curHistogramClass++) { + header += Form("Seed\tSum\tVeto\t"); + } + header += Form("bin noise [ADU]\tbin noise [e]\tnoise"); *fout << header << endl; TString outline; - for(Int_t bin=0;binGetNbinsX();bin++) + for(Int_t bin=1;bin<=cursystemparam.nbins;bin++) { - outline=Form("%.1f\t%.1f\t%.1f\t%.1f\t%.1f\t%.1f\t%.1f", histogram.Seed->GetBinCenter(bin), histogramCalibrated.calibrated?histogramCalibrated.Seed->GetBinCenter(bin):0, histogram.Seed->GetBinContent(bin), histogram.Sum->GetBinContent(bin), histogram.Veto->GetBinContent(bin), histogramthreshold.Seed->GetBinContent(bin), histogramthreshold.Sum->GetBinContent(bin)); - if (bin < systemparamcur.nbinsnoise) - outline+=Form("\t%.1f\t%.1f\t%.1f", histogram.Noise->GetBinCenter(bin),histogramCalibrated.calibrated?histogramCalibrated.Noise->GetBinCenter(bin):0,histogram.Noise->GetBinContent(bin)); + outline=Form("%.1f\t%.1f\t", histogram->Seed->GetBinCenter(bin), histogram->iscalibrated?histogram->calibrated->Seed->GetBinCenter(bin):0); + for (vector::iterator curHistogramClass = HistogramClassVector.begin(); curHistogramClass != HistogramClassVector.end(); curHistogramClass++) { + outline+=Form("%.1f\t%.1f\t%.1f\t", (*curHistogramClass)->Seed->GetBinContent(bin), (*curHistogramClass)->Sum->GetBinContent(bin), (*curHistogramClass)->Veto->GetBinContent(bin)); + } + if (bin <= cursystemparam.nbinsnoise) + outline+=Form("\t%.1f\t%.1f\t%.1f", histogram->Noise->GetBinCenter(bin),histogram->iscalibrated?histogram->calibrated->Noise->GetBinCenter(bin):0,histogram->Noise->GetBinContent(bin)); *fout<close(); @@ -2044,49 +1602,6 @@ void Run::MakeGnuplotFile() //system("okular "+ ergebnisfile + ".png &"); } -void Run::initHistograms(histogramstruct* histogramstructpointer, TString suffix) -{ - histogramstructpointer->histogramdescription = suffix; - initHistogram(histogramstructpointer->Seed, "Seed " + suffix); - initHistogram(histogramstructpointer->Sum, "Sum " + suffix); - initHistogram(histogramstructpointer->Veto, "Veto " + suffix); - - TString prefix = "Noise " + suffix; - - if (suffix.Contains("threshold")) - { - histogramstructpointer->thresholdcluster = 1; - histogramstructpointer->Noise = histogram.Noise; - } - else - { - histogramstructpointer->Noise=new TH1F(prefix.Data(), Form("%s, %s", prefix.Data(), humanreadablestr.Data()), systemparamcur.nbinsnoise, 0, systemparamcur.maxbinnoise); - histogramstructpointer->Noise->SetLineStyle(rootlinestyle[plotStyle]); - histogramstructpointer->Noise->SetLineColor(rootcolors[plotStyle]); - histogramstructpointer->Noise->SetLineWidth(3); - histogramstructpointer->Noise->GetXaxis()->SetTitle("Q_coll [ADU]"); - histogramstructpointer->Noise->GetYaxis()->SetTitle(Form("Entries [1/%.1f ADU]",histogramstructpointer->Noise->GetBinWidth(1))); - histogramstructpointer->Noise->GetXaxis()->CenterTitle(); - histogramstructpointer->Noise->GetYaxis()->CenterTitle(); - } - plothistogrampointer = &histogram.Seed; - plothistogramstructpointer = &histogram; -} - -void Run::initHistogram(TH1F* &histogrampointer, TString prefix) -{ - histogrampointer=new TH1F(prefix.Data(), Form("%s, %s", prefix.Data(), humanreadablestr.Data()), systemparamcur.nbins, 0, systemparamcur.maxbin); - histogrampointer->SetLineStyle(rootlinestyle[plotStyle]); - histogrampointer->SetLineColor(rootcolors[plotStyle]); - histogrampointer->SetStats(kTRUE); - histogrampointer->SetStats(111111111); - histogrampointer->SetLineWidth(3); - histogrampointer->GetXaxis()->SetTitle("Q_coll [ADU]"); - histogrampointer->GetYaxis()->SetTitle(Form("Entries [1/%.1f ADU]",histogrampointer->GetBinWidth(1))); - histogrampointer->GetXaxis()->CenterTitle(); - histogrampointer->GetYaxis()->CenterTitle(); -} - void Run::initClusters(Float_t xcoord_min_step, Float_t xcoord_abs_min, Float_t xcoord_abs_max, Float_t ycoord_min_step, Float_t ycoord_abs_min, Float_t ycoord_abs_max) { if (false) // set to true to debug @@ -2099,22 +1614,11 @@ void Run::initClusters(Float_t xcoord_min_step, Float_t xcoord_abs_min, Float_t cout << "sizearrrotx: " << sizearrrotx << " sizearrroty: " << sizearrroty << endl; } - initCluster(&histogram, "", xcoord_min_step, xcoord_abs_min, xcoord_abs_max, ycoord_min_step, ycoord_abs_min, ycoord_abs_max); - initCluster(&histogramthreshold, "_threshold", xcoord_min_step, xcoord_abs_min, xcoord_abs_max, ycoord_min_step, ycoord_abs_min, ycoord_abs_max); - initCluster(&histogramfixedthreshold, "_fixed_threshold", xcoord_min_step, xcoord_abs_min, xcoord_abs_max, ycoord_min_step, ycoord_abs_min, ycoord_abs_max); + for (vector::iterator curHistogramClass = HistogramClassVector.begin(); curHistogramClass != HistogramClassVector.end(); curHistogramClass++) { + (*curHistogramClass)->initCluster((*curHistogramClass)->histogramdescription, xcoord_min_step, xcoord_abs_min, xcoord_abs_max, ycoord_min_step, ycoord_abs_min, ycoord_abs_max); + } } -void Run::initCluster(histogramstruct* histogramstructpointer, TString suffix, Float_t xcoord_min_step, Float_t xcoord_abs_min, Float_t xcoord_abs_max, Float_t ycoord_min_step, Float_t ycoord_abs_min, Float_t ycoord_abs_max) -{ - TString prefix = "Cluster" + suffix; - Int_t sizearrrotx= (int)((xcoord_abs_max-xcoord_abs_min)/xcoord_min_step+0.5); - Int_t sizearrroty = (int)((ycoord_abs_max-ycoord_abs_min)/ycoord_min_step+0.5); -// cout << "xcoord_abs_min: " << xcoord_abs_min << ", xcoord_abs_max: " << xcoord_abs_max << ", xcoord_min_step: " << xcoord_min_step << endl; -// cout << "ycoord_abs_min: " << ycoord_abs_min << ", ycoord_abs_max: " << ycoord_abs_max << ", ycoord_min_step: " << ycoord_min_step << endl; - histogramstructpointer->histAvgCluster = new TH2F(prefix.Data(), Form("Avg. cluster distribution; x; y %s", Form("%s, %s", prefix.Data(), humanreadablestr.Data())),sizearrrotx, xcoord_abs_min,xcoord_abs_max,sizearrroty,ycoord_abs_min,ycoord_abs_max); - // you can increase the resolution, especially usefull if you observe an rotated (originally staggered) pixel cluster - // histogramstructpointer->histAvgCluster = new TH2F(prefix.Data(), Form("Avg. cluster distribution; x; y %s", humanreadablestr.Data()),2*sizearrrotx, xcoord_abs_min,xcoord_abs_max,sizearrroty*2,ycoord_abs_min,ycoord_abs_max); -} void Run::initRootParameters() { diff --git a/MABS_run_analyzer/Run.h b/MABS_run_analyzer/Run.h index 226286c..d2f8a3b 100644 --- a/MABS_run_analyzer/Run.h +++ b/MABS_run_analyzer/Run.h @@ -36,6 +36,7 @@ #define SERVERPWD "mimosa88" #include "MAPS.c" +#include "HistogramType.c" class MAPS; /** @@ -86,20 +87,7 @@ private: * @brief takes a float value and a precision and appends it to an SQL statement used to update the database */ void constructUpdateString(string*, const string, const Double_t, int); - - /** - * @brief A structure to hold information about a given sensor, like number of #columns and #rows - */ - struct sensorinfostruct - { - /// number of rows read out in the system the run was taken with (USB or PXI) - Int_t rows; - /// number of columns read out in the system the run was taken with (USB or PXI) - Int_t columns; - sensorinfostruct():rows(0),columns(0){}; - sensorinfostruct( Int_t r, Int_t c ) : rows( r ), columns( c ) {} - }; - + /** * @brief writes values back to the SQL database */ @@ -109,11 +97,7 @@ private: * @brief removes all numbers from given string */ string removeNumbers(std::string x); - - /** - * @brief Lorenzian Peak function */ - static Double_t lorentzianPeak(Double_t *x, Double_t *par); - + /** * @brief removes all punctuation from given string */ @@ -161,7 +145,7 @@ private: /** * @brief fills noise histograms in #Run::histogramstruct */ - Bool_t binNoise(); + Bool_t binNoise(HistogramType* oneHistogramClass); /** * @brief fills Seedm, Sum and Veto #Run::histogramstruct */ Bool_t binSeedSumVeto(); @@ -182,23 +166,7 @@ private: * */ void plotVerticalLine(TH1F* histogrampointer, Float_t xVal); - - /** - * @brief rescales all histograms from ADU to electrons */ - Bool_t rescaleHistograms(); - - /** - * @brief rescales one specific histogram from ADU to electrons */ - void rescaleHistogram(TH1F* &histogrampointernew, TH1F* &histogrampointerold ); - - /** - * @brief rescales one specific histogram bin content from ADU to electrons */ - void rescale2DHistogramCounts(TH2F* &histogrampointernew, TH2F* &histogrampointerold ); - - /** - * @brief calculates Charge Collection Efficiency if Fe55 run */ - Bool_t calculteCCE(); - + /** * @brief Checks if a file exists */ @@ -218,50 +186,9 @@ private: }; pixelinfo curpixelinfo; - /** @brief system specific constants - */ - struct systemparam - { - Int_t maxbin; - Int_t nbins; - Int_t vetothreshold; - Int_t maxbinnoise; - Int_t nbinsnoise; - }; - systemparam systemparamUSB { - 2800, // maxbin; - 2800/4,// nbins; - 25, //vetothreshold - 10, - 100 - }; - systemparam systemparamPegasus { - 2800, // maxbin; - 2800/2,// nbins; - 20, //vetothreshold - 10, - 100 - }; - systemparam systemparamPXI { - 800*16, // maxbin - 800,// nbins; - 75, //vetothreshold - 150, - 150 - }; - systemparam systemparamFSBB { - 2800, // maxbin; - 2800/4,// nbins; - 25, //vetothreshold - 10, - 100 - }; Int_t* rootcolors; Int_t* rootlinestyle; - - /// init one specific histogram - void initHistogram(TH1F* &histogrampointer, TString prefix); - + /** @brief calls #initCluster for each structure of type @c Run::histogram * * For the initialization parameters calculated by #findRotatedClusterDimension() are used. @@ -284,13 +211,7 @@ private: /** @brief set this variable to true, to use the common mode noise filter in MAPS:analyzeRun() */ Bool_t commonModeFilter = 1; - const TString colorwhite = "\033[1;29m"; - const TString colorred = "\033[1;31m"; - const TString coloryellow = "\033[1;33m"; - const TString colorgreen = "\033[1;32m"; - const TString colorcyan = "\033[1;36m"; - const TString colorreset = "\033[0m"; - const TString endlr = "\033[0m\n"; + Bool_t rescaleHistogramClasses(); public: /** @brief empty constructor */ @@ -302,18 +223,18 @@ public: /** @brief writes values back to the database, destroys the #MAPS object */ ~Run (void); - systemparam systemparamcur; + systemparam cursystemparam; void setPlotStyle(Int_t); /** @brief Compare histograms * * Compare histograms given in the vector - * #Run::plothistogramstructpointerarray which holds structures + * #Run::plothistogramtypepointerarray which holds structures * of type #Run::histogram into one canvas * * * Use a line such as - * @code runs[runi]->plothistogramstructpointerarray.push_back(&runs[runi]->histogramfixedthreshold); + * @code runs[runi]->plothistogramtypepointerarray.push_back(&runs[runi]->histogramfixedthreshold); * @endcode * to add histograms to compare, after adding a few, execute * @code runs[runi]->plotCompareHistograms(); @@ -322,22 +243,7 @@ public: * */ Bool_t plotCompareHistograms( ); - - /** @brief Plot all histograms from @c Run::histogram into one canvas */ - Bool_t plotAllHistograms(); - - /** @brief Plot all histograms from @c Run::histogram into one canvas */ - Bool_t plotAllHistogramsCalibrated(); - - /** @brief Plot all histograms from @c Run::histogramthreshold into one canvas */ - Bool_t plotAllHistogramsThresholdCluster(); - - /** @brief Plot all histograms from @c Run::histogramthreshold into one canvas */ - Bool_t plotAllHistogramsThresholdClusterCalibrated(); - - /** @brief Plot average cluster distribution from @c Run::histogram into one canvas */ - Bool_t plotClusterDistribution(); - + /** @brief Writes a given 2D-histogram into a file */ Bool_t writeHistogramToFile(TH1F* onehistogram); @@ -378,14 +284,7 @@ public: Int_t getClustersize(); Bool_t runAllreadyAnalyzed(); - - Bool_t plotNoise(); - Bool_t plotSeed(); - Bool_t plotSeedThreshold(); - Bool_t plotSeedThresholdCalibrated(); - Bool_t plotSum(); - Bool_t plotVeto(); - + /** @brief Prints all data gathered from the database */ Bool_t debugDBreadout(); @@ -490,67 +389,12 @@ public: /** @brief calls #Run::setSubMatrixBorders() with the correct boundaries for FSBB submatrix A0 or A1 */ void selectSubMatrixFSBB(TString matrixname); - - /** @brief A structure to hold the histogram data for the run - * - * seed, sum and veto histograms are stored here - * - */ - struct histogramstruct - { - /// Seed spectrum, only the seed pixel is considered when this histogram is build - TH1F* Seed; - /// Sum spectrum, the charge of whole cluster is summed up in binned into this TH1F histogram - TH1F* Sum; - /// Veto spectrum, used to find better calibration peak, only entries where Sum over not seed pixels is below @c Run::vetothreshold are binned into this histogram - TH1F* Veto; - /// Noise histogram - TH1F* Noise; - - /// fitted position of the most probable value of the seed spectrum - Float_t posSeed=0; - /// fitted position of the most probable value in the over cluster summed spectrum - Float_t posSum=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; - - /// Average Noise value - Float_t avgNoise = 0; - /// Positive Noise sigma - Float_t avgNoisePlus = 0; - /// Negative Noise sigma - Float_t avgNoiseMinus = 0; - - - /// threshold for integrating the Sr90 spectrum, is set dynamically in @c integrateSr90Spectra() - Float_t noisethresholdborder = -1; - /// Integral value, after integrating from #noisethresholdborder to maxbin. - Double_t sr90IntegralVal = -1; - Double_t sr90IntegralErr = -1; - - /// set to true, if bins are in electrons, otherwise in ADU - Bool_t calibrated = false; - - /// set to true, if threshold clusters are used - Bool_t thresholdcluster = false; - - /// histogram with cluster data, rotated if staggered - TH2F* histAvgCluster; - - /// number of hits/clusters used to generate all distributions - Double_t numberofhits = 0; - - /// type in here what the histogram is intended for or how it is calculated, will be added to filenames - TString histogramdescription = ""; - }; - histogramstruct histogram; - histogramstruct histogramCalibrated; - histogramstruct histogramthreshold; - histogramstruct histogramthresholdCalibrated; - histogramstruct histogramfixedthreshold; - histogramstruct histogramfixedthresholdCalibrated; - histogramstruct* plothistogramstructpointer; - TH1F** plothistogrampointer; + + HistogramType* histogram; + HistogramType* histogramthreshold; + HistogramType* histogramfixedthreshold; + HistogramType* plothistogramclasspointer; + TH1F* plothistogrampointer; /// /** @brief A vector able to hold pointer of type #Run::histogramstruct * @@ -560,59 +404,30 @@ public: * See Run::plotCompareHistograms( ) for more information and usage. * */ - vector plothistogramstructpointerarray; + vector compareHistogramClassVector; + + + // holds all used HistogramType classes + vector HistogramClassVector; + - /** @brief Plots all histograms from #Run::histogramstruct into one canvas + /** @brief Plots all histograms from #Run::HistogramType into one canvas * - * @param histogramstruct pointer to a histogram structure of type #Run::histogramstruct + * @param HistogramType pointer to a histogram structure of type #Run::histogramstruct */ - Bool_t plotAllHistograms(histogramstruct*); + Bool_t plotAllHistograms(HistogramType*); /** @brief Plots average cluster charge distribution * * Creates a graphical view of the average charge distribution in a * cluster * - * @param histogramstruct Pointer to a structure of type #Run::histogramstruct + * @param HistogramType Pointer to a structure of type #Run::histogramstruct * */ - Bool_t plotClusterDistribution(histogramstruct*); + Bool_t plotClusterDistribution(HistogramType*); - /** - * @brief intern function to calculate and plot fit curve to given histogram - * - * @param histogrampointer histogram pointer to calculate fit to - * @return peak position of the fit - * - */ - Float_t FitPerform(histogramstruct*, TH1F*, TString fitFuncType="landau", Bool_t verbose=true, Double_t* parameters=NULL); - - - TCanvas* plot1DHistogram(histogramstruct* histogramstructpointer, TH1F* onehistogram, TString fitFuncType = "landau", TString titlestr = "", TString legendstr = ""); - - /// - /** @brief initializes all histograms, set binwidth, bin amount and names - * * - * @param histogramstruct Pointer to a structure of type #Run::histogramstruct - * @param suffix A string which will be appended to the generated title of each histograms - */ - void initHistograms(histogramstruct*, TString suffix = ""); - - /** @brief initializes the TH2F cluster for the binning for a specific structure of type #Run::histogramstruct one points to*/ - void initCluster(histogramstruct*, TString suffix, Float_t xcoord_min_step, Float_t xcoord_abs_min, Float_t xcoord_abs_max, Float_t ycoord_min_step, Float_t ycoord_abs_min, Float_t ycoord_abs_max); - - - /** - * @brief find the border between the noise and the signal in Sr90 runs - * - * writes the found threshold into the private variable @c Run::posNoiseThreshold - * - * @param histogrampointer pointer to the histogram structure, threshold will be searched in seed spectra - * @param thresholdborder value from which the integral will be calculated, if not set, the algorithms tries to find the best value - * @return true if succesfull - * - */ - Bool_t integrateSr90Spectra(histogramstruct* histogramstructpointer, TH1F* histogrampointer, Float_t thresholdborder=-1, Bool_t verbose=true); + TCanvas* plot1DHistogram(HistogramType* histogramtypepointer, TH1F* onehistogram, TString fitFuncType = "landau", TString titlestr = "", TString legendstr = ""); pixelinfo pixelinfoMi34[32]; @@ -651,18 +466,7 @@ public: * Fe55 calibration peak. ADU can be converted to e- */ OtherRun Fe55run; - /// to convert from ADU to values in electron one has to divide by this value, it is given in [1/e] - Float_t gain=0; - /// Charge collection efficciency of the cluster in percent - Float_t CCE_in_Perc_25=0; - /// Charge collection efficciency of the seed pixel in percent - Float_t CCE_in_Perc_1=0; - /// threshold for integrating the Sr90 spectrum, is set dynamically in @c integrateSr90Spectra() - Float_t noisethresholdborder = -1; - /// Integral value, after integrating from #noisethresholdborder to maxbin. - Double_t sr90IntegralVal = -1; - Double_t sr90IntegralErr = -1; - + /** @brief specific string wich encodes the database information of the run */ TString runcode=""; @@ -679,7 +483,7 @@ public: MAPS *processed; /// sensor information to use in analysis, is the system read out by USB or PXI? Number of rows differ - sensorinfostruct sensorinfocur; + sensorinfostruct cursensorinfo; }; diff --git a/MABS_run_analyzer/help.h b/MABS_run_analyzer/help.h index 4ab1196..16dc6a9 100644 --- a/MABS_run_analyzer/help.h +++ b/MABS_run_analyzer/help.h @@ -117,6 +117,49 @@ struct frameInfo{ Float_t p [25][MAXHITS]; }; +/** @brief system specific constants + * + * Set in #HistogramType::setSystemSpecificParameters() + */ +struct systemparam { + /// last bin value to use in a histogram + Int_t maxbin; + /// number of bins to use in a histogram + Int_t nbins; + /// threshold to use when binning into the veto spectrum hisogram + Int_t vetothreshold; + /// maximum bin value for the noise histogram + Int_t maxbinnoise; + /// number of bins to use in the noise histogram + Int_t nbinsnoise; + /// default setter + systemparam():maxbin(0),nbins(0),vetothreshold(0),maxbinnoise(0),nbinsnoise(0){}; + /// setter + systemparam( Int_t mb, Int_t nb, Int_t vt, Int_t mbn, Int_t nbn ) : + maxbin(mb), nbins(nb), vetothreshold(vt), maxbinnoise(mbn), nbinsnoise(nbn){} +}; + + +/** + * @brief A structure to hold information about a given sensor, like number of #columns and #rows + */ +struct sensorinfostruct { + /// number of rows read out in the system the run was taken with (USB or PXI) + Int_t rows; + /// number of columns read out in the system the run was taken with (USB or PXI) + Int_t columns; + sensorinfostruct():rows(0),columns(0){}; + sensorinfostruct( Int_t r, Int_t c ) : rows( r ), columns( c ) {} +}; +/// sensor information to use in analysis, is the system read out by USB or PXI? Number of rows differ + +const TString colorwhite = "\033[1;29m"; +const TString colorred = "\033[1;31m"; +const TString coloryellow = "\033[1;33m"; +const TString colorgreen = "\033[1;32m"; +const TString colorcyan = "\033[1;36m"; +const TString colorreset = "\033[0m"; +const TString endlr = "\033[0m\n"; void preparecanvas() { // TCanvas* canvas = new TCanvas("Canvas1", 700, 500); -- 2.43.0