HcalAlgos/src/PadeTableODE.cc

0001 #include <cassert>
0002 #include "FWCore/Utilities/interface/Exception.h"
0003
0004 #include "CalibCalorimetry/HcalAlgos/interface/PadeTableODE.h"
0005
0006 PadeTableODE::PadeTableODE(const unsigned padeRow, const unsigned padeColumn) : row_(padeRow), col_(padeColumn) {
0007   if (row_ > 2U)
0008     throw cms::Exception("In PadeTableODE constructor: Pade table row number out of range");
0009   if (col_ > 3U)
0010     throw cms::Exception("In PadeTableODE constructor: Pade table column number out of range");
0011 }
0012
0013 void PadeTableODE::calculate(const double tau,
0014                              const double currentIn,
0015                              const double dIdt,
0016                              const double d2Id2t,
0017                              const double* x,
0018                              const unsigned lenX,
0019                              const unsigned firstNode,
0020                              double* derivative) const {
0021   // Check input sanity
0022   if (lenX < firstNode + col_)
0023     throw cms::Exception("In PadeTableODE::calculate: insufficient number of variables");
0024   if (tau <= 0.0)
0025     throw cms::Exception("In PadeTableODE::calculate: delay time is not positive");
0026   if (col_)
0027     assert(x);
0028   assert(derivative);
0029
0030   switch (col_) {
0031     case 0U:
0032       // Special case: no ODE to solve
0033       derivative[firstNode] = 0.0;
0034       switch (row_) {
0035         case 2U:
0036           derivative[firstNode] += 0.5 * tau * tau * d2Id2t;
0037           [[fallthrough]];
0038         case 1U:
0039           derivative[firstNode] -= tau * dIdt;
0040           [[fallthrough]];
0041         case 0U:
0042           derivative[firstNode] += currentIn;
0043           break;
0044
0045         default:
0046           assert(0);
0047       }
0048       break;
0049
0050     case 1U:
0051       // First order ODE to solve
0052       switch (row_) {
0053         case 0U:
0054           derivative[firstNode] = (currentIn - x[firstNode]) / tau;
0055           break;
0056
0057         case 1U:
0058           derivative[firstNode] = 2.0 * (currentIn - x[firstNode]) / tau - dIdt;
0059           break;
0060
0061         case 2U:
0062           derivative[firstNode] = 3.0 * (currentIn - x[firstNode]) / tau - 2.0 * dIdt + 0.5 * tau * d2Id2t;
0063           break;
0064
0065         default:
0066           assert(0);
0067       }
0068       break;
0069
0070     case 2U:
0071       // Second order ODE to solve
0072       derivative[firstNode] = x[firstNode + 1];
0073       switch (row_) {
0074         case 0U:
0075           derivative[firstNode + 1] = 2.0 * (currentIn - x[firstNode] - tau * x[firstNode + 1]) / tau / tau;
0076           break;
0077
0078         case 1U:
0079           derivative[firstNode + 1] =
0080               (6.0 * (currentIn - x[firstNode]) - 2.0 * tau * dIdt - 4.0 * tau * x[firstNode + 1]) / tau / tau;
0081           break;
0082
0083         case 2U:
0084           derivative[firstNode + 1] =
0085               12.0 * (currentIn - x[firstNode]) / tau / tau - 6.0 * (x[firstNode + 1] + dIdt) / tau + d2Id2t;
0086           break;
0087
0088         default:
0089           assert(0);
0090       }
0091       break;
0092
0093     case 3U:
0094       // Third order ODE to solve
0095       derivative[firstNode] = x[firstNode + 1];
0096       derivative[firstNode + 1] = x[firstNode + 2];
0097       switch (row_) {
0098         case 0U:
0099           derivative[firstNode + 2] =
0100               6.0 * (currentIn - x[firstNode] - tau * x[firstNode + 1] - 0.5 * tau * tau * x[firstNode + 2]) / tau /
0101               tau / tau;
0102           break;
0103
0104         case 1U:
0105           derivative[firstNode + 2] = 24.0 / tau / tau / tau *
0106                                       (currentIn - x[firstNode] - 0.25 * tau * dIdt - 0.75 * tau * x[firstNode + 1] -
0107                                        0.25 * tau * tau * x[firstNode + 2]);
0108           break;
0109
0110         case 2U:
0111           derivative[firstNode + 2] = 60.0 / tau / tau / tau *
0112                                       (currentIn - x[firstNode] - 0.4 * tau * dIdt + 0.05 * tau * tau * d2Id2t -
0113                                        0.6 * tau * x[firstNode + 1] - 0.15 * tau * tau * x[firstNode + 2]);
0114           break;
0115
0116         default:
0117           assert(0);
0118       }
0119       break;
0120
0121     default:
0122       //
0123       // In principle, it is possible to proceed a bit further, but
0124       // we will soon encounter difficulties. For example, row 0 and
0125       // column 4 is going to generate a 4th order differential
0126       // equation for which all roots of the characteristic equation
0127       // still have negative real parts. The most "inconvenient" pair
0128       // of roots there is (-0.270556 +- 2.50478 I) which leads
0129       // to oscillations with damping. The characteristic equation
0130       // of 5th and higher order ODEs are going to have roots with
0131       // positive real parts. Unless additional damping is
0132       // purposefully introduced into the system, numerical
0133       // solutions of such equations will just blow up.
0134       //
0135       assert(0);
0136   }
0137 }
0138
0139 void PadeTableODE::setParameters(const double* /* pars */, const unsigned nPars) { assert(nPars == 0U); }