HcalAlgos/interface/HcalPulseShapes.h

0001 #ifndef CALIBCALORIMETRY_HCALALGOS_HCALPULSESHAPES_H
0002 #define CALIBCALORIMETRY_HCALALGOS_HCALPULSESHAPES_H 1
0003
0004 #include <map>
0005 #include <vector>
0006 #include <cmath>
0007 #include "CalibCalorimetry/HcalAlgos/interface/HcalPulseShape.h"
0008 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0009 #include "FWCore/Framework/interface/Frameworkfwd.h"
0010 #include "FWCore/Utilities/interface/isFinite.h"
0011 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
0012 #include "CalibFormats/HcalObjects/interface/HcalDbRecord.h"
0013
0014 /** \class HcalPulseShapes
0015   *
0016   * \author J. Mans - Minnesota
0017   */
0018
0019 namespace CLHEP {
0020   class HepRandomEngine;
0021 }
0022
0023 class HcalPulseShapes {
0024 public:
0025   typedef HcalPulseShape Shape;
0026   // Default constructor is for callers that do not call beginRun(EventSetup)
0027   HcalPulseShapes();
0028   explicit HcalPulseShapes(edm::ConsumesCollector iC);
0029   ~HcalPulseShapes();
0030   // only needed if you'll be getting shapes by DetId
0031   void beginRun(edm::EventSetup const& es);
0032   void beginRun(const HcalDbService* conditions);
0033
0034   const Shape& hbShape() const { return hpdShape_; }
0035   const Shape& heShape() const { return hpdShape_; }
0036   const Shape& hfShape() const { return hfShape_; }
0037   const Shape& hoShape(bool sipm = false) const { return sipm ? siPMShapeHO_ : hpdShape_; }
0038   //  return Shape for given shapeType.
0039   const Shape& getShape(int shapeType) const;
0040   /// automatically figures out which shape to return
0041   const Shape& shape(const HcalDetId& detId) const;
0042   const Shape& shapeForReco(const HcalDetId& detId) const;
0043   /// in case of conditions problems
0044   const Shape& defaultShape(const HcalDetId& detId) const;
0045   //public static helpers
0046   static const int nBinsSiPM_ = 250;
0047   static constexpr float deltaTSiPM_ = 0.5;
0048   static constexpr float invDeltaTSiPM_ = 2.0;
0049   static double analyticPulseShapeSiPMHO(double t);
0050   static double analyticPulseShapeSiPMHE(double t);
0051   static constexpr float Y11RANGE_ = nBinsSiPM_;
0052   static constexpr float Y11MAX203_ = 0.04;
0053   static constexpr float Y11MAX206_ = 0.08;
0054   static double Y11203(double t);
0055   static double Y11206(double t);
0056   static double generatePhotonTime(CLHEP::HepRandomEngine* engine, unsigned int signalShape);
0057   static double generatePhotonTime203(CLHEP::HepRandomEngine* engine);
0058   static double generatePhotonTime206(CLHEP::HepRandomEngine* engine);
0059   //this function can take function pointers *or* functors!
0060   template <class F1, class F2>
0061   static std::vector<double> convolve(unsigned nbin, F1 f1, F2 f2) {
0062     std::vector<double> result(2 * nbin - 1, 0.);
0063     for (unsigned i = 0; i < 2 * nbin - 1; ++i) {
0064       for (unsigned j = 0; j < std::min(i + 1, nbin); ++j) {
0065         double tmp = f1(j) * f2(i - j);
0066         if (edm::isNotFinite(tmp))
0067           continue;
0068         result[i] += tmp;
0069       }
0070     }
0071     return result;
0072   }
0073   static std::vector<double> normalize(std::vector<double> nt, unsigned nbin) {
0074     //skip first bin, always 0
0075     double norm = 0.;
0076     for (unsigned int j = 1; j <= nbin; ++j) {
0077       norm += (nt[j] > 0) ? nt[j] : 0.;
0078     }
0079
0080     double normInv = 1. / norm;
0081     for (unsigned int j = 1; j <= nbin; ++j) {
0082       nt[j] *= normInv;
0083     }
0084
0085     return nt;
0086   }
0087   static std::vector<double> normalizeShift(std::vector<double> nt, unsigned nbin, int shift) {
0088     //skip first bin, always 0
0089     double norm = 0.;
0090     for (unsigned int j = std::max(1, -1 * shift); j <= nbin; j++) {
0091       norm += std::max(0., nt[j - shift]);
0092     }
0093     double normInv = 1. / norm;
0094     std::vector<double> nt2(nt.size(), 0);
0095     for (int j = 1; j <= (int)nbin; j++) {
0096       if (j - shift >= 0) {
0097         nt2[j] = nt[j - shift] * normInv;
0098       }
0099     }
0100     return nt2;
0101   }
0102
0103   std::map<int, Shape const*> const& get_all_shapes() const { return theShapes; }
0104
0105 private:
0106   void computeHPDShape(float, float, float, float, float, float, float, float, Shape&);
0107   void computeHFShape();
0108   void computeSiPMShapeHO();
0109   const HcalPulseShape& computeSiPMShapeHE203();
0110   const HcalPulseShape& computeSiPMShapeHE206();
0111   void computeSiPMShapeData2017();
0112   void computeSiPMShapeData2018();
0113   void computeSiPMShapeMCRecoRun3();
0114   Shape hpdShape_, hfShape_, siPMShapeHO_;
0115   Shape siPMShapeData2017_, siPMShapeData2018_, siPMShapeMCRecoRun3_;
0116   Shape hpdShape_v2, hpdShapeMC_v2;
0117   Shape hpdShape_v3, hpdShapeMC_v3;
0118   Shape hpdBV30Shape_v2, hpdBV30ShapeMC_v2;
0119   edm::ESGetToken<HcalDbService, HcalDbRecord> theDbServiceToken;
0120   const HcalDbService* theDbService;
0121   typedef std::map<int, const Shape*> ShapeMap;
0122   ShapeMap theShapes;
0123 };
0124 #endif