From 8a2aeaf200e58ccd8d663435766f0285724883f2 Mon Sep 17 00:00:00 2001 From: Benjamin Linnik Date: Tue, 8 Dec 2015 14:35:18 +0100 Subject: [PATCH] Run analyzer: Added possibility to define noise threshold cut manually and overwrite therefore the automatic threshold finding, fine tuning of automatic threshold finding mechanism - finds cut if noise and signal peak are more close to each other in Sr90 spectra --- MABS_run_analyzer/ChargeSpektrum.c | 1 + MABS_run_analyzer/HistogramType.c | 142 +++++++++++++++-------------- MABS_run_analyzer/HistogramType.h | 2 +- MABS_run_analyzer/Run.c | 38 ++++++-- MABS_run_analyzer/Run.h | 16 ++++ 5 files changed, 122 insertions(+), 77 deletions(-) diff --git a/MABS_run_analyzer/ChargeSpektrum.c b/MABS_run_analyzer/ChargeSpektrum.c index 5b470d0..40c1d83 100644 --- a/MABS_run_analyzer/ChargeSpektrum.c +++ b/MABS_run_analyzer/ChargeSpektrum.c @@ -132,6 +132,7 @@ void ChargeSpektrum(TString runnumber = "") runs[runi]->useDynamicalNoise(true); runs[runi]->useCommonModeFilter(false); runs[runi]->setFixedThresholdValueElectrons(1440); + runs[runi]->setNoisethresholdborderADU(40); // runs[runi]->analyzeFrame(57826); // runs[runi]->analyzeFrame(57827); diff --git a/MABS_run_analyzer/HistogramType.c b/MABS_run_analyzer/HistogramType.c index 78a8f70..148b1ab 100644 --- a/MABS_run_analyzer/HistogramType.c +++ b/MABS_run_analyzer/HistogramType.c @@ -381,78 +381,82 @@ Bool_t HistogramType::calculteCCE(Bool_t verbose) { return 1; } -Bool_t HistogramType::FindNoisethresholdborder(TH1F* histogrampointer, Bool_t verbose) { - Float_t posMaxValHist = histogrampointer->GetXaxis()->GetXmax(); - - TH1F* smoothedcurce = (TH1F*)histogrampointer->Clone(); - // smoothedcurce->Smooth(4); // not working method, rescaling is easier and more reliable - Int_t rebinningfactor = 4; - smoothedcurce->RebinX(rebinningfactor); - smoothedcurce->SetAxisRange(0,histogrampointer->GetBinCenter(histogrampointer->GetNbinsX())); - - if (verbose) { - cout << colorwhite << histogrampointer->GetName() << endlr; - TString canvastitle = Form("%s Noisethreshold border", histogrampointer->GetName()); - TString canvasname = Form("%s Noisethreshold border", histogrampointer->GetName()); - TCanvas* canvas = new TCanvas(canvasname, canvastitle, 900, 700); - smoothedcurce->SetLineColor(color+1); - smoothedcurce->Draw(); - histogrampointer->Draw("SAME"); - canvas->Update(); - } - - // Int_t bini =smoothedcurce->GetMaximumBin(); - // thresholdbincurcandidate = bini; - // if (verbose) - // cout << "GetMaximumBin: smoothedcurce->GetBinContent(" << bini << "): " << smoothedcurce->GetBinContent(bini) << endl; - Int_t bini = 0; - Int_t falling = 0; - Int_t noisepeakcandidate = 0; - - do { - Float_t curval=smoothedcurce->GetBinContent(bini++); - if (curval*0.95 <= smoothedcurce->GetBinContent(bini)) - { - falling = 0; - noisepeakcandidate = bini; - } - else - { - falling++; - if (verbose) - cout << "falling at " << smoothedcurce->GetXaxis()->GetBinCenter(bini) << " as " << curval << " < " << smoothedcurce->GetBinContent(bini) << endl; // debug +Bool_t HistogramType::FindNoisethresholdborder(TH1F* histogrampointer, Bool_t overwrite, Bool_t verbose) { + if (overwrite || noisethresholdborder == -1) + { + Float_t posMaxValHist = histogrampointer->GetXaxis()->GetXmax(); + + TH1F* smoothedcurce = (TH1F*)histogrampointer->Clone(); + // smoothedcurce->Smooth(4); // not working method, rescaling is easier and more reliable + Int_t rebinningfactor = 2; + smoothedcurce->RebinX(rebinningfactor); + smoothedcurce->SetAxisRange(0,histogrampointer->GetBinCenter(histogrampointer->GetNbinsX())); + + if (verbose) { + cout << colorwhite << histogrampointer->GetName() << endlr; + TString canvastitle = Form("%s Noisethreshold border", histogrampointer->GetName()); + TString canvasname = Form("%s Noisethreshold border", histogrampointer->GetName()); + TCanvas* canvas = new TCanvas(canvasname, canvastitle, 900, 700); + smoothedcurce->SetLineColor(color+1); + smoothedcurce->Draw(); + histogrampointer->Draw("SAME"); + canvas->Update(); } - } while (falling < 2 && biniGetNbinsX()); - if (verbose) { - cout << "Noise peak at: smoothedcurce->GetBinContent(" << noisepeakcandidate << "): " << smoothedcurce->GetBinContent(noisepeakcandidate) << endl; - } - Int_t thresholdbincurcandidate = 0; - bini = noisepeakcandidate; - - Int_t rising = 0; - do { - Float_t curval=smoothedcurce->GetBinContent(bini++); - if (curval*0.95 <= smoothedcurce->GetBinContent(bini)) - { - rising++; - if (verbose) - cout << "rising at bin " << bini << ":" << smoothedcurce->GetXaxis()->GetBinCenter(bini) << " as " << curval << " < " << smoothedcurce->GetBinContent(bini) << endl; // debug + + // Int_t bini =smoothedcurce->GetMaximumBin(); + // thresholdbincurcandidate = bini; + // if (verbose) + // cout << "GetMaximumBin: smoothedcurce->GetBinContent(" << bini << "): " << smoothedcurce->GetBinContent(bini) << endl; + Int_t bini = 0; + Int_t falling = 0; + Int_t noisepeakcandidate = 0; + + do { + Float_t curval=smoothedcurce->GetBinContent(bini++); + if (curval*0.95 <= smoothedcurce->GetBinContent(bini)) + { + falling = 0; + noisepeakcandidate = bini; + } + else + { + falling++; + if (verbose) + cout << "falling at " << smoothedcurce->GetXaxis()->GetBinCenter(bini) << " as " << curval << " < " << smoothedcurce->GetBinContent(bini) << endl; // debug + } + } while (falling < 2 && biniGetNbinsX()); + if (verbose) { + cout << "Noise peak at: smoothedcurce->GetBinContent(" << noisepeakcandidate << "): " << smoothedcurce->GetBinContent(noisepeakcandidate) << endl; } - else - { - rising = 0; - thresholdbincurcandidate = bini; - if (verbose) - cout << "falling at bin " << bini << ":" << smoothedcurce->GetXaxis()->GetBinCenter(bini) << " as " << curval << " < " << smoothedcurce->GetBinContent(bini) << endl; // debug + Int_t thresholdbincurcandidate = 0; + bini = noisepeakcandidate; + + Int_t rising = 0; + do { + Float_t curval=smoothedcurce->GetBinContent(bini++); + if (curval*0.95 <= smoothedcurce->GetBinContent(bini)) + { + rising++; + if (verbose) + cout << "rising at bin " << bini << ":" << smoothedcurce->GetXaxis()->GetBinCenter(bini) << " as " << curval << " < " << smoothedcurce->GetBinContent(bini) << endl; // debug + } + else + { + rising = 0; + thresholdbincurcandidate = bini; + if (verbose) + cout << "falling at bin " << bini << ":" << smoothedcurce->GetXaxis()->GetBinCenter(bini) << " as " << curval << " < " << smoothedcurce->GetBinContent(bini) << endl; // debug + } + } while (rising < 3 && biniGetNbinsX()); + thresholdbincurcandidate *= rebinningfactor; + + noisethresholdborder = histogrampointer->GetXaxis()->GetBinUpEdge(thresholdbincurcandidate); + if (verbose) { + cout << " Noise threshold at " << noisethresholdborder << endl; } - } while (rising < 3 && biniGetNbinsX()); - thresholdbincurcandidate *= rebinningfactor; - - noisethresholdborder = histogrampointer->GetXaxis()->GetBinUpEdge(thresholdbincurcandidate); - if (verbose) { - cout << " Noise threshold at " << noisethresholdborder << endl; + return 0; } - return 0; + return 1; } Bool_t HistogramType::integrateSr90Spectra(TH1F* histogrampointer, Int_t frames_found, Float_t thresholdborder, Bool_t verbose) { @@ -461,7 +465,7 @@ Bool_t HistogramType::integrateSr90Spectra(TH1F* histogrampointer, Int_t frames_ if (thresholdborder < 0) { if (noisethresholdborder < 0) { - FindNoisethresholdborder(histogrampointer, verbose); + FindNoisethresholdborder(histogrampointer, true, verbose); } thresholdbincurcandidate = histogrampointer->GetXaxis()->FindBin(noisethresholdborder); } diff --git a/MABS_run_analyzer/HistogramType.h b/MABS_run_analyzer/HistogramType.h index 76c7a17..43e310b 100644 --- a/MABS_run_analyzer/HistogramType.h +++ b/MABS_run_analyzer/HistogramType.h @@ -216,7 +216,7 @@ public: * @return true if succesfull * */ - Bool_t FindNoisethresholdborder(TH1F* histogrampointer, Bool_t verbose=false); + Bool_t FindNoisethresholdborder(TH1F* histogrampointer, Bool_t overwrite=false, Bool_t verbose=false); /** * @brief calculates StoN ratio, if beta source was used diff --git a/MABS_run_analyzer/Run.c b/MABS_run_analyzer/Run.c index f1d67cb..f9e121b 100644 --- a/MABS_run_analyzer/Run.c +++ b/MABS_run_analyzer/Run.c @@ -380,7 +380,7 @@ Bool_t Run::analyzeRun(Bool_t force) } else { - cout << "Skipped analysis of run " << labbook.runnumber << ", I think the root file for this run allready exists." << endl; + cout << endl <<"Skipped analysis of run " << labbook.runnumber << ", I think the root file for this run allready exists." << endl; cout << colorwhite << "initOldRootFile():" << endlr; } if (processed->initOldRootFile()) return 1; @@ -394,7 +394,16 @@ Bool_t Run::analyzeRun(Bool_t force) cout << "---------- loop over all classes --------" << endl; for (vector::iterator curHistogramClass = HistogramClassVector.begin(); curHistogramClass != HistogramClassVector.end(); curHistogramClass++) { cout << "Processing histograms in class: <" << colorwhite << (*curHistogramClass)->histogramdescription << colorreset << " >" << endlr; - binNoise((*curHistogramClass)); + binNoise((*curHistogramClass)); + if (noisethresholdborderADU > -1) + (*curHistogramClass)->noisethresholdborder = noisethresholdborderADU; + if (noisethresholdborderE > -1) { + if (labbook.gainDB > 0) { + (*curHistogramClass)->noisethresholdborder = noisethresholdborderE / labbook.gainDB; + } else { + cout << colorred << "Could not set noise threshold border in units of electrons to " << noisethresholdborderE << "e, no calibration done yet. Please rerun after first analysis." <calculteCCE(true); @@ -407,10 +416,11 @@ Bool_t Run::analyzeRun(Bool_t force) } } cout << "---------------------------------" << endl; - cout << colorwhite << "integrateSr90Spectra():" << endlr; -// histogramthreshold->noisethresholdborder = 35.0; // TODO remove the overwriting of automatic noise threshold setting - histogramthreshold->integrateSr90Spectra(histogramthreshold->Seed, frames_found, -1, true); - histogramfixedthreshold->integrateSr90Spectra(histogramfixedthreshold->Seed, frames_found, 0); + if (labbook.source.Contains("Sr90")) { + cout << colorwhite << "integrateSr90Spectra():" << endlr; + histogramthreshold->integrateSr90Spectra(histogramthreshold->Seed, frames_found, -1, true); + histogramfixedthreshold->integrateSr90Spectra(histogramfixedthreshold->Seed, frames_found, 0); + } rescaleHistogramClasses(); cout << colorwhite << "updateDatabase():" << endlr; updateDatabase(); @@ -422,6 +432,20 @@ Bool_t Run::analyzeRun(Bool_t force) return 1; } +Bool_t Run::setNoisethresholdborderADU(Float_t noisethresholdborder) +{ + noisethresholdborderADU = noisethresholdborder; + cout << " Noise threshold border: " << colorwhite << noisethresholdborder << " ADU"<< endlr; + return 0; +} + +Bool_t Run::setNoisethresholdborderE(Float_t noisethresholdborder) +{ + noisethresholdborderADU = noisethresholdborder; + cout << " Noise threshold border: " << colorwhite << noisethresholdborder << " e"<< colorreset << " <-- only used if GAIN is allready saved in database, rerun after first analysis to make sure" << endl; + return 0; +} + Bool_t Run::rescaleHistogramClasses() { float_t vetopeakposition = -1; @@ -845,7 +869,7 @@ Bool_t Run::binSeedSumVeto() for (vector::iterator curHistogramClass = HistogramClassVector.begin(); curHistogramClass != HistogramClassVector.end(); curHistogramClass++) { Double_t* parameters = (Double_t *)calloc(3, sizeof(Double_t)); - (*curHistogramClass)->FindNoisethresholdborder((*curHistogramClass)->Seed, false); + (*curHistogramClass)->FindNoisethresholdborder((*curHistogramClass)->Seed, false, false); if (labbook.source.Contains("Fe")||labbook.source.Contains("Cd")) parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->Veto, "gaus", true); (*curHistogramClass)->posVeto = parameters[1]; diff --git a/MABS_run_analyzer/Run.h b/MABS_run_analyzer/Run.h index cba89f5..f02ccd0 100644 --- a/MABS_run_analyzer/Run.h +++ b/MABS_run_analyzer/Run.h @@ -205,6 +205,12 @@ private: Bool_t dynamicalNoise = 1; /** @brief set this variable to true, to use the common mode noise filter in MAPS:analyzeRun() */ Bool_t commonModeFilter = 1; + + /** @brief overwrite value from which the integral should be taken, overwrite noise cut in electrons */ + Float_t noisethresholdborderE = -1; + + /** @brief overwrite value from which the integral should be taken, overwrite noise cut in ADU */ + Float_t noisethresholdborderADU = -1; Bool_t rescaleHistogramClasses(); @@ -305,6 +311,16 @@ public: */ Bool_t setDynamicalNoiseMode(TString mode); + /** + * @brief set noise threshold border in ADU + */ + Bool_t setNoisethresholdborderADU(Float_t noisethresholdborder); + + /** + * @brief set noise threshold border in e + */ + Bool_t setNoisethresholdborderE(Float_t noisethresholdborder); + /// stores information from the SQL database of a given run labbooksctruct labbook; -- 2.43.0