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                      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   //hard coded in initializer
0161   static constexpr int pedestalBX_ = 100;
0162
0163   // used to restrict returned time value to a 25 ns window centered
0164   // on the nominal arrival time
0165   static constexpr float timeLimit_ = 12.5f;
0166
0167   // Python-configurables
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   //for pulse shapes
0194   int currentPulseShapeId_ = INT_MIN;
0195   int cntsetPulseShape_ = 0;
0196   FitterFuncs::PulseShapeFunctor* psfPtr_ = nullptr;
0197   std::vector<ShapeWithId> knownPulseShapes_;
0198 };
0199
0200 #endif