From c25fb0e8824f2ca547e270201a51ae853c490458 Mon Sep 17 00:00:00 2001 From: Benjamin Linnik Date: Wed, 13 Dec 2017 18:46:46 +0100 Subject: [PATCH] Analyzer: calibration is done for each cut type seperately, better fitting of sum and veto peak: tested with Mi19 --- MABS_run_analyzer/ChargeSpektrum.c | 82 ++- MABS_run_analyzer/ChargeSpektrumFunctions.c | 174 ++--- MABS_run_analyzer/HistogramType.c | 540 ++++++++-------- MABS_run_analyzer/Run.c | 674 ++++++++++---------- MABS_run_analyzer/Run.h | 4 +- MABS_run_analyzer/help.h | 12 +- 6 files changed, 754 insertions(+), 732 deletions(-) diff --git a/MABS_run_analyzer/ChargeSpektrum.c b/MABS_run_analyzer/ChargeSpektrum.c index 1d8801b..497ff54 100644 --- a/MABS_run_analyzer/ChargeSpektrum.c +++ b/MABS_run_analyzer/ChargeSpektrum.c @@ -78,6 +78,7 @@ void ChargeSpektrum(TString runnumber = "") // variables for RTS analysis vector compareHistogramVectorRTSSeed, compareHistogramVectorRTSSum; + vector compareHistogramClassVectorRTS; // variables for cluster multiplicity vector compareHistogramVectorClusterMultiplicity; @@ -151,13 +152,27 @@ void ChargeSpektrum(TString runnumber = "") if (!isBatch) gROOT->SetBatch(kFALSE); + // **************************************************************** + // + // ADD own anaylsis below + // + // ***************************************************************** // use only lower case when matching anaylsis types + // **************************************************************** + + // default histogram type to be used in the anaylsis below + HistogramType* HistogramTypeDefaultPt = runs[runi]->histogram; + if (analysisType.Contains("rts")) { + HistogramTypeDefaultPt = runs[runi]->histogramwoRTS; + } + + // ################################### if (analysisType.Contains("clas")) { // classic analysis compareHistogramClassVectorClassic.push_back(runs[runi]->histogram->normalized); if (runi+1 == numberRuns) compareHistogramClassVectorVector.push_back(compareHistogramClassVectorClassic); - } + } if (analysisType.Contains("clusthresh")) { // cluster charge > 2 * noise in cluster compareHistogramClassVectorClusterThreshold.push_back(runs[runi]->histogramthreshold->normalized); @@ -174,7 +189,8 @@ void ChargeSpektrum(TString runnumber = "") } if (analysisType.Contains("cali")) { // calibrated analysis - compareHistogramClassVectorCalibrated.push_back(runs[runi]->histogram->normalized->calibrated); + runs[runi]->plot1DHistogram( HistogramTypeDefaultPt->Veto, "GaussTail", true, false, false, HistogramTypeDefaultPt->noisethresholdborder); + compareHistogramClassVectorCalibrated.push_back(HistogramTypeDefaultPt->normalized->calibrated); if (runi+1 == numberRuns) compareHistogramClassVectorVector.push_back(compareHistogramClassVectorCalibrated); } @@ -217,12 +233,12 @@ void ChargeSpektrum(TString runnumber = "") TNtupleO *numpixelincl; numpixelincl =new TNtupleO("Average charge collected", "data", "x:y:xerr:yerr"); - numpixelincl->itsHistogramType = runs[runi]->histogram; + numpixelincl->itsHistogramType = HistogramTypeDefaultPt; numpixelincl->xaxisTitle = "Number of pixel in cluster"; numpixelincl->yaxisTitle = "Average charge collected [e]"; for (UInt_t sumi = 0; sumi < (*runs[runi]->histogramdynamicalthreshold->a_Sum).size(); sumi++) { if (!(sumi%2)) - compareHistogramVectorDepletionSum.push_back((*runs[runi]->histogram->normalized->a_Sum)[sumi]); + compareHistogramVectorDepletionSum.push_back((*HistogramTypeDefaultPt->normalized->a_Sum)[sumi]); numpixelincl->Fill(sumi+1,AverageCharge[sumi],0,0); } ntupleVectorDepltetion3.push_back(numpixelincl); @@ -233,7 +249,11 @@ void ChargeSpektrum(TString runnumber = "") } } - if (analysisType.Contains("rts")) { // rts studies + if (analysisType.Contains("rts")) { // rts studies + compareHistogramClassVectorRTS.push_back(runs[runi]->histogramwoRTS->normalized); + if (runi+1 == numberRuns) + compareHistogramClassVectorVector.push_back(compareHistogramClassVectorRTS); + vector compareHistogramClassVectorRTS; compareHistogramClassVectorRTS.push_back(runs[runi]->histogramwoRTS); compareHistogramClassVectorRTS.push_back(runs[runi]->histogram); @@ -244,7 +264,7 @@ void ChargeSpektrum(TString runnumber = "") compareHistogramVectorRTSSeed.push_back(runs[runi]->histogramwoRTS->normalized->Seed); compareHistogramVectorRTSSum.push_back(*(runs[runi]->histogramwoRTS->normalized->Sum)); - runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS->pixeltimefiredDistrib); + runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS->pixeltimefiredDistrib, "gaus", false, true); runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS->LeakageCurrentInPixelSorted); if (runi+1 == numberRuns) { @@ -254,19 +274,19 @@ void ChargeSpektrum(TString runnumber = "") } if (analysisType.Contains("seedf")) { // seedfit: seed integral anaylsis - runs[runi]->plot1DHistogram( runs[runi]->histogram->normalized->Seed, "gaus", true, false, false, runs[runi]->histogram->fixedThresholdValue); - runs[runi]->plot1DHistogram( runs[runi]->histogramdynamicalthreshold->normalized->Seed, "gaus", true, false, false, runs[runi]->histogram->fixedThresholdValue); + runs[runi]->plot1DHistogram( HistogramTypeDefaultPt->normalized->Seed, "gaus", true, false, false, HistogramTypeDefaultPt->fixedThresholdValue); + runs[runi]->plot1DHistogram( runs[runi]->histogramdynamicalthreshold->normalized->Seed, "gaus", true, false, false, HistogramTypeDefaultPt->fixedThresholdValue); vector compareHistogramVectorSeedFit; - compareHistogramVectorSeedFit.push_back(runs[runi]->histogram->Seed); + compareHistogramVectorSeedFit.push_back(HistogramTypeDefaultPt->Seed); compareHistogramVectorSeedFit.push_back(runs[runi]->histogramdynamicalthreshold->Seed); CompareHistograms(&compareHistogramVectorSeedFit); compareHistogramVectorSeedFit.clear(); - datainputIntSeed->itsHistogramType = runs[runi]->histogram->normalized; + datainputIntSeed->itsHistogramType = HistogramTypeDefaultPt->normalized; datainputIntSeed->xaxisTitle = "Voltage"; datainputIntSeed->yaxisTitle = "Integral of Seed"; - datainputIntSeed->Fill(runs[runi]->labbook.depletionV,runs[runi]->histogram->integralSeed,0,runs[runi]->histogram->integralSeedErr); + datainputIntSeed->Fill(runs[runi]->labbook.depletionV,HistogramTypeDefaultPt->integralSeed,0,HistogramTypeDefaultPt->integralSeedErr); if (runi+1 == numberRuns) { ntupleVectorSeed.push_back(datainputIntSeed); ntupleVectorVector.push_back(ntupleVectorSeed); @@ -274,19 +294,19 @@ void ChargeSpektrum(TString runnumber = "") } if (analysisType.Contains("sumf")) { // seedfit: seed integral anaylsis - runs[runi]->plot1DHistogram(*(runs[runi]->histogram->normalized->Sum), "gaus", true, false, false, runs[runi]->histogram->fixedThresholdValue); - runs[runi]->plot1DHistogram(*(runs[runi]->histogramdynamicalthreshold->normalized->Sum), "gaus", true, false, false, runs[runi]->histogram->fixedThresholdValue); + runs[runi]->plot1DHistogram(*(HistogramTypeDefaultPt->normalized->Sum), "gaus", true, false, false, HistogramTypeDefaultPt->fixedThresholdValue); + runs[runi]->plot1DHistogram(*(runs[runi]->histogramdynamicalthreshold->normalized->Sum), "gaus", true, false, false, HistogramTypeDefaultPt->fixedThresholdValue); vector compareHistogramVectorSumFit; - compareHistogramVectorSumFit.push_back(*(runs[runi]->histogram->Sum)); + compareHistogramVectorSumFit.push_back(*(HistogramTypeDefaultPt->Sum)); compareHistogramVectorSumFit.push_back(*(runs[runi]->histogramdynamicalthreshold->Sum)); CompareHistograms(&compareHistogramVectorSumFit); compareHistogramVectorSumFit.clear(); - datainputIntSum->itsHistogramType = runs[runi]->histogram->normalized; + datainputIntSum->itsHistogramType = HistogramTypeDefaultPt->normalized; datainputIntSum->xaxisTitle = "Voltage"; datainputIntSum->yaxisTitle = "Integral of Sum"; - datainputIntSum->Fill(runs[runi]->labbook.depletionV,runs[runi]->histogram->integralSum,0,runs[runi]->histogram->integralSumErr); + datainputIntSum->Fill(runs[runi]->labbook.depletionV,HistogramTypeDefaultPt->integralSum,0,HistogramTypeDefaultPt->integralSumErr); if (runi+1 == numberRuns) { ntupleVectorSum.push_back(datainputIntSum); ntupleVectorVector.push_back(ntupleVectorSeed); @@ -295,22 +315,22 @@ void ChargeSpektrum(TString runnumber = "") if (analysisType.Contains("cluster")) { // analyze cluster formation // show cluster multiplicity - compareHistogramVectorClusterMultiplicity.push_back(runs[runi]->histogram->normalized->clustermultiplicity); + compareHistogramVectorClusterMultiplicity.push_back(HistogramTypeDefaultPt->normalized->clustermultiplicity); // runs[runi]->plot1DHistogram( runs[runi]->histogram->clustermultiplicity); vector compareHistogramVectorCluster; - compareHistogramVectorCluster.push_back(runs[runi]->histogram->Sum9); - compareHistogramVectorCluster.push_back(runs[runi]->histogram->Sum25); + compareHistogramVectorCluster.push_back(HistogramTypeDefaultPt->Sum9); + compareHistogramVectorCluster.push_back(HistogramTypeDefaultPt->Sum25); CompareHistograms(&compareHistogramVectorCluster); compareHistogramVectorCluster.clear(); vector compareHistogramVectorCluster2; - compareHistogramVectorCluster2.push_back(runs[runi]->histogram->Seed); - compareHistogramVectorCluster2.push_back(runs[runi]->histogram->a_Sum9[1]); - compareHistogramVectorCluster2.push_back(runs[runi]->histogram->Sum9); - compareHistogramVectorCluster2.push_back(runs[runi]->histogram->a_Sum25[1]); - compareHistogramVectorCluster2.push_back(runs[runi]->histogram->Sum25); + compareHistogramVectorCluster2.push_back(HistogramTypeDefaultPt->Seed); + compareHistogramVectorCluster2.push_back(HistogramTypeDefaultPt->a_Sum9[1]); + compareHistogramVectorCluster2.push_back(HistogramTypeDefaultPt->Sum9); + compareHistogramVectorCluster2.push_back(HistogramTypeDefaultPt->a_Sum25[1]); + compareHistogramVectorCluster2.push_back(HistogramTypeDefaultPt->Sum25); CompareHistograms(&compareHistogramVectorCluster2); compareHistogramVectorCluster2.clear(); @@ -328,7 +348,7 @@ void ChargeSpektrum(TString runnumber = "") TNtupleO *numpixelincl2; numpixelincl2 =new TNtupleO("Average charge collected", "data", "x:y:xerr:yerr"); - numpixelincl2->itsHistogramType = runs[runi]->histogram; + numpixelincl2->itsHistogramType = HistogramTypeDefaultPt; numpixelincl2->Title = "5x5 cluster"; numpixelincl2->xaxisTitle = "Number of pixel in cluster"; numpixelincl2->yaxisTitle = "Average charge collected [ADU]"; @@ -351,7 +371,7 @@ void ChargeSpektrum(TString runnumber = "") TNtupleO *numpixelincl; numpixelincl =new TNtupleO("Average charge collected", "data", "x:y:xerr:yerr"); - numpixelincl->itsHistogramType = runs[runi]->histogram; + numpixelincl->itsHistogramType = HistogramTypeDefaultPt; numpixelincl->Title = "3x3 cluster"; numpixelincl->xaxisTitle = "Number of pixel in cluster"; numpixelincl->yaxisTitle = "Average charge collected [ADU]"; @@ -429,13 +449,13 @@ void ChargeSpektrum(TString runnumber = "") if (analysisType.Contains("leak")) { // analyze cluster formation // Plot leakage current - compareHistogramLeakageCurrent.push_back(runs[runi]->histogramwoRTS->LeakageCurrentDistrib); + compareHistogramLeakageCurrent.push_back(HistogramTypeDefaultPt->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); + runs[runi]->plot1DHistogram(HistogramTypeDefaultPt->LeakageCurrentInPixelSorted, "", false, true); } @@ -444,6 +464,12 @@ void ChargeSpektrum(TString runnumber = "") runs[runi]->writeAllHistogramsToFile(); } + + // DEBUGGING + // runs[runi]->histogramwoRTS->FindNoisethresholdborder(runs[runi]->histogramwoRTS->Veto, false, true); + // runs[runi]->plot1DHistogram( runs[runi]->histogramwoRTS->Veto, "GaussTail", true, false, false, runs[runi]->histogramwoRTS->noisethresholdborder); +// runs[runi]->plot1DHistogram( runs[runi]->histogramwoRTS->a_Sum25[22], "gaus", true, false, false, runs[runi]->histogramwoRTS->noisethresholdborder); + //cout << runs[runi]->histogram } } diff --git a/MABS_run_analyzer/ChargeSpektrumFunctions.c b/MABS_run_analyzer/ChargeSpektrumFunctions.c index 09c711f..f9aa303 100644 --- a/MABS_run_analyzer/ChargeSpektrumFunctions.c +++ b/MABS_run_analyzer/ChargeSpektrumFunctions.c @@ -132,7 +132,7 @@ void ChargeSpektrum(Int_t runnumber, Int_t runnumber2) void Init() { -// FindGoodTitle(); + // FindGoodTitle(); TTimeStamp timestamp = TTimeStamp(); savepathresults = Form("./results/%d%06d %s", (int)timestamp.GetDate(kFALSE), (int)timestamp.GetTime(kFALSE), analysisType.Data()); @@ -241,7 +241,7 @@ void InterpreteUserInput(TString runnumber, std::vector *runList, std::vect if (runnumber.Contains(";")) { *analysisis = runnumber(runnumber.First(";")+1,runnumber.Length()); analysisis->ToLower(); -// cout << colorcyan << "*analysisis: " << *analysisis << endlr; + // cout << colorcyan << "*analysisis: " << *analysisis << endlr; } } else @@ -321,7 +321,7 @@ plot1DData(std::vector* ntuplepvectorpointer, Bool_t verbose, Bool_t leg1->SetTextFont(132); leg1->SetFillColor(0); leg1->SetBorderSize(0); - + TString canvastitle = Form("%s", ntuplepointer->GetName()); TTimeStamp* time = new TTimeStamp(); TString canvasname = Form("%s_%d",ntuplepointer->GetName(),(int)time->GetNanoSec()/100000); @@ -370,12 +370,12 @@ plot1DData(std::vector* ntuplepvectorpointer, Bool_t verbose, Bool_t owntitle->Draw("SAME"); } canvas->Update(); - // canvas -> SaveAs( savepathresults + "/" + runcode + " " + onehistogram->GetName() + ".eps"); + // canvas -> SaveAs( savepathresults + "/" + runcode + " " + onehistogram->GetName() + ".eps"); TImageDump *img = new TImageDump(savepathresults + folderadd + "/" + ntuplepointer->GetName() + ".png"); canvas->Paint(); img->Close(); - + writeNTupleVectorToFile(ntuplepvectorpointer); return canvas; } @@ -541,9 +541,9 @@ Bool_t writeObservableToFile(vector* ptCompareHistogramVector) fout->close(); -// for (vector::iterator curHistogram = ptCompareHistogramVector->begin(); curHistogram != ptCompareHistogramVector->end(); curHistogram++) { -// writeHistogramToFile(curHistogram); -// } + // for (vector::iterator curHistogram = ptCompareHistogramVector->begin(); curHistogram != ptCompareHistogramVector->end(); curHistogram++) { + // writeHistogramToFile(curHistogram); + // } return 0; } @@ -556,12 +556,12 @@ Bool_t CompareHistograms(vector* ptCompareHistogramVector, TString title { FindGoodTitle(ptCompareHistogramVector); cout << "Returned from find good title" << endlr; - + gROOT->SetStyle("RadHard_NoTitle"); // legend entries Float_t height = ptCompareHistogramVector->size() * 0.03; -// TLegend* leg1 = new TLegend(0.3,0.89-height,0.95,0.89);//(0.6,0.7,0.89,0.89); + // TLegend* leg1 = new TLegend(0.3,0.89-height,0.95,0.89);//(0.6,0.7,0.89,0.89); TLegend* leg1 = new TLegend(0.7,0.89-height,0.95,0.89);//(0.6,0.7,0.89,0.89); leg1->SetTextSize(0.025); leg1->SetFillStyle(1001); @@ -620,13 +620,13 @@ Bool_t CompareHistograms(vector* ptCompareHistogramVector, TString title leg1->Draw("SAME"); curhistogramclone->SetAxisRange(0,lastbin1*1.1,"X"); -// curhistogramclone->GetYaxis()->UnZoom(); + // curhistogramclone->GetYaxis()->UnZoom(); curhistogramclone->GetYaxis()->SetRangeUser(1,heighestval1*1.4); if (curhistogramclone->GetNbinsX() < 20) { curhistogramclone->GetXaxis()->UnZoom(); curhistogramclone->GetYaxis()->UnZoom(); } -// canvas->Update(); + // canvas->Update(); owntitle->Clear(); owntitle->AddText(trimRunnumberAtBegin(curhistogramclone->GetName())); @@ -637,11 +637,11 @@ Bool_t CompareHistograms(vector* ptCompareHistogramVector, TString title } owntitle->Draw("SAME"); -// gPad->SetLogy(1); -// curhistogramclone->GetYaxis()->UnZoom(); + // gPad->SetLogy(1); + // curhistogramclone->GetYaxis()->UnZoom(); } canvas->Update(); -// gPad->Modified(); gPad->Update(); + // gPad->Modified(); gPad->Update(); MSaveBigPNG(canvas, savepathresults + folderadd + "/" + canvastitle + ".png"); TImageDump *img = new TImageDump(savepathresults + folderadd + "/" + canvastitle + ".png"); @@ -946,7 +946,7 @@ Bool_t plotAllRuns(vector* ptCompareHistogramClassVector) curhistogramclone->SetAxisRange(0,lastbin1*1.1,"X"); curhistogramclone->SetAxisRange(0,heighestval1*1.1,"Y"); curhistogramclone->GetYaxis()->SetRangeUser(1,heighestval1*1.2); -// gPad->SetLogy(1); + // gPad->SetLogy(1); gPad->SetGrid(1); gPad->SetGridy(1); owntitle->Clear(); @@ -972,7 +972,7 @@ Bool_t plotAllRuns(vector* ptCompareHistogramClassVector) // heighestval2 = (curhistogramclone->GetMaximum()>heighestval2?curhistogramclone->GetMaximum():heighestval2); curhistogramclone->GetYaxis()->SetRangeUser(1,heighestval2*1.2); curhistogramclone->Draw("SAME"); -// gPad->SetLogy(1); + // gPad->SetLogy(1); gPad->SetGrid(1); gPad->SetGridy(1); @@ -1106,12 +1106,12 @@ Bool_t FindGoodTitle(vector* ptCompareHistogramClassVector, vect if (ptCompareHistogramVector->size() == ptCompareHistogramClassVector->size()) { mayBeSameRun = kTRUE; // no different classes were provided, only some TH1F histograms. Classes were constructed afterwards -// same_RunNumber = kTRUE; + // same_RunNumber = kTRUE; } -// cout << "Trying to find good title" << endlr; - - headerStringsVector.clear(); + // cout << "Trying to find good title" << endlr; + + headerStringsVector.clear(); legendStringsVector.clear(); if (ptCompareHistogramClassVector->size() == 0 ) @@ -1139,7 +1139,7 @@ Bool_t FindGoodTitle(vector* ptCompareHistogramClassVector, vect headerStringsVector.push_back(((TObjString *)(humanreadablestrings->At(0)))->String()); headerStringsVector.push_back(((TObjString *)(humanreadablestrings->At(1)))->String()); headerStringsVector.push_back(title2); - + } else { if (mayBeSameRun) firsthistogramp = ptCompareHistogramVector->at(0); @@ -1176,8 +1176,8 @@ Bool_t FindGoodTitle(vector* ptCompareHistogramClassVector, vect same_Clock = kFALSE; if (curlabbook.depletionV != firstlabbook.depletionV) same_Depletion= kFALSE; - // if (curlabbook.system != firstlabbook.system) - // same_System = kFALSE; + // if (curlabbook.system != firstlabbook.system) + // same_System = kFALSE; if (curlabbook.runnumber != firstlabbook.runnumber) same_RunNumber= kFALSE; @@ -1186,19 +1186,19 @@ Bool_t FindGoodTitle(vector* ptCompareHistogramClassVector, vect same_HistType = kFALSE; } -// // DEBUG -// cout << colorcyan << "same_Source: " << same_Source << endlr; -// cout << colorcyan << "same_NonIonRad: " << same_NonIonRad << endlr; -// cout << colorcyan << "same_IonRad: " << same_IonRad << endlr; -// cout << colorcyan << "same_ChipNum: " << same_ChipNum << endlr; -// cout << colorcyan << "same_ChipGen: " << same_ChipGen << endlr; -// cout << colorcyan << "same_Matrix: " << same_Matrix << endlr; -// cout << colorcyan << "same_Temp: " << same_Temp << endlr; -// cout << colorcyan << "same_Clock: " << same_Clock << endlr; -// cout << colorcyan << "same_Depletion: " << same_Depletion << endlr; -// cout << colorcyan << "same_HistType: " << same_HistType << endlr; -// cout << colorcyan << "same_RunNumber: " << same_RunNumber << endlr; -// cout << colorcyan << "mayBeSameRun: " << mayBeSameRun << endlr; + // // DEBUG + // cout << colorcyan << "same_Source: " << same_Source << endlr; + // cout << colorcyan << "same_NonIonRad: " << same_NonIonRad << endlr; + // cout << colorcyan << "same_IonRad: " << same_IonRad << endlr; + // cout << colorcyan << "same_ChipNum: " << same_ChipNum << endlr; + // cout << colorcyan << "same_ChipGen: " << same_ChipGen << endlr; + // cout << colorcyan << "same_Matrix: " << same_Matrix << endlr; + // cout << colorcyan << "same_Temp: " << same_Temp << endlr; + // cout << colorcyan << "same_Clock: " << same_Clock << endlr; + // cout << colorcyan << "same_Depletion: " << same_Depletion << endlr; + // cout << colorcyan << "same_HistType: " << same_HistType << endlr; + // cout << colorcyan << "same_RunNumber: " << same_RunNumber << endlr; + // cout << colorcyan << "mayBeSameRun: " << mayBeSameRun << endlr; // construct header string TString title1 = ""; @@ -1224,11 +1224,11 @@ Bool_t FindGoodTitle(vector* ptCompareHistogramClassVector, vect if (same_Depletion && firstlabbook.depletionV >= 0) title1.Append(Form(", U_{dep}=%.1f V", firstlabbook.depletionV)); - // cout << colorred << title1 << endlr; + // cout << colorred << title1 << endlr; if (title1.Length() > 3) title1 = title1.Remove(0,2); // remove trailing " ," - headerStringsVector.push_back(title1); - // TObjArray *humanreadablestrings = trimRunnumberAtBegin(title1).Tokenize("\n"); + headerStringsVector.push_back(title1); + // TObjArray *humanreadablestrings = trimRunnumberAtBegin(title1).Tokenize("\n"); TObjArray *humanreadablestrings = title1.Tokenize("\n"); headerStringsVector.push_back(((TObjString *)(humanreadablestrings->At(0)))->String()); headerStringsVector.push_back(((TObjString *)(humanreadablestrings->At(1)))->String()); @@ -1243,34 +1243,34 @@ Bool_t FindGoodTitle(vector* ptCompareHistogramClassVector, vect if (mayBeSameRun) curhistogramp = ptCompareHistogramVector->at(histogrami); -// if (!same_RunNumber) { if (legendstr.Length()) legendstr.Append(", "); -// legendstr.Append(Form("%d", curhistogramclassp->labbook->runnumber)); } + // if (!same_RunNumber) { if (legendstr.Length()) legendstr.Append(", "); + // legendstr.Append(Form("%d", curhistogramclassp->labbook->runnumber)); } if (!same_HistType && mayBeSameRun) { if (legendstr.Length()) legendstr.Append(", "); legendstr.Append(Form("%s", trimRunnumberAtBegin(curhistogramp->GetName()).Data())); } - if (!same_Source) { if (legendstr.Length()) legendstr.Append(", "); - legendstr.Append(Form("%s", curlabbook.source.Data())); } - if (!same_ChipGen) { if (legendstr.Length()) legendstr.Append(", "); - legendstr.Append(Form("%s", curlabbook.chipGen.Data())); } - if (!same_ChipNum) { if (legendstr.Length()) legendstr.Append(", "); - legendstr.Append(Form("Chip# %s", curlabbook.chip.Data())); } - if (!same_Matrix) { if (legendstr.Length()) legendstr.Append(", "); - legendstr.Append(Form("%s, %.0fx%.0f #mum^{2} pitch, %s", curlabbook.matrix.Data(), curpixelinfo.pitchX, curpixelinfo.pitchY, curpixelinfo.comment.Data())); } - if (!same_Temp) { if (legendstr.Length()) legendstr.Append(", "); - legendstr.Append(Form("T=%.0f {}^{o}C", curlabbook.tempSens)); } - if (!same_IonRad && curlabbook.radDoseIon != 0) { if (legendstr.Length()) legendstr.Append(", "); - legendstr.Append(Form("%.1f MRad", firstlabbook.radDoseIon)); } - if (!same_NonIonRad && curlabbook.radDoseNonIon != 0) { if (legendstr.Length()) legendstr.Append(", "); - legendstr.Append(Form("%.1f*10^{13} n_{eq}/cm^{2}", curlabbook.radDoseNonIon)); } - if (!same_Clock) { if (legendstr.Length()) legendstr.Append(", "); - legendstr.Append(Form("%.2f Mhz", curlabbook.clock)); } - if (!same_Depletion && curlabbook.depletionV >= 0) { if (legendstr.Length()) legendstr.Append(", "); - legendstr.Append(Form("U_{dep}=%.1f V", curlabbook.depletionV)); } - legendStringsVector.push_back(legendstr); - -// cout << colorred << legendstr << endlr; + if (!same_Source) { if (legendstr.Length()) legendstr.Append(", "); + legendstr.Append(Form("%s", curlabbook.source.Data())); } + if (!same_ChipGen) { if (legendstr.Length()) legendstr.Append(", "); + legendstr.Append(Form("%s", curlabbook.chipGen.Data())); } + if (!same_ChipNum) { if (legendstr.Length()) legendstr.Append(", "); + legendstr.Append(Form("Chip# %s", curlabbook.chip.Data())); } + if (!same_Matrix) { if (legendstr.Length()) legendstr.Append(", "); + legendstr.Append(Form("%s, %.0fx%.0f #mum^{2} pitch, %s", curlabbook.matrix.Data(), curpixelinfo.pitchX, curpixelinfo.pitchY, curpixelinfo.comment.Data())); } + if (!same_Temp) { if (legendstr.Length()) legendstr.Append(", "); + legendstr.Append(Form("T=%.0f {}^{o}C", curlabbook.tempSens)); } + if (!same_IonRad && curlabbook.radDoseIon != 0) { if (legendstr.Length()) legendstr.Append(", "); + legendstr.Append(Form("%.1f MRad", firstlabbook.radDoseIon)); } + if (!same_NonIonRad && curlabbook.radDoseNonIon != 0) { if (legendstr.Length()) legendstr.Append(", "); + legendstr.Append(Form("%.1f*10^{13} n_{eq}/cm^{2}", curlabbook.radDoseNonIon)); } + if (!same_Clock) { if (legendstr.Length()) legendstr.Append(", "); + legendstr.Append(Form("%.2f Mhz", curlabbook.clock)); } + if (!same_Depletion && curlabbook.depletionV >= 0) { if (legendstr.Length()) legendstr.Append(", "); + legendstr.Append(Form("U_{dep}=%.1f V", curlabbook.depletionV)); } + legendStringsVector.push_back(legendstr); + + // cout << colorred << legendstr << endlr; } } // end else ptCompareHistogramClassVector->size() == 1 - + // Folder name suffix folderadd = ""; if (same_RunNumber) @@ -1286,29 +1286,29 @@ Bool_t FindGoodTitle(vector* ptCompareHistogramClassVector, vect if (!same_Matrix) for (UInt_t histogrami=0; histogrami < ptCompareHistogramClassVector->size(); histogrami++) folderadd.Append(Form(" -%s", ptCompareHistogramClassVector->at(histogrami)->labbook->matrix.Data())); - if (same_Depletion && firstlabbook.depletionV >= 0) - folderadd.Append(Form(" %.1fV", firstlabbook.depletionV)); - if (!same_Depletion && firstlabbook.depletionV >= 0) - for (UInt_t histogrami=0; histogrami < ptCompareHistogramClassVector->size(); histogrami++) - folderadd.Append(Form(" %.1fV", ptCompareHistogramClassVector->at(histogrami)->labbook->depletionV)); - if (same_Temp) - folderadd.Append(Form(" %.0fC", firstlabbook.tempSens)); - if (same_IonRad && firstlabbook.radDoseIon != 0) - folderadd.Append(Form("%.1fMRad", firstlabbook.radDoseIon)); - if (same_NonIonRad && firstlabbook.radDoseNonIon != 0) - folderadd.Append(Form(" %.1fe13neq", firstlabbook.radDoseNonIon)); - if (same_Clock) - folderadd.Append(Form(" %.2fMhz", firstlabbook.clock)); - if (!same_Clock) - for (UInt_t histogrami=0; histogrami < ptCompareHistogramClassVector->size(); histogrami++) - folderadd.Append(Form(" %.2fMhz", ptCompareHistogramClassVector->at(histogrami)->labbook->clock)); - - removeForbiddenChar(&folderadd); - - system("mkdir \""+ savepathresults + folderadd + "/\"" + " -p"); -// cout << colorred << savepathresults << folderadd << endlr; - - return false; + if (same_Depletion && firstlabbook.depletionV >= 0) + folderadd.Append(Form(" %.1fV", firstlabbook.depletionV)); + if (!same_Depletion && firstlabbook.depletionV >= 0) + for (UInt_t histogrami=0; histogrami < ptCompareHistogramClassVector->size(); histogrami++) + folderadd.Append(Form(" %.1fV", ptCompareHistogramClassVector->at(histogrami)->labbook->depletionV)); + if (same_Temp) + folderadd.Append(Form(" %.0fC", firstlabbook.tempSens)); + if (same_IonRad && firstlabbook.radDoseIon != 0) + folderadd.Append(Form("%.1fMRad", firstlabbook.radDoseIon)); + if (same_NonIonRad && firstlabbook.radDoseNonIon != 0) + folderadd.Append(Form(" %.1fe13neq", firstlabbook.radDoseNonIon)); + if (same_Clock) + folderadd.Append(Form(" %.2fMhz", firstlabbook.clock)); + if (!same_Clock) + for (UInt_t histogrami=0; histogrami < ptCompareHistogramClassVector->size(); histogrami++) + folderadd.Append(Form(" %.2fMhz", ptCompareHistogramClassVector->at(histogrami)->labbook->clock)); + + removeForbiddenChar(&folderadd); + + system("mkdir \""+ savepathresults + folderadd + "/\"" + " -p"); + // cout << colorred << savepathresults << folderadd << endlr; + + return false; } diff --git a/MABS_run_analyzer/HistogramType.c b/MABS_run_analyzer/HistogramType.c index 265ef6c..7fba882 100644 --- a/MABS_run_analyzer/HistogramType.c +++ b/MABS_run_analyzer/HistogramType.c @@ -58,13 +58,13 @@ void HistogramType::initHistograms(Int_t gotcolor, Int_t gotstyle) { initHistogramCustom(LeakageCurrentInPixel, "Leakage current per pixel" + histogramdescription, color, style, 0, cursensorinfo->columns*cursensorinfo->rows, cursensorinfo->columns*cursensorinfo->rows, "Pixel index", "Average CDS"); initHistogramCustom(LeakageCurrentInPixelSorted, "Leakage current per pixel sorted" + histogramdescription, color, style, 0, cursensorinfo->columns*cursensorinfo->rows, cursensorinfo->columns*cursensorinfo->rows, "Pixel index", "Average CDS"); initHistogramCustom(LeakageCurrentDistrib, "Average Leakage current per pixel" + histogramdescription, color, style, 0, 20, 200, "Average CDS", "Count"); -// initHistogram(LeakageCurrentInPixelSorted, "Leakage current" + histogramdescription, color, style); + // initHistogram(LeakageCurrentInPixelSorted, "Leakage current" + histogramdescription, color, style); initHistogramCustom(pixelHadGoodVeto, "Number of times pixel had good veto" + histogramdescription, color, style, 0, cursensorinfo->columns*cursensorinfo->rows, cursensorinfo->columns*cursensorinfo->rows, "Pixel index", "Counter good veto"); SeedPerc->GetXaxis()->SetRangeUser(0,50); Noise->SetBins(cursystempar->nbinsnoise, 0, cursystempar->maxbinnoise); NoiseEnd->SetBins(cursystempar->nbinsnoise, 0, cursystempar->maxbinnoise); -// LeakageCurrentInPixelSorted->SetBins(cursystempar->nbinsnoise*10, 0, cursystempar->maxbinnoise*30); + // LeakageCurrentInPixelSorted->SetBins(cursystempar->nbinsnoise*10, 0, cursystempar->maxbinnoise*30); if( labbook->chipGen.EqualTo("FSBB") || labbook->chipGen.EqualTo("Pipper2") || labbook->chipGen.EqualTo("Pegasus")) { Sum = &Sum9; @@ -127,7 +127,7 @@ void HistogramType::initHistogramArray(vector* histogramarraypointer, in histogrampointer->GetXaxis()->CenterTitle(); histogrampointer->GetYaxis()->CenterTitle(); histogrampointer->itsHistogramType = this; - + histogramarraypointer->push_back(histogrampointer); } } @@ -237,7 +237,7 @@ void HistogramType::calibrateHistogram(TH1FO* &histogrampointernew, TH1FO* &hist scalefactor = gain; for(int i=0; i <= nbins; i++){ new_bins[i] = histogrampointernew->GetBinLowEdge(i)*scalefactor; -// cout << "histogrampointernew->GetBinCenter(i): " << histogrampointernew->GetBinLowEdge(i) << " * gain " << gain << " = " << histogrampointernew->GetBinLowEdge(i)*gain << endl; + // cout << "histogrampointernew->GetBinCenter(i): " << histogrampointernew->GetBinLowEdge(i) << " * gain " << gain << " = " << histogrampointernew->GetBinLowEdge(i)*gain << endl; } histogrampointernew->SetBins(nbins, new_bins); histogrampointernew->GetYaxis()->SetTitle(Form("Entries [1/%.1f e]",histogrampointernew->GetBinWidth(1))); @@ -311,7 +311,7 @@ Bool_t HistogramType::normalizeHistograms( Int_t got_frames_found ) { if (normalized->LeakageCurrentDistrib != 0) normalized->LeakageCurrentDistrib = (TH1FO*)LeakageCurrentDistrib->Clone(); if (normalized->LeakageCurrentInPixel != 0)normalized->LeakageCurrentInPixel = (TH1FO*)LeakageCurrentInPixel->Clone(); if (normalized->LeakageCurrentInPixelSorted != 0)normalized->LeakageCurrentInPixelSorted = (TH1FO*)LeakageCurrentInPixelSorted->Clone(); - + normalized->frames_found = frames_found; normalized->posSeed = posSeed; normalized->posSum = posSum; @@ -386,7 +386,7 @@ void HistogramType::normalizeHistogramXAxis(TH1FO* &histogrampointernew, TH1FO* } histogrampointernew->SetBins(nbins, new_bins); histogrampointernew->itsHistogramType = histogrampointerold->itsHistogramType; -// histogrampointernew->GetYaxis()->SetTitle(Form("%s\b/(%d/1000000)]",histogrampointernew->GetYaxis()->GetTitle(), frames_found)); + // histogrampointernew->GetYaxis()->SetTitle(Form("%s\b/(%d/1000000)]",histogrampointernew->GetYaxis()->GetTitle(), frames_found)); } @@ -406,163 +406,125 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType else if (noisethresholdborder>=0 ) noiseborder = noisethresholdborder; -// cout << colorred << " noiseborder " << noiseborder << " : " << endlr; + // cout << colorred << " noiseborder " << noiseborder << " : " << endlr; -// cout << colorred << " " << histogrampointer->GetName() << " : " << endlr; + // cout << colorred << " " << histogrampointer->GetName() << " : " << endlr; if (fitFuncType.EqualTo("gaus")) { parameters = (Double_t *)calloc(10, sizeof(Double_t)); - 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()); -// cout << colorcyan << "Begin search for fit at: " << histogrampointer->GetXaxis()->FindBin(noiseborder) << endlr; -// cout << colorcyan << "End search for fit at: " << histogrampointer->GetXaxis()->FindBin(posMaxValHist) << endlr; + histogrampointer->GetXaxis()->UnZoom(); // resets range + Int_t binPosMaxHist = histogrampointer->FindLastBinAbove(1,1); + Int_t binPosMaxPeak = 0; + TString histogramname = histogrampointer->GetName(); + posMaxValHist = (histogrampointer->GetBinCenter(binPosMaxHist)); + Bool_t do_loop = kTRUE; + Int_t loopi = 0; + do { + histogrampointer->GetXaxis()->SetRange(histogrampointer->GetXaxis()->FindBin(noiseborder*(++loopi)),histogrampointer->GetXaxis()->FindBin(posMaxValHist)); // look only for maxima with x greater than noiseborder, cut away noise + binPosMaxPeak = histogrampointer->GetMaximumBin(); + // if Sum peak, than the peak must be in the right region of the whole histogram, so search for it + if (!histogramname.Contains("Sum") || ( binPosMaxPeak > binPosMaxHist/2) || (loopi > 3) ) { +// cout << "Exit loop, loopi: " << loopi << " binPosMaxPeak > binPosMaxHist/2: " << binPosMaxPeak << " > " << binPosMaxHist/2 << endlr; + do_loop = kFALSE; + noiseborder = noiseborder*loopi; + } + } while ( do_loop ); + if (histogramname.Contains("Sum")) { + noiseborder = (histogrampointer->GetBinCenter(binPosMaxPeak/2)); // if Sum peak is fitted, otherwise fits fails too often + } + + if (verbose) { + cout << colorcyan << "maxvalhist: " << posMaxValHist << endlr; + cout << colorcyan << "noiseborder: " << noiseborder << endlr; + } + TF1* fitFunc = new TF1("fitFunc",fitFuncType,noiseborder,posMaxValHist); - if (TString(histogrampointer->GetName()).Contains("Veto")) // not used any more - { -// if (verbose) -// cout << "Range: " << histogrampointer->GetXaxis()->FindBin(noiseborder) << " to " << histogrampointer->GetXaxis()->FindBin(posMaxValHist) << endl; - Float_t peak1 = histogrampointer->GetMaximumBin(); - Float_t peak2; -// if (verbose) -// cout << "histogrampointer->GetBinCenter(peak1) " << histogrampointer->GetBinCenter(peak1) << endl; - if (histogrampointer->GetBinCenter(peak1)<3.3*noiseborder) // vermutlich ist veto peak nicht am höchsten - { -// if (verbose) -// cout << "histogrampointer->GetBinCenter(peak1) " << histogrampointer->GetBinCenter(peak1) << " < 3 * noiseborder " << noiseborder << endl; - 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 + const Double_t def_amplitude=histogrampointer->GetBinContent(binPosMaxPeak); + const Double_t def_peakcenter=histogrampointer->GetBinCenter(binPosMaxPeak); + const Double_t def_gausssig=10; + + if (verbose) { + cout << colorcyan << "def_amplitude: " << def_amplitude << endlr; + cout << colorcyan << "def_peakcenter: " << def_peakcenter << endlr; + } + + fitFunc->SetLineWidth(4); + fitFunc->SetLineColor(kGreen); + // set start values for some parameters + fitFunc->SetParName(0,"amplitude of peak"); + fitFunc->SetParameter(0,def_amplitude); + fitFunc->SetParError(0, def_amplitude*0.01); + fitFunc->SetParName(1,"peak centroid"); + fitFunc->SetParameter(1,def_peakcenter); + fitFunc->SetParLimits(1,def_peakcenter*0.8,def_peakcenter*1.2); + fitFunc->SetParError(1, def_peakcenter*0.01); + fitFunc->SetParName(2,"Gaussian sigma"); + fitFunc->SetParameter(2,def_gausssig); + fitFunc->SetParLimits(2,0.0,100.0); + fitFunc->SetParError(2, def_gausssig*0.01); + + int fittries = 0; + Bool_t failed = false; + do { + failed = false; + fit_result_ptr = histogrampointer->Fit(fitFunc, "NMWQS", "", noiseborder, posMaxValHist); + if (gMinuit == nullptr) { + failed = true; + } 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 << colorred << " " << histogrampointer->GetName() << " : " << endlr; -// cout << coloryellow << "noisethresholdborder " << noisethresholdborder << endlr; -// 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; -// cout << coloryellow << "lastbin " << histogrampointer->FindLastBinAbove(0) << endl; -// cout << coloryellow << "posMaxValHist " << posMaxValHist << endl; -// } -// cout << colorred << gMinuit->fCstatu.Data() << endlr; -// cout << colorred << gMinuit->GetStatus() << endlr; -// cout << colorred << fitFunc->GetChisquare() << endlr; - //histogrampointer->Fit(fitFunc, "N,Q,W", "", 70, 500); - int fittries = 0; - do { -// cout << fittries << ": New range: " << histogrampointer->GetXaxis()->GetBinCenter(min+fittries*min/20) << " to " << posMaxValHist << endl; - fit_result_ptr = histogrampointer->Fit(fitFunc, "NQWS", "", histogrampointer->GetXaxis()->GetBinCenter(min+fittries*min/20), posMaxValHist); -// cout << colorred << gMinuit->fCstatu.Data() << endlr; - } while (gMinuit->fCstatu.Contains("FAILED") && fittries++ < 5); - posMax = histogrampointer->GetXaxis()->GetBinCenter(histogrampointer->GetMaximumBin()); // Methode 1 - integralPeak = histogrampointer->IntegralAndError(min+fittries*min/20, histogrampointer->GetXaxis()->FindBin(posMaxValHist), integralPeakError, "width"); -// if (verbose) -// { -// cout << coloryellow << "min " << histogrampointer->GetXaxis()->GetBinCenter( min ) << endlr; -// cout << coloryellow << "posMax " << posMax << endlr; -// } - histogrampointer->GetXaxis()->SetRange(histogrampointer->GetXaxis()->FindBin(noiseborder),histogrampointer->GetXaxis()->FindBin(posMaxValHist)); // look only for maxima with x greater than noiseborder, cut away noise -// TCanvas* canvas = new TCanvas("2121", "212121212", 900, 700); - // histogrampointer->Draw();if (verbose) - fitFunc->SetLineStyle(1); // normal for the following fits - if (verbose) - fitFunc->DrawCopy("same"); - } else { - // modified the gaussian approuch to be more powerful, escpecially the - // integral is calculated in a +- 3 sigma region. - - const Double_t def_amplitude=histogrampointer->GetBinContent(histogrampointer->GetMaximumBin()); - const Double_t def_peakcenter=histogrampointer->GetBinCenter(histogrampointer->GetMaximumBin()); -// cout << colorcyan << "def_amplitude: " << def_amplitude << endlr; -// cout << colorcyan << "def_peakcenter: " << def_peakcenter << endlr; - const Double_t def_gausssig=10; - // set start values for some parameters - fitFunc->SetParName(0,"amplitude of peak"); - fitFunc->SetParameter(0,def_amplitude); - fitFunc->SetParName(1,"peak centroid"); - fitFunc->SetParameter(1,def_peakcenter); - fitFunc->SetParLimits(1,def_peakcenter*0.8,def_peakcenter*1.2); - fitFunc->SetParName(2,"Gaussian sigma"); - fitFunc->SetParameter(2,def_gausssig); - fitFunc->SetParLimits(2,0.0,100.0); - fitFunc->SetLineWidth(4); - fitFunc->SetLineColor(kGreen); - - int fittries = 0; - Bool_t failed = false; - do { - failed = false; - fit_result_ptr = histogrampointer->Fit(fitFunc, "NMWQS", "", noiseborder, posMaxValHist); - if (gMinuit == nullptr) { + if (gMinuit->fCstatu == "NOT POSDEF" || gMinuit->fCstatu.Contains("FAILED")) { failed = true; - } else - { - if (gMinuit->fCstatu == "NOT POSDEF" || gMinuit->fCstatu.Contains("FAILED")) { - failed = true; - } - } - if (failed) - { - fitFunc->SetParameter(0,def_amplitude*(1.0-0.1*++fittries)); - fitFunc->SetParameter(1,def_peakcenter); - fitFunc->SetParameter(2,def_gausssig); - } else fittries = 100; - } while (fittries < 6); - - // get parameters - for (Int_t pari=0; pari<3; pari++) { - //cout << colorcyan << fitFunc->GetParameter(pari) << endlr; - parameters[pari]=fitFunc->GetParameter(pari); - } - if (failed || ( parameters[2]/parameters[0] > 100 ) ) + } + } + if (failed) { - fitFunc->SetParameter(0,def_amplitude); + fitFunc->SetParameter(0,def_amplitude*(1.0-0.1*++fittries)); fitFunc->SetParameter(1,def_peakcenter); fitFunc->SetParameter(2,def_gausssig); - histogrampointer->GetXaxis()->SetRange(histogrampointer->GetXaxis()->FindBin(def_peakcenter*0.7),histogrampointer->GetXaxis()->FindBin(def_peakcenter*1.3)); // look only for maxima with x greater than noiseborder, cut away noise - fit_result_ptr = histogrampointer->Fit(fitFunc, "NMWQS", "", def_peakcenter*0.7, def_peakcenter*1.3); - } - if (failed) - { - cout << colorred << " Could not find " << histogrampointer->GetName() << " peak" << endlr; - parameters = (Double_t *)calloc(8, sizeof(Double_t)); - return parameters; - } - - // get parameters - for (Int_t pari=0; pari<3; pari++) { - //cout << colorcyan << fitFunc->GetParameter(pari) << endlr; - parameters[pari]=fitFunc->GetParameter(pari); - } - // https://suchideas.com/articles/maths/applied/histogram-errors/#The_Sensible_Way - // make an error estimate, in case of rare events one can use the poisson distribution - // error bars become +- sqrt(n) for each bin, where n is the bin content. - - parameters[7] = parameters[1] - 2*parameters[2] ; // starting point of histogram integration - parameters[8] = parameters[1] + 2*parameters[2] ; // end point of histogram integration - integralPeak = histogrampointer->IntegralAndError(histogrampointer->GetXaxis()->FindBin( parameters[7] ), histogrampointer->GetXaxis()->FindBin( parameters[8]), integralPeakError, "width"); // integral value of histogram (NOT fit), normaliezed with bin size! - posMax = fitFunc->GetMaximumX(); // Methode 2 - fitFunc->SetLineStyle(1); // normal for the following fits - - if (verbose) - fitFunc->Draw("SAME"); + } else fittries = 100; + } while (fittries < 6); + + // get parameters + for (Int_t pari=0; pari<3; pari++) { + //cout << colorcyan << fitFunc->GetParameter(pari) << endlr; + parameters[pari]=fitFunc->GetParameter(pari); + } + if (failed || ( parameters[2]/parameters[0] > 100 ) ) + { + fitFunc->SetParameter(0,def_amplitude); + fitFunc->SetParameter(1,def_peakcenter); + fitFunc->SetParameter(2,def_gausssig); + histogrampointer->GetXaxis()->SetRange(histogrampointer->GetXaxis()->FindBin(def_peakcenter*0.7),histogrampointer->GetXaxis()->FindBin(def_peakcenter*1.3)); // look only for maxima with x greater than noiseborder, cut away noise + fit_result_ptr = histogrampointer->Fit(fitFunc, "NMWQS", "", def_peakcenter*0.7, def_peakcenter*1.3); + } + if (failed) + { + cout << colorred << " Could not find " << histogrampointer->GetName() << " peak" << endlr; + parameters = (Double_t *)calloc(8, sizeof(Double_t)); + return parameters; } - + + // get parameters + for (Int_t pari=0; pari<3; pari++) { + //cout << colorcyan << fitFunc->GetParameter(pari) << endlr; + parameters[pari]=fitFunc->GetParameter(pari); + } + // https://suchideas.com/articles/maths/applied/histogram-errors/#The_Sensible_Way + // make an error estimate, in case of rare events one can use the poisson distribution + // error bars become +- sqrt(n) for each bin, where n is the bin content. + + parameters[7] = parameters[1] - 2*parameters[2] ; // starting point of histogram integration + parameters[8] = parameters[1] + 2*parameters[2] ; // end point of histogram integration + integralPeak = histogrampointer->IntegralAndError(histogrampointer->GetXaxis()->FindBin( parameters[7] ), histogrampointer->GetXaxis()->FindBin( parameters[8]), integralPeakError, "width"); // integral value of histogram (NOT fit), normaliezed with bin size! + posMax = fitFunc->GetMaximumX(); // Methode 2 + fitFunc->SetLineStyle(1); // normal for the following fits + + if (verbose) + fitFunc->Draw("SAME"); + // parameters are // parameters[0]: amplitude of peak @@ -580,65 +542,65 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType parameters[9]=integralPeakError; Float_t sigma = fitFunc->GetParameter(2); posMax2 = fitFunc->GetMaximumX(); // Methode 2 -// if (verbose) -// cout << coloryellow << "posMax2 " << posMax2 << endlr; + // if (verbose) + // cout << coloryellow << "posMax2 " << posMax2 << endlr; 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; - } -// } + // 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; + } + // } parameters[1] = -1; return parameters; } 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; - } -// } + // 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; + } + // } } if (verbose) { - cout << colorcyan << "Depletion voltage: " << labbook->depletionV << endlr; + cout << colorcyan << "Depletion voltage: " << labbook->depletionV << endlr; cout << colorcyan << "Integral from bin : " << histogrampointer->FindBin(parameters[7]) << " to " << histogrampointer->GetXaxis()->FindBin(parameters[8]) << endlr; cout << colorcyan << "Integral from val : " << parameters[7] << " to " << parameters[8] << endlr; cout << colorcyan << "Integral value without bg: " << parameters[6] << endlr; cout << colorcyan << "Integral error: " << parameters[9] << endlr; - cout << colorcyan << "Gauss Sigma: " << parameters[2] << endlr; + cout << colorcyan << "Gauss Sigma: " << parameters[2] << endlr; cout << colorcyan << "Number of events : " << frames_found << endlr; (fit_result_ptr.Get())->Print("V"); -// cout << fitFunc->GetParError(0) << endlr; // retrieve the error for the parameter 0 + // cout << fitFunc->GetParError(0) << endlr; // retrieve the error for the parameter 0 } -// 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); + // 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=="landauSIMPLE") { @@ -652,7 +614,8 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType if (verbose) fitFunc->DrawCopy("same"); - } else if (fitFuncType=="landau") + } + else if (fitFuncType=="landau") { fitFuncType = "landau"; if (posVeto > 0) @@ -674,21 +637,21 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType fitFunc->SetLineStyle(1); // normal for the following fits if (verbose) fitFunc->DrawCopy("same"); -// posMax = fitFunc->GetParameter(1); // new method, just use MPV of fit, as noise border is well known. Compare with old method anyway and warn if something suspecious. -// parameters[1]=fitFunc->GetMaximumX(); + // posMax = fitFunc->GetParameter(1); // new method, just use MPV of fit, as noise border is well known. Compare with old method anyway and warn if something suspecious. + // parameters[1]=fitFunc->GetMaximumX(); posMax = fitFunc->GetMaximumX(); // new method, just use MPV of fit, as noise border is well known. Compare with old method anyway and warn if something suspecious. -// cout << "fitFunc->GetParameter(0): " << fitFunc->GetParameter(0) << endl; -// cout << "fitFunc->GetParameter(1): " << fitFunc->GetParameter(1) << endl; -// cout << "fitFunc->GetParameter(2): " << fitFunc->GetParameter(2) << endl; -// cout << "fitFunc->GetMaximumX(): " << fitFunc->GetMaximumX() << endl; -// cout << colorred << "noiseborder: " << noiseborder << endlr; + // cout << "fitFunc->GetParameter(0): " << fitFunc->GetParameter(0) << endl; + // cout << "fitFunc->GetParameter(1): " << fitFunc->GetParameter(1) << endl; + // cout << "fitFunc->GetParameter(2): " << fitFunc->GetParameter(2) << endl; + // cout << "fitFunc->GetMaximumX(): " << fitFunc->GetMaximumX() << endl; + // cout << colorred << "noiseborder: " << noiseborder << endlr; fitMax1 = fitFunc->GetMaximumX(); -// fitFunc->DrawCopy("same"); + // fitFunc->DrawCopy("same"); fit_result_ptr = histogrampointer->Fit(fitFunc, "N,Q,W,S", "", noiseborder, fitMax1*1.1); fitMax2 = fitFunc->GetMaximumX(); fitFunc->SetLineColor(kBlue); fitFunc->SetLineStyle(2); // dashed -// fitFunc->DrawCopy("same"); + // fitFunc->DrawCopy("same"); fit_result_ptr = histogrampointer->Fit(fitFunc, "N,Q,W,S", "", fitMax1*0.9, posMaxValHist); // histogrampointer->Fit(fitFunc, "N,Q,W", "", fitMax1, histogrampointer->GetBinCenter(bini)); fitMax3 = fitFunc->GetMaximumX(); @@ -708,7 +671,8 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType //fitLandauErrorLeft.push_back(posMax - minFitMax); //fitLandauErrorRight.push_back(maxFitMax - posMax); - } else if (fitFuncType=="lorentz") + } + 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); @@ -739,40 +703,52 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType parameters = (Double_t *)calloc(4, sizeof(Double_t)); 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; -// } - } else if (fitFuncType=="GaussTail") // this option is used when fitting the calibration peak + // 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; + // } + } + else if (fitFuncType=="GaussTail") // this option is used when fitting the calibration peak { + parameters = (Double_t *)calloc(10, sizeof(Double_t)); // 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; - } - + // 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))); - + Int_t binPosMaxHist = histogrampointer->FindLastBinAbove(1,1); + Int_t binPosMaxPeak = 0; + TString histogramname = histogrampointer->GetName(); + posMaxValHist = (histogrampointer->GetBinCenter(binPosMaxHist)); + Bool_t do_loop = kTRUE; + Int_t loopi = 0; + do { + histogrampointer->GetXaxis()->SetRange(histogrampointer->GetXaxis()->FindBin(noiseborder*++loopi),histogrampointer->GetXaxis()->FindBin(posMaxValHist)); // look only for maxima with x greater than noiseborder, cut away noise + binPosMaxPeak = histogrampointer->GetMaximumBin(); + // if Sum peak, than the peak must be in the right region of the whole histogram, so search for it + if (( binPosMaxPeak > binPosMaxHist/2) || (loopi > 3) ) { + do_loop = kFALSE; + noiseborder = noiseborder*loopi; + } + } while ( do_loop ); if (verbose) { cout << colorcyan << "maxvalhist: " << posMaxValHist << endlr; cout << colorcyan << "noiseborder: " << noiseborder << endlr; } - histogrampointer->GetXaxis()->SetRange(histogrampointer->GetXaxis()->FindBin(noiseborder),histogrampointer->GetXaxis()->FindBin(posMaxValHist)); // look only for maxima with x greater than noiseborder, cut away noise 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()); + const Double_t def_amplitude=histogrampointer->GetBinContent(binPosMaxPeak); + const Double_t def_peakcenter=histogrampointer->GetBinCenter(binPosMaxPeak); const Double_t def_gausssig=8.68052; const Double_t def_distgauss=20.4; const Double_t def_bgslope=0; const Double_t def_bgoffs=histogrampointer->GetBinContent(histogrampointer->FindBin((noiseborder+def_peakcenter)/2)); -// cout << colorcyan << "histogrampointer->FindBin((noiseborder+def_peakcenter)/2): " << histogrampointer->FindBin((noiseborder+def_peakcenter)/2) << endlr; + // cout << colorcyan << "histogrampointer->FindBin((noiseborder+def_peakcenter)/2): " << histogrampointer->FindBin((noiseborder+def_peakcenter)/2) << endlr; if (verbose) { cout << colorcyan << "def_amplitude: " << def_amplitude << endlr; @@ -783,32 +759,38 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType // set start values for some parameters fitFunc->SetParName(0,"amplitude of peak"); fitFunc->SetParameter(0,def_amplitude); + fitFunc->SetParError(0, def_amplitude*0.01); fitFunc->SetParName(1,"peak centroid"); fitFunc->SetParameter(1,def_peakcenter); -// fitFunc->SetParLimits(1,def_peakcenter*0.8,def_peakcenter*1.2); + fitFunc->SetParError(1, def_peakcenter*0.01); + // fitFunc->SetParLimits(1,def_peakcenter*0.8,def_peakcenter*1.2); fitFunc->SetParName(2,"Gaussian sigma"); fitFunc->SetParameter(2,def_gausssig); fitFunc->SetParLimits(2,0,150); + fitFunc->SetParError(2, def_gausssig*0.01); fitFunc->SetParName(3,"Distance from Gauss"); fitFunc->SetParameter(3,def_distgauss); fitFunc->SetParLimits(3,0,def_distgauss*4); + fitFunc->SetParError(3, def_distgauss*0.01); fitFunc->SetParName(4,"background slope"); fitFunc->SetParameter(4,def_bgslope); -// fitFunc->SetParLimits(4,def_bgslope-0.1,def_bgslope+0.1); + fitFunc->SetParError(4, def_bgslope*0.01); + // fitFunc->SetParLimits(4,def_bgslope-0.1,def_bgslope+0.1); fitFunc->SetParName(5,"background offset"); fitFunc->SetParameter(5,def_bgoffs); fitFunc->SetParLimits(5,0,def_bgoffs*4); + fitFunc->SetParError(5, def_bgoffs*0.01); // TODO: This fix disables the background -// fitFunc->FixParameter(4,0); -// fitFunc->FixParameter(5,0); + // fitFunc->FixParameter(4,0); + // fitFunc->FixParameter(5,0); int fittries = 0; Bool_t failed = false; do { failed = false; fit_result_ptr = histogrampointer->Fit(fitFunc, "N,M,W,Q,S", "", noiseborder, posMaxValHist); -// cout << colorcyan << " AFTER fit " << endlr; + // cout << colorcyan << " AFTER fit " << endlr; if (gMinuit == nullptr) { failed = true; } else @@ -828,7 +810,7 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType } else fittries = 100; -// cout << colorcyan << gMinuit->fCstatu.Data() << " fit tries: " << fittries << endlr; + // cout << colorcyan << gMinuit->fCstatu.Data() << " fit tries: " << fittries << endlr; } while (fittries < 6); fittries = 0; if (failed) @@ -869,7 +851,7 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType fitFunc->SetParameter(3,def_distgauss); fitFunc->SetParameter(4,def_bgslope); fitFunc->SetParameter(5,def_bgoffs); - + fit_result_ptr = histogrampointer->Fit(fitFunc, "N,M,W,Q,S", "", noiseborder, posMaxValHist); if (gMinuit == nullptr) { failed = true; @@ -902,16 +884,16 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType } } } - + if (failed) { cout << colorred << " Could not find " << histogrampointer->GetName() << " calibration peak" << endlr; parameters = (Double_t *)calloc(10, sizeof(Double_t)); return parameters; } -// fitFunc->FixParameter(1,parameters[1]+histogrampointer->GetBinWidth(0)); -// histogrampointer->Fit(fitFunc, "N,M,W,Q", "", noiseborder, posMaxValHist); - + // fitFunc->FixParameter(1,parameters[1]+histogrampointer->GetBinWidth(0)); + // histogrampointer->Fit(fitFunc, "N,M,W,Q", "", noiseborder, posMaxValHist); + // parameters are // parameters[0]: amplitude of peak // parameters[1]: x position of peak @@ -928,7 +910,7 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType parameters[7] = parameters[1] - 3*parameters[2] ; // starting point of histogram integration parameters[8] = parameters[1] + 3*parameters[2] ; // end point of histogram integration integralPeak = histogrampointer->IntegralAndError(histogrampointer->GetXaxis()->FindBin( parameters[7] ), histogrampointer->GetXaxis()->FindBin( parameters[8]), integralPeakError, "width"); // integral value of histogram (NOT fit), - + parameters[6]=integralPeak; parameters[9]=integralPeakError; @@ -938,46 +920,46 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType parameters[6] -= integralbg; if (verbose) { - cout << colorcyan << "Depletion voltage: " << labbook->depletionV << endlr; + cout << colorcyan << "Depletion voltage: " << labbook->depletionV << endlr; cout << colorcyan << "Integral from bin : " << histogrampointer->FindBin(parameters[7]) << " to " << histogrampointer->GetXaxis()->FindBin(parameters[8]) << endlr; cout << colorcyan << "Integral from val : " << parameters[7] << " to " << parameters[8] << endlr; cout << colorcyan << "Integral bg: " << integralbg << endlr; cout << colorcyan << "Integral value without bg: " << parameters[6] << endlr; - cout << colorcyan << "Integral value with bg: " << parameters[6] + integralbg << endlr; - cout << colorcyan << "Gauss Sigma: " << parameters[2] << endlr; + cout << colorcyan << "Integral value with bg: " << parameters[6] + integralbg << endlr; + cout << colorcyan << "Gauss Sigma: " << parameters[2] << endlr; cout << colorcyan << "Number of events : " << frames_found << endlr; - (fit_result_ptr.Get())->Print("V"); - + (fit_result_ptr.Get())->Print("V"); + } - -// DEBUG -// fitFunc->Draw("SAME"); - + + // DEBUG + // fitFunc->Draw("SAME"); + 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"); } } - - + + fit_result_giveback = fit_result_ptr.Get(); -// cout << "Correct address:" << endl; -// cout << fit_result_ptr.Get() << endl; -// -// fit_result_ptr_giveback = new TFitResultPtr(); -// -// fit_result_ptr_giveback = &fit_result_ptr; -// cout << "Gave pointer address:" << endl; -// cout << fit_result_ptr_giveback->Get() << endl; + // cout << "Correct address:" << endl; + // cout << fit_result_ptr.Get() << endl; + // + // fit_result_ptr_giveback = new TFitResultPtr(); + // + // fit_result_ptr_giveback = &fit_result_ptr; + // cout << "Gave pointer address:" << endl; + // cout << fit_result_ptr_giveback->Get() << endl; -// TFitResult resultnew = new TFitResult(fit_result_ptr.Get()); -// cout << "saved address:" << endl; -// cout << &resultnew << endl; + // TFitResult resultnew = new TFitResult(fit_result_ptr.Get()); + // cout << "saved address:" << endl; + // cout << &resultnew << endl; return parameters; } @@ -1001,12 +983,15 @@ Float_t HistogramType::FindNoisethresholdborder(TH1FO* histogrampointer, Bool_t { savetoclass = kTRUE; } -// cout << colorcyan << "Find border" << endl; + // cout << colorcyan << "Find border" << endl; Float_t posMaxValHist = histogrampointer->GetXaxis()->GetXmax(); TH1FO* smoothedcurce = (TH1FO*)histogrampointer->Clone(); // smoothedcurce->Smooth(4); // not working method, rescaling is easier and more reliable Int_t rebinningfactor = 2; +// if (labbook->source.Contains("Cd") && labbook->chipGen.Contains("Mi19")) { +// rebinningfactor *= 2.0; +// } smoothedcurce->RebinX(rebinningfactor); smoothedcurce->SetAxisRange(0,histogrampointer->GetBinCenter(histogrampointer->GetNbinsX())); @@ -1066,6 +1051,9 @@ Float_t HistogramType::FindNoisethresholdborder(TH1FO* histogrampointer, Bool_t cout << "falling at bin " << bini << ":" << smoothedcurce->GetXaxis()->GetBinCenter(bini) << " as " << curval << " < " << smoothedcurce->GetBinContent(bini) << endl; // debug } } while (rising < 3 && biniGetNbinsX()); +// if (labbook->source.Contains("Cd") && thresholdbincurcandidate < 40) { +// thresholdbincurcandidate *= 2.0; +// } thresholdbincurcandidate *= rebinningfactor; Float_t temp_noisethresholdborder = histogrampointer->GetXaxis()->GetBinUpEdge(thresholdbincurcandidate); @@ -1094,16 +1082,16 @@ Double_t HistogramType::FindBorderToPeak(TH1FO* histogrampointer, Double_t start if (verbose) { - cout << colorcyan << "max value: " << peakval << endlr; - cout << colorwhite << histogrampointer->GetName() << endlr; - TString canvastitle = Form("%s border search", histogrampointer->GetName()); - TString canvasname = Form("%s border search", histogrampointer->GetName()); -// TCanvas* canvas = new TCanvas(canvasname, canvastitle, 900, 700); -// smoothedcurce->SetLineColor(color+1); -// smoothedcurce->Draw("SAME"); -// histogrampointer->Draw("SAME"); -// canvas->Update(); - } + cout << colorcyan << "max value: " << peakval << endlr; + cout << colorwhite << histogrampointer->GetName() << endlr; + TString canvastitle = Form("%s border search", histogrampointer->GetName()); + TString canvasname = Form("%s border search", histogrampointer->GetName()); + // TCanvas* canvas = new TCanvas(canvasname, canvastitle, 900, 700); + // smoothedcurce->SetLineColor(color+1); + // smoothedcurce->Draw("SAME"); + // histogrampointer->Draw("SAME"); + // canvas->Update(); + } Int_t bini = smoothedcurce->FindBin(startVal); Int_t rising = 0; Int_t bincurcandidate = bini; @@ -1133,18 +1121,18 @@ Double_t HistogramType::FindBorderToPeak(TH1FO* histogrampointer, Double_t start } Bool_t HistogramType::integrateSr90Spectra(TH1FO* histogrampointer, Int_t frames_found, Float_t thresholdborder, Bool_t verbose) { - + Int_t thresholdbincurcandidate = 0; if (thresholdborder < 0) { if (noisethresholdborder < 0) { -// cout << "noisethresholdborder is not set, searching" << endl; + // cout << "noisethresholdborder is not set, searching" << endl; FindNoisethresholdborder(histogrampointer, true, verbose); } thresholdbincurcandidate = histogrampointer->GetXaxis()->FindBin(noisethresholdborder); -// cout << colorcyan << "noisethresholdborder: " << noisethresholdborder << endl; -// cout << colorcyan << "thresholdbincurcandidate: " << thresholdbincurcandidate << endl; + // cout << colorcyan << "noisethresholdborder: " << noisethresholdborder << endl; + // cout << colorcyan << "thresholdbincurcandidate: " << thresholdbincurcandidate << endl; } else { @@ -1174,11 +1162,11 @@ Bool_t HistogramType::integrateSr90Spectra(TH1FO* histogrampointer, Int_t frames sr90IntegralVal = histogrampointer->IntegralAndError(thresholdbincurcandidate,histogrampointer->GetNbinsX()+1, sr90IntegralErr, "width"); -// histogrampointer->Integral(min+fittries*min/20, histogrampointer->GetXaxis()->FindBin(posMaxValHist), "width"); -// for ( -// cout << "Integrate from bin " << thresholdbincurcandidate << " to " << histogrampointer->GetNbinsX() << endl; -// sr90IntegralErr /= sr90IntegralVal/100; -// cout << " Unscaled integral is " << Form("%e",sr90IntegralVal) << ", error " << sr90IntegralErr << "%" << endl; + // histogrampointer->Integral(min+fittries*min/20, histogrampointer->GetXaxis()->FindBin(posMaxValHist), "width"); + // for ( + // cout << "Integrate from bin " << thresholdbincurcandidate << " to " << histogrampointer->GetNbinsX() << endl; + // sr90IntegralErr /= sr90IntegralVal/100; + // cout << " Unscaled integral is " << Form("%e",sr90IntegralVal) << ", error " << sr90IntegralErr << "%" << endl; // cout << " "; if (verbose) cout << " "; @@ -1191,7 +1179,7 @@ Bool_t HistogramType::integrateSr90Spectra(TH1FO* histogrampointer, Int_t frames cout << "Scaled "; } if (verbose) { -// cout << "Integral is " << colorgreen << Form("%e",sr90IntegralVal) << colorreset << ", error " << sr90IntegralErr << "%" << endl; + // cout << "Integral is " << colorgreen << Form("%e",sr90IntegralVal) << colorreset << ", error " << sr90IntegralErr << "%" << endl; cout << "Integral is " << colorgreen << Form("%e",sr90IntegralVal)<< "+-" << Form("%e",sr90IntegralErr) << endlr; } @@ -1229,7 +1217,7 @@ Bool_t HistogramType::integrateSr90Spectras(Int_t frames_found, Float_t threshol } Bool_t HistogramType::calculteStoN(Bool_t verbose) { -// cout << "posSeed: " << posSeed << endl << "avgNoise: " << avgNoise << endlr; + // cout << "posSeed: " << posSeed << endl << "avgNoise: " << avgNoise << endlr; if (posSeed > 0 && avgNoise > 0) { StoN = posSeed / avgNoise; if (verbose) diff --git a/MABS_run_analyzer/Run.c b/MABS_run_analyzer/Run.c index 8fec386..6b7a33f 100644 --- a/MABS_run_analyzer/Run.c +++ b/MABS_run_analyzer/Run.c @@ -63,30 +63,30 @@ Run::Run(Int_t runnumber, Int_t loopi, TString customtitle) if (labbook.storepath.Length() > 0) { storepathRAWLinux = labbook.storepath; -// storepathRAWLinux = storepathRAWLinux.ReplaceAll("H:","/jspc53_H"); -// storepathRAWLinux = storepathRAWLinux.ReplaceAll("h:","/jspc53_H"); + // storepathRAWLinux = storepathRAWLinux.ReplaceAll("H:","/jspc53_H"); + // storepathRAWLinux = storepathRAWLinux.ReplaceAll("h:","/jspc53_H"); storepathRAWLinux = storepathRAWLinux.ReplaceAll("H:","/d/garlic"); storepathRAWLinux = storepathRAWLinux.ReplaceAll("h:","/d/garlic"); -// if (TString(hostname).EqualTo("jspc48")) -// { - storepathRAWLinux = storepathRAWLinux.ReplaceAll("U:","/d/garlic"); - storepathRAWLinux = storepathRAWLinux.ReplaceAll("u:","/d/garlic"); -// } -// else -// { -// storepathRAWLinux = storepathRAWLinux.ReplaceAll("U:","/jspc53_U"); -// storepathRAWLinux = storepathRAWLinux.ReplaceAll("u:","/jspc53_U"); -// } + // if (TString(hostname).EqualTo("jspc48")) + // { + storepathRAWLinux = storepathRAWLinux.ReplaceAll("U:","/d/garlic"); + storepathRAWLinux = storepathRAWLinux.ReplaceAll("u:","/d/garlic"); + // } + // else + // { + // storepathRAWLinux = storepathRAWLinux.ReplaceAll("U:","/jspc53_U"); + // storepathRAWLinux = storepathRAWLinux.ReplaceAll("u:","/jspc53_U"); + // } storepathRAWLinux = storepathRAWLinux.ReplaceAll("F:","/jspc12_F"); storepathRAWLinux = storepathRAWLinux.ReplaceAll("f:","/jspc12_F"); storepathRAWLinux = storepathRAWLinux.ReplaceAll("O:","/d/garlic"); storepathRAWLinux = storepathRAWLinux.ReplaceAll("o:","/d/garlic"); storepathRAWLinux = storepathRAWLinux.ReplaceAll("S:","/d/garlic"); storepathRAWLinux = storepathRAWLinux.ReplaceAll("s:","/d/garlic"); -// storepathRAWLinux = storepathRAWLinux.ReplaceAll("O:","/home/blinnik/garlic/local"); -// storepathRAWLinux = storepathRAWLinux.ReplaceAll("o:","/home/blinnik/garlic/local"); -// storepathRAWLinux = storepathRAWLinux.ReplaceAll("S:","/home/blinnik/garlic/local"); -// storepathRAWLinux = storepathRAWLinux.ReplaceAll("s:","/home/blinnik/garlic/local"); + // storepathRAWLinux = storepathRAWLinux.ReplaceAll("O:","/home/blinnik/garlic/local"); + // storepathRAWLinux = storepathRAWLinux.ReplaceAll("o:","/home/blinnik/garlic/local"); + // storepathRAWLinux = storepathRAWLinux.ReplaceAll("S:","/home/blinnik/garlic/local"); + // storepathRAWLinux = storepathRAWLinux.ReplaceAll("s:","/home/blinnik/garlic/local"); storepathRAWLinux = storepathRAWLinux.ReplaceAll("\\","/"); } else { storepathRAWLinux=Form("/d/garlic/%s/%d",labbook.chipGen.Data(),runnumber); // default empty path @@ -108,7 +108,7 @@ Run::Run(Int_t runnumber, Int_t loopi, TString customtitle) labbook.StoN = (rowsql->GetField(22) != NULL)?atof(rowsql->GetField(22)):-1; labbook.comment = (rowsql->GetField(23) != NULL)?std::string(rowsql->GetField(23)):""; labbook.depletionV = (rowsql->GetField(24) != NULL)?atof(rowsql->GetField(24)):-1; - // labbook.frames_Analyzed = (rowsql->GetField(20) != NULL)?atoi(rowsql->GetField(20)):-1; + // labbook.frames_Analyzed = (rowsql->GetField(20) != NULL)?atoi(rowsql->GetField(20)):-1; delete res; if (labbook.chipGen.Length() > 0) { @@ -144,8 +144,8 @@ Run::Run(Int_t runnumber, Int_t loopi, TString customtitle) } } } -// if (!(labbook.posVetoDB > 0) && (labbook.source != "Fe55")) // no veto peak position found for this run - if (labbook.source != "Fe55" || labbook.chipGen == "Mi19") // no veto peak position found for this run + // if (!(labbook.posVetoDB > 0) && (labbook.source != "Fe55")) // no veto peak position found for this run + if ((labbook.source != "Fe55" && labbook.source != "Cd109")|| labbook.chipGen == "Mi19") // no veto peak position found for this run { cout << colorwhite << "getVetoPeakPositionFromFe55Run():" << endlr; getVetoPeakPositionFromFe55Run(); @@ -169,7 +169,7 @@ Run::Run(Int_t runnumber, Int_t loopi, TString customtitle) // dynamical threshold cut histogramdynamicalthreshold = new HistogramType(" dynamical Threshold", &cursystemparam, &cursensorinfo, humanreadablestr, &labbook, this, rootcolors[plotStyle], rootlinestyle[plotStyle]); HistogramClassVector.push_back(histogramdynamicalthreshold); - // fixed threshold cut + // fixed threshold cut histogramfixedthreshold = new HistogramType(" fixed Threshold", &cursystemparam, &cursensorinfo, humanreadablestr, &labbook, this, rootcolors[plotStyle], rootlinestyle[plotStyle]); HistogramClassVector.push_back(histogramfixedthreshold); // RTS pixel removed histograms @@ -181,19 +181,19 @@ Run::Run(Int_t runnumber, Int_t loopi, TString customtitle) HistogramClassVector.push_back(histogramwoRTSthreshold); - histogramwoRTSAggresive = new HistogramType(" more RTS cleaned", &cursystemparam, &cursensorinfo, humanreadablestr, &labbook, this, rootcolors[plotStyle], rootlinestyle[plotStyle]); - histogramwoRTSAggresive->maskRTSpixel = true; - histogramwoRTSAggresive->RTSthreshold = 4; - HistogramClassVector.push_back(histogramwoRTSAggresive); + histogramwoRTSAggresive = new HistogramType(" more RTS cleaned", &cursystemparam, &cursensorinfo, humanreadablestr, &labbook, this, rootcolors[plotStyle], rootlinestyle[plotStyle]); + histogramwoRTSAggresive->maskRTSpixel = true; + histogramwoRTSAggresive->RTSthreshold = 4; + HistogramClassVector.push_back(histogramwoRTSAggresive); // histogramwoRTSAggresivethreshold = new HistogramType(" Threshold, more RTS cleaned", &cursystemparam, &cursensorinfo, humanreadablestr, &labbook, this, rootcolors[plotStyle], rootlinestyle[plotStyle] ); -// histogramwoRTSAggresivethreshold->maskRTSpixel = true; -// histogramwoRTSAggresivethreshold->RTSthreshold = 1.0; -// HistogramClassVector.push_back(histogramwoRTSAggresivethreshold); + // histogramwoRTSAggresivethreshold->maskRTSpixel = true; + // histogramwoRTSAggresivethreshold->RTSthreshold = 1.0; + // HistogramClassVector.push_back(histogramwoRTSAggresivethreshold); // histogram with pixel, which have a good veto spectrum // histogramGoodVeto = new HistogramType(" good Veto", &cursystemparam, &cursensorinfo, humanreadablestr, &labbook, this, rootcolors[plotStyle], rootlinestyle[plotStyle] ); -// HistogramClassVector.push_back(histogramGoodVeto); + // HistogramClassVector.push_back(histogramGoodVeto); debugDBreadout(); } @@ -254,7 +254,7 @@ Bool_t Run::debugDBreadout() void Run::setSystemSpecificParameters() { -// systemparam systemparamUSB (2800/*maxbin*/,2800/4/*nbins*/, 25/*vetothreshold*/, 10/*maxbinnoise*/, 100/*nbinsnoise*/); // TODO: uncomment and add again old binning + // systemparam systemparamUSB (2800/*maxbin*/,2800/4/*nbins*/, 25/*vetothreshold*/, 10/*maxbinnoise*/, 100/*nbinsnoise*/); // TODO: uncomment and add again old binning systemparam systemparamUSB (1400/*maxbin*/,1400/2/*nbins*/, 25/*vetothreshold*/, 10/*maxbinnoise*/, 100/*nbinsnoise*/); systemparam systemparamPegasus (2800,2800/2,20,10,100); systemparam systemparamPXI (800*16,800,75,150,150); @@ -327,7 +327,7 @@ Bool_t Run::useDynamicalNoise(Bool_t var) Bool_t Run::fillTTreeWithMoreThenInitialNoiseValues(Bool_t var) { - cout<<" Time noise analysis : " << colorwhite << (var?"1":"0") << colorreset << " <-- only used if RAW files are analyzed, force analysis to make sure" << endl; + cout<<" Time noise analysis : " << colorwhite << (var?"1":"0") << colorreset << " <-- only used if RAW files are analyzed, force analysis to make sure" << endl; processed->RefillNoiseBranch = var; return 0; } @@ -361,7 +361,7 @@ Bool_t Run::setFixedThresholdValueADU(Float_t thresholdValue) Bool_t Run::setFixedThresholdValueADU(HistogramType* HistogramTypepointer, Float_t thresholdValue) { cout<<" Fixed threshold value : " << colorwhite << thresholdValue << " ADU" << colorreset << endl; -// processed->fFixedThresholdValueInADU = thresholdValue; + // processed->fFixedThresholdValueInADU = thresholdValue; HistogramTypepointer->fixedThresholdValue = thresholdValue; return 0; } @@ -425,23 +425,23 @@ Bool_t Run::analyzeRun(Bool_t force) processed->InitialDynNoise(); int start = 0; int nframes = processed->GetNumberFrames(); -// for(int i=0; i<1000;i++) // TODO remove 100000 run 342272 + // for(int i=0; i<1000;i++) // TODO remove 100000 run 342272 for(int i=0; igetFrame(i); if (commonModeFilter) processed->filterCommonMode(); -// cout << "hitana " << i << endl; + // cout << "hitana " << i << endl; processed->hitana(); -// cout << "regetDynNoise " << endl; + // cout << "regetDynNoise " << endl; if (dynamicalNoise) processed->regetDynNoise(); progress = (Int_t)(((i-start)*100)/(nframes-1)*10); if (progress!=progress_tmp) { print_progress( (((i-start)*100.)/(nframes-1)) ); progress_tmp=progress;} } -// cout << processed->plus << " vs. " << processed->minus << " : " << (processed->plus*1.0)/processed->minus<< endl; - + // cout << processed->plus << " vs. " << processed->minus << " : " << (processed->plus*1.0)/processed->minus<< endl; + // print a dummy file to indicate, that a root file was created once if (!processed->used_default_config) { @@ -458,18 +458,18 @@ Bool_t Run::analyzeRun(Bool_t force) } if (processed->initOldRootFile()) return 1; frames_found = processed->GetNumberFrames(); -// cout << colorwhite << "plotPixSignal():"<< flush << endlr; -// processed->plotPixSignal(0,10000000,351); - + // cout << colorwhite << "plotPixSignal():"<< flush << endlr; + // processed->plotPixSignal(0,10000000,351); + // do preperations before binning seed sum and veto cout << "---------- 1. 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; if ((*curHistogramClass)->maskRTSpixel) { -// if (labbook.chipGen=="Mi19") { - cout << colorwhite << " FindRTSPixelToMask():" << endlr; - FindRTSPixelToMask(*curHistogramClass); -// } + // if (labbook.chipGen=="Mi19") { + cout << colorwhite << " FindRTSPixelToMask():" << endlr; + FindRTSPixelToMask(*curHistogramClass); + // } } } cout << "--------------------------------------------" << endl; @@ -477,34 +477,40 @@ Bool_t Run::analyzeRun(Bool_t force) binF0(); cout << colorwhite << "binSeedSumVeto():" << endlr; binSeedSumVeto(); -// cout << colorwhite << "binCluster():" << endlr; -// binCluster(); + // cout << colorwhite << "binCluster():" << endlr; + // binCluster(); cout << "---------- 2. loop over all classes --------" << endl; + HistogramType* curHistogramClassPt; for (vector::iterator curHistogramClass = HistogramClassVector.begin(); curHistogramClass != HistogramClassVector.end(); curHistogramClass++) { - cout << "Processing histograms in class: <" << colorwhite << (*curHistogramClass)->histogramdescription << colorreset << " >" << endlr; - binNoise((*curHistogramClass)); - binLeakageCurrent((*curHistogramClass)); + curHistogramClassPt = (*curHistogramClass); + cout << "Processing histograms in class: <" << colorwhite << curHistogramClassPt->histogramdescription << colorreset << " >" << endlr; + binNoise(curHistogramClassPt); + binLeakageCurrent(curHistogramClassPt); if (noisethresholdborderADU > -1) - (*curHistogramClass)->noisethresholdborder = noisethresholdborderADU; + curHistogramClassPt->noisethresholdborder = noisethresholdborderADU; if (noisethresholdborderE > -1) { if (labbook.gainDB > 0) { - (*curHistogramClass)->noisethresholdborder = noisethresholdborderE / labbook.gainDB; - cout << colorcyan << "Set noisethresholdborder to : " << (*curHistogramClass)->noisethresholdborder << " (" << noisethresholdborderE << " e)" << endlr; + curHistogramClassPt->noisethresholdborder = noisethresholdborderE / labbook.gainDB; + cout << colorcyan << "Set noisethresholdborder to : " << curHistogramClassPt->noisethresholdborder << " (" << noisethresholdborderE << " e)" << endlr; } 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(false); + // cout << colorwhite << " calculateCCE():" << endlr; + curHistogramClassPt->calculteCCE(false); } if (labbook.source.Contains("Sr90")) { //cout << colorwhite << " integrateSr90Spectra():" << endlr; - (*curHistogramClass)->integrateSr90Spectra((*curHistogramClass)->Seed, frames_found, -1, false); - (*curHistogramClass)->integrateSr90Spectras(frames_found, 58, false); - cout << colorwhite << " calculteStoN():" << endlr; - (*curHistogramClass)->calculteStoN(true); - } + curHistogramClassPt->integrateSr90Spectra(curHistogramClassPt->Seed, frames_found, -1, false); + curHistogramClassPt->integrateSr90Spectras(frames_found, 58, false); + // cout << colorwhite << " calculteStoN():" << endlr; + curHistogramClassPt->calculteStoN(true); + } + // cout << colorwhite << "normalizeHistogramClasses():" << endlr; + normalizeHistogramClasses(curHistogramClassPt); + // cout << colorwhite << "rescaleHistogramClasses():" << endlr; + rescaleHistogramClasses(curHistogramClassPt); } cout << "--------------------------------------------" << endl; if (labbook.source.Contains("Sr90")) { @@ -516,14 +522,10 @@ Bool_t Run::analyzeRun(Bool_t force) histogramdynamicalthreshold->integrateSr90Spectras(frames_found, -1, false); } } - cout << colorwhite << "normalizeHistogramClasses():" << endlr; - normalizeHistogramClasses(); - cout << colorwhite << "rescaleHistogramClasses():" << endlr; - rescaleHistogramClasses(); cout << colorwhite << "updateDatabase():" << endlr; updateDatabase(); -// cout << colorwhite << "delete MAPS class:" << endlr; + // cout << colorwhite << "delete MAPS class:" << endlr; //delete processed; return 0; } @@ -550,8 +552,8 @@ Bool_t Run::FindRTSPixelToMask(HistogramType* HistogramTypepointer, Bool_t verbo for(Int_t hiti=0; (unsigned int)hitifFrameInfo.hits;hiti++) { // loop over hits in a frame pixelnumber = processed->fFrameInfo.pixel[hiti]; binnumber = HistogramTypepointer->pixeltimefired->Fill(pixelnumber); -// cout << "binnumber: " << binnumber; // Pixel #10 will be saved in bin number 11 (!) -// cout << colorcyan << " Pixel " << pixelnumber << " fired in frame " << framei << " hiti " << hiti << endlr; + // cout << "binnumber: " << binnumber; // Pixel #10 will be saved in bin number 11 (!) + // cout << colorcyan << " Pixel " << pixelnumber << " fired in frame " << framei << " hiti " << hiti << endlr; } } @@ -562,7 +564,7 @@ Bool_t Run::FindRTSPixelToMask(HistogramType* HistogramTypepointer, Bool_t verbo HistogramTypepointer->pixeltimefiredDistrib->Fill(HistogramTypepointer->pixeltimefired->GetBinContent(pixeli)); meanpixeltimesfired += HistogramTypepointer->pixeltimefired->GetBinContent(pixeli); numberofactivepixel++; -// HistogramTypepointer->pixeltimefiredsorted->Fill(2000); + // HistogramTypepointer->pixeltimefiredsorted->Fill(2000); } } meanpixeltimesfired /= numberofactivepixel; // Very rough estimate of a mean firing time @@ -621,7 +623,7 @@ Bool_t Run::FindRTSPixelToMask(HistogramType* HistogramTypepointer, Bool_t verbo meanpixeltimesfired = meanpixeltimesfired2; } - // Estimate new standard deviation + // Estimate new standard deviation RTSpixel =false; numberofconsideredpixel=0; poscounter = 0; @@ -686,23 +688,23 @@ Bool_t Run::FindRTSPixelToMask(HistogramType* HistogramTypepointer, Bool_t verbo // cout << colorcyan << "Mean value own: " << meanpixeltimesfired << endlr; } -// for (Int_t pixeli=0; pixeli < HistogramTypepointer->pixeltimefired->GetNbinsX(); pixeli++) { -// cout << colorcyan << HistogramTypepointer->pixeltimefired->GetBinContent(pixeli)<< endlr; -// } - //if (processed->fFrameInfo.pixel[hiti] == histogramwoRTS->RTSpixel[RTSpixeli]) - - + // for (Int_t pixeli=0; pixeli < HistogramTypepointer->pixeltimefired->GetNbinsX(); pixeli++) { + // cout << colorcyan << HistogramTypepointer->pixeltimefired->GetBinContent(pixeli)<< endlr; + // } + //if (processed->fFrameInfo.pixel[hiti] == histogramwoRTS->RTSpixel[RTSpixeli]) + + HistogramTypepointer->percentageofRTSpixel = (HistogramTypepointer->RTSpixel.size()*1.0)/(numberofactivepixel*1.0)*100.0; - cout << " cutted away evertyhing with more then " << std::right<< Form("%.1f",HistogramTypepointer->RTSthreshold * stdeviation+meanpixeltimesfired2) << " hits" << endlr; + cout << " cutted away evertyhing with more then " << std::right<< Form("%.1f",HistogramTypepointer->RTSthreshold * stdeviation+meanpixeltimesfired2) << " hits" << endlr; if (HistogramTypepointer->percentageofRTSpixel > 0.1) cout << coloryellow << "Pixel with RTS noise: " << std::right<< HistogramTypepointer->RTSpixel.size() << " out of " << numberofactivepixel << std::cout.width(10) << " (" << HistogramTypepointer->percentageofRTSpixel << " %)" << endlr; else - cout << " Pixel with RTS noise: " << std::right<< HistogramTypepointer->RTSpixel.size() << " out of " << numberofactivepixel << std::cout.width(10) << " (" << HistogramTypepointer->percentageofRTSpixel << " %)" << endlr; + cout << " Pixel with RTS noise: " << std::right<< HistogramTypepointer->RTSpixel.size() << " out of " << numberofactivepixel << std::cout.width(10) << " (" << HistogramTypepointer->percentageofRTSpixel << " %)" << endlr; if ((RTSpixelHits/totalHits * 100.0) > 5) - cout << coloryellow << " Account for: " << std::right<< RTSpixelHits << " out of " << totalHits << " (" << (RTSpixelHits/totalHits * 100.0) << " %) of Hits" << endlr; + cout << coloryellow << " Account for: " << std::right<< RTSpixelHits << " out of " << totalHits << " (" << (RTSpixelHits/totalHits * 100.0) << " %) of Hits" << endlr; else - cout << " Account for: " << std::right<< RTSpixelHits << " out of " << totalHits << " (" << (RTSpixelHits/totalHits * 100.0) << " %) of Hits" << endlr; -// HistogramTypepointer->pixeltimefired->Draw(); + cout << " Account for: " << std::right<< RTSpixelHits << " out of " << totalHits << " (" << (RTSpixelHits/totalHits * 100.0) << " %) of Hits" << endlr; + // HistogramTypepointer->pixeltimefired->Draw(); return 0; } @@ -720,49 +722,49 @@ Bool_t Run::setNoisethresholdborderE(Float_t noisethresholdborder) return 0; } -Bool_t Run::rescaleHistogramClasses() +Bool_t Run::rescaleHistogramClasses(HistogramType* curHistogramClassPt) { float_t vetopeakposition = -1; - if ( histogramdynamicalthreshold->posVeto > 0 ) + Float_t gain; + + if ( curHistogramClassPt->posVeto > 0 ) { - vetopeakposition = histogramdynamicalthreshold->posVeto; - cout << " Use calibration obtained from this run, position veto: " << vetopeakposition << endlr; + vetopeakposition = curHistogramClassPt->posVeto; + cout << " Use calibration obtained from this run, position veto: " << vetopeakposition << endlr; } else if ( labbook.posVetoDB > 0 ) { vetopeakposition = labbook.posVetoDB; - cout << " Use calibration obtained from database value of this run, position veto: " << vetopeakposition << endlr; + cout << " Use calibration obtained from database value of this run, position veto: " << vetopeakposition << endlr; } else if ( Fe55run.posVeto > 0 ) { vetopeakposition = Fe55run.posVeto; - cout << " Use calibration obtained from Fe55 run in database, position veto: " << vetopeakposition << endlr; - cout << " Run number: " << Fe55run.posVetorunnumber << endl; - cout << " Temperature: " << Fe55run.temperature << endl; + cout << " Use calibration obtained from Fe55 run in database, position veto: " << vetopeakposition << endlr; + cout << " Run number: " << Fe55run.posVetorunnumber << endl; + cout << " Temperature: " << Fe55run.temperature << endl; } else { - cout << coloryellow << "Cannot rescale run from [ADU] to [e] units, no calibration peak found." << endlr; + cout << coloryellow << " Cannot rescale run from [ADU] to [e] units, no calibration peak found." << endlr; return 0; } - Float_t gain = 1; + gain = 1; if (labbook.source.Contains("Fe")) gain = 1640.0/vetopeakposition; else if (labbook.source.Contains("Cd")) - gain = 6140.0/vetopeakposition; - for (vector::iterator curHistogramClass = HistogramClassVector.begin(); curHistogramClass != HistogramClassVector.end(); curHistogramClass++) { - (*curHistogramClass)->calibrateHistograms(gain); - (*curHistogramClass)->normalized->calibrateHistograms(gain); - } + gain = 6140.0/vetopeakposition; + curHistogramClassPt->calibrateHistograms(gain); + curHistogramClassPt->normalized->calibrateHistograms(gain); + return 1; } -Bool_t Run::normalizeHistogramClasses() +Bool_t Run::normalizeHistogramClasses(HistogramType* curHistogramClassPt) { - cout << " For normalization use: " << frames_found << " frames" << endl; - for (vector::iterator curHistogramClass = HistogramClassVector.begin(); curHistogramClass != HistogramClassVector.end(); curHistogramClass++) { - (*curHistogramClass)->normalizeHistograms(frames_found); - } + if (curHistogramClassPt->histogramdescription.EqualTo("")) { + cout << " For normalization use: " << frames_found << " frames" << endl; } + curHistogramClassPt->normalizeHistograms(frames_found); return 1; } @@ -787,7 +789,7 @@ Bool_t Run::analyzeFrame(Int_t frame) if (dynamicalNoise) processed->regetDynNoise(); processed->plotFrame(frame); -// delete processed; + // delete processed; return 0; } else @@ -795,7 +797,7 @@ Bool_t Run::analyzeFrame(Int_t frame) cout << "\033[1;33mFrame number too big, max frame: " << entries << "\033[0m" << endl; } } -// delete processed; + // delete processed; return 1; } @@ -815,7 +817,7 @@ Bool_t Run::setLabel(TString newlabel) { if (newlabel.Length() == 0) newlabel = humanreadablestr; -// histogramthreshold->Seed->SetTitle(newlabel); + // histogramthreshold->Seed->SetTitle(newlabel); for (vector::iterator curHistogramClass = HistogramClassVector.begin(); curHistogramClass != HistogramClassVector.end(); curHistogramClass++) { (*curHistogramClass)->Seed->SetTitle(newlabel); (*curHistogramClass)->Veto->SetTitle(newlabel); @@ -841,8 +843,8 @@ Bool_t Run::setLabel(TString newlabel) Bool_t Run::generateReadableRunCode() { humanreadablestr = Form("%s, %s, %s, %.0fx%.0f #mum^{2} pitch, %s, T=%.0f {}^{o}C, %.1f MRad, %.1f*10^{13} n_{eq}/cm^{2}, %.2f Mhz, U_{dep}=%.1f V", labbook.source.Data(), - labbook.chipGen.Data(), labbook.matrix.Data(), curpixelinfo.pitchX, curpixelinfo.pitchY, curpixelinfo.comment.Data(), - labbook.tempSens, labbook.radDoseIon, labbook.radDoseNonIon, labbook.clock, labbook.depletionV); + labbook.chipGen.Data(), labbook.matrix.Data(), curpixelinfo.pitchX, curpixelinfo.pitchY, curpixelinfo.comment.Data(), + labbook.tempSens, labbook.radDoseIon, labbook.radDoseNonIon, labbook.clock, labbook.depletionV); //humanreadablestr = Form("%s, %s, chip %s, %s, %sT=%.1f", labbook.source.Data(), labbook.chipGen.Data(), labbook.chip.Data(), labbook.matrix.Data(), humanreadablesuffix.Data(), labbook.temp); cout << colorwhite << " " << colorgreen << humanreadablestr << colorwhite <<" " << endlr; @@ -921,7 +923,7 @@ void Run::setSubMatrixBorders(u_int x_start, u_int x_end, u_int y_start, u_int y void Run::selectSubMatrixFSBB(TString matrixname) { - + if (matrixname.EqualTo("A0")) { setSubMatrixBorders(0, cursensorinfo.columns/2-2, 0, cursensorinfo.rows, false); runcodesuffix += "_A0_"; @@ -941,7 +943,7 @@ string Run::prepareSQLStatement( T statement ) { stringstream ss; ss << statement; - cout << ss.str() << endl << endl; // for debugging +// cout << ss.str() << endl << endl; // for debugging return ss.str(); } @@ -951,7 +953,7 @@ void Run::getVetoPeakPositionFromFe55Run() { // query database and print results // ignore temperature -// TString query=prepareSQLStatement("select COALESCE(VetoPeak,-1) as VetoPeak, runnumber from radhard.labbook WHERE ChipNum='" + numberToString<>(labbook.chip) + "' AND MABS_comment IS NULL AND RadiationSource='Fe55' AND Matrix='" + numberToString<>(labbook.matrix) + "' AND TempCooling=" + to_str_w_prec(labbook.temp,3) + " AND System='" + numberToString<>(labbook.system) + "'"); + // TString query=prepareSQLStatement("select COALESCE(VetoPeak,-1) as VetoPeak, runnumber from radhard.labbook WHERE ChipNum='" + numberToString<>(labbook.chip) + "' AND MABS_comment IS NULL AND RadiationSource='Fe55' AND Matrix='" + numberToString<>(labbook.matrix) + "' AND TempCooling=" + to_str_w_prec(labbook.temp,3) + " AND System='" + numberToString<>(labbook.system) + "'"); TString query=prepareSQLStatement("select COALESCE(VetoPeak,-1) as VetoPeak, runnumber, TempCooling from radhard.labbook WHERE ChipNum='" + numberToString<>(labbook.chip) + "' AND MABS_comment IS NULL AND RadiationSource='Fe55' AND Matrix='" + numberToString<>(labbook.matrix) + "' AND System='" + numberToString<>(labbook.system) + "' ORDER BY TempCooling ASC"); res = db->Query(query.Data()); nrows = res->GetRowCount(); @@ -1005,37 +1007,37 @@ void Run::getVetoPeakPositionFromFe55Run() void Run::constructUpdateString(string *sqlupdatequery, const string databasevaluename, const Double_t value, const int precision=3, const Double_t minval=-1e100, const Double_t maxval = 1e100, const TString defaultval = "NULL") { -// cout << colorred << databasevaluename << " : " << value << endlr; + // cout << colorred << databasevaluename << " : " << value << endlr; if (std::isnan(value)) { cout << colorred << "Error in updating SQL value for '" << databasevaluename << ", it is 'NaN'." << endlr; return; } - if (std::isinf(value)) { - cout << colorred << "Error in updating SQL value for '" << databasevaluename << ", it is 'NaN'." << endlr; - return; } - if (!(abs(value)>=0.0)) { - cout << colorred << "Error in updating SQL value for '" << databasevaluename << ", it is 'NaN'." << endlr; - return; } - if (value > maxval) { - cout << colorred << "Error in updating SQL value for '" << databasevaluename << ", it is out of range: " << value << " > " << maxval << "." << endlr; - return; } - if (value < minval) { - cout << colorred << "Error in updating SQL value for '" << databasevaluename << ", it is out of range: " << value << " < " << maxval << "." << endlr; - return; } - if ((*sqlupdatequery).length() > 0) - *sqlupdatequery+= ", "; - *sqlupdatequery += "`" + databasevaluename + "`="+ to_str_w_prec(value, precision); + if (std::isinf(value)) { + cout << colorred << "Error in updating SQL value for '" << databasevaluename << ", it is 'NaN'." << endlr; + return; } + if (!(abs(value)>=0.0)) { + cout << colorred << "Error in updating SQL value for '" << databasevaluename << ", it is 'NaN'." << endlr; + return; } + if (value > maxval) { + cout << colorred << "Error in updating SQL value for '" << databasevaluename << ", it is out of range: " << value << " > " << maxval << "." << endlr; + return; } + if (value < minval) { + cout << colorred << "Error in updating SQL value for '" << databasevaluename << ", it is out of range: " << value << " < " << maxval << "." << endlr; + return; } + if ((*sqlupdatequery).length() > 0) + *sqlupdatequery+= ", "; + *sqlupdatequery += "`" + databasevaluename + "`="+ to_str_w_prec(value, precision); } void Run::updateDatabase() { - HistogramType* histogramclassToUseForDB = 0; - // Add code to set the pointer histogramclassToUseForDB to class to use for database values - // histogramclassToUseForDB = histogramthreshold; + HistogramType* histogramclassToUseForDB = 0; + // Add code to set the pointer histogramclassToUseForDB to class to use for database values + // histogramclassToUseForDB = histogramthreshold; if (labbook.chipGen.EqualTo("Pipper2")) { histogramclassToUseForDB = histogramdynamicalthreshold; - } else { - histogramclassToUseForDB = histogramwoRTSthreshold; - } - + } else { + histogramclassToUseForDB = histogramwoRTSthreshold; + } + string sqlupdatequery = ""; constructUpdateString(&sqlupdatequery, "Gain", histogramclassToUseForDB->normalized->gain, 3, 0, 100); constructUpdateString(&sqlupdatequery, "SumPeak", histogramclassToUseForDB->normalized->posSum, 4, 0, 1000); @@ -1077,24 +1079,24 @@ void Run::updateDatabase() { if (histogramclassToUseForDB->normalized->calibrated != 0) if (labbook.source.Contains("Sr") && histogramclassToUseForDB->normalized->calibrated->sr90IntegralVal > 0) constructUpdateString(&sqlupdatequery, "Sr90IntegralVal", histogramclassToUseForDB->normalized->calibrated->sr90IntegralVal, 1000000000, 0,1e7); - - if (sqlupdatequery.length()>0) - { - try + + if (sqlupdatequery.length()>0) { - sqlupdatequery = prepareSQLStatement(" UPDATE `radhard`.`labbook` SET " + sqlupdatequery + " WHERE `runnumber`=" + numberToString<>(labbook.runnumber)); - Bool_t sucess = db->Exec(sqlupdatequery.c_str()); - if (!sucess) + try + { + sqlupdatequery = prepareSQLStatement(" UPDATE `radhard`.`labbook` SET " + sqlupdatequery + " WHERE `runnumber`=" + numberToString<>(labbook.runnumber)); + Bool_t sucess = db->Exec(sqlupdatequery.c_str()); + if (!sucess) + { + cout << "\033[1;31mError while writing laboratory book TO SQL database (run number " << labbook.runnumber << ")!\033[0m\n"; + } + } + catch(...) { cout << "\033[1;31mError while writing laboratory book TO SQL database (run number " << labbook.runnumber << ")!\033[0m\n"; + //exit(EXIT_FAILURE); } } - catch(...) - { - cout << "\033[1;31mError while writing laboratory book TO SQL database (run number " << labbook.runnumber << ")!\033[0m\n"; - //exit(EXIT_FAILURE); - } - } } string Run::to_str_w_prec(const Float_t a_value, int precision = 3) @@ -1108,7 +1110,7 @@ Bool_t Run::binLeakageCurrent(HistogramType* oneHistogramClass) { Double_t const probabilities[] = {0.158655254, 0.5, 1-0.158655254}; // sigma/2 from gaus to the left and to the right //{0.17, 0.5, 1-0.17}; -// Double_t const probabilities[] = {0.25, 0.5, 0.75}; // sigma/2 from gaus to the left and to the right, an area of 68% must be within the interesting range //{0.17, 0.5, 1-0.17}; + // Double_t const probabilities[] = {0.25, 0.5, 0.75}; // sigma/2 from gaus to the left and to the right, an area of 68% must be within the interesting range //{0.17, 0.5, 1-0.17}; Double_t pedestals [cursensorinfo.columns*cursensorinfo.rows]; for (Int_t pixeli=0; pixeliLeakageCurrentDistrib->Fill(oneHistogramClass->LeakageCurrentInPixel->GetBinContent(pixeli)); sortedHistogram.push_back(oneHistogramClass->LeakageCurrentInPixel->GetBinContent(pixeli)); } -// cout << oneHistogramClass->LeakageCurrentInPixel->GetBinContent(pixeli) << " " << endl; + // cout << oneHistogramClass->LeakageCurrentInPixel->GetBinContent(pixeli) << " " << endl; } std::sort(sortedHistogram.begin(),sortedHistogram.end()); -// oneHistogramClass->LeakageCurrentInPixelSorted->SetBins(sortedHistogram.size(), 0, sortedHistogram.size()); + // oneHistogramClass->LeakageCurrentInPixelSorted->SetBins(sortedHistogram.size(), 0, sortedHistogram.size()); for (UInt_t pixeli=0; pixeli < sortedHistogram.size(); pixeli++) { oneHistogramClass->LeakageCurrentInPixelSorted->SetBinContent(pixeli+1, sortedHistogram.at(pixeli)); // be aware that, histograms in root start at 1 and vectors start at 0 ! } oneHistogramClass->avgLeakageCurrentInChip /= numberofconsideredpixel; -// cout << "avgLeakageCurrentInChip: " << oneHistogramClass->avgLeakageCurrentInChip << endl; + // cout << "avgLeakageCurrentInChip: " << oneHistogramClass->avgLeakageCurrentInChip << endl; -// oneHistogramClass->LeakageCurrentInPixelSorted->GetQuantiles( 3, leakagequantiles, probabilities); + // oneHistogramClass->LeakageCurrentInPixelSorted->GetQuantiles( 3, leakagequantiles, probabilities); Int_t lastbinabovezero = oneHistogramClass->LeakageCurrentInPixelSorted->FindLastBinAbove(0); -// cout << "lastbinabovezero: " << lastbinabovezero << endl; + // cout << "lastbinabovezero: " << lastbinabovezero << endl; oneHistogramClass->medianLeakageCurrent = oneHistogramClass->LeakageCurrentInPixelSorted->GetBinContent(lastbinabovezero*probabilities[1]); oneHistogramClass->medianLeakageCurrentPlus = oneHistogramClass->LeakageCurrentInPixelSorted->GetBinContent(lastbinabovezero*probabilities[2]); oneHistogramClass->medianLeakageCurrentPlus = oneHistogramClass->medianLeakageCurrentPlus - oneHistogramClass->medianLeakageCurrent; oneHistogramClass->medianLeakageCurrentMinus = oneHistogramClass->LeakageCurrentInPixelSorted->GetBinContent(lastbinabovezero*probabilities[0]); oneHistogramClass->medianLeakageCurrentMinus = oneHistogramClass->medianLeakageCurrentMinus - oneHistogramClass->medianLeakageCurrent; -// TCanvas *c1 = new TCanvas("c1","c1",600,400); -// oneHistogramClass->LeakageCurrentInPixelSorted->Draw(""); -// c1->Update(); -// TCanvas *c2 = new TCanvas("c2","c2",600,400); -// oneHistogramClass->LeakageCurrentDistrib->Draw(""); -// c2->Update(); -// cout << "medianLeakageCurrent: " << oneHistogramClass->medianLeakageCurrent << endl; -// cout << "medianLeakageCurrentPlus: " << oneHistogramClass->medianLeakageCurrentPlus << endl; -// cout << "medianLeakageCurrentMinus: " << oneHistogramClass->medianLeakageCurrentMinus << endl; -// -// -// oneHistogramClass->medianLeakageCurrent = leakagequantiles[1]; -// oneHistogramClass->medianLeakageCurrentPlus = leakagequantiles[2]; -// oneHistogramClass->medianLeakageCurrentMinus = leakagequantiles[0]; -// cout << "medianLeakageCurrent: " << oneHistogramClass->medianLeakageCurrent << endl; -// cout << "medianLeakageCurrentPlus: " << oneHistogramClass->medianLeakageCurrentPlus << endl; -// cout << "medianLeakageCurrentMinus: " << oneHistogramClass->medianLeakageCurrentMinus << endl; - -// cout << colorcyan << "oneHistogramClass->avgLeakageCurrentInChip: " << oneHistogramClass->avgLeakageCurrentInChip << endlr; -// cout << colorcyan << "oneHistogramClass->medianLeakageCurrent: " << oneHistogramClass->medianLeakageCurrent << endlr; - -// TH2F *hcandle = new TH2F("hcandle","Option CANDLE6 example ",40,-4,4,40,-20,20); -// Float_t px, py; -// for (Int_t pixeli=0; pixeli < oneHistogramClass->LeakageCurrentInPixel->GetNbinsX(); pixeli++) { -// hcandle->Fill(1,oneHistogramClass->LeakageCurrentInPixel->GetBinContent(pixeli)); -// } -// hcandle->SetMarkerSize(0.5); -// hcandle->SetBarWidth(1.0); -// hcandle->Draw("CANDLEX5"); -// gPad->BuildLegend(0.6,0.7,0.7,0.8); + // TCanvas *c1 = new TCanvas("c1","c1",600,400); + // oneHistogramClass->LeakageCurrentInPixelSorted->Draw(""); + // c1->Update(); + // TCanvas *c2 = new TCanvas("c2","c2",600,400); + // oneHistogramClass->LeakageCurrentDistrib->Draw(""); + // c2->Update(); + // cout << "medianLeakageCurrent: " << oneHistogramClass->medianLeakageCurrent << endl; + // cout << "medianLeakageCurrentPlus: " << oneHistogramClass->medianLeakageCurrentPlus << endl; + // cout << "medianLeakageCurrentMinus: " << oneHistogramClass->medianLeakageCurrentMinus << endl; + // + // + // oneHistogramClass->medianLeakageCurrent = leakagequantiles[1]; + // oneHistogramClass->medianLeakageCurrentPlus = leakagequantiles[2]; + // oneHistogramClass->medianLeakageCurrentMinus = leakagequantiles[0]; + // cout << "medianLeakageCurrent: " << oneHistogramClass->medianLeakageCurrent << endl; + // cout << "medianLeakageCurrentPlus: " << oneHistogramClass->medianLeakageCurrentPlus << endl; + // cout << "medianLeakageCurrentMinus: " << oneHistogramClass->medianLeakageCurrentMinus << endl; + + // cout << colorcyan << "oneHistogramClass->avgLeakageCurrentInChip: " << oneHistogramClass->avgLeakageCurrentInChip << endlr; + // cout << colorcyan << "oneHistogramClass->medianLeakageCurrent: " << oneHistogramClass->medianLeakageCurrent << endlr; + + // TH2F *hcandle = new TH2F("hcandle","Option CANDLE6 example ",40,-4,4,40,-20,20); + // Float_t px, py; + // for (Int_t pixeli=0; pixeli < oneHistogramClass->LeakageCurrentInPixel->GetNbinsX(); pixeli++) { + // hcandle->Fill(1,oneHistogramClass->LeakageCurrentInPixel->GetBinContent(pixeli)); + // } + // hcandle->SetMarkerSize(0.5); + // hcandle->SetBarWidth(1.0); + // hcandle->Draw("CANDLEX5"); + // gPad->BuildLegend(0.6,0.7,0.7,0.8); return 0; } @@ -1211,7 +1213,7 @@ Bool_t Run::binLeakageCurrent(HistogramType* oneHistogramClass) Bool_t Run::binNoise(HistogramType* oneHistogramClass) { 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}; -// cout << "processed->fNoiseTree->GetEntries: " << processed->fNoiseTree->GetEntries() << endl; + // cout << "processed->fNoiseTree->GetEntries: " << processed->fNoiseTree->GetEntries() << endl; Bool_t RTSpixel =false; u_int numberofconsideredpixel=0; @@ -1230,7 +1232,7 @@ Bool_t Run::binNoise(HistogramType* oneHistogramClass) numberofconsideredpixel++; } } - + oneHistogramClass->Noise->GetQuantiles( 3, noisequantiles, probabilities); oneHistogramClass->avgNoise = noisequantiles[1]; oneHistogramClass->avgNoisePlus = noisequantiles[2] - noisequantiles[1]; @@ -1268,7 +1270,7 @@ Bool_t Run::binF0() Double_t avgF0InCurFrame = processed->fNoiseInfo.AvgF0; averageF0Hist->SetBinContent(framei, avgF0InCurFrame); labbook.averageF0+=avgF0InCurFrame; -// labbook.averageF0 + = avgF0InCurFrame; + // labbook.averageF0 + = avgF0InCurFrame; } labbook.averageF0 /= processed->fHitTree->GetEntries(); for (Int_t framei=0; frameifHitTree->GetEntries(); framei++) // loop over all frames @@ -1327,7 +1329,7 @@ Bool_t Run::binF0() averageF1Distr->SetLineColor(4); averageF1Distr->itsHistogramType = histogram; - + Double_t const probabilities[] = {0.158655254, 0.5, 1-0.158655254}; // sigma/2 from gaus to the left and to the right //{0.17, 0.5, 1-0.17}; Double_t avgF0ForPixel [cursensorinfo.columns*cursensorinfo.rows]; Double_t avgF1ForPixel [cursensorinfo.columns*cursensorinfo.rows]; @@ -1366,7 +1368,7 @@ Bool_t Run::binF0() avgF0ForPixel[pixeli]=0; avgF1ForPixel[pixeli]=0; } - + for (Int_t framei=processed->fHitTree->GetEntries()-1000; frameifHitTree->GetEntries(); framei++) { // loop over all frames processed->getFrame(framei); for (Int_t pixeli=0; pixeliDraw(""); -// averageF0DistrEnd->Draw("SAME"); -// c1->Update(); -// -// TCanvas *c2 = new TCanvas("c2","c2",600,400); -// averageF0PixelwiseStart->Draw(""); -// averageF0PixelwiseEnd->Draw("SAME"); -// c2->Update(); + // TCanvas *c1 = new TCanvas("c1","c1",600,400); + // averageF0DistrStart->Draw(""); + // averageF0DistrEnd->Draw("SAME"); + // c1->Update(); + // + // TCanvas *c2 = new TCanvas("c2","c2",600,400); + // averageF0PixelwiseStart->Draw(""); + // averageF0PixelwiseEnd->Draw("SAME"); + // c2->Update(); return 0; @@ -1486,15 +1488,15 @@ Bool_t Run::binSeedSumVeto() } } } - - // DEBUG - // for (Int_t clusteri=0; clusterifFrameInfo.p[clusteri*5+clusterj][hiti]; - // cout << endl; - // } - // cout << "......." << endl; - + + // DEBUG + // for (Int_t clusteri=0; clusterifFrameInfo.p[clusteri*5+clusterj][hiti]; + // cout << endl; + // } + // cout << "......." << endl; + Float_t clusterArray9[9]; // temp variable clusterArray necessary, because Sort only accepts 1-dim arrays clusterArray9[0]=clusterArray25[6]; clusterArray9[1]=clusterArray25[7]; @@ -1573,72 +1575,72 @@ Bool_t Run::binSeedSumVeto() // veto spectrum 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++; - fillAHistogramsinclass(histogramthreshold, hiti, completeclustersize, pixelSum, notSeedSum, &a_pixelSum[0], &a_notSeedSum[0], pixelSum25, notSeedSum25, &a_pixelSum25[0], &a_notSeedSum25[0], pixelSum9, notSeedSum9, &a_pixelSum9[0], &a_notSeedSum9[0]); + + + if (processed->fFrameInfo.pixelthreshold[hiti]>0) + { + histogramthreshold->numberofhits++; + fillAHistogramsinclass(histogramthreshold, hiti, completeclustersize, pixelSum, notSeedSum, &a_pixelSum[0], &a_notSeedSum[0], pixelSum25, notSeedSum25, &a_pixelSum25[0], &a_notSeedSum25[0], pixelSum9, notSeedSum9, &a_pixelSum9[0], &a_notSeedSum9[0]); + + // bin the RTS cleaned histogram class + if (histogramwoRTSthreshold != 0) { + RTSpixel = false; + for (u_int RTSpixeli=0; RTSpixeli < histogramwoRTSthreshold->RTSpixel.size(); RTSpixeli++) { + if (pixelno == histogramwoRTSthreshold->RTSpixel[RTSpixeli]) + { + RTSpixel = true; + } + } + if (!RTSpixel) { + fillAHistogramsinclass(histogramwoRTSthreshold, hiti, completeclustersize, pixelSum, notSeedSum, &a_pixelSum[0], &a_notSeedSum[0], pixelSum25, notSeedSum25, &a_pixelSum25[0], &a_notSeedSum25[0], pixelSum9, notSeedSum9, &a_pixelSum9[0], &a_notSeedSum9[0]); + } + } + // bin the more agressive RTS histogram class + if (histogramwoRTSAggresivethreshold != 0) { + RTSpixel = false; + for (u_int RTSpixeli=0; RTSpixeli < histogramwoRTSAggresivethreshold->RTSpixel.size(); RTSpixeli++) { + if (pixelno == histogramwoRTSAggresivethreshold->RTSpixel[RTSpixeli]) + { + RTSpixel = true; + // cout << "not binned! RTS pixel #" << pixelno << endl; + } + } + if (!RTSpixel) { + fillAHistogramsinclass(histogramwoRTSAggresivethreshold, hiti, completeclustersize, pixelSum, notSeedSum, &a_pixelSum[0], &a_notSeedSum[0], pixelSum25, notSeedSum25, &a_pixelSum25[0], &a_notSeedSum25[0], pixelSum9, notSeedSum9, &a_pixelSum9[0], &a_notSeedSum9[0]); + } + } + } // bin the RTS cleaned histogram class - if (histogramwoRTSthreshold != 0) { + if (histogramwoRTS != 0) { RTSpixel = false; - for (u_int RTSpixeli=0; RTSpixeli < histogramwoRTSthreshold->RTSpixel.size(); RTSpixeli++) { - if (pixelno == histogramwoRTSthreshold->RTSpixel[RTSpixeli]) + for (u_int RTSpixeli=0; RTSpixeli < histogramwoRTS->RTSpixel.size(); RTSpixeli++) { + if (pixelno == histogramwoRTS->RTSpixel[RTSpixeli]) { RTSpixel = true; + break; } + } - if (!RTSpixel) { - fillAHistogramsinclass(histogramwoRTSthreshold, hiti, completeclustersize, pixelSum, notSeedSum, &a_pixelSum[0], &a_notSeedSum[0], pixelSum25, notSeedSum25, &a_pixelSum25[0], &a_notSeedSum25[0], pixelSum9, notSeedSum9, &a_pixelSum9[0], &a_notSeedSum9[0]); - } - } + if (!RTSpixel) { + fillAHistogramsinclass(histogramwoRTS, hiti, completeclustersize, pixelSum, notSeedSum, &a_pixelSum[0], &a_notSeedSum[0], pixelSum25, notSeedSum25, &a_pixelSum25[0], &a_notSeedSum25[0], pixelSum9, notSeedSum9, &a_pixelSum9[0], &a_notSeedSum9[0]); + } + } // bin the more agressive RTS histogram class - if (histogramwoRTSAggresivethreshold != 0) { + if (histogramwoRTSAggresive != 0) { RTSpixel = false; - for (u_int RTSpixeli=0; RTSpixeli < histogramwoRTSAggresivethreshold->RTSpixel.size(); RTSpixeli++) { - if (pixelno == histogramwoRTSAggresivethreshold->RTSpixel[RTSpixeli]) + for (u_int RTSpixeli=0; RTSpixeli < histogramwoRTSAggresive->RTSpixel.size(); RTSpixeli++) { + if (pixelno == histogramwoRTSAggresive->RTSpixel[RTSpixeli]) { RTSpixel = true; - // cout << "not binned! RTS pixel #" << pixelno << endl; + // cout << "not binned! RTS pixel #" << processed->fFrameInfo.pixel[hiti] << " with charge: " << processed->fFrameInfo.p[12][hiti] << endl; + // cout << "not binned! RTS pixel #" << processed->fFrameInfo.pixel[hiti] << " with charge: " << processed->fFrameInfo.p[12][hiti] << endl; } } if (!RTSpixel) { - fillAHistogramsinclass(histogramwoRTSAggresivethreshold, hiti, completeclustersize, pixelSum, notSeedSum, &a_pixelSum[0], &a_notSeedSum[0], pixelSum25, notSeedSum25, &a_pixelSum25[0], &a_notSeedSum25[0], pixelSum9, notSeedSum9, &a_pixelSum9[0], &a_notSeedSum9[0]); - } - } - - } - // bin the RTS cleaned histogram class - if (histogramwoRTS != 0) { - RTSpixel = false; - for (u_int RTSpixeli=0; RTSpixeli < histogramwoRTS->RTSpixel.size(); RTSpixeli++) { - if (pixelno == histogramwoRTS->RTSpixel[RTSpixeli]) - { - RTSpixel = true; - break; - } - + fillAHistogramsinclass(histogramwoRTSAggresive, hiti, completeclustersize, pixelSum, notSeedSum, &a_pixelSum[0], &a_notSeedSum[0], pixelSum25, notSeedSum25, &a_pixelSum25[0], &a_notSeedSum25[0], pixelSum9, notSeedSum9, &a_pixelSum9[0], &a_notSeedSum9[0]); + } } - if (!RTSpixel) { - fillAHistogramsinclass(histogramwoRTS, hiti, completeclustersize, pixelSum, notSeedSum, &a_pixelSum[0], &a_notSeedSum[0], pixelSum25, notSeedSum25, &a_pixelSum25[0], &a_notSeedSum25[0], pixelSum9, notSeedSum9, &a_pixelSum9[0], &a_notSeedSum9[0]); - } - } - // bin the more agressive RTS histogram class - if (histogramwoRTSAggresive != 0) { - RTSpixel = false; - for (u_int RTSpixeli=0; RTSpixeli < histogramwoRTSAggresive->RTSpixel.size(); RTSpixeli++) { - if (pixelno == histogramwoRTSAggresive->RTSpixel[RTSpixeli]) - { - RTSpixel = true; - // cout << "not binned! RTS pixel #" << processed->fFrameInfo.pixel[hiti] << " with charge: " << processed->fFrameInfo.p[12][hiti] << endl; - // cout << "not binned! RTS pixel #" << processed->fFrameInfo.pixel[hiti] << " with charge: " << processed->fFrameInfo.p[12][hiti] << endl; - } - } - if (!RTSpixel) { - fillAHistogramsinclass(histogramwoRTSAggresive, hiti, completeclustersize, pixelSum, notSeedSum, &a_pixelSum[0], &a_notSeedSum[0], pixelSum25, notSeedSum25, &a_pixelSum25[0], &a_notSeedSum25[0], pixelSum9, notSeedSum9, &a_pixelSum9[0], &a_notSeedSum9[0]); - } - } } else if (filli == 1) { // second fill round, some histograms are allready filled // histogramdynamicalthreshold relies on "histogram", therefore 'histogram' has to be completly filled before a threshold can be Set @@ -1657,7 +1659,7 @@ Bool_t Run::binSeedSumVeto() } } - // histogramfixedthreshold relies on "histogram", therefore 'histogram' has to be completly filled before a threshold can be Set + // histogramfixedthreshold relies on "histogram", therefore 'histogram' has to be completly filled before a threshold can be Set // bin the fixed threshold for charge histogram class if (histogramfixedthreshold != 0 && histogram != 0) { if (histogramfixedthreshold->fixedThresholdValue <= 0) { @@ -1687,7 +1689,7 @@ Bool_t Run::binSeedSumVeto() }// if good veto } } - + } else if (filli == 2) { // third fill round, some histograms are allready filled if (histogramGoodVeto != 0 && histogram != 0) {// histogramGoodVeto not used any more @@ -1712,26 +1714,28 @@ Bool_t Run::binSeedSumVeto() fillAHistogramsinclass(histogramGoodVeto, hiti, completeclustersize, pixelSum, notSeedSum, &a_pixelSum[0], &a_notSeedSum[0], pixelSum25, notSeedSum25, &a_pixelSum25[0], &a_notSeedSum25[0], pixelSum9, notSeedSum9, &a_pixelSum9[0], &a_notSeedSum9[0]); } } - + } // if fill round == 2 (third round) } // end if in range for submatrix analysis }// end loop over hits in frame } } // end loop over all frames } // end of loop of fill round counter - - + + // gROOT->SetBatch(kTRUE); for (vector::iterator curHistogramClass = HistogramClassVector.begin(); curHistogramClass != HistogramClassVector.end(); curHistogramClass++) { + cout << "Processing histograms in class: <" << colorwhite << (*curHistogramClass)->histogramdescription << colorreset << " >" << endlr; 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; - (*curHistogramClass)->FindNoisethresholdborder((*curHistogramClass)->Seed, false, false); + // cout << colorcyan << (*curHistogramClass)->histogramdescription << endlr; if (labbook.source.Contains("Fe")||labbook.source.Contains("Cd")) { + (*curHistogramClass)->FindNoisethresholdborder((*curHistogramClass)->Veto); parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->Veto, "GaussTail"); (*curHistogramClass)->posVeto = parameters[1]; (*curHistogramClass)->integralVeto = parameters[6]; } + (*curHistogramClass)->FindNoisethresholdborder((*curHistogramClass)->Seed); if (labbook.chipGen.EqualTo("Pipper2") || labbook.chipGen.EqualTo("Pegasus")) parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->Seed, "gaus", 0, false, (*curHistogramClass)->fixedThresholdValue); else @@ -1747,31 +1751,31 @@ Bool_t Run::binSeedSumVeto() for (Int_t bini = 1; bini <= (*curHistogramClass)->clustermultiplicity->GetNbinsX(); bini++) { (*curHistogramClass)->clustermultiplicity->SetBinContent(bini, (*curHistogramClass)->clustermultiplicity->GetBinContent(bini)/(*curHistogramClass)->numberofhits); } - - for (u_int sumi = 0; sumi < (*curHistogramClass)->a_Sum9.size(); sumi++) { - // Float_t integratestart = (*curHistogramClass)->FindNoisethresholdborder((*curHistogramClass)->a_Sum[sumi], false, false); - // parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->a_Sum[sumi], "gaus", 0, false, integratestart); //TODO change back to gauss - parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->a_Sum9[sumi], "gaus", 0, false, (*curHistogramClass)->fixedThresholdValue); //TODO change back to gauss - (*curHistogramClass)->a_posSum9.push_back(parameters[1]); - (*curHistogramClass)->a_integralSum9.push_back(parameters[6]); - (*curHistogramClass)->a_integralSumErr9.push_back(parameters[9]); - } - for (u_int sumi = 0; sumi < (*curHistogramClass)->a_Sum25.size(); sumi++) { - // Float_t integratestart = (*curHistogramClass)->FindNoisethresholdborder((*curHistogramClass)->a_Sum[sumi], false, false); -// parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->a_Sum[sumi], "gaus", 0, false, integratestart); //TODO change back to gauss - parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->a_Sum25[sumi], "gaus", 0, false, (*curHistogramClass)->fixedThresholdValue); //TODO change back to gauss - (*curHistogramClass)->a_posSum25.push_back(parameters[1]); - (*curHistogramClass)->a_integralSum25.push_back(parameters[6]); - (*curHistogramClass)->a_integralSumErr25.push_back(parameters[9]); - } - - parameters =(*curHistogramClass)->FitPerform((*curHistogramClass)->SeedPerc, "landau"); - (*curHistogramClass)->posSeedPerc = parameters[1]; - (*curHistogramClass)->sigmaSeedPerc = parameters[2]; - - for (Int_t bini=1; bini <= 100; bini++) { - (*curHistogramClass)->SeedPerc->SetBinContent(bini,(*curHistogramClass)->SeedPerc->GetBinContent(bini));///((*curHistogramClass)->SeedPerc->GetEntries())); - } + + for (u_int sumi = 0; sumi < (*curHistogramClass)->a_Sum9.size(); sumi++) { + // Float_t integratestart = (*curHistogramClass)->FindNoisethresholdborder((*curHistogramClass)->a_Sum[sumi], false, false); + // parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->a_Sum[sumi], "gaus", 0, false, integratestart); //TODO change back to gauss + parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->a_Sum9[sumi], "gaus", 0, false, (*curHistogramClass)->fixedThresholdValue); //TODO change back to gauss + (*curHistogramClass)->a_posSum9.push_back(parameters[1]); + (*curHistogramClass)->a_integralSum9.push_back(parameters[6]); + (*curHistogramClass)->a_integralSumErr9.push_back(parameters[9]); + } + for (u_int sumi = 0; sumi < (*curHistogramClass)->a_Sum25.size(); sumi++) { + // Float_t integratestart = (*curHistogramClass)->FindNoisethresholdborder((*curHistogramClass)->a_Sum[sumi], false, false); + // parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->a_Sum[sumi], "gaus", 0, false, integratestart); //TODO change back to gauss + parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->a_Sum25[sumi], "gaus", 0, false, (*curHistogramClass)->fixedThresholdValue); //TODO change back to gauss + (*curHistogramClass)->a_posSum25.push_back(parameters[1]); + (*curHistogramClass)->a_integralSum25.push_back(parameters[6]); + (*curHistogramClass)->a_integralSumErr25.push_back(parameters[9]); + } + + parameters =(*curHistogramClass)->FitPerform((*curHistogramClass)->SeedPerc, "landau"); + (*curHistogramClass)->posSeedPerc = parameters[1]; + (*curHistogramClass)->sigmaSeedPerc = parameters[2]; + + for (Int_t bini=1; bini <= 100; bini++) { + (*curHistogramClass)->SeedPerc->SetBinContent(bini,(*curHistogramClass)->SeedPerc->GetBinContent(bini));///((*curHistogramClass)->SeedPerc->GetEntries())); + } } if (histogramdynamicalthreshold != 0) { @@ -1795,7 +1799,7 @@ void Run::fillAHistogramsinclass(HistogramType* histogramtypepointer, Int_t hiti if (TMath::Abs(notSeedSum) < cursystemparam.vetothreshold && (labbook.source.Contains("Fe")||labbook.source.Contains("Cd"))) { histogramtypepointer->Veto->Fill(processed->fFrameInfo.p[12][hiti]); // histogram with the single pixel } - + if (pt_a_pixelSum[5]-pt_a_pixelSum[0]-pt_a_pixelSum[1]-pt_a_pixelSum[2]-pt_a_pixelSum[3]-pt_a_pixelSum[4] > labbook.NoiseAvgADC_DB*5) histogramtypepointer->clustermultiplicity->Fill(6); else if (pt_a_pixelSum[4]-pt_a_pixelSum[0]-pt_a_pixelSum[1]-pt_a_pixelSum[2]-pt_a_pixelSum[3] > labbook.NoiseAvgADC_DB*5) @@ -1868,7 +1872,7 @@ Bool_t Run::binCluster() xcoord = xcoord*TMath::Cos(rotateangle/180*TMath::Pi())-ycoord*TMath::Sin(rotateangle/180*TMath::Pi()); ycoord = xcoord_old*TMath::Sin(rotateangle/180*TMath::Pi())+ycoord*TMath::Cos(rotateangle/180*TMath::Pi()); } - // cout << colorwhite << "oben: xcoord: " << xcoord << ", ycoord: " << ycoord << endl; + // cout << colorwhite << "oben: xcoord: " << xcoord << ", ycoord: " << ycoord << endl; histogram->histAvgCluster->Fill(xcoord,ycoord,processed->fFrameInfo.p[clusteri][hiti]); if (processed->fFrameInfo.pixelthreshold[hiti]>0) @@ -2196,7 +2200,7 @@ Bool_t Run::plotAllHistograms(HistogramType* HistogramTypepointer) owntitle4->AddText(trimRunnumberAtBegin(HistogramTypepointer->Noise->GetName())); owntitle4->Draw("SAME"); -// canvas -> Print( savepathresults + "/" + canvastitle + ".eps"); + // canvas -> Print( savepathresults + "/" + canvastitle + ".eps"); canvas->Update(); TImageDump *img = new TImageDump(savepathresults + "/" + canvastitle + ".png"); @@ -2211,8 +2215,8 @@ Bool_t Run::plotAllHistograms(HistogramType* HistogramTypepointer) f->Append(img); f->Write(); -// gStyle->SetPaperSize(10.,10.); -// canvas->Print(savepathresults + "/" + canvastitle + ".tex"); + // gStyle->SetPaperSize(10.,10.); + // canvas->Print(savepathresults + "/" + canvastitle + ".tex"); return 0; } @@ -2299,7 +2303,7 @@ Bool_t Run::plotClusterDistribution(HistogramType* HistogramTypepointer) Canvas_1->Modified(); Canvas_1->cd(); Canvas_1->SetSelected(Canvas_1); - + TImageDump *img = new TImageDump(savepathresults + "/" + canvastitle + ".png"); Canvas_1->Paint(); img->Close(); @@ -2326,7 +2330,7 @@ Bool_t Run::plotClusterDistribution(HistogramType* HistogramTypepointer) 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 (HistogramTypepointer->iscalibrated) @@ -2351,13 +2355,13 @@ Bool_t Run::plotClusterDistribution(HistogramType* HistogramTypepointer) canvas->Draw(); -// TImageDump *img2 = new TImageDump(savepathresults + "/" + canvastitle + ".png"); -// canvas->Paint(); -// img2->Close(); -// -// f->Append(canvas); -// f->Append(img2); -// f->Write(); + // TImageDump *img2 = new TImageDump(savepathresults + "/" + canvastitle + ".png"); + // canvas->Paint(); + // img2->Close(); + // + // f->Append(canvas); + // f->Append(img2); + // f->Write(); // .Data(), Form("Avg. cluster distribution; x; y %s", humanreadablestr.Data()),sizearrrotx, xcoord_abs_min,xcoord_abs_max,sizearrroty,ycoord_abs_min,ycoord_abs_max); return 0; @@ -2372,11 +2376,14 @@ TCanvas* Run::plot1DHistogram(HistogramType* HistogramTypepointer, TH1FO* onehis { if (onehistogram != nullptr) { + onehistogram->GetYaxis()->UnZoom(); + onehistogram->GetXaxis()->UnZoom(); + Int_t random = random1->Rndm()*100000; TFitResult* fit_result; -// TFitResultPtr* fit_result_ptr; + // TFitResultPtr* fit_result_ptr; fit_result = new TFitResult(); -// fit_result_ptr = new TFitResultPtr(); + // fit_result_ptr = new TFitResultPtr(); TString canvastitle = Form("%s %s", onehistogram->GetName(), runcode.Data()); TString canvasname = Form("%s %s %d", onehistogram->GetName(), runcode.Data(), random); TCanvas* canvas = new TCanvas(canvasname, canvastitle, 900, 700); @@ -2399,33 +2406,34 @@ TCanvas* Run::plot1DHistogram(HistogramType* HistogramTypepointer, TH1FO* onehis } Float_t startvalue = onehistogram->GetBinLowEdge(onehistogram->FindFirstBinAbove(0)); - Float_t endrange = onehistogram->GetBinCenter(onehistogram->FindLastBinAbove(1)); + Float_t endrange = onehistogram->GetBinCenter(onehistogram->FindLastBinAbove(0)); + onehistogram->GetXaxis()->SetRange(fitstart, endrange); Float_t endrangeY = onehistogram->GetMaximum(); - onehistogram->GetXaxis()->SetRangeUser(startvalue,endrange); + onehistogram->GetXaxis()->SetRangeUser(fitstart,endrange*1.1); if (logscale) - onehistogram->GetYaxis()->SetRangeUser(0.1,endrangeY*1.4); + onehistogram->GetYaxis()->SetRangeUser(0.1,endrangeY*5); else - onehistogram->GetYaxis()->SetRangeUser(0,endrangeY*1.4); + onehistogram->GetYaxis()->SetRangeUser(0,endrangeY*1.25); TLegend* leg = new TLegend(0.2,0.85,0.89,0.89);//(0.6,0.7,0.89,0.89); leg->SetFillColor(0); leg->SetBorderSize(0); if (legendstr.Length() < 1) legendstr = Form("%s",onehistogram->GetTitle()); leg->AddEntry(onehistogram, legendstr, "l"); - + leg->SetTextSize(0.025); leg->Draw(); -// canvas->Draw(); -// canvas->Update(); - + // canvas->Draw(); + // canvas->Update(); + if (fitFuncType=="GaussTail" || fitFuncType=="gaus") { Float_t integralstart = parameters[7]; Float_t integralend = parameters[8]; Float_t integralval = parameters[6]; TString integrallbl = Form("Integral: %.0f", integralval); -// cout << colorcyan << " " << HistogramTypepointer->histogramdescription << endlr; -// cout << colorcyan << ": " << parameters[7] << endlr; + // cout << colorcyan << " " << HistogramTypepointer->histogramdescription << endlr; + // cout << colorcyan << ": " << parameters[7] << endlr; plotVerticalLineHeigher(onehistogram, integralstart); plotVerticalLineHeigher(onehistogram, integralend, "Integral to:"); plotTextAtVal(onehistogram, maxValue, integrallbl); @@ -2438,9 +2446,9 @@ TCanvas* Run::plot1DHistogram(HistogramType* HistogramTypepointer, TH1FO* onehis cout << "background slope : " << parameters[4] << endl; cout << "background offset: " << parameters[5] << endl; cout << "FWHM : " << parameters[2]*sqrt(8*log(2)) << endl; -// cout << "got address2:" << endl; -// cout << fit_result << endl; -// fit_result->Print("V"); + // cout << "got address2:" << endl; + // cout << fit_result << endl; + // fit_result->Print("V"); //TF1 *f1 = new TF1("f1","[0] +[1]*x +gaus(2)",0,5); TF1 *f1 = new TF1("f1","[0] +[1]*x",0,1000); @@ -2448,19 +2456,19 @@ TCanvas* Run::plot1DHistogram(HistogramType* HistogramTypepointer, TH1FO* onehis f1->Draw("SAME"); } } -// if (fitFuncType=="gaus") -// { -// Float_t integralstart = parameters[4]; -// Float_t integralval = parameters[3]; -// TString integrallbl = Form("Integral: %.0f", integralval); -// // cout << colorcyan << " " << HistogramTypepointer->histogramdescription << endlr; -// // cout << colorcyan << ": " << parameters[4] << endlr; -// plotVerticalLineHeigher(onehistogram, integralstart); -// plotTextAtVal(onehistogram, maxValue, integrallbl); -// } - -// canvas -> SaveAs( savepathresults + "/" + runcode + " " + onehistogram->GetName() + ".eps"); - + // if (fitFuncType=="gaus") + // { + // Float_t integralstart = parameters[4]; + // Float_t integralval = parameters[3]; + // TString integrallbl = Form("Integral: %.0f", integralval); + // // cout << colorcyan << " " << HistogramTypepointer->histogramdescription << endlr; + // // cout << colorcyan << ": " << parameters[4] << endlr; + // plotVerticalLineHeigher(onehistogram, integralstart); + // plotTextAtVal(onehistogram, maxValue, integrallbl); + // } + + // canvas -> SaveAs( savepathresults + "/" + runcode + " " + onehistogram->GetName() + ".eps"); + TImageDump *img = new TImageDump(savepathresults + "/" + runcode + " " + onehistogram->GetName() + ".png"); canvas->Paint(); img->Close(); @@ -2540,11 +2548,11 @@ Bool_t Run::writeAllHistogramsToFile() } header += Form("\n"); -// header += Form("#"); -// for (vector::iterator curHistogram = compareHistogramVector.begin(); curHistogram != compareHistogramVector.end(); curHistogram++) { -// header += Form("bin\t%s\t",(*curHistogram)->GetName()); -// } -// *fout << header << endl; + // header += Form("#"); + // for (vector::iterator curHistogram = compareHistogramVector.begin(); curHistogram != compareHistogramVector.end(); curHistogram++) { + // header += Form("bin\t%s\t",(*curHistogram)->GetName()); + // } + // *fout << header << endl; header += Form("#posSeed, run: %.1f, DB: %.1f\n", histogram->posSeed, labbook.posSeedDB); header += Form("#posSum, run: %.1f, DB: %.1f\n", histogram->posSum, labbook.posSumDB); diff --git a/MABS_run_analyzer/Run.h b/MABS_run_analyzer/Run.h index b1d811f..9f4c9e3 100644 --- a/MABS_run_analyzer/Run.h +++ b/MABS_run_analyzer/Run.h @@ -212,9 +212,9 @@ private: /** @brief overwrite value from which the integral should be taken, overwrite noise cut in ADU */ Float_t noisethresholdborderADU = -1; /** @brief Uses a found calibration peak value to rescale all the different cuts and histograms of type #HistogramType from ADU to e */ - Bool_t rescaleHistogramClasses(); + Bool_t rescaleHistogramClasses(HistogramType* oneHistogramClass); /** @brief Uses the number of frames found (#frames_found) for calculation to normalize the histograms of type #HistogramType */ - Bool_t normalizeHistogramClasses(); + Bool_t normalizeHistogramClasses(HistogramType* oneHistogramClass); /** @brief finds and and masks RTS pixel if maskRTSpixel isset to true */ Bool_t FindRTSPixelToMask(HistogramType* histogramtypepointer, Bool_t verbose = kFALSE); diff --git a/MABS_run_analyzer/help.h b/MABS_run_analyzer/help.h index 910038d..bdfd649 100644 --- a/MABS_run_analyzer/help.h +++ b/MABS_run_analyzer/help.h @@ -196,14 +196,14 @@ struct sensorinfostruct { }; /// 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 colorwhite = "\033[1;29m"; +const TString colorred = "\033[1;31m"; const TString coloryellow = "\033[1;33m"; -const TString colorgreen = "\033[1;32m"; +const TString colorgreen = "\033[1;32m"; const TString colorpurple = "\033[1;35m"; -const TString colorcyan = "\033[1;36m"; -const TString colorreset = "\033[0m"; -const TString endlr = "\033[0m\n"; +const TString colorcyan = "\033[1;36m"; +const TString colorreset = "\033[0m"; +const TString endlr = "\033[0m\n"; Int_t* rootcolors = new Int_t[13]{1, 2, 4, 6, 8, 13, 46, 28, 32, 33, 12, 20, 40}; /** -- 2.43.0