00001
00002
00003
00004
00005
00006
00007
00008
00009 #include <vector>
00010 #include <functional>
00011
00012
00013 #include <Math/VectorUtil.h>
00014
00015
00016 #include "CommonTools/Utils/interface/StringToEnumValue.h"
00017
00018 #include "RecoEgamma/EgammaIsolationAlgos/plugins/EgammaRecHitExtractor.h"
00019 #include "DataFormats/DetId/interface/DetId.h"
00020 #include "DataFormats/EcalDetId/interface/EcalSubdetector.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 "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00026 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
00027 #include "RecoCaloTools/MetaCollections/interface/CaloRecHitMetaCollections.h"
00028 #include "DataFormats/Math/interface/deltaPhi.h"
00029 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
00030 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
00031 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
00032 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
00033
00034 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
00035
00036 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00037
00038 using namespace std;
00039 using namespace egammaisolation;
00040 using namespace reco::isodeposit;
00041
00042 EgammaRecHitExtractor::EgammaRecHitExtractor(const edm::ParameterSet& par) :
00043 etMin_(par.getParameter<double>("etMin")),
00044 energyMin_(par.getParameter<double>("energyMin")),
00045 extRadius_(par.getParameter<double>("extRadius")),
00046 intRadius_(par.getParameter<double>("intRadius")),
00047 intStrip_(par.getParameter<double>("intStrip")),
00048 barrelEcalHitsTag_(par.getParameter<edm::InputTag>("barrelEcalHits")),
00049 endcapEcalHitsTag_(par.getParameter<edm::InputTag>("endcapEcalHits")),
00050 fakeNegativeDeposit_(par.getParameter<bool>("subtractSuperClusterEnergy")),
00051 tryBoth_(par.getParameter<bool>("tryBoth")),
00052 vetoClustered_(par.getParameter<bool>("vetoClustered")),
00053 sameTag_(false),
00054 severityLevelCut_(par.getParameter<int>("severityLevelCut"))
00055
00056
00057
00058
00059 {
00060 const std::vector<std::string> flagnames =
00061 par.getParameter<std::vector<std::string> >("recHitFlagsToBeExcluded");
00062
00063 v_chstatus_= StringToEnumValue<EcalRecHit::Flags>(flagnames);
00064
00065
00066
00067
00068
00069
00070 if ((intRadius_ != 0.0) && (fakeNegativeDeposit_)) {
00071 throw cms::Exception("Configuration Error") << "EgammaRecHitExtractor: " <<
00072 "If you use 'subtractSuperClusterEnergy', you *must* set 'intRadius' to ZERO; it does not make sense, otherwise.";
00073 }
00074 std::string isoVariable = par.getParameter<std::string>("isolationVariable");
00075 if (isoVariable == "et") {
00076 useEt_ = true;
00077 } else if (isoVariable == "energy") {
00078 useEt_ = false;
00079 } else {
00080 throw cms::Exception("Configuration Error") << "EgammaRecHitExtractor: isolationVariable '" << isoVariable << "' not known. "
00081 << " Supported values are 'et', 'energy'. ";
00082 }
00083 if (endcapEcalHitsTag_.encode() == barrelEcalHitsTag_.encode()) {
00084 sameTag_ = true;
00085 if (tryBoth_) {
00086 edm::LogWarning("EgammaRecHitExtractor") << "If you have configured 'barrelRecHits' == 'endcapRecHits', so I'm switching 'tryBoth' to FALSE.";
00087 tryBoth_ = false;
00088 }
00089 }
00090
00091 }
00092
00093 EgammaRecHitExtractor::~EgammaRecHitExtractor() { }
00094
00095 reco::IsoDeposit EgammaRecHitExtractor::deposit(const edm::Event & iEvent,
00096 const edm::EventSetup & iSetup, const reco::Candidate &emObject ) const {
00097 edm::ESHandle<CaloGeometry> pG;
00098 iSetup.get<CaloGeometryRecord>().get(pG);
00099
00100
00101 edm::ESHandle<EcalChannelStatus> chStatus;
00102 iSetup.get<EcalChannelStatusRcd>().get(chStatus);
00103
00104 edm::ESHandle<EcalSeverityLevelAlgo> sevlv;
00105 iSetup.get<EcalSeverityLevelAlgoRcd>().get(sevlv);
00106
00107 const CaloGeometry* caloGeom = pG.product();
00108 const CaloSubdetectorGeometry* barrelgeom = caloGeom->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
00109 const CaloSubdetectorGeometry* endcapgeom = caloGeom->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);
00110
00111 static std::string metname = "EgammaIsolationAlgos|EgammaRecHitExtractor";
00112
00113 std::auto_ptr<const CaloRecHitMetaCollectionV> barrelRecHits(0), endcapRecHits(0);
00114
00115
00116 edm::Handle<EcalRecHitCollection> barrelEcalRecHitsH;
00117 iEvent.getByLabel(barrelEcalHitsTag_, barrelEcalRecHitsH);
00118
00119
00120 edm::Handle<EcalRecHitCollection> endcapEcalRecHitsH;
00121 iEvent.getByLabel(endcapEcalHitsTag_, endcapEcalRecHitsH);
00122
00123
00124 reco::SuperClusterRef sc = emObject.get<reco::SuperClusterRef>();
00125 math::XYZPoint caloPosition = sc->position();
00126
00127 Direction candDir(caloPosition.eta(), caloPosition.phi());
00128 reco::IsoDeposit deposit( candDir );
00129 deposit.setVeto( reco::IsoDeposit::Veto(candDir, intRadius_) );
00130 double sinTheta = sin(2*atan(exp(-sc->eta())));
00131 deposit.addCandEnergy(sc->energy() * (useEt_ ? sinTheta : 1.0)) ;
00132
00133
00134 double fakeEnergy = -sc->rawEnergy();
00135 if (fakeNegativeDeposit_) {
00136 deposit.addDeposit(candDir, fakeEnergy * (useEt_ ? sinTheta : 1.0));
00137 }
00138
00139
00140 bool inBarrel = sameTag_ || ( abs(sc->eta()) < 1.479 );
00141 if (inBarrel || tryBoth_) {
00142 collect(deposit, sc, barrelgeom, caloGeom, *barrelEcalRecHitsH, chStatus.product(), sevlv.product(), true);
00143 }
00144 if ((!inBarrel) || tryBoth_) {
00145 collect(deposit, sc, endcapgeom, caloGeom, *endcapEcalRecHitsH, chStatus.product(), sevlv.product(), false);
00146 }
00147
00148 return deposit;
00149 }
00150
00151 void EgammaRecHitExtractor::collect(reco::IsoDeposit &deposit,
00152 const reco::SuperClusterRef& sc, const CaloSubdetectorGeometry* subdet,
00153 const CaloGeometry* caloGeom,
00154 const EcalRecHitCollection &hits,
00155 const EcalChannelStatus* chStatus,
00156 const EcalSeverityLevelAlgo* sevLevel,
00157 bool barrel) const
00158 {
00159
00160 GlobalPoint caloPosition(sc->position().x(), sc->position().y() , sc->position().z());
00161 CaloSubdetectorGeometry::DetIdSet chosen = subdet->getCells(caloPosition,extRadius_);
00162 EcalRecHitCollection::const_iterator j=hits.end();
00163 double caloeta=caloPosition.eta();
00164 double calophi=caloPosition.phi();
00165 double r2 = intRadius_*intRadius_;
00166
00167 std::vector< std::pair<DetId, float> >::const_iterator rhIt;
00168
00169
00170 for (CaloSubdetectorGeometry::DetIdSet::const_iterator i = chosen.begin(), end = chosen.end() ; i != end; ++i) {
00171 j=hits.find(*i);
00172 if(j != hits.end()){
00173 const GlobalPoint & position = caloGeom->getPosition(*i);
00174 double eta = position.eta();
00175 double phi = position.phi();
00176 double energy = j->energy();
00177 double et = energy*position.perp()/position.mag();
00178 double phiDiff= reco::deltaPhi(phi,calophi);
00179
00180
00181 if(vetoClustered_) {
00182
00183
00184 bool isClustered = false;
00185 for( reco::CaloCluster_iterator bcIt = sc->clustersBegin();bcIt != sc->clustersEnd(); ++bcIt) {
00186 for(rhIt = (*bcIt)->hitsAndFractions().begin();rhIt != (*bcIt)->hitsAndFractions().end(); ++rhIt) {
00187 if( rhIt->first == *i ) isClustered = true;
00188 if( isClustered ) break;
00189 }
00190 if( isClustered ) break;
00191 }
00192
00193 if(isClustered) continue;
00194 }
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205 if(barrel && sevLevel->severityLevel(EBDetId(j->id()), hits) >= severityLevelCut_)
00206 continue;
00207
00208
00209
00210
00211
00212
00213
00214
00215 std::vector<int>::const_iterator vit = std::find( v_chstatus_.begin(), v_chstatus_.end(), ((EcalRecHit*)(&*j))->recoFlag() );
00216 if ( vit != v_chstatus_.end() ) continue;
00217
00218
00219 if(et > etMin_
00220 && energy > energyMin_
00221 && fabs(eta-caloeta) > intStrip_
00222 && (eta-caloeta)*(eta-caloeta) + phiDiff*phiDiff >r2 ) {
00223
00224 deposit.addDeposit( Direction(eta, phi), (useEt_ ? et : energy) );
00225
00226 }
00227 }
00228 }
00229 }
00230
00231