QIE8Simulator.cc

Go to the documentation of this file.

#include <cmath>
#include <cfloat>

#include "CalibCalorimetry/HcalAlgos/interface/QIE8Simulator.h"

QIE8Simulator::QIE8Simulator()
    : preampOutputCut_(-1.0e100),
      inputGain_(1.0),
      outputGain_(1.0),
      runCount_(0),
      chargeNode_(AbsElectronicODERHS::invalidNode),
      integrateToGetCharge_(false),
      useCubic_(false) {}

QIE8Simulator::QIE8Simulator(const AbsElectronicODERHS& model,
                             const unsigned chargeNode,
                             const bool interpolateCubic,
                             const double preampOutputCut,
                             const double inputGain,
                             const double outputGain)
    : solver_(model),
      preampOutputCut_(preampOutputCut),
      inputGain_(inputGain),
      outputGain_(outputGain),
      runCount_(0),
      chargeNode_(chargeNode),
      useCubic_(interpolateCubic) {
  if (chargeNode >= AbsElectronicODERHS::invalidNode)
    throw cms::Exception("In QIE8Simulator constructor: invalid charge collection node");
  integrateToGetCharge_ = chargeNode == model.outputNode();
  validateGain();
  zeroInitialConditions();
}

unsigned QIE8Simulator::run(
    const double dt, const double tstop, const double tDigitize, double* TS, const unsigned lenTS) {
  if (chargeNode_ >= AbsElectronicODERHS::invalidNode)
    throw cms::Exception("In QIE8Simulator::run: preamp model is not set");

  // Check arguments for validity
  if (dt <= 0.0)
    throw cms::Exception("In QIE8Simulator::run: invalid time step");

  if (tstop < dt)
    throw cms::Exception("In QIE8Simulator::run: invalid stopping time");

  if (lenTS && tDigitize < 0.0)
    throw cms::Exception("In QIE8Simulator::run: can't digitize at t < 0");

  // Determine the number of time steps
  const double dsteps = tstop / dt;
  unsigned n = dsteps;
  if (dsteps != static_cast<double>(n))
    ++n;
  if (n >= maxlen)
    n = maxlen - 1;

  // Run the simulation
  AbsElectronicODERHS& rhs = modifiableRHS();
  const unsigned numNodes = rhs.numberOfNodes();
  if (numNodes) {
    assert(initialConditions_.size() == numNodes);
    solver_.run(&initialConditions_[0], numNodes, dt, n);
  } else {
    // Special situation: the simulation does not
    // need to solve any ODE. Instead, it will fill
    // out the output directly.
    if (!(integrateToGetCharge_ && chargeNode_ == 0U))
      throw cms::Exception(
          "In QIE8Simulator::run: "
          "invalid mode of operation");
    const unsigned runLen = n + 1U;
    if (historyBuffer_.size() < runLen)
      historyBuffer_.resize(runLen);
    double* hbuf = &historyBuffer_[0];
    for (unsigned istep = 0; istep < runLen; ++istep, ++hbuf)
      rhs.calc(istep * dt, nullptr, 0U, hbuf);
    solver_.setHistory(dt, &historyBuffer_[0], 1U, runLen);
  }

  // Truncate the preamp output if this will affect subsequent results
  if (integrateToGetCharge_ && preampOutputCut_ > -0.9e100)
    solver_.truncateCoordinate(chargeNode_, preampOutputCut_, DBL_MAX);

  // Digitize the accumulated charge
  unsigned filled = 0;
  if (lenTS) {
    assert(TS);
    const double lastTStop = solver_.lastMaxT();
    const double tsWidth = this->adcTSWidth();
    double oldCharge = getCharge(tDigitize);
    for (unsigned its = 0; its < lenTS; ++its) {
      const double t0 = tDigitize + its * tsWidth;
      if (t0 < lastTStop) {
        double t1 = t0 + tsWidth;
        if (t1 > lastTStop)
          t1 = lastTStop;
        else
          ++filled;
        const double q = getCharge(t1);
        TS[its] = (q - oldCharge) * outputGain_;
        oldCharge = q;
      } else
        TS[its] = 0.0;
    }
  }
  ++runCount_;
  return filled;
}

double QIE8Simulator::preampOutput(const double t) const {
  if (!runCount_)
    throw cms::Exception("In QIE8Simulator::preampOutput: please run the simulation first");
  const unsigned preampNode = getRHS().outputNode();
  if (preampNode >= AbsElectronicODERHS::invalidNode)
    return 0.0;
  else
    return outputGain_ * solver_.interpolateCoordinate(preampNode, t, useCubic_);
}

double QIE8Simulator::controlOutput(const double t) const {
  if (!runCount_)
    throw cms::Exception("In QIE8Simulator::controlOutput: please run the simulation first");
  const unsigned controlNode = getRHS().controlNode();
  if (controlNode >= AbsElectronicODERHS::invalidNode)
    return 0.0;
  else
    return solver_.interpolateCoordinate(controlNode, t, useCubic_);
}

double QIE8Simulator::preampPeakTime() const {
  if (!runCount_)
    throw cms::Exception("In QIE8Simulator::preampPeakTime: please run the simulation first");
  const unsigned preampNode = getRHS().outputNode();
  if (preampNode >= AbsElectronicODERHS::invalidNode)
    throw cms::Exception("In QIE8Simulator::preampPeakTime: no preamp node in the circuit");
  return solver_.getPeakTime(preampNode);
}

void QIE8Simulator::setInitialConditions(const double* values, const unsigned sz) {
  const unsigned nExpected = getRHS().numberOfNodes();
  if (sz != nExpected)
    throw cms::Exception(
        "In QIE8Simulator::setInitialConditions: unexpected number "
        "of initial conditions");
  assert(sz == initialConditions_.size());
  if (sz) {
    double* c = &initialConditions_[0];
    for (unsigned i = 0; i < sz; ++i)
      *c++ = *values++;
  }
}

void QIE8Simulator::setRHS(const AbsElectronicODERHS& rhs,
                           const unsigned chargeNode,
                           const bool useCubicInterpolation) {
  if (chargeNode >= AbsElectronicODERHS::invalidNode)
    throw cms::Exception("In QIE8Simulator::setRHS invalid charge collection node");
  solver_.setRHS(rhs);
  chargeNode_ = chargeNode;
  integrateToGetCharge_ = chargeNode == rhs.outputNode();
  useCubic_ = useCubicInterpolation;
  zeroInitialConditions();
}

void QIE8Simulator::zeroInitialConditions() {
  const unsigned sz = getRHS().numberOfNodes();
  if (initialConditions_.size() != sz)
    initialConditions_.resize(sz);
  if (sz) {
    double* c = &initialConditions_[0];
    for (unsigned i = 0; i < sz; ++i)
      *c++ = 0.0;
  }
}

void QIE8Simulator::validateGain() const {
  if (inputGain_ <= 0.0)
    throw cms::Exception("In QIE8Simulator::validateGain: invalid input gain");
}

double QIE8Simulator::lastStopTime() const {
  if (!runCount_)
    throw cms::Exception("In QIE8Simulator::lastStopTime: please run the simulation first");
  return solver_.lastMaxT();
}

double QIE8Simulator::totalIntegratedCharge(const double t) const {
  if (!runCount_)
    throw cms::Exception(
        "In QIE8Simulator::totalIntegratedCharge: "
        "please run the simulation first");
  return outputGain_ * getCharge(t);
}