CMS 3D CMS Logo

PhotonIsolationCalculator.cc
Go to the documentation of this file.
1 
23 
29 //#include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
30 
37 
40 
41 #include <string>
42 #include <TMath.h>
43 
45  std::vector<int> const& flagsEB,
46  std::vector<int> const& flagsEE,
47  std::vector<int> const& severitiesEB,
48  std::vector<int> const& severitiesEE,
50  trackInputTag_ = iC.consumes<reco::TrackCollection>(conf.getParameter<edm::InputTag>("trackProducer"));
51  beamSpotProducerTag_ = iC.consumes<reco::BeamSpot>(conf.getParameter<edm::InputTag>("beamSpotProducer"));
53  iC.consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("barrelEcalRecHitCollection"));
55  iC.consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("endcapEcalRecHitCollection"));
56  auto hcRHC = conf.getParameter<edm::InputTag>("HcalRecHitCollection");
57  if (not hcRHC.label().empty())
58  hcalCollection_ = iC.consumes<CaloTowerCollection>(hcRHC);
59 
60  // gsfRecoInputTag_ = conf.getParameter<edm::InputTag>("GsfRecoCollection");
61  modulePhiBoundary_ = conf.getParameter<double>("modulePhiBoundary");
62  moduleEtaBoundary_ = conf.getParameter<std::vector<double> >("moduleEtaBoundary");
63  //
64  vetoClusteredEcalHits_ = conf.getParameter<bool>("vetoClustered");
65  useNumCrystals_ = conf.getParameter<bool>("useNumCrystals");
66 
68  int i = 0;
69  trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("TrackConeOuterRadiusA_Barrel"));
70  trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("TrackConeInnerRadiusA_Barrel"));
71  trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("isolationtrackThresholdA_Barrel"));
72  trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("longImpactParameterA_Barrel"));
73  trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("transImpactParameterA_Barrel"));
74  trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("isolationtrackEtaSliceA_Barrel"));
75 
76  i = 0;
77  ecalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitInnerRadiusA_Barrel"));
78  ecalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitOuterRadiusA_Barrel"));
79  ecalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitEtaSliceA_Barrel"));
80  ecalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitThreshEA_Barrel"));
81  ecalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitThreshEtA_Barrel"));
82 
83  i = 0;
84  hcalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("HcalTowerInnerRadiusA_Barrel"));
85  hcalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("HcalTowerOuterRadiusA_Barrel"));
86  hcalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("HcalTowerThreshEA_Barrel"));
87  hcalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("HcalDepth1TowerInnerRadiusA_Barrel"));
88  hcalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("HcalDepth1TowerOuterRadiusA_Barrel"));
89  hcalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("HcalDepth1TowerThreshEA_Barrel"));
90  hcalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("HcalDepth2TowerInnerRadiusA_Barrel"));
91  hcalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("HcalDepth2TowerOuterRadiusA_Barrel"));
92  hcalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("HcalDepth2TowerThreshEA_Barrel"));
93 
94  i = 0;
95  trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("TrackConeOuterRadiusB_Barrel"));
96  trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("TrackConeInnerRadiusB_Barrel"));
97  trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("isolationtrackThresholdB_Barrel"));
98  trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("longImpactParameterB_Barrel"));
99  trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("transImpactParameterB_Barrel"));
100  trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("isolationtrackEtaSliceB_Barrel"));
101 
102  i = 0;
103  ecalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitInnerRadiusB_Barrel"));
104  ecalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitOuterRadiusB_Barrel"));
105  ecalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitEtaSliceB_Barrel"));
106  ecalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitThreshEB_Barrel"));
107  ecalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitThreshEtB_Barrel"));
108 
109  i = 0;
110  hcalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("HcalTowerInnerRadiusB_Barrel"));
111  hcalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("HcalTowerOuterRadiusB_Barrel"));
112  hcalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("HcalTowerThreshEB_Barrel"));
113  hcalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("HcalDepth1TowerInnerRadiusB_Barrel"));
114  hcalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("HcalDepth1TowerOuterRadiusB_Barrel"));
115  hcalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("HcalDepth1TowerThreshEB_Barrel"));
116  hcalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("HcalDepth2TowerInnerRadiusB_Barrel"));
117  hcalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("HcalDepth2TowerOuterRadiusB_Barrel"));
118  hcalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("HcalDepth2TowerThreshEB_Barrel"));
119 
121  i = 0;
122  trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("TrackConeOuterRadiusA_Endcap"));
123  trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("TrackConeInnerRadiusA_Endcap"));
124  trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("isolationtrackThresholdA_Endcap"));
125  trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("longImpactParameterA_Endcap"));
126  trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("transImpactParameterA_Endcap"));
127  trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("isolationtrackEtaSliceA_Endcap"));
128 
129  i = 0;
130  ecalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitInnerRadiusA_Endcap"));
131  ecalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitOuterRadiusA_Endcap"));
132  ecalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitEtaSliceA_Endcap"));
133  ecalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitThreshEA_Endcap"));
134  ecalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitThreshEtA_Endcap"));
135 
136  i = 0;
137  hcalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("HcalTowerInnerRadiusA_Endcap"));
138  hcalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("HcalTowerOuterRadiusA_Endcap"));
139  hcalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("HcalTowerThreshEA_Endcap"));
140  hcalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("HcalDepth1TowerInnerRadiusA_Endcap"));
141  hcalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("HcalDepth1TowerOuterRadiusA_Endcap"));
142  hcalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("HcalDepth1TowerThreshEA_Endcap"));
143  hcalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("HcalDepth2TowerInnerRadiusA_Endcap"));
144  hcalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("HcalDepth2TowerOuterRadiusA_Endcap"));
145  hcalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("HcalDepth2TowerThreshEA_Endcap"));
146 
147  i = 0;
148  trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("TrackConeOuterRadiusB_Endcap"));
149  trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("TrackConeInnerRadiusB_Endcap"));
150  trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("isolationtrackThresholdB_Endcap"));
151  trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("longImpactParameterB_Endcap"));
152  trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("transImpactParameterB_Endcap"));
153  trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("isolationtrackEtaSliceB_Endcap"));
154 
155  i = 0;
156  ecalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitInnerRadiusB_Endcap"));
157  ecalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitOuterRadiusB_Endcap"));
158  ecalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitEtaSliceB_Endcap"));
159  ecalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitThreshEB_Endcap"));
160  ecalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitThreshEtB_Endcap"));
161 
162  i = 0;
163  hcalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("HcalTowerInnerRadiusB_Endcap"));
164  hcalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("HcalTowerOuterRadiusB_Endcap"));
165  hcalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("HcalTowerThreshEB_Endcap"));
166  hcalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("HcalDepth1TowerInnerRadiusB_Endcap"));
167  hcalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("HcalDepth1TowerOuterRadiusB_Endcap"));
168  hcalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("HcalDepth1TowerThreshEB_Endcap"));
169  hcalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("HcalDepth2TowerInnerRadiusB_Endcap"));
170  hcalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("HcalDepth2TowerOuterRadiusB_Endcap"));
171  hcalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("HcalDepth2TowerThreshEB_Endcap"));
172 
173  //Pick up the variables for the spike removal
174  flagsEB_ = flagsEB;
175  flagsEE_ = flagsEE;
176  severityExclEB_ = severitiesEB;
177  severityExclEE_ = severitiesEE;
178 }
179 
181  const edm::Event& e,
182  const edm::EventSetup& es,
185  reco::Photon::IsolationVariables& phoisolR2) const {
186  //Get fiducial flags. This does not really belong here
187  bool isEBPho = false;
188  bool isEEPho = false;
189  bool isEBEtaGap = false;
190  bool isEBPhiGap = false;
191  bool isEERingGap = false;
192  bool isEEDeeGap = false;
193  bool isEBEEGap = false;
194  classify(pho, isEBPho, isEEPho, isEBEtaGap, isEBPhiGap, isEERingGap, isEEDeeGap, isEBEEGap);
195  phofid.isEB = isEBPho;
196  phofid.isEE = isEEPho;
197  phofid.isEBEtaGap = isEBEtaGap;
198  phofid.isEBPhiGap = isEBPhiGap;
199  phofid.isEERingGap = isEERingGap;
200  phofid.isEEDeeGap = isEEDeeGap;
201  phofid.isEBEEGap = isEBEEGap;
202 
203  // Calculate isolation variables. cone sizes and thresholds
204  // are set for Barrel and endcap separately
205 
206  reco::SuperClusterRef scRef = pho->superCluster();
207  const reco::BasicCluster& seedCluster = *(scRef->seed());
208  DetId seedXtalId = seedCluster.hitsAndFractions()[0].first;
209  int detector = seedXtalId.subdetId();
210 
211  //Isolation parameters variables
212  double photonEcalRecHitConeInnerRadiusA_;
213  double photonEcalRecHitConeOuterRadiusA_;
214  double photonEcalRecHitEtaSliceA_;
215  double photonEcalRecHitThreshEA_;
216  double photonEcalRecHitThreshEtA_;
217  double photonHcalTowerConeInnerRadiusA_;
218  double photonHcalTowerConeOuterRadiusA_;
219  double photonHcalTowerThreshEA_;
220  double photonHcalDepth1TowerConeInnerRadiusA_;
221  double photonHcalDepth1TowerConeOuterRadiusA_;
222  double photonHcalDepth1TowerThreshEA_;
223  double photonHcalDepth2TowerConeInnerRadiusA_;
224  double photonHcalDepth2TowerConeOuterRadiusA_;
225  double photonHcalDepth2TowerThreshEA_;
226  double trackConeOuterRadiusA_;
227  double trackConeInnerRadiusA_;
228  double isolationtrackThresholdA_;
229  double isolationtrackEtaSliceA_;
230  double trackLipRadiusA_;
231  double trackD0RadiusA_;
232  double photonEcalRecHitConeInnerRadiusB_;
233  double photonEcalRecHitConeOuterRadiusB_;
234  double photonEcalRecHitEtaSliceB_;
235  double photonEcalRecHitThreshEB_;
236  double photonEcalRecHitThreshEtB_;
237  double photonHcalTowerConeInnerRadiusB_;
238  double photonHcalTowerConeOuterRadiusB_;
239  double photonHcalTowerThreshEB_;
240  double photonHcalDepth1TowerConeInnerRadiusB_;
241  double photonHcalDepth1TowerConeOuterRadiusB_;
242  double photonHcalDepth1TowerThreshEB_;
243  double photonHcalDepth2TowerConeInnerRadiusB_;
244  double photonHcalDepth2TowerConeOuterRadiusB_;
245  double photonHcalDepth2TowerThreshEB_;
246  double trackConeOuterRadiusB_;
247  double trackConeInnerRadiusB_;
248  double isolationtrackThresholdB_;
249  double isolationtrackEtaSliceB_;
250  double trackLipRadiusB_;
251  double trackD0RadiusB_;
252 
253  if (detector == EcalBarrel) {
254  trackConeOuterRadiusA_ = trkIsoBarrelRadiusA_[0];
255  trackConeInnerRadiusA_ = trkIsoBarrelRadiusA_[1];
256  isolationtrackThresholdA_ = trkIsoBarrelRadiusA_[2];
257  trackLipRadiusA_ = trkIsoBarrelRadiusA_[3];
258  trackD0RadiusA_ = trkIsoBarrelRadiusA_[4];
259  isolationtrackEtaSliceA_ = trkIsoBarrelRadiusA_[5];
260 
261  photonEcalRecHitConeInnerRadiusA_ = ecalIsoBarrelRadiusA_[0];
262  photonEcalRecHitConeOuterRadiusA_ = ecalIsoBarrelRadiusA_[1];
263  photonEcalRecHitEtaSliceA_ = ecalIsoBarrelRadiusA_[2];
264  photonEcalRecHitThreshEA_ = ecalIsoBarrelRadiusA_[3];
265  photonEcalRecHitThreshEtA_ = ecalIsoBarrelRadiusA_[4];
266 
267  photonHcalTowerConeInnerRadiusA_ = hcalIsoBarrelRadiusA_[0];
268  photonHcalTowerConeOuterRadiusA_ = hcalIsoBarrelRadiusA_[1];
269  photonHcalTowerThreshEA_ = hcalIsoBarrelRadiusA_[2];
270  photonHcalDepth1TowerConeInnerRadiusA_ = hcalIsoBarrelRadiusA_[3];
271  photonHcalDepth1TowerConeOuterRadiusA_ = hcalIsoBarrelRadiusA_[4];
272  photonHcalDepth1TowerThreshEA_ = hcalIsoBarrelRadiusA_[5];
273  photonHcalDepth2TowerConeInnerRadiusA_ = hcalIsoBarrelRadiusA_[6];
274  photonHcalDepth2TowerConeOuterRadiusA_ = hcalIsoBarrelRadiusA_[7];
275  photonHcalDepth2TowerThreshEA_ = hcalIsoBarrelRadiusA_[8];
276 
277  trackConeOuterRadiusB_ = trkIsoBarrelRadiusB_[0];
278  trackConeInnerRadiusB_ = trkIsoBarrelRadiusB_[1];
279  isolationtrackThresholdB_ = trkIsoBarrelRadiusB_[2];
280  trackLipRadiusB_ = trkIsoBarrelRadiusB_[3];
281  trackD0RadiusB_ = trkIsoBarrelRadiusB_[4];
282  isolationtrackEtaSliceB_ = trkIsoBarrelRadiusB_[5];
283 
284  photonEcalRecHitConeInnerRadiusB_ = ecalIsoBarrelRadiusB_[0];
285  photonEcalRecHitConeOuterRadiusB_ = ecalIsoBarrelRadiusB_[1];
286  photonEcalRecHitEtaSliceB_ = ecalIsoBarrelRadiusB_[2];
287  photonEcalRecHitThreshEB_ = ecalIsoBarrelRadiusB_[3];
288  photonEcalRecHitThreshEtB_ = ecalIsoBarrelRadiusB_[4];
289 
290  photonHcalTowerConeInnerRadiusB_ = hcalIsoBarrelRadiusB_[0];
291  photonHcalTowerConeOuterRadiusB_ = hcalIsoBarrelRadiusB_[1];
292  photonHcalTowerThreshEB_ = hcalIsoBarrelRadiusB_[2];
293  photonHcalDepth1TowerConeInnerRadiusB_ = hcalIsoBarrelRadiusB_[3];
294  photonHcalDepth1TowerConeOuterRadiusB_ = hcalIsoBarrelRadiusB_[4];
295  photonHcalDepth1TowerThreshEB_ = hcalIsoBarrelRadiusB_[5];
296  photonHcalDepth2TowerConeInnerRadiusB_ = hcalIsoBarrelRadiusB_[6];
297  photonHcalDepth2TowerConeOuterRadiusB_ = hcalIsoBarrelRadiusB_[7];
298  photonHcalDepth2TowerThreshEB_ = hcalIsoBarrelRadiusB_[8];
299 
300  } else {
301  // detector==EcalEndcap
302 
303  trackConeOuterRadiusA_ = trkIsoEndcapRadiusA_[0];
304  trackConeInnerRadiusA_ = trkIsoEndcapRadiusA_[1];
305  isolationtrackThresholdA_ = trkIsoEndcapRadiusA_[2];
306  trackLipRadiusA_ = trkIsoEndcapRadiusA_[3];
307  trackD0RadiusA_ = trkIsoEndcapRadiusA_[4];
308  isolationtrackEtaSliceA_ = trkIsoEndcapRadiusA_[5];
309 
310  photonEcalRecHitConeInnerRadiusA_ = ecalIsoEndcapRadiusA_[0];
311  photonEcalRecHitConeOuterRadiusA_ = ecalIsoEndcapRadiusA_[1];
312  photonEcalRecHitEtaSliceA_ = ecalIsoEndcapRadiusA_[2];
313  photonEcalRecHitThreshEA_ = ecalIsoEndcapRadiusA_[3];
314  photonEcalRecHitThreshEtA_ = ecalIsoEndcapRadiusA_[4];
315 
316  photonHcalTowerConeInnerRadiusA_ = hcalIsoEndcapRadiusA_[0];
317  photonHcalTowerConeOuterRadiusA_ = hcalIsoEndcapRadiusA_[1];
318  photonHcalTowerThreshEA_ = hcalIsoEndcapRadiusA_[2];
319  photonHcalDepth1TowerConeInnerRadiusA_ = hcalIsoEndcapRadiusA_[3];
320  photonHcalDepth1TowerConeOuterRadiusA_ = hcalIsoEndcapRadiusA_[4];
321  photonHcalDepth1TowerThreshEA_ = hcalIsoEndcapRadiusA_[5];
322  photonHcalDepth2TowerConeInnerRadiusA_ = hcalIsoEndcapRadiusA_[6];
323  photonHcalDepth2TowerConeOuterRadiusA_ = hcalIsoEndcapRadiusA_[7];
324  photonHcalDepth2TowerThreshEA_ = hcalIsoEndcapRadiusA_[8];
325 
326  trackConeOuterRadiusB_ = trkIsoEndcapRadiusB_[0];
327  trackConeInnerRadiusB_ = trkIsoEndcapRadiusB_[1];
328  isolationtrackThresholdB_ = trkIsoEndcapRadiusB_[2];
329  trackLipRadiusB_ = trkIsoEndcapRadiusB_[3];
330  trackD0RadiusB_ = trkIsoEndcapRadiusB_[4];
331  isolationtrackEtaSliceB_ = trkIsoEndcapRadiusB_[5];
332 
333  photonEcalRecHitConeInnerRadiusB_ = ecalIsoEndcapRadiusB_[0];
334  photonEcalRecHitConeOuterRadiusB_ = ecalIsoEndcapRadiusB_[1];
335  photonEcalRecHitEtaSliceB_ = ecalIsoEndcapRadiusB_[2];
336  photonEcalRecHitThreshEB_ = ecalIsoEndcapRadiusB_[3];
337  photonEcalRecHitThreshEtB_ = ecalIsoEndcapRadiusB_[4];
338 
339  photonHcalTowerConeInnerRadiusB_ = hcalIsoEndcapRadiusB_[0];
340  photonHcalTowerConeOuterRadiusB_ = hcalIsoEndcapRadiusB_[1];
341  photonHcalTowerThreshEB_ = hcalIsoEndcapRadiusB_[2];
342  photonHcalDepth1TowerConeInnerRadiusB_ = hcalIsoEndcapRadiusB_[3];
343  photonHcalDepth1TowerConeOuterRadiusB_ = hcalIsoEndcapRadiusB_[4];
344  photonHcalDepth1TowerThreshEB_ = hcalIsoEndcapRadiusB_[5];
345  photonHcalDepth2TowerConeInnerRadiusB_ = hcalIsoEndcapRadiusB_[6];
346  photonHcalDepth2TowerConeOuterRadiusB_ = hcalIsoEndcapRadiusB_[7];
347  photonHcalDepth2TowerThreshEB_ = hcalIsoEndcapRadiusB_[8];
348  }
349 
350  //Calculate hollow cone track isolation, CONE A
351  int ntrkA = 0;
352  double trkisoA = 0;
353  calculateTrackIso(pho,
354  e,
355  trkisoA,
356  ntrkA,
357  isolationtrackThresholdA_,
358  trackConeOuterRadiusA_,
359  trackConeInnerRadiusA_,
360  isolationtrackEtaSliceA_,
361  trackLipRadiusA_,
362  trackD0RadiusA_);
363 
364  //Calculate solid cone track isolation, CONE A
365  int sntrkA = 0;
366  double strkisoA = 0;
367  calculateTrackIso(pho,
368  e,
369  strkisoA,
370  sntrkA,
371  isolationtrackThresholdA_,
372  trackConeOuterRadiusA_,
373  0.,
374  isolationtrackEtaSliceA_,
375  trackLipRadiusA_,
376  trackD0RadiusA_);
377 
378  phoisolR1.nTrkHollowCone = ntrkA;
379  phoisolR1.trkSumPtHollowCone = trkisoA;
380  phoisolR1.nTrkSolidCone = sntrkA;
381  phoisolR1.trkSumPtSolidCone = strkisoA;
382 
383  //Calculate hollow cone track isolation, CONE B
384  int ntrkB = 0;
385  double trkisoB = 0;
386  calculateTrackIso(pho,
387  e,
388  trkisoB,
389  ntrkB,
390  isolationtrackThresholdB_,
391  trackConeOuterRadiusB_,
392  trackConeInnerRadiusB_,
393  isolationtrackEtaSliceB_,
394  trackLipRadiusB_,
395  trackD0RadiusB_);
396 
397  //Calculate solid cone track isolation, CONE B
398  int sntrkB = 0;
399  double strkisoB = 0;
400  calculateTrackIso(pho,
401  e,
402  strkisoB,
403  sntrkB,
404  isolationtrackThresholdB_,
405  trackConeOuterRadiusB_,
406  0.,
407  isolationtrackEtaSliceB_,
408  trackLipRadiusB_,
409  trackD0RadiusB_);
410 
411  phoisolR2.nTrkHollowCone = ntrkB;
412  phoisolR2.trkSumPtHollowCone = trkisoB;
413  phoisolR2.nTrkSolidCone = sntrkB;
414  phoisolR2.trkSumPtSolidCone = strkisoB;
415 
416  // std::cout << "Output from solid cone track isolation: ";
417  // std::cout << " Sum pT: " << strkiso << " ntrk: " << sntrk << std::endl;
418 
419  double EcalRecHitIsoA = calculateEcalRecHitIso(pho,
420  e,
421  es,
422  photonEcalRecHitConeOuterRadiusA_,
423  photonEcalRecHitConeInnerRadiusA_,
424  photonEcalRecHitEtaSliceA_,
425  photonEcalRecHitThreshEA_,
426  photonEcalRecHitThreshEtA_,
429  phoisolR1.ecalRecHitSumEt = EcalRecHitIsoA;
430 
431  double EcalRecHitIsoB = calculateEcalRecHitIso(pho,
432  e,
433  es,
434  photonEcalRecHitConeOuterRadiusB_,
435  photonEcalRecHitConeInnerRadiusB_,
436  photonEcalRecHitEtaSliceB_,
437  photonEcalRecHitThreshEB_,
438  photonEcalRecHitThreshEtB_,
441  phoisolR2.ecalRecHitSumEt = EcalRecHitIsoB;
442 
443  double HcalTowerIsoA = calculateHcalTowerIso(
444  pho, e, es, photonHcalTowerConeOuterRadiusA_, photonHcalTowerConeInnerRadiusA_, photonHcalTowerThreshEA_, -1);
445  phoisolR1.hcalTowerSumEt = HcalTowerIsoA;
446 
447  double HcalTowerIsoB = calculateHcalTowerIso(
448  pho, e, es, photonHcalTowerConeOuterRadiusB_, photonHcalTowerConeInnerRadiusB_, photonHcalTowerThreshEB_, -1);
449  phoisolR2.hcalTowerSumEt = HcalTowerIsoB;
450 
452 
453  double HcalDepth1TowerIsoA = calculateHcalTowerIso(pho,
454  e,
455  es,
456  photonHcalDepth1TowerConeOuterRadiusA_,
457  photonHcalDepth1TowerConeInnerRadiusA_,
458  photonHcalDepth1TowerThreshEA_,
459  1);
460  phoisolR1.hcalDepth1TowerSumEt = HcalDepth1TowerIsoA;
461 
462  double HcalDepth1TowerIsoB = calculateHcalTowerIso(pho,
463  e,
464  es,
465  photonHcalDepth1TowerConeOuterRadiusB_,
466  photonHcalDepth1TowerConeInnerRadiusB_,
467  photonHcalDepth1TowerThreshEB_,
468  1);
469  phoisolR2.hcalDepth1TowerSumEt = HcalDepth1TowerIsoB;
470 
472 
473  double HcalDepth2TowerIsoA = calculateHcalTowerIso(pho,
474  e,
475  es,
476  photonHcalDepth2TowerConeOuterRadiusA_,
477  photonHcalDepth2TowerConeInnerRadiusA_,
478  photonHcalDepth2TowerThreshEA_,
479  2);
480  phoisolR1.hcalDepth2TowerSumEt = HcalDepth2TowerIsoA;
481 
482  double HcalDepth2TowerIsoB = calculateHcalTowerIso(pho,
483  e,
484  es,
485  photonHcalDepth2TowerConeOuterRadiusB_,
486  photonHcalDepth2TowerConeInnerRadiusB_,
487  photonHcalDepth2TowerThreshEB_,
488  2);
489  phoisolR2.hcalDepth2TowerSumEt = HcalDepth2TowerIsoB;
490 
491  // New Hcal isolation based on the new H/E definition (towers behind the BCs in the SC are used to evaluated H)
492  double HcalTowerBcIsoA =
493  calculateHcalTowerIso(pho, e, es, photonHcalTowerConeOuterRadiusA_, photonHcalTowerThreshEA_, -1);
494  phoisolR1.hcalTowerSumEtBc = HcalTowerBcIsoA;
495 
496  double HcalTowerBcIsoB =
497  calculateHcalTowerIso(pho, e, es, photonHcalTowerConeOuterRadiusB_, photonHcalTowerThreshEB_, -1);
498  phoisolR2.hcalTowerSumEtBc = HcalTowerBcIsoB;
499 
501 
502  double HcalDepth1TowerBcIsoA =
503  calculateHcalTowerIso(pho, e, es, photonHcalDepth1TowerConeOuterRadiusA_, photonHcalDepth1TowerThreshEA_, 1);
504  phoisolR1.hcalDepth1TowerSumEtBc = HcalDepth1TowerBcIsoA;
505 
506  double HcalDepth1TowerBcIsoB =
507  calculateHcalTowerIso(pho, e, es, photonHcalDepth1TowerConeOuterRadiusB_, photonHcalDepth1TowerThreshEB_, 1);
508  phoisolR2.hcalDepth1TowerSumEtBc = HcalDepth1TowerBcIsoB;
509 
511 
512  double HcalDepth2TowerBcIsoA =
513  calculateHcalTowerIso(pho, e, es, photonHcalDepth2TowerConeOuterRadiusA_, photonHcalDepth2TowerThreshEA_, 2);
514  phoisolR1.hcalDepth2TowerSumEtBc = HcalDepth2TowerBcIsoA;
515 
516  double HcalDepth2TowerBcIsoB =
517  calculateHcalTowerIso(pho, e, es, photonHcalDepth2TowerConeOuterRadiusB_, photonHcalDepth2TowerThreshEB_, 2);
518  phoisolR2.hcalDepth2TowerSumEtBc = HcalDepth2TowerBcIsoB;
519 }
520 
522  bool& isEBPho,
523  bool& isEEPho,
524  bool& isEBEtaGap,
525  bool& isEBPhiGap,
526  bool& isEERingGap,
527  bool& isEEDeeGap,
528  bool& isEBEEGap) {
529  const reco::CaloCluster& seedCluster = *(photon->superCluster()->seed());
530  // following line is temporary until the D. Evans or F. Ferri submit the new tag for ClusterAlgos
531  // DEfinitive will be
532  // DetId seedXtalId = scRef->seed()->seed();
533  DetId seedXtalId = seedCluster.hitsAndFractions()[0].first;
534  int detector = seedXtalId.subdetId();
535 
536  //Set fiducial flags for this photon.
537  double eta = photon->superCluster()->position().eta();
538  double feta = fabs(eta);
539 
540  if (detector == EcalBarrel) {
541  isEBPho = true;
542  if (EBDetId::isNextToEtaBoundary(EBDetId(seedXtalId))) {
543  if (fabs(feta - 1.479) < .1)
544  isEBEEGap = true;
545  else
546  isEBEtaGap = true;
547  }
548 
549  if (EBDetId::isNextToPhiBoundary(EBDetId(seedXtalId)))
550  isEBPhiGap = true;
551 
552  } else if (detector == EcalEndcap) {
553  isEEPho = true;
554  if (EEDetId::isNextToRingBoundary(EEDetId(seedXtalId))) {
555  if (fabs(feta - 1.479) < .1)
556  isEBEEGap = true;
557  else
558  isEERingGap = true;
559  }
560 
561  if (EEDetId::isNextToDBoundary(EEDetId(seedXtalId)))
562  isEEDeeGap = true;
563  }
564 }
565 
567  const edm::Event& e,
568  double& trkCone,
569  int& ntrkCone,
570  double pTThresh,
571  double RCone,
572  double RinnerCone,
573  double etaSlice,
574  double lip,
575  double d0) const {
576  ntrkCone = 0;
577  trkCone = 0;
578  //get the tracks
580  e.getByToken(trackInputTag_, tracks);
581  if (!tracks.isValid()) {
582  return;
583  }
585  //Photon Eta and Phi. Hope these are correct.
586  reco::BeamSpot vertexBeamSpot;
587  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
588  e.getByToken(beamSpotProducerTag_, recoBeamSpotHandle);
589  vertexBeamSpot = *recoBeamSpotHandle;
590 
591  PhotonTkIsolation phoIso(RCone,
592  RinnerCone,
593  etaSlice,
594  pTThresh,
595  lip,
596  d0,
597  trackCollection,
598  math::XYZPoint(vertexBeamSpot.x0(), vertexBeamSpot.y0(), vertexBeamSpot.z0()));
599 
600  std::pair<int, double> res = phoIso.getIso(photon);
601  ntrkCone = res.first;
602  trkCone = res.second;
603 }
604 
606  const edm::Event& iEvent,
607  const edm::EventSetup& iSetup,
608  double RCone,
609  double RConeInner,
610  double etaSlice,
611  double eMin,
612  double etMin,
613  bool vetoClusteredHits,
614  bool useNumXtals) const {
615  edm::Handle<EcalRecHitCollection> ecalhitsCollEB;
616  edm::Handle<EcalRecHitCollection> ecalhitsCollEE;
617  iEvent.getByToken(endcapecalCollection_, ecalhitsCollEE);
618 
619  iEvent.getByToken(barrelecalCollection_, ecalhitsCollEB);
620 
621  const EcalRecHitCollection* rechitsCollectionEE_ = ecalhitsCollEE.product();
622  const EcalRecHitCollection* rechitsCollectionEB_ = ecalhitsCollEB.product();
623 
625  iSetup.get<EcalSeverityLevelAlgoRcd>().get(sevlv);
626  const EcalSeverityLevelAlgo* sevLevel = sevlv.product();
627 
628  edm::ESHandle<CaloGeometry> geoHandle;
629  iSetup.get<CaloGeometryRecord>().get(geoHandle);
630 
631  EgammaRecHitIsolation phoIsoEB(
632  RCone, RConeInner, etaSlice, etMin, eMin, geoHandle, *rechitsCollectionEB_, sevLevel, DetId::Ecal);
633 
634  phoIsoEB.setVetoClustered(vetoClusteredHits);
635  phoIsoEB.setUseNumCrystals(useNumXtals);
636  phoIsoEB.doSeverityChecks(ecalhitsCollEB.product(), severityExclEB_);
637  phoIsoEB.doFlagChecks(flagsEB_);
638  double ecalIsolEB = phoIsoEB.getEtSum(photon);
639 
640  EgammaRecHitIsolation phoIsoEE(
641  RCone, RConeInner, etaSlice, etMin, eMin, geoHandle, *rechitsCollectionEE_, sevLevel, DetId::Ecal);
642 
643  phoIsoEE.setVetoClustered(vetoClusteredHits);
644  phoIsoEE.setUseNumCrystals(useNumXtals);
645  phoIsoEE.doSeverityChecks(ecalhitsCollEE.product(), severityExclEE_);
646  phoIsoEE.doFlagChecks(flagsEE_);
647 
648  double ecalIsolEE = phoIsoEE.getEtSum(photon);
649  // delete phoIso;
650  double ecalIsol = ecalIsolEB + ecalIsolEE;
651 
652  return ecalIsol;
653 }
654 
656  const edm::Event& iEvent,
657  const edm::EventSetup& iSetup,
658  double RCone,
659  double RConeInner,
660  double eMin,
661  signed int depth) const {
662  double hcalIsol = 0.;
663 
664  if (not hcalCollection_.isUninitialized()) {
665  edm::Handle<CaloTowerCollection> hcalhitsCollH;
666  iEvent.getByToken(hcalCollection_, hcalhitsCollH);
667 
668  const CaloTowerCollection* toww = hcalhitsCollH.product();
669 
670  //std::cout << "before iso call" << std::endl;
671  EgammaTowerIsolation phoIso(RCone, RConeInner, eMin, depth, toww);
672  hcalIsol = phoIso.getTowerEtSum(photon);
673  // delete phoIso;
674  //std::cout << "after call" << std::endl;
675  }
676 
677  return hcalIsol;
678 }
679 
681  const edm::Event& iEvent,
682  const edm::EventSetup& iSetup,
683  double RCone,
684  double eMin,
685  signed int depth) const {
686  double hcalIsol = 0.;
687 
688  if (not hcalCollection_.isUninitialized()) {
689  edm::Handle<CaloTowerCollection> hcalhitsCollH;
690  iEvent.getByToken(hcalCollection_, hcalhitsCollH);
691 
692  const CaloTowerCollection* toww = hcalhitsCollH.product();
693 
694  //std::cout << "before iso call" << std::endl;
695  EgammaTowerIsolation phoIso(RCone, 0., eMin, depth, toww);
696  hcalIsol = phoIso.getTowerEtSum(photon, &(photon->hcalTowersBehindClusters()));
697  // delete phoIso;
698  //std::cout << "after call" << std::endl;
699  }
700 
701  return hcalIsol;
702 }
T getParameter(std::string const &) const
double z0() const
z coordinate
Definition: BeamSpot.h:65
static bool isNextToEtaBoundary(EBDetId id)
Definition: EBDetId.cc:108
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
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:525
double getEtSum(const reco::Candidate *emObject) const
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:14
reco::SuperClusterRef superCluster() const override
Ref to SuperCluster.
Definition: Electron.h:6
void doFlagChecks(const std::vector< int > &v)
double calculateHcalTowerIso(const reco::Photon *photon, const edm::Event &iEvent, const edm::EventSetup &iSetup, double RCone, double RConeInner, double eMin, signed int depth) const
void setUseNumCrystals(bool b=true)
int iEvent
Definition: GenABIO.cc:224
const std::vector< std::pair< DetId, float > > & hitsAndFractions() const
Definition: CaloCluster.h:209
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)
static bool isNextToPhiBoundary(EBDetId id)
Definition: EBDetId.cc:113
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:48
static bool isNextToRingBoundary(EEDetId id)
Definition: EEDetId.cc:284
bool isUninitialized() const
Definition: EDGetToken.h:46
static bool isNextToDBoundary(EEDetId id)
Definition: EEDetId.cc:279
bool isValid() const
Definition: HandleBase.h:70
std::pair< int, float > getIso(const reco::Candidate *) const
void doSeverityChecks(const EcalRecHitCollection *const recHits, const std::vector< int > &v)
Definition: DetId.h:17
T const * product() const
Definition: Handle.h:69
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
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
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
std::vector< double > moduleEtaBoundary_
void setVetoClustered(bool b=true)
T get() const
Definition: EventSetup.h:73
double getTowerEtSum(const reco::Candidate *cand, const std::vector< CaloTowerDetId > *detIdToExclude=0) const
double y0() const
y coordinate
Definition: BeamSpot.h:63
const std::vector< CaloTowerDetId > & hcalTowersBehindClusters() const
Definition: Photon.h:224
static void classify(const reco::Photon *photon, bool &isEBPho, bool &isEEPho, bool &isEBEtaGap, bool &isEBPhiGap, bool &isEERingGap, bool &isEEDeeGap, bool &isEBEEGap)
T const * product() const
Definition: ESHandle.h:86
double x0() const
x coordinate
Definition: BeamSpot.h:61