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