HcalSimAlgos/src/HcalAmplifier.cc

0001 #include "SimCalorimetry/HcalSimAlgos/interface/HcalAmplifier.h"
0002 #include "SimCalorimetry/HcalSimAlgos/interface/HcalSimParameters.h"
0003 #include "SimCalorimetry/HcalSimAlgos/interface/HPDIonFeedbackSim.h"
0004 #include "SimCalorimetry/HcalSimAlgos/interface/HcalTimeSlewSim.h"
0005 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVSimParameterMap.h"
0006 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVNoiseSignalGenerator.h"
0007 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
0008 #include "CondFormats/HcalObjects/interface/HcalPedestal.h"
0009 #include "CondFormats/HcalObjects/interface/HcalGain.h"
0010 #include "CondFormats/HcalObjects/interface/HcalPedestalWidth.h"
0011 #include "CondFormats/HcalObjects/interface/HcalGainWidth.h"
0012 #include "CalibFormats/CaloObjects/interface/CaloSamples.h"
0013 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016 #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h"
0017
0018 #include "CLHEP/Random/RandGaussQ.h"
0019 #include "CLHEP/Random/RandFlat.h"
0020
0021 #include <cmath>
0022 #include <cmath>
0023
0024 HcalAmplifier::HcalAmplifier(const CaloVSimParameterMap* parameters, bool addNoise, bool PreMix1, bool PreMix2)
0025     : theDbService(nullptr),
0026       theParameterMap(parameters),
0027       theNoiseSignalGenerator(nullptr),
0028       theIonFeedbackSim(nullptr),
0029       theTimeSlewSim(nullptr),
0030       theStartingCapId(0),
0031       addNoise_(addNoise),
0032       preMixDigi_(PreMix1),
0033       preMixAdd_(PreMix2) {}
0034
0035 void HcalAmplifier::setDbService(const HcalDbService* service) {
0036   theDbService = service;
0037   if (theIonFeedbackSim)
0038     theIonFeedbackSim->setDbService(service);
0039 }
0040
0041 void HcalAmplifier::amplify(CaloSamples& frame, CLHEP::HepRandomEngine* engine) const {
0042   if (theIonFeedbackSim) {
0043     theIonFeedbackSim->addThermalNoise(frame, engine);
0044   }
0045   pe2fC(frame);
0046
0047   if (frame.id().det() == DetId::Hcal && ((frame.id().subdetId() == HcalGenericDetId::HcalGenBarrel) ||
0048                                           (frame.id().subdetId() == HcalGenericDetId::HcalGenEndcap))) {
0049     const HcalSimParameters& params = static_cast<const HcalSimParameters&>(theParameterMap->simParameters(frame.id()));
0050     if (params.delayQIE() > 0)
0051       applyQIEdelay(frame, params.delayQIE());
0052   }
0053
0054   // don't bother for blank signals
0055   if (theTimeSlewSim && frame.size() > 4 && frame[4] > 1.e-6) {
0056     theTimeSlewSim->delay(frame, engine, theTimeSlew);
0057   }
0058
0059   // if we are combining pre-mixed digis, we need noise and peds
0060   if (theNoiseSignalGenerator == nullptr || preMixAdd_ || !theNoiseSignalGenerator->contains(frame.id())) {
0061     addPedestals(frame, engine);
0062   }
0063   LogDebug("HcalAmplifier") << frame;
0064 }
0065
0066 void HcalAmplifier::pe2fC(CaloSamples& frame) const {
0067   const CaloSimParameters& parameters = theParameterMap->simParameters(frame.id());
0068   frame *= parameters.photoelectronsToAnalog(frame.id());
0069 }
0070
0071 void HcalAmplifier::applyQIEdelay(CaloSamples& cs, int delayQIE) const {
0072   DetId detId(cs.id());
0073   int maxbin = cs.size();
0074   int precisebin = cs.preciseSize();
0075   CaloSamples data(detId, maxbin, precisebin);  // make a temporary copy
0076   data = cs;
0077   data.setBlank();
0078   data.resetPrecise();
0079
0080   for (int i = 0; i < precisebin; i++) {
0081     if (i < 2 * delayQIE)
0082       data.preciseAtMod(i) += 0.;
0083     else
0084       data.preciseAtMod(i) += cs.preciseAt(i - 2 * delayQIE);
0085     int samplebin = (int)i * maxbin / precisebin;
0086     data[samplebin] += data.preciseAt(i);
0087   }
0088
0089   cs = data;  // update the sample
0090 }
0091
0092 void HcalAmplifier::addPedestals(CaloSamples& frame, CLHEP::HepRandomEngine* engine) const {
0093   assert(theDbService != nullptr);
0094   HcalGenericDetId hcalGenDetId(frame.id());
0095   HcalGenericDetId::HcalGenericSubdetector hcalSubDet = hcalGenDetId.genericSubdet();
0096
0097   if (!((frame.id().subdetId() == HcalGenericDetId::HcalGenBarrel) ||
0098         (frame.id().subdetId() == HcalGenericDetId::HcalGenEndcap) ||
0099         (frame.id().subdetId() == HcalGenericDetId::HcalGenForward) ||
0100         (frame.id().subdetId() == HcalGenericDetId::HcalGenOuter)))
0101     return;
0102
0103   if (hcalGenDetId.isHcalCastorDetId())
0104     return;
0105   if (hcalGenDetId.isHcalZDCDetId())
0106     return;
0107
0108   const HcalCalibrationWidths& calibWidths = theDbService->getHcalCalibrationWidths(hcalGenDetId);
0109   const HcalCalibrations& calibs = theDbService->getHcalCalibrations(hcalGenDetId);
0110
0111   double noise[32] = {0.};  //big enough
0112   if (addNoise_) {
0113     double gauss[32];  //big enough
0114     for (int i = 0; i < frame.size(); i++)
0115       gauss[i] = CLHEP::RandGaussQ::shoot(engine, 0., 1.);
0116     makeNoise(hcalSubDet, calibWidths, frame.size(), gauss, noise);
0117   }
0118
0119   if (!preMixDigi_) {  // if we are doing initial premix, no pedestals
0120     for (int tbin = 0; tbin < frame.size(); ++tbin) {
0121       int capId = (theStartingCapId + tbin) % 4;
0122       double pedestal = calibs.pedestal(capId) + noise[tbin];
0123       frame[tbin] += pedestal;
0124     }
0125   }
0126 }
0127
0128 void HcalAmplifier::makeNoise(HcalGenericDetId::HcalGenericSubdetector hcalSubDet,
0129                               const HcalCalibrationWidths& width,
0130                               int fFrames,
0131                               double* fGauss,
0132                               double* fNoise) const {
0133   // This is a simplified noise generation scheme using only the diagonal elements
0134   // (proposed by Salavat Abduline).
0135   // This is direct adaptation of the code in HcalPedestalWidth.cc
0136
0137   // average over capId's
0138   double s_xx_mean = 0.25 * (width.pedestal(0) * width.pedestal(0) + width.pedestal(1) * width.pedestal(1) +
0139                              width.pedestal(2) * width.pedestal(2) + width.pedestal(3) * width.pedestal(3));
0140
0141   // Off-diagonal element approximation
0142   // In principle should come from averaging the values of elements (0.1), (1,2), (2,3), (3,0)
0143   // For now use the definition below (but keep structure of the code structure for development)
0144   double s_xy_mean = -0.5 * s_xx_mean;
0145   // Use different parameter for HF to reproduce the noise rate after zero suppression.
0146   // Steven Won/Jim Hirschauer/Radek Ofierzynski 18.03.2010
0147   if (hcalSubDet == HcalGenericDetId::HcalGenForward)
0148     s_xy_mean = 0.08 * s_xx_mean;
0149
0150   double term = s_xx_mean * s_xx_mean - 2. * s_xy_mean * s_xy_mean;
0151
0152   if (term < 0.)
0153     term = 1.e-50;
0154   double sigma = sqrt(0.5 * (s_xx_mean + sqrt(term)));
0155   double corr = sigma == 0. ? 0. : 0.5 * s_xy_mean / sigma;
0156
0157   for (int i = 0; i < fFrames; i++) {
0158     fNoise[i] = fGauss[i] * sigma;
0159     if (i > 0)
0160       fNoise[i] += fGauss[i - 1] * corr;
0161     if (i < fFrames - 1)
0162       fNoise[i] += fGauss[i + 1] * corr;
0163   }
0164 }