File indexing completed on 2023-03-17 10:41:45
0001 #ifndef CalibCalorimetry_HcalAlgos_ConstantStepOdeSolver_h_
0002 #define CalibCalorimetry_HcalAlgos_ConstantStepOdeSolver_h_
0003
0004 #include <vector>
0005 #include <iostream>
0006 #include "FWCore/Utilities/interface/Exception.h"
0007
0008 #include "CalibCalorimetry/HcalAlgos/interface/AbsODERHS.h"
0009
0010
0011
0012
0013
0014 class ConstantStepOdeSolver {
0015 public:
0016 inline ConstantStepOdeSolver() : rhs_(nullptr), dt_(0.0), dim_(0), runLen_(0), lastIntegrated_(0) {}
0017
0018 inline ConstantStepOdeSolver(const AbsODERHS& rhs)
0019 : rhs_(nullptr), dt_(0.0), dim_(0), runLen_(0), lastIntegrated_(0) {
0020 rhs_ = rhs.clone();
0021 }
0022
0023
0024 ConstantStepOdeSolver(const ConstantStepOdeSolver& r);
0025 ConstantStepOdeSolver& operator=(const ConstantStepOdeSolver& r);
0026
0027 inline virtual ~ConstantStepOdeSolver() { delete rhs_; }
0028
0029
0030 inline void setRHS(const AbsODERHS& rhs) {
0031 delete rhs_;
0032 rhs_ = rhs.clone();
0033 }
0034 inline const AbsODERHS* getRHS() const { return rhs_; }
0035 inline AbsODERHS* getRHS() { return rhs_; }
0036
0037
0038 inline unsigned lastDim() const { return dim_; }
0039 inline unsigned lastRunLength() const { return runLen_; }
0040 inline double lastDeltaT() const { return dt_; }
0041 inline double lastMaxT() const { return runLen_ ? dt_ * (runLen_ - 1U) : 0.0; }
0042
0043 inline double getTime(const unsigned idx) const {
0044 if (idx >= runLen_)
0045 throw cms::Exception("In ConstantStepOdeSolver::getTime: index out of range");
0046 return idx * dt_;
0047 }
0048
0049 inline double getCoordinate(const unsigned which, const unsigned idx) const {
0050 if (which >= dim_ || idx >= runLen_)
0051 throw cms::Exception("In ConstantStepOdeSolver::getCoordinate: index out of range");
0052 return historyBuffer_[dim_ * idx + which];
0053 }
0054
0055
0056
0057 double getIntegrated(unsigned which, unsigned idx) const;
0058
0059
0060 void truncateCoordinate(unsigned which, double minValue, double maxValue);
0061
0062
0063
0064 double interpolateCoordinate(unsigned which, double t, bool cubic = false) const;
0065
0066
0067 double interpolateIntegrated(unsigned which, double t, bool cubic = false) const;
0068
0069
0070 double getPeakTime(unsigned which) const;
0071
0072
0073 void run(const double* initialConditions, unsigned lenConditions, double dt, unsigned nSteps);
0074
0075
0076
0077 void setHistory(double dt, const double* data, unsigned dim, unsigned runLen);
0078
0079
0080 void writeHistory(std::ostream& os, double dt, bool cubic = false) const;
0081
0082
0083 void writeIntegrated(std::ostream& os, unsigned which, double dt, bool cubic = false) const;
0084
0085
0086 virtual const char* methodName() const = 0;
0087
0088 protected:
0089 AbsODERHS* rhs_;
0090
0091 private:
0092
0093 virtual void step(double t, double dt, const double* x, unsigned lenX, double* coordIncrement) const = 0;
0094
0095
0096
0097 void integrateCoordinate(const unsigned which);
0098
0099 double dt_;
0100 unsigned dim_;
0101 unsigned runLen_;
0102 unsigned lastIntegrated_;
0103
0104 std::vector<double> historyBuffer_;
0105 std::vector<double> chargeBuffer_;
0106 };
0107
0108
0109 inline std::ostream& operator<<(std::ostream& os, const ConstantStepOdeSolver& s) {
0110 s.writeHistory(os, s.lastDeltaT());
0111 return os;
0112 }
0113
0114
0115 class EulerOdeSolver : public ConstantStepOdeSolver {
0116 public:
0117 inline EulerOdeSolver() : ConstantStepOdeSolver() {}
0118
0119 inline explicit EulerOdeSolver(const AbsODERHS& rhs) : ConstantStepOdeSolver(rhs) {}
0120
0121 inline const char* methodName() const override { return "Euler"; }
0122
0123 private:
0124 void step(double t, double dt, const double* x, unsigned lenX, double* coordIncrement) const override;
0125 };
0126
0127 class RK2 : public ConstantStepOdeSolver {
0128 public:
0129 inline RK2() : ConstantStepOdeSolver() {}
0130
0131 inline explicit RK2(const AbsODERHS& rhs) : ConstantStepOdeSolver(rhs) {}
0132
0133 inline const char* methodName() const override { return "2nd order Runge-Kutta"; }
0134
0135 private:
0136 void step(double t, double dt, const double* x, unsigned lenX, double* coordIncrement) const override;
0137
0138 mutable std::vector<double> buf_;
0139 };
0140
0141 class RK4 : public ConstantStepOdeSolver {
0142 public:
0143 inline RK4() : ConstantStepOdeSolver() {}
0144
0145 inline explicit RK4(const AbsODERHS& rhs) : ConstantStepOdeSolver(rhs) {}
0146
0147 inline const char* methodName() const override { return "4th order Runge-Kutta"; }
0148
0149 private:
0150 void step(double t, double dt, const double* x, unsigned lenX, double* coordIncrement) const override;
0151
0152 mutable std::vector<double> buf_;
0153 };
0154
0155 #endif