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