HcalAlgos/src/HBHERecalibration.cc

0001 #include "CalibCalorimetry/HcalAlgos/interface/HBHERecalibration.h"
0002
0003 #include <vector>
0004 #include <cmath>
0005
0006 //reuse parsing function to read mean energy table

0007 HBHERecalibration::HBHERecalibration(float intlumi, float cutoff, std::string meanenergies)
0008     : intlumi_(intlumi),
0009       cutoff_(cutoff),
0010       ieta_shift_(0),
0011       max_depth_(0),
0012       meanenergies_(HBHEDarkening::readDoseMap(meanenergies)),
0013       darkening_(nullptr) {}
0014
0015 HBHERecalibration::~HBHERecalibration() {}
0016
0017 void HBHERecalibration::setup(const std::vector<std::vector<int>>& m_segmentation, const HBHEDarkening* darkening) {
0018   darkening_ = darkening;
0019   ieta_shift_ = darkening_->get_ieta_shift();
0020
0021   //infer eta bounds

0022   int min_ieta = ieta_shift_ - 1;
0023   int max_ieta = min_ieta + meanenergies_.size();
0024   dsegm_.reserve(max_ieta - min_ieta);
0025   for (int ieta = min_ieta; ieta < max_ieta; ++ieta) {
0026     dsegm_.push_back(m_segmentation[ieta]);
0027     //find maximum

0028     for (unsigned lay = 0; lay < dsegm_.back().size(); ++lay) {
0029       if (lay >= meanenergies_[0].size())
0030         break;
0031       int depth = dsegm_.back()[lay];
0032       if (depth > max_depth_)
0033         max_depth_ = depth;
0034     }
0035   }
0036
0037   initialize();
0038 }
0039
0040 float HBHERecalibration::getCorr(int ieta, int depth) const {
0041   ieta = abs(ieta);
0042   //shift ieta tower index to act as array index

0043   ieta -= ieta_shift_;
0044
0045   //shift depth index to act as array index (depth = 0 - not used!)

0046   depth -= 1;
0047
0048   //bounds check

0049   if (ieta < 0 or ieta >= int(corr_.size()))
0050     return 1.0;
0051   if (depth < 0 or depth >= int(corr_[ieta].size()))
0052     return 1.0;
0053
0054   if (cutoff_ > 1 and corr_[ieta][depth] > cutoff_)
0055     return cutoff_;
0056   else
0057     return corr_[ieta][depth];
0058 }
0059
0060 void HBHERecalibration::initialize() {
0061   std::vector<std::vector<float>> vtmp(dsegm_.size(), std::vector<float>(max_depth_, 0.0));
0062   auto dval =
0063       vtmp;  //conversion of meanenergies into depths-averaged values - denominator (including degradation for intlumi)

0064   auto nval = vtmp;  // conversion of meanenergies into depths-averaged values - numerator (no degradation)

0065   corr_ = vtmp;
0066
0067   //converting energy values from layers into depths

0068   for (unsigned int ieta = 0; ieta < dsegm_.size(); ++ieta) {
0069     //fill sum(means(layer,0)) and sum(means(layer,lumi)) for each depth

0070     for (unsigned int ilay = 0; ilay < std::min(meanenergies_[ieta].size(), dsegm_[ieta].size()); ++ilay) {
0071       int depth = dsegm_[ieta][ilay] - 1;  // depth = 0 - not used!

0072       nval[ieta][depth] += meanenergies_[ieta][ilay];
0073       dval[ieta][depth] +=
0074           meanenergies_[ieta][ilay] *
0075           darkening_->degradation(intlumi_, ieta + ieta_shift_, ilay + 1);  //be careful of eta and layer numbering

0076     }
0077
0078     //compute factors, w/ safety checks

0079     for (int depth = 0; depth < max_depth_; ++depth) {
0080       if (dval[ieta][depth] > 0)
0081         corr_[ieta][depth] = nval[ieta][depth] / dval[ieta][depth];
0082       else
0083         corr_[ieta][depth] = 1.0;
0084
0085       if (corr_[ieta][depth] < 1.0)
0086         corr_[ieta][depth] = 1.0;
0087     }
0088   }
0089 }