1 | |
---|
2 | /** \class ParticlePropagator |
---|
3 | * |
---|
4 | * Propagates charged and neutral particles |
---|
5 | * from a given vertex to a cylinder defined by its radius, |
---|
6 | * its half-length, centered at (0,0,0) and with its axis |
---|
7 | * oriented along the z-axis. |
---|
8 | * |
---|
9 | * $Date: 2013-02-12 15:13:59 +0100 (Tue, 12 Feb 2013) $ |
---|
10 | * $Revision: 907 $ |
---|
11 | * |
---|
12 | * |
---|
13 | * \author P. Demin - UCL, Louvain-la-Neuve |
---|
14 | * |
---|
15 | */ |
---|
16 | |
---|
17 | #include "modules/ParticlePropagator.h" |
---|
18 | |
---|
19 | #include "classes/DelphesClasses.h" |
---|
20 | #include "classes/DelphesFactory.h" |
---|
21 | #include "classes/DelphesFormula.h" |
---|
22 | |
---|
23 | #include "ExRootAnalysis/ExRootResult.h" |
---|
24 | #include "ExRootAnalysis/ExRootFilter.h" |
---|
25 | #include "ExRootAnalysis/ExRootClassifier.h" |
---|
26 | |
---|
27 | #include "TMath.h" |
---|
28 | #include "TString.h" |
---|
29 | #include "TFormula.h" |
---|
30 | #include "TRandom3.h" |
---|
31 | #include "TObjArray.h" |
---|
32 | #include "TDatabasePDG.h" |
---|
33 | #include "TLorentzVector.h" |
---|
34 | |
---|
35 | #include <algorithm> |
---|
36 | #include <stdexcept> |
---|
37 | #include <iostream> |
---|
38 | #include <sstream> |
---|
39 | |
---|
40 | using namespace std; |
---|
41 | |
---|
42 | //------------------------------------------------------------------------------ |
---|
43 | |
---|
44 | ParticlePropagator::ParticlePropagator() : |
---|
45 | fItInputArray(0) |
---|
46 | { |
---|
47 | } |
---|
48 | |
---|
49 | //------------------------------------------------------------------------------ |
---|
50 | |
---|
51 | ParticlePropagator::~ParticlePropagator() |
---|
52 | { |
---|
53 | } |
---|
54 | |
---|
55 | //------------------------------------------------------------------------------ |
---|
56 | |
---|
57 | void ParticlePropagator::Init() |
---|
58 | { |
---|
59 | fRadius = GetDouble("Radius", 1.0); |
---|
60 | fRadius2 = fRadius*fRadius; |
---|
61 | fHalfLength = GetDouble("HalfLength", 3.0); |
---|
62 | fBz = GetDouble("Bz", 0.0); |
---|
63 | if(fRadius < 1.0E-2) |
---|
64 | { |
---|
65 | cout << "ERROR: magnetic field radius is too low\n"; |
---|
66 | return; |
---|
67 | } |
---|
68 | if(fHalfLength < 1.0E-2) |
---|
69 | { |
---|
70 | cout << "ERROR: magnetic field length is too low\n"; |
---|
71 | return; |
---|
72 | } |
---|
73 | |
---|
74 | // import array with output from filter/classifier module |
---|
75 | |
---|
76 | fInputArray = ImportArray(GetString("InputArray", "Delphes/stableParticles")); |
---|
77 | fItInputArray = fInputArray->MakeIterator(); |
---|
78 | |
---|
79 | // create output arrays |
---|
80 | |
---|
81 | fOutputArray = ExportArray(GetString("OutputArray", "stableParticles")); |
---|
82 | fChargedHadronOutputArray = ExportArray(GetString("ChargedHadronOutputArray", "chargedHadrons")); |
---|
83 | fElectronOutputArray = ExportArray(GetString("ElectronOutputArray", "electrons")); |
---|
84 | fMuonOutputArray = ExportArray(GetString("MuonOutputArray", "muons")); |
---|
85 | } |
---|
86 | |
---|
87 | //------------------------------------------------------------------------------ |
---|
88 | |
---|
89 | void ParticlePropagator::Finish() |
---|
90 | { |
---|
91 | if(fItInputArray) delete fItInputArray; |
---|
92 | } |
---|
93 | |
---|
94 | //------------------------------------------------------------------------------ |
---|
95 | |
---|
96 | void ParticlePropagator::Process() |
---|
97 | { |
---|
98 | Candidate *candidate, *mother; |
---|
99 | TLorentzVector candidatePosition, candidateMomentum; |
---|
100 | Double_t px, py, pz, pt, pt2, e, q; |
---|
101 | Double_t x, y, z, t, r, phi; |
---|
102 | Double_t x_c, y_c, r_c, phi_c, phi_0; |
---|
103 | Double_t x_t, y_t, z_t, r_t; |
---|
104 | Double_t t1, t2, t3, t4, t5, t6; |
---|
105 | Double_t t_z, t_r, t_ra, t_rb; |
---|
106 | Double_t tmp, discr, discr2; |
---|
107 | Double_t delta, gammam, omega, asinrho; |
---|
108 | |
---|
109 | const Double_t c_light = 2.99792458E8; |
---|
110 | |
---|
111 | fItInputArray->Reset(); |
---|
112 | while((candidate = static_cast<Candidate*>(fItInputArray->Next()))) |
---|
113 | { |
---|
114 | candidatePosition = candidate->Position; |
---|
115 | candidateMomentum = candidate->Momentum; |
---|
116 | x = candidatePosition.X()*1.0E-3; |
---|
117 | y = candidatePosition.Y()*1.0E-3; |
---|
118 | z = candidatePosition.Z()*1.0E-3; |
---|
119 | q = candidate->Charge; |
---|
120 | |
---|
121 | // check that particle position is inside the cylinder |
---|
122 | if(TMath::Hypot(x, y) > fRadius || TMath::Abs(z) > fHalfLength) |
---|
123 | { |
---|
124 | continue; |
---|
125 | } |
---|
126 | |
---|
127 | px = candidateMomentum.Px(); |
---|
128 | py = candidateMomentum.Py(); |
---|
129 | pz = candidateMomentum.Pz(); |
---|
130 | pt = candidateMomentum.Pt(); |
---|
131 | pt2 = candidateMomentum.Perp2(); |
---|
132 | e = candidateMomentum.E(); |
---|
133 | |
---|
134 | if(pt2 < 1.0E-9) |
---|
135 | { |
---|
136 | continue; |
---|
137 | } |
---|
138 | |
---|
139 | if(TMath::Abs(q) < 1.0E-9 || TMath::Abs(fBz) < 1.0E-9) |
---|
140 | { |
---|
141 | // solve pt2*t^2 + 2*(px*x + py*y)*t + (fRadius2 - x*x - y*y) = 0 |
---|
142 | tmp = px*y - py*x; |
---|
143 | discr2 = pt2*fRadius2 - tmp*tmp; |
---|
144 | |
---|
145 | if(discr2 < 0) |
---|
146 | { |
---|
147 | // no solutions |
---|
148 | continue; |
---|
149 | } |
---|
150 | |
---|
151 | tmp = px*x + py*y; |
---|
152 | discr = TMath::Sqrt(discr2); |
---|
153 | t1 = (-tmp + discr)/pt2; |
---|
154 | t2 = (-tmp - discr)/pt2; |
---|
155 | t = (t1 < 0) ? t2 : t1; |
---|
156 | |
---|
157 | z_t = z + pz*t; |
---|
158 | if(TMath::Abs(z_t) > fHalfLength) |
---|
159 | { |
---|
160 | t3 = (+fHalfLength - z) / pz; |
---|
161 | t4 = (-fHalfLength - z) / pz; |
---|
162 | t = (t3 < 0) ? t4 : t3; |
---|
163 | } |
---|
164 | |
---|
165 | x_t = x + px*t; |
---|
166 | y_t = y + py*t; |
---|
167 | z_t = z + pz*t; |
---|
168 | |
---|
169 | mother = candidate; |
---|
170 | candidate = static_cast<Candidate*>(candidate->Clone()); |
---|
171 | |
---|
172 | candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, 0.0); |
---|
173 | |
---|
174 | candidate->Momentum = candidateMomentum; |
---|
175 | candidate->AddCandidate(mother); |
---|
176 | |
---|
177 | fOutputArray->Add(candidate); |
---|
178 | if(TMath::Abs(q) > 1.0E-9) |
---|
179 | { |
---|
180 | switch(TMath::Abs(candidate->PID)) |
---|
181 | { |
---|
182 | case 11: |
---|
183 | fElectronOutputArray->Add(candidate); |
---|
184 | break; |
---|
185 | case 13: |
---|
186 | fMuonOutputArray->Add(candidate); |
---|
187 | break; |
---|
188 | default: |
---|
189 | fChargedHadronOutputArray->Add(candidate); |
---|
190 | } |
---|
191 | } |
---|
192 | } |
---|
193 | else |
---|
194 | { |
---|
195 | |
---|
196 | // 1. initial transverse momentum p_{T0} : Part->pt |
---|
197 | // initial transverse momentum direction \phi_0 = -atan(p_X0/p_Y0) |
---|
198 | // relativistic gamma : gamma = E/mc² ; gammam = gamma \times m |
---|
199 | // giration frequency \omega = q/(gamma m) fBz |
---|
200 | // helix radius r = p_T0 / (omega gamma m) |
---|
201 | |
---|
202 | gammam = e*1.0E9 / (c_light*c_light); // gammam in [eV/c²] |
---|
203 | omega = q * fBz / (gammam); // omega is here in [ 89875518 / s] |
---|
204 | r = pt / (q * fBz) * 1.0E9/c_light; // in [m] |
---|
205 | |
---|
206 | phi_0 = TMath::ATan2(py, px); // [rad] in [-pi; pi] |
---|
207 | |
---|
208 | // 2. helix axis coordinates |
---|
209 | x_c = x + r*TMath::Sin(phi_0); |
---|
210 | y_c = y - r*TMath::Cos(phi_0); |
---|
211 | r_c = TMath::Hypot(x_c, y_c); |
---|
212 | phi_c = TMath::ATan2(y_c, x_c); |
---|
213 | phi = phi_c; |
---|
214 | if(x_c < 0.0) phi += TMath::Pi(); |
---|
215 | |
---|
216 | // 3. time evaluation t = TMath::Min(t_r, t_z) |
---|
217 | // t_r : time to exit from the sides |
---|
218 | // t_z : time to exit from the front or the back |
---|
219 | t_r = 0.0; // in [ns] |
---|
220 | int sign_pz = (pz > 0.0) ? 1 : -1; |
---|
221 | if(pz == 0.0) t_z = 1.0E99; |
---|
222 | else t_z = gammam / (pz*1.0E9/c_light) * (-z + fHalfLength*sign_pz); |
---|
223 | |
---|
224 | if(r_c + TMath::Abs(r) < fRadius) |
---|
225 | { |
---|
226 | // helix does not cross the cylinder sides |
---|
227 | t = t_z; |
---|
228 | } |
---|
229 | else |
---|
230 | { |
---|
231 | asinrho = TMath::ASin( (fRadius*fRadius - r_c*r_c - r*r) / (2*TMath::Abs(r)*r_c) ); |
---|
232 | delta = phi_0 - phi; |
---|
233 | if(delta <-TMath::Pi()) delta += 2*TMath::Pi(); |
---|
234 | if(delta > TMath::Pi()) delta -= 2*TMath::Pi(); |
---|
235 | t1 = (delta + asinrho) / omega; |
---|
236 | t2 = (delta + TMath::Pi() - asinrho) / omega; |
---|
237 | t3 = (delta + TMath::Pi() + asinrho) / omega; |
---|
238 | t4 = (delta - asinrho) / omega; |
---|
239 | t5 = (delta - TMath::Pi() - asinrho) / omega; |
---|
240 | t6 = (delta - TMath::Pi() + asinrho) / omega; |
---|
241 | |
---|
242 | if(t1 < 0) t1 = 1.0E99; |
---|
243 | if(t2 < 0) t2 = 1.0E99; |
---|
244 | if(t3 < 0) t3 = 1.0E99; |
---|
245 | if(t4 < 0) t4 = 1.0E99; |
---|
246 | if(t5 < 0) t5 = 1.0E99; |
---|
247 | if(t6 < 0) t6 = 1.0E99; |
---|
248 | |
---|
249 | t_ra = TMath::Min(t1, TMath::Min(t2, t3)); |
---|
250 | t_rb = TMath::Min(t4, TMath::Min(t5, t6)); |
---|
251 | t_r = TMath::Min(t_ra, t_rb); |
---|
252 | t = TMath::Min(t_r, t_z); |
---|
253 | } |
---|
254 | |
---|
255 | // 4. position in terms of x(t), y(t), z(t) |
---|
256 | x_t = x_c + r * TMath::Sin(omega * t - phi_0); |
---|
257 | y_t = y_c + r * TMath::Cos(omega * t - phi_0); |
---|
258 | z_t = z + pz*1.0E9 / c_light / gammam * t; |
---|
259 | r_t = TMath::Hypot(x_t, y_t); |
---|
260 | |
---|
261 | if(r_t > 0.0) |
---|
262 | { |
---|
263 | mother = candidate; |
---|
264 | candidate = static_cast<Candidate*>(candidate->Clone()); |
---|
265 | |
---|
266 | candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, 0.0); |
---|
267 | |
---|
268 | candidate->Momentum = candidateMomentum; |
---|
269 | candidate->AddCandidate(mother); |
---|
270 | |
---|
271 | fOutputArray->Add(candidate); |
---|
272 | switch(TMath::Abs(candidate->PID)) |
---|
273 | { |
---|
274 | case 11: |
---|
275 | fElectronOutputArray->Add(candidate); |
---|
276 | break; |
---|
277 | case 13: |
---|
278 | fMuonOutputArray->Add(candidate); |
---|
279 | break; |
---|
280 | default: |
---|
281 | fChargedHadronOutputArray->Add(candidate); |
---|
282 | } |
---|
283 | } |
---|
284 | } |
---|
285 | } |
---|
286 | } |
---|
287 | |
---|
288 | //------------------------------------------------------------------------------ |
---|