Belle II Software  release-05-02-19
CDCDedx1DCellAlgorithm.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2016 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: jikumar, jvbennett *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <algorithm>
12 #include <iostream>
13 
14 #include <TCanvas.h>
15 #include <TH1D.h>
16 #include <TLine.h>
17 #include <TTree.h>
18 #include <TMath.h>
19 
20 #include <reconstruction/calibration/CDCDedx1DCellAlgorithm.h>
21 
22 using namespace Belle2;
23 
24 
25 //-----------------------------------------------------------------
26 // Implementation
27 //-----------------------------------------------------------------
28 
30  CalibrationAlgorithm("CDCDedxElectronCollector"),
31  fnEntaBinG(128),
32  fnEntaBinL(64),
33  feaLE(-TMath::Pi() / 2),
34  feaUE(+TMath::Pi() / 2),
35  fSetPrefix("_it0"),
36  IsLocalBin(true),
37  IsMakePlots(false),
38  IsRS(true)
39 {
40  // Set module properties
41  setDescription("A calibration algorithm for the CDC dE/dx entrance angle cleanup correction");
42 }
43 
44 //-----------------------------------------------------------------
45 // Run the calibration
46 //-----------------------------------------------------------------
47 
49 {
50 
51  //reading electron collector TREE
52  auto ttree = getObjectPtr<TTree>("tree");
53  if (ttree->GetEntries() < 100)return c_NotEnoughData;
54 
55  std::vector<double>* dedxhit = 0, *enta = 0;
56  std::vector<int>* layer = 0;
57 
58  ttree->SetBranchAddress("dedxhit", &dedxhit);
59  ttree->SetBranchAddress("layer", &layer);
60  ttree->SetBranchAddress("entaRS", &enta);
61 
62  // Setting up bins for entra angle
63  feaBS = (feaUE - feaLE) / fnEntaBinG;
64 
65  std::vector<int> globalbins;
66  for (int ibin = 0; ibin < fnEntaBinG; ibin++)globalbins.push_back(ibin);
67 
68  if (!IsLocalBin) {
69  fEntaBinNums = globalbins;
71  } else {
73  fnEntaBinL = fEntaBinNums.at(fEntaBinNums.size() - 1) + 1;
74  }
75 
76  //enta dedx distributions for inner and outer layer
77  std::vector<TH1D*> hILdEdxhitInEntaBin(fnEntaBinL, 0);
78  std::vector<TH1D*> hOLdEdxhitInEntaBin(fnEntaBinL, 0);
79 
80  Double_t ifeaLE = 0, ifeaUE = 0;
81 
82  for (int iea = 0; iea < fnEntaBinL; iea++) {
83 
84  if (IsLocalBin) {
85  ifeaLE = fEntaBinValues.at(iea);
86  ifeaUE = fEntaBinValues.at(iea + 1);
87  } else {
88  ifeaLE = iea * feaBS - feaUE; //- because of -ive range shifting
89  ifeaUE = ifeaLE + feaBS;
90  }
91 
92  hILdEdxhitInEntaBin[iea] = new TH1D(Form("hILdEdxhitInEntaBin%d", iea), "bla-bla", 250, 0., 5.);
93  hILdEdxhitInEntaBin[iea]->SetTitle(Form("IL: dedxhit in EntA = (%0.03f to %0.03f)", ifeaLE, ifeaUE));
94  hILdEdxhitInEntaBin[iea]->GetXaxis()->SetTitle("dedxhits in Inner Layer");
95  hILdEdxhitInEntaBin[iea]->GetYaxis()->SetTitle("Entries");
96 
97  hOLdEdxhitInEntaBin[iea] = new TH1D(Form("hOLdEdxhitInEntaBin%d", iea), "bla-bla", 250, 0., 5.);
98  hOLdEdxhitInEntaBin[iea]->SetTitle(Form("OL: dedxhit in EntA = (%0.03f to %0.03f)", ifeaLE, ifeaUE));
99  hOLdEdxhitInEntaBin[iea]->GetXaxis()->SetTitle("dedxhits in Outer Layer");
100  hOLdEdxhitInEntaBin[iea]->GetYaxis()->SetTitle("Entries");
101  }
102 
103  // Enta stats
104  TH1D* hILEntaG = new TH1D("hILEntaG", "EntA: Inner Layer", fnEntaBinG, feaLE, feaUE);
105  hILEntaG->GetXaxis()->SetTitle("Entrance angle (#theta)");
106  hILEntaG->GetYaxis()->SetTitle("Entries");
107 
108  TH1D* hOLEntaG = new TH1D("hOLEntaG", "EntA: Outer Layer", fnEntaBinG, feaLE, feaUE);
109  hOLEntaG->GetXaxis()->SetTitle("Entrance angle (#theta)");
110  hOLEntaG->GetYaxis()->SetTitle("Entries");
111 
112  //rebinned histogram
113  Double_t* RmapEntaValue = &fEntaBinValues[0];
114 
115  TH1D* hILEntaL = new TH1D("hILEntaL", "EntA: Inner Layer (assym bin)", fnEntaBinL, RmapEntaValue);
116  hILEntaL->GetXaxis()->SetTitle("Entrance angle (#theta)");
117  hILEntaL->GetYaxis()->SetTitle("Entries");
118 
119  TH1D* hOLEntaL = new TH1D("hOLEntaL", "EntA: Outer Layer (assym bin)", fnEntaBinL, RmapEntaValue);
120  hOLEntaL->GetYaxis()->SetTitle("Entries");
121  hOLEntaL->GetXaxis()->SetTitle("Entrance angle (#theta)");
122 
123  TH1D* hILdEdx_all = new TH1D("hILdEdx_all", "", 250, 0., 5.);
124  TH1D* hOLdEdx_all = new TH1D("hOLdEdx_all", "", 250, 0., 5.);
125 
126  Int_t ibinEA = 0;
127  for (int i = 0; i < ttree->GetEntries(); ++i) {
128 
129  ttree->GetEvent(i);
130 
131  for (unsigned int j = 0; j < dedxhit->size(); ++j) {
132 
133  if (dedxhit->at(j) == 0) continue;
134 
135  Double_t ieaHit = enta->at(j);
136  if (ieaHit < -TMath::Pi() / 2.0) ieaHit += TMath::Pi() / 2.0;
137  else if (ieaHit > TMath::Pi() / 2.0) ieaHit -= TMath::Pi() / 2.0;
138  if (abs(ieaHit) > TMath::Pi() / 2.0) continue;
139 
140  //Bin corresponds to enta
141  ibinEA = (ieaHit - feaLE) / feaBS ; //from 0
142  if (ibinEA >= fnEntaBinG) continue; //bin stats from 0
143 
144  if (IsLocalBin) ibinEA = fEntaBinNums.at(ibinEA);
145 
146  if (layer->at(j) < 8) {
147  hILEntaG->Fill(ieaHit);
148  if (IsLocalBin)hILEntaL->Fill(ieaHit);
149  hILdEdx_all->Fill(dedxhit->at(j));
150  hILdEdxhitInEntaBin[ibinEA]->Fill(dedxhit->at(j));
151  } else {
152  hOLEntaG->Fill(ieaHit);
153  if (IsLocalBin)hOLEntaL->Fill(ieaHit);
154  hOLdEdx_all->Fill(dedxhit->at(j));
155  hOLdEdxhitInEntaBin[ibinEA]->Fill(dedxhit->at(j));
156  }
157  }
158  }
159 
160  if (IsMakePlots) {
161  TCanvas* ctmpde = new TCanvas("hEntaDist", "Enta distributions", 400, 400);
162  if (IsLocalBin) {
163  ctmpde->SetCanvasSize(800, 400);
164  ctmpde->Divide(2, 1);
165  ctmpde->cd(1); gPad->SetLogy(); hOLEntaG->SetMarkerColor(kBlue); hOLEntaG->Draw(""); hILEntaG->Draw("same");
166  ctmpde->cd(2); gPad->SetLogy(); hOLEntaL->SetMarkerColor(kBlue); hOLEntaL->Draw(""); hILEntaL->Draw("same");
167  } else {
168  hOLEntaG->Draw(""); hILEntaG->Draw("same");
169  }
170  ctmpde->SaveAs(Form("hEntaDistributions_%s.pdf", fSetPrefix.data()));
171 
172  }
173 
174  double InsumPer5 = 0.0, InsumPer75 = 0.0;
175  double OutsumPer5 = 0.0, OutsumPer75 = 0.0;
176  double InLayInt = hILdEdx_all->Integral();
177  double OutLayInt = hOLdEdx_all->Integral();
178 
179  Int_t lBinInLayer = 1, hBinInLayer = 1;
180  Int_t lBinOutLayer = 1, hBinOutLayer = 1;
181 
182  for (int ibin = 1; ibin <= hILdEdx_all->GetNbinsX(); ibin++) {
183 
184  if (InsumPer5 <= 0.05 * InLayInt) {
185  InsumPer5 += hILdEdx_all->GetBinContent(ibin);
186  lBinInLayer = ibin;
187  }
188 
189  if (InsumPer75 <= 0.75 * InLayInt) {
190  InsumPer75 += hILdEdx_all->GetBinContent(ibin);
191  hBinInLayer = ibin;
192  }
193 
194  if (OutsumPer5 <= 0.05 * OutLayInt) {
195  OutsumPer5 += hOLdEdx_all->GetBinContent(ibin);
196  lBinOutLayer = ibin;
197  }
198 
199  if (OutsumPer75 <= 0.75 * OutLayInt) {
200  OutsumPer75 += hOLdEdx_all->GetBinContent(ibin);
201  hBinOutLayer = ibin;
202  }
203  }
204 
205  if (IsMakePlots) {
206 
207  TCanvas* ctem = new TCanvas("Layerhisto", "Inner and Outer Histo", 1200, 600);
208  ctem->Divide(2, 1);
209  ctem->cd(1);
210  hILdEdx_all->SetMarkerColor(kBlue);
211  hILdEdx_all->Draw("histo");
212 
213  TLine* tlInF = new TLine();
214  tlInF->SetLineColor(kBlue);
215  tlInF->SetX1(InsumPer5); tlInF->SetX2(InsumPer5);
216  tlInF->SetY1(0); tlInF->SetY2(hILdEdx_all->GetMaximum());
217  tlInF->DrawClone("same");
218 
219  TLine* tlInL = new TLine();
220  tlInL->SetLineColor(kBlue);
221  tlInL->SetX1(InsumPer75); tlInL->SetX2(InsumPer75);
222  tlInL->SetY1(0); tlInL->SetY2(hILdEdx_all->GetMaximum());
223  tlInL->DrawClone("same");
224 
225  ctem->cd(2);
226  hOLdEdx_all->Draw("histo");
227  hOLdEdx_all->SetMarkerColor(kRed);
228 
229  TLine* tlOutF = new TLine();
230  tlOutF->SetLineColor(kBlue);
231  tlOutF->SetX1(OutsumPer5); tlOutF->SetX2(OutsumPer5);
232  tlOutF->SetY1(0); tlOutF->SetY2(hOLdEdx_all->GetMaximum());
233  tlOutF->DrawClone("same");
234 
235  TLine* tlOutL = new TLine();
236  tlOutL->SetLineColor(kBlue);
237  tlOutL->SetX1(OutsumPer75); tlOutL->SetX2(OutsumPer75);
238  tlOutL->SetY1(0); tlOutL->SetY2(hOLdEdx_all->GetMaximum());
239  tlOutL->DrawClone("same");
240  ctem->SaveAs(Form("Layerhistodedxhit_OneDCorr_%s.pdf", fSetPrefix.data()));
241  }
242 
243  //short version = 0;
244  TCanvas* ctmp = new TCanvas("tmp", "tmp", 1200, 1200);
245  ctmp->Divide(4, 4);
246  std::stringstream psname; psname << Form("dedx_1dcell_%s.pdf[", fSetPrefix.data());
247  TLine* tl = new TLine();
248  tl->SetLineColor(kRed);
249 
250  TBox* tb = new TBox();
251 
252  if (IsMakePlots) {
253  ctmp->Print(psname.str().c_str());
254  psname.str(""); psname << Form("dedx_1dcell_%s.pdf", fSetPrefix.data());
255  }
256 
257  TH1D* htemp = 0x0;
258  std::vector<std::vector<double>> onedcors; // prev->std::vector<std::vector<double>> ones;
259  std::vector<double> onedcorIorOL, onedcorIorOLtemp;
260 
261  TH1D* hILEntaConst = new TH1D("hILEntaConst", "EntA: Outer Layer", fnEntaBinG, feaLE, feaUE);
262  hILEntaConst->GetXaxis()->SetTitle("Entrance angle (#theta)");
263  hILEntaConst->GetYaxis()->SetTitle("Constant");
264 
265  TH1D* hOLEntaConst = new TH1D("hOLEntaConst", "EntA: Outer Layer", fnEntaBinG, feaLE, feaUE);
266  hOLEntaConst->GetXaxis()->SetTitle("Entrance angle (#theta)");
267  hOLEntaConst->GetYaxis()->SetTitle("Constant");
268 
269  for (int iIOLayer = 0; iIOLayer <= 1; iIOLayer++) {
270 
271  int startfrom = 1, endat = 1;
272 
273  if (iIOLayer == 0) {
274  startfrom = lBinInLayer; endat = hBinInLayer;
275  } else if (iIOLayer == 1) {
276  startfrom = lBinOutLayer; endat = hBinOutLayer;
277  } else continue;
278 
279  for (int iea = 1; iea <= fnEntaBinL; iea++) {
280 
281  Int_t ieaprime = iea; //rotation symmtery for 1<->3 and 4<->2
282 
283  if (IsRS)ieaprime = GetRotationSymmericBin(fnEntaBinL, iea);
284 
285  if (iIOLayer == 0)htemp = (TH1D*)hILdEdxhitInEntaBin[ieaprime - 1]->Clone(Form("hL%d_Ea%d", iIOLayer, iea));
286  else if (iIOLayer == 1)htemp = (TH1D*)hOLdEdxhitInEntaBin[ieaprime - 1]->Clone(Form("hL%d_Ea%d", iIOLayer, iea));
287  else continue;
288 
289  double truncMean = 1.0;
290  if (htemp->Integral() < 250) truncMean = 1.0; //low stats
291  else {
292  double binweights = 0.0;
293  int sumofbc = 0;
294  for (int ibin = startfrom; ibin <= endat; ibin++) {
295 
296  if (htemp->GetBinContent(ibin) > 0) {
297  binweights += (htemp->GetBinContent(ibin) * htemp->GetBinCenter(ibin));
298  sumofbc += htemp->GetBinContent(ibin);
299  }
300  }
301  if (sumofbc > 0)truncMean = (double)(binweights / sumofbc);
302  else truncMean = 1.0;
303  }
304 
305  if (truncMean <= 0)truncMean = 1.0; //protection only
306  onedcorIorOLtemp.push_back(truncMean);
307 
308  if (IsMakePlots) {
309  ctmp->cd(((iea - 1) % 16) + 1);
310  htemp->SetFillColorAlpha(kGreen, 0.30);
311  htemp->SetTitle(Form("%s, #mean = %0.2f", htemp->GetTitle(), truncMean));
312  if (truncMean >= 1.02 || truncMean <= 0.98)htemp->SetFillColor(kYellow);
313  if (truncMean >= 1.05 || truncMean <= 0.95)htemp->SetFillColor(kOrange);
314  if (truncMean >= 1.10 || truncMean <= 0.90)htemp->SetFillColor(kRed);
315  htemp->DrawClone("hist"); //clone is nessesory for pointer survival
316  tl->SetLineColor(kRed);
317  tl->SetX1(truncMean); tl->SetX2(truncMean);
318  tl->SetY1(0); tl->SetY2(htemp->GetMaximum());
319  tl->DrawClone("same");
320 
321  tb->SetLineColor(kPink);
322  tb->SetLineStyle(6);
323  tb->SetFillColorAlpha(kPink, 0.15);
324  tb->SetX1(htemp->GetBinLowEdge(startfrom));
325  tb->SetY1(0);
326  tb->SetX2(htemp->GetBinLowEdge(endat));
327  tb->SetY2(htemp->GetMaximum());
328  tb->DrawClone("same");
329  if (iea % 16 == 0)ctmp->Print(psname.str().c_str());
330  }
331 
332  htemp->Reset();
333  }
334 
335  ibinEA = 0;
336  for (int iea = 0; iea < fnEntaBinG; iea++) {
337  ibinEA = iea;
338  if (IsLocalBin)ibinEA = fEntaBinNums.at(iea);
339  onedcorIorOL.push_back(onedcorIorOLtemp.at(ibinEA));
340  if (iIOLayer == 0)hILEntaConst->SetBinContent(iea + 1, onedcorIorOLtemp.at(ibinEA));
341  else if (iIOLayer == 1)hOLEntaConst->SetBinContent(iea + 1, onedcorIorOLtemp.at(ibinEA));
342  }
343 
344  onedcors.push_back(onedcorIorOL);
345  onedcorIorOL.clear();
346  onedcorIorOLtemp.clear();
347  }
348 
349 
350  if (IsMakePlots) {
351  psname.str(""); psname << Form("dedx_1dcell_%s.pdf]", fSetPrefix.data());
352  ctmp->Print(psname.str().c_str());
353 
354  // //Drawing final constants
355  TCanvas* cconst = new TCanvas("FinalConstantHistoMap", "Inner and Outer Histo", 800, 400);
356  cconst->Divide(2, 1);
357  cconst->cd(1); hILEntaConst->Draw();
358  cconst->cd(2); hOLEntaConst->Draw();
359  cconst->SaveAs(Form("FinalOneDConstantMap_%s.pdf", fSetPrefix.data()));
360  }
361 
362  B2INFO("dE/dx Calibration done for 1D cleanup correction");
363  CDCDedx1DCell* gain = new CDCDedx1DCell(0, onedcors);
364  saveCalibration(gain, "CDCDedx1DCell");
365 
366  delete htemp;
367  delete ctmp;
368  delete tl;
369  delete tb;
370  return c_OK;
371 
372 }
Belle2::CDCDedx1DCellAlgorithm::fnEntaBinL
int fnEntaBinL
etna angle bins, Local
Definition: CDCDedx1DCellAlgorithm.h:162
Belle2::CDCDedx1DCellAlgorithm::IsLocalBin
bool IsLocalBin
if local variable bins requested
Definition: CDCDedx1DCellAlgorithm.h:172
Belle2::CDCDedx1DCellAlgorithm::feaBS
Double_t feaBS
Binwidth edge of enta angle.
Definition: CDCDedx1DCellAlgorithm.h:166
Belle2::CalibrationAlgorithm::saveCalibration
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
Definition: CalibrationAlgorithm.cc:290
Belle2::CDCDedx1DCellAlgorithm::fEntaBinNums
std::vector< int > fEntaBinNums
Vector for enta asym bin values.
Definition: CDCDedx1DCellAlgorithm.h:169
Belle2::CDCDedx1DCell
dE/dx wire gain calibration constants
Definition: CDCDedx1DCell.h:37
Belle2::CalibrationAlgorithm::c_OK
@ c_OK
Finished successfuly =0 in Python.
Definition: CalibrationAlgorithm.h:51
Belle2::CalibrationAlgorithm::setDescription
void setDescription(const std::string &description)
Set algorithm description (in constructor)
Definition: CalibrationAlgorithm.h:331
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::CDCDedx1DCellAlgorithm::IsMakePlots
bool IsMakePlots
produce plots for status
Definition: CDCDedx1DCellAlgorithm.h:173
Belle2::CDCDedx1DCellAlgorithm::fEntaBinValues
std::vector< double > fEntaBinValues
Vector for doca asym bin values.
Definition: CDCDedx1DCellAlgorithm.h:170
Belle2::CalibrationAlgorithm::EResult
EResult
The result of calibration.
Definition: CalibrationAlgorithm.h:50
Belle2::CDCDedx1DCellAlgorithm::CDCDedx1DCellAlgorithm
CDCDedx1DCellAlgorithm()
Constructor: Sets the description, the properties and the parameters of the algorithm.
Definition: CDCDedx1DCellAlgorithm.cc:29
Belle2::CDCDedx1DCellAlgorithm::GetVariableBin
void GetVariableBin(int nbin, std::vector< int > &nBinEnta0to100Per)
funtion to set variable bins
Definition: CDCDedx1DCellAlgorithm.h:100
Belle2::CDCDedx1DCellAlgorithm::IsRS
bool IsRS
if rotation symmtery requested
Definition: CDCDedx1DCellAlgorithm.h:174
Belle2::CDCDedx1DCellAlgorithm::feaUE
Double_t feaUE
Upper edge of enta angle.
Definition: CDCDedx1DCellAlgorithm.h:165
Belle2::CDCDedx1DCellAlgorithm::fSetPrefix
std::string fSetPrefix
suffix to filename
Definition: CDCDedx1DCellAlgorithm.h:168
Belle2::CalibrationAlgorithm::c_NotEnoughData
@ c_NotEnoughData
Needs more data =2 in Python.
Definition: CalibrationAlgorithm.h:53
Belle2::CDCDedx1DCellAlgorithm::feaLE
Double_t feaLE
Lower edge of enta angle.
Definition: CDCDedx1DCellAlgorithm.h:164
Belle2::CDCDedx1DCellAlgorithm::fnEntaBinG
int fnEntaBinG
Save arithmetic and truncated mean for the 'dedx' values.
Definition: CDCDedx1DCellAlgorithm.h:161
Belle2::CalibrationAlgorithm
Base class for calibration algorithms.
Definition: CalibrationAlgorithm.h:47
Belle2::CDCDedx1DCellAlgorithm::calibrate
virtual EResult calibrate() override
1D cell algorithm
Definition: CDCDedx1DCellAlgorithm.cc:48
Belle2::CDCDedx1DCellAlgorithm::GetRotationSymmericBin
int GetRotationSymmericBin(int nbin, int ibin)
funtion to set rotation symmetry
Definition: CDCDedx1DCellAlgorithm.h:85