Belle II Software  release-08-01-10
EDepInGas.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/modules/cdcDigitizer/EDepInGas.h>
10 #include <framework/gearbox/Unit.h>
11 #include <framework/logging/Logger.h>
12 
13 #include <TDatabasePDG.h>
14 #include <TRandom.h>
15 //#include <iostream>
16 
17 using namespace std;
18 using namespace Belle2;
19 using namespace CDC;
20 
21 EDepInGas* EDepInGas::m_pntr = 0;
22 
23 EDepInGas& EDepInGas::getInstance()
24 {
25  if (!m_pntr) {
26  m_pntr = new EDepInGas();
27  m_pntr->initialize();
28  }
29  return *m_pntr;
30 }
31 
32 EDepInGas::EDepInGas() {}
33 
34 EDepInGas::~EDepInGas() {}
35 
36 void EDepInGas::initialize()
37 {
38 // gas
39  const double zHe = 2;
40  const double aHe = 4.003;
41  double fHe = 0.1108;
42 
43  const double zC = 6;
44  const double aC = 12.011;
45  double fC = 0.2223;
46 
47  const double zH = 1;
48  const double aH = 1.008;
49  double fH = 0.6669;
50 
51  double f = fHe + fC + fH;
52  fHe /= f;
53  fC /= f;
54  fH /= f;
55  m_z1 = fHe * zHe + fC * zC + fH * zH;
56  m_a1 = fHe * aHe + fC * aC + fH * aH;
57  m_i1 = 49.229;
58  m_rho1 = 0.703e-3;
59  // B2DEBUG(m_dbg, m_z1 <<" "<< m_a1 <<" "<< m_i1 <<" "<< m_rho1*m_z1/m_a1);
60 
61  // wire
62  const double zW = 74;
63  const double aW = 183.842;
64  double fW = 0.0198;
65 
66  const double zAl = 13;
67  const double aAl = 26.982;
68  double fAl = 0.9802;
69 
70  f = fW + fAl;
71  fW /= f;
72  fAl /= f;
73  m_z2 = fW * zW + fAl * zAl;
74  m_a2 = fW * aW + fAl * aAl;
75  m_i2 = 193.301;
76  m_rho2 = 0.419e-3;
77  // B2DEBUG(m_dbg, m_z2 <<" "<< m_a2 <<" "<< m_i2 <<" "<< m_rho2*m_z2/m_a2);
78 
79  fHe = 0.01031;
80  fC = 0.02068;
81  fH = 0.06204;
82  fW = 0.0014;
83  fAl = 0.0684;
84  f = fHe + fC + fH + fW + fAl;
85  fHe /= f;
86  fC /= f;
87  fH /= f;
88  fW /= f;
89  fAl /= f;
90  m_z0 = fHe * zHe + fC * zC + fH * zH + fW * zW + fAl * zAl;
91  m_a0 = fHe * aHe + fC * aC + fH * aH + fW * aW + fAl * aAl;
92  m_i0 = 76.693;
93  // m_rho0 = m_rho1 + m_rho2;
94  m_rho0 = 1.121e-3;
95 
96  //Define coefficients of density correction term
97  m_ak1[0][0] = 10.065;
98  m_ak1[0][1] = -13.63;
99  m_ak1[1][0] = 4.1334;
100  m_ak1[1][1] = -12.4;
101  m_ak1[2][0] = 0.8135;
102  m_ak1[2][1] = -11.33;
103  m_ak1[3][0] = 0.00402;
104  m_ak1[3][1] = -10.36;
105  m_ak1[4][0] = - 1.7564;
106  m_ak1[4][1] = - 9.611;
107  m_ak1[5][0] = - 1.5065;
108  m_ak1[5][1] = - 8.6753;
109 
110  m_bk1[0][0] = - 8.5619;
111  m_bk1[0][1] = 4.6;
112  m_bk1[1][0] = - 4.5857;
113  m_bk1[1][1] = 4.6;
114  m_bk1[2][0] = - 2.1447;
115  m_bk1[2][1] = 4.6;
116  m_bk1[3][0] = - 1.1582;
117  m_bk1[3][1] = 4.6;
118  m_bk1[4][0] = 0.67256;
119  m_bk1[4][1] = 4.6;
120  m_bk1[5][0] = 1.0156;
121  m_bk1[5][1] = 4.6;
122 
123  m_ck1[0][0] = 1.82804;
124  m_ck1[0][1] = 0.0;
125  m_ck1[1][0] = 1.27579;
126  m_ck1[1][1] = 0.0;
127  m_ck1[2][0] = 0.93676;
128  m_ck1[2][1] = 0.0;
129  m_ck1[3][0] = 0.79975;
130  m_ck1[3][1] = 0.0;
131  m_ck1[4][0] = 0.49093;
132  m_ck1[4][1] = 0.0;
133  m_ck1[5][0] = 0.44805;
134  m_ck1[5][1] = 0.0;
135 }
136 
137 double EDepInGas::getEDepInGas(int mode, int pdg, double mom, double dx, double e3) const
138 {
139  double xi1 = m_rho1 * (m_z1 / m_a1);
140  double xi2 = m_rho2 * (m_z2 / m_a2);
141  // B2DEBUG(m_dbg, (xi1+xi2)/xi1 <<" "<< e3);
142  if (mode == 0) return e3 * xi1 / (xi1 + xi2);
143 
144  static TDatabasePDG* ptrTDatabasePDG = TDatabasePDG::Instance();
145  const double mass = ptrTDatabasePDG->GetParticle(pdg)->Mass();
146  const double chrg = fabs(ptrTDatabasePDG->GetParticle(pdg)->Charge() / 3);
147 
148  const double me1 = getMostProbabEDep(mom, mass, chrg, dx, m_z1, m_a1, m_i1, m_rho1);
149  if (mode == 2) {
150  const double me0 = getMostProbabEDep(mom, mass, chrg, dx, m_z0, m_a0, m_i0, m_rho0);
151  return e3 * me1 / me0;
152  }
153  const double me2 = getMostProbabEDep(mom, mass, chrg, dx, m_z2, m_a2, m_i2, m_rho2);
154  if (mode == 1) return e3 * me1 / (me1 + me2);
155 
156  // extract energy deposit (e1) based on probab., L(e1)*L(e3-e1),
157  // where L is the analytic approx. of Landau dist.
158  // find lm1 which maximizes probability
159  const double mom2 = mom * mom;
160  const double beta2 = mom2 / (mass * mass + mom2);
161  const double buf = 0.1535e-3 * dx * chrg * chrg / beta2;
162  xi1 *= buf;
163  xi2 *= buf;
164  const double xi1i = 1 / xi1;
165  const double xi2i = 1 / xi2;
166  const double lm1min = -me1 * xi1i;
167  const double lm1max = (e3 - me1) * xi1i;
168  double lm1 = (e3 * me1 / (me1 + me2) - me1) * xi1i; // initial guess
169  lm1 = min(lm1, lm1max);
170  lm1 = max(lm1, lm1min);
171  // double edm = 1;
172  int itry = 0;
173  // B2DEBUG(m_dbg, "e3,me1,me2,xi1i=" << e3 <<" "<< me1 <<" "<< me2 <<" "<< xi1i);
174  // B2DEBUG(m_dbg, "lm1min,lm1max,lm1=" << lm1min <<" "<< lm1max <<" "<< lm1);
175  const double lm2b = (e3 - me1 - me2) * xi2i;
176  const double r = xi1 * xi2i;
177  double deltal = 0;
178  double lm1old = lm1;
179  while (itry <= 2 || fabs(lm1 - lm1old) > 0.001) {
180  ++itry;
181  lm1old = lm1;
182  if (itry > 50) {
183  // B2DEBUG(m_dbg, "itry > 50" << "itry,lm1min,lm1max,lm1,edm=" << itry <<" "<< lm1min <<" "<< lm1max <<" "<< lm1);
184  break;
185  }
186  lm1 -= deltal;
187  lm1 = max(lm1min, lm1);
188  lm1 = min(lm1max, lm1);
189  double lm2 = lm2b - r * lm1;
190  double emlm1 = exp(-lm1);
191  double emlm2 = exp(-lm2);
192  double dldlm1 = 1 - emlm1 - r + r * emlm2;
193  double dl2dlm12 = emlm1 + r * r * emlm2;
194  deltal = dldlm1 / dl2dlm12;
195  // edm = dldlm1 * deltal; //drop 0.5* for speed up
196  }
197 
198  if (mode == 3) return xi1 * lm1 + me1;
199 
200  double lm2 = lm2b - r * lm1;
201  double lmin = 0.5 * (lm1 + lm2 + exp(-lm1) + exp(-lm2));
202  const double ymax = exp(-lmin);
203  // if (ymax < 1.e-10) return e3*me1/(me1+me2);
204 
205  // generate e1 according to probab.
206  const double lm1w = lm1max - lm1min;
207  double y = 10;
208  double p = 1;
209  int jtry = 0;
210  double urndm[2];
211  while (y > p) {
212  ++jtry;
213  if (jtry > 500) {
214  // B2DEBUG(m_dbg, "jtry > 500" <<"ymax" << ymax);
215  lm1 = (e3 * me1 / (me1 + me2) - me1) * xi1i;
216  break;
217  }
218  gRandom->RndmArray(2, urndm);
219  lm1 = lm1min + lm1w * urndm[0];
220  lm2 = lm2b - r * lm1;
221  double l = 0.5 * (lm1 + lm2 + exp(-lm1) + exp(-lm2));
222  p = exp(-l);
223  y = ymax * urndm[1];
224  }
225 
226  return xi1 * lm1 + me1;
227 }
228 
229 double EDepInGas::getMostProbabEDep(double p, double mass, double zi, double dx, double Z,
230  double A, double I, double rho) const
231 {
232 // Calculate most probable energy deposit in GeV.
233 // Ref. J.Va'vra et al., NIM 203 (1982) 109.
234 // ----------------------------------------------------------------------
235 // <Input>
236 // p : total momentum (GeV)
237 // mass : mass of particle (GeV)
238 // zi : charge of particle (e)
239 // dx : thickness of gas (cm)
240 // Z : atomic no. of gas
241 // A : atomic weight of gas
242 // I : mean ionization potential of gas (eV)
243 // rho : density of gas (g/cm**3)
244 // ----------------------------------------------------------------------
245  const double psq = p * p;
246  const double masq = mass * mass;
247  const double tote = psq + masq;
248  const double beta2 = psq / tote;
249  const double bgam2 = psq / masq;
250  const double gamma = sqrt(tote) / mass;
251  const double alft = 0.1535e6 * rho * dx * Z / A;
252  const double k = log(m_massE * alft / I / I) + 0.891;
253  const double k1 = k - log(dx);
254  const double alngam = log(gamma);
255  const double alggam = log10(gamma);
256  // calculate density correction term
257  int irgn = 0;
258  double aa, bb, cc;
259  if (k1 < 6.) {
260  if (alggam > 3.6) irgn = 1;
261  aa = (m_ak1[1][irgn] - m_ak1[0][irgn]) * (6. - k1) + m_ak1[0][irgn];
262  bb = (m_bk1[1][irgn] - m_bk1[0][irgn]) * (6. - k1) + m_bk1[0][irgn];
263  cc = (m_ck1[1][irgn] - m_ck1[0][irgn]) * (6. - k1) + m_ck1[0][irgn];
264  } else if (k1 < 7.) {
265  if (alggam > 3.6) irgn = 1;
266  aa = (m_ak1[1][irgn] - m_ak1[0][irgn]) * (k1 - 6.) + m_ak1[0][irgn];
267  bb = (m_bk1[1][irgn] - m_bk1[0][irgn]) * (k1 - 6.) + m_bk1[0][irgn];
268  cc = (m_ck1[1][irgn] - m_ck1[0][irgn]) * (k1 - 6.) + m_ck1[0][irgn];
269  } else if (k1 < 8.) {
270  if (alggam > 3.6) irgn = 1;
271  aa = (m_ak1[2][irgn] - m_ak1[1][irgn]) * (k1 - 7.) + m_ak1[1][irgn];
272  bb = (m_bk1[2][irgn] - m_bk1[1][irgn]) * (k1 - 7.) + m_bk1[1][irgn];
273  cc = (m_ck1[2][irgn] - m_ck1[1][irgn]) * (k1 - 7.) + m_ck1[1][irgn];
274  } else if (k1 < 9.) {
275  if (alggam > 3.6) irgn = 1;
276  aa = (m_ak1[3][irgn] - m_ak1[2][irgn]) * (k1 - 8.) + m_ak1[2][irgn];
277  bb = (m_bk1[3][irgn] - m_bk1[2][irgn]) * (k1 - 8.) + m_bk1[2][irgn];
278  cc = (m_ck1[3][irgn] - m_ck1[2][irgn]) * (k1 - 8.) + m_ck1[2][irgn];
279  } else if (k1 < 10.) {
280  if (alggam > 3.6) irgn = 1;
281  aa = (m_ak1[4][irgn] - m_ak1[3][irgn]) * (k1 - 9.) + m_ak1[3][irgn];
282  bb = (m_bk1[4][irgn] - m_bk1[3][irgn]) * (k1 - 9.) + m_bk1[3][irgn];
283  cc = (m_ck1[4][irgn] - m_ck1[3][irgn]) * (k1 - 9.) + m_ck1[3][irgn];
284  } else if (k1 < 11.) {
285  if (alggam > 4.0) irgn = 1;
286  aa = (m_ak1[5][irgn] - m_ak1[4][irgn]) * (k1 - 10.) + m_ak1[4][irgn];
287  bb = (m_bk1[5][irgn] - m_bk1[4][irgn]) * (k1 - 10.) + m_bk1[4][irgn];
288  cc = (m_ck1[5][irgn] - m_ck1[4][irgn]) * (k1 - 10.) + m_ck1[4][irgn];
289  } else {
290  if (alggam > 4.0) irgn = 1;
291  aa = (m_ak1[5][irgn] - m_ak1[4][irgn]) * (k1 - 11.) + m_ak1[5][irgn];
292  bb = (m_bk1[5][irgn] - m_bk1[4][irgn]) * (k1 - 11.) + m_bk1[5][irgn];
293  cc = (m_ck1[5][irgn] - m_ck1[4][irgn]) * (k1 - 11.) + m_ck1[5][irgn];
294  }
295 
296  double delta = aa + alggam * (bb + cc * alggam);
297  delta = max(delta, 0.);
298  delta = min(delta, fabs(log(bgam2)));
299 
300  double val = zi * zi * alft * (k + 2 * alngam - beta2 - delta) / beta2;
301  val = max(val, 0.);
302  val *= Unit::eV;
303  return val;
304 }
Belle2::CDC::EDepInGas
The Class for Energy deposit in the gas.
Definition: EDepInGas.h:20
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17