HcalAlgos/src/QIE8Simulator.cc

0001 #include <cmath>
0002 #include <cfloat>
0003
0004 #include "CalibCalorimetry/HcalAlgos/interface/QIE8Simulator.h"
0005
0006 QIE8Simulator::QIE8Simulator()
0007     : preampOutputCut_(-1.0e100),
0008       inputGain_(1.0),
0009       outputGain_(1.0),
0010       runCount_(0),
0011       chargeNode_(AbsElectronicODERHS::invalidNode),
0012       integrateToGetCharge_(false),
0013       useCubic_(false) {}
0014
0015 QIE8Simulator::QIE8Simulator(const AbsElectronicODERHS& model,
0016                              const unsigned chargeNode,
0017                              const bool interpolateCubic,
0018                              const double preampOutputCut,
0019                              const double inputGain,
0020                              const double outputGain)
0021     : solver_(model),
0022       preampOutputCut_(preampOutputCut),
0023       inputGain_(inputGain),
0024       outputGain_(outputGain),
0025       runCount_(0),
0026       chargeNode_(chargeNode),
0027       useCubic_(interpolateCubic) {
0028   if (chargeNode >= AbsElectronicODERHS::invalidNode)
0029     throw cms::Exception("In QIE8Simulator constructor: invalid charge collection node");
0030   integrateToGetCharge_ = chargeNode == model.outputNode();
0031   validateGain();
0032   zeroInitialConditions();
0033 }
0034
0035 unsigned QIE8Simulator::run(
0036     const double dt, const double tstop, const double tDigitize, double* TS, const unsigned lenTS) {
0037   if (chargeNode_ >= AbsElectronicODERHS::invalidNode)
0038     throw cms::Exception("In QIE8Simulator::run: preamp model is not set");
0039
0040   // Check arguments for validity
0041   if (dt <= 0.0)
0042     throw cms::Exception("In QIE8Simulator::run: invalid time step");
0043
0044   if (tstop < dt)
0045     throw cms::Exception("In QIE8Simulator::run: invalid stopping time");
0046
0047   if (lenTS && tDigitize < 0.0)
0048     throw cms::Exception("In QIE8Simulator::run: can't digitize at t < 0");
0049
0050   // Determine the number of time steps
0051   const double dsteps = tstop / dt;
0052   unsigned n = dsteps;
0053   if (dsteps != static_cast<double>(n))
0054     ++n;
0055   if (n >= maxlen)
0056     n = maxlen - 1;
0057
0058   // Run the simulation
0059   AbsElectronicODERHS& rhs = modifiableRHS();
0060   const unsigned numNodes = rhs.numberOfNodes();
0061   if (numNodes) {
0062     assert(initialConditions_.size() == numNodes);
0063     solver_.run(&initialConditions_[0], numNodes, dt, n);
0064   } else {
0065     // Special situation: the simulation does not
0066     // need to solve any ODE. Instead, it will fill
0067     // out the output directly.
0068     if (!(integrateToGetCharge_ && chargeNode_ == 0U))
0069       throw cms::Exception(
0070           "In QIE8Simulator::run: "
0071           "invalid mode of operation");
0072     const unsigned runLen = n + 1U;
0073     if (historyBuffer_.size() < runLen)
0074       historyBuffer_.resize(runLen);
0075     double* hbuf = &historyBuffer_[0];
0076     for (unsigned istep = 0; istep < runLen; ++istep, ++hbuf)
0077       rhs.calc(istep * dt, nullptr, 0U, hbuf);
0078     solver_.setHistory(dt, &historyBuffer_[0], 1U, runLen);
0079   }
0080
0081   // Truncate the preamp output if this will affect subsequent results
0082   if (integrateToGetCharge_ && preampOutputCut_ > -0.9e100)
0083     solver_.truncateCoordinate(chargeNode_, preampOutputCut_, DBL_MAX);
0084
0085   // Digitize the accumulated charge
0086   unsigned filled = 0;
0087   if (lenTS) {
0088     assert(TS);
0089     const double lastTStop = solver_.lastMaxT();
0090     const double tsWidth = this->adcTSWidth();
0091     double oldCharge = getCharge(tDigitize);
0092     for (unsigned its = 0; its < lenTS; ++its) {
0093       const double t0 = tDigitize + its * tsWidth;
0094       if (t0 < lastTStop) {
0095         double t1 = t0 + tsWidth;
0096         if (t1 > lastTStop)
0097           t1 = lastTStop;
0098         else
0099           ++filled;
0100         const double q = getCharge(t1);
0101         TS[its] = (q - oldCharge) * outputGain_;
0102         oldCharge = q;
0103       } else
0104         TS[its] = 0.0;
0105     }
0106   }
0107   ++runCount_;
0108   return filled;
0109 }
0110
0111 double QIE8Simulator::preampOutput(const double t) const {
0112   if (!runCount_)
0113     throw cms::Exception("In QIE8Simulator::preampOutput: please run the simulation first");
0114   const unsigned preampNode = getRHS().outputNode();
0115   if (preampNode >= AbsElectronicODERHS::invalidNode)
0116     return 0.0;
0117   else
0118     return outputGain_ * solver_.interpolateCoordinate(preampNode, t, useCubic_);
0119 }
0120
0121 double QIE8Simulator::controlOutput(const double t) const {
0122   if (!runCount_)
0123     throw cms::Exception("In QIE8Simulator::controlOutput: please run the simulation first");
0124   const unsigned controlNode = getRHS().controlNode();
0125   if (controlNode >= AbsElectronicODERHS::invalidNode)
0126     return 0.0;
0127   else
0128     return solver_.interpolateCoordinate(controlNode, t, useCubic_);
0129 }
0130
0131 double QIE8Simulator::preampPeakTime() const {
0132   if (!runCount_)
0133     throw cms::Exception("In QIE8Simulator::preampPeakTime: please run the simulation first");
0134   const unsigned preampNode = getRHS().outputNode();
0135   if (preampNode >= AbsElectronicODERHS::invalidNode)
0136     throw cms::Exception("In QIE8Simulator::preampPeakTime: no preamp node in the circuit");
0137   return solver_.getPeakTime(preampNode);
0138 }
0139
0140 void QIE8Simulator::setInitialConditions(const double* values, const unsigned sz) {
0141   const unsigned nExpected = getRHS().numberOfNodes();
0142   if (sz != nExpected)
0143     throw cms::Exception(
0144         "In QIE8Simulator::setInitialConditions: unexpected number "
0145         "of initial conditions");
0146   assert(sz == initialConditions_.size());
0147   if (sz) {
0148     double* c = &initialConditions_[0];
0149     for (unsigned i = 0; i < sz; ++i)
0150       *c++ = *values++;
0151   }
0152 }
0153
0154 void QIE8Simulator::setRHS(const AbsElectronicODERHS& rhs,
0155                            const unsigned chargeNode,
0156                            const bool useCubicInterpolation) {
0157   if (chargeNode >= AbsElectronicODERHS::invalidNode)
0158     throw cms::Exception("In QIE8Simulator::setRHS invalid charge collection node");
0159   solver_.setRHS(rhs);
0160   chargeNode_ = chargeNode;
0161   integrateToGetCharge_ = chargeNode == rhs.outputNode();
0162   useCubic_ = useCubicInterpolation;
0163   zeroInitialConditions();
0164 }
0165
0166 void QIE8Simulator::zeroInitialConditions() {
0167   const unsigned sz = getRHS().numberOfNodes();
0168   if (initialConditions_.size() != sz)
0169     initialConditions_.resize(sz);
0170   if (sz) {
0171     double* c = &initialConditions_[0];
0172     for (unsigned i = 0; i < sz; ++i)
0173       *c++ = 0.0;
0174   }
0175 }
0176
0177 void QIE8Simulator::validateGain() const {
0178   if (inputGain_ <= 0.0)
0179     throw cms::Exception("In QIE8Simulator::validateGain: invalid input gain");
0180 }
0181
0182 double QIE8Simulator::lastStopTime() const {
0183   if (!runCount_)
0184     throw cms::Exception("In QIE8Simulator::lastStopTime: please run the simulation first");
0185   return solver_.lastMaxT();
0186 }
0187
0188 double QIE8Simulator::totalIntegratedCharge(const double t) const {
0189   if (!runCount_)
0190     throw cms::Exception(
0191         "In QIE8Simulator::totalIntegratedCharge: "
0192         "please run the simulation first");
0193   return outputGain_ * getCharge(t);
0194 }