HcalRecAlgos/src/HFStripFilter.cc

0001 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "DataFormats/METReco/interface/HcalPhase1FlagLabels.h"
0004
0005 #include "RecoLocalCalo/HcalRecAlgos/interface/HFStripFilter.h"
0006 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0007
0008 #include <cmath>
0009
0010 HFStripFilter::HFStripFilter(const double stripThreshold,
0011                              const double maxThreshold,
0012                              const double timeMax,
0013                              const double maxStripTime,
0014                              const double wedgeCut,
0015                              const int seedHitIetaMax,
0016                              const int gap,
0017                              const int lstrips,
0018                              const int verboseLevel)
0019     : stripThreshold_(stripThreshold),
0020       maxThreshold_(maxThreshold),
0021       timeMax_(timeMax),
0022       maxStripTime_(maxStripTime),
0023       wedgeCut_(wedgeCut),
0024       seedHitIetaMax_(seedHitIetaMax),
0025       gap_(gap),
0026       lstrips_(lstrips),
0027       verboseLevel_(verboseLevel) {
0028   // For the description of CMSSW message logging, see
0029   // https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideMessageLogger
0030   if (verboseLevel_ >= 20)
0031     edm::LogInfo("HFStripFilter") << "constructor called";
0032 }
0033
0034 HFStripFilter::~HFStripFilter() {
0035   if (verboseLevel_ >= 20)
0036     edm::LogInfo("HFStripFilter") << "destructor called";
0037 }
0038
0039 void HFStripFilter::runFilter(HFRecHitCollection& rec, const HcalChannelQuality* myqual) const {
0040   if (verboseLevel_ >= 20)
0041     edm::LogInfo("HFStripFilter") << "runFilter called";
0042
0043   std::vector<HFRecHit> d1strip;
0044   std::vector<HFRecHit> d2strip;
0045
0046   HFRecHit d1max;
0047   HFRecHit d2max;
0048
0049   d1max.setEnergy(-10);
0050   // find d1 and d2 seed hits with max energy and time < timeMax_
0051   for (auto& it : rec) {
0052     if (it.time() > timeMax_ || it.time() < 0 || it.energy() < stripThreshold_)
0053       continue;
0054     // find HF hit with maximum signal in depth = 1
0055     if (it.id().depth() == 1) {
0056       if (it.energy() > d1max.energy() && std::abs(it.id().ieta()) < seedHitIetaMax_) {
0057         d1max = it;
0058       }
0059     }
0060   }
0061
0062   bool signStripIeta = false;
0063   int stripIphiMax = 0;
0064   int stripIetaMax = 0;
0065   int stripDepthMax = 0;
0066
0067   if (d1max.energy() > 0) {
0068     d1strip.push_back(d1max);
0069     signStripIeta = std::signbit(d1max.id().ieta());
0070     stripIphiMax = d1max.id().iphi();
0071     stripIetaMax = d1max.id().ieta();
0072     stripDepthMax = d1max.id().depth();
0073   }
0074
0075   d2max.setEnergy(-10);
0076   for (auto& it : rec) {
0077     if (it.time() > timeMax_ || it.time() < 0 || it.energy() < stripThreshold_)
0078       continue;
0079     // find HFhit with maximum signal in depth = 2
0080     if (it.id().depth() == 2) {
0081       if (it.energy() > d2max.energy() && std::abs(it.id().ieta()) < seedHitIetaMax_) {
0082         if (d1max.energy() > 0) {
0083           bool signIeta = std::signbit(it.id().ieta());
0084           if (it.id().iphi() == stripIphiMax && signIeta == signStripIeta) {
0085             d2max = it;
0086           }
0087         } else {
0088           d2max = it;
0089         }
0090       }
0091     }
0092   }
0093
0094   if (d2max.energy() > 0)
0095     d2strip.push_back(d2max);
0096
0097   // check if at least one of possible seed hits has enough energy
0098   if (d1max.energy() < maxThreshold_ && d2max.energy() < maxThreshold_)
0099     return;
0100
0101   if (d1max.energy() <= 0 && d2max.energy() > 0) {
0102     signStripIeta = std::signbit(d2max.id().ieta());
0103     stripIphiMax = d2max.id().iphi();
0104     stripIetaMax = d2max.id().ieta();
0105     stripDepthMax = d2max.id().depth();
0106   }
0107
0108   std::stringstream ss;
0109   if (verboseLevel_ >= 30) {
0110     if (d1max.energy() > 0) {
0111       ss << "  MaxHit in Depth 1: ieta = " << d1max.id().ieta() << " iphi = " << stripIphiMax
0112          << " energy = " << d1max.energy() << " time = " << d1max.time() << std::endl;
0113     }
0114     if (d2max.energy() > 0) {
0115       ss << "  MaxHit in Depth 2: ieta = " << d2max.id().ieta() << " iphi = " << stripIphiMax
0116          << " energy = " << d2max.energy() << " time = " << d2max.time() << std::endl;
0117     }
0118     ss << "  stripThreshold_ = " << stripThreshold_ << std::endl;
0119   }
0120
0121   // prepare the strips: all hits along given ieta in one wedge (d1strip and d2strip)
0122   //---------------------------------------------------------------------------------
0123   for (auto& it : rec) {
0124     if (it.energy() < stripThreshold_)
0125       continue;
0126     bool signIeta = std::signbit(it.id().ieta());
0127
0128     if (verboseLevel_ >= 30) {
0129       ss << " HF hit: ieta = " << it.id().ieta() << "\t iphi = " << it.id().iphi() << "\t depth = " << it.id().depth()
0130          << "\t time = " << it.time() << "\t energy = " << it.energy() << "\t flags = " << it.flags() << std::endl;
0131     }
0132
0133     // collect hits with the same iphi but different ieta into strips
0134     if (it.id().iphi() == stripIphiMax && signIeta == signStripIeta && it.time() < maxStripTime_) {
0135       if (it.id().depth() == 1) {
0136         // check if hit = it is already in d1strip
0137         bool pass = false;
0138         if (d1strip.empty()) {
0139           if (std::abs(it.id().iphi() - stripIetaMax) <= gap_) {
0140             d1strip.push_back(it);
0141             pass = true;
0142           }
0143         } else {
0144           for (auto& it1 : d1strip) {
0145             if (it.id().ieta() == it1.id().ieta() && it.energy() == it1.energy()) {
0146               pass = true;
0147               break;
0148             }
0149           }
0150         }
0151         if (pass)
0152           continue;
0153         // add hit = it to d1strip if distance to the closest hit in d1strip <= gap_
0154         for (auto& it1 : d1strip) {
0155           // check distance along Ieta to the closest hit
0156           if (std::abs(it1.id().ieta() - it.id().ieta()) <= gap_) {
0157             d1strip.push_back(it);
0158             break;
0159           }
0160         }
0161       } else if (it.id().depth() == 2) {
0162         // check if hit = it is already in d2strip
0163         bool pass = false;
0164         if (d2strip.empty()) {
0165           if (std::abs(it.id().ieta() - stripIetaMax) <= gap_) {
0166             d2strip.push_back(it);
0167             pass = true;
0168           }
0169         } else {
0170           for (auto& it1 : d2strip) {
0171             if (it.id().ieta() == it1.id().ieta() && it.energy() == it1.energy()) {
0172               pass = true;
0173               break;
0174             }
0175           }
0176         }
0177         if (pass)
0178           continue;
0179
0180         // add hit = it to d2strip if distance to the closest hit in d1strip <= gap_
0181         for (auto& it1 : d2strip) {
0182           // check distance along Ieta to the closest hit
0183           if (std::abs(it1.id().ieta() - it.id().ieta()) <= gap_) {
0184             d2strip.push_back(it);
0185             break;
0186           }
0187         }
0188       }
0189     }
0190   }
0191
0192   if (verboseLevel_ >= 30) {
0193     ss << " Lstrip1 = " << (int)d1strip.size() << " (iphi = " << stripIphiMax << ")  Lstrip2 = " << (int)d2strip.size()
0194        << std::endl
0195        << " Strip1: ";
0196     for (auto& it1 : d1strip) {
0197       ss << it1.energy() << " (" << it1.id().ieta() << ") ";
0198     }
0199     ss << std::endl << " Strip2: ";
0200     for (auto& it1 : d2strip) {
0201       ss << it1.energy() << " (" << it1.id().ieta() << ") ";
0202     }
0203     ss << std::endl;
0204   }
0205
0206   // check if at least one of strips in depth1 or depth2 is not less than lstrips_
0207   if ((int)d1strip.size() < lstrips_ && (int)d2strip.size() < lstrips_) {
0208     if (verboseLevel_ >= 30) {
0209       LogDebug("HFSFilter") << ss.str();
0210     }
0211     return;
0212   }
0213
0214   // define range of strips in ieta
0215   int ietaMin1 = 1000;  // for d1strip
0216   int ietaMax1 = -1000;
0217   for (auto& it1 : d1strip) {
0218     if (it1.id().ieta() < ietaMin1)
0219       ietaMin1 = it1.id().ieta();
0220     if (it1.id().ieta() > ietaMax1)
0221       ietaMax1 = it1.id().ieta();
0222   }
0223   int ietaMin2 = 1000;  // for d2strip
0224   int ietaMax2 = -1000;
0225   for (auto& it1 : d2strip) {
0226     if (it1.id().ieta() < ietaMin2)
0227       ietaMin2 = it1.id().ieta();
0228     if (it1.id().ieta() > ietaMax2)
0229       ietaMax2 = it1.id().ieta();
0230   }
0231
0232   // define ietamin and ietamax - common area for d1strip and d2strip
0233   int ietaMin = ietaMin1;
0234   int ietaMax = ietaMax1;
0235
0236   if (ietaMin2 >= ietaMin1 && ietaMin2 <= ietaMax1) {
0237     if (ietaMax2 > ietaMax1)
0238       ietaMax = ietaMax2;
0239   } else if (ietaMin2 <= ietaMin1 && ietaMax2 >= ietaMin1) {
0240     if (ietaMin2 < ietaMin1)
0241       ietaMin = ietaMin2;
0242     if (ietaMax2 > ietaMax1)
0243       ietaMax = ietaMax2;
0244   }
0245
0246   // calculate the total energy in strips
0247   double eStrip = 0;
0248   for (auto& it1 : d1strip) {
0249     if (it1.id().ieta() < ietaMin || it1.id().ieta() > ietaMax)
0250       continue;
0251     eStrip += it1.energy();
0252   }
0253   for (auto& it1 : d2strip) {
0254     if (it1.id().ieta() < ietaMin || it1.id().ieta() > ietaMax)
0255       continue;
0256     eStrip += it1.energy();
0257   }
0258
0259   if (verboseLevel_ >= 30) {
0260     ss << " ietaMin1 = " << ietaMin1 << "  ietaMax1 = " << ietaMax1 << std::endl
0261        << " ietaMin2 = " << ietaMin2 << "  ietaMax2 = " << ietaMax2 << std::endl
0262        << " Common strip:  ietaMin = " << ietaMin << "  ietaMax = " << ietaMax << std::endl;
0263   }
0264
0265   // energies in the neighboring wedges
0266   double energyIphi1 = 0;
0267   double energyIphi2 = 0;
0268   int iphi1 = 0, iphi2 = 0;
0269
0270   // get information about the neighboring wedges from HcalTopology
0271   HcalDetId neighbour;
0272
0273   for (auto& it : rec) {
0274     if (it.energy() < stripThreshold_)
0275       continue;
0276     if (it.id().ieta() < ietaMin || it.id().ieta() > ietaMax)
0277       continue;
0278     HcalDetId id(HcalForward, it.id().ieta(), stripIphiMax, stripDepthMax);
0279     bool neigh = myqual->topo()->decIPhi(id, neighbour);
0280     iphi1 = (neigh) ? neighbour.iphi() : 0;
0281     neigh = myqual->topo()->incIPhi(id, neighbour);
0282     iphi2 = (neigh) ? neighbour.iphi() : 0;
0283     if (it.id().iphi() == iphi1)
0284       energyIphi1 += it.energy();  // Energy iphi1
0285     else if (it.id().iphi() == iphi2)
0286       energyIphi2 += it.energy();  // Energy iphi2
0287   }
0288
0289   double ratio1 = eStrip > 0 ? energyIphi1 / eStrip : 0;
0290   double ratio2 = eStrip > 0 ? energyIphi2 / eStrip : 0;
0291
0292   if (verboseLevel_ >= 30) {
0293     ss << "  iphi = " << stripIphiMax << "  iphi1 = " << iphi1 << "  iphi2 = " << iphi2 << std::endl
0294        << "  Estrip = " << eStrip << "  EnergyIphi1 = " << energyIphi1 << "  Ratio = " << ratio1 << std::endl
0295        << "                  "
0296        << "  EnergyIphi2 = " << energyIphi2 << "  Ratio = " << ratio2 << std::endl;
0297   }
0298
0299   // check if our wedge does not have substantial leak into adjacent wedges
0300   if (ratio1 < wedgeCut_ && ratio2 < wedgeCut_) {  // noise event with strips (d1 and/or d2)
0301
0302     if (verboseLevel_ >= 30) {
0303       ss << "  stripPass = false" << std::endl << "  mark hits in strips now " << std::endl;
0304     }
0305
0306     // Figure out which rechits need to be tagged (d1strip)
0307     for (auto& it1 : d1strip) {
0308       if (it1.id().ieta() < ietaMin || it1.id().ieta() > ietaMax)
0309         continue;
0310       HFRecHitCollection::iterator hit = rec.find(it1.id());
0311       if (hit != rec.end()) {
0312         // tag a rechit with the anomalous hit flag
0313         if (verboseLevel_ >= 30) {
0314           ss << " d1strip: marked hit = " << (*hit) << std::endl;
0315         }
0316         hit->setFlagField(1U, HcalPhase1FlagLabels::HFAnomalousHit);
0317       }
0318     }
0319
0320     // Figure out which rechits need to be tagged (dstrip)
0321     for (auto& it1 : d2strip) {
0322       if (it1.id().ieta() < ietaMin || it1.id().ieta() > ietaMax)
0323         continue;
0324       HFRecHitCollection::iterator hit = rec.find(it1.id());
0325       if (hit != rec.end()) {
0326         // tag a rechit with the anomalous hit flag
0327         if (verboseLevel_ >= 30) {
0328           ss << " d2strip: marked hit = " << (*hit) << std::endl;
0329         }
0330         hit->setFlagField(1U, HcalPhase1FlagLabels::HFAnomalousHit);
0331       }
0332     }
0333     if (verboseLevel_ >= 30) {
0334       LogDebug("HFSFilter") << ss.str();
0335     }
0336   } else {
0337     if (verboseLevel_ >= 30) {
0338       ss << "  stripPass = true" << std::endl;
0339       LogDebug("HFSFilter") << ss.str();
0340     }
0341   }
0342 }
0343
0344 std::unique_ptr<HFStripFilter> HFStripFilter::parseParameterSet(const edm::ParameterSet& ps) {
0345   return std::make_unique<HFStripFilter>(ps.getParameter<double>("stripThreshold"),
0346                                          ps.getParameter<double>("maxThreshold"),
0347                                          ps.getParameter<double>("timeMax"),
0348                                          ps.getParameter<double>("maxStripTime"),
0349                                          ps.getParameter<double>("wedgeCut"),
0350                                          ps.getParameter<int>("seedHitIetaMax"),
0351                                          ps.getParameter<int>("gap"),
0352                                          ps.getParameter<int>("lstrips"),
0353                                          ps.getUntrackedParameter<int>("verboseLevel"));
0354 }
0355
0356 edm::ParameterSetDescription HFStripFilter::fillDescription() {
0357   edm::ParameterSetDescription desc;
0358
0359   desc.add<double>("stripThreshold", 40.0)->setComment(" threshold to include hits into strips");
0360   desc.add<double>("maxThreshold", 100.0)->setComment(" threshold for seed hits in the strips (depth1 and depth2)");
0361   desc.add<double>("timeMax", 6.0)->setComment(" seed hits should have time < timeMax");
0362   desc.add<double>("maxStripTime", 10.0)->setComment(" maximum time for hits in the strips");
0363   desc.add<double>("wedgeCut", 0.05)->setComment(" the possible level of energy leak into adjacent wedges");
0364   desc.add<int>("seedHitIetaMax", 35)->setComment(" maximum possible Ieta value for seed hit");
0365   desc.add<int>("gap", 2)->setComment(" maximum distance between hits in the strip (along Ieta direction)");
0366   desc.add<int>("lstrips", 2)->setComment(" at least one of strips in depth1 or depth2 is not less than lstrips");
0367   desc.addUntracked<int>("verboseLevel", 0)
0368       ->setComment(" verboseLevel for debugging printouts, should be > 20 to get output");
0369
0370   return desc;
0371 }