From 66b52be98dc6a62a06a149c2efe1009b698ceb9e Mon Sep 17 00:00:00 2001 From: Benjamin Linnik Date: Thu, 31 Aug 2017 15:49:40 +0200 Subject: [PATCH] Analyzer: RTS routine bug fix, update database improvements --- MABS_run_analyzer/ChargeSpektrum.c | 147 ++------ MABS_run_analyzer/ChargeSpektrumFunctions.c | 112 +++++- MABS_run_analyzer/HistogramType.c | 135 +++++-- MABS_run_analyzer/HistogramType.h | 74 ++-- MABS_run_analyzer/MAPS.c | 54 ++- MABS_run_analyzer/Run.c | 397 +++++++++++--------- MABS_run_analyzer/Run.h | 18 +- MABS_run_analyzer/help.h | 1 + 8 files changed, 523 insertions(+), 415 deletions(-) diff --git a/MABS_run_analyzer/ChargeSpektrum.c b/MABS_run_analyzer/ChargeSpektrum.c index fbb844f..36350d0 100644 --- a/MABS_run_analyzer/ChargeSpektrum.c +++ b/MABS_run_analyzer/ChargeSpektrum.c @@ -59,17 +59,18 @@ void ChargeSpektrum(TString runnumber = "") // runs[runi]->setSubMatrixBorders(0,100,0,100); // if (runi == 2) // runs[runi]->setSubMatrixBorders(100,120,100,120); - //runs[runi]->setFixedThresholdValueElectrons(1440); - //runs[runi]->setNoisethresholdborderADU(30); + //runs[runi]->setFixedThresholdValueElectrons(1440); + runs[runi]->setFixedThresholdValueADU(runs[runi]->labbook.posSeedDB*0.7); + //runs[runi]->setNoisethresholdborderADU(800); //runs[runi]->setNoisethresholdborderE(220); // runs[runi]->analyzeFrame(6384); // runs[runi]->analyzeFrame(6605); -// runs[runi]->analyzeFrame(675); +// runs[runi]->analyzeFrame(1994914); // runs[runi]->analyzeFrame(680); // runs[runi]->analyzeFrame(685); // runs[runi]->analyzeFrame(690); -// continue; +// continue; // creates or opens .root file, can analyze the RAW data @@ -102,134 +103,32 @@ void ChargeSpektrum(TString runnumber = "") gROOT->SetBatch(kFALSE); // Uncomment below to do classical analysis withour RTS -// compareHistogramClassVector.push_back(runs[runi]->histogramwoRTS->normalized); + compareHistogramClassVector.push_back(runs[runi]->histogramwoRTS->normalized); // compareHistogramClassVector2.push_back(runs[runi]->histogramwoRTS); // compareHistogramClassVector.push_back(runs[runi]->histogram->normalized); // compareHistogramClassVector2.push_back(runs[runi]->histogramwoRTS->normalized); - compareHistogramClassVector.push_back(runs[runi]->histogramwoRTS->normalized); -// runs[runi]->plot1DHistogram(runs[runi]->histogram, runs[runi]->histogramwoRTS->pixeltimefired, "", 0); - //compareHistogramClassVector.push_back(runs[runi]->histogramwoRTS); -// compareHistogramClassVector2.push_back(runs[runi]->histogramthreshold); - //compareHistogramVector.push_back(runs[runi]->histogram->normalized->calibrated->Seed); - //compareHistogramVector2.push_back(runs[runi]->histogram->normalized->Seed); - -// runs[runi]->plot1DHistogram(runs[runi]->histogram->normalized->calibrated, runs[runi]->processed->fNoiseInfo.fPedestals, "", 1); -// runs[runi]->plot1DHistogram(runs[runi]->histogram->normalized, runs[runi]->histogramwoRTS->normalized->calibrated->Seed, "landau", 1, 1, 0, 200); -// runs[runi]->plot1DHistogram(runs[runi]->histogram, runs[runi]->histogramwoRTS->pixeltimefiredsorted, "", 0); - -// Integralberechnung vom Veto Peak -// runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS, runs[runi]->histogramwoRTS->Seed, "GaussTail", true); -// runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS, runs[runi]->histogramwoRTS->Sum, "GaussTail", true); - - // compareHistogramVector.push_back(runs[runi]->histogram->normalized->Seed); -// compareHistogramVector.push_back(runs[runi]->histogram->normalized->Seed); -// compareHistogramVector.push_back(runs[runi]->histogramthreshold->normalized->Seed); -// compareHistogramVector.push_back(runs[runi]->histogramwoRTS->normalized->Seed); -// compareHistogramVector.push_back(runs[runi]->histogramwoRTS->normalized->Sum); -// compareHistogramVector.push_back(runs[runi]->histogramwoRTSAggresive->normalized->Seed); - - runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS, runs[runi]->histogramwoRTS->Seed, "GaussTail", true, true, false, 500); - runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS, runs[runi]->histogramwoRTS->Sum, "GaussTail", true, true, false, 500); - -// compareHistogramVector2.push_back(runs[runi]->histogramwoRTS->normalized->Seed); -// compareHistogramVector3.push_back(runs[runi]->histogramwoRTS->normalized->Sum); - -// compareHistogramVector3.push_back(runs[runi]->histogram->pixeltimefired); -// compareHistogramVector3.push_back(runs[runi]->histogramwoRTS->pixeltimefired); -// compareHistogramVector3.push_back(runs[runi]->histogramwoRTSAggresive->pixeltimefired); - // - // compareHistogramVector5.push_back(runs[runi]->histogram->pixeltimefiredsorted); - -// compareHistogramVector2.push_back(runs[runi]->histogramwoRTS->Seed); -// compareHistogramVector3.push_back(runs[runi]->histogramwoRTS->Sum); - -// runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS->normalized->calibrated, runs[runi]->histogramwoRTS->normalized->calibrated->Sum, "gaus", true); - -// // Uncomment below to do analysis without RTS pixel -// compareHistogramClassVector2.push_back(runs[runi]->histogramwoRTSthreshold); -// // compareHistogramVector2.push_back(runs[runi]->histogramwoRTSthreshold->Seed); -// compareHistogramClassVector3.push_back(runs[runi]->histogramwoRTSthreshold->calibrated); -// compareHistogramVector4.push_back(runs[runi]->histogramwoRTSthreshold->calibrated->Seed); -// -// // RTS analysis -// runs[runi]->compareHistogramClassVector.push_back(runs[runi]->histogramthreshold); -// runs[runi]->compareHistogramClassVector.push_back(runs[runi]->histogramwoRTSthreshold); -// // runs[runi]->compareHistogramClassVector.push_back(runs[runi]->histogramwoRTSAggresivethreshold); -// plotAllRuns(&runs[runi]->compareHistogramClassVector); -// runs[runi]->compareHistogramClassVector.clear(); -// -// // RTS Hit analysis -// runs[runi]->plot1DHistogram(runs[runi]->histogram, runs[runi]->histogramwoRTSthreshold->pixeltimefired, "", 0); -// runs[runi]->plot1DHistogram(runs[runi]->histogram, runs[runi]->histogramwoRTSthreshold->pixeltimefiredsorted, "", 0); -// // compareHistogramVector3.push_back(runs[runi]->histogramwoRTSthreshold->pixeltimefiredsorted); -// compareHistogramVector6.push_back(runs[runi]->histogramwoRTSthreshold->calibrated->pixeltimefiredsorted); - - // Leakage current -// runs[runi]->plot1DHistogram(runs[runi]->histogram, runs[runi]->histogram->calibrated->LeakageCurrentInPixel, "", 0); -// runs[runi]->plot1DHistogram(runs[runi]->histogram, runs[runi]->histogramthreshold->LeakageCurrentInPixelSorted, "", 0); -// compareHistogramVector5.push_back(runs[runi]->histogram->calibrated->LeakageCurrentInPixelSorted); -// compareHistogramClassVector4.push_back(runs[runi]->histogramwoRTS->calibrated); -// compareHistogramClassVector5.push_back(runs[runi]->histogramwoRTS); -// compareHistogramVector7.push_back(runs[runi]->histogram->calibrated->LeakageCurrentInPixelSorted); - - +// compareHistogramClassVector.push_back(runs[runi]->histogramwoRTS->normalized); +// compareHistogramClassVector.push_back(runs[runi]->histogramfixedthreshold->normalized); + // compareHistogramClassVector.push_back(runs[runi]->histogramfixedthreshold->normalized); +// runs[runi]->plot1DHistogram(runs[runi]->histogram->Seed); +// runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS->Seed); +// compareHistogramVector.push_back(runs[runi]->histogram->Seed); +// compareHistogramVector.push_back(runs[runi]->histogramwoRTS->Seed); +// compareHistogramVector.push_back(runs[runi]->histogramwoRTSAggresive->Seed); +// compareHistogramVector.push_back(runs[runi]->histogramfixedthreshold->Seed); +// compareHistogramVector.push_back(runs[runi]->histogramfixedthreshold->normalized->Seed); +// runs[runi]->plot1DHistogram( runs[runi]->histogramfixedthreshold->normalized->Seed, "gaus", true, false, false, 500); + runs[runi]->plot1DHistogram( runs[runi]->histogramwoRTS->normalized->Seed, "GaussTail"); +// runs[runi]->plot1DHistogram( runs[runi]->histogramfixedthreshold->normalized->Seed, "gaus", true, false, false, runs[runi]->histogramfixedthreshold->normalized->fixedThresholdValue); +// runs[runi]->plot1DHistogram( runs[runi]->histogramfixedthreshold->normalized->Sum, "gaus", true, false, false, 500); -// runs[runi]->plot1DHistogram(runs[runi]->histogram,runs[runi]->histogram->Seed,"landau","Test"); -// runs[runi]->compareHistogramClassVector.push_back(runs[runi]->histogram->calibrated); -// runs[runi]->compareHistogramClassVector.push_back(runs[runi]->histogramthreshold->calibrated); -// runs[runi]->compareHistogramClassVector.push_back(runs[runi]->histogramfixedthreshold); -// runs[runi]->plotCompareHistograms(); -// runs[runi]->plotAllHistograms(runs[runi]->histogram); -// runs[runi]->plotAllHistograms(runs[runi]->histogramwoRTSthreshold); -// -// compareHistogramClassVector.push_back(runs[runi]->histogramwoRTSAggresivethreshold); -// compareHistogramClassVector.push_back(runs[runi]->histogramwoRTSthreshold); -// compareHistogramClassVector.push_back(runs[runi]->histogramwoRTSAggresivethreshold); -// compareHistogramClassVector.push_back(runs[runi]->histogramGoodVeto); + // plot RTS pixel histogram +// runs[runi]->plot1DHistogram( runs[runi]->histogramwoRTS->pixeltimefiredDistrib, "", 0); -// compareHistogramClassVector2.push_back(runs[runi]->histogramGoodVeto); -// -// compareHistogramVector.push_back(runs[runi]->histogram->Seed); -// compareHistogramVector.push_back(runs[runi]->histogramwoRTSthreshold->Seed); -// compareHistogramVector.push_back(runs[runi]->histogramwoRTSAggresivethreshold->Seed); -// compareHistogramVector2.push_back(runs[runi]->histogramwoRTSthreshold->Seed); -// -// // compareHistogramClassVector.push_back(runs[runi]->histogram); -// // compareHistogramVector2.push_back(runs[runi]->histogram->LeakageCurrentInPixel); -// -// compareHistogramVector5.push_back(runs[runi]->histogram->pixeltimefiredsorted); -// compareHistogramVector6.push_back(runs[runi]->histogramwoRTS->pixeltimefiredsorted); + -// compareHistogramVector.push_back(runs[runi]->histogramwoRTSthreshold->normalized->Seed); -// compareHistogramVector2.push_back(runs[runi]->histogramwoRTS->normalized->Noise); -// compareHistogramVector3.push_back(runs[runi]->histogramwoRTS->pixeltimefiredsorted); -// compareHistogramVector.push_back(runs[runi]->histogramwoRTS->normalized->Seed); -// runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS, runs[runi]->histogramwoRTS->pixeltimefired, "", 0); -// runs[runi]->plot1DHistogram(runs[runi]->histogram, runs[runi]->histogram->Seed, "landau", 1); -// runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS, runs[runi]->histogramwoRTS->pixeltimefiredsorted, "", 0, true, true); -// runs[runi]->plot1DHistogram(runs[runi]->histogram, runs[runi]->histogram->pixeltimefired, "", 0, true, true); -// runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS, runs[runi]->histogramwoRTS->pixeltimefired, "", 0, true, true); -// runs[runi]->plot1DHistogram(runs[runi]->histogram, runs[runi]->histogram->pixeltimefired, "", 0, true, true); -// runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS, runs[runi]->histogramwoRTS->LeakageCurrentInPixel, "", 0); -// runs[runi]->plot1DHistogram(runs[runi]->histogram, runs[runi]->histogram->LeakageCurrentInPixel, "", 0); -// runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS, runs[runi]->histogramwoRTS->LeakageCurrentInPixelSorted, "", 0); -// runs[runi]->plot1DHistogram(runs[runi]->histogram, runs[runi]->histogram->LeakageCurrentInPixelSorted, "", 0); - -// runs[runi]->plot1DHistogram(runs[runi]->histogram, runs[runi]->histogram->pixeltimefired, "gaus", 1); -// compareHistogramVector2.push_back(runs[runi]->histogramthreshold->normalized->Seed); - //compareHistogramClassVector.push_back(runs[runi]->histogram); - // runs[runi]->plot1DHistogram(runs[runi]->histogramthreshold, runs[runi]->histogramthreshold->Veto, "gaus"); -// runs[runi]->plot1DHistogram(runs[runi]->histogram, runs[runi]->histogram->Sum, "gaus", 1); -// runs[runi]->plot1DHistogram(runs[runi]->histogram, runs[runi]->histogram->normalized->Seed, "landau", 1); - //runs[runi]->plot1DHistogram(runs[runi]->histogram, runs[runi]->histogram->normalized->Seed, "GaussTail", 1, 0,0,70); - // compareHistogramVector.push_back(runs[runi]->histogramthreshold->calibrated->normalized->Veto); -// runs[runi]->plot1DHistogram(runs[runi]->histogram, runs[runi]->histogram->Seed, "landau", 1); runs[runi]->writeAllHistogramsToFile(); -// CompareHistograms(&compareHistogramVector); -// compareHistogramVector.clear(); - //plotAllRuns(&compareHistogramClassVector); -// compareHistogramClassVector.clear(); } //cout << runs[runi]->histogram } diff --git a/MABS_run_analyzer/ChargeSpektrumFunctions.c b/MABS_run_analyzer/ChargeSpektrumFunctions.c index 909ad0e..45fc2d6 100644 --- a/MABS_run_analyzer/ChargeSpektrumFunctions.c +++ b/MABS_run_analyzer/ChargeSpektrumFunctions.c @@ -1,7 +1,7 @@ /** * @file ChargeSpektrumFunctions.c * @brief This file contains functions used in ChargeSpectrum - * + * * It is used, to make #ChargeSpectrum.c shorter and easier to maintain * Generally speaking, all type of helper functions are here. For example for comparing plots of different Runs or write ourt summary tables. * @@ -40,7 +40,8 @@ TString folderadd; * If findind a common nominator fails it return a false value * */ -Bool_t FindGoodTitle(vector* ptCompareHistogramClassVector); +Bool_t FindGoodTitle(vector* ptCompareHistogramClassVector, vector* ptCompareHistogramVector = 0); +Bool_t FindGoodTitle(vector* ptCompareHistogramVector); Bool_t FindGoodTitle(); @@ -60,7 +61,7 @@ template string printTableElement(varType t1, varType t2, cons * You must push TH1F pointers to #compareHistogramVector to be able to use this * */ -Bool_t CompareHistograms(vector* ptCompareHistogramVector, TString titlestr = "", TString YAxisTitle = ""); +Bool_t CompareHistograms(vector* ptCompareHistogramVector, TString titlestr = "", TString YAxisTitle = ""); Bool_t writeOneHistogramTypeToFile(vector*); @@ -86,7 +87,7 @@ Bool_t testifMixingCalibration(vector*); /** @brief Turns a value into a string with fixed precision */ string to_str_w_prec(const Float_t a_value, int precision = 1); vector compareHistogramClassVector, compareHistogramClassVector2, compareHistogramClassVector3, compareHistogramClassVector4, compareHistogramClassVector5, compareHistogramClassVector6, compareHistogramClassVector7, compareHistogramClassVector8, compareHistogramClassVectorMABSBot; -vector compareHistogramVector, compareHistogramVector2, compareHistogramVector3, compareHistogramVector4, compareHistogramVector5, compareHistogramVector6, compareHistogramVector7, compareHistogramVector8; +vector compareHistogramVector, compareHistogramVector2, compareHistogramVector3, compareHistogramVector4, compareHistogramVector5, compareHistogramVector6, compareHistogramVector7, compareHistogramVector8; TString ownpath = ""; TString savepathresults; @@ -320,7 +321,7 @@ Bool_t writeObservableToFile() filename+=".dat"; fstream* fout = new fstream(filename,ios::out); header += Form("#"); - for (vector::iterator curHistogram = compareHistogramVector.begin(); curHistogram != compareHistogramVector.end(); curHistogram++) { + for (vector::iterator curHistogram = compareHistogramVector.begin(); curHistogram != compareHistogramVector.end(); curHistogram++) { header += Form("bin\t%s\t",(*curHistogram)->GetName()); } *fout << header << endl; @@ -328,7 +329,7 @@ Bool_t writeObservableToFile() for(Int_t bin=1;bin<=runs[0]->cursystemparam.nbins;bin++) { outline =""; - for (vector::iterator curHistogram = compareHistogramVector.begin(); curHistogram != compareHistogramVector.end(); curHistogram++) { + for (vector::iterator curHistogram = compareHistogramVector.begin(); curHistogram != compareHistogramVector.end(); curHistogram++) { outline+=Form("%.1f\t%.1f\t", (*curHistogram)->GetBinCenter(bin), (*curHistogram)->GetBinContent(bin)); } *fout<* ptCompareHistogramVector, TString titlestr, TString YAxisTitle) +Bool_t CompareHistograms(vector* ptCompareHistogramVector, TString titlestr, TString YAxisTitle) { if (ptCompareHistogramVector->size() > 0) { + FindGoodTitle(ptCompareHistogramVector); + gROOT->SetStyle("RadHard_NoTitle"); // legend entries @@ -395,17 +398,28 @@ Bool_t CompareHistograms(vector* ptCompareHistogramVector, TString titles // legendEntry = Form("%s %s", getRunnumberAtBegin(curhistogramclone->GetName()).Data(), curhistogramclone->GetTitle()); legendEntry = Form("%s", curhistogramclone->GetTitle()); + if (legendStringsVector.size() > 0) + legendEntry = Form("%s", legendStringsVector.at(histogrami).Data()); leg1->AddEntry(curhistogramclone, legendEntry, "l"); leg1->Draw("SAME"); curhistogramclone->SetAxisRange(0,lastbin1*1.1,"X"); // curhistogramclone->GetYaxis()->SetRangeUser(5,heighestval1*4); - gPad->SetLogy(1); + owntitle->Clear(); owntitle->AddText(trimRunnumberAtBegin(curhistogramclone->GetName())); + if (headerStringsVector.size() > 0) { + owntitle->Clear(); + owntitle->AddText(Form("%s", headerStringsVector.at(1).Data())); + owntitle->AddText(headerStringsVector.at(2)); + } owntitle->Draw("SAME"); + + gPad->SetLogy(1); + curhistogramclone->GetYaxis()->UnZoom(); } canvas->Update(); +// gPad->Modified(); gPad->Update(); MSaveBigPNG(canvas, savepathresults + "/" + canvastitle + ".png"); TImageDump *img = new TImageDump(savepathresults + "/" + canvastitle + ".png"); @@ -700,12 +714,12 @@ 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.0); - //gPad->SetLogy(1); + gPad->SetLogy(1); gPad->SetGrid(1); gPad->SetGridy(1); owntitle->Clear(); if (headerStringsVector.size() > 0) { - owntitle->AddText(Form("%s, %s",trimRunnumberAtBegin(curhistogramclone->GetName()).Data(), headerStringsVector.at(1).Data())); + owntitle->AddText(Form("%s, %s", trimRunnumberAtBegin(curhistogramclone->GetName()).Data(), headerStringsVector.at(1).Data())); owntitle->AddText(headerStringsVector.at(2)); } owntitle->Draw("SAME"); @@ -726,7 +740,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); @@ -797,6 +811,9 @@ Bool_t plotAllRuns(vector* ptCompareHistogramClassVector) f->Write(); // gROOT->SetStyle("RadHard_AutoTitle"); + + headerStringsVector.clear(); + legendStringsVector.clear(); return 0; } return 1; @@ -817,8 +834,18 @@ Bool_t FindGoodTitle() { return FindGoodTitle(&compareHistogramClassVector); } +Bool_t FindGoodTitle(vector* ptCompareHistogramVector) { + vector compareHistogramClassVector; + for (UInt_t histogrami=0; histogrami < ptCompareHistogramVector->size(); histogrami++) + { + TH1FO* curHistogram = ptCompareHistogramVector->at(histogrami); + if (curHistogram->itsHistogramType != 0) + compareHistogramClassVector.push_back(curHistogram->itsHistogramType); + } + return FindGoodTitle(&compareHistogramClassVector, ptCompareHistogramVector); +} -Bool_t FindGoodTitle(vector* ptCompareHistogramClassVector) { +Bool_t FindGoodTitle(vector* ptCompareHistogramClassVector, vector* ptCompareHistogramVector) { // Title with everything included: // @@ -836,20 +863,40 @@ Bool_t FindGoodTitle(vector* ptCompareHistogramClassVector) { Bool_t same_Temp = kTRUE; Bool_t same_Clock = kTRUE; Bool_t same_Depletion = kTRUE; + Bool_t same_RunNumber = kTRUE; + Bool_t same_HistType = kTRUE; //Bool_t same_Annealing = kTRUE; //Bool_t same_System = kTRUE; + Bool_t mayBeSameRun = kFALSE; + if (ptCompareHistogramVector != 0) + if (ptCompareHistogramVector->size() == ptCompareHistogramClassVector->size()) + { + mayBeSameRun = kTRUE; +// same_RunNumber = kTRUE; + } + + cout << "Trying to find good title" << endlr; + headerStringsVector.clear(); legendStringsVector.clear(); - if (ptCompareHistogramClassVector->size() == 1) { + if (ptCompareHistogramClassVector->size() == 0 ) + return 1; + + if (ptCompareHistogramClassVector->at(0) == 0 ) + return 1; + + if (ptCompareHistogramClassVector->size() == 1) { HistogramType* curhistogramclassp = ptCompareHistogramClassVector->at(0); labbooksctruct curlabbook = curhistogramclassp->run->labbook; pixelinfo curpixelinfo = curhistogramclassp->run->curpixelinfo; - TString title1 = Form("%s, %s, %s, %.0fx%.0f #mum^{2} pitch, %s\nT=%.0f {}^{o}C, %.1f MRad, %.1f*10^{13} n_{eq}/cm^{2}, %.2f Mhz, U_{dep}=%.1f V", curlabbook.source.Data(), + TString title1 = Form("%s, %s, %s, %.0fx%.0f #mum^{2} pitch, %s\nT=%.0f {}^{o}C, %.1f MRad, %.1f*10^{13} n_{eq}/cm^{2}, %.2f Mhz", curlabbook.source.Data(), curlabbook.chipGen.Data(), curlabbook.matrix.Data(), curpixelinfo.pitchX, curpixelinfo.pitchY, curpixelinfo.comment.Data(), - curlabbook.tempSens, curlabbook.radDoseIon, curlabbook.radDoseNonIon, curlabbook.clock, curlabbook.depletionV); + curlabbook.tempSens, curlabbook.radDoseIon, curlabbook.radDoseNonIon, curlabbook.clock); + if (curlabbook.depletionV > 0) + title1 += Form(", U_{dep}=%.1f V", curlabbook.depletionV); TString title2 = Form("%s", curlabbook.comment.Data()); @@ -865,11 +912,17 @@ Bool_t FindGoodTitle(vector* ptCompareHistogramClassVector) { HistogramType* firsthistogramclassp = ptCompareHistogramClassVector->at(0); labbooksctruct firstlabbook = firsthistogramclassp->run->labbook; pixelinfo firstpixelinfo = firsthistogramclassp->run->curpixelinfo; + TH1FO* firsthistogramp = 0; + if (mayBeSameRun) + firsthistogramp = ptCompareHistogramVector->at(0); for (UInt_t histogrami=1; histogrami < ptCompareHistogramClassVector->size(); histogrami++) { HistogramType* curhistogramclassp = ptCompareHistogramClassVector->at(histogrami); labbooksctruct curlabbook = curhistogramclassp->run->labbook; pixelinfo curpixelinfo = curhistogramclassp->run->curpixelinfo; + TH1FO* curhistogramp = 0; + if (mayBeSameRun) + curhistogramp = ptCompareHistogramVector->at(histogrami); if (!curlabbook.source.EqualTo(firstlabbook.source)) same_Source = kFALSE; @@ -896,11 +949,21 @@ Bool_t FindGoodTitle(vector* ptCompareHistogramClassVector) { if (curlabbook.depletionV != firstlabbook.depletionV) same_Depletion= kFALSE; // if (curlabbook.system != firstlabbook.system) -// same_System = kFALSE; +// same_System = kFALSE; + if (curlabbook.runnumber != firstlabbook.runnumber) + same_RunNumber= kFALSE; + + if ( curhistogramp != 0 ) + if (!(trimRunnumberAtBegin(curhistogramp->GetName()).EqualTo(trimRunnumberAtBegin(firsthistogramp->GetName())))) + same_HistType = kFALSE; } // construct header string TString title1 = ""; + if (same_RunNumber) + title1.Append(Form(", %d", firstlabbook.runnumber)); + if (same_HistType && mayBeSameRun) + title1.Append(Form(", %s", trimRunnumberAtBegin(firsthistogramp->GetName()).Data())); if (same_Source) title1.Append(Form(", %s", firstlabbook.source.Data())); if (same_ChipGen) @@ -919,8 +982,11 @@ Bool_t FindGoodTitle(vector* ptCompareHistogramClassVector) { if (same_Depletion && firstlabbook.depletionV >= 0) title1.Append(Form(", U_{dep}=%.1f V", firstlabbook.depletionV)); +// cout << colorred << title1 << endlr; + title1 = title1.Remove(0,2); // remove trailing " ," headerStringsVector.push_back(title1); - TObjArray *humanreadablestrings = trimRunnumberAtBegin(title1).Tokenize("\n"); +// 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()); @@ -930,6 +996,14 @@ Bool_t FindGoodTitle(vector* ptCompareHistogramClassVector) { HistogramType* curhistogramclassp = ptCompareHistogramClassVector->at(histogrami); labbooksctruct curlabbook = curhistogramclassp->run->labbook; pixelinfo curpixelinfo = curhistogramclassp->run->curpixelinfo; + TH1FO* curhistogramp = 0; + if (mayBeSameRun) + curhistogramp = ptCompareHistogramVector->at(histogrami); + + 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(", "); @@ -947,8 +1021,10 @@ Bool_t FindGoodTitle(vector* ptCompareHistogramClassVector) { 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)); } + legendstr.Append(Form("U_{dep}=%.1f V", curlabbook.depletionV)); } legendStringsVector.push_back(legendstr); + +// cout << colorred << legendstr << endlr; } // Folder name suffix diff --git a/MABS_run_analyzer/HistogramType.c b/MABS_run_analyzer/HistogramType.c index 93f7ae5..566596c 100644 --- a/MABS_run_analyzer/HistogramType.c +++ b/MABS_run_analyzer/HistogramType.c @@ -42,14 +42,16 @@ HistogramType::HistogramType(TString suffix, systemparam* gotsystempar, sensorin void HistogramType::initHistograms(Int_t gotcolor, Int_t gotstyle) { color = gotcolor; style = gotstyle; + //Seed2=new TH1F(Form("%d %s",labbook->runnumber, ("Seed" + histogramdescription).Data()), Form("%s, %s",("Seed" + histogramdescription).Data(), humanreadablestr.Data()), cursystempar->nbins, 0, cursystempar->maxbin); initHistogram(Seed, "Seed" + histogramdescription, color, style); initHistogram(Sum, "Sum" + histogramdescription, color, style); + initHistogramArray(&a_Sum, 25, "Sum" + histogramdescription, color, style); initHistogram(Veto, "Veto" + histogramdescription, color, style); initHistogram(Noise, "Noise" + histogramdescription, color, style); initHistogram(NoiseEnd, "Noise at end" + histogramdescription, 16, style); initHistogramCustom(SeedPerc, "Seed Percentage" + histogramdescription, color, style, 0, 120, 121, Form("Entries [1/%%]"), "Q_coll [%]"); initHistogramCustom(pixeltimefired, "Pixel fired, used for " + histogramdescription, color, style, 0, cursensorinfo->columns*cursensorinfo->rows-1, cursensorinfo->columns*cursensorinfo->rows, "Pixel index", "Counts"); - initHistogramCustom(pixeltimefiredsorted, "Fire counter" + histogramdescription, color, style, 0, labbook->frames_foundDB/100-1, labbook->frames_foundDB/100, "times fired", "Counts"); + initHistogramCustom(pixeltimefiredDistrib, "Fire counter" + histogramdescription, color, style, 0, labbook->frames_foundDB/100-1, labbook->frames_foundDB/100, "times fired", "Counts"); initHistogramCustom(LeakageCurrentInPixel, "Leakage current per pixel" + histogramdescription, color, style, 0, cursensorinfo->columns*cursensorinfo->rows-1, cursensorinfo->columns*cursensorinfo->rows, "Pixel index", "Average CDS"); initHistogramCustom(LeakageCurrentInPixelSorted, "Leakage current per pixel sorted" + histogramdescription, color, style, 0, cursensorinfo->columns*cursensorinfo->rows-1, cursensorinfo->columns*cursensorinfo->rows, "Pixel index", "Average CDS"); // initHistogram(LeakageCurrentInPixelSorted, "Leakage current" + histogramdescription, color, style); @@ -61,8 +63,8 @@ void HistogramType::initHistograms(Int_t gotcolor, Int_t gotstyle) { // LeakageCurrentInPixelSorted->SetBins(cursystempar->nbinsnoise*10, 0, cursystempar->maxbinnoise*30); } -void HistogramType::initHistogram(TH1F* &histogrampointer, TString prefix, Int_t color, Int_t style) { - histogrampointer=new TH1F(Form("%d %s",labbook->runnumber, prefix.Data()), Form("%s, %s",prefix.Data(), humanreadablestr.Data()), cursystempar->nbins, 0, cursystempar->maxbin); +void HistogramType::initHistogram(TH1FO* &histogrampointer, TString prefix, Int_t color, Int_t style) { + histogrampointer=(TH1FO*)new TH1F(Form("%d %s",labbook->runnumber, prefix.Data()), Form("%s, %s",prefix.Data(), humanreadablestr.Data()), cursystempar->nbins, 0, cursystempar->maxbin); histogrampointer->SetLineStyle(style); histogrampointer->SetLineColor(color); histogrampointer->SetStats(kTRUE); @@ -73,10 +75,11 @@ void HistogramType::initHistogram(TH1F* &histogrampointer, TString prefix, Int_t histogrampointer->GetYaxis()->SetTitle(Form("Yield [1/%.1f ADU]",histogrampointer->GetBinWidth(1))); histogrampointer->GetXaxis()->CenterTitle(); histogrampointer->GetYaxis()->CenterTitle(); + histogrampointer->itsHistogramType = this; } -void HistogramType::initHistogramCustom(TH1F* &histogrampointer, TString prefix, Int_t color, Int_t style, Double_t xLow, Double_t xUp, Int_t nBins, TString xaxistitle, TString yaxistitle) { - histogrampointer=new TH1F(Form("%d %s",labbook->runnumber, prefix.Data()), Form("%s, %s",prefix.Data(), humanreadablestr.Data()), nBins, xLow, xUp); +void HistogramType::initHistogramCustom(TH1FO* &histogrampointer, TString prefix, Int_t color, Int_t style, Double_t xLow, Double_t xUp, Int_t nBins, TString xaxistitle, TString yaxistitle) { + histogrampointer=(TH1FO*)new TH1F(Form("%d %s",labbook->runnumber, prefix.Data()), Form("%s, %s",prefix.Data(), humanreadablestr.Data()), nBins, xLow, xUp); histogrampointer->SetLineStyle(style); histogrampointer->SetLineColor(color); histogrampointer->SetStats(kTRUE); @@ -86,6 +89,27 @@ void HistogramType::initHistogramCustom(TH1F* &histogrampointer, TString prefix, histogrampointer->GetYaxis()->SetTitle(Form("%s [1/%.1f]",yaxistitle.Data(), histogrampointer->GetBinWidth(1))); histogrampointer->GetXaxis()->CenterTitle(); histogrampointer->GetYaxis()->CenterTitle(); + histogrampointer->itsHistogramType = this; +} + + +void HistogramType::initHistogramArray(vector* histogramarraypointer, int numberofhisto, TString prefix, Int_t color, Int_t style) { + for (int i=0; i < numberofhisto; ++i) { + TH1FO* histogrampointer = (TH1FO*)new TH1F(Form("%d %s %d",labbook->runnumber, prefix.Data(), i+1), Form("%s %d, %s",prefix.Data(), i+1, humanreadablestr.Data()), cursystempar->nbins, 0, cursystempar->maxbin); + histogrampointer->SetLineStyle(style); + histogrampointer->SetLineColor(color); + histogrampointer->SetStats(kTRUE); + histogrampointer->SetStats(111111111); + histogrampointer->SetLineWidth(2); // TODO: set to 3 again + histogrampointer->GetXaxis()->SetTitle("Collected charge [ADU]"); + //histogrampointer->GetYaxis()->SetTitle(Form("Entries [1/%.1f ADU]",histogrampointer->GetBinWidth(1))); + histogrampointer->GetYaxis()->SetTitle(Form("Yield [1/%.1f ADU]",histogrampointer->GetBinWidth(1))); + histogrampointer->GetXaxis()->CenterTitle(); + histogrampointer->GetYaxis()->CenterTitle(); + histogrampointer->itsHistogramType = this; + + histogramarraypointer->push_back(histogrampointer); + } } void HistogramType::initCluster(TString suffix, Float_t xcoord_min_step, Float_t xcoord_abs_min, Float_t xcoord_abs_max, Float_t ycoord_min_step, Float_t ycoord_abs_min, Float_t ycoord_abs_max) { @@ -113,10 +137,12 @@ Bool_t HistogramType::calibrateHistograms( Float_t gotgain ) { calibrated = new HistogramType(histogramdescription+" calibrated", cursystempar, cursensorinfo, humanreadablestr, labbook, run, color, style); if (Seed != 0) calibrateHistogram(calibrated->Seed, Seed); if (Sum != 0) calibrateHistogram(calibrated->Sum, Sum); + for (unsigned int i=0; i < a_Sum.size(); ++i) + normalizeHistogram(calibrated->a_Sum[i], a_Sum[i]); if (Veto != 0) calibrateHistogram(calibrated->Veto, Veto); if (Noise != 0) calibrateHistogram(calibrated->Noise, Noise); if (pixeltimefired != 0) calibrated->pixeltimefired = pixeltimefired; - if (pixeltimefiredsorted != 0) calibrated->pixeltimefiredsorted = (TH1F*)pixeltimefiredsorted->Clone(); + if (pixeltimefiredDistrib != 0) calibrated->pixeltimefiredDistrib = (TH1FO*)pixeltimefiredDistrib->Clone(); if (SeedPerc != 0) calibrated->SeedPerc = SeedPerc; if (histAvgCluster != 0) calibrate2DHistogramCounts(calibrated->histAvgCluster, histAvgCluster); if (LeakageCurrentInPixel != 0) calibrateHistogramYAxis(calibrated->LeakageCurrentInPixel, LeakageCurrentInPixel); @@ -140,6 +166,7 @@ Bool_t HistogramType::calibrateHistograms( Float_t gotgain ) { calibrated->CCE_in_Perc_1 = CCE_in_Perc_1; calibrated->CCE_in_Perc_25 = CCE_in_Perc_25; calibrated->noisethresholdborder = noisethresholdborder; + calibrated->fixedThresholdValue = fixedThresholdValue * gain; calibrated->maskRTSpixel = maskRTSpixel; calibrated->RTSthreshold = RTSthreshold; calibrated->RTSpixel = RTSpixel; @@ -154,8 +181,8 @@ Bool_t HistogramType::calibrateHistograms( Float_t gotgain ) { return 1; } -void HistogramType::calibrateHistogram(TH1F* &histogrampointernew, TH1F* &histogrampointerold) { - histogrampointernew = (TH1F*)histogrampointerold->Clone(); +void HistogramType::calibrateHistogram(TH1FO* &histogrampointernew, TH1FO* &histogrampointerold) { + histogrampointernew = (TH1FO*)histogrampointerold->Clone(); histogrampointernew->SetName(Form("%s Calibr.", histogrampointerold->GetName())); histogrampointernew->SetTitle(Form("%s Calibr.", histogrampointerold->GetTitle())); histogrampointernew->GetXaxis()->SetTitle("Collected charge [e]"); @@ -167,11 +194,12 @@ void HistogramType::calibrateHistogram(TH1F* &histogrampointernew, TH1F* &histog } histogrampointernew->SetBins(nbins, new_bins); histogrampointernew->GetYaxis()->SetTitle(Form("Entries [1/%.1f e]",histogrampointernew->GetBinWidth(1))); + histogrampointernew->itsHistogramType = histogrampointerold->itsHistogramType; } -void HistogramType::calibrateHistogramYAxis(TH1F* &histogrampointernew, TH1F* &histogrampointerold) { - histogrampointernew = (TH1F*)histogrampointerold->Clone(); +void HistogramType::calibrateHistogramYAxis(TH1FO* &histogrampointernew, TH1FO* &histogrampointerold) { + histogrampointernew = (TH1FO*)histogrampointerold->Clone(); histogrampointernew->SetName(Form("%s Calibr.", histogrampointerold->GetName())); histogrampointernew->SetTitle(Form("%s Calibr.", histogrampointerold->GetTitle())); histogrampointernew->GetYaxis()->SetTitle("Collected charge [e]"); @@ -180,6 +208,7 @@ void HistogramType::calibrateHistogramYAxis(TH1F* &histogrampointernew, TH1F* &h histogrampointernew->SetBinContent(x,histogrampointerold->GetBinContent(x)*gain); } histogrampointernew->GetYaxis()->SetTitle(Form("Entries [1/%.1f e]",histogrampointernew->GetBinWidth(1))); + histogrampointernew->itsHistogramType = histogrampointerold->itsHistogramType; } void HistogramType::calibrate2DHistogramCounts(TH2F* &histogrampointernew, TH2F* &histogrampointerold ) { @@ -214,14 +243,17 @@ Bool_t HistogramType::normalizeHistograms( Int_t got_frames_found ) { normalized = new HistogramType(histogramdescription+" normalized", cursystempar, cursensorinfo, humanreadablestr, labbook, run, color, style); if (Seed != 0) normalizeHistogram(normalized->Seed, Seed); if (Sum != 0) normalizeHistogram(normalized->Sum, Sum); + for (unsigned int i=0; i < a_Sum.size(); ++i){ + normalizeHistogram(normalized->a_Sum[i], a_Sum[i]); + } if (Veto != 0) normalizeHistogram(normalized->Veto, Veto); if (pixeltimefired != 0) normalizeHistogram(normalized->pixeltimefired, pixeltimefired); - if (pixeltimefiredsorted != 0) normalizeHistogramXAxis(normalized->pixeltimefiredsorted, pixeltimefiredsorted); - if (SeedPerc != 0) normalized->SeedPerc = SeedPerc; - if (normalized->Noise != 0)normalized->Noise = (TH1F*)Noise->Clone(); + if (pixeltimefiredDistrib != 0) normalizeHistogramXAxis(normalized->pixeltimefiredDistrib, pixeltimefiredDistrib); + if (SeedPerc != 0) normalizeHistogram(normalized->SeedPerc, SeedPerc); + if (normalized->Noise != 0)normalized->Noise = (TH1FO*)Noise->Clone(); if (normalized->histAvgCluster != 0) normalized->histAvgCluster = (TH2F*)histAvgCluster->Clone(); - if (normalized->LeakageCurrentInPixel != 0)normalized->LeakageCurrentInPixel = (TH1F*)LeakageCurrentInPixel->Clone(); - if (normalized->LeakageCurrentInPixelSorted != 0)normalized->LeakageCurrentInPixelSorted = (TH1F*)LeakageCurrentInPixelSorted->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; @@ -235,6 +267,7 @@ Bool_t HistogramType::normalizeHistograms( Int_t got_frames_found ) { normalized->posSeedPerc = posSeedPerc; normalized->sigmaSeedPerc = sigmaSeedPerc; normalized->noisethresholdborder = noisethresholdborder; + normalized->fixedThresholdValue = fixedThresholdValue; normalized->avgNoise = avgNoise; normalized->avgNoisePlus = avgNoisePlus; normalized->avgNoiseMinus = avgNoiseMinus; @@ -257,19 +290,20 @@ Bool_t HistogramType::normalizeHistograms( Int_t got_frames_found ) { return 1; } -void HistogramType::normalizeHistogram(TH1F* &histogrampointernew, TH1F* &histogrampointerold) { - histogrampointernew = (TH1F*)histogrampointerold->Clone(); +void HistogramType::normalizeHistogram(TH1FO* &histogrampointernew, TH1FO* &histogrampointerold) { + histogrampointernew = (TH1FO*)histogrampointerold->Clone(); histogrampointernew->SetName(Form("%s norm", histogrampointerold->GetName())); histogrampointernew->SetTitle(Form("%s, norm.", histogrampointerold->GetTitle())); int nbins = histogrampointernew->GetXaxis()->GetNbins(); for(int x=0; x <= nbins; x++){ histogrampointernew->SetBinContent(x, (histogrampointerold->GetBinContent(x)*1000000.0)/frames_found); } + histogrampointernew->itsHistogramType = histogrampointerold->itsHistogramType; //histogrampointernew->GetYaxis()->SetTitle(Form("%s\b/(%d/1000000)]",histogrampointernew->GetYaxis()->GetTitle(), frames_found)); } -void HistogramType::normalizeHistogramXAxis(TH1F* &histogrampointernew, TH1F* &histogrampointerold) { - histogrampointernew = (TH1F*)histogrampointerold->Clone(); +void HistogramType::normalizeHistogramXAxis(TH1FO* &histogrampointernew, TH1FO* &histogrampointerold) { + histogrampointernew = (TH1FO*)histogrampointerold->Clone(); histogrampointernew->SetName(Form("%s norm", histogrampointerold->GetName())); histogrampointernew->SetTitle(Form("%s, norm.", histogrampointerold->GetTitle())); int nbins = histogrampointernew->GetXaxis()->GetNbins(); @@ -278,12 +312,13 @@ void HistogramType::normalizeHistogramXAxis(TH1F* &histogrampointernew, TH1F* &h new_bins[i] = histogrampointernew->GetBinLowEdge(i)*1000000/frames_found; } histogrampointernew->SetBins(nbins, new_bins); + histogrampointernew->itsHistogramType = histogrampointerold->itsHistogramType; // histogrampointernew->GetYaxis()->SetTitle(Form("%s\b/(%d/1000000)]",histogrampointernew->GetYaxis()->GetTitle(), frames_found)); } -//Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType, TFitResult* fit_result, Bool_t verbose, Float_t fitstart) { -Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType, TFitResult* fit_result_giveback, Bool_t verbose, Float_t fitstart) { +//Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType, TFitResult* fit_result, Bool_t verbose, Float_t fitstart) { +Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType, TFitResult* fit_result_giveback, Bool_t verbose, Float_t fitstart) { TFitResultPtr fit_result_ptr; Double_t* parameters = 0; Float_t posMax = 0; @@ -297,7 +332,6 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType, noiseborder = fitstart; else if (noisethresholdborder>=0 ) noiseborder = noisethresholdborder; -// cout << "noisethresholdborder: " << noisethresholdborder << endl; } // cout << colorred << " " << histogrampointer->GetName() << " : " << endlr; @@ -308,6 +342,8 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType, 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; TF1* fitFunc = new TF1("fitFunc",fitFuncType,noiseborder,posMaxValHist); if (TString(histogrampointer->GetName()).Contains("Veto")) // not used any more @@ -374,11 +410,13 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType, fitFunc->DrawCopy("same"); } else { // modified the gaussian approuch to be more powerful, escpecially the - // integral is calculated in a +- 2 sigma region. + // integral is calculated in a +- 3 sigma region. - const Double_t def_amplitude=15; + const Double_t def_amplitude=histogrampointer->GetBinContent(histogrampointer->GetMaximumBin()); const Double_t def_peakcenter=histogrampointer->GetBinCenter(histogrampointer->GetMaximumBin()); - const Double_t def_gausssig=21; +// cout << colorcyan << "def_amplitude: " << def_amplitude << endlr; +// cout << colorcyan << "def_peakcenter: " << def_peakcenter << endlr; + const Double_t def_gausssig=12; // set start values for some parameters fitFunc->SetParName(0,"amplitude of peak"); fitFunc->SetParameter(0,def_amplitude); @@ -422,6 +460,7 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType, 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) @@ -440,8 +479,8 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType, // 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 + 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), normaliezed with bin size! posMax = fitFunc->GetMaximumX(); // Methode 2 fitFunc->SetLineStyle(1); // normal for the following fits @@ -631,7 +670,7 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType, // histogrampointer->GetXaxis()->SetRangeUser(noiseborder,posMaxValHist); TF1* fitFunc = new TF1("fitFunc",GaussTail,noiseborder,posMaxValHist,6); parameters = (Double_t *)calloc(10, sizeof(Double_t)); - const Double_t def_amplitude=459.951; + const Double_t def_amplitude=histogrampointer->GetBinContent(histogrampointer->GetMaximumBin()); const Double_t def_peakcenter=histogrampointer->GetBinCenter(histogrampointer->GetMaximumBin()); // cout << colorcyan << "def_peakcenter: " << def_peakcenter << endlr; // cout << colorcyan << "histogrampointer->GetMaximumBin(): " << histogrampointer->GetMaximumBin() << endlr; @@ -690,10 +729,6 @@ Double_t* HistogramType::FitPerform(TH1F* histogrampointer, TString fitFuncType, fitFunc->SetParameter(2,def_gausssig); fitFunc->SetParameter(4,def_bgslope); fitFunc->SetParameter(5,def_bgoffs); - - // TODO: This fix disables the background - fitFunc->FixParameter(4,0); -// fitFunc->FixParameter(5,0); } else fittries = 100; @@ -709,12 +744,30 @@ Double_t* HistogramType::FitPerform(TH1F* 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; + } else + { + if (gMinuit->fCstatu == "NOT POSDEF" || gMinuit->fCstatu.Contains("FAILED")) { + failed = true; + } + } + } + if (failed) + { + failed = false; + 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 + fitFunc->SetParameter(0,def_amplitude); + fitFunc->SetParameter(1,def_peakcenter); + fitFunc->SetParameter(2,def_gausssig); + fitFunc->SetParameter(3,def_distgauss); + fitFunc->SetParameter(4,def_bgslope); + fitFunc->SetParameter(5,def_bgoffs); - // TODO: This fix disables the background - fitFunc->FixParameter(4,0); -// fitFunc->FixParameter(5,0); - fit_result_ptr = histogrampointer->Fit(fitFunc, "N,M,W,Q,S", "", noiseborder, posMaxValHist); + fit_result_ptr = histogrampointer->Fit(fitFunc, "N,M,W,Q,S", "", def_peakcenter*0.7, def_peakcenter*1.3); if (gMinuit == nullptr) { failed = true; } else @@ -815,13 +868,13 @@ Bool_t HistogramType::calculteCCE(Bool_t verbose) { return 1; } -Bool_t HistogramType::FindNoisethresholdborder(TH1F* histogrampointer, Bool_t overwrite, Bool_t verbose) { +Bool_t HistogramType::FindNoisethresholdborder(TH1FO* histogrampointer, Bool_t overwrite, Bool_t verbose) { if (overwrite || noisethresholdborder == -1) { // cout << colorcyan << "Find border" << endl; Float_t posMaxValHist = histogrampointer->GetXaxis()->GetXmax(); - TH1F* smoothedcurce = (TH1F*)histogrampointer->Clone(); + TH1FO* smoothedcurce = (TH1FO*)histogrampointer->Clone(); // smoothedcurce->Smooth(4); // not working method, rescaling is easier and more reliable Int_t rebinningfactor = 2; smoothedcurce->RebinX(rebinningfactor); @@ -894,10 +947,10 @@ Bool_t HistogramType::FindNoisethresholdborder(TH1F* histogrampointer, Bool_t ov return 1; } -Double_t HistogramType::FindBorderToPeak(TH1F* histogrampointer, Double_t startVal, Double_t endVal, Bool_t verbose) { +Double_t HistogramType::FindBorderToPeak(TH1FO* histogrampointer, Double_t startVal, Double_t endVal, Bool_t verbose) { Double_t borderval=0; Double_t peakval=0; - TH1F* smoothedcurce = (TH1F*)histogrampointer->Clone(); + TH1FO* smoothedcurce = (TH1FO*)histogrampointer->Clone(); Int_t rebinningfactor = 2; smoothedcurce->RebinX(rebinningfactor); @@ -948,7 +1001,7 @@ Double_t HistogramType::FindBorderToPeak(TH1F* histogrampointer, Double_t startV return borderval; } -Bool_t HistogramType::integrateSr90Spectra(TH1F* histogrampointer, Int_t frames_found, Float_t thresholdborder, Bool_t verbose) { +Bool_t HistogramType::integrateSr90Spectra(TH1FO* histogrampointer, Int_t frames_found, Float_t thresholdborder, Bool_t verbose) { Int_t thresholdbincurcandidate = 0; if (thresholdborder < 0) diff --git a/MABS_run_analyzer/HistogramType.h b/MABS_run_analyzer/HistogramType.h index 8a512cd..c1fbe66 100644 --- a/MABS_run_analyzer/HistogramType.h +++ b/MABS_run_analyzer/HistogramType.h @@ -16,6 +16,22 @@ #include "help.h" +class HistogramType; + +class TH1FO: public TH1F { +public: + HistogramType* itsHistogramType = 0; + + // Overload = operator + TH1FO operator=(const TH1FO &h1) { + TH1FO h2; + h2 = h1; + h2.itsHistogramType = h1.itsHistogramType; + return h2; + }; + TH1FO& operator=(const TH1F &h1); +}; + /** * @file HistogramType.h * @brief A class to store the histograms used in Run @@ -46,10 +62,11 @@ private: */ void initHistograms(Int_t gotcolor, Int_t gotstyle); /// init one specific histogram - void initHistogram(TH1F* &histogrampointer, TString prefix, Int_t color, Int_t style); + void initHistogram(TH1FO* &histogrampointer, TString prefix, Int_t color, Int_t style); /// init one custom histogram - void initHistogramCustom(TH1F* &histogrampointer, TString prefix, Int_t color, Int_t style, Double_t xLow, Double_t xUp, Int_t nBins, TString xaxistitle, TString yaxistitle); - + void initHistogramCustom(TH1FO* &histogrampointer, TString prefix, Int_t color, Int_t style, Double_t xLow, Double_t xUp, Int_t nBins, TString xaxistitle, TString yaxistitle); + /// init a histogram pointer array + void initHistogramArray(vector* histogramarraypointer, int numberofhisto, TString prefix, Int_t color, Int_t style); //***************** // GENERAL HISTOGRAM PROPERTIES //***************** @@ -70,19 +87,19 @@ private: //***************** /** * @brief rescales one specific histogram from ADU to electrons along the x axis */ - void calibrateHistogram(TH1F* &histogrampointernew, TH1F* &histogrampointerold); + void calibrateHistogram(TH1FO* &histogrampointernew, TH1FO* &histogrampointerold); /** * @brief rescales one specific histogram from ADU to electrons along the y axis */ - void calibrateHistogramYAxis(TH1F* &histogrampointernew, TH1F* &histogrampointerold); + void calibrateHistogramYAxis(TH1FO* &histogrampointernew, TH1FO* &histogrampointerold); /** * @brief normalizes one specific histogram with #frames_found along the y axis, so it can be compared graphically */ - void normalizeHistogram(TH1F* &histogrampointernew, TH1F* &histogrampointerold); + void normalizeHistogram(TH1FO* &histogrampointernew, TH1FO* &histogrampointerold); /** * @brief normalizes one specific histogram with #frames_found along the x axis, so it can be compared graphically */ - void normalizeHistogramXAxis(TH1F* &histogrampointernew, TH1F* &histogrampointerold); + void normalizeHistogramXAxis(TH1FO* &histogrampointernew, TH1FO* &histogrampointerold); //***************** // OTHER @@ -122,30 +139,32 @@ public: //***************** /// Seed spectrum, only the seed pixel is considered when this histogram is build - TH1F* Seed = 0; + TH1FO* Seed = 0; /// Sum spectrum, the charge of whole cluster is summed up in binned into this TH1F histogram - TH1F* Sum = 0; + TH1FO* Sum = 0; + /// Sum spectrum array, the charge of whole cluster is summed up in binned into this TH1F histogram + vector a_Sum; /// Veto spectrum, used to find better calibration peak, only entries where Sum over not seed pixels is below @c Run::vetothreshold are binned into this histogram - TH1F* Veto = 0; + TH1FO* Veto = 0; /// Noise histogram - TH1F* Noise = 0; + TH1FO* Noise = 0; /// Noise histogram at the end of a run - TH1F* NoiseEnd = 0; + TH1FO* NoiseEnd = 0; /// histogram with cluster data, rotated if staggered TH2F* histAvgCluster = 0; /// Percentage of charge gathered by seed pixel (Seed to Sum ratio) - TH1F* SeedPerc = 0; + TH1FO* SeedPerc = 0; /// Pixel numbers fired, the ixel number is on the x axis - TH1F* pixeltimefired = 0; + TH1FO* pixeltimefired = 0; /// Pixel numbers fired, times fired is on the x axis - TH1F* pixeltimefiredsorted = 0; + TH1FO* pixeltimefiredDistrib = 0; /// Pixel numbers fired, times fired is on the x axis - TH1F* LeakageCurrentInPixel = 0; + TH1FO* LeakageCurrentInPixel = 0; /// Pixel numbers fired, times fired is on the x axis - TH1F* LeakageCurrentInPixelSorted = 0; + TH1FO* LeakageCurrentInPixelSorted = 0; /// shows the time a specific pixel number had a "good" veto value, in the range of the veto peak - TH1F* pixelHadGoodVeto = 0; + TH1FO* pixelHadGoodVeto = 0; //***************** // FITTING RESULTS OF TH1F HISTOGRAMS @@ -191,14 +210,19 @@ public: Float_t sigmaSeedPerc = 0; /// Average leakage current in whole chip - Float_t avgLeakageCurrentInChip = 0; - + Float_t avgLeakageCurrentInChip = 0; /// threshold for integrating the Sr90 spectrum, is set dynamically in @c integrateSr90Spectra() - Float_t noisethresholdborder = -1; + Float_t noisethresholdborder = -1; + + /// threshold border used by "fixed Threshold" HistogramType instance. Only hits where charge in cluster is higher then this value are accounted + Float_t fixedThresholdValue = -1; /// percentage of pixel showing RTS signature Float_t percentageofRTSpixel = -1; + + /// depleted Volume, not used yet + Float_t depletedVolume = -1; /** @brief what is the threshold for a pixel to be considered as a RTS pixel * * It is set as a multiply of the standard deviation @@ -273,7 +297,7 @@ public: * @return peak position of the fit * */ - Double_t* FitPerform(TH1F* histogrampointer, TString fitFuncType, TFitResult* r = 0, Bool_t verbose = 0, Float_t fitstart = -1); + Double_t* FitPerform(TH1FO* histogrampointer, TString fitFuncType, TFitResult* r = 0, Bool_t verbose = 0, Float_t fitstart = -1); //Double_t* FitPerform(TH1F* histogrampointer, TString fitFuncType, TFitResultPtr* r = 0, Bool_t verbose = 0, Float_t fitstart = -1); /** @brief calculates Charge Collection Efficiency if Fe55 run */ Bool_t calculteCCE(Bool_t verbose=false); @@ -287,7 +311,7 @@ public: * @return true if succesfull * */ - Bool_t integrateSr90Spectra(TH1F* histogrampointer, Int_t frames_found=-1, Float_t thresholdborder=-1, Bool_t verbose=true); + Bool_t integrateSr90Spectra(TH1FO* histogrampointer, Int_t frames_found=-1, Float_t thresholdborder=-1, Bool_t verbose=true); /** * @brief Tries to find a border between the noise and the signal @@ -297,7 +321,7 @@ public: * @return true if succesfull * */ - Bool_t FindNoisethresholdborder(TH1F* histogrampointer, Bool_t overwrite=false, Bool_t verbose=false); + Bool_t FindNoisethresholdborder(TH1FO* histogrampointer, Bool_t overwrite=false, Bool_t verbose=false); /** @@ -306,7 +330,7 @@ public: * @return the value found * */ - Double_t FindBorderToPeak(TH1F* histogrampointer, Double_t startVal=0, Double_t endVal=0, Bool_t verbose=false); + Double_t FindBorderToPeak(TH1FO* histogrampointer, Double_t startVal=0, Double_t endVal=0, Bool_t verbose=false); /** * @brief calculates StoN ratio, if beta source was used diff --git a/MABS_run_analyzer/MAPS.c b/MABS_run_analyzer/MAPS.c index 99a22bc..7d77efa 100644 --- a/MABS_run_analyzer/MAPS.c +++ b/MABS_run_analyzer/MAPS.c @@ -1191,9 +1191,18 @@ void MAPS::hitana() { // debugStream<>(fPedestals, fPixels, fColumns, 2, 20); for(Int_t i=0; i (50.) ) + // DEBUG find out hit candidates +// A = (i)/fColumns; // A: row of seed +// B = (i)%fColumns; // B: column of seed +// if ( B == 53 && A == 6 ) { +// cout << coloryellow << "fCdsmatrix[i]: " << fCdsmatrix[i] << " fPedestals[i] " << fNoiseInfo.fPedestals[i] << " fNoise[i] " << fNoiseInfo.fNoise[i] << endlr; +// } if( (float)(1.*fCdsmatrix[i]-fNoiseInfo.fPedestals[i]) > (5.*fNoiseInfo.fNoise[i]) ) { HITNR++; @@ -1214,7 +1223,7 @@ void MAPS::hitana() { Hitlist= new Int_t[HITNR]; for(Int_t i=0; i pixelchargeincluster[12]) ) { + if ( (i!=12) && (pixelchargeincluster[i] > pixelchargeincluster[12]) ) { CHANCE=0; // cout << "kill cluster: " << Hitlist[hit] << " with ADC: " << pixelchargeincluster[12] << endl; break; @@ -1514,11 +1524,11 @@ void MAPS::hitana() { { if(fOrderCode=="FSBB") { - if (B < 2 || B > fColumns-3 || A < 1 || A > fRows-2) + if (B < 2 || B > fColumns-3 || A < 1 || A > fRows-2) // allow seed pixel to be in the second row from the edge in x-direction bordercluster = true; else bordercluster = false; - } else if(fOrderCode=="Pipper2") // collects very efficient in one seed + } else if(fOrderCode=="Pipper2") // collects very efficient in one seed, therefore smaller clusters are ok. // allow seed pixel to be in the second row from the edge in x-direction { if (B < 2 || B > fColumns-3 || A < 1 || A > fRows-2) bordercluster = true; @@ -1547,7 +1557,7 @@ void MAPS::hitana() { else if ( (column==3) && (B==fColumns-1)) { } else if ( (column==4) && (B>fColumns-3)) { } else { - fHittedPixel[Hitlist[hit]+(row-2)*fColumns+(column-2)] = bordercluster?-1:1; + fHittedPixel[Hitlist[hit]+(row-2)*fColumns+(column-2)] = bordercluster?-1:1; // mark cluster pixel which are on the edge with a "-1", good to know that a pixel belongs to the border. } } } @@ -1738,6 +1748,8 @@ void MAPS::plotFrame(Int_t FrameNumber) { { cdsmatrixCorrPed[i]=(float)(1.*fCdsmatrix[i]-fNoiseInfo.fPedestals[i]); } + cout<<"Noise:"<(fNoiseInfo.fNoise, fPixels, fColumns, 0, 5); cout<<"CDS matrix minus pedestals:"<(cdsmatrixCorrPed, fPixels, fColumns, 0, 10); // cout<Divide(2,3); gStyle->SetOptFit(1011); - gStyle->SetPalette(1); + gStyle->SetPalette(kDarkBodyRadiator); gStyle->SetCanvasColor(0); gStyle->SetFrameFillColor(10); gStyle->SetOptStat(0); - TH2F *h1 = new TH2F("CDS matrix", "CDS matrix", fColumns, 0, fColumns, fRows, 0, fRows); - TH2F *h2 = new TH2F("Frame 0 matrix", "Frame 0 matrix", fColumns, 0, fColumns, fRows, 0, fRows); - TH2F *h3 = new TH2F("Frame 1 matrix", "Frame 1 matrix", fColumns, 0, fColumns, fRows, 0, fRows); - TH2F *h6 = new TH2F("CDS matrix - Pedestial", "CDS matrix - Pedestial", fColumns, 0, fColumns, fRows, 0, fRows); + TH2F *h1 = new TH2F("CDS matrix", "CDS matrix", fColumns, 0, fColumns, fRows, -1*fRows, 0); + TH2F *h2 = new TH2F("Frame 0 matrix", "Frame 0 matrix", fColumns, 0, fColumns, fRows, -1*fRows, 0); + TH2F *h3 = new TH2F("Frame 1 matrix", "Frame 1 matrix", fColumns, 0, fColumns, fRows, -1*fRows, 0); + TH2F *h6 = new TH2F("CDS matrix - Pedestial", "CDS matrix - Pedestial", fColumns, 0, fColumns, fRows, -1*fRows, 0); TH1F *h4 = new TH1F("Frame 0 histo", "Frame 0 histo", 2*16384, -16384, 16384); TH1F *h5 = new TH1F("Frame 1 histo", "Frame 1 histo", 2*16384, -16384, 16384); @@ -1765,8 +1777,8 @@ void MAPS::plotFrame(Int_t FrameNumber) { for(Int_t i=0; iUpdate(); cout<<"\rPIXELMATRIX plotted! "<Draw(); + + + TString canvastitle = Form("%s %s", "CDS - Ped", ""); + TString canvasname = Form("%s %s", "CDS - Ped", ""); + TCanvas* canvas = new TCanvas(canvasname, canvastitle, 900, 300); + gStyle->SetOptFit(1011); + gStyle->SetPalette(70); //kBlackBody + gStyle->SetCanvasColor(0); + gStyle->SetFrameFillColor(10); + gStyle->SetOptStat(0); + h6->Draw("colz"); + canvas->Update(); cout<<"-------------------"<maskRTSpixel = true; @@ -342,10 +342,16 @@ Bool_t Run::setFixedThresholdValueElectrons(Float_t thresholdValueIne) Bool_t Run::setFixedThresholdValueADU(Float_t thresholdValue) { - cout<<" Fixed threshold value : " << colorwhite << thresholdValue << " ADU" << colorreset << " <-- only used if RAW files are analyzed, force analysis to make sure" << endl; - processed->fFixedThresholdValueInADU = thresholdValue; + return setFixedThresholdValueADU(histogramfixedthreshold, thresholdValue); +} + +Bool_t Run::setFixedThresholdValueADU(HistogramType* HistogramTypepointer, Float_t thresholdValue) +{ + cout<<" Fixed threshold value : " << colorwhite << thresholdValue << " ADU" << colorreset << endl; +// processed->fFixedThresholdValueInADU = thresholdValue; + HistogramTypepointer->fixedThresholdValue = thresholdValue; return 0; -} +} Bool_t Run::setResultsPath(TString path) { @@ -507,8 +513,8 @@ Bool_t Run::analyzeRun(Bool_t force) return 1; } -Bool_t Run::FindRTSPixelToMask(HistogramType* HistogramTypepointer) -{ +Bool_t Run::FindRTSPixelToMask(HistogramType* HistogramTypepointer, Bool_t verbose) +{ Double_t meanpixeltimesfired = 0; Double_t RTSpixelHits = 0; @@ -520,19 +526,22 @@ Bool_t Run::FindRTSPixelToMask(HistogramType* HistogramTypepointer) // first run, bin number of times pixel fired, sort by pixelindex uint pixelnumber = 0; + uint binnumber = 0; for (Int_t framei=0; frameifHitTree->GetEntries(); framei++) { // loop over all frames processed->fHitTree->GetEntry(framei); for(Int_t hiti=0; (unsigned int)hitifFrameInfo.hits;hiti++) { // loop over hits in a frame pixelnumber = processed->fFrameInfo.pixel[hiti]; - HistogramTypepointer->pixeltimefired->Fill(pixelnumber); + 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; } } // calculate the mean firing times and bin fireing times in a propability histogram, first approach - for (Int_t pixeli=0; pixeli < HistogramTypepointer->pixeltimefired->GetNbinsX(); pixeli++) { + for (Int_t pixeli=1; pixeli <= HistogramTypepointer->pixeltimefired->GetNbinsX(); pixeli++) { if (HistogramTypepointer->pixeltimefired->GetBinContent(pixeli) > 0) { - HistogramTypepointer->pixeltimefiredsorted->Fill(HistogramTypepointer->pixeltimefired->GetBinContent(pixeli)); + HistogramTypepointer->pixeltimefiredDistrib->Fill(HistogramTypepointer->pixeltimefired->GetBinContent(pixeli)); meanpixeltimesfired += HistogramTypepointer->pixeltimefired->GetBinContent(pixeli); numberofactivepixel++; // HistogramTypepointer->pixeltimefiredsorted->Fill(2000); @@ -541,16 +550,18 @@ Bool_t Run::FindRTSPixelToMask(HistogramType* HistogramTypepointer) meanpixeltimesfired /= numberofactivepixel; // Very rough estimate of a mean firing time // vrey rough estimate on standard deviation - for (Int_t pixeli=0; pixeli < HistogramTypepointer->pixeltimefired->GetNbinsX(); pixeli++) { + for (Int_t pixeli=1; pixeli <= HistogramTypepointer->pixeltimefired->GetNbinsX(); pixeli++) { if (HistogramTypepointer->pixeltimefired->GetBinContent(pixeli) > 0) stdeviation += pow(HistogramTypepointer->pixeltimefired->GetBinContent(pixeli)-meanpixeltimesfired,2); } stdeviation /= numberofactivepixel; stdeviation = sqrt(stdeviation); - cout << colorcyan << "Mean firing rate first approximation: " << meanpixeltimesfired << endlr; - cout << colorcyan << "Deviation value first approximation: " << stdeviation << endlr; - cout << colorcyan << "number of considered pixel: " << numberofactivepixel << endlr; - cout << colorcyan << endlr; + if (verbose) { + cout << colorcyan << "Mean firing rate first approximation: " << meanpixeltimesfired << endlr; + cout << colorcyan << "Deviation value first approximation: " << stdeviation << endlr; + cout << colorcyan << "number of considered pixel: " << numberofactivepixel << endlr; + cout << colorcyan << endlr; + } // better estimate on average firing time Double_t meanpixeltimesfired2, stdeviation2; @@ -560,19 +571,24 @@ Bool_t Run::FindRTSPixelToMask(HistogramType* HistogramTypepointer) u_int poscounter; - for (UInt_t i=0; i < 3; ++i) // 3 times remove most fired pixel + for (UInt_t i=0; i < 2; ++i) // 2 times remove most fired pixel { + RTSpixelHits = 0; + totalHits = 0; meanpixeltimesfired2 = 0; numberofactivepixel2=0; HistogramTypepointer->RTSpixel.clear(); - cout << "i " << i << endl; - cout << "HistogramTypepointer->pixeltimefired->GetBinContent(0) " << HistogramTypepointer->pixeltimefired->GetBinContent(0) << endl; - cout << "meanpixeltimesfired " << meanpixeltimesfired << endl; - cout << "HistogramTypepointer->RTSthreshold " << HistogramTypepointer->RTSthreshold << endl; - cout << "stdeviation " << stdeviation << endl; - for (Int_t pixeli=0; pixeli < HistogramTypepointer->pixeltimefired->GetNbinsX(); pixeli++) { - if (abs(HistogramTypepointer->pixeltimefired->GetBinContent(pixeli)-meanpixeltimesfired) > HistogramTypepointer->RTSthreshold * stdeviation) { - HistogramTypepointer->RTSpixel.push_back(pixeli); + if (verbose) { + cout << "i " << i << endl; + cout << "HistogramTypepointer->pixeltimefired->GetBinContent(1) " << HistogramTypepointer->pixeltimefired->GetBinContent(1) << endl; // first pixel + cout << "meanpixeltimesfired " << meanpixeltimesfired << endl; + cout << "HistogramTypepointer->RTSthreshold " << HistogramTypepointer->RTSthreshold << endl; + cout << "stdeviation " << stdeviation << endl; + } + for (Int_t pixeli=1; pixeli <= HistogramTypepointer->pixeltimefired->GetNbinsX(); pixeli++) { + totalHits += HistogramTypepointer->pixeltimefired->GetBinContent(pixeli); + if ( abs(HistogramTypepointer->pixeltimefired->GetBinContent(pixeli)-meanpixeltimesfired) > HistogramTypepointer->RTSthreshold * stdeviation ) { + HistogramTypepointer->RTSpixel.push_back(pixeli-1); // Pixel no 0 is saved in bin 1 (!) RTSpixelHits+=HistogramTypepointer->pixeltimefired->GetBinContent(pixeli); } else @@ -586,47 +602,17 @@ Bool_t Run::FindRTSPixelToMask(HistogramType* HistogramTypepointer) meanpixeltimesfired2 /= numberofactivepixel2; meanpixeltimesfired = meanpixeltimesfired2; } - -// Double_t const probabilities[] = {1-erf(HistogramTypepointer->RTSthreshold/sqrt(2)), 0.5, erf(HistogramTypepointer->RTSthreshold/sqrt(2))}; // sigma/2 from gaus to the left and to the right //{0.17, 0.5, 1-0.17}; -// HistogramTypepointer->pixeltimefiredsorted->GetQuantiles( 3, leakagequantiles, probabilities); - -// HistogramTypepointer->medianLeakageCurrent = leakagequantiles[1]; -// HistogramTypepointer->medianLeakageCurrentPlus = leakagequantiles[2] - leakagequantiles[1]; -// HistogramTypepointer->medianLeakageCurrentMinus = leakagequantiles[1] - leakagequantiles[0]; -// cout << "Median is at: " << HistogramTypepointer->medianLeakageCurrent << endl; -// cout << "Upper limit is at: " << leakagequantiles[2] << endl; -// cout << "Lower limit is at: " << leakagequantiles[0] << endl; -// cout << "Plus: " << HistogramTypepointer->medianLeakageCurrentPlus << endl; -// cout << "minus: " << HistogramTypepointer->medianLeakageCurrentMinus << endl; - -// cout << "last value above 1: " << HistogramTypepointer->pixeltimefiredsorted->GetBinCenter(HistogramTypepointer->pixeltimefiredsorted->FindLastBinAbove(4)) << endl; -// cout << "last value above 1 in overflow TH1F: " << pixeltimefiredWithOverflow->GetBinCenter(pixeltimefiredWithOverflow->FindLastBinAbove(4)) << endl; - // cout << colorcyan << "Deviation value after one more loop: " << stdeviation << endlr; - -// HistogramTypepointer->RTSpixel.clear(); -// for (Int_t pixeli=0; pixeli < HistogramTypepointer->pixeltimefired->GetNbinsX(); pixeli++) { -// if (HistogramTypepointer->pixeltimefired->GetBinContent(pixeli) > leakagequantiles[2]) { -// // cout << coloryellow << numberToString(5) << " mal drüber, Böser pixel #" << pixeli << " Wert: " << HistogramTypepointer->pixeltimefired->GetBinContent(pixeli) << endlr; -// HistogramTypepointer->RTSpixel.push_back(pixeli); -// } -// } -// cout << colorred << " number of RTS pixel in quantile run: " << HistogramTypepointer->RTSpixel.size() << endlr; - - -// meanpixeltimesfired /= HistogramTypepointer->pixeltimefired->GetNbinsX(); // Very rough estimate of a mean firing time -// meanpixeltimesfired = leakagequantiles[1]; -// cout << colorcyan << "Mean value own: " << meanpixeltimesfired << endlr; // Estimate new standard deviation RTSpixel =false; numberofconsideredpixel=0; poscounter = 0; stdeviation2 = 0; - for (UInt_t pixeli=0; pixeli < (UInt_t)HistogramTypepointer->pixeltimefired->GetNbinsX(); pixeli++) { + for (UInt_t pixeli=1; pixeli <= (UInt_t)HistogramTypepointer->pixeltimefired->GetNbinsX(); pixeli++) { if (HistogramTypepointer->pixeltimefired->GetBinContent(pixeli) > 0) { RTSpixel = false; for (u_int RTSpixeli=0; RTSpixeli < (u_int) HistogramTypepointer->RTSpixel.size(); RTSpixeli++) { - if (pixeli == HistogramTypepointer->RTSpixel[RTSpixeli]) { + if (pixeli-1 == HistogramTypepointer->RTSpixel[RTSpixeli]) { poscounter = RTSpixeli; RTSpixel = true; } @@ -642,10 +628,12 @@ Bool_t Run::FindRTSPixelToMask(HistogramType* HistogramTypepointer) stdeviation2 = sqrt(stdeviation); stdeviation = stdeviation2; } - cout << colorred << " number of RTS pixel in " << i << " run: " << HistogramTypepointer->RTSpixel.size() << endlr; - cout << colorred << " number of considered pixel: " << numberofconsideredpixel << endlr; - cout << colorcyan << " mean value: " << meanpixeltimesfired << endl; - cout << colorcyan << " Deviation value: " << stdeviation << endlr; + if (verbose) { + cout << colorred << " number of RTS pixel in " << i << " run: " << HistogramTypepointer->RTSpixel.size() << endlr; + cout << colorred << " number of considered pixel: " << numberofconsideredpixel << endlr; + cout << colorcyan << " mean value: " << meanpixeltimesfired << endl; + cout << colorcyan << " Deviation value: " << stdeviation << endlr; + } @@ -680,6 +668,10 @@ Bool_t Run::FindRTSPixelToMask(HistogramType* HistogramTypepointer) // 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]) HistogramTypepointer->percentageofRTSpixel = (HistogramTypepointer->RTSpixel.size()*1.0)/(numberofactivepixel*1.0)*100.0; @@ -761,11 +753,11 @@ Bool_t Run::analyzeFrame(Int_t frame) if (!error) { processed = new MAPS(this); - processed->initNewRootFile(); + processed->initOldRootFile(); int entries = processed->GetNumberFrames(); if (frame < entries) { - processed->InitialDynNoise(); + processed->InitialDynNoise(frame-100); processed->getFrame(frame); if (commonModeFilter) processed->filterCommonMode(); @@ -842,7 +834,9 @@ Bool_t Run::generateReadableRunCode() if (labbook.radDoseNonIon > 0) runcode+= Form("Rn%.0f",labbook.radDoseNonIon*1000); else - runcode+= Form("Rn0%.0f",labbook.radDoseNonIon); + runcode+= Form("Rn0%.0f",labbook.radDoseNonIon); + if (labbook.depletionV > 0) + runcode+= Form("Ud%.1f",labbook.depletionV); runcode+= labbook.matrix; runcode+= runcodesuffix; if (labbook.resistivity > 0) @@ -982,18 +976,27 @@ void Run::getVetoPeakPositionFromFe55Run() } } -void Run::constructUpdateString(string *sqlupdatequery, const string databasevaluename, const Double_t value, const int precision=3) +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; - if (!std::isinf(value)) - { - if (abs(value)>0.0) - { - if ((*sqlupdatequery).length() > 0) - *sqlupdatequery+= ", "; - *sqlupdatequery += "`" + databasevaluename + "`="+ to_str_w_prec(value, precision); - } - } + 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); } void Run::updateDatabase() { @@ -1005,36 +1008,32 @@ void Run::updateDatabase() { } else { histogramclassToUseForDB = histogramwoRTSthreshold; } - + string sqlupdatequery = ""; - constructUpdateString(&sqlupdatequery, "Gain", histogramclassToUseForDB->normalized->gain); - constructUpdateString(&sqlupdatequery, "SumPeak", histogramclassToUseForDB->normalized->posSum, 4); - constructUpdateString(&sqlupdatequery, "SeedPeak", histogramclassToUseForDB->normalized->posSeed, 4); - constructUpdateString(&sqlupdatequery, "AvgF0", labbook.averageF0, 6); - constructUpdateString(&sqlupdatequery, "SigmaF0", labbook.sigmaF0, 6); - constructUpdateString(&sqlupdatequery, "SeedIntegralMinusBG2Sigma", histogramclassToUseForDB->normalized->integralSeed, 10); - constructUpdateString(&sqlupdatequery, "SeedIntegralMinusBG2SigmaErr", histogramclassToUseForDB->normalized->integralSeedErr, 5); - constructUpdateString(&sqlupdatequery, "VetoPeak", histogramclassToUseForDB->normalized->posVeto, 4); - constructUpdateString(&sqlupdatequery, "VetoIntegral", histogramclassToUseForDB->normalized->integralVeto, 10); - constructUpdateString(&sqlupdatequery, "SumIntegral", histogramclassToUseForDB->normalized->integralSum, 10); - constructUpdateString(&sqlupdatequery, "SumIntegralErr", histogramclassToUseForDB->normalized->integralSumErr, 5); + constructUpdateString(&sqlupdatequery, "Gain", histogramclassToUseForDB->normalized->gain, 3, 0, 100); + constructUpdateString(&sqlupdatequery, "SumPeak", histogramclassToUseForDB->normalized->posSum, 4, 0, 1000); + constructUpdateString(&sqlupdatequery, "SeedPeak", histogramclassToUseForDB->normalized->posSeed, 4, 0, 1000); + constructUpdateString(&sqlupdatequery, "AvgF0", labbook.averageF0, 6, 0, 10000); + constructUpdateString(&sqlupdatequery, "SigmaF0", labbook.sigmaF0, 6, 0, 10000); + constructUpdateString(&sqlupdatequery, "SeedIntegralMinusBG2Sigma", histogramclassToUseForDB->normalized->integralSeed, 10, 0, 1e7); + constructUpdateString(&sqlupdatequery, "SeedIntegralMinusBG2SigmaErr", histogramclassToUseForDB->normalized->integralSeedErr, 5, 0, 1e7); + constructUpdateString(&sqlupdatequery, "VetoPeak", histogramclassToUseForDB->normalized->posVeto, 4, 0, 1000); + constructUpdateString(&sqlupdatequery, "VetoIntegral", histogramclassToUseForDB->normalized->integralVeto, 10, 0, 1e7); + constructUpdateString(&sqlupdatequery, "SumIntegral", histogramclassToUseForDB->normalized->integralSum, 10, 0, 1e7); + constructUpdateString(&sqlupdatequery, "SumIntegralErr", histogramclassToUseForDB->normalized->integralSumErr, 5, 0, 1e7); if (histogramclassToUseForDB->normalized->calibrated != 0) -// if (histogramthreshold->normalized->calibrated->avgNoise < 100) - constructUpdateString(&sqlupdatequery, "Avg.Noise", histogramclassToUseForDB->normalized->calibrated->avgNoise); + constructUpdateString(&sqlupdatequery, "Avg.Noise", histogramclassToUseForDB->normalized->calibrated->avgNoise, 3, 0, 1000); if (histogramclassToUseForDB->normalized->calibrated != 0) -// if (histogramthreshold->normalized->calibrated->avgNoisePlus < 100) - constructUpdateString(&sqlupdatequery, "Avg.Noise+", histogramclassToUseForDB->normalized->calibrated->avgNoisePlus, 2); + constructUpdateString(&sqlupdatequery, "Avg.Noise+", histogramclassToUseForDB->normalized->calibrated->avgNoisePlus, 3, 0, 1000); if (histogramclassToUseForDB->normalized->calibrated != 0) -// if (histogramthreshold->normalized->calibrated->avgNoiseMinus < 100) - constructUpdateString(&sqlupdatequery, "Avg.Noise-", histogramclassToUseForDB->normalized->calibrated->avgNoiseMinus, 2); - constructUpdateString(&sqlupdatequery, "CCE_1", histogramclassToUseForDB->normalized->CCE_in_Perc_1); - constructUpdateString(&sqlupdatequery, "CCE_25", histogramclassToUseForDB->normalized->CCE_in_Perc_25); - constructUpdateString(&sqlupdatequery, "StoN", histogramclassToUseForDB->normalized->StoN, 3); -// if (histogramclassToUseForDB->normalized->avgNoise < 100) - constructUpdateString(&sqlupdatequery, "Avg.NoiseADC", histogramclassToUseForDB->normalized->avgNoise); - constructUpdateString(&sqlupdatequery, "Avg.NoiseADC+", histogramclassToUseForDB->normalized->avgNoisePlus); - constructUpdateString(&sqlupdatequery, "Avg.NoiseADC-", histogramclassToUseForDB->normalized->avgNoiseMinus); - constructUpdateString(&sqlupdatequery, "Frames_found", frames_found,100000000); + constructUpdateString(&sqlupdatequery, "Avg.Noise-", histogramclassToUseForDB->normalized->calibrated->avgNoiseMinus, 3, 0, 1000); + constructUpdateString(&sqlupdatequery, "CCE_1", histogramclassToUseForDB->normalized->CCE_in_Perc_1, 2, 0 , 100); + constructUpdateString(&sqlupdatequery, "CCE_25", histogramclassToUseForDB->normalized->CCE_in_Perc_25, 2, 0 , 100); + constructUpdateString(&sqlupdatequery, "StoN", histogramclassToUseForDB->normalized->StoN, 3, 0 , 100); + constructUpdateString(&sqlupdatequery, "Avg.NoiseADC", histogramclassToUseForDB->normalized->avgNoise, 2, 0 , 1000); + constructUpdateString(&sqlupdatequery, "Avg.NoiseADC+", histogramclassToUseForDB->normalized->avgNoisePlus, 2, 0 , 1000); + constructUpdateString(&sqlupdatequery, "Avg.NoiseADC-", histogramclassToUseForDB->normalized->avgNoiseMinus, 2, 0 , 1000); + constructUpdateString(&sqlupdatequery, "Frames_found", frames_found, 100000000); if (histogramwoRTS != 0) { constructUpdateString(&sqlupdatequery, "RTSpixel", histogramwoRTS->normalized->RTSpixel.size()); constructUpdateString(&sqlupdatequery, "RTSpixel_percentage", histogramwoRTS->normalized->percentageofRTSpixel); @@ -1045,12 +1044,12 @@ void Run::updateDatabase() { constructUpdateString(&sqlupdatequery, "LeakageCurfA", histogramwoRTS->calibrated->medianLeakageCurrent); constructUpdateString(&sqlupdatequery, "LeakageCurfA+", histogramwoRTS->calibrated->medianLeakageCurrentPlus); constructUpdateString(&sqlupdatequery, "LeakageCurfA-", histogramwoRTS->calibrated->medianLeakageCurrentMinus); - constructUpdateString(&sqlupdatequery, "CalibrationPeak", Fe55run.posVeto); + constructUpdateString(&sqlupdatequery, "CalibrationPeak", Fe55run.posVeto, 4, 0, 1000); } } if (histogramclassToUseForDB->normalized->calibrated != 0) if (labbook.source.Contains("Sr") && histogramclassToUseForDB->normalized->calibrated->sr90IntegralVal > 0) - constructUpdateString(&sqlupdatequery, "Sr90IntegralVal", histogramclassToUseForDB->normalized->calibrated->sr90IntegralVal,1000000000); + constructUpdateString(&sqlupdatequery, "Sr90IntegralVal", histogramclassToUseForDB->normalized->calibrated->sr90IntegralVal, 1000000000, 0,1e7); if (sqlupdatequery.length()>0) { @@ -1254,19 +1253,29 @@ Bool_t Run::binF0() Bool_t Run::binSeedSumVeto() { /// collected charge in cluster + Int_t completeclustersize = processed->clustersize*processed->clustersize; + Float_t a_pixelSum[completeclustersize]; + Float_t a_notSeedSum[completeclustersize]; Float_t pixelSum = 0; Float_t notSeedSum = 0; + for (Int_t clusteri=0; clusterifHitTree->GetEntries(); framei++) // loop over all frames { processed->fHitTree->GetEntry(framei); // account only frames with less then 10 hits -// cout << colorcyan << processed->fFrameInfo.hits << endlr; +// cout << colorcyan << processed->fFrameInfo.hits << endlr; // if (processed->fFrameInfo.hits<(unsigned int)10) { for(Int_t hiti=0; (unsigned int)hitifFrameInfo.hits;hiti++) { +// if (!hiti && framei==1994914) +// cout << coloryellow << " " << framei << endlr; uint pixel_column_x = processed->fFrameInfo.pixel[hiti]%cursensorinfo.columns; // column of seed uint pixel_row_y = processed->fFrameInfo.pixel[hiti]/cursensorinfo.columns; // row of seed if (false) { @@ -1282,54 +1291,85 @@ Bool_t Run::binSeedSumVeto() histogram->Seed->Fill(processed->fFrameInfo.p[12][hiti]); // cout << colorcyan << "filled seed" << endlr; - // sum histogram pixelSum = 0; - notSeedSum = 0; - //if(labbook.chipGen=="FSBB" || labbook.chipGen=="Mi19") - if(labbook.chipGen=="FSBB") - { - Float_t clusterArray[processed->clustersize*processed->clustersize];// temp variable clusterArray necessary, because Sort only accepts 1-dim arrays - Int_t index[processed->clustersize*processed->clustersize]; - for (Int_t clusteri=0; clustericlustersize*processed->clustersize; clusteri++) - { - clusterArray[clusteri] = processed->fFrameInfo.p[clusteri][hiti]; - if (clusteri != 12) - notSeedSum += processed->fFrameInfo.p[clusteri][hiti]; - } - TMath::Sort(processed->clustersize*processed->clustersize,clusterArray,index,1); - for (Int_t clusteri=0; clusteri<4; clusteri++) - { - pixelSum += clusterArray[index[clusteri]]; - - } + notSeedSum = 0; + // sum histogram + for (Int_t clusteri=0; clusterifFrameInfo.p[clusteri][hiti] = 0; +// for (Int_t clusteri=completeclustersize-sqrt(completeclustersize); clusterifFrameInfo.p[clusteri][hiti] = 0; +// for (Int_t clusteri=0; clusterifFrameInfo.p[clusteri][hiti] = 0; +// for (Int_t clusteri=4; clusterifFrameInfo.p[clusteri][hiti] = 0; +// } + + // DEBUG +// for (Int_t clusteri=0; clusterifFrameInfo.p[clusteri*5+clusterj][hiti]; +// cout << endl; +// } +// cout << "......." << endl; + + Float_t clusterArray[completeclustersize];// temp variable clusterArray necessary, because Sort only accepts 1-dim arrays + for (Int_t clusteri=0; clusterifFrameInfo.p[clusteri][hiti]; + } + Int_t index[completeclustersize]; + TMath::Sort(completeclustersize,clusterArray,index,1); + for (Int_t clusteri=0; clustericlustersize*processed->clustersize; clusteri++) - { - pixelSum += processed->fFrameInfo.p[clusteri][hiti]; - if (clusteri != 12) - notSeedSum += processed->fFrameInfo.p[clusteri][hiti]; + pixelSum += clusterArray[index[clusteri]]; + for (Int_t clusterj=clusteri; clusterjfFrameInfo.p[clusteri][hiti]; + for (Int_t clusterj=clusteri; clusterjSum->Fill(pixelSum); + for (Int_t clusterj=0; clusterja_Sum[clusterj]->Fill(a_pixelSum[clusterj]); + } // seed percentage spectrum - histogram->SeedPerc->Fill(processed->fFrameInfo.p[12][hiti]/pixelSum*100); + histogram->SeedPerc->Fill(processed->fFrameInfo.p[12][hiti]/pixelSum*100); +// if (framei==1994914) { +// cout << "pixel_column_x: " << pixel_column_x << ", pixel_row_y: " << pixel_row_y << " pixelcharge: " << processed->fFrameInfo.p[12][hiti] << endl; +// cout << coloryellow << processed->fFrameInfo.p[12][hiti]/pixelSum*100 << endlr; +// } // 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 + + // bin the fixed threshold for charge histogram class + if (histogramfixedthreshold != 0) { + if (pixelSum > histogramfixedthreshold->fixedThresholdValue) // charge is more then histogramfixedthreshold->fixedThresholdValue in whole cluster + { + histogramfixedthreshold->numberofhits++; + fillAHistogramsinclass(histogramfixedthreshold, hiti, completeclustersize, pixelSum, notSeedSum, &a_pixelSum[0], &a_notSeedSum[0]); + } + } + if (processed->fFrameInfo.pixelthreshold[hiti]>0) { - histogramthreshold->numberofhits++; - - histogramthreshold->Seed->Fill(processed->fFrameInfo.p[12][hiti]); - histogramthreshold->Sum->Fill(pixelSum); - histogramthreshold->SeedPerc->Fill(processed->fFrameInfo.p[12][hiti]/pixelSum*100); - if (TMath::Abs(notSeedSum) < cursystemparam.vetothreshold && (labbook.source.Contains("Fe")||labbook.source.Contains("Cd"))) - histogramthreshold->Veto->Fill(processed->fFrameInfo.p[12][hiti]); // histogram with the single pixel + histogramthreshold->numberofhits++; + fillAHistogramsinclass(histogramthreshold, hiti, completeclustersize, pixelSum, notSeedSum, &a_pixelSum[0], &a_notSeedSum[0]); // bin the RTS cleaned histogram class if (histogramwoRTSthreshold != 0) { @@ -1341,12 +1381,8 @@ Bool_t Run::binSeedSumVeto() // cout << "not binned! RTS pixel #" << processed->fFrameInfo.pixel[hiti] << endl; } } - if (!RTSpixel) { - histogramwoRTSthreshold->Seed->Fill(processed->fFrameInfo.p[12][hiti]); - histogramwoRTSthreshold->Sum->Fill(pixelSum); - histogramwoRTSthreshold->SeedPerc->Fill(processed->fFrameInfo.p[12][hiti]/pixelSum*100); - if (TMath::Abs(notSeedSum) < cursystemparam.vetothreshold && (labbook.source.Contains("Fe")||labbook.source.Contains("Cd"))) - histogramwoRTSthreshold->Veto->Fill(processed->fFrameInfo.p[12][hiti]); // histogram with the single pixel + if (!RTSpixel) { + fillAHistogramsinclass(histogramwoRTSthreshold, hiti, completeclustersize, pixelSum, notSeedSum, &a_pixelSum[0], &a_notSeedSum[0]); } } // bin the more agressive RTS histogram class @@ -1360,46 +1396,30 @@ Bool_t Run::binSeedSumVeto() } } if (!RTSpixel) { - histogramwoRTSAggresivethreshold->Seed->Fill(processed->fFrameInfo.p[12][hiti]); - histogramwoRTSAggresivethreshold->Sum->Fill(pixelSum); - histogramwoRTSAggresivethreshold->SeedPerc->Fill(processed->fFrameInfo.p[12][hiti]/pixelSum*100); - if (TMath::Abs(notSeedSum) < cursystemparam.vetothreshold && (labbook.source.Contains("Fe")||labbook.source.Contains("Cd"))) - histogramwoRTSAggresivethreshold->Veto->Fill(processed->fFrameInfo.p[12][hiti]); // histogram with the single pixel + fillAHistogramsinclass(histogramwoRTSAggresivethreshold, hiti, completeclustersize, pixelSum, notSeedSum, &a_pixelSum[0], &a_notSeedSum[0]); } } } - - // bin the fixed threshold for charge histogram class - if (histogramfixedthreshold != 0) { - if (processed->fFrameInfo.pixelfixedthreshold[hiti]>0) - { - histogramfixedthreshold->numberofhits++; - - histogramfixedthreshold->Seed->Fill(processed->fFrameInfo.p[12][hiti]); - histogramfixedthreshold->Sum->Fill(pixelSum); - histogramfixedthreshold->SeedPerc->Fill(processed->fFrameInfo.p[12][hiti]/pixelSum*100); - if (TMath::Abs(notSeedSum) < cursystemparam.vetothreshold && (labbook.source.Contains("Fe")||labbook.source.Contains("Cd"))) - histogramfixedthreshold->Veto->Fill(processed->fFrameInfo.p[12][hiti]); // histogram with the single pixel - } - } // bin the RTS cleaned histogram class if (histogramwoRTS != 0) { +// cout << "is " << processed->fFrameInfo.pixel[hiti] << " an RTS pixel? size of rts pixel matrix: " << histogramwoRTS->RTSpixel.size(); RTSpixel = false; for (u_int RTSpixeli=0; RTSpixeli < histogramwoRTS->RTSpixel.size(); RTSpixeli++) { if (processed->fFrameInfo.pixel[hiti] == histogramwoRTS->RTSpixel[RTSpixeli]) { + // if (RTSpixeli == 0) RTSpixel = true; -// cout << "not binned! RTS pixel #" << processed->fFrameInfo.pixel[hiti] << " with charge: " << processed->fFrameInfo.p[12][hiti] << endl; + break; } + } +// if (RTSpixel) +// cout << " yes"; +// cout << endlr; if (!RTSpixel) { - histogramwoRTS->Seed->Fill(processed->fFrameInfo.p[12][hiti]); - histogramwoRTS->Sum->Fill(pixelSum); - histogramwoRTS->SeedPerc->Fill(processed->fFrameInfo.p[12][hiti]/pixelSum*100); - if (TMath::Abs(notSeedSum) < cursystemparam.vetothreshold && (labbook.source.Contains("Fe")||labbook.source.Contains("Cd"))) - histogramwoRTS->Veto->Fill(processed->fFrameInfo.p[12][hiti]); // histogram with the single pixel - } + fillAHistogramsinclass(histogramwoRTS, hiti, completeclustersize, pixelSum, notSeedSum, &a_pixelSum[0], &a_notSeedSum[0]); + } } // bin the more agressive RTS histogram class if (histogramwoRTSAggresive != 0) { @@ -1413,19 +1433,13 @@ Bool_t Run::binSeedSumVeto() } } if (!RTSpixel) { -// cout << "binned! pixel #" << processed->fFrameInfo.pixel[hiti] << " with charge: " << processed->fFrameInfo.p[12][hiti] << endl; - histogramwoRTSAggresive->Seed->Fill(processed->fFrameInfo.p[12][hiti]); - histogramwoRTSAggresive->Sum->Fill(pixelSum); - histogramwoRTSAggresive->SeedPerc->Fill(processed->fFrameInfo.p[12][hiti]/pixelSum*100); - if (TMath::Abs(notSeedSum) < cursystemparam.vetothreshold && (labbook.source.Contains("Fe")||labbook.source.Contains("Cd"))) - histogramwoRTSAggresive->Veto->Fill(processed->fFrameInfo.p[12][hiti]); // histogram with the single pixel + fillAHistogramsinclass(histogramwoRTSAggresive, hiti, completeclustersize, pixelSum, notSeedSum, &a_pixelSum[0], &a_notSeedSum[0]); } } } // end if in range for submatrix analysis } // end loop over hits in frame } } // end loop over all frames - cout << colorred << " histogramwoRTSAggresive->Seed size: " << histogramwoRTSAggresive->Seed->GetEntries() << endlr; // gROOT->SetBatch(kTRUE); for (vector::iterator curHistogramClass = HistogramClassVector.begin(); curHistogramClass != HistogramClassVector.end(); curHistogramClass++) { Double_t* parameters = (Double_t *)calloc(11, sizeof(Double_t)); // allocate 11 parameters for safety, maximum 10 used at the moment @@ -1449,13 +1463,13 @@ Bool_t Run::binSeedSumVeto() (*curHistogramClass)->integralVeto = parameters[6]; } if (labbook.chipGen.EqualTo("Pipper2")) - parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->Seed, "GaussTail", 0, false, 500); + parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->Seed, "GaussTail", 0, false, (*curHistogramClass)->fixedThresholdValue); else parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->Seed, "landau"); (*curHistogramClass)->integralSeed = parameters[6]; (*curHistogramClass)->posSeed = parameters[1]; (*curHistogramClass)->integralSeedErr = parameters[9]; - parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->Sum, "GaussTail", 0, false, 500); //TODO change back to gauss + parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->Sum, "GaussTail", 0, false, (*curHistogramClass)->fixedThresholdValue); //TODO change back to gauss (*curHistogramClass)->posSum = parameters[1]; (*curHistogramClass)->integralSum = parameters[6]; (*curHistogramClass)->integralSumErr = parameters[9]; @@ -1464,12 +1478,8 @@ Bool_t Run::binSeedSumVeto() (*curHistogramClass)->posSeedPerc = parameters[1]; (*curHistogramClass)->sigmaSeedPerc = parameters[2]; - for (Int_t bini=1; bini <= cursystemparam.nbins; bini++) { // TODO: rescaled histogram to number of frames found, remove the inner for loop later -// (*curHistogramClass)->Seed->SetBinContent(bini,(*curHistogramClass)->Seed->GetBinContent(bini)/(frames_found*1.0)); -// (*curHistogramClass)->Sum->SetBinContent(bini,(*curHistogramClass)->Sum->GetBinContent(bini)/(frames_found*1.0)); - } for (Int_t bini=1; bini <= 100; bini++) { - (*curHistogramClass)->SeedPerc->SetBinContent(bini,(*curHistogramClass)->SeedPerc->GetBinContent(bini)/0.01/((*curHistogramClass)->SeedPerc->GetEntries())); + (*curHistogramClass)->SeedPerc->SetBinContent(bini,(*curHistogramClass)->SeedPerc->GetBinContent(bini));//((*curHistogramClass)->SeedPerc->GetEntries())); } } if (histogramGoodVeto != 0 && histogram != 0) { @@ -1600,6 +1610,18 @@ Bool_t Run::binSeedSumVeto() return 0; } +void Run::fillAHistogramsinclass(HistogramType* histogramtypepointer, Int_t hiti, Int_t completeclustersize, Float_t pixelSum, Float_t notSeedSum, Float_t* pt_a_pixelSum, Float_t* pt_a_notSeedSum) { + histogramtypepointer->Seed->Fill(processed->fFrameInfo.p[12][hiti]); + histogramtypepointer->Sum->Fill(pixelSum); + for (Int_t clusterj=0; clusterja_Sum[clusterj]->Fill(pt_a_pixelSum[clusterj]); + } + histogramtypepointer->SeedPerc->Fill(processed->fFrameInfo.p[12][hiti]/pixelSum*100); + 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 +} + + Bool_t Run::binCluster() { Float_t rotateangle = 0; @@ -2104,13 +2126,13 @@ Bool_t Run::plotClusterDistribution(HistogramType* HistogramTypepointer) // canvas->Print(savepathresults + "/" + canvastitle + ".tex"); // create lorentz fit of a slice - TH1F* xslicetroughcluster = new TH1F(Form("%s_xslice", HistogramTypepointer->histAvgCluster->GetTitle()),"X slice",HistogramTypepointer->histAvgCluster->GetXaxis()->GetNbins(),HistogramTypepointer->histAvgCluster->GetXaxis()->GetBinLowEdge(1),HistogramTypepointer->histAvgCluster->GetXaxis()->GetBinUpEdge(HistogramTypepointer->histAvgCluster->GetXaxis()->GetNbins())); + TH1FO* xslicetroughcluster = (TH1FO*) new TH1F(Form("%s_xslice", HistogramTypepointer->histAvgCluster->GetTitle()),"X slice",HistogramTypepointer->histAvgCluster->GetXaxis()->GetNbins(),HistogramTypepointer->histAvgCluster->GetXaxis()->GetBinLowEdge(1),HistogramTypepointer->histAvgCluster->GetXaxis()->GetBinUpEdge(HistogramTypepointer->histAvgCluster->GetXaxis()->GetNbins())); Double_t* xslicetroughclusterpar = (Double_t *)calloc(11, sizeof(Double_t)); // allocate 11 parameters for safety, maximum 10 used at the moment Int_t middlebin = (int)(HistogramTypepointer->histAvgCluster->GetXaxis()->GetNbins()/2.0 + 0.5); for (Int_t bini=1; bini <= HistogramTypepointer->histAvgCluster->GetXaxis()->GetNbins(); bini++) xslicetroughcluster->SetBinContent(bini,HistogramTypepointer->histAvgCluster->GetBinContent(bini,middlebin)); - TH1F* yslicetroughcluster = new TH1F(Form("%s_yslice", HistogramTypepointer->histAvgCluster->GetTitle()),"Y slice",HistogramTypepointer->histAvgCluster->GetYaxis()->GetNbins(),HistogramTypepointer->histAvgCluster->GetYaxis()->GetBinLowEdge(1),HistogramTypepointer->histAvgCluster->GetYaxis()->GetBinUpEdge(HistogramTypepointer->histAvgCluster->GetYaxis()->GetNbins())); + TH1FO* yslicetroughcluster = (TH1FO*) new TH1F(Form("%s_yslice", HistogramTypepointer->histAvgCluster->GetTitle()),"Y slice",HistogramTypepointer->histAvgCluster->GetYaxis()->GetNbins(),HistogramTypepointer->histAvgCluster->GetYaxis()->GetBinLowEdge(1),HistogramTypepointer->histAvgCluster->GetYaxis()->GetBinUpEdge(HistogramTypepointer->histAvgCluster->GetYaxis()->GetNbins())); Double_t* yslicetroughclusterpar = (Double_t *)calloc(11, sizeof(Double_t)); // allocate 11 parameters for safety, maximum 10 used at the moment middlebin = (int)(HistogramTypepointer->histAvgCluster->GetYaxis()->GetNbins()/2.0 + 0.5); for (Int_t bini=1; bini <= HistogramTypepointer->histAvgCluster->GetYaxis()->GetNbins(); bini++) @@ -2152,9 +2174,12 @@ Bool_t Run::plotClusterDistribution(HistogramType* HistogramTypepointer) return 0; } +TCanvas* Run::plot1DHistogram(TH1FO* onehistogram, TString fitFuncType, Bool_t verbose, Bool_t logscale, Bool_t withOverflow, Float_t fitstart, TString titlestr, TString legendstr ) { + HistogramType* HistogramTypepointer = onehistogram->itsHistogramType; + return plot1DHistogram(HistogramTypepointer, onehistogram, fitFuncType, verbose, logscale, withOverflow, fitstart, titlestr, legendstr ); +} - -TCanvas* Run::plot1DHistogram(HistogramType* HistogramTypepointer, TH1F* onehistogram, TString fitFuncType, Bool_t verbose, Bool_t logscale, Bool_t withOverflow, Float_t fitstart, TString titlestr, TString legendstr ) +TCanvas* Run::plot1DHistogram(HistogramType* HistogramTypepointer, TH1FO* onehistogram, TString fitFuncType, Bool_t verbose, Bool_t logscale, Bool_t withOverflow, Float_t fitstart, TString titlestr, TString legendstr ) { if (onehistogram != nullptr) { @@ -2172,7 +2197,7 @@ TCanvas* Run::plot1DHistogram(HistogramType* HistogramTypepointer, TH1F* onehist gPad->SetLogy(1); } if (withOverflow) - onehistogram=ShowOverflow(onehistogram); + onehistogram=(TH1FO*)ShowOverflow((TH1F*)onehistogram); onehistogram->Draw(); Double_t* parameters = (Double_t *)calloc(11, sizeof(Double_t)); // allocate 10 parameters for safety, maximum 9 used at the moment diff --git a/MABS_run_analyzer/Run.h b/MABS_run_analyzer/Run.h index 1703891..7457e8c 100644 --- a/MABS_run_analyzer/Run.h +++ b/MABS_run_analyzer/Run.h @@ -58,8 +58,8 @@ class Run; * to the databases and stores the analyzed data as histograms. * */ -class Run { +class Run { private: /// Object to connect to the database located at jspc29 @@ -94,7 +94,7 @@ private: /** * @brief takes a float value and a precision and appends it to an SQL statement used to update the database */ - void constructUpdateString(string*, const string, const Double_t, int); + void constructUpdateString(string *sqlupdatequery, const string databasevaluename, const Double_t value, const int precision, const Double_t minval, const Double_t maxval, const TString defaultval); /** * @brief writes values back to the SQL database @@ -195,6 +195,13 @@ private: Bool_t findRotatedClusterDimension(Float_t rotateangle = 45); void initRootParameters(); + + /** + * @brief Fills typical histograms in HistomTyoe + * + * + */ + void fillAHistogramsinclass(HistogramType* histogramtypepointer, Int_t hiti, Int_t completeclustersize, Float_t pixelSum, Float_t notSeedSum, Float_t* pt_a_pixelSum, Float_t* pt_a_notSeedSum); /** @brief set this variable to true, to use the dynamical noise calculation in MAPS:analyzeRun() */ Bool_t dynamicalNoise = 1; @@ -212,7 +219,7 @@ private: Bool_t normalizeHistogramClasses(); /** @brief finds and and masks RTS pixel if maskRTSpixel isset to true */ - Bool_t FindRTSPixelToMask(HistogramType* histogramtypepointer); + Bool_t FindRTSPixelToMask(HistogramType* histogramtypepointer, Bool_t verbose = kFALSE); public: /** @brief empty constructor */ @@ -311,6 +318,7 @@ public: Bool_t setFixedThresholdValueElectrons(Float_t); /** @brief set a threshold value in ADU to use in the histogram Run::histogramfixedthreshold, this value is used in @c MAPS::initMapsRun() */ + Bool_t setFixedThresholdValueADU(HistogramType*, Float_t); Bool_t setFixedThresholdValueADU(Float_t); /** @@ -416,7 +424,8 @@ public: */ Bool_t plotClusterDistribution(HistogramType*); - TCanvas* plot1DHistogram(HistogramType* histogramtypepointer, TH1F* onehistogram, TString fitFuncType = "landau", Bool_t verbose = false, Bool_t logscale =false, Bool_t withOverflow=false, Float_t fitstart=-1, TString titlestr = "", TString legendstr = ""); + TCanvas* plot1DHistogram(HistogramType* histogramtypepointer, TH1FO* onehistogram, TString fitFuncType = "landau", Bool_t verbose = false, Bool_t logscale =false, Bool_t withOverflow=false, Float_t fitstart=-1, TString titlestr = "", TString legendstr = ""); + TCanvas* plot1DHistogram(TH1FO* onehistogram, TString fitFuncType = "landau", Bool_t verbose = false, Bool_t logscale =false, Bool_t withOverflow=false, Float_t fitstart=-1, TString titlestr = "", TString legendstr = ""); pixelinfo pixelinfoMi34[32]; @@ -471,4 +480,5 @@ public: /// sensor information to use in analysis, is the system read out by USB or PXI? Number of rows differ sensorinfostruct cursensorinfo; }; + #endif diff --git a/MABS_run_analyzer/help.h b/MABS_run_analyzer/help.h index a396ae7..d0aae5c 100644 --- a/MABS_run_analyzer/help.h +++ b/MABS_run_analyzer/help.h @@ -24,6 +24,7 @@ #define MAXHITS 1000000 #define MAXPIXELS 100000 + //#################################################################### /** * @file help.h -- 2.43.0