MahiFit.h

Go to the documentation of this file.

#ifndef RecoLocalCalo_HcalRecAlgos_MahiFit_HH
#define RecoLocalCalo_HcalRecAlgos_MahiFit_HH
 
#include <climits>
#include <utility>
#include <memory>
 
#include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
#include "CalibFormats/HcalObjects/interface/HcalCoder.h"
#include "RecoLocalCalo/HcalRecAlgos/interface/EigenMatrixTypes.h"
#include "DataFormats/HcalRecHit/interface/HBHEChannelInfo.h"
 
#include "CalibCalorimetry/HcalAlgos/interface/HcalPulseShapes.h"
#include "CalibCalorimetry/HcalAlgos/interface/HcalTimeSlew.h"
#include "RecoLocalCalo/HcalRecAlgos/interface/PulseShapeFunctor.h"
 
struct MahiNnlsWorkspace {
  unsigned int nPulseTot;
  unsigned int tsSize = 0U;
  unsigned int tsOffset;
  int bxOffset;
  int maxoffset;
  float dt;
 
  //holds active bunch crossings
  BXVector bxs;
 
  //holds data samples
  SampleVector amplitudes;
 
  //holds diagonal noise terms
  SampleVector noiseTerms;
 
  //holds diagonal pedestal noise terms
  SampleVector pedVals;
 
  //holds flat pedestal uncertainty
  float pedVal;
 
  float noisecorr;
 
  //holds full covariance matrix for a pulse shape
  //varied in time
  std::array<SampleMatrix, MaxPVSize> pulseCovArray;
 
  //holds matrix of pulse shape templates for each BX
  SamplePulseMatrix pulseMat;
 
  //holds matrix of pulse shape derivatives for each BX
  SamplePulseMatrix pulseDerivMat;
 
  //for FNNLS algorithm
  unsigned int nP;
  PulseVector ampVec;
 
  SamplePulseMatrix invcovp;
  PulseMatrix aTaMat;  // A-transpose A (matrix)
  PulseVector aTbVec;  // A-transpose b (vector)
 
  SampleDecompLLT covDecomp;
};
 
struct MahiDebugInfo {
  int nSamples;
  int soi;
 
  bool use3;
 
  float inTimeConst;
  float inDarkCurrent;
  float inPedAvg;
  float inGain;
 
  float inNoiseADC[MaxSVSize];
  float inNoiseDC[MaxSVSize];
  float inNoisePhoto[MaxSVSize];
  float inPedestal[MaxSVSize];
 
  float totalUCNoise[MaxSVSize];
 
  float mahiEnergy;
  float chiSq;
  float arrivalTime;
  float pedEnergy;
  float ootEnergy[7];
 
  float count[MaxSVSize];
  float inputTS[MaxSVSize];
  int inputTDC[MaxSVSize];
  float itPulse[MaxSVSize];
  float ootPulse[7][MaxSVSize];
};
 
class MahiFit {
public:
  MahiFit();
  ~MahiFit(){};
 
  void setParameters(bool iDynamicPed,
                     double iTS4Thresh,
                     double chiSqSwitch,
                     bool iApplyTimeSlew,
                     HcalTimeSlew::BiasSetting slewFlavor,
                     bool iCalculateArrivalTime,
                     double iMeanTime,
                     double iTimeSigmaHPD,
                     double iTimeSigmaSiPM,
                     const std::vector<int>& iActiveBXs,
                     int iNMaxItersMin,
                     int iNMaxItersNNLS,
                     double iDeltaChiSqThresh,
                     double iNnlsThresh);
 
  void phase1Apply(const HBHEChannelInfo& channelData,
                   float& reconstructedEnergy,
                   float& reconstructedTime,
                   bool& useTriple,
                   float& chi2) const;
 
  void phase1Debug(const HBHEChannelInfo& channelData, MahiDebugInfo& mdi) const;
 
  void doFit(std::array<float, 3>& correctedOutput, const int nbx) const;
 
  void setPulseShapeTemplate(int pulseShapeId,
                             const HcalPulseShapes& ps,
                             bool hasTimeInfo,
                             const HcalTimeSlew* hcalTimeSlewDelay,
                             unsigned int nSamples);
 
  typedef BXVector::Index Index;
  const HcalTimeSlew* hcalTimeSlewDelay_ = nullptr;
 
private:
  typedef std::pair<int, std::shared_ptr<FitterFuncs::PulseShapeFunctor> > ShapeWithId;
 
  const float minimize() const;
  void onePulseMinimize() const;
  void updateCov(const SampleMatrix& invCovMat) const;
  void resetPulseShapeTemplate(int pulseShapeId, const HcalPulseShapes& ps, unsigned int nSamples);
  void updatePulseShape(const float itQ,
                        FullSampleVector& pulseShape,
                        FullSampleVector& pulseDeriv,
                        FullSampleMatrix& pulseCov) const;
 
  float calculateArrivalTime(const unsigned int iBX) const;
  float calculateChiSq() const;
  void nnls() const;
  void resetWorkspace() const;
 
  void nnlsUnconstrainParameter(Index idxp) const;
  void nnlsConstrainParameter(Index minratioidx) const;
 
  void solveSubmatrix(PulseMatrix& mat, PulseVector& invec, PulseVector& outvec, unsigned nP) const;
 
  mutable MahiNnlsWorkspace nnlsWork_;
 
  //hard coded in initializer
  static constexpr int pedestalBX_ = 100;
 
  // used to restrict returned time value to a 25 ns window centered
  // on the nominal arrival time
  static constexpr float timeLimit_ = 12.5f;
 
  // Python-configurables
  bool dynamicPed_;
  float ts4Thresh_;
  float chiSqSwitch_;
 
  bool applyTimeSlew_;
  HcalTimeSlew::BiasSetting slewFlavor_;
  float tsDelay1GeV_ = 0.f;
  float norm_ = (1.f / std::sqrt(12));
 
  bool calculateArrivalTime_;
  float meanTime_;
  float timeSigmaHPD_;
  float timeSigmaSiPM_;
 
  std::vector<int> activeBXs_;
 
  int nMaxItersMin_;
  int nMaxItersNNLS_;
 
  float deltaChiSqThresh_;
  float nnlsThresh_;
 
  unsigned int bxSizeConf_;
  int bxOffsetConf_;
 
  //for pulse shapes
  int currentPulseShapeId_ = INT_MIN;
  int cntsetPulseShape_ = 0;
  FitterFuncs::PulseShapeFunctor* psfPtr_ = nullptr;
  std::vector<ShapeWithId> knownPulseShapes_;
};
 
#endif