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