From: Benjamin Linnik Date: Tue, 15 Sep 2015 12:28:36 +0000 (+0200) Subject: Run analyzer: Added option to integrate histograms from a certain threshold X-Git-Url: https://jspc29.x-matter.uni-frankfurt.de/git/?a=commitdiff_plain;h=2dc27f1d0189b0b70870a25abbfe852c968497c7;p=radhard.git Run analyzer: Added option to integrate histograms from a certain threshold --- diff --git a/MABS_run_analyzer/MAPS.c b/MABS_run_analyzer/MAPS.c index 0757651..ca788a0 100644 --- a/MABS_run_analyzer/MAPS.c +++ b/MABS_run_analyzer/MAPS.c @@ -1644,7 +1644,7 @@ void MAPS::hitana() { else fFrameInfo.pixelthreshold[fHits] = 0; - if (1.0*chargesumincluster4 > 150) + if (1.0*chargesumincluster4 > fFixedThresholdValueInADU) fFrameInfo.pixelfixedthreshold[fHits] = Hitlist[hit]; else fFrameInfo.pixelfixedthreshold[fHits] = 0; @@ -1656,8 +1656,10 @@ void MAPS::hitana() { else fFrameInfo.pixelthreshold[fHits] = 0; - if (1.0*chargesumincluster > 150) + if (1.0*chargesumincluster > fFixedThresholdValueInADU) + { fFrameInfo.pixelfixedthreshold[fHits] = Hitlist[hit]; + } else fFrameInfo.pixelfixedthreshold[fHits] = 0; } diff --git a/MABS_run_analyzer/MAPS.h b/MABS_run_analyzer/MAPS.h index 8a372de..67a46ac 100644 --- a/MABS_run_analyzer/MAPS.h +++ b/MABS_run_analyzer/MAPS.h @@ -76,7 +76,7 @@ private: /// number of columns in the sensor, set and passed initMapsRun() to in the constructorm, if Mimosa34 then it is equal to 64 Int_t fColumns; /// Variable to be able to subdivide the Matrix according to fOrderCode - TString fMatrix; + TString fMatrix; UInt_t fEventsSum; Int_t fFile; @@ -445,6 +445,9 @@ public: * @param Noise Float number to set the pedestial of each pixel */ bool setPedestals (Float_t Pedestals); + /// Value to use when filling up #MAPS::fFrameInfo.pixelfixedthreshold, set in Run::setFixedThresholdValueElectrons() and Run::setFixedThresholdValueADU() + Float_t fFixedThresholdValueInADU = 150.0; + /// stores infomration about the current frame, like cluster frameInfo fFrameInfo; // pointer to TTree with hit, frame and cluster information diff --git a/MABS_run_analyzer/Run.c b/MABS_run_analyzer/Run.c index c8a086a..35c5c44 100644 --- a/MABS_run_analyzer/Run.c +++ b/MABS_run_analyzer/Run.c @@ -88,15 +88,15 @@ Run::Run(Int_t runnumber, Int_t loopi) storepathRAWLinux=Form("/d/garlic/%s/%d",labbook.chipGen.Data(),runnumber); // default empty path } savepathresults = storepathOutLinux; - labbook.posVetoDB = (rowsql->GetField(10) != NULL)?atoi(rowsql->GetField(10)):-1; - labbook.posSeedDB = (rowsql->GetField(11) != NULL)?atoi(rowsql->GetField(11)):-1; - labbook.posSumDB = (rowsql->GetField(12) != NULL)?atoi(rowsql->GetField(12)):-1; - labbook.gainDB = (rowsql->GetField(13) != NULL)?atoi(rowsql->GetField(13)):-1; - labbook.NoiseAvgDB = (rowsql->GetField(14) != NULL)?atoi(rowsql->GetField(14)):-1; - labbook.NoiseAvgPlusDB = (rowsql->GetField(15) != NULL)?atoi(rowsql->GetField(15)):-1; - labbook.NoiseAvgMinusDB = (rowsql->GetField(16) != NULL)?atoi(rowsql->GetField(16)):-1; - labbook.CCE_in_Perc_1DB = (rowsql->GetField(17) != NULL)?atoi(rowsql->GetField(17)):-1; - labbook.CCE_in_Perc_25DB = (rowsql->GetField(18) != NULL)?atoi(rowsql->GetField(18)):-1; + labbook.posVetoDB = (rowsql->GetField(10) != NULL)?atof(rowsql->GetField(10)):-1; + labbook.posSeedDB = (rowsql->GetField(11) != NULL)?atof(rowsql->GetField(11)):-1; + labbook.posSumDB = (rowsql->GetField(12) != NULL)?atof(rowsql->GetField(12)):-1; + labbook.gainDB = (rowsql->GetField(13) != NULL)?atof(rowsql->GetField(13)):-1; + labbook.NoiseAvgDB = (rowsql->GetField(14) != NULL)?atof(rowsql->GetField(14)):-1; + labbook.NoiseAvgPlusDB = (rowsql->GetField(15) != NULL)?atof(rowsql->GetField(15)):-1; + labbook.NoiseAvgMinusDB = (rowsql->GetField(16) != NULL)?atof(rowsql->GetField(16)):-1; + labbook.CCE_in_Perc_1DB = (rowsql->GetField(17) != NULL)?atof(rowsql->GetField(17)):-1; + labbook.CCE_in_Perc_25DB = (rowsql->GetField(18) != NULL)?atof(rowsql->GetField(18)):-1; labbook.frames_foundDB = (rowsql->GetField(19) != NULL)?atoi(rowsql->GetField(19)):-1; // labbook.frames_Analyzed = (rowsql->GetField(20) != NULL)?atoi(rowsql->GetField(20)):-1; delete res; @@ -196,6 +196,7 @@ Bool_t Run::debugDBreadout() cout << "| NoiseAvgMinusDB: " << std::right << colorwhite << labbook.NoiseAvgMinusDB << endlr; cout << "| CCE_in_Perc_25DB:" << std::right << colorwhite << labbook.CCE_in_Perc_25DB << endlr; cout << "| CCE_in_Perc_1DB: " << std::right << colorwhite << labbook.CCE_in_Perc_1DB << endlr; + cout << "| frames_foundDB : " << std::right << colorwhite << labbook.frames_foundDB << endlr; cout << "|______________________________ " << endlr; cout << endlr; return 0; @@ -215,6 +216,27 @@ Bool_t Run::useCommonModeFilter(Bool_t var) return 0; } +Bool_t Run::setFixedThresholdValueElectrons(Float_t thresholdValueIne) +{ + if (labbook.gainDB > 0) { + cout << " Fixed threshold value : " << colorwhite << thresholdValueIne << " e" << colorreset << " <-- only used if RAW files are analyzed, force analysis to make sure" << endl; + processed->fFixedThresholdValueInADU = thresholdValueIne / labbook.gainDB; + return 0; + } + else + { + cout << colorred << "Could not set value for fixed threshold to " << thresholdValueIne << ", no calibration done yet. Please rerun after first analysis." <fFixedThresholdValueInADU = thresholdValue; + return 0; +} + Bool_t Run::setResultsPath(TString path) { cout << "Changing save directory to '" << path << "'" << endl; @@ -322,7 +344,7 @@ Bool_t Run::analyzeRun(Bool_t force) calculteCCE(); if (labbook.source.Contains("Sr90")) { cout << colorwhite << "integrateSr90Spectra():" << endlr; - integrateSr90Spectra(histogramthresholdCalibrated.Seed); + integrateSr90Spectra(&histogramthresholdCalibrated, histogramthresholdCalibrated.Seed); } cout << colorwhite << "updateDatabase():" << endlr; updateDatabase(); @@ -738,7 +760,7 @@ void Run::updateDatabase() { constructUpdateString(&sqlupdatequery, "CCE_25", CCE_in_Perc_25); constructUpdateString(&sqlupdatequery, "Avg.NoiseADC", histogram.avgNoise); constructUpdateString(&sqlupdatequery, "Frames_found", frames_found,100000000); - constructUpdateString(&sqlupdatequery, "Sr90IntegralVal", sr90IntegralVal,1000000000); + constructUpdateString(&sqlupdatequery, "Sr90IntegralVal", histogramthresholdCalibrated.sr90IntegralVal,1000000000); if (sqlupdatequery.length()>0) { @@ -886,17 +908,17 @@ Bool_t Run::binSeedSumVeto() } // gROOT->SetBatch(kTRUE); if (labbook.source.Contains("Fe")||labbook.source.Contains("Cd")) - histogram.posVeto=FitPerform(histogram.Veto, "gaus", true); - histogram.posSeed=FitPerform(histogram.Seed, "landau", true); - histogram.posSum=FitPerform(histogram.Sum, "gaus", true); + 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.Veto, "gaus", false); - histogramthreshold.posSeed=FitPerform(histogram.Seed, "landau", false); - histogramthreshold.posSum=FitPerform(histogram.Sum, "gaus", false); + 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.Veto, "gaus", false); - histogramfixedthreshold.posSeed=FitPerform(histogramfixedthreshold.Seed, "landau", false); - histogramfixedthreshold.posSum=FitPerform(histogramfixedthreshold.Sum, "gaus", false); + histogramfixedthreshold.posVeto=FitPerform(&histogramfixedthreshold, histogramfixedthreshold.Veto, "gaus", false); + histogramfixedthreshold.posSeed=FitPerform(&histogramfixedthreshold, histogramfixedthreshold.Seed, "landau", false); + histogramfixedthreshold.posSum=FitPerform(&histogramfixedthreshold, histogramfixedthreshold.Sum, "gaus", false); // gROOT->SetBatch(kFALSE); return 0; } @@ -1131,7 +1153,7 @@ Bool_t Run::plotNoise() if (!error) { TString legendEntry = Form("Noise: %.2f + %.2f - %.2f", histogram.avgNoise, histogram.avgNoisePlus, histogram.avgNoiseMinus ); - noisecanvas = plot1DHistogram(histogram.Noise, "landau", "", legendEntry); + noisecanvas = plot1DHistogram(&histogram, histogram.Noise, "landau", "", legendEntry); return 0; } return 1; @@ -1141,7 +1163,7 @@ Bool_t Run::plotSeed() { if (!error) { - plot1DHistogram(histogram.Seed, "landau"); + plot1DHistogram(&histogram, histogram.Seed, "landau"); return 0; } return 1; @@ -1151,7 +1173,7 @@ Bool_t Run::plotSeedThreshold() { if (!error) { - plot1DHistogram(histogramthreshold.Seed, "landau"); + plot1DHistogram(&histogramthreshold, histogramthreshold.Seed, "landau"); return 0; } return 1; @@ -1161,7 +1183,7 @@ Bool_t Run::plotSeedThresholdCalibrated() { if (!error) { - plot1DHistogram(histogramthresholdCalibrated.Seed, "landau"); + plot1DHistogram(&histogramthresholdCalibrated, histogramthresholdCalibrated.Seed, "landau"); return 0; } return 1; @@ -1171,7 +1193,7 @@ Bool_t Run::plotSum() { if (!error) { - plot1DHistogram(histogram.Sum, "gaus"); + plot1DHistogram(&histogram, histogram.Sum, "gaus"); return 0; } return 1; @@ -1181,7 +1203,7 @@ Bool_t Run::plotVeto() { if (!error) { - plot1DHistogram(histogram.Veto, "gaus"); + plot1DHistogram(&histogram, histogram.Veto, "gaus"); return 0; } return 1; @@ -1460,13 +1482,13 @@ Bool_t Run::plotClusterDistribution(histogramstruct* histogramstructpointer) canvas->cd(1); TPad* firstpart = (TPad*)canvas->GetPad(1); firstpart->SetGrid(); - FitPerform(xslicetroughcluster,"lorentz", false, &xslicetroughclusterpar[0]); + FitPerform(&histogram, xslicetroughcluster,"lorentz", false, &xslicetroughclusterpar[0]); canvas->cd(2); TPad* secondpart = (TPad*)canvas->GetPad(2); secondpart->SetGrid(); - FitPerform(yslicetroughcluster,"lorentz", false, &yslicetroughclusterpar[0]); + FitPerform(&histogram, yslicetroughcluster,"lorentz", false, &yslicetroughclusterpar[0]); canvas->Draw(); @@ -1534,32 +1556,45 @@ Bool_t Run::plotAllHistogramsThresholdClusterCalibrated() return 1; } -Bool_t Run::integrateSr90Spectra(TH1F* histogrampointer, Bool_t verbose) +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 rising = 0; - Int_t bini =smoothedcurce->GetMaximumBin(); - Float_t curval = smoothedcurce->GetXaxis()->GetBinCenter(bini); - Float_t thresholdbincurcandidate = 0; +// 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(); - 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()); - + 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)); @@ -1567,9 +1602,10 @@ Bool_t Run::integrateSr90Spectra(TH1F* histogrampointer, Bool_t verbose) // histogrampointer->GetXaxis()->SetRange(posNoiseMax,histogrampointer->GetXaxis()->FindBin(posMaxValHist)); // Float_t posChargeMax= histogrampointer->GetMaximum(); - noisethresholdborder = histogrampointer->GetXaxis()->GetBinUpEdge(thresholdbincurcandidate); + 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]")) @@ -1587,25 +1623,32 @@ Bool_t Run::integrateSr90Spectra(TH1F* histogrampointer, Bool_t verbose) } } } - sr90IntegralVal = histogrampointer->IntegralAndError(thresholdbincurcandidate,histogrampointer->GetMaximumBin(), sr90IntegralErr); - sr90IntegralErr /= sr90IntegralVal/100; + 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) { - sr90IntegralVal/=(labbook.frames_foundDB>0)?labbook.frames_foundDB:frames_found; + histogramstructpointer->sr90IntegralVal/=(labbook.frames_foundDB>0)?labbook.frames_foundDB:frames_found; cout << "Scaled "; } - cout << "Integral is " << Form("%e",sr90IntegralVal) << ", error " << sr90IntegralErr << "%" << endl; + cout << "Integral is " << Form("%e",histogramstructpointer->sr90IntegralVal) << ", error " << histogramstructpointer->sr90IntegralErr << "%" << endl; return 0; } -Float_t Run::FitPerform(TH1F* histogrampointer, TString fitFuncType, Bool_t verbose, Double_t* parameters) +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(); - Float_t noiseborder = (noisethresholdborder>0)?noisethresholdborder:90; // posMaxValHist/10 for USB system, the value is 90 + // 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) { @@ -1800,7 +1843,7 @@ void Run::plotVerticalLine(TH1F* histogrampointer, Float_t xVal) { } } -TCanvas* Run::plot1DHistogram(TH1F* onehistogram, TString fitFuncType, TString titlestr, TString legendstr) +TCanvas* Run::plot1DHistogram(histogramstruct* histogramstructpointer, TH1F* onehistogram, TString fitFuncType, TString titlestr, TString legendstr) { if (onehistogram != nullptr) { @@ -1812,7 +1855,7 @@ TCanvas* Run::plot1DHistogram(TH1F* onehistogram, TString fitFuncType, TString t TCanvas* canvas = new TCanvas(canvasname, canvastitle, 900, 700); onehistogram->SetTitle(titlestr); onehistogram->Draw(); - Float_t maxValue = FitPerform(onehistogram, fitFuncType); + Float_t maxValue = FitPerform(histogramstructpointer, 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); diff --git a/MABS_run_analyzer/Run.h b/MABS_run_analyzer/Run.h index b0b35c6..f396971 100644 --- a/MABS_run_analyzer/Run.h +++ b/MABS_run_analyzer/Run.h @@ -171,7 +171,6 @@ private: /// noise quantiles: mean value, sigma in postive direction and sigma in negative direction Double_t noisequantiles[3]; - TCanvas* plot1DHistogram(TH1F* onehistogram, TString fitFuncType = "landau", TString titlestr = "", TString legendstr = ""); /** @@ -182,28 +181,7 @@ private: * */ void plotVerticalLine(TH1F* histogrampointer, Float_t xVal); - - /** - * @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*, TString fitFuncType="landau", Bool_t verbose=true, Double_t* parameters=NULL); - - /** - * @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 - * @return true if succesfull - * - */ - Bool_t integrateSr90Spectra(TH1F* histogrampointer, Bool_t verbose=true); - - + /** * @brief rescales all histograms from ADU to electrons */ Bool_t rescaleHistograms(); @@ -370,7 +348,7 @@ public: /** @brief Makes a gnuplot file to plot the histogram data created in @c Run::writeAllHistogramsToFile() */ void MakeGnuplotFile(); - + /// is set to true if an error occured Bool_t error = 0; @@ -396,7 +374,6 @@ public: */ Bool_t analyzeFrame(Int_t frame); - Int_t getClustersize(); Bool_t runAllreadyAnalyzed(); @@ -416,6 +393,12 @@ public: /** @brief Turn on or off the use of common mode noise filter @c MAPS::regetDynNoise() */ Bool_t useCommonModeFilter(Bool_t); + + /** @brief set a threshold value in electrons to use in the histogram Run::histogramfixedthreshold, this value is used in @c MAPS::initMapsRun() */ + Bool_t setFixedThresholdValueElectrons(Float_t); + + /** @brief set a threshold value in ADU to use in the histogram Run::histogramfixedthreshold, this value is used in @c MAPS::initMapsRun() */ + Bool_t setFixedThresholdValueADU(Float_t); /** * @brief set pathwhere images and data to save to @@ -537,6 +520,13 @@ public: /// 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; @@ -587,6 +577,18 @@ public: */ Bool_t plotClusterDistribution(histogramstruct*); + /** + * @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 * * @@ -598,6 +600,19 @@ public: /** @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); + pixelinfo pixelinfoMi34[32]; /// analyze only half of matrix