#include <HPDIonFeedbackSim.h>

Inheritance diagram for HPDIonFeedbackSim:

Public Member Functions
void	addThermalNoise (CaloSamples &samples)
virtual double	correctPE (const DetId &detId, double npe) const
double	getIonFeedback (DetId detId, double signal, double pedWidth, bool doThermal, bool isInGeV)
	HPDIonFeedbackSim (const edm::ParameterSet &)
void	setDbService (const HcalDbService *service)
void	setRandomEngine (CLHEP::HepRandomEngine &engine)
void	setShape (const CaloVShape *shape)
	need a shaper in order to set thermal noise
	~HPDIonFeedbackSim ()
Private Member Functions
double	fCtoGeV (const DetId &detId) const
Private Attributes
const HcalDbService *	theDbService
CLHEP::RandBinomial *	theRandBinomial
CLHEP::RandFlat *	theRandFlat
CLHEP::RandGaussQ *	theRandGauss
CLHEP::RandPoissonQ *	theRandPoissonQ
const CaloVShape *	theShape

Detailed Description

Definition at line 28 of file HPDIonFeedbackSim.h.

Constructor & Destructor Documentation

HPDIonFeedbackSim::HPDIonFeedbackSim ( const edm::ParameterSet & iConfig ) [explicit]

Definition at line 31 of file HPDIonFeedbackSim.cc.

: theDbService(0), theShape(0),
 theRandBinomial(0), theRandFlat(0), theRandGauss(0), theRandPoissonQ(0)
{
}

HPDIonFeedbackSim::~HPDIonFeedbackSim ( )

Definition at line 37 of file HPDIonFeedbackSim.cc.

References theRandBinomial, theRandFlat, theRandGauss, and theRandPoissonQ.

{
    if(theRandBinomial) delete theRandBinomial;
    if (theRandFlat) delete theRandFlat;
    if (theRandGauss) delete theRandGauss;
    if (theRandPoissonQ) delete theRandPoissonQ;
}

Member Function Documentation

void HPDIonFeedbackSim::addThermalNoise ( CaloSamples & samples )

Definition at line 121 of file HPDIonFeedbackSim.cc.

References correctPE(), i, CaloSamples::id(), j, CaloSamples::size(), and theRandPoissonQ.

Referenced by HcalAmplifier::amplify().

{
  // make some chance to add a PE (with a chance of feedback)
  // for each time sample
  double meanPE = 0.02;
  DetId detId(samples.id());
  int nSamples = samples.size();
  for(int i = 0; i < nSamples; ++i) 
  {
    double npe = theRandPoissonQ->fire(meanPE);
    // TODOprobably should time-smear these
    if(npe > 0.)
    {
      // chance of feedback
      npe = correctPE(detId, npe);
      for(int j = i; j < nSamples; ++j)
      {
        double timeFromPE = (j-i) * 25.;
        samples[j] += (*theShape)(timeFromPE) * npe;
      }
    }
  }
}

double HPDIonFeedbackSim::correctPE	(	const DetId &	detId,
		double	npe
	)		const `[virtual]`

Implements CaloVPECorrection.

Definition at line 84 of file HPDIonFeedbackSim.cc.

References j, max(), p4, pe2Charge, theRandBinomial, and theRandGauss.

Referenced by addThermalNoise(), and getIonFeedback().

{
    double rateInTail = 0.000211988;//read this from XML file
    double rateInSecondTail = 4.61579e-06;//read this from XML file

    // three gauss fit is applied to data to get ion feedback distribution
    // parameters (in fC)
    // first gaussian
    // double p0 = 9.53192e+05;
    // double p1 = -3.13653e-01;
    // double p2 = 2.78350e+00;

    // second gaussian
    // double p3 = 2.41611e+03;
    double p4 = 2.06117e+01;
    double p5 = 1.09239e+01;

    // third gaussian
    // double p6 = 3.42793e+01;
    double p7 = 5.45548e+01;
    double p8 = 1.59696e+01;

    double noise = 0.;          // fC
    int nFirst  = (int)(theRandBinomial->fire(npe, rateInTail));
    int nSecond = (int)(theRandBinomial->fire(npe, rateInSecondTail));

    for (int j = 0; j < nFirst; ++j) {
      noise += theRandGauss->fire(p4, p5);
    }
    for (int j = 0; j < nSecond; ++j) {
      noise += theRandGauss->fire(p7, p8);
    }

    return npe + std::max(noise/pe2Charge, 0.);
}

double HPDIonFeedbackSim::fCtoGeV ( const DetId & detId ) const [private]

Definition at line 146 of file HPDIonFeedbackSim.cc.

References HcalDbService::getGain(), HcalDbService::getGainWidth(), HcalGain::getValue(), query::result, and theDbService.

Referenced by getIonFeedback().

{
  assert(theDbService != 0);
  HcalGenericDetId hcalGenDetId(detId);
  const HcalGain* gains = theDbService->getGain(hcalGenDetId);
  const HcalGainWidth* gwidths = theDbService->getGainWidth(hcalGenDetId);
  if (!gains || !gwidths )
  {
    edm::LogError("HcalAmplifier") << "Could not fetch HCAL conditions for channel " << hcalGenDetId;
  }
  // only one gain will be recorded per channel, so just use capID 0 for now
  double result = gains->getValue(0);
  return result;
}

double HPDIonFeedbackSim::getIonFeedback	(	DetId	detId,
		double	signal,
		double	pedWidth,
		bool	doThermal,
		bool	isInGeV
	)

Definition at line 59 of file HPDIonFeedbackSim.cc.

References DeDxDiscriminatorTools::charge(), correctPE(), fCtoGeV(), pe2Charge, and theRandPoissonQ.

{
    
    HcalDetId id = detId;
    
    double GeVperfC = 1.;
    if(isInGeV) GeVperfC = 1./fCtoGeV(detId);
    
    double charge = signal / GeVperfC;
    
    double noise = 0.;          // fC
    if (charge > 3. * pedWidth) {    // 3 sigma away from pedestal mean
        int npe = int (charge / pe2Charge);
        if(doThermal) {
          double electronEmission = 0.08;
          npe += theRandPoissonQ->fire(electronEmission);
        }

        noise = correctPE(detId, npe) - npe;
    }
    return (noise * GeVperfC);

}

void HPDIonFeedbackSim::setDbService ( const HcalDbService * service ) [inline]

Definition at line 35 of file HPDIonFeedbackSim.h.

References theDbService.

Referenced by HcalAmplifier::setDbService().

{theDbService = service;}

void HPDIonFeedbackSim::setRandomEngine ( CLHEP::HepRandomEngine & engine )

Definition at line 46 of file HPDIonFeedbackSim.cc.

References theRandBinomial, theRandFlat, theRandGauss, and theRandPoissonQ.

Referenced by HcalDigitizer::HcalDigitizer(), and HcalAmplifier::setRandomEngine().

{
    if(theRandBinomial) delete theRandBinomial;
    if (theRandFlat) delete theRandFlat;
    if (theRandGauss) delete theRandGauss;
    if (theRandPoissonQ) delete theRandPoissonQ;

    theRandBinomial = new CLHEP::RandBinomial(engine);
    theRandFlat = new CLHEP::RandFlat(engine);
    theRandGauss = new CLHEP::RandGaussQ(engine);
    theRandPoissonQ = new CLHEP::RandPoissonQ(engine);
}