From 6b433dd1da6acfce55b863cdc4ee35bb9122a115 Mon Sep 17 00:00:00 2001 From: Benjamin Linnik Date: Mon, 18 Sep 2017 13:03:26 +0200 Subject: [PATCH] Analyzer: adjustments for Ali --- MABS_run_analyzer/ChargeSpektrum.c | 140 +++++++++++- MABS_run_analyzer/ChargeSpektrumFunctions.c | 209 +++++++++++------ MABS_run_analyzer/HistogramType.c | 238 +++++++++++++------- MABS_run_analyzer/HistogramType.h | 23 +- MABS_run_analyzer/Run.c | 49 +++- MABS_run_analyzer/help.h | 15 +- 6 files changed, 495 insertions(+), 179 deletions(-) diff --git a/MABS_run_analyzer/ChargeSpektrum.c b/MABS_run_analyzer/ChargeSpektrum.c index b3c77ec..19e523b 100644 --- a/MABS_run_analyzer/ChargeSpektrum.c +++ b/MABS_run_analyzer/ChargeSpektrum.c @@ -40,6 +40,12 @@ void ChargeSpektrum(TString runnumber = "") // set file path and other variables to be set Init(); cout << "Found " << numberRuns << " run(s) in 'runlist.txt' or given as parameters." << endl; + TNtupleO *datainput2; + datainput2 =new TNtupleO("Position Veto", "data", "x:y:xerr:yerr"); + TNtupleO *datainput3; + datainput3 =new TNtupleO("ntuple3", "data", "x:y:xerr:yerr"); + TNtupleO *datainput4; + datainput4 =new TNtupleO("ntuple4", "data", "x:y:xerr:yerr"); for(Int_t runi=0; runi0) @@ -95,21 +101,35 @@ void ChargeSpektrum(TString runnumber = "") runs[runi]->setLabel(runListCustomTitle[runi]); } } - compareHistogramClassVectorMABSBot.push_back(runs[runi]->histogram); + compareHistogramClassVectorMABSBot.push_back(runs[runi]->histogramwoRTS); // gROOT->SetBatch(kTRUE); if (!isBatch) gROOT->SetBatch(kFALSE); // Uncomment below to do classical analysis withour RTS + compareHistogramClassVector.push_back(runs[runi]->histogramwoRTS->normalized); compareHistogramClassVector.push_back(runs[runi]->histogram->normalized); - compareHistogramClassVector.push_back(runs[runi]->histogramfixedthreshold->normalized); + compareHistogramClassVector2.push_back(runs[runi]->histogramwoRTS->normalized); - compareHistogramClassVector2.push_back(runs[runi]->histogram->normalized); - compareHistogramVector.push_back(runs[runi]->histogram->normalized->Seed); - compareHistogramVector.push_back(runs[runi]->histogram->normalized->Sum); +// compareHistogramVector.push_back(runs[runi]->histogramfixedthreshold->SeedPerc); + + compareHistogramVector2.push_back(runs[runi]->histogramwoRTS->normalized->Seed); + compareHistogramVector3.push_back(runs[runi]->histogramwoRTS->normalized->Sum); + + + runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS->pixeltimefiredDistrib); + runs[runi]->plot1DHistogram(runs[runi]->histogramwoRTS->LeakageCurrentInPixelSorted); + +// compareHistogramClassVector2.push_back(runs[runi]->histogram->normalized); +// +// compareHistogramVector3.push_back(runs[runi]->histogram->normalized->Seed); +// compareHistogramVector4.push_back(runs[runi]->histogram->normalized->Sum); +// +// +// compareHistogramVector5.push_back(runs[runi]->histogramfixedthreshold->calibrated->Seed); // compareHistogramClassVector2.push_back(runs[runi]->histogramwoRTS); // compareHistogramClassVector2.push_back(runs[runi]->histogram->normalized); @@ -120,10 +140,90 @@ void ChargeSpektrum(TString runnumber = "") // 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); +// runs[runi]->plot1DHistogram( runs[runi]->histogram->Seed, "gaus", true, false, false, runs[runi]->histogram->fixedThresholdValue); + - ntuple->Fill(runs[runi]->labbook.depletionV,runs[runi]->histogramfixedthreshold->posSeed,0,0); +// Double_t AverageCharge[9]; +// Double_t TotalHits[9]; +// for (UInt_t sumi = 0; sumi < 9; sumi++) { +// AverageCharge[sumi] = 0; +// TotalHits[sumi] = 0; +// for (UInt_t bini = 1; bini <= runs[runi]->histogramfixedthreshold->a_Sum[sumi]->GetXaxis()->GetNbins();++bini) { +// AverageCharge[sumi] += runs[runi]->histogramfixedthreshold->a_Sum[sumi]->GetBinContent(bini)*runs[runi]->histogramfixedthreshold->a_Sum[sumi]->GetBinCenter(bini); +// TotalHits[sumi] += runs[runi]->histogramfixedthreshold->a_Sum[sumi]->GetBinContent(bini); +// } +// AverageCharge[sumi] /= TotalHits[sumi]/runs[runi]->histogramfixedthreshold->gain; +// } + + + + +// TNtupleO *datainput; +// datainput =new TNtupleO("Position Seed", "data", "x:y:xerr:yerr"); +// datainput->itsHistogramType = runs[runi]->histogramfixedthreshold; +// datainput->xaxisTitle = "Number of pixel in cluster"; +// datainput->yaxisTitle = "Position Seed"; +// for (UInt_t sumi = 0; sumi < 9; sumi++) { +// // if (!(sumi%2)) +// // compareHistogramVector2.push_back(runs[runi]->histogram->normalized->a_Sum[sumi]); +// datainput->Fill(sumi+1,runs[runi]->histogramfixedthreshold->a_posSum[sumi],0,0); +// // runs[runi]->plot1DHistogram( runs[runi]->histogram->a_Sum[sumi], "gaus", true, false, false, runs[runi]->histogram->fixedThresholdValue); +// } +// ntupleVector.push_back(datainput); +// +// +// datainput2->itsHistogramType = runs[runi]->histogramfixedthreshold; +// datainput2->xaxisTitle = "Voltage"; +// datainput2->yaxisTitle = "Position Veto"; +// datainput2->Fill(runs[runi]->labbook.depletionV,runs[runi]->histogramfixedthreshold->posVeto,0,0); +// +// +// +// datainput3->itsHistogramType = runs[runi]->histogramfixedthreshold; +// datainput3->xaxisTitle = "Voltage"; +// datainput3->yaxisTitle = "Position Sum"; +// // datainput3->Fill(runs[runi]->labbook.depletionV,runs[runi]->histogramfixedthreshold->posSum,0,0); +// +// +// +// +// datainput4->itsHistogramType = runs[runi]->histogramfixedthreshold; +// datainput4->xaxisTitle = "Voltage"; +// datainput4->yaxisTitle = "Cluster multiplicity"; +// datainput4->Fill(runs[runi]->labbook.depletionV,runs[runi]->histogramfixedthreshold->clustermultiplicity,0,0); + + +// TNtupleO *datainput; +// datainput =new TNtupleO("ntuple", "data", "x:y:xerr:yerr"); +// datainput->itsHistogramType = runs[runi]->histogram; +// datainput->xaxisTitle = "Number of pixel in cluster"; +// datainput->yaxisTitle = "Average charge collected [e]"; +// for (UInt_t sumi = 0; sumi < 9; sumi++) { +// if (!(sumi%2)) +// compareHistogramVector2.push_back(runs[runi]->histogram->normalized->a_Sum[sumi]); +// datainput->Fill(sumi+1,AverageCharge[sumi],0,0); +// // runs[runi]->plot1DHistogram( runs[runi]->histogram->a_Sum[sumi], "gaus", true, false, false, runs[runi]->histogram->fixedThresholdValue); +// } +// ntupleVector.push_back(datainput); + + + +// TNtupleO *datainput; +// datainput =new TNtupleO("ntuple", "data", "x:y:xerr:yerr"); +// datainput->itsHistogramType = runs[runi]->histogram; +// datainput->xaxisTitle = "Number of pixel in cluster"; +// datainput->yaxisTitle = "Average charge collected [ADU]"; +// for (UInt_t sumi = 0; sumi < 9; sumi++) { +// if (!(sumi%2)) +// compareHistogramVector2.push_back(runs[runi]->histogram->normalized->a_Sum[sumi]); +// datainput->Fill(sumi+1,runs[runi]->histogram->a_integralSum[sumi],0,runs[runi]->histogram->a_integralSumErr[sumi]); +// runs[runi]->plot1DHistogram( runs[runi]->histogram->a_Sum[sumi], "gaus", true, false, false, runs[runi]->histogram->fixedThresholdValue); +// } +// ntupleVector.push_back(datainput); + +// ntuple.Clear(); + // ntuple->Fill(2,2,1,1); // ntuple->Fill(3,3,1,1); // ntuple->Fill(4,4,1,1); @@ -136,7 +236,17 @@ void ChargeSpektrum(TString runnumber = "") // 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); +// Float_t integratestart = 0; +// runs[runi]->plot1DHistogram( runs[runi]->histogram->normalized->Seed, "landau", true, false, false); +// integratestart = runs[runi]->histogram->FindNoisethresholdborder(runs[runi]->histogram->Sum, false, false); +// cout << "integratestart: " << integratestart << endlr; +// runs[runi]->plot1DHistogram( runs[runi]->histogram->normalized->Sum, "landau", true, false, false, integratestart); +// integratestart = runs[runi]->histogram->FindNoisethresholdborder(runs[runi]->histogram->a_Sum[2], false, false); +// cout << "integratestart: " << integratestart << endlr; +// runs[runi]->plot1DHistogram( runs[runi]->histogram->normalized->a_Sum[2], "landau", true, false, false, integratestart); +// integratestart = runs[runi]->histogram->FindNoisethresholdborder(runs[runi]->histogram->a_Sum[5], false, false); +// cout << "integratestart: " << integratestart << endlr; +// runs[runi]->plot1DHistogram( runs[runi]->histogram->normalized->a_Sum[5], "landau", true, false, false, integratestart); //runs[runi]->plot1DHistogram( runs[runi]->histogramwoRTS->normalized->Seed, "GaussTail"); //runs[runi]->plot1DHistogram( runs[runi]->histogramwoRTS->a_Sum[0]); // compareHistogramVector.push_back(runs[runi]->histogramwoRTS->a_Sum[0]); // Seed @@ -165,11 +275,19 @@ void ChargeSpektrum(TString runnumber = "") // plot and clear compareHistogramClassVector vector, delete these two lines if you want to compare among different runs plotAllRuns(&compareHistogramClassVector); compareHistogramClassVector.clear(); + + CompareHistograms(&compareHistogramVector); + compareHistogramVector.clear(); } //cout << runs[runi]->histogram } } - } + } // loop over all runs end + + ntupleVector2.push_back(datainput2); + ntupleVector2.push_back(datainput3); + + ntupleVector3.push_back(datainput4); @@ -199,7 +317,9 @@ plotAllRuns(&compareHistogramClassVector2); plotAllRuns(&compareHistogramClassVector3); plotAllRuns(&compareHistogramClassVector4); plotAllRuns(&compareHistogramClassVector5); -plot1DData(ntuple); +plot1DData(&ntupleVector); +plot1DData(&ntupleVector2); +plot1DData(&ntupleVector3); CompareLeageCurrent(&compareHistogramClassVector4); CompareLeageCurrent(&compareHistogramClassVector5); writeObservableToFile(); diff --git a/MABS_run_analyzer/ChargeSpektrumFunctions.c b/MABS_run_analyzer/ChargeSpektrumFunctions.c index c7ef5e1..ff13df5 100644 --- a/MABS_run_analyzer/ChargeSpektrumFunctions.c +++ b/MABS_run_analyzer/ChargeSpektrumFunctions.c @@ -16,10 +16,10 @@ #include #include #include +#include "help.h" #include "HistogramType.h" #include "Run.h" -#include "help.h" using namespace std; Int_t* ReadRunList(); @@ -45,9 +45,6 @@ Bool_t FindGoodTitle(vector* ptCompareHistogramClassVector, vect Bool_t FindGoodTitle(vector* ptCompareHistogramVector); Bool_t FindGoodTitle(); -/** @brief A function to write data from TNtuple to a file. File is save with the prefix "Oberservable_" in the result folder - * */ -Bool_t writeNTupleToFile(TNtuple* ntuplepointer); /** @brief A function to plot TH1F histograms of different runs into one file @@ -80,7 +77,11 @@ void Init(); /** @brief A function to plot ntuplöe data into a canvas * */ -TCanvas* plot1DData(TNtuple* ntuplepointer, Bool_t verbose=0, Bool_t logscale=0, TString titlestr="", TString legendstr="" ); +TCanvas* plot1DData(std::vector* ntuplepvectorpointer, Bool_t verbose=0, Bool_t logscale=0, TString titlestr="", TString legendstr="" ); + +/** @brief A function to write data from TNtuple to a file. File is save with the prefix "Oberservable_" in the result folder + * */ +Bool_t writeNTupleVectorToFile(std::vector* ntuplepvectorpointer); Bool_t writeObservableToFile(); /** @brief A vector able to hold pointer of type #Run::histogramstruct @@ -99,7 +100,7 @@ Bool_t testifMixingCalibration(vector*); 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; -TNtuple *ntuple, *ntuple2, *ntuple3, *ntuple4; +vector ntupleVector, ntupleVector2, ntupleVector3, ntupleVector4; TString ownpath = ""; TString savepathresults; @@ -127,10 +128,6 @@ void Init() //FindGoodTitle(); TTimeStamp timestamp = TTimeStamp(); savepathresults = Form("./results/%d%06d/", (int)timestamp.GetDate(kFALSE), (int)timestamp.GetTime(kFALSE)); - ntuple = new TNtuple("ntuple", "data", "x:y:xerr:yerr"); - ntuple2 = new TNtuple("ntuple2","data2","x:y:xerr:yerr"); - ntuple3 = new TNtuple("ntuple3","data3","x:y:xerr:yerr"); - ntuple4 = new TNtuple("ntuple4","data4","x:y:xerr:yerr"); } @@ -267,67 +264,111 @@ Int_t ReadRunList(std::vector* runlist) -TCanvas* plot1DData(TNtuple* ntuplepointer, Bool_t verbose, Bool_t logscale, TString titlestr, TString legendstr ) +TCanvas* +plot1DData(std::vector* ntuplepvectorpointer, Bool_t verbose, Bool_t logscale, TString titlestr, TString legendstr ) { - if (ntuplepointer != nullptr) + if (ntuplepvectorpointer != nullptr) { - if (ntuplepointer->GetEntries() > 0) { - TString canvastitle = Form("%s", ntuplepointer->GetName()); - TString canvasname = Form("%s", ntuplepointer->GetName()); - TCanvas* canvas = new TCanvas(canvasname, canvastitle, 900, 700); - canvas->SetGrid(); - - if (titlestr.Length()>0) - ntuplepointer->SetTitle(titlestr); - if (logscale) { - gPad->SetLogy(1); + if (ntuplepvectorpointer->size() > 0) { + TNtupleO* ntuplepointer = ntuplepvectorpointer->at(0); + if (ntuplepointer->GetEntries() > 0) { + // TNtuple* ntuplepointer = &(ntuplepvectorpointer->at(0)); + TString canvastitle = Form("%s", ntuplepointer->GetName()); + TString canvasname = Form("%s", ntuplepointer->GetName()); + // titlestr = ntuplepointer->GetTitle(); + + if (titlestr.Length()>0) + ntuplepointer->SetTitle(titlestr); + if (logscale) { + gPad->SetLogy(1); + } + + TCanvas* canvas = new TCanvas(canvasname, canvastitle, 900, 700); + canvas->SetGrid(); + + for(UInt_t runi=0;runisize();runi++) // loop over runs or other stacked elements in vector + { + // TNtuple* ntuplepointer = &(ntuplepvectorpointer->at(runi)); + TNtupleO* ntuplepointer = ntuplepvectorpointer->at(runi); + ntuplepointer->Draw("x:y:xerr:yerr","","goff"); + + TGraphErrors *gr = new TGraphErrors(ntuplepointer->GetEntries(),ntuplepointer->GetV1(),ntuplepointer->GetV2(),ntuplepointer->GetV3(),ntuplepointer->GetV4()); + if ( ntuplepointer->Title.Length() > 0 ) { gr->SetTitle(ntuplepointer->Title); } else { gr->SetTitle("");} + if ( ntuplepointer->xaxisTitle.Length() > 0 ) { gr->GetXaxis()->SetTitle(ntuplepointer->xaxisTitle); } else { gr->GetXaxis()->SetTitle("");} + if ( ntuplepointer->yaxisTitle.Length() > 0 ) { gr->GetYaxis()->SetTitle(ntuplepointer->yaxisTitle); } else { gr->GetYaxis()->SetTitle("");} + if (ntuplepointer->itsHistogramType != 0) gr->SetMarkerColor(ntuplepointer->itsHistogramType->color); else gr->SetMarkerColor(runi); gr->SetMarkerStyle(21); + if (!runi) gr->Draw("ALP"); else gr->Draw("LP"); + } + + TPaveText *owntitle = new TPaveText(0.15,0.9,0.930401,0.995,"nbNDC"); + owntitle->SetFillStyle(0); + owntitle->SetBorderSize(0); + owntitle->SetTextAlign(22); + owntitle->SetTextSize(0.05); + owntitle->SetTextFont(42); + + owntitle->Clear(); + if (titlestr.Length()>0) { + owntitle->AddText(titlestr); + owntitle->Draw("SAME"); + } + + canvas->Update(); + // canvas -> SaveAs( savepathresults + "/" + runcode + " " + onehistogram->GetName() + ".eps"); + + TImageDump *img = new TImageDump(savepathresults + ntuplepointer->GetName() + ".png"); + canvas->Paint(); + img->Close(); + + writeNTupleVectorToFile(ntuplepvectorpointer); + return canvas; } - - ntuplepointer->Draw("x:y:xerr:yerr","","goff"); - - TGraphErrors *gr = new TGraphErrors(ntuplepointer->GetEntries(),ntuplepointer->GetV1(),ntuplepointer->GetV2(),ntuplepointer->GetV3(),ntuplepointer->GetV4()); - gr->SetTitle("TGraphErrors Example"); - gr->SetMarkerColor(4); - gr->SetMarkerStyle(21); - gr->Draw("ALP"); - - canvas->Update(); - // canvas -> SaveAs( savepathresults + "/" + runcode + " " + onehistogram->GetName() + ".eps"); - - TImageDump *img = new TImageDump(savepathresults + ntuplepointer->GetName() + ".png"); - canvas->Paint(); - img->Close(); - - writeNTupleToFile(ntuplepointer); - return canvas; } } return nullptr; } -Bool_t writeNTupleToFile(TNtuple* ntuplepointer) +Bool_t writeNTupleVectorToFile(std::vector* ntuplepvectorpointer) { - if (ntuplepointer != nullptr) + if (ntuplepvectorpointer != nullptr) { - system("mkdir "+ savepathresults + " -p"); - - TString filename = savepathresults + Form("Observable_%s.dat", ntuplepointer->GetName()); - - TString header = "x\ty\txerr\tyerr"; - fstream* fout = new fstream(filename,ios::out); - - *fout << header << endl; - TString outline; - for(int i=0;iGetEntries();i++) { - ntuplepointer->GetEntry(i); - Float_t *p = ntuplepointer->GetArgs(); - outline =""; - outline+=Form("%.2f\t%.2f\t%.2f\t%.2f", p[0], p[1], p[2], p[3]); - *fout<size() > 0) { + TNtupleO* ntuplepointer = ntuplepvectorpointer->at(0); + system("mkdir "+ savepathresults + " -p"); + TString filename = savepathresults + Form("Observable_%s.dat", ntuplepointer->GetName()); + + TString header = ""; + for(UInt_t runi=0;runisize();runi++) // loop over runs or other stacked elements in vector + { + header += "x\ty\txerr\tyerr\t"; + } + fstream* fout = new fstream(filename,ios::out); + *fout << header << endl; + + // find max entries + Int_t maxEntries = 0; + for(UInt_t runi=0;runisize();runi++) { + ntuplepointer = ntuplepvectorpointer->at(runi); + maxEntries = (ntuplepointer->GetEntries() > maxEntries)?ntuplepointer->GetEntries():maxEntries; + } + + + TString outline=""; + for(int i=0;i < maxEntries;i++) { + for(UInt_t runi=0;runisize();runi++) // loop over runs or other stacked elements in vector + { + ntuplepointer = ntuplepvectorpointer->at(runi); + ntuplepointer->GetEntry(i); + Float_t *p = ntuplepointer->GetArgs(); + outline+=Form("%.2f\t%.2f\t%.2f\t%.2f\t", p[0], p[1], p[2], p[3]); + } + *fout<close(); + return 0; } - fout->close(); - return 0; } return 1; } @@ -436,7 +477,8 @@ Bool_t CompareHistograms(vector* ptCompareHistogramVector, TString title // legend entries Float_t height = ptCompareHistogramVector->size() * 0.03; - TLegend* leg1 = new TLegend(0.3,0.89-height,0.95,0.89);//(0.6,0.7,0.89,0.89); +// TLegend* leg1 = new TLegend(0.3,0.89-height,0.95,0.89);//(0.6,0.7,0.89,0.89); + TLegend* leg1 = new TLegend(0.3,0.89-height,0.6,0.89);//(0.6,0.7,0.89,0.89); leg1->SetTextSize(0.025); leg1->SetFillStyle(1001); leg1->SetTextFont(132); @@ -462,7 +504,7 @@ Bool_t CompareHistograms(vector* ptCompareHistogramVector, TString title { TH1F* curhistogramclone = (TH1F*) ptCompareHistogramVector->at(histogrami)->Clone(); heighestval1 = (curhistogramclone->GetMaximum()>heighestval1?curhistogramclone->GetMaximum():heighestval1); - lastbin1 = (curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,1))>lastbin1)?curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,1)):lastbin1; + lastbin1 = (curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,2))>lastbin1)?curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,2)):lastbin1; canvastitle+= Form("_%s",getRunnumberAtBegin(curhistogramclone->GetName()).Data()); } TCanvas* canvas = new TCanvas(canvasname, canvastitle, 1200, 800); @@ -489,7 +531,7 @@ Bool_t CompareHistograms(vector* ptCompareHistogramVector, TString title leg1->Draw("SAME"); curhistogramclone->SetAxisRange(0,lastbin1*1.1,"X"); - // curhistogramclone->GetYaxis()->SetRangeUser(5,heighestval1*4); + curhistogramclone->GetYaxis()->SetRangeUser(1,heighestval1*1.4); owntitle->Clear(); owntitle->AddText(trimRunnumberAtBegin(curhistogramclone->GetName())); @@ -500,8 +542,8 @@ Bool_t CompareHistograms(vector* ptCompareHistogramVector, TString title } owntitle->Draw("SAME"); - gPad->SetLogy(1); - curhistogramclone->GetYaxis()->UnZoom(); +// gPad->SetLogy(1); +// curhistogramclone->GetYaxis()->UnZoom(); } canvas->Update(); // gPad->Modified(); gPad->Update(); @@ -759,7 +801,7 @@ Bool_t plotAllRuns(vector* ptCompareHistogramClassVector) curhistogramclone->GetXaxis()->UnZoom(); curhistogramclone->GetXaxis()->SetRange(curhistogramclone->GetXaxis()->FindBin(curhistogramclassp->noisethresholdborder),curhistogramclone->GetXaxis()->FindBin(posMaxValHist)); // look only for maxima with x greater than noiseborder, cut away noise heighestval1 = (curhistogramclone->GetMaximum()>heighestval1?curhistogramclone->GetMaximum():heighestval1); - lastbin1 = (curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,1))>lastbin1)?curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,1)):lastbin1; + lastbin1 = (curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,2))>lastbin1)?curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,2)):lastbin1; // cout << "Last bin: " << colorcyan << lastbin1 << endlr; // cout << "curhistogramclone->GetMaximum(): " << colorcyan << curhistogramclone->GetMaximum() << endlr; @@ -770,14 +812,14 @@ Bool_t plotAllRuns(vector* ptCompareHistogramClassVector) heighestval2 = (curhistogramclone->GetMaximum()>heighestval2?curhistogramclone->GetMaximum():heighestval2); - lastbin2 = (curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,1))>lastbin2)?curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,1)):lastbin2; + lastbin2 = (curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,2))>lastbin2)?curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,2)):lastbin2; curhistogramclone = (TH1F*) curhistogramclassp->Veto->Clone(); posMaxValHist = curhistogramclone->GetXaxis()->GetXmax(); curhistogramclone->GetXaxis()->UnZoom(); curhistogramclone->GetXaxis()->SetRange(curhistogramclone->GetXaxis()->FindBin(curhistogramclassp->noisethresholdborder),curhistogramclone->GetXaxis()->FindBin(posMaxValHist)); // look only for maxima with x greater than noiseborder, cut away noise heighestval3 = (curhistogramclone->GetMaximum()>heighestval3?curhistogramclone->GetMaximum():heighestval3); - lastbin3 = (curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,1))>lastbin3)?curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,1)):lastbin3; + lastbin3 = (curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,2))>lastbin3)?curhistogramclone->GetBinCenter(curhistogramclone->FindLastBinAbove(2,2)):lastbin3; // cout << "heighestval3: " << colorcyan << heighestval3 << endlr; // cout << "curhistogramclone->GetMaximum(): " << colorcyan << curhistogramclone->GetMaximum() << endlr; @@ -986,7 +1028,7 @@ Bool_t FindGoodTitle(vector* ptCompareHistogramClassVector, vect 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", 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); if (curlabbook.depletionV > 0) @@ -1121,8 +1163,35 @@ Bool_t FindGoodTitle(vector* ptCompareHistogramClassVector, vect // cout << colorred << legendstr << endlr; } - // Folder name suffix - TString folderadd = ""; +// // Folder name suffix +// TString folderadd = ""; +// if (same_RunNumber) +// folderadd.Append(Form(", %d", firstlabbook.runnumber)); +// if (same_HistType && mayBeSameRun) +// folderadd.Append(Form(", %s", trimRunnumberAtBegin(firsthistogramp->GetName()).Data())); +// if (same_Source) +// folderadd.Append(Form(", %s", firstlabbook.source.Data())); +// if (same_ChipGen) +// folderadd.Append(Form(", %s", firstlabbook.chipGen.Data())); +// if (same_Matrix) +// folderadd.Append(Form(", %s, %.0fx%.0f #mum^{2} pitch, %s", firstlabbook.matrix.Data(), firstpixelinfo.pitchX, firstpixelinfo.pitchY, firstpixelinfo.comment.Data())); +// title1.Append(Form("\n")); +// if (same_Temp) +// title1.Append(Form("T=%.0f {}^{o}C", firstlabbook.tempSens)); +// if (same_IonRad && firstlabbook.radDoseIon != 0) +// title1.Append(Form(", %.1f MRad", firstlabbook.radDoseIon)); +// if (same_NonIonRad && firstlabbook.radDoseNonIon != 0) +// title1.Append(Form(", %.1f*10^{13} n_{eq}/cm^{2}", firstlabbook.radDoseNonIon)); +// if (same_Clock) +// title1.Append(Form(", %.2f Mhz", firstlabbook.clock)); +// if (same_Depletion && firstlabbook.depletionV >= 0) +// title1.Append(Form(", U_{dep}=%.1f V", firstlabbook.depletionV)); +// title1 = title1.Remove(0,2); // remove trailing " ," + + + + + HistogramType* curhistogramclassp = ptCompareHistogramClassVector->at(0); labbooksctruct curlabbook = curhistogramclassp->run->labbook; pixelinfo curpixelinfo = curhistogramclassp->run->curpixelinfo; diff --git a/MABS_run_analyzer/HistogramType.c b/MABS_run_analyzer/HistogramType.c index 24de8dd..4f3c132 100644 --- a/MABS_run_analyzer/HistogramType.c +++ b/MABS_run_analyzer/HistogramType.c @@ -137,8 +137,10 @@ 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]); + for (unsigned int i=0; i < a_Sum.size(); ++i) { + calibrated->a_Sum.push_back(a_Sum[i]); + calibrateHistogram(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; @@ -149,19 +151,27 @@ Bool_t HistogramType::calibrateHistograms( Float_t gotgain ) { if (LeakageCurrentInPixelSorted != 0) calibrateHistogramYAxis(calibrated->LeakageCurrentInPixelSorted, LeakageCurrentInPixelSorted); calibrated->posSeed = posSeed * gain; - calibrated->posSum = posSum * gain; + calibrated->posSum = posSum * gain; + for (unsigned int i=0; i < a_posSum.size(); ++i){ + calibrated->a_posSum.push_back(a_posSum[i] * gain); + } calibrated->posVeto = posVeto * gain; calibrated->integralSeed = integralSeed; calibrated->integralSeedErr = integralSeedErr; calibrated->integralVeto = integralVeto; calibrated->integralSum = integralVeto; - calibrated->integralSumErr = integralSumErr; + calibrated->integralSumErr = integralSumErr; + for (unsigned int i=0; i < a_integralSum.size(); ++i){ + calibrated->a_integralSum.push_back( a_integralSum[i]); + calibrated->a_integralSumErr.push_back( a_integralSumErr[i] ); + } calibrated->posSeedPerc = posSeedPerc; calibrated->sigmaSeedPerc = sigmaSeedPerc; calibrated->avgNoise = avgNoise * gain; calibrated->avgNoisePlus = avgNoisePlus * gain; calibrated->avgNoiseMinus = avgNoiseMinus * gain; calibrated->sr90IntegralVal = sr90IntegralVal; + calibrated->sr90IntegralErr = sr90IntegralErr; calibrated->StoN = StoN; calibrated->CCE_in_Perc_1 = CCE_in_Perc_1; calibrated->CCE_in_Perc_25 = CCE_in_Perc_25; @@ -244,6 +254,7 @@ Bool_t HistogramType::normalizeHistograms( Int_t got_frames_found ) { if (Seed != 0) normalizeHistogram(normalized->Seed, Seed); if (Sum != 0) normalizeHistogram(normalized->Sum, Sum); for (unsigned int i=0; i < a_Sum.size(); ++i){ + normalized->a_Sum.push_back(a_Sum[i]); normalizeHistogram(normalized->a_Sum[i], a_Sum[i]); } if (Veto != 0) normalizeHistogram(normalized->Veto, Veto); @@ -257,12 +268,19 @@ Bool_t HistogramType::normalizeHistograms( Int_t got_frames_found ) { normalized->frames_found = frames_found; normalized->posSeed = posSeed; normalized->posSum = posSum; + for (unsigned int i=0; i < a_posSum.size(); ++i){ + normalized->a_posSum.push_back(a_posSum[i]); + } normalized->posVeto = posVeto; normalized->integralSeed = integralSeed/frames_found*pow10(6); normalized->integralSeedErr = integralSeedErr * sqrt(pow10(6)/frames_found); normalized->integralVeto= integralVeto/frames_found*pow10(6); normalized->integralVeto= integralVeto/frames_found*pow10(6); normalized->integralSum = integralSum/frames_found*pow10(6); + for (unsigned int i=0; i < a_integralSum.size(); ++i){ + normalized->a_integralSum.push_back(a_integralSum[i]/frames_found*pow10(6)); + normalized->a_integralSumErr.push_back(a_integralSumErr[i]* sqrt(pow10(6)/frames_found)); + } normalized->integralSumErr = integralSumErr * sqrt(pow10(6)/frames_found); normalized->posSeedPerc = posSeedPerc; normalized->sigmaSeedPerc = sigmaSeedPerc; @@ -271,7 +289,8 @@ Bool_t HistogramType::normalizeHistograms( Int_t got_frames_found ) { normalized->avgNoise = avgNoise; normalized->avgNoisePlus = avgNoisePlus; normalized->avgNoiseMinus = avgNoiseMinus; - normalized->sr90IntegralVal = sr90IntegralVal; + normalized->sr90IntegralVal = sr90IntegralVal/frames_found*pow10(6);; + normalized->sr90IntegralErr = sr90IntegralErr * sqrt(pow10(6)/frames_found); normalized->StoN = StoN; normalized->CCE_in_Perc_1 = CCE_in_Perc_1; normalized->CCE_in_Perc_25 = CCE_in_Perc_25; @@ -333,6 +352,7 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType else if (noisethresholdborder>=0 ) noiseborder = noisethresholdborder; +// cout << colorred << " noiseborder " << noiseborder << " : " << endlr; // cout << colorred << " " << histogrampointer->GetName() << " : " << endlr; @@ -565,8 +585,21 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType // leg->Draw(); // } - else if (fitFuncType=="landau") + else if (fitFuncType=="landauSIMPLE") { + TF1* fitFunc = new TF1("fitFunc",fitFuncType,noiseborder,posMaxValHist); + fit_result_ptr = histogrampointer->Fit(fitFunc, "N,Q,W,S", "", noiseborder, posMaxValHist); + parameters = (Double_t *)calloc(3, sizeof(Double_t)); + for (Int_t pari=0; pari<3; pari++) + parameters[pari]=fitFunc->GetParameter(pari); + fitFunc->SetLineColor(kGreen); + fitFunc->SetLineStyle(1); // normal for the following fits + if (verbose) + fitFunc->DrawCopy("same"); + + } else if (fitFuncType=="landau") + { + fitFuncType = "landau"; if (posVeto > 0) posMaxValHist = posVeto-posVeto/8; histogrampointer->GetXaxis()->SetRange(histogrampointer->GetXaxis()->FindBin(noiseborder),histogrampointer->GetXaxis()->FindBin(posMaxValHist)); // look only for maxima with x greater than noiseborder, cut away noise @@ -581,6 +614,11 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType fit_result_ptr = histogrampointer->Fit(fitFunc, "N,Q,W,S", "", noiseborder, posMaxValHist); for (Int_t pari=0; pari<3; pari++) parameters[pari]=fitFunc->GetParameter(pari); + + fitFunc->SetLineColor(kGreen); + fitFunc->SetLineStyle(1); // normal for the following fits + if (verbose) + fitFunc->DrawCopy("same"); // posMax = fitFunc->GetParameter(1); // new method, just use MPV of fit, as noise border is well known. Compare with old method anyway and warn if something suspecious. // parameters[1]=fitFunc->GetMaximumX(); posMax = fitFunc->GetMaximumX(); // new method, just use MPV of fit, as noise border is well known. Compare with old method anyway and warn if something suspecious. @@ -599,10 +637,8 @@ Double_t* HistogramType::FitPerform(TH1FO* histogrampointer, TString fitFuncType fit_result_ptr = histogrampointer->Fit(fitFunc, "N,Q,W,S", "", fitMax1*0.9, posMaxValHist); // histogrampointer->Fit(fitFunc, "N,Q,W", "", fitMax1, histogrampointer->GetBinCenter(bini)); fitMax3 = fitFunc->GetMaximumX(); - fitFunc->SetLineColor(kGreen); + fitFunc->SetLineColor(kAzure); fitFunc->SetLineStyle(1); // normal for the following fits - if (verbose) - fitFunc->DrawCopy("same"); //Sort the three fits and save error estimation minFitMax = TMath::Min(TMath::Min(fitMax1,fitMax2),fitMax3); @@ -868,83 +904,87 @@ Bool_t HistogramType::calculteCCE(Bool_t verbose) { return 1; } -Bool_t HistogramType::FindNoisethresholdborder(TH1FO* histogrampointer, Bool_t overwrite, Bool_t verbose) { +Float_t HistogramType::FindNoisethresholdborder(TH1FO* histogrampointer, Bool_t overwrite, Bool_t verbose) { + Bool_t savetoclass = kFALSE; if (overwrite || noisethresholdborder == -1) { + savetoclass = kTRUE; + } // cout << colorcyan << "Find border" << endl; - Float_t posMaxValHist = histogrampointer->GetXaxis()->GetXmax(); - - TH1FO* smoothedcurce = (TH1FO*)histogrampointer->Clone(); - // smoothedcurce->Smooth(4); // not working method, rescaling is easier and more reliable - Int_t rebinningfactor = 2; - smoothedcurce->RebinX(rebinningfactor); - smoothedcurce->SetAxisRange(0,histogrampointer->GetBinCenter(histogrampointer->GetNbinsX())); - - if (verbose) { - cout << colorwhite << histogrampointer->GetName() << endlr; - TString canvastitle = Form("%s Noisethreshold border", histogrampointer->GetName()); - TString canvasname = Form("%s Noisethreshold border", histogrampointer->GetName()); - TCanvas* canvas = new TCanvas(canvasname, canvastitle, 900, 700); - smoothedcurce->SetLineColor(color+1); - smoothedcurce->Draw(); - histogrampointer->Draw("SAME"); - canvas->Update(); + Float_t posMaxValHist = histogrampointer->GetXaxis()->GetXmax(); + + TH1FO* smoothedcurce = (TH1FO*)histogrampointer->Clone(); + // smoothedcurce->Smooth(4); // not working method, rescaling is easier and more reliable + Int_t rebinningfactor = 2; + smoothedcurce->RebinX(rebinningfactor); + smoothedcurce->SetAxisRange(0,histogrampointer->GetBinCenter(histogrampointer->GetNbinsX())); + + if (verbose) { + cout << colorwhite << histogrampointer->GetName() << endlr; + TString canvastitle = Form("%s Noisethreshold border", histogrampointer->GetName()); + TString canvasname = Form("%s Noisethreshold border", histogrampointer->GetName()); + TCanvas* canvas = new TCanvas(canvasname, canvastitle, 900, 700); + smoothedcurce->SetLineColor(color+1); + smoothedcurce->Draw(); + histogrampointer->Draw("SAME"); + canvas->Update(); + } + + // Int_t bini =smoothedcurce->GetMaximumBin(); + // thresholdbincurcandidate = bini; + // if (verbose) + // cout << "GetMaximumBin: smoothedcurce->GetBinContent(" << bini << "): " << smoothedcurce->GetBinContent(bini) << endl; + Int_t bini = 0; + Int_t falling = 0; + Int_t noisepeakcandidate = 0; + + do { + Float_t curval=smoothedcurce->GetBinContent(bini++); + if (curval*1.05 <= smoothedcurce->GetBinContent(bini)) + { + falling = 0; + noisepeakcandidate = bini; } - - // Int_t bini =smoothedcurce->GetMaximumBin(); - // thresholdbincurcandidate = bini; - // if (verbose) - // cout << "GetMaximumBin: smoothedcurce->GetBinContent(" << bini << "): " << smoothedcurce->GetBinContent(bini) << endl; - Int_t bini = 0; - Int_t falling = 0; - Int_t noisepeakcandidate = 0; - - do { - Float_t curval=smoothedcurce->GetBinContent(bini++); - if (curval*1.05 <= smoothedcurce->GetBinContent(bini)) - { - falling = 0; - noisepeakcandidate = bini; - } - else - { - falling++; - if (verbose) - cout << "falling at " << smoothedcurce->GetXaxis()->GetBinCenter(bini) << " as " << curval << " < " << smoothedcurce->GetBinContent(bini) << endl; // debug - } - } while (falling < 2 && biniGetNbinsX()); - if (verbose) { - cout << "Noise peak at: smoothedcurce->GetBinContent(" << noisepeakcandidate << "): " << smoothedcurce->GetBinContent(noisepeakcandidate) << endl; + else + { + falling++; + if (verbose) + cout << "falling at " << smoothedcurce->GetXaxis()->GetBinCenter(bini) << " as " << curval << " < " << smoothedcurce->GetBinContent(bini) << endl; // debug } - Int_t thresholdbincurcandidate = 0; - bini = noisepeakcandidate; - - Int_t rising = 0; - do { - Float_t curval=smoothedcurce->GetBinContent(bini++); - if (curval*0.95 <= smoothedcurce->GetBinContent(bini)) - { - rising++; - if (verbose) - cout << "rising at bin " << bini << ":" << smoothedcurce->GetXaxis()->GetBinCenter(bini) << " as " << curval << " < " << smoothedcurce->GetBinContent(bini) << endl; // debug - } - else - { - rising = 0; - thresholdbincurcandidate = bini; - if (verbose) - cout << "falling at bin " << bini << ":" << smoothedcurce->GetXaxis()->GetBinCenter(bini) << " as " << curval << " < " << smoothedcurce->GetBinContent(bini) << endl; // debug - } - } while (rising < 3 && biniGetNbinsX()); - thresholdbincurcandidate *= rebinningfactor; - - noisethresholdborder = histogrampointer->GetXaxis()->GetBinUpEdge(thresholdbincurcandidate); - if (verbose) { - cout << " Noise threshold at " << noisethresholdborder << endl; + } while (falling < 2 && biniGetNbinsX()); + if (verbose) { + cout << "Noise peak at: smoothedcurce->GetBinContent(" << noisepeakcandidate << "): " << smoothedcurce->GetBinContent(noisepeakcandidate) << endl; + } + Int_t thresholdbincurcandidate = 0; + bini = noisepeakcandidate; + + Int_t rising = 0; + do { + Float_t curval=smoothedcurce->GetBinContent(bini++); + if (curval*0.95 <= smoothedcurce->GetBinContent(bini)) + { + rising++; + if (verbose) + cout << "rising at bin " << bini << ":" << smoothedcurce->GetXaxis()->GetBinCenter(bini) << " as " << curval << " < " << smoothedcurce->GetBinContent(bini) << endl; // debug } - return 0; + else + { + rising = 0; + thresholdbincurcandidate = bini; + if (verbose) + cout << "falling at bin " << bini << ":" << smoothedcurce->GetXaxis()->GetBinCenter(bini) << " as " << curval << " < " << smoothedcurce->GetBinContent(bini) << endl; // debug + } + } while (rising < 3 && biniGetNbinsX()); + thresholdbincurcandidate *= rebinningfactor; + + Float_t temp_noisethresholdborder = histogrampointer->GetXaxis()->GetBinUpEdge(thresholdbincurcandidate); + if (savetoclass) + noisethresholdborder = temp_noisethresholdborder; + + if (verbose) { + cout << " Noise threshold at " << temp_noisethresholdborder << endl; } - return 1; + return temp_noisethresholdborder; } Double_t HistogramType::FindBorderToPeak(TH1FO* histogrampointer, Double_t startVal, Double_t endVal, Bool_t verbose) { @@ -1004,6 +1044,7 @@ Double_t HistogramType::FindBorderToPeak(TH1FO* histogrampointer, Double_t start Bool_t HistogramType::integrateSr90Spectra(TH1FO* histogrampointer, Int_t frames_found, Float_t thresholdborder, Bool_t verbose) { Int_t thresholdbincurcandidate = 0; + if (thresholdborder < 0) { if (noisethresholdborder < 0) { @@ -1039,7 +1080,9 @@ Bool_t HistogramType::integrateSr90Spectra(TH1FO* histogrampointer, Int_t frames } } } - sr90IntegralVal = histogrampointer->Integral(thresholdbincurcandidate,histogrampointer->GetNbinsX(), "width"); + + sr90IntegralVal = histogrampointer->IntegralAndError(thresholdbincurcandidate,histogrampointer->GetNbinsX()+1, sr90IntegralErr, "width"); + // histogrampointer->Integral(min+fittries*min/20, histogrampointer->GetXaxis()->FindBin(posMaxValHist), "width"); // for ( // cout << "Integrate from bin " << thresholdbincurcandidate << " to " << histogrampointer->GetNbinsX() << endl; @@ -1050,15 +1093,46 @@ Bool_t HistogramType::integrateSr90Spectra(TH1FO* histogrampointer, Int_t frames cout << " "; if (frames_found>0) { + cout << "Integral is " << colorgreen << Form("%e",sr90IntegralVal) << "+-" << Form("%e",sr90IntegralErr) << colorreset << endl; sr90IntegralVal/=frames_found; - cout << "Integral is " << colorgreen << Form("%e",sr90IntegralVal) << colorreset << endl; + sr90IntegralErr/=frames_found; if (verbose) cout << "Scaled "; } if (verbose) { // cout << "Integral is " << colorgreen << Form("%e",sr90IntegralVal) << colorreset << ", error " << sr90IntegralErr << "%" << endl; - cout << "Integral is " << colorgreen << Form("%e",sr90IntegralVal) << endlr; + cout << "Integral is " << colorgreen << Form("%e",sr90IntegralVal)<< "+-" << Form("%e",sr90IntegralErr) << endlr; + } + + return 0; +} + + +Bool_t HistogramType::integrateSr90Spectras(Int_t frames_found, Float_t thresholdborder, Bool_t verbose) { + + Int_t thresholdbincurcandidate = 0; + Double_t cur_sr90IntegralVal = 0; + Double_t cur_sr90IntegralErr = 0; + + if (thresholdborder < 0) + { + if (noisethresholdborder < 0) { + FindNoisethresholdborder(Seed, true, verbose); + } + thresholdbincurcandidate = Seed->GetXaxis()->FindBin(noisethresholdborder); } + else + { + thresholdbincurcandidate = Seed->GetXaxis()->FindBin(thresholdborder); + noisethresholdborder = thresholdborder; + } + + for(UInt_t sumi=0;sumiIntegralAndError(thresholdbincurcandidate,a_Sum[sumi]->GetNbinsX()+1, cur_sr90IntegralErr, "width"); + a_sr90IntegralVal.push_back(cur_sr90IntegralVal); + a_sr90IntegralErr.push_back(cur_sr90IntegralErr); + } return 0; } diff --git a/MABS_run_analyzer/HistogramType.h b/MABS_run_analyzer/HistogramType.h index 51ecb57..e56e10c 100644 --- a/MABS_run_analyzer/HistogramType.h +++ b/MABS_run_analyzer/HistogramType.h @@ -71,10 +71,6 @@ private: // GENERAL HISTOGRAM PROPERTIES //***************** - Int_t color; - Int_t style; - - //***************** // GENERAL SYSTEM PROPERTIES //***************** @@ -133,6 +129,12 @@ public: /// Pointer to the run class wich initialized this histogram type Run* run = 0; + //***************** + // GENERAL HISTOGRAM PROPERTIES + //***************** + + Int_t color; + Int_t style; //***************** // TH HISTOGRAMS @@ -142,6 +144,8 @@ public: TH1FO* Seed = 0; /// Sum spectrum, the charge of whole cluster is summed up in binned into this TH1F histogram TH1FO* Sum = 0; + /// Sum spectrum, the charge of whole cluster is summed up in binned into this TH1F histogram +// TH1FO* Sum9 = 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 @@ -178,14 +182,20 @@ public: Double_t integralSeedErr=0; /// fitted position of the most probable value in the over cluster summed spectrum Float_t posSum=0; + vector a_posSum; /// fintegral over the sum peak, normalized to number of events Double_t integralSum=0; + vector a_integralSum; /// fintegral error over the Seed peak, poisson dist. assumed -> sqrt(N) for each bin Double_t integralSumErr=0; + vector a_integralSumErr; /// fitted position of the calibration peak of Fe55-beta-photons in the seed spectrum, from a run best suited to the current Float_t posVeto=0; /// fintegral over the veto peak, normalized to number of events Double_t integralVeto=0; + + /// what is the propability that a second pixel is over the threshold + Double_t clustermultiplicity=0; /// Average Noise value Float_t avgNoise = 0; @@ -204,6 +214,8 @@ public: /// Integral value, after integrating from #noisethresholdborder to maxbin. Double_t sr90IntegralVal = -1; Double_t sr90IntegralErr = -1; + vector a_sr90IntegralVal; + vector a_sr90IntegralErr; /// Peak and sigma of Seed to Sum ratio Float_t posSeedPerc = 0; @@ -312,6 +324,7 @@ public: * */ Bool_t integrateSr90Spectra(TH1FO* histogrampointer, Int_t frames_found=-1, Float_t thresholdborder=-1, Bool_t verbose=true); + Bool_t integrateSr90Spectras(Int_t frames_found, Float_t thresholdborder, Bool_t verbose); /** * @brief Tries to find a border between the noise and the signal @@ -321,7 +334,7 @@ public: * @return true if succesfull * */ - Bool_t FindNoisethresholdborder(TH1FO* histogrampointer, Bool_t overwrite=false, Bool_t verbose=false); + Float_t FindNoisethresholdborder(TH1FO* histogrampointer, Bool_t overwrite=false, Bool_t verbose=false); /** diff --git a/MABS_run_analyzer/Run.c b/MABS_run_analyzer/Run.c index 98058c1..26c5074 100644 --- a/MABS_run_analyzer/Run.c +++ b/MABS_run_analyzer/Run.c @@ -42,7 +42,7 @@ Run::Run(Int_t runnumber, Int_t loopi, TString customtitle) //db = TSQLServer::Connect("mysql://jspc29.x-matter.uni-frankfurt.de","radhard","mimosa88"); try { - string selectquery = prepareSQLStatement("select `System`, `TempCooling`, COALESCE(`TempChipStart`,-10000) as `TempChipStart`, COALESCE(`TempChipEnd`,-10000) as `TempChipEnd`, `ChipNum`, `RadiationSource`, `Matrix`, `Clock`, `StorePath`, `ChipGen`,COALESCE(`VetoPeak`,-1) as `VetoPeak`,COALESCE(`SeedPeak`,-1) as `SeedPeak`,COALESCE(`SumPeak`,-1) as `SumPeak`,COALESCE(`Gain`,-1) as `Gain`,COALESCE(`Avg.Noise`,-1) as `Avg.Noise`,COALESCE(`Avg.Noise+`,-1) as `Avg.Noise+`,COALESCE(`Avg.Noise-`,-1) as `Avg.Noise-`,COALESCE(`CCE_1`,-1) as `CCE_1`,COALESCE(`CCE_25`,-1) as `CCE_25`,COALESCE(`Frames_found`,-1) as `Frames_found`,COALESCE(`Sr90IntegralVal`,-1) as `Sr90IntegralVal`,COALESCE(`StoN`,-1) as `StoN`, `Comment`, `DepletionVoltage` from `radhard`.`labbook` WHERE `runnumber`=" + numberToString<>(labbook.runnumber)); + string selectquery = prepareSQLStatement("select `System`, `TempCooling`, COALESCE(`TempChipStart`,-10000) as `TempChipStart`, COALESCE(`TempChipEnd`,-10000) as `TempChipEnd`, `ChipNum`, `RadiationSource`, `Matrix`, `Clock`, `StorePath`, `ChipGen`,COALESCE(`VetoPeak`,-1) as `VetoPeak`,COALESCE(`SeedPeak`,-1) as `SeedPeak`,COALESCE(`SumPeak`,-1) as `SumPeak`,COALESCE(`Gain`,-1) as `Gain`,COALESCE(`Avg.Noise`,-1) as `Avg.Noise`,COALESCE(`Avg.NoiseADC`,-1) as `Avg.NoiseADC`,COALESCE(`Avg.Noise+`,-1) as `Avg.Noise+`,COALESCE(`Avg.Noise-`,-1) as `Avg.Noise-`,COALESCE(`CCE_1`,-1) as `CCE_1`,COALESCE(`CCE_25`,-1) as `CCE_25`,COALESCE(`Frames_found`,-1) as `Frames_found`,COALESCE(`Sr90IntegralVal`,-1) as `Sr90IntegralVal`,COALESCE(`StoN`,-1) as `StoN`, `Comment`, `DepletionVoltage` from `radhard`.`labbook` WHERE `runnumber`=" + numberToString<>(labbook.runnumber)); res = db->Query(selectquery.c_str()); nrows = res->GetRowCount(); if (nrows > 0) @@ -98,15 +98,16 @@ Run::Run(Int_t runnumber, Int_t loopi, TString customtitle) labbook.posSumDB = (rowsql->GetField(12) != NULL)?atof(rowsql->GetField(12)):-1; labbook.gainDB = (rowsql->GetField(13) != NULL)?atof(rowsql->GetField(13)):-1; labbook.NoiseAvgDB = (rowsql->GetField(14) != NULL)?atof(rowsql->GetField(14)):-1; - labbook.NoiseAvgPlusDB = (rowsql->GetField(15) != NULL)?atof(rowsql->GetField(15)):-1; - labbook.NoiseAvgMinusDB = (rowsql->GetField(16) != NULL)?atof(rowsql->GetField(16)):-1; - labbook.CCE_in_Perc_1DB = (rowsql->GetField(17) != NULL)?atof(rowsql->GetField(17)):-1; - labbook.CCE_in_Perc_25DB = (rowsql->GetField(18) != NULL)?atof(rowsql->GetField(18)):-1; - labbook.frames_foundDB = (rowsql->GetField(19) != NULL)?atoi(rowsql->GetField(19)):-1; - labbook.Sr90IntegralVal = (rowsql->GetField(20) != NULL)?atof(rowsql->GetField(20)):-1; - labbook.StoN = (rowsql->GetField(21) != NULL)?atof(rowsql->GetField(21)):-1; - labbook.comment = (rowsql->GetField(22) != NULL)?std::string(rowsql->GetField(22)):""; - labbook.depletionV = (rowsql->GetField(23) != NULL)?atof(rowsql->GetField(23)):-1; + labbook.NoiseAvgADC_DB = (rowsql->GetField(15) != NULL)?atof(rowsql->GetField(15)):-1; + labbook.NoiseAvgPlusDB = (rowsql->GetField(16) != NULL)?atof(rowsql->GetField(16)):-1; + labbook.NoiseAvgMinusDB = (rowsql->GetField(17) != NULL)?atof(rowsql->GetField(17)):-1; + labbook.CCE_in_Perc_1DB = (rowsql->GetField(18) != NULL)?atof(rowsql->GetField(18)):-1; + labbook.CCE_in_Perc_25DB = (rowsql->GetField(19) != NULL)?atof(rowsql->GetField(19)):-1; + labbook.frames_foundDB = (rowsql->GetField(20) != NULL)?atoi(rowsql->GetField(20)):-1; + labbook.Sr90IntegralVal = (rowsql->GetField(21) != NULL)?atof(rowsql->GetField(21)):-1; + labbook.StoN = (rowsql->GetField(22) != NULL)?atof(rowsql->GetField(22)):-1; + labbook.comment = (rowsql->GetField(23) != NULL)?std::string(rowsql->GetField(23)):""; + labbook.depletionV = (rowsql->GetField(24) != NULL)?atof(rowsql->GetField(24)):-1; // labbook.frames_Analyzed = (rowsql->GetField(20) != NULL)?atoi(rowsql->GetField(20)):-1; delete res; if (labbook.chipGen.Length() > 0) @@ -233,6 +234,7 @@ Bool_t Run::debugDBreadout() cout << "| posSumDB: " << std::right << colorwhite << labbook.posSumDB << endlr; cout << "| posVetoDB: " << std::right << colorwhite << labbook.posVetoDB << endlr; cout << "| gainDB: " << std::right << colorwhite << labbook.gainDB << endlr; + cout << "| NoiseAvgADC_DB " << std::right << colorwhite << labbook.NoiseAvgADC_DB << endlr; cout << "| NoiseAvgDB: " << std::right << colorwhite << labbook.NoiseAvgDB << endlr; cout << "| NoiseAvgPlusDB: " << std::right << colorwhite << labbook.NoiseAvgPlusDB << endlr; cout << "| NoiseAvgMinusDB: " << std::right << colorwhite << labbook.NoiseAvgMinusDB << endlr; @@ -256,6 +258,7 @@ void Run::setSystemSpecificParameters() systemparam systemparamFSBB (2800,2800/4,25,10,100); systemparam systemparamUSBMi19 (400/*maxbin*/,400/1/*nbins*/, 25/*vetothreshold*/, 10/*maxbinnoise*/, 100/*nbinsnoise*/); systemparam systemparamPipper2 (1000,500,25,10,100); // change bin size to 200 if Cd109 + systemparam systemparamPipper2_Sr (3000,1200,25,10,100); // change bin size to 200 if Cd109 systemparam systemparamPipper2_Cd (1000,200,25,10,100); // change bin size to 200 if Cd109 if (labbook.system.EqualTo("USB") && labbook.chipGen.EqualTo("Mi34") ) cursystemparam = systemparamUSB; @@ -270,6 +273,8 @@ void Run::setSystemSpecificParameters() else if (labbook.system.EqualTo("USB") && labbook.chipGen.EqualTo("Pipper2") ) { if (labbook.source.Contains("Cd")) cursystemparam = systemparamPipper2_Cd; + else if (labbook.source.Contains("Sr")) + cursystemparam = systemparamPipper2_Sr; else cursystemparam = systemparamPipper2; } @@ -493,6 +498,7 @@ Bool_t Run::analyzeRun(Bool_t force) if (labbook.source.Contains("Sr90")) { //cout << colorwhite << " integrateSr90Spectra():" << endlr; (*curHistogramClass)->integrateSr90Spectra((*curHistogramClass)->Seed, frames_found, -1, false); + (*curHistogramClass)->integrateSr90Spectras(frames_found, 58, false); cout << colorwhite << " calculteStoN():" << endlr; (*curHistogramClass)->calculteStoN(true); } @@ -501,8 +507,11 @@ Bool_t Run::analyzeRun(Bool_t force) if (labbook.source.Contains("Sr90")) { cout << colorwhite << "integrateSr90Spectra():" << endlr; histogramthreshold->integrateSr90Spectra(histogramthreshold->Seed, frames_found, -1, true); - if (histogramfixedthreshold != 0) + histogramthreshold->integrateSr90Spectras(frames_found, -1, false); + if (histogramfixedthreshold != 0) { histogramfixedthreshold->integrateSr90Spectra(histogramfixedthreshold->Seed, frames_found, 0); + histogramfixedthreshold->integrateSr90Spectras(frames_found, -1, false); + } } cout << colorwhite << "normalizeHistogramClasses():" << endlr; normalizeHistogramClasses(); @@ -1280,6 +1289,7 @@ Bool_t Run::binSeedSumVeto() Double_t gausssigma = 0; + for (Int_t filli=0; filli<3; filli++) { // first and second fill loop for (Int_t framei=0; frameifHitTree->GetEntries(); framei++) // loop over all frames { @@ -1358,6 +1368,9 @@ Bool_t Run::binSeedSumVeto() histogram->a_Sum[clusterj]->Fill(a_pixelSum[clusterj]); } + if (a_pixelSum[1]-a_pixelSum[0] > labbook.NoiseAvgADC_DB*5) + histogram->clustermultiplicity++; + // seed percentage spectrum histogram->SeedPerc->Fill(processed->fFrameInfo.p[12][hiti]/pixelSum*100); @@ -1530,6 +1543,17 @@ Bool_t Run::binSeedSumVeto() (*curHistogramClass)->integralSum = parameters[6]; (*curHistogramClass)->integralSumErr = parameters[9]; + (*curHistogramClass)->clustermultiplicity /=(*curHistogramClass)->numberofhits; + + for (u_int sumi = 0; sumi < (*curHistogramClass)->a_Sum.size(); sumi++) { + // Float_t integratestart = (*curHistogramClass)->FindNoisethresholdborder((*curHistogramClass)->a_Sum[sumi], false, false); +// parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->a_Sum[sumi], "gaus", 0, false, integratestart); //TODO change back to gauss + parameters = (*curHistogramClass)->FitPerform((*curHistogramClass)->a_Sum[sumi], "gaus", 0, false, (*curHistogramClass)->fixedThresholdValue); //TODO change back to gauss + (*curHistogramClass)->a_posSum.push_back(parameters[1]); + (*curHistogramClass)->a_integralSum.push_back(parameters[6]); + (*curHistogramClass)->a_integralSumErr.push_back(parameters[9]); + } + parameters =(*curHistogramClass)->FitPerform((*curHistogramClass)->SeedPerc, "landau"); (*curHistogramClass)->posSeedPerc = parameters[1]; (*curHistogramClass)->sigmaSeedPerc = parameters[2]; @@ -1555,6 +1579,9 @@ void Run::fillAHistogramsinclass(HistogramType* histogramtypepointer, Int_t hiti 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 + + if (pt_a_pixelSum[1]-pt_a_pixelSum[0] > labbook.NoiseAvgADC_DB*5) + histogramtypepointer->clustermultiplicity++; } diff --git a/MABS_run_analyzer/help.h b/MABS_run_analyzer/help.h index 880334e..4fc72e7 100644 --- a/MABS_run_analyzer/help.h +++ b/MABS_run_analyzer/help.h @@ -16,6 +16,7 @@ #include #include +#include #include #include #include @@ -25,6 +26,9 @@ #define MAXHITS 1000000 #define MAXPIXELS 100000 + + + //#################################################################### /** * @file help.h @@ -235,7 +239,9 @@ struct labbooksctruct Float_t posVetoDB = -1; /// gain found in db Float_t gainDB = -1; - /// average noise found in db + /// average noise found in db in ADU + Float_t NoiseAvgADC_DB = -1; + /// average noise found in db in e Float_t NoiseAvgDB = -1; /// Positive noise error found in db Float_t NoiseAvgPlusDB = -1; @@ -509,6 +515,13 @@ Float_t stringToNumber ( string Text ) return ss >> result ? result : 0; } +class HistogramType; +class TNtupleO: public TNtuple { +public: + HistogramType* itsHistogramType = 0; + TString xaxisTitle = "", yaxisTitle = "", Title = ""; + using TNtuple::TNtuple; +}; /** -- 2.43.0