MuonIsolation/plugins/PixelTrackExtractor.cc

0001 #include "PixelTrackExtractor.h"
0002
0003 #include "RecoMuon/MuonIsolation/interface/Range.h"
0004 #include "DataFormats/RecoCandidate/interface/IsoDepositDirection.h"
0005 #include "TrackSelector.h"
0006 #include "DataFormats/Common/interface/Handle.h"
0007 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0010
0011 using namespace edm;
0012 using namespace std;
0013 using namespace reco;
0014 using namespace muonisolation;
0015 using reco::isodeposit::Direction;
0016
0017 PixelTrackExtractor::PixelTrackExtractor(const ParameterSet& par, edm::ConsumesCollector&& iC)
0018     : theFieldToken(iC.esConsumes()),
0019       theTrackCollectionToken(iC.consumes<TrackCollection>(par.getParameter<edm::InputTag>("inputTrackCollection"))),
0020       theDepositLabel(par.getUntrackedParameter<string>("DepositLabel")),
0021       theDiff_r(par.getParameter<double>("Diff_r")),
0022       theDiff_z(par.getParameter<double>("Diff_z")),
0023       theDR_Max(par.getParameter<double>("DR_Max")),
0024       theDR_Veto(par.getParameter<double>("DR_Veto")),
0025       theBeamlineOption(par.getParameter<string>("BeamlineOption")),
0026       theBeamSpotToken(iC.mayConsume<BeamSpot>(par.getParameter<edm::InputTag>("BeamSpotLabel"))),
0027       theNHits_Min(par.getParameter<unsigned int>("NHits_Min")),
0028       theChi2Ndof_Max(par.getParameter<double>("Chi2Ndof_Max")),
0029       theChi2Prob_Min(par.getParameter<double>("Chi2Prob_Min")),
0030       thePt_Min(par.getParameter<double>("Pt_Min")),
0031       thePropagateTracksToRadius(par.getParameter<bool>("PropagateTracksToRadius")),
0032       theReferenceRadius(par.getParameter<double>("ReferenceRadius")),
0033       theVetoLeadingTrack(par.getParameter<bool>("VetoLeadingTrack")),  //! will veto leading track if
0034       thePtVeto_Min(par.getParameter<double>("PtVeto_Min")),            //! .. it is above this threshold
0035       theDR_VetoPt(par.getParameter<double>("DR_VetoPt"))               //!.. and is inside this cone
0036 {}
0037
0038 reco::IsoDeposit::Vetos PixelTrackExtractor::vetos(const edm::Event& ev,
0039                                                    const edm::EventSetup& evSetup,
0040                                                    const reco::Track& track) const {
0041   Direction dir(track.eta(), track.phi());
0042   return reco::IsoDeposit::Vetos(1, veto(dir));
0043 }
0044
0045 reco::IsoDeposit::Veto PixelTrackExtractor::veto(const Direction& dir) const {
0046   reco::IsoDeposit::Veto result;
0047   result.vetoDir = dir;
0048   result.dR = theDR_Veto;
0049   return result;
0050 }
0051
0052 Direction PixelTrackExtractor::directionAtPresetRadius(const Track& tk, double bz) const {
0053   if (!thePropagateTracksToRadius) {
0054     return Direction(tk.eta(), tk.phi());
0055   }
0056
0057   // this should represent a cylinder in global frame at R=refRadius cm, roughly where mid-layer of pixels is
0058   double psRadius = theReferenceRadius;
0059   double tkDxy = tk.dxy();
0060   double s2D = fabs(tk.dxy()) < psRadius ? sqrt(psRadius * psRadius - tkDxy * tkDxy) : 0;
0061
0062   // the field we get from the caller is already in units of GeV/cm
0063   double dPhi = -s2D * tk.charge() * bz / tk.pt();
0064
0065   return Direction(tk.eta(), tk.phi() + dPhi);
0066 }
0067
0068 IsoDeposit PixelTrackExtractor::deposit(const Event& event, const EventSetup& eventSetup, const Track& muon) const {
0069   static const std::string metname = "MuonIsolation|PixelTrackExtractor";
0070
0071   auto const& bField = eventSetup.getData(theFieldToken);
0072   double bz = bField.inInverseGeV(GlobalPoint(0., 0., 0.)).z();
0073
0074   Direction muonDir(directionAtPresetRadius(muon, bz));
0075   IsoDeposit deposit(muonDir);
0076   //! Note, this can be reset below if theVetoLeadingTrack is set and conditions are met
0077   deposit.setVeto(veto(muonDir));
0078
0079   deposit.addCandEnergy(muon.pt());
0080
0081   Handle<TrackCollection> tracksH;
0082   event.getByToken(theTrackCollectionToken, tracksH);
0083   //  const TrackCollection tracks = *(tracksH.product());
0084   LogTrace(metname) << "***** TRACK COLLECTION SIZE: " << tracksH->size();
0085
0086   double vtx_z = muon.vz();
0087   LogTrace(metname) << "***** Muon vz: " << vtx_z;
0088   reco::TrackBase::Point beamPoint(0, 0, 0);
0089
0090   if (theBeamlineOption == "BeamSpotFromEvent") {
0091     //pick beamSpot
0092     reco::BeamSpot beamSpot;
0093     edm::Handle<reco::BeamSpot> beamSpotH;
0094
0095     event.getByToken(theBeamSpotToken, beamSpotH);
0096
0097     if (beamSpotH.isValid()) {
0098       beamPoint = beamSpotH->position();
0099       LogTrace(metname) << "Extracted beam point at " << beamPoint << std::endl;
0100     }
0101   }
0102
0103   LogTrace(metname) << "Using beam point at " << beamPoint << std::endl;
0104
0105   TrackSelector::Parameters pars(
0106       TrackSelector::Range(vtx_z - theDiff_z, vtx_z + theDiff_z), theDiff_r, muonDir, theDR_Max, beamPoint);
0107
0108   pars.nHitsMin = theNHits_Min;
0109   pars.chi2NdofMax = theChi2Ndof_Max;
0110   pars.chi2ProbMin = theChi2Prob_Min;
0111   pars.ptMin = thePt_Min;
0112
0113   TrackSelector selection(pars);
0114   TrackSelector::result_type sel_tracks = selection(*tracksH);
0115   LogTrace(metname) << "all tracks: " << tracksH->size() << " selected: " << sel_tracks.size();
0116
0117   double maxPt = -1;
0118   Direction maxPtDir;
0119   TrackSelector::result_type::const_iterator tkI = sel_tracks.begin();
0120   for (; tkI != sel_tracks.end(); ++tkI) {
0121     const reco::Track* tk = *tkI;
0122     LogTrace(metname) << "This track has: pt= " << tk->pt() << ", eta= " << tk->eta() << ", phi= " << tk->phi();
0123     Direction dirTrk(directionAtPresetRadius(*tk, bz));
0124     deposit.addDeposit(dirTrk, tk->pt());
0125     double tkDr = (muonDir - dirTrk).deltaR;
0126     double tkPt = tk->pt();
0127     if (theVetoLeadingTrack && tkPt > thePtVeto_Min && tkDr < theDR_VetoPt && maxPt < tkPt) {
0128       maxPt = tkPt;
0129       maxPtDir = dirTrk;
0130     }
0131   }
0132   if (maxPt > 0) {
0133     deposit.setVeto(veto(maxPtDir));
0134     LogTrace(metname) << " Set track veto the leading track with pt " << maxPt << " in direction  (eta,phi) "
0135                       << maxPtDir.eta() << ", " << maxPtDir.phi();
0136   }
0137
0138   return deposit;
0139 }