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