EcalShapeBase.cc

Go to the documentation of this file.

#include "SimCalorimetry/EcalSimAlgos/interface/EcalShapeBase.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
 
#include <cmath>
#include <algorithm>
 
#include <iostream>
 
EcalShapeBase::~EcalShapeBase() {}
 
EcalShapeBase::EcalShapeBase(bool useDBShape)
    : m_useDBShape(useDBShape),
      m_firstIndexOverThreshold(0),
      m_firstTimeOverThreshold(0.0),
      m_indexOfMax(0),
      m_timeOfMax(0.0),
      m_thresh(0.0) {}
 
void EcalShapeBase::setEventSetup(const edm::EventSetup& evtSetup) { buildMe(&evtSetup); }
 
double EcalShapeBase::timeOfThr() const { return m_firstTimeOverThreshold; }
 
double EcalShapeBase::timeOfMax() const { return m_timeOfMax; }
 
double EcalShapeBase::timeToRise() const { return timeOfMax() - timeOfThr(); }
 
double EcalShapeBase::threshold() const { return m_thresh; }
 
void EcalShapeBase::buildMe(const edm::EventSetup* evtSetup) {
  DVec shapeArray;
 
  float time_interval = 0;
  fillShape(time_interval,
            m_thresh,
            shapeArray,
            evtSetup);              // pure virtual function, implementation may vary for EB/EE/APD ...
  m_arraySize = shapeArray.size();  // original data
 
  m_denseArraySize = 10 * m_arraySize;  // dense array with interpolation between data
  m_kNBinsPerNSec =
      (unsigned int)(10 /
                     time_interval);  // used to be an unsigned int = 10 in  < CMSSW10X, should work for time intervals ~0.1, 0.2, 0.5, 1
  m_qNSecPerBin = time_interval / 10.;
 
  m_deriv.resize(m_denseArraySize);
  m_shape.resize(m_denseArraySize);
 
  const double maxel(*max_element(shapeArray.begin(), shapeArray.end()));
 
  const double maxelt(1.e-5 < maxel ? maxel : 1);
 
  for (unsigned int i(0); i != shapeArray.size(); ++i) {
    shapeArray[i] = shapeArray[i] / maxelt;
  }
 
  const double thresh(threshold() / maxelt);
 
  const double delta(m_qNSecPerBin / 2.);
 
  for (unsigned int denseIndex(0); denseIndex != m_denseArraySize; ++denseIndex) {
    const double xb((denseIndex + 0.5) * m_qNSecPerBin);
 
    const unsigned int ibin(denseIndex / 10);
 
    double value = 0.0;
    double deriv = 0.0;
 
    if (0 == ibin || shapeArray.size() == 1 + ibin)  // cannot do quadratic interpolation at ends
    {
      value = shapeArray[ibin];
      deriv = 0.0;
    } else {
      const double x(xb - (ibin + 0.5) * time_interval);
      const double f1(shapeArray[ibin - 1]);
      const double f2(shapeArray[ibin]);
      const double f3(shapeArray[ibin + 1]);
      const double a(f2);
      const double b((f3 - f1) / (2. * time_interval));
      const double c(((f1 + f3) / 2. - f2) / (time_interval * time_interval));
      value = a + b * x + c * x * x;
      deriv = (b + 2 * c * x) / delta;
    }
 
    m_shape[denseIndex] = value;
    m_deriv[denseIndex] = deriv;
 
    if (0 < denseIndex && thresh < value && 0 == m_firstIndexOverThreshold) {
      m_firstIndexOverThreshold = denseIndex - 1;
      m_firstTimeOverThreshold = m_firstIndexOverThreshold * m_qNSecPerBin;
    }
 
    if (m_shape[m_indexOfMax] < value) {
      m_indexOfMax = denseIndex;
    }
  }
  m_timeOfMax = m_indexOfMax * m_qNSecPerBin;
}
 
unsigned int EcalShapeBase::timeIndex(double aTime) const {
  const int index(m_firstIndexOverThreshold + (unsigned int)(aTime * m_kNBinsPerNSec + 0.5));
 
  const bool bad((int)m_firstIndexOverThreshold > index || (int)m_denseArraySize <= index);
 
  if ((int)m_denseArraySize <= index) {
    LogDebug("EcalShapeBase") << " ECAL MGPA shape requested for out of range time " << aTime;
  }
  return (bad ? m_denseArraySize : (unsigned int)index);
}
 
double EcalShapeBase::operator()(double aTime) const {
  // return pulse amplitude for request time in ns
 
  const unsigned int index(timeIndex(aTime));
  return (m_denseArraySize == index ? 0 : m_shape[index]);
}
 
double EcalShapeBase::derivative(double aTime) const {
  const unsigned int index(timeIndex(aTime));
  return (m_denseArraySize == index ? 0 : m_deriv[index]);
}
 
void EcalShapeBase::m_shape_print(const char* fileName) {
  std::ofstream fs;
  fs.open(fileName);
  for (auto i : m_shape)
    fs << "vec.push_back(" << i << ");\n";
  fs.close();
}