PFTrackProducer.cc

Go to the documentation of this file.

#include <memory>
#include "RecoParticleFlow/PFTracking/interface/PFTrackProducer.h"
#include "RecoParticleFlow/PFTracking/interface/PFTrackTransformer.h"
#include "DataFormats/ParticleFlowReco/interface/PFRecTrack.h"
#include "DataFormats/ParticleFlowReco/interface/PFRecTrackFwd.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
#include "DataFormats/EgammaReco/interface/ElectronSeed.h"
#include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h"
#include "TrackingTools/PatternTools/interface/Trajectory.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "MagneticField/Engine/interface/MagneticField.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "DataFormats/MuonReco/interface/Muon.h"
#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
#include "TrackingTools/Records/interface/TransientTrackRecord.h"
#include "TrackingTools/IPTools/interface/IPTools.h"
#include "DataFormats/VertexReco/interface/Vertex.h"

using namespace std;
using namespace edm;
using namespace reco;
PFTrackProducer::PFTrackProducer(const ParameterSet& iConfig):
  pfTransformer_()
{
  produces<reco::PFRecTrackCollection>();
  
  
  std::vector<InputTag> tags=iConfig.getParameter< vector < InputTag > >("TkColList");
  trajinev_ = iConfig.getParameter<bool>("TrajInEvents");
  tracksContainers_.reserve(tags.size());
  if(trajinev_) { trajContainers_.reserve(tags.size()); }
  for( auto const& tag: tags) {
    tracksContainers_.push_back(consumes<reco::TrackCollection>(tag));
    if(trajinev_) {
      trajContainers_.push_back(consumes<std::vector<Trajectory> >(tag));
    }
  }

  useQuality_   = iConfig.getParameter<bool>("UseQuality");
  
  gsfinev_ = iConfig.getParameter<bool>("GsfTracksInEvents");
  if(gsfinev_) {
    gsfTrackLabel_ = consumes<reco::GsfTrackCollection>(iConfig.getParameter<InputTag>
                                                        ("GsfTrackModuleLabel"));  
  }

  trackQuality_=reco::TrackBase::qualityByName(iConfig.getParameter<std::string>("TrackQuality"));
  
  muonColl_ = consumes<reco::MuonCollection>(iConfig.getParameter< InputTag >("MuColl"));
  
  
  vtx_h=consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("PrimaryVertexLabel"));
  
}

void
PFTrackProducer::produce(Event& iEvent, const EventSetup& iSetup)
{
  
  //create the empty collections 
  auto PfTrColl = std::make_unique<reco::PFRecTrackCollection>();
  
  //read track collection
  Handle<GsfTrackCollection> gsftrackcoll;
  bool foundgsf = false;
  if(gsfinev_) {
    foundgsf = iEvent.getByToken(gsfTrackLabel_,gsftrackcoll);
  }

  //Get PV for STIP calculation, if there is none then take the dummy  
  Handle<reco::VertexCollection> vertex;
  iEvent.getByToken(vtx_h, vertex);
  reco::Vertex dummy;
  const reco::Vertex* pv=&dummy;  
  if (vertex.isValid()) 
    {
      pv = &*vertex->begin();    
    } 
  else 
    { // create a dummy PV
      reco::Vertex::Error e;
      e(0, 0) = 0.0015 * 0.0015;
      e(1, 1) = 0.0015 * 0.0015;
      e(2, 2) = 15. * 15.;
      reco::Vertex::Point p(0, 0, 0);
      dummy = reco::Vertex(p, e, 0, 0, 0);   
    } 
  
  edm::ESHandle<TransientTrackBuilder> builder;
  iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder", builder);
  TransientTrackBuilder thebuilder = *(builder.product());
  
  // read muon collection
  Handle< reco::MuonCollection > recMuons;
  iEvent.getByToken(muonColl_, recMuons);
  
  //default value for when trajinev_ is false
  const vector<Trajectory> dummyTj(0);
  
  for (unsigned int istr=0; istr<tracksContainers_.size();istr++){
    
    //Track collection
    Handle<reco::TrackCollection> tkRefCollection;
    iEvent.getByToken(tracksContainers_[istr], tkRefCollection);
    reco::TrackCollection  Tk=*(tkRefCollection.product());
    
    //Use a pointer to aoid unnecessary copying of the collection
    const vector<Trajectory>* Tj = &dummyTj;
    if(trajinev_) {
      //Trajectory collection
      Handle<vector<Trajectory> > tjCollection;
      iEvent.getByToken(trajContainers_[istr], tjCollection);
    
      Tj =tjCollection.product();
    }
    
    
    for(unsigned int i=0;i<Tk.size();i++){
      
      reco::TrackRef trackRef(tkRefCollection, i);
    
      if (useQuality_ &&
      (!(Tk[i].quality(trackQuality_)))){
    
    bool isMuCandidate = false;
    
    //TrackRef trackRef(tkRefCollection, i);
    
    if(recMuons.isValid() ) {
      for(unsigned j=0;j<recMuons->size(); j++) {
        reco::MuonRef muonref( recMuons, j );
        if (muonref->track().isNonnull()) 
          if( muonref->track() == trackRef && muonref->isGlobalMuon()){
        isMuCandidate=true;
        //cout<<" SAVING TRACK "<<endl;
        break;
          }
      }
    }
    if(!isMuCandidate)
      {
        continue;     
      }
    
      }
           
      // find the pre-id kf track
      bool preId = false;
      if(foundgsf) {
        //NOTE: foundgsf is only true if gsftrackcoll is valid
    for (auto const& gsfTrack: *gsftrackcoll) {
      if (gsfTrack.seedRef().isNull()) continue;
          auto const& seed = *(gsfTrack.extra()->seedRef());
          auto const& ElSeed = dynamic_cast<ElectronSeed const&>(seed);
      if (ElSeed.ctfTrack().isNonnull()) {
        if(ElSeed.ctfTrack() == trackRef) {
              preId = true;
              break;
            }
      }
    }
      }
      if(preId) {
    // Set PFRecTrack of type KF_ElCAND
    reco::PFRecTrack pftrack( trackRef->charge(), 
                  reco::PFRecTrack::KF_ELCAND, 
                  i, trackRef );
    
    bool valid = false;
    if(trajinev_) {
      valid = pfTransformer_->addPoints( pftrack, *trackRef, (*Tj)[i]);
    }
    else {
      Trajectory FakeTraj;
      valid = pfTransformer_->addPoints( pftrack, *trackRef, FakeTraj);
    }
    if(valid) {
      //calculate STIP
      double stip=-999;
      const reco::PFTrajectoryPoint& atECAL=pftrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance);
      if(atECAL.isValid()) //if track extrapolates to ECAL
        {
          GlobalVector direction(pftrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance).position().x(),
                     pftrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance).position().y(), 
                     pftrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance).position().z());
          stip = IPTools::signedTransverseImpactParameter(thebuilder.build(*trackRef), direction, *pv).second.significance();
        }
      pftrack.setSTIP(stip);
      PfTrColl->push_back(pftrack);
    }       
      }
      else {
    reco::PFRecTrack pftrack( trackRef->charge(), 
                  reco::PFRecTrack::KF, 
                  i, trackRef );
    bool valid = false;
    if(trajinev_) {
      valid = pfTransformer_->addPoints( pftrack, *trackRef, (*Tj)[i]);
    }
    else {
      Trajectory FakeTraj;
      valid = pfTransformer_->addPoints( pftrack, *trackRef, FakeTraj);
    }
    
    if(valid) {
      double stip=-999;
      const reco::PFTrajectoryPoint& atECAL=pftrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance);
      if(atECAL.isValid())
        {
          GlobalVector direction(pftrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance).position().x(),
                     pftrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance).position().y(), 
                     pftrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance).position().z());
          stip = IPTools::signedTransverseImpactParameter(thebuilder.build(*trackRef), direction, *pv).second.significance();
        }
      pftrack.setSTIP(stip);
      PfTrColl->push_back(pftrack);
    }
      }
    }
  }
  iEvent.put(std::move(PfTrColl));
}

// ------------ method called once each job just before starting event loop  ------------
void 
PFTrackProducer::beginRun(const edm::Run& run,
              const EventSetup& iSetup)
{
  ESHandle<MagneticField> magneticField;
  iSetup.get<IdealMagneticFieldRecord>().get(magneticField);
  pfTransformer_.reset( new PFTrackTransformer(math::XYZVector(magneticField->inTesla(GlobalPoint(0,0,0)))) );
  if(!trajinev_)
    pfTransformer_->OnlyProp();
}

// ------------ method called once each job just after ending the event loop  ------------
void 
PFTrackProducer::endRun(const edm::Run& run,
            const EventSetup& iSetup) {
  pfTransformer_.reset();
}