EcalSimAlgos/src/EcalShapeBase.cc

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);              // pure virtual function, implementation may vary for EB/EE/APD ...
0037   m_arraySize = shapeArray.size();  // original data
0038
0039   m_denseArraySize = 10 * m_arraySize;  // dense array with interpolation between data
0040   m_kNBinsPerNSec =
0041       (unsigned int)(10 /
0042                      time_interval);  // used to be an unsigned int = 10 in  < CMSSW10X, should work for time intervals ~0.1, 0.2, 0.5, 1
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)  // cannot do quadratic interpolation at ends
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   // return pulse amplitude for request time in ns
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 }