File indexing completed on 2023-04-02 22:48:16
0001 #include "SimCalorimetry/EcalSimAlgos/interface/EcalShapeBase.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 #include <cmath>
0005 #include <algorithm>
0006
0007 #include <iostream>
0008
0009 EcalShapeBase::~EcalShapeBase() {}
0010
0011 EcalShapeBase::EcalShapeBase(bool useDBShape)
0012 : m_useDBShape(useDBShape),
0013 m_firstIndexOverThreshold(0),
0014 m_firstTimeOverThreshold(0.0),
0015 m_indexOfMax(0),
0016 m_timeOfMax(0.0),
0017 m_thresh(0.0) {}
0018
0019 void EcalShapeBase::setEventSetup(const edm::EventSetup& evtSetup, bool normalize) { buildMe(&evtSetup, normalize); }
0020
0021 double EcalShapeBase::timeOfThr() const { return m_firstTimeOverThreshold; }
0022
0023 double EcalShapeBase::timeOfMax() const { return m_timeOfMax; }
0024
0025 double EcalShapeBase::timeToRise() const { return timeOfMax() - timeOfThr(); }
0026
0027 double EcalShapeBase::threshold() const { return m_thresh; }
0028
0029 void EcalShapeBase::buildMe(const edm::EventSetup* evtSetup, bool normalize) {
0030 DVec shapeArray;
0031
0032 float time_interval = 0;
0033 fillShape(time_interval,
0034 m_thresh,
0035 shapeArray,
0036 evtSetup);
0037 m_arraySize = shapeArray.size();
0038
0039 m_denseArraySize = 10 * m_arraySize;
0040 m_kNBinsPerNSec =
0041 (unsigned int)(10 /
0042 time_interval);
0043 m_qNSecPerBin = time_interval / 10.;
0044
0045 m_deriv.resize(m_denseArraySize);
0046 m_shape.resize(m_denseArraySize);
0047
0048 const double maxel(*max_element(shapeArray.begin(), shapeArray.end()));
0049
0050 const double maxelt(1.e-5 < maxel ? maxel : 1);
0051
0052 if (normalize) {
0053 for (unsigned int i(0); i != shapeArray.size(); ++i) {
0054 shapeArray[i] = shapeArray[i] / maxelt;
0055 }
0056 }
0057
0058 const double thresh(threshold() / maxelt);
0059
0060 const double delta(m_qNSecPerBin / 2.);
0061
0062 for (unsigned int denseIndex(0); denseIndex != m_denseArraySize; ++denseIndex) {
0063 const double xb((denseIndex + 0.5) * m_qNSecPerBin);
0064
0065 const unsigned int ibin(denseIndex / 10);
0066
0067 double value = 0.0;
0068 double deriv = 0.0;
0069
0070 if (0 == ibin || shapeArray.size() == 1 + ibin)
0071 {
0072 value = shapeArray[ibin];
0073 deriv = 0.0;
0074 } else {
0075 const double x(xb - (ibin + 0.5) * time_interval);
0076 const double f1(shapeArray[ibin - 1]);
0077 const double f2(shapeArray[ibin]);
0078 const double f3(shapeArray[ibin + 1]);
0079 const double a(f2);
0080 const double b((f3 - f1) / (2. * time_interval));
0081 const double c(((f1 + f3) / 2. - f2) / (time_interval * time_interval));
0082 value = a + b * x + c * x * x;
0083 deriv = (b + 2 * c * x) / delta;
0084 }
0085
0086 m_shape[denseIndex] = value;
0087 m_deriv[denseIndex] = deriv;
0088
0089 if (0 < denseIndex && thresh < value && 0 == m_firstIndexOverThreshold) {
0090 m_firstIndexOverThreshold = denseIndex - 1;
0091 m_firstTimeOverThreshold = m_firstIndexOverThreshold * m_qNSecPerBin;
0092 }
0093
0094 if (m_shape[m_indexOfMax] < value) {
0095 m_indexOfMax = denseIndex;
0096 }
0097 }
0098 m_timeOfMax = m_indexOfMax * m_qNSecPerBin;
0099 }
0100
0101 unsigned int EcalShapeBase::timeIndex(double aTime) const {
0102 const int index(m_firstIndexOverThreshold + (unsigned int)(aTime * m_kNBinsPerNSec + 0.5));
0103
0104 const bool bad((int)m_firstIndexOverThreshold > index || (int)m_denseArraySize <= index);
0105
0106 if ((int)m_denseArraySize <= index) {
0107 LogDebug("EcalShapeBase") << " ECAL MGPA shape requested for out of range time " << aTime;
0108 }
0109 return (bad ? m_denseArraySize : (unsigned int)index);
0110 }
0111
0112 double EcalShapeBase::operator()(double aTime) const {
0113
0114
0115 const unsigned int index(timeIndex(aTime));
0116 return (m_denseArraySize == index ? 0 : m_shape[index]);
0117 }
0118
0119 double EcalShapeBase::derivative(double aTime) const {
0120 const unsigned int index(timeIndex(aTime));
0121 return (m_denseArraySize == index ? 0 : m_deriv[index]);
0122 }
0123
0124 void EcalShapeBase::m_shape_print(const char* fileName) const {
0125 std::ofstream fs;
0126 fs.open(fileName);
0127 fs << "{\n";
0128 for (auto i : m_shape)
0129 fs << "vec.push_back(" << i << ");\n";
0130 fs << "}\n";
0131 fs.close();
0132 }