CMS 3D CMS Logo

HBHEDarkening.cc
Go to the documentation of this file.
2 
4 
5 #include <vector>
6 #include <string>
7 #include <map>
8 #include <algorithm>
9 #include <fstream>
10 #include <iostream>
11 #include <string>
12 #include <sstream>
13 #include <cmath>
14 #include <cassert>
15 
17  float drdA,
18  float drdB,
19  const std::map<int, std::vector<std::vector<float>>>& dosemaps,
20  const std::vector<LumiYear>& years)
21  : ieta_shift_(ieta_shift), drdA_(drdA), drdB_(drdB), dosemaps_(dosemaps), years_(years) {
22  //finish initializing years
23  std::sort(years_.begin(), years_.end());
24  //sum up int lumi
25  float sumlumi = 0.0;
26  for (auto& year : years_) {
27  sumlumi += year.intlumi_;
28  year.sumlumi_ = sumlumi;
29  }
30 }
31 
32 std::vector<std::vector<float>> HBHEDarkening::readDoseMap(const std::string& fullpath) {
33  std::ifstream infile(fullpath.c_str());
34  if (!infile.is_open()) {
35  throw cms::Exception("FileNotFound") << "Unable to open '" << fullpath << "'" << std::endl;
36  }
38  std::vector<std::vector<float>> result;
39  while (getline(infile, line)) {
40  //space-separated
41  std::stringstream linestream(line);
42  std::vector<float> lineresult;
43  float doseval;
44  while (linestream >> doseval)
45  lineresult.push_back(doseval);
46  result.push_back(lineresult);
47  }
48  return result;
49 }
50 
51 float HBHEDarkening::dose(int ieta, int lay, int energy) const {
52  //existence check
53  const auto dosemapIt = dosemaps_.find(energy);
54  if (dosemapIt == dosemaps_.end())
55  return 0.0;
56 
57  //bounds check
58  const auto& dosemap = dosemapIt->second;
59  if (ieta < 0 or ieta >= int(dosemap.size()))
60  return 0.0;
61 
62  //bounds check
63  const auto& doserow = dosemap[ieta];
64  if (lay < 0 or lay >= int(doserow.size()))
65  return 0.0;
66 
67  return doserow[lay];
68 }
69 
71  //compare based on sum lumi value
72  auto lb = std::lower_bound(years_.begin(), years_.end(), intlumi, LumiYearComp());
73  if (lb == years_.end() or lb->sumlumi_ < intlumi) {
74  throw cms::Exception("ValueError") << "HBHEDarkening: insufficient LHC run information provided to simulate "
75  << intlumi << "/fb - check the python config" << std::endl;
76  }
77  return lb->year_;
78 }
79 
80 float HBHEDarkening::degradationYear(const LumiYear& year, float intlumi, int ieta, int lay) const {
81  float doseToUse = dose(ieta, lay, year.energy_);
82  if (doseToUse == 0.0)
83  return 1.0;
84 
85  //apply dose rate dependence model to the provided year
86  //get krad/hr from Mrad/fb-1 and fb-1/hr
87  float decayConst = drdA_ * std::pow(1000 * doseToUse * year.lumirate_, drdB_);
88 
89  //determine if this is a partial year
90  float intlumiToUse = year.intlumi_;
91  if (intlumi < year.sumlumi_)
92  intlumiToUse = intlumi - (year.sumlumi_ - year.intlumi_);
93 
94  //calculate degradation
95  return std::exp(-(intlumiToUse * doseToUse) / decayConst);
96 }
97 
98 float HBHEDarkening::degradation(float intlumi, int ieta, int lay) const {
99  ieta = abs(ieta);
100  //shift ieta tower index to act as array index
101  ieta -= ieta_shift_;
102  //shift layer index by 1 to act as array index
103  lay -= 1;
104 
105  //accumulate degradation over years
106  float response = 1.0;
107  std::string yearForLumi = getYearForLumi(intlumi);
108  assert(!yearForLumi.empty());
109 
110  for (const auto& year : years_) {
111  response *= degradationYear(year, intlumi, ieta, lay);
112  if (year.year_ == yearForLumi)
113  break;
114  }
115 
116  return response;
117 }
HBHEDarkening(int ieta_shift, float drdA, float drdB, const std::map< int, std::vector< std::vector< float >>> &dosemaps, const std::vector< LumiYear > &years)
float degradation(float intlumi, int ieta, int lay) const
std::string getYearForLumi(float intlumi) const
constexpr int pow(int x)
Definition: conifer.h:24
assert(be >=bs)
float degradationYear(const LumiYear &year, float intlumi, int ieta, int lay) const
std::map< int, std::vector< std::vector< float > > > dosemaps_
Definition: HBHEDarkening.h:63
static std::vector< std::vector< float > > readDoseMap(const std::string &fullpath)
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e< void, edm::EventID const &, edm::Timestamp const & > We also list in braces which AR_WATCH_USING_METHOD_ is used for those or
Definition: Activities.doc:12
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< LumiYear > years_
Definition: HBHEDarkening.h:64
float dose(int ieta, int lay, int energy) const