RecoTau/plugins/TauDiscriminationAgainstCaloMuon.cc

0001
0002 /*
0003  * class TauDiscriminationAgainstCaloMuon
0004  * created : Nov 20 2010
0005  * revised :
0006  * Authors : Christian Veelken (UC Davis)
0007  */
0008
0009 #include "RecoTauTag/RecoTau/interface/TauDiscriminationProducerBase.h"
0010
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/Utilities/interface/InputTag.h"
0013
0014 #include <FWCore/ParameterSet/interface/ConfigurationDescriptions.h>
0015 #include <FWCore/ParameterSet/interface/ParameterSetDescription.h>
0016
0017 #include "DataFormats/MuonReco/interface/Muon.h"
0018 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0019
0020 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0021 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0022 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0023 #include "DataFormats/HcalRecHit/interface/HBHERecHit.h"
0024 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0025 #include "DataFormats/CaloTowers/interface/CaloTower.h"
0026 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0027
0028 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0029 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0030 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0031 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateClosestToPoint.h"
0032
0033 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0034 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0035 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0036 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0037 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0038
0039 #include "DataFormats/VertexReco/interface/Vertex.h"
0040 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0041
0042 #include "DataFormats/Math/interface/deltaR.h"
0043
0044 #include <TVector3.h>
0045 #include <TMath.h>
0046
0047 #include <string>
0048
0049 namespace {
0050
0051   using namespace reco;
0052
0053   // define acess to lead. track for CaloTaus
0054   template <typename T>
0055   class TauLeadTrackExtractor {
0056   public:
0057     reco::TrackRef getLeadTrack(const T& tau) const { return tau.leadTrack(); }
0058     double getTrackPtSum(const T& tau) const {
0059       double trackPtSum = 0.;
0060       for (TrackRefVector::const_iterator signalTrack = tau.signalTracks().begin();
0061            signalTrack != tau.signalTracks().end();
0062            ++signalTrack) {
0063         trackPtSum += (*signalTrack)->pt();
0064       }
0065       return trackPtSum;
0066     }
0067   };
0068
0069   // define acess to lead. track for PFTaus
0070   template <>
0071   class TauLeadTrackExtractor<reco::PFTau> {
0072   public:
0073     reco::TrackRef getLeadTrack(const reco::PFTau& tau) const { return tau.leadPFChargedHadrCand()->trackRef(); }
0074     double getTrackPtSum(const reco::PFTau& tau) const {
0075       double trackPtSum = 0.;
0076       for (std::vector<CandidatePtr>::const_iterator signalTrack = tau.signalChargedHadrCands().begin();
0077            signalTrack != tau.signalChargedHadrCands().end();
0078            ++signalTrack) {
0079         trackPtSum += (*signalTrack)->pt();
0080       }
0081       return trackPtSum;
0082     }
0083   };
0084
0085   template <class TauType, class TauDiscriminator>
0086   class TauDiscriminationAgainstCaloMuon final : public TauDiscriminationProducerBase<TauType, TauDiscriminator> {
0087   public:
0088     // setup framework types for this tautype
0089     typedef std::vector<TauType> TauCollection;
0090     typedef edm::Ref<TauCollection> TauRef;
0091
0092     explicit TauDiscriminationAgainstCaloMuon(const edm::ParameterSet&);
0093     ~TauDiscriminationAgainstCaloMuon() override {}
0094
0095     // called at the beginning of every event
0096     void beginEvent(const edm::Event&, const edm::EventSetup&) override;
0097
0098     double discriminate(const TauRef&) const override;
0099
0100     static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0101
0102   private:
0103     edm::InputTag srcEcalRecHitsBarrel_;
0104     edm::Handle<EcalRecHitCollection> ebRecHits_;
0105     edm::InputTag srcEcalRecHitsEndcap_;
0106     edm::Handle<EcalRecHitCollection> eeRecHits_;
0107     edm::InputTag srcHcalRecHits_;
0108     edm::Handle<HBHERecHitCollection> hbheRecHits_;
0109
0110     edm::InputTag srcVertex_;
0111     GlobalPoint eventVertexPosition_;
0112
0113     const TransientTrackBuilder* trackBuilder_;
0114     const CaloGeometry* caloGeometry_;
0115     const edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> trackBuilderToken_;
0116     const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0117
0118     TauLeadTrackExtractor<TauType> leadTrackExtractor_;
0119
0120     double minLeadTrackPt_;
0121     double minLeadTrackPtFraction_;
0122
0123     double drEcal_;
0124     double drHcal_;
0125
0126     double maxEnEcal_;
0127     double maxEnHcal_;
0128
0129     double maxEnToTrackRatio_;
0130   };
0131
0132   template <class TauType, class TauDiscriminator>
0133   TauDiscriminationAgainstCaloMuon<TauType, TauDiscriminator>::TauDiscriminationAgainstCaloMuon(
0134       const edm::ParameterSet& cfg)
0135       : TauDiscriminationProducerBase<TauType, TauDiscriminator>(cfg),
0136         trackBuilderToken_(this->esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))),
0137         caloGeometryToken_(this->esConsumes()) {
0138     srcEcalRecHitsBarrel_ = cfg.getParameter<edm::InputTag>("srcEcalRecHitsBarrel");
0139     srcEcalRecHitsEndcap_ = cfg.getParameter<edm::InputTag>("srcEcalRecHitsEndcap");
0140     srcHcalRecHits_ = cfg.getParameter<edm::InputTag>("srcHcalRecHits");
0141
0142     srcVertex_ = cfg.getParameter<edm::InputTag>("srcVertex");
0143
0144     minLeadTrackPt_ = cfg.getParameter<double>("minLeadTrackPt");
0145     minLeadTrackPtFraction_ = cfg.getParameter<double>("minLeadTrackPtFraction");
0146
0147     drEcal_ = cfg.getParameter<double>("dRecal");
0148     drHcal_ = cfg.getParameter<double>("dRhcal");
0149
0150     maxEnEcal_ = cfg.getParameter<double>("maxEnEcal");
0151     maxEnHcal_ = cfg.getParameter<double>("maxEnHcal");
0152
0153     maxEnToTrackRatio_ = cfg.getParameter<double>("maxEnToTrackRatio");
0154   }
0155
0156   template <class TauType, class TauDiscriminator>
0157   void TauDiscriminationAgainstCaloMuon<TauType, TauDiscriminator>::beginEvent(const edm::Event& evt,
0158                                                                                const edm::EventSetup& evtSetup) {
0159     evt.getByLabel(srcEcalRecHitsBarrel_, ebRecHits_);
0160     evt.getByLabel(srcEcalRecHitsEndcap_, eeRecHits_);
0161     evt.getByLabel(srcHcalRecHits_, hbheRecHits_);
0162
0163     trackBuilder_ = &evtSetup.getData(trackBuilderToken_);
0164     if (!trackBuilder_) {
0165       edm::LogError("TauDiscriminationAgainstCaloMuon::discriminate") << " Failed to access TransientTrackBuilder !!";
0166     }
0167
0168     caloGeometry_ = &evtSetup.getData(caloGeometryToken_);
0169     if (!caloGeometry_) {
0170       edm::LogError("TauDiscriminationAgainstCaloMuon::discriminate") << " Failed to access CaloGeometry !!";
0171     }
0172
0173     edm::Handle<reco::VertexCollection> vertices;
0174     evt.getByLabel(srcVertex_, vertices);
0175     eventVertexPosition_ = GlobalPoint(0., 0., 0.);
0176     if (!vertices->empty()) {
0177       const reco::Vertex& thePrimaryEventVertex = (*vertices->begin());
0178       eventVertexPosition_ =
0179           GlobalPoint(thePrimaryEventVertex.x(), thePrimaryEventVertex.y(), thePrimaryEventVertex.z());
0180     }
0181   }
0182
0183   double compEcalEnergySum(const EcalRecHitCollection& ecalRecHits,
0184                            const CaloSubdetectorGeometry* detGeometry,
0185                            const reco::TransientTrack& transientTrack,
0186                            double dR,
0187                            const GlobalPoint& eventVertexPosition) {
0188     double ecalEnergySum = 0.;
0189     for (EcalRecHitCollection::const_iterator ecalRecHit = ecalRecHits.begin(); ecalRecHit != ecalRecHits.end();
0190          ++ecalRecHit) {
0191       auto cellGeometry = detGeometry->getGeometry(ecalRecHit->detid());
0192
0193       if (!cellGeometry) {
0194         edm::LogError("compEcalEnergySum")
0195             << " Failed to access ECAL geometry for detId = " << ecalRecHit->detid().rawId() << " --> skipping !!";
0196         continue;
0197       }
0198
0199       const GlobalPoint& cellPosition = cellGeometry->getPosition();
0200
0201       //--- CV: speed up computation by requiring eta-phi distance
0202       //        between cell position and track direction to be dR < 0.5
0203       Vector3DBase<float, GlobalTag> cellPositionRelVertex = (cellPosition)-eventVertexPosition;
0204       if (deltaR(cellPositionRelVertex.eta(),
0205                  cellPositionRelVertex.phi(),
0206                  transientTrack.track().eta(),
0207                  transientTrack.track().phi()) > 0.5)
0208         continue;
0209
0210       TrajectoryStateClosestToPoint dcaPosition = transientTrack.trajectoryStateClosestToPoint(cellPosition);
0211
0212       Vector3DBase<float, GlobalTag> d = (cellPosition - dcaPosition.position());
0213
0214       TVector3 d3(d.x(), d.y(), d.z());
0215       TVector3 dir(transientTrack.track().px(), transientTrack.track().py(), transientTrack.track().pz());
0216
0217       double dPerp = d3.Cross(dir.Unit()).Mag();
0218       double dParl = TVector3(cellPosition.x(), cellPosition.y(), cellPosition.z()).Dot(dir.Unit());
0219
0220       if (dPerp < dR && dParl > 100.) {
0221         ecalEnergySum += ecalRecHit->energy();
0222       }
0223     }
0224
0225     return ecalEnergySum;
0226   }
0227
0228   double compHcalEnergySum(const HBHERecHitCollection& hcalRecHits,
0229                            const HcalGeometry* hcGeometry,
0230                            const reco::TransientTrack& transientTrack,
0231                            double dR,
0232                            const GlobalPoint& eventVertexPosition) {
0233     double hcalEnergySum = 0.;
0234     for (HBHERecHitCollection::const_iterator hcalRecHit = hcalRecHits.begin(); hcalRecHit != hcalRecHits.end();
0235          ++hcalRecHit) {
0236       const GlobalPoint cellPosition = hcGeometry->getPosition(hcalRecHit->detid());
0237
0238       //--- CV: speed up computation by requiring eta-phi distance
0239       //        between cell position and track direction to be dR < 0.5
0240       Vector3DBase<float, GlobalTag> cellPositionRelVertex = (cellPosition)-eventVertexPosition;
0241       if (deltaR(cellPositionRelVertex.eta(),
0242                  cellPositionRelVertex.phi(),
0243                  transientTrack.track().eta(),
0244                  transientTrack.track().phi()) > 0.5)
0245         continue;
0246
0247       TrajectoryStateClosestToPoint dcaPosition = transientTrack.trajectoryStateClosestToPoint(cellPosition);
0248
0249       Vector3DBase<float, GlobalTag> d = ((cellPosition)-dcaPosition.position());
0250
0251       TVector3 d3(d.x(), d.y(), d.z());
0252       TVector3 dir(transientTrack.track().px(), transientTrack.track().py(), transientTrack.track().pz());
0253
0254       double dPerp = d3.Cross(dir.Unit()).Mag();
0255       double dParl = TVector3(cellPosition.x(), cellPosition.y(), cellPosition.z()).Dot(dir.Unit());
0256
0257       if (dPerp < dR && dParl > 100.) {
0258         hcalEnergySum += hcalRecHit->energy();
0259       }
0260     }
0261
0262     return hcalEnergySum;
0263   }
0264
0265   template <class TauType, class TauDiscriminator>
0266   double TauDiscriminationAgainstCaloMuon<TauType, TauDiscriminator>::discriminate(const TauRef& tau) const {
0267     if (!(trackBuilder_ && caloGeometry_))
0268       return 0.;
0269
0270     const CaloSubdetectorGeometry* ebGeometry = caloGeometry_->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0271     const CaloSubdetectorGeometry* eeGeometry = caloGeometry_->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0272     const HcalGeometry* hcGeometry =
0273         static_cast<const HcalGeometry*>(caloGeometry_->getSubdetectorGeometry(DetId::Hcal, HcalBarrel));
0274
0275     TrackRef leadTrackRef = leadTrackExtractor_.getLeadTrack(*tau);
0276
0277     if ((leadTrackRef.isAvailable() || leadTrackRef.isTransient()) && leadTrackRef.isNonnull()) {
0278       double leadTrackPt = leadTrackRef->pt();
0279       double trackPtSum = leadTrackExtractor_.getTrackPtSum(*tau);
0280
0281       double leadTrackPtFraction = (trackPtSum > 0.) ? (leadTrackPt / trackPtSum) : -1.;
0282
0283       if (leadTrackPt > minLeadTrackPt_ && leadTrackPtFraction > minLeadTrackPtFraction_) {
0284         reco::TransientTrack transientTrack = trackBuilder_->build(leadTrackRef);
0285
0286         double ebEnergySum = compEcalEnergySum(*ebRecHits_, ebGeometry, transientTrack, drEcal_, eventVertexPosition_);
0287         double eeEnergySum = compEcalEnergySum(*eeRecHits_, eeGeometry, transientTrack, drEcal_, eventVertexPosition_);
0288         double ecalEnergySum = ebEnergySum + eeEnergySum;
0289
0290         double hbheEnergySum =
0291             compHcalEnergySum(*hbheRecHits_, hcGeometry, transientTrack, drHcal_, eventVertexPosition_);
0292
0293         double caloEnergySum = ecalEnergySum + hbheEnergySum;
0294
0295         if (ecalEnergySum < maxEnEcal_ && hbheEnergySum < maxEnHcal_ &&
0296             caloEnergySum < (maxEnToTrackRatio_ * leadTrackPt))
0297           return 0.;
0298       }
0299     }
0300
0301     return 1.;
0302   }
0303
0304 }  // namespace
0305
0306 #include "FWCore/Framework/interface/MakerMacros.h"
0307
0308 typedef TauDiscriminationAgainstCaloMuon<PFTau, PFTauDiscriminator> PFRecoTauDiscriminationAgainstCaloMuon;
0309
0310 // accordingly method for specific class
0311 template <>
0312 void TauDiscriminationAgainstCaloMuon<PFTau, PFTauDiscriminator>::fillDescriptions(
0313     edm::ConfigurationDescriptions& descriptions) {
0314   // pfRecoTauDiscriminationAgainstCaloMuon
0315   edm::ParameterSetDescription desc;
0316   desc.add<edm::InputTag>("srcHcalRecHits", edm::InputTag("hbhereco"));
0317   desc.add<double>("minLeadTrackPt", 15.0);
0318   desc.add<double>("maxEnToTrackRatio", 0.25);
0319   desc.add<edm::InputTag>("srcVertex", edm::InputTag("offlinePrimaryVertices"));
0320   desc.add<edm::InputTag>("PFTauProducer", edm::InputTag("pfRecoTauProducer"));
0321   desc.add<edm::InputTag>("srcEcalRecHitsBarrel", edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0322   desc.add<double>("dRhcal", 25.0);
0323   {
0324     edm::ParameterSetDescription psd0;
0325     psd0.add<std::string>("BooleanOperator", "and");
0326     {
0327       edm::ParameterSetDescription psd1;
0328       psd1.add<double>("cut");
0329       psd1.add<edm::InputTag>("Producer");
0330       psd0.addOptional<edm::ParameterSetDescription>("leadTrack", psd1);
0331     }
0332     desc.add<edm::ParameterSetDescription>("Prediscriminants", psd0);
0333   }
0334   desc.add<double>("maxEnHcal", 8.0);
0335   desc.add<double>("dRecal", 15.0);
0336   desc.add<edm::InputTag>("srcEcalRecHitsEndcap", edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0337   desc.add<double>("minLeadTrackPtFraction", 0.8);
0338   desc.add<double>("maxEnEcal", 3.0);
0339   descriptions.add("pfRecoTauDiscriminationAgainstCaloMuon", desc);
0340 }
0341
0342 DEFINE_FWK_MODULE(PFRecoTauDiscriminationAgainstCaloMuon);