#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 >, CaloRecHitMetaCollectionV , const EcalSeverityLevelAlgo , DetId::Detector detector)
double	getEnergySum (const reco::SuperCluster *emObject) const
double	getEnergySum (const reco::Candidate *emObject) const
double	getEtSum (const reco::SuperCluster *emObject) const
double	getEtSum (const reco::Candidate *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
CaloRecHitMetaCollectionV *	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 27 of file EgammaRecHitIsolation.h.

Constructor & Destructor Documentation

EgammaRecHitIsolation::EgammaRecHitIsolation	(	double	extRadius,
		double	intRadius,
		double	etaSlice,
		double	etLow,
		double	eLow,
		edm::ESHandle< CaloGeometry >	theCaloGeom,
		CaloRecHitMetaCollectionV *	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.

{}

Member Function Documentation

void EgammaRecHitIsolation::doFlagChecks ( const std::vector< int > v ) [inline]

Definition at line 57 of file EgammaRecHitIsolation.h.

References flags_, and python::multivaluedict::sort().

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 49 of file EgammaRecHitIsolation.h.

References ecalBarHits_, severitiesexcl_, and python::multivaluedict::sort().

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::SuperCluster * emObject ) const [inline]

Definition at line 45 of file EgammaRecHitIsolation.h.

References getSum_().

{ return  getSum_(emObject,false);}

double EgammaRecHitIsolation::getEnergySum ( const reco::Candidate * emObject ) const [inline]

Definition at line 42 of file EgammaRecHitIsolation.h.

References getSum_().

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

{ return  getSum_(emObject,false);}

double EgammaRecHitIsolation::getEtSum ( const reco::Candidate * emObject ) const [inline]

Definition at line 41 of file EgammaRecHitIsolation.h.

References getSum_().

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

{return getSum_(emObject,true);}

double EgammaRecHitIsolation::getEtSum ( const reco::SuperCluster * emObject ) const [inline]

Definition at line 44 of file EgammaRecHitIsolation.h.

References getSum_().

{return getSum_(emObject,true);}

double EgammaRecHitIsolation::getSum_	(	const reco::SuperCluster *	sc,
		bool	returnEt
	)		const `[private]`

Definition at line 183 of file EgammaRecHitIsolation.cc.

References caloHits_, reco::SuperCluster::clustersBegin(), reco::SuperCluster::clustersEnd(), SiPixelRawToDigiRegional_cfi::deltaPhi, ecalBarHits_, eLow_, CaloRecHitMetaCollectionV::end(), relval_parameters_module::energy, CastorDataFrameFilter_impl::energySum(), eta, PV3DBase< T, PVType, FrameType >::eta(), etaSlice_, etLow_, extRadius_, CaloRecHitMetaCollectionV::find(), spr::find(), flags_, CaloSubdetectorGeometry::getCells(), i, intRadius_, j, EcalRecHit::kGood, PV3DBase< T, PVType, FrameType >::mag(), PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), 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_){
    //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
      CaloSubdetectorGeometry::DetIdSet chosen = subdet_[subdetnr]->getCells(pclu,extRadius_);// select cells around cluster
      CaloRecHitMetaCollectionV::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;
}

double EgammaRecHitIsolation::getSum_	(	const reco::Candidate *	emObject,
		bool	returnEt
	)		const `[private]`

Definition at line 69 of file EgammaRecHitIsolation.cc.

References caloHits_, SiPixelRawToDigiRegional_cfi::deltaPhi, cond::rpcobgas::detid, ecalBarHits_, eLow_, CaloRecHitMetaCollectionV::end(), relval_parameters_module::energy, CastorDataFrameFilter_impl::energySum(), eta, PV3DBase< T, PVType, FrameType >::eta(), etaSlice_, etLow_, extRadius_, CaloRecHitMetaCollectionV::find(), spr::find(), flags_, reco::Candidate::get(), edm::Ref< C, T, F >::get(), CaloSubdetectorGeometry::getCells(), i, intRadius_, j, EcalRecHit::kGood, phi, edm::ESHandle< T >::product(), diffTwoXMLs::r2, severitiesexcl_, EcalSeverityLevelAlgo::severityLevel(), sevLevel_, mathSSE::sqrt(), subdet_, theCaloGeom_, useNumCrystals_, and vetoClustered_.

Referenced by getEnergySum(), and getEtSum().

                                                                                       {
  
  double energySum = 0.;
  if (caloHits_){
    //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
      CaloSubdetectorGeometry::DetIdSet chosen = subdet_[subdetnr]->getCells(pclu,extRadius_);// select cells around cluster
      CaloRecHitMetaCollectionV::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;
}