CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions | Private Member Functions | Private Attributes
HcalHaloAlgo Class Reference

#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 CaloGeometrygeo
 
float HBRecHitEnergyThreshold
 
float HERecHitEnergyThreshold
 
int NHitsThreshold
 
float SumEnergyThreshold
 

Detailed Description

Definition at line 47 of file HcalHaloAlgo.h.

Constructor & Destructor Documentation

HcalHaloAlgo::HcalHaloAlgo ( )

Definition at line 26 of file HcalHaloAlgo.cc.

27 {
30  SumEnergyThreshold = 0.;
31  NHitsThreshold = 0;
32 
33  geo = 0;
34 }
float HERecHitEnergyThreshold
Definition: HcalHaloAlgo.h:86
float HBRecHitEnergyThreshold
Definition: HcalHaloAlgo.h:85
float SumEnergyThreshold
Definition: HcalHaloAlgo.h:89
int NHitsThreshold
Definition: HcalHaloAlgo.h:90
const CaloGeometry * geo
Definition: HcalHaloAlgo.h:92
HcalHaloAlgo::~HcalHaloAlgo ( )
inline

Definition at line 52 of file HcalHaloAlgo.h.

52 {}

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 41 of file HcalHaloAlgo.cc.

References funct::abs(), Abs(), CompareTime(), CompareTowers(), CaloTower::emEt(), F(), edm::EventSetup::get(), reco::HcalHaloData::GetPhiWedges(), reco::HcalHaloData::getProblematicStrips(), CaloTower::hadEt(), HcalBarrel, HcalEndcap, i, hit::id, CaloTower::id(), HcalDetId::ieta(), CaloTower::ieta(), CaloTower::iphi(), j, bookConverter::max, CaloTower::numProblematicHcalCells(), edm::ESHandle< class >::product(), reco::HcalHaloData::setHaloClusterCandidatesHB(), reco::HcalHaloData::setHaloClusterCandidatesHE(), and reco::PhiWedge::SetPlusZOriginConfidence().

Referenced by reco::HcalHaloDataProducer::produce().

42 {
43 
44  HcalHaloData TheHcalHaloData;
45 
46  // Store Energy sum of rechits as a function of iPhi (iPhi goes from 1 to 72)
47  float SumE[73];
48  // Store Number of rechits as a function of iPhi
49  int NumHits[73];
50  // Store minimum time of rechit as a function of iPhi
51  float MinTimeHits[73];
52  // Store maximum time of rechit as a function of iPhi
53  float MaxTimeHits[73];
54  for(unsigned int i = 0 ; i < 73 ; i++ )
55  {
56  SumE[i] = 0;
57  NumHits[i]= 0;
58  MinTimeHits[i] = 0.;
59  MaxTimeHits[i] = 0.;
60  }
61 
62  for( HBHERecHitCollection::const_iterator hit = TheHBHERecHits->begin() ; hit != TheHBHERecHits->end() ; hit++ )
63  {
64  HcalDetId id = HcalDetId(hit->id());
65  switch ( id.subdet() )
66  {
67  case HcalBarrel:
68  if(hit->energy() < HBRecHitEnergyThreshold )continue;
69  break;
70  case HcalEndcap:
71  if(hit->energy() < HERecHitEnergyThreshold ) continue;
72  break;
73  default:
74  continue;
75  }
76 
77  int iEta = id.ieta();
78  int iPhi = id.iphi();
79  if(iPhi < 73 && TMath::Abs(iEta) < 23 )
80  {
81  SumE[iPhi]+= hit->energy();
82  NumHits[iPhi] ++;
83 
84  float time = hit->time();
85  MinTimeHits[iPhi] = time < MinTimeHits[iPhi] ? time : MinTimeHits[iPhi];
86  MaxTimeHits[iPhi] = time > MaxTimeHits[iPhi] ? time : MaxTimeHits[iPhi];
87  }
88  }
89 
90  for( int iPhi = 1 ; iPhi < 73 ; iPhi++ )
91  {
92  if( SumE[iPhi] >= SumEnergyThreshold && NumHits[iPhi] > NHitsThreshold )
93  {
94  // Build PhiWedge and store to HcalHaloData if energy or #hits pass thresholds
95  PhiWedge wedge(SumE[iPhi], iPhi, NumHits[iPhi], MinTimeHits[iPhi], MaxTimeHits[iPhi]);
96 
97  // Loop over rechits again to calculate direction based on timing info
98  std::vector<const HBHERecHit*> Hits;
99  for( HBHERecHitCollection::const_iterator hit = TheHBHERecHits->begin() ; hit != TheHBHERecHits->end() ; hit++ )
100  {
101 
102  HcalDetId id = HcalDetId(hit->id());
103  if( id.iphi() != iPhi ) continue;
104  if( TMath::Abs(id.ieta() ) > 22 ) continue; // has to overlap geometrically w/ HB
105  switch ( id.subdet() )
106  {
107  case HcalBarrel:
108  if(hit->energy() < HBRecHitEnergyThreshold )continue;
109  break;
110  case HcalEndcap:
111  if(hit->energy() < HERecHitEnergyThreshold ) continue;
112  break;
113  default:
114  continue;
115  }
116  Hits.push_back(&(*hit));
117  }
118 
119  std::sort( Hits.begin() , Hits.end() , CompareTime);
120  float MinusToPlus = 0.;
121  float PlusToMinus = 0.;
122  for( unsigned int i = 0 ; i < Hits.size() ; i++ )
123  {
124  HcalDetId id_i = HcalDetId(Hits[i]->id() );
125  int ieta_i = id_i.ieta();
126  for( unsigned int j = (i+1) ; j < Hits.size() ; j++ )
127  {
128  HcalDetId id_j = HcalDetId(Hits[j]->id() );
129  int ieta_j = id_j.ieta();
130  if( ieta_i > ieta_j ) PlusToMinus += TMath::Abs(ieta_i - ieta_j ) ;
131  else MinusToPlus += TMath::Abs(ieta_i - ieta_j);
132  }
133  }
134  float PlusZOriginConfidence = (PlusToMinus + MinusToPlus )? PlusToMinus / ( PlusToMinus + MinusToPlus ) : -1. ;
135  wedge.SetPlusZOriginConfidence( PlusZOriginConfidence );
136  TheHcalHaloData.GetPhiWedges().push_back( wedge );
137  }
138  }
139 
140 
141  // Don't use HF.
142  int maxAbsIEta = 29;
143 
144 
145  std::map<int, float> iPhiHadEtMap;
146  std::vector<const CaloTower*> sortedCaloTowers;
147  for(CaloTowerCollection::const_iterator tower = TheCaloTowers->begin(); tower != TheCaloTowers->end(); tower++) {
148  if(abs(tower->ieta()) > maxAbsIEta) continue;
149 
150  int iPhi = tower->iphi();
151  if(!iPhiHadEtMap.count(iPhi)) iPhiHadEtMap[iPhi] = 0.0;
152  iPhiHadEtMap[iPhi] += tower->hadEt();
153 
154  if(tower->numProblematicHcalCells() > 0) sortedCaloTowers.push_back(&(*tower));
155 
156  }
157 
158 
159  // Sort towers such that lowest iphi and ieta are first, highest last, and towers
160  // with same iphi value are consecutive. Then we can do everything else in one loop.
161  std::sort(sortedCaloTowers.begin(), sortedCaloTowers.end(), CompareTowers);
162 
163  HaloTowerStrip strip;
164 
165 
166  int prevIEta = -99, prevIPhi = -99;
167  float prevHadEt = 0.;
168  float prevEmEt = 0.;
169  std::pair<uint8_t, CaloTowerDetId> prevPair, towerPair;
170  bool wasContiguous = true;
171 
172  // Loop through and store a vector of pairs (problematicCells, DetId) for each contiguous strip we find
173  for(unsigned int i = 0; i < sortedCaloTowers.size(); i++) {
174  const CaloTower* tower = sortedCaloTowers[i];
175 
176  towerPair = std::make_pair((uint8_t)tower->numProblematicHcalCells(), tower->id());
177 
178  bool newIPhi = tower->iphi() != prevIPhi;
179  bool isContiguous = tower->ieta() == 1 ? tower->ieta() - 2 == prevIEta : tower->ieta() - 1 == prevIEta;
180 
181  isContiguous = isContiguous || (tower->ieta() == -maxAbsIEta);
182  if(newIPhi) isContiguous = false;
183 
184  if(!wasContiguous && isContiguous) {
185  strip.cellTowerIds.push_back(prevPair);
186  strip.cellTowerIds.push_back(towerPair);
187  strip.hadEt += prevHadEt + tower->hadEt();
188  strip.emEt += prevEmEt + tower->emEt();
189  }
190 
191  if(wasContiguous && isContiguous) {
192  strip.cellTowerIds.push_back(towerPair);
193  strip.hadEt += tower->hadEt();
194  strip.emEt += tower->emEt();
195  }
196 
197  if((wasContiguous && !isContiguous) || i == sortedCaloTowers.size()-1) { //ended the strip, so flush it
198 
199  if(strip.cellTowerIds.size() > 3) {
200 
201  int iPhi = strip.cellTowerIds.at(0).second.iphi();
202  int iPhiLower = (iPhi == 1) ? 72 : iPhi - 1;
203  int iPhiUpper = (iPhi == 72) ? 1 : iPhi + 1;
204 
205  float energyRatio = 0.0;
206  if(iPhiHadEtMap.count(iPhiLower)) energyRatio += iPhiHadEtMap[iPhiLower];
207  if(iPhiHadEtMap.count(iPhiUpper)) energyRatio += iPhiHadEtMap[iPhiUpper];
208  iPhiHadEtMap[iPhi] = max(iPhiHadEtMap[iPhi], 0.001F);
209 
210  energyRatio /= iPhiHadEtMap[iPhi];
211  strip.energyRatio = energyRatio;
212 
213  TheHcalHaloData.getProblematicStrips().push_back( strip );
214 
215  }
216  strip = HaloTowerStrip();
217  }
218 
219  wasContiguous = isContiguous;
220  prevPair = towerPair;
221  prevEmEt = tower->emEt();
222  prevIPhi = tower->iphi();
223  prevIEta = tower->ieta();
224  prevHadEt = tower->hadEt();
225  }
226 
227  geo = 0;
229  TheSetup.get<CaloGeometryRecord>().get(pGeo);
230  geo = pGeo.product();
231 
232  //Halo cluster building:
233  //Various clusters are built, depending on the subdetector.
234  //In barrel, one looks for deposits narrow in phi.
235  //In endcaps, one looks for localized deposits (dr condition in EE where r =sqrt(dphi*dphi+deta*deta)
236  //E/H condition is also applied.
237  //The halo cluster building step targets a large efficiency (ideally >99%) for beam halo deposits.
238  //These clusters are used as input for the halo pattern finding methods in HcalHaloAlgo and for the CSC-calo matching methods in GlobalHaloAlgo.
239 
240  //Et threshold hardcoded for now. Might one to get it from config
241 
242  std::vector<HaloClusterCandidateHCAL> haloclustercands_HB;
243  haloclustercands_HB= GetHaloClusterCandidateHB(TheEBRecHits , TheHBHERecHits, 5);
244 
245  std::vector<HaloClusterCandidateHCAL> haloclustercands_HE;
246  haloclustercands_HE= GetHaloClusterCandidateHE(TheEERecHits , TheHBHERecHits, 10);
247 
248  TheHcalHaloData.setHaloClusterCandidatesHB(haloclustercands_HB);
249  TheHcalHaloData.setHaloClusterCandidatesHE(haloclustercands_HE);
250 
251 
252  return TheHcalHaloData;
253 
254 }
int i
Definition: DBlmapReader.cc:9
float HERecHitEnergyThreshold
Definition: HcalHaloAlgo.h:86
void setHaloClusterCandidatesHB(const std::vector< HaloClusterCandidateHCAL > &x)
Definition: HcalHaloData.h:61
int ieta() const
Definition: CaloTower.h:182
double hadEt() const
Definition: CaloTower.h:113
std::vector< HBHERecHit >::const_iterator const_iterator
int iphi() const
Definition: CaloTower.h:184
std::vector< reco::HaloClusterCandidateHCAL > GetHaloClusterCandidateHE(edm::Handle< EcalRecHitCollection > &eerechitcoll, edm::Handle< HBHERecHitCollection > &hbherechitcoll, float et_thresh_seedrh)
T Abs(T a)
Definition: MathUtil.h:49
int ieta() const
get the cell ieta
Definition: HcalDetId.h:56
void setHaloClusterCandidatesHE(const std::vector< HaloClusterCandidateHCAL > &x)
Definition: HcalHaloData.h:62
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
unsigned int id
CaloTowerDetId id() const
Definition: CaloTower.h:100
unsigned int numProblematicHcalCells() const
Definition: CaloTower.h:199
const T & get() const
Definition: EventSetup.h:56
bool CompareTowers(const CaloTower *x, const CaloTower *y)
Definition: HcalHaloAlgo.cc:22
T const * product() const
Definition: ESHandle.h:86
float HBRecHitEnergyThreshold
Definition: HcalHaloAlgo.h:85
float SumEnergyThreshold
Definition: HcalHaloAlgo.h:89
const std::vector< HaloTowerStrip > & getProblematicStrips() const
Definition: HcalHaloData.h:52
bool CompareTime(const EcalRecHit *x, const EcalRecHit *y)
Definition: EcalHaloAlgo.cc:17
const std::vector< PhiWedge > & GetPhiWedges() const
Definition: HcalHaloData.h:48
int NHitsThreshold
Definition: HcalHaloAlgo.h:90
static uInt32 F(BLOWFISH_CTX *ctx, uInt32 x)
Definition: blowfish.cc:281
double emEt() const
Definition: CaloTower.h:112
const CaloGeometry * geo
Definition: HcalHaloAlgo.h:92
std::vector< reco::HaloClusterCandidateHCAL > GetHaloClusterCandidateHB(edm::Handle< EcalRecHitCollection > &ebrechitcoll, edm::Handle< HBHERecHitCollection > &hbherechitcoll, float et_thresh_seedrh)
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 36 of file HcalHaloAlgo.cc.

36  {
38  return Calculate(TheCaloGeometry, TheHBHERecHits, TheCaloTowers,TheEBRecHits,TheEERecHits,TheSetup);
39 }
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)
Definition: HcalHaloAlgo.cc:41
std::vector< HaloClusterCandidateHCAL > HcalHaloAlgo::GetHaloClusterCandidateHB ( edm::Handle< EcalRecHitCollection > &  ebrechitcoll,
edm::Handle< HBHERecHitCollection > &  hbherechitcoll,
float  et_thresh_seedrh 
)

Definition at line 263 of file HcalHaloAlgo.cc.

References funct::abs(), SiPixelRawToDigiRegional_cfi::deltaPhi, CaloRecHit::energy(), EcalRecHit::energy(), eta, HBHERecHit::id(), EcalRecHit::id(), phi(), edm::RefVector< C, T, F >::push_back(), reco::HaloClusterCandidateHCAL::setBeamHaloRecHitsCandidates(), reco::HaloClusterCandidateHCAL::setClusterEt(), reco::HaloClusterCandidateHCAL::setEoverH(), reco::HaloClusterCandidateHCAL::setEtStripPhiSeedMinus1(), reco::HaloClusterCandidateHCAL::setEtStripPhiSeedPlus1(), reco::HaloClusterCandidateHCAL::setIsHaloFromPattern(), reco::HaloClusterCandidateHCAL::setIsHaloFromPattern_HLT(), reco::HaloClusterCandidateHCAL::setNbTowersInEta(), reco::HaloClusterCandidateHCAL::setSeedEt(), reco::HaloClusterCandidateHCAL::setSeedEta(), reco::HaloClusterCandidateHCAL::setSeedPhi(), reco::HaloClusterCandidateHCAL::setSeedR(), reco::HaloClusterCandidateHCAL::setSeedTime(), reco::HaloClusterCandidateHCAL::setSeedZ(), reco::HaloClusterCandidateHCAL::setTimeDiscriminatorITBH(), reco::HaloClusterCandidateHCAL::setTimeDiscriminatorOTBH(), mathSSE::sqrt(), and CaloRecHit::time().

263  {
264 
265  std::vector<HaloClusterCandidateHCAL> TheHaloClusterCandsHB;
266 
267  reco::Vertex::Point vtx(0,0,0);
268 
269  for(size_t ihit = 0; ihit<hbherechitcoll->size(); ++ ihit){
270  HaloClusterCandidateHCAL clustercand;
271 
272  const HBHERecHit & rechit = (*hbherechitcoll)[ ihit ];
273  math::XYZPoint rhpos = getPosition(rechit.id(),vtx);
274  //Et condition
275  double rhet = rechit.energy()* sqrt(rhpos.perp2()/rhpos.mag2());
276  if(rhet<et_thresh_seedrh) continue;
277  if(std::abs(rhpos.z())>zseparation_HBHE) continue;
278  double eta = rhpos.eta();
279  double phi = rhpos.phi();
280 
281  bool isiso = true;
282  double etcluster(0);
283  int nbtowerssameeta(0);
284  double timediscriminatorITBH(0),timediscriminatorOTBH(0);
285  double etstrip_phiseedplus1(0), etstrip_phiseedminus1(0);
286 
287  //Building the cluster
289  for(size_t jhit = 0; jhit<hbherechitcoll->size(); ++ jhit){
290  const HBHERecHit & rechitj = (*hbherechitcoll)[ jhit ];
291  HBHERecHitRef rhRef(hbherechitcoll,jhit);
292  math::XYZPoint rhposj = getPosition(rechitj.id(),vtx);
293  double rhetj = rechitj.energy()* sqrt(rhposj.perp2()/rhposj.mag2());
294  if(rhetj<2) continue;
295  if(std::abs(rhposj.z())>zseparation_HBHE) continue;
296  double etaj = rhposj.eta();
297  double phij = rhposj.phi();
298  double deta = eta - etaj;
299  double dphi = deltaPhi(phi,phij);
300  if(std::abs(deta)>0.4) continue;//This means +/-4 towers in eta
301  if(std::abs(dphi)>0.2) continue;//This means +/-2 towers in phi
302  if(std::abs(dphi)>0.1&&std::abs(deta)<0.2){isiso=false;break;}//The strip should be isolated
303  if(std::abs(dphi)>0.1)continue;
304  if(std::abs(dphi)<0.05) nbtowerssameeta++;
305  if(dphi>0.05) etstrip_phiseedplus1+=rhetj;
306  if(dphi<-0.05) etstrip_phiseedminus1+=rhetj;
307 
308  etcluster+=rhetj;
309  //Timing discriminator
310  //We assign a weight to the rechit defined as:
311  //Log10(Et)*f(T,R,Z)
312  //where f(T,R,Z) is the separation curve between halo-like and IP-like times.
313  //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:
314  //dt= ( - sqrt(R^2+z^2) + |z| )/c
315  // For OT beam halo, the time difference is:
316  //dt= ( 25 + sqrt(R^2+z^2) + |z| )/c
317  //only consider the central part of HB as things get hard at large z.
318  //The best fitted value for R leads to 240 cm (IT) and 330 cm (OT)
319  double rhtj = rechitj.time();
320  timediscriminatorITBH+= std::log10( rhetj )* ( rhtj +0.5*(sqrt(240.*240.+rhposj.z()*rhposj.z()) -std::abs(rhposj.z()))/c_cm_per_ns);
321  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) );
322  bhrhcandidates.push_back(rhRef);
323  }
324  //Isolation conditions
325  if(!isiso) continue;
326  if(etstrip_phiseedplus1/etcluster>0.2&& etstrip_phiseedminus1/etcluster>0.2) continue;
327 
328  //Calculate E/H
329  double eoh(0);
330  for(size_t jhit = 0; jhit<ecalrechitcoll->size(); ++ jhit){
331  const EcalRecHit & rechitj = (*ecalrechitcoll)[ jhit ];
332  math::XYZPoint rhposj = getPosition(rechitj.id(),vtx);
333  double rhetj = rechitj.energy()* sqrt(rhposj.perp2()/rhposj.mag2());
334  if(rhetj<2) continue;
335  double etaj = rhposj.eta();
336  double phij = rhposj.phi();
337  if(std::abs(eta-etaj)>0.2) continue;
338  if(std::abs(deltaPhi(phi,phij))>0.2) continue;
339  eoh+=rhetj/etcluster;
340  }
341  //E/H condition
342  if(eoh>0.1) continue;
343 
344 
345  clustercand.setClusterEt(etcluster);
346  clustercand.setSeedEt(rhet);
347  clustercand.setSeedEta(eta);
348  clustercand.setSeedPhi(phi);
349  clustercand.setSeedZ(rhpos.Z());
350  clustercand.setSeedR(sqrt(rhpos.perp2()));
351  clustercand.setSeedTime(rechit.time());
352  clustercand.setEoverH(eoh);
353  clustercand.setNbTowersInEta(nbtowerssameeta);
354  clustercand.setEtStripPhiSeedPlus1(etstrip_phiseedplus1);
355  clustercand.setEtStripPhiSeedMinus1(etstrip_phiseedminus1);
356  clustercand.setTimeDiscriminatorITBH(timediscriminatorITBH);
357  clustercand.setTimeDiscriminatorOTBH(timediscriminatorOTBH);
358  clustercand.setBeamHaloRecHitsCandidates(bhrhcandidates);
359 
360  bool isbeamhalofrompattern = HBClusterShapeandTimeStudy(clustercand,false);
361  clustercand.setIsHaloFromPattern(isbeamhalofrompattern);
362  bool isbeamhalofrompattern_hlt = HBClusterShapeandTimeStudy(clustercand,true);
363  clustercand.setIsHaloFromPattern_HLT(isbeamhalofrompattern_hlt);
364 
365 
366  TheHaloClusterCandsHB.push_back(clustercand);
367  }
368 
369  return TheHaloClusterCandsHB;
370 }
HcalDetId id() const
get the id
Definition: HBHERecHit.h:23
float time() const
Definition: CaloRecHit.h:19
bool HBClusterShapeandTimeStudy(reco::HaloClusterCandidateHCAL hcand, bool ishlt)
float energy() const
Definition: CaloRecHit.h:17
T sqrt(T t)
Definition: SSEVec.h:18
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
float energy() const
Definition: EcalRecHit.h:68
math::XYZPoint Point
point in the space
Definition: Vertex.h:39
DetId id() const
get the id
Definition: EcalRecHit.h:76
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
math::XYZPoint getPosition(const DetId &id, reco::Vertex::Point vtx)
void setBeamHaloRecHitsCandidates(edm::RefVector< HBHERecHitCollection > x)
void push_back(value_type const &ref)
Add a Ref&lt;C, T&gt; to the RefVector.
Definition: RefVector.h:62
std::vector< HaloClusterCandidateHCAL > HcalHaloAlgo::GetHaloClusterCandidateHE ( edm::Handle< EcalRecHitCollection > &  eerechitcoll,
edm::Handle< HBHERecHitCollection > &  hbherechitcoll,
float  et_thresh_seedrh 
)

Definition at line 373 of file HcalHaloAlgo.cc.

References funct::abs(), SiPixelRawToDigiRegional_cfi::deltaPhi, CaloRecHit::energy(), EcalRecHit::energy(), eta, HBHERecHit::id(), EcalRecHit::id(), phi(), edm::RefVector< C, T, F >::push_back(), reco::HaloClusterCandidateHCAL::setBeamHaloRecHitsCandidates(), reco::HaloClusterCandidateHCAL::setClusterEt(), reco::HaloClusterCandidateHCAL::setClusterSize(), reco::HaloClusterCandidateHCAL::setEoverH(), reco::HaloClusterCandidateHCAL::setEtStripPhiSeedMinus1(), reco::HaloClusterCandidateHCAL::setEtStripPhiSeedPlus1(), reco::HaloClusterCandidateHCAL::setH1overH123(), reco::HaloClusterCandidateHCAL::setIsHaloFromPattern(), reco::HaloClusterCandidateHCAL::setIsHaloFromPattern_HLT(), reco::HaloClusterCandidateHCAL::setSeedEt(), reco::HaloClusterCandidateHCAL::setSeedEta(), reco::HaloClusterCandidateHCAL::setSeedPhi(), reco::HaloClusterCandidateHCAL::setSeedR(), reco::HaloClusterCandidateHCAL::setSeedTime(), reco::HaloClusterCandidateHCAL::setSeedZ(), reco::HaloClusterCandidateHCAL::setTimeDiscriminator(), mathSSE::sqrt(), and CaloRecHit::time().

373  {
374 
375  std::vector<HaloClusterCandidateHCAL> TheHaloClusterCandsHE;
376 
377  reco::Vertex::Point vtx(0,0,0);
378 
379  for(size_t ihit = 0; ihit<hbherechitcoll->size(); ++ ihit){
380  HaloClusterCandidateHCAL clustercand;
381 
382  const HBHERecHit & rechit = (*hbherechitcoll)[ ihit ];
383  math::XYZPoint rhpos = getPosition(rechit.id(),vtx);
384  //Et condition
385  double rhet = rechit.energy()* sqrt(rhpos.perp2()/rhpos.mag2());
386  if(rhet<et_thresh_seedrh) continue;
387  if(std::abs(rhpos.z())<zseparation_HBHE) continue;
388  double eta = rhpos.eta();
389  double phi = rhpos.phi();
390  double rhr = sqrt(rhpos.perp2());
391  bool isiso = true;
392  double etcluster(0),hdepth1(0);
393  int clustersize(0);
394  double etstrip_phiseedplus1(0), etstrip_phiseedminus1(0);
395 
396  //Building the cluster
398  for(size_t jhit = 0; jhit<hbherechitcoll->size(); ++ jhit){
399  const HBHERecHit & rechitj = (*hbherechitcoll)[ jhit ];
400  HBHERecHitRef rhRef(hbherechitcoll,jhit);
401  math::XYZPoint rhposj = getPosition(rechitj.id(),vtx);
402  double rhetj = rechitj.energy()* sqrt(rhposj.perp2()/rhposj.mag2());
403  if(rhetj<2) continue;
404  if(std::abs(rhposj.z())<zseparation_HBHE) continue;
405  if(rhpos.z()*rhposj.z()<0) continue;
406  double phij = rhposj.phi();
407  double dphi = deltaPhi(phi,phij);
408  if(std::abs(dphi)>0.4) continue;
409  double rhrj = sqrt(rhposj.perp2());
410  if(std::abs( rhr-rhrj )>50) continue;
411  if(std::abs(dphi)>0.2 ||std::abs( rhr-rhrj )>20 ){isiso=false;break;}//The deposit should be isolated
412  if(dphi>0.05) etstrip_phiseedplus1+=rhetj;
413  if(dphi<-0.05) etstrip_phiseedminus1+=rhetj;
414  clustersize++;
415  etcluster+=rhetj;
416  if(std::abs( rhposj.z())<405 )hdepth1+=rhetj;
417  //No timing condition for now in HE
418  bhrhcandidates.push_back(rhRef);
419  }
420  //Isolation conditions
421  if(!isiso) continue;
422  if(etstrip_phiseedplus1/etcluster>0.1&& etstrip_phiseedminus1/etcluster>0.1) continue;
423 
424  //Calculate E/H
425  double eoh(0);
426  for(size_t jhit = 0; jhit<ecalrechitcoll->size(); ++ jhit){
427  const EcalRecHit & rechitj = (*ecalrechitcoll)[ jhit ];
428  math::XYZPoint rhposj = getPosition(rechitj.id(),vtx);
429  double rhetj = rechitj.energy()* sqrt(rhposj.perp2()/rhposj.mag2());
430  if(rhetj<2) continue;
431  if(rhpos.z()*rhposj.z()<0) continue;
432  double etaj = rhposj.eta();
433  double phij = rhposj.phi();
434  double dr = sqrt( (eta-etaj)*(eta-etaj)+deltaPhi(phi,phij)*deltaPhi(phi,phij));
435  if(dr>0.3) continue;
436 
437  eoh+=rhetj/etcluster;
438  }
439  //E/H condition
440  if(eoh>0.1) continue;
441 
442 
443  clustercand.setClusterEt(etcluster);
444  clustercand.setSeedEt(rhet);
445  clustercand.setSeedEta(eta);
446  clustercand.setSeedPhi(phi);
447  clustercand.setSeedZ(rhpos.Z());
448  clustercand.setSeedR(sqrt(rhpos.perp2()));
449  clustercand.setSeedTime(rechit.time());
450  clustercand.setEoverH(eoh);
451  clustercand.setH1overH123(hdepth1/etcluster);
452  clustercand.setClusterSize(clustersize);
453  clustercand.setEtStripPhiSeedPlus1(etstrip_phiseedplus1);
454  clustercand.setEtStripPhiSeedMinus1(etstrip_phiseedminus1);
455  clustercand.setTimeDiscriminator(0);
456  clustercand.setBeamHaloRecHitsCandidates(bhrhcandidates);
457 
458  bool isbeamhalofrompattern = HEClusterShapeandTimeStudy(clustercand,false);
459  clustercand.setIsHaloFromPattern(isbeamhalofrompattern);
460  bool isbeamhalofrompattern_hlt = HEClusterShapeandTimeStudy(clustercand,true);
461  clustercand.setIsHaloFromPattern_HLT(isbeamhalofrompattern_hlt);
462 
463 
464  TheHaloClusterCandsHE.push_back(clustercand);
465  }
466 
467  return TheHaloClusterCandsHE;
468 }
HcalDetId id() const
get the id
Definition: HBHERecHit.h:23
float time() const
Definition: CaloRecHit.h:19
float energy() const
Definition: CaloRecHit.h:17
T sqrt(T t)
Definition: SSEVec.h:18
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
float energy() const
Definition: EcalRecHit.h:68
math::XYZPoint Point
point in the space
Definition: Vertex.h:39
DetId id() const
get the id
Definition: EcalRecHit.h:76
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
math::XYZPoint getPosition(const DetId &id, reco::Vertex::Point vtx)
void setBeamHaloRecHitsCandidates(edm::RefVector< HBHERecHitCollection > x)
void push_back(value_type const &ref)
Add a Ref&lt;C, T&gt; to the RefVector.
Definition: RefVector.h:62
bool HEClusterShapeandTimeStudy(reco::HaloClusterCandidateHCAL hcand, bool ishlt)
float HcalHaloAlgo::GetHBRecHitEnergyThreshold ( )
inline

Definition at line 68 of file HcalHaloAlgo.h.

References HBRecHitEnergyThreshold.

68 { return HBRecHitEnergyThreshold;}
float HBRecHitEnergyThreshold
Definition: HcalHaloAlgo.h:85
float HcalHaloAlgo::GetHERecHitEnergyThreshold ( )
inline

Definition at line 69 of file HcalHaloAlgo.h.

References HERecHitEnergyThreshold.

69 { return HERecHitEnergyThreshold;}
float HERecHitEnergyThreshold
Definition: HcalHaloAlgo.h:86
float HcalHaloAlgo::GetPhiWedgeEnergyThreshold ( )
inline

Definition at line 72 of file HcalHaloAlgo.h.

References SumEnergyThreshold.

72 { return SumEnergyThreshold;}
float SumEnergyThreshold
Definition: HcalHaloAlgo.h:89
int HcalHaloAlgo::GetPhiWedgeNHitsThreshold ( )
inline

Definition at line 73 of file HcalHaloAlgo.h.

References NHitsThreshold.

73 { return NHitsThreshold;}
int NHitsThreshold
Definition: HcalHaloAlgo.h:90
math::XYZPoint HcalHaloAlgo::getPosition ( const DetId id,
reco::Vertex::Point  vtx 
)
private

Definition at line 521 of file HcalHaloAlgo.cc.

References PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

521  {
522 
523  const GlobalPoint& pos=geo->getPosition(id);
524  math::XYZPoint posV(pos.x() - vtx.x(),pos.y() - vtx.y(),pos.z() - vtx.z());
525  return posV;
526 }
T y() const
Definition: PV3DBase.h:63
T z() const
Definition: PV3DBase.h:64
const GlobalPoint & getPosition(const DetId &id) const
Get the position of a given detector id.
Definition: CaloGeometry.cc:68
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
T x() const
Definition: PV3DBase.h:62
const CaloGeometry * geo
Definition: HcalHaloAlgo.h:92
bool HcalHaloAlgo::HBClusterShapeandTimeStudy ( reco::HaloClusterCandidateHCAL  hcand,
bool  ishlt 
)

Definition at line 472 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().

472  {
473  //Conditions on the central strip size in eta.
474  //For low size, extra conditions on seed et, isolation and cluster timing
475  //Here we target both IT and OT beam halo. Two separate discriminators were built for the two cases.
476 
477  if(hcand.getSeedEt()<10)return false;
478 
479  if(hcand.getNbTowersInEta()<3) return false;
480  //Isolation criteria for very short eta strips
481  if(hcand.getNbTowersInEta()==3 && (hcand.getEtStripPhiSeedPlus1()>0.1 || hcand.getEtStripPhiSeedMinus1()>0.1) ) return false;
482  if(hcand.getNbTowersInEta()<=5 && (hcand.getEtStripPhiSeedPlus1()>0.1 && hcand.getEtStripPhiSeedMinus1()>0.1) ) return false;
483 
484  //Timing conditions for short eta strips
485  if(hcand.getNbTowersInEta()==3 && hcand.getTimeDiscriminatorITBH()>=0.) return false;
486  if(hcand.getNbTowersInEta()<=6 && hcand.getTimeDiscriminatorITBH()>=5. &&hcand.getTimeDiscriminatorOTBH()<0.) return false;
487 
488  //For HLT, only use conditions without timing
489  if(ishlt && hcand.getNbTowersInEta()<7) return false;
490 
491  hcand.setIsHaloFromPattern(true);
492 
493  return true;
494 }
bool HcalHaloAlgo::HEClusterShapeandTimeStudy ( reco::HaloClusterCandidateHCAL  hcand,
bool  ishlt 
)

Definition at line 498 of file HcalHaloAlgo.cc.

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

498  {
499  //Conditions on H1/H123 to spot halo interacting only in one HCAL layer.
500  //For R> about 170cm, HE has only one layer and this condition cannot be applied
501  //Note that for R>170 cm, the halo is in CSC acceptance and will most likely be spotted by the CSC-calo matching method
502  //A method to identify halos interacting in both H1 and H2/H3 at low R is still missing.
503 
504  if(hcand.getSeedEt()<20)return false;
505  if(hcand.getSeedR()>170) return false;
506 
507  if(hcand.getH1overH123()>0.02 &&hcand.getH1overH123()<0.98) return false;
508 
509  //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.
510  //At HLT, one only cares about large deposits from BH that would lead to a MET/SinglePhoton trigger to be fired.
511  //Rising the seed Et threshold at HLT has therefore little impact on the HLT performances but ensures that possible controversial events are still recorded.
512  if(ishlt && hcand.getSeedEt()<50)return false;
513 
514  hcand.setIsHaloFromPattern(true);
515 
516  return true;
517 
518 }
void HcalHaloAlgo::SetPhiWedgeEnergyThreshold ( float  SumE)
inline

Definition at line 63 of file HcalHaloAlgo.h.

References SumEnergyThreshold.

63 { SumEnergyThreshold = SumE ;}
float SumEnergyThreshold
Definition: HcalHaloAlgo.h:89
void HcalHaloAlgo::SetPhiWedgeNHitsThreshold ( int  nhits)
inline

Definition at line 64 of file HcalHaloAlgo.h.

References nhits, and NHitsThreshold.

void HcalHaloAlgo::SetPhiWedgeThresholds ( float  SumE,
int  nhits 
)
inline

Definition at line 65 of file HcalHaloAlgo.h.

References nhits, NHitsThreshold, and SumEnergyThreshold.

Referenced by reco::HcalHaloDataProducer::produce().

void HcalHaloAlgo::SetRecHitEnergyThresholds ( float  HB,
float  HE 
)
inline

Member Data Documentation

const CaloGeometry* HcalHaloAlgo::geo
private

Definition at line 92 of file HcalHaloAlgo.h.

float HcalHaloAlgo::HBRecHitEnergyThreshold
private

Definition at line 85 of file HcalHaloAlgo.h.

Referenced by GetHBRecHitEnergyThreshold(), and SetRecHitEnergyThresholds().

float HcalHaloAlgo::HERecHitEnergyThreshold
private

Definition at line 86 of file HcalHaloAlgo.h.

Referenced by GetHERecHitEnergyThreshold(), and SetRecHitEnergyThresholds().

int HcalHaloAlgo::NHitsThreshold
private
float HcalHaloAlgo::SumEnergyThreshold
private