From ae499b85c30b9557939eb3267d20a4a2a964d17b Mon Sep 17 00:00:00 2001 From: Benjamin Linnik Date: Tue, 28 Nov 2017 14:53:52 +0100 Subject: [PATCH] Analyzer: Bugfixes, Added 'leak' command line flag for leakage current analysis --- MABS_run_analyzer/ChargeSpektrum.c | 21 +++++- MABS_run_analyzer/ChargeSpektrumFunctions.c | 36 +++++++++- MABS_run_analyzer/HistogramType.c | 76 +++++++++++---------- MABS_run_analyzer/Run.c | 23 ++----- 4 files changed, 99 insertions(+), 57 deletions(-) diff --git a/MABS_run_analyzer/ChargeSpektrum.c b/MABS_run_analyzer/ChargeSpektrum.c index b24821b..ac756af 100644 --- a/MABS_run_analyzer/ChargeSpektrum.c +++ b/MABS_run_analyzer/ChargeSpektrum.c @@ -85,9 +85,13 @@ void ChargeSpektrum(TString runnumber = "") // f0 analysis variables vector compareHistogramAllRunsVectorF0; + + // leakage current analysis variables + vector compareHistogramLeakageCurrent; for(Int_t runi=0; runi0) { runs[runi] = new Run(runList[runi], runi); @@ -148,7 +152,7 @@ void ChargeSpektrum(TString runnumber = "") gROOT->SetBatch(kFALSE); // use only lower case when matching anaylsis types - + if (analysisType.Contains("clas")) { // classic analysis compareHistogramClassVectorClassic.push_back(runs[runi]->histogram->normalized); if (runi+1 == numberRuns) @@ -410,6 +414,21 @@ void ChargeSpektrum(TString runnumber = "") compareHistogramAllRunsVectorF0.push_back(runs[runi]->averageF0Distr); runs[runi]->plot1DHistogram( runs[runi]->averageF0Distr, "gaus", true); } + + + + if (analysisType.Contains("leak")) { // analyze cluster formation + // Plot leakage current + compareHistogramLeakageCurrent.push_back(runs[runi]->histogramwoRTS->LeakageCurrentDistrib); + if (runi+1 == numberRuns) { + compareHistogramVectorVector.push_back(compareHistogramLeakageCurrent); + } + + runs[runi]->plot1DHistogram(runs[runi]->histogramdynamicalthreshold->Veto, "GaussTail", true); + runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS->LeakageCurrentInPixelSorted, "", false, true); + } + + // write all histograms to a dat file, can be imported by ORIGIN or Excel. Leave this line ;) runs[runi]->writeAllHistogramsToFile(); diff --git a/MABS_run_analyzer/ChargeSpektrumFunctions.c b/MABS_run_analyzer/ChargeSpektrumFunctions.c index f7ed149..09c711f 100644 --- a/MABS_run_analyzer/ChargeSpektrumFunctions.c +++ b/MABS_run_analyzer/ChargeSpektrumFunctions.c @@ -83,6 +83,9 @@ TCanvas* plot1DData(std::vector* ntuplepvectorpointer, Bool_t verbose Bool_t writeNTupleVectorToFile(std::vector* ntuplepvectorpointer); Bool_t writeObservableToFile(vector* ptCompareHistogramVector); + +// Write one histogram to its own file +Bool_t writeHistogramToFile(TH1FO* onehistogram); /** @brief A vector able to hold pointer of type #Run::histogramstruct * * You push in pointer of the type Run::histogramstruct in ChargeSpectrum.C and then @@ -478,13 +481,35 @@ Bool_t writeOneHistogramTypeToFile(vector* ptCompareHistogramCla } +Bool_t writeHistogramToFile(TH1FO* onehistogram) +{ + Run* runp = onehistogram->itsHistogramType->run; + system("mkdir \""+ runp->savepathresults + "\" -p"); + TString filename= runp->savepathresults + "/" + runp->runcode + " " + onehistogram->GetName() + " hist.dat"; + fstream* fout = new fstream(filename,ios::out); + + TString header = Form("#bin [ADU]\tCounts"); + *fout << header << endl; + + TString outline; + for(Int_t bin=1;bin<=onehistogram->GetNbinsX();bin++) + { + outline=Form("%.3f\t%.3f", onehistogram->GetBinCenter(bin), onehistogram->GetBinContent(bin)); + *fout<close(); + return 0; +} + + Bool_t writeObservableToFile(vector* ptCompareHistogramVector) { if (numberRuns>0) { - system("mkdir \""+ savepathresults + folderadd + "/\"" + " -p"); - TString filename = savepathresults + folderadd + "/Spectra"; + TH1FO* histogrampfirst = ptCompareHistogramVector->at(0); + system("mkdir \""+ savepathresults + folderadd + "/\"" + " -p"); + TString filename = savepathresults + folderadd + "/Spectra " + histogrampfirst->GetName(); TString header = ""; for(Int_t runi=0;runi* ptCompareHistogramVector) } *fout << header << endl; TString outline; - for(Int_t bin=1;bin<=runs[0]->cursystemparam.nbins;bin++) + for(Int_t bin=1;bin<=histogrampfirst->GetNbinsX();bin++) { outline =""; for (vector::iterator curHistogram = ptCompareHistogramVector->begin(); curHistogram != ptCompareHistogramVector->end(); curHistogram++) { @@ -514,6 +539,11 @@ Bool_t writeObservableToFile(vector* ptCompareHistogramVector) *fout<close(); + + +// for (vector::iterator curHistogram = ptCompareHistogramVector->begin(); curHistogram != ptCompareHistogramVector->end(); curHistogram++) { +// writeHistogramToFile(curHistogram); +// } return 0; } diff --git a/MABS_run_analyzer/HistogramType.c b/MABS_run_analyzer/HistogramType.c index 0a5c06f..9002c0c 100644 --- a/MABS_run_analyzer/HistogramType.c +++ b/MABS_run_analyzer/HistogramType.c @@ -747,25 +747,28 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType cout << "peak: " << fitFcn->GetParameter(2) << endl; cout << "y-shift:" << fitFcn->GetParameter(3) << endl; // } - } else if (fitFuncType=="GaussTail") - { - histogrampointer->GetXaxis()->UnZoom(); - posMaxValHist = (histogrampointer->GetBinCenter(histogrampointer->FindLastBinAbove(1,1))); + } else if (fitFuncType=="GaussTail") // this option is used when fitting the calibration peak + { + // for Cadmium, the Veto spectrum is not so clean, set the noise border to higher values + // this way the other peaks which come not from the K-Alpha line will be ignored. + if (labbook->source.Contains("Cd")) { + noiseborder = 100; + cout << "Set noiseborder to: " << noiseborder << endlr; + } + + histogrampointer->GetXaxis()->UnZoom(); // resets range + posMaxValHist = (histogrampointer->GetBinCenter(histogrampointer->FindLastBinAbove(1,1))); - if (verbose) { + if (verbose) { cout << colorcyan << "maxvalhist: " << posMaxValHist << endlr; cout << colorcyan << "noiseborder: " << noiseborder << endlr; - } -// histogrampointer->GetXaxis()->SetRange(noiseborder,histogrampointer->FindLastBinAbove(0)); // look only for maxima with x greater than noiseborder, 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 -// histogrampointer->GetXaxis()->SetRangeUser(noiseborder,posMaxValHist); TF1* fitFunc = new TF1("fitFunc",GaussTail,noiseborder,posMaxValHist,6); parameters = (Double_t *)calloc(10, sizeof(Double_t)); const Double_t def_amplitude=histogrampointer->GetBinContent(histogrampointer->GetMaximumBin()); const Double_t def_peakcenter=histogrampointer->GetBinCenter(histogrampointer->GetMaximumBin()); -// cout << colorcyan << "def_peakcenter: " << def_peakcenter << endlr; -// cout << colorcyan << "histogrampointer->GetMaximumBin(): " << histogrampointer->GetMaximumBin() << endlr; const Double_t def_gausssig=8.68052; const Double_t def_distgauss=20.4; const Double_t def_bgslope=0; @@ -829,34 +832,37 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType // cout << colorcyan << gMinuit->fCstatu.Data() << " fit tries: " << fittries << endlr; } while (fittries < 6); fittries = 0; - do { - failed = false; - fit_result_ptr = histogrampointer->Fit(fitFunc, "N,M,W,Q,S", "", posMaxValHist-posMaxValHist/4*(++fittries), posMaxValHist); - // cout << colorcyan << " AFTER fit " << endlr; - if (gMinuit == nullptr) { - failed = true; - } else - { - if (gMinuit->fCstatu == "NOT POSDEF" || gMinuit->fCstatu.Contains("FAILED")) { + if (failed) + { + do { + failed = false; + fit_result_ptr = histogrampointer->Fit(fitFunc, "N,M,W,Q,S", "", posMaxValHist-posMaxValHist/4*(++fittries), posMaxValHist); + // cout << colorcyan << " AFTER fit " << endlr; + if (gMinuit == nullptr) { failed = true; + } else + { + if (gMinuit->fCstatu == "NOT POSDEF" || gMinuit->fCstatu.Contains("FAILED")) { + failed = true; + } } - } - if (failed) - { - fitFunc->SetParameter(0,def_amplitude); - fitFunc->SetParameter(1,posMaxValHist-posMaxValHist/4*(fittries)); - fitFunc->SetParameter(2,def_gausssig); - fitFunc->SetParameter(3,def_distgauss); - fitFunc->SetParameter(4,def_bgslope); - fitFunc->SetParameter(5,def_bgoffs); - } - else - fittries = 100; - } while (fittries < 4); + if (failed) + { + fitFunc->SetParameter(0,def_amplitude); + fitFunc->SetParameter(1,posMaxValHist-posMaxValHist/4*(fittries)); + fitFunc->SetParameter(2,def_gausssig); + fitFunc->SetParameter(3,def_distgauss); + fitFunc->SetParameter(4,def_bgslope); + fitFunc->SetParameter(5,def_bgoffs); + } + else + fittries = 100; + } while (fittries < 4); + } if (failed) { failed = false; - cout << "Most emphasized estimate" << endlr; + cout << "More emphasized estimate" << endlr; histogrampointer->GetXaxis()->SetRange(0,posMaxValHist); // look only for maxima with x greater than noiseborder, cut away noise fitFunc->SetParameter(0,def_amplitude); fitFunc->SetParameter(1,def_peakcenter); @@ -952,8 +958,8 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType fitFunc->SetLineStyle(1); // normal for the following fits if (verbose) { - TCanvas* canvas = new TCanvas(histogrampointer->GetName(), histogrampointer->GetName(), 900, 700); - histogrampointer->Draw(); +// TCanvas* canvas = new TCanvas(histogrampointer->GetName(), histogrampointer->GetName(), 900, 700); +// histogrampointer->Draw(); fitFunc->Draw("SAME"); } diff --git a/MABS_run_analyzer/Run.c b/MABS_run_analyzer/Run.c index b357019..bfa8e9b 100644 --- a/MABS_run_analyzer/Run.c +++ b/MABS_run_analyzer/Run.c @@ -1150,9 +1150,10 @@ Bool_t Run::binLeakageCurrent(HistogramType* oneHistogramClass) std::vector sortedHistogram; // for (Int_t pixeli=1; pixeli <= oneHistogramClass->LeakageCurrentInPixel->GetNbinsX(); pixeli++) { - oneHistogramClass->LeakageCurrentDistrib->Fill(oneHistogramClass->LeakageCurrentInPixel->GetBinContent(pixeli)); - if (oneHistogramClass->LeakageCurrentInPixel->GetBinContent(pixeli) != 0) // != 0 // not RTS pixel + if (oneHistogramClass->LeakageCurrentInPixel->GetBinContent(pixeli) != 0) { // != 0 // not RTS pixel + oneHistogramClass->LeakageCurrentDistrib->Fill(oneHistogramClass->LeakageCurrentInPixel->GetBinContent(pixeli)); sortedHistogram.push_back(oneHistogramClass->LeakageCurrentInPixel->GetBinContent(pixeli)); + } // cout << oneHistogramClass->LeakageCurrentInPixel->GetBinContent(pixeli) << " " << endl; } std::sort(sortedHistogram.begin(),sortedHistogram.end()); @@ -1698,25 +1699,11 @@ Bool_t Run::binSeedSumVeto() for (vector::iterator curHistogramClass = HistogramClassVector.begin(); curHistogramClass != HistogramClassVector.end(); curHistogramClass++) { Double_t* parameters = (Double_t *)calloc(11, sizeof(Double_t)); // allocate 11 parameters for safety, maximum 10 used at the moment - cout << colorcyan << (*curHistogramClass)->histogramdescription << endlr; +// cout << colorcyan << (*curHistogramClass)->histogramdescription << endlr; (*curHistogramClass)->FindNoisethresholdborder((*curHistogramClass)->Seed, false, false); - if (labbook.runnumber == 343056) - (*curHistogramClass)->noisethresholdborder = 34; - if (labbook.source.Contains("Fe")||labbook.source.Contains("Cd")) { - + if (labbook.source.Contains("Fe")||labbook.source.Contains("Cd")) { parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->Veto, "GaussTail"); (*curHistogramClass)->posVeto = parameters[1]; - if (labbook.runnumber == 343056) - (*curHistogramClass)->posVeto = 234; - if (labbook.runnumber == 345045) - (*curHistogramClass)->posVeto = 425; - if (labbook.runnumber == 342824) - (*curHistogramClass)->posVeto = 425; - if (labbook.runnumber == 342265) - (*curHistogramClass)->posVeto = 417; - if (labbook.runnumber == 345054) - (*curHistogramClass)->posVeto = 417; - // cout << (*curHistogramClass)->histogramdescription << ": " << colorcyan << parameters[1] << endlr; (*curHistogramClass)->integralVeto = parameters[6]; } if (labbook.chipGen.EqualTo("Pipper2") || labbook.chipGen.EqualTo("Pegasus")) -- 2.43.0