Fork me on GitHub

source: git/modules/ParticlePropagator.cc@ df35033

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

remove svn tags and fix formatting

  • 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 ang_mom, 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)
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) ? 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) ? 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² ; gammam = gamma \times m
214 // giration 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²]
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 ang_mom = (xd*py - yd*px);
243 dxy = ang_mom/pt;
244
245
246 // 3. time evaluation t = TMath::Min(t_r, t_z)
247 // t_r : time to exit from the sides
248 // t_z : time to exit from the front or the back
249 t_r = 0.0; // in [ns]
250 int sign_pz = (pz > 0.0) ? 1 : -1;
251 if(pz == 0.0) t_z = 1.0E99;
252 else t_z = gammam / (pz*1.0E9/c_light) * (-z + fHalfLength*sign_pz);
253
254 if(r_c + TMath::Abs(r) < fRadius)
255 {
256 // helix does not cross the cylinder sides
257 t = t_z;
258 }
259 else
260 {
261 asinrho = TMath::ASin( (fRadius*fRadius - r_c*r_c - r*r) / (2*TMath::Abs(r)*r_c) );
262 delta = phi_0 - phi;
263 if(delta <-TMath::Pi()) delta += 2*TMath::Pi();
264 if(delta > TMath::Pi()) delta -= 2*TMath::Pi();
265 t1 = (delta + asinrho) / omega;
266 t2 = (delta + TMath::Pi() - asinrho) / omega;
267 t3 = (delta + TMath::Pi() + asinrho) / omega;
268 t4 = (delta - asinrho) / omega;
269 t5 = (delta - TMath::Pi() - asinrho) / omega;
270 t6 = (delta - TMath::Pi() + asinrho) / omega;
271
272 if(t1 < 0) t1 = 1.0E99;
273 if(t2 < 0) t2 = 1.0E99;
274 if(t3 < 0) t3 = 1.0E99;
275 if(t4 < 0) t4 = 1.0E99;
276 if(t5 < 0) t5 = 1.0E99;
277 if(t6 < 0) t6 = 1.0E99;
278
279 t_ra = TMath::Min(t1, TMath::Min(t2, t3));
280 t_rb = TMath::Min(t4, TMath::Min(t5, t6));
281 t_r = TMath::Min(t_ra, t_rb);
282 t = TMath::Min(t_r, t_z);
283 }
284
285 // 4. position in terms of x(t), y(t), z(t)
286 x_t = x_c + r * TMath::Sin(omega * t - phi_0);
287 y_t = y_c + r * TMath::Cos(omega * t - phi_0);
288 z_t = z + pz*1.0E9 / c_light / gammam * t;
289 r_t = TMath::Hypot(x_t, y_t);
290
291 if(r_t > 0.0)
292 {
293 mother = candidate;
294 candidate = static_cast<Candidate*>(candidate->Clone());
295
296 candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t*c_light*1.0E3);
297
298 candidate->Momentum = candidateMomentum;
299 candidate->Xd = xd*1.0E3;
300 candidate->Yd = yd*1.0E3;
301 candidate->Zd = zd*1.0E3;
302
303 candidate->AddCandidate(mother);
304
305 fOutputArray->Add(candidate);
306 switch(TMath::Abs(candidate->PID))
307 {
308 case 11:
309 fElectronOutputArray->Add(candidate);
310 break;
311 case 13:
312 fMuonOutputArray->Add(candidate);
313 break;
314 default:
315 fChargedHadronOutputArray->Add(candidate);
316 }
317 }
318 }
319 }
320}
321
322//------------------------------------------------------------------------------
323
Note: See TracBrowser for help on using the repository browser.