text: #include "SimCalorimetry/HcalSimAlgos/interface/HcalAmplifier.h"
  #include "SimCalorimetry/HcalSimAlgos/interface/HcalSimParameters.h"
  #include "SimCalorimetry/HcalSimAlgos/interface/HPDIonFeedbackSim.h"
  #include "SimCalorimetry/HcalSimAlgos/interface/HcalTimeSlewSim.h"
  #include "SimCalorimetry/CaloSimAlgos/interface/CaloVSimParameterMap.h"
  #include "SimCalorimetry/CaloSimAlgos/interface/CaloVNoiseSignalGenerator.h"
  #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
  #include "CondFormats/HcalObjects/interface/HcalPedestal.h"
  #include "CondFormats/HcalObjects/interface/HcalGain.h"
  #include "CondFormats/HcalObjects/interface/HcalPedestalWidth.h"
  #include "CondFormats/HcalObjects/interface/HcalGainWidth.h"
  #include "CalibFormats/CaloObjects/interface/CaloSamples.h"
  #include "DataFormats/HcalDetId/interface/HcalDetId.h"
  #include "FWCore/MessageLogger/interface/MessageLogger.h"
  
  #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
  
  #include "CLHEP/Random/RandGaussQ.h"
  #include "CLHEP/Random/RandFlat.h"
  
  #include <cmath>
  #include <cmath>
  
  HcalAmplifier::HcalAmplifier(const CaloVSimParameterMap* parameters, bool addNoise, bool PreMix1, bool PreMix2)
      : theDbService(nullptr),
        theParameterMap(parameters),
        theNoiseSignalGenerator(nullptr),
        theIonFeedbackSim(nullptr),
        theTimeSlewSim(nullptr),
        theStartingCapId(0),
        addNoise_(addNoise),
        preMixDigi_(PreMix1),
        preMixAdd_(PreMix2) {}
  
  void HcalAmplifier::setDbService(const HcalDbService* service) {
    theDbService = service;
    if (theIonFeedbackSim)
      theIonFeedbackSim->setDbService(service);
  }
  
  void HcalAmplifier::amplify(CaloSamples& frame, CLHEP::HepRandomEngine* engine) const {
    if (theIonFeedbackSim) {
      theIonFeedbackSim->addThermalNoise(frame, engine);
    }
    pe2fC(frame);
  
    if (frame.id().det() == DetId::Hcal && ((frame.id().subdetId() == HcalGenericDetId::HcalGenBarrel) ||
                                            (frame.id().subdetId() == HcalGenericDetId::HcalGenEndcap))) {
      const HcalSimParameters& params = static_cast<const HcalSimParameters&>(theParameterMap->simParameters(frame.id()));
      if (params.delayQIE() > 0)
        applyQIEdelay(frame, params.delayQIE());
    }
  
    // don't bother for blank signals
    if (theTimeSlewSim && frame.size() > 4 && frame[4] > 1.e-6) {
      theTimeSlewSim->delay(frame, engine, theTimeSlew);
    }
  
    // if we are combining pre-mixed digis, we need noise and peds
    if (theNoiseSignalGenerator == nullptr || preMixAdd_ || !theNoiseSignalGenerator->contains(frame.id())) {
      addPedestals(frame, engine);
    }
    LogDebug("HcalAmplifier") << frame;
  }
  
  void HcalAmplifier::pe2fC(CaloSamples& frame) const {
    const CaloSimParameters& parameters = theParameterMap->simParameters(frame.id());
    frame *= parameters.photoelectronsToAnalog(frame.id());
  }
  
  void HcalAmplifier::applyQIEdelay(CaloSamples& cs, int delayQIE) const {
    DetId detId(cs.id());
    int maxbin = cs.size();
    int precisebin = cs.preciseSize();
    CaloSamples data(detId, maxbin, precisebin);  // make a temporary copy
    data = cs;
    data.setBlank();
    data.resetPrecise();
  
    for (int i = 0; i < precisebin; i++) {
      if (i < 2 * delayQIE)
        data.preciseAtMod(i) += 0.;
      else
        data.preciseAtMod(i) += cs.preciseAt(i - 2 * delayQIE);
      int samplebin = (int)i * maxbin / precisebin;
      data[samplebin] += data.preciseAt(i);
    }
  
    cs = data;  // update the sample
  }
  
  void HcalAmplifier::addPedestals(CaloSamples& frame, CLHEP::HepRandomEngine* engine) const {
    assert(theDbService != nullptr);
    HcalGenericDetId hcalGenDetId(frame.id());
    HcalGenericDetId::HcalGenericSubdetector hcalSubDet = hcalGenDetId.genericSubdet();
  
    if (!((frame.id().subdetId() == HcalGenericDetId::HcalGenBarrel) ||
          (frame.id().subdetId() == HcalGenericDetId::HcalGenEndcap) ||
          (frame.id().subdetId() == HcalGenericDetId::HcalGenForward) ||
          (frame.id().subdetId() == HcalGenericDetId::HcalGenOuter)))
      return;
  
    if (hcalGenDetId.isHcalCastorDetId())
      return;
    if (hcalGenDetId.isHcalZDCDetId())
      return;
  
    const HcalCalibrationWidths& calibWidths = theDbService->getHcalCalibrationWidths(hcalGenDetId);
    const HcalCalibrations& calibs = theDbService->getHcalCalibrations(hcalGenDetId);
  
    double noise[32] = {0.};  //big enough
    if (addNoise_) {
      double gauss[32];  //big enough
      for (int i = 0; i < frame.size(); i++)
        gauss[i] = CLHEP::RandGaussQ::shoot(engine, 0., 1.);
      makeNoise(hcalSubDet, calibWidths, frame.size(), gauss, noise);
    }
  
    if (!preMixDigi_) {  // if we are doing initial premix, no pedestals
      for (int tbin = 0; tbin < frame.size(); ++tbin) {
        int capId = (theStartingCapId + tbin) % 4;
        double pedestal = calibs.pedestal(capId) + noise[tbin];
        frame[tbin] += pedestal;
      }
    }
  }
  
  void HcalAmplifier::makeNoise(HcalGenericDetId::HcalGenericSubdetector hcalSubDet,
                                const HcalCalibrationWidths& width,
                                int fFrames,
                                double* fGauss,
                                double* fNoise) const {
    // This is a simplified noise generation scheme using only the diagonal elements
    // (proposed by Salavat Abduline).
    // This is direct adaptation of the code in HcalPedestalWidth.cc
  
    // average over capId's
    double s_xx_mean = 0.25 * (width.pedestal(0) * width.pedestal(0) + width.pedestal(1) * width.pedestal(1) +
                               width.pedestal(2) * width.pedestal(2) + width.pedestal(3) * width.pedestal(3));
  
    // Off-diagonal element approximation
    // In principle should come from averaging the values of elements (0.1), (1,2), (2,3), (3,0)
    // For now use the definition below (but keep structure of the code structure for development)
    double s_xy_mean = -0.5 * s_xx_mean;
    // Use different parameter for HF to reproduce the noise rate after zero suppression.
    // Steven Won/Jim Hirschauer/Radek Ofierzynski 18.03.2010
    if (hcalSubDet == HcalGenericDetId::HcalGenForward)
      s_xy_mean = 0.08 * s_xx_mean;
  
    double term = s_xx_mean * s_xx_mean - 2. * s_xy_mean * s_xy_mean;
  
    if (term < 0.)
      term = 1.e-50;
    double sigma = sqrt(0.5 * (s_xx_mean + sqrt(term)));
    double corr = sigma == 0. ? 0. : 0.5 * s_xy_mean / sigma;
  
    for (int i = 0; i < fFrames; i++) {
      fNoise[i] = fGauss[i] * sigma;
      if (i > 0)
        fNoise[i] += fGauss[i - 1] * corr;
      if (i < fFrames - 1)
        fNoise[i] += fGauss[i + 1] * corr;
    }
  }
  