Reference Manual

Go to the documentation of this file.
00001 //*****************************************************************************
00002 // File:      EgammaRecHitIsolation.cc
00003 // ----------------------------------------------------------------------------
00004 // OrigAuth:  Matthias Mozer, hacked by Sam Harper (ie the ugly stuff is mine)
00005 // Institute: IIHE-VUB, RAL
00006 //=============================================================================
00007 //*****************************************************************************
00008 //C++ includes
00009 #include <vector>
00010 #include <functional>
00011 
00012 //ROOT includes
00013 #include <Math/VectorUtil.h>
00014 
00015 //CMSSW includes
00016 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaRecHitIsolation.h"
00017 #include "DataFormats/DetId/interface/DetId.h"
00018 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
00019 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00020 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
00021 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
00022 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
00023 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
00024 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
00025 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00026 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
00027 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
00028 #include "DataFormats/Math/interface/deltaPhi.h"
00029 
00030 using namespace std;
00031 
00032 EgammaRecHitIsolation::EgammaRecHitIsolation (double extRadius,
00033                                               double intRadius,
00034                                               double etaSlice,
00035                                               double etLow,
00036                                               double eLow,
00037                                               edm::ESHandle<CaloGeometry> theCaloGeom,
00038                                               CaloRecHitMetaCollectionV* caloHits,
00039                                               const EcalSeverityLevelAlgo* sl,
00040                                               DetId::Detector detector):  // not used anymore, kept for compatibility
00041     extRadius_(extRadius),
00042     intRadius_(intRadius),
00043     etaSlice_(etaSlice),
00044     etLow_(etLow),
00045     eLow_(eLow),
00046     theCaloGeom_(theCaloGeom) ,  
00047     caloHits_(caloHits),
00048     sevLevel_(sl),
00049     useNumCrystals_(false),
00050     vetoClustered_(false),
00051     ecalBarHits_(0),
00052     chStatus_(0),
00053     severityLevelCut_(-1),
00054     //severityRecHitThreshold_(0),
00055     //spId_(EcalSeverityLevelAlgo::kSwissCross),
00056     //spIdThreshold_(0),
00057     v_chstatus_(0)
00058 {
00059     //set up the geometry and selector
00060     const CaloGeometry* caloGeom = theCaloGeom_.product();
00061     subdet_[0] = caloGeom->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
00062     subdet_[1] = caloGeom->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);
00063 
00064 }
00065 
00066 EgammaRecHitIsolation::~EgammaRecHitIsolation ()
00067 {
00068 }
00069 
00070 double EgammaRecHitIsolation::getSum_(const reco::Candidate* emObject,bool returnEt) const
00071 {
00072 
00073     double energySum = 0.;
00074     if (caloHits_){
00075         //Take the SC position
00076         reco::SuperClusterRef sc = emObject->get<reco::SuperClusterRef>();
00077         math::XYZPoint theCaloPosition = sc.get()->position();
00078         GlobalPoint pclu (theCaloPosition.x () ,
00079                 theCaloPosition.y () ,
00080                 theCaloPosition.z () );
00081         double etaclus = pclu.eta();
00082         double phiclus = pclu.phi();
00083         double r2 = intRadius_*intRadius_;
00084 
00085 
00086 
00087         std::vector< std::pair<DetId, float> >::const_iterator rhIt;
00088 
00089 
00090         for(int subdetnr=0; subdetnr<=1 ; subdetnr++){  // look in barrel and endcap
00091             CaloSubdetectorGeometry::DetIdSet chosen = subdet_[subdetnr]->getCells(pclu,extRadius_);// select cells around cluster
00092             CaloRecHitMetaCollectionV::const_iterator j=caloHits_->end();
00093             for (CaloSubdetectorGeometry::DetIdSet::const_iterator  i = chosen.begin ();i!= chosen.end ();++i){//loop selected cells
00094 
00095                 j=caloHits_->find(*i); // find selected cell among rechits
00096                 if( j!=caloHits_->end()){ // add rechit only if available 
00097                     const  GlobalPoint & position = theCaloGeom_.product()->getPosition(*i);
00098                     double eta = position.eta();
00099                     double phi = position.phi();
00100                     double etaDiff = eta - etaclus;
00101                     double phiDiff= reco::deltaPhi(phi,phiclus);
00102                     double energy = j->energy();
00103 
00104                     if(useNumCrystals_) {
00105                         if( fabs(etaclus) < 1.479 ) {  // Barrel num crystals, crystal width = 0.0174
00106                             if ( fabs(etaDiff) < 0.0174*etaSlice_) continue;  
00107                             if ( sqrt(etaDiff*etaDiff + phiDiff*phiDiff) < 0.0174*intRadius_) continue; 
00108                         } else {                       // Endcap num crystals, crystal width = 0.00864*fabs(sinh(eta))
00109                             if ( fabs(etaDiff) < 0.00864*fabs(sinh(eta))*etaSlice_) continue;  
00110                             if ( sqrt(etaDiff*etaDiff + phiDiff*phiDiff) < 0.00864*fabs(sinh(eta))*intRadius_) continue; 
00111                         }
00112                     } else {
00113                         if ( fabs(etaDiff) < etaSlice_) continue;  // jurassic strip cut
00114                         if ( etaDiff*etaDiff + phiDiff*phiDiff < r2) continue; // jurassic exclusion cone cut
00115                     }
00116 
00117                     //Check if RecHit is in SC
00118                     if(vetoClustered_) {
00119 
00120                         //Loop over basic clusters:
00121                         bool isClustered = false;
00122                         for(reco::CaloCluster_iterator bcIt = sc->clustersBegin();bcIt != sc->clustersEnd(); ++bcIt) {
00123                             for(rhIt = (*bcIt)->hitsAndFractions().begin();rhIt != (*bcIt)->hitsAndFractions().end(); ++rhIt) {
00124                                 if( rhIt->first == *i ) isClustered = true;
00125                                 if( isClustered ) break;
00126                             }
00127                             if( isClustered ) break;
00128                         } //end loop over basic clusters
00129 
00130                         if(isClustered) continue;
00131                     }  //end if removeClustered
00132 
00133                     //Severity level check
00134                     //make sure we have a barrel rechit                                     
00135                     //call the severity level method                                        
00136                     //passing the EBDetId                                                   
00137                     //the rechit collection in order to calculate the swiss crss            
00138                     //and the EcalChannelRecHitRcd                                          
00139                     //only consider rechits with ET >                                       
00140                     //the SpikeId method (currently kE1OverE9 or kSwissCross)               
00141                     //cut value for above                                                   
00142                     //then if the severity level is too high, we continue to the next rechit
00143 
00144                     if( severityLevelCut_!=-1 && ecalBarHits_ && 
00145                         sevLevel_->severityLevel(EBDetId(j->detid()), *ecalBarHits_) >= severityLevelCut_) 
00146                       continue;                    
00147                     //                            *chStatus_,                           
00148                     //        severityRecHitThreshold_,             
00149                     //        spId_,                                
00150                     //        spIdThreshold_                        
00151                     //    ) >= severityLevelCut_) continue;         
00152 
00153 
00154 
00155 
00156                     //Check based on flags to protect from recovered channels from non-read towers
00157                     //Assumption is that v_chstatus_ is empty unless doFlagChecks() has been called
00158                     std::vector<int>::const_iterator vit = std::find( v_chstatus_.begin(), v_chstatus_.end(),  ((EcalRecHit*)(&*j))->recoFlag() );
00159                     if ( vit != v_chstatus_.end() ) continue; // the recHit has to be excluded from the iso sum
00160 
00161 
00162                     double et = energy*position.perp()/position.mag();
00163                     if ( fabs(et) > etLow_ && fabs(energy) > eLow_){ //Changed energy --> fabs(energy)
00164                         if(returnEt) energySum+=et;
00165                         else energySum+=energy;
00166                     }
00167 
00168                 }                //End if not end of list
00169             }                    //End loop over rechits
00170         }                        //End loop over barrel/endcap
00171     }                            //End if caloHits_
00172     return energySum;
00173 }
00174 
00175 
00176 
00177 double EgammaRecHitIsolation::getSum_(const reco::SuperCluster* sc,bool returnEt) const
00178 {
00179 
00180     double energySum = 0.;
00181     if (caloHits_){
00182         //Take the SC position
00183  
00184         math::XYZPoint theCaloPosition = sc->position();
00185         GlobalPoint pclu (theCaloPosition.x () ,
00186                 theCaloPosition.y () ,
00187                 theCaloPosition.z () );
00188         double etaclus = pclu.eta();
00189         double phiclus = pclu.phi();
00190         double r2 = intRadius_*intRadius_;
00191 
00192         std::vector< std::pair<DetId, float> >::const_iterator rhIt;
00193 
00194 
00195         for(int subdetnr=0; subdetnr<=1 ; subdetnr++){  // look in barrel and endcap
00196             CaloSubdetectorGeometry::DetIdSet chosen = subdet_[subdetnr]->getCells(pclu,extRadius_);// select cells around cluster
00197             CaloRecHitMetaCollectionV::const_iterator j=caloHits_->end();
00198             for (CaloSubdetectorGeometry::DetIdSet::const_iterator  i = chosen.begin ();i!= chosen.end ();++i){//loop selected cells
00199 
00200                 j=caloHits_->find(*i); // find selected cell among rechits
00201                 if( j!=caloHits_->end()){ // add rechit only if available 
00202                     const  GlobalPoint & position = theCaloGeom_.product()->getPosition(*i);
00203                     double eta = position.eta();
00204                     double phi = position.phi();
00205                     double etaDiff = eta - etaclus;
00206                     double phiDiff= reco::deltaPhi(phi,phiclus);
00207                     double energy = j->energy();
00208 
00209                     if(useNumCrystals_) {
00210                         if( fabs(etaclus) < 1.479 ) {  // Barrel num crystals, crystal width = 0.0174
00211                             if ( fabs(etaDiff) < 0.0174*etaSlice_) continue;  
00212                             if ( sqrt(etaDiff*etaDiff + phiDiff*phiDiff) < 0.0174*intRadius_) continue; 
00213                         } else {                       // Endcap num crystals, crystal width = 0.00864*fabs(sinh(eta))
00214                             if ( fabs(etaDiff) < 0.00864*fabs(sinh(eta))*etaSlice_) continue;  
00215                             if ( sqrt(etaDiff*etaDiff + phiDiff*phiDiff) < 0.00864*fabs(sinh(eta))*intRadius_) continue; 
00216                         }
00217                     } else {
00218                         if ( fabs(etaDiff) < etaSlice_) continue;  // jurassic strip cut
00219                         if ( etaDiff*etaDiff + phiDiff*phiDiff < r2) continue; // jurassic exclusion cone cut
00220                     }
00221 
00222                     //Check if RecHit is in SC
00223                     if(vetoClustered_) {
00224 
00225                         //Loop over basic clusters:
00226                         bool isClustered = false;
00227                         for(reco::CaloCluster_iterator bcIt = sc->clustersBegin();bcIt != sc->clustersEnd(); ++bcIt) {
00228                             for(rhIt = (*bcIt)->hitsAndFractions().begin();rhIt != (*bcIt)->hitsAndFractions().end(); ++rhIt) {
00229                                 if( rhIt->first == *i ) isClustered = true;
00230                                 if( isClustered ) break;
00231                             }
00232                             if( isClustered ) break;
00233                         } //end loop over basic clusters
00234 
00235                         if(isClustered) continue;
00236                     }  //end if removeClustered
00237 
00238                     //Severity level check
00239                     //make sure we have a barrel rechit                                     
00240                     //call the severity level method                                        
00241                     //passing the EBDetId                                                   
00242                     //the rechit collection in order to calculate the swiss crss            
00243                     //and the EcalChannelRecHitRcd                                          
00244                     //only consider rechits with ET >                                       
00245                     //the SpikeId method (currently kE1OverE9 or kSwissCross)               
00246                     //cut value for above                                                   
00247                     //then if the severity level is too high, we continue to the next rechit
00248 
00249                     if( severityLevelCut_!=-1 && ecalBarHits_ && 
00250                         sevLevel_->severityLevel(EBDetId(j->detid()), *ecalBarHits_) >= severityLevelCut_) 
00251                       continue;                    
00252                     //                            *chStatus_,                           
00253                     //        severityRecHitThreshold_,             
00254                     //        spId_,                                
00255                     //        spIdThreshold_                        
00256                     //    ) >= severityLevelCut_) continue;         
00257 
00258 
00259 
00260 
00261                     //Check based on flags to protect from recovered channels from non-read towers
00262                     //Assumption is that v_chstatus_ is empty unless doFlagChecks() has been called
00263                     std::vector<int>::const_iterator vit = std::find( v_chstatus_.begin(), v_chstatus_.end(),  ((EcalRecHit*)(&*j))->recoFlag() );
00264                     if ( vit != v_chstatus_.end() ) continue; // the recHit has to be excluded from the iso sum
00265 
00266 
00267                     double et = energy*position.perp()/position.mag();
00268                     if ( fabs(et) > etLow_ && fabs(energy) > eLow_){ //Changed energy --> fabs(energy)
00269                         if(returnEt) energySum+=et;
00270                         else energySum+=energy;
00271                     }
00272 
00273                 }                //End if not end of list
00274             }                    //End loop over rechits
00275         }                        //End loop over barrel/endcap
00276     }                            //End if caloHits_
00277     return energySum;
00278 }
00279
/data/refman/pasoursint/CMSSW_5_2_7_hltpatch1/src/RecoEgamma/EgammaIsolationAlgos/src/EgammaRecHitIsolation.cc
