HcalRecAlgos/src/HBHEStatusBitSetter.cc

0001 #include <cassert>
0002 #include "RecoLocalCalo/HcalRecAlgos/interface/HBHEStatusBitSetter.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 HBHEStatusBitSetter::HBHEStatusBitSetter() : frontEndMap_(nullptr) {
0006   nominalPedestal_ = 3.0;
0007   hitEnergyMinimum_ = 2.0;
0008   hitMultiplicityThreshold_ = 17;
0009 }
0010
0011 HBHEStatusBitSetter::HBHEStatusBitSetter(double nominalPedestal,
0012                                          double hitEnergyMinimum,
0013                                          int hitMultiplicityThreshold,
0014                                          const std::vector<edm::ParameterSet>& pulseShapeParameterSets)
0015     : hitEnergyMinimum_(hitEnergyMinimum),
0016       hitMultiplicityThreshold_(hitMultiplicityThreshold),
0017       nominalPedestal_(nominalPedestal),
0018       frontEndMap_(nullptr) {
0019   const unsigned sz = pulseShapeParameterSets.size();
0020   pulseShapeParameters_.reserve(sz);
0021   for (unsigned iPSet = 0; iPSet < sz; iPSet++) {
0022     const edm::ParameterSet& pset = pulseShapeParameterSets[iPSet];
0023     const std::vector<double>& params = pset.getParameter<std::vector<double> >("pulseShapeParameters");
0024     pulseShapeParameters_.push_back(params);
0025   }
0026 }
0027
0028 HBHEStatusBitSetter::~HBHEStatusBitSetter() {}
0029
0030 void HBHEStatusBitSetter::SetFrontEndMap(const HcalFrontEndMap* m) {
0031   frontEndMap_ = m;
0032   hpdMultiplicity_.clear();
0033   if (frontEndMap_) {
0034     const int sz = frontEndMap_->maxRMIndex();
0035     hpdMultiplicity_.reserve(sz);
0036     for (int iRm = 0; iRm < sz; iRm++) {
0037       hpdMultiplicity_.push_back(0);
0038     }
0039   }
0040 }
0041
0042 void HBHEStatusBitSetter::Clear() {
0043   const unsigned sz = hpdMultiplicity_.size();
0044   for (unsigned i = 0; i < sz; i++)
0045     hpdMultiplicity_[i] = 0;
0046 }
0047
0048 void HBHEStatusBitSetter::rememberHit(const HBHERecHit& hbhe) {
0049   if (frontEndMap_ == nullptr) {
0050     edm::LogError("HBHEStatusBitSetter") << "No HcalFrontEndMap in rememberHit";
0051     return;
0052   }
0053   //increment hit multiplicity
0054   if (hbhe.energy() > hitEnergyMinimum_) {
0055     const int index = frontEndMap_->lookupRMIndex(hbhe.detid());
0056     hpdMultiplicity_.at(index)++;
0057   }
0058 }
0059
0060 void HBHEStatusBitSetter::SetFlagsFromDigi(HBHERecHit& hbhe,
0061                                            const HBHEDataFrame& digi,
0062                                            const HcalCoder& coder,
0063                                            const HcalCalibrations& calib) {
0064   if (frontEndMap_ == nullptr) {
0065     edm::LogError("HBHEStatusBitSetter") << "No HcalFrontEndMap in SetFlagsFromDigi";
0066     return;
0067   }
0068   rememberHit(hbhe);
0069
0070   //set pulse shape bits
0071   // Shuichi's algorithm uses the "correct" charge & pedestals, while Ted's uses "nominal" values.
0072   // Perhaps we should correct Ted's algorithm in the future, though that will mean re-tuning thresholds for his cuts. -- Jeff, 28 May 2010
0073   //double shuichi_charge_total=0.0;
0074   double nominal_charge_total = 0.0;
0075   double charge_max3 = -100.0;
0076   double charge_late3 = -100.0;
0077   unsigned int slice_max3 = 0;
0078   unsigned int size = digi.size();
0079
0080   CaloSamples tool;
0081   coder.adc2fC(digi, tool);
0082
0083   //  int capid=-1;
0084   for (unsigned int iSlice = 0; iSlice < size; iSlice++) {
0085     //      capid  = digi.sample(iSlice).capid();
0086     //shuichi_charge_total+=tool[iSlice]-calib.pedestal(capid);
0087     nominal_charge_total += digi[iSlice].nominal_fC() - nominalPedestal_;
0088
0089     if (iSlice < 2)
0090       continue;
0091     // digi[i].nominal_fC() could be replaced by tool[iSlice], I think...  -- Jeff, 11 April 2011
0092     double qsum3 = digi[iSlice].nominal_fC() + digi[iSlice - 1].nominal_fC() + digi[iSlice - 2].nominal_fC() -
0093                    3 * nominalPedestal_;
0094     if (qsum3 > charge_max3) {
0095       charge_max3 = qsum3;
0096       slice_max3 = iSlice;
0097     }
0098   }
0099
0100   if ((4 + slice_max3) > size)
0101     return;
0102   charge_late3 = digi[slice_max3 + 1].nominal_fC() + digi[slice_max3 + 2].nominal_fC() +
0103                  digi[slice_max3 + 3].nominal_fC() - 3 * nominalPedestal_;
0104
0105   for (unsigned int iCut = 0; iCut < pulseShapeParameters_.size(); iCut++) {
0106     if (pulseShapeParameters_[iCut].size() != 6)
0107       continue;
0108     if (nominal_charge_total < pulseShapeParameters_[iCut].at(0) ||
0109         nominal_charge_total >= pulseShapeParameters_[iCut].at(1))
0110       continue;
0111     if (charge_late3 < (pulseShapeParameters_[iCut].at(2) + nominal_charge_total * pulseShapeParameters_[iCut].at(3)))
0112       continue;
0113     if (charge_late3 >= (pulseShapeParameters_[iCut].at(4) + nominal_charge_total * pulseShapeParameters_[iCut].at(5)))
0114       continue;
0115     hbhe.setFlagField(1, HcalCaloFlagLabels::HBHEPulseShape);
0116     return;
0117   }
0118 }
0119
0120 void HBHEStatusBitSetter::SetFlagsFromRecHits(HBHERecHitCollection& rec) {
0121   if (frontEndMap_ == nullptr) {
0122     edm::LogError("HBHEStatusBitSetter") << "No HcalFrontEndMap in SetFlagsFromRecHits";
0123     return;
0124   }
0125
0126   for (HBHERecHitCollection::iterator iHBHE = rec.begin(); iHBHE != rec.end(); ++iHBHE) {
0127     const int index = frontEndMap_->lookupRMIndex(iHBHE->detid());
0128     if (hpdMultiplicity_.at(index) < hitMultiplicityThreshold_)
0129       continue;
0130     iHBHE->setFlagField(1, HcalCaloFlagLabels::HBHEHpdHitMultiplicity);
0131   }
0132 }