Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-08-23 03:25:28

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 }