plugins/iso_deposit_extractors/EgammaTowerExtractor.cc

0001 //*****************************************************************************
0002 // File:      EgammaTowerExtractor.cc
0003 // ----------------------------------------------------------------------------
0004 // OrigAuth:  Matthias Mozer
0005 // Institute: IIHE-VUB
0006 //=============================================================================
0007 //*****************************************************************************
0008
0009 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0010 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0011 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0012 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0013 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0014 #include "DataFormats/CaloTowers/interface/CaloTower.h"
0015 #include "DataFormats/Candidate/interface/Candidate.h"
0016 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "FWCore/Framework/interface/EventSetup.h"
0019 #include "FWCore/Framework/interface/ConsumesCollector.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/Utilities/interface/InputTag.h"
0022 #include "DataFormats/TrackReco/interface/Track.h"
0023 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0024 #include "PhysicsTools/IsolationAlgos/interface/IsoDepositExtractor.h"
0025 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
0026
0027 #include <Math/VectorUtil.h>
0028
0029 #include <vector>
0030 #include <functional>
0031 #include <cmath>
0032
0033 namespace egammaisolation {
0034
0035   class EgammaTowerExtractor : public reco::isodeposit::IsoDepositExtractor {
0036   public:
0037     enum HcalDepth { AllDepths = -1, Undefined = 0, Depth1 = 1, Depth2 = 2 };
0038
0039   public:
0040     EgammaTowerExtractor(const edm::ParameterSet& par, edm::ConsumesCollector&& iC) : EgammaTowerExtractor(par, iC) {}
0041     EgammaTowerExtractor(const edm::ParameterSet& par, edm::ConsumesCollector& iC)
0042         : extRadius2_(par.getParameter<double>("extRadius")),
0043           intRadius_(par.getParameter<double>("intRadius")),
0044           etLow_(par.getParameter<double>("etMin")),
0045           caloTowerToken(iC.consumes<CaloTowerCollection>(par.getParameter<edm::InputTag>("caloTowers"))),
0046           depth_(par.getParameter<int>("hcalDepth")) {
0047       extRadius2_ *= extRadius2_;
0048       //lets just check we have a valid depth
0049       //should we throw an exception or just warn and then fail gracefully later?
0050       if (depth_ != AllDepths && depth_ != Depth1 && depth_ != Depth2) {
0051         throw cms::Exception("Configuration Error")
0052             << "hcalDepth passed to EgammaTowerExtractor is invalid " << std::endl;
0053       }
0054     }
0055
0056     ~EgammaTowerExtractor() override;
0057
0058     void fillVetos(const edm::Event& ev, const edm::EventSetup& evSetup, const reco::TrackCollection& tracks) override {
0059     }
0060     reco::IsoDeposit deposit(const edm::Event& ev,
0061                              const edm::EventSetup& evSetup,
0062                              const reco::Track& track) const override {
0063       throw cms::Exception("Configuration Error")
0064           << "This extractor " << (typeid(this).name()) << " is not made for tracks";
0065     }
0066     reco::IsoDeposit deposit(const edm::Event& ev,
0067                              const edm::EventSetup& evSetup,
0068                              const reco::Candidate& c) const override;
0069
0070   private:
0071     double extRadius2_;
0072     double intRadius_;
0073     double etLow_;
0074
0075     edm::EDGetTokenT<CaloTowerCollection> caloTowerToken;
0076     int depth_;
0077     //const CaloTowerCollection *towercollection_ ;
0078   };
0079 }  // namespace egammaisolation
0080
0081 #include "FWCore/Framework/interface/MakerMacros.h"
0082 #include "PhysicsTools/IsolationAlgos/interface/IsoDepositExtractorFactory.h"
0083 DEFINE_EDM_PLUGIN(IsoDepositExtractorFactory, egammaisolation::EgammaTowerExtractor, "EgammaTowerExtractor");
0084
0085 using namespace ROOT::Math::VectorUtil;
0086
0087 using namespace egammaisolation;
0088 using namespace reco::isodeposit;
0089
0090 EgammaTowerExtractor::~EgammaTowerExtractor() {}
0091
0092 reco::IsoDeposit EgammaTowerExtractor::deposit(const edm::Event& iEvent,
0093                                                const edm::EventSetup& iSetup,
0094                                                const reco::Candidate& emObject) const {
0095   //Take the SC position
0096   reco::SuperClusterRef sc = emObject.get<reco::SuperClusterRef>();
0097   math::XYZPoint caloPosition = sc->position();
0098
0099   Direction candDir(caloPosition.eta(), caloPosition.phi());
0100   reco::IsoDeposit deposit(candDir);
0101   deposit.setVeto(reco::IsoDeposit::Veto(candDir, intRadius_));
0102   deposit.addCandEnergy(sc->energy() * sin(2 * atan(exp(-sc->eta()))));
0103
0104   //loop over tracks
0105   for (auto const& trItr : iEvent.get(caloTowerToken)) {
0106     double depEt = 0;
0107     //the hcal can be seperated into different depths
0108     //currently it is setup to check that the depth is valid in constructor
0109     //if the depth is not valid it fails gracefully
0110     //small bug fix, hadEnergyHeInnerLater returns zero for towers which are only depth 1
0111     //but we want Depth1 isolation to include these so we have to manually check for this
0112     if (depth_ == AllDepths)
0113       depEt = trItr.hadEt();
0114     else if (depth_ == Depth1)
0115       depEt = trItr.ietaAbs() < 18 || trItr.ietaAbs() > 29 ? trItr.hadEt()
0116                                                            : trItr.hadEnergyHeInnerLayer() * sin(trItr.p4().theta());
0117     else if (depth_ == Depth2)
0118       depEt = trItr.hadEnergyHeOuterLayer() * sin(trItr.p4().theta());
0119
0120     if (depEt < etLow_)
0121       continue;
0122
0123     Direction towerDir(trItr.eta(), trItr.phi());
0124     double dR2 = candDir.deltaR2(towerDir);
0125
0126     if (dR2 < extRadius2_) {
0127       deposit.addDeposit(towerDir, depEt);
0128     }
0129
0130   }  //end loop over tracks
0131
0132   return deposit;
0133 }