Belle II Software  release-05-02-19
ECLShowerCorrectorModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2016 - Belle II Collaboration *
4  * *
5  * This module performs the energy correction for EM shower (mainly *
6  * longitudinal leakage): corr = raw * correctionfactor. *
7  * *
8  * Author: The Belle II Collaboration *
9  * Contributors: Torben Ferber (ferber@physics.ubc.ca) (TF) *
10  * Alon Hershenhorn (hershen@physics.ubc.ca) (AH) *
11  * Suman Koirala (Suman Koirala <suman@ntu.edu.tw>) (SK) *
12  * *
13  * This software is provided "as is" without any warranty. *
14  **************************************************************************/
15 
16 // THIS MODULE
17 #include <ecl/modules/eclShowerCorrection/ECLShowerCorrectorModule.h>
18 
19 // FRAMEWORK
20 #include <framework/logging/Logger.h>
21 
22 // ROOT
23 #include <TMath.h>
24 
25 // ECL
26 #include <ecl/dataobjects/ECLShower.h>
27 #include <ecl/dbobjects/ECLShowerCorrectorLeakageCorrection.h>
28 #include <ecl/dbobjects/ECLShowerEnergyCorrectionTemporary.h>
29 
30 // MDST
31 #include <mdst/dataobjects/EventLevelClusteringInfo.h>
32 
33 using namespace Belle2;
34 
35 //-----------------------------------------------------------------
36 // Register the Module
37 //-----------------------------------------------------------------
38 REG_MODULE(ECLShowerCorrector)
39 REG_MODULE(ECLShowerCorrectorPureCsI)
40 
41 //-----------------------------------------------------------------
42 // Implementation
43 //-----------------------------------------------------------------
44 
46  m_leakageCorrectionPtr_bgx0("ecl_shower_corrector_leakage_corrections"),
47  m_leakageCorrectionPtr_thetaGeo_phase2bgx1("ECLLeakageCorrection_thetaGeometry_phase2"),
48  m_leakageCorrectionPtr_thetaGeo_phase3bgx1("ECLLeakageCorrection_thetaGeometry_phase3"),
49  m_leakageCorrectionPtr_phiGeo_phase2bgx1("ECLLeakageCorrection_phiGeometry_phase2"),
50  m_leakageCorrectionPtr_phiGeo_phase3bgx1("ECLLeakageCorrection_phiGeometry_phase3"),
51  m_leakageCorrectionPtr_thetaEn_phase2bgx1("ECLLeakageCorrection_thetaEnergy_phase2"),
52  m_leakageCorrectionPtr_thetaEn_phase3bgx1("ECLLeakageCorrection_thetaEnergy_phase3"),
53  m_leakageCorrectionPtr_phiEn_phase2bgx1("ECLLeakageCorrection_phiEnergy_phase2"),
54  m_leakageCorrectionPtr_phiEn_phase3bgx1("ECLLeakageCorrection_phiEnergy_phase3"),
55  m_eclShowers(eclShowerArrayName()),
56  m_eventLevelClusteringInfo(eventLevelClusteringInfoName())
57 {
58 
59  // Set description
60  setDescription("ECLShowerCorrectorModule: Corrects for MC truth to reconstruction shower and highest crystal energy differences");
61  setPropertyFlags(c_ParallelProcessingCertified);
62 
63 }
64 
66 {
67  ;
68 }
69 
71 {
72  B2DEBUG(28, "ECLShowerCorrectorModule::initialize()");
73 
74  // Register in datastore
75  m_eclShowers.registerInDataStore(eclShowerArrayName());
77 }
78 
80 {
81  // TODO: callback!
83 }
84 
86 {
87 
88  // Get the event background level.
89  const int bkgdcount = m_eventLevelClusteringInfo->getNECLCalDigitsOutOfTime();
90  double backgroundLevel = 0.0; // from out of time digit counting
91  if (m_fullBkgdCount > 0) {
92  backgroundLevel = static_cast<double>(bkgdcount) / m_fullBkgdCount;
93  }
94 
95  // Loop over all ECLShowers.
96  for (auto& eclShower : m_eclShowers) {
97 
98  // Only correct EM showers! Other showers keep the raw energy!
99  if (eclShower.getHypothesisId() == ECLShower::c_nPhotons) {
100 
101  const double energy = eclShower.getEnergy();
102  const double energyHighest = eclShower.getEnergyHighestCrystal();
103  const double theta = eclShower.getTheta() * TMath::RadToDeg();
104  const double phi = eclShower.getPhi() * TMath::RadToDeg();
105 
106  // Get the correction
107  double correctionFactor = 1.0;
108 
109  if (backgroundLevel < 0.1) correctionFactor = getLeakageCorrection(theta, phi, energy, backgroundLevel); // "Sumans corrections"
110  else correctionFactor = getLeakageCorrectionTemporary(theta, phi, energy, backgroundLevel); // "Elisas and Claudias corrections"
111 
112  B2DEBUG(28, "theta=" << theta << ", phi=" << phi << ", E=" << energy << ", BG=" << backgroundLevel << " --> correction factor=" <<
113  correctionFactor);
114 
115  if (correctionFactor < 1.e-5 or correctionFactor > 10.) {
116  B2ERROR("theta=" << theta << ", phi=" << phi << ", E=" << energy << " --> Correction factor=" << correctionFactor <<
117  " is very small/too large! Resetting to 1.0.");
118  correctionFactor = 1.0;
119  }
120 
121  const double correctedEnergy = energy * correctionFactor;
122  const double correctedEnergyHighest = energyHighest * correctionFactor;
123  B2DEBUG(28, "Correction factor=" << correctionFactor << ", correctedEnergy=" << correctedEnergy << ", correctedEnergyHighest=" <<
124  correctedEnergyHighest);
125 
126  // Set the correction
127  eclShower.setEnergy(correctedEnergy);
128  eclShower.setEnergyHighestCrystal(correctedEnergyHighest);
129 
130  } // end correction
131  } // end loop over all shower
132 
133 }
134 
136 {
137  ;
138 }
139 
141 {
142  ;
143 }
144 
146 {
147  //thetaGeometry phase2
149  m_leakage_bgx1_limits[0].resize(4);
154 
155  //thetaGeometry phase3
157  m_leakage_bgx1_limits[1].resize(4);
162  B2DEBUG(28, "Preparing leakage corrections, thetaGeometry phase3, theta and energy boundaries : theta_min=" <<
163  m_leakage_bgx1_limits[1][0] <<
164  ", thata_max=" << m_leakage_bgx1_limits[1][1]
165  << ", E_min=" << m_leakage_bgx1_limits[1][2] << " E_max= " << m_leakage_bgx1_limits[1][3]);
166 
167  //phiGeometry phase2
169  m_leakage_bgx1_limits[2].resize(4);
174 
175  //phiGeometry phase3
177  m_leakage_bgx1_limits[3].resize(4);
182  B2DEBUG(28, "Preparing leakage corrections, phiGeometry phase3, phi and energy boundaries : phi_min=" <<
183  m_leakage_bgx1_limits[3][0] <<
184  ", phi_max=" << m_leakage_bgx1_limits[3][1]
185  << ", E_min=" << m_leakage_bgx1_limits[3][2] << " E_max= " << m_leakage_bgx1_limits[3][3]);
186 
187  //thetaEnergy phase2
189  m_leakage_bgx1_limits[4].resize(4);
194 
195  //thetaEnergy phase3
197  m_leakage_bgx1_limits[5].resize(4);
202  B2DEBUG(28, "Preparing leakage corrections, thetaEnergy phase3, theta and energy boundaries : theta_min=" <<
203  m_leakage_bgx1_limits[5][0] <<
204  ", thata_max=" << m_leakage_bgx1_limits[5][1]
205  << ", E_min=" << m_leakage_bgx1_limits[5][2] << " E_max= " << m_leakage_bgx1_limits[5][3]);
206 
207 
208  //phiEnergy phase2
210  m_leakage_bgx1_limits[6].resize(4);
215 
216  //phiEnergy phase3
218  m_leakage_bgx1_limits[7].resize(4);
223  B2DEBUG(28, "Preparing leakage corrections, phiEnergy phase3, phi and energy boundaries : phi_min=" << m_leakage_bgx1_limits[7][0]
224  <<
225  ", phi_max=" << m_leakage_bgx1_limits[7][1]
226  << ", E_min=" << m_leakage_bgx1_limits[7][2] << " E_max= " << m_leakage_bgx1_limits[7][3]);
227 
228  // Prepare energy correction constants taken from the database to be used in an interpolation correction
229  // get all information from the payload
230  m_numOfBfBins = m_leakageCorrectionPtr_bgx0->getNumOfBfBins()[0];
231  m_numOfEnergyBins = m_leakageCorrectionPtr_bgx0->getNumOfEnergyBins()[0];
232  m_numOfPhiBins = m_leakageCorrectionPtr_bgx0->getNumOfPhiBins()[0];
233  m_numOfReg1ThetaBins = m_leakageCorrectionPtr_bgx0->getNumOfReg1ThetaBins()[0];
234  m_numOfReg2ThetaBins = m_leakageCorrectionPtr_bgx0->getNumOfReg2ThetaBins()[0];
235  m_numOfReg3ThetaBins = m_leakageCorrectionPtr_bgx0->getNumOfReg3ThetaBins()[0];
236 
237  // Resize the multidimensional vectors
241 
242  for (int iBfBin = 0; iBfBin < m_numOfBfBins; iBfBin++) {
246 
247  for (int iEBin = 0; iEBin < m_numOfEnergyBins; iEBin++) {
248  m_reg1CorrFactorArrays[iBfBin][iEBin].resize(m_numOfPhiBins);
249  m_reg2CorrFactorArrays[iBfBin][iEBin].resize(m_numOfPhiBins);
250  m_reg3CorrFactorArrays[iBfBin][iEBin].resize(m_numOfPhiBins);
251 
252  for (int iPhiBin = 0; iPhiBin < m_numOfPhiBins; iPhiBin++) {
253  m_reg1CorrFactorArrays[iBfBin][iEBin][iPhiBin].resize(m_numOfReg1ThetaBins);
254  m_reg2CorrFactorArrays[iBfBin][iEBin][iPhiBin].resize(m_numOfReg2ThetaBins);
255  m_reg3CorrFactorArrays[iBfBin][iEBin][iPhiBin].resize(m_numOfReg3ThetaBins);
256  }
257  }
258  }
259 
260  m_avgRecEn = m_leakageCorrectionPtr_bgx0->getAvgRecEn();
261 
262  // Fill the multidimensional vectors
263  m_correctionFactor = m_leakageCorrectionPtr_bgx0->getCorrectionFactor();
264  m_bgFractionBinNum = m_leakageCorrectionPtr_bgx0->getBgFractionBinNum();
265  m_regNum = m_leakageCorrectionPtr_bgx0->getRegNum();
266  m_phiBinNum = m_leakageCorrectionPtr_bgx0->getPhiBinNum();
267  m_thetaBinNum = m_leakageCorrectionPtr_bgx0->getThetaBinNum();
268  m_energyBinNum = m_leakageCorrectionPtr_bgx0->getEnergyBinNum();
269  m_correctionFactor = m_leakageCorrectionPtr_bgx0->getCorrectionFactor();
270 
271  for (unsigned i = 0; i < m_correctionFactor.size(); ++i) {
272  if (m_regNum[i] == 1) {
274  } else if (m_regNum[i] == 2) {
276  } else if (m_regNum[i] == 3) {
278  }
279  }
280 
281  m_phiPeriodicity = m_leakageCorrectionPtr_bgx0->getPhiPeriodicity()[0];
282  m_lReg1Theta = m_leakageCorrectionPtr_bgx0->getLReg1Theta()[0];
283  m_hReg1Theta = m_leakageCorrectionPtr_bgx0->getHReg1Theta()[0];
284  m_lReg2Theta = m_leakageCorrectionPtr_bgx0->getLReg2Theta()[0];
285  m_hReg2Theta = m_leakageCorrectionPtr_bgx0->getHReg2Theta()[0];
286  m_lReg3Theta = m_leakageCorrectionPtr_bgx0->getLReg3Theta()[0];
287  m_hReg3Theta = m_leakageCorrectionPtr_bgx0->getHReg3Theta()[0];
288 
289 }
290 
292  const double phi,
293  const double energy,
294  const double background)
295 {
296  // Corrections are available for Phase2BG15x1 and Phase3BG15x1.
297  int add_type = 0;
298  if (background > 0.75) add_type = 1;
299 
300  // m_leakage_bgx1_limits are ordered this way:
301  // 0 : theta_geo ph2, 1 : theta_geo ph3, 2 : phi_geo ph2, 3 : phi_geo ph3,
302  // 4 : theta_en ph2, 5 : theta_en ph3, 6 : phi_en ph2, 7 : phi_en ph3
303  const int type_theta_geo = 0 + add_type;
304  const int type_phi_geo = 2 + add_type;
305  const int type_theta_en = 4 + add_type ;
306  const int type_phi_en = 6 + add_type;
307 
308  B2DEBUG(28, "Index to select payloads: type_theta_geo = " << type_theta_geo << " , type_phi_geo = " << type_phi_geo
309  << " type_theta_en = " << type_theta_en << " , type_phi_en = " << type_phi_en);
310 
311  double theta_clip = theta;
312  double phi_clip = phi;
313  double energy_clip = energy;
314 
315  B2DEBUG(28, "Cluster info to compute leakage corrections: theta=" << theta_clip << ", phi=" << phi_clip << ", E=" << energy_clip);
316 
317  // check and clip boundaries since TGraph2D returns zero outside of them
318  if (theta_clip < m_leakage_bgx1_limits[type_theta_geo][0]) theta_clip = m_leakage_bgx1_limits[type_theta_geo][0] + 1e-4;
319  if (theta_clip > m_leakage_bgx1_limits[type_theta_geo][1]) theta_clip = m_leakage_bgx1_limits[type_theta_geo][1] - 1e-4;
320  if (energy_clip < m_leakage_bgx1_limits[type_theta_geo][2]) energy_clip = m_leakage_bgx1_limits[type_theta_geo][2] + 1e-5;
321  if (energy_clip > m_leakage_bgx1_limits[type_theta_geo][3]) energy_clip = m_leakage_bgx1_limits[type_theta_geo][3] - 1e-5;
322  const double corr_theta_geo = m_leakage_bgx1[type_theta_geo].Interpolate(theta_clip, energy_clip);
323  B2DEBUG(28, "After corr_theta_geo computation : theta=" << theta_clip << ", phi=" << phi_clip << ", E=" << energy_clip <<
324  " , corr_theta_geo : " << corr_theta_geo);
325  theta_clip = theta; energy_clip = energy; // re-initialize for the computation of the remaining corrections
326 
327  if (phi_clip < m_leakage_bgx1_limits[type_phi_geo][0]) phi_clip = m_leakage_bgx1_limits[type_phi_geo][0] + 1e-4;
328  if (phi_clip > m_leakage_bgx1_limits[type_phi_geo][1]) phi_clip = m_leakage_bgx1_limits[type_phi_geo][1] - 1e-4;
329  if (energy_clip < m_leakage_bgx1_limits[type_phi_geo][2]) energy_clip = m_leakage_bgx1_limits[type_phi_geo][2] + 1e-5;
330  if (energy_clip > m_leakage_bgx1_limits[type_phi_geo][3]) energy_clip = m_leakage_bgx1_limits[type_phi_geo][3] - 1e-5;
331  const double corr_phi_geo = m_leakage_bgx1[type_phi_geo].Interpolate(phi_clip, energy_clip);
332  B2DEBUG(28, "After corr_phi_geo computation : theta=" << theta_clip << ", phi=" << phi_clip << ", E=" << energy_clip <<
333  " , corr_phi_geo : " << corr_phi_geo);
334  phi_clip = phi; energy_clip = energy;
335 
336  if (theta_clip < m_leakage_bgx1_limits[type_theta_en][0]) theta_clip = m_leakage_bgx1_limits[type_theta_en][0] + 1e-4;
337  if (theta_clip > m_leakage_bgx1_limits[type_theta_en][1]) theta_clip = m_leakage_bgx1_limits[type_theta_en][1] - 1e-4;
338  if (energy_clip < m_leakage_bgx1_limits[type_theta_en][2]) energy_clip = m_leakage_bgx1_limits[type_theta_en][2] + 1e-5;
339  if (energy_clip > m_leakage_bgx1_limits[type_theta_en][3]) energy_clip = m_leakage_bgx1_limits[type_theta_en][3] - 1e-5;
340  const double corr_theta_en = m_leakage_bgx1[type_theta_en].Interpolate(theta_clip, energy_clip);
341  B2DEBUG(28, "After corr_theta_en computation : theta=" << theta_clip << ", phi=" << phi_clip << ", E=" << energy_clip <<
342  " , corr_theta_en : " << corr_theta_en);
343  theta_clip = theta; energy_clip = energy;
344 
345  if (phi_clip < m_leakage_bgx1_limits[type_phi_en][0]) phi_clip = m_leakage_bgx1_limits[type_phi_en][0] + 1e-4;
346  if (phi_clip > m_leakage_bgx1_limits[type_phi_en][1]) phi_clip = m_leakage_bgx1_limits[type_phi_en][1] - 1e-4;
347  if (energy_clip < m_leakage_bgx1_limits[type_phi_en][2]) energy_clip = m_leakage_bgx1_limits[type_phi_en][2] + 1e-5;
348  if (energy_clip > m_leakage_bgx1_limits[type_phi_en][3]) energy_clip = m_leakage_bgx1_limits[type_phi_en][3] - 1e-5;
349  const double corr_phi_en = m_leakage_bgx1[type_phi_en].Interpolate(phi_clip, energy_clip);
350  B2DEBUG(28, "After corr_phi_en computation : theta=" << theta_clip << ", phi=" << phi_clip << ", E=" << energy_clip <<
351  " , corr_phi_en : " << corr_phi_en);
352 
353 
354  // The final correction formula is:
355  // correctedEnergy = (rawEnergy + corr_theta_geo+corr_phi_geo ) *corr_theta_en*corr_phi_en
356  // to use "corr" as multiplicative factor in ECLShowerCorrectorModule::event,
357  // the following form should be returned
358  const double corr = (1 + (corr_theta_geo + corr_phi_geo) / energy_clip) * corr_theta_en * corr_phi_en;
359 
360  B2DEBUG(28, "Geometrical correction factors: corr_theta_geo=" << corr_theta_geo << ", corr_phi_geo=" << corr_phi_geo);
361  B2DEBUG(28, "Energy correction factors: corr_theta_en=" << corr_theta_en << ", corr_phi_en=" << corr_phi_en);
362  B2DEBUG(28, "Final correction=" << corr);
363 
364  return corr;
365 
366 }
367 
368 
369 
371  const double phi,
372  const double energy,
373  const double background) const
374 {
375  // Corrections are available for BGx0 and BGx1.0 (12th Campaign) only.
376  int bf = 0;
377  if (background > 0.1) bf = 1; // If users choose values other than BGx0 and BGx1 corrections will not be optimal.
378 
379  // Correction factor (multiplicative).
380  double result = 1.0;
381  double x0 = 0.0;
382  double x1 = 0.0;
383  double xd = 0.0;
384 
385  // Convert -180..180 to 0..360
386  const double phiMod = phi + 180;
387 
388  int x0Bin = 0;
389  int x1Bin = m_numOfPhiBins - 1;
390  int phiGap = int(phiMod / (360.0 / (m_phiPeriodicity * 1.0)));
391  double x = phiMod - phiGap * (360.0 / (m_phiPeriodicity * 1.0));
392 
393  const double phiBinWidth = 360.0 / (1.0 * m_numOfPhiBins * m_phiPeriodicity);
394 
395  if (x <= phiBinWidth / 2.0) {
396  x0 = 0.0;
397  x1 = phiBinWidth / 2.0;
398  x0Bin = 0;
399  x1Bin = m_numOfPhiBins - 1;
400  }
401 
402  for (int ii = 0; ii < (m_numOfPhiBins - 1); ii++) {
403  if (x >= (ii * phiBinWidth + phiBinWidth / 2.0) and x < ((ii + 1)*phiBinWidth + phiBinWidth / 2.0)) {
404  x0 = (ii * phiBinWidth + phiBinWidth / 2.0);
405  x1 = ((ii + 1) * phiBinWidth + phiBinWidth / 2.0);
406  x0Bin = ii;
407  x1Bin = ii + 1;
408  }
409  }
410 
411  if (x >= phiBinWidth * m_numOfPhiBins - phiBinWidth / 2.0) {
412  x0 = phiBinWidth * m_numOfPhiBins - phiBinWidth / 2.0;
413  x1 = phiBinWidth * m_numOfPhiBins;
414  x0Bin = m_numOfPhiBins - 1;
415  x1Bin = 0;
416  }
417 
418  B2DEBUG(28, "m_numOfPhiBins=" << m_numOfPhiBins << " x0Bin=" << x0Bin << " x1Bin=" << x1Bin);
419 
420  double y = theta;
421  double y0 = 0.;
422  double y1 = 0.;
423  double yd = 0.;
424  int y0Bin = 0;;
425  int y1Bin = 0;
426 
427  if (theta >= m_lReg1Theta and theta <= m_hReg1Theta) {
428  int reg1thetaBin = 0;
429  reg1thetaBin = int((theta - m_lReg1Theta) * m_numOfReg1ThetaBins / (m_hReg1Theta - m_lReg1Theta));
430  y = theta;
431  if (y < m_lReg1Theta + 0.5 * ((m_hReg1Theta - m_lReg1Theta) / (1.0 * m_numOfReg1ThetaBins))) {
432  y0 = m_lReg1Theta;
433  y1 = m_lReg1Theta + 0.5 * ((m_hReg1Theta - m_lReg1Theta) / (1.0 * m_numOfReg1ThetaBins));
434  y0Bin = 0;
435  y1Bin = 0;
436  } else if (y >= m_lReg1Theta + 0.5 * ((m_hReg1Theta - m_lReg1Theta) / (1.0 * m_numOfReg1ThetaBins))
437  and y < m_hReg1Theta - 0.5 * ((m_hReg1Theta - m_lReg1Theta) / (1.0 * m_numOfReg1ThetaBins))) {
438  y0 = m_lReg1Theta + 0.5 * ((m_hReg1Theta - m_lReg1Theta) / (1.0 * m_numOfReg1ThetaBins)) + reg1thetaBin * ((
440  y1 = m_lReg1Theta + 0.5 * ((m_hReg1Theta - m_lReg1Theta) / (1.0 * m_numOfReg1ThetaBins)) + (reg1thetaBin + 1) * ((
442  y0Bin = reg1thetaBin;
443  y1Bin = reg1thetaBin + 1;
444  if (y0Bin == m_numOfReg1ThetaBins - 1) y1Bin = m_numOfReg1ThetaBins - 1;
445  } else if (y >= m_hReg1Theta - 0.5 * ((m_hReg1Theta - m_lReg1Theta) / (1.0 * m_numOfReg1ThetaBins))) {
446  y0 = m_hReg1Theta - 0.5 * ((m_hReg1Theta - m_lReg1Theta) / (1.0 * m_numOfReg1ThetaBins));
447  y1 = m_hReg1Theta;
448  y0Bin = m_numOfReg1ThetaBins - 1;
449  y1Bin = m_numOfReg1ThetaBins - 1;
450  } else {
451  y1 = m_hReg1Theta;
452  y0Bin = m_numOfReg1ThetaBins - 1;
453  y1Bin = m_numOfReg1ThetaBins - 1;
454  }
455  } //End of region 1
456  else if (theta >= m_lReg2Theta && theta <= m_hReg2Theta) {
457  int reg2thetaBin = 0;
458  reg2thetaBin = int((theta - m_lReg2Theta) * m_numOfReg2ThetaBins / (m_hReg2Theta - m_lReg2Theta));
459  y = theta;
460  if (y < m_lReg2Theta + 0.5 * ((m_hReg2Theta - m_lReg2Theta) / (1.0 * m_numOfReg2ThetaBins))) {
461  y0 = m_lReg2Theta;
462  y1 = m_lReg2Theta + 0.5 * ((m_hReg2Theta - m_lReg2Theta) / (1.0 * m_numOfReg2ThetaBins));
463  y0Bin = 0;
464  y1Bin = 0;
465  } else if (y >= m_lReg2Theta + 0.5 * ((m_hReg2Theta - m_lReg2Theta) / (1.0 * m_numOfReg2ThetaBins))
466  and y < m_hReg2Theta - 0.5 * ((m_hReg2Theta - m_lReg2Theta) / (1.0 * m_numOfReg2ThetaBins))) {
467  y0 = m_lReg2Theta + 0.5 * ((m_hReg2Theta - m_lReg2Theta) / (1.0 * m_numOfReg2ThetaBins)) + reg2thetaBin * ((
469  y1 = m_lReg2Theta + 0.5 * ((m_hReg2Theta - m_lReg2Theta) / (1.0 * m_numOfReg2ThetaBins)) + (reg2thetaBin + 1) * ((
471  y0Bin = reg2thetaBin;
472  y1Bin = reg2thetaBin + 1;
473  if (y0Bin == m_numOfReg2ThetaBins - 1) y1Bin = m_numOfReg2ThetaBins - 1;
474  } else if (y >= m_hReg2Theta - 0.5 * ((m_hReg2Theta - m_lReg2Theta) / (1.0 * m_numOfReg2ThetaBins))) {
475  y0 = m_hReg2Theta - 0.5 * ((m_hReg2Theta - m_lReg2Theta) / (1.0 * m_numOfReg2ThetaBins)); y1 = m_hReg2Theta;
476  y0Bin = m_numOfReg2ThetaBins - 1; y1Bin = m_numOfReg2ThetaBins - 1;
477  } else {y1 = m_hReg2Theta; y0Bin = m_numOfReg2ThetaBins - 1; y1Bin = m_numOfReg2ThetaBins - 1;}
478  } //End of region 2
479  else if (theta >= m_lReg3Theta && theta <= m_hReg3Theta) {
480  int reg3thetaBin = 0;
481  reg3thetaBin = int((theta - m_lReg3Theta) * m_numOfReg3ThetaBins / (m_hReg3Theta - m_lReg3Theta));
482  y = theta;
483  if (y < m_lReg3Theta + 0.5 * ((m_hReg3Theta - m_lReg3Theta) / (1.0 * m_numOfReg3ThetaBins))) {
484  y0 = m_lReg3Theta;
485  y1 = m_lReg3Theta + 0.5 * ((m_hReg3Theta - m_lReg3Theta) / (1.0 * m_numOfReg3ThetaBins));
486  y0Bin = 0;
487  y1Bin = 0;
488  } else if (y >= m_lReg3Theta + 0.5 * ((m_hReg3Theta - m_lReg3Theta) / (1.0 * m_numOfReg3ThetaBins))
489  and y < m_hReg3Theta - 0.5 * ((m_hReg3Theta - m_lReg3Theta) / (1.0 * m_numOfReg3ThetaBins))) {
490  y0 = m_lReg3Theta + 0.5 * ((m_hReg3Theta - m_lReg3Theta) / (1.0 * m_numOfReg3ThetaBins)) + reg3thetaBin * ((
492  y1 = m_lReg3Theta + 0.5 * ((m_hReg3Theta - m_lReg3Theta) / (1.0 * m_numOfReg3ThetaBins)) + (reg3thetaBin + 1) * ((
494  y0Bin = reg3thetaBin;
495  y1Bin = reg3thetaBin + 1;
496  if (y0Bin == m_numOfReg3ThetaBins - 1) y1Bin = m_numOfReg3ThetaBins - 1;
497  } else if (y >= m_hReg3Theta - 0.5 * ((m_hReg3Theta - m_lReg3Theta) / (1.0 * m_numOfReg3ThetaBins))) {
498  y0 = m_hReg3Theta - 0.5 * ((m_hReg3Theta - m_lReg3Theta) / (1.0 * m_numOfReg3ThetaBins));
499  y1 = m_hReg3Theta;
500  y0Bin = m_numOfReg3ThetaBins - 1;
501  y1Bin = m_numOfReg3ThetaBins - 1;
502  } else {
503  y1 = m_hReg3Theta;
504  y0Bin = m_numOfReg3ThetaBins - 1;
505  y1Bin = m_numOfReg3ThetaBins - 1;
506  }
507  } //End of region 3
508  else {
509  B2DEBUG(28, "No valid theta region found (theta=" << theta << " deg), return correction factor of 1.0.");
510  return 1.0; // We have not found a theta region that is valid, return no correction at all (i.e. 1.0)
511  }
512 
513  B2DEBUG(28, "y0Bin=" << y0Bin << " y1Bin=" << y1Bin);
514 
515  //int energyBin = 0;
516  double z = energy;
517  double z0 = 0.0;
518  double z1 = 0.0;
519  double zd = 0.0;
520  int z0Bin = 0;
521  int z1Bin = 1;
522 
523  if (z >= 0.0 and z < m_avgRecEn[0]) {
524  z0 = 0.01;
525  z1 = m_avgRecEn[0];
526  z0Bin = 0;
527  z1Bin = 0;
528  } else if (z >= m_avgRecEn[0] and z < m_avgRecEn[m_numOfEnergyBins - 1]) {
529  for (int ii = 0; ii < (m_numOfEnergyBins - 1); ii++) {
530  if (z >= m_avgRecEn[ii] and z < m_avgRecEn[ii + 1]) {
531  z0 = m_avgRecEn[ii];
532  z1 = m_avgRecEn[ii + 1];
533  z0Bin = ii;
534  z1Bin = ii + 1;
535  }
536  }
537  } else if (z >= m_avgRecEn[m_numOfEnergyBins - 1] and z < 9.0) {
538  z0 = m_avgRecEn[m_numOfEnergyBins - 1];
539  z1 = 9.0;
540  z0Bin = (m_numOfEnergyBins - 1);
541  z1Bin = (m_numOfEnergyBins - 1);
542  } else {
543  z0 = m_avgRecEn[m_numOfEnergyBins - 1];
544  z1 = 9.0;
545  z0Bin = (m_numOfEnergyBins - 1);
546  z1Bin = (m_numOfEnergyBins - 1);
547  }
548 
549  B2DEBUG(28, "z0Bin=" << z0Bin << " z1Bin=" << z1Bin);
550 
551  xd = (x - x0) / (x1 - x0);
552  yd = (y - y0) / (y1 - y0);
553  zd = (z - z0) / (z1 - z0);
554 
555  if (theta >= m_lReg2Theta && theta <= m_hReg2Theta) {
556  double c00, c01, c10, c11, c0, c1;
557  c00 = m_reg2CorrFactorArrays[bf][z0Bin][x0Bin][y0Bin] * (1 - xd) + m_reg2CorrFactorArrays[bf][z0Bin][x1Bin][y0Bin] * xd;
558  c01 = m_reg2CorrFactorArrays[bf][z1Bin][x0Bin][y0Bin] * (1 - xd) + m_reg2CorrFactorArrays[bf][z1Bin][x1Bin][y0Bin] * xd;
559  c10 = m_reg2CorrFactorArrays[bf][z0Bin][x0Bin][y1Bin] * (1 - xd) + m_reg2CorrFactorArrays[bf][z0Bin][x1Bin][y1Bin] * xd;
560  c11 = m_reg2CorrFactorArrays[bf][z1Bin][x0Bin][y1Bin] * (1 - xd) + m_reg2CorrFactorArrays[bf][z1Bin][x1Bin][y1Bin] * xd;
561 
562  c0 = c00 * (1 - yd) + c10 * yd;
563  c1 = c01 * (1 - yd) + c11 * yd;
564 
565  result = c0 * (1 - zd) + c1 * zd;
566  } else if (theta >= m_lReg1Theta && theta <= m_hReg1Theta) {
567  double c00, c01, c10, c11, c0, c1;
568  c00 = m_reg1CorrFactorArrays[bf][z0Bin][x0Bin][y0Bin] * (1 - xd) + m_reg1CorrFactorArrays[bf][z0Bin][x1Bin][y0Bin] * xd;
569  c01 = m_reg1CorrFactorArrays[bf][z1Bin][x0Bin][y0Bin] * (1 - xd) + m_reg1CorrFactorArrays[bf][z1Bin][x1Bin][y0Bin] * xd;
570  c10 = m_reg1CorrFactorArrays[bf][z0Bin][x0Bin][y1Bin] * (1 - xd) + m_reg1CorrFactorArrays[bf][z0Bin][x1Bin][y1Bin] * xd;
571  c11 = m_reg1CorrFactorArrays[bf][z1Bin][x0Bin][y1Bin] * (1 - xd) + m_reg1CorrFactorArrays[bf][z1Bin][x1Bin][y1Bin] * xd;
572 
573  c0 = c00 * (1 - yd) + c10 * yd;
574  c1 = c01 * (1 - yd) + c11 * yd;
575 
576  result = c0 * (1 - zd) + c1 * zd;
577  } else if (theta >= m_lReg3Theta && theta <= m_hReg3Theta) {
578  float c00, c01, c10, c11, c0, c1;
579  c00 = m_reg3CorrFactorArrays[bf][z0Bin][x0Bin][y0Bin] * (1 - xd) + m_reg3CorrFactorArrays[bf][z0Bin][x1Bin][y0Bin] * xd;
580  c01 = m_reg3CorrFactorArrays[bf][z1Bin][x0Bin][y0Bin] * (1 - xd) + m_reg3CorrFactorArrays[bf][z1Bin][x1Bin][y0Bin] * xd;
581  c10 = m_reg3CorrFactorArrays[bf][z0Bin][x0Bin][y1Bin] * (1 - xd) + m_reg3CorrFactorArrays[bf][z0Bin][x1Bin][y1Bin] * xd;
582  c11 = m_reg3CorrFactorArrays[bf][z1Bin][x0Bin][y1Bin] * (1 - xd) + m_reg3CorrFactorArrays[bf][z1Bin][x1Bin][y1Bin] * xd;
583 
584  c0 = c00 * (1 - yd) + c10 * yd;
585  c1 = c01 * (1 - yd) + c11 * yd;
586 
587  result = c0 * (1 - zd) + c1 * zd;
588  } else {
589  result = (m_reg2CorrFactorArrays[bf][z0Bin][0][0] + m_reg2CorrFactorArrays[bf][z1Bin][0][0]) / 2.0;
590  }
591 
592  return result;
593 }
Belle2::ECLShowerCorrectorModule::m_avgRecEn
std::vector< float > m_avgRecEn
averages of the energy bins
Definition: ECLShowerCorrectorModule.h:136
Belle2::ECLShowerCorrectorModule::m_eventLevelClusteringInfo
StoreObjPtr< EventLevelClusteringInfo > m_eventLevelClusteringInfo
Store object pointer: EventLevelClusteringInfo.
Definition: ECLShowerCorrectorModule.h:154
Belle2::ECLShowerCorrectorModule::beginRun
virtual void beginRun() override
Begin run.
Definition: ECLShowerCorrectorModule.cc:79
Belle2::ECLShowerCorrectorModule::m_reg2CorrFactorArrays
std::vector< std::vector< std::vector< std::vector< float > > > > m_reg2CorrFactorArrays
region 2 corrections
Definition: ECLShowerCorrectorModule.h:147
Belle2::ECLShowerCorrectorModule::m_phiPeriodicity
int m_phiPeriodicity
repeating pattern in phi direction, for barrel it is 72
Definition: ECLShowerCorrectorModule.h:127
Belle2::ECLShowerCorrectorModule::getLeakageCorrectionTemporary
double getLeakageCorrectionTemporary(const double theta, const double phi, const double energy, const double background)
Get correction for BGx1 (temporary)
Definition: ECLShowerCorrectorModule.cc:291
Belle2::ECLShowerCorrectorModule::m_hReg3Theta
float m_hReg3Theta
upper boundary of the region 3 theta
Definition: ECLShowerCorrectorModule.h:133
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::ECLShowerCorrectorModule::endRun
virtual void endRun() override
End run.
Definition: ECLShowerCorrectorModule.cc:135
Belle2::ECLShowerCorrectorModule::m_numOfBfBins
int m_numOfBfBins
number of background fraction bins; currently only two
Definition: ECLShowerCorrectorModule.h:121
Belle2::ECLShowerCorrectorModule::m_numOfReg1ThetaBins
int m_numOfReg1ThetaBins
number of region 1 theta bins
Definition: ECLShowerCorrectorModule.h:124
Belle2::ECLShowerCorrectorModule::m_reg1CorrFactorArrays
std::vector< std::vector< std::vector< std::vector< float > > > > m_reg1CorrFactorArrays
region 1 corrections
Definition: ECLShowerCorrectorModule.h:146
Belle2::ECLShowerCorrectorModule::m_lReg2Theta
float m_lReg2Theta
lower boundary of the region 2 theta
Definition: ECLShowerCorrectorModule.h:130
Belle2::ECLShowerCorrectorModule::m_leakageCorrectionPtr_phiGeo_phase2bgx1
DBObjPtr< ECLShowerEnergyCorrectionTemporary > m_leakageCorrectionPtr_phiGeo_phase2bgx1
Leakage corrections from DB for Phase2 BG15x1.0, geometry correction as a function of phi.
Definition: ECLShowerCorrectorModule.h:98
Belle2::ECLShowerCorrectorModule::m_leakageCorrectionPtr_bgx0
DBObjPtr< ECLShowerCorrectorLeakageCorrection > m_leakageCorrectionPtr_bgx0
Leakage corrections from DB for BGx0.
Definition: ECLShowerCorrectorModule.h:91
Belle2::ECLShowerCorrectorModule::m_leakage_bgx1
TGraph2D m_leakage_bgx1[8]
the leakage in BGx1
Definition: ECLShowerCorrectorModule.h:115
Belle2::ECLShowerCorrectorModule::m_regNum
std::vector< int > m_regNum
region bin
Definition: ECLShowerCorrectorModule.h:140
Belle2::ECLShowerCorrectorModule::event
virtual void event() override
Event.
Definition: ECLShowerCorrectorModule.cc:85
Belle2::ECLShowerCorrectorModule::terminate
virtual void terminate() override
Terminate.
Definition: ECLShowerCorrectorModule.cc:140
Belle2::ECLShowerCorrectorModule::m_phiBinNum
std::vector< int > m_phiBinNum
phi bin
Definition: ECLShowerCorrectorModule.h:141
Belle2::ECLShowerCorrectorModule::m_numOfReg3ThetaBins
int m_numOfReg3ThetaBins
number of region 3 theta bins
Definition: ECLShowerCorrectorModule.h:126
Belle2::ECLShowerCorrectorModule::~ECLShowerCorrectorModule
~ECLShowerCorrectorModule()
Destructor.
Definition: ECLShowerCorrectorModule.cc:65
Belle2::ECLShowerCorrectorModule::m_reg3CorrFactorArrays
std::vector< std::vector< std::vector< std::vector< float > > > > m_reg3CorrFactorArrays
region 3 corrections
Definition: ECLShowerCorrectorModule.h:148
Belle2::ECLShowerCorrectorModule::m_hReg1Theta
float m_hReg1Theta
upper boundary of the region 1 theta
Definition: ECLShowerCorrectorModule.h:129
Belle2::ECLShowerCorrectorModule::m_lReg1Theta
float m_lReg1Theta
lower boundary of the region 1 theta
Definition: ECLShowerCorrectorModule.h:128
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::ECLShowerCorrectorModule::m_bgFractionBinNum
std::vector< int > m_bgFractionBinNum
BG fraction bin.
Definition: ECLShowerCorrectorModule.h:139
Belle2::ECLShowerCorrectorModule::m_energyBinNum
std::vector< int > m_energyBinNum
energu bin
Definition: ECLShowerCorrectorModule.h:143
Belle2::ECLShowerCorrectorModule::eventLevelClusteringInfoName
virtual const char * eventLevelClusteringInfoName() const
Name to be used for default option: EventLevelClusteringInfo.
Definition: ECLShowerCorrectorModule.h:164
Belle2::ECLShowerCorrectorModule::eclShowerArrayName
virtual const char * eclShowerArrayName() const
We need names for the data objects to differentiate between PureCsI and default.
Definition: ECLShowerCorrectorModule.h:160
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::ECLShowerCorrectorModule::m_leakageCorrectionPtr_thetaGeo_phase3bgx1
DBObjPtr< ECLShowerEnergyCorrectionTemporary > m_leakageCorrectionPtr_thetaGeo_phase3bgx1
Leakage corrections from DB for Phase3 BG15x1.0, geometry correction as a function of theta.
Definition: ECLShowerCorrectorModule.h:96
Belle2::ECLShowerCorrectorModule
Class to perform the shower correction.
Definition: ECLShowerCorrectorModule.h:57
Belle2::ECLShowerCorrectorModule::m_leakageCorrectionPtr_phiEn_phase2bgx1
DBObjPtr< ECLShowerEnergyCorrectionTemporary > m_leakageCorrectionPtr_phiEn_phase2bgx1
Leakage corrections from DB for Phase2 BG15x1.0, energy correction as a function of phi.
Definition: ECLShowerCorrectorModule.h:106
Belle2::ECLShowerCorrectorModule::initialize
virtual void initialize() override
Initialize.
Definition: ECLShowerCorrectorModule.cc:70
Belle2::ECLShower::c_nPhotons
@ c_nPhotons
CR is split into n photons (N1)
Definition: ECLShower.h:54
Belle2::ECLShowerCorrectorModule::m_eclShowers
StoreArray< ECLShower > m_eclShowers
Store array: ECLShower.
Definition: ECLShowerCorrectorModule.h:151
Belle2::ECLShowerCorrectorModule::m_leakageCorrectionPtr_phiEn_phase3bgx1
DBObjPtr< ECLShowerEnergyCorrectionTemporary > m_leakageCorrectionPtr_phiEn_phase3bgx1
Leakage corrections from DB for Phase3 BG15x1.0, energy correction as a function of phi.
Definition: ECLShowerCorrectorModule.h:108
Belle2::ECLShowerCorrectorModule::m_lReg3Theta
float m_lReg3Theta
lower boundary of the region 3 theta
Definition: ECLShowerCorrectorModule.h:132
Belle2::ECLShowerCorrectorModule::prepareLeakageCorrections
void prepareLeakageCorrections()
Prepare correction.
Definition: ECLShowerCorrectorModule.cc:145
Belle2::ECLShowerCorrectorModule::m_leakageCorrectionPtr_thetaEn_phase3bgx1
DBObjPtr< ECLShowerEnergyCorrectionTemporary > m_leakageCorrectionPtr_thetaEn_phase3bgx1
Leakage corrections from DB for Phase3 BG15x1.0, energy correction as a function of theta.
Definition: ECLShowerCorrectorModule.h:104
Belle2::ECLShowerCorrectorModule::m_thetaBinNum
std::vector< int > m_thetaBinNum
theta bin
Definition: ECLShowerCorrectorModule.h:142
Belle2::ECLShowerCorrectorModule::m_leakageCorrectionPtr_thetaEn_phase2bgx1
DBObjPtr< ECLShowerEnergyCorrectionTemporary > m_leakageCorrectionPtr_thetaEn_phase2bgx1
Leakage corrections from DB for Phase2 BG15x1.0, energy correction as a function of theta.
Definition: ECLShowerCorrectorModule.h:102
Belle2::ECLShowerCorrectorModule::m_numOfEnergyBins
int m_numOfEnergyBins
number of energy bins
Definition: ECLShowerCorrectorModule.h:122
Belle2::ECLShowerCorrectorModule::getLeakageCorrection
double getLeakageCorrection(const double theta, const double phi, const double energy, const double background) const
Get correction for BGx0.
Definition: ECLShowerCorrectorModule.cc:370
Belle2::ECLShowerCorrectorModule::m_leakageCorrectionPtr_thetaGeo_phase2bgx1
DBObjPtr< ECLShowerEnergyCorrectionTemporary > m_leakageCorrectionPtr_thetaGeo_phase2bgx1
Leakage corrections from DB for Phase2 BG15x1.0, geometry correction as a function of theta.
Definition: ECLShowerCorrectorModule.h:94
Belle2::ECLShowerCorrectorModule::m_correctionFactor
std::vector< float > m_correctionFactor
correction value
Definition: ECLShowerCorrectorModule.h:144
Belle2::ECLShowerCorrectorModule::m_hReg2Theta
float m_hReg2Theta
upper boundary of the region 2 theta
Definition: ECLShowerCorrectorModule.h:131
Belle2::ECLShowerCorrectorModule::m_leakageCorrectionPtr_phiGeo_phase3bgx1
DBObjPtr< ECLShowerEnergyCorrectionTemporary > m_leakageCorrectionPtr_phiGeo_phase3bgx1
Leakage corrections from DB for Phase3 BG15x1.0, geometry correction as a function of phi.
Definition: ECLShowerCorrectorModule.h:100
Belle2::ECLShowerCorrectorModule::m_numOfPhiBins
int m_numOfPhiBins
number of phi bins
Definition: ECLShowerCorrectorModule.h:123
Belle2::ECLShowerCorrectorModule::m_leakage_bgx1_limits
std::vector< double > m_leakage_bgx1_limits[8]
limits for the leakage in BGx1
Definition: ECLShowerCorrectorModule.h:117
Belle2::ECLShowerCorrectorModule::m_fullBkgdCount
const double m_fullBkgdCount
Nominal Background at BGx1.0 (MC12)
Definition: ECLShowerCorrectorModule.h:111
Belle2::ECLShowerCorrectorModule::m_numOfReg2ThetaBins
int m_numOfReg2ThetaBins
number of region 2 theta bins
Definition: ECLShowerCorrectorModule.h:125