PFSimProducer/plugins/PFSimParticleProducer.cc

0001 /**\class PFSimParticleProducer
0002 \brief Producer for PFRecTracks and PFSimParticles
0003
0004 \todo Remove the PFRecTrack part, which is now handled by PFTracking
0005 \author Colin Bernet
0006 \date   April 2007
0007 */
0008
0009 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0010 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0011 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0012 #include "DataFormats/GeometrySurface/interface/BoundCylinder.h"
0013 #include "DataFormats/GeometrySurface/interface/SimpleCylinderBounds.h"
0014 #include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"
0015 #include "DataFormats/GeometrySurface/interface/TkRotation.h"
0016 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0017 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0018 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0019 #include "DataFormats/ParticleFlowReco/interface/PFLayer.h"
0020 #include "DataFormats/ParticleFlowReco/interface/PFRecTrack.h"
0021 #include "DataFormats/ParticleFlowReco/interface/PFRecTrackFwd.h"
0022 #include "DataFormats/ParticleFlowReco/interface/PFSimParticle.h"
0023 #include "DataFormats/ParticleFlowReco/interface/PFSimParticleFwd.h"
0024 #include "DataFormats/TrackReco/interface/Track.h"
0025 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0026 #include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHit.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/EventSetup.h"
0029 #include "FWCore/Framework/interface/stream/EDProducer.h"
0030 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0031 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0032 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0033 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0034 #include "FWCore/Utilities/interface/EDGetToken.h"
0035 #include "FWCore/Utilities/interface/Exception.h"
0036 #include "FastSimulation/Event/interface/FSimEvent.h"
0037 #include "FastSimulation/Event/interface/FSimTrack.h"
0038 #include "FastSimulation/Event/interface/FSimVertex.h"
0039 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0040 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0041 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0042 #include "MagneticField/Engine/interface/MagneticField.h"
0043 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0044 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h"
0045 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyResolution.h"
0046 #include "RecoParticleFlow/PFProducer/interface/PFBlockAlgo.h"
0047 #include "RecoParticleFlow/PFTracking/interface/PFGeometry.h"
0048 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0049 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0050 #include "SimDataFormats/Track/interface/SimTrack.h"
0051 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0052 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0053 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0054 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0055 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0056 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0057 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
0058 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0059
0060 #include <memory>
0061 #include <set>
0062 #include <sstream>
0063 #include <string>
0064
0065 class PFSimParticleProducer : public edm::stream::EDProducer<> {
0066 public:
0067   explicit PFSimParticleProducer(const edm::ParameterSet&);
0068
0069   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0070
0071   void produce(edm::Event&, const edm::EventSetup&) override;
0072
0073   typedef edm::Handle<reco::PFRecTrackCollection> TrackHandle;
0074   void getSimIDs(const TrackHandle& trackh, std::vector<unsigned>& recTrackSimID);
0075
0076 private:
0077   /// module label for retrieving input simtrack and simvertex
0078   edm::InputTag inputTagSim_;
0079   edm::EDGetTokenT<std::vector<SimTrack> > tokenSim_;
0080   edm::EDGetTokenT<std::vector<SimVertex> > tokenSimVertices_;
0081
0082   //MC Truth Matching
0083   //modif-beg
0084   bool mctruthMatchingInfo_;
0085   edm::InputTag inputTagFastSimProducer_;
0086   edm::EDGetTokenT<edm::PCaloHitContainer> tokenFastSimProducer_;
0087   //modif-end
0088
0089   edm::InputTag inputTagRecTracks_;
0090   edm::EDGetTokenT<reco::PFRecTrackCollection> tokenRecTracks_;
0091   edm::InputTag inputTagEcalRecHitsEB_;
0092   edm::EDGetTokenT<EcalRecHitCollection> tokenEcalRecHitsEB_;
0093   edm::InputTag inputTagEcalRecHitsEE_;
0094   edm::EDGetTokenT<EcalRecHitCollection> tokenEcalRecHitsEE_;
0095
0096   edm::ESGetToken<HepPDT::ParticleDataTable, edm::DefaultRecord> pdtToken_;
0097
0098   // parameters for retrieving true particles information --
0099
0100   edm::ParameterSet particleFilter_;
0101   std::unique_ptr<FSimEvent> mySimEvent;
0102
0103   // flags for the various tasks ---------------------------
0104
0105   /// process particles on/off
0106   bool processParticles_;
0107
0108   /// verbose ?
0109   bool verbose_;
0110 };
0111
0112 #include "FWCore/Framework/interface/MakerMacros.h"
0113 DEFINE_FWK_MODULE(PFSimParticleProducer);
0114
0115 void PFSimParticleProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0116   // particleFlowSimParticle
0117   edm::ParameterSetDescription desc;
0118   desc.addUntracked<edm::InputTag>("fastSimProducer", edm::InputTag("fastSimProducer", "EcalHitsEB"));
0119   desc.addUntracked<bool>("MCTruthMatchingInfo", false);
0120   desc.add<edm::InputTag>("RecTracks", edm::InputTag("trackerDrivenElectronSeeds"));
0121   desc.add<std::string>("Fitter", "KFFittingSmoother");
0122   desc.add<edm::InputTag>("ecalRecHitsEE", edm::InputTag("caloRecHits", "EcalRecHitsEE"));
0123   desc.add<edm::InputTag>("ecalRecHitsEB", edm::InputTag("caloRecHits", "EcalRecHitsEB"));
0124   desc.addUntracked<bool>("process_RecTracks", false);
0125   {
0126     edm::ParameterSetDescription psd0;
0127     psd0.setUnknown();
0128     desc.add<edm::ParameterSetDescription>("ParticleFilter", psd0);
0129   }
0130   desc.add<std::string>("TTRHBuilder", "WithTrackAngle");
0131   desc.addUntracked<bool>("process_Particles", true);
0132   desc.add<std::string>("Propagator", "PropagatorWithMaterial");
0133   desc.add<edm::InputTag>("sim", edm::InputTag("g4SimHits"));
0134   desc.addUntracked<bool>("verbose", false);
0135   descriptions.add("particleFlowSimParticle", desc);
0136 }
0137
0138 using namespace std;
0139 using namespace edm;
0140
0141 PFSimParticleProducer::PFSimParticleProducer(const edm::ParameterSet& iConfig) {
0142   processParticles_ = iConfig.getUntrackedParameter<bool>("process_Particles", true);
0143
0144   inputTagSim_ = iConfig.getParameter<InputTag>("sim");
0145   tokenSim_ = consumes<std::vector<SimTrack> >(inputTagSim_);
0146   tokenSimVertices_ = consumes<std::vector<SimVertex> >(inputTagSim_);
0147
0148   //retrieving collections for MC Truth Matching
0149
0150   //modif-beg
0151   inputTagFastSimProducer_ = iConfig.getUntrackedParameter<InputTag>("fastSimProducer");
0152   tokenFastSimProducer_ = consumes<edm::PCaloHitContainer>(inputTagFastSimProducer_);
0153   mctruthMatchingInfo_ = iConfig.getUntrackedParameter<bool>("MCTruthMatchingInfo", false);
0154   //modif-end
0155
0156   inputTagRecTracks_ = iConfig.getParameter<InputTag>("RecTracks");
0157   tokenRecTracks_ = consumes<reco::PFRecTrackCollection>(inputTagRecTracks_);
0158
0159   inputTagEcalRecHitsEB_ = iConfig.getParameter<InputTag>("ecalRecHitsEB");
0160   tokenEcalRecHitsEB_ = consumes<EcalRecHitCollection>(inputTagEcalRecHitsEB_);
0161   inputTagEcalRecHitsEE_ = iConfig.getParameter<InputTag>("ecalRecHitsEE");
0162   tokenEcalRecHitsEE_ = consumes<EcalRecHitCollection>(inputTagEcalRecHitsEE_);
0163
0164   pdtToken_ = esConsumes();
0165
0166   verbose_ = iConfig.getUntrackedParameter<bool>("verbose", false);
0167
0168   // register products
0169   produces<reco::PFSimParticleCollection>();
0170
0171   particleFilter_ = iConfig.getParameter<ParameterSet>("ParticleFilter");
0172
0173   mySimEvent = std::make_unique<FSimEvent>(particleFilter_);
0174 }
0175
0176 void PFSimParticleProducer::produce(Event& iEvent, const EventSetup& iSetup) {
0177   // init Particle data table (from Pythia)
0178   mySimEvent->initializePdt(&iSetup.getData(pdtToken_));
0179
0180   LogDebug("PFSimParticleProducer") << "START event: " << iEvent.id().event() << " in run " << iEvent.id().run()
0181                                     << endl;
0182
0183   //MC Truth Matching only with Famos and UnFoldedMode option to true!!
0184
0185   //vector to store the trackIDs of rectracks corresponding
0186   //to the simulated particle.
0187   std::vector<unsigned> recTrackSimID;
0188
0189   //In order to know which simparticule contribute to
0190   //a given Ecal RecHit energy, we need to access
0191   //the PCAloHit from FastSim.
0192
0193   typedef std::pair<double, unsigned> hitSimID;
0194   std::vector<std::list<hitSimID> > caloHitsEBID(62000);
0195   std::vector<double> caloHitsEBTotE(62000, 0.0);
0196
0197   if (mctruthMatchingInfo_) {
0198     //getting the PCAloHit
0199     auto pcalohits = iEvent.getHandle(tokenFastSimProducer_);
0200
0201     if (!pcalohits) {
0202       ostringstream err;
0203       err << "could not find pcaloHit "
0204           << "fastSimProducer:EcalHitsEB";
0205       LogError("PFSimParticleProducer") << err.str() << endl;
0206
0207       throw cms::Exception("MissingProduct", err.str());
0208     } else {
0209       assert(pcalohits.isValid());
0210
0211       edm::PCaloHitContainer::const_iterator it = pcalohits.product()->begin();
0212       edm::PCaloHitContainer::const_iterator itend = pcalohits.product()->end();
0213
0214       //loop on the PCaloHit from FastSim Calorimetry
0215       for (; it != itend; ++it) {
0216         EBDetId detid(it->id());
0217
0218         if (it->energy() > 0.0) {
0219           std::pair<double, unsigned> phitsimid = make_pair(it->energy(), it->geantTrackId());
0220           caloHitsEBID[detid.hashedIndex()].push_back(phitsimid);
0221           caloHitsEBTotE[detid.hashedIndex()] += it->energy();  //summing pcalhit energy
0222         }                                                       //energy > 0
0223
0224       }  //loop PcaloHits
0225     }    //pcalohit handle access
0226
0227     //Retrieving the PFRecTrack collection for
0228     //Monte Carlo Truth Matching tool
0229     Handle<reco::PFRecTrackCollection> recTracks;
0230     try {
0231       LogDebug("PFSimParticleProducer") << "getting PFRecTracks" << endl;
0232       iEvent.getByToken(tokenRecTracks_, recTracks);
0233
0234     } catch (cms::Exception& err) {
0235       LogError("PFSimParticleProducer") << err << " cannot get collection "
0236                                         << "particleFlowBlock"
0237                                         << ":"
0238                                         << "" << endl;
0239     }  //pfrectrack handle access
0240
0241     //getting the simID corresponding to
0242     //each of the PFRecTracks
0243     getSimIDs(recTracks, recTrackSimID);
0244
0245   }  //mctruthMatchingInfo_ //modif
0246
0247   // deal with true particles
0248   if (processParticles_) {
0249     auto pOutputPFSimParticleCollection = std::make_unique<reco::PFSimParticleCollection>();
0250
0251     Handle<vector<SimTrack> > simTracks;
0252     bool found = iEvent.getByToken(tokenSim_, simTracks);
0253     if (!found) {
0254       ostringstream err;
0255       err << "cannot find sim tracks: " << inputTagSim_;
0256       LogError("PFSimParticleProducer") << err.str() << endl;
0257
0258       throw cms::Exception("MissingProduct", err.str());
0259     }
0260
0261     Handle<vector<SimVertex> > simVertices;
0262     found = iEvent.getByToken(tokenSimVertices_, simVertices);
0263     if (!found) {
0264       LogError("PFSimParticleProducer") << "cannot find sim vertices: " << inputTagSim_ << endl;
0265       return;
0266     }
0267
0268     mySimEvent->fill(*simTracks, *simVertices);
0269
0270     if (verbose_)
0271       mySimEvent->print();
0272
0273     for (unsigned i = 0; i < mySimEvent->nTracks(); i++) {
0274       const FSimTrack& fst = mySimEvent->track(i);
0275
0276       int motherId = -1;
0277       if (!fst.noMother())  // a mother exist
0278         motherId = fst.mother().id();
0279
0280       //This is finding out the simID corresponding
0281       //to the recTrack
0282
0283       //GETTING THE TRACK ID
0284       unsigned recTrackID = 99999;
0285       vector<unsigned> recHitContrib;    //modif
0286       vector<double> recHitContribFrac;  //modif
0287
0288       if (mctruthMatchingInfo_) {  //modif
0289
0290         for (unsigned lo = 0; lo < recTrackSimID.size(); lo++) {
0291           if (i == recTrackSimID[lo]) {
0292             recTrackID = lo;
0293           }  //match track
0294         }    //loop rectrack
0295
0296         // get the ecalBarrel rechits for MC truth matching tool
0297         edm::Handle<EcalRecHitCollection> rhcHandle;
0298         bool found = iEvent.getByToken(tokenEcalRecHitsEB_, rhcHandle);
0299         if (!found) {
0300           ostringstream err;
0301           err << "could not find rechits " << inputTagEcalRecHitsEB_;
0302           LogError("PFSimParticleProducer") << err.str() << endl;
0303
0304           throw cms::Exception("MissingProduct", err.str());
0305         } else {
0306           assert(rhcHandle.isValid());
0307
0308           EBRecHitCollection::const_iterator it_rh = rhcHandle.product()->begin();
0309           EBRecHitCollection::const_iterator itend_rh = rhcHandle.product()->end();
0310
0311           for (; it_rh != itend_rh; ++it_rh) {
0312             unsigned rhit_hi = EBDetId(it_rh->id()).hashedIndex();
0313             EBDetId detid(it_rh->id());
0314
0315             auto it_phit = caloHitsEBID[rhit_hi].begin();
0316             auto itend_phit = caloHitsEBID[rhit_hi].end();
0317             for (; it_phit != itend_phit; ++it_phit) {
0318               if (i == it_phit->second) {
0319                 //Alex (08/10/08) TO BE REMOVED, eliminating
0320                 //duplicated rechits
0321                 bool alreadyin = false;
0322                 for (unsigned ihit = 0; ihit < recHitContrib.size(); ++ihit)
0323                   if (detid.rawId() == recHitContrib[ihit])
0324                     alreadyin = true;
0325
0326                 if (!alreadyin) {
0327                   double pcalofraction = 0.0;
0328                   if (caloHitsEBTotE[rhit_hi] != 0.0)
0329                     pcalofraction = (it_phit->first / caloHitsEBTotE[rhit_hi]) * 100.0;
0330
0331                   //store info
0332                   recHitContrib.push_back(it_rh->id());
0333                   recHitContribFrac.push_back(pcalofraction);
0334                 }  //selected rechits
0335               }    //matching
0336             }      //loop pcalohit
0337
0338           }  //loop rechits
0339
0340         }  //getting the rechits
0341
0342       }  //mctruthMatchingInfo_ //modif
0343
0344       reco::PFSimParticle particle(
0345           fst.charge(), fst.type(), fst.id(), motherId, fst.daughters(), recTrackID, recHitContrib, recHitContribFrac);
0346
0347       const FSimVertex& originVtx = fst.vertex();
0348
0349       math::XYZPoint posOrig(originVtx.position().x(), originVtx.position().y(), originVtx.position().z());
0350
0351       math::XYZTLorentzVector momOrig(
0352           fst.momentum().px(), fst.momentum().py(), fst.momentum().pz(), fst.momentum().e());
0353       reco::PFTrajectoryPoint pointOrig(-1, reco::PFTrajectoryPoint::ClosestApproach, posOrig, momOrig);
0354
0355       // point 0 is origin vertex
0356       particle.addPoint(pointOrig);
0357
0358       if (!fst.noEndVertex()) {
0359         const FSimVertex& endVtx = fst.endVertex();
0360
0361         math::XYZPoint posEnd(endVtx.position().x(), endVtx.position().y(), endVtx.position().z());
0362
0363         math::XYZTLorentzVector momEnd;
0364
0365         reco::PFTrajectoryPoint pointEnd(-1, reco::PFTrajectoryPoint::BeamPipeOrEndVertex, posEnd, momEnd);
0366
0367         particle.addPoint(pointEnd);
0368       } else {  // add a dummy point
0369         reco::PFTrajectoryPoint dummy;
0370         particle.addPoint(dummy);
0371       }
0372
0373       if (fst.onLayer1()) {  // PS layer1
0374         const RawParticle& rp = fst.layer1Entrance();
0375
0376         math::XYZPoint posLayer1(rp.x(), rp.y(), rp.z());
0377         math::XYZTLorentzVector momLayer1(rp.px(), rp.py(), rp.pz(), rp.e());
0378         reco::PFTrajectoryPoint layer1Pt(-1, reco::PFTrajectoryPoint::PS1, posLayer1, momLayer1);
0379
0380         particle.addPoint(layer1Pt);
0381
0382         // extrapolate to cluster depth
0383       } else {  // add a dummy point
0384         reco::PFTrajectoryPoint dummy;
0385         particle.addPoint(dummy);
0386       }
0387
0388       if (fst.onLayer2()) {  // PS layer2
0389         const RawParticle& rp = fst.layer2Entrance();
0390
0391         math::XYZPoint posLayer2(rp.x(), rp.y(), rp.z());
0392         math::XYZTLorentzVector momLayer2(rp.px(), rp.py(), rp.pz(), rp.e());
0393         reco::PFTrajectoryPoint layer2Pt(-1, reco::PFTrajectoryPoint::PS2, posLayer2, momLayer2);
0394
0395         particle.addPoint(layer2Pt);
0396
0397         // extrapolate to cluster depth
0398       } else {  // add a dummy point
0399         reco::PFTrajectoryPoint dummy;
0400         particle.addPoint(dummy);
0401       }
0402
0403       if (fst.onEcal()) {
0404         const RawParticle& rp = fst.ecalEntrance();
0405
0406         math::XYZPoint posECAL(rp.x(), rp.y(), rp.z());
0407         math::XYZTLorentzVector momECAL(rp.px(), rp.py(), rp.pz(), rp.e());
0408         reco::PFTrajectoryPoint ecalPt(-1, reco::PFTrajectoryPoint::ECALEntrance, posECAL, momECAL);
0409
0410         particle.addPoint(ecalPt);
0411
0412         // extrapolate to cluster depth
0413       } else {  // add a dummy point
0414         reco::PFTrajectoryPoint dummy;
0415         particle.addPoint(dummy);
0416       }
0417
0418       // add a dummy point for ECAL Shower max
0419       reco::PFTrajectoryPoint dummy;
0420       particle.addPoint(dummy);
0421
0422       if (fst.onHcal()) {
0423         const RawParticle& rpin = fst.hcalEntrance();
0424
0425         math::XYZPoint posHCALin(rpin.x(), rpin.y(), rpin.z());
0426         math::XYZTLorentzVector momHCALin(rpin.px(), rpin.py(), rpin.pz(), rpin.e());
0427         reco::PFTrajectoryPoint hcalPtin(-1, reco::PFTrajectoryPoint::HCALEntrance, posHCALin, momHCALin);
0428
0429         particle.addPoint(hcalPtin);
0430
0431         const RawParticle& rpout = fst.hcalExit();
0432
0433         math::XYZPoint posHCALout(rpout.x(), rpout.y(), rpout.z());
0434         math::XYZTLorentzVector momHCALout(rpout.px(), rpout.py(), rpout.pz(), rpout.e());
0435         reco::PFTrajectoryPoint hcalPtout(-1, reco::PFTrajectoryPoint::HCALExit, posHCALout, momHCALout);
0436
0437         particle.addPoint(hcalPtout);
0438
0439         const RawParticle& rpho = fst.hoEntrance();
0440
0441         math::XYZPoint posHOEntrance(rpho.x(), rpho.y(), rpho.z());
0442         math::XYZTLorentzVector momHOEntrance(rpho.px(), rpho.py(), rpho.pz(), rpho.e());
0443         reco::PFTrajectoryPoint hoPtin(-1, reco::PFTrajectoryPoint::HOLayer, posHOEntrance, momHOEntrance);
0444
0445         particle.addPoint(hoPtin);
0446
0447       } else {  // add a dummy point
0448         reco::PFTrajectoryPoint dummy;
0449         particle.addPoint(dummy);
0450       }
0451
0452       pOutputPFSimParticleCollection->push_back(particle);
0453     }
0454
0455     iEvent.put(std::move(pOutputPFSimParticleCollection));
0456   }
0457
0458   LogDebug("PFSimParticleProducer") << "STOP event: " << iEvent.id().event() << " in run " << iEvent.id().run() << endl;
0459 }
0460
0461 void PFSimParticleProducer::getSimIDs(const TrackHandle& trackh, std::vector<unsigned>& recTrackSimID) {
0462   if (trackh.isValid()) {
0463     for (unsigned i = 0; i < trackh->size(); i++) {
0464       const reco::PFRecTrackRef ref(trackh, i);
0465
0466       for (auto const& hit : ref->trackRef()->recHits()) {
0467         if (hit->isValid()) {
0468           auto rechit = dynamic_cast<const FastTrackerRecHit*>(hit);
0469
0470           for (unsigned int st_index = 0; st_index < rechit->nSimTrackIds(); ++st_index) {
0471             recTrackSimID.push_back(rechit->simTrackId(st_index));
0472           }
0473           break;
0474         }
0475       }  //loop track rechit
0476     }    //loop recTracks
0477   }      //track handle valid
0478 }