PFTracking/plugins/PFTrackProducer.cc

0001 /// \brief Abstract
0002 /*!
0003 \author Daniele Benedetti
0004 \date November 2010
0005   PFTrackTransformer transforms all the tracks in the PFRecTracks.
0006   NOTE the PFRecTrack collection is transient.
0007 */
0008
0009 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
0010 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0011 #include "DataFormats/MuonReco/interface/Muon.h"
0012 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0013 #include "DataFormats/ParticleFlowReco/interface/PFRecTrack.h"
0014 #include "DataFormats/TrackReco/interface/Track.h"
0015 #include "DataFormats/VertexReco/interface/Vertex.h"
0016 #include "FWCore/Framework/interface/ESHandle.h"
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "FWCore/Framework/interface/stream/EDProducer.h"
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020 #include "MagneticField/Engine/interface/MagneticField.h"
0021 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0022 #include "RecoParticleFlow/PFTracking/interface/PFTrackTransformer.h"
0023 #include "TrackingTools/IPTools/interface/IPTools.h"
0024 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0025 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0026 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0027
0028 class PFTrackProducer : public edm::stream::EDProducer<> {
0029 public:
0030   ///Constructor
0031   explicit PFTrackProducer(const edm::ParameterSet&);
0032
0033   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0034
0035 private:
0036   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0037   void endRun(const edm::Run&, const edm::EventSetup&) override;
0038
0039   ///Produce the PFRecTrack collection
0040   void produce(edm::Event&, const edm::EventSetup&) override;
0041
0042   const edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> transientTrackToken_;
0043   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldToken_;
0044
0045   ///PFTrackTransformer
0046   std::unique_ptr<PFTrackTransformer> pfTransformer_;
0047   std::vector<edm::EDGetTokenT<reco::TrackCollection>> tracksContainers_;
0048   std::vector<edm::EDGetTokenT<std::vector<Trajectory>>> trajContainers_;
0049   edm::EDGetTokenT<reco::GsfTrackCollection> gsfTrackLabel_;
0050   edm::EDGetTokenT<reco::MuonCollection> muonColl_;
0051   edm::EDGetTokenT<reco::VertexCollection> vtx_h;
0052   ///TRACK QUALITY
0053   bool useQuality_;
0054   reco::TrackBase::TrackQuality trackQuality_;
0055   bool trajinev_;
0056   bool gsfinev_;
0057 };
0058
0059 #include "FWCore/Framework/interface/MakerMacros.h"
0060 DEFINE_FWK_MODULE(PFTrackProducer);
0061
0062 void PFTrackProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0063   edm::ParameterSetDescription desc;
0064   desc.add<std::string>("TrackQuality", "highPurity");
0065   desc.add<bool>("UseQuality", true);
0066   desc.add<edm::InputTag>("GsfTrackModuleLabel", {"electronGsfTracks"});
0067   desc.add<std::vector<edm::InputTag>>("TkColList", {edm::InputTag("generalTracks")});
0068   desc.add<edm::InputTag>("PrimaryVertexLabel", {"offlinePrimaryVertices"});
0069   desc.add<edm::InputTag>("MuColl", {"muons1stStep"});
0070   desc.add<bool>("TrajInEvents", false);
0071   desc.add<bool>("GsfTracksInEvents", true);
0072   descriptions.add("pfTrack", desc);
0073 }
0074
0075 using namespace std;
0076 using namespace edm;
0077 using namespace reco;
0078 PFTrackProducer::PFTrackProducer(const ParameterSet& iConfig)
0079     : transientTrackToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0080       magneticFieldToken_(esConsumes<edm::Transition::BeginRun>()),
0081       pfTransformer_() {
0082   produces<reco::PFRecTrackCollection>();
0083
0084   std::vector<InputTag> tags = iConfig.getParameter<vector<InputTag>>("TkColList");
0085   trajinev_ = iConfig.getParameter<bool>("TrajInEvents");
0086   tracksContainers_.reserve(tags.size());
0087   if (trajinev_) {
0088     trajContainers_.reserve(tags.size());
0089   }
0090   for (auto const& tag : tags) {
0091     tracksContainers_.push_back(consumes<reco::TrackCollection>(tag));
0092     if (trajinev_) {
0093       trajContainers_.push_back(consumes<std::vector<Trajectory>>(tag));
0094     }
0095   }
0096
0097   useQuality_ = iConfig.getParameter<bool>("UseQuality");
0098
0099   gsfinev_ = iConfig.getParameter<bool>("GsfTracksInEvents");
0100   if (gsfinev_) {
0101     gsfTrackLabel_ = consumes<reco::GsfTrackCollection>(iConfig.getParameter<InputTag>("GsfTrackModuleLabel"));
0102   }
0103
0104   trackQuality_ = reco::TrackBase::qualityByName(iConfig.getParameter<std::string>("TrackQuality"));
0105
0106   muonColl_ = consumes<reco::MuonCollection>(iConfig.getParameter<InputTag>("MuColl"));
0107
0108   vtx_h = consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("PrimaryVertexLabel"));
0109 }
0110
0111 void PFTrackProducer::produce(Event& iEvent, const EventSetup& iSetup) {
0112   //create the empty collections
0113   auto PfTrColl = std::make_unique<reco::PFRecTrackCollection>();
0114
0115   //read track collection
0116   Handle<GsfTrackCollection> gsftrackcoll;
0117   bool foundgsf = false;
0118   if (gsfinev_) {
0119     foundgsf = iEvent.getByToken(gsfTrackLabel_, gsftrackcoll);
0120   }
0121
0122   //Get PV for STIP calculation, if there is none then take the dummy
0123   Handle<reco::VertexCollection> vertex;
0124   iEvent.getByToken(vtx_h, vertex);
0125   reco::Vertex dummy;
0126   const reco::Vertex* pv = &dummy;
0127   if (vertex.isValid()) {
0128     pv = &*vertex->begin();
0129   } else {  // create a dummy PV
0130     reco::Vertex::Error e;
0131     e(0, 0) = 0.0015 * 0.0015;
0132     e(1, 1) = 0.0015 * 0.0015;
0133     e(2, 2) = 15. * 15.;
0134     reco::Vertex::Point p(0, 0, 0);
0135     dummy = reco::Vertex(p, e, 0, 0, 0);
0136   }
0137
0138   //setup transient track builder
0139   TransientTrackBuilder const& thebuilder = iSetup.getData(transientTrackToken_);
0140
0141   // read muon collection
0142   Handle<reco::MuonCollection> recMuons;
0143   iEvent.getByToken(muonColl_, recMuons);
0144
0145   //default value for when trajinev_ is false
0146   const vector<Trajectory> dummyTj(0);
0147
0148   for (unsigned int istr = 0; istr < tracksContainers_.size(); istr++) {
0149     //Track collection
0150     Handle<reco::TrackCollection> tkRefCollection;
0151     iEvent.getByToken(tracksContainers_[istr], tkRefCollection);
0152     reco::TrackCollection Tk = *(tkRefCollection.product());
0153
0154     //Use a pointer to aoid unnecessary copying of the collection
0155     const vector<Trajectory>* Tj = &dummyTj;
0156     if (trajinev_) {
0157       //Trajectory collection
0158       Handle<vector<Trajectory>> tjCollection;
0159       iEvent.getByToken(trajContainers_[istr], tjCollection);
0160
0161       Tj = tjCollection.product();
0162     }
0163
0164     for (unsigned int i = 0; i < Tk.size(); i++) {
0165       reco::TrackRef trackRef(tkRefCollection, i);
0166
0167       if (useQuality_ && (!(Tk[i].quality(trackQuality_)))) {
0168         bool isMuCandidate = false;
0169
0170         //TrackRef trackRef(tkRefCollection, i);
0171
0172         if (recMuons.isValid()) {
0173           for (unsigned j = 0; j < recMuons->size(); j++) {
0174             reco::MuonRef muonref(recMuons, j);
0175             if (muonref->track().isNonnull())
0176               if (muonref->track() == trackRef && muonref->isGlobalMuon()) {
0177                 isMuCandidate = true;
0178                 //cout<<" SAVING TRACK "<<endl;
0179                 break;
0180               }
0181           }
0182         }
0183         if (!isMuCandidate) {
0184           continue;
0185         }
0186       }
0187
0188       // find the pre-id kf track
0189       bool preId = false;
0190       if (foundgsf) {
0191         //NOTE: foundgsf is only true if gsftrackcoll is valid
0192         for (auto const& gsfTrack : *gsftrackcoll) {
0193           if (gsfTrack.seedRef().isNull())
0194             continue;
0195           auto const& seed = *(gsfTrack.extra()->seedRef());
0196           auto const& ElSeed = dynamic_cast<ElectronSeed const&>(seed);
0197           if (ElSeed.ctfTrack().isNonnull()) {
0198             if (ElSeed.ctfTrack() == trackRef) {
0199               preId = true;
0200               break;
0201             }
0202           }
0203         }
0204       }
0205       if (preId) {
0206         // Set PFRecTrack of type KF_ElCAND
0207         reco::PFRecTrack pftrack(trackRef->charge(), reco::PFRecTrack::KF_ELCAND, i, trackRef);
0208
0209         bool valid = false;
0210         if (trajinev_) {
0211           valid = pfTransformer_->addPoints(pftrack, *trackRef, (*Tj)[i]);
0212         } else {
0213           Trajectory FakeTraj;
0214           valid = pfTransformer_->addPoints(pftrack, *trackRef, FakeTraj);
0215         }
0216         if (valid) {
0217           //calculate STIP
0218           double stip = -999;
0219           const reco::PFTrajectoryPoint& atECAL = pftrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance);
0220           if (atECAL.isValid())  //if track extrapolates to ECAL
0221           {
0222             GlobalVector direction(pftrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance).position().x(),
0223                                    pftrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance).position().y(),
0224                                    pftrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance).position().z());
0225             stip = IPTools::signedTransverseImpactParameter(thebuilder.build(*trackRef), direction, *pv)
0226                        .second.significance();
0227           }
0228           pftrack.setSTIP(stip);
0229           PfTrColl->push_back(pftrack);
0230         }
0231       } else {
0232         reco::PFRecTrack pftrack(trackRef->charge(), reco::PFRecTrack::KF, i, trackRef);
0233         bool valid = false;
0234         if (trajinev_) {
0235           valid = pfTransformer_->addPoints(pftrack, *trackRef, (*Tj)[i]);
0236         } else {
0237           Trajectory FakeTraj;
0238           valid = pfTransformer_->addPoints(pftrack, *trackRef, FakeTraj);
0239         }
0240
0241         if (valid) {
0242           double stip = -999;
0243           const reco::PFTrajectoryPoint& atECAL = pftrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance);
0244           if (atECAL.isValid()) {
0245             GlobalVector direction(pftrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance).position().x(),
0246                                    pftrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance).position().y(),
0247                                    pftrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance).position().z());
0248             stip = IPTools::signedTransverseImpactParameter(thebuilder.build(*trackRef), direction, *pv)
0249                        .second.significance();
0250           }
0251           pftrack.setSTIP(stip);
0252           PfTrColl->push_back(pftrack);
0253         }
0254       }
0255     }
0256   }
0257   iEvent.put(std::move(PfTrColl));
0258 }
0259
0260 // ------------ method called once each job just before starting event loop  ------------
0261 void PFTrackProducer::beginRun(const edm::Run& run, const EventSetup& iSetup) {
0262   auto const& magneticField = iSetup.getData(magneticFieldToken_);
0263   pfTransformer_ = std::make_unique<PFTrackTransformer>(math::XYZVector(magneticField.inTesla(GlobalPoint(0, 0, 0))));
0264   if (!trajinev_)
0265     pfTransformer_->OnlyProp();
0266 }
0267
0268 // ------------ method called once each job just after ending the event loop  ------------
0269 void PFTrackProducer::endRun(const edm::Run& run, const EventSetup& iSetup) { pfTransformer_.reset(); }