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   desc.add<edm::ParameterSetDescription>("ParticleFilter", {});
0126   desc.add<std::string>("TTRHBuilder", "WithTrackAngle");
0127   desc.addUntracked<bool>("process_Particles", true);
0128   desc.add<std::string>("Propagator", "PropagatorWithMaterial");
0129   desc.add<edm::InputTag>("sim", edm::InputTag("g4SimHits"));
0130   desc.addUntracked<bool>("verbose", false);
0131   descriptions.add("particleFlowSimParticle", desc);
0132 }
0133
0134 using namespace std;
0135 using namespace edm;
0136
0137 PFSimParticleProducer::PFSimParticleProducer(const edm::ParameterSet& iConfig) {
0138   processParticles_ = iConfig.getUntrackedParameter<bool>("process_Particles", true);
0139
0140   inputTagSim_ = iConfig.getParameter<InputTag>("sim");
0141   tokenSim_ = consumes<std::vector<SimTrack> >(inputTagSim_);
0142   tokenSimVertices_ = consumes<std::vector<SimVertex> >(inputTagSim_);
0143
0144   //retrieving collections for MC Truth Matching
0145
0146   //modif-beg
0147   inputTagFastSimProducer_ = iConfig.getUntrackedParameter<InputTag>("fastSimProducer");
0148   tokenFastSimProducer_ = consumes<edm::PCaloHitContainer>(inputTagFastSimProducer_);
0149   mctruthMatchingInfo_ = iConfig.getUntrackedParameter<bool>("MCTruthMatchingInfo", false);
0150   //modif-end
0151
0152   inputTagRecTracks_ = iConfig.getParameter<InputTag>("RecTracks");
0153   tokenRecTracks_ = consumes<reco::PFRecTrackCollection>(inputTagRecTracks_);
0154
0155   inputTagEcalRecHitsEB_ = iConfig.getParameter<InputTag>("ecalRecHitsEB");
0156   tokenEcalRecHitsEB_ = consumes<EcalRecHitCollection>(inputTagEcalRecHitsEB_);
0157   inputTagEcalRecHitsEE_ = iConfig.getParameter<InputTag>("ecalRecHitsEE");
0158   tokenEcalRecHitsEE_ = consumes<EcalRecHitCollection>(inputTagEcalRecHitsEE_);
0159
0160   pdtToken_ = esConsumes();
0161
0162   verbose_ = iConfig.getUntrackedParameter<bool>("verbose", false);
0163
0164   // register products
0165   produces<reco::PFSimParticleCollection>();
0166
0167   particleFilter_ = iConfig.getParameter<ParameterSet>("ParticleFilter");
0168
0169   mySimEvent = std::make_unique<FSimEvent>(particleFilter_);
0170 }
0171
0172 void PFSimParticleProducer::produce(Event& iEvent, const EventSetup& iSetup) {
0173   // init Particle data table (from Pythia)
0174   mySimEvent->initializePdt(&iSetup.getData(pdtToken_));
0175
0176   LogDebug("PFSimParticleProducer") << "START event: " << iEvent.id().event() << " in run " << iEvent.id().run()
0177                                     << endl;
0178
0179   //MC Truth Matching only with Famos and UnFoldedMode option to true!!
0180
0181   //vector to store the trackIDs of rectracks corresponding
0182   //to the simulated particle.
0183   std::vector<unsigned> recTrackSimID;
0184
0185   //In order to know which simparticule contribute to
0186   //a given Ecal RecHit energy, we need to access
0187   //the PCAloHit from FastSim.
0188
0189   typedef std::pair<double, unsigned> hitSimID;
0190   std::vector<std::list<hitSimID> > caloHitsEBID(62000);
0191   std::vector<double> caloHitsEBTotE(62000, 0.0);
0192
0193   if (mctruthMatchingInfo_) {
0194     //getting the PCAloHit
0195     auto pcalohits = iEvent.getHandle(tokenFastSimProducer_);
0196
0197     if (!pcalohits) {
0198       ostringstream err;
0199       err << "could not find pcaloHit "
0200           << "fastSimProducer:EcalHitsEB";
0201       LogError("PFSimParticleProducer") << err.str() << endl;
0202
0203       throw cms::Exception("MissingProduct", err.str());
0204     } else {
0205       assert(pcalohits.isValid());
0206
0207       edm::PCaloHitContainer::const_iterator it = pcalohits.product()->begin();
0208       edm::PCaloHitContainer::const_iterator itend = pcalohits.product()->end();
0209
0210       //loop on the PCaloHit from FastSim Calorimetry
0211       for (; it != itend; ++it) {
0212         EBDetId detid(it->id());
0213
0214         if (it->energy() > 0.0) {
0215           std::pair<double, unsigned> phitsimid = make_pair(it->energy(), it->geantTrackId());
0216           caloHitsEBID[detid.hashedIndex()].push_back(phitsimid);
0217           caloHitsEBTotE[detid.hashedIndex()] += it->energy();  //summing pcalhit energy
0218         }                                                       //energy > 0
0219
0220       }  //loop PcaloHits
0221     }    //pcalohit handle access
0222
0223     //Retrieving the PFRecTrack collection for
0224     //Monte Carlo Truth Matching tool
0225     Handle<reco::PFRecTrackCollection> recTracks;
0226     try {
0227       LogDebug("PFSimParticleProducer") << "getting PFRecTracks" << endl;
0228       iEvent.getByToken(tokenRecTracks_, recTracks);
0229
0230     } catch (cms::Exception& err) {
0231       LogError("PFSimParticleProducer") << err << " cannot get collection "
0232                                         << "particleFlowBlock"
0233                                         << ":"
0234                                         << "" << endl;
0235     }  //pfrectrack handle access
0236
0237     //getting the simID corresponding to
0238     //each of the PFRecTracks
0239     getSimIDs(recTracks, recTrackSimID);
0240
0241   }  //mctruthMatchingInfo_ //modif
0242
0243   // deal with true particles
0244   if (processParticles_) {
0245     auto pOutputPFSimParticleCollection = std::make_unique<reco::PFSimParticleCollection>();
0246
0247     Handle<vector<SimTrack> > simTracks;
0248     bool found = iEvent.getByToken(tokenSim_, simTracks);
0249     if (!found) {
0250       ostringstream err;
0251       err << "cannot find sim tracks: " << inputTagSim_;
0252       LogError("PFSimParticleProducer") << err.str() << endl;
0253
0254       throw cms::Exception("MissingProduct", err.str());
0255     }
0256
0257     Handle<vector<SimVertex> > simVertices;
0258     found = iEvent.getByToken(tokenSimVertices_, simVertices);
0259     if (!found) {
0260       LogError("PFSimParticleProducer") << "cannot find sim vertices: " << inputTagSim_ << endl;
0261       return;
0262     }
0263
0264     mySimEvent->fill(*simTracks, *simVertices);
0265
0266     if (verbose_)
0267       mySimEvent->print();
0268
0269     for (unsigned i = 0; i < mySimEvent->nTracks(); i++) {
0270       const FSimTrack& fst = mySimEvent->track(i);
0271
0272       int motherId = -1;
0273       if (!fst.noMother())  // a mother exist
0274         motherId = fst.mother().id();
0275
0276       //This is finding out the simID corresponding
0277       //to the recTrack
0278
0279       //GETTING THE TRACK ID
0280       unsigned recTrackID = 99999;
0281       vector<unsigned> recHitContrib;    //modif
0282       vector<double> recHitContribFrac;  //modif
0283
0284       if (mctruthMatchingInfo_) {  //modif
0285
0286         for (unsigned lo = 0; lo < recTrackSimID.size(); lo++) {
0287           if (i == recTrackSimID[lo]) {
0288             recTrackID = lo;
0289           }  //match track
0290         }    //loop rectrack
0291
0292         // get the ecalBarrel rechits for MC truth matching tool
0293         edm::Handle<EcalRecHitCollection> rhcHandle;
0294         bool found = iEvent.getByToken(tokenEcalRecHitsEB_, rhcHandle);
0295         if (!found) {
0296           ostringstream err;
0297           err << "could not find rechits " << inputTagEcalRecHitsEB_;
0298           LogError("PFSimParticleProducer") << err.str() << endl;
0299
0300           throw cms::Exception("MissingProduct", err.str());
0301         } else {
0302           assert(rhcHandle.isValid());
0303
0304           EBRecHitCollection::const_iterator it_rh = rhcHandle.product()->begin();
0305           EBRecHitCollection::const_iterator itend_rh = rhcHandle.product()->end();
0306
0307           for (; it_rh != itend_rh; ++it_rh) {
0308             unsigned rhit_hi = EBDetId(it_rh->id()).hashedIndex();
0309             EBDetId detid(it_rh->id());
0310
0311             auto it_phit = caloHitsEBID[rhit_hi].begin();
0312             auto itend_phit = caloHitsEBID[rhit_hi].end();
0313             for (; it_phit != itend_phit; ++it_phit) {
0314               if (i == it_phit->second) {
0315                 //Alex (08/10/08) TO BE REMOVED, eliminating
0316                 //duplicated rechits
0317                 bool alreadyin = false;
0318                 for (unsigned ihit = 0; ihit < recHitContrib.size(); ++ihit)
0319                   if (detid.rawId() == recHitContrib[ihit])
0320                     alreadyin = true;
0321
0322                 if (!alreadyin) {
0323                   double pcalofraction = 0.0;
0324                   if (caloHitsEBTotE[rhit_hi] != 0.0)
0325                     pcalofraction = (it_phit->first / caloHitsEBTotE[rhit_hi]) * 100.0;
0326
0327                   //store info
0328                   recHitContrib.push_back(it_rh->id());
0329                   recHitContribFrac.push_back(pcalofraction);
0330                 }  //selected rechits
0331               }    //matching
0332             }      //loop pcalohit
0333
0334           }  //loop rechits
0335
0336         }  //getting the rechits
0337
0338       }  //mctruthMatchingInfo_ //modif
0339
0340       reco::PFSimParticle particle(
0341           fst.charge(), fst.type(), fst.id(), motherId, fst.daughters(), recTrackID, recHitContrib, recHitContribFrac);
0342
0343       const FSimVertex& originVtx = fst.vertex();
0344
0345       math::XYZPoint posOrig(originVtx.position().x(), originVtx.position().y(), originVtx.position().z());
0346
0347       math::XYZTLorentzVector momOrig(
0348           fst.momentum().px(), fst.momentum().py(), fst.momentum().pz(), fst.momentum().e());
0349       reco::PFTrajectoryPoint pointOrig(-1, reco::PFTrajectoryPoint::ClosestApproach, posOrig, momOrig);
0350
0351       // point 0 is origin vertex
0352       particle.addPoint(pointOrig);
0353
0354       if (!fst.noEndVertex()) {
0355         const FSimVertex& endVtx = fst.endVertex();
0356
0357         math::XYZPoint posEnd(endVtx.position().x(), endVtx.position().y(), endVtx.position().z());
0358
0359         math::XYZTLorentzVector momEnd;
0360
0361         reco::PFTrajectoryPoint pointEnd(-1, reco::PFTrajectoryPoint::BeamPipeOrEndVertex, posEnd, momEnd);
0362
0363         particle.addPoint(pointEnd);
0364       } else {  // add a dummy point
0365         reco::PFTrajectoryPoint dummy;
0366         particle.addPoint(dummy);
0367       }
0368
0369       if (fst.onLayer1()) {  // PS layer1
0370         const RawParticle& rp = fst.layer1Entrance();
0371
0372         math::XYZPoint posLayer1(rp.x(), rp.y(), rp.z());
0373         math::XYZTLorentzVector momLayer1(rp.px(), rp.py(), rp.pz(), rp.e());
0374         reco::PFTrajectoryPoint layer1Pt(-1, reco::PFTrajectoryPoint::PS1, posLayer1, momLayer1);
0375
0376         particle.addPoint(layer1Pt);
0377
0378         // extrapolate to cluster depth
0379       } else {  // add a dummy point
0380         reco::PFTrajectoryPoint dummy;
0381         particle.addPoint(dummy);
0382       }
0383
0384       if (fst.onLayer2()) {  // PS layer2
0385         const RawParticle& rp = fst.layer2Entrance();
0386
0387         math::XYZPoint posLayer2(rp.x(), rp.y(), rp.z());
0388         math::XYZTLorentzVector momLayer2(rp.px(), rp.py(), rp.pz(), rp.e());
0389         reco::PFTrajectoryPoint layer2Pt(-1, reco::PFTrajectoryPoint::PS2, posLayer2, momLayer2);
0390
0391         particle.addPoint(layer2Pt);
0392
0393         // extrapolate to cluster depth
0394       } else {  // add a dummy point
0395         reco::PFTrajectoryPoint dummy;
0396         particle.addPoint(dummy);
0397       }
0398
0399       if (fst.onEcal()) {
0400         const RawParticle& rp = fst.ecalEntrance();
0401
0402         math::XYZPoint posECAL(rp.x(), rp.y(), rp.z());
0403         math::XYZTLorentzVector momECAL(rp.px(), rp.py(), rp.pz(), rp.e());
0404         reco::PFTrajectoryPoint ecalPt(-1, reco::PFTrajectoryPoint::ECALEntrance, posECAL, momECAL);
0405
0406         particle.addPoint(ecalPt);
0407
0408         // extrapolate to cluster depth
0409       } else {  // add a dummy point
0410         reco::PFTrajectoryPoint dummy;
0411         particle.addPoint(dummy);
0412       }
0413
0414       // add a dummy point for ECAL Shower max
0415       reco::PFTrajectoryPoint dummy;
0416       particle.addPoint(dummy);
0417
0418       if (fst.onHcal()) {
0419         const RawParticle& rpin = fst.hcalEntrance();
0420
0421         math::XYZPoint posHCALin(rpin.x(), rpin.y(), rpin.z());
0422         math::XYZTLorentzVector momHCALin(rpin.px(), rpin.py(), rpin.pz(), rpin.e());
0423         reco::PFTrajectoryPoint hcalPtin(-1, reco::PFTrajectoryPoint::HCALEntrance, posHCALin, momHCALin);
0424
0425         particle.addPoint(hcalPtin);
0426
0427         const RawParticle& rpout = fst.hcalExit();
0428
0429         math::XYZPoint posHCALout(rpout.x(), rpout.y(), rpout.z());
0430         math::XYZTLorentzVector momHCALout(rpout.px(), rpout.py(), rpout.pz(), rpout.e());
0431         reco::PFTrajectoryPoint hcalPtout(-1, reco::PFTrajectoryPoint::HCALExit, posHCALout, momHCALout);
0432
0433         particle.addPoint(hcalPtout);
0434
0435         const RawParticle& rpho = fst.hoEntrance();
0436
0437         math::XYZPoint posHOEntrance(rpho.x(), rpho.y(), rpho.z());
0438         math::XYZTLorentzVector momHOEntrance(rpho.px(), rpho.py(), rpho.pz(), rpho.e());
0439         reco::PFTrajectoryPoint hoPtin(-1, reco::PFTrajectoryPoint::HOLayer, posHOEntrance, momHOEntrance);
0440
0441         particle.addPoint(hoPtin);
0442
0443       } else {  // add a dummy point
0444         reco::PFTrajectoryPoint dummy;
0445         particle.addPoint(dummy);
0446       }
0447
0448       pOutputPFSimParticleCollection->push_back(particle);
0449     }
0450
0451     iEvent.put(std::move(pOutputPFSimParticleCollection));
0452   }
0453
0454   LogDebug("PFSimParticleProducer") << "STOP event: " << iEvent.id().event() << " in run " << iEvent.id().run() << endl;
0455 }
0456
0457 void PFSimParticleProducer::getSimIDs(const TrackHandle& trackh, std::vector<unsigned>& recTrackSimID) {
0458   if (trackh.isValid()) {
0459     for (unsigned i = 0; i < trackh->size(); i++) {
0460       const reco::PFRecTrackRef ref(trackh, i);
0461
0462       for (auto const& hit : ref->trackRef()->recHits()) {
0463         if (hit->isValid()) {
0464           auto rechit = dynamic_cast<const FastTrackerRecHit*>(hit);
0465
0466           for (unsigned int st_index = 0; st_index < rechit->nSimTrackIds(); ++st_index) {
0467             recTrackSimID.push_back(rechit->simTrackId(st_index));
0468           }
0469           break;
0470         }
0471       }  //loop track rechit
0472     }    //loop recTracks
0473   }      //track handle valid
0474 }