Belle II Software  release-05-02-19
ECLCovarianceMatrixModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2016 - Belle II Collaboration *
4  * *
5  * This module calculates the covariance matrix for photon showers. *
6  * The matrix depends on the shower region (FWD, Bartel, BWD) *
7  * *
8  * Author: The Belle II Collaboration *
9  * Contributors: Torben Ferber (ferber@physics.ubc.ca) *
10  * *
11  * This software is provided "as is" without any warranty. *
12  **************************************************************************/
13 
14 // THIS MODULE
15 #include <ecl/modules/eclCovarianceMatrix/ECLCovarianceMatrixModule.h>
16 
17 // ROOT
18 #include "TMath.h"
19 
20 // MDST
21 #include <mdst/dataobjects/EventLevelClusteringInfo.h>
22 
23 // ECL
24 #include <ecl/dataobjects/ECLShower.h>
25 
26 // NAMESPACES
27 using namespace Belle2;
28 
29 //-----------------------------------------------------------------
30 // Register the Module
31 //-----------------------------------------------------------------
32 REG_MODULE(ECLCovarianceMatrix)
33 REG_MODULE(ECLCovarianceMatrixPureCsI)
34 
35 //-----------------------------------------------------------------
36 // Implementation
37 //-----------------------------------------------------------------
39  m_eclShowers(eclShowerArrayName()),
40  m_eventLevelClusteringInfo(eventLevelClusteringInfoName())
41 {
42  // Set description
43  setDescription("ECLCovarianceMatrix: Sets the ECL photon shower covariance matrix.");
44  setPropertyFlags(c_ParallelProcessingCertified);
45 
46 }
47 
49 {
50 }
51 
53 {
54  // Register in datastore
55  m_eclShowers.registerInDataStore(eclShowerArrayName());
57 }
58 
60 {
61  // TODO: callback
62  ;
63 }
64 
66 {
67 
68  // Get the event background level
69  const int bkgdcount = m_eventLevelClusteringInfo->getNECLCalDigitsOutOfTime();
70  double background = 0.0; // from out of time digit counting
71  if (m_fullBkgdCount > 0) {
72  background = static_cast<double>(bkgdcount) / m_fullBkgdCount;
73  }
74 
75  // loop over all ECLShowers
76  for (auto& eclShower : m_eclShowers) {
77 
78  // Only calculate for photon showers
79  if (eclShower.getHypothesisId() == ECLShower::c_nPhotons) {
80 
81  const double energy = eclShower.getEnergy();
82 
83  // 1/energy
84  double invEnergy = 0;
85  if (energy > 0.) invEnergy = 1. / energy;
86 
87  //1/energy^0.5
88  double invRoot2Energy = 0;
89  if (energy > 0.) invRoot2Energy = 1. / TMath::Power(energy, 0.5);
90 
91  //1/energy^0.25
92  double invRoot4Energy = 0;
93  if (energy > 0.) invRoot4Energy = 1. / TMath::Power(energy, 0.25);
94 
95  int detregion = eclShower.getDetectorRegion(); // FWD: 1, Barrel: 2, BWD: 3
96  if (detregion == 11 or detregion == 13) detregion = 2;
97 
98  double sigmaEnergy = 0.;
99  double sigmaTheta = 0.;
100  double sigmaPhi = 0.;
101 
102  // three background levels, three detector regions, energy, theta and phi (needs to be revisited soon!)
103  if (detregion == 1 and background < 0.05) {
104  if (energy <= 0.022) {sigmaEnergy = energy * (0.0449188); }
105  else {sigmaEnergy = energy * (-0.0912379 * invEnergy + 1.91849 * invRoot2Energy + -2.82169 * invRoot4Energy + 3.03119) / 100.;}
106  } else if (detregion == 1 and background <= 0.3) {
107  if (energy <= 0.022) {sigmaEnergy = energy * (0.100769); }
108  else {sigmaEnergy = energy * (-0.0725928 * invEnergy + 3.4685 * invRoot2Energy + -5.44127 * invRoot4Energy + 4.12045) / 100.;}
109  } else if (detregion == 1 and background > 0.3) {
110  if (energy <= 0.022) {sigmaEnergy = energy * (0.148462); }
111  else {sigmaEnergy = energy * (-0.239931 * invEnergy + 6.94958 * invRoot2Energy + -10.4085 * invRoot4Energy + 5.92412) / 100.;}
112  }
113 
114  if (detregion == 1 and background < 0.05) {
115  if (energy <= 0.0289609) { sigmaTheta = 1e-3 * (6.4008);}
116  else {sigmaTheta = 1e-3 * (-0.110397 * invEnergy + 0.753603 * invRoot2Energy + 2.63652 * invRoot4Energy + -0.606703);}
117  } else if (detregion == 1 and background <= 0.3) {
118  if (energy <= 0.0274562) { sigmaTheta = 1e-3 * (7.39868);}
119  else {sigmaTheta = 1e-3 * (-0.207278 * invEnergy + 2.68616 * invRoot2Energy + -0.905487 * invRoot4Energy + 0.961485);}
120  } else if (detregion == 1 and background > 0.3) {
121  if (energy <= 0.022) { sigmaTheta = 1e-3 * (8.83505);}
122  else {sigmaTheta = 1e-3 * (-0.160921 * invEnergy + 2.35311 * invRoot2Energy + -0.107975 * invRoot4Energy + 0.565367);}
123  }
124 
125  if (detregion == 1 and background < 0.05) {
126  if (energy <= 0.0391279) {sigmaPhi = 1e-3 * (15.4147);}
127  else {sigmaPhi = 1e-3 * (-0.35934 * invEnergy + 2.01807 * invRoot2Energy + 7.26313 * invRoot4Energy + -1.93438);}
128  } else if (detregion == 1 and background <= 0.3) {
129  if (energy <= 0.022) {sigmaPhi = 1e-3 * (18.1096);}
130  else {sigmaPhi = 1e-3 * (-0.178258 * invEnergy + 1.39155 * invRoot2Energy + 6.97462 * invRoot4Energy + -1.2795);}
131  } else if (detregion == 1 and background > 0.3) {
132  if (energy <= 0.0224565) {sigmaPhi = 1e-3 * (20.1605);}
133  else {sigmaPhi = 1e-3 * (-0.344617 * invEnergy + 3.80823 * invRoot2Energy + 4.08729 * invRoot4Energy + -0.464714);}
134  }
135 
136  if (detregion == 2 and background < 0.05) {
137  if (energy <= 0.022) {sigmaEnergy = energy * (0.0466552); }
138  else {sigmaEnergy = energy * (-0.0473491 * invEnergy + 1.02015 * invRoot2Energy + -0.705164 * invRoot4Energy + 1.77086) / 100.;}
139  } else if (detregion == 2 and background <= 0.3) {
140  if (energy <= 0.022) {sigmaEnergy = energy * (0.112299); }
141  else {sigmaEnergy = energy * (0.0469588 * invEnergy + 1.81959 * invRoot2Energy + -2.13666 * invRoot4Energy + 2.37568) / 100.;}
142  } else if (detregion == 2 and background > 0.3) {
143  if (energy <= 0.022) {sigmaEnergy = energy * (0.164881); }
144  else {sigmaEnergy = energy * (-0.267069 * invEnergy + 7.6176 * invRoot2Energy + -11.0341 * invRoot4Energy + 5.9203) / 100.;}
145  }
146 
147  if (detregion == 2 and background < 0.05) {
148  if (energy <= 0.0309077) { sigmaTheta = 1e-3 * (9.31098);}
149  else {sigmaTheta = 1e-3 * (-0.221073 * invEnergy + 2.19876 * invRoot2Energy + 1.50833 * invRoot4Energy + 0.359609);}
150  } else if (detregion == 2 and background <= 0.3) {
151  if (energy <= 0.022) { sigmaTheta = 1e-3 * (11.6166);}
152  else {sigmaTheta = 1e-3 * (-0.126603 * invEnergy + 1.81898 * invRoot2Energy + 1.85132 * invRoot4Energy + 0.300728);}
153  } else if (detregion == 2 and background > 0.3) {
154  if (energy <= 0.022) { sigmaTheta = 1e-3 * (13.159);}
155  else {sigmaTheta = 1e-3 * (-0.316561 * invEnergy + 4.94686 * invRoot2Energy + -3.12199 * invRoot4Energy + 2.30275);}
156  }
157 
158  if (detregion == 2 and background < 0.05) {
159  if (energy <= 0.0381178) {sigmaPhi = 1e-3 * (10.7661);}
160  else {sigmaPhi = 1e-3 * (-0.21842 * invEnergy + 0.648976 * invRoot2Energy + 7.1901 * invRoot4Energy + -3.10025);}
161  } else if (detregion == 2 and background <= 0.3) {
162  if (energy <= 0.022237) {sigmaPhi = 1e-3 * (12.6804);}
163  else {sigmaPhi = 1e-3 * (-0.141507 * invEnergy + 0.488293 * invRoot2Energy + 7.30055 * invRoot4Energy + -3.13591);}
164  } else if (detregion == 2 and background > 0.3) {
165  if (energy <= 0.024905) {sigmaPhi = 1e-3 * (14.8726);}
166  else {sigmaPhi = 1e-3 * (-0.323752 * invEnergy + 3.60933 * invRoot2Energy + 2.48526 * invRoot4Energy + -1.25493);}
167  }
168 
169  if (detregion == 3 and background < 0.05) {
170  if (energy <= 0.022) {sigmaEnergy = energy * (0.0466662); }
171  else {sigmaEnergy = energy * (-0.106729 * invEnergy + 2.38163 * invRoot2Energy + -4.11299 * invRoot4Energy + 4.14056) / 100.;}
172  } else if (detregion == 3 and background <= 0.3) {
173  if (energy <= 0.022) {sigmaEnergy = energy * (0.128851); }
174  else {sigmaEnergy = energy * (-0.204362 * invEnergy + 6.43 * invRoot2Energy + -10.8673 * invRoot4Energy + 7.04059) / 100.;}
175  } else if (detregion == 3 and background > 0.3) {
176  if (energy <= 0.022) {sigmaEnergy = energy * (0.255153); }
177  else {sigmaEnergy = energy * (-0.316245 * invEnergy + 11.5211 * invRoot2Energy + -18.3626 * invRoot4Energy + 9.89422) / 100.;}
178  }
179 
180  if (detregion == 3 and background < 0.05) {
181  if (energy <= 0.0318079) { sigmaTheta = 1e-3 * (10.4446);}
182  else {sigmaTheta = 1e-3 * (-0.295827 * invEnergy + 3.37748 * invRoot2Energy + -0.284446 * invRoot4Energy + 1.48098);}
183  } else if (detregion == 3 and background <= 0.3) {
184  if (energy <= 0.022) { sigmaTheta = 1e-3 * (13.828);}
185  else {sigmaTheta = 1e-3 * (-0.252479 * invEnergy + 4.18833 * invRoot2Energy + -2.04822 * invRoot4Energy + 2.38486);}
186  } else if (detregion == 3 and background > 0.3) {
187  if (energy <= 0.022) { sigmaTheta = 1e-3 * (17.0354);}
188  else {sigmaTheta = 1e-3 * (-0.420215 * invEnergy + 7.08904 * invRoot2Energy + -5.88269 * invRoot4Energy + 3.61643);}
189  }
190 
191  if (detregion == 3 and background < 0.05) {
192  if (energy <= 0.022) {sigmaPhi = 1e-3 * (15.129);}
193  else {sigmaPhi = 1e-3 * (-0.0534441 * invEnergy + -1.72466 * invRoot2Energy + 12.8625 * invRoot4Energy + -4.21203);}
194  } else if (detregion == 3 and background <= 0.3) {
195  if (energy <= 0.022) {sigmaPhi = 1e-3 * (20.1408);}
196  else {sigmaPhi = 1e-3 * (0.0869293 * invEnergy + -1.93352 * invRoot2Energy + 12.8105 * invRoot4Energy + -4.03756);}
197  } else if (detregion == 3 and background > 0.3) {
198  if (energy <= 0.022) {sigmaPhi = 1e-3 * (23.2771);}
199  else {sigmaPhi = 1e-3 * (-0.439611 * invEnergy + 5.74196 * invRoot2Energy + 1.76693 * invRoot4Energy + -0.0407866);}
200  }
201 
202  B2DEBUG(175, "energy=" << energy << ", detector region=" << detregion);
203  B2DEBUG(175, "sigmaEnergy=" << sigmaEnergy << ", sigmaPhi=" << sigmaPhi << ", sigmaTheta=" << sigmaTheta);
204 
205  double covMatrix[6] = {sigmaEnergy * sigmaEnergy, 0.0, sigmaPhi * sigmaPhi, 0.0, 0.0, sigmaTheta * sigmaTheta};
206  eclShower.setCovarianceMatrix(covMatrix);
207  }
208  }
209 }
210 
212 {
213  ;
214 }
215 
217 {
218  ;
219 }
Belle2::ECLCovarianceMatrixModule::m_eventLevelClusteringInfo
StoreObjPtr< EventLevelClusteringInfo > m_eventLevelClusteringInfo
Store object pointer: EventLevelClusteringInfo.
Definition: ECLCovarianceMatrixModule.h:77
Belle2::ECLCovarianceMatrixModule
Class to perform the shower correction.
Definition: ECLCovarianceMatrixModule.h:45
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::ECLCovarianceMatrixModule::endRun
virtual void endRun() override
End run.
Definition: ECLCovarianceMatrixModule.cc:211
Belle2::ECLCovarianceMatrixModule::initialize
virtual void initialize() override
Initialize.
Definition: ECLCovarianceMatrixModule.cc:52
Belle2::ECLCovarianceMatrixModule::eclShowerArrayName
virtual const char * eclShowerArrayName() const
Default name ECLShowers.
Definition: ECLCovarianceMatrixModule.h:81
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::ECLCovarianceMatrixModule::~ECLCovarianceMatrixModule
~ECLCovarianceMatrixModule()
Destructor.
Definition: ECLCovarianceMatrixModule.cc:48
Belle2::ECLCovarianceMatrixModule::m_fullBkgdCount
const double m_fullBkgdCount
Nominal Background at BGx1.0 (MC12)
Definition: ECLCovarianceMatrixModule.h:71
Belle2::ECLShower::c_nPhotons
@ c_nPhotons
CR is split into n photons (N1)
Definition: ECLShower.h:54
Belle2::ECLCovarianceMatrixModule::terminate
virtual void terminate() override
Terminate.
Definition: ECLCovarianceMatrixModule.cc:216
Belle2::ECLCovarianceMatrixModule::event
virtual void event() override
Event.
Definition: ECLCovarianceMatrixModule.cc:65
Belle2::ECLCovarianceMatrixModule::m_eclShowers
StoreArray< ECLShower > m_eclShowers
Store array: ECLShower.
Definition: ECLCovarianceMatrixModule.h:74
Belle2::ECLCovarianceMatrixModule::beginRun
virtual void beginRun() override
Begin run.
Definition: ECLCovarianceMatrixModule.cc:59
Belle2::ECLCovarianceMatrixModule::eventLevelClusteringInfoName
virtual const char * eventLevelClusteringInfoName() const
Name to be used for default option: EventLevelClusteringInfo.
Definition: ECLCovarianceMatrixModule.h:85