Fork me on GitHub

source: git/modules/ParticlePropagator.cc@ c614dd7

ImprovedOutputFile Timing
Last change on this file since c614dd7 was 341014c, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

apply .clang-format to all .h, .cc and .cpp files

  • Property mode set to 100644
File size: 12.3 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, phi;
128 Double_t x_c, y_c, r_c, phi_c, phi_0;
129 Double_t x_t, y_t, z_t, r_t;
130 Double_t t1, t2, t3, t4, t5, t6;
131 Double_t t_z, t_r, t_ra, t_rb;
132 Double_t tmp, discr, discr2;
133 Double_t delta, gammam, omega, asinrho;
134 Double_t rcu, rc2, xd, yd, zd;
135 Double_t l, d0, dz, p, ctgTheta, phip, etap, alpha;
136 Double_t bsx, bsy, bsz;
137
138 const Double_t c_light = 2.99792458E8;
139
140 if(!fBeamSpotInputArray || fBeamSpotInputArray->GetSize() == 0)
141 beamSpotPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
142 else
143 {
144 Candidate &beamSpotCandidate = *((Candidate *)fBeamSpotInputArray->At(0));
145 beamSpotPosition = beamSpotCandidate.Position;
146 }
147
148 fItInputArray->Reset();
149 while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
150 {
151 if(candidate->GetCandidates()->GetEntriesFast() == 0)
152 {
153 particle = candidate;
154 }
155 else
156 {
157 particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
158 }
159
160 particlePosition = particle->Position;
161 particleMomentum = particle->Momentum;
162 x = particlePosition.X() * 1.0E-3;
163 y = particlePosition.Y() * 1.0E-3;
164 z = particlePosition.Z() * 1.0E-3;
165
166 bsx = beamSpotPosition.X() * 1.0E-3;
167 bsy = beamSpotPosition.Y() * 1.0E-3;
168 bsz = beamSpotPosition.Z() * 1.0E-3;
169
170 q = particle->Charge;
171
172 // check that particle position is inside the cylinder
173 if(TMath::Hypot(x, y) > fRadiusMax || TMath::Abs(z) > fHalfLengthMax)
174 {
175 continue;
176 }
177
178 px = particleMomentum.Px();
179 py = particleMomentum.Py();
180 pz = particleMomentum.Pz();
181 pt = particleMomentum.Pt();
182 pt2 = particleMomentum.Perp2();
183 e = particleMomentum.E();
184
185 if(pt2 < 1.0E-9)
186 {
187 continue;
188 }
189
190 if(TMath::Hypot(x, y) > fRadius || TMath::Abs(z) > fHalfLength)
191 {
192 mother = candidate;
193 candidate = static_cast<Candidate *>(candidate->Clone());
194
195 candidate->InitialPosition = particlePosition;
196 candidate->Position = particlePosition;
197 candidate->L = 0.0;
198
199 candidate->Momentum = particleMomentum;
200 candidate->AddCandidate(mother);
201
202 fOutputArray->Add(candidate);
203 }
204 else if(TMath::Abs(q) < 1.0E-9 || TMath::Abs(fBz) < 1.0E-9)
205 {
206 // solve pt2*t^2 + 2*(px*x + py*y)*t - (fRadius2 - x*x - y*y) = 0
207 tmp = px * y - py * x;
208 discr2 = pt2 * fRadius2 - tmp * tmp;
209
210 if(discr2 < 0.0)
211 {
212 // no solutions
213 continue;
214 }
215
216 tmp = px * x + py * y;
217 discr = TMath::Sqrt(discr2);
218 t1 = (-tmp + discr) / pt2;
219 t2 = (-tmp - discr) / pt2;
220 t = (t1 < 0.0) ? t2 : t1;
221
222 z_t = z + pz * t;
223 if(TMath::Abs(z_t) > fHalfLength)
224 {
225 t3 = (+fHalfLength - z) / pz;
226 t4 = (-fHalfLength - z) / pz;
227 t = (t3 < 0.0) ? t4 : t3;
228 }
229
230 x_t = x + px * t;
231 y_t = y + py * t;
232 z_t = z + pz * t;
233
234 l = TMath::Sqrt((x_t - x) * (x_t - x) + (y_t - y) * (y_t - y) + (z_t - z) * (z_t - z));
235
236 mother = candidate;
237 candidate = static_cast<Candidate *>(candidate->Clone());
238
239 candidate->InitialPosition = particlePosition;
240 candidate->Position.SetXYZT(x_t * 1.0E3, y_t * 1.0E3, z_t * 1.0E3, particlePosition.T() + t * e * 1.0E3);
241 candidate->L = l * 1.0E3;
242
243 candidate->Momentum = particleMomentum;
244 candidate->AddCandidate(mother);
245
246 fOutputArray->Add(candidate);
247 if(TMath::Abs(q) > 1.0E-9)
248 {
249 switch(TMath::Abs(candidate->PID))
250 {
251 case 11:
252 fElectronOutputArray->Add(candidate);
253 break;
254 case 13:
255 fMuonOutputArray->Add(candidate);
256 break;
257 default:
258 fChargedHadronOutputArray->Add(candidate);
259 }
260 }
261 else
262 {
263 fNeutralOutputArray->Add(candidate);
264 }
265 }
266 else
267 {
268
269 // 1. initial transverse momentum p_{T0}: Part->pt
270 // initial transverse momentum direction phi_0 = -atan(p_X0/p_Y0)
271 // relativistic gamma: gamma = E/mc^2; gammam = gamma * m
272 // gyration frequency omega = q/(gamma m) fBz
273 // helix radius r = p_{T0} / (omega gamma m)
274
275 gammam = e * 1.0E9 / (c_light * c_light); // gammam in [eV/c^2]
276 omega = q * fBz / (gammam); // omega is here in [89875518/s]
277 r = pt / (q * fBz) * 1.0E9 / c_light; // in [m]
278
279 phi_0 = TMath::ATan2(py, px); // [rad] in [-pi, pi]
280
281 // 2. helix axis coordinates
282 x_c = x + r * TMath::Sin(phi_0);
283 y_c = y - r * TMath::Cos(phi_0);
284 r_c = TMath::Hypot(x_c, y_c);
285 phi_c = TMath::ATan2(y_c, x_c);
286 phi = phi_c;
287 if(x_c < 0.0) phi += TMath::Pi();
288
289 rcu = TMath::Abs(r);
290 rc2 = r_c * r_c;
291
292 // calculate coordinates of closest approach to track circle in transverse plane xd, yd, zd
293 xd = x_c * x_c * x_c - x_c * rcu * r_c + x_c * y_c * y_c;
294 xd = (rc2 > 0.0) ? xd / rc2 : -999;
295 yd = y_c * (-rcu * r_c + rc2);
296 yd = (rc2 > 0.0) ? yd / rc2 : -999;
297 zd = z + (TMath::Sqrt(xd * xd + yd * yd) - TMath::Sqrt(x * x + y * y)) * pz / pt;
298
299 // use perigee momentum rather than original particle
300 // momentum, since the orignal particle momentum isn't known
301
302 px = TMath::Sign(1.0, r) * pt * (-y_c / r_c);
303 py = TMath::Sign(1.0, r) * pt * (x_c / r_c);
304 etap = particleMomentum.Eta();
305 phip = TMath::ATan2(py, px);
306
307 particleMomentum.SetPtEtaPhiE(pt, etap, phip, particleMomentum.E());
308
309 // calculate additional track parameters (correct for beamspot position)
310
311 d0 = ((x - bsx) * py - (y - bsy) * px) / pt;
312 dz = z - ((x - bsx) * px + (y - bsy) * py) / pt * (pz / pt);
313 p = particleMomentum.P();
314 ctgTheta = 1.0 / TMath::Tan(particleMomentum.Theta());
315
316 // 3. time evaluation t = TMath::Min(t_r, t_z)
317 // t_r : time to exit from the sides
318 // t_z : time to exit from the front or the back
319 t_r = 0.0; // in [ns]
320 int sign_pz = (pz > 0.0) ? 1 : -1;
321 if(pz == 0.0)
322 t_z = 1.0E99;
323 else
324 t_z = gammam / (pz * 1.0E9 / c_light) * (-z + fHalfLength * sign_pz);
325
326 if(r_c + TMath::Abs(r) < fRadius)
327 {
328 // helix does not cross the cylinder sides
329 t = t_z;
330 }
331 else
332 {
333 asinrho = TMath::ASin((fRadius * fRadius - r_c * r_c - r * r) / (2 * TMath::Abs(r) * r_c));
334 delta = phi_0 - phi;
335 if(delta < -TMath::Pi()) delta += 2 * TMath::Pi();
336 if(delta > TMath::Pi()) delta -= 2 * TMath::Pi();
337 t1 = (delta + asinrho) / omega;
338 t2 = (delta + TMath::Pi() - asinrho) / omega;
339 t3 = (delta + TMath::Pi() + asinrho) / omega;
340 t4 = (delta - asinrho) / omega;
341 t5 = (delta - TMath::Pi() - asinrho) / omega;
342 t6 = (delta - TMath::Pi() + asinrho) / omega;
343
344 if(t1 < 0.0) t1 = 1.0E99;
345 if(t2 < 0.0) t2 = 1.0E99;
346 if(t3 < 0.0) t3 = 1.0E99;
347 if(t4 < 0.0) t4 = 1.0E99;
348 if(t5 < 0.0) t5 = 1.0E99;
349 if(t6 < 0.0) t6 = 1.0E99;
350
351 t_ra = TMath::Min(t1, TMath::Min(t2, t3));
352 t_rb = TMath::Min(t4, TMath::Min(t5, t6));
353 t_r = TMath::Min(t_ra, t_rb);
354 t = TMath::Min(t_r, t_z);
355 }
356
357 // 4. position in terms of x(t), y(t), z(t)
358 x_t = x_c + r * TMath::Sin(omega * t - phi_0);
359 y_t = y_c + r * TMath::Cos(omega * t - phi_0);
360 z_t = z + pz * 1.0E9 / c_light / gammam * t;
361 r_t = TMath::Hypot(x_t, y_t);
362
363 // compute path length for an helix
364
365 alpha = pz * 1.0E9 / c_light / gammam;
366 l = t * TMath::Sqrt(alpha * alpha + r * r * omega * omega);
367
368 if(r_t > 0.0)
369 {
370
371 // store these variables before cloning
372 if(particle == candidate)
373 {
374 particle->D0 = d0 * 1.0E3;
375 particle->DZ = dz * 1.0E3;
376 particle->P = p;
377 particle->PT = pt;
378 particle->CtgTheta = ctgTheta;
379 particle->Phi = phip;
380 }
381
382 mother = candidate;
383 candidate = static_cast<Candidate *>(candidate->Clone());
384
385 candidate->InitialPosition = particlePosition;
386 candidate->Position.SetXYZT(x_t * 1.0E3, y_t * 1.0E3, z_t * 1.0E3, particlePosition.T() + t * c_light * 1.0E3);
387
388 candidate->Momentum = particleMomentum;
389
390 candidate->L = l * 1.0E3;
391
392 candidate->Xd = xd * 1.0E3;
393 candidate->Yd = yd * 1.0E3;
394 candidate->Zd = zd * 1.0E3;
395
396 candidate->AddCandidate(mother);
397
398 fOutputArray->Add(candidate);
399 switch(TMath::Abs(candidate->PID))
400 {
401 case 11:
402 fElectronOutputArray->Add(candidate);
403 break;
404 case 13:
405 fMuonOutputArray->Add(candidate);
406 break;
407 default:
408 fChargedHadronOutputArray->Add(candidate);
409 }
410 }
411 }
412 }
413}
414
415//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.