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 void HcalHFStatusBitFromDigis::fillHFDigiTimeParamsDesc(edm::ParameterSetDescription& desc) {
0054   desc.add<int>("HFdigiflagFirstSample", 1);
0055   desc.add<int>("HFdigiflagSamplesToAdd", 1);
0056   desc.add<int>("HFdigiflagExpectedPeak", 2);
0057   desc.add<double>("HFdigiflagMinEthreshold", {40});
0058   desc.add<std::vector<double> >("HFdigiflagCoef", {0.93, -0.38275, -0.012667});
0059 }
0060
0061 void HcalHFStatusBitFromDigis::fillHFTimeInWindowParamsDesc(edm::ParameterSetDescription& desc) {
0062   desc.add<std::vector<double> >("hflongMinWindowTime", {-10.});
0063   desc.add<std::vector<double> >("hflongMaxWindowTime", {10.});
0064   desc.add<double>("hflongEthresh", 40.);
0065   desc.add<std::vector<double> >("hfshortMinWindowTime", {-12.});
0066   desc.add<std::vector<double> >("hfshortMaxWindowTime", {10.});
0067   desc.add<double>("hfshortEthresh", 40.);
0068 }
0069
0070 HcalHFStatusBitFromDigis::~HcalHFStatusBitFromDigis() {}
0071
0072 void HcalHFStatusBitFromDigis::resetParamsFromDB(
0073     int firstSample, int samplesToAdd, int expectedPeak, double minthreshold, const std::vector<double>& coef) {
0074   // Resets values based on values in database.
0075   firstSample_ = firstSample;
0076   samplesToAdd_ = samplesToAdd;
0077   expectedPeak_ = expectedPeak;
0078   minthreshold_ = minthreshold;
0079   coef_ = coef;
0080 }
0081
0082 void HcalHFStatusBitFromDigis::resetFlagTimeSamples(int firstSample, int samplesToAdd, int expectedPeak) {
0083   // This resets the time samples used in the HF flag.  These samples are not necessarily the same
0084   // as the flags used by the energy reconstruction
0085   firstSample_ = firstSample;
0086   samplesToAdd_ = samplesToAdd;
0087   expectedPeak_ = expectedPeak;
0088 }  // void HcalHFStatusBitFromDigis
0089
0090 void HcalHFStatusBitFromDigis::hfSetFlagFromDigi(HFRecHit& hf,
0091                                                  const HFDataFrame& digi,
0092                                                  const HcalCoder& coder,
0093                                                  const HcalCalibrations& calib) {
0094   // The following 3 values are computed by Igor's algorithm
0095   //only in the window [firstSample_, firstSample_ + samplesToAdd_),
0096   //which may not be the same as the default reco window.
0097
0098   double totalCharge = 0;
0099   double peakCharge = 0;
0100   double RecomputedEnergy = 0;
0101
0102   CaloSamples tool;
0103   coder.adc2fC(digi, tool);
0104
0105   // Compute quantities needed for HFDigiTime, Fraction2TS FLAGS
0106   for (int i = 0; i < digi.size(); ++i) {
0107     int capid = digi.sample(i).capid();
0108     double value = tool[i] - calib.pedestal(capid);
0109
0110     // Sum all charge within flagging window, find charge in expected peak time slice
0111     if (i >= firstSample_ && i < firstSample_ + samplesToAdd_) {
0112       totalCharge += value;
0113       RecomputedEnergy += value * calib.respcorrgain(capid);
0114       if (i == expectedPeak_)
0115         peakCharge = value;
0116     }
0117   }  // for (int i=0;i<digi.size();++i)
0118
0119   // FLAG:  HcalCaloLabels::HFDigiTime
0120   // Igor's algorithm:  compare charge in peak to total charge in window
0121   if (RecomputedEnergy >= minthreshold_)  // don't set noise flags for cells below a given threshold
0122   {
0123     // Calculate allowed minimum value of (TS4/TS3+4+5+6):
0124     double cutoff = 0;  // no arguments specified; no cutoff applied
0125     if (!coef_.empty())
0126       cutoff = coef_[0];
0127     // default cutoff takes the form:
0128     // cutoff = coef_[0] - exp(coef_[1]+coef_[2]*E+coef_[3]*E^2+...)
0129     double powRE = 1;
0130     double expo_arg = 0;
0131     for (unsigned int zz = 1; zz < coef_.size(); ++zz) {
0132       expo_arg += coef_[zz] * powRE;
0133       powRE *= RecomputedEnergy;
0134     }
0135     cutoff -= exp(expo_arg);
0136
0137     if (peakCharge / totalCharge < cutoff)
0138       hf.setFlagField(1, HcalCaloFlagLabels::HFDigiTime);
0139   }
0140
0141   // FLAG:  HcalCaloLabels:: HFInTimeWindow
0142   // Timing algorithm
0143   if (hf.id().depth() == 1) {
0144     if (hf.energy() >= HFlongwindowEthresh_) {
0145       float mult = 1. / hf.energy();
0146       float enPow = 1.;
0147       float mintime = 0;
0148       float maxtime = 0;
0149       for (unsigned int i = 0; i < HFlongwindowMinTime_.size(); ++i) {
0150         mintime += HFlongwindowMinTime_[i] * enPow;
0151         maxtime += HFlongwindowMaxTime_[i] * enPow;
0152         enPow *= mult;
0153       }
0154       if (hf.time() < mintime || hf.time() > maxtime)
0155         hf.setFlagField(1, HcalCaloFlagLabels::HFInTimeWindow);
0156     }
0157   } else if (hf.id().depth() == 2) {
0158     if (hf.energy() >= HFshortwindowEthresh_) {
0159       float mult = 1. / hf.energy();
0160       float enPow = 1.;
0161       float mintime = 0;
0162       float maxtime = 0;
0163       for (unsigned int i = 0; i < HFshortwindowMinTime_.size(); ++i) {
0164         mintime += HFshortwindowMinTime_[i] * enPow;
0165         maxtime += HFshortwindowMaxTime_[i] * enPow;
0166         enPow *= mult;
0167       }
0168       if (hf.time() < mintime || hf.time() > maxtime)
0169         hf.setFlagField(1, HcalCaloFlagLabels::HFInTimeWindow);
0170     }
0171   }
0172
0173   return;
0174 }