HFClusterAlgo.cc

Go to the documentation of this file.

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/METReco/interface/HcalCaloFlagLabels.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
 
#include "HFClusterAlgo.h"
 
using namespace std;
using namespace reco;
HFClusterAlgo::HFClusterAlgo() {
  m_isMC = true;  // safest
}
 
class CompareHFCompleteHitET {
public:
  bool operator()(const HFClusterAlgo::HFCompleteHit& h1, const HFClusterAlgo::HFCompleteHit& h2) const {
    return h1.et > h2.et;
  }
};
 
class CompareHFCore {
public:
  bool operator()(const double c1, const double c2) const { return c1 > c2; }
};
 
static int indexByEta(HcalDetId id) { return (id.zside() > 0) ? (id.ietaAbs() - 29 + 13) : (41 - id.ietaAbs()); }
static const double MCMaterialCorrections_3XX[] = {1.000, 1.000, 1.105, 0.970, 0.965, 0.975, 0.956, 0.958, 0.981,
                                                   1.005, 0.986, 1.086, 1.000, 1.000, 1.086, 0.986, 1.005, 0.981,
                                                   0.958, 0.956, 0.975, 0.965, 0.970, 1.105, 1.000, 1.000};
 
void HFClusterAlgo::setup(double minTowerEnergy,
                          double seedThreshold,
                          double maximumSL,
                          double maximumRenergy,
                          bool usePMTflag,
                          bool usePulseflag,
                          bool forcePulseFlagMC,
                          int correctionSet) {
  m_seedThreshold = seedThreshold;
  m_minTowerEnergy = minTowerEnergy;
  m_maximumSL = maximumSL;
  m_usePMTFlag = usePMTflag;
  m_usePulseFlag = usePulseflag;
  m_forcePulseFlagMC = forcePulseFlagMC;
  m_maximumRenergy = maximumRenergy;
  m_correctionSet = correctionSet;
 
  for (int ii = 0; ii < 13; ii++) {
    m_cutByEta.push_back(-1);
    m_seedmnEta.push_back(99);
    m_seedMXeta.push_back(-1);
  }
  for (int ii = 0; ii < 73; ii++) {
    double minphi = 0.0872664 * (ii - 1);
    double maxphi = 0.0872664 * (ii + 1);
    while (minphi < -M_PI)
      minphi += 2 * M_PI;
    while (minphi > M_PI)
      minphi -= 2 * M_PI;
    while (maxphi < -M_PI)
      maxphi += 2 * M_PI;
    while (maxphi > M_PI)
      maxphi -= 2 * M_PI;
    if (ii == 37)
      minphi = -3.1415904;
    m_seedmnPhi.push_back(minphi);
    m_seedMXphi.push_back(maxphi);
  }
 
  // always set all the corrections to one...
  for (int ii = 0; ii < 13 * 2; ii++)
    m_correctionByEta.push_back(1.0);
  if (m_correctionSet == 1) {  // corrections for material from MC
    for (int ii = 0; ii < 13 * 2; ii++)
      m_correctionByEta[ii] = MCMaterialCorrections_3XX[ii];
  }
}
 
void HFClusterAlgo::clusterize(const HFRecHitCollection& hf,
                               const CaloGeometry& geomO,
                               HFEMClusterShapeCollection& clusterShapes,
                               SuperClusterCollection& SuperClusters) {
  const CaloGeometry* geom(&geomO);
  std::vector<HFCompleteHit> protoseeds, seeds;
  HFRecHitCollection::const_iterator j, j2;
  std::vector<HFCompleteHit>::iterator i;
  std::vector<HFCompleteHit>::iterator k;
  int dP, dE, PWrap;
  bool isok = true;
  HFEMClusterShape clusShp;
 
  SuperCluster Sclus;
  bool doCluster = false;
 
  for (j = hf.begin(); j != hf.end(); j++) {
    const int aieta = j->id().ietaAbs();
    int iz = (aieta - 29);
    // only long fibers and not 29,40,41 allowed to be considered as seeds
    if (j->id().depth() != 1)
      continue;
    if (aieta == 40 || aieta == 41 || aieta == 29)
      continue;
 
    if (iz < 0 || iz > 12) {
      edm::LogWarning("HFClusterAlgo") << "Strange invalid HF hit: " << j->id();
      continue;
    }
 
    if (m_cutByEta[iz] < 0) {
      double eta = geom->getPosition(j->id()).eta();
      m_cutByEta[iz] = m_seedThreshold * cosh(eta);  // convert ET to E for this ring
      auto ccg = geom->getGeometry(j->id());
      const CaloCellGeometry::CornersVec& CellCorners = ccg->getCorners();
      for (size_t sc = 0; sc < CellCorners.size(); sc++) {
        if (fabs(CellCorners[sc].z()) < 1200) {
          if (fabs(CellCorners[sc].eta()) < m_seedmnEta[iz])
            m_seedmnEta[iz] = fabs(CellCorners[sc].eta());
          if (fabs(CellCorners[sc].eta()) > m_seedMXeta[iz])
            m_seedMXeta[iz] = fabs(CellCorners[sc].eta());
        }
      }
    }
    double elong = j->energy() * m_correctionByEta[indexByEta(j->id())];
    if (elong > m_cutByEta[iz]) {
      j2 = hf.find(HcalDetId(HcalForward, j->id().ieta(), j->id().iphi(), 2));
      double eshort = (j2 == hf.end()) ? (0) : (j2->energy());
      if (j2 != hf.end())
        eshort *= m_correctionByEta[indexByEta(j2->id())];
      if (((elong - eshort) / (elong + eshort)) > m_maximumSL)
        continue;
      //if ((m_usePMTFlag)&&(j->flagField(4,1))) continue;
      //if ((m_usePulseFlag)&&(j->flagField(1,1))) continue;
      if (isPMTHit(*j))
        continue;
 
      HFCompleteHit ahit;
      double eta = geom->getPosition(j->id()).eta();
      ahit.id = j->id();
      ahit.energy = elong;
      ahit.et = ahit.energy / cosh(eta);
      protoseeds.push_back(ahit);
    }
  }
 
  if (!protoseeds.empty()) {
    std::sort(protoseeds.begin(), protoseeds.end(), CompareHFCompleteHitET());
    for (i = protoseeds.begin(); i != protoseeds.end(); i++) {
      isok = true;
      doCluster = false;
 
      if ((i == protoseeds.begin()) && (isok)) {
        doCluster = true;
      } else {
        // check for overlap with existing clusters
        for (k = seeds.begin(); isok && k != seeds.end(); k++) {  //i->hits, k->seeds
 
          for (dE = -2; dE <= 2; dE++)
            for (dP = -4; dP <= 4; dP += 2) {
              PWrap = k->id.iphi() + dP;
              if (PWrap < 0)
                PWrap += 72;
              if (PWrap > 72)
                PWrap -= 72;
 
              if ((i->id.iphi() == PWrap) && (i->id.ieta() == k->id.ieta() + dE))
                isok = false;
            }
        }
        if (isok) {
          doCluster = true;
        }
      }
      if (doCluster) {
        seeds.push_back(*i);
 
        bool clusterOk = makeCluster(i->id(), hf, geom, clusShp, Sclus);
        if (clusterOk) {  // cluster is _not_ ok if seed is rejected due to other cuts
          clusterShapes.push_back(clusShp);
          SuperClusters.push_back(Sclus);
        }
      }
    }  //end protoseed loop
  }    //end if seeCount
}
 
bool HFClusterAlgo::makeCluster(const HcalDetId& seedid,
                                const HFRecHitCollection& hf,
                                const CaloGeometry* geom,
                                HFEMClusterShape& clusShp,
                                SuperCluster& Sclus) {
  double w = 0;  //sum over all log E's
  double wgt = 0;
  double w_e = 0;  //sum over ieat*energy
  double w_x = 0;
  double w_y = 0;
  double w_z = 0;
  double wp_e = 0;  //sum over iphi*energy
  double e_e = 0;   //nonwieghted eta sum
  double e_ep = 0;  //nonweighted phi sum
 
  double l_3 = 0;  //sum for enenergy in 3x3 long fibers etc.
  double s_3 = 0;
  double l_5 = 0;
  double s_5 = 0;
  double l_1 = 0;
  double s_1 = 0;
  int de, dp, phiWrap;
  double l_1e = 0;
  const GlobalPoint& sp = geom->getPosition(seedid);
  std::vector<double> coreCanid;
  std::vector<double>::const_iterator ci;
  HFRecHitCollection::const_iterator is, il;
  std::vector<DetId> usedHits;
 
  HFRecHitCollection::const_iterator si;
  HcalDetId sid(HcalForward, seedid.ieta(), seedid.iphi(), 1);
  si = hf.find(sid);
 
  bool clusterOk = true;  // assume the best to start...
 
  // lots happens here
  // edge type 1 has 40/41 in 3x3 and 5x5
  bool edge_type1 = seedid.ietaAbs() == 39 && (seedid.iphi() % 4) == 3;
 
  double e_seed = si->energy() * m_correctionByEta[indexByEta(si->id())];
 
  for (de = -2; de <= 2; de++)
    for (dp = -4; dp <= 4; dp += 2) {
      phiWrap = seedid.iphi() + dp;
      if (phiWrap < 0)
        phiWrap += 72;
      if (phiWrap > 72)
        phiWrap -= 72;
 
      /* Handling of phi-width change problems */
      if (edge_type1 && de == seedid.zside()) {
        if (dp == -2) {  // we want it in the 3x3
          phiWrap -= 2;
          if (phiWrap < 0)
            phiWrap += 72;
        } else if (dp == -4) {
          continue;  // but not double counted in 5x5
        }
      }
 
      HcalDetId idl(HcalForward, seedid.ieta() + de, phiWrap, 1);
      HcalDetId ids(HcalForward, seedid.ieta() + de, phiWrap, 2);
 
      il = hf.find(idl);
      is = hf.find(ids);
 
      double e_long = 1.0;
      double e_short = 0.0;
      if (il != hf.end())
        e_long = il->energy() * m_correctionByEta[indexByEta(il->id())];
      if (e_long <= m_minTowerEnergy)
        e_long = 0.0;
      if (is != hf.end())
        e_short = is->energy() * m_correctionByEta[indexByEta(is->id())];
      if (e_short <= m_minTowerEnergy)
        e_short = 0.0;
      double eRatio = (e_long - e_short) / std::max(1.0, (e_long + e_short));
 
      // require S/L > a minimum amount for inclusion
      if ((abs(eRatio) > m_maximumSL) && (std::max(e_long, e_short) > m_maximumRenergy)) {
        if (dp == 0 && de == 0)
          clusterOk = false;  // somehow, the seed is hosed
        continue;
      }
 
      if ((il != hf.end()) && (isPMTHit(*il))) {
        if (dp == 0 && de == 0)
          clusterOk = false;  // somehow, the seed is hosed
        continue;             //continue to next hit, do not include this one in cluster
      }
 
      if (e_long > m_minTowerEnergy && il != hf.end()) {
        // record usage
        usedHits.push_back(idl.rawId());
        // always in the 5x5
        l_5 += e_long;
        // maybe in the 3x3
        if ((de > -2) && (de < 2) && (dp > -4) && (dp < 4)) {
          l_3 += e_long;
 
          // sometimes in the 1x1
          if ((dp == 0) && (de == 0)) {
            l_1 = e_long;
          }
 
          // maybe in the core?
          if ((de > -2) && (de < 2) && (dp > -4) && (dp < 4) && (e_long > (.5 * e_seed))) {
            coreCanid.push_back(e_long);
          }
 
          // position calculation
          const GlobalPoint& p = geom->getPosition(idl);
 
          double d_p = p.phi() - sp.phi();
          while (d_p < -M_PI)
            d_p += 2 * M_PI;
          while (d_p > M_PI)
            d_p -= 2 * M_PI;
          double d_e = p.eta() - sp.eta();
 
          wgt = log((e_long));
          if (wgt > 0) {
            w += wgt;
            w_e += (d_e)*wgt;
            wp_e += (d_p)*wgt;
            e_e += d_e;
            e_ep += d_p;
 
            w_x += (p.x()) * wgt;  //(p.x()-sp.x())*wgt;
            w_y += (p.y()) * wgt;
            w_z += (p.z()) * wgt;
          }
        }
      } else {
        if (dp == 0 && de == 0)
          clusterOk = false;  // somehow, the seed is hosed
      }
 
      if (e_short > m_minTowerEnergy && is != hf.end()) {
        // record usage
        usedHits.push_back(ids.rawId());
        // always in the 5x5
        s_5 += e_short;
        // maybe in the 3x3
        if ((de > -2) && (de < 2) && (dp > -4) && (dp < 4)) {
          s_3 += e_short;
        }
        // sometimes in the 1x1
        if ((dp == 0) && (de == 0)) {
          s_1 = e_short;
        }
      }
    }
 
  if (!clusterOk)
    return false;
 
  //Core sorting done here
  std::sort(coreCanid.begin(), coreCanid.end(), CompareHFCore());
  for (ci = coreCanid.begin(); ci != coreCanid.end(); ci++) {
    if (ci == coreCanid.begin()) {
      l_1e = *ci;
    } else if (*ci > .5 * l_1e) {
      l_1e += *ci;
    }
  }  //core sorting end
 
  double z_ = w_z / w;
  double x_ = w_x / w;
  double y_ = w_y / w;
  math::XYZPoint xyzclus(x_, y_, z_);
  //calcualte position, final
  double eta = xyzclus.eta();  //w_e/w+sp.eta();
 
  double phi = xyzclus.phi();  //(wp_e/w)+sp.phi();
 
  while (phi < -M_PI)
    phi += 2 * M_PI;
  while (phi > M_PI)
    phi -= 2 * M_PI;
 
  //calculate cell phi and cell eta
  int idx = fabs(seedid.ieta()) - 29;
  int ipx = seedid.iphi();
  double Cphi = (phi - m_seedmnPhi[ipx]) / (m_seedMXphi[ipx] - m_seedmnPhi[ipx]);
  double Ceta = (fabs(eta) - m_seedmnEta[idx]) / (m_seedMXeta[idx] - m_seedmnEta[idx]);
 
  //return  HFEMClusterShape, SuperCluster
  HFEMClusterShape myClusShp(l_1, s_1, l_3, s_3, l_5, s_5, l_1e, Ceta, Cphi, seedid);
  clusShp = myClusShp;
 
  SuperCluster MySclus(l_3, xyzclus);
  Sclus = MySclus;
 
  return clusterOk;
}
bool HFClusterAlgo::isPMTHit(const HFRecHit& hfr) {
  bool pmthit = false;
 
  if ((hfr.flagField(HcalCaloFlagLabels::HFLongShort)) && (m_usePMTFlag))
    pmthit = true;
  if (!(m_isMC && !m_forcePulseFlagMC))
    if ((hfr.flagField(HcalCaloFlagLabels::HFDigiTime)) && (m_usePulseFlag))
      pmthit = true;
 
  return pmthit;
}
void HFClusterAlgo::resetForRun() {
  edm::LogInfo("HFClusterAlgo") << "Resetting for Run!";
  for (int ii = 0; ii < 13; ii++) {
    m_cutByEta.push_back(-1);
    m_seedmnEta.push_back(99);
    m_seedMXeta.push_back(-1);
  }
}