CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PFSimParticleProducer.cc
Go to the documentation of this file.
2 
5 
8 
11 
12 // include files used for reconstructed tracks
13 // #include "DataFormats/TrackReco/interface/Track.h"
14 // #include "DataFormats/TrackReco/interface/TrackFwd.h"
15 // #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
16 // #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
17 // #include "DataFormats/TrackCandidate/interface/TrackCandidate.h"
18 // #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
19 // #include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
20 
37 // #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
39 
43 
49 
59 
60 #include <set>
61 #include <sstream>
62 
63 using namespace std;
64 using namespace edm;
65 
67 {
68 
69 
70  processParticles_ =
71  iConfig.getUntrackedParameter<bool>("process_Particles",true);
72 
73 
74  inputTagSim_
75  = iConfig.getParameter<InputTag>("sim");
76 
77  //retrieving collections for MC Truth Matching
78 
79  //modif-beg
80  inputTagFamosSimHits_
81  = iConfig.getUntrackedParameter<InputTag>("famosSimHits");
82  mctruthMatchingInfo_ =
83  iConfig.getUntrackedParameter<bool>("MCTruthMatchingInfo",false);
84  //modif-end
85 
86  inputTagRecTracks_
87  = iConfig.getParameter<InputTag>("RecTracks");
88  inputTagEcalRecHitsEB_
89  = iConfig.getParameter<InputTag>("ecalRecHitsEB");
90  inputTagEcalRecHitsEE_
91  = iConfig.getParameter<InputTag>("ecalRecHitsEE");
92 
93  verbose_ =
94  iConfig.getUntrackedParameter<bool>("verbose",false);
95 
96 
97  // register products
98  produces<reco::PFSimParticleCollection>();
99 
100  particleFilter_ = iConfig.getParameter<ParameterSet>
101  ( "ParticleFilter" );
102 
103  mySimEvent = new FSimEvent( particleFilter_ );
104 }
105 
106 
107 
109 {
110  delete mySimEvent;
111 }
112 
113 
114 void
116  const edm::EventSetup & es)
117 {
118  // init Particle data table (from Pythia)
120  // edm::ESHandle < DefaultConfig::ParticleDataTable > pdt;
121  es.getData(pdt);
122  mySimEvent->initializePdt(&(*pdt));
123 }
124 
125 
127  const EventSetup& iSetup)
128 {
129  ParticleTable::Sentry ptable(mySimEvent->theTable());
130  LogDebug("PFSimParticleProducer")<<"START event: "<<iEvent.id().event()
131  <<" in run "<<iEvent.id().run()<<endl;
132 
133  //MC Truth Matching only with Famos and UnFoldedMode option to true!!
134 
135  //vector to store the trackIDs of rectracks corresponding
136  //to the simulated particle.
137  std::vector<unsigned> recTrackSimID;
138 
139  //In order to know which simparticule contribute to
140  //a given Ecal RecHit energy, we need to access
141  //the PCAloHit from FastSim.
142 
143  typedef std::pair<double, unsigned> hitSimID;
144  typedef std::list< std::pair<double, unsigned> >::iterator ITM;
145  std::vector< std::list <hitSimID> > caloHitsEBID(62000);
146  std::vector<double> caloHitsEBTotE(62000,0.0);
147 
148  if(mctruthMatchingInfo_){
149 
150  //getting the PCAloHit
152  // bool found_phit
153  // = iEvent.getByLabel("famosSimHits","EcalHitsEB",
154  // pcalohits);
155  //modif-beg
156  bool found_phit
157  = iEvent.getByLabel(inputTagFamosSimHits_,pcalohits);
158  //modif-end
159 
160  if(!found_phit) {
161  ostringstream err;
162  err<<"could not find pcaloHit "<<"famosSimHits:EcalHitsEB";
163  LogError("PFSimParticleProducer")<<err.str()<<endl;
164 
165  throw cms::Exception( "MissingProduct", err.str());
166  }
167  else {
168  assert( pcalohits.isValid() );
169 
170  // cout << "PFSimParticleProducer: number of pcalohits="
171  // << pcalohits->size() << endl;
172 
173  edm::PCaloHitContainer::const_iterator it
174  = pcalohits.product()->begin();
175  edm::PCaloHitContainer::const_iterator itend
176  = pcalohits.product()->end();
177 
178  //loop on the PCaloHit from FastSim Calorimetry
179  for(;it!=itend;++it)
180  {
181  EBDetId detid(it->id());
182 
183  // cout << detid << " " << detid.rawId()
184  // << " " << detid.hashedIndex()
185  // << " " << it->energy()
186  // << " " << it->id()
187  // << " trackId="
188  // << it->geantTrackId()
189  // << endl;
190 
191  if(it->energy() > 0.0) {
192  std::pair<double, unsigned> phitsimid
193  = make_pair(it->energy(),it->geantTrackId());
194  caloHitsEBID[detid.hashedIndex()].push_back(phitsimid);
195  caloHitsEBTotE[detid.hashedIndex()]
196  += it->energy(); //summing pcalhit energy
197  }//energy > 0
198 
199  }//loop PcaloHits
200  }//pcalohit handle access
201 
202  //Retrieving the PFRecTrack collection for
203  //Monte Carlo Truth Matching tool
205  try{
206  LogDebug("PFSimParticleProducer")<<"getting PFRecTracks"<<endl;
207  iEvent.getByLabel(inputTagRecTracks_, recTracks);
208 
209  } catch (cms::Exception& err) {
210  LogError("PFSimParticleProducer")<<err
211  <<" cannot get collection "
212  <<"particleFlowBlock"<<":"
213  <<""
214  <<endl;
215  }//pfrectrack handle access
216 
217  //getting the simID corresponding to
218  //each of the PFRecTracks
219  getSimIDs( recTracks, recTrackSimID );
220 
221  }//mctruthMatchingInfo_ //modif
222 
223  // deal with true particles
224  if( processParticles_) {
225 
226  auto_ptr< reco::PFSimParticleCollection >
227  pOutputPFSimParticleCollection(new reco::PFSimParticleCollection );
228 
229  Handle<vector<SimTrack> > simTracks;
230  bool found = iEvent.getByLabel(inputTagSim_,simTracks);
231  if(!found) {
232 
233  ostringstream err;
234  err<<"cannot find sim tracks: "<<inputTagSim_;
235  LogError("PFSimParticleProducer")<<err.str()<<endl;
236 
237  throw cms::Exception( "MissingProduct", err.str());
238  }
239 
240 
241 
242  Handle<vector<SimVertex> > simVertices;
243  found = iEvent.getByLabel(inputTagSim_,simVertices);
244  if(!found) {
245  LogError("PFSimParticleProducer")
246  <<"cannot find sim vertices: "<<inputTagSim_<<endl;
247  return;
248  }
249 
250 
251 
252  // for(unsigned it = 0; it<simTracks->size(); it++ ) {
253  // cout<<"\t track "<< (*simTracks)[it]<<" "
254  // <<(*simTracks)[it].momentum().vect().perp()<<" "
255  // <<(*simTracks)[it].momentum().e()<<endl;
256  // }
257 
258  mySimEvent->fill( *simTracks, *simVertices );
259 
260  if(verbose_)
261  mySimEvent->print();
262 
263  // const std::vector<FSimTrack>& fsimTracks = *(mySimEvent->tracks() );
264  for(unsigned i=0; i<mySimEvent->nTracks(); i++) {
265 
266  const FSimTrack& fst = mySimEvent->track(i);
267 
268  int motherId = -1;
269  if( ! fst.noMother() ) // a mother exist
270  motherId = fst.mother().id();
271 
272  //This is finding out the simID corresponding
273  //to the recTrack
274 // cout << "Particle " << i
275 // << " " << fst.genpartIndex()
276 // << " -------------------------------------" << endl;
277 
278  //GETTING THE TRACK ID
279  unsigned recTrackID = 99999;
280  vector<unsigned> recHitContrib; //modif
281  vector<double> recHitContribFrac; //modif
282 
283  if(mctruthMatchingInfo_){ //modif
284 
285  for(unsigned lo=0; lo<recTrackSimID.size();
286  lo++) {
287  if( i == recTrackSimID[lo] ) {
288 // cout << "Corresponding Rec Track "
289 // << lo << endl;
290  recTrackID = lo;
291  }//match track
292  }//loop rectrack
293 // if( recTrackID == 99999 )
294 // cout << "Sim Track not reconstructed pT=" <<
295 // fst.momentum().pt() << endl;
296 
297  // get the ecalBarrel rechits for MC truth matching tool
299  bool found = iEvent.getByLabel(inputTagEcalRecHitsEB_,
300  rhcHandle);
301  if(!found) {
302  ostringstream err;
303  err<<"could not find rechits "<< inputTagEcalRecHitsEB_;
304  LogError("PFSimParticleProducer")<<err.str()<<endl;
305 
306  throw cms::Exception( "MissingProduct", err.str());
307  }
308  else {
309  assert( rhcHandle.isValid() );
310 // cout << "PFSimParticleProducer: number of rechits="
311 // << rhcHandle->size() << endl;
312 
314  = rhcHandle.product()->begin();
316  = rhcHandle.product()->end();
317 
318  for(;it_rh!=itend_rh;++it_rh)
319  {
320  unsigned rhit_hi
321  = EBDetId(it_rh->id()).hashedIndex();
322  EBDetId detid(it_rh->id());
323 // cout << detid << " " << detid.rawId()
324 // << " " << detid.hashedIndex()
325 // << " " << it_rh->energy() << endl;
326 
327  ITM it_phit = caloHitsEBID[rhit_hi].begin();
328  ITM itend_phit = caloHitsEBID[rhit_hi].end();
329  for(;it_phit!=itend_phit;++it_phit)
330  {
331  if(i == it_phit->second)
332  {
333  //Alex (08/10/08) TO BE REMOVED, eliminating
334  //duplicated rechits
335  bool alreadyin = false;
336  for( unsigned ihit = 0; ihit < recHitContrib.size();
337  ++ihit )
338  if(detid.rawId() == recHitContrib[ihit])
339  alreadyin = true;
340 
341  if(!alreadyin){
342  double pcalofraction = 0.0;
343  if(caloHitsEBTotE[rhit_hi] != 0.0)
344  pcalofraction
345  = (it_phit->first/caloHitsEBTotE[rhit_hi])*100.0;
346 
347  //store info
348  recHitContrib.push_back(it_rh->id());
349  recHitContribFrac.push_back(pcalofraction);
350  }//selected rechits
351  }//matching
352  }//loop pcalohit
353 
354  }//loop rechits
355 
356  }//getting the rechits
357 
358  }//mctruthMatchingInfo_ //modif
359 
360 // cout << "This particule has " << recHitContrib.size()
361 // << " rechit contribution" << endl;
362 // for( unsigned ih = 0; ih < recHitContrib.size(); ++ih )
363 // cout << recHitContrib[ih]
364 // << " f=" << recHitContribFrac[ih] << " ";
365 // cout << endl;
366 
367  reco::PFSimParticle particle( fst.charge(),
368  fst.type(),
369  fst.id(),
370  motherId,
371  fst.daughters(),
372  recTrackID,
373  recHitContrib,
374  recHitContribFrac);
375 
376 
377  const FSimVertex& originVtx = fst.vertex();
378 
379  math::XYZPoint posOrig( originVtx.position().x(),
380  originVtx.position().y(),
381  originVtx.position().z() );
382 
383  math::XYZTLorentzVector momOrig( fst.momentum().px(),
384  fst.momentum().py(),
385  fst.momentum().pz(),
386  fst.momentum().e() );
388  pointOrig(-1,
390  posOrig, momOrig);
391 
392  // point 0 is origin vertex
393  particle.addPoint(pointOrig);
394 
395 
396  if( ! fst.noEndVertex() ) {
397  const FSimVertex& endVtx = fst.endVertex();
398 
399  math::XYZPoint posEnd( endVtx.position().x(),
400  endVtx.position().y(),
401  endVtx.position().z() );
402  // cout<<"end vertex : "
403  // <<endVtx.position().x()<<" "
404  // <<endVtx.position().y()<<endl;
405 
407 
409  pointEnd( -1,
411  posEnd, momEnd);
412 
413  particle.addPoint(pointEnd);
414  }
415  else { // add a dummy point
417  particle.addPoint(dummy);
418  }
419 
420 
421  if( fst.onLayer1() ) { // PS layer1
422  const RawParticle& rp = fst.layer1Entrance();
423 
424  math::XYZPoint posLayer1( rp.x(), rp.y(), rp.z() );
425  math::XYZTLorentzVector momLayer1( rp.px(), rp.py(), rp.pz(), rp.e() );
427  posLayer1, momLayer1);
428 
429  particle.addPoint( layer1Pt );
430 
431  // extrapolate to cluster depth
432  }
433  else { // add a dummy point
435  particle.addPoint(dummy);
436  }
437 
438  if( fst.onLayer2() ) { // PS layer2
439  const RawParticle& rp = fst.layer2Entrance();
440 
441  math::XYZPoint posLayer2( rp.x(), rp.y(), rp.z() );
442  math::XYZTLorentzVector momLayer2( rp.px(), rp.py(), rp.pz(), rp.e() );
444  posLayer2, momLayer2);
445 
446  particle.addPoint( layer2Pt );
447 
448  // extrapolate to cluster depth
449  }
450  else { // add a dummy point
452  particle.addPoint(dummy);
453  }
454 
455  if( fst.onEcal() ) {
456  const RawParticle& rp = fst.ecalEntrance();
457 
458  math::XYZPoint posECAL( rp.x(), rp.y(), rp.z() );
459  math::XYZTLorentzVector momECAL( rp.px(), rp.py(), rp.pz(), rp.e() );
460  reco::PFTrajectoryPoint ecalPt(-1,
462  posECAL, momECAL);
463 
464  particle.addPoint( ecalPt );
465 
466  // extrapolate to cluster depth
467  }
468  else { // add a dummy point
470  particle.addPoint(dummy);
471  }
472 
473  // add a dummy point for ECAL Shower max
475  particle.addPoint(dummy);
476 
477  if( fst.onHcal() ) {
478 
479  const RawParticle& rpin = fst.hcalEntrance();
480 
481  math::XYZPoint posHCALin( rpin.x(), rpin.y(), rpin.z() );
482  math::XYZTLorentzVector momHCALin( rpin.px(), rpin.py(), rpin.pz(),
483  rpin.e() );
484  reco::PFTrajectoryPoint hcalPtin(-1,
486  posHCALin, momHCALin);
487 
488  particle.addPoint( hcalPtin );
489 
490  const RawParticle& rpout = fst.hcalExit();
491 
492  math::XYZPoint posHCALout( rpout.x(), rpout.y(), rpout.z() );
493  math::XYZTLorentzVector momHCALout( rpout.px(), rpout.py(), rpout.pz(),
494  rpout.e() );
496  hcalPtout(-1, reco::PFTrajectoryPoint::HCALExit,
497  posHCALout, momHCALout);
498 
499  particle.addPoint( hcalPtout );
500 
501  const RawParticle& rpho = fst.hoEntrance();
502 
503  math::XYZPoint posHOEntrance( rpho.x(), rpho.y(), rpho.z() );
504  math::XYZTLorentzVector momHOEntrance( rpho.px(), rpho.py(), rpho.pz(),
505  rpho.e() );
508  posHOEntrance, momHOEntrance);
509 
510  particle.addPoint( hoPtin );
511 
512 
513 
514 
515  }
516  else { // add a dummy point
518  particle.addPoint(dummy);
519  }
520 
521  pOutputPFSimParticleCollection->push_back( particle );
522  }
523 
524  iEvent.put(pOutputPFSimParticleCollection);
525  }
526 
527 
528  LogDebug("PFSimParticleProducer")<<"STOP event: "<<iEvent.id().event()
529  <<" in run "<<iEvent.id().run()<<endl;
530 }
531 
533  std::vector<unsigned>& recTrackSimID )
534 {
535 
536  if( trackh.isValid() ) {
537 // cout << "Size=" << trackh->size() << endl;
538 
539  for(unsigned i=0;i<trackh->size(); i++) {
540 
541  reco::PFRecTrackRef ref( trackh,i );
542  const reco::PFRecTrack& PFT = *ref;
543  const reco::TrackRef trackref = PFT.trackRef();
544 
545 // double Pt = trackref->pt();
546 // double DPt = trackref->ptError();
547 // cout << " PFBlockProducer: PFrecTrack->Track Pt= "
548 // << Pt << " DPt = " << DPt << endl;
549 
550  trackingRecHit_iterator rhitbeg
551  = trackref->recHitsBegin();
552  trackingRecHit_iterator rhitend
553  = trackref->recHitsEnd();
554  for (trackingRecHit_iterator it = rhitbeg;
555  it != rhitend; it++){
556 
557  if( (*it)->isValid() ){
558 
559  const FastTrackerRecHit * rechit
560  = (const FastTrackerRecHit*) (*it);
561 
562 // cout << "rechit"
563 // << " corresponding simId "
564 // << rechit->simtrackId()
565 // << endl;
566 
567  for(unsigned int st_index=0;st_index<rechit->nSimTrackIds();++st_index){
568  recTrackSimID.push_back(rechit->simTrackId(st_index));
569  }
570  break;
571  }
572  }//loop track rechit
573  }//loop recTracks
574  }//track handle valid
575 }
PFSimParticleProducer(const edm::ParameterSet &)
#define LogDebug(id)
RunNumber_t run() const
Definition: EventID.h:39
T getParameter(std::string const &) const
EventNumber_t event() const
Definition: EventID.h:41
T getUntrackedParameter(std::string const &, T const &) const
reconstructed track used as an input to particle flow
Definition: PFRecTrack.h:22
int i
Definition: DBlmapReader.cc:9
bool noEndVertex() const
no end vertex
float charge() const
charge
Definition: FSimTrack.h:47
const RawParticle & hoEntrance() const
The particle at HCAL exir.
Definition: FSimTrack.h:145
int onLayer2() const
Definition: FSimTrack.h:96
const FSimVertex & endVertex() const
end vertex
const RawParticle & hcalExit() const
The particle at HCAL exir.
Definition: FSimTrack.h:142
double y() const
y of vertex
Definition: RawParticle.h:271
double z() const
z of vertex
Definition: RawParticle.h:272
virtual void beginRun(const edm::Run &r, const edm::EventSetup &c) override
const XYZTLorentzVector & momentum() const
Temporary (until move of SimTrack to Mathcore) - No! Actually very useful.
Definition: FSimTrack.h:190
assert(m_qm.get())
const RawParticle & layer1Entrance() const
The particle at Preshower Layer 1.
Definition: FSimTrack.h:127
std::vector< EcalRecHit >::const_iterator const_iterator
std::vector< PFSimParticle > PFSimParticleCollection
collection of PFSimParticle objects
const std::vector< int > & daughters() const
Vector of daughter indices.
double x() const
x of vertex
Definition: RawParticle.h:270
const reco::TrackRef & trackRef() const
Definition: PFRecTrack.h:54
void getData(T &iHolder) const
Definition: EventSetup.h:79
int onEcal() const
Definition: FSimTrack.h:101
virtual void produce(edm::Event &, const edm::EventSetup &) override
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
int hashedIndex(int ieta, int iphi)
Definition: EcalPyUtils.cc:42
int iEvent
Definition: GenABIO.cc:230
Point of closest approach from beam axis (initial point in the case of PFSimParticle) ...
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:120
const math::XYZTLorentzVector & position() const
Temporary (until CMSSW moves to Mathcore) - No ! Actually very useful.
Definition: FSimVertex.h:49
virtual int32_t simTrackId(size_t i) const
const RawParticle & ecalEntrance() const
The particle at ECAL entrance.
Definition: FSimTrack.h:133
const FSimVertex & vertex() const
Origin vertex.
bool isValid() const
Definition: HandleBase.h:75
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:420
true particle for particle flow
Definition: PFSimParticle.h:19
T const * product() const
Definition: Handle.h:81
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
virtual size_t nSimTrackIds() const
edm::EventID id() const
Definition: EventBase.h:60
int type() const
particle type (HEP PDT convension)
Definition: CoreSimTrack.h:25
int id() const
the index in FBaseSimEvent and other vectors
Definition: FSimTrack.h:86
const RawParticle & layer2Entrance() const
The particle at Preshower Layer 2.
Definition: FSimTrack.h:130
int onLayer1() const
Definition: FSimTrack.h:91
int onHcal() const
Definition: FSimTrack.h:106
bool noMother() const
no mother particle
const RawParticle & hcalEntrance() const
The particle at HCAL entrance.
Definition: FSimTrack.h:136
A PFTrack holds several trajectory points, which basically contain the position and momentum of a tra...
void getSimIDs(const TrackHandle &trackh, std::vector< unsigned > &recTrackSimID)
const FSimTrack & mother() const
mother
TrackingRecHitCollection::base::const_iterator trackingRecHit_iterator
iterator over a vector of reference to TrackingRecHit in the same collection
Definition: Run.h:43