File indexing completed on 2023-03-17 10:47:00
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 const std::map<int, std::vector<std::vector<float>>>& dosemaps,
0020 const std::vector<LumiYear>& years)
0021 : ieta_shift_(ieta_shift), drdA_(drdA), drdB_(drdB), dosemaps_(dosemaps), years_(years) {
0022
0023 std::sort(years_.begin(), years_.end());
0024
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
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
0053 const auto dosemapIt = dosemaps_.find(energy);
0054 if (dosemapIt == dosemaps_.end())
0055 return 0.0;
0056
0057
0058 const auto& dosemap = dosemapIt->second;
0059 if (ieta < 0 or ieta >= int(dosemap.size()))
0060 return 0.0;
0061
0062
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
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
0086
0087 float decayConst = drdA_ * std::pow(1000 * doseToUse * year.lumirate_, drdB_);
0088
0089
0090 float intlumiToUse = year.intlumi_;
0091 if (intlumi < year.sumlumi_)
0092 intlumiToUse = intlumi - (year.sumlumi_ - year.intlumi_);
0093
0094
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
0101 ieta -= ieta_shift_;
0102
0103 lay -= 1;
0104
0105
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 }