CMS 3D CMS Logo

ParticleManager.cc
Go to the documentation of this file.
2 
3 #include "HepMC/GenEvent.h"
4 #include "HepMC/Units.h"
5 #include "HepPDT/ParticleDataTable.hh"
6 
11 
16 
18  const HepMC::GenEvent & genEvent,
19  const HepPDT::ParticleDataTable & particleDataTable,
20  double beamPipeRadius,
21  double deltaRchargedMother,
22  const fastsim::ParticleFilter & particleFilter,
23  std::vector<SimTrack> & simTracks,
24  std::vector<SimVertex> & simVertices)
25  : genEvent_(&genEvent)
26  , genParticleIterator_(genEvent_->particles_begin())
27  , genParticleEnd_(genEvent_->particles_end())
28  , genParticleIndex_(1)
29  , particleDataTable_(&particleDataTable)
30  , beamPipeRadius2_(beamPipeRadius*beamPipeRadius)
31  , deltaRchargedMother_(deltaRchargedMother)
32  , particleFilter_(&particleFilter)
33  , simTracks_(&simTracks)
34  , simVertices_(&simVertices)
35  // prepare unit convsersions
36  // --------------------------------------------
37  // | | hepmc | cms |
38  // --------------------------------------------
39  // | length | genEvent_->length_unit | cm |
40  // | momentum | genEvent_->momentum_unit | GeV |
41  // | time | length unit (t*c) | ns |
42  // --------------------------------------------
43  , momentumUnitConversionFactor_(conversion_factor( genEvent_->momentum_unit(), HepMC::Units::GEV ))
44  , lengthUnitConversionFactor_(conversion_factor(genEvent_->length_unit(),HepMC::Units::LengthUnit::CM))
45  , lengthUnitConversionFactor2_(lengthUnitConversionFactor_*lengthUnitConversionFactor_)
46  , timeUnitConversionFactor_(lengthUnitConversionFactor_/fastsim::Constants::speedOfLight)
47 
48 {
49 
50  // add the main vertex from the signal event to the simvertex collection
51  if(genEvent.vertices_begin() != genEvent_->vertices_end())
52  {
53  const HepMC::FourVector & position = (*genEvent.vertices_begin())->position();
55  position.y()*lengthUnitConversionFactor_,
56  position.z()*lengthUnitConversionFactor_,
57  position.t()*timeUnitConversionFactor_)
58  ,-1);
59  }
60 }
61 
63 
65 {
66  std::unique_ptr<fastsim::Particle> particle;
67 
68  // retrieve particle from buffer
69  if(!particleBuffer_.empty())
70  {
71  particle = std::move(particleBuffer_.back());
72  particleBuffer_.pop_back();
73  }
74  // or from genParticle list
75  else
76  {
77  particle = nextGenParticle();
78  if(!particle) return nullptr;
79  }
80 
81  // if filter does not accept, skip particle
82  if(!particleFilter_->accepts(*particle))
83  {
84  return nextParticle(random);
85  }
86 
87  // lifetime or charge of particle are not yet set
88  if(!particle->remainingProperLifeTimeIsSet() || !particle->chargeIsSet())
89  {
90  // retrieve the particle data
91  const HepPDT::ParticleData * particleData = particleDataTable_->particle(HepPDT::ParticleID(particle->pdgId()));
92  if(!particleData)
93  {
94  throw cms::Exception("fastsim::ParticleManager") << "unknown pdg id: " << particle->pdgId() << std::endl;
95  }
96 
97  // set lifetime
98  if(!particle->remainingProperLifeTimeIsSet())
99  {
100  // The lifetime is 0. in the Pythia Particle Data Table! Calculate from width instead (ct=hbar/width).
101  // ct=particleData->lifetime().value();
102  double width = particleData->totalWidth().value();
103  if(width > 1.0e-35)
104  {
105  particle->setRemainingProperLifeTimeC(-log(random.flatShoot())*6.582119e-25/width/10.); // ct in cm
106  }
107  else
108  {
109  particle->setStable();
110  }
111  }
112 
113  // set charge
114  if(!particle->chargeIsSet())
115  {
116  particle->setCharge(particleData->charge());
117  }
118  }
119 
120  // add corresponding simTrack to simTrack collection
121  unsigned simTrackIndex = addSimTrack(particle.get());
122  particle->setSimTrackIndex(simTrackIndex);
123 
124  // and return
125  return particle;
126 }
127 
128 
130  const math::XYZTLorentzVector & vertexPosition,
131  int parentSimTrackIndex,
132  std::vector<std::unique_ptr<Particle> > & secondaries,
133  const SimplifiedGeometry * layer)
134 {
135 
136  // vertex must be within the accepted volume
137  if(!particleFilter_->acceptsVtx(vertexPosition))
138  {
139  return;
140  }
141 
142  // no need to create vertex in case no particles are produced
143  if(secondaries.empty()){
144  return;
145  }
146 
147  // add simVertex
148  unsigned simVertexIndex = addSimVertex(vertexPosition,parentSimTrackIndex);
149 
150  // closest charged daughter continues the track of the mother particle
151  // simplified tracking algorithm for fastSim
152  double distMin = 99999.;
153  int idx = -1;
154  int idxMin = -1;
155  for(auto & secondary : secondaries)
156  {
157  idx++;
158  if(secondary->getMotherDeltaR() != -1){
159  if(secondary->getMotherDeltaR() > deltaRchargedMother_){
160  // larger than max requirement on deltaR
161  secondary->resetMother();
162  }else{
163  if(secondary->getMotherDeltaR() < distMin){
164  distMin = secondary->getMotherDeltaR();
165  idxMin = idx;
166  }
167  }
168  }
169  }
170 
171  // add secondaries to buffer
172  idx = -1;
173  for(auto & secondary : secondaries)
174  {
175  idx++;
176  if(idxMin != -1){
177  // reset all but the particle with the lowest deltaR (which is at idxMin)
178  if(secondary->getMotherDeltaR() != -1 && idx != idxMin){
179  secondary->resetMother();
180  }
181  }
182 
183  // set origin vertex
184  secondary->setSimVertexIndex(simVertexIndex);
185  //
186  if(layer)
187  {
188  secondary->setOnLayer(layer->isForward(), layer->index());
189  }
190  // ...and add particle to buffer
191  particleBuffer_.push_back(std::move(secondary));
192  }
193 
194 }
195 
197 {
198  return this->addSimVertex(particle->position(), particle->simTrackIndex());
199 }
200 
203  int parentSimTrackIndex)
204 {
205  int simVertexIndex = simVertices_->size();
206  simVertices_->emplace_back(position.Vect(),
207  position.T(),
208  parentSimTrackIndex,
209  simVertexIndex);
210  return simVertexIndex;
211 }
212 
214 {
215  int simTrackIndex = simTracks_->size();
216  simTracks_->emplace_back(particle->pdgId(),
217  particle->momentum(),
218  particle->simVertexIndex(),
219  particle->genParticleIndex());
220  simTracks_->back().setTrackId(simTrackIndex);
221  return simTrackIndex;
222 }
223 
224 std::unique_ptr<fastsim::Particle> fastsim::ParticleManager::nextGenParticle()
225 {
226  // only consider particles that start in the beam pipe and end outside the beam pipe
227  // try to get the decay time from pythia
228  // use hepmc units
229  // make the link simtrack to simvertex
230  // try not to change the simvertex structure
231 
232  // loop over gen particles
234  {
235  // some handy pointers and references
236  const HepMC::GenParticle & particle = **genParticleIterator_;
237  const HepMC::GenVertex * productionVertex = particle.production_vertex();
238  const HepMC::GenVertex * endVertex = particle.end_vertex();
239 
240  // skip incoming particles
241  if(!productionVertex){
242  continue;
243  }
244 
245  // particle must be produced within the beampipe
246  if(productionVertex->position().perp2()*lengthUnitConversionFactor2_ > beamPipeRadius2_)
247  {
248  continue;
249  }
250 
251  // particle must not decay before it reaches the beam pipe
252  if(endVertex && endVertex->position().perp2()*lengthUnitConversionFactor2_ < beamPipeRadius2_)
253  {
254  continue;
255  }
256 
257  // make the particle
258  std::unique_ptr<Particle> newParticle(
259  new Particle(particle.pdg_id(),
260  math::XYZTLorentzVector(productionVertex->position().x()*lengthUnitConversionFactor_,
261  productionVertex->position().y()*lengthUnitConversionFactor_,
262  productionVertex->position().z()*lengthUnitConversionFactor_,
263  productionVertex->position().t()*timeUnitConversionFactor_),
265  particle.momentum().y()*momentumUnitConversionFactor_,
266  particle.momentum().z()*momentumUnitConversionFactor_,
267  particle.momentum().e()*momentumUnitConversionFactor_)));
268  newParticle->setGenParticleIndex(genParticleIndex_);
269 
270  // try to get the life time of the particle from the genEvent
271  if(endVertex)
272  {
273  double labFrameLifeTime = (endVertex->position().t() - productionVertex->position().t())*timeUnitConversionFactor_;
274  newParticle->setRemainingProperLifeTimeC(labFrameLifeTime / newParticle->gamma() * fastsim::Constants::speedOfLight);
275  }
276 
277  // Find production vertex if it already exists. Otherwise create new vertex
278  // Possible to recreate the whole GenEvent using SimTracks/SimVertices (see FBaseSimEvent::fill(..))
279  bool foundVtx = false;
280  for(const auto& simVtx : *simVertices_){
281  if(std::abs(simVtx.position().x() - newParticle->position().x()) < 1E-3 && std::abs(simVtx.position().y() - newParticle->position().y()) < 1E-3 && std::abs(simVtx.position().z() - newParticle->position().z()) < 1E-3){
282  newParticle->setSimVertexIndex(simVtx.vertexId());
283  foundVtx = true;
284  break;
285  }
286  }
287  if(!foundVtx) newParticle->setSimVertexIndex(addSimVertex(newParticle->position(), -1));
288 
289  // iterator/index has to be increased in case of return (is not done by the loop then)
291  // and return
292  return newParticle;
293  }
294 
295  return std::unique_ptr<Particle>();
296 }
Implementation of a generic detector layer (base class for forward/barrel layers).
const math::XYZTLorentzVector & position() const
Return position of the particle.
Definition: Particle.h:142
double flatShoot(double xmin=0.0, double xmax=1.0) const
HepPDT::ParticleDataTable ParticleDataTable
const double deltaRchargedMother_
For FastSim (cheat) tracking: cut on the angle between a charged mother and charged daughter...
bool accepts(const Particle &particle) const
Check all if all criteria are fullfilled.
const ParticleFilter *const particleFilter_
(Kinematic) cuts on the particles that have to be propagated.
int genParticleIndex() const
Return index of the particle in the genParticle vector.
Definition: Particle.h:167
double momentumUnitConversionFactor_
Convert pythia units to GeV (FastSim standard)
std::vector< std::unique_ptr< Particle > > particleBuffer_
The vector of all secondaries that are not yet propagated in the event.
static double constexpr speedOfLight
Speed of light [cm / ns].
Definition: Constants.h:16
virtual bool isForward() const =0
Returns false/true depending if the object is a (non-abstract) barrel/forward layer.
HepMC::GenEvent::particle_const_iterator genParticleIterator_
Iterator to keep track on which GenParticles where already considered.
int pdgId() const
Return pdgId of the particle.
Definition: Particle.h:136
TRandom random
Definition: MVATrainer.cc:138
const HepMC::GenEvent *const genEvent_
The GenEvent.
const HepMC::GenEvent::particle_const_iterator genParticleEnd_
The last particle of the GenEvent.
double lengthUnitConversionFactor2_
Convert pythia unis to cm^2 (FastSim standard)
unsigned addSimTrack(const Particle *particle)
Add a simTrack (simTrack contains some basic info about the particle, e.g. pdgId).
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
int simVertexIndex() const
Return index of the origin vertex.
Definition: Particle.h:161
const double beamPipeRadius2_
(Radius of the beampipe)^2
ParticleManager(const HepMC::GenEvent &genEvent, const HepPDT::ParticleDataTable &particleDataTable, double beamPipeRadius, double deltaRchargedMother, const ParticleFilter &particleFilter, std::vector< SimTrack > &simTracks, std::vector< SimVertex > &simVertices)
Constructor.
~ParticleManager()
Default destructor.
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
void addSecondaries(const math::XYZTLorentzVector &vertexPosition, int motherSimTrackId, std::vector< std::unique_ptr< Particle > > &secondaries, const SimplifiedGeometry *layer=nullptr)
Adds secondaries that are produced by any of the interactions (or particle decay) to the buffer...
std::unique_ptr< Particle > nextParticle(const RandomEngineAndDistribution &random)
Returns the next particle that has to be propagated (secondary or genParticle).
HepPDT::ParticleData ParticleData
int simTrackIndex() const
Return index of the SimTrack.
Definition: Particle.h:155
std::vector< SimTrack > * simTracks_
The generated SimTrack of this event.
std::unique_ptr< Particle > nextGenParticle()
Returns next particle from the GenEvent that has to be propagated.
double lengthUnitConversionFactor_
Convert pythia unis to cm (FastSim standard)
const HepPDT::ParticleDataTable *const particleDataTable_
Necessary to get information like lifetime and charge of a particle if unknown.
double timeUnitConversionFactor_
Convert pythia unis to ns (FastSim standard)
int index() const
Return index of this layer (layers are numbered according to their position in the detector)...
static int position[264][3]
Definition: ReadPGInfo.cc:509
std::vector< SimVertex > * simVertices_
The generated SimVertices of this event.
const math::XYZTLorentzVector & momentum() const
Return momentum of the particle.
Definition: Particle.h:145
(Kinematic) cuts on the particles that are propagated.
Definition of a generic FastSim Particle which can be propagated through the detector (formerly Parti...
Definition: Particle.h:19
int genParticleIndex_
Index of particle in the GenEvent (if it is a GenParticle)
unsigned addSimVertex(const math::XYZTLorentzVector &position, int motherIndex)
Add a simVertex (simVertex contains information about the track it was produced). ...
unsigned addEndVertex(const Particle *particle)
Necessary to add an end vertex to a particle.
def move(src, dest)
Definition: eostools.py:510
bool acceptsVtx(const math::XYZTLorentzVector &originVertexPosition) const
Vertex within tracker volume.