Fork me on GitHub

source: git/modules/ParticlePropagator.cc@ 308ab1a

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 308ab1a was 764f12a0, checked in by Pavel Demin <pavel.demin@…>, 10 years ago

set Dxy in ParticlePropagator (close #330)

  • Property mode set to 100644
File size: 9.4 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
71void ParticlePropagator::Init()
72{
73 fRadius = GetDouble("Radius", 1.0);
74 fRadius2 = fRadius*fRadius;
75 fHalfLength = GetDouble("HalfLength", 3.0);
76 fBz = GetDouble("Bz", 0.0);
77 if(fRadius < 1.0E-2)
78 {
79 cout << "ERROR: magnetic field radius is too low\n";
80 return;
81 }
82 if(fHalfLength < 1.0E-2)
83 {
84 cout << "ERROR: magnetic field length is too low\n";
85 return;
86 }
87
88 // import array with output from filter/classifier module
89
90 fInputArray = ImportArray(GetString("InputArray", "Delphes/stableParticles"));
91 fItInputArray = fInputArray->MakeIterator();
92
93 // create output arrays
94
95 fOutputArray = ExportArray(GetString("OutputArray", "stableParticles"));
96 fChargedHadronOutputArray = ExportArray(GetString("ChargedHadronOutputArray", "chargedHadrons"));
97 fElectronOutputArray = ExportArray(GetString("ElectronOutputArray", "electrons"));
98 fMuonOutputArray = ExportArray(GetString("MuonOutputArray", "muons"));
99}
100
101//------------------------------------------------------------------------------
102
103void ParticlePropagator::Finish()
104{
105 if(fItInputArray) delete fItInputArray;
106}
107
108//------------------------------------------------------------------------------
109
110void ParticlePropagator::Process()
111{
112 Candidate *candidate, *mother;
113 TLorentzVector candidatePosition, candidateMomentum;
114 Double_t px, py, pz, pt, pt2, e, q;
115 Double_t x, y, z, t, r, phi;
116 Double_t x_c, y_c, r_c, phi_c, phi_0;
117 Double_t x_t, y_t, z_t, r_t;
118 Double_t t1, t2, t3, t4, t5, t6;
119 Double_t t_z, t_r, t_ra, t_rb;
120 Double_t tmp, discr, discr2;
121 Double_t delta, gammam, omega, asinrho;
122 Double_t rcu, rc2, dxy, xd, yd, zd;
123
124 const Double_t c_light = 2.99792458E8;
125
126 fItInputArray->Reset();
127 while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
128 {
129 candidatePosition = candidate->Position;
130 candidateMomentum = candidate->Momentum;
131 x = candidatePosition.X()*1.0E-3;
132 y = candidatePosition.Y()*1.0E-3;
133 z = candidatePosition.Z()*1.0E-3;
134 q = candidate->Charge;
135
136 // check that particle position is inside the cylinder
137 if(TMath::Hypot(x, y) > fRadius || TMath::Abs(z) > fHalfLength)
138 {
139 continue;
140 }
141
142 px = candidateMomentum.Px();
143 py = candidateMomentum.Py();
144 pz = candidateMomentum.Pz();
145 pt = candidateMomentum.Pt();
146 pt2 = candidateMomentum.Perp2();
147 e = candidateMomentum.E();
148
149 if(pt2 < 1.0E-9)
150 {
151 continue;
152 }
153
154 if(TMath::Abs(q) < 1.0E-9 || TMath::Abs(fBz) < 1.0E-9)
155 {
156 // solve pt2*t^2 + 2*(px*x + py*y)*t - (fRadius2 - x*x - y*y) = 0
157 tmp = px*y - py*x;
158 discr2 = pt2*fRadius2 - tmp*tmp;
159
160 if(discr2 < 0.0)
161 {
162 // no solutions
163 continue;
164 }
165
166 tmp = px*x + py*y;
167 discr = TMath::Sqrt(discr2);
168 t1 = (-tmp + discr)/pt2;
169 t2 = (-tmp - discr)/pt2;
170 t = (t1 < 0.0) ? t2 : t1;
171
172 z_t = z + pz*t;
173 if(TMath::Abs(z_t) > fHalfLength)
174 {
175 t3 = (+fHalfLength - z) / pz;
176 t4 = (-fHalfLength - z) / pz;
177 t = (t3 < 0.0) ? t4 : t3;
178 }
179
180 x_t = x + px*t;
181 y_t = y + py*t;
182 z_t = z + pz*t;
183
184 mother = candidate;
185 candidate = static_cast<Candidate*>(candidate->Clone());
186
187 candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t*e*1.0E3);
188
189 candidate->Momentum = candidateMomentum;
190 candidate->AddCandidate(mother);
191
192 fOutputArray->Add(candidate);
193 if(TMath::Abs(q) > 1.0E-9)
194 {
195 switch(TMath::Abs(candidate->PID))
196 {
197 case 11:
198 fElectronOutputArray->Add(candidate);
199 break;
200 case 13:
201 fMuonOutputArray->Add(candidate);
202 break;
203 default:
204 fChargedHadronOutputArray->Add(candidate);
205 }
206 }
207 }
208 else
209 {
210
211 // 1. initial transverse momentum p_{T0}: Part->pt
212 // initial transverse momentum direction phi_0 = -atan(p_X0/p_Y0)
213 // relativistic gamma: gamma = E/mc^2; gammam = gamma * m
214 // gyration frequency omega = q/(gamma m) fBz
215 // helix radius r = p_{T0} / (omega gamma m)
216
217 gammam = e*1.0E9 / (c_light*c_light); // gammam in [eV/c^2]
218 omega = q * fBz / (gammam); // omega is here in [89875518/s]
219 r = pt / (q * fBz) * 1.0E9/c_light; // in [m]
220
221 phi_0 = TMath::ATan2(py, px); // [rad] in [-pi, pi]
222
223 // 2. helix axis coordinates
224 x_c = x + r*TMath::Sin(phi_0);
225 y_c = y - r*TMath::Cos(phi_0);
226 r_c = TMath::Hypot(x_c, y_c);
227 phi_c = TMath::ATan2(y_c, x_c);
228 phi = phi_c;
229 if(x_c < 0.0) phi += TMath::Pi();
230
231 rcu = TMath::Abs(r);
232 rc2 = r_c*r_c;
233
234 // calculate coordinates of closest approach to track circle in transverse plane xd, yd, zd
235 xd = x_c*x_c*x_c - x_c*rcu*r_c + x_c*y_c*y_c;
236 xd = (rc2 > 0.0) ? xd / rc2 : -999;
237 yd = y_c*(-rcu*r_c + rc2);
238 yd = (rc2 > 0.0) ? yd / rc2 : -999;
239 zd = z + (TMath::Sqrt(xd*xd + yd*yd) - TMath::Sqrt(x*x + y*y))*pz/pt;
240
241 // calculate impact paramater
242 dxy = (xd*py - yd*px)/pt;
243
244 // 3. time evaluation t = TMath::Min(t_r, t_z)
245 // t_r : time to exit from the sides
246 // t_z : time to exit from the front or the back
247 t_r = 0.0; // in [ns]
248 int sign_pz = (pz > 0.0) ? 1 : -1;
249 if(pz == 0.0) t_z = 1.0E99;
250 else t_z = gammam / (pz*1.0E9/c_light) * (-z + fHalfLength*sign_pz);
251
252 if(r_c + TMath::Abs(r) < fRadius)
253 {
254 // helix does not cross the cylinder sides
255 t = t_z;
256 }
257 else
258 {
259 asinrho = TMath::ASin( (fRadius*fRadius - r_c*r_c - r*r) / (2*TMath::Abs(r)*r_c) );
260 delta = phi_0 - phi;
261 if(delta <-TMath::Pi()) delta += 2*TMath::Pi();
262 if(delta > TMath::Pi()) delta -= 2*TMath::Pi();
263 t1 = (delta + asinrho) / omega;
264 t2 = (delta + TMath::Pi() - asinrho) / omega;
265 t3 = (delta + TMath::Pi() + asinrho) / omega;
266 t4 = (delta - asinrho) / omega;
267 t5 = (delta - TMath::Pi() - asinrho) / omega;
268 t6 = (delta - TMath::Pi() + asinrho) / omega;
269
270 if(t1 < 0.0) t1 = 1.0E99;
271 if(t2 < 0.0) t2 = 1.0E99;
272 if(t3 < 0.0) t3 = 1.0E99;
273 if(t4 < 0.0) t4 = 1.0E99;
274 if(t5 < 0.0) t5 = 1.0E99;
275 if(t6 < 0.0) t6 = 1.0E99;
276
277 t_ra = TMath::Min(t1, TMath::Min(t2, t3));
278 t_rb = TMath::Min(t4, TMath::Min(t5, t6));
279 t_r = TMath::Min(t_ra, t_rb);
280 t = TMath::Min(t_r, t_z);
281 }
282
283 // 4. position in terms of x(t), y(t), z(t)
284 x_t = x_c + r * TMath::Sin(omega * t - phi_0);
285 y_t = y_c + r * TMath::Cos(omega * t - phi_0);
286 z_t = z + pz*1.0E9 / c_light / gammam * t;
287 r_t = TMath::Hypot(x_t, y_t);
288
289 if(r_t > 0.0)
290 {
291 mother = candidate;
292 candidate = static_cast<Candidate*>(candidate->Clone());
293
294 candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t*c_light*1.0E3);
295
296 candidate->Momentum = candidateMomentum;
297 candidate->Dxy = dxy*1.0E3;
298 candidate->Xd = xd*1.0E3;
299 candidate->Yd = yd*1.0E3;
300 candidate->Zd = zd*1.0E3;
301
302 candidate->AddCandidate(mother);
303
304 fOutputArray->Add(candidate);
305 switch(TMath::Abs(candidate->PID))
306 {
307 case 11:
308 fElectronOutputArray->Add(candidate);
309 break;
310 case 13:
311 fMuonOutputArray->Add(candidate);
312 break;
313 default:
314 fChargedHadronOutputArray->Add(candidate);
315 }
316 }
317 }
318 }
319}
320
321//------------------------------------------------------------------------------
322
Note: See TracBrowser for help on using the repository browser.