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