Reference Manual

Go to the documentation of this file.
00001 
00007 #include "RecoEgamma/PhotonIdentification/interface/PhotonIsolationCalculator.h"
00008 
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 
00048 
00049   trackInputTag_ = conf.getParameter<edm::InputTag>("trackProducer");
00050   beamSpotProducerTag_ = conf.getParameter<edm::InputTag>("beamSpotProducer");
00051   barrelecalCollection_ = conf.getParameter<edm::InputTag>("barrelEcalRecHitCollection");
00052   endcapecalCollection_ = conf.getParameter<edm::InputTag>("endcapEcalRecHitCollection");
00053   hcalCollection_ = conf.getParameter<edm::InputTag>("HcalRecHitCollection");
00054 
00055   //  gsfRecoInputTag_ = conf.getParameter<edm::InputTag>("GsfRecoCollection");
00056   modulePhiBoundary_ = conf.getParameter<double>("modulePhiBoundary");
00057   moduleEtaBoundary_ = conf.getParameter<std::vector<double> >("moduleEtaBoundary");
00058   //
00059   vetoClusteredEcalHits_ = conf.getParameter<bool>("vetoClustered");
00060   useNumCrystals_ = conf.getParameter<bool>("useNumCrystals");
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   severityLevelCut_        = conf.getParameter<int>("severityLevelCut");
00161   //severityRecHitThreshold_ = conf.getParameter<double>("severityRecHitThreshold");
00162   //spikeIdThreshold_        = conf.getParameter<double>("spikeIdThreshold");
00163 
00164   const std::vector<int> &tempV(conf.getParameter<std::vector<int> >("recHitFlagsToBeExcluded"));
00165   v_chstatus_.clear();
00166   v_chstatus_.insert(v_chstatus_.begin(),tempV.begin(),tempV.end());
00167   std::sort( v_chstatus_.begin(), v_chstatus_.end() );
00168 
00169 
00170   //Need to figure out which algo to use
00171   //if(!conf.getParameter<std::string>("spikeIdString").compare("kE1OverE9") )   {
00172   //  spId_ = EcalSeverityLevelAlgo::kE1OverE9;
00173   //} else if(!conf.getParameter<std::string>("spikeIdString").compare("kSwissCross") ) {
00174   //  spId_ = EcalSeverityLevelAlgo::kSwissCross;
00175   //} else if(!conf.getParameter<std::string>("spikeIdString").compare("kSwissCrossBordersIncluded") ) {
00176   //  spId_ = EcalSeverityLevelAlgo::kSwissCrossBordersIncluded;
00177   //} else {
00178   //  spId_ = EcalSeverityLevelAlgo::kSwissCrossBordersIncluded;
00179   //   edm::LogWarning("PhotonIsolationCalculator|SpikeRemovalForIsolation")
00180   //     << "Cannot find the requested method. kSwissCross set instead.";
00181   //}
00182 }
00183 
00184 
00185 void PhotonIsolationCalculator::calculate(const reco::Photon* pho,
00186                                      const edm::Event& e,
00187                                      const edm::EventSetup& es,
00188                                      reco::Photon::FiducialFlags& phofid, 
00189                                      reco::Photon::IsolationVariables& phoisolR1, 
00190                                      reco::Photon::IsolationVariables& phoisolR2){
00191 
00192 
00193   //Get fiducial flags. This does not really belong here
00194   bool isEBPho     = false;
00195   bool isEEPho     = false;
00196   bool isEBEtaGap  = false;
00197   bool isEBPhiGap  = false;
00198   bool isEERingGap = false;
00199   bool isEEDeeGap  = false;
00200   bool isEBEEGap   = false;
00201   classify(pho, isEBPho, isEEPho, isEBEtaGap, isEBPhiGap, isEERingGap, isEEDeeGap, isEBEEGap);
00202   phofid.isEB = isEBPho;
00203   phofid.isEE = isEEPho;
00204   phofid.isEBEtaGap  = isEBEtaGap;
00205   phofid.isEBPhiGap  = isEBPhiGap;
00206   phofid.isEERingGap = isEERingGap;
00207   phofid.isEEDeeGap  = isEEDeeGap;
00208   phofid.isEBEEGap   = isEBEEGap;
00209   
00210   // Calculate isolation variables. cone sizes and thresholds
00211   // are set for Barrel and endcap separately 
00212 
00213   reco::SuperClusterRef scRef=pho->superCluster();
00214   const reco::BasicCluster & seedCluster = *(scRef->seed()) ;
00215   DetId seedXtalId = seedCluster.hitsAndFractions()[0].first ;
00216   int detector = seedXtalId.subdetId() ;
00217 
00218   if (detector==EcalBarrel) {
00219 
00220     trackConeOuterRadiusA_    = trkIsoBarrelRadiusA_[0];
00221     trackConeInnerRadiusA_    = trkIsoBarrelRadiusA_[1];
00222     isolationtrackThresholdA_ = trkIsoBarrelRadiusA_[2];
00223     trackLipRadiusA_          = trkIsoBarrelRadiusA_[3];
00224     trackD0RadiusA_           = trkIsoBarrelRadiusA_[4];
00225     isolationtrackEtaSliceA_  = trkIsoBarrelRadiusA_[5];
00226 
00227     photonEcalRecHitConeInnerRadiusA_  = ecalIsoBarrelRadiusA_[0];
00228     photonEcalRecHitConeOuterRadiusA_  = ecalIsoBarrelRadiusA_[1];
00229     photonEcalRecHitEtaSliceA_         = ecalIsoBarrelRadiusA_[2];
00230     photonEcalRecHitThreshEA_          = ecalIsoBarrelRadiusA_[3];
00231     photonEcalRecHitThreshEtA_         = ecalIsoBarrelRadiusA_[4];
00232 
00233     photonHcalTowerConeInnerRadiusA_       = hcalIsoBarrelRadiusA_[0];
00234     photonHcalTowerConeOuterRadiusA_       = hcalIsoBarrelRadiusA_[1];
00235     photonHcalTowerThreshEA_               = hcalIsoBarrelRadiusA_[2];
00236     photonHcalDepth1TowerConeInnerRadiusA_ = hcalIsoBarrelRadiusA_[3];
00237     photonHcalDepth1TowerConeOuterRadiusA_ = hcalIsoBarrelRadiusA_[4];
00238     photonHcalDepth1TowerThreshEA_         = hcalIsoBarrelRadiusA_[5];
00239     photonHcalDepth2TowerConeInnerRadiusA_ = hcalIsoBarrelRadiusA_[6];
00240     photonHcalDepth2TowerConeOuterRadiusA_ = hcalIsoBarrelRadiusA_[7];
00241     photonHcalDepth2TowerThreshEA_         = hcalIsoBarrelRadiusA_[8];
00242 
00243 
00244     trackConeOuterRadiusB_    = trkIsoBarrelRadiusB_[0];
00245     trackConeInnerRadiusB_    = trkIsoBarrelRadiusB_[1];
00246     isolationtrackThresholdB_ = trkIsoBarrelRadiusB_[2];
00247     trackLipRadiusB_          = trkIsoBarrelRadiusB_[3];
00248     trackD0RadiusB_           = trkIsoBarrelRadiusB_[4];
00249     isolationtrackEtaSliceB_  = trkIsoBarrelRadiusB_[5];
00250 
00251     photonEcalRecHitConeInnerRadiusB_  = ecalIsoBarrelRadiusB_[0];
00252     photonEcalRecHitConeOuterRadiusB_  = ecalIsoBarrelRadiusB_[1];
00253     photonEcalRecHitEtaSliceB_         = ecalIsoBarrelRadiusB_[2];
00254     photonEcalRecHitThreshEB_          = ecalIsoBarrelRadiusB_[3];
00255     photonEcalRecHitThreshEtB_         = ecalIsoBarrelRadiusB_[4];
00256 
00257     photonHcalTowerConeInnerRadiusB_       = hcalIsoBarrelRadiusB_[0];
00258     photonHcalTowerConeOuterRadiusB_       = hcalIsoBarrelRadiusB_[1];
00259     photonHcalTowerThreshEB_               = hcalIsoBarrelRadiusB_[2];
00260     photonHcalDepth1TowerConeInnerRadiusB_ = hcalIsoBarrelRadiusB_[3];
00261     photonHcalDepth1TowerConeOuterRadiusB_ = hcalIsoBarrelRadiusB_[4];
00262     photonHcalDepth1TowerThreshEB_         = hcalIsoBarrelRadiusB_[5];
00263     photonHcalDepth2TowerConeInnerRadiusB_ = hcalIsoBarrelRadiusB_[6];
00264     photonHcalDepth2TowerConeOuterRadiusB_ = hcalIsoBarrelRadiusB_[7];
00265     photonHcalDepth2TowerThreshEB_         = hcalIsoBarrelRadiusB_[8];
00266 
00267   } else if 
00268     (detector==EcalEndcap) {
00269 
00270     trackConeOuterRadiusA_    = trkIsoEndcapRadiusA_[0];
00271     trackConeInnerRadiusA_    = trkIsoEndcapRadiusA_[1];
00272     isolationtrackThresholdA_ = trkIsoEndcapRadiusA_[2];
00273     trackLipRadiusA_          = trkIsoEndcapRadiusA_[3];
00274     trackD0RadiusA_           = trkIsoEndcapRadiusA_[4];
00275     isolationtrackEtaSliceA_  = trkIsoEndcapRadiusA_[5];
00276 
00277     photonEcalRecHitConeInnerRadiusA_  = ecalIsoEndcapRadiusA_[0];
00278     photonEcalRecHitConeOuterRadiusA_  = ecalIsoEndcapRadiusA_[1];
00279     photonEcalRecHitEtaSliceA_         = ecalIsoEndcapRadiusA_[2];
00280     photonEcalRecHitThreshEA_          = ecalIsoEndcapRadiusA_[3];
00281     photonEcalRecHitThreshEtA_         = ecalIsoEndcapRadiusA_[4];
00282 
00283     photonHcalTowerConeInnerRadiusA_       = hcalIsoEndcapRadiusA_[0];
00284     photonHcalTowerConeOuterRadiusA_       = hcalIsoEndcapRadiusA_[1];
00285     photonHcalTowerThreshEA_               = hcalIsoEndcapRadiusA_[2];
00286     photonHcalDepth1TowerConeInnerRadiusA_ = hcalIsoEndcapRadiusA_[3];
00287     photonHcalDepth1TowerConeOuterRadiusA_ = hcalIsoEndcapRadiusA_[4];
00288     photonHcalDepth1TowerThreshEA_         = hcalIsoEndcapRadiusA_[5];
00289     photonHcalDepth2TowerConeInnerRadiusA_ = hcalIsoEndcapRadiusA_[6];
00290     photonHcalDepth2TowerConeOuterRadiusA_ = hcalIsoEndcapRadiusA_[7];
00291     photonHcalDepth2TowerThreshEA_         = hcalIsoEndcapRadiusA_[8];
00292 
00293 
00294     trackConeOuterRadiusB_    = trkIsoEndcapRadiusB_[0];
00295     trackConeInnerRadiusB_    = trkIsoEndcapRadiusB_[1];
00296     isolationtrackThresholdB_ = trkIsoEndcapRadiusB_[2];
00297     trackLipRadiusB_          = trkIsoEndcapRadiusB_[3];
00298     trackD0RadiusB_           = trkIsoEndcapRadiusB_[4];
00299     isolationtrackEtaSliceB_  = trkIsoEndcapRadiusB_[5];
00300 
00301 
00302     photonEcalRecHitConeInnerRadiusB_  = ecalIsoEndcapRadiusB_[0];
00303     photonEcalRecHitConeOuterRadiusB_  = ecalIsoEndcapRadiusB_[1];
00304     photonEcalRecHitEtaSliceB_         = ecalIsoEndcapRadiusB_[2];
00305     photonEcalRecHitThreshEB_          = ecalIsoEndcapRadiusB_[3];
00306     photonEcalRecHitThreshEtB_         = ecalIsoEndcapRadiusB_[4];
00307 
00308     photonHcalTowerConeInnerRadiusB_       = hcalIsoEndcapRadiusB_[0];
00309     photonHcalTowerConeOuterRadiusB_       = hcalIsoEndcapRadiusB_[1];
00310     photonHcalTowerThreshEB_               = hcalIsoEndcapRadiusB_[2];
00311     photonHcalDepth1TowerConeInnerRadiusB_ = hcalIsoEndcapRadiusB_[3];
00312     photonHcalDepth1TowerConeOuterRadiusB_ = hcalIsoEndcapRadiusB_[4];
00313     photonHcalDepth1TowerThreshEB_         = hcalIsoEndcapRadiusB_[5];
00314     photonHcalDepth2TowerConeInnerRadiusB_ = hcalIsoEndcapRadiusB_[6];
00315     photonHcalDepth2TowerConeOuterRadiusB_ = hcalIsoEndcapRadiusB_[7];
00316     photonHcalDepth2TowerThreshEB_         = hcalIsoEndcapRadiusB_[8];
00317 
00318   }
00319 
00320  
00321   //Calculate hollow cone track isolation, CONE A
00322   int ntrkA=0;
00323   double trkisoA=0;
00324   calculateTrackIso(pho, e, trkisoA, ntrkA, isolationtrackThresholdA_,    
00325                     trackConeOuterRadiusA_, trackConeInnerRadiusA_, isolationtrackEtaSliceA_, trackLipRadiusA_, trackD0RadiusA_);
00326 
00327   //Calculate solid cone track isolation, CONE A
00328   int sntrkA=0;
00329   double strkisoA=0;
00330   calculateTrackIso(pho, e, strkisoA, sntrkA, isolationtrackThresholdA_,    
00331                     trackConeOuterRadiusA_, 0., isolationtrackEtaSliceA_, trackLipRadiusA_, trackD0RadiusA_ );
00332 
00333   phoisolR1.nTrkHollowCone = ntrkA;
00334   phoisolR1.trkSumPtHollowCone = trkisoA;
00335   phoisolR1.nTrkSolidCone = sntrkA;
00336   phoisolR1.trkSumPtSolidCone = strkisoA;
00337 
00338   //Calculate hollow cone track isolation, CONE B
00339   int ntrkB=0;
00340   double trkisoB=0;
00341   calculateTrackIso(pho, e, trkisoB, ntrkB, isolationtrackThresholdB_,    
00342                     trackConeOuterRadiusB_, trackConeInnerRadiusB_,  isolationtrackEtaSliceB_, trackLipRadiusB_, trackD0RadiusB_ );
00343 
00344   //Calculate solid cone track isolation, CONE B
00345   int sntrkB=0;
00346   double strkisoB=0;
00347   calculateTrackIso(pho, e, strkisoB, sntrkB, isolationtrackThresholdB_,    
00348                     trackConeOuterRadiusB_, 0., isolationtrackEtaSliceB_, trackLipRadiusB_, trackD0RadiusB_);
00349 
00350   phoisolR2.nTrkHollowCone = ntrkB;
00351   phoisolR2.trkSumPtHollowCone = trkisoB;
00352   phoisolR2.nTrkSolidCone = sntrkB;
00353   phoisolR2.trkSumPtSolidCone = strkisoB;
00354 
00355 //   std::cout << "Output from solid cone track isolation: ";
00356 //   std::cout << " Sum pT: " << strkiso << " ntrk: " << sntrk << std::endl;
00357   
00358   double EcalRecHitIsoA = calculateEcalRecHitIso(pho, e, es,
00359                                                  photonEcalRecHitConeOuterRadiusA_,
00360                                                  photonEcalRecHitConeInnerRadiusA_,
00361                                                  photonEcalRecHitEtaSliceA_,
00362                                                  photonEcalRecHitThreshEA_,
00363                                                  photonEcalRecHitThreshEtA_,
00364                                                  vetoClusteredEcalHits_,
00365                                                  useNumCrystals_);
00366   phoisolR1.ecalRecHitSumEt = EcalRecHitIsoA;
00367 
00368   double EcalRecHitIsoB = calculateEcalRecHitIso(pho, e, es,
00369                                                  photonEcalRecHitConeOuterRadiusB_,
00370                                                  photonEcalRecHitConeInnerRadiusB_,
00371                                                  photonEcalRecHitEtaSliceB_,
00372                                                  photonEcalRecHitThreshEB_,
00373                                                  photonEcalRecHitThreshEtB_,
00374                                                  vetoClusteredEcalHits_,
00375                                                  useNumCrystals_);
00376   phoisolR2.ecalRecHitSumEt = EcalRecHitIsoB;
00377 
00378   double HcalTowerIsoA = calculateHcalTowerIso(pho, e, es, photonHcalTowerConeOuterRadiusA_,
00379                                               photonHcalTowerConeInnerRadiusA_,
00380                                               photonHcalTowerThreshEA_, -1 );
00381   phoisolR1.hcalTowerSumEt = HcalTowerIsoA;
00382 
00383 
00384   double HcalTowerIsoB = calculateHcalTowerIso(pho, e, es, photonHcalTowerConeOuterRadiusB_,
00385                                               photonHcalTowerConeInnerRadiusB_,
00386                                               photonHcalTowerThreshEB_, -1 );
00387   phoisolR2.hcalTowerSumEt = HcalTowerIsoB;
00388 
00390 
00391   double HcalDepth1TowerIsoA = calculateHcalTowerIso(pho, e, es, photonHcalDepth1TowerConeOuterRadiusA_,
00392                                               photonHcalDepth1TowerConeInnerRadiusA_,
00393                                               photonHcalDepth1TowerThreshEA_, 1 );
00394   phoisolR1.hcalDepth1TowerSumEt = HcalDepth1TowerIsoA;
00395 
00396 
00397   double HcalDepth1TowerIsoB = calculateHcalTowerIso(pho, e, es, photonHcalDepth1TowerConeOuterRadiusB_,
00398                                               photonHcalDepth1TowerConeInnerRadiusB_,
00399                                               photonHcalDepth1TowerThreshEB_, 1 );
00400   phoisolR2.hcalDepth1TowerSumEt = HcalDepth1TowerIsoB;
00401 
00402 
00403 
00405 
00406   double HcalDepth2TowerIsoA = calculateHcalTowerIso(pho, e, es, photonHcalDepth2TowerConeOuterRadiusA_,
00407                                               photonHcalDepth2TowerConeInnerRadiusA_,
00408                                               photonHcalDepth2TowerThreshEA_, 2 );
00409   phoisolR1.hcalDepth2TowerSumEt = HcalDepth2TowerIsoA;
00410 
00411 
00412   double HcalDepth2TowerIsoB = calculateHcalTowerIso(pho, e, es, photonHcalDepth2TowerConeOuterRadiusB_,
00413                                               photonHcalDepth2TowerConeInnerRadiusB_,
00414                                               photonHcalDepth2TowerThreshEB_, 2 );
00415   phoisolR2.hcalDepth2TowerSumEt = HcalDepth2TowerIsoB;
00416 
00417 
00418 
00419 
00420 
00421 
00422 }
00423 
00424 
00425 void PhotonIsolationCalculator::classify(const reco::Photon* photon, 
00426                                          bool &isEBPho,
00427                                          bool &isEEPho,
00428                                          bool &isEBEtaGap,
00429                                          bool &isEBPhiGap,
00430                                          bool &isEERingGap,
00431                                          bool &isEEDeeGap,
00432                                          bool &isEBEEGap){
00433 
00434 
00435   const reco::CaloCluster & seedCluster = *(photon->superCluster()->seed()) ;
00436   // following line is temporary until the D. Evans or F. Ferri submit the new tag for ClusterAlgos 
00437   // DEfinitive will be 
00438   // DetId seedXtalId = scRef->seed()->seed();
00439   DetId seedXtalId = seedCluster.hitsAndFractions()[0].first ;
00440   int detector = seedXtalId.subdetId() ;
00441 
00442   //Set fiducial flags for this photon.
00443   double eta = photon->superCluster()->position().eta();
00444   double feta = fabs(eta);
00445 
00446   if (detector==EcalBarrel) {
00447 
00448 
00449     isEBPho = true;
00450     if (EBDetId::isNextToEtaBoundary(EBDetId(seedXtalId))) {
00451       if (fabs(feta-1.479)<.1) isEBEEGap = true ; 
00452       else
00453         isEBEtaGap = true ; 
00454     }
00455 
00456     if (EBDetId::isNextToPhiBoundary(EBDetId(seedXtalId)))
00457       isEBPhiGap = true ; 
00458 
00459 
00460 
00461   } else if (detector==EcalEndcap) {
00462 
00463     isEEPho = true;
00464     if (EEDetId::isNextToRingBoundary(EEDetId(seedXtalId))) {
00465       if (fabs(feta-1.479)<.1) isEBEEGap = true ; 
00466       else
00467         isEERingGap = true ; 
00468     }
00469 
00470     if (EEDetId::isNextToDBoundary(EEDetId(seedXtalId)))
00471       isEEDeeGap = true ; 
00472   }
00473 
00474 
00475 }
00476 
00477 void PhotonIsolationCalculator::calculateTrackIso(const reco::Photon* photon,
00478                                                   const edm::Event& e,
00479                                                   double &trkCone,
00480                                                   int &ntrkCone,
00481                                                   double pTThresh,
00482                                                   double RCone,
00483                                                   double RinnerCone,
00484                                                   double etaSlice, 
00485                                                   double lip, 
00486                                                   double d0){
00487   int counter  =0;
00488   double ptSum =0.;
00489   
00490   
00491   //get the tracks
00492   edm::Handle<reco::TrackCollection> tracks;
00493   e.getByLabel(trackInputTag_,tracks);
00494   if(!tracks.isValid()) {
00495     return;
00496   }
00497   const reco::TrackCollection* trackCollection = tracks.product();
00498   //Photon Eta and Phi.  Hope these are correct.
00499   reco::BeamSpot vertexBeamSpot;
00500   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
00501   e.getByLabel(beamSpotProducerTag_,recoBeamSpotHandle);
00502   vertexBeamSpot = *recoBeamSpotHandle;
00503   
00504   PhotonTkIsolation phoIso(RCone, 
00505                            RinnerCone, 
00506                            etaSlice,  
00507                            pTThresh, 
00508                            lip , 
00509                            d0, 
00510                            trackCollection, 
00511                            math::XYZPoint(vertexBeamSpot.x0(),vertexBeamSpot.y0(),vertexBeamSpot.z0()));
00512 
00513   counter = phoIso.getNumberTracks(photon);
00514   ptSum = phoIso.getPtTracks(photon);
00515   ntrkCone = counter;
00516   trkCone = ptSum;
00517 }
00518 
00519 
00520 
00521 double PhotonIsolationCalculator::calculateEcalRecHitIso(const reco::Photon* photon,
00522                                                          const edm::Event& iEvent,
00523                                                          const edm::EventSetup& iSetup,
00524                                                          double RCone,
00525                                                          double RConeInner,
00526                                                          double etaSlice,
00527                                                          double eMin,
00528                                                          double etMin,
00529                                                          bool vetoClusteredHits,
00530                                                          bool useNumXtals)
00531 {
00532   
00533   edm::Handle<EcalRecHitCollection> ecalhitsCollEB;
00534   edm::Handle<EcalRecHitCollection> ecalhitsCollEE;
00535 
00536   iEvent.getByLabel(endcapecalCollection_, ecalhitsCollEE);
00537   
00538   iEvent.getByLabel(barrelecalCollection_, ecalhitsCollEB);
00539  
00540   const EcalRecHitCollection* rechitsCollectionEE_ = ecalhitsCollEE.product();
00541   const EcalRecHitCollection* rechitsCollectionEB_ = ecalhitsCollEB.product();
00542 
00543   //Get the channel status from the db
00544   edm::ESHandle<EcalChannelStatus> chStatus;
00545   iSetup.get<EcalChannelStatusRcd>().get(chStatus);
00546 
00547   edm::ESHandle<EcalSeverityLevelAlgo> sevlv;
00548   iSetup.get<EcalSeverityLevelAlgoRcd>().get(sevlv);
00549   const EcalSeverityLevelAlgo* sevLevel = sevlv.product();
00550 
00551   std::auto_ptr<CaloRecHitMetaCollectionV> RecHitsEE(0); 
00552   RecHitsEE = std::auto_ptr<CaloRecHitMetaCollectionV>(new EcalRecHitMetaCollection(*rechitsCollectionEE_));
00553  
00554   std::auto_ptr<CaloRecHitMetaCollectionV> RecHitsEB(0); 
00555   RecHitsEB = std::auto_ptr<CaloRecHitMetaCollectionV>(new EcalRecHitMetaCollection(*rechitsCollectionEB_));
00556 
00557   edm::ESHandle<CaloGeometry> geoHandle;
00558   iSetup.get<CaloGeometryRecord>().get(geoHandle);
00559 
00560   EgammaRecHitIsolation phoIsoEB(RCone,
00561                                  RConeInner,
00562                                  etaSlice,
00563                                  etMin,
00564                                  eMin,
00565                                  geoHandle,
00566                                  &(*RecHitsEB),
00567                                  sevLevel,
00568                                  DetId::Ecal);
00569 
00570   phoIsoEB.setVetoClustered(vetoClusteredHits);
00571   phoIsoEB.setUseNumCrystals(useNumXtals);
00572   phoIsoEB.doSpikeRemoval(ecalhitsCollEB.product(),chStatus.product(),severityLevelCut_);//,severityRecHitThreshold_,spId_,spikeIdThreshold_);
00573   phoIsoEB.doFlagChecks(v_chstatus_);
00574   double ecalIsolEB = phoIsoEB.getEtSum(photon);
00575   
00576   EgammaRecHitIsolation phoIsoEE(RCone,
00577                                  RConeInner,
00578                                  etaSlice,
00579                                  etMin,
00580                                  eMin,
00581                                  geoHandle,
00582                                  &(*RecHitsEE),
00583                                  sevLevel,
00584                                  DetId::Ecal);
00585   
00586   phoIsoEE.setVetoClustered(vetoClusteredHits);
00587   phoIsoEE.setUseNumCrystals(useNumXtals);
00588   phoIsoEE.doFlagChecks(v_chstatus_);
00589 
00590   double ecalIsolEE = phoIsoEE.getEtSum(photon);
00591   //  delete phoIso;
00592   double ecalIsol = ecalIsolEB + ecalIsolEE;
00593   
00594   return ecalIsol;
00595   
00596 
00597 }
00598 
00599 double PhotonIsolationCalculator::calculateHcalTowerIso(const reco::Photon* photon,
00600                                                         const edm::Event& iEvent,
00601                                                         const edm::EventSetup& iSetup,
00602                                                         double RCone,
00603                                                         double RConeInner,
00604                                                         double eMin,
00605                                                         signed int depth )
00606 {
00607 
00608   edm::Handle<CaloTowerCollection> hcalhitsCollH;
00609  
00610   iEvent.getByLabel(hcalCollection_, hcalhitsCollH);
00611   
00612   const CaloTowerCollection *toww = hcalhitsCollH.product();
00613 
00614   double ecalIsol=0.;
00615   
00616   //std::cout << "before iso call" << std::endl;
00617   EgammaTowerIsolation phoIso(RCone,
00618                               RConeInner,
00619                               eMin,depth,
00620                               toww);
00621   ecalIsol = phoIso.getTowerEtSum(photon);
00622   //  delete phoIso;
00623   //std::cout << "after call" << std::endl;
00624   return ecalIsol;
00625   
00626 
00627 }
00628 
00629 
00630 
00631
/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/RecoEgamma/PhotonIdentification/src/PhotonIsolationCalculator.cc
