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