Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-10-08 23:09:55 

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 PhotonIsolationCalculator::PhotonIsolationCalculator(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   ecalPFRechitThresholdsToken_ = iC.esConsumes();
0053   modulePhiBoundary_ = conf.getParameter<double>("modulePhiBoundary");
0054   moduleEtaBoundary_ = conf.getParameter<std::vector<double>>("moduleEtaBoundary");
0055   //
0056   vetoClusteredEcalHits_ = conf.getParameter<bool>("vetoClustered");
0057   useNumCrystals_ = conf.getParameter<bool>("useNumCrystals");
0058 
0059   /// Isolation parameters for barrel and for two different cone sizes
0060   int i = 0;
0061   trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("TrackConeOuterRadiusA_Barrel"));
0062   trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("TrackConeInnerRadiusA_Barrel"));
0063   trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("isolationtrackThresholdA_Barrel"));
0064   trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("longImpactParameterA_Barrel"));
0065   trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("transImpactParameterA_Barrel"));
0066   trkIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("isolationtrackEtaSliceA_Barrel"));
0067 
0068   i = 0;
0069   ecalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitInnerRadiusA_Barrel"));
0070   ecalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitOuterRadiusA_Barrel"));
0071   ecalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitEtaSliceA_Barrel"));
0072   ecalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitThreshEA_Barrel"));
0073   ecalIsoBarrelRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitThreshEtA_Barrel"));
0074 
0075   hcalIsoInnerRadAEB_ = conf.getParameter<std::array<double, 7>>("HcalRecHitInnerRadiusA_Barrel");
0076   hcalIsoOuterRadAEB_ = conf.getParameter<std::array<double, 7>>("HcalRecHitOuterRadiusA_Barrel");
0077 
0078   i = 0;
0079   trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("TrackConeOuterRadiusB_Barrel"));
0080   trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("TrackConeInnerRadiusB_Barrel"));
0081   trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("isolationtrackThresholdB_Barrel"));
0082   trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("longImpactParameterB_Barrel"));
0083   trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("transImpactParameterB_Barrel"));
0084   trkIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("isolationtrackEtaSliceB_Barrel"));
0085 
0086   i = 0;
0087   ecalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitInnerRadiusB_Barrel"));
0088   ecalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitOuterRadiusB_Barrel"));
0089   ecalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitEtaSliceB_Barrel"));
0090   ecalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitThreshEB_Barrel"));
0091   ecalIsoBarrelRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitThreshEtB_Barrel"));
0092 
0093   hcalIsoInnerRadBEB_ = conf.getParameter<std::array<double, 7>>("HcalRecHitInnerRadiusB_Barrel");
0094   hcalIsoOuterRadBEB_ = conf.getParameter<std::array<double, 7>>("HcalRecHitOuterRadiusB_Barrel");
0095 
0096   /// Isolation parameters for Endcap and for two different cone sizes
0097   i = 0;
0098   trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("TrackConeOuterRadiusA_Endcap"));
0099   trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("TrackConeInnerRadiusA_Endcap"));
0100   trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("isolationtrackThresholdA_Endcap"));
0101   trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("longImpactParameterA_Endcap"));
0102   trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("transImpactParameterA_Endcap"));
0103   trkIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("isolationtrackEtaSliceA_Endcap"));
0104 
0105   i = 0;
0106   ecalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitInnerRadiusA_Endcap"));
0107   ecalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitOuterRadiusA_Endcap"));
0108   ecalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitEtaSliceA_Endcap"));
0109   ecalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitThreshEA_Endcap"));
0110   ecalIsoEndcapRadiusA_[i++] = (conf.getParameter<double>("EcalRecHitThreshEtA_Endcap"));
0111 
0112   hcalIsoInnerRadAEE_ = conf.getParameter<std::array<double, 7>>("HcalRecHitInnerRadiusA_Endcap");
0113   hcalIsoOuterRadAEE_ = conf.getParameter<std::array<double, 7>>("HcalRecHitOuterRadiusA_Endcap");
0114 
0115   i = 0;
0116   trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("TrackConeOuterRadiusB_Endcap"));
0117   trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("TrackConeInnerRadiusB_Endcap"));
0118   trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("isolationtrackThresholdB_Endcap"));
0119   trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("longImpactParameterB_Endcap"));
0120   trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("transImpactParameterB_Endcap"));
0121   trkIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("isolationtrackEtaSliceB_Endcap"));
0122 
0123   i = 0;
0124   ecalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitInnerRadiusB_Endcap"));
0125   ecalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitOuterRadiusB_Endcap"));
0126   ecalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitEtaSliceB_Endcap"));
0127   ecalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitThreshEB_Endcap"));
0128   ecalIsoEndcapRadiusB_[i++] = (conf.getParameter<double>("EcalRecHitThreshEtB_Endcap"));
0129 
0130   hcalIsoInnerRadBEE_ = conf.getParameter<std::array<double, 7>>("HcalRecHitInnerRadiusB_Endcap");
0131   hcalIsoOuterRadBEE_ = conf.getParameter<std::array<double, 7>>("HcalRecHitOuterRadiusB_Endcap");
0132 
0133   //Pick up the variables for the spike removal
0134   flagsEB_ = flagsEB;
0135   flagsEE_ = flagsEE;
0136   severityExclEB_ = severitiesEB;
0137   severityExclEE_ = severitiesEE;
0138 
0139   hcalIsoEThresHB_ = conf.getParameter<EgammaHcalIsolation::arrayHB>("recHitEThresholdHB");
0140   hcalIsoEThresHE_ = conf.getParameter<EgammaHcalIsolation::arrayHE>("recHitEThresholdHE");
0141   maxHcalSeverity_ = conf.getParameter<int>("maxHcalRecHitSeverity");
0142 }
0143 
0144 void PhotonIsolationCalculator::calculate(const reco::Photon* pho,
0145                                           const edm::Event& e,
0146                                           const edm::EventSetup& es,
0147                                           reco::Photon::FiducialFlags& phofid,
0148                                           reco::Photon::IsolationVariables& phoisolR1,
0149                                           reco::Photon::IsolationVariables& phoisolR2,
0150                                           const HcalPFCuts* hcalCuts) 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,
0366                   hcalCuts](double outer, double inner, int depth) {
0367       return calculateHcalRecHitIso<false>(
0368           pho, caloGeometry, hcalTopo, hcalQual, hcalSev, towerMap, hbheRecHits, outer, inner, depth, hcalCuts);
0369     };
0370 
0371     auto fbc = [this,
0372                 pho,
0373                 &caloGeometry,
0374                 &hcalTopo = *hcalTopology,
0375                 &hcalQual = *hcalChannelQuality,
0376                 &hcalSev = *hcalSevLvlComputer,
0377                 &towerMap,
0378                 &hbheRecHits,
0379                 hcalCuts](double outer, int depth) {
0380       return calculateHcalRecHitIso<true>(
0381           pho, caloGeometry, hcalTopo, hcalQual, hcalSev, towerMap, hbheRecHits, outer, 0., depth, hcalCuts);
0382     };
0383 
0384     for (size_t id = 0; id < phoisolR1.hcalRecHitSumEt.size(); ++id) {
0385       const auto& outerA = detector == EcalBarrel ? hcalIsoOuterRadAEB_[id] : hcalIsoOuterRadAEE_[id];
0386       const auto& outerB = detector == EcalBarrel ? hcalIsoOuterRadBEB_[id] : hcalIsoOuterRadBEE_[id];
0387       const auto& innerA = detector == EcalBarrel ? hcalIsoInnerRadAEB_[id] : hcalIsoInnerRadAEE_[id];
0388       const auto& innerB = detector == EcalBarrel ? hcalIsoInnerRadBEB_[id] : hcalIsoInnerRadBEE_[id];
0389 
0390       phoisolR1.hcalRecHitSumEt[id] = fcone(outerA, innerA, id + 1);
0391       phoisolR2.hcalRecHitSumEt[id] = fcone(outerB, innerB, id + 1);
0392 
0393       phoisolR1.hcalRecHitSumEtBc[id] = fbc(outerA, id + 1);
0394       phoisolR2.hcalRecHitSumEtBc[id] = fbc(outerB, id + 1);
0395     }
0396   }
0397 
0398   phoisolR1.pre7DepthHcal = false;
0399   phoisolR2.pre7DepthHcal = false;
0400 }
0401 
0402 void PhotonIsolationCalculator::classify(const reco::Photon* photon,
0403                                          bool& isEBPho,
0404                                          bool& isEEPho,
0405                                          bool& isEBEtaGap,
0406                                          bool& isEBPhiGap,
0407                                          bool& isEERingGap,
0408                                          bool& isEEDeeGap,
0409                                          bool& isEBEEGap) {
0410   const reco::CaloCluster& seedCluster = *(photon->superCluster()->seed());
0411   // following line is temporary until the D. Evans or F. Ferri submit the new tag for ClusterAlgos
0412   // DEfinitive will be
0413   // DetId seedXtalId = scRef->seed()->seed();
0414   DetId seedXtalId = seedCluster.hitsAndFractions()[0].first;
0415   int detector = seedXtalId.subdetId();
0416 
0417   //Set fiducial flags for this photon.
0418   double eta = photon->superCluster()->position().eta();
0419   double feta = fabs(eta);
0420 
0421   if (detector == EcalBarrel) {
0422     isEBPho = true;
0423     if (EBDetId::isNextToEtaBoundary(EBDetId(seedXtalId))) {
0424       if (fabs(feta - 1.479) < .1)
0425         isEBEEGap = true;
0426       else
0427         isEBEtaGap = true;
0428     }
0429 
0430     if (EBDetId::isNextToPhiBoundary(EBDetId(seedXtalId)))
0431       isEBPhiGap = true;
0432 
0433   } else if (detector == EcalEndcap) {
0434     isEEPho = true;
0435     if (EEDetId::isNextToRingBoundary(EEDetId(seedXtalId))) {
0436       if (fabs(feta - 1.479) < .1)
0437         isEBEEGap = true;
0438       else
0439         isEERingGap = true;
0440     }
0441 
0442     if (EEDetId::isNextToDBoundary(EEDetId(seedXtalId)))
0443       isEEDeeGap = true;
0444   }
0445 }
0446 
0447 void PhotonIsolationCalculator::calculateTrackIso(const reco::Photon* photon,
0448                                                   const edm::Event& e,
0449                                                   double& trkCone,
0450                                                   int& ntrkCone,
0451                                                   double pTThresh,
0452                                                   double RCone,
0453                                                   double RinnerCone,
0454                                                   double etaSlice,
0455                                                   double lip,
0456                                                   double d0) const {
0457   ntrkCone = 0;
0458   trkCone = 0;
0459   //get the tracks
0460   edm::Handle<reco::TrackCollection> tracks;
0461   e.getByToken(trackInputTag_, tracks);
0462   if (!tracks.isValid()) {
0463     return;
0464   }
0465   const reco::TrackCollection* trackCollection = tracks.product();
0466   //Photon Eta and Phi.  Hope these are correct.
0467   reco::BeamSpot vertexBeamSpot;
0468   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0469   e.getByToken(beamSpotProducerTag_, recoBeamSpotHandle);
0470   vertexBeamSpot = *recoBeamSpotHandle;
0471 
0472   PhotonTkIsolation phoIso(RCone,
0473                            RinnerCone,
0474                            etaSlice,
0475                            pTThresh,
0476                            lip,
0477                            d0,
0478                            trackCollection,
0479                            math::XYZPoint(vertexBeamSpot.x0(), vertexBeamSpot.y0(), vertexBeamSpot.z0()));
0480 
0481   std::pair<int, double> res = phoIso.getIso(photon);
0482   ntrkCone = res.first;
0483   trkCone = res.second;
0484 }
0485 
0486 double PhotonIsolationCalculator::calculateEcalRecHitIso(const reco::Photon* photon,
0487                                                          const edm::Event& iEvent,
0488                                                          const edm::EventSetup& iSetup,
0489                                                          double RCone,
0490                                                          double RConeInner,
0491                                                          double etaSlice,
0492                                                          double eMin,
0493                                                          double etMin,
0494                                                          bool vetoClusteredHits,
0495                                                          bool useNumXtals) const {
0496   edm::Handle<EcalRecHitCollection> ecalhitsCollEB;
0497   edm::Handle<EcalRecHitCollection> ecalhitsCollEE;
0498   iEvent.getByToken(endcapecalCollection_, ecalhitsCollEE);
0499 
0500   iEvent.getByToken(barrelecalCollection_, ecalhitsCollEB);
0501 
0502   auto const& thresholds = iSetup.getData(ecalPFRechitThresholdsToken_);
0503 
0504   const EcalRecHitCollection* rechitsCollectionEE_ = ecalhitsCollEE.product();
0505   const EcalRecHitCollection* rechitsCollectionEB_ = ecalhitsCollEB.product();
0506 
0507   const EcalSeverityLevelAlgo* sevLevel = &iSetup.getData(ecalSevLvlToken_);
0508 
0509   edm::ESHandle<CaloGeometry> geoHandle = iSetup.getHandle(caloGeometryToken_);
0510 
0511   EgammaRecHitIsolation phoIsoEB(
0512       RCone, RConeInner, etaSlice, etMin, eMin, geoHandle, *rechitsCollectionEB_, sevLevel, DetId::Ecal);
0513 
0514   phoIsoEB.setVetoClustered(vetoClusteredHits);
0515   phoIsoEB.setUseNumCrystals(useNumXtals);
0516   phoIsoEB.doSeverityChecks(ecalhitsCollEB.product(), severityExclEB_);
0517   phoIsoEB.doFlagChecks(flagsEB_);
0518   double ecalIsolEB = phoIsoEB.getEtSum(photon, &thresholds);
0519 
0520   EgammaRecHitIsolation phoIsoEE(
0521       RCone, RConeInner, etaSlice, etMin, eMin, geoHandle, *rechitsCollectionEE_, sevLevel, DetId::Ecal);
0522 
0523   phoIsoEE.setVetoClustered(vetoClusteredHits);
0524   phoIsoEE.setUseNumCrystals(useNumXtals);
0525   phoIsoEE.doSeverityChecks(ecalhitsCollEE.product(), severityExclEE_);
0526   phoIsoEE.doFlagChecks(flagsEE_);
0527 
0528   double ecalIsolEE = phoIsoEE.getEtSum(photon, &thresholds);
0529   //  delete phoIso;
0530   double ecalIsol = ecalIsolEB + ecalIsolEE;
0531 
0532   return ecalIsol;
0533 }
0534 
0535 template <bool isoBC>
0536 double PhotonIsolationCalculator::calculateHcalRecHitIso(const reco::Photon* photon,
0537                                                          const CaloGeometry& geometry,
0538                                                          const HcalTopology& hcalTopology,
0539                                                          const HcalChannelQuality& hcalChStatus,
0540                                                          const HcalSeverityLevelComputer& hcalSevLvlComputer,
0541                                                          const CaloTowerConstituentsMap& towerMap,
0542                                                          const HBHERecHitCollection& hbheRecHits,
0543                                                          double RCone,
0544                                                          double RConeInner,
0545                                                          int depth,
0546                                                          const HcalPFCuts* hcalCuts) const {
0547   const EgammaHcalIsolation::arrayHB e04{{0., 0., 0., 0.}};
0548   const EgammaHcalIsolation::arrayHE e07{{0., 0., 0., 0., 0., 0., 0.}};
0549 
0550   if constexpr (isoBC) {
0551     auto hcaliso = EgammaHcalIsolation(EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
0552                                        RCone,
0553                                        EgammaHcalIsolation::InclusionRule::isBehindClusterSeed,
0554                                        RConeInner,
0555                                        hcalIsoEThresHB_,
0556                                        e04,
0557                                        maxHcalSeverity_,
0558                                        hcalIsoEThresHE_,
0559                                        e07,
0560                                        maxHcalSeverity_,
0561                                        hbheRecHits,
0562                                        geometry,
0563                                        hcalTopology,
0564                                        hcalChStatus,
0565                                        hcalSevLvlComputer,
0566                                        towerMap);
0567 
0568     return hcaliso.getHcalEtSumBc(photon, depth, hcalCuts);
0569   } else {
0570     auto hcaliso = EgammaHcalIsolation(EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
0571                                        RCone,
0572                                        EgammaHcalIsolation::InclusionRule::withinConeAroundCluster,
0573                                        RConeInner,
0574                                        hcalIsoEThresHB_,
0575                                        e04,
0576                                        maxHcalSeverity_,
0577                                        hcalIsoEThresHE_,
0578                                        e07,
0579                                        maxHcalSeverity_,
0580                                        hbheRecHits,
0581                                        geometry,
0582                                        hcalTopology,
0583                                        hcalChStatus,
0584                                        hcalSevLvlComputer,
0585                                        towerMap);
0586 
0587     return hcaliso.getHcalEtSum(photon, depth, hcalCuts);
0588   }
0589 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!