plugins/iso_deposit_extractors/EgammaTrackExtractor.cc

0001 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0002 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0003 #include "DataFormats/Candidate/interface/Candidate.h"
0004 #include "DataFormats/Common/interface/Handle.h"
0005 #include "DataFormats/Common/interface/View.h"
0006 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0007 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
0008 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0009 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0010 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0011 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
0012 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
0013 #include "DataFormats/RecoCandidate/interface/IsoDepositDirection.h"
0014 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0015 #include "DataFormats/TrackReco/interface/Track.h"
0016 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0017 #include "FWCore/Framework/interface/ConsumesCollector.h"
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0020 #include "PhysicsTools/IsolationAlgos/interface/IsoDepositExtractor.h"
0021 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaRecHitIsolation.h"
0022 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaTrackSelector.h"
0023
0024 #include <string>
0025 #include <vector>
0026
0027 namespace egammaisolation {
0028
0029   class EgammaTrackExtractor : public reco::isodeposit::IsoDepositExtractor {
0030   public:
0031     EgammaTrackExtractor(){};
0032     EgammaTrackExtractor(const edm::ParameterSet& par, edm::ConsumesCollector&& iC) : EgammaTrackExtractor(par, iC) {}
0033     EgammaTrackExtractor(const edm::ParameterSet& par, edm::ConsumesCollector& iC);
0034
0035     ~EgammaTrackExtractor() override {}
0036
0037     void fillVetos(const edm::Event& ev, const edm::EventSetup& evSetup, const reco::TrackCollection& track) override {}
0038
0039     virtual reco::IsoDeposit::Vetos vetos(const edm::Event& ev,
0040                                           const edm::EventSetup& evSetup,
0041                                           const reco::Track& track) const;
0042
0043     reco::IsoDeposit deposit(const edm::Event& ev,
0044                              const edm::EventSetup& evSetup,
0045                              const reco::Track& muon) const override {
0046       edm::LogWarning("EgammaIsolationAlgos|EgammaTrackExtractor")
0047           << "This Function is not implemented, bad IsoDeposit Returned";
0048       return reco::IsoDeposit(reco::isodeposit::Direction(1, 1));
0049     }
0050
0051     reco::IsoDeposit deposit(const edm::Event& ev,
0052                              const edm::EventSetup& evSetup,
0053                              const reco::Candidate& muon) const override;
0054
0055   private:
0056     reco::IsoDeposit::Veto veto(const reco::IsoDeposit::Direction& dir) const;
0057
0058   private:
0059     // Parameter set
0060     edm::EDGetTokenT<edm::View<reco::Track> > theTrackCollectionToken;  //! Track Collection Label
0061     std::string theDepositLabel;                                        //! name for deposit
0062     double minCandEt_;                                                  //! minimum candidate et
0063     double theDiff_r;                                                   //! transverse distance to vertex
0064     double theDiff_z;                                                   //! z distance to vertex
0065     double theDR_Max;                                                   //! Maximum cone angle for deposits
0066     double theDR_Veto;                                                  //! Veto cone angle
0067     std::string theBeamlineOption;                                      //! "NONE", "BeamSpotFromEvent"
0068     edm::InputTag barrelEcalHitsTag_;
0069     edm::InputTag endcapEcalHitsTag_;
0070     edm::EDGetTokenT<reco::BeamSpot> theBeamSpotToken;  //! BeamSpot name
0071     unsigned int theNHits_Min;                          //! trk.numberOfValidHits >= theNHits_Min
0072     double theChi2Ndof_Max;                             //! trk.normalizedChi2 < theChi2Ndof_Max
0073     double theChi2Prob_Min;                             //! ChiSquaredProbability(trk.chi2,trk.ndof) > theChi2Prob_Min
0074     double thePt_Min;                                   //! min track pt to include into iso deposit
0075     std::vector<double> paramForIsolBarrel_;  //! Barrel requirements to determine if isolated for selective filling
0076     std::vector<double> paramForIsolEndcap_;  //! Endcap requirements to determine if isolated for selective filling
0077     std::string dzOptionString;
0078     int dzOption;
0079   };
0080
0081 }  // namespace egammaisolation
0082
0083 #include "FWCore/Framework/interface/MakerMacros.h"
0084 #include "PhysicsTools/IsolationAlgos/interface/IsoDepositExtractorFactory.h"
0085 DEFINE_EDM_PLUGIN(IsoDepositExtractorFactory, egammaisolation::EgammaTrackExtractor, "EgammaTrackExtractor");
0086
0087 using namespace edm;
0088 using namespace std;
0089 using namespace reco;
0090 using namespace egammaisolation;
0091 using reco::isodeposit::Direction;
0092
0093 EgammaTrackExtractor::EgammaTrackExtractor(const ParameterSet& par, edm::ConsumesCollector& iC)
0094     : theTrackCollectionToken(iC.consumes<View<Track> >(par.getParameter<edm::InputTag>("inputTrackCollection"))),
0095       theDepositLabel(par.getUntrackedParameter<std::string>("DepositLabel")),
0096       theDiff_r(par.getParameter<double>("Diff_r")),
0097       theDiff_z(par.getParameter<double>("Diff_z")),
0098       theDR_Max(par.getParameter<double>("DR_Max")),
0099       theDR_Veto(par.getParameter<double>("DR_Veto")),
0100       theBeamlineOption(par.getParameter<std::string>("BeamlineOption")),
0101       theBeamSpotToken(iC.mayConsume<reco::BeamSpot>(par.getParameter<edm::InputTag>("BeamSpotLabel"))),
0102       theNHits_Min(par.getParameter<unsigned int>("NHits_Min")),
0103       theChi2Ndof_Max(par.getParameter<double>("Chi2Ndof_Max")),
0104       theChi2Prob_Min(par.getParameter<double>("Chi2Prob_Min")),
0105       thePt_Min(par.getParameter<double>("Pt_Min")),
0106       dzOptionString(par.getParameter<std::string>("dzOption")) {
0107   if (!dzOptionString.compare("dz"))
0108     dzOption = EgammaTrackSelector::dz;
0109   else if (!dzOptionString.compare("vz"))
0110     dzOption = EgammaTrackSelector::vz;
0111   else if (!dzOptionString.compare("bs"))
0112     dzOption = EgammaTrackSelector::bs;
0113   else if (!dzOptionString.compare("vtx"))
0114     dzOption = EgammaTrackSelector::vtx;
0115   else
0116     dzOption = EgammaTrackSelector::dz;
0117 }
0118
0119 reco::IsoDeposit::Vetos EgammaTrackExtractor::vetos(const edm::Event& ev,
0120                                                     const edm::EventSetup& evSetup,
0121                                                     const reco::Track& track) const {
0122   reco::isodeposit::Direction dir(track.eta(), track.phi());
0123   return reco::IsoDeposit::Vetos(1, veto(dir));
0124 }
0125
0126 reco::IsoDeposit::Veto EgammaTrackExtractor::veto(const reco::IsoDeposit::Direction& dir) const {
0127   reco::IsoDeposit::Veto result;
0128   result.vetoDir = dir;
0129   result.dR = theDR_Veto;
0130   return result;
0131 }
0132
0133 IsoDeposit EgammaTrackExtractor::deposit(const Event& event,
0134                                          const EventSetup& eventSetup,
0135                                          const Candidate& candTk) const {
0136   static const std::string metname = "EgammaIsolationAlgos|EgammaTrackExtractor";
0137
0138   reco::isodeposit::Direction candDir;
0139   double dzCut = 0;
0140
0141   reco::TrackBase::Point beamPoint(0, 0, 0);
0142   if (theBeamlineOption == "BeamSpotFromEvent") {
0143     //pick beamSpot
0144     auto beamSpotH = event.getHandle(theBeamSpotToken);
0145
0146     if (beamSpotH.isValid()) {
0147       beamPoint = beamSpotH->position();
0148     }
0149   }
0150
0151   if (candTk.isElectron()) {
0152     const reco::GsfElectron* elec = dynamic_cast<const reco::GsfElectron*>(&candTk);
0153     candDir = reco::isodeposit::Direction(elec->gsfTrack()->eta(), elec->gsfTrack()->phi());
0154   } else {
0155     candDir = reco::isodeposit::Direction(candTk.eta(), candTk.phi());
0156   }
0157
0158   IsoDeposit deposit(candDir);
0159   deposit.setVeto(veto(candDir));
0160   deposit.addCandEnergy(candTk.et());
0161
0162   for (auto const& itrTr : event.get(theTrackCollectionToken)) {
0163     if (candDir.deltaR(reco::isodeposit::Direction(itrTr.eta(), itrTr.phi())) > theDR_Max)
0164       continue;
0165
0166     if (itrTr.normalizedChi2() > theChi2Ndof_Max)
0167       continue;
0168
0169     if (itrTr.pt() < thePt_Min)
0170       continue;
0171
0172     if (theChi2Prob_Min > 0 && ChiSquaredProbability(itrTr.chi2(), itrTr.ndof()) < theChi2Prob_Min)
0173       continue;
0174
0175     if (theNHits_Min > 0 && itrTr.numberOfValidHits() < theNHits_Min)
0176       continue;
0177
0178     if (candTk.isElectron()) {
0179       const reco::GsfElectron* elec = dynamic_cast<const reco::GsfElectron*>(&candTk);
0180       switch (dzOption) {
0181         case EgammaTrackSelector::dz:
0182           dzCut = elec->gsfTrack()->dz() - itrTr.dz();
0183           break;
0184         case EgammaTrackSelector::vz:
0185           dzCut = elec->gsfTrack()->vz() - itrTr.vz();
0186           break;
0187         case EgammaTrackSelector::bs:
0188           dzCut = elec->gsfTrack()->dz(beamPoint) - itrTr.dz(beamPoint);
0189           break;
0190         case EgammaTrackSelector::vtx:
0191           dzCut = itrTr.dz(elec->gsfTrack()->vertex());
0192           break;
0193         default:
0194           dzCut = elec->gsfTrack()->vz() - itrTr.vz();
0195           break;
0196       }
0197     } else {
0198       switch (dzOption) {
0199         case EgammaTrackSelector::dz:
0200           dzCut = itrTr.dz() - candTk.vertex().z();
0201           break;
0202         case EgammaTrackSelector::vz:
0203           dzCut = itrTr.vz() - candTk.vertex().z();
0204           break;
0205         case EgammaTrackSelector::bs:
0206           dzCut = itrTr.dz(beamPoint) - candTk.vertex().z();
0207           break;
0208         case EgammaTrackSelector::vtx:
0209           dzCut = itrTr.dz(candTk.vertex());
0210           break;
0211         default:
0212           dzCut = itrTr.vz() - candTk.vertex().z();
0213           break;
0214       }
0215     }
0216
0217     if (fabs(dzCut) > theDiff_z)
0218       continue;
0219
0220     if (fabs(itrTr.dxy(beamPoint)) > theDiff_r)
0221       continue;
0222
0223     deposit.addDeposit(reco::isodeposit::Direction(itrTr.eta(), itrTr.phi()), itrTr.pt());
0224   }
0225
0226   return deposit;
0227 }