HcalRecAlgos/src/HBHETimingShapedFlag.cc

0001 #include <iostream>
0002 #include "RecoLocalCalo/HcalRecAlgos/interface/HBHETimingShapedFlag.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 /*
0006 v1.0 by Jeff Temple
0007 29 August 2009
0008
0009 This takes the timing envelope algorithms developed by
0010 Phil Dudero and uses them to apply a RecHit Flag.
0011 */
0012
0013 HBHETimingShapedFlagSetter::HBHETimingShapedFlagSetter() {}
0014
0015 HBHETimingShapedFlagSetter::~HBHETimingShapedFlagSetter() {}
0016
0017 HBHETimingShapedFlagSetter::HBHETimingShapedFlagSetter(const std::vector<double>& v_userEnvelope) {
0018   makeTfilterEnvelope(v_userEnvelope);
0019
0020   ignorelowest_ = false;
0021   ignorehighest_ = false;
0022   win_offset_ = 0.;
0023   win_gain_ = 1.;
0024 }
0025
0026 HBHETimingShapedFlagSetter::HBHETimingShapedFlagSetter(const std::vector<double>& v_userEnvelope,
0027                                                        bool ignorelowest,
0028                                                        bool ignorehighest,
0029                                                        double win_offset,
0030                                                        double win_gain) {
0031   makeTfilterEnvelope(v_userEnvelope);
0032
0033   ignorelowest_ = ignorelowest;    // can ignore flagging hits below lowest energy threshold
0034   ignorehighest_ = ignorehighest;  // can ignore flagging hits above highest energy threshold
0035   win_offset_ = win_offset;        // timing offset
0036   win_gain_ = win_gain;            // time gain
0037 }
0038
0039 void HBHETimingShapedFlagSetter::makeTfilterEnvelope(const std::vector<double>& v_userEnvelope) {
0040   // Transform vector of doubles into a map of <energy,lowtime,hitime> triplets
0041   // Add extra protection in case vector of doubles is not a multiple of 3
0042   if (v_userEnvelope.size() % 3)
0043     throw cms::Exception("Invalid tfilterEnvelope definition") << "Must be one energy and two times per point";
0044
0045   for (unsigned int i = 0; i < v_userEnvelope.size(); i += 3) {
0046     int intGeV = (int)(v_userEnvelope[i] + 0.5);
0047     std::pair<double, double> pairOfTimes = std::make_pair(v_userEnvelope[i + 1], v_userEnvelope[i + 2]);
0048     if ((pairOfTimes.first > 0) || (pairOfTimes.second < 0))
0049       throw cms::Exception("Invalid tfilterEnvelope definition")
0050           << "Min and max time values must straddle t=0; use win_offset to shift";
0051
0052     tfilterEnvelope_[intGeV] = pairOfTimes;
0053   }
0054 }
0055
0056 void HBHETimingShapedFlagSetter::dumpInfo() {
0057   std::cout << "Timing Energy/Time parameters are:" << std::endl;
0058   TfilterEnvelope_t::const_iterator it;
0059   for (it = tfilterEnvelope_.begin(); it != tfilterEnvelope_.end(); ++it)
0060     std::cout << "\t" << it->first << " GeV\t" << it->second.first << " ns\t" << it->second.second << " ns"
0061               << std::endl;
0062
0063   std::cout << "ignorelowest  = " << ignorelowest_ << std::endl;
0064   std::cout << "ignorehighest = " << ignorehighest_ << std::endl;
0065   std::cout << "win_offset    = " << win_offset_ << std::endl;
0066   std::cout << "win_gain      = " << win_gain_ << std::endl;
0067 }
0068
0069 int HBHETimingShapedFlagSetter::timingStatus(const HBHERecHit& hbhe) {
0070   int status = 0;  // 3 bits reserved;status can range from 0-7
0071
0072   // tfilterEnvelope stores triples of energy and high/low time limits;
0073   // make sure we're checking over an even number of values
0074   // energies are guaranteed by std::map to appear in increasing order
0075
0076   // need at least two values to make comparison, and must
0077   // always have energy, time pair; otherwise, assume "in time" and don't set bits
0078   if (tfilterEnvelope_.empty())
0079     return 0;
0080
0081   double twinmin, twinmax;  // min, max 'good' time; values outside this range have flag set
0082   double rhtime = hbhe.time();
0083   double energy = hbhe.energy();
0084
0085   if (energy < (double)tfilterEnvelope_.begin()->first)  // less than lowest energy threshold
0086   {
0087     // Can skip timing test on cells below lowest threshold if so desired
0088     if (ignorelowest_)
0089       return 0;
0090     else {
0091       twinmin = tfilterEnvelope_.begin()->second.first;
0092       twinmax = tfilterEnvelope_.begin()->second.second;
0093     }
0094   } else {
0095     // Loop over energies in tfilterEnvelope
0096     TfilterEnvelope_t::const_iterator it;
0097     for (it = tfilterEnvelope_.begin(); it != tfilterEnvelope_.end(); ++it) {
0098       // Identify tfilterEnvelope index for this rechit energy
0099       if (energy < (double)it->first)
0100         break;
0101     }
0102
0103     if (it == tfilterEnvelope_.end()) {
0104       // Skip timing test on cells above highest threshold if so desired
0105       if (ignorehighest_)
0106         return 0;
0107       else {
0108         twinmin = tfilterEnvelope_.rbegin()->second.first;
0109         twinmax = tfilterEnvelope_.rbegin()->second.second;
0110       }
0111     } else {
0112       // Perform linear interpolation between energy boundaries
0113
0114       std::map<int, std::pair<double, double> >::const_iterator prev = it;
0115       prev--;
0116
0117       // twinmax interpolation
0118       double energy1 = prev->first;
0119       double lim1 = prev->second.second;
0120       double energy2 = it->first;
0121       double lim2 = it->second.second;
0122
0123       twinmax = lim1 + ((lim2 - lim1) * (energy - energy1) / (energy2 - energy1));
0124
0125       // twinmin interpolation
0126       lim1 = prev->second.first;
0127       lim2 = it->second.first;
0128
0129       twinmin = lim1 + ((lim2 - lim1) * (energy - energy1) / (energy2 - energy1));
0130     }
0131   }
0132
0133   // Apply offset, gain
0134   twinmin = win_offset_ + twinmin * win_gain_;
0135   twinmax = win_offset_ + twinmax * win_gain_;
0136
0137   // Set status high if time outside expected range
0138   if (rhtime <= twinmin || rhtime >= twinmax)
0139     status = 1;  // set status to 1
0140
0141   return status;
0142 }
0143
0144 void HBHETimingShapedFlagSetter::SetTimingShapedFlags(HBHERecHit& hbhe) {
0145   int status = timingStatus(hbhe);
0146
0147   // Though we're only using a single bit right now, 3 bits are reserved for these cuts
0148   hbhe.setFlagField(status, HcalCaloFlagLabels::HBHETimingShapedCutsBits, 3);
0149
0150   return;
0151 }