Fork me on GitHub

source: git/modules/ParticlePropagator.cc@ ef8a06d

ImprovedOutputFile Timing dual_readout
Last change on this file since ef8a06d was 5b51d33, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

check if candidate is particle before setting impact parameters

  • Property mode set to 100644
File size: 12.2 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
20/** \class ParticlePropagator
21 *
22 * Propagates charged and neutral particles
23 * from a given vertex to a cylinder defined by its radius,
24 * its half-length, centered at (0,0,0) and with its axis
25 * oriented along the z-axis.
26 *
27 * \author P. Demin - UCL, Louvain-la-Neuve
28 *
29 */
30
31#include "modules/ParticlePropagator.h"
32
33#include "classes/DelphesClasses.h"
34#include "classes/DelphesFactory.h"
35#include "classes/DelphesFormula.h"
36
37#include "ExRootAnalysis/ExRootResult.h"
38#include "ExRootAnalysis/ExRootFilter.h"
39#include "ExRootAnalysis/ExRootClassifier.h"
40
41#include "TMath.h"
42#include "TString.h"
43#include "TFormula.h"
44#include "TRandom3.h"
45#include "TObjArray.h"
46#include "TDatabasePDG.h"
47#include "TLorentzVector.h"
48
49#include <algorithm>
50#include <stdexcept>
51#include <iostream>
52#include <sstream>
53
54using namespace std;
55
56//------------------------------------------------------------------------------
57
58ParticlePropagator::ParticlePropagator() :
59 fItInputArray(0)
60{
61}
62
63//------------------------------------------------------------------------------
64
65ParticlePropagator::~ParticlePropagator()
66{
67}
68
69
70//------------------------------------------------------------------------------
71
72void ParticlePropagator::Init()
73{
74 fRadius = GetDouble("Radius", 1.0);
75 fRadius2 = fRadius*fRadius;
76 fHalfLength = GetDouble("HalfLength", 3.0);
77 fBz = GetDouble("Bz", 0.0);
78 if(fRadius < 1.0E-2)
79 {
80 cout << "ERROR: magnetic field radius is too low\n";
81 return;
82 }
83 if(fHalfLength < 1.0E-2)
84 {
85 cout << "ERROR: magnetic field length is too low\n";
86 return;
87 }
88
89 fRadiusMax = GetDouble("RadiusMax", fRadius);
90 fHalfLengthMax = GetDouble("HalfLengthMax", fHalfLength);
91
92 // import array with output from filter/classifier module
93
94 fInputArray = ImportArray(GetString("InputArray", "Delphes/stableParticles"));
95 fItInputArray = fInputArray->MakeIterator();
96
97 // import beamspot
98 try
99 {
100 fBeamSpotInputArray = ImportArray(GetString("BeamSpotInputArray", "BeamSpotFilter/beamSpotParticle"));
101 }
102 catch(runtime_error &e)
103 {
104 fBeamSpotInputArray = 0;
105 }
106 // create output arrays
107
108 fOutputArray = ExportArray(GetString("OutputArray", "stableParticles"));
109 fNeutralOutputArray = ExportArray(GetString("NeutralOutputArray", "neutralParticles"));
110 fChargedHadronOutputArray = ExportArray(GetString("ChargedHadronOutputArray", "chargedHadrons"));
111 fElectronOutputArray = ExportArray(GetString("ElectronOutputArray", "electrons"));
112 fMuonOutputArray = ExportArray(GetString("MuonOutputArray", "muons"));
113}
114
115//------------------------------------------------------------------------------
116
117void ParticlePropagator::Finish()
118{
119 if(fItInputArray) delete fItInputArray;
120}
121
122//------------------------------------------------------------------------------
123
124void ParticlePropagator::Process()
125{
126 Candidate *candidate, *mother, *particle;
127 TLorentzVector particlePosition, particleMomentum, beamSpotPosition;
128 Double_t px, py, pz, pt, pt2, e, q;
129 Double_t x, y, z, t, r, phi;
130 Double_t x_c, y_c, r_c, phi_c, phi_0;
131 Double_t x_t, y_t, z_t, r_t;
132 Double_t t1, t2, t3, t4, t5, t6;
133 Double_t t_z, t_r, t_ra, t_rb;
134 Double_t tmp, discr, discr2;
135 Double_t delta, gammam, omega, asinrho;
136 Double_t rcu, rc2, xd, yd, zd;
137 Double_t l, d0, dz, p, ctgTheta, phip, etap, alpha;
138 Double_t bsx, bsy, bsz;
139
140 const Double_t c_light = 2.99792458E8;
141
142 if (!fBeamSpotInputArray || fBeamSpotInputArray->GetSize () == 0)
143 beamSpotPosition.SetXYZT(0.0, 0.0, 0.0, 0.0);
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 x = particlePosition.X()*1.0E-3;
165 y = particlePosition.Y()*1.0E-3;
166 z = particlePosition.Z()*1.0E-3;
167
168 bsx = beamSpotPosition.X()*1.0E-3;
169 bsy = beamSpotPosition.Y()*1.0E-3;
170 bsz = beamSpotPosition.Z()*1.0E-3;
171
172 q = particle->Charge;
173
174 // check that particle position is inside the cylinder
175 if(TMath::Hypot(x, y) > fRadiusMax || TMath::Abs(z) > fHalfLengthMax)
176 {
177 continue;
178 }
179
180 px = particleMomentum.Px();
181 py = particleMomentum.Py();
182 pz = particleMomentum.Pz();
183 pt = particleMomentum.Pt();
184 pt2 = particleMomentum.Perp2();
185 e = particleMomentum.E();
186
187 if(pt2 < 1.0E-9)
188 {
189 continue;
190 }
191
192 if(TMath::Hypot(x, y) > fRadius || TMath::Abs(z) > fHalfLength)
193 {
194 mother = candidate;
195 candidate = static_cast<Candidate*>(candidate->Clone());
196
197 candidate->InitialPosition = particlePosition;
198 candidate->Position = particlePosition;
199 candidate->L = 0.0;
200
201 candidate->Momentum = particleMomentum;
202 candidate->AddCandidate(mother);
203
204 fOutputArray->Add(candidate);
205 }
206 else if(TMath::Abs(q) < 1.0E-9 || TMath::Abs(fBz) < 1.0E-9)
207 {
208 // solve pt2*t^2 + 2*(px*x + py*y)*t - (fRadius2 - x*x - y*y) = 0
209 tmp = px*y - py*x;
210 discr2 = pt2*fRadius2 - tmp*tmp;
211
212 if(discr2 < 0.0)
213 {
214 // no solutions
215 continue;
216 }
217
218 tmp = px*x + py*y;
219 discr = TMath::Sqrt(discr2);
220 t1 = (-tmp + discr)/pt2;
221 t2 = (-tmp - discr)/pt2;
222 t = (t1 < 0.0) ? t2 : t1;
223
224 z_t = z + pz*t;
225 if(TMath::Abs(z_t) > fHalfLength)
226 {
227 t3 = (+fHalfLength - z) / pz;
228 t4 = (-fHalfLength - z) / pz;
229 t = (t3 < 0.0) ? t4 : t3;
230 }
231
232 x_t = x + px*t;
233 y_t = y + py*t;
234 z_t = z + pz*t;
235
236 l = TMath::Sqrt( (x_t - x)*(x_t - x) + (y_t - y)*(y_t - y) + (z_t - z)*(z_t - z));
237
238 mother = candidate;
239 candidate = static_cast<Candidate*>(candidate->Clone());
240
241 candidate->InitialPosition = particlePosition;
242 candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, particlePosition.T() + t*e*1.0E3);
243 candidate->L = l*1.0E3;
244
245 candidate->Momentum = particleMomentum;
246 candidate->AddCandidate(mother);
247
248 fOutputArray->Add(candidate);
249 if(TMath::Abs(q) > 1.0E-9)
250 {
251 switch(TMath::Abs(candidate->PID))
252 {
253 case 11:
254 fElectronOutputArray->Add(candidate);
255 break;
256 case 13:
257 fMuonOutputArray->Add(candidate);
258 break;
259 default:
260 fChargedHadronOutputArray->Add(candidate);
261 }
262 }
263 else
264 {
265 fNeutralOutputArray->Add(candidate);
266 }
267 }
268 else
269 {
270
271 // 1. initial transverse momentum p_{T0}: Part->pt
272 // initial transverse momentum direction phi_0 = -atan(p_X0/p_Y0)
273 // relativistic gamma: gamma = E/mc^2; gammam = gamma * m
274 // gyration frequency omega = q/(gamma m) fBz
275 // helix radius r = p_{T0} / (omega gamma m)
276
277 gammam = e*1.0E9 / (c_light*c_light); // gammam in [eV/c^2]
278 omega = q * fBz / (gammam); // omega is here in [89875518/s]
279 r = pt / (q * fBz) * 1.0E9/c_light; // in [m]
280
281 phi_0 = TMath::ATan2(py, px); // [rad] in [-pi, pi]
282
283 // 2. helix axis coordinates
284 x_c = x + r*TMath::Sin(phi_0);
285 y_c = y - r*TMath::Cos(phi_0);
286 r_c = TMath::Hypot(x_c, y_c);
287 phi_c = TMath::ATan2(y_c, x_c);
288 phi = phi_c;
289 if(x_c < 0.0) phi += TMath::Pi();
290
291 rcu = TMath::Abs(r);
292 rc2 = r_c*r_c;
293
294 // calculate coordinates of closest approach to track circle in transverse plane xd, yd, zd
295 xd = x_c*x_c*x_c - x_c*rcu*r_c + x_c*y_c*y_c;
296 xd = (rc2 > 0.0) ? xd / rc2 : -999;
297 yd = y_c*(-rcu*r_c + rc2);
298 yd = (rc2 > 0.0) ? yd / rc2 : -999;
299 zd = z + (TMath::Sqrt(xd*xd + yd*yd) - TMath::Sqrt(x*x + y*y))*pz/pt;
300
301 // use perigee momentum rather than original particle
302 // momentum, since the orignal particle momentum isn't known
303
304 px = TMath::Sign(1.0, r) * pt * (-y_c / r_c);
305 py = TMath::Sign(1.0, r) * pt * (x_c / r_c);
306 etap = particleMomentum.Eta();
307 phip = TMath::ATan2(py, px);
308
309 particleMomentum.SetPtEtaPhiE(pt, etap, phip, particleMomentum.E());
310
311 // calculate additional track parameters (correct for beamspot position)
312
313 d0 = ((x - bsx) * py - (y - bsy) * px) / pt;
314 dz = z - ((x - bsx) * px + (y - bsy) * py) / pt * (pz / pt);
315 p = particleMomentum.P();
316 ctgTheta = 1.0 / TMath::Tan (particleMomentum.Theta());
317
318
319 // 3. time evaluation t = TMath::Min(t_r, t_z)
320 // t_r : time to exit from the sides
321 // t_z : time to exit from the front or the back
322 t_r = 0.0; // in [ns]
323 int sign_pz = (pz > 0.0) ? 1 : -1;
324 if(pz == 0.0) t_z = 1.0E99;
325 else t_z = gammam / (pz*1.0E9/c_light) * (-z + fHalfLength*sign_pz);
326
327 if(r_c + TMath::Abs(r) < fRadius)
328 {
329 // helix does not cross the cylinder sides
330 t = t_z;
331 }
332 else
333 {
334 asinrho = TMath::ASin((fRadius*fRadius - r_c*r_c - r*r) / (2*TMath::Abs(r)*r_c));
335 delta = phi_0 - phi;
336 if(delta <-TMath::Pi()) delta += 2*TMath::Pi();
337 if(delta > TMath::Pi()) delta -= 2*TMath::Pi();
338 t1 = (delta + asinrho) / omega;
339 t2 = (delta + TMath::Pi() - asinrho) / omega;
340 t3 = (delta + TMath::Pi() + asinrho) / omega;
341 t4 = (delta - asinrho) / omega;
342 t5 = (delta - TMath::Pi() - asinrho) / omega;
343 t6 = (delta - TMath::Pi() + asinrho) / omega;
344
345 if(t1 < 0.0) t1 = 1.0E99;
346 if(t2 < 0.0) t2 = 1.0E99;
347 if(t3 < 0.0) t3 = 1.0E99;
348 if(t4 < 0.0) t4 = 1.0E99;
349 if(t5 < 0.0) t5 = 1.0E99;
350 if(t6 < 0.0) t6 = 1.0E99;
351
352 t_ra = TMath::Min(t1, TMath::Min(t2, t3));
353 t_rb = TMath::Min(t4, TMath::Min(t5, t6));
354 t_r = TMath::Min(t_ra, t_rb);
355 t = TMath::Min(t_r, t_z);
356 }
357
358 // 4. position in terms of x(t), y(t), z(t)
359 x_t = x_c + r * TMath::Sin(omega * t - phi_0);
360 y_t = y_c + r * TMath::Cos(omega * t - phi_0);
361 z_t = z + pz*1.0E9 / c_light / gammam * t;
362 r_t = TMath::Hypot(x_t, y_t);
363
364
365 // compute path length for an helix
366
367 alpha = pz*1.0E9 / c_light / gammam;
368 l = t * TMath::Sqrt(alpha*alpha + r*r*omega*omega);
369
370 if(r_t > 0.0)
371 {
372
373 // store these variables before cloning
374 if(particle == candidate)
375 {
376 particle->D0 = d0*1.0E3;
377 particle->DZ = dz*1.0E3;
378 particle->P = p;
379 particle->PT = pt;
380 particle->CtgTheta = ctgTheta;
381 particle->Phi = phip;
382 }
383
384 mother = candidate;
385 candidate = static_cast<Candidate*>(candidate->Clone());
386
387 candidate->InitialPosition = particlePosition;
388 candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, particlePosition.T() + t*c_light*1.0E3);
389
390 candidate->Momentum = particleMomentum;
391
392 candidate->L = l*1.0E3;
393
394 candidate->Xd = xd*1.0E3;
395 candidate->Yd = yd*1.0E3;
396 candidate->Zd = zd*1.0E3;
397
398 candidate->AddCandidate(mother);
399
400 fOutputArray->Add(candidate);
401 switch(TMath::Abs(candidate->PID))
402 {
403 case 11:
404 fElectronOutputArray->Add(candidate);
405 break;
406 case 13:
407 fMuonOutputArray->Add(candidate);
408 break;
409 default:
410 fChargedHadronOutputArray->Add(candidate);
411 }
412 }
413 }
414 }
415}
416
417//------------------------------------------------------------------------------
418
Note: See TracBrowser for help on using the repository browser.