Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-09-10 10:03:07

0001 /** \class PhotonIsolationCalculator
0002  *  Determine and Set quality information on Photon Objects
0003  *
0004  *  \author A. Askew, N. Marinelli, M.B. Anderson
0005  */
0006 
0007 #include "RecoEgamma/PhotonIdentification/interface/PhotonIsolationCalculator.h"
0008 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0009 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
0010 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0011 #include "DataFormats/TrackReco/interface/Track.h"
0012 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0013 #include "DataFormats/DetId/interface/DetId.h"
0014 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0015 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0016 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0017 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0018 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0019 
0020 #include "RecoEgamma/EgammaIsolationAlgos/interface/PhotonTkIsolation.h"
0021 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaEcalIsolation.h"
0022 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaRecHitIsolation.h"
0023 
0024 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0025 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0026 
0027 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0028 
0029 void PhotonIsolationCalculator::setup(const edm::ParameterSet& conf,
0030                                       std::vector<int> const& flagsEB,
0031                                       std::vector<int> const& flagsEE,
0032                                       std::vector<int> const& severitiesEB,
0033                                       std::vector<int> const& severitiesEE,
0034                                       edm::ConsumesCollector&& iC) {
0035   trackInputTag_ = iC.consumes<reco::TrackCollection>(conf.getParameter<edm::InputTag>("trackProducer"));
0036   beamSpotProducerTag_ = iC.consumes<reco::BeamSpot>(conf.getParameter<edm::InputTag>("beamSpotProducer"));
0037   barrelecalCollection_ =
0038       iC.consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("barrelEcalRecHitCollection"));
0039   endcapecalCollection_ =
0040       iC.consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("endcapEcalRecHitCollection"));
0041 
0042   auto hbhetag = conf.getParameter<edm::InputTag>("HBHERecHitCollection");
0043   if (not hbhetag.label().empty())
0044     hbheRecHitsTag_ = iC.consumes<HBHERecHitCollection>(hbhetag);
0045 
0046   caloGeometryToken_ = decltype(caloGeometryToken_){iC.esConsumes()};
0047   hcalTopologyToken_ = decltype(hcalTopologyToken_){iC.esConsumes()};
0048   hcalChannelQualityToken_ = decltype(hcalChannelQualityToken_){iC.esConsumes(edm::ESInputTag("", "withTopo"))};
0049   hcalSevLvlComputerToken_ = decltype(hcalSevLvlComputerToken_){iC.esConsumes()};
0050   towerMapToken_ = decltype(towerMapToken_){iC.esConsumes()};
0051   ecalSevLvlToken_ = iC.esConsumes();
0052 
0053   //  gsfRecoInputTag_ = conf.getParameter<edm::InputTag>("GsfRecoCollection");
0054   modulePhiBoundary_ = conf.getParameter<double>("modulePhiBoundary");
0055   moduleEtaBoundary_ = conf.getParameter<std::vector<double>>("moduleEtaBoundary");
0056   //
0057   vetoClusteredEcalHits_ = conf.getParameter<bool>("vetoClustered");
0058   useNumCrystals_ = conf.getParameter<bool>("useNumCrystals");
0059 
0060   /// Isolation parameters for barrel and for two different cone sizes
0061   int i = 0;
0062   trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("TrackConeOuterRadiusA_Barrel"));
0063   trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("TrackConeInnerRadiusA_Barrel"));
0064   trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("isolationtrackThresholdA_Barrel"));
0065   trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("longImpactParameterA_Barrel"));
0066   trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("transImpactParameterA_Barrel"));
0067   trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("isolationtrackEtaSliceA_Barrel"));
0068 
0069   i = 0;
0070   ecalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitInnerRadiusA_Barrel"));
0071   ecalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitOuterRadiusA_Barrel"));
0072   ecalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitEtaSliceA_Barrel"));
0073   ecalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitThreshEA_Barrel"));
0074   ecalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitThreshEtA_Barrel"));
0075 
0076   hcalIsoInnerRadAEB_ = conf.getParameter<std::array<double, 7>>("HcalRecHitInnerRadiusA_Barrel");
0077   hcalIsoOuterRadAEB_ = conf.getParameter<std::array<double, 7>>("HcalRecHitOuterRadiusA_Barrel");
0078 
0079   i = 0;
0080   trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("TrackConeOuterRadiusB_Barrel"));
0081   trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("TrackConeInnerRadiusB_Barrel"));
0082   trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("isolationtrackThresholdB_Barrel"));
0083   trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("longImpactParameterB_Barrel"));
0084   trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("transImpactParameterB_Barrel"));
0085   trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("isolationtrackEtaSliceB_Barrel"));
0086 
0087   i = 0;
0088   ecalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitInnerRadiusB_Barrel"));
0089   ecalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitOuterRadiusB_Barrel"));
0090   ecalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitEtaSliceB_Barrel"));
0091   ecalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitThreshEB_Barrel"));
0092   ecalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitThreshEtB_Barrel"));
0093 
0094   hcalIsoInnerRadBEB_ = conf.getParameter<std::array<double, 7>>("HcalRecHitInnerRadiusB_Barrel");
0095   hcalIsoOuterRadBEB_ = conf.getParameter<std::array<double, 7>>("HcalRecHitOuterRadiusB_Barrel");
0096 
0097   /// Isolation parameters for Endcap and for two different cone sizes
0098   i = 0;
0099   trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("TrackConeOuterRadiusA_Endcap"));
0100   trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("TrackConeInnerRadiusA_Endcap"));
0101   trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("isolationtrackThresholdA_Endcap"));
0102   trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("longImpactParameterA_Endcap"));
0103   trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("transImpactParameterA_Endcap"));
0104   trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("isolationtrackEtaSliceA_Endcap"));
0105 
0106   i = 0;
0107   ecalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitInnerRadiusA_Endcap"));
0108   ecalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitOuterRadiusA_Endcap"));
0109   ecalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitEtaSliceA_Endcap"));
0110   ecalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitThreshEA_Endcap"));
0111   ecalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitThreshEtA_Endcap"));
0112 
0113   hcalIsoInnerRadAEE_ = conf.getParameter<std::array<double, 7>>("HcalRecHitInnerRadiusA_Endcap");
0114   hcalIsoOuterRadAEE_ = conf.getParameter<std::array<double, 7>>("HcalRecHitOuterRadiusA_Endcap");
0115 
0116   i = 0;
0117   trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("TrackConeOuterRadiusB_Endcap"));
0118   trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("TrackConeInnerRadiusB_Endcap"));
0119   trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("isolationtrackThresholdB_Endcap"));
0120   trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("longImpactParameterB_Endcap"));
0121   trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("transImpactParameterB_Endcap"));
0122   trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("isolationtrackEtaSliceB_Endcap"));
0123 
0124   i = 0;
0125   ecalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitInnerRadiusB_Endcap"));
0126   ecalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitOuterRadiusB_Endcap"));
0127   ecalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitEtaSliceB_Endcap"));
0128   ecalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitThreshEB_Endcap"));
0129   ecalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitThreshEtB_Endcap"));
0130 
0131   hcalIsoInnerRadBEE_ = conf.getParameter<std::array<double, 7>>("HcalRecHitInnerRadiusB_Endcap");
0132   hcalIsoOuterRadBEE_ = conf.getParameter<std::array<double, 7>>("HcalRecHitOuterRadiusB_Endcap");
0133 
0134   //Pick up the variables for the spike removal
0135   flagsEB_ = flagsEB;
0136   flagsEE_ = flagsEE;
0137   severityExclEB_ = severitiesEB;
0138   severityExclEE_ = severitiesEE;
0139 
0140   hcalIsoEThresHB_ = conf.getParameter<EgammaHcalIsolation::arrayHB>("recHitEThresholdHB");
0141   hcalIsoEThresHE_ = conf.getParameter<EgammaHcalIsolation::arrayHE>("recHitEThresholdHE");
0142   maxHcalSeverity_ = conf.getParameter<int>("maxHcalRecHitSeverity");
0143 }
0144 
0145 void PhotonIsolationCalculator::calculate(const reco::Photon* pho,
0146                                           const edm::Event& e,
0147                                           const edm::EventSetup& es,
0148                                           reco::Photon::FiducialFlags& phofid,
0149                                           reco::Photon::IsolationVariables& phoisolR1,
0150                                           reco::Photon::IsolationVariables& phoisolR2) const {
0151   //Get fiducial flags. This does not really belong here
0152   bool isEBPho = false;
0153   bool isEEPho = false;
0154   bool isEBEtaGap = false;
0155   bool isEBPhiGap = false;
0156   bool isEERingGap = false;
0157   bool isEEDeeGap = false;
0158   bool isEBEEGap = false;
0159   classify(pho, isEBPho, isEEPho, isEBEtaGap, isEBPhiGap, isEERingGap, isEEDeeGap, isEBEEGap);
0160   phofid.isEB = isEBPho;
0161   phofid.isEE = isEEPho;
0162   phofid.isEBEtaGap = isEBEtaGap;
0163   phofid.isEBPhiGap = isEBPhiGap;
0164   phofid.isEERingGap = isEERingGap;
0165   phofid.isEEDeeGap = isEEDeeGap;
0166   phofid.isEBEEGap = isEBEEGap;
0167 
0168   auto const& caloGeometry = es.getData(caloGeometryToken_);
0169   auto const& hcalTopology = &es.getData(hcalTopologyToken_);
0170   auto const& hcalChannelQuality = &es.getData(hcalChannelQualityToken_);
0171   auto const& hcalSevLvlComputer = &es.getData(hcalSevLvlComputerToken_);
0172   auto const& towerMap = es.getData(towerMapToken_);
0173 
0174   // Calculate isolation variables. cone sizes and thresholds
0175   // are set for Barrel and endcap separately
0176 
0177   reco::SuperClusterRef scRef = pho->superCluster();
0178   const reco::BasicCluster& seedCluster = *(scRef->seed());
0179   DetId seedXtalId = seedCluster.hitsAndFractions()[0].first;
0180   int detector = seedXtalId.subdetId();
0181 
0182   //Isolation parameters variables
0183   double photonEcalRecHitConeInnerRadiusA_;
0184   double photonEcalRecHitConeOuterRadiusA_;
0185   double photonEcalRecHitEtaSliceA_;
0186   double photonEcalRecHitThreshEA_;
0187   double photonEcalRecHitThreshEtA_;
0188   double trackConeOuterRadiusA_;
0189   double trackConeInnerRadiusA_;
0190   double isolationtrackThresholdA_;
0191   double isolationtrackEtaSliceA_;
0192   double trackLipRadiusA_;
0193   double trackD0RadiusA_;
0194   double photonEcalRecHitConeInnerRadiusB_;
0195   double photonEcalRecHitConeOuterRadiusB_;
0196   double photonEcalRecHitEtaSliceB_;
0197   double photonEcalRecHitThreshEB_;
0198   double photonEcalRecHitThreshEtB_;
0199   double trackConeOuterRadiusB_;
0200   double trackConeInnerRadiusB_;
0201   double isolationtrackThresholdB_;
0202   double isolationtrackEtaSliceB_;
0203   double trackLipRadiusB_;
0204   double trackD0RadiusB_;
0205 
0206   if (detector == EcalBarrel) {
0207     trackConeOuterRadiusA_ = trkIsoBarrelRadiusA_[0];
0208     trackConeInnerRadiusA_ = trkIsoBarrelRadiusA_[1];
0209     isolationtrackThresholdA_ = trkIsoBarrelRadiusA_[2];
0210     trackLipRadiusA_ = trkIsoBarrelRadiusA_[3];
0211     trackD0RadiusA_ = trkIsoBarrelRadiusA_[4];
0212     isolationtrackEtaSliceA_ = trkIsoBarrelRadiusA_[5];
0213 
0214     photonEcalRecHitConeInnerRadiusA_ = ecalIsoBarrelRadiusA_[0];
0215     photonEcalRecHitConeOuterRadiusA_ = ecalIsoBarrelRadiusA_[1];
0216     photonEcalRecHitEtaSliceA_ = ecalIsoBarrelRadiusA_[2];
0217     photonEcalRecHitThreshEA_ = ecalIsoBarrelRadiusA_[3];
0218     photonEcalRecHitThreshEtA_ = ecalIsoBarrelRadiusA_[4];
0219 
0220     trackConeOuterRadiusB_ = trkIsoBarrelRadiusB_[0];
0221     trackConeInnerRadiusB_ = trkIsoBarrelRadiusB_[1];
0222     isolationtrackThresholdB_ = trkIsoBarrelRadiusB_[2];
0223     trackLipRadiusB_ = trkIsoBarrelRadiusB_[3];
0224     trackD0RadiusB_ = trkIsoBarrelRadiusB_[4];
0225     isolationtrackEtaSliceB_ = trkIsoBarrelRadiusB_[5];
0226 
0227     photonEcalRecHitConeInnerRadiusB_ = ecalIsoBarrelRadiusB_[0];
0228     photonEcalRecHitConeOuterRadiusB_ = ecalIsoBarrelRadiusB_[1];
0229     photonEcalRecHitEtaSliceB_ = ecalIsoBarrelRadiusB_[2];
0230     photonEcalRecHitThreshEB_ = ecalIsoBarrelRadiusB_[3];
0231     photonEcalRecHitThreshEtB_ = ecalIsoBarrelRadiusB_[4];
0232   } else {
0233     // detector==EcalEndcap
0234 
0235     trackConeOuterRadiusA_ = trkIsoEndcapRadiusA_[0];
0236     trackConeInnerRadiusA_ = trkIsoEndcapRadiusA_[1];
0237     isolationtrackThresholdA_ = trkIsoEndcapRadiusA_[2];
0238     trackLipRadiusA_ = trkIsoEndcapRadiusA_[3];
0239     trackD0RadiusA_ = trkIsoEndcapRadiusA_[4];
0240     isolationtrackEtaSliceA_ = trkIsoEndcapRadiusA_[5];
0241 
0242     photonEcalRecHitConeInnerRadiusA_ = ecalIsoEndcapRadiusA_[0];
0243     photonEcalRecHitConeOuterRadiusA_ = ecalIsoEndcapRadiusA_[1];
0244     photonEcalRecHitEtaSliceA_ = ecalIsoEndcapRadiusA_[2];
0245     photonEcalRecHitThreshEA_ = ecalIsoEndcapRadiusA_[3];
0246     photonEcalRecHitThreshEtA_ = ecalIsoEndcapRadiusA_[4];
0247 
0248     trackConeOuterRadiusB_ = trkIsoEndcapRadiusB_[0];
0249     trackConeInnerRadiusB_ = trkIsoEndcapRadiusB_[1];
0250     isolationtrackThresholdB_ = trkIsoEndcapRadiusB_[2];
0251     trackLipRadiusB_ = trkIsoEndcapRadiusB_[3];
0252     trackD0RadiusB_ = trkIsoEndcapRadiusB_[4];
0253     isolationtrackEtaSliceB_ = trkIsoEndcapRadiusB_[5];
0254 
0255     photonEcalRecHitConeInnerRadiusB_ = ecalIsoEndcapRadiusB_[0];
0256     photonEcalRecHitConeOuterRadiusB_ = ecalIsoEndcapRadiusB_[1];
0257     photonEcalRecHitEtaSliceB_ = ecalIsoEndcapRadiusB_[2];
0258     photonEcalRecHitThreshEB_ = ecalIsoEndcapRadiusB_[3];
0259     photonEcalRecHitThreshEtB_ = ecalIsoEndcapRadiusB_[4];
0260   }
0261 
0262   //Calculate hollow cone track isolation, CONE A
0263   int ntrkA = 0;
0264   double trkisoA = 0;
0265   calculateTrackIso(pho,
0266                     e,
0267                     trkisoA,
0268                     ntrkA,
0269                     isolationtrackThresholdA_,
0270                     trackConeOuterRadiusA_,
0271                     trackConeInnerRadiusA_,
0272                     isolationtrackEtaSliceA_,
0273                     trackLipRadiusA_,
0274                     trackD0RadiusA_);
0275 
0276   //Calculate solid cone track isolation, CONE A
0277   int sntrkA = 0;
0278   double strkisoA = 0;
0279   calculateTrackIso(pho,
0280                     e,
0281                     strkisoA,
0282                     sntrkA,
0283                     isolationtrackThresholdA_,
0284                     trackConeOuterRadiusA_,
0285                     0.,
0286                     isolationtrackEtaSliceA_,
0287                     trackLipRadiusA_,
0288                     trackD0RadiusA_);
0289 
0290   phoisolR1.nTrkHollowCone = ntrkA;
0291   phoisolR1.trkSumPtHollowCone = trkisoA;
0292   phoisolR1.nTrkSolidCone = sntrkA;
0293   phoisolR1.trkSumPtSolidCone = strkisoA;
0294 
0295   //Calculate hollow cone track isolation, CONE B
0296   int ntrkB = 0;
0297   double trkisoB = 0;
0298   calculateTrackIso(pho,
0299                     e,
0300                     trkisoB,
0301                     ntrkB,
0302                     isolationtrackThresholdB_,
0303                     trackConeOuterRadiusB_,
0304                     trackConeInnerRadiusB_,
0305                     isolationtrackEtaSliceB_,
0306                     trackLipRadiusB_,
0307                     trackD0RadiusB_);
0308 
0309   //Calculate solid cone track isolation, CONE B
0310   int sntrkB = 0;
0311   double strkisoB = 0;
0312   calculateTrackIso(pho,
0313                     e,
0314                     strkisoB,
0315                     sntrkB,
0316                     isolationtrackThresholdB_,
0317                     trackConeOuterRadiusB_,
0318                     0.,
0319                     isolationtrackEtaSliceB_,
0320                     trackLipRadiusB_,
0321                     trackD0RadiusB_);
0322 
0323   phoisolR2.nTrkHollowCone = ntrkB;
0324   phoisolR2.trkSumPtHollowCone = trkisoB;
0325   phoisolR2.nTrkSolidCone = sntrkB;
0326   phoisolR2.trkSumPtSolidCone = strkisoB;
0327 
0328   //   std::cout << "Output from solid cone track isolation: ";
0329   //   std::cout << " Sum pT: " << strkiso << " ntrk: " << sntrk << std::endl;
0330 
0331   double EcalRecHitIsoA = calculateEcalRecHitIso(pho,
0332                                                  e,
0333                                                  es,
0334                                                  photonEcalRecHitConeOuterRadiusA_,
0335                                                  photonEcalRecHitConeInnerRadiusA_,
0336                                                  photonEcalRecHitEtaSliceA_,
0337                                                  photonEcalRecHitThreshEA_,
0338                                                  photonEcalRecHitThreshEtA_,
0339                                                  vetoClusteredEcalHits_,
0340                                                  useNumCrystals_);
0341   phoisolR1.ecalRecHitSumEt = EcalRecHitIsoA;
0342 
0343   double EcalRecHitIsoB = calculateEcalRecHitIso(pho,
0344                                                  e,
0345                                                  es,
0346                                                  photonEcalRecHitConeOuterRadiusB_,
0347                                                  photonEcalRecHitConeInnerRadiusB_,
0348                                                  photonEcalRecHitEtaSliceB_,
0349                                                  photonEcalRecHitThreshEB_,
0350                                                  photonEcalRecHitThreshEtB_,
0351                                                  vetoClusteredEcalHits_,
0352                                                  useNumCrystals_);
0353   phoisolR2.ecalRecHitSumEt = EcalRecHitIsoB;
0354 
0355   if (not hbheRecHitsTag_.isUninitialized()) {
0356     auto const& hbheRecHits = e.get(hbheRecHitsTag_);
0357 
0358     auto fcone = [this,
0359                   pho,
0360                   &caloGeometry,
0361                   &hcalTopo = *hcalTopology,
0362                   &hcalQual = *hcalChannelQuality,
0363                   &hcalSev = *hcalSevLvlComputer,
0364                   &towerMap,
0365                   &hbheRecHits](double outer, double inner, int depth) {
0366       return calculateHcalRecHitIso<false>(
0367           pho, caloGeometry, hcalTopo, hcalQual, hcalSev, towerMap, hbheRecHits, outer, inner, depth);
0368     };
0369 
0370     auto fbc = [this,
0371                 pho,
0372                 &caloGeometry,
0373                 &hcalTopo = *hcalTopology,
0374                 &hcalQual = *hcalChannelQuality,
0375                 &hcalSev = *hcalSevLvlComputer,
0376                 &towerMap,
0377                 &hbheRecHits](double outer, int depth) {
0378       return calculateHcalRecHitIso<true>(
0379           pho, caloGeometry, hcalTopo, hcalQual, hcalSev, towerMap, hbheRecHits, outer, 0., depth);
0380     };
0381 
0382     for (size_t id = 0; id < phoisolR1.hcalRecHitSumEt.size(); ++id) {
0383       const auto& outerA = detector == EcalBarrel ? hcalIsoOuterRadAEB_[id] : hcalIsoOuterRadAEE_[id];
0384       const auto& outerB = detector == EcalBarrel ? hcalIsoOuterRadBEB_[id] : hcalIsoOuterRadBEE_[id];
0385       const auto& innerA = detector == EcalBarrel ? hcalIsoInnerRadAEB_[id] : hcalIsoInnerRadAEE_[id];
0386       const auto& innerB = detector == EcalBarrel ? hcalIsoInnerRadBEB_[id] : hcalIsoInnerRadBEE_[id];
0387 
0388       phoisolR1.hcalRecHitSumEt[id] = fcone(outerA, innerA, id + 1);
0389       phoisolR2.hcalRecHitSumEt[id] = fcone(outerB, innerB, id + 1);
0390 
0391       phoisolR1.hcalRecHitSumEtBc[id] = fbc(outerA, id + 1);
0392       phoisolR2.hcalRecHitSumEtBc[id] = fbc(outerB, id + 1);
0393     }
0394   }
0395 
0396   phoisolR1.pre7DepthHcal = false;
0397   phoisolR2.pre7DepthHcal = false;
0398 }
0399 
0400 void PhotonIsolationCalculator::classify(const reco::Photon* photon,
0401                                          bool& isEBPho,
0402                                          bool& isEEPho,
0403                                          bool& isEBEtaGap,
0404                                          bool& isEBPhiGap,
0405                                          bool& isEERingGap,
0406                                          bool& isEEDeeGap,
0407                                          bool& isEBEEGap) {
0408   const reco::CaloCluster& seedCluster = *(photon->superCluster()->seed());
0409   // following line is temporary until the D. Evans or F. Ferri submit the new tag for ClusterAlgos
0410   // DEfinitive will be
0411   // DetId seedXtalId = scRef->seed()->seed();
0412   DetId seedXtalId = seedCluster.hitsAndFractions()[0].first;
0413   int detector = seedXtalId.subdetId();
0414 
0415   //Set fiducial flags for this photon.
0416   double eta = photon->superCluster()->position().eta();
0417   double feta = fabs(eta);
0418 
0419   if (detector == EcalBarrel) {
0420     isEBPho = true;
0421     if (EBDetId::isNextToEtaBoundary(EBDetId(seedXtalId))) {
0422       if (fabs(feta - 1.479) < .1)
0423         isEBEEGap = true;
0424       else
0425         isEBEtaGap = true;
0426     }
0427 
0428     if (EBDetId::isNextToPhiBoundary(EBDetId(seedXtalId)))
0429       isEBPhiGap = true;
0430 
0431   } else if (detector == EcalEndcap) {
0432     isEEPho = true;
0433     if (EEDetId::isNextToRingBoundary(EEDetId(seedXtalId))) {
0434       if (fabs(feta - 1.479) < .1)
0435         isEBEEGap = true;
0436       else
0437         isEERingGap = true;
0438     }
0439 
0440     if (EEDetId::isNextToDBoundary(EEDetId(seedXtalId)))
0441       isEEDeeGap = true;
0442   }
0443 }
0444 
0445 void PhotonIsolationCalculator::calculateTrackIso(const reco::Photon* photon,
0446                                                   const edm::Event& e,
0447                                                   double& trkCone,
0448                                                   int& ntrkCone,
0449                                                   double pTThresh,
0450                                                   double RCone,
0451                                                   double RinnerCone,
0452                                                   double etaSlice,
0453                                                   double lip,
0454                                                   double d0) const {
0455   ntrkCone = 0;
0456   trkCone = 0;
0457   //get the tracks
0458   edm::Handle<reco::TrackCollection> tracks;
0459   e.getByToken(trackInputTag_, tracks);
0460   if (!tracks.isValid()) {
0461     return;
0462   }
0463   const reco::TrackCollection* trackCollection = tracks.product();
0464   //Photon Eta and Phi.  Hope these are correct.
0465   reco::BeamSpot vertexBeamSpot;
0466   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0467   e.getByToken(beamSpotProducerTag_, recoBeamSpotHandle);
0468   vertexBeamSpot = *recoBeamSpotHandle;
0469 
0470   PhotonTkIsolation phoIso(RCone,
0471                            RinnerCone,
0472                            etaSlice,
0473                            pTThresh,
0474                            lip,
0475                            d0,
0476                            trackCollection,
0477                            math::XYZPoint(vertexBeamSpot.x0(), vertexBeamSpot.y0(), vertexBeamSpot.z0()));
0478 
0479   std::pair<int, double> res = phoIso.getIso(photon);
0480   ntrkCone = res.first;
0481   trkCone = res.second;
0482 }
0483 
0484 double PhotonIsolationCalculator::calculateEcalRecHitIso(const reco::Photon* photon,
0485                                                          const edm::Event& iEvent,
0486                                                          const edm::EventSetup& iSetup,
0487                                                          double RCone,
0488                                                          double RConeInner,
0489                                                          double etaSlice,
0490                                                          double eMin,
0491                                                          double etMin,
0492                                                          bool vetoClusteredHits,
0493                                                          bool useNumXtals) const {
0494   edm::Handle<EcalRecHitCollection> ecalhitsCollEB;
0495   edm::Handle<EcalRecHitCollection> ecalhitsCollEE;
0496   iEvent.getByToken(endcapecalCollection_, ecalhitsCollEE);
0497 
0498   iEvent.getByToken(barrelecalCollection_, ecalhitsCollEB);
0499 
0500   const EcalRecHitCollection* rechitsCollectionEE_ = ecalhitsCollEE.product();
0501   const EcalRecHitCollection* rechitsCollectionEB_ = ecalhitsCollEB.product();
0502 
0503   const EcalSeverityLevelAlgo* sevLevel = &iSetup.getData(ecalSevLvlToken_);
0504 
0505   edm::ESHandle<CaloGeometry> geoHandle = iSetup.getHandle(caloGeometryToken_);
0506 
0507   EgammaRecHitIsolation phoIsoEB(
0508       RCone, RConeInner, etaSlice, etMin, eMin, geoHandle, *rechitsCollectionEB_, sevLevel, DetId::Ecal);
0509 
0510   phoIsoEB.setVetoClustered(vetoClusteredHits);
0511   phoIsoEB.setUseNumCrystals(useNumXtals);
0512   phoIsoEB.doSeverityChecks(ecalhitsCollEB.product(), severityExclEB_);
0513   phoIsoEB.doFlagChecks(flagsEB_);
0514   double ecalIsolEB = phoIsoEB.getEtSum(photon);
0515 
0516   EgammaRecHitIsolation phoIsoEE(
0517       RCone, RConeInner, etaSlice, etMin, eMin, geoHandle, *rechitsCollectionEE_, sevLevel, DetId::Ecal);
0518 
0519   phoIsoEE.setVetoClustered(vetoClusteredHits);
0520   phoIsoEE.setUseNumCrystals(useNumXtals);
0521   phoIsoEE.doSeverityChecks(ecalhitsCollEE.product(), severityExclEE_);
0522   phoIsoEE.doFlagChecks(flagsEE_);
0523 
0524   double ecalIsolEE = phoIsoEE.getEtSum(photon);
0525   //  delete phoIso;
0526   double ecalIsol = ecalIsolEB + ecalIsolEE;
0527 
0528   return ecalIsol;
0529 }
0530 
0531 template <bool isoBC>
0532 double PhotonIsolationCalculator::calculateHcalRecHitIso(const reco::Photon* photon,
0533                                                          const CaloGeometry& geometry,
0534                                                          const HcalTopology& hcalTopology,
0535                                                          const HcalChannelQuality& hcalChStatus,
0536                                                          const HcalSeverityLevelComputer& hcalSevLvlComputer,
0537                                                          const CaloTowerConstituentsMap& towerMap,
0538                                                          const HBHERecHitCollection& hbheRecHits,
0539                                                          double RCone,
0540                                                          double RConeInner,
0541                                                          int depth) const {
0542   const EgammaHcalIsolation::arrayHB e04{{0., 0., 0., 0.}};
0543   const EgammaHcalIsolation::arrayHE e07{{0., 0., 0., 0., 0., 0., 0.}};
0544 
0545   if constexpr (isoBC) {
0546     auto hcaliso = EgammaHcalIsolation(EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
0547                                        RCone,
0548                                        EgammaHcalIsolation::InclusionRule::isBehindClusterSeed,
0549                                        RConeInner,
0550                                        hcalIsoEThresHB_,
0551                                        e04,
0552                                        maxHcalSeverity_,
0553                                        hcalIsoEThresHE_,
0554                                        e07,
0555                                        maxHcalSeverity_,
0556                                        hbheRecHits,
0557                                        geometry,
0558                                        hcalTopology,
0559                                        hcalChStatus,
0560                                        hcalSevLvlComputer,
0561                                        towerMap);
0562 
0563     return hcaliso.getHcalEtSumBc(photon, depth);
0564   } else {
0565     auto hcaliso = EgammaHcalIsolation(EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
0566                                        RCone,
0567                                        EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
0568                                        RConeInner,
0569                                        hcalIsoEThresHB_,
0570                                        e04,
0571                                        maxHcalSeverity_,
0572                                        hcalIsoEThresHE_,
0573                                        e07,
0574                                        maxHcalSeverity_,
0575                                        hbheRecHits,
0576                                        geometry,
0577                                        hcalTopology,
0578                                        hcalChStatus,
0579                                        hcalSevLvlComputer,
0580                                        towerMap);
0581 
0582     return hcaliso.getHcalEtSum(photon, depth);
0583   }
0584 }