File indexing completed on 2021-02-14 14:28:39
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
0055 if (theTimeSlewSim && frame.size() > 4 && frame[4] > 1.e-6) {
0056 theTimeSlewSim->delay(frame, engine, theTimeSlew);
0057 }
0058
0059
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);
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;
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.};
0112 if (addNoise_) {
0113 double gauss[32];
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_) {
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
0134
0135
0136
0137
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
0142
0143
0144 double s_xy_mean = -0.5 * s_xx_mean;
0145
0146
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 }