CMS 3D CMS Logo

PFSimParticleProducer.cc
Go to the documentation of this file.
1 
59 
60 #include <memory>
61 #include <set>
62 #include <sstream>
63 #include <string>
64 
66 public:
68 
69  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
70 
71  void produce(edm::Event&, const edm::EventSetup&) override;
72 
74  void getSimIDs(const TrackHandle& trackh, std::vector<unsigned>& recTrackSimID);
75 
76 private:
81 
82  //MC Truth Matching
83  //modif-beg
87  //modif-end
88 
95 
97 
98  // parameters for retrieving true particles information --
99 
101  std::unique_ptr<FSimEvent> mySimEvent;
102 
103  // flags for the various tasks ---------------------------
104 
107 
109  bool verbose_;
110 };
111 
114 
116  // particleFlowSimParticle
118  desc.addUntracked<edm::InputTag>("fastSimProducer", edm::InputTag("fastSimProducer", "EcalHitsEB"));
119  desc.addUntracked<bool>("MCTruthMatchingInfo", false);
120  desc.add<edm::InputTag>("RecTracks", edm::InputTag("trackerDrivenElectronSeeds"));
121  desc.add<std::string>("Fitter", "KFFittingSmoother");
122  desc.add<edm::InputTag>("ecalRecHitsEE", edm::InputTag("caloRecHits", "EcalRecHitsEE"));
123  desc.add<edm::InputTag>("ecalRecHitsEB", edm::InputTag("caloRecHits", "EcalRecHitsEB"));
124  desc.addUntracked<bool>("process_RecTracks", false);
125  desc.add<edm::ParameterSetDescription>("ParticleFilter", {});
126  desc.add<std::string>("TTRHBuilder", "WithTrackAngle");
127  desc.addUntracked<bool>("process_Particles", true);
128  desc.add<std::string>("Propagator", "PropagatorWithMaterial");
129  desc.add<edm::InputTag>("sim", edm::InputTag("g4SimHits"));
130  desc.addUntracked<bool>("verbose", false);
131  descriptions.add("particleFlowSimParticle", desc);
132 }
133 
134 using namespace std;
135 using namespace edm;
136 
138  processParticles_ = iConfig.getUntrackedParameter<bool>("process_Particles", true);
139 
140  inputTagSim_ = iConfig.getParameter<InputTag>("sim");
141  tokenSim_ = consumes<std::vector<SimTrack> >(inputTagSim_);
142  tokenSimVertices_ = consumes<std::vector<SimVertex> >(inputTagSim_);
143 
144  //retrieving collections for MC Truth Matching
145 
146  //modif-beg
147  inputTagFastSimProducer_ = iConfig.getUntrackedParameter<InputTag>("fastSimProducer");
148  tokenFastSimProducer_ = consumes<edm::PCaloHitContainer>(inputTagFastSimProducer_);
149  mctruthMatchingInfo_ = iConfig.getUntrackedParameter<bool>("MCTruthMatchingInfo", false);
150  //modif-end
151 
152  inputTagRecTracks_ = iConfig.getParameter<InputTag>("RecTracks");
153  tokenRecTracks_ = consumes<reco::PFRecTrackCollection>(inputTagRecTracks_);
154 
155  inputTagEcalRecHitsEB_ = iConfig.getParameter<InputTag>("ecalRecHitsEB");
156  tokenEcalRecHitsEB_ = consumes<EcalRecHitCollection>(inputTagEcalRecHitsEB_);
157  inputTagEcalRecHitsEE_ = iConfig.getParameter<InputTag>("ecalRecHitsEE");
158  tokenEcalRecHitsEE_ = consumes<EcalRecHitCollection>(inputTagEcalRecHitsEE_);
159 
160  pdtToken_ = esConsumes();
161 
162  verbose_ = iConfig.getUntrackedParameter<bool>("verbose", false);
163 
164  // register products
165  produces<reco::PFSimParticleCollection>();
166 
167  particleFilter_ = iConfig.getParameter<ParameterSet>("ParticleFilter");
168 
169  mySimEvent = std::make_unique<FSimEvent>(particleFilter_);
170 }
171 
173  // init Particle data table (from Pythia)
174  mySimEvent->initializePdt(&iSetup.getData(pdtToken_));
175 
176  LogDebug("PFSimParticleProducer") << "START event: " << iEvent.id().event() << " in run " << iEvent.id().run()
177  << endl;
178 
179  //MC Truth Matching only with Famos and UnFoldedMode option to true!!
180 
181  //vector to store the trackIDs of rectracks corresponding
182  //to the simulated particle.
183  std::vector<unsigned> recTrackSimID;
184 
185  //In order to know which simparticule contribute to
186  //a given Ecal RecHit energy, we need to access
187  //the PCAloHit from FastSim.
188 
189  typedef std::pair<double, unsigned> hitSimID;
190  std::vector<std::list<hitSimID> > caloHitsEBID(62000);
191  std::vector<double> caloHitsEBTotE(62000, 0.0);
192 
193  if (mctruthMatchingInfo_) {
194  //getting the PCAloHit
195  auto pcalohits = iEvent.getHandle(tokenFastSimProducer_);
196 
197  if (!pcalohits) {
198  ostringstream err;
199  err << "could not find pcaloHit "
200  << "fastSimProducer:EcalHitsEB";
201  LogError("PFSimParticleProducer") << err.str() << endl;
202 
203  throw cms::Exception("MissingProduct", err.str());
204  } else {
205  assert(pcalohits.isValid());
206 
207  edm::PCaloHitContainer::const_iterator it = pcalohits.product()->begin();
208  edm::PCaloHitContainer::const_iterator itend = pcalohits.product()->end();
209 
210  //loop on the PCaloHit from FastSim Calorimetry
211  for (; it != itend; ++it) {
212  EBDetId detid(it->id());
213 
214  if (it->energy() > 0.0) {
215  std::pair<double, unsigned> phitsimid = make_pair(it->energy(), it->geantTrackId());
216  caloHitsEBID[detid.hashedIndex()].push_back(phitsimid);
217  caloHitsEBTotE[detid.hashedIndex()] += it->energy(); //summing pcalhit energy
218  } //energy > 0
219 
220  } //loop PcaloHits
221  } //pcalohit handle access
222 
223  //Retrieving the PFRecTrack collection for
224  //Monte Carlo Truth Matching tool
226  try {
227  LogDebug("PFSimParticleProducer") << "getting PFRecTracks" << endl;
228  iEvent.getByToken(tokenRecTracks_, recTracks);
229 
230  } catch (cms::Exception& err) {
231  LogError("PFSimParticleProducer") << err << " cannot get collection "
232  << "particleFlowBlock"
233  << ":"
234  << "" << endl;
235  } //pfrectrack handle access
236 
237  //getting the simID corresponding to
238  //each of the PFRecTracks
239  getSimIDs(recTracks, recTrackSimID);
240 
241  } //mctruthMatchingInfo_ //modif
242 
243  // deal with true particles
244  if (processParticles_) {
245  auto pOutputPFSimParticleCollection = std::make_unique<reco::PFSimParticleCollection>();
246 
248  bool found = iEvent.getByToken(tokenSim_, simTracks);
249  if (!found) {
250  ostringstream err;
251  err << "cannot find sim tracks: " << inputTagSim_;
252  LogError("PFSimParticleProducer") << err.str() << endl;
253 
254  throw cms::Exception("MissingProduct", err.str());
255  }
256 
258  found = iEvent.getByToken(tokenSimVertices_, simVertices);
259  if (!found) {
260  LogError("PFSimParticleProducer") << "cannot find sim vertices: " << inputTagSim_ << endl;
261  return;
262  }
263 
264  mySimEvent->fill(*simTracks, *simVertices);
265 
266  if (verbose_)
267  mySimEvent->print();
268 
269  for (unsigned i = 0; i < mySimEvent->nTracks(); i++) {
270  const FSimTrack& fst = mySimEvent->track(i);
271 
272  int motherId = -1;
273  if (!fst.noMother()) // a mother exist
274  motherId = fst.mother().id();
275 
276  //This is finding out the simID corresponding
277  //to the recTrack
278 
279  //GETTING THE TRACK ID
280  unsigned recTrackID = 99999;
281  vector<unsigned> recHitContrib; //modif
282  vector<double> recHitContribFrac; //modif
283 
284  if (mctruthMatchingInfo_) { //modif
285 
286  for (unsigned lo = 0; lo < recTrackSimID.size(); lo++) {
287  if (i == recTrackSimID[lo]) {
288  recTrackID = lo;
289  } //match track
290  } //loop rectrack
291 
292  // get the ecalBarrel rechits for MC truth matching tool
294  bool found = iEvent.getByToken(tokenEcalRecHitsEB_, rhcHandle);
295  if (!found) {
296  ostringstream err;
297  err << "could not find rechits " << inputTagEcalRecHitsEB_;
298  LogError("PFSimParticleProducer") << err.str() << endl;
299 
300  throw cms::Exception("MissingProduct", err.str());
301  } else {
302  assert(rhcHandle.isValid());
303 
304  EBRecHitCollection::const_iterator it_rh = rhcHandle.product()->begin();
305  EBRecHitCollection::const_iterator itend_rh = rhcHandle.product()->end();
306 
307  for (; it_rh != itend_rh; ++it_rh) {
308  unsigned rhit_hi = EBDetId(it_rh->id()).hashedIndex();
309  EBDetId detid(it_rh->id());
310 
311  auto it_phit = caloHitsEBID[rhit_hi].begin();
312  auto itend_phit = caloHitsEBID[rhit_hi].end();
313  for (; it_phit != itend_phit; ++it_phit) {
314  if (i == it_phit->second) {
315  //Alex (08/10/08) TO BE REMOVED, eliminating
316  //duplicated rechits
317  bool alreadyin = false;
318  for (unsigned ihit = 0; ihit < recHitContrib.size(); ++ihit)
319  if (detid.rawId() == recHitContrib[ihit])
320  alreadyin = true;
321 
322  if (!alreadyin) {
323  double pcalofraction = 0.0;
324  if (caloHitsEBTotE[rhit_hi] != 0.0)
325  pcalofraction = (it_phit->first / caloHitsEBTotE[rhit_hi]) * 100.0;
326 
327  //store info
328  recHitContrib.push_back(it_rh->id());
329  recHitContribFrac.push_back(pcalofraction);
330  } //selected rechits
331  } //matching
332  } //loop pcalohit
333 
334  } //loop rechits
335 
336  } //getting the rechits
337 
338  } //mctruthMatchingInfo_ //modif
339 
340  reco::PFSimParticle particle(
341  fst.charge(), fst.type(), fst.id(), motherId, fst.daughters(), recTrackID, recHitContrib, recHitContribFrac);
342 
343  const FSimVertex& originVtx = fst.vertex();
344 
345  math::XYZPoint posOrig(originVtx.position().x(), originVtx.position().y(), originVtx.position().z());
346 
347  math::XYZTLorentzVector momOrig(
348  fst.momentum().px(), fst.momentum().py(), fst.momentum().pz(), fst.momentum().e());
349  reco::PFTrajectoryPoint pointOrig(-1, reco::PFTrajectoryPoint::ClosestApproach, posOrig, momOrig);
350 
351  // point 0 is origin vertex
352  particle.addPoint(pointOrig);
353 
354  if (!fst.noEndVertex()) {
355  const FSimVertex& endVtx = fst.endVertex();
356 
357  math::XYZPoint posEnd(endVtx.position().x(), endVtx.position().y(), endVtx.position().z());
358 
360 
362 
363  particle.addPoint(pointEnd);
364  } else { // add a dummy point
366  particle.addPoint(dummy);
367  }
368 
369  if (fst.onLayer1()) { // PS layer1
370  const RawParticle& rp = fst.layer1Entrance();
371 
372  math::XYZPoint posLayer1(rp.x(), rp.y(), rp.z());
373  math::XYZTLorentzVector momLayer1(rp.px(), rp.py(), rp.pz(), rp.e());
374  reco::PFTrajectoryPoint layer1Pt(-1, reco::PFTrajectoryPoint::PS1, posLayer1, momLayer1);
375 
376  particle.addPoint(layer1Pt);
377 
378  // extrapolate to cluster depth
379  } else { // add a dummy point
381  particle.addPoint(dummy);
382  }
383 
384  if (fst.onLayer2()) { // PS layer2
385  const RawParticle& rp = fst.layer2Entrance();
386 
387  math::XYZPoint posLayer2(rp.x(), rp.y(), rp.z());
388  math::XYZTLorentzVector momLayer2(rp.px(), rp.py(), rp.pz(), rp.e());
389  reco::PFTrajectoryPoint layer2Pt(-1, reco::PFTrajectoryPoint::PS2, posLayer2, momLayer2);
390 
391  particle.addPoint(layer2Pt);
392 
393  // extrapolate to cluster depth
394  } else { // add a dummy point
396  particle.addPoint(dummy);
397  }
398 
399  if (fst.onEcal()) {
400  const RawParticle& rp = fst.ecalEntrance();
401 
402  math::XYZPoint posECAL(rp.x(), rp.y(), rp.z());
403  math::XYZTLorentzVector momECAL(rp.px(), rp.py(), rp.pz(), rp.e());
405 
406  particle.addPoint(ecalPt);
407 
408  // extrapolate to cluster depth
409  } else { // add a dummy point
411  particle.addPoint(dummy);
412  }
413 
414  // add a dummy point for ECAL Shower max
416  particle.addPoint(dummy);
417 
418  if (fst.onHcal()) {
419  const RawParticle& rpin = fst.hcalEntrance();
420 
421  math::XYZPoint posHCALin(rpin.x(), rpin.y(), rpin.z());
422  math::XYZTLorentzVector momHCALin(rpin.px(), rpin.py(), rpin.pz(), rpin.e());
423  reco::PFTrajectoryPoint hcalPtin(-1, reco::PFTrajectoryPoint::HCALEntrance, posHCALin, momHCALin);
424 
425  particle.addPoint(hcalPtin);
426 
427  const RawParticle& rpout = fst.hcalExit();
428 
429  math::XYZPoint posHCALout(rpout.x(), rpout.y(), rpout.z());
430  math::XYZTLorentzVector momHCALout(rpout.px(), rpout.py(), rpout.pz(), rpout.e());
431  reco::PFTrajectoryPoint hcalPtout(-1, reco::PFTrajectoryPoint::HCALExit, posHCALout, momHCALout);
432 
433  particle.addPoint(hcalPtout);
434 
435  const RawParticle& rpho = fst.hoEntrance();
436 
437  math::XYZPoint posHOEntrance(rpho.x(), rpho.y(), rpho.z());
438  math::XYZTLorentzVector momHOEntrance(rpho.px(), rpho.py(), rpho.pz(), rpho.e());
439  reco::PFTrajectoryPoint hoPtin(-1, reco::PFTrajectoryPoint::HOLayer, posHOEntrance, momHOEntrance);
440 
441  particle.addPoint(hoPtin);
442 
443  } else { // add a dummy point
445  particle.addPoint(dummy);
446  }
447 
448  pOutputPFSimParticleCollection->push_back(particle);
449  }
450 
451  iEvent.put(std::move(pOutputPFSimParticleCollection));
452  }
453 
454  LogDebug("PFSimParticleProducer") << "STOP event: " << iEvent.id().event() << " in run " << iEvent.id().run() << endl;
455 }
456 
457 void PFSimParticleProducer::getSimIDs(const TrackHandle& trackh, std::vector<unsigned>& recTrackSimID) {
458  if (trackh.isValid()) {
459  for (unsigned i = 0; i < trackh->size(); i++) {
460  const reco::PFRecTrackRef ref(trackh, i);
461 
462  for (auto const& hit : ref->trackRef()->recHits()) {
463  if (hit->isValid()) {
464  auto rechit = dynamic_cast<const FastTrackerRecHit*>(hit);
465 
466  for (unsigned int st_index = 0; st_index < rechit->nSimTrackIds(); ++st_index) {
467  recTrackSimID.push_back(rechit->simTrackId(st_index));
468  }
469  break;
470  }
471  } //loop track rechit
472  } //loop recTracks
473  } //track handle valid
474 }
PFSimParticleProducer(const edm::ParameterSet &)
bool noMother() const
no mother particle
std::unique_ptr< FSimEvent > mySimEvent
bool processParticles_
process particles on/off
edm::InputTag inputTagEcalRecHitsEE_
ESGetTokenH3DDVariant esConsumes(std::string const &Record, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
double x() const
x of vertex
Definition: RawParticle.h:282
edm::EDGetTokenT< EcalRecHitCollection > tokenEcalRecHitsEE_
T const & getData(const ESGetToken< T, R > &iToken) const noexcept(false)
Definition: EventSetup.h:119
const RawParticle & hcalEntrance() const
The particle at HCAL entrance.
Definition: FSimTrack.h:146
const RawParticle & layer2Entrance() const
The particle at Preshower Layer 2.
Definition: FSimTrack.h:140
edm::Handle< reco::PFRecTrackCollection > TrackHandle
const RawParticle & ecalEntrance() const
The particle at ECAL entrance.
Definition: FSimTrack.h:143
edm::ParameterSet particleFilter_
T const * product() const
Definition: Handle.h:70
std::vector< EcalRecHit >::const_iterator const_iterator
edm::InputTag inputTagEcalRecHitsEB_
int onHcal() const
Definition: FSimTrack.h:116
const RawParticle & hoEntrance() const
The particle at HCAL exir.
Definition: FSimTrack.h:155
Log< level::Error, false > LogError
bool noEndVertex() const
no end vertex
assert(be >=bs)
Producer for PFRecTracks and PFSimParticles.
const FSimVertex vertex() const
Origin vertex.
edm::EDGetTokenT< reco::PFRecTrackCollection > tokenRecTracks_
double px() const
x of the momentum
Definition: RawParticle.h:296
const FSimVertex & endVertex() const
end vertex
T getUntrackedParameter(std::string const &, T const &) const
void produce(edm::Event &, const edm::EventSetup &) override
int type() const
particle type (HEP PDT convension)
Definition: CoreSimTrack.h:22
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
int hashedIndex(int ieta, int iphi)
Definition: EcalPyUtils.cc:36
int iEvent
Definition: GenABIO.cc:224
const XYZTLorentzVector & momentum() const
Temporary (until move of SimTrack to Mathcore) - No! Actually very useful.
Definition: FSimTrack.h:209
edm::InputTag inputTagSim_
module label for retrieving input simtrack and simvertex
const RawParticle & hcalExit() const
The particle at HCAL exir.
Definition: FSimTrack.h:152
Point of closest approach from beam axis (initial point in the case of PFSimParticle) ...
edm::EDGetTokenT< std::vector< SimTrack > > tokenSim_
float charge() const
charge
Definition: FSimTrack.h:56
edm::ESGetToken< HepPDT::ParticleDataTable, edm::DefaultRecord > pdtToken_
double z() const
z of vertex
Definition: RawParticle.h:284
const std::vector< int > & daughters() const
Vector of daughter indices.
int onLayer2() const
Definition: FSimTrack.h:106
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
double py() const
y of the momentum
Definition: RawParticle.h:299
int onLayer1() const
Definition: FSimTrack.h:101
const_iterator begin() const
double y() const
y of vertex
Definition: RawParticle.h:283
true particle for particle flow
Definition: PFSimParticle.h:19
int onEcal() const
Definition: FSimTrack.h:111
double pz() const
z of the momentum
Definition: RawParticle.h:302
edm::EDGetTokenT< edm::PCaloHitContainer > tokenFastSimProducer_
const_iterator end() const
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
edm::EDGetTokenT< EcalRecHitCollection > tokenEcalRecHitsEB_
const math::XYZTLorentzVector & position() const
Temporary (until CMSSW moves to Mathcore) - No ! Actually very useful.
Definition: FSimVertex.h:48
void add(std::string const &label, ParameterSetDescription const &psetDescription)
bool isValid() const
Definition: HandleBase.h:70
HLT enums.
A PFTrack holds several trajectory points, which basically contain the position and momentum of a tra...
const RawParticle & layer1Entrance() const
The particle at Preshower Layer 1.
Definition: FSimTrack.h:137
double e() const
energy of the momentum
Definition: RawParticle.h:305
void getSimIDs(const TrackHandle &trackh, std::vector< unsigned > &recTrackSimID)
def move(src, dest)
Definition: eostools.py:511
int id() const
the index in FBaseSimEvent and other vectors
Definition: FSimTrack.h:96
edm::EDGetTokenT< std::vector< SimVertex > > tokenSimVertices_
const FSimTrack & mother() const
mother
edm::InputTag inputTagFastSimProducer_
#define LogDebug(id)