Belle II Software  release-05-01-25
ECLCrystalData.cc
1 
2 //This module
3 #include <ecl/modules/eclBackgroundStudy/ECLCrystalData.h>
4 
5 //Framework
6 #include <framework/logging/Logger.h>
7 
8 //STL
9 #include <cmath>
10 #define PI 3.14159265358979323846
11 #define DEGTORAD 1.74532925199432955e-02
12 
13 using namespace Belle2;
14 
15 //ECL has 16-fold symmetry in phi. Multiply this by 16 to get the total number of crystals at each theta location.
16 int const ECLCrystalData::Ring[69] = {
17  3, 3, 4, 4, 4, 6, 6, 6, 6, 6, 6, 9, 9, //forward calorimeter
18  9, 9, 9, 9, 9 , 9, 9, 9, 9, 9 , 9, 9, 9, 9, 9 , 9, 9, 9, 9, 9, //barrel
19  9, 9, 9, 9, 9 , 9, 9, 9, 9, 9 , 9, 9, 9, 9, 9 , 9, 9, 9, 9, 9, //barrel
20  9, 9, 9, 9, 9 , 9, //barrel
21  9, 9, 6, 6, 6, 6, 6, 4, 4, 4 //backward calorimeter
22 };
23 
24 int const ECLCrystalData::sumPrevious[24] = {0, 3, 6, 10, 14, 18, 24, 30, 36, 42, 48, 54, 63, 72, 81, 90, 96, 102, 108, 114, 120, 124, 128, 132};
25 
26 //endcap crystal masses
27 float const ECLCrystalData::CrystalMassEndcap[132] = {4.99, 5.01, 4.99, 5.82, 5.94, 5.83, 4.75, 4.76, 4.77, 4.74, 5.14, 5.16, 5.16,
28  5.14, 5.76, 5.8, 5.79, 5.78, 4.14, 4.28, 4.15, 4.16, 4.28, 4.14, 4.32, 4.34,
29  4.34, 4.34, 4.34, 4.32, 4.85, 4.87, 4.87, 4.87, 4.87, 4.85, 4.99, 5.02, 5.02,
30  5.02, 5.02, 4.99, 5.37, 5.39, 5.39, 5.39, 5.39, 5.37, 5.89, 5.92, 5.92, 5.92,
31  5.92, 5.9, 4.03, 4.12, 4.05, 4.05, 4.12, 4.05, 4.05, 4.12, 4.03, 4.16, 4.18,
32  4.18, 4.18, 4.18, 4.18, 4.18, 4.18, 4.16, 4.47, 4.49, 4.49, 4.49, 4.49, 4.49,
33  4.49, 4.49, 4.47, 4.26, 4.34, 4.28, 4.28, 4.34, 4.28, 4.28, 4.34, 4.26, 5.88,
34  5.91, 5.91, 5.91, 5.91, 5.88, 5.46, 5.48, 5.48, 5.48, 5.48, 5.46, 5.02, 5.04,
35  5.04, 5.04, 5.04, 5.02, 4.63, 4.65, 4.65, 4.65, 4.65, 4.63, 4.4, 4.53, 4.42,
36  4.42, 4.53, 4.4, 5.74, 5.77, 5.76, 5.74, 5.09, 5.11, 5.11, 5.09, 4.56, 4.58,
37  4.58, 4.56
38  };
39 
40 //Barrel crystal masses
41 float const ECLCrystalData::CrystalMassBarrel[46] = {4.599, 4.619, 4.639, 4.659, 4.680, 4.701, 4.722, 4.744, 4.765,
42  4.786, 4.807, 4.828, 4.848, 4.868, 4.887, 4.906, 4.923, 4.940,
43  4.955, 4.969, 4.981, 4.992, 5.001, 5.008, 5.013, 5.016, 5.017,
44  5.016, 4.455, 4.455, 5.016, 5.017, 5.016, 5.013, 5.008, 5.001,
45  4.992, 4.981, 4.969, 4.955, 4.940, 4.923, 4.906, 4.887, 4.868,
46  4.848
47  };
48 
49 
50 //Theta ID to theta(deg)
51 double const ECLCrystalData::theta[69] = {13.9486666667, 15.4862, 17.110575, 18.757625, 20.255425, 21.8637333333,
52  23.4259666667, 24.96125, 26.4777833333, 27.9407, 29.36275, 30.8382555556,
53  32.2618444444, 33.667039, 35.062886, 36.509337, 38.007345, 39.557769,
54  41.161373, 42.818633, 44.530027, 46.295838, 48.115876, 49.990235, 51.91882,
55  53.900865, 55.935415, 58.021324, 60.157158, 62.3412, 64.571442, 66.8455,
56  69.160698, 71.513983, 73.902011, 76.321253, 78.767888, 81.237718, 83.726351,
57  86.229301, 88.741891, 90, 90, 91.258109, 93.770699, 96.273649, 98.762282,
58  101.232112, 103.678747, 106.097989, 108.486017, 110.839302, 113.154501,
59  115.428558, 117.658801, 119.842842, 121.978676, 124.064585, 126.099135,
60  128.2902222222, 130.4424444444, 132.6918333333, 134.9551666667, 137.3811666667,
61  139.9081666667, 142.5951666667, 145.374, 148.0185, 150.8355
62  };
63 
65 double const ECLCrystalData::EndcapRadius[132] = {
66 // 3,3,4,4,4,6,6,6,6,6,6,9,9, //forward calorimeter
67  217.25, 217.25, 217.25,
68  218.72, 218.72, 218.72,
69  220.29, 220.29, 220.29, 220.29,
70  222.00, 222.00, 222.00, 222.00,
71  223.96, 223.96, 223.96, 223.96,
72  226.01, 226.01, 226.01, 226.01, 226.01, 226.01,
73  228.21, 228.21, 228.21, 228.21, 228.21, 228.21,
74  230.67, 230.67, 230.67, 230.67, 230.67, 230.67,
75  233.20, 233.20, 233.20, 233.20, 233.20, 233.20,
76  235.89, 235.89, 235.89, 235.89, 235.89, 235.89,
77  238.89, 238.89, 238.89, 238.89, 238.89, 238.89,
78  241.94, 241.94, 241.94, 241.94, 241.94, 241.94, 241.94, 241.94, 241.94,
79  245.16, 245.16, 245.16, 245.16, 245.16, 245.16, 245.16, 245.16, 245.16,
80 // 9,9,6,6,6,6,6,4,4,4 //backward calorimeter
81  172.06, 172.06, 172.06, 172.06, 172.06, 172.06, 172.06, 172.06, 172.06,
82  165.69, 165.69, 165.69, 165.69, 165.69, 165.69, 165.69, 165.69, 165.69,
83  159.78, 159.78, 159.78, 159.78, 159.78, 159.78,
84  154.36, 154.36, 154.36, 154.36, 154.36, 154.36,
85  149.37, 149.37, 149.37, 149.37, 149.37, 149.37,
86  144.82, 144.82, 144.82, 144.82, 144.82, 144.82,
87  140.67, 140.67, 140.67, 140.67, 140.67, 140.67,
88  136.89, 136.89, 136.89, 136.89,
89  133.49, 133.49, 133.49, 133.49,
90  133.49, 133.49, 133.49, 133.49
91 };
92 
93 
95 double const ECLCrystalData::BarrelZ[46] = {210.15, 200.14, 190.44, 181.05, 171.95, 163.03, 154.45, 146.10, 137.98, 129.99, 122.27, 114.72, 107.34, 100.05, 92.96, 86.00, 79.15, 72.36, 65.72, 59.16, 52.69, 46.22, 39.87, 33.56, 27.29, 21.01, 14.81, 8.62, 2.78, -2.78, -14.81, -21.01, -27.29, -33.56, -39.87, -46.22, -52.69, -59.16, -65.72, -72.36, -79.15, -86.00, -92.95, -100.05, -107 - 28};
96 
97 
98 
99 
101 ECLCrystalData::ECLCrystalData()
102 {
103  m_cell_ID = -1;
104  m_phi_ID = -1;
105  m_theta_ID = -1;
106  m_phi_idx = -1;
107  m_theta_idx = -1;
108 
109  m_mass = -1;
110  m_volume = -1;
111  m_PosR = -1;
112  m_PosP = -1;
113  m_PosX = -1;
114  m_PosY = -1;
115  m_PosZ = -1;
116 }
117 
118 ECLCrystalData::ECLCrystalData(int cid)
119 {
120  m_cell_ID = cid;
121  Mapping(cid);
122 
123  Eval();
124 
125 }
126 
127 ECLCrystalData::ECLCrystalData(int tid, int pid)
128 {
129  m_phi_idx = -1;
130  m_theta_idx = -1;
131  m_theta_ID = tid;
132  m_phi_ID = pid;
133  m_cell_ID = GetCellID(tid, pid);
134 
135  Eval();
136 }
137 
138 ECLCrystalData::~ECLCrystalData()
139 {
140 }
141 
142 
144 void ECLCrystalData::Eval()
145 {
146  //Crystal properties
147  m_mass = EvalMass();
148  m_volume = m_mass / 4.51;
149 
150  //Crystal Location
151  m_PosP = EvalPhi();
152  m_PosR = EvalR();
153 
154  m_PosZ = EvalZ();
155  m_PosX = EvalX();
156  m_PosY = EvalY();
157 
158 
159 }
160 
161 double ECLCrystalData::EvalX()
162 {
163  return m_PosR * cos(DEGTORAD * m_PosP);
164 }
165 
166 double ECLCrystalData::EvalY()
167 {
168  return m_PosR * sin(DEGTORAD * m_PosP);
169 }
170 
171 
172 //Determine mass of crystal
173 float ECLCrystalData::EvalMass()
174 {
175  if (IsEndCap()) {
176  return CrystalMassEndcap[GetCrystalIndex()];
177 
178  } else {
179  return CrystalMassBarrel[GetCrystalIndex()];
180  }
181 }
182 
183 int ECLCrystalData::GetCrystalIndex()
184 {
185  if (m_theta_ID <= 12) { //forward ECL crystal
186  return m_phi_ID % Ring[m_theta_ID] + sumPrevious[m_theta_ID];
187 
188  } else if (m_theta_ID >= 59) { //backward ECL crystal
189  return m_phi_ID % Ring[m_theta_ID] + sumPrevious[m_theta_ID - 46];
190 
191  } else if (12 < m_theta_ID && m_theta_ID < 59) { //barrel ECL crystal
192  return m_theta_ID - 13;
193  }
194 
195  B2ERROR("ERROR: invalid theta_ID=" << m_theta_ID);
196  return -1;
197 }
198 
199 inline bool ECLCrystalData::IsEndCap()
200 {
201  return m_theta_ID < 13 || m_theta_ID > 58 ;
202 }
203 
204 inline bool ECLCrystalData::IsBarrel()
205 {
206  return !IsEndCap();
207 }
208 
209 
210 inline int ECLCrystalData::GetCellID()
211 {
212  return m_cell_ID;
213 }
214 
215 double ECLCrystalData::EvalPhi()
216 {
217  return 22.5 * m_phi_ID / Ring[m_theta_ID];
218 }
219 
220 double ECLCrystalData::EvalZ()
221 {
222  if (IsEndCap()) {
223  return EndcapRadius[GetCrystalIndex()] * cos(DEGTORAD * theta[m_theta_ID]);
224  }
225 
226  return BarrelZ[GetCrystalIndex()];
227 }
228 
229 double ECLCrystalData::EvalR()
230 {
231  if (IsEndCap()) {
232  return EndcapRadius[GetCrystalIndex()] * sin(DEGTORAD * theta[m_theta_ID]);
233  }
234 
235  return 147.7;
236 }
237 
238 
239 //Get CellID from theta and phi IDs
240 int ECLCrystalData::GetCellID(int ThetaId, int PhiId)
241 {
245  if (ThetaId < 13) {
246  int forwRing[13] = {0, 3, 6, 10, 14, 18, 24, 30, 36, 42, 48, 54, 63 };
247  return forwRing[ThetaId] * 16 + PhiId;
248 
249  } else if (ThetaId > 58) {
250  int backRing[10] = {0, 9, 18, 24, 30, 36, 42, 48, 52, 56} ;
251  return 7776 + backRing[ThetaId - 59] * 16 + PhiId;
252 
253  } else if (ThetaId > 12 && ThetaId < 59) {
254  return 1152 + 144 * (ThetaId - 13) + PhiId;
255 
256  } else
257  B2ERROR("ERROR (ECLCrystalData): CellID int ThetaId " << ThetaId << " int PhiId " << PhiId << ". Out of range.");
258 
259  return -1;
260 }
261 
262 void ECLCrystalData::Mapping(int cid)
263 {
264  if (cid < 0) {
265  B2ERROR("ECL ECLGeometryPar Mapping " << cid << ". Out of range.");
266 
267  } else if (cid < 3 * 16) { //Forkward start
268  m_theta_ID = 0;
269  m_phi_ID = cid;
270  m_phi_idx = m_phi_ID % 3;
271  m_theta_idx = m_phi_ID / 3;
272  } else if (cid < 6 * 16) {
273  m_theta_ID = 1;
274  m_phi_ID = cid - 3 * 16;
275  m_phi_idx = m_phi_ID % 3 + 3;
276  m_theta_idx = m_phi_ID / 3;
277  } else if (cid < 10 * 16) {
278  m_theta_ID = 2;
279  m_phi_ID = cid - 6 * 16 ;
280  m_phi_idx = m_phi_ID % 4 + 6;
281  m_theta_idx = m_phi_ID / 4;
282  } else if (cid < 14 * 16) {
283  m_theta_ID = 3;
284  m_phi_ID = cid - 10 * 16 ;
285  m_phi_idx = m_phi_ID % 4 + 10;
286  m_theta_idx = m_phi_ID / 4;
287  } else if (cid < 18 * 16) {
288  m_theta_ID = 4;
289  m_phi_ID = cid - 14 * 16 ;
290  m_phi_idx = m_phi_ID % 4 + 14;
291  m_theta_idx = m_phi_ID / 4;
292  } else if (cid < 24 * 16) {
293  m_theta_ID = 5;
294  m_phi_ID = cid - 18 * 16 ;
295  m_phi_idx = m_phi_ID % 6 + 18;
296  m_theta_idx = m_phi_ID / 6;
297  } else if (cid < 30 * 16) {
298  m_theta_ID = 6;
299  m_phi_ID = cid - 24 * 16 ;
300  m_phi_idx = m_phi_ID % 6 + 24;
301  m_theta_idx = m_phi_ID / 6;
302  } else if (cid < 36 * 16) {
303  m_theta_ID = 7;
304  m_phi_ID = cid - 30 * 16 ;
305  m_phi_idx = m_phi_ID % 6 + 30;
306  m_theta_idx = m_phi_ID / 6;
307  } else if (cid < 42 * 16) {
308  m_theta_ID = 8;
309  m_phi_ID = cid - 36 * 16 ;
310  m_phi_idx = m_phi_ID % 6 + 36;
311  m_theta_idx = m_phi_ID / 6;
312  } else if (cid < 48 * 16) {
313  m_theta_ID = 9;
314  m_phi_ID = cid - 42 * 16 ;
315  m_phi_idx = m_phi_ID % 6 + 42;
316  m_theta_idx = m_phi_ID / 6;
317  } else if (cid < 54 * 16) {
318  m_theta_ID = 10;
319  m_phi_ID = cid - 48 * 16 ;
320  m_phi_idx = m_phi_ID % 6 + 48;
321  m_theta_idx = m_phi_ID / 6;
322  } else if (cid < 63 * 16) {
323  m_theta_ID = 11;
324  m_phi_ID = cid - 54 * 16 ;
325  m_phi_idx = m_phi_ID % 9 + 54;
326  m_theta_idx = m_phi_ID / 9;
327  } else if (cid < 72 * 16) {
328  m_theta_ID = 12;
329  m_phi_ID = cid - 63 * 16 ;
330  m_phi_idx = m_phi_ID % 9 + 63;
331  m_theta_idx = m_phi_ID / 9;
332  } else if (cid < 7776) {//Barrel start
333  m_phi_ID = (cid - 1152) % 144;
334  m_theta_ID = (cid - 1152) / 144 + 13;
335  m_theta_idx = (cid - 1152) / 144;
336  m_phi_idx = m_phi_ID;
337  } else if (cid < 7776 + 9 * 16) { //Backward start
338  m_theta_ID = 59;
339  m_phi_ID = cid - 7776 ;
340  m_phi_idx = m_phi_ID % 9 + 72;
341  m_theta_idx = m_phi_ID / 9;
342  } else if (cid < 7776 + 18 * 16) {
343  m_theta_ID = 60;
344  m_phi_ID = cid - 7776 - 9 * 16 ;
345  m_phi_idx = m_phi_ID % 9 + 81;
346  m_theta_idx = m_phi_ID / 9;
347  } else if (cid < 7776 + 24 * 16) {
348  m_theta_ID = 61;
349  m_phi_ID = cid - 7776 - 18 * 16 ;
350  m_phi_idx = m_phi_ID % 6 + 90;
351  m_theta_idx = m_phi_ID / 6;
352  } else if (cid < 7776 + 30 * 16) {
353  m_theta_ID = 62;
354  m_phi_ID = cid - 7776 - 24 * 16 ;
355  m_phi_idx = m_phi_ID % 6 + 96;
356  m_theta_idx = m_phi_ID / 6;
357  } else if (cid < 7776 + 36 * 16) {
358  m_theta_ID = 63;
359  m_phi_ID = cid - 7776 - 30 * 16 ;
360  m_phi_idx = m_phi_ID % 6 + 102;
361  m_theta_idx = m_phi_ID / 6;
362  } else if (cid < 7776 + 42 * 16) {
363  m_theta_ID = 64;
364  m_phi_ID = cid - 7776 - 36 * 16 ;
365  m_phi_idx = m_phi_ID % 6 + 108;
366  m_theta_idx = m_phi_ID / 6;
367  } else if (cid < 7776 + 48 * 16) {
368  m_theta_ID = 65;
369  m_phi_ID = cid - 7776 - 42 * 16 ;
370  m_phi_idx = m_phi_ID % 6 + 114;
371  m_theta_idx = m_phi_ID / 6;
372  } else if (cid < 7776 + 52 * 16) {
373  m_theta_ID = 66;
374  m_phi_ID = cid - 7776 - 48 * 16 ;
375  m_phi_idx = m_phi_ID % 4 + 120;
376  m_theta_idx = m_phi_ID / 4;
377  } else if (cid < 7776 + 56 * 16) {
378  m_theta_ID = 67;
379  m_phi_ID = cid - 7776 - 52 * 16 ;
380  m_phi_idx = m_phi_ID % 4 + 124;
381  m_theta_idx = m_phi_ID / 4;
382  } else if (cid < 7776 + 60 * 16) {
383  m_theta_ID = 68;
384  m_phi_ID = cid - 7776 - 56 * 16 ;
385  m_phi_idx = m_phi_ID % 4 + 128;
386  m_theta_idx = m_phi_ID / 4;
387  } else {
388  B2ERROR("ECL ECLCrystalData Mapping " << cid << ". Out of range.");
389  }
390 }
Belle2::ECLCrystalData::EvalZ
double EvalZ()
Set Z value of crystal.
Definition: ECLCrystalData.cc:220
Belle2::ECLCrystalData::m_PosY
float m_PosY
Crystal Position - Y.
Definition: ECLCrystalData.h:64
Belle2::ECLCrystalData::Eval
void Eval()
Set values for the crystal.
Definition: ECLCrystalData.cc:144
Belle2::ECLCrystalData::m_phi_ID
int m_phi_ID
phi ID of crystal
Definition: ECLCrystalData.h:38
Belle2::ECLCrystalData::BarrelZ
static const double BarrelZ[46]
BarrelCrystalID to z-coordinate.
Definition: ECLCrystalData.h:109
Belle2::ECLCrystalData::EndcapRadius
static const double EndcapRadius[132]
EndCapCrystalID to radius (spherical)
Definition: ECLCrystalData.h:106
Belle2::ECLCrystalData::Mapping
void Mapping(int cid)
set theta and phi value of crystal
Definition: ECLCrystalData.cc:262
Belle2::ECLCrystalData::~ECLCrystalData
~ECLCrystalData()
destructor
Definition: ECLCrystalData.cc:138
Belle2::ECLCrystalData::IsBarrel
bool IsBarrel()
True if crystal is in barrel.
Definition: ECLCrystalData.cc:204
Belle2::ECLCrystalData::sumPrevious
static const int sumPrevious[24]
Number of crystals in preceding theta rings (endcaps only)
Definition: ECLCrystalData.h:88
Belle2::ECLCrystalData::EvalR
double EvalR()
Set R value of crystal.
Definition: ECLCrystalData.cc:229
Belle2::ECLCrystalData::Ring
static const int Ring[69]
ECL has 16-fold symmetry in phi.
Definition: ECLCrystalData.h:86
Belle2::ECLCrystalData::ECLCrystalData
ECLCrystalData()
Blank Constructor.
Definition: ECLCrystalData.cc:101
Belle2::ECLCrystalData::m_PosX
float m_PosX
Crystal Position - X.
Definition: ECLCrystalData.h:62
Belle2::ECLCrystalData::m_theta_ID
int m_theta_ID
theta ID of crystal
Definition: ECLCrystalData.h:40
Belle2::ECLCrystalData::GetCellID
int GetCellID()
return cell ID
Definition: ECLCrystalData.cc:210
Belle2::ECLCrystalData::m_volume
float m_volume
Crystal Volume.
Definition: ECLCrystalData.h:52
Belle2::ECLCrystalData::m_cell_ID
int m_cell_ID
Cell ID of crystal.
Definition: ECLCrystalData.h:36
Belle2::ECLCrystalData::m_PosP
float m_PosP
Crystal Position - P.
Definition: ECLCrystalData.h:60
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::ECLCrystalData::theta
static const double theta[69]
Theta ID to theta(deg)
Definition: ECLCrystalData.h:103
Belle2::ECLCrystalData::GetCrystalIndex
int GetCrystalIndex()
return crystal index
Definition: ECLCrystalData.cc:183
Belle2::ECLCrystalData::CrystalMassBarrel
static const float CrystalMassBarrel[46]
Endcap crystal mass.
Definition: ECLCrystalData.h:99
Belle2::ECLCrystalData::CrystalMassEndcap
static const float CrystalMassEndcap[132]
Barrel crystal mass.
Definition: ECLCrystalData.h:97
Belle2::ECLCrystalData::EvalX
double EvalX()
Set X value of crystal.
Definition: ECLCrystalData.cc:161
Belle2::ECLCrystalData::EvalPhi
double EvalPhi()
Set phi value of crystal.
Definition: ECLCrystalData.cc:215
Belle2::ECLCrystalData::m_mass
float m_mass
Crystal Mass.
Definition: ECLCrystalData.h:50
Belle2::ECLCrystalData::m_phi_idx
int m_phi_idx
Index.
Definition: ECLCrystalData.h:44
Belle2::ECLCrystalData::m_PosZ
float m_PosZ
Crystal Position - Z.
Definition: ECLCrystalData.h:66
Belle2::ECLCrystalData::EvalY
double EvalY()
Set Y value of crystal.
Definition: ECLCrystalData.cc:166
Belle2::ECLCrystalData::EvalMass
float EvalMass()
Evaluate mass of crystal.
Definition: ECLCrystalData.cc:173
Belle2::ECLCrystalData::IsEndCap
bool IsEndCap()
True if crystal is in endcaps.
Definition: ECLCrystalData.cc:199
Belle2::ECLCrystalData::m_PosR
float m_PosR
Crystal Position - R.
Definition: ECLCrystalData.h:58
Belle2::ECLCrystalData::m_theta_idx
int m_theta_idx
Index.
Definition: ECLCrystalData.h:46