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