HcalRecAlgos/interface/MahiFit.h

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   //holds active bunch crossings
0028   BXVector bxs;
0029
0030   //holds data samples
0031   SampleVector amplitudes;
0032
0033   //holds diagonal noise terms
0034   SampleVector noiseTerms;
0035
0036   //holds diagonal pedestal noise terms
0037   SampleVector pedVals;
0038
0039   //holds flat pedestal uncertainty
0040   float pedVal;
0041
0042   float noisecorr;
0043
0044   //holds full covariance matrix for a pulse shape
0045   //varied in time
0046   std::array<SampleMatrix, MaxPVSize> pulseCovArray;
0047
0048   //holds matrix of pulse shape templates for each BX
0049   SamplePulseMatrix pulseMat;
0050
0051   //holds matrix of pulse shape derivatives for each BX
0052   SamplePulseMatrix pulseDerivMat;
0053
0054   //for FNNLS algorithm
0055   unsigned int nP;
0056   PulseVector ampVec;
0057
0058   SamplePulseMatrix invcovp;
0059   PulseMatrix aTaMat;  // A-transpose A (matrix)
0060   PulseVector aTbVec;  // A-transpose b (vector)
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                      int iTimeAlgo,
0108                      double iThEnergeticPulses,
0109                      double iMeanTime,
0110                      double iTimeSigmaHPD,
0111                      double iTimeSigmaSiPM,
0112                      const std::vector<int>& iActiveBXs,
0113                      int iNMaxItersMin,
0114                      int iNMaxItersNNLS,
0115                      double iDeltaChiSqThresh,
0116                      double iNnlsThresh);
0117
0118   void phase1Apply(const HBHEChannelInfo& channelData,
0119                    float& reconstructedEnergy,
0120                    float& soiPlusOneEnergy,
0121                    float& reconstructedTime,
0122                    bool& useTriple,
0123                    float& chi2) const;
0124
0125   void phase1Debug(const HBHEChannelInfo& channelData, MahiDebugInfo& mdi) const;
0126
0127   void doFit(std::array<float, 4>& correctedOutput, const int nbx) const;
0128
0129   void setPulseShapeTemplate(int pulseShapeId,
0130                              const HcalPulseShapes& ps,
0131                              bool hasTimeInfo,
0132                              const HcalTimeSlew* hcalTimeSlewDelay,
0133                              unsigned int nSamples,
0134                              const float gain);
0135
0136   typedef BXVector::Index Index;
0137   const HcalTimeSlew* hcalTimeSlewDelay_ = nullptr;
0138
0139   float thEnergeticPulses_;
0140   float thEnergeticPulsesFC_;
0141
0142 private:
0143   typedef std::pair<int, std::shared_ptr<FitterFuncs::PulseShapeFunctor> > ShapeWithId;
0144
0145   const float minimize() const;
0146   void onePulseMinimize() const;
0147   void updateCov(const SampleMatrix& invCovMat) const;
0148   void resetPulseShapeTemplate(int pulseShapeId, const HcalPulseShapes& ps, unsigned int nSamples);
0149
0150   float ccTime(const float itQ) const;
0151   void updatePulseShape(const float itQ,
0152                         FullSampleVector& pulseShape,
0153                         FullSampleVector& pulseDeriv,
0154                         FullSampleMatrix& pulseCov) const;
0155
0156   float calculateArrivalTime(const unsigned int iBX) const;
0157   float calculateChiSq() const;
0158   void nnls() const;
0159   void resetWorkspace() const;
0160
0161   void nnlsUnconstrainParameter(Index idxp) const;
0162   void nnlsConstrainParameter(Index minratioidx) const;
0163
0164   void solveSubmatrix(PulseMatrix& mat, PulseVector& invec, PulseVector& outvec, unsigned nP) const;
0165
0166   mutable MahiNnlsWorkspace nnlsWork_;
0167
0168   //hard coded in initializer
0169   static constexpr int pedestalBX_ = 100;
0170
0171   // used to restrict returned time value to a 25 ns window centered
0172   // on the nominal arrival time
0173   static constexpr float timeLimit_ = 12.5f;
0174
0175   // Python-configurables
0176   int timeAlgo_;
0177
0178   bool dynamicPed_;
0179   float ts4Thresh_;
0180   float chiSqSwitch_;
0181
0182   bool applyTimeSlew_;
0183   HcalTimeSlew::BiasSetting slewFlavor_;
0184   float tsDelay1GeV_ = 0.f;
0185   float norm_ = (1.f / std::sqrt(12));
0186
0187   bool calculateArrivalTime_;
0188   float meanTime_;
0189   float timeSigmaHPD_;
0190   float timeSigmaSiPM_;
0191
0192   std::vector<int> activeBXs_;
0193
0194   int nMaxItersMin_;
0195   int nMaxItersNNLS_;
0196
0197   float deltaChiSqThresh_;
0198   float nnlsThresh_;
0199
0200   unsigned int bxSizeConf_;
0201   int bxOffsetConf_;
0202
0203   //for pulse shapes
0204   int currentPulseShapeId_ = INT_MIN;
0205   int cntsetPulseShape_ = 0;
0206   FitterFuncs::PulseShapeFunctor* psfPtr_ = nullptr;
0207   std::vector<ShapeWithId> knownPulseShapes_;
0208 };
0209
0210 #endif