HcalSimAlgos/src/HcalSimParameters.cc

0001 #include "SimCalorimetry/HcalSimAlgos/interface/HcalSimParameters.h"
0002 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0003 #include "DataFormats/HcalDetId/interface/HcalGenericDetId.h"
0004 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "CondFormats/HcalObjects/interface/HcalGain.h"
0007 #include "CondFormats/HcalObjects/interface/HcalGainWidth.h"
0008 #include "CalibFormats/HcalObjects/interface/HcalSiPMType.h"
0009 #include "CLHEP/Random/RandGaussQ.h"
0010 using namespace std;
0011
0012 HcalSimParameters::HcalSimParameters(double simHitToPhotoelectrons,
0013                                      double samplingFactor,
0014                                      double timePhase,
0015                                      int readoutFrameSize,
0016                                      int binOfMaximum,
0017                                      bool doPhotostatistics,
0018                                      bool syncPhase,
0019                                      int firstRing,
0020                                      const std::vector<double>& samplingFactors,
0021                                      double sipmTau)
0022     : CaloSimParameters(simHitToPhotoelectrons,
0023                         0.0,
0024                         samplingFactor,
0025                         timePhase,
0026                         readoutFrameSize,
0027                         binOfMaximum,
0028                         doPhotostatistics,
0029                         syncPhase),
0030       theDbService(nullptr),
0031       theSiPMcharacteristics(nullptr),
0032       theFirstRing(firstRing),
0033       theSamplingFactors(samplingFactors),
0034       theSiPMSmearing(false),
0035       doTimeSmear_(true),
0036       theSiPMTau(sipmTau) {
0037   defaultTimeSmearing();
0038
0039   edm::LogInfo("HcalSimParameters:") << " doSiPMsmearing    = " << theSiPMSmearing;
0040 }
0041
0042 HcalSimParameters::HcalSimParameters(const edm::ParameterSet& p)
0043     : CaloSimParameters(p, true),
0044       theDbService(nullptr),
0045       theFirstRing(p.getParameter<int>("firstRing")),
0046       theSamplingFactors(p.getParameter<std::vector<double> >("samplingFactors")),
0047       theSiPMSmearing(p.getParameter<bool>("doSiPMSmearing")),
0048       doTimeSmear_(p.getParameter<bool>("timeSmearing")),
0049       theSiPMTau(p.getParameter<double>("sipmTau")),
0050       threshold_currentTDC_(p.getParameter<double>("threshold_currentTDC")),
0051       delayQIE_(p.getParameter<int>("delayQIE")) {
0052   defaultTimeSmearing();
0053
0054   edm::LogInfo("HcalSimParameters:") << " doSiPMsmearing    = " << theSiPMSmearing;
0055 }
0056
0057 void HcalSimParameters::setDbService(const HcalDbService* service) {
0058   assert(service);
0059   theDbService = service;
0060   theSiPMcharacteristics = service->getHcalSiPMCharacteristics();
0061   assert(theSiPMcharacteristics);
0062 }
0063
0064 double HcalSimParameters::simHitToPhotoelectrons(const DetId& detId) const {
0065   // the gain is in units of GeV/fC.  We want a constant with pe/dGeV
0066   // pe/dGeV = (GeV/dGeV) / (GeV/fC) / (fC/pe)
0067   double result = CaloSimParameters::simHitToPhotoelectrons(detId);
0068   if (HcalGenericDetId(detId).genericSubdet() != HcalGenericDetId::HcalGenForward ||
0069       HcalGenericDetId(detId).genericSubdet() != HcalGenericDetId::HcalGenZDC) {
0070     result = samplingFactor(detId) / fCtoGeV(detId) / photoelectronsToAnalog(detId);
0071   }
0072   return result;
0073 }
0074
0075 double HcalSimParameters::fCtoGeV(const DetId& detId) const {
0076   assert(theDbService != nullptr);
0077   HcalGenericDetId hcalGenDetId(detId);
0078   const HcalGain* gains = theDbService->getGain(hcalGenDetId);
0079   const HcalGainWidth* gwidths = theDbService->getGainWidth(hcalGenDetId);
0080   double result = 0.0;
0081   if (!gains || !gwidths) {
0082     edm::LogError("HcalAmplifier") << "Could not fetch HCAL conditions for channel " << hcalGenDetId;
0083   } else {
0084     // only one gain will be recorded per channel, so just use capID 0 for now
0085     result = gains->getValue(0);
0086     //  if(doNoise_)
0087     ///  {
0088     //    result += CLHEP::RandGaussQ::shoot(0.,  gwidths->getValue(0));
0089     //  }
0090   }
0091   return result;
0092 }
0093
0094 double HcalSimParameters::samplingFactor(const DetId& detId) const {
0095   HcalDetId hcalDetId(detId);
0096   return theSamplingFactors.at(hcalDetId.ietaAbs() - theFirstRing);
0097 }
0098
0099 double HcalSimParameters::photoelectronsToAnalog(const DetId& detId) const {
0100   //now always taken from database for HPDs or SiPMs (HB, HE, HO)
0101   assert(theDbService);
0102   return theDbService->getHcalSiPMParameter(detId)->getFCByPE();
0103 }
0104
0105 //static const double GeV2fC = 1.0/0.145;
0106 static const double GeV2fC = 1.0 / 0.4;
0107
0108 void HcalSimParameters::defaultTimeSmearing() {
0109   // GeV->ampl (fC), time (ns)
0110   theSmearSettings.emplace_back(4.00 * GeV2fC, 4.050);
0111   theSmearSettings.emplace_back(20.00 * GeV2fC, 3.300);
0112   theSmearSettings.emplace_back(25.00 * GeV2fC, 2.925);
0113   theSmearSettings.emplace_back(30.00 * GeV2fC, 2.714);
0114   theSmearSettings.emplace_back(37.00 * GeV2fC, 2.496);
0115   theSmearSettings.emplace_back(44.50 * GeV2fC, 2.278);
0116   theSmearSettings.emplace_back(56.00 * GeV2fC, 2.138);
0117   theSmearSettings.emplace_back(63.50 * GeV2fC, 2.022);
0118   theSmearSettings.emplace_back(81.00 * GeV2fC, 1.788);
0119   theSmearSettings.emplace_back(88.50 * GeV2fC, 1.695);
0120   theSmearSettings.emplace_back(114.50 * GeV2fC, 1.716);
0121   theSmearSettings.emplace_back(175.50 * GeV2fC, 1.070);
0122   theSmearSettings.emplace_back(350.00 * GeV2fC, 1.564);
0123   theSmearSettings.emplace_back(99999.00 * GeV2fC, 1.564);
0124 }
0125
0126 double HcalSimParameters::timeSmearRMS(double ampl) const {
0127   HcalTimeSmearSettings::size_type i;
0128   double smearsigma = 0;
0129
0130   for (i = 0; i < theSmearSettings.size(); i++)
0131     if (theSmearSettings[i].first > ampl)
0132       break;
0133
0134   // Smearing occurs only within the envelope definitions.
0135   if (i != 0 && (i < theSmearSettings.size())) {
0136     double energy1 = theSmearSettings[i - 1].first;
0137     double sigma1 = theSmearSettings[i - 1].second;
0138     double energy2 = theSmearSettings[i].first;
0139     double sigma2 = theSmearSettings[i].second;
0140
0141     if (energy2 != energy1)
0142       smearsigma = sigma1 + ((sigma2 - sigma1) * (ampl - energy1) / (energy2 - energy1));
0143     else
0144       smearsigma = (sigma2 + sigma1) / 2.;
0145   }
0146
0147   return smearsigma;
0148 }
0149
0150 int HcalSimParameters::pixels(const DetId& detId) const {
0151   assert(theDbService);
0152   int type = theDbService->getHcalSiPMParameter(detId)->getType();
0153   return theSiPMcharacteristics->getPixels(type);
0154 }
0155
0156 double HcalSimParameters::sipmDarkCurrentuA(const DetId& detId) const {
0157   assert(theDbService);
0158   return theDbService->getHcalSiPMParameter(detId)->getDarkCurrent();
0159 }
0160
0161 double HcalSimParameters::sipmCrossTalk(const DetId& detId) const {
0162   assert(theDbService);
0163   int type = theDbService->getHcalSiPMParameter(detId)->getType();
0164   return theSiPMcharacteristics->getCrossTalk(type);
0165 }
0166 std::vector<float> HcalSimParameters::sipmNonlinearity(const DetId& detId) const {
0167   assert(theDbService);
0168   int type = theDbService->getHcalSiPMParameter(detId)->getType();
0169   return theSiPMcharacteristics->getNonLinearities(type);
0170 }
0171
0172 unsigned int HcalSimParameters::signalShape(const DetId& detId) const {
0173   assert(theDbService);
0174   return theDbService->getHcalMCParam(detId)->signalShape();
0175 }