CMSSW Reference Manual

00001 //*****************************************************************************
00002 // File:      EgammaRecHitExtractor.cc
00003 // ----------------------------------------------------------------------------
00004 // OrigAuth:  Matthias Mozer, hacked by Sam Harper (ie the ugly stuff is mine)
00005 // Institute: IIHE-VUB, RAL
00006 //=============================================================================
00007 //*****************************************************************************
00008 //C++ includes
00009 #include <vector>
00010 #include <functional>
00011 
00012 //ROOT includes
00013 #include <Math/VectorUtil.h>
00014 
00015 //CMSSW includes
00016 #include "RecoEgamma/EgammaIsolationAlgos/plugins/EgammaRecHitExtractor.h"
00017 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaRecHitIsolation.h"
00018 #include "DataFormats/DetId/interface/DetId.h"
00019 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
00020 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
00021 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00022 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00023 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00024 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
00025 #include "RecoCaloTools/Selectors/interface/CaloConeSelector.h"
00026 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00027 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
00028 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
00029 #include "RecoCaloTools/MetaCollections/interface/CaloRecHitMetaCollections.h"
00030 
00031 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00032 
00033 using namespace std;
00034 using namespace egammaisolation;
00035 using namespace reco::isodeposit;
00036 
00037 EgammaRecHitExtractor::EgammaRecHitExtractor(const edm::ParameterSet& par) : 
00038     etMin_(par.getParameter<double>("etMin")),
00039     energyMin_(par.getParameter<double>("energyMin")),
00040     minCandEt_(par.getParameter<double>("minCandEt")),
00041     extRadius_(par.getParameter<double>("extRadius")),
00042     intRadius_(par.getParameter<double>("intRadius")),
00043     intStrip_(par.getParameter<double>("intStrip")),
00044     barrelRecHitsTag_(par.getParameter<edm::InputTag>("barrelRecHits")), 
00045     endcapRecHitsTag_(par.getParameter<edm::InputTag>("endcapRecHits")),
00046     barrelEcalHitsTag_(par.getParameter<edm::InputTag>("barrelEcalHits")), 
00047     endcapEcalHitsTag_(par.getParameter<edm::InputTag>("endcapEcalHits")),
00048     fakeNegativeDeposit_(par.getParameter<bool>("subtractSuperClusterEnergy")),
00049     tryBoth_(par.getParameter<bool>("tryBoth")),
00050     sameTag_(false)
00051 { 
00052     if ((intRadius_ != 0.0) && (fakeNegativeDeposit_)) {
00053         throw cms::Exception("Configuration Error") << "EgammaRecHitExtractor: " << 
00054             "If you use 'subtractSuperClusterEnergy', you *must* set 'intRadius' to ZERO; it does not make sense, otherwise.";
00055     }
00056     std::string isoVariable = par.getParameter<std::string>("isolationVariable");
00057     if (isoVariable == "et") {
00058         useEt_ = true;
00059     } else if (isoVariable == "energy") {
00060         useEt_ = false;
00061     } else {
00062         throw cms::Exception("Configuration Error") << "EgammaRecHitExtractor: isolationVariable '" << isoVariable << "' not known. " 
00063             << " Supported values are 'et', 'energy'. ";
00064     }
00065     std::string detector = par.getParameter<std::string>("detector");
00066     if (detector == "Ecal") {
00067         detector_ = DetId::Ecal;
00068     } else if (detector == "Hcal") {
00069         detector_ = DetId::Hcal;
00070     } else {
00071         throw cms::Exception("Configuration Error") << "EgammaRecHitExtractor: detector '" << detector << "' not known. " 
00072             << " Supported values are 'Ecal', 'Hcal'. ";
00073     }
00074 
00075     if (endcapRecHitsTag_.encode() ==  barrelRecHitsTag_.encode()) {
00076         sameTag_ = true;
00077         if (tryBoth_) {
00078             edm::LogWarning("EgammaRecHitExtractor") << "If you have configured 'barrelRecHits' == 'endcapRecHits', so I'm switching 'tryBoth' to FALSE.";
00079             tryBoth_ = false;
00080         }
00081     }
00082 
00083     paramForIsolBarrel_.push_back(par.getParameter<double>("checkIsoExtRBarrel")); 
00084     paramForIsolBarrel_.push_back(par.getParameter<double>("checkIsoInnRBarrel"));
00085     paramForIsolBarrel_.push_back(par.getParameter<double>("checkIsoEtaStripBarrel"));
00086     paramForIsolBarrel_.push_back(par.getParameter<double>("checkIsoEtRecHitBarrel"));     
00087     paramForIsolBarrel_.push_back(par.getParameter<double>("checkIsoEtCutBarrel"));     
00088 
00089     paramForIsolEndcap_.push_back(par.getParameter<double>("checkIsoExtREndcap")); 
00090     paramForIsolEndcap_.push_back(par.getParameter<double>("checkIsoInnREndcap"));
00091     paramForIsolEndcap_.push_back(par.getParameter<double>("checkIsoEtaStripEndcap"));
00092     paramForIsolEndcap_.push_back(par.getParameter<double>("checkIsoEtRecHitEndcap"));     
00093     paramForIsolEndcap_.push_back(par.getParameter<double>("checkIsoEtCutEndcap"));     
00094 
00095 }
00096 
00097 EgammaRecHitExtractor::~EgammaRecHitExtractor() { }
00098 
00099 reco::IsoDeposit EgammaRecHitExtractor::deposit(const edm::Event & iEvent, 
00100         const edm::EventSetup & iSetup, const reco::Candidate &emObject ) const {
00101     edm::ESHandle<CaloGeometry> pG;
00102     iSetup.get<CaloGeometryRecord>().get(pG);
00103     const CaloGeometry* caloGeom = pG.product(); 
00104 
00105     static std::string metname = "EgammaIsolationAlgos|EgammaRecHitExtractor";
00106 
00107     std::auto_ptr<const CaloRecHitMetaCollectionV> barrelRecHits(0), endcapRecHits(0);
00108 
00109     //Get barrel ECAL RecHits for Iso checking
00110     edm::Handle<EcalRecHitCollection> barrelEcalRecHitsH;
00111     iEvent.getByLabel(barrelEcalHitsTag_, barrelEcalRecHitsH);
00112 
00113     //Get endcap ECAL RecHits for Iso checking
00114     edm::Handle<EcalRecHitCollection> endcapEcalRecHitsH;
00115     iEvent.getByLabel(endcapEcalHitsTag_, endcapEcalRecHitsH);
00116 
00117     if (detector_ == DetId::Ecal) {
00118         barrelRecHits = std::auto_ptr<const CaloRecHitMetaCollectionV>(new EcalRecHitMetaCollection(*barrelEcalRecHitsH));
00119         if (!sameTag_) {
00120             endcapRecHits = std::auto_ptr<const CaloRecHitMetaCollectionV>(new EcalRecHitMetaCollection(*endcapEcalRecHitsH));
00121         }
00122     } else if (detector_ == DetId::Hcal) {
00123         edm::Handle<HBHERecHitCollection> barrelHcalRecHitsH, endcapHcalRecHitsH;
00124         iEvent.getByLabel(barrelRecHitsTag_, barrelHcalRecHitsH);
00125         barrelRecHits = std::auto_ptr<const CaloRecHitMetaCollectionV>(new HBHERecHitMetaCollection(*barrelHcalRecHitsH));
00126         if (!sameTag_) {
00127             iEvent.getByLabel(endcapRecHitsTag_, endcapHcalRecHitsH);
00128             endcapRecHits = std::auto_ptr<const CaloRecHitMetaCollectionV>(new HBHERecHitMetaCollection(*endcapHcalRecHitsH));
00129         }
00130     }
00131 
00132     reco::SuperClusterRef sc = emObject.get<reco::SuperClusterRef>();
00133     math::XYZPoint caloPosition = sc->position();
00134     GlobalPoint point(caloPosition.x(), caloPosition.y() , caloPosition.z());
00135 
00136     Direction candDir(caloPosition.eta(), caloPosition.phi());
00137     reco::IsoDeposit deposit( candDir );
00138     deposit.setVeto( reco::IsoDeposit::Veto(candDir, intRadius_) ); 
00139     double sinTheta = sin(2*atan(exp(-sc->eta())));
00140     deposit.addCandEnergy(sc->energy() * (useEt_ ? sinTheta : 1.0)) ;
00141 
00142     //avoid slow preshower
00143     int bid,eid;
00144     if(detector_==DetId::Ecal){
00145         bid=EcalBarrel;
00146         eid=EcalEndcap;
00147     }else{
00148         bid=0;eid=0;
00149     }
00150 
00151     double fakeEnergy = -sc->rawEnergy();
00152     if (fakeNegativeDeposit_) {
00153         deposit.addDeposit(candDir, fakeEnergy * (useEt_ ?  sinTheta : 1.0)); // not exactly clean...
00154     }
00155 
00156     std::auto_ptr<CaloRecHitMetaCollectionV> ecalRecHits(0); 
00157     double extRadius, innRadius, etaStrip, minEtRecHit, isolEtCut;
00158     if( abs(sc->eta()) < 1.5 ) {
00159         extRadius   = paramForIsolBarrel_[0];
00160         innRadius   = paramForIsolBarrel_[1];
00161         etaStrip    = paramForIsolBarrel_[2];
00162         minEtRecHit = paramForIsolBarrel_[3];
00163         isolEtCut   = paramForIsolBarrel_[4];
00164         ecalRecHits = std::auto_ptr<CaloRecHitMetaCollectionV>(new EcalRecHitMetaCollection(*barrelEcalRecHitsH));
00165     } else {
00166         extRadius   = paramForIsolEndcap_[0];
00167         innRadius   = paramForIsolEndcap_[1];
00168         etaStrip    = paramForIsolEndcap_[2];
00169         minEtRecHit = paramForIsolEndcap_[3];
00170         isolEtCut   = paramForIsolEndcap_[4];
00171         ecalRecHits = std::auto_ptr<CaloRecHitMetaCollectionV>(new EcalRecHitMetaCollection(*endcapEcalRecHitsH));
00172     }
00173 
00174     EgammaRecHitIsolation candIso(extRadius,innRadius,etaStrip,minEtRecHit,pG,&(*ecalRecHits),DetId::Ecal);
00175     if ( sc->energy()*sinTheta < minCandEt_ || candIso.getEtSum(&emObject) > isolEtCut ) {
00176         deposit.addDeposit( Direction(caloPosition.eta(), caloPosition.phi()+0.15), 10000 );
00177         deposit.addDeposit( Direction(caloPosition.eta(), caloPosition.phi()+0.25), 100000 );
00178     } else {
00179         CaloDualConeSelector doubleConeSelBarrel(intRadius_ ,extRadius_, caloGeom, detector_,bid);
00180         CaloDualConeSelector doubleConeSelEndcap(intRadius_ ,extRadius_, caloGeom, detector_,eid);
00181 
00182         // 3 possible options
00183         bool inBarrel = sameTag_ || ( abs(sc->eta()) < 1.5 );
00184         if (inBarrel || tryBoth_) {
00185             collect(deposit, point, doubleConeSelBarrel, caloGeom, *barrelRecHits);
00186         } 
00187         if ((!inBarrel) || tryBoth_) {
00188             collect(deposit, point, doubleConeSelEndcap, caloGeom, *endcapRecHits);
00189         }
00190     }
00191     return deposit;
00192 }
00193 
00194 void EgammaRecHitExtractor::collect(reco::IsoDeposit &deposit, 
00195         const GlobalPoint &caloPosition, CaloDualConeSelector &cone, 
00196         const CaloGeometry* caloGeom,
00197         const CaloRecHitMetaCollectionV &hits) const 
00198 {
00199     std::auto_ptr<CaloRecHitMetaCollectionV> chosen = cone.select(caloPosition, hits);
00200     for (CaloRecHitMetaCollectionV::const_iterator i = chosen->begin(), end = chosen->end() ; i != end;  ++i)  {
00201         const  GlobalPoint & position = caloGeom->getPosition(i->detid());
00202         double eta = position.eta();
00203         double energy = i->energy();
00204         double et = energy*sin(2*atan(exp(-eta)));
00205         if ( et > etMin_ && energy > energyMin_ && fabs(eta-caloPosition.eta()) > intStrip_ ){
00206             deposit.addDeposit( Direction(eta, position.phi()), (useEt_ ? et : energy) );
00207         }
00208     }
00209 } 
00210 
00211