File indexing completed on 2023-03-17 11:23:55
0001 #include "SimCalorimetry/HGCalSimAlgos/interface/HGCalRadiationMap.h"
0002 #include "FWCore/ParameterSet/interface/FileInPath.h"
0003 #include <fstream>
0004
0005
0006 HGCalRadiationMap::HGCalRadiationMap() : fluenceSFlog10_(0.) {}
0007
0008
0009 void HGCalRadiationMap::setDoseMap(const std::string& fullpath, const unsigned int algo) {
0010 doseMap_ = readDosePars(fullpath);
0011 algo_ = algo;
0012 }
0013
0014
0015 void HGCalRadiationMap::setGeometry(const CaloSubdetectorGeometry* geom) {
0016 hgcalGeom_ = static_cast<const HGCalGeometry*>(geom);
0017 hgcalTopology_ = &(hgcalGeom_->topology());
0018 hgcalDDD_ = &(hgcalTopology_->dddConstants());
0019 }
0020
0021
0022 double HGCalRadiationMap::computeRadius(const HGCScintillatorDetId& cellId) {
0023 GlobalPoint global = geom()->getPosition(cellId);
0024 return std::sqrt(std::pow(global.x(), 2) + std::pow(global.y(), 2));
0025 }
0026
0027
0028 double HGCalRadiationMap::getDoseValue(const int subdet, const int layer, const double radius, bool logVal) {
0029 std::pair<int, int> key(subdet, layer);
0030
0031 if (doseMap_.find(key) == doseMap_.end()) {
0032 return logVal ? -10. : 0.;
0033 }
0034
0035 double r(radius - doseMap_[key].doff_);
0036 double r2(r * r);
0037 double r3(r2 * r);
0038 double r4(r3 * r);
0039
0040 double cellDoseLog10 =
0041 doseMap_[key].a_ + doseMap_[key].b_ * r + doseMap_[key].c_ * r2 + doseMap_[key].d_ * r3 + doseMap_[key].e_ * r4;
0042
0043 return logVal ? cellDoseLog10 * M_LN10 + log(grayToKrad_) : std::pow(10, cellDoseLog10) * grayToKrad_;
0044 }
0045
0046
0047 double HGCalRadiationMap::getFluenceValue(const int subdet, const int layer, const double radius, bool logVal) {
0048 std::pair<int, int> key(subdet, layer);
0049
0050 double r(radius - doseMap_[key].foff_);
0051 double r2(r * r);
0052 double r3(r2 * r);
0053 double r4(r3 * r);
0054
0055 double cellFluenceLog10 =
0056 doseMap_[key].f_ + doseMap_[key].g_ * r + doseMap_[key].h_ * r2 + doseMap_[key].i_ * r3 + doseMap_[key].j_ * r4;
0057 cellFluenceLog10 += fluenceSFlog10_;
0058
0059 return logVal ? cellFluenceLog10 * M_LN10 : std::pow(10, cellFluenceLog10);
0060 }
0061
0062
0063 std::map<std::pair<int, int>, HGCalRadiationMap::DoseParameters> HGCalRadiationMap::readDosePars(
0064 const std::string& fullpath) {
0065 doseParametersMap result;
0066
0067
0068 if (fullpath.empty())
0069 return result;
0070
0071 edm::FileInPath fp(fullpath);
0072 std::ifstream infile(fp.fullPath());
0073 if (!infile.is_open()) {
0074 throw cms::Exception("FileNotFound") << "Unable to open '" << fullpath << "'" << std::endl;
0075 }
0076 std::string line;
0077 while (getline(infile, line)) {
0078 int subdet;
0079 int layer;
0080 DoseParameters dosePars;
0081
0082
0083 std::stringstream linestream(line);
0084 linestream >> subdet >> layer >> dosePars.a_ >> dosePars.b_ >> dosePars.c_ >> dosePars.d_ >> dosePars.e_ >>
0085 dosePars.doff_ >> dosePars.f_ >> dosePars.g_ >> dosePars.h_ >> dosePars.i_ >> dosePars.j_ >> dosePars.foff_;
0086
0087 std::pair<int, int> key(subdet, layer);
0088 result[key] = dosePars;
0089 }
0090 return result;
0091 }