File indexing completed on 2021-11-28 23:48:59
0001 #ifndef RecoLocalCalo_HcalRecAlgos_MahiFit_HH
0002 #define RecoLocalCalo_HcalRecAlgos_MahiFit_HH
0003
0004 #include <climits>
0005 #include <utility>
0006 #include <memory>
0007
0008 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0009 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0010 #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
0011 #include "CalibFormats/HcalObjects/interface/HcalCoder.h"
0012 #include "RecoLocalCalo/HcalRecAlgos/interface/EigenMatrixTypes.h"
0013 #include "DataFormats/HcalRecHit/interface/HBHEChannelInfo.h"
0014
0015 #include "CalibCalorimetry/HcalAlgos/interface/HcalPulseShapes.h"
0016 #include "CalibCalorimetry/HcalAlgos/interface/HcalTimeSlew.h"
0017 #include "RecoLocalCalo/HcalRecAlgos/interface/PulseShapeFunctor.h"
0018
0019 struct MahiNnlsWorkspace {
0020 unsigned int nPulseTot;
0021 unsigned int tsSize = 0U;
0022 unsigned int tsOffset;
0023 int bxOffset;
0024 int maxoffset;
0025 float dt;
0026
0027
0028 BXVector bxs;
0029
0030
0031 SampleVector amplitudes;
0032
0033
0034 SampleVector noiseTerms;
0035
0036
0037 SampleVector pedVals;
0038
0039
0040 float pedVal;
0041
0042 float noisecorr;
0043
0044
0045
0046 std::array<SampleMatrix, MaxPVSize> pulseCovArray;
0047
0048
0049 SamplePulseMatrix pulseMat;
0050
0051
0052 SamplePulseMatrix pulseDerivMat;
0053
0054
0055 unsigned int nP;
0056 PulseVector ampVec;
0057
0058 SamplePulseMatrix invcovp;
0059 PulseMatrix aTaMat;
0060 PulseVector aTbVec;
0061
0062 SampleDecompLLT covDecomp;
0063 };
0064
0065 struct MahiDebugInfo {
0066 int nSamples;
0067 int soi;
0068
0069 bool use3;
0070
0071 float inTimeConst;
0072 float inDarkCurrent;
0073 float inPedAvg;
0074 float inGain;
0075
0076 float inNoiseADC[MaxSVSize];
0077 float inNoiseDC[MaxSVSize];
0078 float inNoisePhoto[MaxSVSize];
0079 float inPedestal[MaxSVSize];
0080
0081 float totalUCNoise[MaxSVSize];
0082
0083 float mahiEnergy;
0084 float chiSq;
0085 float arrivalTime;
0086 float pedEnergy;
0087 float ootEnergy[7];
0088
0089 float count[MaxSVSize];
0090 float inputTS[MaxSVSize];
0091 int inputTDC[MaxSVSize];
0092 float itPulse[MaxSVSize];
0093 float ootPulse[7][MaxSVSize];
0094 };
0095
0096 class MahiFit {
0097 public:
0098 MahiFit();
0099 ~MahiFit(){};
0100
0101 void setParameters(bool iDynamicPed,
0102 double iTS4Thresh,
0103 double chiSqSwitch,
0104 bool iApplyTimeSlew,
0105 HcalTimeSlew::BiasSetting slewFlavor,
0106 bool iCalculateArrivalTime,
0107 double iMeanTime,
0108 double iTimeSigmaHPD,
0109 double iTimeSigmaSiPM,
0110 const std::vector<int>& iActiveBXs,
0111 int iNMaxItersMin,
0112 int iNMaxItersNNLS,
0113 double iDeltaChiSqThresh,
0114 double iNnlsThresh);
0115
0116 void phase1Apply(const HBHEChannelInfo& channelData,
0117 float& reconstructedEnergy,
0118 float& soiPlusOneEnergy,
0119 float& reconstructedTime,
0120 bool& useTriple,
0121 float& chi2) const;
0122
0123 void phase1Debug(const HBHEChannelInfo& channelData, MahiDebugInfo& mdi) const;
0124
0125 void doFit(std::array<float, 4>& correctedOutput, const int nbx) const;
0126
0127 void setPulseShapeTemplate(int pulseShapeId,
0128 const HcalPulseShapes& ps,
0129 bool hasTimeInfo,
0130 const HcalTimeSlew* hcalTimeSlewDelay,
0131 unsigned int nSamples);
0132
0133 typedef BXVector::Index Index;
0134 const HcalTimeSlew* hcalTimeSlewDelay_ = nullptr;
0135
0136 private:
0137 typedef std::pair<int, std::shared_ptr<FitterFuncs::PulseShapeFunctor> > ShapeWithId;
0138
0139 const float minimize() const;
0140 void onePulseMinimize() const;
0141 void updateCov(const SampleMatrix& invCovMat) const;
0142 void resetPulseShapeTemplate(int pulseShapeId, const HcalPulseShapes& ps, unsigned int nSamples);
0143 void updatePulseShape(const float itQ,
0144 FullSampleVector& pulseShape,
0145 FullSampleVector& pulseDeriv,
0146 FullSampleMatrix& pulseCov) const;
0147
0148 float calculateArrivalTime(const unsigned int iBX) const;
0149 float calculateChiSq() const;
0150 void nnls() const;
0151 void resetWorkspace() const;
0152
0153 void nnlsUnconstrainParameter(Index idxp) const;
0154 void nnlsConstrainParameter(Index minratioidx) const;
0155
0156 void solveSubmatrix(PulseMatrix& mat, PulseVector& invec, PulseVector& outvec, unsigned nP) const;
0157
0158 mutable MahiNnlsWorkspace nnlsWork_;
0159
0160
0161 static constexpr int pedestalBX_ = 100;
0162
0163
0164
0165 static constexpr float timeLimit_ = 12.5f;
0166
0167
0168 bool dynamicPed_;
0169 float ts4Thresh_;
0170 float chiSqSwitch_;
0171
0172 bool applyTimeSlew_;
0173 HcalTimeSlew::BiasSetting slewFlavor_;
0174 float tsDelay1GeV_ = 0.f;
0175 float norm_ = (1.f / std::sqrt(12));
0176
0177 bool calculateArrivalTime_;
0178 float meanTime_;
0179 float timeSigmaHPD_;
0180 float timeSigmaSiPM_;
0181
0182 std::vector<int> activeBXs_;
0183
0184 int nMaxItersMin_;
0185 int nMaxItersNNLS_;
0186
0187 float deltaChiSqThresh_;
0188 float nnlsThresh_;
0189
0190 unsigned int bxSizeConf_;
0191 int bxOffsetConf_;
0192
0193
0194 int currentPulseShapeId_ = INT_MIN;
0195 int cntsetPulseShape_ = 0;
0196 FitterFuncs::PulseShapeFunctor* psfPtr_ = nullptr;
0197 std::vector<ShapeWithId> knownPulseShapes_;
0198 };
0199
0200 #endif