Fork me on GitHub

source: git/modules/ParticlePropagator.cc

Last change on this file was ae93700, checked in by Pavel Demin <pavel.demin@…>, 3 years ago

improve ParticlePropagator code

  • Property mode set to 100644
File size: 10.8 KB
Line 
1/*
2 * Delphes: a framework for fast simulation of a generic collider experiment
3 * Copyright (C) 2012-2014 Universite catholique de Louvain (UCL), Belgium
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19/** \class ParticlePropagator
20 *
21 * Propagates charged and neutral particles
22 * from a given vertex to a cylinder defined by its radius,
23 * its half-length, centered at (0,0,0) and with its axis
24 * oriented along the z-axis.
25 *
26 * \author P. Demin - UCL, Louvain-la-Neuve
27 *
28 */
29
30#include "modules/ParticlePropagator.h"
31
32#include "classes/DelphesClasses.h"
33#include "classes/DelphesFactory.h"
34#include "classes/DelphesFormula.h"
35
36#include "ExRootAnalysis/ExRootClassifier.h"
37#include "ExRootAnalysis/ExRootFilter.h"
38#include "ExRootAnalysis/ExRootResult.h"
39
40#include "TDatabasePDG.h"
41#include "TFormula.h"
42#include "TLorentzVector.h"
43#include "TMath.h"
44#include "TObjArray.h"
45#include "TRandom3.h"
46#include "TString.h"
47
48#include <algorithm>
49#include <iostream>
50#include <sstream>
51#include <stdexcept>
52
53using namespace std;
54
55//------------------------------------------------------------------------------
56
57ParticlePropagator::ParticlePropagator() :
58 fItInputArray(0)
59{
60}
61
62//------------------------------------------------------------------------------
63
64ParticlePropagator::~ParticlePropagator()
65{
66}
67
68//------------------------------------------------------------------------------
69
70void ParticlePropagator::Init()
71{
72 fRadius = GetDouble("Radius", 1.0);
73 fRadius2 = fRadius * fRadius;
74 fHalfLength = GetDouble("HalfLength", 3.0);
75 fBz = GetDouble("Bz", 0.0);
76 if(fRadius < 1.0E-2)
77 {
78 cout << "ERROR: magnetic field radius is too low\n";
79 return;
80 }
81 if(fHalfLength < 1.0E-2)
82 {
83 cout << "ERROR: magnetic field length is too low\n";
84 return;
85 }
86
87 fRadiusMax = GetDouble("RadiusMax", fRadius);
88 fHalfLengthMax = GetDouble("HalfLengthMax", fHalfLength);
89
90 // import array with output from filter/classifier module
91
92 fInputArray = ImportArray(GetString("InputArray", "Delphes/stableParticles"));
93 fItInputArray = fInputArray->MakeIterator();
94
95 // import beamspot
96 try
97 {
98 fBeamSpotInputArray = ImportArray(GetString("BeamSpotInputArray", "BeamSpotFilter/beamSpotParticle"));
99 }
100 catch(runtime_error &e)
101 {
102 fBeamSpotInputArray = 0;
103 }
104 // create output arrays
105
106 fOutputArray = ExportArray(GetString("OutputArray", "stableParticles"));
107 fNeutralOutputArray = ExportArray(GetString("NeutralOutputArray", "neutralParticles"));
108 fChargedHadronOutputArray = ExportArray(GetString("ChargedHadronOutputArray", "chargedHadrons"));
109 fElectronOutputArray = ExportArray(GetString("ElectronOutputArray", "electrons"));
110 fMuonOutputArray = ExportArray(GetString("MuonOutputArray", "muons"));
111}
112
113//------------------------------------------------------------------------------
114
115void ParticlePropagator::Finish()
116{
117 if(fItInputArray) delete fItInputArray;
118}
119
120//------------------------------------------------------------------------------
121
122void ParticlePropagator::Process()
123{
124 Candidate *candidate, *mother, *particle;
125 TLorentzVector particlePosition, particleMomentum, beamSpotPosition;
126 Double_t px, py, pz, pt, pt2, e, q;
127 Double_t x, y, z, t, r;
128 Double_t x_c, y_c, r_c, phi_0;
129 Double_t x_t, y_t, z_t, r_t, phi_t;
130 Double_t t_r, t_z;
131 Double_t tmp;
132 Double_t gammam, omega;
133 Double_t xd, yd, zd;
134 Double_t l, d0, dz, ctgTheta, alpha;
135 Double_t bsx, bsy, bsz;
136 Double_t td, pio, phid, vz;
137
138 const Double_t c_light = 2.99792458E8;
139
140 if(!fBeamSpotInputArray || fBeamSpotInputArray->GetSize() == 0)
141 {
142 beamSpotPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
143 }
144 else
145 {
146 Candidate &beamSpotCandidate = *((Candidate *)fBeamSpotInputArray->At(0));
147 beamSpotPosition = beamSpotCandidate.Position;
148 }
149
150 fItInputArray->Reset();
151 while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
152 {
153 if(candidate->GetCandidates()->GetEntriesFast() == 0)
154 {
155 particle = candidate;
156 }
157 else
158 {
159 particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
160 }
161
162 particlePosition = particle->Position;
163 particleMomentum = particle->Momentum;
164
165 x = particlePosition.X() * 1.0E-3;
166 y = particlePosition.Y() * 1.0E-3;
167 z = particlePosition.Z() * 1.0E-3;
168
169 bsx = beamSpotPosition.X() * 1.0E-3;
170 bsy = beamSpotPosition.Y() * 1.0E-3;
171 bsz = beamSpotPosition.Z() * 1.0E-3;
172
173 q = particle->Charge;
174
175 // check that particle position is inside the cylinder
176 if(TMath::Hypot(x, y) > fRadiusMax || TMath::Abs(z) > fHalfLengthMax)
177 {
178 continue;
179 }
180
181 px = particleMomentum.Px();
182 py = particleMomentum.Py();
183 pz = particleMomentum.Pz();
184 pt = particleMomentum.Pt();
185 pt2 = particleMomentum.Perp2();
186 e = particleMomentum.E();
187
188 if(pt2 < 1.0E-9)
189 {
190 continue;
191 }
192
193 if(TMath::Hypot(x, y) > fRadius || TMath::Abs(z) > fHalfLength)
194 {
195 mother = candidate;
196 candidate = static_cast<Candidate *>(candidate->Clone());
197
198 candidate->InitialPosition = particlePosition;
199 candidate->Position = particlePosition;
200 candidate->L = 0.0;
201
202 candidate->Momentum = particleMomentum;
203 candidate->AddCandidate(mother);
204
205 fOutputArray->Add(candidate);
206 }
207 else if(TMath::Abs(q) < 1.0E-9 || TMath::Abs(fBz) < 1.0E-9)
208 {
209 // solve pt2*t^2 + 2*(px*x + py*y)*t - (fRadius2 - x*x - y*y) = 0
210 tmp = px * y - py * x;
211 t_r = (TMath::Sqrt(pt2 * fRadius2 - tmp * tmp) - px * x - py * y) / pt2;
212
213 t_z = (TMath::Sign(fHalfLength, pz) - z) / pz;
214
215 t = TMath::Min(t_r, t_z);
216
217 x_t = x + px * t;
218 y_t = y + py * t;
219 z_t = z + pz * t;
220
221 l = TMath::Sqrt((x_t - x) * (x_t - x) + (y_t - y) * (y_t - y) + (z_t - z) * (z_t - z));
222
223 mother = candidate;
224 candidate = static_cast<Candidate *>(candidate->Clone());
225
226 candidate->InitialPosition = particlePosition;
227 candidate->Position.SetXYZT(x_t * 1.0E3, y_t * 1.0E3, z_t * 1.0E3, particlePosition.T() + t * e * 1.0E3);
228 candidate->L = l * 1.0E3;
229
230 candidate->Momentum = particleMomentum;
231 candidate->AddCandidate(mother);
232
233 fOutputArray->Add(candidate);
234
235 if(TMath::Abs(q) > 1.0E-9)
236 {
237 switch(TMath::Abs(candidate->PID))
238 {
239 case 11:
240 fElectronOutputArray->Add(candidate);
241 break;
242 case 13:
243 fMuonOutputArray->Add(candidate);
244 break;
245 default:
246 fChargedHadronOutputArray->Add(candidate);
247 }
248 }
249 else
250 {
251 fNeutralOutputArray->Add(candidate);
252 }
253 }
254 else
255 {
256
257 // 1. initial transverse momentum p_{T0}: Part->pt
258 // initial transverse momentum direction phi_0 = -atan(p_{X0} / p_{Y0})
259 // relativistic gamma: gamma = E / mc^2; gammam = gamma * m
260 // gyration frequency omega = q * Bz / (gammam)
261 // helix radius r = p_{T0} / (omega * gammam)
262
263 gammam = e * 1.0E9 / (c_light * c_light); // gammam in [eV/c^2]
264 omega = q * fBz / gammam; // omega is here in [89875518/s]
265 r = pt / (q * fBz) * 1.0E9 / c_light; // in [m]
266
267 phi_0 = TMath::ATan2(py, px); // [rad] in [-pi, pi]
268
269 // 2. helix axis coordinates
270 x_c = x + r * TMath::Sin(phi_0);
271 y_c = y - r * TMath::Cos(phi_0);
272 r_c = TMath::Hypot(x_c, y_c);
273
274 // time of closest approach
275 td = (phi_0 + TMath::ATan2(x_c, y_c)) / omega;
276
277 // remove all the modulo pi that might have come from the atan
278 pio = TMath::Abs(TMath::Pi() / omega);
279 while(TMath::Abs(td) > 0.5 * pio)
280 {
281 td -= TMath::Sign(1.0, td) * pio;
282 }
283
284 vz = pz * c_light / e;
285
286 // calculate coordinates of closest approach to z axis
287 phid = phi_0 - omega * td;
288 xd = x_c - r * TMath::Sin(phid);
289 yd = y_c + r * TMath::Cos(phid);
290 zd = z + vz * td;
291
292 // momentum at closest approach
293 px = pt * TMath::Cos(phid);
294 py = pt * TMath::Sin(phid);
295
296 particleMomentum.SetPtEtaPhiE(pt, particleMomentum.Eta(), phid, particleMomentum.E());
297
298 // calculate additional track parameters (correct for beamspot position)
299 d0 = ((xd - bsx) * py - (yd - bsy) * px) / pt;
300 dz = zd - bsz;
301 ctgTheta = 1.0 / TMath::Tan(particleMomentum.Theta());
302
303 // 3. time evaluation t = TMath::Min(t_r, t_z)
304 // t_r : time to exit from the sides
305 // t_z : time to exit from the front or the back
306 t_z = (vz == 0.0) ? 1.0E99 : (TMath::Sign(fHalfLength, pz) - z) / vz;
307
308 if(r_c + TMath::Abs(r) < fRadius)
309 {
310 // helix does not cross the cylinder sides
311 t = t_z;
312 }
313 else
314 {
315 alpha = TMath::ACos((r * r + r_c * r_c - fRadius * fRadius) / (2 * TMath::Abs(r) * r_c));
316 t_r = td + TMath::Abs(alpha / omega);
317
318 t = TMath::Min(t_r, t_z);
319 }
320
321 // 4. position in terms of x(t), y(t), z(t)
322 phi_t = phi_0 - omega * t;
323 x_t = x_c - r * TMath::Sin(phi_t);
324 y_t = y_c + r * TMath::Cos(phi_t);
325 z_t = z + vz * t;
326 r_t = TMath::Hypot(x_t, y_t);
327
328 // lenght of the path from production to tracker
329 l = t * TMath::Hypot(vz, r * omega);
330
331 if(r_t > 0.0)
332 {
333 // store these variables before cloning
334 if(particle == candidate)
335 {
336 particle->D0 = d0 * 1.0E3;
337 particle->DZ = dz * 1.0E3;
338 particle->P = particleMomentum.P();
339 particle->PT = pt;
340 particle->CtgTheta = ctgTheta;
341 particle->Phi = particleMomentum.Phi();
342 }
343
344 mother = candidate;
345 candidate = static_cast<Candidate *>(candidate->Clone());
346
347 candidate->InitialPosition = particlePosition;
348 candidate->Position.SetXYZT(x_t * 1.0E3, y_t * 1.0E3, z_t * 1.0E3, particlePosition.T() + t * c_light * 1.0E3);
349
350 candidate->Momentum = particleMomentum;
351
352 candidate->L = l * 1.0E3;
353
354 candidate->Xd = xd * 1.0E3;
355 candidate->Yd = yd * 1.0E3;
356 candidate->Zd = zd * 1.0E3;
357
358 candidate->AddCandidate(mother);
359
360 fOutputArray->Add(candidate);
361 switch(TMath::Abs(candidate->PID))
362 {
363 case 11:
364 fElectronOutputArray->Add(candidate);
365 break;
366 case 13:
367 fMuonOutputArray->Add(candidate);
368 break;
369 default:
370 fChargedHadronOutputArray->Add(candidate);
371 }
372 }
373 }
374 }
375}
376
377//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.