CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
PhotonIsolationCalculator.cc
Go to the documentation of this file.
1 
7 #include "RecoEgamma/PhotonIdentification/interface/PhotonIsolationCalculator.h"
8 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
9 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
10 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
11 #include "DataFormats/TrackReco/interface/Track.h"
12 #include "DataFormats/TrackReco/interface/TrackFwd.h"
13 #include "DataFormats/DetId/interface/DetId.h"
14 #include "DataFormats/EcalDetId/interface/EEDetId.h"
15 #include "DataFormats/EcalDetId/interface/EBDetId.h"
16 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
17 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
18 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
19 
20 #include "RecoEgamma/EgammaIsolationAlgos/interface/PhotonTkIsolation.h"
21 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaEcalIsolation.h"
22 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaRecHitIsolation.h"
23 
24 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
25 #include "Geometry/Records/interface/CaloGeometryRecord.h"
26 
27 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
28 
29 void PhotonIsolationCalculator::setup(const edm::ParameterSet& conf,
30  std::vector<int> const& flagsEB,
31  std::vector<int> const& flagsEE,
32  std::vector<int> const& severitiesEB,
33  std::vector<int> const& severitiesEE,
34  edm::ConsumesCollector&& iC) {
35  trackInputTag_ = iC.consumes<reco::TrackCollection>(conf.getParameter<edm::InputTag>("trackProducer"));
36  beamSpotProducerTag_ = iC.consumes<reco::BeamSpot>(conf.getParameter<edm::InputTag>("beamSpotProducer"));
37  barrelecalCollection_ =
38  iC.consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("barrelEcalRecHitCollection"));
39  endcapecalCollection_ =
40  iC.consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("endcapEcalRecHitCollection"));
41 
42  auto hbhetag = conf.getParameter<edm::InputTag>("HBHERecHitCollection");
43  if (not hbhetag.label().empty())
44  hbheRecHitsTag_ = iC.consumes<HBHERecHitCollection>(hbhetag);
45 
46  caloGeometryToken_ = decltype(caloGeometryToken_){iC.esConsumes()};
47  hcalTopologyToken_ = decltype(hcalTopologyToken_){iC.esConsumes()};
48  hcalChannelQualityToken_ = decltype(hcalChannelQualityToken_){iC.esConsumes(edm::ESInputTag("", "withTopo"))};
49  hcalSevLvlComputerToken_ = decltype(hcalSevLvlComputerToken_){iC.esConsumes()};
50  towerMapToken_ = decltype(towerMapToken_){iC.esConsumes()};
51  ecalSevLvlToken_ = iC.esConsumes();
52 
53  // gsfRecoInputTag_ = conf.getParameter<edm::InputTag>("GsfRecoCollection");
54  modulePhiBoundary_ = conf.getParameter<double>("modulePhiBoundary");
55  moduleEtaBoundary_ = conf.getParameter<std::vector<double>>("moduleEtaBoundary");
56  //
57  vetoClusteredEcalHits_ = conf.getParameter<bool>("vetoClustered");
58  useNumCrystals_ = conf.getParameter<bool>("useNumCrystals");
59 
61  int i = 0;
62  trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("TrackConeOuterRadiusA_Barrel"));
63  trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("TrackConeInnerRadiusA_Barrel"));
64  trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("isolationtrackThresholdA_Barrel"));
65  trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("longImpactParameterA_Barrel"));
66  trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("transImpactParameterA_Barrel"));
67  trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("isolationtrackEtaSliceA_Barrel"));
68 
69  i = 0;
70  ecalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitInnerRadiusA_Barrel"));
71  ecalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitOuterRadiusA_Barrel"));
72  ecalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitEtaSliceA_Barrel"));
73  ecalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitThreshEA_Barrel"));
74  ecalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitThreshEtA_Barrel"));
75 
76  hcalIsoInnerRadAEB_ = conf.getParameter<std::array<double, 7>>("HcalRecHitInnerRadiusA_Barrel");
77  hcalIsoOuterRadAEB_ = conf.getParameter<std::array<double, 7>>("HcalRecHitOuterRadiusA_Barrel");
78 
79  i = 0;
80  trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("TrackConeOuterRadiusB_Barrel"));
81  trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("TrackConeInnerRadiusB_Barrel"));
82  trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("isolationtrackThresholdB_Barrel"));
83  trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("longImpactParameterB_Barrel"));
84  trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("transImpactParameterB_Barrel"));
85  trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("isolationtrackEtaSliceB_Barrel"));
86 
87  i = 0;
88  ecalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitInnerRadiusB_Barrel"));
89  ecalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitOuterRadiusB_Barrel"));
90  ecalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitEtaSliceB_Barrel"));
91  ecalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitThreshEB_Barrel"));
92  ecalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitThreshEtB_Barrel"));
93 
94  hcalIsoInnerRadBEB_ = conf.getParameter<std::array<double, 7>>("HcalRecHitInnerRadiusB_Barrel");
95  hcalIsoOuterRadBEB_ = conf.getParameter<std::array<double, 7>>("HcalRecHitOuterRadiusB_Barrel");
96 
98  i = 0;
99  trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("TrackConeOuterRadiusA_Endcap"));
100  trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("TrackConeInnerRadiusA_Endcap"));
101  trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("isolationtrackThresholdA_Endcap"));
102  trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("longImpactParameterA_Endcap"));
103  trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("transImpactParameterA_Endcap"));
104  trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("isolationtrackEtaSliceA_Endcap"));
105 
106  i = 0;
107  ecalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitInnerRadiusA_Endcap"));
108  ecalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitOuterRadiusA_Endcap"));
109  ecalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitEtaSliceA_Endcap"));
110  ecalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitThreshEA_Endcap"));
111  ecalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitThreshEtA_Endcap"));
112 
113  hcalIsoInnerRadAEE_ = conf.getParameter<std::array<double, 7>>("HcalRecHitInnerRadiusA_Endcap");
114  hcalIsoOuterRadAEE_ = conf.getParameter<std::array<double, 7>>("HcalRecHitOuterRadiusA_Endcap");
115 
116  i = 0;
117  trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("TrackConeOuterRadiusB_Endcap"));
118  trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("TrackConeInnerRadiusB_Endcap"));
119  trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("isolationtrackThresholdB_Endcap"));
120  trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("longImpactParameterB_Endcap"));
121  trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("transImpactParameterB_Endcap"));
122  trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("isolationtrackEtaSliceB_Endcap"));
123 
124  i = 0;
125  ecalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitInnerRadiusB_Endcap"));
126  ecalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitOuterRadiusB_Endcap"));
127  ecalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitEtaSliceB_Endcap"));
128  ecalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitThreshEB_Endcap"));
129  ecalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitThreshEtB_Endcap"));
130 
131  hcalIsoInnerRadBEE_ = conf.getParameter<std::array<double, 7>>("HcalRecHitInnerRadiusB_Endcap");
132  hcalIsoOuterRadBEE_ = conf.getParameter<std::array<double, 7>>("HcalRecHitOuterRadiusB_Endcap");
133 
134  //Pick up the variables for the spike removal
135  flagsEB_ = flagsEB;
136  flagsEE_ = flagsEE;
137  severityExclEB_ = severitiesEB;
138  severityExclEE_ = severitiesEE;
139 
140  hcalIsoEThresHB_ = conf.getParameter<EgammaHcalIsolation::arrayHB>("recHitEThresholdHB");
141  hcalIsoEThresHE_ = conf.getParameter<EgammaHcalIsolation::arrayHE>("recHitEThresholdHE");
142  maxHcalSeverity_ = conf.getParameter<int>("maxHcalRecHitSeverity");
143 }
144 
145 void PhotonIsolationCalculator::calculate(const reco::Photon* pho,
146  const edm::Event& e,
147  const edm::EventSetup& es,
148  reco::Photon::FiducialFlags& phofid,
149  reco::Photon::IsolationVariables& phoisolR1,
150  reco::Photon::IsolationVariables& phoisolR2) const {
151  //Get fiducial flags. This does not really belong here
152  bool isEBPho = false;
153  bool isEEPho = false;
154  bool isEBEtaGap = false;
155  bool isEBPhiGap = false;
156  bool isEERingGap = false;
157  bool isEEDeeGap = false;
158  bool isEBEEGap = false;
159  classify(pho, isEBPho, isEEPho, isEBEtaGap, isEBPhiGap, isEERingGap, isEEDeeGap, isEBEEGap);
160  phofid.isEB = isEBPho;
161  phofid.isEE = isEEPho;
162  phofid.isEBEtaGap = isEBEtaGap;
163  phofid.isEBPhiGap = isEBPhiGap;
164  phofid.isEERingGap = isEERingGap;
165  phofid.isEEDeeGap = isEEDeeGap;
166  phofid.isEBEEGap = isEBEEGap;
167 
168  auto const& caloGeometry = es.getData(caloGeometryToken_);
169  auto const& hcalTopology = &es.getData(hcalTopologyToken_);
170  auto const& hcalChannelQuality = &es.getData(hcalChannelQualityToken_);
171  auto const& hcalSevLvlComputer = &es.getData(hcalSevLvlComputerToken_);
172  auto const& towerMap = es.getData(towerMapToken_);
173 
174  // Calculate isolation variables. cone sizes and thresholds
175  // are set for Barrel and endcap separately
176 
177  reco::SuperClusterRef scRef = pho->superCluster();
178  const reco::BasicCluster& seedCluster = *(scRef->seed());
179  DetId seedXtalId = seedCluster.hitsAndFractions()[0].first;
180  int detector = seedXtalId.subdetId();
181 
182  //Isolation parameters variables
183  double photonEcalRecHitConeInnerRadiusA_;
184  double photonEcalRecHitConeOuterRadiusA_;
185  double photonEcalRecHitEtaSliceA_;
186  double photonEcalRecHitThreshEA_;
187  double photonEcalRecHitThreshEtA_;
188  double trackConeOuterRadiusA_;
189  double trackConeInnerRadiusA_;
190  double isolationtrackThresholdA_;
191  double isolationtrackEtaSliceA_;
192  double trackLipRadiusA_;
193  double trackD0RadiusA_;
194  double photonEcalRecHitConeInnerRadiusB_;
195  double photonEcalRecHitConeOuterRadiusB_;
196  double photonEcalRecHitEtaSliceB_;
197  double photonEcalRecHitThreshEB_;
198  double photonEcalRecHitThreshEtB_;
199  double trackConeOuterRadiusB_;
200  double trackConeInnerRadiusB_;
201  double isolationtrackThresholdB_;
202  double isolationtrackEtaSliceB_;
203  double trackLipRadiusB_;
204  double trackD0RadiusB_;
205 
206  if (detector == EcalBarrel) {
207  trackConeOuterRadiusA_ = trkIsoBarrelRadiusA_[0];
208  trackConeInnerRadiusA_ = trkIsoBarrelRadiusA_[1];
209  isolationtrackThresholdA_ = trkIsoBarrelRadiusA_[2];
210  trackLipRadiusA_ = trkIsoBarrelRadiusA_[3];
211  trackD0RadiusA_ = trkIsoBarrelRadiusA_[4];
212  isolationtrackEtaSliceA_ = trkIsoBarrelRadiusA_[5];
213 
214  photonEcalRecHitConeInnerRadiusA_ = ecalIsoBarrelRadiusA_[0];
215  photonEcalRecHitConeOuterRadiusA_ = ecalIsoBarrelRadiusA_[1];
216  photonEcalRecHitEtaSliceA_ = ecalIsoBarrelRadiusA_[2];
217  photonEcalRecHitThreshEA_ = ecalIsoBarrelRadiusA_[3];
218  photonEcalRecHitThreshEtA_ = ecalIsoBarrelRadiusA_[4];
219 
220  trackConeOuterRadiusB_ = trkIsoBarrelRadiusB_[0];
221  trackConeInnerRadiusB_ = trkIsoBarrelRadiusB_[1];
222  isolationtrackThresholdB_ = trkIsoBarrelRadiusB_[2];
223  trackLipRadiusB_ = trkIsoBarrelRadiusB_[3];
224  trackD0RadiusB_ = trkIsoBarrelRadiusB_[4];
225  isolationtrackEtaSliceB_ = trkIsoBarrelRadiusB_[5];
226 
227  photonEcalRecHitConeInnerRadiusB_ = ecalIsoBarrelRadiusB_[0];
228  photonEcalRecHitConeOuterRadiusB_ = ecalIsoBarrelRadiusB_[1];
229  photonEcalRecHitEtaSliceB_ = ecalIsoBarrelRadiusB_[2];
230  photonEcalRecHitThreshEB_ = ecalIsoBarrelRadiusB_[3];
231  photonEcalRecHitThreshEtB_ = ecalIsoBarrelRadiusB_[4];
232  } else {
233  // detector==EcalEndcap
234 
235  trackConeOuterRadiusA_ = trkIsoEndcapRadiusA_[0];
236  trackConeInnerRadiusA_ = trkIsoEndcapRadiusA_[1];
237  isolationtrackThresholdA_ = trkIsoEndcapRadiusA_[2];
238  trackLipRadiusA_ = trkIsoEndcapRadiusA_[3];
239  trackD0RadiusA_ = trkIsoEndcapRadiusA_[4];
240  isolationtrackEtaSliceA_ = trkIsoEndcapRadiusA_[5];
241 
242  photonEcalRecHitConeInnerRadiusA_ = ecalIsoEndcapRadiusA_[0];
243  photonEcalRecHitConeOuterRadiusA_ = ecalIsoEndcapRadiusA_[1];
244  photonEcalRecHitEtaSliceA_ = ecalIsoEndcapRadiusA_[2];
245  photonEcalRecHitThreshEA_ = ecalIsoEndcapRadiusA_[3];
246  photonEcalRecHitThreshEtA_ = ecalIsoEndcapRadiusA_[4];
247 
248  trackConeOuterRadiusB_ = trkIsoEndcapRadiusB_[0];
249  trackConeInnerRadiusB_ = trkIsoEndcapRadiusB_[1];
250  isolationtrackThresholdB_ = trkIsoEndcapRadiusB_[2];
251  trackLipRadiusB_ = trkIsoEndcapRadiusB_[3];
252  trackD0RadiusB_ = trkIsoEndcapRadiusB_[4];
253  isolationtrackEtaSliceB_ = trkIsoEndcapRadiusB_[5];
254 
255  photonEcalRecHitConeInnerRadiusB_ = ecalIsoEndcapRadiusB_[0];
256  photonEcalRecHitConeOuterRadiusB_ = ecalIsoEndcapRadiusB_[1];
257  photonEcalRecHitEtaSliceB_ = ecalIsoEndcapRadiusB_[2];
258  photonEcalRecHitThreshEB_ = ecalIsoEndcapRadiusB_[3];
259  photonEcalRecHitThreshEtB_ = ecalIsoEndcapRadiusB_[4];
260  }
261 
262  //Calculate hollow cone track isolation, CONE A
263  int ntrkA = 0;
264  double trkisoA = 0;
265  calculateTrackIso(pho,
266  e,
267  trkisoA,
268  ntrkA,
269  isolationtrackThresholdA_,
270  trackConeOuterRadiusA_,
271  trackConeInnerRadiusA_,
272  isolationtrackEtaSliceA_,
273  trackLipRadiusA_,
274  trackD0RadiusA_);
275 
276  //Calculate solid cone track isolation, CONE A
277  int sntrkA = 0;
278  double strkisoA = 0;
279  calculateTrackIso(pho,
280  e,
281  strkisoA,
282  sntrkA,
283  isolationtrackThresholdA_,
284  trackConeOuterRadiusA_,
285  0.,
286  isolationtrackEtaSliceA_,
287  trackLipRadiusA_,
288  trackD0RadiusA_);
289 
290  phoisolR1.nTrkHollowCone = ntrkA;
291  phoisolR1.trkSumPtHollowCone = trkisoA;
292  phoisolR1.nTrkSolidCone = sntrkA;
293  phoisolR1.trkSumPtSolidCone = strkisoA;
294 
295  //Calculate hollow cone track isolation, CONE B
296  int ntrkB = 0;
297  double trkisoB = 0;
298  calculateTrackIso(pho,
299  e,
300  trkisoB,
301  ntrkB,
302  isolationtrackThresholdB_,
303  trackConeOuterRadiusB_,
304  trackConeInnerRadiusB_,
305  isolationtrackEtaSliceB_,
306  trackLipRadiusB_,
307  trackD0RadiusB_);
308 
309  //Calculate solid cone track isolation, CONE B
310  int sntrkB = 0;
311  double strkisoB = 0;
312  calculateTrackIso(pho,
313  e,
314  strkisoB,
315  sntrkB,
316  isolationtrackThresholdB_,
317  trackConeOuterRadiusB_,
318  0.,
319  isolationtrackEtaSliceB_,
320  trackLipRadiusB_,
321  trackD0RadiusB_);
322 
323  phoisolR2.nTrkHollowCone = ntrkB;
324  phoisolR2.trkSumPtHollowCone = trkisoB;
325  phoisolR2.nTrkSolidCone = sntrkB;
326  phoisolR2.trkSumPtSolidCone = strkisoB;
327 
328  // std::cout << "Output from solid cone track isolation: ";
329  // std::cout << " Sum pT: " << strkiso << " ntrk: " << sntrk << std::endl;
330 
331  double EcalRecHitIsoA = calculateEcalRecHitIso(pho,
332  e,
333  es,
334  photonEcalRecHitConeOuterRadiusA_,
335  photonEcalRecHitConeInnerRadiusA_,
336  photonEcalRecHitEtaSliceA_,
337  photonEcalRecHitThreshEA_,
338  photonEcalRecHitThreshEtA_,
339  vetoClusteredEcalHits_,
340  useNumCrystals_);
341  phoisolR1.ecalRecHitSumEt = EcalRecHitIsoA;
342 
343  double EcalRecHitIsoB = calculateEcalRecHitIso(pho,
344  e,
345  es,
346  photonEcalRecHitConeOuterRadiusB_,
347  photonEcalRecHitConeInnerRadiusB_,
348  photonEcalRecHitEtaSliceB_,
349  photonEcalRecHitThreshEB_,
350  photonEcalRecHitThreshEtB_,
351  vetoClusteredEcalHits_,
352  useNumCrystals_);
353  phoisolR2.ecalRecHitSumEt = EcalRecHitIsoB;
354 
355  if (not hbheRecHitsTag_.isUninitialized()) {
356  auto const& hbheRecHits = e.get(hbheRecHitsTag_);
357 
358  auto fcone = [this,
359  pho,
360  &caloGeometry,
361  &hcalTopo = *hcalTopology,
362  &hcalQual = *hcalChannelQuality,
363  &hcalSev = *hcalSevLvlComputer,
364  &towerMap,
365  &hbheRecHits](double outer, double inner, int depth) {
366  return calculateHcalRecHitIso<false>(
367  pho, caloGeometry, hcalTopo, hcalQual, hcalSev, towerMap, hbheRecHits, outer, inner, depth);
368  };
369 
370  auto fbc = [this,
371  pho,
372  &caloGeometry,
373  &hcalTopo = *hcalTopology,
374  &hcalQual = *hcalChannelQuality,
375  &hcalSev = *hcalSevLvlComputer,
376  &towerMap,
377  &hbheRecHits](double outer, int depth) {
378  return calculateHcalRecHitIso<true>(
379  pho, caloGeometry, hcalTopo, hcalQual, hcalSev, towerMap, hbheRecHits, outer, 0., depth);
380  };
381 
382  for (size_t id = 0; id < phoisolR1.hcalRecHitSumEt.size(); ++id) {
383  const auto& outerA = detector == EcalBarrel ? hcalIsoOuterRadAEB_[id] : hcalIsoOuterRadAEE_[id];
384  const auto& outerB = detector == EcalBarrel ? hcalIsoOuterRadBEB_[id] : hcalIsoOuterRadBEE_[id];
385  const auto& innerA = detector == EcalBarrel ? hcalIsoInnerRadAEB_[id] : hcalIsoInnerRadAEE_[id];
386  const auto& innerB = detector == EcalBarrel ? hcalIsoInnerRadBEB_[id] : hcalIsoInnerRadBEE_[id];
387 
388  phoisolR1.hcalRecHitSumEt[id] = fcone(outerA, innerA, id + 1);
389  phoisolR2.hcalRecHitSumEt[id] = fcone(outerB, innerB, id + 1);
390 
391  phoisolR1.hcalRecHitSumEtBc[id] = fbc(outerA, id + 1);
392  phoisolR2.hcalRecHitSumEtBc[id] = fbc(outerB, id + 1);
393  }
394  }
395 
396  phoisolR1.pre7DepthHcal = false;
397  phoisolR2.pre7DepthHcal = false;
398 }
399 
400 void PhotonIsolationCalculator::classify(const reco::Photon* photon,
401  bool& isEBPho,
402  bool& isEEPho,
403  bool& isEBEtaGap,
404  bool& isEBPhiGap,
405  bool& isEERingGap,
406  bool& isEEDeeGap,
407  bool& isEBEEGap) {
408  const reco::CaloCluster& seedCluster = *(photon->superCluster()->seed());
409  // following line is temporary until the D. Evans or F. Ferri submit the new tag for ClusterAlgos
410  // DEfinitive will be
411  // DetId seedXtalId = scRef->seed()->seed();
412  DetId seedXtalId = seedCluster.hitsAndFractions()[0].first;
413  int detector = seedXtalId.subdetId();
414 
415  //Set fiducial flags for this photon.
416  double eta = photon->superCluster()->position().eta();
417  double feta = fabs(eta);
418 
419  if (detector == EcalBarrel) {
420  isEBPho = true;
421  if (EBDetId::isNextToEtaBoundary(EBDetId(seedXtalId))) {
422  if (fabs(feta - 1.479) < .1)
423  isEBEEGap = true;
424  else
425  isEBEtaGap = true;
426  }
427 
428  if (EBDetId::isNextToPhiBoundary(EBDetId(seedXtalId)))
429  isEBPhiGap = true;
430 
431  } else if (detector == EcalEndcap) {
432  isEEPho = true;
433  if (EEDetId::isNextToRingBoundary(EEDetId(seedXtalId))) {
434  if (fabs(feta - 1.479) < .1)
435  isEBEEGap = true;
436  else
437  isEERingGap = true;
438  }
439 
440  if (EEDetId::isNextToDBoundary(EEDetId(seedXtalId)))
441  isEEDeeGap = true;
442  }
443 }
444 
445 void PhotonIsolationCalculator::calculateTrackIso(const reco::Photon* photon,
446  const edm::Event& e,
447  double& trkCone,
448  int& ntrkCone,
449  double pTThresh,
450  double RCone,
451  double RinnerCone,
452  double etaSlice,
453  double lip,
454  double d0) const {
455  ntrkCone = 0;
456  trkCone = 0;
457  //get the tracks
458  edm::Handle<reco::TrackCollection> tracks;
459  e.getByToken(trackInputTag_, tracks);
460  if (!tracks.isValid()) {
461  return;
462  }
463  const reco::TrackCollection* trackCollection = tracks.product();
464  //Photon Eta and Phi. Hope these are correct.
465  reco::BeamSpot vertexBeamSpot;
466  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
467  e.getByToken(beamSpotProducerTag_, recoBeamSpotHandle);
468  vertexBeamSpot = *recoBeamSpotHandle;
469 
470  PhotonTkIsolation phoIso(RCone,
471  RinnerCone,
472  etaSlice,
473  pTThresh,
474  lip,
475  d0,
476  trackCollection,
477  math::XYZPoint(vertexBeamSpot.x0(), vertexBeamSpot.y0(), vertexBeamSpot.z0()));
478 
479  std::pair<int, double> res = phoIso.getIso(photon);
480  ntrkCone = res.first;
481  trkCone = res.second;
482 }
483 
484 double PhotonIsolationCalculator::calculateEcalRecHitIso(const reco::Photon* photon,
485  const edm::Event& iEvent,
486  const edm::EventSetup& iSetup,
487  double RCone,
488  double RConeInner,
489  double etaSlice,
490  double eMin,
491  double etMin,
492  bool vetoClusteredHits,
493  bool useNumXtals) const {
494  edm::Handle<EcalRecHitCollection> ecalhitsCollEB;
495  edm::Handle<EcalRecHitCollection> ecalhitsCollEE;
496  iEvent.getByToken(endcapecalCollection_, ecalhitsCollEE);
497 
498  iEvent.getByToken(barrelecalCollection_, ecalhitsCollEB);
499 
500  const EcalRecHitCollection* rechitsCollectionEE_ = ecalhitsCollEE.product();
501  const EcalRecHitCollection* rechitsCollectionEB_ = ecalhitsCollEB.product();
502 
503  const EcalSeverityLevelAlgo* sevLevel = &iSetup.getData(ecalSevLvlToken_);
504 
505  edm::ESHandle<CaloGeometry> geoHandle = iSetup.getHandle(caloGeometryToken_);
506 
507  EgammaRecHitIsolation phoIsoEB(
508  RCone, RConeInner, etaSlice, etMin, eMin, geoHandle, *rechitsCollectionEB_, sevLevel, DetId::Ecal);
509 
510  phoIsoEB.setVetoClustered(vetoClusteredHits);
511  phoIsoEB.setUseNumCrystals(useNumXtals);
512  phoIsoEB.doSeverityChecks(ecalhitsCollEB.product(), severityExclEB_);
513  phoIsoEB.doFlagChecks(flagsEB_);
514  double ecalIsolEB = phoIsoEB.getEtSum(photon);
515 
516  EgammaRecHitIsolation phoIsoEE(
517  RCone, RConeInner, etaSlice, etMin, eMin, geoHandle, *rechitsCollectionEE_, sevLevel, DetId::Ecal);
518 
519  phoIsoEE.setVetoClustered(vetoClusteredHits);
520  phoIsoEE.setUseNumCrystals(useNumXtals);
521  phoIsoEE.doSeverityChecks(ecalhitsCollEE.product(), severityExclEE_);
522  phoIsoEE.doFlagChecks(flagsEE_);
523 
524  double ecalIsolEE = phoIsoEE.getEtSum(photon);
525  // delete phoIso;
526  double ecalIsol = ecalIsolEB + ecalIsolEE;
527 
528  return ecalIsol;
529 }
530 
531 template <bool isoBC>
532 double PhotonIsolationCalculator::calculateHcalRecHitIso(const reco::Photon* photon,
533  const CaloGeometry& geometry,
534  const HcalTopology& hcalTopology,
535  const HcalChannelQuality& hcalChStatus,
536  const HcalSeverityLevelComputer& hcalSevLvlComputer,
537  const CaloTowerConstituentsMap& towerMap,
538  const HBHERecHitCollection& hbheRecHits,
539  double RCone,
540  double RConeInner,
541  int depth) const {
542  const EgammaHcalIsolation::arrayHB e04{{0., 0., 0., 0.}};
543  const EgammaHcalIsolation::arrayHE e07{{0., 0., 0., 0., 0., 0., 0.}};
544 
545  if constexpr (isoBC) {
546  auto hcaliso = EgammaHcalIsolation(EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
547  RCone,
548  EgammaHcalIsolation::InclusionRule::isBehindClusterSeed,
549  RConeInner,
550  hcalIsoEThresHB_,
551  e04,
552  maxHcalSeverity_,
553  hcalIsoEThresHE_,
554  e07,
555  maxHcalSeverity_,
556  hbheRecHits,
557  geometry,
558  hcalTopology,
559  hcalChStatus,
560  hcalSevLvlComputer,
561  towerMap);
562 
563  return hcaliso.getHcalEtSumBc(photon, depth);
564  } else {
565  auto hcaliso = EgammaHcalIsolation(EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
566  RCone,
567  EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
568  RConeInner,
569  hcalIsoEThresHB_,
570  e04,
571  maxHcalSeverity_,
572  hcalIsoEThresHE_,
573  e07,
574  maxHcalSeverity_,
575  hbheRecHits,
576  geometry,
577  hcalTopology,
578  hcalChStatus,
579  hcalSevLvlComputer,
580  towerMap);
581 
582  return hcaliso.getHcalEtSum(photon, depth);
583  }
584 }
PhotonIsolationCalculator::hcalIsoInnerRadBEB_
std::array< double, 7 > hcalIsoInnerRadBEB_
Definition: PhotonIsolationCalculator.h:127
reco::BeamSpot::z0
double z0() const
z coordinate
Definition: BeamSpot.h:65
PhotonIsolationCalculator::trkIsoBarrelRadiusA_
std::array< double, 6 > trkIsoBarrelRadiusA_
Definition: PhotonIsolationCalculator.h:114
EBDetId::isNextToEtaBoundary
static bool isNextToEtaBoundary(EBDetId id)
Definition: EBDetId.cc:108
PhotonIsolationCalculator::calculate
void calculate(const reco::Photon *, const edm::Event &, const edm::EventSetup &es, reco::Photon::FiducialFlags &phofid, reco::Photon::IsolationVariables &phoisolR03, reco::Photon::IsolationVariables &phoisolR04) const
Definition: PhotonIsolationCalculator.cc:145
reco::Photon::IsolationVariables::hcalRecHitSumEtBc
std::array< float, 7 > hcalRecHitSumEtBc
Definition: Photon.h:420
mps_fire.i
i
Definition: mps_fire.py:428
HLT_FULL_cff.lip
tuple lip
Definition: HLT_FULL_cff.py:39859
HLT_FULL_cff.etMin
tuple etMin
Definition: HLT_FULL_cff.py:8483
gpuClustering::id
uint16_t *__restrict__ id
Definition: gpuClusterChargeCut.h:21
PhotonIsolationCalculator::hcalIsoEThresHE_
EgammaHcalIsolation::arrayHE hcalIsoEThresHE_
Definition: PhotonIsolationCalculator.h:137
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:539
EgammaRecHitIsolation::getEtSum
double getEtSum(const reco::Candidate *emObject) const
Definition: EgammaRecHitIsolation.h:39
ootPhotons_cff.hbheRecHits
string hbheRecHits
Definition: ootPhotons_cff.py:41
PhotonIsolationCalculator::ecalSevLvlToken_
edm::ESGetToken< EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd > ecalSevLvlToken_
Definition: PhotonIsolationCalculator.h:105
reco::TrackCollection
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:14
PhotonIsolationCalculator::caloGeometryToken_
edm::ESGetToken< CaloGeometry, CaloGeometryRecord > caloGeometryToken_
Definition: PhotonIsolationCalculator.h:100
tracks
auto const & tracks
cannot be loose
Definition: CAHitNtupletGeneratorKernelsImpl.h:148
EgammaRecHitIsolation::doFlagChecks
void doFlagChecks(const std::vector< int > &v)
Definition: EgammaRecHitIsolation.h:54
edm::EventSetup::getData
bool getData(T &iHolder) const
Definition: EventSetup.h:122
PhotonIsolationCalculator::hbheRecHitsTag_
edm::EDGetTokenT< HBHERecHitCollection > hbheRecHitsTag_
Definition: PhotonIsolationCalculator.h:98
EgammaRecHitIsolation::setUseNumCrystals
void setUseNumCrystals(bool b=true)
Definition: EgammaRecHitIsolation.h:45
iEvent
int iEvent
Definition: GenABIO.cc:224
reco::CaloCluster::hitsAndFractions
const std::vector< std::pair< DetId, float > > & hitsAndFractions() const
Definition: CaloCluster.h:210
reco::Photon::superCluster
reco::SuperClusterRef superCluster() const override
Ref to SuperCluster.
PhotonIsolationCalculator::hcalIsoOuterRadAEE_
std::array< double, 7 > hcalIsoOuterRadAEE_
Definition: PhotonIsolationCalculator.h:131
PhotonIsolationCalculator::setup
void setup(const edm::ParameterSet &conf, std::vector< int > const &flagsEB_, std::vector< int > const &flagsEE_, std::vector< int > const &severitiesEB_, std::vector< int > const &severitiesEE_, edm::ConsumesCollector &&iC)
Definition: PhotonIsolationCalculator.cc:29
PhotonIsolationCalculator::hcalIsoInnerRadAEE_
std::array< double, 7 > hcalIsoInnerRadAEE_
Definition: PhotonIsolationCalculator.h:130
EBDetId::isNextToPhiBoundary
static bool isNextToPhiBoundary(EBDetId id)
Definition: EBDetId.cc:113
PhotonIsolationCalculator::hcalIsoOuterRadBEE_
std::array< double, 7 > hcalIsoOuterRadBEE_
Definition: PhotonIsolationCalculator.h:134
DetId::subdetId
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:48
EEDetId::isNextToRingBoundary
static bool isNextToRingBoundary(EEDetId id)
Definition: EEDetId.cc:284
edm::Event::get
bool get(ProductID const &oid, Handle< PROD > &result) const
Definition: Event.h:346
PhotonIsolationCalculator::ecalIsoEndcapRadiusB_
std::array< double, 5 > ecalIsoEndcapRadiusB_
Definition: PhotonIsolationCalculator.h:122
PhotonIsolationCalculator::hcalIsoOuterRadAEB_
std::array< double, 7 > hcalIsoOuterRadAEB_
Definition: PhotonIsolationCalculator.h:125
EEDetId::isNextToDBoundary
static bool isNextToDBoundary(EEDetId id)
Definition: EEDetId.cc:279
edm::HandleBase::isValid
bool isValid() const
Definition: HandleBase.h:70
PhotonIsolationCalculator::trkIsoEndcapRadiusB_
std::array< double, 6 > trkIsoEndcapRadiusB_
Definition: PhotonIsolationCalculator.h:121
PhotonTkIsolation::getIso
std::pair< int, float > getIso(const reco::Candidate *) const
Definition: PhotonTkIsolation.cc:52
EgammaRecHitIsolation::doSeverityChecks
void doSeverityChecks(const EcalRecHitCollection *const recHits, const std::vector< int > &v)
Definition: EgammaRecHitIsolation.h:47
HLT_FULL_cff.depth
tuple depth
Definition: HLT_FULL_cff.py:12859
DetId
Definition: DetId.h:17
PhotonIsolationCalculator::hcalIsoEThresHB_
EgammaHcalIsolation::arrayHB hcalIsoEThresHB_
Definition: PhotonIsolationCalculator.h:136
fbc
cc *****************************************************cc the common blocks pinput and cwdidth are for input parameters cc these parameters needed to be interfaced to other program common pinput fbc
Definition: inclcon.h:4
d0
static constexpr float d0
Definition: L1EGammaCrystalsEmulatorProducer.cc:85
PhotonIsolationCalculator::hcalChannelQualityToken_
edm::ESGetToken< HcalChannelQuality, HcalChannelQualityRcd > hcalChannelQualityToken_
Definition: PhotonIsolationCalculator.h:102
HLT_FULL_cff.trackCollection
tuple trackCollection
Definition: HLT_FULL_cff.py:66975
PhotonIsolationCalculator::ecalIsoBarrelRadiusB_
std::array< double, 5 > ecalIsoBarrelRadiusB_
Definition: PhotonIsolationCalculator.h:117
PhotonIsolationCalculator::ecalIsoEndcapRadiusA_
std::array< double, 5 > ecalIsoEndcapRadiusA_
Definition: PhotonIsolationCalculator.h:120
PhotonIsolationCalculator::towerMapToken_
edm::ESGetToken< CaloTowerConstituentsMap, CaloGeometryRecord > towerMapToken_
Definition: PhotonIsolationCalculator.h:104
PhotonIsolationCalculator::trkIsoBarrelRadiusB_
std::array< double, 6 > trkIsoBarrelRadiusB_
Definition: PhotonIsolationCalculator.h:116
edm::Handle::product
T const * product() const
Definition: Handle.h:70
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
reco::Photon::IsolationVariables::hcalRecHitSumEt
std::array< float, 7 > hcalRecHitSumEt
Definition: Photon.h:419
PhotonIsolationCalculator::ecalIsoBarrelRadiusA_
std::array< double, 5 > ecalIsoBarrelRadiusA_
Definition: PhotonIsolationCalculator.h:115
PhotonIsolationCalculator::calculateEcalRecHitIso
double calculateEcalRecHitIso(const reco::Photon *photon, const edm::Event &iEvent, const edm::EventSetup &iSetup, double RCone, double RConeInner, double etaSlice, double eMin, double etMin, bool vetoClusteredHits, bool useNumCrystals) const
Definition: PhotonIsolationCalculator.cc:484
PhotonIsolationCalculator::calculateTrackIso
void calculateTrackIso(const reco::Photon *photon, const edm::Event &e, double &trkCone, int &ntrkCone, double pTThresh=0, double RCone=.4, double RinnerCone=.1, double etaSlice=0.015, double lip=0.2, double d0=0.1) const
Definition: PhotonIsolationCalculator.cc:445
PhotonIsolationCalculator::hcalTopologyToken_
edm::ESGetToken< HcalTopology, HcalRecNumberingRecord > hcalTopologyToken_
Definition: PhotonIsolationCalculator.h:101
PhotonIsolationCalculator::calculateHcalRecHitIso
double calculateHcalRecHitIso(const reco::Photon *photon, const CaloGeometry &geometry, const HcalTopology &hcalTopology, const HcalChannelQuality &hcalChStatus, const HcalSeverityLevelComputer &hcalSevLvlComputer, const CaloTowerConstituentsMap &towerMap, const HBHERecHitCollection &hbheRecHits, double RCone, double RConeInner, int depth) const
Definition: PhotonIsolationCalculator.cc:532
PhotonIsolationCalculator::moduleEtaBoundary_
std::vector< double > moduleEtaBoundary_
Definition: PhotonIsolationCalculator.h:110
EgammaRecHitIsolation::setVetoClustered
void setVetoClustered(bool b=true)
Definition: EgammaRecHitIsolation.h:46
PhotonIsolationCalculator::hcalSevLvlComputerToken_
edm::ESGetToken< HcalSeverityLevelComputer, HcalSeverityLevelComputerRcd > hcalSevLvlComputerToken_
Definition: PhotonIsolationCalculator.h:103
reco::BeamSpot::y0
double y0() const
y coordinate
Definition: BeamSpot.h:63
PhotonIsolationCalculator::hcalIsoInnerRadBEE_
std::array< double, 7 > hcalIsoInnerRadBEE_
Definition: PhotonIsolationCalculator.h:133
edm::EventSetup::getHandle
ESHandle< T > getHandle(const ESGetToken< T, R > &iToken) const
Definition: EventSetup.h:151
EgammaHcalIsolation::arrayHB
std::array< double, 4 > arrayHB
Definition: EgammaHcalIsolation.h:41
PhotonIsolationCalculator::classify
static void classify(const reco::Photon *photon, bool &isEBPho, bool &isEEPho, bool &isEBEtaGap, bool &isEBPhiGap, bool &isEERingGap, bool &isEEDeeGap, bool &isEBEEGap)
Definition: PhotonIsolationCalculator.cc:400
PhotonIsolationCalculator::hcalIsoOuterRadBEB_
std::array< double, 7 > hcalIsoOuterRadBEB_
Definition: PhotonIsolationCalculator.h:128
PhotonIsolationCalculator::hcalIsoInnerRadAEB_
std::array< double, 7 > hcalIsoInnerRadAEB_
Definition: PhotonIsolationCalculator.h:124
PhotonIsolationCalculator::trkIsoEndcapRadiusA_
std::array< double, 6 > trkIsoEndcapRadiusA_
Definition: PhotonIsolationCalculator.h:119
reco::BeamSpot::x0
double x0() const
x coordinate
Definition: BeamSpot.h:61
EgammaHcalIsolation::arrayHE
std::array< double, 7 > arrayHE
Definition: EgammaHcalIsolation.h:42
Generated for CMSSW Reference Manual by doxygen 1.8.5