SiStripObjects/src/SiStripApvSimulationParameters.cc

0001 #include "CondFormats/SiStripObjects/interface/SiStripApvSimulationParameters.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "CLHEP/Random/RandFlat.h"
0004
0005 namespace {
0006   PhysicsTools::Calibration::HistogramF2D calculateXInt(const SiStripApvSimulationParameters::LayerParameters& params) {
0007     auto hXInt = (params.hasEquidistantBinsY()
0008                       ? (params.hasEquidistantBinsZ()
0009                              ? PhysicsTools::Calibration::HistogramF2D(
0010                                    params.numberOfBinsY(), params.rangeY(), params.numberOfBinsZ(), params.rangeZ())
0011                              : PhysicsTools::Calibration::HistogramF2D(
0012                                    params.numberOfBinsY(), params.rangeY(), params.upperLimitsZ()))
0013                       : (params.hasEquidistantBinsZ()
0014                              ? PhysicsTools::Calibration::HistogramF2D(
0015                                    params.upperLimitsY(), params.numberOfBinsZ(), params.rangeZ())
0016                              : PhysicsTools::Calibration::HistogramF2D(params.upperLimitsY(), params.upperLimitsZ())));
0017     for (int i{0}; i != params.numberOfBinsY() + 2; ++i) {
0018       for (int j{0}; j != params.numberOfBinsZ() + 2; ++j) {
0019         float xInt = 0.;
0020         for (int k{0}; k != params.numberOfBinsX() + 2; ++k) {
0021           xInt += params.binContent(k, i, j);
0022         }
0023         hXInt.setBinContent(i, j, xInt);
0024       }
0025     }
0026     return hXInt;
0027   }
0028
0029   float xBinPos(const SiStripApvSimulationParameters::LayerParameters& hist, int iBin, float pos = 0.5) {
0030     // NOTE: does not work for under- and overflow bins (iBin = 0 and iBIn == hist.numberOfBinsX()+1)
0031     if (hist.hasEquidistantBinsX()) {
0032       const auto range = hist.rangeX();
0033       const auto binWidth = (range.max - range.min) / hist.numberOfBinsX();
0034       return range.min + (iBin - 1 + pos) * binWidth;
0035     } else {
0036       return (1. - pos) * hist.upperLimitsX()[iBin - 1] + pos * hist.upperLimitsX()[iBin];
0037     }
0038   }
0039 }  // namespace
0040
0041 void SiStripApvSimulationParameters::calculateIntegrals() {
0042   if (m_barrelParam.size() != m_barrelParam_xInt.size()) {
0043     m_barrelParam_xInt.resize(m_barrelParam.size());
0044     for (unsigned int i{0}; i != m_barrelParam.size(); ++i) {
0045       m_barrelParam_xInt[i] = calculateXInt(m_barrelParam[i]);
0046     }
0047   }
0048   if (m_endcapParam.size() != m_endcapParam_xInt.size()) {
0049     m_endcapParam_xInt.resize(m_endcapParam.size());
0050     for (unsigned int i{0}; i != m_endcapParam.size(); ++i) {
0051       m_endcapParam_xInt[i] = calculateXInt(m_endcapParam[i]);
0052     }
0053   }
0054 }
0055
0056 bool SiStripApvSimulationParameters::putTIB(SiStripApvSimulationParameters::layerid layer,
0057                                             SiStripApvSimulationParameters::LayerParameters&& params) {
0058   if ((layer > m_nTIB) || (layer < 1)) {
0059     edm::LogError("SiStripApvSimulationParameters")
0060         << "[" << __PRETTY_FUNCTION__ << "] layer index " << layer << " out of range [1," << m_nTIB << "]";
0061     return false;
0062   }
0063   m_barrelParam[layer - 1] = params;
0064   m_barrelParam_xInt[layer - 1] = calculateXInt(params);
0065   return true;
0066 }
0067
0068 bool SiStripApvSimulationParameters::putTOB(SiStripApvSimulationParameters::layerid layer,
0069                                             SiStripApvSimulationParameters::LayerParameters&& params) {
0070   if ((layer > m_nTOB) || (layer < 1)) {
0071     edm::LogError("SiStripApvSimulationParameters")
0072         << "[" << __PRETTY_FUNCTION__ << "] layer index " << layer << " out of range [1," << m_nTOB << ")";
0073     return false;
0074   }
0075   m_barrelParam[m_nTIB + layer - 1] = params;
0076   m_barrelParam_xInt[m_nTIB + layer - 1] = calculateXInt(params);
0077   return true;
0078 }
0079
0080 bool SiStripApvSimulationParameters::putTID(SiStripApvSimulationParameters::layerid wheel,
0081                                             SiStripApvSimulationParameters::LayerParameters&& params) {
0082   if ((wheel > m_nTID) || (wheel < 1)) {
0083     edm::LogError("SiStripApvSimulationParameters")
0084         << "[" << __PRETTY_FUNCTION__ << "] wheel index " << wheel << " out of range [1," << m_nTID << "]";
0085     return false;
0086   }
0087   m_endcapParam[wheel - 1] = params;
0088   m_endcapParam_xInt[wheel - 1] = calculateXInt(params);
0089   return true;
0090 }
0091
0092 bool SiStripApvSimulationParameters::putTEC(SiStripApvSimulationParameters::layerid wheel,
0093                                             SiStripApvSimulationParameters::LayerParameters&& params) {
0094   if ((wheel > m_nTEC) || (wheel < 1)) {
0095     edm::LogError("SiStripApvSimulationParameters")
0096         << "[" << __PRETTY_FUNCTION__ << "] wheel index " << wheel << " out of range [1," << m_nTEC << ")";
0097     return false;
0098   }
0099   m_endcapParam[m_nTID + wheel - 1] = params;
0100   m_endcapParam_xInt[m_nTID + wheel - 1] = calculateXInt(params);
0101   return true;
0102 }
0103
0104 float SiStripApvSimulationParameters::sampleBarrel(layerid layerIdx,
0105                                                    float z,
0106                                                    float pu,
0107                                                    CLHEP::HepRandomEngine* engine) const {
0108   if (m_barrelParam.size() != m_barrelParam_xInt.size()) {
0109     throw cms::Exception("LogicError") << "x-integrals of 3D histograms have not been calculated";
0110   }
0111   const auto layerParam = m_barrelParam[layerIdx];
0112   const int ip = layerParam.findBinY(pu);
0113   const int iz = layerParam.findBinZ(z);
0114   const float norm = m_barrelParam_xInt[layerIdx].binContent(ip, iz);
0115   const auto val = CLHEP::RandFlat::shoot(engine) * norm;
0116   if (val < layerParam.binContent(0, ip, iz)) {  // underflow
0117     return layerParam.rangeX().min;
0118   } else if (norm - val < layerParam.binContent(layerParam.numberOfBinsX() + 1, ip, iz)) {  // overflow
0119     return layerParam.rangeX().max;
0120   } else {  // loop over bins, return center of our bin
0121     float sum = layerParam.binContent(0, ip, iz);
0122     for (int i{1}; i != layerParam.numberOfBinsX() + 1; ++i) {
0123       sum += layerParam.binContent(i, ip, iz);
0124       if (sum > val) {
0125         return xBinPos(layerParam, i, (sum - val) / layerParam.binContent(i, ip, iz));
0126       }
0127     }
0128   }
0129   throw cms::Exception("LogicError") << "Problem drawing a random number from the distribution";
0130 }
0131
0132 float SiStripApvSimulationParameters::sampleEndcap(layerid wheelIdx,
0133                                                    float r,
0134                                                    float pu,
0135                                                    CLHEP::HepRandomEngine* engine) const {
0136   if (m_endcapParam.size() != m_endcapParam_xInt.size()) {
0137     throw cms::Exception("LogicError") << "x-integrals of 3D histograms have not been calculated";
0138   }
0139   const auto layerParam = m_endcapParam[wheelIdx];
0140   const int ip = layerParam.findBinY(pu);
0141   const int ir = layerParam.findBinZ(r);
0142   const float norm = m_endcapParam_xInt[wheelIdx].binContent(ip, ir);
0143   const auto val = CLHEP::RandFlat::shoot(engine) * norm;
0144   if (val < layerParam.binContent(0, ip, ir)) {  // underflow
0145     return layerParam.rangeX().min;
0146   } else if (norm - val < layerParam.binContent(layerParam.numberOfBinsX() + 1, ip, ir)) {  // overflow
0147     return layerParam.rangeX().max;
0148   } else {  // loop over bins, return center of our bin
0149     float sum = layerParam.binContent(0, ip, ir);
0150     for (int i{1}; i != layerParam.numberOfBinsX() + 1; ++i) {
0151       sum += layerParam.binContent(i, ip, ir);
0152       if (sum > val) {
0153         return xBinPos(layerParam, i, (sum - val) / layerParam.binContent(i, ip, ir));
0154       }
0155     }
0156   }
0157   throw cms::Exception("LogicError") << "Problem drawing a random number from the distribution";
0158 }