HGCalRecAlgos/interface/HGCalUncalibRecHitRecWeightsAlgo.h

0001 #ifndef RecoLocalCalo_HGCalRecAlgos_HGCalUncalibRecHitRecWeightsAlgo_HH
0002 #define RecoLocalCalo_HGCalRecAlgos_HGCalUncalibRecHitRecWeightsAlgo_HH
0003
0004 /** \class HGalUncalibRecHitRecWeightsAlgo
0005   *  compute amplitude, pedestal, time jitter, chi2 of a pulse
0006   *  using a weights method, a la Ecal
0007   *
0008   *  \author Valeri Andreev
0009   *
0010   *
0011   */
0012
0013 #include "RecoLocalCalo/HGCalRecAlgos/interface/HGCalUncalibRecHitRecAbsAlgo.h"
0014 #include <vector>
0015 #include <cmath>
0016
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0019
0020 template <class C>
0021 class HGCalUncalibRecHitRecWeightsAlgo {
0022 public:
0023   // destructor
0024   virtual ~HGCalUncalibRecHitRecWeightsAlgo() {}
0025
0026   void set_isSiFESim(const bool isSiFE) { isSiFESim_ = isSiFE; }
0027   bool isSiFESim() const { return isSiFESim_; }
0028
0029   void set_ADCLSB(const double adclsb) { adcLSB_ = adclsb; }
0030   void set_TDCLSB(const double tdclsb) { tdcLSB_ = tdclsb; }
0031
0032   void set_toaLSBToNS(const double lsb2ns) { toaLSBToNS_ = lsb2ns; }
0033
0034   void set_tdcOnsetfC(const double tdcOnset) { tdcOnsetfC_ = tdcOnset; }
0035
0036   void set_tofDelay(const double tofDelay) { tofDelay_ = tofDelay; }
0037
0038   void set_fCPerMIP(const std::vector<double>& fCPerMIP) {
0039     if (std::any_of(fCPerMIP.cbegin(), fCPerMIP.cend(), [](double conv) { return conv <= 0.0; })) {
0040       throw cms::Exception("BadConversionFactor") << "At least one of fCPerMIP is zero!" << std::endl;
0041     }
0042     fCPerMIP_ = fCPerMIP;
0043   }
0044
0045   bool setGeometry(const edm::ESHandle<HGCalGeometry>& geom) {
0046     geom_ = geom.isValid() ? geom.product() : nullptr;
0047     if (isSiFESim_ && geom_)
0048       ddd_ = &(geom_->topology().dddConstants());
0049     else
0050       ddd_ = nullptr;
0051     return geom_ != nullptr;
0052   }
0053
0054   /// Compute HGCUncalibratedRecHit from DataFrame
0055   virtual HGCUncalibratedRecHit makeRecHit(const C& dataFrame, const bool computeLocalTime) {
0056     double amplitude_(-1.), pedestal_(-1.), jitter_(-99.), chi2_(-1.);
0057     uint32_t flag = 0;
0058
0059     constexpr int iSample = 2;  //only in-time sample
0060     const auto& sample = dataFrame.sample(iSample);
0061
0062     // Were digis done w/ the complete digitization and using signal shape?
0063     // (originally done only for the silicon, while for scitillator it was trivial. Fomr 11_ also scinti uses shape)
0064     if (isSiFESim_) {
0065       // mode == true: TDC readout was activated and amplitude comes from TimeOverThreshold
0066       if (sample.mode()) {
0067         flag = !sample.threshold();  // raise flag if busy cell
0068         // LG (23/06/2015):
0069         // to get a continuous energy spectrum we must add here the maximum value in fC ever
0070         // reported by the ADC. Namely: floor(tdcOnset/adcLSB_) * adcLSB_
0071         // need to increment by one so TDC doesn't overlap with ADC last bin
0072         // LG (11/04/2016):
0073         // offset the TDC upwards to reflect the bin center
0074         amplitude_ = (std::floor(tdcOnsetfC_ / adcLSB_) + 1.0) * adcLSB_ + (double(sample.data()) + 0.5) * tdcLSB_;
0075       } else {
0076         amplitude_ = double(sample.data()) * adcLSB_;  // why do we not have +0.5 here ?
0077       }  //isSiFESim_
0078     }  //mode()
0079
0080     // trivial digitization, i.e. no signal shape
0081     else {
0082       amplitude_ = double(sample.data()) * adcLSB_;
0083     }
0084
0085     if (sample.getToAValid()) {
0086       const auto& dist2center = geom_ ? geom_->getPosition(dataFrame.id()).mag() : 0;
0087       jitter_ = computeLocalTime ? double(sample.toa()) * toaLSBToNS_ - tofDelay_
0088                                  : double(sample.toa()) * toaLSBToNS_ - dist2center / c_cm_ns - tofDelay_;
0089     }
0090
0091     int thickness = (ddd_ != nullptr) ? ddd_->waferType(dataFrame.id(), false) : 0;
0092     amplitude_ = amplitude_ / fCPerMIP_[thickness];
0093
0094     LogDebug("HGCUncalibratedRecHit") << "isSiFESim_: " << isSiFESim_ << " ADC+: set the charge to: " << amplitude_
0095                                       << ' ' << sample.data() << ' ' << adcLSB_ << ' '
0096                                       << "               TDC+: set the ToA to: " << jitter_ << ' ' << sample.toa()
0097                                       << ' ' << toaLSBToNS_ << ' ' << tdcLSB_
0098                                       << "               getToAValid(): " << sample.getToAValid()
0099                                       << " mode(): " << sample.mode() << std::endl;
0100     LogDebug("HGCUncalibratedRecHit") << "Final uncalibrated amplitude : " << amplitude_ << std::endl;
0101
0102     return HGCUncalibratedRecHit(dataFrame.id(), amplitude_, pedestal_, jitter_, chi2_, flag);
0103   }
0104
0105 private:
0106   static constexpr float c_cm_ns = CLHEP::c_light * CLHEP::ns / CLHEP::cm;
0107   double adcLSB_, tdcLSB_, toaLSBToNS_, tdcOnsetfC_, tofDelay_;
0108   bool isSiFESim_;
0109   std::vector<double> fCPerMIP_;
0110   const HGCalDDDConstants* ddd_;
0111   const HGCalGeometry* geom_;
0112 };
0113 #endif