Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:24:54 

0001 //*****************************************************************************
0002 // File:      EgammaRecHitExtractor.cc
0003 // ----------------------------------------------------------------------------
0004 // OrigAuth:  Matthias Mozer, adapted from EgammaHcalExtractor by S. Harper
0005 // Institute: IIHE-VUB, RAL
0006 //=============================================================================
0007 //*****************************************************************************
0008 
0009 #include "CommonTools/Utils/interface/StringToEnumValue.h"
0010 #include "DataFormats/Candidate/interface/Candidate.h"
0011 #include "DataFormats/DetId/interface/DetId.h"
0012 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0013 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0014 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
0015 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0016 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0017 #include "DataFormats/Math/interface/deltaPhi.h"
0018 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
0019 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0020 #include "DataFormats/TrackReco/interface/Track.h"
0021 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0022 #include "FWCore/Framework/interface/ConsumesCollector.h"
0023 #include "FWCore/Framework/interface/Event.h"
0024 #include "FWCore/Framework/interface/EventSetup.h"
0025 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0026 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0027 #include "FWCore/Utilities/interface/InputTag.h"
0028 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0029 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0030 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0031 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0032 #include "PhysicsTools/IsolationAlgos/interface/IsoDepositExtractor.h"
0033 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0034 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0035 
0036 #include <Math/VectorUtil.h>
0037 
0038 #include <vector>
0039 #include <functional>
0040 
0041 namespace egammaisolation {
0042 
0043   class EgammaRecHitExtractor : public reco::isodeposit::IsoDepositExtractor {
0044   public:
0045     EgammaRecHitExtractor(const edm::ParameterSet& par, edm::ConsumesCollector&& iC) : EgammaRecHitExtractor(par, iC) {}
0046     EgammaRecHitExtractor(const edm::ParameterSet& par, edm::ConsumesCollector& iC);
0047     ~EgammaRecHitExtractor() override;
0048     void fillVetos(const edm::Event& ev, const edm::EventSetup& evSetup, const reco::TrackCollection& tracks) override {
0049     }
0050     reco::IsoDeposit deposit(const edm::Event& ev,
0051                              const edm::EventSetup& evSetup,
0052                              const reco::Track& track) const override {
0053       throw cms::Exception("Configuration Error")
0054           << "This extractor " << (typeid(this).name()) << " is not made for tracks";
0055     }
0056     reco::IsoDeposit deposit(const edm::Event& ev,
0057                              const edm::EventSetup& evSetup,
0058                              const reco::Candidate& c) const override;
0059 
0060   private:
0061     void collect(reco::IsoDeposit& deposit,
0062                  const reco::SuperClusterRef& sc,
0063                  const CaloSubdetectorGeometry* subdet,
0064                  const CaloGeometry* caloGeom,
0065                  const EcalRecHitCollection& hits,
0066                  //const EcalChannelStatus* chStatus,
0067                  const EcalSeverityLevelAlgo* sevLevel,
0068                  bool barrel) const;
0069 
0070     double etMin_;
0071     double energyMin_;
0072     double extRadius_;
0073     double intRadius_;
0074     double intStrip_;
0075     edm::InputTag barrelEcalHitsTag_;
0076     edm::InputTag endcapEcalHitsTag_;
0077     edm::EDGetTokenT<EcalRecHitCollection> barrelEcalHitsToken_;
0078     edm::EDGetTokenT<EcalRecHitCollection> endcapEcalHitsToken_;
0079     edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geometryToken_;
0080     edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> sevlvToken_;
0081     bool fakeNegativeDeposit_;
0082     bool tryBoth_;
0083     bool useEt_;
0084     bool vetoClustered_;
0085     bool sameTag_;
0086     //int   severityLevelCut_;
0087     //float severityRecHitThreshold_;
0088     //std::string spIdString_;
0089     //float spIdThreshold_;
0090     //EcalSeverityLevelAlgo::SpikeId spId_;
0091     //std::vector<int> v_chstatus_;
0092     std::vector<int> severitiesexclEB_;
0093     std::vector<int> severitiesexclEE_;
0094     std::vector<int> flagsexclEB_;
0095     std::vector<int> flagsexclEE_;
0096   };
0097 }  // namespace egammaisolation
0098 
0099 #include "FWCore/Framework/interface/MakerMacros.h"
0100 #include "PhysicsTools/IsolationAlgos/interface/IsoDepositExtractorFactory.h"
0101 DEFINE_EDM_PLUGIN(IsoDepositExtractorFactory, egammaisolation::EgammaRecHitExtractor, "EgammaRecHitExtractor");
0102 
0103 using namespace std;
0104 using namespace egammaisolation;
0105 using namespace reco::isodeposit;
0106 
0107 EgammaRecHitExtractor::EgammaRecHitExtractor(const edm::ParameterSet& par, edm::ConsumesCollector& iC)
0108     : etMin_(par.getParameter<double>("etMin")),
0109       energyMin_(par.getParameter<double>("energyMin")),
0110       extRadius_(par.getParameter<double>("extRadius")),
0111       intRadius_(par.getParameter<double>("intRadius")),
0112       intStrip_(par.getParameter<double>("intStrip")),
0113       barrelEcalHitsTag_(par.getParameter<edm::InputTag>("barrelEcalHits")),
0114       endcapEcalHitsTag_(par.getParameter<edm::InputTag>("endcapEcalHits")),
0115       barrelEcalHitsToken_(iC.consumes<EcalRecHitCollection>(barrelEcalHitsTag_)),
0116       endcapEcalHitsToken_(iC.consumes<EcalRecHitCollection>(endcapEcalHitsTag_)),
0117       geometryToken_(iC.esConsumes()),
0118       sevlvToken_(iC.esConsumes()),
0119       fakeNegativeDeposit_(par.getParameter<bool>("subtractSuperClusterEnergy")),
0120       tryBoth_(par.getParameter<bool>("tryBoth")),
0121       vetoClustered_(par.getParameter<bool>("vetoClustered")),
0122       sameTag_(false)
0123 //severityLevelCut_(par.getParameter<int>("severityLevelCut"))
0124 //severityRecHitThreshold_(par.getParameter<double>("severityRecHitThreshold")),
0125 //spIdString_(par.getParameter<std::string>("spikeIdString")),
0126 //spIdThreshold_(par.getParameter<double>("spikeIdThreshold")),
0127 {
0128   const std::vector<std::string> flagnamesEB = par.getParameter<std::vector<std::string> >("RecHitFlagToBeExcludedEB");
0129 
0130   const std::vector<std::string> flagnamesEE = par.getParameter<std::vector<std::string> >("RecHitFlagToBeExcludedEE");
0131 
0132   flagsexclEB_ = StringToEnumValue<EcalRecHit::Flags>(flagnamesEB);
0133 
0134   flagsexclEE_ = StringToEnumValue<EcalRecHit::Flags>(flagnamesEE);
0135 
0136   const std::vector<std::string> severitynamesEB =
0137       par.getParameter<std::vector<std::string> >("RecHitSeverityToBeExcludedEB");
0138 
0139   severitiesexclEB_ = StringToEnumValue<EcalSeverityLevel::SeverityLevel>(severitynamesEB);
0140 
0141   const std::vector<std::string> severitynamesEE =
0142       par.getParameter<std::vector<std::string> >("RecHitSeverityToBeExcludedEE");
0143 
0144   severitiesexclEE_ = StringToEnumValue<EcalSeverityLevel::SeverityLevel>(severitynamesEE);
0145 
0146   if ((intRadius_ != 0.0) && (fakeNegativeDeposit_)) {
0147     throw cms::Exception("Configuration Error") << "EgammaRecHitExtractor: "
0148                                                 << "If you use 'subtractSuperClusterEnergy', you *must* set "
0149                                                    "'intRadius' to ZERO; it does not make sense, otherwise.";
0150   }
0151   std::string isoVariable = par.getParameter<std::string>("isolationVariable");
0152   if (isoVariable == "et") {
0153     useEt_ = true;
0154   } else if (isoVariable == "energy") {
0155     useEt_ = false;
0156   } else {
0157     throw cms::Exception("Configuration Error")
0158         << "EgammaRecHitExtractor: isolationVariable '" << isoVariable << "' not known. "
0159         << " Supported values are 'et', 'energy'. ";
0160   }
0161   if (endcapEcalHitsTag_.encode() == barrelEcalHitsTag_.encode()) {
0162     sameTag_ = true;
0163     if (tryBoth_) {
0164       edm::LogWarning("EgammaRecHitExtractor")
0165           << "If you have configured 'barrelRecHits' == 'endcapRecHits', so I'm switching 'tryBoth' to FALSE.";
0166       tryBoth_ = false;
0167     }
0168   }
0169 }
0170 
0171 EgammaRecHitExtractor::~EgammaRecHitExtractor() {}
0172 
0173 reco::IsoDeposit EgammaRecHitExtractor::deposit(const edm::Event& iEvent,
0174                                                 const edm::EventSetup& iSetup,
0175                                                 const reco::Candidate& emObject) const {
0176   //Get the channel status from the db
0177   //edm::ESHandle<EcalChannelStatus> chStatus;
0178   //iSetup.get<EcalChannelStatusRcd>().get(chStatus);
0179 
0180   const EcalSeverityLevelAlgo* sevLevel = &iSetup.getData(sevlvToken_);
0181 
0182   const CaloGeometry* caloGeom = &iSetup.getData(geometryToken_);
0183   const CaloSubdetectorGeometry* barrelgeom = caloGeom->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0184   const CaloSubdetectorGeometry* endcapgeom = caloGeom->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0185 
0186   static const std::string metname = "EgammaIsolationAlgos|EgammaRecHitExtractor";
0187 
0188   //define isodeposit starting from candidate
0189   reco::SuperClusterRef sc = emObject.get<reco::SuperClusterRef>();
0190   math::XYZPoint caloPosition = sc->position();
0191 
0192   Direction candDir(caloPosition.eta(), caloPosition.phi());
0193   reco::IsoDeposit deposit(candDir);
0194   deposit.setVeto(reco::IsoDeposit::Veto(candDir, intRadius_));
0195   double sinTheta = sin(2 * atan(exp(-sc->eta())));
0196   deposit.addCandEnergy(sc->energy() * (useEt_ ? sinTheta : 1.0));
0197 
0198   // subtract supercluster if desired
0199   double fakeEnergy = -sc->rawEnergy();
0200   if (fakeNegativeDeposit_) {
0201     deposit.addDeposit(candDir, fakeEnergy * (useEt_ ? sinTheta : 1.0));  // not exactly clean...
0202   }
0203 
0204   // fill rechits
0205   bool inBarrel = sameTag_ || (abs(sc->eta()) < 1.479);  //check for barrel. If only one collection is used, use barrel
0206   if (inBarrel || tryBoth_) {
0207     collect(deposit, sc, barrelgeom, caloGeom, iEvent.get(barrelEcalHitsToken_), sevLevel, true);
0208   }
0209 
0210   if ((!inBarrel) || tryBoth_) {
0211     collect(deposit, sc, endcapgeom, caloGeom, iEvent.get(endcapEcalHitsToken_), sevLevel, false);
0212   }
0213 
0214   return deposit;
0215 }
0216 
0217 void EgammaRecHitExtractor::collect(reco::IsoDeposit& deposit,
0218                                     const reco::SuperClusterRef& sc,
0219                                     const CaloSubdetectorGeometry* subdet,
0220                                     const CaloGeometry* caloGeom,
0221                                     const EcalRecHitCollection& hits,
0222                                     //const EcalChannelStatus* chStatus,
0223                                     const EcalSeverityLevelAlgo* sevLevel,
0224                                     bool barrel) const {
0225   GlobalPoint caloPosition(sc->position().x(), sc->position().y(), sc->position().z());
0226   CaloSubdetectorGeometry::DetIdSet chosen = subdet->getCells(caloPosition, extRadius_);
0227   EcalRecHitCollection::const_iterator j = hits.end();
0228   double caloeta = caloPosition.eta();
0229   double calophi = caloPosition.phi();
0230   double r2 = intRadius_ * intRadius_;
0231 
0232   std::vector<std::pair<DetId, float> >::const_iterator rhIt;
0233 
0234   for (CaloSubdetectorGeometry::DetIdSet::const_iterator i = chosen.begin(), end = chosen.end(); i != end; ++i) {
0235     j = hits.find(*i);
0236     if (j != hits.end()) {
0237       const GlobalPoint& position = caloGeom->getPosition(*i);
0238       double eta = position.eta();
0239       double phi = position.phi();
0240       double energy = j->energy();
0241       double et = energy * position.perp() / position.mag();
0242       double phiDiff = reco::deltaPhi(phi, calophi);
0243 
0244       //check if we are supposed to veto clustered and then do so
0245       if (vetoClustered_) {
0246         //Loop over basic clusters:
0247         bool isClustered = false;
0248         for (auto bcIt = sc->clustersBegin(); bcIt != sc->clustersEnd(); ++bcIt) {
0249           for (rhIt = (*bcIt)->hitsAndFractions().begin(); rhIt != (*bcIt)->hitsAndFractions().end(); ++rhIt) {
0250             if (rhIt->first == *i)
0251               isClustered = true;
0252             if (isClustered)
0253               break;
0254           }
0255           if (isClustered)
0256             break;
0257         }  //end loop over basic clusters
0258 
0259         if (isClustered)
0260           continue;
0261       }  //end if removeClustered
0262 
0263       std::vector<int>::const_iterator sit;
0264       int severityFlag = sevLevel->severityLevel(j->detid(), hits);
0265       if (barrel) {
0266         sit = std::find(severitiesexclEB_.begin(), severitiesexclEB_.end(), severityFlag);
0267         if (sit != severitiesexclEB_.end())
0268           continue;
0269       } else {
0270         sit = std::find(severitiesexclEE_.begin(), severitiesexclEE_.end(), severityFlag);
0271         if (sit != severitiesexclEE_.end())
0272           continue;
0273       }
0274 
0275       if (barrel) {
0276         // new rechit flag checks
0277         if (!j->checkFlag(EcalRecHit::kGood)) {
0278           if (j->checkFlags(flagsexclEB_)) {
0279             continue;
0280           }
0281         }
0282       } else {
0283         // new rechit flag checks
0284         if (!j->checkFlag(EcalRecHit::kGood)) {
0285           if (j->checkFlags(flagsexclEE_)) {
0286             continue;
0287           }
0288         }
0289       }
0290 
0291       if (et > etMin_ && energy > energyMin_  //Changed to fabs - then changed back to energy
0292           && fabs(eta - caloeta) > intStrip_ && (eta - caloeta) * (eta - caloeta) + phiDiff * phiDiff > r2) {
0293         deposit.addDeposit(Direction(eta, phi), (useEt_ ? et : energy));
0294       }
0295     }
0296   }
0297 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!