HcalObjects/src/HBHEDarkening.cc

0001 #include "CondFormats/HcalObjects/interface/HBHEDarkening.h"
0002
0003 #include "FWCore/Utilities/interface/Exception.h"
0004
0005 #include <vector>
0006 #include <string>
0007 #include <map>
0008 #include <algorithm>
0009 #include <fstream>
0010 #include <iostream>
0011 #include <string>
0012 #include <sstream>
0013 #include <cmath>
0014 #include <cassert>
0015
0016 HBHEDarkening::HBHEDarkening(int ieta_shift,
0017                              float drdA,
0018                              float drdB,
0019                              std::map<int, std::vector<std::vector<float>>> dosemaps,
0020                              std::vector<LumiYear> years)
0021     : ieta_shift_(ieta_shift), drdA_(drdA), drdB_(drdB), dosemaps_(std::move(dosemaps)), years_(std::move(years)) {
0022   //finish initializing years
0023   std::sort(years_.begin(), years_.end());
0024   //sum up int lumi
0025   float sumlumi = 0.0;
0026   for (auto& year : years_) {
0027     sumlumi += year.intlumi_;
0028     year.sumlumi_ = sumlumi;
0029   }
0030 }
0031
0032 std::vector<std::vector<float>> HBHEDarkening::readDoseMap(const std::string& fullpath) {
0033   std::ifstream infile(fullpath.c_str());
0034   if (!infile.is_open()) {
0035     throw cms::Exception("FileNotFound") << "Unable to open '" << fullpath << "'" << std::endl;
0036   }
0037   std::string line;
0038   std::vector<std::vector<float>> result;
0039   while (getline(infile, line)) {
0040     //space-separated
0041     std::stringstream linestream(line);
0042     std::vector<float> lineresult;
0043     float doseval;
0044     while (linestream >> doseval)
0045       lineresult.push_back(doseval);
0046     result.push_back(lineresult);
0047   }
0048   return result;
0049 }
0050
0051 float HBHEDarkening::dose(int ieta, int lay, int energy) const {
0052   //existence check
0053   const auto dosemapIt = dosemaps_.find(energy);
0054   if (dosemapIt == dosemaps_.end())
0055     return 0.0;
0056
0057   //bounds check
0058   const auto& dosemap = dosemapIt->second;
0059   if (ieta < 0 or ieta >= int(dosemap.size()))
0060     return 0.0;
0061
0062   //bounds check
0063   const auto& doserow = dosemap[ieta];
0064   if (lay < 0 or lay >= int(doserow.size()))
0065     return 0.0;
0066
0067   return doserow[lay];
0068 }
0069
0070 std::string HBHEDarkening::getYearForLumi(float intlumi) const {
0071   //compare based on sum lumi value
0072   auto lb = std::lower_bound(years_.begin(), years_.end(), intlumi, LumiYearComp());
0073   if (lb == years_.end() or lb->sumlumi_ < intlumi) {
0074     throw cms::Exception("ValueError") << "HBHEDarkening: insufficient LHC run information provided to simulate "
0075                                        << intlumi << "/fb - check the python config" << std::endl;
0076   }
0077   return lb->year_;
0078 }
0079
0080 float HBHEDarkening::degradationYear(const LumiYear& year, float intlumi, int ieta, int lay) const {
0081   float doseToUse = dose(ieta, lay, year.energy_);
0082   if (doseToUse == 0.0)
0083     return 1.0;
0084
0085   //apply dose rate dependence model to the provided year
0086   //get krad/hr from Mrad/fb-1 and fb-1/hr
0087   float decayConst = drdA_ * std::pow(1000 * doseToUse * year.lumirate_, drdB_);
0088
0089   //determine if this is a partial year
0090   float intlumiToUse = year.intlumi_;
0091   if (intlumi < year.sumlumi_)
0092     intlumiToUse = intlumi - (year.sumlumi_ - year.intlumi_);
0093
0094   //calculate degradation
0095   return std::exp(-(intlumiToUse * doseToUse) / decayConst);
0096 }
0097
0098 float HBHEDarkening::degradation(float intlumi, int ieta, int lay) const {
0099   ieta = abs(ieta);
0100   //shift ieta tower index to act as array index
0101   ieta -= ieta_shift_;
0102   //shift layer index by 1 to act as array index
0103   lay -= 1;
0104
0105   //accumulate degradation over years
0106   float response = 1.0;
0107   std::string yearForLumi = getYearForLumi(intlumi);
0108   assert(!yearForLumi.empty());
0109
0110   for (const auto& year : years_) {
0111     response *= degradationYear(year, intlumi, ieta, lay);
0112     if (year.year_ == yearForLumi)
0113       break;
0114   }
0115
0116   return response;
0117 }