HcalRecAlgos/src/HcalHFStatusBitFromDigis.cc

0001 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalHFStatusBitFromDigis.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 #include <algorithm>  // for "max"
0005 #include <cmath>
0006 #include <iostream>
0007
0008 HcalHFStatusBitFromDigis::HcalHFStatusBitFromDigis() {
0009   // use simple values in default constructor
0010   minthreshold_ = 40;  // minimum energy threshold (in GeV)
0011
0012   firstSample_ =
0013       1;  // these are the firstSample, samplesToAdd value of Igor's algorithm -- not necessarily the same as the reco first, toadd values (which are supplied individually for each hit)
0014   samplesToAdd_ = 3;
0015   expectedPeak_ = 2;
0016
0017   // Based on Igor V's algorithm:
0018   //TS4/(TS3+TS4+TS5+TS6) > 0.93 - exp(-0.38275-0.012667*E)
0019   coef_.push_back(0.93);
0020   coef_.push_back(-0.38275);
0021   coef_.push_back(-0.012667);
0022   // Minimum energy of 10 GeV required
0023   HFlongwindowEthresh_ = 40;
0024   HFlongwindowMinTime_.clear();
0025   HFlongwindowMinTime_.push_back(-10);
0026   HFlongwindowMaxTime_.clear();
0027   HFlongwindowMaxTime_.push_back(8);
0028   HFshortwindowEthresh_ = 40;
0029   HFshortwindowMinTime_.clear();
0030   HFshortwindowMinTime_.push_back(-10);
0031   HFshortwindowMaxTime_.clear();
0032   HFshortwindowMaxTime_.push_back(8);
0033 }
0034
0035 HcalHFStatusBitFromDigis::HcalHFStatusBitFromDigis(const edm::ParameterSet& HFDigiTimeParams,
0036                                                    const edm::ParameterSet& HFTimeInWindowParams) {
0037   // Specify parameters used in forming the HFDigiTime flag
0038   firstSample_ = HFDigiTimeParams.getParameter<int>("HFdigiflagFirstSample");
0039   samplesToAdd_ = HFDigiTimeParams.getParameter<int>("HFdigiflagSamplesToAdd");
0040   expectedPeak_ = HFDigiTimeParams.getParameter<int>("HFdigiflagExpectedPeak");
0041   minthreshold_ = HFDigiTimeParams.getParameter<double>("HFdigiflagMinEthreshold");
0042   coef_ = HFDigiTimeParams.getParameter<std::vector<double> >("HFdigiflagCoef");
0043
0044   // Specify parameters used in forming HFInTimeWindow flag
0045   HFlongwindowMinTime_ = HFTimeInWindowParams.getParameter<std::vector<double> >("hflongMinWindowTime");
0046   HFlongwindowMaxTime_ = HFTimeInWindowParams.getParameter<std::vector<double> >("hflongMaxWindowTime");
0047   HFlongwindowEthresh_ = HFTimeInWindowParams.getParameter<double>("hflongEthresh");
0048   HFshortwindowMinTime_ = HFTimeInWindowParams.getParameter<std::vector<double> >("hfshortMinWindowTime");
0049   HFshortwindowMaxTime_ = HFTimeInWindowParams.getParameter<std::vector<double> >("hfshortMaxWindowTime");
0050   HFshortwindowEthresh_ = HFTimeInWindowParams.getParameter<double>("hfshortEthresh");
0051 }
0052
0053 HcalHFStatusBitFromDigis::~HcalHFStatusBitFromDigis() {}
0054
0055 void HcalHFStatusBitFromDigis::resetParamsFromDB(
0056     int firstSample, int samplesToAdd, int expectedPeak, double minthreshold, const std::vector<double>& coef) {
0057   // Resets values based on values in database.
0058   firstSample_ = firstSample;
0059   samplesToAdd_ = samplesToAdd;
0060   expectedPeak_ = expectedPeak;
0061   minthreshold_ = minthreshold;
0062   coef_ = coef;
0063 }
0064
0065 void HcalHFStatusBitFromDigis::resetFlagTimeSamples(int firstSample, int samplesToAdd, int expectedPeak) {
0066   // This resets the time samples used in the HF flag.  These samples are not necessarily the same
0067   // as the flags used by the energy reconstruction
0068   firstSample_ = firstSample;
0069   samplesToAdd_ = samplesToAdd;
0070   expectedPeak_ = expectedPeak;
0071 }  // void HcalHFStatusBitFromDigis
0072
0073 void HcalHFStatusBitFromDigis::hfSetFlagFromDigi(HFRecHit& hf,
0074                                                  const HFDataFrame& digi,
0075                                                  const HcalCoder& coder,
0076                                                  const HcalCalibrations& calib) {
0077   // The following 3 values are computed by Igor's algorithm
0078   //only in the window [firstSample_, firstSample_ + samplesToAdd_),
0079   //which may not be the same as the default reco window.
0080
0081   double totalCharge = 0;
0082   double peakCharge = 0;
0083   double RecomputedEnergy = 0;
0084
0085   CaloSamples tool;
0086   coder.adc2fC(digi, tool);
0087
0088   // Compute quantities needed for HFDigiTime, Fraction2TS FLAGS
0089   for (int i = 0; i < digi.size(); ++i) {
0090     int capid = digi.sample(i).capid();
0091     double value = tool[i] - calib.pedestal(capid);
0092
0093     // Sum all charge within flagging window, find charge in expected peak time slice
0094     if (i >= firstSample_ && i < firstSample_ + samplesToAdd_) {
0095       totalCharge += value;
0096       RecomputedEnergy += value * calib.respcorrgain(capid);
0097       if (i == expectedPeak_)
0098         peakCharge = value;
0099     }
0100   }  // for (int i=0;i<digi.size();++i)
0101
0102   // FLAG:  HcalCaloLabels::HFDigiTime
0103   // Igor's algorithm:  compare charge in peak to total charge in window
0104   if (RecomputedEnergy >= minthreshold_)  // don't set noise flags for cells below a given threshold
0105   {
0106     // Calculate allowed minimum value of (TS4/TS3+4+5+6):
0107     double cutoff = 0;  // no arguments specified; no cutoff applied
0108     if (!coef_.empty())
0109       cutoff = coef_[0];
0110     // default cutoff takes the form:
0111     // cutoff = coef_[0] - exp(coef_[1]+coef_[2]*E+coef_[3]*E^2+...)
0112     double powRE = 1;
0113     double expo_arg = 0;
0114     for (unsigned int zz = 1; zz < coef_.size(); ++zz) {
0115       expo_arg += coef_[zz] * powRE;
0116       powRE *= RecomputedEnergy;
0117     }
0118     cutoff -= exp(expo_arg);
0119
0120     if (peakCharge / totalCharge < cutoff)
0121       hf.setFlagField(1, HcalCaloFlagLabels::HFDigiTime);
0122   }
0123
0124   // FLAG:  HcalCaloLabels:: HFInTimeWindow
0125   // Timing algorithm
0126   if (hf.id().depth() == 1) {
0127     if (hf.energy() >= HFlongwindowEthresh_) {
0128       float mult = 1. / hf.energy();
0129       float enPow = 1.;
0130       float mintime = 0;
0131       float maxtime = 0;
0132       for (unsigned int i = 0; i < HFlongwindowMinTime_.size(); ++i) {
0133         mintime += HFlongwindowMinTime_[i] * enPow;
0134         maxtime += HFlongwindowMaxTime_[i] * enPow;
0135         enPow *= mult;
0136       }
0137       if (hf.time() < mintime || hf.time() > maxtime)
0138         hf.setFlagField(1, HcalCaloFlagLabels::HFInTimeWindow);
0139     }
0140   } else if (hf.id().depth() == 2) {
0141     if (hf.energy() >= HFshortwindowEthresh_) {
0142       float mult = 1. / hf.energy();
0143       float enPow = 1.;
0144       float mintime = 0;
0145       float maxtime = 0;
0146       for (unsigned int i = 0; i < HFshortwindowMinTime_.size(); ++i) {
0147         mintime += HFshortwindowMinTime_[i] * enPow;
0148         maxtime += HFshortwindowMaxTime_[i] * enPow;
0149         enPow *= mult;
0150       }
0151       if (hf.time() < mintime || hf.time() > maxtime)
0152         hf.setFlagField(1, HcalCaloFlagLabels::HFInTimeWindow);
0153     }
0154   }
0155
0156   return;
0157 }