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     severitiesexcl_(0),
00054     //severityRecHitThreshold_(0),
00055     //spId_(EcalSeverityLevelAlgo::kSwissCross),
00056     //spIdThreshold_(0),
00057     flags_(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 double EgammaRecHitIsolation::getSum_(const reco::Candidate* emObject,bool returnEt) const {
00070   
00071   double energySum = 0.;
00072   if (caloHits_){
00073     //Take the SC position
00074     reco::SuperClusterRef sc = emObject->get<reco::SuperClusterRef>();
00075     math::XYZPoint const & theCaloPosition = sc.get()->position();
00076     GlobalPoint pclu (theCaloPosition.x () ,
00077                       theCaloPosition.y () ,
00078                       theCaloPosition.z () );
00079     float etaclus = pclu.eta();
00080     float phiclus = pclu.phi();
00081     float r2 = intRadius_*intRadius_;
00082     
00083     std::vector< std::pair<DetId, float> >::const_iterator rhIt;
00084     
00085     for(int subdetnr=0; subdetnr<=1 ; subdetnr++){  // look in barrel and endcap
00086       CaloSubdetectorGeometry::DetIdSet chosen = subdet_[subdetnr]->getCells(pclu,extRadius_);// select cells around cluster
00087       CaloRecHitMetaCollectionV::const_iterator j = caloHits_->end();
00088       for (CaloSubdetectorGeometry::DetIdSet::const_iterator  i = chosen.begin ();i != chosen.end (); ++i){ //loop selected cells
00089         j = caloHits_->find(*i); // find selected cell among rechits
00090         if(j != caloHits_->end()) { // add rechit only if available 
00091           auto const cell  = theCaloGeom_.product()->getGeometry(*i);
00092           float eta = cell->etaPos();
00093           float phi = cell->phiPos();
00094           float etaDiff = eta - etaclus;
00095           float phiDiff= reco::deltaPhi(phi,phiclus);
00096           float energy = j->energy();
00097 
00098           if(useNumCrystals_) {
00099             if(fabs(etaclus) < 1.479) {  // Barrel num crystals, crystal width = 0.0174
00100               if (fabs(etaDiff) < 0.0174*etaSlice_) 
00101                 continue;  
00102               //if (sqrt(etaDiff*etaDiff + phiDiff*phiDiff) < 0.0174*intRadius_) 
00103               //continue; 
00104               if ((etaDiff*etaDiff + phiDiff*phiDiff) < 0.00030276*r2) 
00105                 continue; 
00106             } else {                       // Endcap num crystals, crystal width = 0.00864*fabs(sinh(eta))
00107               if (fabs(etaDiff) < 0.00864*fabs(sinh(eta))*etaSlice_) 
00108                 continue;  
00109               //if (sqrt(etaDiff*etaDiff + phiDiff*phiDiff) < 0.00864*fabs(sinh(eta))*intRadius_) 
00110               //        continue; 
00111               if ((etaDiff*etaDiff + phiDiff*phiDiff) < (0.000037325*(cosh(2*eta)-1)*r2))
00112                 continue; 
00113             }
00114           } else {
00115             if (fabs(etaDiff) < etaSlice_) 
00116               continue;  // jurassic strip cut
00117             if (etaDiff*etaDiff + phiDiff*phiDiff < r2) 
00118               continue; // jurassic exclusion cone cut
00119           }
00120           //Check if RecHit is in SC
00121           if(vetoClustered_) {
00122             
00123             //Loop over basic clusters:
00124             bool isClustered = false;
00125             for(reco::CaloCluster_iterator bcIt = sc->clustersBegin();bcIt != sc->clustersEnd(); ++bcIt) {
00126               for(rhIt = (*bcIt)->hitsAndFractions().begin();rhIt != (*bcIt)->hitsAndFractions().end(); ++rhIt) {
00127                 if(rhIt->first == *i)
00128                   isClustered = true;
00129                 if(isClustered) 
00130                   break;
00131               }
00132               
00133               if(isClustered) 
00134                 break;
00135             } //end loop over basic clusters
00136             
00137             if(isClustered) 
00138               continue;
00139           }  //end if removeClustered
00140           
00141           
00142           
00143           
00144           //std::cout << "detid " << ((EcalRecHit*)(&*j))->detid() << std::endl;
00145           int severityFlag = ecalBarHits_ == 0 ? -1 : sevLevel_->severityLevel(((const EcalRecHit*)(&*j))->detid(), *ecalBarHits_);
00146           std::vector<int>::const_iterator sit = std::find(severitiesexcl_.begin(), 
00147                                                            severitiesexcl_.end(), 
00148                                                            severityFlag);
00149 
00150           if (sit!= severitiesexcl_.end())
00151             continue;
00152           
00153           // new rechit flag checks
00154           //std::vector<int>::const_iterator vit = std::find(flags_.begin(), 
00155           //                                               flags_.end(),
00156           //                                               ((const EcalRecHit*)(&*j))->recoFlag());
00157           //if (vit != flags_.end()) 
00158           //  continue;
00159           if (!((EcalRecHit*)(&*j))->checkFlag(EcalRecHit::kGood)) {
00160             if (((EcalRecHit*)(&*j))->checkFlags(flags_)) {                
00161               continue;
00162             }
00163           }
00164           
00165           float et = energy*std::sqrt(cell->getPosition().perp2()/cell->getPosition().mag2());
00166           if ( et > etLow_ && energy > eLow_) { //Changed energy --> fabs(energy) - now changed back to energy
00167             if(returnEt) 
00168               energySum += et;
00169             else 
00170               energySum += energy;
00171           }
00172           
00173         }                //End if not end of list
00174       }                  //End loop over rechits
00175     }                    //End loop over barrel/endcap
00176   }                        //End if caloHits_
00177   
00178   return energySum;
00179 }
00180 
00181 
00182 
00183 double EgammaRecHitIsolation::getSum_(const reco::SuperCluster* sc, bool returnEt) const {
00184 
00185   double energySum = 0.;
00186   if (caloHits_){
00187     //Take the SC position
00188  
00189     math::XYZPoint theCaloPosition = sc->position();
00190     GlobalPoint pclu (theCaloPosition.x () ,
00191                       theCaloPosition.y () ,
00192                       theCaloPosition.z () );
00193     double etaclus = pclu.eta();
00194     double phiclus = pclu.phi();
00195     double r2 = intRadius_*intRadius_;
00196     
00197     std::vector< std::pair<DetId, float> >::const_iterator rhIt;
00198     
00199     
00200     for(int subdetnr=0; subdetnr<=1 ; subdetnr++){  // look in barrel and endcap
00201       CaloSubdetectorGeometry::DetIdSet chosen = subdet_[subdetnr]->getCells(pclu,extRadius_);// select cells around cluster
00202       CaloRecHitMetaCollectionV::const_iterator j=caloHits_->end();
00203       for (CaloSubdetectorGeometry::DetIdSet::const_iterator  i = chosen.begin ();i!= chosen.end ();++i){//loop selected cells
00204         
00205         j=caloHits_->find(*i); // find selected cell among rechits
00206         if( j!=caloHits_->end()){ // add rechit only if available 
00207           const  GlobalPoint & position = theCaloGeom_.product()->getPosition(*i);
00208           double eta = position.eta();
00209           double phi = position.phi();
00210           double etaDiff = eta - etaclus;
00211           double phiDiff= reco::deltaPhi(phi,phiclus);
00212           double energy = j->energy();
00213           
00214           if(useNumCrystals_) {
00215             if( fabs(etaclus) < 1.479 ) {  // Barrel num crystals, crystal width = 0.0174
00216               if ( fabs(etaDiff) < 0.0174*etaSlice_) continue;  
00217               // if ( sqrt(etaDiff*etaDiff + phiDiff*phiDiff) < 0.0174*intRadius_) continue; 
00218               if ((etaDiff*etaDiff + phiDiff*phiDiff) < 0.00030276*r2) continue;
00219             } else {                       // Endcap num crystals, crystal width = 0.00864*fabs(sinh(eta))
00220               if ( fabs(etaDiff) < 0.00864*fabs(sinh(eta))*etaSlice_) continue;  
00221               // if ( sqrt(etaDiff*etaDiff + phiDiff*phiDiff) < 0.00864*fabs(sinh(eta))*intRadius_) continue; 
00222               if ((etaDiff*etaDiff + phiDiff*phiDiff) < (0.000037325*(cosh(2*eta)-1)*r2)) continue;
00223             }
00224           } else {
00225             if ( fabs(etaDiff) < etaSlice_) continue;  // jurassic strip cut
00226             if ( etaDiff*etaDiff + phiDiff*phiDiff < r2) continue; // jurassic exclusion cone cut
00227           }
00228           
00229           //Check if RecHit is in SC
00230           if(vetoClustered_) {
00231             
00232             //Loop over basic clusters:
00233             bool isClustered = false;
00234             for(reco::CaloCluster_iterator bcIt = sc->clustersBegin();bcIt != sc->clustersEnd(); ++bcIt) {
00235               for(rhIt = (*bcIt)->hitsAndFractions().begin();rhIt != (*bcIt)->hitsAndFractions().end(); ++rhIt) {
00236                 if( rhIt->first == *i ) isClustered = true;
00237                 if( isClustered ) break;
00238               }
00239               if( isClustered ) break;
00240             } //end loop over basic clusters
00241             
00242             if(isClustered) continue;
00243           }  //end if removeClustered
00244           
00245           
00246           int severityFlag = sevLevel_->severityLevel(j->detid(), *ecalBarHits_);
00247           std::vector<int>::const_iterator sit = std::find(severitiesexcl_.begin(), 
00248                                                            severitiesexcl_.end(), 
00249                                                            severityFlag);
00250           
00251           if (sit!= severitiesexcl_.end())
00252             continue;
00253           
00254           // new rechit flag checks
00255           //std::vector<int>::const_iterator vit = std::find(flags_.begin(), 
00256           //                                               flags_.end(),
00257           //                                               ((EcalRecHit*)(&*j))->recoFlag());
00258           //if (vit != flags_.end()) 
00259           //  continue;
00260           if (!((EcalRecHit*)(&*j))->checkFlag(EcalRecHit::kGood)) {
00261             if (((EcalRecHit*)(&*j))->checkFlags(flags_)) {                
00262               continue;
00263             }
00264           }
00265           
00266           
00267           double et = energy*position.perp()/position.mag();
00268           if ( et > etLow_ && energy > eLow_){ //Changed energy --> fabs(energy) -- then changed into 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 
00278   return energySum;
00279 }
00280
/data/refman/pasoursint/CMSSW_6_1_2_SLHC2/src/RecoEgamma/EgammaIsolationAlgos/src/EgammaRecHitIsolation.cc
