HcalSimAlgos/src/HPDIonFeedbackSim.cc

0001
0002 // --------------------------------------------------------
0003 // A class to simulated HPD ion feedback noise.
0004 // The deliverable of the class is the ion feedback noise
0005 // for an HcalDetId units of fC or GeV
0006 //
0007 // Project: HPD ion feedback
0008 // Author: T.Yetkin University of Iowa, Feb. 16, 2010
0009 // --------------------------------------------------------
0010
0011 #include "SimCalorimetry/HcalSimAlgos/interface/HPDIonFeedbackSim.h"
0012 #include "SimCalorimetry/HcalSimAlgos/interface/HcalShapes.h"
0013 #include "SimCalorimetry/CaloSimAlgos/interface/CaloVShape.h"
0014 #include "CondFormats/HcalObjects/interface/HcalGain.h"
0015 #include "CondFormats/HcalObjects/interface/HcalGainWidth.h"
0016 #include "DataFormats/HcalDetId/interface/HcalGenericDetId.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019 #include "FWCore/PluginManager/interface/PluginManager.h"
0020 #include "FWCore/PluginManager/interface/standard.h"
0021
0022 #include "CLHEP/Random/RandBinomial.h"
0023 #include "CLHEP/Random/RandGaussQ.h"
0024 #include "CLHEP/Random/RandPoissonQ.h"
0025
0026 using namespace edm;
0027 using namespace std;
0028
0029 // constants for simulation/parameterization
0030 static const double pe2Charge = 0.333333;  // fC/p.e.
0031
0032 HPDIonFeedbackSim::HPDIonFeedbackSim(const edm::ParameterSet& iConfig, const CaloShapes* shapes)
0033     : theDbService(nullptr), theShapes(shapes) {}
0034
0035 HPDIonFeedbackSim::~HPDIonFeedbackSim() {}
0036
0037 double HPDIonFeedbackSim::getIonFeedback(
0038     DetId detId, double signal, double pedWidth, bool doThermal, bool isInGeV, CLHEP::HepRandomEngine* engine) {
0039   //    HcalDetId id = detId;
0040
0041   double GeVperfC = 1.;
0042   if (isInGeV)
0043     GeVperfC = 1. / fCtoGeV(detId);
0044
0045   double charge = signal / GeVperfC;
0046
0047   double noise = 0.;             // fC
0048   if (charge > 3. * pedWidth) {  // 3 sigma away from pedestal mean
0049     int npe = int(charge / pe2Charge);
0050     if (doThermal) {
0051       double electronEmission = 0.08;
0052       CLHEP::RandPoissonQ theRandPoissonQ(*engine, electronEmission);
0053       npe += theRandPoissonQ.fire();
0054     }
0055
0056     noise = correctPE(detId, npe, engine) - npe;
0057   }
0058   return (noise * GeVperfC);
0059 }
0060
0061 double HPDIonFeedbackSim::correctPE(const DetId& detId, double npe, CLHEP::HepRandomEngine* engine) const {
0062   double rateInTail = 0.000211988;        //read this from XML file
0063   double rateInSecondTail = 4.61579e-06;  //read this from XML file
0064
0065   // three gauss fit is applied to data to get ion feedback distribution
0066   // parameters (in fC)
0067   // first gaussian
0068   // double p0 = 9.53192e+05;
0069   // double p1 = -3.13653e-01;
0070   // double p2 = 2.78350e+00;
0071
0072   // second gaussian
0073   // double p3 = 2.41611e+03;
0074   double p4 = 2.06117e+01;
0075   double p5 = 1.09239e+01;
0076
0077   // third gaussian
0078   // double p6 = 3.42793e+01;
0079   double p7 = 5.45548e+01;
0080   double p8 = 1.59696e+01;
0081
0082   double noise = 0.;  // fC
0083   int nFirst = (int)(CLHEP::RandBinomial::shoot(engine, npe, rateInTail));
0084   int nSecond = (int)(CLHEP::RandBinomial::shoot(engine, npe, rateInSecondTail));
0085
0086   for (int j = 0; j < nFirst; ++j) {
0087     noise += CLHEP::RandGaussQ::shoot(engine, p4, p5);
0088   }
0089   for (int j = 0; j < nSecond; ++j) {
0090     noise += CLHEP::RandGaussQ::shoot(engine, p7, p8);
0091   }
0092
0093   return npe + std::max(noise / pe2Charge, 0.);
0094 }
0095
0096 void HPDIonFeedbackSim::addThermalNoise(CaloSamples& samples, CLHEP::HepRandomEngine* engine) {
0097   // make some chance to add a PE (with a chance of feedback)
0098   // for each time sample
0099   double meanPE = 0.02;
0100   DetId detId(samples.id());
0101   int nSamples = samples.size();
0102   const CaloVShape* shape = theShapes->shape(detId);
0103   for (int i = 0; i < nSamples; ++i) {
0104     CLHEP::RandPoissonQ theRandPoissonQ(*engine, meanPE);
0105     double npe = theRandPoissonQ.fire();
0106     // TODOprobably should time-smear these
0107     if (npe > 0.) {
0108       // chance of feedback
0109       npe = correctPE(detId, npe, engine);
0110       for (int j = i; j < nSamples; ++j) {
0111         double timeFromPE = (j - i) * 25.;
0112         samples[j] += (*shape)(timeFromPE)*npe;
0113       }
0114     }
0115   }
0116 }
0117
0118 double HPDIonFeedbackSim::fCtoGeV(const DetId& detId) const {
0119   assert(theDbService != nullptr);
0120   HcalGenericDetId hcalGenDetId(detId);
0121   const HcalGain* gains = theDbService->getGain(hcalGenDetId);
0122   const HcalGainWidth* gwidths = theDbService->getGainWidth(hcalGenDetId);
0123   double result = 0.0;
0124   if (!gains || !gwidths) {
0125     edm::LogError("HcalAmplifier") << "Could not fetch HCAL conditions for channel " << hcalGenDetId;
0126   } else {
0127     // only one gain will be recorded per channel, so just use capID 0 for now
0128     result = gains->getValue(0);
0129   }
0130   return result;
0131 }