#include <EgammaRecHitExtractor.h>

Inheritance diagram for egammaisolation::EgammaRecHitExtractor:

Public Member Functions
virtual reco::IsoDeposit	deposit (const edm::Event &ev, const edm::EventSetup &evSetup, const reco::Track &track) const
virtual reco::IsoDeposit	deposit (const edm::Event &ev, const edm::EventSetup &evSetup, const reco::Candidate &c) const
	EgammaRecHitExtractor (const edm::ParameterSet &par)
virtual void	fillVetos (const edm::Event &ev, const edm::EventSetup &evSetup, const reco::TrackCollection &tracks)
virtual	~EgammaRecHitExtractor ()
Private Member Functions
void	collect (reco::IsoDeposit &deposit, const reco::SuperClusterRef &sc, const CaloSubdetectorGeometry subdet, const CaloGeometry caloGeom, const EcalRecHitCollection &hits, const EcalChannelStatus *chStatus, bool barrel) const
Private Attributes
edm::InputTag	barrelEcalHitsTag_
edm::InputTag	endcapEcalHitsTag_
double	energyMin_
double	etMin_
double	extRadius_
bool	fakeNegativeDeposit_
double	intRadius_
double	intStrip_
bool	sameTag_
int	severityLevelCut_
float	severityRecHitThreshold_
EcalSeverityLevelAlgo::SpikeId	spId_
std::string	spIdString_
float	spIdThreshold_
bool	tryBoth_
bool	useEt_
std::vector< int >	v_chstatus_
bool	vetoClustered_

Detailed Description

Definition at line 39 of file EgammaRecHitExtractor.h.

Constructor & Destructor Documentation

EgammaRecHitExtractor::EgammaRecHitExtractor ( const edm::ParameterSet & par )

Definition at line 38 of file EgammaRecHitExtractor.cc.

References barrelEcalHitsTag_, edm::InputTag::encode(), endcapEcalHitsTag_, Exception, fakeNegativeDeposit_, edm::ParameterSet::getParameter(), intRadius_, EcalSeverityLevelAlgo::kE1OverE9, EcalSeverityLevelAlgo::kSwissCross, EcalSeverityLevelAlgo::kSwissCrossBordersIncluded, sameTag_, spId_, spIdString_, tryBoth_, and useEt_.

                                                                       : 
    etMin_(par.getParameter<double>("etMin")),
    energyMin_(par.getParameter<double>("energyMin")),
    extRadius_(par.getParameter<double>("extRadius")),
    intRadius_(par.getParameter<double>("intRadius")),
    intStrip_(par.getParameter<double>("intStrip")),
    barrelEcalHitsTag_(par.getParameter<edm::InputTag>("barrelEcalHits")), 
    endcapEcalHitsTag_(par.getParameter<edm::InputTag>("endcapEcalHits")),
    fakeNegativeDeposit_(par.getParameter<bool>("subtractSuperClusterEnergy")),
    tryBoth_(par.getParameter<bool>("tryBoth")),
    vetoClustered_(par.getParameter<bool>("vetoClustered")),
    sameTag_(false),
    severityLevelCut_(par.getParameter<int>("severityLevelCut")),
    severityRecHitThreshold_(par.getParameter<double>("severityRecHitThreshold")),
    spIdString_(par.getParameter<std::string>("spikeIdString")),
    spIdThreshold_(par.getParameter<double>("spikeIdThreshold")),
    v_chstatus_(par.getParameter<std::vector<int> >("recHitFlagsToBeExcluded"))
{ 

    if     ( !spIdString_.compare("kE1OverE9") )                  spId_ = EcalSeverityLevelAlgo::kE1OverE9;
    else if( !spIdString_.compare("kSwissCross") )                spId_ = EcalSeverityLevelAlgo::kSwissCross;
    else if( !spIdString_.compare("kSwissCrossBordersIncluded") ) spId_ = EcalSeverityLevelAlgo::kSwissCrossBordersIncluded;
    else                                                          spId_ = EcalSeverityLevelAlgo::kSwissCross;

    if ((intRadius_ != 0.0) && (fakeNegativeDeposit_)) {
        throw cms::Exception("Configuration Error") << "EgammaRecHitExtractor: " << 
            "If you use 'subtractSuperClusterEnergy', you *must* set 'intRadius' to ZERO; it does not make sense, otherwise.";
    }
    std::string isoVariable = par.getParameter<std::string>("isolationVariable");
    if (isoVariable == "et") {
        useEt_ = true;
    } else if (isoVariable == "energy") {
        useEt_ = false;
    } else {
        throw cms::Exception("Configuration Error") << "EgammaRecHitExtractor: isolationVariable '" << isoVariable << "' not known. " 
            << " Supported values are 'et', 'energy'. ";
    }
    if (endcapEcalHitsTag_.encode() ==  barrelEcalHitsTag_.encode()) {
        sameTag_ = true;
        if (tryBoth_) {
            edm::LogWarning("EgammaRecHitExtractor") << "If you have configured 'barrelRecHits' == 'endcapRecHits', so I'm switching 'tryBoth' to FALSE.";
            tryBoth_ = false;
        }
    }

}

EgammaRecHitExtractor::~EgammaRecHitExtractor ( ) [virtual]

Definition at line 85 of file EgammaRecHitExtractor.cc.

{ }

Member Function Documentation

void EgammaRecHitExtractor::collect	(	reco::IsoDeposit &	deposit,
		const reco::SuperClusterRef &	sc,
		const CaloSubdetectorGeometry *	subdet,
		const CaloGeometry *	caloGeom,
		const EcalRecHitCollection &	hits,
		const EcalChannelStatus *	chStatus,
		bool	barrel
	)		const `[private]`

Definition at line 141 of file EgammaRecHitExtractor.cc.

References reco::IsoDeposit::addDeposit(), Geom::deltaPhi(), end, edm::SortedCollection< T, SORT >::end(), relval_parameters_module::energy, energyMin_, ExpressReco_HICollisions_FallBack::et, eta(), PV3DBase< T, PVType, FrameType >::eta(), etMin_, extRadius_, edm::SortedCollection< T, SORT >::find(), spr::find(), CaloSubdetectorGeometry::getCells(), CaloGeometry::getPosition(), i, intRadius_, intStrip_, j, PV3DBase< T, PVType, FrameType >::mag(), PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), phi, position, EcalSeverityLevelAlgo::severityLevel(), severityLevelCut_, severityRecHitThreshold_, spId_, spIdThreshold_, useEt_, v_chstatus_, and vetoClustered_.

Referenced by deposit().

{

    GlobalPoint caloPosition(sc->position().x(), sc->position().y() , sc->position().z());
    CaloSubdetectorGeometry::DetIdSet chosen = subdet->getCells(caloPosition,extRadius_);
    EcalRecHitCollection::const_iterator j=hits.end();
    double caloeta=caloPosition.eta();
    double calophi=caloPosition.phi();
    double r2 = intRadius_*intRadius_;

    std::vector< std::pair<DetId, float> >::const_iterator rhIt;


    for (CaloSubdetectorGeometry::DetIdSet::const_iterator i = chosen.begin(), end = chosen.end() ; i != end;  ++i)  {  
        j=hits.find(*i);
        if(j != hits.end()){
            const  GlobalPoint & position = caloGeom->getPosition(*i);
            double eta = position.eta();
            double phi = position.phi();
            double energy = j->energy();
            double et = energy*position.perp()/position.mag();
            double phiDiff= reco::deltaPhi(phi,calophi);

            //check if we are supposed to veto clustered and then do so
            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


            if(barrel &&                                 //make sure we have a barrel rechit
//             if(//j->id().subdetId() == 1 &&                //make sure we have a barrel rechit
               EcalSeverityLevelAlgo::severityLevel(     //call the severity level method
                   EBDetId(j->id()),                     //passing the EBDetId
                   hits,                                 //the rechit collection in order to calculate the swiss crss
                   *chStatus,                            //and the EcalChannelRecHitRcd
                   severityRecHitThreshold_,             //only consider rechits with ET >
                   spId_,                                //the SpikeId method (currently kE1OverE9 or kSwissCross)
                   spIdThreshold_                        //cut value for above
               ) >= severityLevelCut_) continue;         //then if the severity level is too high, we continue to the next rechit

            //Check based on flags to protect from recovered channels from non-read towers
            //Assumption is that v_chstatus_ is empty unless doFlagChecks() has been called
            std::vector<int>::const_iterator vit = std::find( v_chstatus_.begin(), v_chstatus_.end(),  ((EcalRecHit*)(&*j))->recoFlag() );
            if ( vit != v_chstatus_.end() ) continue; // the recHit has to be excluded from the iso sum


            if(   fabs(et) > etMin_ 
                    && fabs(energy) > energyMin_  //Changed to fabs
                    && fabs(eta-caloeta) > intStrip_ 
                    && (eta-caloeta)*(eta-caloeta) + phiDiff*phiDiff >r2 ) {

                deposit.addDeposit( Direction(eta, phi), (useEt_ ? et : energy) );

            }
        }
    }
}

virtual reco::IsoDeposit egammaisolation::EgammaRecHitExtractor::deposit	(	const edm::Event &	ev,
		const edm::EventSetup &	evSetup,
		const reco::Track &	track
	)		const `[inline, virtual]`

make single IsoDeposit based on track as input purely virtual: have to implement in concrete implementations

Implements reco::isodeposit::IsoDepositExtractor.

Definition at line 44 of file EgammaRecHitExtractor.h.

References Exception, and AlCaRecoCosmics_cfg::name.

Referenced by deposit().

                                                                                                                           {
            throw cms::Exception("Configuration Error") << "This extractor " << (typeid(this).name()) << " is not made for tracks";
         }

reco::IsoDeposit EgammaRecHitExtractor::deposit	(	const edm::Event &	ev,
		const edm::EventSetup &	evSetup,
		const reco::Candidate &	track
	)		const `[virtual]`

make single IsoDeposit based on a candidate as input purely virtual: have to implement in concrete implementations

Reimplemented from reco::isodeposit::IsoDepositExtractor.

Definition at line 87 of file EgammaRecHitExtractor.cc.

References abs, barrelEcalHitsTag_, collect(), deposit(), DetId::Ecal, EcalBarrel, EcalEndcap, endcapEcalHitsTag_, funct::exp(), fakeNegativeDeposit_, edm::EventSetup::get(), reco::Candidate::get(), edm::Event::getByLabel(), CaloGeometry::getSubdetectorGeometry(), intRadius_, metname, edm::ESHandle< T >::product(), sameTag_, funct::sin(), tryBoth_, and useEt_.

                                                                            {
    edm::ESHandle<CaloGeometry> pG;
    iSetup.get<CaloGeometryRecord>().get(pG);

    //Get the channel status from the db
    edm::ESHandle<EcalChannelStatus> chStatus;
    iSetup.get<EcalChannelStatusRcd>().get(chStatus);

    const CaloGeometry* caloGeom = pG.product(); 
    const CaloSubdetectorGeometry* barrelgeom = caloGeom->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
    const CaloSubdetectorGeometry* endcapgeom = caloGeom->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);

    static std::string metname = "EgammaIsolationAlgos|EgammaRecHitExtractor";

    std::auto_ptr<const CaloRecHitMetaCollectionV> barrelRecHits(0), endcapRecHits(0);

    //Get barrel ECAL RecHits 
    edm::Handle<EcalRecHitCollection> barrelEcalRecHitsH;
    iEvent.getByLabel(barrelEcalHitsTag_, barrelEcalRecHitsH);

    //Get endcap ECAL RecHits 
    edm::Handle<EcalRecHitCollection> endcapEcalRecHitsH;
    iEvent.getByLabel(endcapEcalHitsTag_, endcapEcalRecHitsH);


    //define isodeposit starting from candidate
    reco::SuperClusterRef sc = emObject.get<reco::SuperClusterRef>();
    math::XYZPoint caloPosition = sc->position();

    Direction candDir(caloPosition.eta(), caloPosition.phi());
    reco::IsoDeposit deposit( candDir );
    deposit.setVeto( reco::IsoDeposit::Veto(candDir, intRadius_) ); 
    double sinTheta = sin(2*atan(exp(-sc->eta())));
    deposit.addCandEnergy(sc->energy() * (useEt_ ? sinTheta : 1.0)) ;

    // subtract supercluster if desired
    double fakeEnergy = -sc->rawEnergy();
    if (fakeNegativeDeposit_) {
        deposit.addDeposit(candDir, fakeEnergy * (useEt_ ?  sinTheta : 1.0)); // not exactly clean...
    }

    // fill rechits
    bool inBarrel = sameTag_ || ( abs(sc->eta()) < 1.479 ); //check for barrel. If only one collection is used, use barrel
    if (inBarrel || tryBoth_) {
        collect(deposit, sc, barrelgeom, caloGeom, *barrelEcalRecHitsH, chStatus.product(), true);
    } 
    if ((!inBarrel) || tryBoth_) {
        collect(deposit, sc, endcapgeom, caloGeom, *endcapEcalRecHitsH, chStatus.product(), false);
    }

    return deposit;
}

virtual void egammaisolation::EgammaRecHitExtractor::fillVetos	(	const edm::Event &	ev,
		const edm::EventSetup &	evSetup,
		const reco::TrackCollection &	tracks
	)		`[inline, virtual]`