Reference Manual

Go to the documentation of this file.
00001 
00007 #include "RecoEgamma/PhotonIdentification/interface/PhotonIsolationCalculator.h"
00008 #include "CommonTools/Utils/interface/StringToEnumValue.h"
00009 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
00010 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
00011 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00012 #include "DataFormats/TrackReco/interface/Track.h"
00013 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00014 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
00015 #include "DataFormats/EgammaCandidates/interface/GsfElectronFwd.h"
00016 #include "DataFormats/DetId/interface/DetId.h"
00017 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00018 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00019 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
00020 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00021 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
00022 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00023 
00024 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaHcalIsolation.h"
00025 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaRecHitIsolation.h"
00026 #include "RecoEgamma/EgammaIsolationAlgos/interface/PhotonTkIsolation.h"
00027 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaEcalIsolation.h"
00028 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaTowerIsolation.h"
00029 #include "RecoCaloTools/MetaCollections/interface/CaloRecHitMetaCollections.h"
00030 //#include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
00031 
00032 #include "Geometry/CaloTopology/interface/CaloTopology.h"
00033 #include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
00034 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00035 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
00036 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
00037 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00038 
00039 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
00040 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
00041 
00042 #include <string>
00043 #include <TMath.h>
00044 
00045 
00046 void PhotonIsolationCalculator::setup(const edm::ParameterSet& conf, std::vector<int> flagsEB, std::vector<int> flagsEE, std::vector<int> severitiesEB, std::vector<int> severitiesEE) {
00047 
00048   trackInputTag_ = conf.getParameter<edm::InputTag>("trackProducer");
00049   beamSpotProducerTag_ = conf.getParameter<edm::InputTag>("beamSpotProducer");
00050   barrelecalCollection_ = conf.getParameter<edm::InputTag>("barrelEcalRecHitCollection");
00051   endcapecalCollection_ = conf.getParameter<edm::InputTag>("endcapEcalRecHitCollection");
00052   hcalCollection_ = conf.getParameter<edm::InputTag>("HcalRecHitCollection");
00053 
00054   //  gsfRecoInputTag_ = conf.getParameter<edm::InputTag>("GsfRecoCollection");
00055   modulePhiBoundary_ = conf.getParameter<double>("modulePhiBoundary");
00056   moduleEtaBoundary_ = conf.getParameter<std::vector<double> >("moduleEtaBoundary");
00057   //
00058   vetoClusteredEcalHits_ = conf.getParameter<bool>("vetoClustered");
00059   useNumCrystals_ = conf.getParameter<bool>("useNumCrystals");
00060 
00062   trkIsoBarrelRadiusA_.push_back(  conf.getParameter<double>("TrackConeOuterRadiusA_Barrel") );
00063   trkIsoBarrelRadiusA_.push_back(  conf.getParameter<double>("TrackConeInnerRadiusA_Barrel") ) ;
00064   trkIsoBarrelRadiusA_.push_back(  conf.getParameter<double>("isolationtrackThresholdA_Barrel") );
00065   trkIsoBarrelRadiusA_.push_back(  conf.getParameter<double>("longImpactParameterA_Barrel") );
00066   trkIsoBarrelRadiusA_.push_back(  conf.getParameter<double>("transImpactParameterA_Barrel") );
00067   trkIsoBarrelRadiusA_.push_back(  conf.getParameter<double>("isolationtrackEtaSliceA_Barrel") );
00068 
00069   ecalIsoBarrelRadiusA_.push_back( conf.getParameter<double>("EcalRecHitInnerRadiusA_Barrel") );
00070   ecalIsoBarrelRadiusA_.push_back( conf.getParameter<double>("EcalRecHitOuterRadiusA_Barrel") );
00071   ecalIsoBarrelRadiusA_.push_back( conf.getParameter<double>("EcalRecHitEtaSliceA_Barrel") );
00072   ecalIsoBarrelRadiusA_.push_back( conf.getParameter<double>("EcalRecHitThreshEA_Barrel") );
00073   ecalIsoBarrelRadiusA_.push_back( conf.getParameter<double>("EcalRecHitThreshEtA_Barrel") );
00074 
00075   hcalIsoBarrelRadiusA_.push_back( conf.getParameter<double>("HcalTowerInnerRadiusA_Barrel") );
00076   hcalIsoBarrelRadiusA_.push_back( conf.getParameter<double>("HcalTowerOuterRadiusA_Barrel") );
00077   hcalIsoBarrelRadiusA_.push_back( conf.getParameter<double>("HcalTowerThreshEA_Barrel") );
00078   hcalIsoBarrelRadiusA_.push_back( conf.getParameter<double>("HcalDepth1TowerInnerRadiusA_Barrel") );
00079   hcalIsoBarrelRadiusA_.push_back( conf.getParameter<double>("HcalDepth1TowerOuterRadiusA_Barrel") );
00080   hcalIsoBarrelRadiusA_.push_back( conf.getParameter<double>("HcalDepth1TowerThreshEA_Barrel") );
00081   hcalIsoBarrelRadiusA_.push_back( conf.getParameter<double>("HcalDepth2TowerInnerRadiusA_Barrel") );
00082   hcalIsoBarrelRadiusA_.push_back( conf.getParameter<double>("HcalDepth2TowerOuterRadiusA_Barrel") );
00083   hcalIsoBarrelRadiusA_.push_back( conf.getParameter<double>("HcalDepth2TowerThreshEA_Barrel") );
00084 
00085 
00086   trkIsoBarrelRadiusB_.push_back(  conf.getParameter<double>("TrackConeOuterRadiusB_Barrel") );
00087   trkIsoBarrelRadiusB_.push_back(  conf.getParameter<double>("TrackConeInnerRadiusB_Barrel") ) ;
00088   trkIsoBarrelRadiusB_.push_back(  conf.getParameter<double>("isolationtrackThresholdB_Barrel") );
00089   trkIsoBarrelRadiusB_.push_back(  conf.getParameter<double>("longImpactParameterB_Barrel") );
00090   trkIsoBarrelRadiusB_.push_back(  conf.getParameter<double>("transImpactParameterB_Barrel") );
00091   trkIsoBarrelRadiusB_.push_back(  conf.getParameter<double>("isolationtrackEtaSliceB_Barrel") );
00092 
00093   ecalIsoBarrelRadiusB_.push_back( conf.getParameter<double>("EcalRecHitInnerRadiusB_Barrel") );
00094   ecalIsoBarrelRadiusB_.push_back( conf.getParameter<double>("EcalRecHitOuterRadiusB_Barrel") );
00095   ecalIsoBarrelRadiusB_.push_back( conf.getParameter<double>("EcalRecHitEtaSliceB_Barrel") );
00096   ecalIsoBarrelRadiusB_.push_back( conf.getParameter<double>("EcalRecHitThreshEB_Barrel") );
00097   ecalIsoBarrelRadiusB_.push_back( conf.getParameter<double>("EcalRecHitThreshEtB_Barrel") );
00098 
00099   hcalIsoBarrelRadiusB_.push_back(  conf.getParameter<double>("HcalTowerInnerRadiusB_Barrel") );
00100   hcalIsoBarrelRadiusB_.push_back( conf.getParameter<double>("HcalTowerOuterRadiusB_Barrel") );
00101   hcalIsoBarrelRadiusB_.push_back( conf.getParameter<double>("HcalTowerThreshEB_Barrel") );
00102   hcalIsoBarrelRadiusB_.push_back( conf.getParameter<double>("HcalDepth1TowerInnerRadiusB_Barrel") );
00103   hcalIsoBarrelRadiusB_.push_back( conf.getParameter<double>("HcalDepth1TowerOuterRadiusB_Barrel") );
00104   hcalIsoBarrelRadiusB_.push_back( conf.getParameter<double>("HcalDepth1TowerThreshEB_Barrel") );
00105   hcalIsoBarrelRadiusB_.push_back( conf.getParameter<double>("HcalDepth2TowerInnerRadiusB_Barrel") );
00106   hcalIsoBarrelRadiusB_.push_back( conf.getParameter<double>("HcalDepth2TowerOuterRadiusB_Barrel") );
00107   hcalIsoBarrelRadiusB_.push_back( conf.getParameter<double>("HcalDepth2TowerThreshEB_Barrel") );
00108 
00110   trkIsoEndcapRadiusA_.push_back(  conf.getParameter<double>("TrackConeOuterRadiusA_Endcap") );
00111   trkIsoEndcapRadiusA_.push_back(  conf.getParameter<double>("TrackConeInnerRadiusA_Endcap") ) ;
00112   trkIsoEndcapRadiusA_.push_back(  conf.getParameter<double>("isolationtrackThresholdA_Endcap") );
00113   trkIsoEndcapRadiusA_.push_back(  conf.getParameter<double>("longImpactParameterA_Endcap") );
00114   trkIsoEndcapRadiusA_.push_back(  conf.getParameter<double>("transImpactParameterA_Endcap") );
00115   trkIsoEndcapRadiusA_.push_back(  conf.getParameter<double>("isolationtrackEtaSliceA_Endcap") );
00116 
00117 
00118   ecalIsoEndcapRadiusA_.push_back( conf.getParameter<double>("EcalRecHitInnerRadiusA_Endcap") );
00119   ecalIsoEndcapRadiusA_.push_back( conf.getParameter<double>("EcalRecHitOuterRadiusA_Endcap") );
00120   ecalIsoEndcapRadiusA_.push_back( conf.getParameter<double>("EcalRecHitEtaSliceA_Endcap") );
00121   ecalIsoEndcapRadiusA_.push_back( conf.getParameter<double>("EcalRecHitThreshEA_Endcap") );
00122   ecalIsoEndcapRadiusA_.push_back( conf.getParameter<double>("EcalRecHitThreshEtA_Endcap") );
00123 
00124 
00125   hcalIsoEndcapRadiusA_.push_back(  conf.getParameter<double>("HcalTowerInnerRadiusA_Endcap") );
00126   hcalIsoEndcapRadiusA_.push_back( conf.getParameter<double>("HcalTowerOuterRadiusA_Endcap") );
00127   hcalIsoEndcapRadiusA_.push_back( conf.getParameter<double>("HcalTowerThreshEA_Endcap") );
00128   hcalIsoEndcapRadiusA_.push_back( conf.getParameter<double>("HcalDepth1TowerInnerRadiusA_Endcap") );
00129   hcalIsoEndcapRadiusA_.push_back( conf.getParameter<double>("HcalDepth1TowerOuterRadiusA_Endcap") );
00130   hcalIsoEndcapRadiusA_.push_back( conf.getParameter<double>("HcalDepth1TowerThreshEA_Endcap") );
00131   hcalIsoEndcapRadiusA_.push_back( conf.getParameter<double>("HcalDepth2TowerInnerRadiusA_Endcap") );
00132   hcalIsoEndcapRadiusA_.push_back( conf.getParameter<double>("HcalDepth2TowerOuterRadiusA_Endcap") );
00133   hcalIsoEndcapRadiusA_.push_back( conf.getParameter<double>("HcalDepth2TowerThreshEA_Endcap") );
00134 
00135 
00136   trkIsoEndcapRadiusB_.push_back(  conf.getParameter<double>("TrackConeOuterRadiusB_Endcap") );
00137   trkIsoEndcapRadiusB_.push_back(  conf.getParameter<double>("TrackConeInnerRadiusB_Endcap") ) ;
00138   trkIsoEndcapRadiusB_.push_back(  conf.getParameter<double>("isolationtrackThresholdB_Endcap") );
00139   trkIsoEndcapRadiusB_.push_back(  conf.getParameter<double>("longImpactParameterB_Endcap") );
00140   trkIsoEndcapRadiusB_.push_back(  conf.getParameter<double>("transImpactParameterB_Endcap") );
00141   trkIsoEndcapRadiusB_.push_back(  conf.getParameter<double>("isolationtrackEtaSliceB_Endcap") );
00142 
00143   ecalIsoEndcapRadiusB_.push_back( conf.getParameter<double>("EcalRecHitInnerRadiusB_Endcap") );
00144   ecalIsoEndcapRadiusB_.push_back( conf.getParameter<double>("EcalRecHitOuterRadiusB_Endcap") );
00145   ecalIsoEndcapRadiusB_.push_back( conf.getParameter<double>("EcalRecHitEtaSliceB_Endcap") );
00146   ecalIsoEndcapRadiusB_.push_back( conf.getParameter<double>("EcalRecHitThreshEB_Endcap") );
00147   ecalIsoEndcapRadiusB_.push_back( conf.getParameter<double>("EcalRecHitThreshEtB_Endcap") );
00148 
00149   hcalIsoEndcapRadiusB_.push_back(  conf.getParameter<double>("HcalTowerInnerRadiusB_Endcap") );
00150   hcalIsoEndcapRadiusB_.push_back( conf.getParameter<double>("HcalTowerOuterRadiusB_Endcap") );
00151   hcalIsoEndcapRadiusB_.push_back( conf.getParameter<double>("HcalTowerThreshEB_Endcap") );
00152   hcalIsoEndcapRadiusB_.push_back( conf.getParameter<double>("HcalDepth1TowerInnerRadiusB_Endcap") );
00153   hcalIsoEndcapRadiusB_.push_back( conf.getParameter<double>("HcalDepth1TowerOuterRadiusB_Endcap") );
00154   hcalIsoEndcapRadiusB_.push_back( conf.getParameter<double>("HcalDepth1TowerThreshEB_Endcap") );
00155   hcalIsoEndcapRadiusB_.push_back( conf.getParameter<double>("HcalDepth2TowerInnerRadiusB_Endcap") );
00156   hcalIsoEndcapRadiusB_.push_back( conf.getParameter<double>("HcalDepth2TowerOuterRadiusB_Endcap") );
00157   hcalIsoEndcapRadiusB_.push_back( conf.getParameter<double>("HcalDepth2TowerThreshEB_Endcap") );
00158 
00159   //Pick up the variables for the spike removal
00160   flagsEB_      = flagsEB;
00161   flagsEE_      = flagsEE;
00162   severityExclEB_ = severitiesEB;
00163   severityExclEE_ = severitiesEE;
00164 }
00165 
00166 
00167 void PhotonIsolationCalculator::calculate(const reco::Photon* pho,
00168                                           const edm::Event& e,
00169                                           const edm::EventSetup& es,
00170                                           reco::Photon::FiducialFlags& phofid, 
00171                                           reco::Photon::IsolationVariables& phoisolR1, 
00172                                           reco::Photon::IsolationVariables& phoisolR2) {
00173 
00174   //Get fiducial flags. This does not really belong here
00175   bool isEBPho     = false;
00176   bool isEEPho     = false;
00177   bool isEBEtaGap  = false;
00178   bool isEBPhiGap  = false;
00179   bool isEERingGap = false;
00180   bool isEEDeeGap  = false;
00181   bool isEBEEGap   = false;
00182   classify(pho, isEBPho, isEEPho, isEBEtaGap, isEBPhiGap, isEERingGap, isEEDeeGap, isEBEEGap);
00183   phofid.isEB = isEBPho;
00184   phofid.isEE = isEEPho;
00185   phofid.isEBEtaGap  = isEBEtaGap;
00186   phofid.isEBPhiGap  = isEBPhiGap;
00187   phofid.isEERingGap = isEERingGap;
00188   phofid.isEEDeeGap  = isEEDeeGap;
00189   phofid.isEBEEGap   = isEBEEGap;
00190   
00191   // Calculate isolation variables. cone sizes and thresholds
00192   // are set for Barrel and endcap separately 
00193 
00194   reco::SuperClusterRef scRef=pho->superCluster();
00195   const reco::BasicCluster & seedCluster = *(scRef->seed()) ;
00196   DetId seedXtalId = seedCluster.hitsAndFractions()[0].first ;
00197   int detector = seedXtalId.subdetId() ;
00198 
00199   if (detector==EcalBarrel) {
00200 
00201     trackConeOuterRadiusA_    = trkIsoBarrelRadiusA_[0];
00202     trackConeInnerRadiusA_    = trkIsoBarrelRadiusA_[1];
00203     isolationtrackThresholdA_ = trkIsoBarrelRadiusA_[2];
00204     trackLipRadiusA_          = trkIsoBarrelRadiusA_[3];
00205     trackD0RadiusA_           = trkIsoBarrelRadiusA_[4];
00206     isolationtrackEtaSliceA_  = trkIsoBarrelRadiusA_[5];
00207 
00208     photonEcalRecHitConeInnerRadiusA_  = ecalIsoBarrelRadiusA_[0];
00209     photonEcalRecHitConeOuterRadiusA_  = ecalIsoBarrelRadiusA_[1];
00210     photonEcalRecHitEtaSliceA_         = ecalIsoBarrelRadiusA_[2];
00211     photonEcalRecHitThreshEA_          = ecalIsoBarrelRadiusA_[3];
00212     photonEcalRecHitThreshEtA_         = ecalIsoBarrelRadiusA_[4];
00213 
00214     photonHcalTowerConeInnerRadiusA_       = hcalIsoBarrelRadiusA_[0];
00215     photonHcalTowerConeOuterRadiusA_       = hcalIsoBarrelRadiusA_[1];
00216     photonHcalTowerThreshEA_               = hcalIsoBarrelRadiusA_[2];
00217     photonHcalDepth1TowerConeInnerRadiusA_ = hcalIsoBarrelRadiusA_[3];
00218     photonHcalDepth1TowerConeOuterRadiusA_ = hcalIsoBarrelRadiusA_[4];
00219     photonHcalDepth1TowerThreshEA_         = hcalIsoBarrelRadiusA_[5];
00220     photonHcalDepth2TowerConeInnerRadiusA_ = hcalIsoBarrelRadiusA_[6];
00221     photonHcalDepth2TowerConeOuterRadiusA_ = hcalIsoBarrelRadiusA_[7];
00222     photonHcalDepth2TowerThreshEA_         = hcalIsoBarrelRadiusA_[8];
00223 
00224 
00225     trackConeOuterRadiusB_    = trkIsoBarrelRadiusB_[0];
00226     trackConeInnerRadiusB_    = trkIsoBarrelRadiusB_[1];
00227     isolationtrackThresholdB_ = trkIsoBarrelRadiusB_[2];
00228     trackLipRadiusB_          = trkIsoBarrelRadiusB_[3];
00229     trackD0RadiusB_           = trkIsoBarrelRadiusB_[4];
00230     isolationtrackEtaSliceB_  = trkIsoBarrelRadiusB_[5];
00231 
00232     photonEcalRecHitConeInnerRadiusB_  = ecalIsoBarrelRadiusB_[0];
00233     photonEcalRecHitConeOuterRadiusB_  = ecalIsoBarrelRadiusB_[1];
00234     photonEcalRecHitEtaSliceB_         = ecalIsoBarrelRadiusB_[2];
00235     photonEcalRecHitThreshEB_          = ecalIsoBarrelRadiusB_[3];
00236     photonEcalRecHitThreshEtB_         = ecalIsoBarrelRadiusB_[4];
00237 
00238     photonHcalTowerConeInnerRadiusB_       = hcalIsoBarrelRadiusB_[0];
00239     photonHcalTowerConeOuterRadiusB_       = hcalIsoBarrelRadiusB_[1];
00240     photonHcalTowerThreshEB_               = hcalIsoBarrelRadiusB_[2];
00241     photonHcalDepth1TowerConeInnerRadiusB_ = hcalIsoBarrelRadiusB_[3];
00242     photonHcalDepth1TowerConeOuterRadiusB_ = hcalIsoBarrelRadiusB_[4];
00243     photonHcalDepth1TowerThreshEB_         = hcalIsoBarrelRadiusB_[5];
00244     photonHcalDepth2TowerConeInnerRadiusB_ = hcalIsoBarrelRadiusB_[6];
00245     photonHcalDepth2TowerConeOuterRadiusB_ = hcalIsoBarrelRadiusB_[7];
00246     photonHcalDepth2TowerThreshEB_         = hcalIsoBarrelRadiusB_[8];
00247 
00248   } else if 
00249     (detector==EcalEndcap) {
00250 
00251     trackConeOuterRadiusA_    = trkIsoEndcapRadiusA_[0];
00252     trackConeInnerRadiusA_    = trkIsoEndcapRadiusA_[1];
00253     isolationtrackThresholdA_ = trkIsoEndcapRadiusA_[2];
00254     trackLipRadiusA_          = trkIsoEndcapRadiusA_[3];
00255     trackD0RadiusA_           = trkIsoEndcapRadiusA_[4];
00256     isolationtrackEtaSliceA_  = trkIsoEndcapRadiusA_[5];
00257 
00258     photonEcalRecHitConeInnerRadiusA_  = ecalIsoEndcapRadiusA_[0];
00259     photonEcalRecHitConeOuterRadiusA_  = ecalIsoEndcapRadiusA_[1];
00260     photonEcalRecHitEtaSliceA_         = ecalIsoEndcapRadiusA_[2];
00261     photonEcalRecHitThreshEA_          = ecalIsoEndcapRadiusA_[3];
00262     photonEcalRecHitThreshEtA_         = ecalIsoEndcapRadiusA_[4];
00263 
00264     photonHcalTowerConeInnerRadiusA_       = hcalIsoEndcapRadiusA_[0];
00265     photonHcalTowerConeOuterRadiusA_       = hcalIsoEndcapRadiusA_[1];
00266     photonHcalTowerThreshEA_               = hcalIsoEndcapRadiusA_[2];
00267     photonHcalDepth1TowerConeInnerRadiusA_ = hcalIsoEndcapRadiusA_[3];
00268     photonHcalDepth1TowerConeOuterRadiusA_ = hcalIsoEndcapRadiusA_[4];
00269     photonHcalDepth1TowerThreshEA_         = hcalIsoEndcapRadiusA_[5];
00270     photonHcalDepth2TowerConeInnerRadiusA_ = hcalIsoEndcapRadiusA_[6];
00271     photonHcalDepth2TowerConeOuterRadiusA_ = hcalIsoEndcapRadiusA_[7];
00272     photonHcalDepth2TowerThreshEA_         = hcalIsoEndcapRadiusA_[8];
00273 
00274 
00275     trackConeOuterRadiusB_    = trkIsoEndcapRadiusB_[0];
00276     trackConeInnerRadiusB_    = trkIsoEndcapRadiusB_[1];
00277     isolationtrackThresholdB_ = trkIsoEndcapRadiusB_[2];
00278     trackLipRadiusB_          = trkIsoEndcapRadiusB_[3];
00279     trackD0RadiusB_           = trkIsoEndcapRadiusB_[4];
00280     isolationtrackEtaSliceB_  = trkIsoEndcapRadiusB_[5];
00281 
00282 
00283     photonEcalRecHitConeInnerRadiusB_  = ecalIsoEndcapRadiusB_[0];
00284     photonEcalRecHitConeOuterRadiusB_  = ecalIsoEndcapRadiusB_[1];
00285     photonEcalRecHitEtaSliceB_         = ecalIsoEndcapRadiusB_[2];
00286     photonEcalRecHitThreshEB_          = ecalIsoEndcapRadiusB_[3];
00287     photonEcalRecHitThreshEtB_         = ecalIsoEndcapRadiusB_[4];
00288 
00289     photonHcalTowerConeInnerRadiusB_       = hcalIsoEndcapRadiusB_[0];
00290     photonHcalTowerConeOuterRadiusB_       = hcalIsoEndcapRadiusB_[1];
00291     photonHcalTowerThreshEB_               = hcalIsoEndcapRadiusB_[2];
00292     photonHcalDepth1TowerConeInnerRadiusB_ = hcalIsoEndcapRadiusB_[3];
00293     photonHcalDepth1TowerConeOuterRadiusB_ = hcalIsoEndcapRadiusB_[4];
00294     photonHcalDepth1TowerThreshEB_         = hcalIsoEndcapRadiusB_[5];
00295     photonHcalDepth2TowerConeInnerRadiusB_ = hcalIsoEndcapRadiusB_[6];
00296     photonHcalDepth2TowerConeOuterRadiusB_ = hcalIsoEndcapRadiusB_[7];
00297     photonHcalDepth2TowerThreshEB_         = hcalIsoEndcapRadiusB_[8];
00298 
00299   }
00300 
00301  
00302   //Calculate hollow cone track isolation, CONE A
00303   int ntrkA=0;
00304   double trkisoA=0;
00305   calculateTrackIso(pho, e, trkisoA, ntrkA, isolationtrackThresholdA_,    
00306                     trackConeOuterRadiusA_, trackConeInnerRadiusA_, isolationtrackEtaSliceA_, trackLipRadiusA_, trackD0RadiusA_);
00307 
00308   //Calculate solid cone track isolation, CONE A
00309   int sntrkA=0;
00310   double strkisoA=0;
00311   calculateTrackIso(pho, e, strkisoA, sntrkA, isolationtrackThresholdA_,    
00312                     trackConeOuterRadiusA_, 0., isolationtrackEtaSliceA_, trackLipRadiusA_, trackD0RadiusA_ );
00313 
00314   phoisolR1.nTrkHollowCone = ntrkA;
00315   phoisolR1.trkSumPtHollowCone = trkisoA;
00316   phoisolR1.nTrkSolidCone = sntrkA;
00317   phoisolR1.trkSumPtSolidCone = strkisoA;
00318 
00319   //Calculate hollow cone track isolation, CONE B
00320   int ntrkB=0;
00321   double trkisoB=0;
00322   calculateTrackIso(pho, e, trkisoB, ntrkB, isolationtrackThresholdB_,    
00323                     trackConeOuterRadiusB_, trackConeInnerRadiusB_,  isolationtrackEtaSliceB_, trackLipRadiusB_, trackD0RadiusB_ );
00324 
00325   //Calculate solid cone track isolation, CONE B
00326   int sntrkB=0;
00327   double strkisoB=0;
00328   calculateTrackIso(pho, e, strkisoB, sntrkB, isolationtrackThresholdB_,    
00329                     trackConeOuterRadiusB_, 0., isolationtrackEtaSliceB_, trackLipRadiusB_, trackD0RadiusB_);
00330 
00331   phoisolR2.nTrkHollowCone = ntrkB;
00332   phoisolR2.trkSumPtHollowCone = trkisoB;
00333   phoisolR2.nTrkSolidCone = sntrkB;
00334   phoisolR2.trkSumPtSolidCone = strkisoB;
00335 
00336 //   std::cout << "Output from solid cone track isolation: ";
00337 //   std::cout << " Sum pT: " << strkiso << " ntrk: " << sntrk << std::endl;
00338   
00339   double EcalRecHitIsoA = calculateEcalRecHitIso(pho, e, es,
00340                                                  photonEcalRecHitConeOuterRadiusA_,
00341                                                  photonEcalRecHitConeInnerRadiusA_,
00342                                                  photonEcalRecHitEtaSliceA_,
00343                                                  photonEcalRecHitThreshEA_,
00344                                                  photonEcalRecHitThreshEtA_,
00345                                                  vetoClusteredEcalHits_,
00346                                                  useNumCrystals_);
00347   phoisolR1.ecalRecHitSumEt = EcalRecHitIsoA;
00348 
00349   double EcalRecHitIsoB = calculateEcalRecHitIso(pho, e, es,
00350                                                  photonEcalRecHitConeOuterRadiusB_,
00351                                                  photonEcalRecHitConeInnerRadiusB_,
00352                                                  photonEcalRecHitEtaSliceB_,
00353                                                  photonEcalRecHitThreshEB_,
00354                                                  photonEcalRecHitThreshEtB_,
00355                                                  vetoClusteredEcalHits_,
00356                                                  useNumCrystals_);
00357   phoisolR2.ecalRecHitSumEt = EcalRecHitIsoB;
00358 
00359 
00360   double HcalTowerIsoA = calculateHcalTowerIso(pho, e, es, photonHcalTowerConeOuterRadiusA_,
00361                                               photonHcalTowerConeInnerRadiusA_,
00362                                               photonHcalTowerThreshEA_, -1 );
00363   phoisolR1.hcalTowerSumEt = HcalTowerIsoA;
00364 
00365 
00366   double HcalTowerIsoB = calculateHcalTowerIso(pho, e, es, photonHcalTowerConeOuterRadiusB_,
00367                                               photonHcalTowerConeInnerRadiusB_,
00368                                               photonHcalTowerThreshEB_, -1 );
00369   phoisolR2.hcalTowerSumEt = HcalTowerIsoB;
00370 
00372 
00373   double HcalDepth1TowerIsoA = calculateHcalTowerIso(pho, e, es, photonHcalDepth1TowerConeOuterRadiusA_,
00374                                               photonHcalDepth1TowerConeInnerRadiusA_,
00375                                               photonHcalDepth1TowerThreshEA_, 1 );
00376   phoisolR1.hcalDepth1TowerSumEt = HcalDepth1TowerIsoA;
00377 
00378 
00379   double HcalDepth1TowerIsoB = calculateHcalTowerIso(pho, e, es, photonHcalDepth1TowerConeOuterRadiusB_,
00380                                               photonHcalDepth1TowerConeInnerRadiusB_,
00381                                               photonHcalDepth1TowerThreshEB_, 1 );
00382   phoisolR2.hcalDepth1TowerSumEt = HcalDepth1TowerIsoB;
00383 
00384 
00385 
00387 
00388   double HcalDepth2TowerIsoA = calculateHcalTowerIso(pho, e, es, photonHcalDepth2TowerConeOuterRadiusA_,
00389                                                        photonHcalDepth2TowerConeInnerRadiusA_,
00390                                                        photonHcalDepth2TowerThreshEA_, 2 );
00391   phoisolR1.hcalDepth2TowerSumEt = HcalDepth2TowerIsoA;
00392 
00393 
00394   double HcalDepth2TowerIsoB = calculateHcalTowerIso(pho, e, es, photonHcalDepth2TowerConeOuterRadiusB_,
00395                                               photonHcalDepth2TowerConeInnerRadiusB_,
00396                                               photonHcalDepth2TowerThreshEB_, 2 );
00397   phoisolR2.hcalDepth2TowerSumEt = HcalDepth2TowerIsoB;
00398 
00399 
00400   // New Hcal isolation based on the new H/E definition (towers behind the BCs in the SC are used to evaluated H)
00401   double HcalTowerBcIsoA = calculateHcalTowerIso(pho, e, es, photonHcalTowerConeOuterRadiusA_,
00402                                                  photonHcalTowerThreshEA_, -1 );
00403   phoisolR1.hcalTowerSumEtBc = HcalTowerBcIsoA;
00404 
00405 
00406   double HcalTowerBcIsoB = calculateHcalTowerIso(pho, e, es, photonHcalTowerConeOuterRadiusB_,
00407                                                  photonHcalTowerThreshEB_, -1 );
00408   phoisolR2.hcalTowerSumEtBc = HcalTowerBcIsoB;
00409 
00411 
00412   double HcalDepth1TowerBcIsoA = calculateHcalTowerIso(pho, e, es, photonHcalDepth1TowerConeOuterRadiusA_,
00413                                                        photonHcalDepth1TowerThreshEA_, 1 );
00414   phoisolR1.hcalDepth1TowerSumEtBc = HcalDepth1TowerBcIsoA;
00415 
00416 
00417   double HcalDepth1TowerBcIsoB = calculateHcalTowerIso(pho, e, es, photonHcalDepth1TowerConeOuterRadiusB_,
00418                                                        photonHcalDepth1TowerThreshEB_, 1 );
00419   phoisolR2.hcalDepth1TowerSumEtBc = HcalDepth1TowerBcIsoB;
00420 
00421 
00422 
00424 
00425   double HcalDepth2TowerBcIsoA = calculateHcalTowerIso(pho, e, es, photonHcalDepth2TowerConeOuterRadiusA_,
00426                                                        photonHcalDepth2TowerThreshEA_, 2 );
00427   phoisolR1.hcalDepth2TowerSumEtBc = HcalDepth2TowerBcIsoA;
00428 
00429 
00430   double HcalDepth2TowerBcIsoB = calculateHcalTowerIso(pho, e, es, photonHcalDepth2TowerConeOuterRadiusB_,
00431                                                        photonHcalDepth2TowerThreshEB_, 2 );
00432   phoisolR2.hcalDepth2TowerSumEtBc = HcalDepth2TowerBcIsoB;
00433 
00434 
00435 
00436 }
00437 
00438 
00439 void PhotonIsolationCalculator::classify(const reco::Photon* photon, 
00440                                          bool &isEBPho,
00441                                          bool &isEEPho,
00442                                          bool &isEBEtaGap,
00443                                          bool &isEBPhiGap,
00444                                          bool &isEERingGap,
00445                                          bool &isEEDeeGap,
00446                                          bool &isEBEEGap){
00447 
00448 
00449   const reco::CaloCluster & seedCluster = *(photon->superCluster()->seed()) ;
00450   // following line is temporary until the D. Evans or F. Ferri submit the new tag for ClusterAlgos 
00451   // DEfinitive will be 
00452   // DetId seedXtalId = scRef->seed()->seed();
00453   DetId seedXtalId = seedCluster.hitsAndFractions()[0].first ;
00454   int detector = seedXtalId.subdetId() ;
00455 
00456   //Set fiducial flags for this photon.
00457   double eta = photon->superCluster()->position().eta();
00458   double feta = fabs(eta);
00459 
00460   if (detector==EcalBarrel) {
00461 
00462 
00463     isEBPho = true;
00464     if (EBDetId::isNextToEtaBoundary(EBDetId(seedXtalId))) {
00465       if (fabs(feta-1.479)<.1) isEBEEGap = true ; 
00466       else
00467         isEBEtaGap = true ; 
00468     }
00469 
00470     if (EBDetId::isNextToPhiBoundary(EBDetId(seedXtalId)))
00471       isEBPhiGap = true ; 
00472 
00473 
00474 
00475   } else if (detector==EcalEndcap) {
00476 
00477     isEEPho = true;
00478     if (EEDetId::isNextToRingBoundary(EEDetId(seedXtalId))) {
00479       if (fabs(feta-1.479)<.1) isEBEEGap = true ; 
00480       else
00481         isEERingGap = true ; 
00482     }
00483 
00484     if (EEDetId::isNextToDBoundary(EEDetId(seedXtalId)))
00485       isEEDeeGap = true ; 
00486   }
00487 
00488 
00489 }
00490 
00491 void PhotonIsolationCalculator::calculateTrackIso(const reco::Photon* photon,
00492                                                   const edm::Event& e,
00493                                                   double &trkCone,
00494                                                   int &ntrkCone,
00495                                                   double pTThresh,
00496                                                   double RCone,
00497                                                   double RinnerCone,
00498                                                   double etaSlice, 
00499                                                   double lip, 
00500                                                   double d0){
00501   int counter  =0;
00502   double ptSum =0.;
00503   
00504   
00505   //get the tracks
00506   edm::Handle<reco::TrackCollection> tracks;
00507   e.getByLabel(trackInputTag_,tracks);
00508   if(!tracks.isValid()) {
00509     return;
00510   }
00511   const reco::TrackCollection* trackCollection = tracks.product();
00512   //Photon Eta and Phi.  Hope these are correct.
00513   reco::BeamSpot vertexBeamSpot;
00514   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
00515   e.getByLabel(beamSpotProducerTag_,recoBeamSpotHandle);
00516   vertexBeamSpot = *recoBeamSpotHandle;
00517   
00518   PhotonTkIsolation phoIso(RCone, 
00519                            RinnerCone, 
00520                            etaSlice,  
00521                            pTThresh, 
00522                            lip , 
00523                            d0, 
00524                            trackCollection, 
00525                            math::XYZPoint(vertexBeamSpot.x0(),vertexBeamSpot.y0(),vertexBeamSpot.z0()));
00526 
00527   counter = phoIso.getNumberTracks(photon);
00528   ptSum = phoIso.getPtTracks(photon);
00529   ntrkCone = counter;
00530   trkCone = ptSum;
00531 }
00532 
00533 
00534 
00535 double PhotonIsolationCalculator::calculateEcalRecHitIso(const reco::Photon* photon,
00536                                                          const edm::Event& iEvent,
00537                                                          const edm::EventSetup& iSetup,
00538                                                          double RCone,
00539                                                          double RConeInner,
00540                                                          double etaSlice,
00541                                                          double eMin,
00542                                                          double etMin,
00543                                                          bool vetoClusteredHits,
00544                                                          bool useNumXtals)
00545 {
00546   
00547   edm::Handle<EcalRecHitCollection> ecalhitsCollEB;
00548   edm::Handle<EcalRecHitCollection> ecalhitsCollEE;
00549   iEvent.getByLabel(endcapecalCollection_, ecalhitsCollEE);
00550  
00551   iEvent.getByLabel(barrelecalCollection_, ecalhitsCollEB);
00552  
00553   const EcalRecHitCollection* rechitsCollectionEE_ = ecalhitsCollEE.product();
00554   const EcalRecHitCollection* rechitsCollectionEB_ = ecalhitsCollEB.product();
00555 
00556   edm::ESHandle<EcalSeverityLevelAlgo> sevlv;
00557   iSetup.get<EcalSeverityLevelAlgoRcd>().get(sevlv);
00558   const EcalSeverityLevelAlgo* sevLevel = sevlv.product();
00559 
00560   std::auto_ptr<CaloRecHitMetaCollectionV> RecHitsEE(0); 
00561   RecHitsEE = std::auto_ptr<CaloRecHitMetaCollectionV>(new EcalRecHitMetaCollection(*rechitsCollectionEE_));
00562  
00563   std::auto_ptr<CaloRecHitMetaCollectionV> RecHitsEB(0); 
00564   RecHitsEB = std::auto_ptr<CaloRecHitMetaCollectionV>(new EcalRecHitMetaCollection(*rechitsCollectionEB_));
00565 
00566   edm::ESHandle<CaloGeometry> geoHandle;
00567   iSetup.get<CaloGeometryRecord>().get(geoHandle);
00568 
00569   EgammaRecHitIsolation phoIsoEB(RCone,
00570                                  RConeInner,
00571                                  etaSlice,
00572                                  etMin,
00573                                  eMin,
00574                                  geoHandle,
00575                                  &(*RecHitsEB),
00576                                  sevLevel,
00577                                  DetId::Ecal);
00578 
00579   phoIsoEB.setVetoClustered(vetoClusteredHits);
00580   phoIsoEB.setUseNumCrystals(useNumXtals);
00581   phoIsoEB.doSeverityChecks(ecalhitsCollEB.product(), severityExclEB_);
00582   phoIsoEB.doFlagChecks(flagsEB_);
00583   double ecalIsolEB = phoIsoEB.getEtSum(photon);
00584   
00585   EgammaRecHitIsolation phoIsoEE(RCone,
00586                                  RConeInner,
00587                                  etaSlice,
00588                                  etMin,
00589                                  eMin,
00590                                  geoHandle,
00591                                  &(*RecHitsEE),
00592                                  sevLevel,
00593                                  DetId::Ecal);
00594   
00595   phoIsoEE.setVetoClustered(vetoClusteredHits);
00596   phoIsoEE.setUseNumCrystals(useNumXtals);
00597   phoIsoEE.doSeverityChecks(ecalhitsCollEE.product(), severityExclEE_); 
00598   phoIsoEE.doFlagChecks(flagsEE_);
00599 
00600   double ecalIsolEE = phoIsoEE.getEtSum(photon);
00601   //  delete phoIso;
00602   double ecalIsol = ecalIsolEB + ecalIsolEE;
00603   
00604   return ecalIsol;
00605 }
00606 
00607 double PhotonIsolationCalculator::calculateHcalTowerIso(const reco::Photon* photon,
00608                                                         const edm::Event& iEvent,
00609                                                         const edm::EventSetup& iSetup,
00610                                                         double RCone,
00611                                                         double RConeInner,
00612                                                         double eMin,
00613                                                         signed int depth )
00614 {
00615 
00616   edm::Handle<CaloTowerCollection> hcalhitsCollH;
00617  
00618   iEvent.getByLabel(hcalCollection_, hcalhitsCollH);
00619   
00620   const CaloTowerCollection *toww = hcalhitsCollH.product();
00621 
00622   double hcalIsol=0.;
00623   
00624   //std::cout << "before iso call" << std::endl;
00625   EgammaTowerIsolation phoIso(RCone,
00626                               RConeInner,
00627                               eMin,depth,
00628                               toww);
00629   hcalIsol = phoIso.getTowerEtSum(photon);
00630   //  delete phoIso;
00631   //std::cout << "after call" << std::endl;
00632   return hcalIsol;
00633   
00634 
00635 }
00636 
00637 
00638 
00639 double PhotonIsolationCalculator::calculateHcalTowerIso(const reco::Photon* photon,
00640                                                         const edm::Event& iEvent,
00641                                                         const edm::EventSetup& iSetup,
00642                                                         double RCone,
00643                                                         double eMin,
00644                                                         signed int depth )
00645 {
00646 
00647   edm::Handle<CaloTowerCollection> hcalhitsCollH;
00648  
00649   iEvent.getByLabel(hcalCollection_, hcalhitsCollH);
00650   
00651   const CaloTowerCollection *toww = hcalhitsCollH.product();
00652 
00653   double hcalIsol=0.;
00654   
00655   //std::cout << "before iso call" << std::endl;
00656   EgammaTowerIsolation phoIso(RCone,
00657                               0.,
00658                               eMin,depth,
00659                               toww);
00660   hcalIsol = phoIso.getTowerEtSum(photon, &(photon->hcalTowersBehindClusters()) );
00661   //  delete phoIso;
00662   //std::cout << "after call" << std::endl;
00663   return hcalIsol;
00664   
00665 
00666 }
00667 
00668 
00669 
00670
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/RecoEgamma/PhotonIdentification/src/PhotonIsolationCalculator.cc
