Fork me on GitHub

source: git/modules/ParticlePropagator.cc@ d108fdc

ImprovedOutputFile Timing dual_readout llp
Last change on this file since d108fdc was 9330b7b, checked in by Pavel Demin <pavel-demin@…>, 8 years ago

add RadiusMax and HalfLengthMax to ParticlePropagator

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