HgcrocEmulator.h

 
#ifndef TOOLS_HGCROCEMULATOR_H
#define TOOLS_HGCROCEMULATOR_H
 
#include "Conditions/SimpleTableCondition.h"
#include "Framework/Configure/Parameters.h"
#include "Recon/Event/HgcrocDigiCollection.h"
#include "SimCore/Event/SimCalorimeterHit.h"
#include "Tools/NoiseGenerator.h"
 
//----------//
//   ROOT   //
//----------//
#include "TF1.h"
#include "TRandom3.h"
 
namespace ldmx {
 
class HgcrocEmulator {
 public:
  HgcrocEmulator(const framework::config::Parameters& ps);
 
  ~HgcrocEmulator() {}
 
  bool hasSeed() const { return noiseInjector_.get() != nullptr; }
 
  void seedGenerator(uint64_t seed);
 
  void condition(const conditions::DoubleTableCondition& table) {
    // reset cache of column numbers if table changes
    if (&table != chipConditions_) conditionNamesToIndex_.clear();
    chipConditions_ = &table;
  }
 
  bool digitize(
      const int& channelID,
      std::vector<std::pair<double, double>>& arriving_pulses,
      std::vector<ldmx::HgcrocDigiCollection::Sample>& digiToAdd) const;
 
  std::vector<ldmx::HgcrocDigiCollection::Sample> noiseDigi(
      const int& channel, const double& soi_amplitude = 0) const;
 
  double noise(const int& channelID) const {
    return noiseInjector_->Gaus(
        0, getCondition(channelID, "NOISE") * gain(channelID));
  };
 
  double gain(const int& channelID) const {
    return getCondition(channelID, "GAIN");
  }
 
  double pedestal(const int& id) const { return getCondition(id, "PEDESTAL"); }
 
  double readoutThreshold(const int& id) const {
    return getCondition(id, "READOUT_THRESHOLD");
  }
 
 private:
  double getCondition(int id, const std::string& name) const {
    // check if emulator has been passed a table of conditions
    if (!chipConditions_) {
      EXCEPTION_RAISE("HgcrocCond",
                      "HGC ROC Emulator was not given a conditions table.");
    }
 
    // cache column index for the input name
    if (conditionNamesToIndex_.count(name) == 0)
      conditionNamesToIndex_[name] = chipConditions_->getColumnNumber(name);
 
    // get condition
    return chipConditions_->get(id, conditionNamesToIndex_.at(name));
  }
 
 private:
  class CompositePulse {
   public:
    CompositePulse(TF1& func, const double& g, const double& p)
        : pulseFunc_{func}, gain_{g}, pedestal_{p} {}
 
    void addOrMerge(const std::pair<double, double>& hit, double hit_merge_ns) {
      auto imerge{hits_.begin()};
      for (; imerge != hits_.end(); imerge++)
        if (fabs(imerge->second - hit.second) < hit_merge_ns) break;
      if (imerge == hits_.end()) {  // didn't find a match, add to the list
        hits_.push_back(hit);
      } else {  // merge hits, shifting time to average
        imerge->second =
            (imerge->second * imerge->first + hit.first * hit.second);
        imerge->first += hit.first;
        imerge->second /= imerge->first;
      }
    }
 
    double findCrossing(double low, double high, double level,
                        double prec = 0.01) {
      // use midpoint algorithm, assumes low is below and high is above
      double step = high - low;
      double pt = (high + low) / 2;
      while (step > prec) {
        double vmid = at(pt);
        if (vmid < level) {
          low = pt;
        } else {
          high = pt;
        }
        step = high - low;
        pt = (high + low) / 2;
      }
      return pt;
    }
 
    void setGainPedestal(double gain, double pedestal) {
      gain_ = gain;
      pedestal_ = pedestal;
    }
 
    double operator()(double time) const { return at(time); }
 
    double at(double time) const {
      double signal = gain_ * pedestal_;
      for (auto hit : hits_)
        signal += hit.first * pulseFunc_.Eval(time - hit.second);
      return signal;
    };
 
    const std::vector<std::pair<double, double>>& hits() const { return hits_; }
 
   private:
    std::vector<std::pair<double, double>> hits_;
 
    TF1& pulseFunc_;
 
    double gain_;
 
    double pedestal_;
 
  };  // CompositePulse
 
 private:
  /**************************************************************************************
   * Parameters Identical for all Chips
   *************************************************************************************/
 
  bool noise_{true};
 
  int nADCs_;
 
  int iSOI_;
 
  double clockCycle_;
 
  double timingJitter_;
 
  double ns_;
 
  double rateUpSlope_;
 
  double timeUpSlope_;
 
  double rateDnSlope_;
 
  double timeDnSlope_;
 
  double timePeak_;
 
  double hit_merge_ns_;
 
  /**************************************************************************************
   * Chip-Dependent Parameters (Conditions)
   *************************************************************************************/
 
  const conditions::DoubleTableCondition* chipConditions_{nullptr};
 
  mutable std::map<std::string, int> conditionNamesToIndex_;
 
  /**************************************************************************************
   * Helpful Member Objects
   *************************************************************************************/
 
  std::unique_ptr<TRandom3> noiseInjector_;
 
  mutable TF1 pulseFunc_;
 
};  // HgcrocEmulator
 
}  // namespace ldmx
 
#endif  // TOOLS_HGCROCEMULATOR_H