Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:28:32

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) { buildMe(&evtSetup); }
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) {
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   for (unsigned int i(0); i != shapeArray.size(); ++i) {
0053     shapeArray[i] = shapeArray[i] / maxelt;
0054   }
0055 
0056   const double thresh(threshold() / maxelt);
0057 
0058   const double delta(m_qNSecPerBin / 2.);
0059 
0060   for (unsigned int denseIndex(0); denseIndex != m_denseArraySize; ++denseIndex) {
0061     const double xb((denseIndex + 0.5) * m_qNSecPerBin);
0062 
0063     const unsigned int ibin(denseIndex / 10);
0064 
0065     double value = 0.0;
0066     double deriv = 0.0;
0067 
0068     if (0 == ibin || shapeArray.size() == 1 + ibin)  // cannot do quadratic interpolation at ends
0069     {
0070       value = shapeArray[ibin];
0071       deriv = 0.0;
0072     } else {
0073       const double x(xb - (ibin + 0.5) * time_interval);
0074       const double f1(shapeArray[ibin - 1]);
0075       const double f2(shapeArray[ibin]);
0076       const double f3(shapeArray[ibin + 1]);
0077       const double a(f2);
0078       const double b((f3 - f1) / (2. * time_interval));
0079       const double c(((f1 + f3) / 2. - f2) / (time_interval * time_interval));
0080       value = a + b * x + c * x * x;
0081       deriv = (b + 2 * c * x) / delta;
0082     }
0083 
0084     m_shape[denseIndex] = value;
0085     m_deriv[denseIndex] = deriv;
0086 
0087     if (0 < denseIndex && thresh < value && 0 == m_firstIndexOverThreshold) {
0088       m_firstIndexOverThreshold = denseIndex - 1;
0089       m_firstTimeOverThreshold = m_firstIndexOverThreshold * m_qNSecPerBin;
0090     }
0091 
0092     if (m_shape[m_indexOfMax] < value) {
0093       m_indexOfMax = denseIndex;
0094     }
0095   }
0096   m_timeOfMax = m_indexOfMax * m_qNSecPerBin;
0097 }
0098 
0099 unsigned int EcalShapeBase::timeIndex(double aTime) const {
0100   const int index(m_firstIndexOverThreshold + (unsigned int)(aTime * m_kNBinsPerNSec + 0.5));
0101 
0102   const bool bad((int)m_firstIndexOverThreshold > index || (int)m_denseArraySize <= index);
0103 
0104   if ((int)m_denseArraySize <= index) {
0105     LogDebug("EcalShapeBase") << " ECAL MGPA shape requested for out of range time " << aTime;
0106   }
0107   return (bad ? m_denseArraySize : (unsigned int)index);
0108 }
0109 
0110 double EcalShapeBase::operator()(double aTime) const {
0111   // return pulse amplitude for request time in ns
0112 
0113   const unsigned int index(timeIndex(aTime));
0114   return (m_denseArraySize == index ? 0 : m_shape[index]);
0115 }
0116 
0117 double EcalShapeBase::derivative(double aTime) const {
0118   const unsigned int index(timeIndex(aTime));
0119   return (m_denseArraySize == index ? 0 : m_deriv[index]);
0120 }
0121 
0122 void EcalShapeBase::m_shape_print(const char* fileName) {
0123   std::ofstream fs;
0124   fs.open(fileName);
0125   for (auto i : m_shape)
0126     fs << "vec.push_back(" << i << ");\n";
0127   fs.close();
0128 }