Belle II Software development
CDCDedx2DCellAlgorithm.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <cdc/calibration/CDCdEdx/CDCDedx2DCellAlgorithm.h>
10
11#include <cdc/dbobjects/CDCDedx2DCell.h>
12
13#include <TBox.h>
14#include <TCanvas.h>
15#include <TH2F.h>
16#include <TLine.h>
17#include <TMath.h>
18#include <TStyle.h>
19
20#include <fstream>
21
22using namespace Belle2;
23
24
25//-----------------------------------------------------------------
26// Implementation
27//-----------------------------------------------------------------
28
30 CalibrationAlgorithm("CDCDedxElectronCollector"),
31 fnEntaBinG(128),
32 fnDocaBinG(64),
33 fnEntaBinL(64),
34 fnDocaBinL(28),
35 feaLE(-TMath::Pi() / 2),
36 feaUE(+TMath::Pi() / 2),
37 fdocaLE(-1.50),
38 fdocaUE(1.50),
39 fSetPrefix("_it0"),
40 IsLocalBin(true),
41 IsMakePlots(false),
42 IsRS(true)
43{
44 // Set module properties
45 setDescription("A calibration algorithm for the CDC dE/dx two dimensional correction");
46}
47
48//-----------------------------------------------------------------
49// Run the calibration
50//-----------------------------------------------------------------
51
53{
54
55 //reading electron collector TREE
56 auto ttree = getObjectPtr<TTree>("tree");
57 if (ttree->GetEntries() < 100)return c_NotEnoughData;
58
59 std::vector<double>* dedxhit = 0, *doca = 0, *enta = 0;
60 std::vector<int>* layer = 0;
61
62 ttree->SetBranchAddress("dedxhit", &dedxhit);
63 ttree->SetBranchAddress("layer", &layer);
64 ttree->SetBranchAddress("docaRS", &doca);
65 ttree->SetBranchAddress("entaRS", &enta);
66
67 // Setting up bins for doca and entra angle
70
71 std::vector<int> globalbinsEnta, globalbinsDoca;
72 for (int ibin = 0; ibin < fnEntaBinG; ibin++)globalbinsEnta.push_back(ibin);
73 for (int ibin = 0; ibin < fnDocaBinG; ibin++)globalbinsDoca.push_back(ibin);
74
75 if (!IsLocalBin) {
76 fEntaBinNums = globalbinsEnta;
78 } else {
80 fnEntaBinL = fEntaBinNums.at(fEntaBinNums.size() - 1) + 1;
81 }
82
83 fDocaBinNums = globalbinsDoca;
85
86 std::vector<std::vector<TH1F*>> hILdEdxhitInEntaDocaBin(fnEntaBinL, std::vector<TH1F*>(fnDocaBinL, 0));
87 std::vector<std::vector<TH1F*>> hOLdEdxhitInEntaDocaBin(fnEntaBinL, std::vector<TH1F*>(fnDocaBinL, 0));
88 Double_t ifeaLE = 0, ifeaUE = 0, ifdocaLE = 0, ifdocaUE = 0;
89
90 for (int iea = 0; iea < fnEntaBinL; iea++) {
91 if (IsLocalBin) {
92 ifeaLE = fEntaBinValues.at(iea);
93 ifeaUE = fEntaBinValues.at(iea + 1);
94 } else {
95 ifeaLE = iea * feaBS - feaUE; //- because of -ive range shifting
96 ifeaUE = ifeaLE + feaBS;
97 }
98
99 for (int idoca = 0; idoca < fnDocaBinL; idoca++) {
100
101 ifdocaLE = idoca * fdocaBS - fdocaUE;
102 ifdocaUE = ifdocaLE + fdocaBS;
103
104 hILdEdxhitInEntaDocaBin[iea][idoca] = new TH1F(Form("hILdEdxhitInEntaBin%dDocaBin%d", iea, idoca), "bla-bla", 500, 0., 5.);
105 hILdEdxhitInEntaDocaBin[iea][idoca]->SetTitle(Form("IL: EntA = (%0.03f to %0.03f) and Doca = (%0.03f to %0.03f)", ifeaLE,
106 ifeaUE, ifdocaLE, ifdocaUE));
107 hILdEdxhitInEntaDocaBin[iea][idoca]->GetXaxis()->SetTitle("dedxhits in Inner Layer");
108 hILdEdxhitInEntaDocaBin[iea][idoca]->GetYaxis()->SetTitle("Entries");
109
110 hOLdEdxhitInEntaDocaBin[iea][idoca] = new TH1F(Form("hOLdEdxhitInEntaBin%dDocaBin%d", iea, idoca), "bla-bla", 500, 0., 5.);
111 hOLdEdxhitInEntaDocaBin[iea][idoca]->SetTitle(Form("OL: EntA = (%0.03f to %0.03f) and Doca = (%0.03f to %0.03f)", ifeaLE,
112 ifeaUE, ifdocaLE, ifdocaUE));
113 hOLdEdxhitInEntaDocaBin[iea][idoca]->GetXaxis()->SetTitle("dedxhits in Outer Layer");
114 hOLdEdxhitInEntaDocaBin[iea][idoca]->GetYaxis()->SetTitle("Entries");
115
116 }
117 }
118
119 //Doca vs Enta stats
120 TH2D* hILDocaEntaG = new TH2D("hILDocaEntaG", "Doca vs EntA: Inner Layer", fnDocaBinG, fdocaLE, fdocaUE, fnEntaBinG, feaLE, feaUE);
121 hILDocaEntaG->GetXaxis()->SetTitle("Doca");
122 hILDocaEntaG->GetYaxis()->SetTitle("Entrance angle (#theta)");
123
124 TH2D* hOLDocaEntaG = new TH2D("hOLDocaEntaG", "Doca vs EntA: Outer Layer", fnDocaBinG, fdocaLE, fdocaUE, fnEntaBinG, feaLE, feaUE);
125 hOLDocaEntaG->GetXaxis()->SetTitle("Doca");
126 hOLDocaEntaG->GetYaxis()->SetTitle("Entrance angle (#theta)");
127
128 //when local enta angle is demanded
129 Double_t* RmapEntaValue = &fEntaBinValues[0];
130 TH2D* hILDocaEntaL = new TH2D("hILDocaEntaL", "Doca vs EntA: Inner Layer (rebin)", fnDocaBinL, fdocaLE, fdocaUE, fnEntaBinL,
131 RmapEntaValue);
132 hILDocaEntaL->GetXaxis()->SetTitle("Doca");
133 hILDocaEntaL->GetYaxis()->SetTitle("Entrance angle (#theta)");
134
135 TH2D* hOLDocaEntaL = new TH2D("hOLDocaEntaL", "Doca vs EntA: Outer Layer (rebin)", fnDocaBinL, fdocaLE, fdocaUE, fnEntaBinL,
136 RmapEntaValue);
137 hOLDocaEntaL->GetXaxis()->SetTitle("Doca");
138 hOLDocaEntaL->GetYaxis()->SetTitle("Entrance angle (#theta)");
139
140 TH1D* hILdEdx_all = new TH1D("hILdEdx_all", "", 500, 0., 5.);
141 TH1D* hOLdEdx_all = new TH1D("hOLdEdx_all", "", 500, 0., 5.);
142
143
144 Int_t ibinEA = 0, ibinDOCA = 0;
145 for (int i = 0; i < ttree->GetEntries(); ++i) {
146
147 ttree->GetEvent(i);
148
149 for (unsigned int j = 0; j < dedxhit->size(); ++j) {
150
151 if (dedxhit->at(j) == 0) continue;
152
153 Double_t ieaHit = enta->at(j);
154 if (ieaHit < -TMath::Pi() / 2.0) ieaHit += TMath::Pi() / 2.0;
155 else if (ieaHit > TMath::Pi() / 2.0) ieaHit -= TMath::Pi() / 2.0;
156 if (abs(ieaHit) > TMath::Pi() / 2.0) continue;
157
158 Double_t idocaHit = doca->at(j);
159 if (abs(idocaHit) > 1.50) continue;
160
161 //Bin corresponds to enta and doca value
162 ibinEA = (ieaHit - feaLE) / feaBS ; //from 0
163
164 ibinDOCA = (idocaHit - fdocaLE) / fdocaBS;
165
166 if (ibinEA >= fnEntaBinG || ibinDOCA >= fnDocaBinG) continue; //bin stats from 0
167
168 if (IsLocalBin) ibinEA = fEntaBinNums.at(ibinEA);
169
170 if (layer->at(j) < 8) {
171 hILDocaEntaG->Fill(idocaHit, ieaHit);
172 if (IsLocalBin)hILDocaEntaL->Fill(idocaHit, ieaHit);
173 hILdEdx_all->Fill(dedxhit->at(j));
174 hILdEdxhitInEntaDocaBin[ibinEA][ibinDOCA]->Fill(dedxhit->at(j));
175 } else {
176 hOLDocaEntaG->Fill(idocaHit, ieaHit);
177 if (IsLocalBin)hOLDocaEntaL->Fill(idocaHit, ieaHit);
178 hOLdEdx_all->Fill(dedxhit->at(j));
179 hOLdEdxhitInEntaDocaBin[ibinEA][ibinDOCA]->Fill(dedxhit->at(j));
180 }
181 }
182 }
183
184
185 if (IsMakePlots) {
186 TCanvas* ctmpde = new TCanvas("DocavsEnta", "Doca vs Enta distributions", 800, 400);
187 if (IsLocalBin) {
188 ctmpde->SetCanvasSize(800, 800);
189 ctmpde->Divide(2, 2);
190 ctmpde->cd(1); hILDocaEntaG->Draw("colz");
191 ctmpde->cd(2); hOLDocaEntaG->Draw("colz");
192 ctmpde->cd(3); hILDocaEntaL->Draw("colz");
193 ctmpde->cd(4); hOLDocaEntaL->Draw("colz");
194 } else {
195 ctmpde->Divide(2, 1);
196 ctmpde->cd(1); hILDocaEntaG->Draw("colz");
197 ctmpde->cd(2); hOLDocaEntaG->Draw("colz");
198 }
199 ctmpde->SaveAs(Form("DocavsEnta_%s.pdf", fSetPrefix.data()));
200 ctmpde->SaveAs(Form("DocavsEnta_%s.root", fSetPrefix.data()));
201
202 TCanvas* ctem = new TCanvas("Layerhisto", "Inner and Outer Histo", 600, 600);
203 hOLdEdx_all->Draw("histo");
204 hILdEdx_all->SetMarkerColor(kRed);
205 hILdEdx_all->Draw("same histo");
206 ctem->SaveAs(Form("Layerhistodedxhit_TwoDCorr_%s.pdf", fSetPrefix.data()));
207 }
208
209
210 //Calculating 5-75% global truncation mean
211 double InsumPer5 = 0.0, InsumPer75 = 0.0;
212 double OutsumPer5 = 0.0, OutsumPer75 = 0.0;
213 double InLayInt = hILdEdx_all->Integral();
214 double OutLayInt = hOLdEdx_all->Integral();
215
216 Int_t lBinInLayer = 1, hBinInLayer = 1;
217 Int_t lBinOutLayer = 1, hBinOutLayer = 1;
218
219 for (int ibin = 1; ibin <= hILdEdx_all->GetNbinsX(); ibin++) {
220
221 if (InsumPer5 <= 0.05 * InLayInt) {
222 InsumPer5 += hILdEdx_all->GetBinContent(ibin);
223 lBinInLayer = ibin;
224 }
225
226 if (InsumPer75 <= 0.75 * InLayInt) {
227 InsumPer75 += hILdEdx_all->GetBinContent(ibin);
228 hBinInLayer = ibin;
229 }
230
231 if (OutsumPer5 <= 0.05 * OutLayInt) {
232 OutsumPer5 += hOLdEdx_all->GetBinContent(ibin);
233 lBinOutLayer = ibin;
234 }
235
236 if (OutsumPer75 <= 0.75 * OutLayInt) {
237 OutsumPer75 += hOLdEdx_all->GetBinContent(ibin);
238 hBinOutLayer = ibin;
239 }
240 }
241
242 short version = 1;
243 TCanvas* ctmp = new TCanvas("tmp", "tmp", 1200, 1200);
244 ctmp->Divide(4, 4);
245 std::stringstream psname; psname << Form("dedx_2dcell_%s.pdf[", fSetPrefix.data());
246 TLine* tl = new TLine();
247 tl->SetLineColor(kRed);
248
249 TBox* tb = new TBox();
250
251 if (IsMakePlots) {
252 ctmp->Print(psname.str().c_str());
253 psname.str(""); psname << Form("dedx_2dcell_%s.pdf", fSetPrefix.data());
254 }
255
256 gStyle->SetOptStat("nemriou");
257 TH1D* htemp = NULL;
258 std::vector<TH2F> twodcors;
259 TH2F tempTwoD = TH2F("tempTwoD", "dE/dx in bins of DOCA/Enta;DOCA;Entrance Angle", fnDocaBinL, fdocaLE, fdocaUE, fnEntaBinL,
260 RmapEntaValue);
261 TH2F twodcor = TH2F("twodcorrection", "dE/dx in bins of DOCA/Enta;DOCA;Entrance Angle", fnDocaBinG, fdocaLE, fdocaUE, fnEntaBinG,
262 feaLE, feaUE);
263
264 std::ofstream file2DILout;
265
266 for (int iIOLayer = 0; iIOLayer <= 1; iIOLayer++) {
267
268 file2DILout.open(Form("f2D_Constants_Local_Layer%d_%s.txt", iIOLayer, fSetPrefix.data()));
269
270 int startfrom = 1, endat = 1;
271 if (iIOLayer == 0) {
272 startfrom = lBinInLayer; endat = hBinInLayer;
273 } else {
274 startfrom = lBinOutLayer; endat = hBinOutLayer;
275 }
276
277 //std::cout << "Layer I/O # = " << iIOLayer << std::endl;
278 for (int iea = 1; iea <= fnEntaBinL; iea++) {
279
280 Int_t ieaprime = iea; //rotation symmtery for 1<->3 and 4<->2
281 if (IsRS)ieaprime = GetRotationSymmericBin(fnEntaBinL, iea);
282
283 for (int idoca = 1; idoca <= fnDocaBinL; idoca++) {
284
285 if (iIOLayer == 0)
286 htemp = (TH1D*)hILdEdxhitInEntaDocaBin[ieaprime - 1][idoca - 1]->Clone(Form("hL%d_Ea%d_Doca%d", iIOLayer, iea, idoca));
287 else if (iIOLayer == 1)
288 htemp = (TH1D*)hOLdEdxhitInEntaDocaBin[ieaprime - 1][idoca - 1]->Clone(Form("hL%d_Ea%d_Doca%d", iIOLayer, iea, idoca));
289 else continue;
290
291 double truncMean = 1.0;
292 if (htemp->GetEntries() < 400) truncMean = 1.0; //low stats
293 else {
294 double binweights = 0.0;
295 int sumofbc = 0;
296 for (int ibin = startfrom; ibin <= endat; ibin++) {
297 //std::cout << " dedxhit bin = " << ibin << ", Entries =" << htemp->GetBinContent(ibin) << std::endl;
298 if (htemp->GetBinContent(ibin) > 0) {
299 binweights += (htemp->GetBinContent(ibin) * htemp->GetBinCenter(ibin));
300 sumofbc += htemp->GetBinContent(ibin);
301 }
302 }
303 if (sumofbc > 0)truncMean = (double)(binweights / sumofbc);
304 else truncMean = 1.0;
305 }
306
307 if (truncMean <= 0)truncMean = 1.0; //protection only
308 tempTwoD.SetBinContent(idoca, iea, truncMean); //binning starts at 1
309
310 if (IsMakePlots) {
311 ctmp->cd(((idoca - 1) % 16) + 1);
312 htemp->SetTitle(Form("%s, mean = %0.4f", htemp->GetTitle(), truncMean));
313 htemp->SetFillColorAlpha(kGreen, 0.30);
314 if (truncMean >= 1.02 || truncMean <= 0.98)htemp->SetFillColor(kYellow);
315 if (truncMean >= 1.05 || truncMean <= 0.95)htemp->SetFillColor(kOrange);
316 if (truncMean >= 1.10 || truncMean <= 0.90)htemp->SetFillColor(kRed);
317 htemp->DrawClone(); //clone is nessesory for pointer survival
318 tl->SetLineColor(kRed);
319 tl->SetX1(truncMean); tl->SetX2(truncMean);
320 tl->SetY1(0); tl->SetY2(htemp->GetMaximum());
321 tl->DrawClone("same");
322
323 tb->SetLineColor(kPink);
324 tb->SetLineStyle(6);
325 tb->SetFillColorAlpha(kPink, 0.15);
326 tb->SetX1(htemp->GetBinLowEdge(startfrom));
327 tb->SetY1(0);
328 tb->SetX2(htemp->GetBinLowEdge(endat));
329 tb->SetY2(htemp->GetMaximum());
330 tb->DrawClone("same");
331 if (idoca % 16 == 0)ctmp->Print(psname.str().c_str());
332 }
333
334 file2DILout << iea << ", \t" << idoca << ", \t" << truncMean << "\n";
335 htemp->Reset();
336 }
337 }
338
339 //conversion from local to global
340 for (int iea = 0; iea < fnEntaBinG; iea++) {
341 ibinEA = iea;
342 if (IsLocalBin)ibinEA = fEntaBinNums.at(iea);
343 for (int idoca = 0; idoca < fnDocaBinG; idoca++) {
344 ibinDOCA = idoca;
345 twodcor.SetBinContent(idoca + 1, iea + 1, tempTwoD.GetBinContent(ibinDOCA + 1, ibinEA + 1));
346 }
347 }
348
349 if (iIOLayer == 0)twodcor.SetTitle("InnerLayer: dE/dx in bins of DOCA/Enta");
350 else if (iIOLayer == 1)twodcor.SetTitle("OuterLayer: dE/dx in bins of DOCA/Enta");
351 twodcors.push_back(twodcor);
352 twodcor.Reset();
353 tempTwoD.Reset();
354 file2DILout.close();
355 }
356
357 if (IsMakePlots) {
358 psname.str(""); psname << Form("dedx_2dcell_%s.pdf]", fSetPrefix.data());
359 ctmp->Print(psname.str().c_str());
360
361 // //Drawing final constants
362 TCanvas* cconst = new TCanvas("FinalConstantHistoMap", "Inner and Outer Histo", 800, 400);
363 cconst->Divide(2, 1);
364 cconst->cd(1); twodcors.at(0).Draw("colz");
365 cconst->cd(2); twodcors.at(1).Draw("colz");
366 cconst->SaveAs(Form("Final2DConstantMap_wGlobalBin_%s.pdf", fSetPrefix.data()));
367 cconst->SaveAs(Form("Final2DConstantMap_wGlobalBin_%s.root", fSetPrefix.data()));
368 }
369
370
371 B2INFO("dE/dx calibration done for 2D correction");
372 CDCDedx2DCell* gain = new CDCDedx2DCell(version, twodcors);
373 saveCalibration(gain, "CDCDedx2DCell");
374
375 delete htemp;
376 delete ctmp;
377 delete tl;
378 delete tb;
379 return c_OK;
380
381}
bool IsRS
if rotation symmtery requested
bool IsLocalBin
if local variable bin requested
std::vector< double > fEntaBinValues
Vector for doca asym bin values.
Double_t fdocaUE
Upper edge of doca.
Double_t fdocaLE
Lower edge of doca.
int fnDocaBinG
doca angle bins, Global
Double_t feaBS
Binwidth edge of enta angle.
std::vector< int > fEntaBinNums
Vector for enta asym bin values.
int fnDocaBinL
doca angle bins, Local
CDCDedx2DCellAlgorithm()
Constructor: Sets the description the properties and the parameters of the algorithm.
int fnEntaBinL
etna angle bins, Local
std::string fSetPrefix
prefix to filename
Double_t fdocaBS
Binwidth edge of doca.
int GetRotationSymmericBin(int nbin, int ibin)
function to set rotation symmetry
bool IsMakePlots
produce plots for status
virtual EResult calibrate() override
2D Cell algorithm algorithm
int fnEntaBinG
Save arithmetic and truncated mean for the 'dedx' values.
std::vector< int > fDocaBinNums
Vector for enta asym bin #.
Double_t feaLE
Lower edge of enta angle.
void GetVariableBin(int nbin, std::vector< int > &nBinEnta0to100Per)
function to set variable bins
Double_t feaUE
Upper edge of enta angle.
dE/dx wire gain calibration constants
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
void setDescription(const std::string &description)
Set algorithm description (in constructor)
EResult
The result of calibration.
@ c_OK
Finished successfully =0 in Python.
@ c_NotEnoughData
Needs more data =2 in Python.
CalibrationAlgorithm(const std::string &collectorModuleName)
Constructor - sets the prefix for collected objects (won't be accesses until execute(....
std::shared_ptr< T > getObjectPtr(const std::string &name, const std::vector< Calibration::ExpRun > &requestedRuns)
Get calibration data object by name and list of runs, the Merge function will be called to generate t...
Abstract base class for different kinds of events.