[HCAL] Mahi for HBHE local reconstruction

HLTrigger/Configuration/python/customizeHLTforCMSSW.py

@@ -17,6 +17,23 @@
 #                     pset.minGoodStripCharge = cms.PSet(refToPSet_ = cms.string('HLTSiStripClusterChargeCutNone'))
 #     return process
 
+# Add mahi to HCAL local reconstruction 
+def customiseFor21664(process):
+    from RecoLocalCalo.HcalRecProducers.HBHEPhase1Reconstructor_cfi import hbheprereco
+    for producer in producers_by_type(process, "HBHEPhase1Reconstructor"):
+        producer.algorithm.useM2     = hbheprereco.algorithm.useM2
+        producer.algorithm.useM3     = hbheprereco.algorithm.useM3
+        producer.algorithm.useMahi     = hbheprereco.algorithm.useMahi
+        producer.algorithm.dynamicPed   = hbheprereco.algorithm.dynamicPed
+        producer.algorithm.ts4Thresh    = hbheprereco.algorithm.ts4Thresh
+        producer.algorithm.chiSqSwitch  = hbheprereco.algorithm.chiSqSwitch
+        producer.algorithm.activeBXs    = hbheprereco.algorithm.activeBXs
+        producer.algorithm.nMaxItersMin = hbheprereco.algorithm.nMaxItersMin
+        producer.algorithm.nMaxItersNNLS = hbheprereco.algorithm.nMaxItersNNLS
+        producer.algorithm.deltaChiSqThresh = hbheprereco.algorithm.deltaChiSqThresh
+        producer.algorithm.nnlsThresh   = hbheprereco.algorithm.nnlsThresh
+    return process
+
 # Add new parameters to RecoTrackRefSelector
 def customiseFor19029(process):
     for producer in producers_by_type(process, "RecoTrackRefSelector"):
@@ -72,6 +89,8 @@ def customizeHLTforCMSSW(process, menuType="GRun"):
         process = customiseFor20269(process)
         process = customiseFor19989(process)
         process = customiseFor20429(process)
-        process = customiseFor21437(process)
+        process = customiseFor21437(process)        
+
+    process = customiseFor21664(process)
 
     return process

RecoLocalCalo/HcalRecAlgos/BuildFile.xml

@@ -1,4 +1,6 @@
 <use   name="boost"/>
+<use   name="eigen"/>
+<use   name="clhep"/>
 <use   name="DataFormats/HcalDigi"/>
 <use   name="DataFormats/HcalRecHit"/>
 <use   name="CalibFormats/HcalObjects"/>
@@ -9,6 +11,8 @@
 <use   name="FWCore/Framework"/>
 <use   name="FWCore/PluginManager"/>
 <use   name="FWCore/ParameterSet"/>
+<use   name="CondFormats/DataRecord"/>
+<use   name="vdt_headers"/>
 <use   name="rootminuit2"/> 
 <export>
   <lib   name="1"/>

RecoLocalCalo/HcalRecAlgos/interface/EigenMatrixTypes.h

@@ -0,0 +1,28 @@
+#ifndef RecoLocalCalo_HcalRecAlgos_EigenMatrixTypes_h
+#define RecoLocalCalo_HcalRecAlgos_EigenMatrixTypes_h
+
+#include <Eigen/Dense>
+
+constexpr int MaxSVSize=10;
+constexpr int MaxFSVSize=19;
+constexpr int MaxPVSize=10;
+
+typedef Eigen::Matrix<double,Eigen::Dynamic,1,0,MaxSVSize,1> SampleVector;
+typedef Eigen::Matrix<double,Eigen::Dynamic,1,0,MaxPVSize,1> PulseVector;
+typedef Eigen::Matrix<int,Eigen::Dynamic,1,0,MaxPVSize,1> BXVector;
+
+typedef Eigen::Matrix<double,MaxFSVSize,1> FullSampleVector;
+typedef Eigen::Matrix<double,MaxFSVSize,MaxFSVSize> FullSampleMatrix;
+
+typedef Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic,0,MaxSVSize,MaxSVSize> SampleMatrix;
+typedef Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic,0,MaxPVSize,MaxPVSize> PulseMatrix;
+typedef Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic,0,MaxSVSize,MaxPVSize> SamplePulseMatrix;
+
+typedef Eigen::LLT<SampleMatrix> SampleDecompLLT;
+typedef Eigen::LLT<PulseMatrix> PulseDecompLLT;
+typedef Eigen::LDLT<PulseMatrix> PulseDecompLDLT;
+
+typedef Eigen::Matrix<double,1,1> SingleMatrix;
+typedef Eigen::Matrix<double,1,1> SingleVector;
+
+#endif

RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h

@@ -0,0 +1,167 @@
+#ifndef RecoLocalCalo_HcalRecAlgos_MahiFit_HH
+#define RecoLocalCalo_HcalRecAlgos_MahiFit_HH
+
+#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"
+
+#include <Math/Functor.h>
+
+struct MahiNnlsWorkspace {
+
+  unsigned int nPulseTot;
+  unsigned int tsSize;
+  unsigned int tsOffset;
+  unsigned int fullTSOffset;
+  int bxOffset;
+  double dt;
+
+  //holds active bunch crossings
+  BXVector bxs;  
+
+  //holds data samples
+  SampleVector amplitudes;
+
+  //holds inverse covariance matrix
+  SampleMatrix invCovMat;
+
+  //holds diagonal noise terms
+  SampleVector noiseTerms;
+
+  //holds flat pedestal uncertainty
+  SampleMatrix pedConstraint;
+
+  //holds full covariance matrix for a pulse shape 
+  //varied in time
+  std::array<FullSampleMatrix, MaxPVSize> pulseCovArray;
+
+  //holds full pulse shape template
+  std::array<FullSampleVector, MaxPVSize> pulseShapeArray;
+
+  //holds full pulse shape derivatives
+  std::array<FullSampleVector, MaxPVSize> pulseDerivArray;
+
+  //holders for calculating pulse shape & covariance matrices
+  std::array<double, MaxSVSize> pulseN;
+  std::array<double, MaxSVSize> pulseM;
+  std::array<double, MaxSVSize> pulseP;
+
+  //holds matrix of pulse shape templates for each BX
+  SamplePulseMatrix pulseMat;
+
+  //holds matrix of pulse shape derivatives for each BX
+  SamplePulseMatrix pulseDerivMat;
+
+  //holds residual vector
+  PulseVector residuals;
+
+  //for FNNLS algorithm
+  unsigned int nP;
+  PulseVector ampVec;
+
+  PulseVector errVec;
+  PulseVector ampvecpermtest;
+
+  SamplePulseMatrix invcovp;
+  PulseMatrix aTaMat; // A-transpose A (matrix)
+  PulseVector aTbVec; // A-transpose b (vector)
+  PulseVector updateWork; // w (vector)
+
+  SampleDecompLLT covDecomp;
+  SampleMatrix covDecompLinv;
+  PulseMatrix topleft_work;
+  PulseDecompLDLT pulseDecomp;
+
+};
+
+class MahiFit
+{
+ public:
+  MahiFit();
+  ~MahiFit() { };
+
+  void setParameters(bool iDynamicPed, double iTS4Thresh, double chiSqSwitch, 
+		     bool iApplyTimeSlew, HcalTimeSlew::BiasSetting slewFlavor,
+		     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 doFit(std::array<float,3> &correctedOutput, int nbx) const;
+
+  void setPulseShapeTemplate  (const HcalPulseShapes::Shape& ps);
+  void resetPulseShapeTemplate(const HcalPulseShapes::Shape& ps);
+
+  typedef BXVector::Index Index;
+  const HcalPulseShapes::Shape* currentPulseShape_=nullptr;
+
+ private:
+
+  double minimize() const;
+  void onePulseMinimize() const;
+  void updateCov() const;
+  void updatePulseShape(double itQ, FullSampleVector &pulseShape, 
+			FullSampleVector &pulseDeriv,
+			FullSampleMatrix &pulseCov) const;
+  double calculateArrivalTime() const;
+  double 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;
+
+  double getSiPMDarkCurrent(double darkCurrent, double fcByPE, double lambda) const;
+
+  mutable MahiNnlsWorkspace nnlsWork_;
+
+  //hard coded in initializer
+  const unsigned int fullTSSize_;
+  const unsigned int fullTSofInterest_;
+
+  static constexpr int pedestalBX_ = 100;
+
+  // Python-configurables
+  bool dynamicPed_;
+  float ts4Thresh_; 
+  float chiSqSwitch_; 
+
+  bool applyTimeSlew_; 
+  HcalTimeSlew::BiasSetting slewFlavor_;
+
+  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 cntsetPulseShape_;
+  std::unique_ptr<FitterFuncs::PulseShapeFunctor> psfPtr_;
+  std::unique_ptr<ROOT::Math::Functor> pfunctor_;
+
+}; 
+#endif

RecoLocalCalo/HcalRecAlgos/interface/PulseShapeFunctor.h

@@ -39,6 +39,10 @@ namespace FitterFuncs{
      double singlePulseShapeFunc( const double *x );
      double doublePulseShapeFunc( const double *x );
      double triplePulseShapeFunc( const double *x );
+
+     void getPulseShape(std::array<double,HcalConst::maxSamples>& fillPulseShape) { 
+       fillPulseShape = pulse_shape_;
+     }
 
    private:
      std::array<float,HcalConst::maxPSshapeBin> pulse_hist;

RecoLocalCalo/HcalRecAlgos/interface/SimpleHBHEPhase1Algo.h

@@ -12,7 +12,7 @@
 
 #include "RecoLocalCalo/HcalRecAlgos/interface/PulseShapeFitOOTPileupCorrection.h"
 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalDeterministicFit.h"
-
+#include "RecoLocalCalo/HcalRecAlgos/interface/MahiFit.h"
 
 class SimpleHBHEPhase1Algo : public AbsHBHEPhase1Algo
 {
@@ -45,7 +45,8 @@ class SimpleHBHEPhase1Algo : public AbsHBHEPhase1Algo
                          float timeShift,
                          bool correctForPhaseContainment,
                          std::unique_ptr<PulseShapeFitOOTPileupCorrection> m2,
-                         std::unique_ptr<HcalDeterministicFit> detFit);
+                         std::unique_ptr<HcalDeterministicFit> detFit,
+			 std::unique_ptr<MahiFit> mahi);
 
     inline ~SimpleHBHEPhase1Algo() override {}
 
@@ -101,6 +102,9 @@ class SimpleHBHEPhase1Algo : public AbsHBHEPhase1Algo
     // "Metod 3" algorithm
     std::unique_ptr<HcalDeterministicFit> hltOOTpuCorr_;
 
+    // Mahi algorithm
+    std::unique_ptr<MahiFit> mahiOOTpuCorr_;
+
     HcalPulseShapes theHcalPulseShapes_;
 };