CMS 3D CMS Logo

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

#include <HFClusterAlgo.h>

Classes

struct  HFCompleteHit
 

Public Member Functions

void clusterize (const HFRecHitCollection &hf, const CaloGeometry &geom, reco::HFEMClusterShapeCollection &clusters, reco::SuperClusterCollection &SuperClusters)
 
 HFClusterAlgo ()
 
void isMC (bool isMC)
 
void resetForRun ()
 
void setup (double minTowerEnergy, double seedThreshold, double maximumSL, double m_maximumRenergy, bool usePMTflag, bool usePulseflag, bool forcePulseFlagMC, int correctionSet)
 

Private Member Functions

bool isPMTHit (const HFRecHit &hfr)
 
bool makeCluster (const HcalDetId &seedid, const HFRecHitCollection &hf, const CaloGeometry *geom, reco::HFEMClusterShape &clusShp, reco::SuperCluster &SClus)
 

Private Attributes

std::vector< double > m_correctionByEta
 
int m_correctionSet
 
std::vector< double > m_cutByEta
 
bool m_forcePulseFlagMC
 
bool m_isMC
 
double m_maximumRenergy
 
double m_maximumSL
 
double m_minTowerEnergy
 
std::vector< double > m_seedmnEta
 
std::vector< double > m_seedmnPhi
 
std::vector< double > m_seedMXeta
 
std::vector< double > m_seedMXphi
 
double m_seedThreshold
 
bool m_usePMTFlag
 
bool m_usePulseFlag
 

Friends

class CompareHFCompleteHitET
 
class CompareHFCore
 

Detailed Description

Author
Kevin Klapoetke (Minnesota)

$Id:version 1.2

Author
K. Klapoetke – Minnesota

Definition at line 18 of file HFClusterAlgo.h.

Constructor & Destructor Documentation

◆ HFClusterAlgo()

HFClusterAlgo::HFClusterAlgo ( )

Definition at line 16 of file HFClusterAlgo.cc.

16  {
17  m_isMC = true; // safest
18 }

Member Function Documentation

◆ clusterize()

void HFClusterAlgo::clusterize ( const HFRecHitCollection hf,
const CaloGeometry geomO,
reco::HFEMClusterShapeCollection clusters,
reco::SuperClusterCollection SuperClusters 
)

Analyze the hits

Definition at line 86 of file HFClusterAlgo.cc.

References HFClusterAlgo::HFCompleteHit::energy, HFClusterAlgo::HFCompleteHit::et, PVValHelper::eta, relativeConstraints::geom, HcalForward, photonIsolationHIProducer_cfi::hf, mps_fire::i, HFClusterAlgo::HFCompleteHit::id, indexByEta(), dqmiolumiharvest::j, dqmdumpme::k, DetachedQuadStep_cff::seeds, EZArrayFL< T >::size(), and jetUpdater_cfi::sort.

Referenced by HFEMClusterProducer::produce().

89  {
90  const CaloGeometry* geom(&geomO);
91  std::vector<HFCompleteHit> protoseeds, seeds;
93  std::vector<HFCompleteHit>::iterator i;
94  std::vector<HFCompleteHit>::iterator k;
95  int dP, dE, PWrap;
96  bool isok = true;
97  HFEMClusterShape clusShp;
98 
99  SuperCluster Sclus;
100  bool doCluster = false;
101 
102  for (j = hf.begin(); j != hf.end(); j++) {
103  const int aieta = j->id().ietaAbs();
104  int iz = (aieta - 29);
105  // only long fibers and not 29,40,41 allowed to be considered as seeds
106  if (j->id().depth() != 1)
107  continue;
108  if (aieta == 40 || aieta == 41 || aieta == 29)
109  continue;
110 
111  if (iz < 0 || iz > 12) {
112  edm::LogWarning("HFClusterAlgo") << "Strange invalid HF hit: " << j->id();
113  continue;
114  }
115 
116  if (m_cutByEta[iz] < 0) {
117  double eta = geom->getPosition(j->id()).eta();
118  m_cutByEta[iz] = m_seedThreshold * cosh(eta); // convert ET to E for this ring
119  auto ccg = geom->getGeometry(j->id());
120  const CaloCellGeometry::CornersVec& CellCorners = ccg->getCorners();
121  for (size_t sc = 0; sc < CellCorners.size(); sc++) {
122  if (fabs(CellCorners[sc].z()) < 1200) {
123  if (fabs(CellCorners[sc].eta()) < m_seedmnEta[iz])
124  m_seedmnEta[iz] = fabs(CellCorners[sc].eta());
125  if (fabs(CellCorners[sc].eta()) > m_seedMXeta[iz])
126  m_seedMXeta[iz] = fabs(CellCorners[sc].eta());
127  }
128  }
129  }
130  double elong = j->energy() * m_correctionByEta[indexByEta(j->id())];
131  if (elong > m_cutByEta[iz]) {
132  j2 = hf.find(HcalDetId(HcalForward, j->id().ieta(), j->id().iphi(), 2));
133  double eshort = (j2 == hf.end()) ? (0) : (j2->energy());
134  if (j2 != hf.end())
135  eshort *= m_correctionByEta[indexByEta(j2->id())];
136  if (((elong - eshort) / (elong + eshort)) > m_maximumSL)
137  continue;
138  //if ((m_usePMTFlag)&&(j->flagField(4,1))) continue;
139  //if ((m_usePulseFlag)&&(j->flagField(1,1))) continue;
140  if (isPMTHit(*j))
141  continue;
142 
143  HFCompleteHit ahit;
144  double eta = geom->getPosition(j->id()).eta();
145  ahit.id = j->id();
146  ahit.energy = elong;
147  ahit.et = ahit.energy / cosh(eta);
148  protoseeds.push_back(ahit);
149  }
150  }
151 
152  if (!protoseeds.empty()) {
153  std::sort(protoseeds.begin(), protoseeds.end(), CompareHFCompleteHitET());
154  for (i = protoseeds.begin(); i != protoseeds.end(); i++) {
155  isok = true;
156  doCluster = false;
157 
158  if ((i == protoseeds.begin()) && (isok)) {
159  doCluster = true;
160  } else {
161  // check for overlap with existing clusters
162  for (k = seeds.begin(); isok && k != seeds.end(); k++) { //i->hits, k->seeds
163 
164  for (dE = -2; dE <= 2; dE++)
165  for (dP = -4; dP <= 4; dP += 2) {
166  PWrap = k->id.iphi() + dP;
167  if (PWrap < 0)
168  PWrap += 72;
169  if (PWrap > 72)
170  PWrap -= 72;
171 
172  if ((i->id.iphi() == PWrap) && (i->id.ieta() == k->id.ieta() + dE))
173  isok = false;
174  }
175  }
176  if (isok) {
177  doCluster = true;
178  }
179  }
180  if (doCluster) {
181  seeds.push_back(*i);
182 
183  bool clusterOk = makeCluster(i->id(), hf, geom, clusShp, Sclus);
184  if (clusterOk) { // cluster is _not_ ok if seed is rejected due to other cuts
185  clusterShapes.push_back(clusShp);
186  SuperClusters.push_back(Sclus);
187  }
188  }
189  } //end protoseed loop
190  } //end if seeCount
191 }
size_type size() const
Definition: EZArrayFL.h:67
std::vector< double > m_cutByEta
Definition: HFClusterAlgo.h:50
double m_maximumSL
Definition: HFClusterAlgo.h:45
std::vector< double > m_correctionByEta
Definition: HFClusterAlgo.h:51
std::vector< T >::const_iterator const_iterator
std::vector< double > m_seedmnEta
Definition: HFClusterAlgo.h:52
friend class CompareHFCompleteHitET
Definition: HFClusterAlgo.h:42
std::vector< double > m_seedMXeta
Definition: HFClusterAlgo.h:53
static int indexByEta(HcalDetId id)
double m_seedThreshold
Definition: HFClusterAlgo.h:45
Log< level::Warning, false > LogWarning
bool makeCluster(const HcalDetId &seedid, const HFRecHitCollection &hf, const CaloGeometry *geom, reco::HFEMClusterShape &clusShp, reco::SuperCluster &SClus)
bool isPMTHit(const HFRecHit &hfr)

◆ isMC()

void HFClusterAlgo::isMC ( bool  isMC)
inline

Definition at line 31 of file HFClusterAlgo.h.

References isMC(), and m_isMC.

Referenced by isMC(), and HFEMClusterProducer::produce().

31 { m_isMC = isMC; }
void isMC(bool isMC)
Definition: HFClusterAlgo.h:31

◆ isPMTHit()

bool HFClusterAlgo::isPMTHit ( const HFRecHit hfr)
private

Definition at line 380 of file HFClusterAlgo.cc.

References CaloRecHit::flagField(), HcalCaloFlagLabels::HFDigiTime, and HcalCaloFlagLabels::HFLongShort.

380  {
381  bool pmthit = false;
382 
384  pmthit = true;
385  if (!(m_isMC && !m_forcePulseFlagMC))
387  pmthit = true;
388 
389  return pmthit;
390 }
bool m_forcePulseFlagMC
Definition: HFClusterAlgo.h:47
constexpr uint32_t flagField(int base, int width=1) const
Definition: CaloRecHit.h:46

◆ makeCluster()

bool HFClusterAlgo::makeCluster ( const HcalDetId seedid,
const HFRecHitCollection hf,
const CaloGeometry geom,
reco::HFEMClusterShape clusShp,
reco::SuperCluster SClus 
)
private

Definition at line 193 of file HFClusterAlgo.cc.

References funct::abs(), Calorimetry_cff::dp, PVValHelper::eta, relativeConstraints::geom, HcalForward, photonIsolationHIProducer_cfi::hf, heavyIonCSV_trainingSettings::idx, HcalDetId::ieta(), HcalDetId::ietaAbs(), indexByEta(), HcalDetId::iphi(), dqm-mbProfile::log, M_PI, SiStripPI::max, AlCaHLTBitMon_ParallelJobs::p, PV3DBase< T, PVType, FrameType >::phi(), jetUpdater_cfi::sort, w(), and HcalDetId::zside().

197  {
198  double w = 0; //sum over all log E's
199  double wgt = 0;
200  double w_x = 0;
201  double w_y = 0;
202  double w_z = 0;
203 
204  double l_3 = 0; //sum for enenergy in 3x3 long fibers etc.
205  double s_3 = 0;
206  double l_5 = 0;
207  double s_5 = 0;
208  double l_1 = 0;
209  double s_1 = 0;
210  int de, dp, phiWrap;
211  double l_1e = 0;
212  const GlobalPoint& sp = geom->getPosition(seedid);
213  std::vector<double> coreCanid;
214  std::vector<double>::const_iterator ci;
216  std::vector<DetId> usedHits;
217 
219  HcalDetId sid(HcalForward, seedid.ieta(), seedid.iphi(), 1);
220  si = hf.find(sid);
221 
222  bool clusterOk = true; // assume the best to start...
223 
224  // lots happens here
225  // edge type 1 has 40/41 in 3x3 and 5x5
226  bool edge_type1 = seedid.ietaAbs() == 39 && (seedid.iphi() % 4) == 3;
227 
228  double e_seed = si->energy() * m_correctionByEta[indexByEta(si->id())];
229 
230  for (de = -2; de <= 2; de++)
231  for (dp = -4; dp <= 4; dp += 2) {
232  phiWrap = seedid.iphi() + dp;
233  if (phiWrap < 0)
234  phiWrap += 72;
235  if (phiWrap > 72)
236  phiWrap -= 72;
237 
238  /* Handling of phi-width change problems */
239  if (edge_type1 && de == seedid.zside()) {
240  if (dp == -2) { // we want it in the 3x3
241  phiWrap -= 2;
242  if (phiWrap < 0)
243  phiWrap += 72;
244  } else if (dp == -4) {
245  continue; // but not double counted in 5x5
246  }
247  }
248 
249  HcalDetId idl(HcalForward, seedid.ieta() + de, phiWrap, 1);
250  HcalDetId ids(HcalForward, seedid.ieta() + de, phiWrap, 2);
251 
252  il = hf.find(idl);
253  is = hf.find(ids);
254 
255  double e_long = 1.0;
256  double e_short = 0.0;
257  if (il != hf.end())
258  e_long = il->energy() * m_correctionByEta[indexByEta(il->id())];
259  if (e_long <= m_minTowerEnergy)
260  e_long = 0.0;
261  if (is != hf.end())
262  e_short = is->energy() * m_correctionByEta[indexByEta(is->id())];
263  if (e_short <= m_minTowerEnergy)
264  e_short = 0.0;
265  double eRatio = (e_long - e_short) / std::max(1.0, (e_long + e_short));
266 
267  // require S/L > a minimum amount for inclusion
268  if ((abs(eRatio) > m_maximumSL) && (std::max(e_long, e_short) > m_maximumRenergy)) {
269  if (dp == 0 && de == 0)
270  clusterOk = false; // somehow, the seed is hosed
271  continue;
272  }
273 
274  if ((il != hf.end()) && (isPMTHit(*il))) {
275  if (dp == 0 && de == 0)
276  clusterOk = false; // somehow, the seed is hosed
277  continue; //continue to next hit, do not include this one in cluster
278  }
279 
280  if (e_long > m_minTowerEnergy && il != hf.end()) {
281  // record usage
282  usedHits.push_back(idl.rawId());
283  // always in the 5x5
284  l_5 += e_long;
285  // maybe in the 3x3
286  if ((de > -2) && (de < 2) && (dp > -4) && (dp < 4)) {
287  l_3 += e_long;
288 
289  // sometimes in the 1x1
290  if ((dp == 0) && (de == 0)) {
291  l_1 = e_long;
292  }
293 
294  // maybe in the core?
295  if ((de > -2) && (de < 2) && (dp > -4) && (dp < 4) && (e_long > (.5 * e_seed))) {
296  coreCanid.push_back(e_long);
297  }
298 
299  // position calculation
300  const GlobalPoint& p = geom->getPosition(idl);
301 
302  double d_p = p.phi() - sp.phi();
303  while (d_p < -M_PI)
304  d_p += 2 * M_PI;
305  while (d_p > M_PI)
306  d_p -= 2 * M_PI;
307 
308  wgt = log((e_long));
309  if (wgt > 0) {
310  w += wgt;
311 
312  w_x += (p.x()) * wgt; //(p.x()-sp.x())*wgt;
313  w_y += (p.y()) * wgt;
314  w_z += (p.z()) * wgt;
315  }
316  }
317  } else {
318  if (dp == 0 && de == 0)
319  clusterOk = false; // somehow, the seed is hosed
320  }
321 
322  if (e_short > m_minTowerEnergy && is != hf.end()) {
323  // record usage
324  usedHits.push_back(ids.rawId());
325  // always in the 5x5
326  s_5 += e_short;
327  // maybe in the 3x3
328  if ((de > -2) && (de < 2) && (dp > -4) && (dp < 4)) {
329  s_3 += e_short;
330  }
331  // sometimes in the 1x1
332  if ((dp == 0) && (de == 0)) {
333  s_1 = e_short;
334  }
335  }
336  }
337 
338  if (!clusterOk)
339  return false;
340 
341  //Core sorting done here
342  std::sort(coreCanid.begin(), coreCanid.end(), CompareHFCore());
343  for (ci = coreCanid.begin(); ci != coreCanid.end(); ci++) {
344  if (ci == coreCanid.begin()) {
345  l_1e = *ci;
346  } else if (*ci > .5 * l_1e) {
347  l_1e += *ci;
348  }
349  } //core sorting end
350 
351  double z_ = w_z / w;
352  double x_ = w_x / w;
353  double y_ = w_y / w;
354  math::XYZPoint xyzclus(x_, y_, z_);
355  //calcualte position, final
356  double eta = xyzclus.eta(); //w_e/w+sp.eta();
357 
358  double phi = xyzclus.phi(); //(wp_e/w)+sp.phi();
359 
360  while (phi < -M_PI)
361  phi += 2 * M_PI;
362  while (phi > M_PI)
363  phi -= 2 * M_PI;
364 
365  //calculate cell phi and cell eta
366  int idx = fabs(seedid.ieta()) - 29;
367  int ipx = seedid.iphi();
368  double Cphi = (phi - m_seedmnPhi[ipx]) / (m_seedMXphi[ipx] - m_seedmnPhi[ipx]);
369  double Ceta = (fabs(eta) - m_seedmnEta[idx]) / (m_seedMXeta[idx] - m_seedmnEta[idx]);
370 
371  //return HFEMClusterShape, SuperCluster
372  HFEMClusterShape myClusShp(l_1, s_1, l_3, s_3, l_5, s_5, l_1e, Ceta, Cphi, seedid);
373  clusShp = myClusShp;
374 
375  SuperCluster MySclus(l_3, xyzclus);
376  Sclus = MySclus;
377 
378  return clusterOk;
379 }
double m_minTowerEnergy
Definition: HFClusterAlgo.h:45
double m_maximumRenergy
Definition: HFClusterAlgo.h:45
constexpr int zside() const
get the z-side of the cell (1/-1)
Definition: HcalDetId.h:141
T w() const
constexpr int ietaAbs() const
get the absolute value of the cell ieta
Definition: HcalDetId.h:148
double m_maximumSL
Definition: HFClusterAlgo.h:45
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
std::vector< double > m_correctionByEta
Definition: HFClusterAlgo.h:51
std::vector< T >::const_iterator const_iterator
std::vector< double > m_seedmnEta
Definition: HFClusterAlgo.h:52
friend class CompareHFCore
Definition: HFClusterAlgo.h:43
std::vector< double > m_seedMXeta
Definition: HFClusterAlgo.h:53
constexpr int ieta() const
get the cell ieta
Definition: HcalDetId.h:155
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< double > m_seedmnPhi
Definition: HFClusterAlgo.h:54
#define M_PI
static int indexByEta(HcalDetId id)
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
std::vector< double > m_seedMXphi
Definition: HFClusterAlgo.h:55
bool isPMTHit(const HFRecHit &hfr)
constexpr int iphi() const
get the cell iphi
Definition: HcalDetId.h:157

◆ resetForRun()

void HFClusterAlgo::resetForRun ( )

Definition at line 391 of file HFClusterAlgo.cc.

References cuy::ii.

Referenced by HFEMClusterProducer::beginRun().

391  {
392  edm::LogInfo("HFClusterAlgo") << "Resetting for Run!";
393  for (int ii = 0; ii < 13; ii++) {
394  m_cutByEta.push_back(-1);
395  m_seedmnEta.push_back(99);
396  m_seedMXeta.push_back(-1);
397  }
398 }
std::vector< double > m_cutByEta
Definition: HFClusterAlgo.h:50
std::vector< double > m_seedmnEta
Definition: HFClusterAlgo.h:52
std::vector< double > m_seedMXeta
Definition: HFClusterAlgo.h:53
ii
Definition: cuy.py:589
Log< level::Info, false > LogInfo

◆ setup()

void HFClusterAlgo::setup ( double  minTowerEnergy,
double  seedThreshold,
double  maximumSL,
double  m_maximumRenergy,
bool  usePMTflag,
bool  usePulseflag,
bool  forcePulseFlagMC,
int  correctionSet 
)

Definition at line 37 of file HFClusterAlgo.cc.

References hfClusterShapes_cfi::forcePulseFlagMC, cuy::ii, M_PI, hfClusterShapes_cfi::maximumRenergy, hfClusterShapes_cfi::maximumSL, MCMaterialCorrections_3XX, and hfClusterShapes_cfi::minTowerEnergy.

44  {
45  m_seedThreshold = seedThreshold;
48  m_usePMTFlag = usePMTflag;
49  m_usePulseFlag = usePulseflag;
52  m_correctionSet = correctionSet;
53 
54  for (int ii = 0; ii < 13; ii++) {
55  m_cutByEta.push_back(-1);
56  m_seedmnEta.push_back(99);
57  m_seedMXeta.push_back(-1);
58  }
59  for (int ii = 0; ii < 73; ii++) {
60  double minphi = 0.0872664 * (ii - 1);
61  double maxphi = 0.0872664 * (ii + 1);
62  while (minphi < -M_PI)
63  minphi += 2 * M_PI;
64  while (minphi > M_PI)
65  minphi -= 2 * M_PI;
66  while (maxphi < -M_PI)
67  maxphi += 2 * M_PI;
68  while (maxphi > M_PI)
69  maxphi -= 2 * M_PI;
70  if (ii == 37)
71  minphi = -3.1415904;
72  m_seedmnPhi.push_back(minphi);
73  m_seedMXphi.push_back(maxphi);
74  }
75 
76  // always set all the corrections to one...
77  for (int ii = 0; ii < 13 * 2; ii++)
78  m_correctionByEta.push_back(1.0);
79  if (m_correctionSet == 1) { // corrections for material from MC
80  for (int ii = 0; ii < 13 * 2; ii++)
82  }
83 }
double m_minTowerEnergy
Definition: HFClusterAlgo.h:45
double m_maximumRenergy
Definition: HFClusterAlgo.h:45
bool m_forcePulseFlagMC
Definition: HFClusterAlgo.h:47
std::vector< double > m_cutByEta
Definition: HFClusterAlgo.h:50
double m_maximumSL
Definition: HFClusterAlgo.h:45
std::vector< double > m_correctionByEta
Definition: HFClusterAlgo.h:51
std::vector< double > m_seedmnEta
Definition: HFClusterAlgo.h:52
static const double MCMaterialCorrections_3XX[]
std::vector< double > m_seedMXeta
Definition: HFClusterAlgo.h:53
std::vector< double > m_seedmnPhi
Definition: HFClusterAlgo.h:54
ii
Definition: cuy.py:589
#define M_PI
std::vector< double > m_seedMXphi
Definition: HFClusterAlgo.h:55
double m_seedThreshold
Definition: HFClusterAlgo.h:45

Friends And Related Function Documentation

◆ CompareHFCompleteHitET

friend class CompareHFCompleteHitET
friend

Definition at line 42 of file HFClusterAlgo.h.

◆ CompareHFCore

friend class CompareHFCore
friend

Definition at line 43 of file HFClusterAlgo.h.

Member Data Documentation

◆ m_correctionByEta

std::vector<double> HFClusterAlgo::m_correctionByEta
private

Definition at line 51 of file HFClusterAlgo.h.

◆ m_correctionSet

int HFClusterAlgo::m_correctionSet
private

Definition at line 49 of file HFClusterAlgo.h.

◆ m_cutByEta

std::vector<double> HFClusterAlgo::m_cutByEta
private

Definition at line 50 of file HFClusterAlgo.h.

◆ m_forcePulseFlagMC

bool HFClusterAlgo::m_forcePulseFlagMC
private

Definition at line 47 of file HFClusterAlgo.h.

◆ m_isMC

bool HFClusterAlgo::m_isMC
private

Definition at line 48 of file HFClusterAlgo.h.

Referenced by isMC().

◆ m_maximumRenergy

double HFClusterAlgo::m_maximumRenergy
private

Definition at line 45 of file HFClusterAlgo.h.

◆ m_maximumSL

double HFClusterAlgo::m_maximumSL
private

Definition at line 45 of file HFClusterAlgo.h.

◆ m_minTowerEnergy

double HFClusterAlgo::m_minTowerEnergy
private

Definition at line 45 of file HFClusterAlgo.h.

◆ m_seedmnEta

std::vector<double> HFClusterAlgo::m_seedmnEta
private

Definition at line 52 of file HFClusterAlgo.h.

◆ m_seedmnPhi

std::vector<double> HFClusterAlgo::m_seedmnPhi
private

Definition at line 54 of file HFClusterAlgo.h.

◆ m_seedMXeta

std::vector<double> HFClusterAlgo::m_seedMXeta
private

Definition at line 53 of file HFClusterAlgo.h.

◆ m_seedMXphi

std::vector<double> HFClusterAlgo::m_seedMXphi
private

Definition at line 55 of file HFClusterAlgo.h.

◆ m_seedThreshold

double HFClusterAlgo::m_seedThreshold
private

Definition at line 45 of file HFClusterAlgo.h.

◆ m_usePMTFlag

bool HFClusterAlgo::m_usePMTFlag
private

Definition at line 46 of file HFClusterAlgo.h.

◆ m_usePulseFlag

bool HFClusterAlgo::m_usePulseFlag
private

Definition at line 47 of file HFClusterAlgo.h.