From: Benjamin Linnik Date: Mon, 4 May 2015 17:50:24 +0000 (+0200) Subject: Analysis: Improved Veto peak finding routine, added user warnings X-Git-Url: https://jspc29.x-matter.uni-frankfurt.de/git/?a=commitdiff_plain;h=3878800506dc8bee578700151fb163da7f2ee77d;p=radhard.git Analysis: Improved Veto peak finding routine, added user warnings --- diff --git a/MABS_run_analyzer/Run.c b/MABS_run_analyzer/Run.c index a4b3e97..77de290 100644 --- a/MABS_run_analyzer/Run.c +++ b/MABS_run_analyzer/Run.c @@ -211,7 +211,7 @@ Bool_t Run::analyzeRun(Bool_t force) } else { - cout << "\033[1;33mSkipped analysis of run " << labbook.runnumber << ", I think the root file for this run allready exists.\033[0m" << endl; + cout << "Skipped analysis of run " << labbook.runnumber << ", I think the root file for this run allready exists." << endl; cout << colorwhite << "initOldRootFile():" << endlr; error = processed->initOldRootFile(); } @@ -762,59 +762,72 @@ Bool_t Run::plotAllHistogramsCalibrated() return 1; } -Float_t Run::FitPerform(TH1F* histogrampointer, TString fitFuncType) +Float_t Run::FitPerform(TH1F* histogrampointer, TString fitFuncType, Bool_t verbose) { Float_t posMax = 0; Float_t posMaxValHist = histogrampointer->GetXaxis()->GetXmax(); + Float_t noiseborder = posMaxValHist/15; if (doFits) { - histogrampointer->GetXaxis()->SetRange(histogrampointer->GetXaxis()->FindBin(posMaxValHist/15),histogrampointer->GetXaxis()->FindBin(posMaxValHist)); // look only for maxima with x greater than posMaxValHist/15, cut away noise + 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,posMaxValHist/15,posMaxValHist); + TF1* fitFunc = new TF1("fitFunc",fitFuncType,noiseborder,posMaxValHist); if (fitFuncType.Contains("gaus")) { if (TString(histogrampointer->GetName()).Contains("Veto")) { -// histogrampointer->GetXaxis()->SetRange(200,histogrampointer->GetXaxis()->FindBin(posMaxValHist)); // look only for maxima with x greater than posMaxValHist/15, cut away noise - - Float_t peak11 = histogrampointer->GetXaxis()->GetBinCenter(histogrampointer->GetMaximumBin()); Float_t peak1 = histogrampointer->GetMaximumBin(); - Float_t peak1val = histogrampointer->GetBinContent(histogrampointer->GetMaximumBin()); - Float_t avg = 0; - for(Int_t bin=histogrampointer->GetXaxis()->FindBin(posMaxValHist/15);binFindLastBinAbove(0);bin++) + Float_t peak2; + if (histogrampointer->GetBinCenter(peak1)<200) { - avg += histogrampointer->GetBinContent(bin); + 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/4); } - avg /= histogrampointer->FindLastBinAbove(0) - histogrampointer->GetXaxis()->FindBin(posMaxValHist/15); - Float_t peak2 = histogrampointer->FindLastBinAbove(avg/4); - Float_t peak21 = histogrampointer->GetXaxis()->GetBinCenter(peak2); - Float_t peak2val = histogrampointer->GetBinContent(peak2); - histogrampointer->GetXaxis()->SetRange(peak1,peak2); // look only for maxima with x greater than posMaxValHist/15, cut away noise - // - cout << coloryellow << "peak1 " << peak11 << endlr; - cout << coloryellow << "peak1val " << peak1val << endlr; - cout << coloryellow << "peak2 " << peak21 << endlr; - cout << coloryellow << "peak2val " << peak1val << endlr; - cout << coloryellow << "avg " << avg << endlr; - + 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->GetXaxis()->SetRange(min,histogrampointer->FindLastBinAbove(0)); // look only for maxima with x greater than posMaxValHist/15, cut away noise - - - + } histogrampointer->Fit(fitFunc, "N,Q,W", "", histogrampointer->GetXaxis()->GetBinCenter(min), posMaxValHist); + 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", "", posMaxValHist/15, posMaxValHist); + histogrampointer->Fit(fitFunc, "N,Q,W", "", noiseborder, posMaxValHist); } Float_t sigma = fitFunc->GetParameter(2); - if (sigma > 50) + if (sigma > 260) { - cout << coloryellow << "sigma suspiciously height when fitting " << histogrampointer->GetName() << " spectrum: " << sigma << endlr; - return 0; + if (verbose) + { + cout << "Sigma suspiciously height when fitting " << histogrampointer->GetName() << " spectrum: " << sigma << endl; + cout << colorred << " Could not find " << histogrampointer->GetName() << " peak" << endlr; + } + return posMax; } + else if (sigma > 40) + { + if (verbose) + { + cout << "Sigma suspiciously height when fitting " << histogrampointer->GetName() << " spectrum: " << sigma << endl; + cout << coloryellow << " Please check " << histogrampointer->GetName() << " peak position: " << colorreset << colorwhite << fitFunc->GetMaximumX() << endlr; + } + } posMax = fitFunc->GetMaximumX(); fitFunc->DrawCopy("same"); TString legendEntry = TString(Form("%s",runcode.Data())); @@ -833,13 +846,10 @@ Float_t Run::FitPerform(TH1F* histogrampointer, TString fitFuncType) Float_t fitMax3 = 1000; Float_t minFitMax = 1000; Float_t maxFitMax = 1000; - - cout << "posMaxValHist/15: " << posMaxValHist/15 << endl; - - histogrampointer->Fit(fitFunc, "N,Q,W", "", posMaxValHist/15, posMaxValHist); + histogrampointer->Fit(fitFunc, "N,Q,W", "", noiseborder, posMaxValHist); fitMax1 = fitFunc->GetMaximumX(); fitFunc->DrawCopy("same"); - histogrampointer->Fit(fitFunc, "N,Q,W", "", posMaxValHist/15, fitMax1*1.1); + histogrampointer->Fit(fitFunc, "N,Q,W", "", noiseborder, fitMax1*1.1); fitMax2 = fitFunc->GetMaximumX(); fitFunc->SetLineColor(kBlue); fitFunc->SetLineStyle(2); // dashed @@ -895,7 +905,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(onehistogram, fitFuncType, false); 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);