CMS 3D CMS Logo

HGCalTriggerNtupleGen.cc
Go to the documentation of this file.
6 
7 
10 
19 
24 
25 // NOTE: most of this code is borrowed by https://github.com/CMS-HGCAL/reco-ntuples
26 // kudos goes to the original authors. Ideally the 2 repos should be merged since they share part of the use case
27 #include <memory>
28 
29 namespace HGCal_helpers {
30 
31  class Coordinates {
32  public:
33  Coordinates() : x(0), y(0), z(0), eta(0), phi(0) {}
34  float x, y, z, eta, phi;
36  };
37 
38 
40  public:
42  : field_(f), prod_(field_, alongMomentum), absz_target_(0) {
43  ROOT::Math::SMatrixIdentity id;
45  const float uncert = 0.001;
46  C *= uncert;
48  }
49  void setPropagationTargetZ(const float &z);
50 
51  bool propagate(const double px, const double py, const double pz, const double x, const double y,
52  const double z, const float charge, Coordinates &coords) const;
53 
54  bool propagate(const math::XYZTLorentzVectorD &momentum, const math::XYZTLorentzVectorD &position,
55  const float charge, Coordinates &coords) const;
56 
57  private:
58  SimpleTrackPropagator() : field_(nullptr), prod_(field_, alongMomentum), absz_target_(0) {}
59  const RKPropagatorInS &RKProp() const { return prod_.propagator; }
64  float absz_target_;
65  };
66 
68  targetPlaneForward_ = Plane::build(Plane::PositionType(0, 0, std::abs(z)), Plane::RotationType());
69  targetPlaneBackward_ =
71  absz_target_ = std::abs(z);
72  }
73  bool SimpleTrackPropagator::propagate(const double px, const double py, const double pz,
74  const double x, const double y, const double z,
75  const float charge, Coordinates &output) const {
76  output = Coordinates();
77 
79  GlobalPoint startingPosition(x, y, z);
80  GlobalVector startingMomentum(px, py, pz);
81  Plane::PlanePointer startingPlane =
83  TSOS startingStateP(
84  GlobalTrajectoryParameters(startingPosition, startingMomentum, charge, field_), err_,
85  *startingPlane);
86 
87  TSOS trackStateP;
88  if (pz > 0) {
89  trackStateP = RKProp().propagate(startingStateP, *targetPlaneForward_);
90  } else {
91  trackStateP = RKProp().propagate(startingStateP, *targetPlaneBackward_);
92  }
93  if (trackStateP.isValid()) {
94  output.x = trackStateP.globalPosition().x();
95  output.y = trackStateP.globalPosition().y();
96  output.z = trackStateP.globalPosition().z();
97  output.phi = trackStateP.globalPosition().phi();
98  output.eta = trackStateP.globalPosition().eta();
99  return true;
100  }
101  return false;
102  }
103 
105  const math::XYZTLorentzVectorD &position, const float charge,
106  Coordinates &output) const {
107  return propagate(momentum.px(), momentum.py(), momentum.pz(), position.x(), position.y(),
108  position.z(), charge, output);
109  }
110 
111 } // HGCal_helpers
112 
113 
115 {
116 
117  public:
119 
120  void initialize(TTree&, const edm::ParameterSet&, edm::ConsumesCollector&&) final;
121  void fill(const edm::Event&, const edm::EventSetup& ) final;
122 
125 
126  enum ReachHGCal {
127  notReach = 0,
128  outsideEESurface = 1,
129  onEESurface = 2
130  };
131 
132  private:
133  void clear() final;
134 
135  edm::EDGetToken gen_token_;
136  edm::EDGetToken gen_PU_token_;
137 
138  int gen_n_;
139  int gen_PUNumInt_;
140  float gen_TrueNumInt_;
141 
142  float vtx_x_;
143  float vtx_y_;
144  float vtx_z_;
145 
146 
148  // GenParticles
149  //
150  std::vector<float> genpart_eta_;
151  std::vector<float> genpart_phi_;
152  std::vector<float> genpart_pt_;
153  std::vector<float> genpart_energy_;
154  std::vector<float> genpart_dvx_;
155  std::vector<float> genpart_dvy_;
156  std::vector<float> genpart_dvz_;
157  std::vector<float> genpart_ovx_;
158  std::vector<float> genpart_ovy_;
159  std::vector<float> genpart_ovz_;
160  std::vector<float> genpart_exx_;
161  std::vector<float> genpart_exy_;
162  std::vector<int> genpart_mother_;
163  std::vector<float> genpart_exphi_;
164  std::vector<float> genpart_exeta_;
165  std::vector<float> genpart_fbrem_;
166  std::vector<int> genpart_pid_;
167  std::vector<int> genpart_gen_;
168  std::vector<int> genpart_reachedEE_;
169  std::vector<bool> genpart_fromBeamPipe_;
170  std::vector<std::vector<float>> genpart_posx_;
171  std::vector<std::vector<float>> genpart_posy_;
172  std::vector<std::vector<float>> genpart_posz_;
173 
175  // reco::GenParticles
176  //
177  std::vector<float> gen_eta_;
178  std::vector<float> gen_phi_;
179  std::vector<float> gen_pt_;
180  std::vector<float> gen_energy_;
181  std::vector<int> gen_charge_;
182  std::vector<int> gen_pdgid_;
183  std::vector<int> gen_status_;
184  std::vector<std::vector<int>> gen_daughters_;
185 
186 
187  // -------convenient tool to deal with simulated tracks
188  std::unique_ptr<FSimEvent> mySimEvent_;
189 
190  // and also the magnetic field
191  const MagneticField *aField_;
192 
193  HGCalTriggerTools triggerTools_;
194 
195  // edm::EDGetTokenT<std::vector<reco::GenParticle> > genParticles_;
196  edm::EDGetToken simTracks_token_;
197  edm::EDGetToken simVertices_token_;
198  edm::EDGetToken hepmcev_token_;
199 
200  edm::ESWatcher<PDTRecord> pdt_watcher_;
202 
203 
204 };
205 
208  "HGCalTriggerNtupleGen" );
209 
210 
211 HGCalTriggerNtupleGen::
212 HGCalTriggerNtupleGen(const edm::ParameterSet& conf):HGCalTriggerNtupleBase(conf)
213 {
214 }
215 
216 
217 void
219 initialize(TTree& tree, const edm::ParameterSet& conf, edm::ConsumesCollector&& collector)
220 {
221 
222  edm::ParameterSet particleFilter_(conf.getParameter<edm::ParameterSet>("particleFilter"));
223  mySimEvent_ = std::make_unique<FSimEvent>(particleFilter_);
224 
225  gen_token_ = collector.consumes<reco::GenParticleCollection>(conf.getParameter<edm::InputTag>("GenParticles"));
226  gen_PU_token_ = collector.consumes<std::vector<PileupSummaryInfo>>(conf.getParameter<edm::InputTag>("GenPU"));
227  tree.Branch("gen_n", &gen_n_, "gen_n/I");
228  tree.Branch("gen_PUNumInt", &gen_PUNumInt_ ,"gen_PUNumInt/I");
229  tree.Branch("gen_TrueNumInt", &gen_TrueNumInt_ ,"gen_TrueNumInt/F");
230 
231  hepmcev_token_ = collector.consumes<edm::HepMCProduct>(conf.getParameter<edm::InputTag>("MCEvent"));
232 
233  simTracks_token_ = collector.consumes<std::vector<SimTrack>>(conf.getParameter<edm::InputTag>("SimTracks"));
234  simVertices_token_ = collector.consumes<std::vector<SimVertex>>(conf.getParameter<edm::InputTag>("SimVertices"));
235 
236  tree.Branch("vtx_x", &vtx_x_);
237  tree.Branch("vtx_y", &vtx_y_);
238  tree.Branch("vtx_z", &vtx_z_);
239 
240 
241  tree.Branch("gen_eta", &gen_eta_);
242  tree.Branch("gen_phi", &gen_phi_);
243  tree.Branch("gen_pt", &gen_pt_);
244  tree.Branch("gen_energy", &gen_energy_);
245  tree.Branch("gen_charge", &gen_charge_);
246  tree.Branch("gen_pdgid", &gen_pdgid_);
247  tree.Branch("gen_status", &gen_status_);
248  tree.Branch("gen_daughters", &gen_daughters_);
249 
250  tree.Branch("genpart_eta", &genpart_eta_);
251  tree.Branch("genpart_phi", &genpart_phi_);
252  tree.Branch("genpart_pt", &genpart_pt_);
253  tree.Branch("genpart_energy", &genpart_energy_);
254  tree.Branch("genpart_dvx", &genpart_dvx_);
255  tree.Branch("genpart_dvy", &genpart_dvy_);
256  tree.Branch("genpart_dvz", &genpart_dvz_);
257  tree.Branch("genpart_ovx", &genpart_ovx_);
258  tree.Branch("genpart_ovy", &genpart_ovy_);
259  tree.Branch("genpart_ovz", &genpart_ovz_);
260  tree.Branch("genpart_mother", &genpart_mother_);
261  tree.Branch("genpart_exphi", &genpart_exphi_);
262  tree.Branch("genpart_exeta", &genpart_exeta_);
263  tree.Branch("genpart_exx", &genpart_exx_);
264  tree.Branch("genpart_exy", &genpart_exy_);
265  tree.Branch("genpart_fbrem", &genpart_fbrem_);
266  tree.Branch("genpart_pid", &genpart_pid_);
267  tree.Branch("genpart_gen", &genpart_gen_);
268  tree.Branch("genpart_reachedEE", &genpart_reachedEE_);
269  tree.Branch("genpart_fromBeamPipe", &genpart_fromBeamPipe_);
270  tree.Branch("genpart_posx", &genpart_posx_);
271  tree.Branch("genpart_posy", &genpart_posy_);
272  tree.Branch("genpart_posz", &genpart_posz_);
273 
274 
275 }
276 
277 void
280 {
281  clear();
282 
284  iEvent.getByToken(gen_PU_token_, PupInfo_h);
285  const std::vector< PileupSummaryInfo >& PupInfo = *PupInfo_h;
286 
287  if(pdt_watcher_.check(es))
288  {
290  es.get<PDTRecord>().get(pdt);
291  mySimEvent_->initializePdt(&(*pdt));
292  }
293 
294  if(magfield_watcher_.check(es))
295  {
297  es.get<IdealMagneticFieldRecord>().get(magfield);
298  aField_ = &(*magfield);
299  }
300 
301  triggerTools_.eventSetup(es);
302 
303  // This balck magic is needed to use the mySimEvent_
304  ParticleTable::Sentry ptable(mySimEvent_->theTable());
306  edm::Handle<std::vector<SimTrack>> simTracksHandle;
307  edm::Handle<std::vector<SimVertex>> simVerticesHandle;
308 
309  iEvent.getByToken(hepmcev_token_, hevH);
310  iEvent.getByToken(simTracks_token_, simTracksHandle);
311  iEvent.getByToken(simVertices_token_, simVerticesHandle);
312  mySimEvent_->fill(*simTracksHandle, *simVerticesHandle);
313 
314  HepMC::GenVertex *primaryVertex = *(hevH)->GetEvent()->vertices_begin();
315  const float mm2cm = 0.1;
316  vtx_x_ = primaryVertex->position().x() * mm2cm; // to put in official units
317  vtx_y_ = primaryVertex->position().y() * mm2cm;
318  vtx_z_ = primaryVertex->position().z() * mm2cm;
319  Point sim_pv(vtx_x_, vtx_y_, vtx_z_);
320 
321 
322  HGCal_helpers::SimpleTrackPropagator toHGCalPropagator(aField_);
323  toHGCalPropagator.setPropagationTargetZ(triggerTools_.getLayerZ(1));
324  std::vector<FSimTrack *> allselectedgentracks;
325  const float eeInnerRadius = 25.;
326  const float eeOuterRadius = 160.;
327  unsigned int npart = mySimEvent_->nTracks();
328  for (unsigned int i = 0; i < npart; ++i) {
329  std::vector<float> xp, yp, zp;
330  FSimTrack &myTrack(mySimEvent_->track(i));
331  math::XYZTLorentzVectorD vtx(0, 0, 0, 0);
332 
333  int reachedEE = ReachHGCal::notReach; // compute the extrapolations for the particles reaching EE
334  // and for the gen particles
335  double fbrem = -1;
336 
337  if (std::abs(myTrack.vertex().position().z()) >= triggerTools_.getLayerZ(1)) continue;
338 
339  const unsigned nlayers = triggerTools_.lastLayerBH();
340  if (myTrack.noEndVertex()) // || myTrack.genpartIndex()>=0)
341  {
342  HGCal_helpers::Coordinates propcoords;
343  bool reachesHGCal = toHGCalPropagator.propagate(
344  myTrack.momentum(), myTrack.vertex().position(), myTrack.charge(), propcoords);
345  vtx = propcoords.toVector();
346 
347  if (reachesHGCal && vtx.Rho() < eeOuterRadius && vtx.Rho() > eeInnerRadius) {
348  reachedEE = ReachHGCal::onEESurface;
349  double dpt = 0;
350 
351  for (int i = 0; i < myTrack.nDaughters(); ++i) dpt += myTrack.daughter(i).momentum().pt();
352  if (abs(myTrack.type()) == 11) fbrem = dpt / myTrack.momentum().pt();
353  } else if (reachesHGCal && vtx.Rho() > eeOuterRadius)
354  reachedEE = ReachHGCal::outsideEESurface;
355 
356  HGCal_helpers::SimpleTrackPropagator indiv_particleProp(aField_);
357  for (unsigned il = 1; il <= nlayers; ++il) {
358  const float charge = myTrack.charge();
359  indiv_particleProp.setPropagationTargetZ(triggerTools_.getLayerZ(il));
360  HGCal_helpers::Coordinates propCoords;
361  indiv_particleProp.propagate(myTrack.momentum(), myTrack.vertex().position(), charge,
362  propCoords);
363 
364  xp.push_back(propCoords.x);
365  yp.push_back(propCoords.y);
366  zp.push_back(propCoords.z);
367  }
368  } else {
369  vtx = myTrack.endVertex().position();
370  }
371  auto orig_vtx = myTrack.vertex().position();
372 
373  allselectedgentracks.push_back(&mySimEvent_->track(i));
374  // fill branches
375  genpart_eta_.push_back(myTrack.momentum().eta());
376  genpart_phi_.push_back(myTrack.momentum().phi());
377  genpart_pt_.push_back(myTrack.momentum().pt());
378  genpart_energy_.push_back(myTrack.momentum().energy());
379  genpart_dvx_.push_back(vtx.x());
380  genpart_dvy_.push_back(vtx.y());
381  genpart_dvz_.push_back(vtx.z());
382 
383  genpart_ovx_.push_back(orig_vtx.x());
384  genpart_ovy_.push_back(orig_vtx.y());
385  genpart_ovz_.push_back(orig_vtx.z());
386 
387  HGCal_helpers::Coordinates hitsHGCal;
388  toHGCalPropagator.propagate(myTrack.momentum(), orig_vtx, myTrack.charge(), hitsHGCal);
389 
390  genpart_exphi_.push_back(hitsHGCal.phi);
391  genpart_exeta_.push_back(hitsHGCal.eta);
392  genpart_exx_.push_back(hitsHGCal.x);
393  genpart_exy_.push_back(hitsHGCal.y);
394 
395  genpart_fbrem_.push_back(fbrem);
396  genpart_pid_.push_back(myTrack.type());
397  genpart_gen_.push_back(myTrack.genpartIndex());
398  genpart_reachedEE_.push_back(reachedEE);
399  genpart_fromBeamPipe_.push_back(true);
400 
401  genpart_posx_.push_back(xp);
402  genpart_posy_.push_back(yp);
403  genpart_posz_.push_back(zp);
404  }
405 
406 
407  edm::Handle<std::vector<reco::GenParticle>> genParticlesHandle;
408  iEvent.getByToken(gen_token_, genParticlesHandle);
409  gen_n_ = genParticlesHandle->size();
410 
411  for (const auto& particle : *genParticlesHandle) {
412  gen_eta_.push_back(particle.eta());
413  gen_phi_.push_back(particle.phi());
414  gen_pt_.push_back(particle.pt());
415  gen_energy_.push_back(particle.energy());
416  gen_charge_.push_back(particle.charge());
417  gen_pdgid_.push_back(particle.pdgId());
418  gen_status_.push_back(particle.status());
419  std::vector<int> daughters(particle.daughterRefVector().size(), 0);
420  for (unsigned j = 0; j < particle.daughterRefVector().size(); ++j) {
421  daughters[j] = static_cast<int>(particle.daughterRefVector().at(j).key());
422  }
423  gen_daughters_.push_back(daughters);
424  }
425 
426 
427  // associate gen particles to mothers
428  genpart_mother_.resize(genpart_posz_.size(), -1);
429  for (size_t i = 0; i < allselectedgentracks.size(); i++) {
430  const auto tracki = allselectedgentracks.at(i);
431 
432  for (size_t j = i + 1; j < allselectedgentracks.size(); j++) {
433  const auto trackj = allselectedgentracks.at(j);
434 
435  if (!tracki->noMother()) {
436  if (&tracki->mother() == trackj) genpart_mother_.at(i) = j;
437  }
438  if (!trackj->noMother()) {
439  if (&trackj->mother() == tracki) genpart_mother_.at(j) = i;
440  }
441  }
442  }
443 
444  for(const auto& PVI : PupInfo)
445  {
446  if(PVI.getBunchCrossing() == 0)
447  {
448  gen_PUNumInt_ = PVI.getPU_NumInteractions();
449  gen_TrueNumInt_ = PVI.getTrueNumInteractions();
450  }
451  }
452 
453 }
454 
455 
456 void
459 {
460 
461  gen_n_ = 0;
462  gen_PUNumInt_ = 0;
463  gen_TrueNumInt_ = 0.;
464 
465  vtx_x_ = 0;
466  vtx_y_ = 0;
467  vtx_z_ = 0;
468 
469 
470  //
471  genpart_eta_.clear();
472  genpart_phi_.clear();
473  genpart_pt_.clear();
474  genpart_energy_.clear();
475  genpart_dvx_.clear();
476  genpart_dvy_.clear();
477  genpart_dvz_.clear();
478  genpart_ovx_.clear();
479  genpart_ovy_.clear();
480  genpart_ovz_.clear();
481  genpart_exx_.clear();
482  genpart_exy_.clear();
483  genpart_mother_.clear();
484  genpart_exphi_.clear();
485  genpart_exeta_.clear();
486  genpart_fbrem_.clear();
487  genpart_pid_.clear();
488  genpart_gen_.clear();
489  genpart_reachedEE_.clear();
490  genpart_fromBeamPipe_.clear();
491  genpart_posx_.clear();
492  genpart_posy_.clear();
493  genpart_posz_.clear();
494 
496  // reco::GenParticles
497  //
498  gen_eta_.clear();
499  gen_phi_.clear();
500  gen_pt_.clear();
501  gen_energy_.clear();
502  gen_charge_.clear();
503  gen_pdgid_.clear();
504  gen_status_.clear();
505  gen_daughters_.clear();
506 }
std::vector< GenParticle > GenParticleCollection
collection of GenParticles
T getParameter(std::string const &) const
static AlgebraicMatrix initialize()
PositionVector3D< Cartesian3D< double >, DefaultCoordinateSystemTag > Point
Definition: Transform3DPJ.h:67
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:14
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:579
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
T y() const
Definition: PV3DBase.h:63
GlobalPoint globalPosition() const
double npart
Definition: HydjetWrapper.h:49
ROOT::Math::SMatrix< double, 5, 5, ROOT::Math::MatRepSym< double, 5 > > AlgebraicSymMatrix55
#define nullptr
void initialize(TTree &, const edm::ParameterSet &, edm::ConsumesCollector &&) final
int iEvent
Definition: GenABIO.cc:230
void clear(CLHEP::HepGenMatrix &m)
Helper function: Reset all elements of a matrix to 0.
Definition: matutil.cc:167
static PlanePointer build(Args &&...args)
Definition: Plane.h:33
T z() const
Definition: PV3DBase.h:64
void fill(const edm::Event &, const edm::EventSetup &) final
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double f[11][100]
const RKPropagatorInS & RKProp() const
ROOT::Math::Transform3DPJ Transform3D
TrajectoryStateOnSurface TSOS
Definition: TestHits.cc:19
bool propagate(const double px, const double py, const double pz, const double x, const double y, const double z, const float charge, Coordinates &coords) const
T eta() const
Definition: PV3DBase.h:76
HLT enums.
math::XYZTLorentzVectorD toVector()
static int position[264][3]
Definition: ReadPGInfo.cc:509
T get() const
Definition: EventSetup.h:63
#define DEFINE_EDM_PLUGIN(factory, type, name)
ROOT::Math::Transform3DPJ::Point Point
Definition: tree.py:1
T x() const
Definition: PV3DBase.h:62