#include <EgammaRecHitIsolation.h>

Public Member Functions
void	doFlagChecks (const std::vector< int > &v)

void	doSeverityChecks (const EcalRecHitCollection *const recHits, const std::vector< int > &v)

	EgammaRecHitIsolation (double extRadius, double intRadius, double etaSlice, double etLow, double eLow, edm::ESHandle< CaloGeometry >, const EcalRecHitCollection &, const EcalSeverityLevelAlgo *, DetId::Detector detector)

double	getEnergySum (const reco::Candidate *emObject) const

double	getEnergySum (const reco::SuperCluster *emObject) const

double	getEtSum (const reco::Candidate *emObject) const

double	getEtSum (const reco::SuperCluster *emObject) const

void	setUseNumCrystals (bool b=true)

void	setVetoClustered (bool b=true)

	~EgammaRecHitIsolation ()

Private Member Functions
double	getSum_ (const reco::Candidate *, bool returnEt) const

double	getSum_ (const reco::SuperCluster *, bool returnEt) const

Private Attributes
const EcalRecHitCollection &	caloHits_

const EcalRecHitCollection *	ecalBarHits_

double	eLow_

double	etaSlice_

double	etLow_

double	extRadius_

std::vector< int >	flags_

double	intRadius_

std::vector< int >	severitiesexcl_

const EcalSeverityLevelAlgo *	sevLevel_

const CaloSubdetectorGeometry *	subdet_ [2]

edm::ESHandle< CaloGeometry >	theCaloGeom_

bool	useNumCrystals_

bool	vetoClustered_

Detailed Description

Definition at line 28 of file EgammaRecHitIsolation.h.

Constructor & Destructor Documentation

EgammaRecHitIsolation::EgammaRecHitIsolation	(	double	extRadius,
		double	intRadius,
		double	etaSlice,
		double	etLow,
		double	eLow,
		edm::ESHandle< CaloGeometry >	theCaloGeom,
		const EcalRecHitCollection &	caloHits,
		const EcalSeverityLevelAlgo *	sl,
		DetId::Detector	detector
	)

Definition at line 32 of file EgammaRecHitIsolation.cc.

References DetId::Ecal, EcalBarrel, EcalEndcap, CaloGeometry::getSubdetectorGeometry(), edm::ESHandle< T >::product(), subdet_, and theCaloGeom_.

                                                                      :  // not used anymore, kept for compatibility
     extRadius_(extRadius),
     intRadius_(intRadius),
     etaSlice_(etaSlice),
     etLow_(etLow),
     eLow_(eLow),
     theCaloGeom_(theCaloGeom) ,  
     caloHits_(caloHits),
     sevLevel_(sl),
     useNumCrystals_(false),
     vetoClustered_(false),
     ecalBarHits_(0),
     //chStatus_(0),
     severitiesexcl_(0),
     //severityRecHitThreshold_(0),
     //spId_(EcalSeverityLevelAlgo::kSwissCross),
     //spIdThreshold_(0),
     flags_(0)
 {
     //set up the geometry and selector
     const CaloGeometry* caloGeom = theCaloGeom_.product();
     subdet_[0] = caloGeom->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
     subdet_[1] = caloGeom->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);
 
 }

EgammaRecHitIsolation::~EgammaRecHitIsolation ( )

Definition at line 66 of file EgammaRecHitIsolation.cc.

Referenced by doFlagChecks().

67 {}

Member Function Documentation

void EgammaRecHitIsolation::doFlagChecks ( const std::vector< int > & v )

inline

Definition at line 58 of file EgammaRecHitIsolation.h.

References flags_, getSum_(), and ~EgammaRecHitIsolation().

Referenced by GsfElectronAlgo::beginEvent(), ConversionTrackCandidateProducer::buildCollections(), and PhotonIsolationCalculator::calculateEcalRecHitIso().

                                              {
     flags_.clear();
     flags_.insert(flags_.begin(), v.begin(), v.end());
     std::sort(flags_.begin(), flags_.end() );
   }

void EgammaRecHitIsolation::doSeverityChecks	(	const EcalRecHitCollection *const	recHits,
		const std::vector< int > &	v
	)

inline

Definition at line 50 of file EgammaRecHitIsolation.h.

References ecalBarHits_, and severitiesexcl_.

Referenced by GsfElectronAlgo::beginEvent(), ConversionTrackCandidateProducer::buildCollections(), and PhotonIsolationCalculator::calculateEcalRecHitIso().

                                      { 
     ecalBarHits_ = recHits; 
     severitiesexcl_.clear();
     severitiesexcl_.insert(severitiesexcl_.begin(), v.begin(), v.end());
     std::sort(severitiesexcl_.begin(), severitiesexcl_.end());
   }

double EgammaRecHitIsolation::getEnergySum ( const reco::Candidate * emObject ) const

inline

Definition at line 43 of file EgammaRecHitIsolation.h.

References getSum_().

Referenced by EgammaHLTEcalRecIsolationProducer::produce(), and EgammaEcalRecHitIsolationProducer::produce().

43 { return getSum_(emObject,false);}

EgammaRecHitIsolation::getSum_

double getSum_(const reco::Candidate *, bool returnEt) const

Definition: EgammaRecHitIsolation.cc:69

double EgammaRecHitIsolation::getEnergySum ( const reco::SuperCluster * emObject ) const

inline

Definition at line 46 of file EgammaRecHitIsolation.h.

References getSum_().

46 { return getSum_(emObject,false);}

EgammaRecHitIsolation::getSum_

double getSum_(const reco::Candidate *, bool returnEt) const

Definition: EgammaRecHitIsolation.cc:69

double EgammaRecHitIsolation::getEtSum ( const reco::Candidate * emObject ) const

inline

Definition at line 42 of file EgammaRecHitIsolation.h.

References getSum_().

Referenced by ConversionTrackCandidateProducer::buildCollections(), PhotonIsolationCalculator::calculateEcalRecHitIso(), GsfElectronAlgo::createElectron(), EgammaHLTEcalRecIsolationProducer::produce(), and EgammaEcalRecHitIsolationProducer::produce().

42 {return getSum_(emObject,true);}

EgammaRecHitIsolation::getSum_

double getSum_(const reco::Candidate *, bool returnEt) const

Definition: EgammaRecHitIsolation.cc:69

double EgammaRecHitIsolation::getEtSum ( const reco::SuperCluster * emObject ) const

inline

Definition at line 45 of file EgammaRecHitIsolation.h.

References getSum_().

45 {return getSum_(emObject,true);}

EgammaRecHitIsolation::getSum_

double getSum_(const reco::Candidate *, bool returnEt) const

Definition: EgammaRecHitIsolation.cc:69

double EgammaRecHitIsolation::getSum_	(	const reco::Candidate *	emObject,
		bool	returnEt
	)		const

private

Definition at line 69 of file EgammaRecHitIsolation.cc.

References caloHits_, reco::deltaPhi(), ecalBarHits_, eLow_, edm::SortedCollection< T, SORT >::empty(), edm::SortedCollection< T, SORT >::end(), CastorDataFrameFilter_impl::energySum(), stringResolutionProvider_cfi::et, PVValHelper::eta, PV3DBase< T, PVType, FrameType >::eta(), etaSlice_, etLow_, extRadius_, spr::find(), edm::SortedCollection< T, SORT >::find(), flags_, reco::Candidate::get(), edm::Ref< C, T, F >::get(), CaloSubdetectorGeometry::getCells(), CaloGeometry::getGeometry(), mps_fire::i, intRadius_, EcalRecHit::kGood, phi, edm::ESHandle< T >::product(), diffTwoXMLs::r2, SimDataFormats::CaloAnalysis::sc, severitiesexcl_, EcalSeverityLevelAlgo::severityLevel(), sevLevel_, mathSSE::sqrt(), subdet_, theCaloGeom_, useNumCrystals_, and vetoClustered_.

Referenced by doFlagChecks(), getEnergySum(), and getEtSum().

                                                                                        {
   
   double energySum = 0.;
   if (! caloHits_.empty()) {
     //Take the SC position
     reco::SuperClusterRef sc = emObject->get<reco::SuperClusterRef>();
     math::XYZPoint const & theCaloPosition = sc.get()->position();
     GlobalPoint pclu (theCaloPosition.x () ,
               theCaloPosition.y () ,
               theCaloPosition.z () );
     float etaclus = pclu.eta();
     float phiclus = pclu.phi();
     float r2 = intRadius_*intRadius_;
     
     std::vector< std::pair<DetId, float> >::const_iterator rhIt;
     
     for(int subdetnr=0; subdetnr<=1 ; subdetnr++){  // look in barrel and endcap
       if( nullptr == subdet_[subdetnr] ) continue;
 
       CaloSubdetectorGeometry::DetIdSet chosen = subdet_[subdetnr]->getCells(pclu,extRadius_);// select cells around cluster
       EcalRecHitCollection::const_iterator j = caloHits_.end();
 
       for (CaloSubdetectorGeometry::DetIdSet::const_iterator  i = chosen.begin ();i != chosen.end (); ++i){ //loop selected cells
     j = caloHits_.find(*i); // find selected cell among rechits
     if(j != caloHits_.end()) { // add rechit only if available 
       auto const cell  = theCaloGeom_.product()->getGeometry(*i);
       float eta = cell->etaPos();
       float phi = cell->phiPos();
       float etaDiff = eta - etaclus;
       float phiDiff= reco::deltaPhi(phi,phiclus);
       float energy = j->energy();
 
       if(useNumCrystals_) {
         if(fabs(etaclus) < 1.479) {  // Barrel num crystals, crystal width = 0.0174
           if (fabs(etaDiff) < 0.0174*etaSlice_) 
         continue;  
           //if (sqrt(etaDiff*etaDiff + phiDiff*phiDiff) < 0.0174*intRadius_) 
           //continue; 
           if ((etaDiff*etaDiff + phiDiff*phiDiff) < 0.00030276*r2) 
         continue; 
         } else {                       // Endcap num crystals, crystal width = 0.00864*fabs(sinh(eta))
           if (fabs(etaDiff) < 0.00864*fabs(sinh(eta))*etaSlice_) 
         continue;  
           //if (sqrt(etaDiff*etaDiff + phiDiff*phiDiff) < 0.00864*fabs(sinh(eta))*intRadius_) 
           //    continue; 
           if ((etaDiff*etaDiff + phiDiff*phiDiff) < (0.000037325*(cosh(2*eta)-1)*r2))
             continue; 
         }
       } else {
         if (fabs(etaDiff) < etaSlice_) 
           continue;  // jurassic strip cut
         if (etaDiff*etaDiff + phiDiff*phiDiff < r2) 
           continue; // jurassic exclusion cone cut
       }
       //Check if RecHit is in SC
       if(vetoClustered_) {
         
         //Loop over basic clusters:
         bool isClustered = false;
         for(reco::CaloCluster_iterator bcIt = sc->clustersBegin();bcIt != sc->clustersEnd(); ++bcIt) {
           for(rhIt = (*bcIt)->hitsAndFractions().begin();rhIt != (*bcIt)->hitsAndFractions().end(); ++rhIt) {
         if(rhIt->first == *i)
           isClustered = true;
         if(isClustered) 
           break;
           }
           
           if(isClustered) 
         break;
         } //end loop over basic clusters
         
         if(isClustered) 
           continue;
       }  //end if removeClustered
       
       
       
       
       //std::cout << "detid " << ((EcalRecHit*)(&*j))->detid() << std::endl;
       int severityFlag = ecalBarHits_ == 0 ? -1 : sevLevel_->severityLevel(((const EcalRecHit*)(&*j))->detid(), *ecalBarHits_);
       std::vector<int>::const_iterator sit = std::find(severitiesexcl_.begin(), 
                                severitiesexcl_.end(), 
                                severityFlag);
 
       if (sit!= severitiesexcl_.end())
         continue;
       
       // new rechit flag checks
       //std::vector<int>::const_iterator vit = std::find(flags_.begin(), 
       //                           flags_.end(),
       //                           ((const EcalRecHit*)(&*j))->recoFlag());
       //if (vit != flags_.end()) 
       //  continue;
       if (!((EcalRecHit*)(&*j))->checkFlag(EcalRecHit::kGood)) {
         if (((EcalRecHit*)(&*j))->checkFlags(flags_)) {                
           continue;
         }
       }
       
       float et = energy*std::sqrt(cell->getPosition().perp2()/cell->getPosition().mag2());
       if ( et > etLow_ && energy > eLow_) { //Changed energy --> fabs(energy) - now changed back to energy
         if(returnEt) 
           energySum += et;
         else 
           energySum += energy;
       }
       
     }                //End if not end of list
       }                  //End loop over rechits
     }                    //End loop over barrel/endcap
   }                        //End if caloHits_
   
   return energySum;
 }

double EgammaRecHitIsolation::getSum_	(	const reco::SuperCluster *	sc,
		bool	returnEt
	)		const

private

Definition at line 186 of file EgammaRecHitIsolation.cc.

References caloHits_, reco::SuperCluster::clustersBegin(), reco::SuperCluster::clustersEnd(), reco::deltaPhi(), ecalBarHits_, eLow_, edm::SortedCollection< T, SORT >::empty(), edm::SortedCollection< T, SORT >::end(), CastorDataFrameFilter_impl::energySum(), stringResolutionProvider_cfi::et, PVValHelper::eta, PV3DBase< T, PVType, FrameType >::eta(), etaSlice_, etLow_, extRadius_, spr::find(), edm::SortedCollection< T, SORT >::find(), flags_, CaloSubdetectorGeometry::getCells(), CaloGeometry::getPosition(), mps_fire::i, intRadius_, EcalRecHit::kGood, PV3DBase< T, PVType, FrameType >::mag(), PV3DBase< T, PVType, FrameType >::perp(), phi, PV3DBase< T, PVType, FrameType >::phi(), reco::CaloCluster::position(), position, edm::ESHandle< T >::product(), diffTwoXMLs::r2, severitiesexcl_, EcalSeverityLevelAlgo::severityLevel(), sevLevel_, subdet_, theCaloGeom_, useNumCrystals_, and vetoClustered_.

                                                                                      {
 
   double energySum = 0.;
   if (! caloHits_.empty()){
     //Take the SC position
  
     math::XYZPoint theCaloPosition = sc->position();
     GlobalPoint pclu (theCaloPosition.x () ,
               theCaloPosition.y () ,
               theCaloPosition.z () );
     double etaclus = pclu.eta();
     double phiclus = pclu.phi();
     double r2 = intRadius_*intRadius_;
     
     std::vector< std::pair<DetId, float> >::const_iterator rhIt;
     
     
     for(int subdetnr=0; subdetnr<=1 ; subdetnr++){  // look in barrel and endcap
       if( nullptr == subdet_[subdetnr] ) continue;
       CaloSubdetectorGeometry::DetIdSet chosen = subdet_[subdetnr]->getCells(pclu,extRadius_);// select cells around cluster
       EcalRecHitCollection::const_iterator j=caloHits_.end();
       for (CaloSubdetectorGeometry::DetIdSet::const_iterator  i = chosen.begin ();i!= chosen.end ();++i){//loop selected cells
     
     j=caloHits_.find(*i); // find selected cell among rechits
     if( j!=caloHits_.end()){ // add rechit only if available 
       const  GlobalPoint & position = theCaloGeom_.product()->getPosition(*i);
       double eta = position.eta();
       double phi = position.phi();
       double etaDiff = eta - etaclus;
       double phiDiff= reco::deltaPhi(phi,phiclus);
       double energy = j->energy();
       
       if(useNumCrystals_) {
         if( fabs(etaclus) < 1.479 ) {  // Barrel num crystals, crystal width = 0.0174
           if ( fabs(etaDiff) < 0.0174*etaSlice_) continue;  
           // if ( sqrt(etaDiff*etaDiff + phiDiff*phiDiff) < 0.0174*intRadius_) continue; 
           if ((etaDiff*etaDiff + phiDiff*phiDiff) < 0.00030276*r2) continue;
         } else {                       // Endcap num crystals, crystal width = 0.00864*fabs(sinh(eta))
           if ( fabs(etaDiff) < 0.00864*fabs(sinh(eta))*etaSlice_) continue;  
           // if ( sqrt(etaDiff*etaDiff + phiDiff*phiDiff) < 0.00864*fabs(sinh(eta))*intRadius_) continue; 
           if ((etaDiff*etaDiff + phiDiff*phiDiff) < (0.000037325*(cosh(2*eta)-1)*r2)) continue;
         }
       } else {
         if ( fabs(etaDiff) < etaSlice_) continue;  // jurassic strip cut
         if ( etaDiff*etaDiff + phiDiff*phiDiff < r2) continue; // jurassic exclusion cone cut
       }
       
       //Check if RecHit is in SC
       if(vetoClustered_) {
         
         //Loop over basic clusters:
         bool isClustered = false;
         for(reco::CaloCluster_iterator bcIt = sc->clustersBegin();bcIt != sc->clustersEnd(); ++bcIt) {
           for(rhIt = (*bcIt)->hitsAndFractions().begin();rhIt != (*bcIt)->hitsAndFractions().end(); ++rhIt) {
         if( rhIt->first == *i ) isClustered = true;
         if( isClustered ) break;
           }
           if( isClustered ) break;
         } //end loop over basic clusters
         
         if(isClustered) continue;
       }  //end if removeClustered
       
       
       int severityFlag = sevLevel_->severityLevel(j->detid(), *ecalBarHits_);
       std::vector<int>::const_iterator sit = std::find(severitiesexcl_.begin(), 
                                severitiesexcl_.end(), 
                                severityFlag);
       
       if (sit!= severitiesexcl_.end())
         continue;
       
       // new rechit flag checks
       //std::vector<int>::const_iterator vit = std::find(flags_.begin(), 
       //                           flags_.end(),
       //                           ((EcalRecHit*)(&*j))->recoFlag());
       //if (vit != flags_.end()) 
       //  continue;
       if (!((EcalRecHit*)(&*j))->checkFlag(EcalRecHit::kGood)) {
         if (((EcalRecHit*)(&*j))->checkFlags(flags_)) {                
           continue;
         }
       }
       
       
       double et = energy*position.perp()/position.mag();
       if ( et > etLow_ && energy > eLow_){ //Changed energy --> fabs(energy) -- then changed into energy
         if(returnEt) energySum+=et;
         else energySum+=energy;
       }
       
     }                //End if not end of list
       }                  //End loop over rechits
     }                    //End loop over barrel/endcap
   }                      //End if caloHits_
 
   return energySum;
 }