Fork me on GitHub

source: svn/trunk/modules/ParticlePropagator.cc

Last change on this file was 1367, checked in by Michele Selvaggi, 10 years ago

track counting btagging
  • Property svn:keywords set to Id Revision Date
File size: 8.7 KB
Line 
1
2/** \class ParticlePropagator
3 *
4 * Propagates charged and neutral particles
5 * from a given vertex to a cylinder defined by its radius,
6 * its half-length, centered at (0,0,0) and with its axis
7 * oriented along the z-axis.
8 *
9 * $Date: 2014-04-16 14:08:33 +0000 (Wed, 16 Apr 2014) $
10 * $Revision: 1367 $
11 *
12 *
13 * \author P. Demin - UCL, Louvain-la-Neuve
14 *
15 */
16
17#include "modules/ParticlePropagator.h"
18
19#include "classes/DelphesClasses.h"
20#include "classes/DelphesFactory.h"
21#include "classes/DelphesFormula.h"
22
23#include "ExRootAnalysis/ExRootResult.h"
24#include "ExRootAnalysis/ExRootFilter.h"
25#include "ExRootAnalysis/ExRootClassifier.h"
26
27#include "TMath.h"
28#include "TString.h"
29#include "TFormula.h"
30#include "TRandom3.h"
31#include "TObjArray.h"
32#include "TDatabasePDG.h"
33#include "TLorentzVector.h"
34
35#include <algorithm>
36#include <stdexcept>
37#include <iostream>
38#include <sstream>
39
40using namespace std;
41
42//------------------------------------------------------------------------------
43
44ParticlePropagator::ParticlePropagator() :
45 fItInputArray(0)
46{
47}
48
49//------------------------------------------------------------------------------
50
51ParticlePropagator::~ParticlePropagator()
52{
53}
54
55//------------------------------------------------------------------------------
56
57void ParticlePropagator::Init()
58{
59 fRadius = GetDouble("Radius", 1.0);
60 fRadius2 = fRadius*fRadius;
61 fHalfLength = GetDouble("HalfLength", 3.0);
62 fBz = GetDouble("Bz", 0.0);
63 if(fRadius < 1.0E-2)
64 {
65 cout << "ERROR: magnetic field radius is too low\n";
66 return;
67 }
68 if(fHalfLength < 1.0E-2)
69 {
70 cout << "ERROR: magnetic field length is too low\n";
71 return;
72 }
73
74 // import array with output from filter/classifier module
75
76 fInputArray = ImportArray(GetString("InputArray", "Delphes/stableParticles"));
77 fItInputArray = fInputArray->MakeIterator();
78
79 // create output arrays
80
81 fOutputArray = ExportArray(GetString("OutputArray", "stableParticles"));
82 fChargedHadronOutputArray = ExportArray(GetString("ChargedHadronOutputArray", "chargedHadrons"));
83 fElectronOutputArray = ExportArray(GetString("ElectronOutputArray", "electrons"));
84 fMuonOutputArray = ExportArray(GetString("MuonOutputArray", "muons"));
85}
86
87//------------------------------------------------------------------------------
88
89void ParticlePropagator::Finish()
90{
91 if(fItInputArray) delete fItInputArray;
92}
93
94//------------------------------------------------------------------------------
95
96void ParticlePropagator::Process()
97{
98 Candidate *candidate, *mother;
99 TLorentzVector candidatePosition, candidateMomentum;
100 Double_t px, py, pz, pt, pt2, e, q;
101 Double_t x, y, z, t, r, phi;
102 Double_t x_c, y_c, r_c, phi_c, phi_0;
103 Double_t x_t, y_t, z_t, r_t;
104 Double_t t1, t2, t3, t4, t5, t6;
105 Double_t t_z, t_r, t_ra, t_rb;
106 Double_t tmp, discr, discr2;
107 Double_t delta, gammam, omega, asinrho;
108 Double_t ang_mom, rcu, rc2, dxy, xd, yd, zd;
109
110 const Double_t c_light = 2.99792458E8;
111
112 fItInputArray->Reset();
113 while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
114 {
115 candidatePosition = candidate->Position;
116 candidateMomentum = candidate->Momentum;
117 x = candidatePosition.X()*1.0E-3;
118 y = candidatePosition.Y()*1.0E-3;
119 z = candidatePosition.Z()*1.0E-3;
120 q = candidate->Charge;
121
122 // check that particle position is inside the cylinder
123 if(TMath::Hypot(x, y) > fRadius || TMath::Abs(z) > fHalfLength)
124 {
125 continue;
126 }
127
128 px = candidateMomentum.Px();
129 py = candidateMomentum.Py();
130 pz = candidateMomentum.Pz();
131 pt = candidateMomentum.Pt();
132 pt2 = candidateMomentum.Perp2();
133 e = candidateMomentum.E();
134
135 if(pt2 < 1.0E-9)
136 {
137 continue;
138 }
139
140 if(TMath::Abs(q) < 1.0E-9 || TMath::Abs(fBz) < 1.0E-9)
141 {
142 // solve pt2*t^2 + 2*(px*x + py*y)*t + (fRadius2 - x*x - y*y) = 0
143 tmp = px*y - py*x;
144 discr2 = pt2*fRadius2 - tmp*tmp;
145
146 if(discr2 < 0)
147 {
148 // no solutions
149 continue;
150 }
151
152 tmp = px*x + py*y;
153 discr = TMath::Sqrt(discr2);
154 t1 = (-tmp + discr)/pt2;
155 t2 = (-tmp - discr)/pt2;
156 t = (t1 < 0) ? t2 : t1;
157
158 z_t = z + pz*t;
159 if(TMath::Abs(z_t) > fHalfLength)
160 {
161 t3 = (+fHalfLength - z) / pz;
162 t4 = (-fHalfLength - z) / pz;
163 t = (t3 < 0) ? t4 : t3;
164 }
165
166 x_t = x + px*t;
167 y_t = y + py*t;
168 z_t = z + pz*t;
169
170 mother = candidate;
171 candidate = static_cast<Candidate*>(candidate->Clone());
172
173 candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t*e*1.0E3);
174
175 candidate->Momentum = candidateMomentum;
176 candidate->AddCandidate(mother);
177
178 fOutputArray->Add(candidate);
179 if(TMath::Abs(q) > 1.0E-9)
180 {
181 switch(TMath::Abs(candidate->PID))
182 {
183 case 11:
184 fElectronOutputArray->Add(candidate);
185 break;
186 case 13:
187 fMuonOutputArray->Add(candidate);
188 break;
189 default:
190 fChargedHadronOutputArray->Add(candidate);
191 }
192 }
193 }
194 else
195 {
196
197 // 1. initial transverse momentum p_{T0} : Part->pt
198 // initial transverse momentum direction \phi_0 = -atan(p_X0/p_Y0)
199 // relativistic gamma : gamma = E/mc² ; gammam = gamma \times m
200 // giration frequency \omega = q/(gamma m) fBz
201 // helix radius r = p_T0 / (omega gamma m)
202
203 gammam = e*1.0E9 / (c_light*c_light); // gammam in [eV/c²]
204 omega = q * fBz / (gammam); // omega is here in [ 89875518 / s]
205 r = pt / (q * fBz) * 1.0E9/c_light; // in [m]
206
207 phi_0 = TMath::ATan2(py, px); // [rad] in [-pi; pi]
208
209 // 2. helix axis coordinates
210 x_c = x + r*TMath::Sin(phi_0);
211 y_c = y - r*TMath::Cos(phi_0);
212 r_c = TMath::Hypot(x_c, y_c);
213 phi_c = TMath::ATan2(y_c, x_c);
214 phi = phi_c;
215 if(x_c < 0.0) phi += TMath::Pi();
216
217 rcu = TMath::Abs(r);
218 rc2 = r_c*r_c;
219
220 // calculate coordinates of closest approach to track circle in transverse plane xd, yd, zd
221 xd = x_c*x_c*x_c - x_c*rcu*r_c + x_c*y_c*y_c;
222 xd = ( rc2 > 0.0 ) ? xd / rc2 : -999;
223 yd = y_c*(-rcu*r_c + rc2);
224 yd = ( rc2 > 0.0 ) ? yd / rc2 : -999;
225 zd = z + (TMath::Sqrt(xd*xd+yd*yd) - TMath::Sqrt(x*x+y*y))*pz/pt;
226
227 // calculate impact paramater
228 ang_mom = (xd*py - yd*px);
229 dxy = ang_mom/pt;
230
231
232 // 3. time evaluation t = TMath::Min(t_r, t_z)
233 // t_r : time to exit from the sides
234 // t_z : time to exit from the front or the back
235 t_r = 0.0; // in [ns]
236 int sign_pz = (pz > 0.0) ? 1 : -1;
237 if(pz == 0.0) t_z = 1.0E99;
238 else t_z = gammam / (pz*1.0E9/c_light) * (-z + fHalfLength*sign_pz);
239
240 if(r_c + TMath::Abs(r) < fRadius)
241 {
242 // helix does not cross the cylinder sides
243 t = t_z;
244 }
245 else
246 {
247 asinrho = TMath::ASin( (fRadius*fRadius - r_c*r_c - r*r) / (2*TMath::Abs(r)*r_c) );
248 delta = phi_0 - phi;
249 if(delta <-TMath::Pi()) delta += 2*TMath::Pi();
250 if(delta > TMath::Pi()) delta -= 2*TMath::Pi();
251 t1 = (delta + asinrho) / omega;
252 t2 = (delta + TMath::Pi() - asinrho) / omega;
253 t3 = (delta + TMath::Pi() + asinrho) / omega;
254 t4 = (delta - asinrho) / omega;
255 t5 = (delta - TMath::Pi() - asinrho) / omega;
256 t6 = (delta - TMath::Pi() + asinrho) / omega;
257
258 if(t1 < 0) t1 = 1.0E99;
259 if(t2 < 0) t2 = 1.0E99;
260 if(t3 < 0) t3 = 1.0E99;
261 if(t4 < 0) t4 = 1.0E99;
262 if(t5 < 0) t5 = 1.0E99;
263 if(t6 < 0) t6 = 1.0E99;
264
265 t_ra = TMath::Min(t1, TMath::Min(t2, t3));
266 t_rb = TMath::Min(t4, TMath::Min(t5, t6));
267 t_r = TMath::Min(t_ra, t_rb);
268 t = TMath::Min(t_r, t_z);
269 }
270
271 // 4. position in terms of x(t), y(t), z(t)
272 x_t = x_c + r * TMath::Sin(omega * t - phi_0);
273 y_t = y_c + r * TMath::Cos(omega * t - phi_0);
274 z_t = z + pz*1.0E9 / c_light / gammam * t;
275 r_t = TMath::Hypot(x_t, y_t);
276
277 if(r_t > 0.0)
278 {
279 mother = candidate;
280 candidate = static_cast<Candidate*>(candidate->Clone());
281
282 candidate->Position.SetXYZT(x_t*1.0E3, y_t*1.0E3, z_t*1.0E3, candidatePosition.T() + t*c_light*1.0E3);
283
284 candidate->Momentum = candidateMomentum;
285 candidate->Xd = xd*1.0E3;
286 candidate->Yd = yd*1.0E3;
287 candidate->Zd = zd*1.0E3;
288
289 candidate->AddCandidate(mother);
290
291 fOutputArray->Add(candidate);
292 switch(TMath::Abs(candidate->PID))
293 {
294 case 11:
295 fElectronOutputArray->Add(candidate);
296 break;
297 case 13:
298 fMuonOutputArray->Add(candidate);
299 break;
300 default:
301 fChargedHadronOutputArray->Add(candidate);
302 }
303 }
304 }
305 }
306}
307
308//------------------------------------------------------------------------------
309
Note: See TracBrowser for help on using the repository browser.