Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:34:06

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_e = 0;  //sum over ieat*energy
0201   double w_x = 0;
0202   double w_y = 0;
0203   double w_z = 0;
0204   double wp_e = 0;  //sum over iphi*energy
0205   double e_e = 0;   //nonwieghted eta sum
0206   double e_ep = 0;  //nonweighted phi sum
0207 
0208   double l_3 = 0;  //sum for enenergy in 3x3 long fibers etc.
0209   double s_3 = 0;
0210   double l_5 = 0;
0211   double s_5 = 0;
0212   double l_1 = 0;
0213   double s_1 = 0;
0214   int de, dp, phiWrap;
0215   double l_1e = 0;
0216   const GlobalPoint& sp = geom->getPosition(seedid);
0217   std::vector<double> coreCanid;
0218   std::vector<double>::const_iterator ci;
0219   HFRecHitCollection::const_iterator is, il;
0220   std::vector<DetId> usedHits;
0221 
0222   HFRecHitCollection::const_iterator si;
0223   HcalDetId sid(HcalForward, seedid.ieta(), seedid.iphi(), 1);
0224   si = hf.find(sid);
0225 
0226   bool clusterOk = true;  // assume the best to start...
0227 
0228   // lots happens here
0229   // edge type 1 has 40/41 in 3x3 and 5x5
0230   bool edge_type1 = seedid.ietaAbs() == 39 && (seedid.iphi() % 4) == 3;
0231 
0232   double e_seed = si->energy() * m_correctionByEta[indexByEta(si->id())];
0233 
0234   for (de = -2; de <= 2; de++)
0235     for (dp = -4; dp <= 4; dp += 2) {
0236       phiWrap = seedid.iphi() + dp;
0237       if (phiWrap < 0)
0238         phiWrap += 72;
0239       if (phiWrap > 72)
0240         phiWrap -= 72;
0241 
0242       /* Handling of phi-width change problems */
0243       if (edge_type1 && de == seedid.zside()) {
0244         if (dp == -2) {  // we want it in the 3x3
0245           phiWrap -= 2;
0246           if (phiWrap < 0)
0247             phiWrap += 72;
0248         } else if (dp == -4) {
0249           continue;  // but not double counted in 5x5
0250         }
0251       }
0252 
0253       HcalDetId idl(HcalForward, seedid.ieta() + de, phiWrap, 1);
0254       HcalDetId ids(HcalForward, seedid.ieta() + de, phiWrap, 2);
0255 
0256       il = hf.find(idl);
0257       is = hf.find(ids);
0258 
0259       double e_long = 1.0;
0260       double e_short = 0.0;
0261       if (il != hf.end())
0262         e_long = il->energy() * m_correctionByEta[indexByEta(il->id())];
0263       if (e_long <= m_minTowerEnergy)
0264         e_long = 0.0;
0265       if (is != hf.end())
0266         e_short = is->energy() * m_correctionByEta[indexByEta(is->id())];
0267       if (e_short <= m_minTowerEnergy)
0268         e_short = 0.0;
0269       double eRatio = (e_long - e_short) / std::max(1.0, (e_long + e_short));
0270 
0271       // require S/L > a minimum amount for inclusion
0272       if ((abs(eRatio) > m_maximumSL) && (std::max(e_long, e_short) > m_maximumRenergy)) {
0273         if (dp == 0 && de == 0)
0274           clusterOk = false;  // somehow, the seed is hosed
0275         continue;
0276       }
0277 
0278       if ((il != hf.end()) && (isPMTHit(*il))) {
0279         if (dp == 0 && de == 0)
0280           clusterOk = false;  // somehow, the seed is hosed
0281         continue;             //continue to next hit, do not include this one in cluster
0282       }
0283 
0284       if (e_long > m_minTowerEnergy && il != hf.end()) {
0285         // record usage
0286         usedHits.push_back(idl.rawId());
0287         // always in the 5x5
0288         l_5 += e_long;
0289         // maybe in the 3x3
0290         if ((de > -2) && (de < 2) && (dp > -4) && (dp < 4)) {
0291           l_3 += e_long;
0292 
0293           // sometimes in the 1x1
0294           if ((dp == 0) && (de == 0)) {
0295             l_1 = e_long;
0296           }
0297 
0298           // maybe in the core?
0299           if ((de > -2) && (de < 2) && (dp > -4) && (dp < 4) && (e_long > (.5 * e_seed))) {
0300             coreCanid.push_back(e_long);
0301           }
0302 
0303           // position calculation
0304           const GlobalPoint& p = geom->getPosition(idl);
0305 
0306           double d_p = p.phi() - sp.phi();
0307           while (d_p < -M_PI)
0308             d_p += 2 * M_PI;
0309           while (d_p > M_PI)
0310             d_p -= 2 * M_PI;
0311           double d_e = p.eta() - sp.eta();
0312 
0313           wgt = log((e_long));
0314           if (wgt > 0) {
0315             w += wgt;
0316             w_e += (d_e)*wgt;
0317             wp_e += (d_p)*wgt;
0318             e_e += d_e;
0319             e_ep += d_p;
0320 
0321             w_x += (p.x()) * wgt;  //(p.x()-sp.x())*wgt;
0322             w_y += (p.y()) * wgt;
0323             w_z += (p.z()) * wgt;
0324           }
0325         }
0326       } else {
0327         if (dp == 0 && de == 0)
0328           clusterOk = false;  // somehow, the seed is hosed
0329       }
0330 
0331       if (e_short > m_minTowerEnergy && is != hf.end()) {
0332         // record usage
0333         usedHits.push_back(ids.rawId());
0334         // always in the 5x5
0335         s_5 += e_short;
0336         // maybe in the 3x3
0337         if ((de > -2) && (de < 2) && (dp > -4) && (dp < 4)) {
0338           s_3 += e_short;
0339         }
0340         // sometimes in the 1x1
0341         if ((dp == 0) && (de == 0)) {
0342           s_1 = e_short;
0343         }
0344       }
0345     }
0346 
0347   if (!clusterOk)
0348     return false;
0349 
0350   //Core sorting done here
0351   std::sort(coreCanid.begin(), coreCanid.end(), CompareHFCore());
0352   for (ci = coreCanid.begin(); ci != coreCanid.end(); ci++) {
0353     if (ci == coreCanid.begin()) {
0354       l_1e = *ci;
0355     } else if (*ci > .5 * l_1e) {
0356       l_1e += *ci;
0357     }
0358   }  //core sorting end
0359 
0360   double z_ = w_z / w;
0361   double x_ = w_x / w;
0362   double y_ = w_y / w;
0363   math::XYZPoint xyzclus(x_, y_, z_);
0364   //calcualte position, final
0365   double eta = xyzclus.eta();  //w_e/w+sp.eta();
0366 
0367   double phi = xyzclus.phi();  //(wp_e/w)+sp.phi();
0368 
0369   while (phi < -M_PI)
0370     phi += 2 * M_PI;
0371   while (phi > M_PI)
0372     phi -= 2 * M_PI;
0373 
0374   //calculate cell phi and cell eta
0375   int idx = fabs(seedid.ieta()) - 29;
0376   int ipx = seedid.iphi();
0377   double Cphi = (phi - m_seedmnPhi[ipx]) / (m_seedMXphi[ipx] - m_seedmnPhi[ipx]);
0378   double Ceta = (fabs(eta) - m_seedmnEta[idx]) / (m_seedMXeta[idx] - m_seedmnEta[idx]);
0379 
0380   //return  HFEMClusterShape, SuperCluster
0381   HFEMClusterShape myClusShp(l_1, s_1, l_3, s_3, l_5, s_5, l_1e, Ceta, Cphi, seedid);
0382   clusShp = myClusShp;
0383 
0384   SuperCluster MySclus(l_3, xyzclus);
0385   Sclus = MySclus;
0386 
0387   return clusterOk;
0388 }
0389 bool HFClusterAlgo::isPMTHit(const HFRecHit& hfr) {
0390   bool pmthit = false;
0391 
0392   if ((hfr.flagField(HcalCaloFlagLabels::HFLongShort)) && (m_usePMTFlag))
0393     pmthit = true;
0394   if (!(m_isMC && !m_forcePulseFlagMC))
0395     if ((hfr.flagField(HcalCaloFlagLabels::HFDigiTime)) && (m_usePulseFlag))
0396       pmthit = true;
0397 
0398   return pmthit;
0399 }
0400 void HFClusterAlgo::resetForRun() {
0401   edm::LogInfo("HFClusterAlgo") << "Resetting for Run!";
0402   for (int ii = 0; ii < 13; ii++) {
0403     m_cutByEta.push_back(-1);
0404     m_seedmnEta.push_back(99);
0405     m_seedMXeta.push_back(-1);
0406   }
0407 }