#include <HcalHaloAlgo.h>

Public Member Functions
reco::HcalHaloData	Calculate (const CaloGeometry &TheCaloGeometry, edm::Handle< HBHERecHitCollection > &TheHBHERecHits, edm::Handle< CaloTowerCollection > &TheCaloTowers, edm::Handle< EBRecHitCollection > &TheEBRecHits, edm::Handle< EERecHitCollection > &TheEERecHits, const edm::EventSetup &TheSetup)

reco::HcalHaloData	Calculate (const CaloGeometry &TheCaloGeometry, edm::Handle< HBHERecHitCollection > &TheHBHERecHits, edm::Handle< EBRecHitCollection > &TheEBRecHits, edm::Handle< EERecHitCollection > &TheEERecHits, const edm::EventSetup &TheSetup)

std::vector< reco::HaloClusterCandidateHCAL >	GetHaloClusterCandidateHB (edm::Handle< EcalRecHitCollection > &ebrechitcoll, edm::Handle< HBHERecHitCollection > &hbherechitcoll, float et_thresh_seedrh)

std::vector< reco::HaloClusterCandidateHCAL >	GetHaloClusterCandidateHE (edm::Handle< EcalRecHitCollection > &eerechitcoll, edm::Handle< HBHERecHitCollection > &hbherechitcoll, float et_thresh_seedrh)

float	GetHBRecHitEnergyThreshold ()

float	GetHERecHitEnergyThreshold ()

float	GetPhiWedgeEnergyThreshold ()

int	GetPhiWedgeNHitsThreshold ()

bool	HBClusterShapeandTimeStudy (reco::HaloClusterCandidateHCAL hcand, bool ishlt)

	HcalHaloAlgo ()

bool	HEClusterShapeandTimeStudy (reco::HaloClusterCandidateHCAL hcand, bool ishlt)

void	SetPhiWedgeEnergyThreshold (float SumE)

void	SetPhiWedgeNHitsThreshold (int nhits)

void	SetPhiWedgeThresholds (float SumE, int nhits)

void	SetRecHitEnergyThresholds (float HB, float HE)

	~HcalHaloAlgo ()

Private Member Functions
math::XYZPoint	getPosition (const DetId &id, reco::Vertex::Point vtx)

Private Attributes
const CaloGeometry *	geo_

float	HBRecHitEnergyThreshold

float	HERecHitEnergyThreshold

const HcalGeometry *	hgeo_

int	NHitsThreshold

float	SumEnergyThreshold

Detailed Description

Definition at line 41 of file HcalHaloAlgo.h.

Constructor & Destructor Documentation

HcalHaloAlgo::HcalHaloAlgo ( )

Definition at line 23 of file HcalHaloAlgo.cc.

References HBRecHitEnergyThreshold, HERecHitEnergyThreshold, NHitsThreshold, and SumEnergyThreshold.

                            : geo_(nullptr), hgeo_(nullptr) {
   HBRecHitEnergyThreshold = 0.;
   HERecHitEnergyThreshold = 0.;
   SumEnergyThreshold = 0.;
   NHitsThreshold = 0;
 }

HcalHaloAlgo::~HcalHaloAlgo ( )

inline

Definition at line 46 of file HcalHaloAlgo.h.

References Calculate().

46 {}

Member Function Documentation

HcalHaloData HcalHaloAlgo::Calculate	(	const CaloGeometry &	TheCaloGeometry,
		edm::Handle< HBHERecHitCollection > &	TheHBHERecHits,
		edm::Handle< CaloTowerCollection > &	TheCaloTowers,
		edm::Handle< EBRecHitCollection > &	TheEBRecHits,
		edm::Handle< EERecHitCollection > &	TheEERecHits,
		const edm::EventSetup &	TheSetup
	)

Definition at line 35 of file HcalHaloAlgo.cc.

References funct::abs(), CompareTime(), CompareTowers(), CaloTower::emEt(), F(), geo_, edm::EventSetup::get(), GetHaloClusterCandidateHB(), GetHaloClusterCandidateHE(), reco::HcalHaloData::GetPhiWedges(), reco::HcalHaloData::getProblematicStrips(), CaloGeometry::getSubdetectorGeometry(), CaloTower::hadEt(), HBRecHitEnergyThreshold, DetId::Hcal, HcalBarrel, HcalEndcap, HERecHitEnergyThreshold, hgeo_, mps_fire::i, hit::id, CaloTower::id(), HcalDetId::ieta(), CaloTower::ieta(), CaloTower::iphi(), SiStripPI::max, NHitsThreshold, CaloTower::numProblematicHcalCells(), edm::ESHandle< T >::product(), reco::HcalHaloData::setHaloClusterCandidatesHB(), reco::HcalHaloData::setHaloClusterCandidatesHE(), reco::PhiWedge::SetPlusZOriginConfidence(), digitizers_cfi::strip, SumEnergyThreshold, and ntuplemaker::time.

Referenced by Calculate(), reco::HcalHaloDataProducer::produce(), and ~HcalHaloAlgo().

 {
 
   HcalHaloData TheHcalHaloData;
   // ieta overlap geometrically w/ HB
   const int iEtaOverlap = 22;
   const int nPhiMax=73;
   // Store Energy sum of rechits as a function of iPhi (iPhi goes from 1 to 72)
   float SumE[nPhiMax];
   // Store Number of rechits as a function of iPhi 
   int NumHits[nPhiMax];
   // Store minimum time of rechit as a function of iPhi
   float MinTimeHits[nPhiMax];
   // Store maximum time of rechit as a function of iPhi
   float MaxTimeHits[nPhiMax];
   for(unsigned int i = 0 ; i < nPhiMax ; i++ ) 
     {
       SumE[i] = 0;
       NumHits[i]= 0;
       MinTimeHits[i] = 0.;
       MaxTimeHits[i] = 0.;
     }
   
   for (const auto & hit : (*TheHBHERecHits)) {
     HcalDetId id = HcalDetId(hit.id());
     switch ( id.subdet() ) {      
     case HcalBarrel:
       if(hit.energy() < HBRecHitEnergyThreshold )continue;
       break;
     case HcalEndcap:
       if(hit.energy() < HERecHitEnergyThreshold ) continue;
       break;
     default:
       continue;
     }
     
     int iEta = id.ieta();
     int iPhi = id.iphi();
     if(iPhi < nPhiMax && std::abs(iEta) <= iEtaOverlap) { 
       SumE[iPhi]+= hit.energy();
       NumHits[iPhi] ++;
       
       float time = hit.time();
       MinTimeHits[iPhi] = time < MinTimeHits[iPhi] ? time : MinTimeHits[iPhi];
       MaxTimeHits[iPhi] = time > MaxTimeHits[iPhi] ? time : MaxTimeHits[iPhi];
     }
   }
   
   for (int iPhi = 1 ; iPhi < nPhiMax ; iPhi++ )   {
     if( SumE[iPhi] >= SumEnergyThreshold && NumHits[iPhi] > NHitsThreshold ) {
       // Build PhiWedge and store to HcalHaloData if energy or #hits pass thresholds
       PhiWedge wedge(SumE[iPhi], iPhi, NumHits[iPhi], MinTimeHits[iPhi], MaxTimeHits[iPhi]);
       
       // Loop over rechits again to calculate direction based on timing info
       std::vector<const HBHERecHit*> Hits;
       for (const auto & hit : (*TheHBHERecHits)) {
     HcalDetId id = HcalDetId(hit.id());
     if( id.iphi() != iPhi ) continue;
     if( std::abs(id.ieta() ) > iEtaOverlap ) continue;  // has to overlap geometrically w/ HB
     switch ( id.subdet() ) {
     case HcalBarrel:
       if(hit.energy() < HBRecHitEnergyThreshold )continue;
       break;
     case HcalEndcap:
       if(hit.energy() < HERecHitEnergyThreshold ) continue;
       break;
     default:
       continue;
     }
     Hits.push_back(&(hit));
       }
           
       std::sort( Hits.begin() , Hits.end() , CompareTime);
       float MinusToPlus = 0.;
       float PlusToMinus = 0.;
       for( unsigned int i = 0 ; i < Hits.size() ; i++ )  {
     HcalDetId id_i = HcalDetId(Hits[i]->id() );
     int ieta_i = id_i.ieta();
     for( unsigned int j = (i+1) ; j < Hits.size() ; j++ ) {
       HcalDetId id_j = HcalDetId(Hits[j]->id() );
       int ieta_j = id_j.ieta();
       if( ieta_i > ieta_j ) PlusToMinus += std::abs(ieta_i - ieta_j ) ;
       else                  MinusToPlus += std::abs(ieta_i - ieta_j);
     }
       }
       float PlusZOriginConfidence = (PlusToMinus + MinusToPlus )? PlusToMinus / ( PlusToMinus + MinusToPlus ) : -1. ;
       wedge.SetPlusZOriginConfidence( PlusZOriginConfidence );
       TheHcalHaloData.GetPhiWedges().push_back( wedge );
     }
   }
 
 
   // Don't use HF.
   int maxAbsIEta = 29;
 
 
   std::map<int, float> iPhiHadEtMap;
   std::vector<const CaloTower*> sortedCaloTowers;
   for (const auto & tower : (*TheCaloTowers)) {
     if(std::abs(tower.ieta()) > maxAbsIEta) continue;
     
     int iPhi = tower.iphi();
     if(!iPhiHadEtMap.count(iPhi)) iPhiHadEtMap[iPhi] = 0.0;
     iPhiHadEtMap[iPhi] += tower.hadEt();
 
     if(tower.numProblematicHcalCells() > 0) sortedCaloTowers.push_back(&(tower));
     
   }
 
 
   // Sort towers such that lowest iphi and ieta are first, highest last, and towers
   // with same iphi value are consecutive. Then we can do everything else in one loop.
   std::sort(sortedCaloTowers.begin(), sortedCaloTowers.end(), CompareTowers);
 
   HaloTowerStrip strip;
 
 
   int prevIEta = -99, prevIPhi = -99;
   float prevHadEt = 0.;
   float prevEmEt = 0.;
   std::pair<uint8_t, CaloTowerDetId> prevPair, towerPair;
   bool wasContiguous = true;
 
   // Loop through and store a vector of pairs (problematicCells, DetId) for each contiguous strip we find
   for(unsigned int i = 0; i < sortedCaloTowers.size(); i++) {
     const CaloTower* tower = sortedCaloTowers[i];
 
     towerPair = std::make_pair((uint8_t)tower->numProblematicHcalCells(), tower->id());
 
     bool newIPhi = tower->iphi() != prevIPhi;
     bool isContiguous = tower->ieta() == 1 ? tower->ieta() - 2 == prevIEta : tower->ieta() - 1 == prevIEta;
 
     isContiguous = isContiguous || (tower->ieta() == -maxAbsIEta);
     if(newIPhi) isContiguous = false;
 
     if(!wasContiguous && isContiguous) {
       strip.cellTowerIds.push_back(prevPair);
       strip.cellTowerIds.push_back(towerPair);
       strip.hadEt += prevHadEt + tower->hadEt();
       strip.emEt += prevEmEt + tower->emEt();
     }
 
     if(wasContiguous && isContiguous) {
       strip.cellTowerIds.push_back(towerPair);
       strip.hadEt += tower->hadEt();
       strip.emEt += tower->emEt();
     }
 
     if((wasContiguous && !isContiguous) || i == sortedCaloTowers.size()-1) { //ended the strip, so flush it
 
       if(strip.cellTowerIds.size() > 3) {
 
         int iPhi = strip.cellTowerIds.at(0).second.iphi();
         int iPhiLower = (iPhi == 1) ? 72 : iPhi - 1;
         int iPhiUpper = (iPhi == 72) ? 1 : iPhi + 1;
 
         float energyRatio = 0.0;
         if(iPhiHadEtMap.count(iPhiLower)) energyRatio += iPhiHadEtMap[iPhiLower];
         if(iPhiHadEtMap.count(iPhiUpper)) energyRatio += iPhiHadEtMap[iPhiUpper];
         iPhiHadEtMap[iPhi] = std::max(iPhiHadEtMap[iPhi], 0.001F);
 
         energyRatio /= iPhiHadEtMap[iPhi];
         strip.energyRatio = energyRatio;
 
         TheHcalHaloData.getProblematicStrips().push_back( strip );
 
       }
       strip = HaloTowerStrip();
     }
 
     wasContiguous = isContiguous;
     prevPair = towerPair;
     prevEmEt = tower->emEt();
     prevIPhi = tower->iphi();
     prevIEta = tower->ieta();
     prevHadEt = tower->hadEt();
   }
 
   edm::ESHandle<CaloGeometry> pGeo;
   TheSetup.get<CaloGeometryRecord>().get(pGeo);
   geo_  = pGeo.product();
   hgeo_ = dynamic_cast<const HcalGeometry*>(geo_->getSubdetectorGeometry(DetId::Hcal, 1));
 
   //Halo cluster building:
   //Various clusters are built, depending on the subdetector. 
   //In barrel, one looks for deposits narrow in phi.
   //In endcaps, one looks for localized deposits (dr condition in EE where r =sqrt(dphi*dphi+deta*deta)
   //E/H condition is also applied.
   //The halo cluster building step targets a large efficiency (ideally >99%) for beam halo deposits.
   //These clusters are used as input for the halo pattern finding methods in HcalHaloAlgo and for the CSC-calo matching methods in GlobalHaloAlgo.
   
   //Et threshold hardcoded for now. Might one to get it from config
 
   std::vector<HaloClusterCandidateHCAL> haloclustercands_HB;
   haloclustercands_HB=  GetHaloClusterCandidateHB(TheEBRecHits , TheHBHERecHits, 5);
 
   std::vector<HaloClusterCandidateHCAL> haloclustercands_HE;
   haloclustercands_HE=  GetHaloClusterCandidateHE(TheEERecHits , TheHBHERecHits, 10);
 
   TheHcalHaloData.setHaloClusterCandidatesHB(haloclustercands_HB);
   TheHcalHaloData.setHaloClusterCandidatesHE(haloclustercands_HE);
 
 
   return TheHcalHaloData;
   
 }

HcalHaloData HcalHaloAlgo::Calculate	(	const CaloGeometry &	TheCaloGeometry,
		edm::Handle< HBHERecHitCollection > &	TheHBHERecHits,
		edm::Handle< EBRecHitCollection > &	TheEBRecHits,
		edm::Handle< EERecHitCollection > &	TheEERecHits,
		const edm::EventSetup &	TheSetup
	)

Definition at line 30 of file HcalHaloAlgo.cc.

References Calculate().

                                                                                                                                                                                                                                                {
     edm::Handle<CaloTowerCollection> TheCaloTowers;
     return Calculate(TheCaloGeometry, TheHBHERecHits, TheCaloTowers,TheEBRecHits,TheEERecHits,TheSetup);
 }

std::vector< HaloClusterCandidateHCAL > HcalHaloAlgo::GetHaloClusterCandidateHB	(	edm::Handle< EcalRecHitCollection > &	ebrechitcoll,
		edm::Handle< HBHERecHitCollection > &	hbherechitcoll,
		float	et_thresh_seedrh
	)

Definition at line 249 of file HcalHaloAlgo.cc.

Referenced by Calculate(), and GetPhiWedgeNHitsThreshold().

                                                                                                                                                                                                   {
 
   std::vector<HaloClusterCandidateHCAL> TheHaloClusterCandsHB;
 
   reco::Vertex::Point vtx(0,0,0);
 
   for(size_t ihit = 0; ihit<hbherechitcoll->size(); ++ ihit){
     HaloClusterCandidateHCAL  clustercand;
     
     const HBHERecHit & rechit = (*hbherechitcoll)[ ihit ];
     math::XYZPoint rhpos = getPosition(rechit.id(),vtx);
     //Et condition
     double rhet = rechit.energy()* sqrt(rhpos.perp2()/rhpos.mag2());
     if(rhet<et_thresh_seedrh) continue;
     if(std::abs(rhpos.z())>zseparation_HBHE) continue;
     double eta = rhpos.eta();
     double phi = rhpos.phi();
     
     bool isiso = true;
     double etcluster(0);
     int nbtowerssameeta(0);
     double timediscriminatorITBH(0),timediscriminatorOTBH(0);
     double etstrip_phiseedplus1(0), etstrip_phiseedminus1(0);
 
     //Building the cluster 
     edm::RefVector<HBHERecHitCollection> bhrhcandidates;
     for(size_t jhit = 0; jhit<hbherechitcoll->size(); ++ jhit){
       const HBHERecHit & rechitj = (*hbherechitcoll)[ jhit ];
       HBHERecHitRef rhRef(hbherechitcoll,jhit);
       math::XYZPoint rhposj = getPosition(rechitj.id(),vtx);
       double rhetj = rechitj.energy()* sqrt(rhposj.perp2()/rhposj.mag2());
       if(rhetj<2) continue;
       if(std::abs(rhposj.z())>zseparation_HBHE) continue;
       double etaj = rhposj.eta();
       double phij = rhposj.phi();
       double deta = eta - etaj;
       double dphi = deltaPhi(phi,phij);
       if(std::abs(deta)>0.4) continue;//This means +/-4 towers in eta 
       if(std::abs(dphi)>0.2) continue;//This means +/-2 towers in phi 
       if(std::abs(dphi)>0.1&&std::abs(deta)<0.2){isiso=false;break;}//The strip should be isolated
       if(std::abs(dphi)>0.1)continue;
       if(std::abs(dphi)<0.05) nbtowerssameeta++;
       if(dphi>0.05) etstrip_phiseedplus1+=rhetj;
       if(dphi<-0.05) etstrip_phiseedminus1+=rhetj;
       
       etcluster+=rhetj;
       //Timing discriminator
       //We assign a weight to the rechit defined as: 
       //Log10(Et)*f(T,R,Z)
       //where f(T,R,Z) is the separation curve between halo-like and IP-like times.
       //The time difference between a deposit from a outgoing IT halo and a deposit coming from a particle emitted at the IP is given by:  
       //dt= ( - sqrt(R^2+z^2) + |z| )/c 
       // For OT beam halo, the time difference is: 
       //dt= ( 25 + sqrt(R^2+z^2) + |z| )/c 
       //only consider the central part of HB as things get hard at large z.
       //The best fitted value for R leads to 240 cm (IT) and 330 cm (OT)
       double rhtj = rechitj.time();
       timediscriminatorITBH+= std::log10( rhetj )* ( rhtj +0.5*(sqrt(240.*240.+rhposj.z()*rhposj.z()) -std::abs(rhposj.z()))/c_cm_per_ns);
       if(std::abs(rhposj.z())<300) timediscriminatorOTBH+= std::log10( rhetj )* ( rhtj -0.5*(25-(sqrt(330.*330.+rhposj.z()*rhposj.z()) +std::abs(rhposj.z()))/c_cm_per_ns) );
       bhrhcandidates.push_back(rhRef);
     }
     //Isolation conditions
     if(!isiso) continue;
     if(etstrip_phiseedplus1/etcluster>0.2&& etstrip_phiseedminus1/etcluster>0.2) continue;
     
     //Calculate E/H
     double eoh(0);
     for(size_t jhit = 0; jhit<ecalrechitcoll->size(); ++ jhit){
       const EcalRecHit & rechitj = (*ecalrechitcoll)[ jhit ];
       math::XYZPoint rhposj = getPosition(rechitj.id(),vtx);
       double rhetj = rechitj.energy()* sqrt(rhposj.perp2()/rhposj.mag2());
       if(rhetj<2) continue;
       double etaj = rhposj.eta();
       double phij = rhposj.phi();
       if(std::abs(eta-etaj)>0.2) continue;
       if(std::abs(deltaPhi(phi,phij))>0.2) continue;
       eoh+=rhetj/etcluster;
     }
     //E/H condition
     if(eoh>0.1) continue; 
         
 
     clustercand.setClusterEt(etcluster); 
     clustercand.setSeedEt(rhet); 
     clustercand.setSeedEta(eta); 
     clustercand.setSeedPhi(phi); 
     clustercand.setSeedZ(rhpos.Z());
     clustercand.setSeedR(sqrt(rhpos.perp2())); 
     clustercand.setSeedTime(rechit.time()); 
     clustercand.setEoverH(eoh);
     clustercand.setNbTowersInEta(nbtowerssameeta);
     clustercand.setEtStripPhiSeedPlus1(etstrip_phiseedplus1);
     clustercand.setEtStripPhiSeedMinus1(etstrip_phiseedminus1);
     clustercand.setTimeDiscriminatorITBH(timediscriminatorITBH);
     clustercand.setTimeDiscriminatorOTBH(timediscriminatorOTBH);
     clustercand.setBeamHaloRecHitsCandidates(bhrhcandidates);
 
     bool isbeamhalofrompattern = HBClusterShapeandTimeStudy(clustercand,false);
     clustercand.setIsHaloFromPattern(isbeamhalofrompattern);
     bool isbeamhalofrompattern_hlt = HBClusterShapeandTimeStudy(clustercand,true);
     clustercand.setIsHaloFromPattern_HLT(isbeamhalofrompattern_hlt);
 
 
     TheHaloClusterCandsHB.push_back(clustercand);
   }
 
   return  TheHaloClusterCandsHB;
 } 

std::vector< HaloClusterCandidateHCAL > HcalHaloAlgo::GetHaloClusterCandidateHE	(	edm::Handle< EcalRecHitCollection > &	eerechitcoll,
		edm::Handle< HBHERecHitCollection > &	hbherechitcoll,
		float	et_thresh_seedrh
	)

Definition at line 359 of file HcalHaloAlgo.cc.

Referenced by Calculate(), and GetPhiWedgeNHitsThreshold().

                                                                                                                                                                                                   {
 
   std::vector<HaloClusterCandidateHCAL> TheHaloClusterCandsHE;
 
   reco::Vertex::Point vtx(0,0,0);
 
   for(size_t ihit = 0; ihit<hbherechitcoll->size(); ++ ihit){
     HaloClusterCandidateHCAL  clustercand;
     
     const HBHERecHit & rechit = (*hbherechitcoll)[ ihit ];
     math::XYZPoint rhpos = getPosition(rechit.id(),vtx);
     //Et condition
     double rhet = rechit.energy()* sqrt(rhpos.perp2()/rhpos.mag2());
     if(rhet<et_thresh_seedrh) continue;
     if(std::abs(rhpos.z())<zseparation_HBHE) continue;
     double eta = rhpos.eta();
     double phi = rhpos.phi();
     double rhr = sqrt(rhpos.perp2());
     bool isiso = true;
     double etcluster(0),hdepth1(0);
     int clustersize(0);
     double etstrip_phiseedplus1(0), etstrip_phiseedminus1(0);
 
     //Building the cluster 
     edm::RefVector<HBHERecHitCollection> bhrhcandidates;
     for(size_t jhit = 0; jhit<hbherechitcoll->size(); ++ jhit){
       const HBHERecHit & rechitj = (*hbherechitcoll)[ jhit ];
       HBHERecHitRef rhRef(hbherechitcoll,jhit);
       math::XYZPoint rhposj = getPosition(rechitj.id(),vtx);
       double rhetj = rechitj.energy()* sqrt(rhposj.perp2()/rhposj.mag2());
       if(rhetj<2) continue;
       if(std::abs(rhposj.z())<zseparation_HBHE) continue;
       if(rhpos.z()*rhposj.z()<0) continue;
       double phij = rhposj.phi();
       double dphi = deltaPhi(phi,phij);
       if(std::abs(dphi)>0.4) continue;
       double rhrj = sqrt(rhposj.perp2()); 
       if(std::abs( rhr-rhrj )>50) continue;
       if(std::abs(dphi)>0.2 ||std::abs( rhr-rhrj )>20 ){isiso=false;break;}//The deposit should be isolated
       if(dphi>0.05) etstrip_phiseedplus1+=rhetj;
       if(dphi<-0.05) etstrip_phiseedminus1+=rhetj;
       clustersize++;
       etcluster+=rhetj;
       if(std::abs( rhposj.z())<405 )hdepth1+=rhetj;
       //No timing condition for now in HE
       bhrhcandidates.push_back(rhRef);
     }
     //Isolation conditions
     if(!isiso) continue;
     if(etstrip_phiseedplus1/etcluster>0.1&& etstrip_phiseedminus1/etcluster>0.1) continue;
     
     //Calculate E/H
     double eoh(0);
     for(size_t jhit = 0; jhit<ecalrechitcoll->size(); ++ jhit){
       const EcalRecHit & rechitj = (*ecalrechitcoll)[ jhit ];
       math::XYZPoint rhposj = getPosition(rechitj.id(),vtx);
       double rhetj = rechitj.energy()* sqrt(rhposj.perp2()/rhposj.mag2());
       if(rhetj<2) continue;
       if(rhpos.z()*rhposj.z()<0) continue;
       double etaj = rhposj.eta();
       double phij = rhposj.phi();
       double dr = sqrt( (eta-etaj)*(eta-etaj)+deltaPhi(phi,phij)*deltaPhi(phi,phij));
       if(dr>0.3) continue;
 
       eoh+=rhetj/etcluster;
     }
     //E/H condition
     if(eoh>0.1) continue; 
         
 
     clustercand.setClusterEt(etcluster); 
     clustercand.setSeedEt(rhet); 
     clustercand.setSeedEta(eta); 
     clustercand.setSeedPhi(phi); 
     clustercand.setSeedZ(rhpos.Z());
     clustercand.setSeedR(sqrt(rhpos.perp2())); 
     clustercand.setSeedTime(rechit.time()); 
     clustercand.setEoverH(eoh);
     clustercand.setH1overH123(hdepth1/etcluster);
     clustercand.setClusterSize(clustersize);
     clustercand.setEtStripPhiSeedPlus1(etstrip_phiseedplus1);
     clustercand.setEtStripPhiSeedMinus1(etstrip_phiseedminus1);
     clustercand.setTimeDiscriminator(0);
     clustercand.setBeamHaloRecHitsCandidates(bhrhcandidates);
 
     bool isbeamhalofrompattern = HEClusterShapeandTimeStudy(clustercand,false);
     clustercand.setIsHaloFromPattern(isbeamhalofrompattern);
     bool isbeamhalofrompattern_hlt = HEClusterShapeandTimeStudy(clustercand,true);
     clustercand.setIsHaloFromPattern_HLT(isbeamhalofrompattern_hlt);
 
 
     TheHaloClusterCandsHE.push_back(clustercand);
   }
 
   return  TheHaloClusterCandsHE;
 } 

float HcalHaloAlgo::GetHBRecHitEnergyThreshold ( )

inline

Definition at line 62 of file HcalHaloAlgo.h.

References HBRecHitEnergyThreshold.

62 { return HBRecHitEnergyThreshold;}

HcalHaloAlgo::HBRecHitEnergyThreshold

float HBRecHitEnergyThreshold

Definition: HcalHaloAlgo.h:79

float HcalHaloAlgo::GetHERecHitEnergyThreshold ( )

inline

Definition at line 63 of file HcalHaloAlgo.h.

References HERecHitEnergyThreshold.

63 { return HERecHitEnergyThreshold;}

HcalHaloAlgo::HERecHitEnergyThreshold

float HERecHitEnergyThreshold

Definition: HcalHaloAlgo.h:80

float HcalHaloAlgo::GetPhiWedgeEnergyThreshold ( )

inline

Definition at line 66 of file HcalHaloAlgo.h.

References SumEnergyThreshold.

66 { return SumEnergyThreshold;}

HcalHaloAlgo::SumEnergyThreshold

float SumEnergyThreshold

Definition: HcalHaloAlgo.h:83

int HcalHaloAlgo::GetPhiWedgeNHitsThreshold ( )

inline

Definition at line 67 of file HcalHaloAlgo.h.

References GetHaloClusterCandidateHB(), GetHaloClusterCandidateHE(), HBClusterShapeandTimeStudy(), HEClusterShapeandTimeStudy(), and NHitsThreshold.

67 { return NHitsThreshold;}

HcalHaloAlgo::NHitsThreshold

int NHitsThreshold

Definition: HcalHaloAlgo.h:84

math::XYZPoint HcalHaloAlgo::getPosition	(	const DetId &	id,
		reco::Vertex::Point	vtx
	)

private

Definition at line 507 of file HcalHaloAlgo.cc.

References geo_, CaloGeometry::getPosition(), HcalGeometry::getPosition(), DetId::Hcal, hgeo_, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by GetHaloClusterCandidateHB(), and GetHaloClusterCandidateHE().

                                                                           {
 
   const GlobalPoint pos = ((id.det() == DetId::Hcal) ? hgeo_->getPosition(id) :
                GlobalPoint(geo_->getPosition(id)));
   math::XYZPoint posV(pos.x() - vtx.x(),pos.y() - vtx.y(),pos.z() - vtx.z());
   return posV;
 }

bool HcalHaloAlgo::HBClusterShapeandTimeStudy	(	reco::HaloClusterCandidateHCAL	hcand,
		bool	ishlt
	)

Definition at line 458 of file HcalHaloAlgo.cc.

References reco::HaloClusterCandidateHCAL::getEtStripPhiSeedMinus1(), reco::HaloClusterCandidateHCAL::getEtStripPhiSeedPlus1(), reco::HaloClusterCandidateHCAL::getNbTowersInEta(), reco::HaloClusterCandidateHCAL::getSeedEt(), reco::HaloClusterCandidateHCAL::getTimeDiscriminatorITBH(), reco::HaloClusterCandidateHCAL::getTimeDiscriminatorOTBH(), and reco::HaloClusterCandidateHCAL::setIsHaloFromPattern().

Referenced by GetHaloClusterCandidateHB(), and GetPhiWedgeNHitsThreshold().

                                                                                         {
   //Conditions on the central strip size in eta.
   //For low size, extra conditions on seed et, isolation and cluster timing 
   //Here we target both IT and OT beam halo. Two separate discriminators were built for the two cases.
   
   if(hcand.getSeedEt()<10)return false;
   
   if(hcand.getNbTowersInEta()<3) return false;
   //Isolation criteria for very short eta strips
   if(hcand.getNbTowersInEta()==3 && (hcand.getEtStripPhiSeedPlus1()>0.1 || hcand.getEtStripPhiSeedMinus1()>0.1) ) return false;
   if(hcand.getNbTowersInEta()<=5 && (hcand.getEtStripPhiSeedPlus1()>0.1 && hcand.getEtStripPhiSeedMinus1()>0.1) ) return false;
   
   //Timing conditions for short eta strips
   if(hcand.getNbTowersInEta()==3 && hcand.getTimeDiscriminatorITBH()>=0.) return false;
   if(hcand.getNbTowersInEta()<=6 && hcand.getTimeDiscriminatorITBH()>=5. &&hcand.getTimeDiscriminatorOTBH()<0.) return false; 
   
   //For HLT, only use conditions without timing 
   if(ishlt && hcand.getNbTowersInEta()<7) return false;
 
   hcand.setIsHaloFromPattern(true);
   
   return true;
 }

bool HcalHaloAlgo::HEClusterShapeandTimeStudy	(	reco::HaloClusterCandidateHCAL	hcand,
		bool	ishlt
	)

Definition at line 484 of file HcalHaloAlgo.cc.

References reco::HaloClusterCandidateHCAL::getH1overH123(), reco::HaloClusterCandidateHCAL::getSeedEt(), reco::HaloClusterCandidateHCAL::getSeedR(), and reco::HaloClusterCandidateHCAL::setIsHaloFromPattern().

Referenced by GetHaloClusterCandidateHE(), and GetPhiWedgeNHitsThreshold().

                                                                                         {
   //Conditions on H1/H123 to spot halo interacting only in one HCAL layer. 
   //For R> about 170cm, HE has only one layer and this condition cannot be applied
   //Note that for R>170 cm, the halo is in CSC acceptance and will most likely be spotted by the CSC-calo matching method
   //A method to identify halos interacting in both H1 and H2/H3 at low R is still missing. 
   
   if(hcand.getSeedEt()<20)return false;
   if(hcand.getSeedR()>170) return false;
   
   if(hcand.getH1overH123()>0.02 &&hcand.getH1overH123()<0.98) return false;
   
   //This method is one of the ones with the highest fake rate: in JetHT dataset, it happens in around 0.1% of the cases that a low pt jet (pt= 20) leaves all of its energy in only one HCAL layer. 
   //At HLT, one only cares about large deposits from BH that would lead to a MET/SinglePhoton trigger to be fired.
   //Rising the seed Et threshold at HLT has therefore little impact on the HLT performances but ensures that possible controversial events are still recorded.
   if(ishlt && hcand.getSeedEt()<50)return false;
   
   hcand.setIsHaloFromPattern(true);
   
   return true;
   
 }