Fork me on GitHub

source: git/modules/TimeOfFlight.cc@ 71eba5e

Last change on this file since 71eba5e was 5cc021e, checked in by michele <michele.selvaggi@…>, 3 years ago

removed trailing comment

  • Property mode set to 100644
File size: 6.7 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 /** \class TimeOfFlight
20 *
21 * Calculates Time-Of-Flight
22 *
23 * \author Michele Selvaggi - CERN
24 *
25 */
26
27#include "modules/TimeOfFlight.h"
28
29#include "classes/DelphesClasses.h"
30#include "classes/DelphesFactory.h"
31#include "classes/DelphesFormula.h"
32
33#include "ExRootAnalysis/ExRootClassifier.h"
34#include "ExRootAnalysis/ExRootFilter.h"
35#include "ExRootAnalysis/ExRootResult.h"
36
37#include "TDatabasePDG.h"
38#include "TFormula.h"
39#include "TLorentzVector.h"
40#include "TMath.h"
41#include "TObjArray.h"
42#include "TRandom3.h"
43#include "TString.h"
44
45#include <algorithm>
46#include <iostream>
47#include <sstream>
48#include <stdexcept>
49
50using namespace std;
51//------------------------------------------------------------------------------
52
53TimeOfFlight::TimeOfFlight() :
54 fItTrackInputArray(0), fItVertexInputArray(0)
55{
56}
57
58//------------------------------------------------------------------------------
59
60TimeOfFlight::~TimeOfFlight()
61{
62}
63
64//------------------------------------------------------------------------------
65
66void TimeOfFlight::Init()
67{
68
69 // method to compute vertex time
70 fVertexTimeMode = GetInt("VertexTimeMode", 0);
71
72 // import track input array
73 fTrackInputArray = ImportArray(GetString("TrackInputArray", "MuonMomentumSmearing/muons"));
74 fItTrackInputArray = fTrackInputArray->MakeIterator();
75
76 // import vertex input array
77 fVertexInputArray = ImportArray(GetString("VertexInputArray", "TruthVertexFinder/vertices"));
78 fItVertexInputArray = fVertexInputArray->MakeIterator();
79
80 // create output array
81 fOutputArray = ExportArray(GetString("OutputArray", "tracks"));
82}
83
84//------------------------------------------------------------------------------
85
86void TimeOfFlight::Finish()
87{
88 if(fItTrackInputArray) delete fItTrackInputArray;
89 if(fItVertexInputArray) delete fItVertexInputArray;
90}
91
92//------------------------------------------------------------------------------
93
94void TimeOfFlight::Process()
95{
96 Candidate *candidate, *particle, *vertex, *constituent, *mother;
97 Double_t ti, t_truth, tf;
98 Double_t l, tof, beta, p, mass;
99
100 const Double_t c_light = 2.99792458E8;
101
102 // first compute momenta of vertices based on reconstructed tracks
103 ComputeVertexMomenta();
104
105 fItTrackInputArray->Reset();
106 while((candidate = static_cast<Candidate *>(fItTrackInputArray->Next())))
107 {
108
109 particle = static_cast<Candidate *>(candidate->GetCandidates()->At(0));
110
111 const TLorentzVector &candidateInitialPosition = particle->Position;
112 const TLorentzVector &candidateInitialPositionSmeared = candidate->InitialPosition;
113 const TLorentzVector &candidateFinalPosition = candidate->Position;
114 const TLorentzVector &candidateMomentum = particle->Momentum;
115
116 // time at vertex from MC truth
117 t_truth = candidateInitialPosition.T() * 1.0E-3 / c_light;
118
119 if (candidate->Position.Vect().Mag() < 5.) continue;
120
121 // various options on how to calculate the vertex time
122 ti=0;
123 switch (fVertexTimeMode)
124 {
125 case 0:
126 {
127 // assume ti from MC truth
128 // most aggressive, we are cheating and assume we can perfectly reconstruct time of primary and secondary vertices
129 ti = t_truth;
130 break;
131 }
132 case 1:
133 {
134 // always assume t=0, most conservative assumption
135 // reasonable assumption for particles originating from PV, if beamSpot has small time spread compared to timing resolution
136 // probably bad assumption for particles from highly displaced vertices (i.e Ks)
137 ti=0;
138 break;
139 }
140 case 2:
141 {
142 // same as 2 but attempt at estimate beta from vertex mass and momentum
143 beta = 1.;
144 fItVertexInputArray->Reset();
145 while((vertex = static_cast<Candidate *>(fItVertexInputArray->Next())))
146 {
147 TIter itGenParts(vertex->GetCandidates());
148 itGenParts.Reset();
149
150 while((constituent = static_cast<Candidate *>(itGenParts.Next())))
151 {
152 if (particle == constituent)
153 {
154 beta = vertex->Momentum.Beta();
155 break;
156 }
157 }
158 } // end vertex loop
159
160 // track displacement to be possibily replaced by vertex fitted position
161 ti = candidateInitialPositionSmeared.Vect().Mag() * 1.0E-3 /(beta*c_light);
162 }
163 break;
164 }
165
166 p = candidateMomentum.P();
167
168 // this quantity has already been smeared by another module
169 tf = candidateFinalPosition.T() * 1.0E-3 / c_light;
170
171 // calculate time-of-flight
172 tof = tf - ti;
173 // path length of the full helix
174 l = candidate->L * 1.0E-3;
175
176 // particle velocity
177 beta = l/(c_light*tof);
178
179 // calculate particle mass (i.e particle ID)
180 mass = 0.;
181 if (beta<1) mass = p* TMath::Sqrt(1/(beta*beta) - 1);
182
183 mother = candidate;
184 candidate = static_cast<Candidate*>(candidate->Clone());
185
186 // update time at vertex based on option
187 candidate->InitialPosition.SetT(ti * 1.0E3 * c_light);
188
189 // update particle mass based on TOF-based PID (commented for now, assume this calculation is done offline)
190 //candidate->Momentum.SetVectM(candidateMomentum.Vect(), mass);
191
192 candidate->AddCandidate(mother);
193 fOutputArray->Add(candidate);
194 }
195}
196
197//------------------------------------------------------------------------------
198
199void TimeOfFlight::ComputeVertexMomenta()
200{
201 Candidate *track, *constituent, *particle, *vertex;
202
203 fItVertexInputArray->Reset();
204 while((vertex = static_cast<Candidate *>(fItVertexInputArray->Next())))
205 {
206
207 TIter itGenParts(vertex->GetCandidates());
208 itGenParts.Reset();
209
210 while((constituent = static_cast<Candidate *>(itGenParts.Next())))
211 {
212 fItTrackInputArray->Reset();
213 while((track = static_cast<Candidate *>(fItTrackInputArray->Next())))
214 {
215 // get gen part that generated track
216 particle = static_cast<Candidate *>(track->GetCandidates()->At(0));
217 if (particle == constituent)
218 {
219 vertex->Momentum += track->Momentum;
220 }
221
222 } // end track loop
223 } // end vertex consitutent loop
224 } // end vertex loop
225}
Note: See TracBrowser for help on using the repository browser.