From 67ebc5835182140a7b64e4a7d7fa19e9102ec278 Mon Sep 17 00:00:00 2001 From: Benjamin Linnik Date: Thu, 20 Aug 2015 17:11:03 +0200 Subject: [PATCH] Run analyzer: Added average cluster calculation routine --- MABS_run_analyzer/ChargeSpektrum.c | 3 +- MABS_run_analyzer/ChargeSpektrumFunctions.c | 3 +- MABS_run_analyzer/Run.c | 468 ++++++++++++++++++-- MABS_run_analyzer/Run.h | 60 ++- 4 files changed, 495 insertions(+), 39 deletions(-) diff --git a/MABS_run_analyzer/ChargeSpektrum.c b/MABS_run_analyzer/ChargeSpektrum.c index 6f94117..2375ddb 100644 --- a/MABS_run_analyzer/ChargeSpektrum.c +++ b/MABS_run_analyzer/ChargeSpektrum.c @@ -160,7 +160,8 @@ void ChargeSpektrum(TString runnumber = "") if (!isBatch) gROOT->SetBatch(kFALSE); // runs[runi]->plotAllHistograms(); - runs[runi]->plotAllHistogramsThresholdCluster(); +// runs[runi]->plotAllHistogramsThresholdCluster(); + runs[runi]->plotClusterDistribution(&runs[runi]->histogramCalibrated); runs[runi]->plotAllHistogramsThresholdClusterCalibrated(); // runs[runi]->plotAllHistogramsCalibrated(); // runs[runi]->writeAllHistogramsToFile(); diff --git a/MABS_run_analyzer/ChargeSpektrumFunctions.c b/MABS_run_analyzer/ChargeSpektrumFunctions.c index 52fe9a3..030273f 100644 --- a/MABS_run_analyzer/ChargeSpektrumFunctions.c +++ b/MABS_run_analyzer/ChargeSpektrumFunctions.c @@ -121,7 +121,7 @@ Bool_t writeObservableToFile() headerInfo+=runs[runi]->runcode+"\t\t\t"; } TH1F* plothistogrampointer = *runs[0]->plothistogrampointer; - TString runnumberListe=Form("%s_",plothistogrampointer->GetName()); + TString runnumberListe=""; TString header=""; for(Int_t runi=0;runiGetName()); TString filename= runs[0]->savepathresults + "/" + runnumberListe + "histograms.dat"; fstream* fout = new fstream(filename,ios::out); *fout << headerInfo << endl; diff --git a/MABS_run_analyzer/Run.c b/MABS_run_analyzer/Run.c index 103c449..ecb94f0 100644 --- a/MABS_run_analyzer/Run.c +++ b/MABS_run_analyzer/Run.c @@ -99,19 +99,37 @@ Run::Run(Int_t runnumber, Int_t loopi) labbook.CCE_in_Perc_25DB = (rowsql->GetField(18) != NULL)?atoi(rowsql->GetField(18)):-1; labbook.frames_foundDB = (rowsql->GetField(19) != NULL)?atoi(rowsql->GetField(19)):-1; delete res; - if (labbook.chip.Length() > 0 && labbook.chipGen.Length() > 0) // versuche infos zum Chip aus der ChipDatenbank zu bekommen + if (labbook.chipGen.Length() > 0) { - selectquery=prepareSQLStatement("select `epi_thickness`, `resistivity`, `ChipRadiation Ion`, `ChipRadiation NonIon` from `radhard`.`chips` WHERE `no`='" + numberToString<>(labbook.chip) + "' AND `chipgen`='" + labbook.chipGen + "'"); - res = db->Query(selectquery.c_str()); - nrows = res->GetRowCount(); - if (nrows > 0) + if (labbook.chip.Length() > 0) // versuche infos zum Chip aus der ChipDatenbank zu bekommen { - rowsql = res->Next(); - labbook.epi_thickness = (rowsql->GetField(0) != NULL)?atol(rowsql->GetField(0)):-1; - labbook.resistivity = (rowsql->GetField(1) != NULL)?atof(rowsql->GetField(1)):-1; - labbook.radDoseIon = (rowsql->GetField(2) != NULL)?atof(rowsql->GetField(2)):-1; - labbook.radDoseNonIon = (rowsql->GetField(3) != NULL)?atof(rowsql->GetField(3)):-1; - delete res; + selectquery=prepareSQLStatement("select `epi_thickness`, `resistivity`, `ChipRadiation Ion`, `ChipRadiation NonIon` from `radhard`.`chips` WHERE `no`='" + numberToString<>(labbook.chip) + "' AND `chipgen`='" + labbook.chipGen + "'"); + res = db->Query(selectquery.c_str()); + nrows = res->GetRowCount(); + if (nrows > 0) + { + rowsql = res->Next(); + labbook.epi_thickness = (rowsql->GetField(0) != NULL)?atol(rowsql->GetField(0)):-1; + labbook.resistivity = (rowsql->GetField(1) != NULL)?atof(rowsql->GetField(1)):-1; + labbook.radDoseIon = (rowsql->GetField(2) != NULL)?atof(rowsql->GetField(2)):-1; + labbook.radDoseNonIon = (rowsql->GetField(3) != NULL)?atof(rowsql->GetField(3)):-1; + delete res; + } + } + if (labbook.matrix.Length() > 0) // versuche infos zum Pixel aus der ChipDatenbank zu bekommen + { + selectquery=prepareSQLStatement("select `pitchX`, `pitchY`, `num_diod`, `staggered` from `radhard`.`pixelinfo` WHERE `matrix`='" + labbook.matrix + "' AND `ChipGen`='" + labbook.chipGen + "'"); + res = db->Query(selectquery.c_str()); + nrows = res->GetRowCount(); + if (nrows > 0) + { + rowsql = res->Next(); + curpixelinfo.pitchX = (rowsql->GetField(0) != NULL)?atof(rowsql->GetField(0)):-1; + curpixelinfo.pitchY = (rowsql->GetField(1) != NULL)?atof(rowsql->GetField(1)):-1; + curpixelinfo.ndiods = (rowsql->GetField(2) != NULL)?atoi(rowsql->GetField(2)):-1; + curpixelinfo.staggered = (rowsql->GetField(3) != NULL)?atoi(rowsql->GetField(3)):0; + delete res; + } } } // if (!(labbook.posVetoDB > 0) && (labbook.source != "Fe55")) // no veto peak position found for this run @@ -155,6 +173,10 @@ Bool_t Run::debugDBreadout() cout << "| chipGen: " << std::right << colorwhite << labbook.chipGen << endlr; cout << "| source: " << std::right << colorwhite << labbook.source << endlr; cout << "| matrix: " << std::right << colorwhite << labbook.matrix << endlr; + cout << "| pitch X: " << std::right << colorwhite << curpixelinfo.pitchX << endlr; + cout << "| pitch Y: " << std::right << colorwhite << curpixelinfo.pitchY << endlr; + cout << "| staggered: " << std::right << colorwhite << (curpixelinfo.staggered?"Yes":"No") << endlr; + cout << "| matrix: " << std::right << colorwhite << labbook.matrix << endlr; cout << "| radDoseIon: " << std::right << colorwhite << labbook.radDoseIon << endlr; cout << "| radDoseNonIon: " << std::right << colorwhite << labbook.radDoseNonIon << endlr; cout << "| resistivity: " << std::right << colorwhite << labbook.resistivity << endlr; @@ -277,7 +299,9 @@ Bool_t Run::analyzeRun(Bool_t force) cout << colorwhite << "binNoise():" << endlr; binNoise(); cout << colorwhite << "binSeedSumVeto():" << endlr; - binSeedSumVeto(); + binSeedSumVeto(); + cout << colorwhite << "binCluster():" << endlr; + binCluster(); cout << colorwhite << "rescaleHistograms():" << endlr; histogramCalibrated.calibrated = rescaleHistograms(); histogramthresholdCalibrated.calibrated = histogramCalibrated.calibrated; @@ -341,6 +365,11 @@ Bool_t Run::rescaleHistograms() rescaleHistogram(histogramCalibrated.Veto, histogram.Veto); rescaleHistogram(histogramCalibrated.Noise, histogram.Noise); + if (histogram.histAvgCluster != nullptr) + { + rescale2DHistogramCounts(histogramCalibrated.histAvgCluster, histogram.histAvgCluster); + } + histogramCalibrated.posSeed = histogram.posSeed * gain; histogramCalibrated.posSum = histogram.posSum * gain; histogramCalibrated.posVeto = histogram.posVeto * gain; @@ -360,6 +389,11 @@ Bool_t Run::rescaleHistograms() histogramthresholdCalibrated.avgNoisePlus = histogramthreshold.avgNoisePlus * gain; histogramthresholdCalibrated.avgNoiseMinus = histogramthreshold.avgNoiseMinus * gain; + if (histogramthreshold.histAvgCluster != nullptr) + { + rescale2DHistogramCounts(histogramthresholdCalibrated.histAvgCluster, histogramthreshold.histAvgCluster); + } + return 1; } @@ -379,6 +413,23 @@ void Run::rescaleHistogram(TH1F* &histogrampointernew, TH1F* &histogrampointerol } +void Run::rescale2DHistogramCounts(TH2F* &histogrampointernew, TH2F* &histogrampointerold ) +{ + histogrampointernew = (TH2F*)histogrampointerold->Clone(); + histogrampointernew->SetName(Form("%sC", histogrampointerold->GetName())); + histogrampointernew->SetTitle(Form("%sC", histogrampointerold->GetTitle())); + histogrampointernew->GetZaxis()->SetTitle("Q_coll [e]"); + int nbinsx = histogrampointernew->GetXaxis()->GetNbins(); + int nbinsy = histogrampointernew->GetYaxis()->GetNbins(); + + double new_bins[nbinsx][nbinsy]; + for(int y=0; y < nbinsy; y++){ + for(int x=0; x < nbinsx; x++){ + histogrampointernew->SetBinContent(x,y,histogrampointerold->GetBinContent(x,y)*gain); + } + } +} + Bool_t Run::analyzeFrame(Int_t frame) { if (!error) @@ -671,7 +722,7 @@ Bool_t Run::binSeedSumVeto() { // /// pixel number of seed pixel, position on sensor // UInt_t seedPixel[10000]; -// /// Array of CLUSTERSIZE * CLUSTERSIZE clusters, seed pixel in the middle +// /// Array of processed->clustersize * processed->clustersize clusters, seed pixel in the middle // Float_t pixelcluster[processed->clustersize*processed->clustersize][10000]; /// collected charge in cluster @@ -686,8 +737,10 @@ Bool_t Run::binSeedSumVeto() { for(Int_t hiti=0; (unsigned int)hitifFrameInfo.hits;hiti++) { + histogram.numberofhits++; + // histogram with the single pixel - histogram.Seed->Fill(processed->fFrameInfo.p[12][hiti]); + histogram.Seed->Fill(processed->fFrameInfo.p[12][hiti]); // sum histogram pixelSum = 0; @@ -706,6 +759,8 @@ Bool_t Run::binSeedSumVeto() if (processed->fFrameInfo.pixelthreshold[hiti]>0) { + histogramthreshold.numberofhits++; + histogramthreshold.Seed->Fill(processed->fFrameInfo.p[12][hiti]); histogramthreshold.Sum->Fill(pixelSum); if (TMath::Abs(notSeedSum) < systemparamcur.vetothreshold && (labbook.source.Contains("Fe")||labbook.source.Contains("Cd"))) @@ -727,22 +782,227 @@ Bool_t Run::binSeedSumVeto() return 0; } -void Run::setPlotStyle(Int_t x){ - plotStyle = x%13; - histogram.Seed->SetLineColor(rootcolors[plotStyle]); - histogram.Sum->SetLineColor(rootcolors[plotStyle]); - histogram.Veto->SetLineColor(rootcolors[plotStyle]); - histogram.Noise->SetLineColor(rootcolors[plotStyle]); +Bool_t Run::binCluster() +{ + // histogram.histAvgCluster->SetBit(TH1::kCanRebin); + // histogramthreshold.histAvgCluster->SetBit(TH1::kCanRebin); + Float_t rotateangle = 0; + if (curpixelinfo.staggered) { + rotateangle = TMath::ATan(curpixelinfo.pitchY/(curpixelinfo.pitchX/2))*180.0/TMath::Pi(); + cout << " Rotate angle: " << rotateangle << endl; + findRotatedClusterDimension(rotateangle); + } + else + { + initClusters(curpixelinfo.pitchY, 0,curpixelinfo.pitchY*5,curpixelinfo.pitchX,0,curpixelinfo.pitchX*5); + // prefix = "Cluster" + suffix; + // histogramstructpointer->histAvgCluster = new TH2F(prefix.Data(), Form("Avg. cluster distribution; x; y %s", humanreadablestr.Data()),5, 0,curpixelinfo.pitchY*5,5,0,curpixelinfo.pitchX*5); // TODO: remove hardcoded 5 + } + + + Float_t xcoord, xcoord_old, ycoord; + + for (Int_t framei=0; frameifHitTree->GetEntries(); framei++) // loop over all frames + { + processed->fHitTree->GetEntry(framei); + // account only frames with less then 10 hits + if (processed->fFrameInfo.hits<(unsigned int)10) + { + for(Int_t hiti=0; (unsigned int)hitifFrameInfo.hits;hiti++) + { + if (processed->fFrameInfo.pixel[hiti]%2==0) // Diode sitzt oben im SeedPixel, da nach PitchY angeordnet + { + for (Int_t clusteri=0; clustericlustersize*processed->clustersize; clusteri++) + { + if ((clusteri%processed->clustersize)%2==0 || rotateangle == 0) + { + // betrachte für Pixel mit Diode links + xcoord = (clusteri%processed->clustersize)*curpixelinfo.pitchY+11; + ycoord = (int)(clusteri/processed->clustersize)*curpixelinfo.pitchX+11; + xcoord_old = xcoord; + xcoord = xcoord*TMath::Cos(rotateangle/180*TMath::Pi())-ycoord*TMath::Sin(rotateangle/180*TMath::Pi()); + ycoord = xcoord_old*TMath::Sin(rotateangle/180*TMath::Pi())+ycoord*TMath::Cos(rotateangle/180*TMath::Pi()); + } + else + { + // betrachte für Pixel mit Diode links + xcoord = (clusteri%processed->clustersize)*curpixelinfo.pitchY+11; + ycoord = (int)(clusteri/processed->clustersize)*curpixelinfo.pitchX+curpixelinfo.pitchX/2+11; + xcoord_old = xcoord; + xcoord = xcoord*TMath::Cos(rotateangle/180*TMath::Pi())-ycoord*TMath::Sin(rotateangle/180*TMath::Pi()); + ycoord = xcoord_old*TMath::Sin(rotateangle/180*TMath::Pi())+ycoord*TMath::Cos(rotateangle/180*TMath::Pi()); + } +// cout << colorwhite << "oben: xcoord: " << xcoord << ", ycoord: " << ycoord << endl; + + histogram.histAvgCluster->Fill(xcoord,ycoord,processed->fFrameInfo.p[clusteri][hiti]); + if (processed->fFrameInfo.pixelthreshold[hiti]>0) + histogramthreshold.histAvgCluster->Fill(xcoord,ycoord,processed->fFrameInfo.p[clusteri][hiti]); + } + } + else // Diode sitzt unten im SeedPixel, da nach PitchY angeordnet + { + for (Int_t clusteri=0; clustericlustersize*processed->clustersize; clusteri++) + { + if ((clusteri%processed->clustersize)%2==1 || rotateangle == 0) + { + // betrachte für Pixel mit Diode links + xcoord = (clusteri%processed->clustersize)*curpixelinfo.pitchY+11; + ycoord = (int)(clusteri/processed->clustersize)*curpixelinfo.pitchX+11; + xcoord_old = xcoord; + xcoord = xcoord*TMath::Cos(rotateangle/180*TMath::Pi())-ycoord*TMath::Sin(rotateangle/180*TMath::Pi()); + ycoord = xcoord_old*TMath::Sin(rotateangle/180*TMath::Pi())+ycoord*TMath::Cos(rotateangle/180*TMath::Pi()); + } + else + { + // betrachte für Pixel mit Diode links + xcoord = (clusteri%processed->clustersize)*curpixelinfo.pitchY+11; + ycoord = (int)(clusteri/processed->clustersize)*curpixelinfo.pitchX+curpixelinfo.pitchX/2+11; + xcoord_old = xcoord; + xcoord = xcoord*TMath::Cos(rotateangle/180*TMath::Pi())-ycoord*TMath::Sin(rotateangle/180*TMath::Pi()); + ycoord = xcoord_old*TMath::Sin(rotateangle/180*TMath::Pi())+ycoord*TMath::Cos(rotateangle/180*TMath::Pi()); + } +// cout << colorwhite << "unten: xcoord: " << xcoord << ", ycoord: " << ycoord << endl; + + histogram.histAvgCluster->Fill(xcoord,ycoord,processed->fFrameInfo.p[clusteri][hiti]); + if (processed->fFrameInfo.pixelthreshold[hiti]>0) + histogramthreshold.histAvgCluster->Fill(xcoord,ycoord,processed->fFrameInfo.p[clusteri][hiti]); + } + } + } + } + } + + for (Int_t clusterx = 1; clusterx <=histogram.histAvgCluster->GetXaxis()->GetNbins();clusterx++) { + for (Int_t clustery = 1; clustery <=histogram.histAvgCluster->GetYaxis()->GetNbins();clustery++) { + histogram.histAvgCluster->SetBinContent(clusterx,clustery,histogram.histAvgCluster->GetBinContent(clusterx,clustery)/histogram.numberofhits); + histogramthreshold.histAvgCluster->SetBinContent(clusterx,clustery,histogramthreshold.histAvgCluster->GetBinContent(clusterx,clustery)/histogramthreshold.numberofhits); + } + } + + return 0; +} + +Bool_t Run:: findRotatedClusterDimension(Float_t rotateangle) +{ + /// x coordinate in micrometer + Float_t xcoord = 0; + /// y coordinate in micrometer + Float_t ycoord = 0; + /// temporal x coordinate + Float_t xcoord_old; + /// after rotation, the minimal coordinate in x direction for one sort of a cluster (seed with diode in the upper part or lower part) [micrometer] + Float_t xcoord_min=999999; //(int)(0/5)*66+66/6*1*cos(rotateangle/180*M_PI)-(0%5)*33+33/3*1*sin(rotateangle/180*M_PI); + /// after rotation, the minimal coordinate in x direction for both, clusters with diode in upper part and lower part [micrometer] + Float_t xcoord_abs_min=999999; //(int)(0/5)*66+66/6*1*cos(rotateangle/180*M_PI)-(0%5)*33+33/3*1*sin(rotateangle/180*M_PI); + /// after rotation, the minimal coordinate in x direction for one sort of a cluster (seed with diode in the upper part or lower part) [micrometer] + Float_t ycoord_min=999999; + /// after rotation, the minimal coordinate in x direction for both: clusters with diode in upper part and lower part [micrometer] + Float_t ycoord_abs_min=999999; + Float_t xcoord_max=0; + /// after rotation, the maximal coordinate in x direction for both, clusters with diode in upper part and lower part [micrometer] + Float_t xcoord_abs_max=0; + /// after rotation, the maximal coordinate in y direction for one clusters with diode in upper part or lower part [micrometer] + Float_t ycoord_max=0; + /// after rotation, the global maximal coordinate in y direction for clusters with diode in upper part and lower part [micrometer] + Float_t ycoord_abs_max=0; + /// after rotation, the step size in x direction from one diode to the next + Float_t xcoord_min_step=999999; + /// after rotation, the step size in y direction from one diode to the next + Float_t ycoord_min_step=999999; + /// after rotation, the maximal number of diodes in x direction + Int_t sizearrrotx = 1; + /// after rotation, the maximal number of diodes in y direction + Int_t sizearrroty = 1; - histogramCalibrated.Seed->SetLineColor(rootcolors[plotStyle]); - histogramCalibrated.Sum->SetLineColor(rootcolors[plotStyle]); - histogramCalibrated.Veto->SetLineColor(rootcolors[plotStyle]); - histogramCalibrated.Noise->SetLineColor(rootcolors[plotStyle]); - histogramthreshold.Seed->SetLineColor(rootcolors[plotStyle]); - histogramthreshold.Sum->SetLineColor(rootcolors[plotStyle]); - histogramthreshold.Veto->SetLineColor(rootcolors[plotStyle]); - histogramthreshold.Noise->SetLineColor(rootcolors[plotStyle]); + //TODO beatify code + // find out size coord_min_step values after rotation and coord_min values + for (Int_t clusteri=0; clustericlustersize*processed->clustersize; clusteri++) + { + if ((clusteri%processed->clustersize)%2==1) + { + // betrachte für Pixel mit Diode links + xcoord = (clusteri%processed->clustersize)*curpixelinfo.pitchY+11; + ycoord = (int)(clusteri/processed->clustersize)*curpixelinfo.pitchX+11; + + /// after rotation, the minimal coordinate in x direction for both, clusters with diode in upper part and lower part [micrometer]LUSTERSIZE)*curpixelinfo.pitchX+11; + xcoord_old = xcoord; + xcoord = xcoord*cos(rotateangle/180*TMath::Pi())-ycoord*sin(rotateangle/180*TMath::Pi()); + ycoord = xcoord_old*sin(rotateangle/180*TMath::Pi())+ycoord*cos(rotateangle/180*TMath::Pi()); + xcoord_min=(xcoordxcoord_max?xcoord:xcoord_max); + xcoord_min_step=(abs(xcoord-xcoord_min)0.1?abs(xcoord-xcoord_min):xcoord_min_step); + // cout << "xcoord: " << xcoord << " xcoord_min: " << xcoord_min << " xcoord-xcoord_min: " << xcoord-xcoord_min << endl; + ycoord_min=(ycoordycoord_max?ycoord:ycoord_max); + ycoord_min_step=(abs(ycoord-ycoord_min)0.1?abs(ycoord-ycoord_min):ycoord_min_step); + // cout << "ycoord: " << ycoord << " ycoord_min: " << ycoord_min << " ycoord-ycoord_min: " << ycoord-ycoord_min << endl; + // cout << "xcoord: " << xcoord << endl; + + // betrachte für Pixel mit Diode rechts, jedoch nur für globale minima maxima + xcoord = (clusteri%processed->clustersize)*curpixelinfo.pitchY+11; + ycoord = (int)(clusteri/processed->clustersize)*curpixelinfo.pitchX+11+curpixelinfo.pitchX/2; + xcoord_old = xcoord; + xcoord = xcoord*cos(rotateangle/180*TMath::Pi())-ycoord*sin(rotateangle/180*TMath::Pi()); + ycoord = xcoord_old*sin(rotateangle/180*TMath::Pi())+ycoord*cos(rotateangle/180*TMath::Pi()); + xcoord_abs_min=(xcoord_minxcoord_abs_max?xcoord_max:xcoord_abs_max); + xcoord_abs_max=(xcoord>xcoord_abs_max?xcoord:xcoord_abs_max); + ycoord_abs_min=(ycoord_minycoord_abs_max?ycoord_max:ycoord_abs_max); + ycoord_abs_max=(ycoord>ycoord_abs_max?ycoord:ycoord_abs_max); + } + else + { + xcoord = (clusteri%processed->clustersize)*curpixelinfo.pitchY+11; + ycoord = (int)(clusteri/processed->clustersize)*curpixelinfo.pitchX+11+curpixelinfo.pitchX/2; + xcoord_old = xcoord; + xcoord = xcoord*cos(rotateangle/180*TMath::Pi())-ycoord*sin(rotateangle/180*TMath::Pi()); + ycoord = xcoord_old*sin(rotateangle/180*TMath::Pi())+ycoord*cos(rotateangle/180*TMath::Pi()); + xcoord_min=(xcoordxcoord_max?xcoord:xcoord_max); + // cout << "xcoord: " << xcoord << " xcoord_min: " << xcoord_min << " xcoord-xcoord_min: " << xcoord-xcoord_min << endl; + xcoord_min_step=(abs(xcoord-xcoord_min)0.1?abs(xcoord-xcoord_min):xcoord_min_step); + ycoord_min=(ycoordycoord_max?ycoord:ycoord_max); + ycoord_min_step=(abs(ycoord-ycoord_min)0.1?abs(ycoord-ycoord_min):ycoord_min_step); + // cout << "ycoord: " << ycoord << " ycoord_min: " << ycoord_min << " ycoord-ycoord_min: " << ycoord-ycoord_min << endl; + + // cout << "xcoord: " << xcoord << endl; + + xcoord = (clusteri%processed->clustersize)*curpixelinfo.pitchY+11; + ycoord = (int)(clusteri/processed->clustersize)*curpixelinfo.pitchX+11; + xcoord_old = xcoord; + xcoord = xcoord*cos(rotateangle/180*TMath::Pi())-ycoord*sin(rotateangle/180*TMath::Pi()); + ycoord = xcoord_old*sin(rotateangle/180*TMath::Pi())+ycoord*cos(rotateangle/180*TMath::Pi()); + xcoord_abs_min=(xcoord_minxcoord_abs_max?xcoord_max:xcoord_abs_max); + xcoord_abs_max=(xcoord>xcoord_abs_max?xcoord:xcoord_abs_max); + ycoord_abs_min=(ycoord_minycoord_abs_max?ycoord_max:ycoord_abs_max); + ycoord_abs_max=(ycoord>ycoord_abs_max?ycoord:ycoord_abs_max); + } + } + // xcoord_min_step=abs(curpixelinfo.pitchX*cos(rotateangle/180*M_PI)-curpixelinfo.pitchY*sin(rotateangle/180*M_PI)); + // ycoord_min_step=abs(curpixelinfo.pitchX*sin(rotateangle/180*M_PI)+curpixelinfo.pitchY*cos(rotateangle/180*M_PI)); +// cout << "xcoord_min: " << xcoord_min << " xcoord_max: " << xcoord_max << endl; +// cout << "ycoord_min: " << ycoord_min << " ycoord_max: " << ycoord_max << endl; +// cout << "xcoord_min_step: " << xcoord_min_step << " ycoord_min_step: " << ycoord_min_step << endl; + // double deleteme = curpixelinfo.pitchX*processed->clustersize*cos(rotateangle/180*M_PI)-curpixelinfo.pitchY*processed->clustersize*sin(rotateangle/180*M_PI); + // double deleteme2 = curpixelinfo.pitchX*processed->clustersize*sin(rotateangle/180*M_PI)+curpixelinfo.pitchY*processed->clustersize*cos(rotateangle/180*M_PI); + // cout << "xmax2: " << deleteme << " ymax2: " << deleteme2 << endl; +// sizearrrotx= (int)((xcoord_max-xcoord_min+1)/xcoord_min_step)+1; +// sizearrroty = (int)((ycoord_max-ycoord_min+1)/ycoord_min_step)+1; +// cout << "sizearrrotx: " << sizearrrotx << " sizearrroty: " << sizearrroty << endl; + + initClusters(xcoord_min_step, xcoord_abs_min, xcoord_abs_max, ycoord_min_step, ycoord_abs_min, ycoord_abs_max); + //sizearrrotx*2, xcoord_abs_min,xcoord_abs_max,sizearrroty*2,ycoord_abs_min,ycoord_abs_max); + + return 0; } Bool_t Run::plotNoise() @@ -869,6 +1129,115 @@ Bool_t Run::plotAllHistograms(histogramstruct* histogramstructpointer) return 1; } +Bool_t Run::plotClusterDistribution(histogramstruct* histogramstructpointer) +{ + Int_t random = random1->Rndm()*100000; + TString canvastitle = Form("%s_Cluster", runcode.Data()); + if (histogramstructpointer->calibrated) + canvastitle += "_el"; + if (histogramstructpointer->thresholdcluster) + canvastitle += "_trsh"; + TString canvasname = Form("%s%d",runcode.Data(),random); + TCanvas* Canvas_1 = new TCanvas(canvasname, canvastitle, 1200, 800); + TView *view = TView::CreateView(1);; + + Canvas_1->SetFillColor(0); + Canvas_1->SetBorderMode(0); + Canvas_1->SetBorderSize(2); + // Canvas_1->SetTheta(64.291); + // Canvas_1->SetPhi(32.46973); + Canvas_1->TAttPad::SetFrameBorderMode(0); + Canvas_1->TAttPad::SetFrameBorderMode(0); + histogramstructpointer->histAvgCluster->SetContour(20); + + TPaletteAxis *palette = new TPaletteAxis(0.8596098,-0.8025656,0.8949913,0.9836416,histogramstructpointer->histAvgCluster); + palette->SetLabelColor(1); + palette->SetLabelFont(42); + palette->SetLabelOffset(0.005); + palette->SetLabelSize(0.025); + palette->SetTitleOffset(1); + palette->SetTitleSize(0.035); + palette->SetFillColor(100); + palette->SetFillStyle(1001); + histogramstructpointer->histAvgCluster->GetListOfFunctions()->Add(palette,"br"); + + TPaveStats *ptstats = new TPaveStats(0.7013721,0.7794677,0.850686,0.9226869,"brNDC"); + ptstats->SetTextFont(42); + ptstats->SetTextSize(0.02534854); + ptstats->SetOptStat(1011); + ptstats->SetOptFit(0); + ptstats->Draw(); + histogramstructpointer->histAvgCluster->GetListOfFunctions()->Add(ptstats); + ptstats->SetParent(histogramstructpointer->histAvgCluster); + + Int_t ci; // for color index setting + ci = TColor::GetColor("#000099"); + histogramstructpointer->histAvgCluster->SetLineColor(ci); + histogramstructpointer->histAvgCluster->GetXaxis()->SetTitle(" x [#mum]"); + histogramstructpointer->histAvgCluster->GetXaxis()->SetLabelFont(42); + histogramstructpointer->histAvgCluster->GetXaxis()->SetLabelOffset(0.004); + histogramstructpointer->histAvgCluster->GetXaxis()->SetLabelSize(0.025); + histogramstructpointer->histAvgCluster->GetXaxis()->SetTitleSize(0.035); + histogramstructpointer->histAvgCluster->GetXaxis()->SetTitleFont(42); + histogramstructpointer->histAvgCluster->GetXaxis()->SetTitleOffset(1.5); + histogramstructpointer->histAvgCluster->GetYaxis()->SetTitle(" y [#mum]"); + histogramstructpointer->histAvgCluster->GetYaxis()->SetLabelFont(42); + histogramstructpointer->histAvgCluster->GetYaxis()->SetLabelOffset(-0.016); + histogramstructpointer->histAvgCluster->GetYaxis()->SetLabelSize(0.025); + histogramstructpointer->histAvgCluster->GetYaxis()->SetTitleSize(0.035); + histogramstructpointer->histAvgCluster->GetYaxis()->SetTitleFont(42); + histogramstructpointer->histAvgCluster->GetYaxis()->SetTitleOffset(1.5); + if (histogramstructpointer->calibrated) + histogramstructpointer->histAvgCluster->GetZaxis()->SetTitle(" Charge [e]"); + else + histogramstructpointer->histAvgCluster->GetZaxis()->SetTitle(" Charge [ADU]"); + histogramstructpointer->histAvgCluster->GetZaxis()->SetTitleOffset(1.5); + histogramstructpointer->histAvgCluster->GetZaxis()->SetLabelFont(42); + histogramstructpointer->histAvgCluster->GetZaxis()->SetLabelSize(0.025); + histogramstructpointer->histAvgCluster->GetZaxis()->SetTitleSize(0.035); + histogramstructpointer->histAvgCluster->GetZaxis()->SetTitleFont(42); + histogramstructpointer->histAvgCluster->GetZaxis()->SetNdivisions(5,5,0); + histogramstructpointer->histAvgCluster->Draw("LEGO2Z "); + + TPaveText *pt = new TPaveText(0.3177154,0.9348119,0.6822846,0.995,"blNDC"); + pt->SetName("title"); + pt->SetBorderSize(0); + pt->SetFillColor(0); + pt->SetFillStyle(0); + pt->SetTextFont(42); + pt->Draw(); + Canvas_1->Modified(); + Canvas_1->cd(); + Canvas_1->SetSelected(Canvas_1); + + + TImageDump *img = new TImageDump(savepathresults + "/" + canvastitle + ".png"); + Canvas_1->Paint(); + img->Close(); + + TFile *f = new TFile(savepathresults + "/" + canvastitle + ".root","RECREATE"); + f->cd(); + f->Append(Canvas_1); + f->Append(img); + f->Write(); + + write2DHistogramToFile(histogramstructpointer->histAvgCluster); + // gStyle->SetPaperSize(10.,10.); + // canvas->Print(savepathresults + "/" + canvastitle + ".tex"); + + return 0; +} + +Bool_t Run::plotClusterDistribution() +{ + if (!error) + { + plotClusterDistribution(&histogram); + return 0; + } + return 1; +} + Bool_t Run::plotAllHistograms() { if (!error) @@ -1185,14 +1554,33 @@ Bool_t Run::writeHistogramToFile(TH1F* onehistogram) *fout << header << endl; TString outline; - for(Int_t bin=0;binGetNbinsX();bin++) + for(Int_t bin=1;bin<=onehistogram->GetNbinsX();bin++) { outline=Form("%.3f\t%.3f", onehistogram->GetBinCenter(bin), onehistogram->GetBinContent(bin)); *fout<close(); - - fout = new fstream(filename,ios::out); + return 0; +} + +Bool_t Run::write2DHistogramToFile(TH2F* onehistogram) +{ + system("mkdir "+ savepathresults + " -p"); + TString filename= savepathresults + "/" + runcode + " " + onehistogram->GetName() + " hist.dat"; + fstream* fout = new fstream(filename,ios::out); + + TString header = Form("#x\ty\tCounts"); + *fout << header << endl; + + TString outline; + for (Int_t biny_i = 1; biny_i<=onehistogram->GetXaxis()->GetNbins(); biny_i++) { + for (Int_t binx_i = 1; binx_i<=onehistogram->GetXaxis()->GetNbins(); binx_i++) { + outline=Form("%.2f\t%.2f\t%.2f", onehistogram->GetXaxis()->GetBinCenter(binx_i), onehistogram->GetYaxis()->GetBinCenter(biny_i), onehistogram->GetBinContent(binx_i, biny_i)); + cout << outline << endl; // debug line + *fout<close(); return 0; } @@ -1369,6 +1757,24 @@ void Run::initHistogram(TH1F* &histogrampointer, TString prefix) histogrampointer->GetYaxis()->CenterTitle(); } +void Run::initClusters(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) +{ + initCluster(&histogram, "", xcoord_min_step, xcoord_abs_min, xcoord_abs_max, ycoord_min_step, ycoord_abs_min, ycoord_abs_max); + initCluster(&histogramthreshold, "_threshold", xcoord_min_step, xcoord_abs_min, xcoord_abs_max, ycoord_min_step, ycoord_abs_min, ycoord_abs_max); +} + +void Run::initCluster(histogramstruct* histogramstructpointer, 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) +{ + TString prefix = "Cluster" + suffix; + TString humanreadablestr = Form("%s, %s cluster, %s, chip %s, %s, T=%.1f",prefix.Data(), labbook.source.Data(), labbook.chipGen.Data(), labbook.chip.Data(), labbook.matrix.Data(), labbook.temp); + Int_t sizearrrotx= (int)((xcoord_abs_max-xcoord_abs_min)/xcoord_min_step+0.5); + Int_t sizearrroty = (int)((ycoord_abs_max-ycoord_abs_min)/ycoord_min_step+0.5); +// cout << "xcoord_abs_min: " << xcoord_abs_min << ", xcoord_abs_max: " << xcoord_abs_max << ", xcoord_min_step: " << xcoord_min_step << endl; +// cout << "ycoord_abs_min: " << ycoord_abs_min << ", ycoord_abs_max: " << ycoord_abs_max << ", ycoord_min_step: " << ycoord_min_step << endl; + histogramstructpointer->histAvgCluster = new TH2F(prefix.Data(), Form("Avg. cluster distribution; x; y %s", humanreadablestr.Data()),sizearrrotx, xcoord_abs_min,xcoord_abs_max,sizearrroty,ycoord_abs_min,ycoord_abs_max); + // you can increase the resolution, especially usefull if you observe an rotated (originally staggered) pixel cluster + // histogramstructpointer->histAvgCluster = new TH2F(prefix.Data(), Form("Avg. cluster distribution; x; y %s", humanreadablestr.Data()),2*sizearrrotx, xcoord_abs_min,xcoord_abs_max,sizearrroty*2,ycoord_abs_min,ycoord_abs_max); +} void Run::initRootParameters() { rootcolors = new Int_t[13]{1, 2, 4, 6, 8, 13, 46, 28, 32, 33, 12, 20, 40}; diff --git a/MABS_run_analyzer/Run.h b/MABS_run_analyzer/Run.h index 5c48ec3..b1133d0 100644 --- a/MABS_run_analyzer/Run.h +++ b/MABS_run_analyzer/Run.h @@ -17,9 +17,14 @@ #include #include #include +#include #include #include +#include +#include +#include + #include "help.h" @@ -126,7 +131,7 @@ private: * The labbook data is used to make a string with useful information about the run */ Bool_t generateReadableRunCode(); - + /// is set to true if an entry was found in the SQL database Bool_t runexistsAsRootFile = 0; @@ -144,9 +149,11 @@ private: /** * @brief fills Seedm, Sum and Veto #histogram */ Bool_t binSeedSumVeto(); + /** + * @brief fills the average cluster distribution into a 2 dimension histogram #histogram*/ + Bool_t binCluster(); /// noise quantiles: mean value, sigma in postive direction and sigma in negative direction - Double_t noisequantiles[3]; - + Double_t noisequantiles[3]; TCanvas* plot1DHistogram(TH1F* onehistogram, TString fitFuncType = "landau", TString titlestr = "", TString legendstr = ""); @@ -189,6 +196,10 @@ private: * @brief rescales one specific histogram from ADU to electrons */ void rescaleHistogram(TH1F* &histogrampointernew, TH1F* &histogrampointerold ); + /** + * @brief rescales one specific histogram bin content from ADU to electrons */ + void rescale2DHistogramCounts(TH2F* &histogrampointernew, TH2F* &histogrampointerold ); + /** * @brief calculates Charge Collection Efficiency if Fe55 run */ Bool_t calculteCCE(); @@ -203,10 +214,14 @@ private: */ struct pixelinfo { - Int_t ndiods; - Float_t diodesize; + Int_t ndiods = -1; + Float_t diodesize = -1; + Float_t pitchX = -1; + Float_t pitchY = -1; + Bool_t staggered = 0; /// TODO: add more }; + pixelinfo curpixelinfo; /** @brief system specific constants */ @@ -246,6 +261,21 @@ private: /// init one specific histogram void initHistogram(TH1F* &histogrampointer, TString prefix); + /** @brief calls #initCluster for each structure of type @c Run::histogram + * + * For the initialization parameters calculated by #findRotatedClusterDimension() are used. + */ + void initClusters(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); + + /** + * @brief Find dimensions and borders of rotated cluster + * + * This function ias needed to find out the + * final cluster size and its dimensions, to properly initialize the + * TH2F cluster bins + */ + Bool_t findRotatedClusterDimension(Float_t rotateangle = 45); + void initRootParameters(); Bool_t dynamicalNoise = 1; @@ -284,8 +314,14 @@ public: /** @brief Plot all histograms from @c Run::histogramthreshold into one canvas */ Bool_t plotAllHistogramsThresholdClusterCalibrated(); - /** @brief Writes a given histogram into a file */ + /** @brief Plot average cluster distribution from @c Run::histogram into one canvas */ + Bool_t plotClusterDistribution(); + + /** @brief Writes a given 2D-histogram into a file */ Bool_t writeHistogramToFile(TH1F* onehistogram); + + /** @brief Writes a given histogram into a file */ + Bool_t write2DHistogramToFile(TH2F* onehistogram); /** @brief Writes all histogram of @c Run::histogramsctruct into a file */ Bool_t writeAllHistogramsToFile(); @@ -452,6 +488,12 @@ public: /// set to true, if threshold clusters are used Bool_t thresholdcluster = false; + + /// histogram with cluster data, rotated if staggered + TH2F* histAvgCluster; + + /// number of hits/clusters used to generate all distributions + Double_t numberofhits = 0; }; histogramstruct histogram; histogramstruct histogramCalibrated; @@ -463,9 +505,15 @@ public: /** @brief Plots all histograms from #histogram into one canvas */ Bool_t plotAllHistograms(histogramstruct*); + /** @brief Plots average cluster charge distribution */ + Bool_t plotClusterDistribution(histogramstruct*); + /// init all histograms, set binwidth, bin amount and names void initHistograms(histogramstruct*, TString suffix = ""); + /** @brief initializes the TH2F cluster for the binning for a specific structure of type #histogramstruct one points to*/ + void initCluster(histogramstruct*, 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); + pixelinfo pixelinfoMi34[32]; /// analyze only half of matrix -- 2.43.0