EgammaHFProducers/plugins/HFClusterAlgo.cc

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "DataFormats/METReco/interface/HcalCaloFlagLabels.h"
0003 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0004
0005 #include "HFClusterAlgo.h"
0006
0007 using namespace std;
0008 using namespace reco;
0009 /** \class HFClusterAlgo
0010  *
0011  *  \author Kevin Klapoetke (Minnesota)
0012  *
0013  * $Id:version 1.2
0014  */
0015
0016 HFClusterAlgo::HFClusterAlgo() {
0017   m_isMC = true;  // safest
0018 }
0019
0020 class CompareHFCompleteHitET {
0021 public:
0022   bool operator()(const HFClusterAlgo::HFCompleteHit& h1, const HFClusterAlgo::HFCompleteHit& h2) const {
0023     return h1.et > h2.et;
0024   }
0025 };
0026
0027 class CompareHFCore {
0028 public:
0029   bool operator()(const double c1, const double c2) const { return c1 > c2; }
0030 };
0031
0032 static int indexByEta(HcalDetId id) { return (id.zside() > 0) ? (id.ietaAbs() - 29 + 13) : (41 - id.ietaAbs()); }
0033 static const double MCMaterialCorrections_3XX[] = {1.000, 1.000, 1.105, 0.970, 0.965, 0.975, 0.956, 0.958, 0.981,
0034                                                    1.005, 0.986, 1.086, 1.000, 1.000, 1.086, 0.986, 1.005, 0.981,
0035                                                    0.958, 0.956, 0.975, 0.965, 0.970, 1.105, 1.000, 1.000};
0036
0037 void HFClusterAlgo::setup(double minTowerEnergy,
0038                           double seedThreshold,
0039                           double maximumSL,
0040                           double maximumRenergy,
0041                           bool usePMTflag,
0042                           bool usePulseflag,
0043                           bool forcePulseFlagMC,
0044                           int correctionSet) {
0045   m_seedThreshold = seedThreshold;
0046   m_minTowerEnergy = minTowerEnergy;
0047   m_maximumSL = maximumSL;
0048   m_usePMTFlag = usePMTflag;
0049   m_usePulseFlag = usePulseflag;
0050   m_forcePulseFlagMC = forcePulseFlagMC;
0051   m_maximumRenergy = maximumRenergy;
0052   m_correctionSet = correctionSet;
0053
0054   for (int ii = 0; ii < 13; ii++) {
0055     m_cutByEta.push_back(-1);
0056     m_seedmnEta.push_back(99);
0057     m_seedMXeta.push_back(-1);
0058   }
0059   for (int ii = 0; ii < 73; ii++) {
0060     double minphi = 0.0872664 * (ii - 1);
0061     double maxphi = 0.0872664 * (ii + 1);
0062     while (minphi < -M_PI)
0063       minphi += 2 * M_PI;
0064     while (minphi > M_PI)
0065       minphi -= 2 * M_PI;
0066     while (maxphi < -M_PI)
0067       maxphi += 2 * M_PI;
0068     while (maxphi > M_PI)
0069       maxphi -= 2 * M_PI;
0070     if (ii == 37)
0071       minphi = -3.1415904;
0072     m_seedmnPhi.push_back(minphi);
0073     m_seedMXphi.push_back(maxphi);
0074   }
0075
0076   // always set all the corrections to one...
0077   for (int ii = 0; ii < 13 * 2; ii++)
0078     m_correctionByEta.push_back(1.0);
0079   if (m_correctionSet == 1) {  // corrections for material from MC
0080     for (int ii = 0; ii < 13 * 2; ii++)
0081       m_correctionByEta[ii] = MCMaterialCorrections_3XX[ii];
0082   }
0083 }
0084
0085 /** Analyze the hits */
0086 void HFClusterAlgo::clusterize(const HFRecHitCollection& hf,
0087                                const CaloGeometry& geomO,
0088                                HFEMClusterShapeCollection& clusterShapes,
0089                                SuperClusterCollection& SuperClusters) {
0090   const CaloGeometry* geom(&geomO);
0091   std::vector<HFCompleteHit> protoseeds, seeds;
0092   HFRecHitCollection::const_iterator j, j2;
0093   std::vector<HFCompleteHit>::iterator i;
0094   std::vector<HFCompleteHit>::iterator k;
0095   int dP, dE, PWrap;
0096   bool isok = true;
0097   HFEMClusterShape clusShp;
0098
0099   SuperCluster Sclus;
0100   bool doCluster = false;
0101
0102   for (j = hf.begin(); j != hf.end(); j++) {
0103     const int aieta = j->id().ietaAbs();
0104     int iz = (aieta - 29);
0105     // only long fibers and not 29,40,41 allowed to be considered as seeds
0106     if (j->id().depth() != 1)
0107       continue;
0108     if (aieta == 40 || aieta == 41 || aieta == 29)
0109       continue;
0110
0111     if (iz < 0 || iz > 12) {
0112       edm::LogWarning("HFClusterAlgo") << "Strange invalid HF hit: " << j->id();
0113       continue;
0114     }
0115
0116     if (m_cutByEta[iz] < 0) {
0117       double eta = geom->getPosition(j->id()).eta();
0118       m_cutByEta[iz] = m_seedThreshold * cosh(eta);  // convert ET to E for this ring
0119       auto ccg = geom->getGeometry(j->id());
0120       const CaloCellGeometry::CornersVec& CellCorners = ccg->getCorners();
0121       for (size_t sc = 0; sc < CellCorners.size(); sc++) {
0122         if (fabs(CellCorners[sc].z()) < 1200) {
0123           if (fabs(CellCorners[sc].eta()) < m_seedmnEta[iz])
0124             m_seedmnEta[iz] = fabs(CellCorners[sc].eta());
0125           if (fabs(CellCorners[sc].eta()) > m_seedMXeta[iz])
0126             m_seedMXeta[iz] = fabs(CellCorners[sc].eta());
0127         }
0128       }
0129     }
0130     double elong = j->energy() * m_correctionByEta[indexByEta(j->id())];
0131     if (elong > m_cutByEta[iz]) {
0132       j2 = hf.find(HcalDetId(HcalForward, j->id().ieta(), j->id().iphi(), 2));
0133       double eshort = (j2 == hf.end()) ? (0) : (j2->energy());
0134       if (j2 != hf.end())
0135         eshort *= m_correctionByEta[indexByEta(j2->id())];
0136       if (((elong - eshort) / (elong + eshort)) > m_maximumSL)
0137         continue;
0138       //if ((m_usePMTFlag)&&(j->flagField(4,1))) continue;
0139       //if ((m_usePulseFlag)&&(j->flagField(1,1))) continue;
0140       if (isPMTHit(*j))
0141         continue;
0142
0143       HFCompleteHit ahit;
0144       double eta = geom->getPosition(j->id()).eta();
0145       ahit.id = j->id();
0146       ahit.energy = elong;
0147       ahit.et = ahit.energy / cosh(eta);
0148       protoseeds.push_back(ahit);
0149     }
0150   }
0151
0152   if (!protoseeds.empty()) {
0153     std::sort(protoseeds.begin(), protoseeds.end(), CompareHFCompleteHitET());
0154     for (i = protoseeds.begin(); i != protoseeds.end(); i++) {
0155       isok = true;
0156       doCluster = false;
0157
0158       if ((i == protoseeds.begin()) && (isok)) {
0159         doCluster = true;
0160       } else {
0161         // check for overlap with existing clusters
0162         for (k = seeds.begin(); isok && k != seeds.end(); k++) {  //i->hits, k->seeds
0163
0164           for (dE = -2; dE <= 2; dE++)
0165             for (dP = -4; dP <= 4; dP += 2) {
0166               PWrap = k->id.iphi() + dP;
0167               if (PWrap < 0)
0168                 PWrap += 72;
0169               if (PWrap > 72)
0170                 PWrap -= 72;
0171
0172               if ((i->id.iphi() == PWrap) && (i->id.ieta() == k->id.ieta() + dE))
0173                 isok = false;
0174             }
0175         }
0176         if (isok) {
0177           doCluster = true;
0178         }
0179       }
0180       if (doCluster) {
0181         seeds.push_back(*i);
0182
0183         bool clusterOk = makeCluster(i->id(), hf, geom, clusShp, Sclus);
0184         if (clusterOk) {  // cluster is _not_ ok if seed is rejected due to other cuts
0185           clusterShapes.push_back(clusShp);
0186           SuperClusters.push_back(Sclus);
0187         }
0188       }
0189     }  //end protoseed loop
0190   }  //end if seeCount
0191 }
0192
0193 bool HFClusterAlgo::makeCluster(const HcalDetId& seedid,
0194                                 const HFRecHitCollection& hf,
0195                                 const CaloGeometry* geom,
0196                                 HFEMClusterShape& clusShp,
0197                                 SuperCluster& Sclus) {
0198   double w = 0;  //sum over all log E's
0199   double wgt = 0;
0200   double w_x = 0;
0201   double w_y = 0;
0202   double w_z = 0;
0203
0204   double l_3 = 0;  //sum for enenergy in 3x3 long fibers etc.
0205   double s_3 = 0;
0206   double l_5 = 0;
0207   double s_5 = 0;
0208   double l_1 = 0;
0209   double s_1 = 0;
0210   int de, dp, phiWrap;
0211   double l_1e = 0;
0212   const GlobalPoint& sp = geom->getPosition(seedid);
0213   std::vector<double> coreCanid;
0214   std::vector<double>::const_iterator ci;
0215   HFRecHitCollection::const_iterator is, il;
0216   std::vector<DetId> usedHits;
0217
0218   HFRecHitCollection::const_iterator si;
0219   HcalDetId sid(HcalForward, seedid.ieta(), seedid.iphi(), 1);
0220   si = hf.find(sid);
0221
0222   bool clusterOk = true;  // assume the best to start...
0223
0224   // lots happens here
0225   // edge type 1 has 40/41 in 3x3 and 5x5
0226   bool edge_type1 = seedid.ietaAbs() == 39 && (seedid.iphi() % 4) == 3;
0227
0228   double e_seed = si->energy() * m_correctionByEta[indexByEta(si->id())];
0229
0230   for (de = -2; de <= 2; de++)
0231     for (dp = -4; dp <= 4; dp += 2) {
0232       phiWrap = seedid.iphi() + dp;
0233       if (phiWrap < 0)
0234         phiWrap += 72;
0235       if (phiWrap > 72)
0236         phiWrap -= 72;
0237
0238       /* Handling of phi-width change problems */
0239       if (edge_type1 && de == seedid.zside()) {
0240         if (dp == -2) {  // we want it in the 3x3
0241           phiWrap -= 2;
0242           if (phiWrap < 0)
0243             phiWrap += 72;
0244         } else if (dp == -4) {
0245           continue;  // but not double counted in 5x5
0246         }
0247       }
0248
0249       HcalDetId idl(HcalForward, seedid.ieta() + de, phiWrap, 1);
0250       HcalDetId ids(HcalForward, seedid.ieta() + de, phiWrap, 2);
0251
0252       il = hf.find(idl);
0253       is = hf.find(ids);
0254
0255       double e_long = 1.0;
0256       double e_short = 0.0;
0257       if (il != hf.end())
0258         e_long = il->energy() * m_correctionByEta[indexByEta(il->id())];
0259       if (e_long <= m_minTowerEnergy)
0260         e_long = 0.0;
0261       if (is != hf.end())
0262         e_short = is->energy() * m_correctionByEta[indexByEta(is->id())];
0263       if (e_short <= m_minTowerEnergy)
0264         e_short = 0.0;
0265       double eRatio = (e_long - e_short) / std::max(1.0, (e_long + e_short));
0266
0267       // require S/L > a minimum amount for inclusion
0268       if ((abs(eRatio) > m_maximumSL) && (std::max(e_long, e_short) > m_maximumRenergy)) {
0269         if (dp == 0 && de == 0)
0270           clusterOk = false;  // somehow, the seed is hosed
0271         continue;
0272       }
0273
0274       if ((il != hf.end()) && (isPMTHit(*il))) {
0275         if (dp == 0 && de == 0)
0276           clusterOk = false;  // somehow, the seed is hosed
0277         continue;             //continue to next hit, do not include this one in cluster
0278       }
0279
0280       if (e_long > m_minTowerEnergy && il != hf.end()) {
0281         // record usage
0282         usedHits.push_back(idl.rawId());
0283         // always in the 5x5
0284         l_5 += e_long;
0285         // maybe in the 3x3
0286         if ((de > -2) && (de < 2) && (dp > -4) && (dp < 4)) {
0287           l_3 += e_long;
0288
0289           // sometimes in the 1x1
0290           if ((dp == 0) && (de == 0)) {
0291             l_1 = e_long;
0292           }
0293
0294           // maybe in the core?
0295           if ((de > -2) && (de < 2) && (dp > -4) && (dp < 4) && (e_long > (.5 * e_seed))) {
0296             coreCanid.push_back(e_long);
0297           }
0298
0299           // position calculation
0300           const GlobalPoint& p = geom->getPosition(idl);
0301
0302           double d_p = p.phi() - sp.phi();
0303           while (d_p < -M_PI)
0304             d_p += 2 * M_PI;
0305           while (d_p > M_PI)
0306             d_p -= 2 * M_PI;
0307
0308           wgt = log((e_long));
0309           if (wgt > 0) {
0310             w += wgt;
0311
0312             w_x += (p.x()) * wgt;  //(p.x()-sp.x())*wgt;
0313             w_y += (p.y()) * wgt;
0314             w_z += (p.z()) * wgt;
0315           }
0316         }
0317       } else {
0318         if (dp == 0 && de == 0)
0319           clusterOk = false;  // somehow, the seed is hosed
0320       }
0321
0322       if (e_short > m_minTowerEnergy && is != hf.end()) {
0323         // record usage
0324         usedHits.push_back(ids.rawId());
0325         // always in the 5x5
0326         s_5 += e_short;
0327         // maybe in the 3x3
0328         if ((de > -2) && (de < 2) && (dp > -4) && (dp < 4)) {
0329           s_3 += e_short;
0330         }
0331         // sometimes in the 1x1
0332         if ((dp == 0) && (de == 0)) {
0333           s_1 = e_short;
0334         }
0335       }
0336     }
0337
0338   if (!clusterOk)
0339     return false;
0340
0341   //Core sorting done here
0342   std::sort(coreCanid.begin(), coreCanid.end(), CompareHFCore());
0343   for (ci = coreCanid.begin(); ci != coreCanid.end(); ci++) {
0344     if (ci == coreCanid.begin()) {
0345       l_1e = *ci;
0346     } else if (*ci > .5 * l_1e) {
0347       l_1e += *ci;
0348     }
0349   }  //core sorting end
0350
0351   double z_ = w_z / w;
0352   double x_ = w_x / w;
0353   double y_ = w_y / w;
0354   math::XYZPoint xyzclus(x_, y_, z_);
0355   //calcualte position, final
0356   double eta = xyzclus.eta();  //w_e/w+sp.eta();
0357
0358   double phi = xyzclus.phi();  //(wp_e/w)+sp.phi();
0359
0360   while (phi < -M_PI)
0361     phi += 2 * M_PI;
0362   while (phi > M_PI)
0363     phi -= 2 * M_PI;
0364
0365   //calculate cell phi and cell eta
0366   int idx = fabs(seedid.ieta()) - 29;
0367   int ipx = seedid.iphi();
0368   double Cphi = (phi - m_seedmnPhi[ipx]) / (m_seedMXphi[ipx] - m_seedmnPhi[ipx]);
0369   double Ceta = (fabs(eta) - m_seedmnEta[idx]) / (m_seedMXeta[idx] - m_seedmnEta[idx]);
0370
0371   //return  HFEMClusterShape, SuperCluster
0372   HFEMClusterShape myClusShp(l_1, s_1, l_3, s_3, l_5, s_5, l_1e, Ceta, Cphi, seedid);
0373   clusShp = myClusShp;
0374
0375   SuperCluster MySclus(l_3, xyzclus);
0376   Sclus = MySclus;
0377
0378   return clusterOk;
0379 }
0380 bool HFClusterAlgo::isPMTHit(const HFRecHit& hfr) {
0381   bool pmthit = false;
0382
0383   if ((hfr.flagField(HcalCaloFlagLabels::HFLongShort)) && (m_usePMTFlag))
0384     pmthit = true;
0385   if (!(m_isMC && !m_forcePulseFlagMC))
0386     if ((hfr.flagField(HcalCaloFlagLabels::HFDigiTime)) && (m_usePulseFlag))
0387       pmthit = true;
0388
0389   return pmthit;
0390 }
0391 void HFClusterAlgo::resetForRun() {
0392   edm::LogInfo("HFClusterAlgo") << "Resetting for Run!";
0393   for (int ii = 0; ii < 13; ii++) {
0394     m_cutByEta.push_back(-1);
0395     m_seedmnEta.push_back(99);
0396     m_seedMXeta.push_back(-1);
0397   }
0398 }