HGCalSimAlgos/src/HGCalRadiationMap.cc

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   //no dose file means no aging
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     //space-separated
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 }