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 (abs(params.delayQIE()) <= 25)
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     int j = i + 2 * delayQIE;
0082     data.preciseAtMod(i) +=
0083         // value positive (signal moves earlier in time)
0084         delayQIE > 0 ? (j < precisebin ? cs.preciseAt(j) : 0.) :
0085                      //  value = 0 or negative (signal gets delayed)
0086             (j < 0 ? 0. : cs.preciseAt(j));
0087
0088     int samplebin = (int)i * maxbin / precisebin;
0089     data[samplebin] += data.preciseAt(i);
0090   }
0091
0092   cs = data;  // update the sample
0093 }
0094
0095 void HcalAmplifier::addPedestals(CaloSamples& frame, CLHEP::HepRandomEngine* engine) const {
0096   assert(theDbService != nullptr);
0097   HcalGenericDetId hcalGenDetId(frame.id());
0098   HcalGenericDetId::HcalGenericSubdetector hcalSubDet = hcalGenDetId.genericSubdet();
0099
0100   if (!((frame.id().subdetId() == HcalGenericDetId::HcalGenBarrel) ||
0101         (frame.id().subdetId() == HcalGenericDetId::HcalGenEndcap) ||
0102         (frame.id().subdetId() == HcalGenericDetId::HcalGenForward) ||
0103         (frame.id().subdetId() == HcalGenericDetId::HcalGenOuter)))
0104     return;
0105
0106   if (hcalGenDetId.isHcalCastorDetId())
0107     return;
0108   if (hcalGenDetId.isHcalZDCDetId())
0109     return;
0110
0111   const HcalCalibrationWidths& calibWidths = theDbService->getHcalCalibrationWidths(hcalGenDetId);
0112   const HcalCalibrations& calibs = theDbService->getHcalCalibrations(hcalGenDetId);
0113
0114   double noise[32] = {0.};  //big enough
0115   if (addNoise_) {
0116     double gauss[32];  //big enough
0117     for (int i = 0; i < frame.size(); i++)
0118       gauss[i] = CLHEP::RandGaussQ::shoot(engine, 0., 1.);
0119     makeNoise(hcalSubDet, calibWidths, frame.size(), gauss, noise);
0120   }
0121
0122   if (!preMixDigi_) {  // if we are doing initial premix, no pedestals
0123     for (int tbin = 0; tbin < frame.size(); ++tbin) {
0124       int capId = (theStartingCapId + tbin) % 4;
0125       double pedestal = calibs.pedestal(capId) + noise[tbin];
0126       frame[tbin] += pedestal;
0127     }
0128   }
0129 }
0130
0131 void HcalAmplifier::makeNoise(HcalGenericDetId::HcalGenericSubdetector hcalSubDet,
0132                               const HcalCalibrationWidths& width,
0133                               int fFrames,
0134                               double* fGauss,
0135                               double* fNoise) const {
0136   // This is a simplified noise generation scheme using only the diagonal elements
0137   // (proposed by Salavat Abduline).
0138   // This is direct adaptation of the code in HcalPedestalWidth.cc
0139
0140   // average over capId's
0141   double s_xx_mean = 0.25 * (width.pedestal(0) * width.pedestal(0) + width.pedestal(1) * width.pedestal(1) +
0142                              width.pedestal(2) * width.pedestal(2) + width.pedestal(3) * width.pedestal(3));
0143
0144   // Off-diagonal element approximation
0145   // In principle should come from averaging the values of elements (0.1), (1,2), (2,3), (3,0)
0146   // For now use the definition below (but keep structure of the code structure for development)
0147   double s_xy_mean = -0.5 * s_xx_mean;
0148   // Use different parameter for HF to reproduce the noise rate after zero suppression.
0149   // Steven Won/Jim Hirschauer/Radek Ofierzynski 18.03.2010
0150   if (hcalSubDet == HcalGenericDetId::HcalGenForward)
0151     s_xy_mean = 0.08 * s_xx_mean;
0152
0153   double term = s_xx_mean * s_xx_mean - 2. * s_xy_mean * s_xy_mean;
0154
0155   if (term < 0.)
0156     term = 1.e-50;
0157   double sigma = sqrt(0.5 * (s_xx_mean + sqrt(term)));
0158   double corr = sigma == 0. ? 0. : 0.5 * s_xy_mean / sigma;
0159
0160   for (int i = 0; i < fFrames; i++) {
0161     fNoise[i] = fGauss[i] * sigma;
0162     if (i > 0)
0163       fNoise[i] += fGauss[i - 1] * corr;
0164     if (i < fFrames - 1)
0165       fNoise[i] += fGauss[i + 1] * corr;
0166   }
0167 }