In This Package:

EsPmtEffectPulseTool.cc

Go to the documentation of this file.

#include "EsPmtEffectPulseTool.h"

#include "Event/SimHit.h"
#include "Event/ElecPulse.h"
#include "Event/ElecRpcPulse.h"
#include "Event/SimTrackReference.h"
#include "Event/SimTrack.h"

#include "Conventions/PulseType.h"

#include "GaudiKernel/SystemOfUnits.h"
#include "GaudiKernel/PhysicalConstants.h"

#include <math.h>


// DWL: After-pulse timing distribution

#define NUM_BINS_DIST 96

static double afterPusleTimingDist[NUM_BINS_DIST] = {0.0000000, 0.0012950, 0.0030194, 0.0064297, 0.0142225,
                                                     0.0247712, 0.0417949, 0.0696709, 0.1063259, 0.1485606,
                                                     0.1957973, 0.2465997, 0.2957179, 0.3401377, 0.3798086,
                                                     0.4135611, 0.4437462, 0.4702155, 0.4928434, 0.5123140,
                                                     0.5289294, 0.5432514, 0.5556360, 0.5664919, 0.5758642,
                                                     0.5843533, 0.5920903, 0.5990121, 0.6052796, 0.6111246,
                                                     0.6165560, 0.6215807, 0.6261028, 0.6304032, 0.6346304,
                                                     0.6387545, 0.6425872, 0.6465577, 0.6500606, 0.6534989,
                                                     0.6568080, 0.6602148, 0.6635239, 0.6667457, 0.6700025, 
                                                     0.6737148, 0.6780693, 0.6829002, 0.6879599, 0.6935379, 
                                                     0.6994667, 0.7058650, 0.7121394, 0.7187122, 0.7255085, 
                                                     0.7324600, 0.7394151, 0.7467314, 0.7546551, 0.7630884, 
                                                     0.7718917, 0.7816427, 0.7914601, 0.8023700, 0.8137512, 
                                                     0.8255861, 0.8375781, 0.8498354, 0.8621276, 0.8741685, 
                                                     0.8855410, 0.8962833, 0.9061268, 0.9156126, 0.9243357, 
                                                     0.9320726, 0.9392022, 0.9459897, 0.9519587, 0.9574180, 
                                                     0.9623468, 0.9667851, 0.9708308, 0.9744296, 0.9777980, 
                                                     0.9807231, 0.9835296, 0.9861040, 0.9885404, 0.9906504, 
                                                     0.9926209, 0.9943923, 0.9959998, 0.9974658, 0.9987678, 
                                                     1.0000000};


using namespace CLHEP;

EsPmtEffectPulseTool::EsPmtEffectPulseTool(const std::string& type,
                                           const std::string& name, 
                                           const IInterface* parent)
    : GaudiTool(type,name,parent)
    , m_cableSvc(0)
    , m_simDataSvc(0)
    , m_simTimeEarliest(0)
    , m_simTimeLatest(0)
{
  // Initialization
  m_prePulsePdf.clear();     m_prePulsePdf.push_back(1);
  m_prePulseEdges.clear();   m_prePulseEdges.push_back(-50.0*ns);  m_prePulseEdges.push_back(0.0*ns);
  //  m_afterPulsePdf.clear();   m_afterPulsePdf.push_back(1);
  //  m_afterPulseEdges.clear(); m_afterPulseEdges.push_back(0.0*ns); m_afterPulseEdges.push_back(10.0e3*ns);
  m_afterPulsePdf.clear();
  m_afterPulseEdges.clear();

  for(int ii=0; ii < NUM_BINS_DIST; ii++) {
    m_afterPulsePdf.push_back(afterPusleTimingDist[ii]);
    m_afterPulseEdges.push_back(ii*100.+500.);
  }

  debug() << "time(ns) Afterpulse PDF (integrated) " <<endreq;
  for(int numLine=0; numLine < NUM_BINS_DIST; numLine++)
    debug() << m_afterPulseEdges[numLine] << ",   " << m_afterPulsePdf[numLine] << endreq;
  //    std::cout << m_afterPulseEdges[numLine] << ",   " << m_afterPulsePdf[numLine] << std::endl;

  declareInterface< IEsPulseTool >(this) ;   
  declareProperty("CableSvcName",m_cableSvcName="StaticCableSvc",
                  "Name of service to map between detector, hardware, and electronic IDs");
  declareProperty("SimDataSvcName",m_simDataSvcName="StaticSimDataSvc",
                  "Name of service to provide pmt properties for simulation");
  declareProperty("EnablePrePulse",m_enablePrePulse=true,"Switch on/off prepulse simulation");
  declareProperty("EnableAfterPulse",m_enableAfterPulse=true,"Switch on/off afterpulse simulation");
  declareProperty("PrePulsePdf",m_prePulsePdf,"User-defined PDF for timing distribution of pre-pulse");
  declareProperty("PrePulsePdfEdges",m_prePulseEdges,"Bin edges of PrePulsePdf"); 
  declareProperty("AfterPulseAmpMode",m_afterPulseAmpMode="SinglePE","Mode for afterpulse amplitude distribution");
  declareProperty("DarkPulseAmpMode",m_darkPulseAmpMode="SinglePE","Mode for dark pulse amplitude distribution");
  declareProperty("AfterPulsePdf",m_afterPulsePdf,"User-defined PDF for timing distribution of after-pulse");
  declareProperty("AfterPulsePdfEdges",m_afterPulseEdges,"Bin edges of AfterPulsePdf");
  declareProperty("AfterPulseTimeInterval", m_timeInterval=20*ns, 
                  "PMT hit counting window size");
  declareProperty("EnableNonlinearAfterpulse", m_enableNonlinearAfterpulse=true,
                  "Enable nonlinear suppression of afterpulsing.");
  declareProperty("LinearAfterpulseThreshold", m_linearAfterpulseThreshold=400,
                  "Upper limit of linear afterpulsing in number of PE.");


  declareProperty("EnableVeryLongTimeSuppression",m_enableVeryLongTimeSuppression=false,
                  "Enable suppression of hit times in the far future");
  declareProperty("VeryLongTime",m_veryLongTime=3e7*s,
                  "Definition of very long time for hit time suppression");
  declareProperty("ExpWeight", m_expWeight=0,"Weight of the exponential contribution to the spe response function");
  declareProperty("SpeExpDecay", m_speExpDecay=0.5,"Decay time of the exponential contribution to the spe response function");
}

EsPmtEffectPulseTool::~EsPmtEffectPulseTool(){}

StatusCode EsPmtEffectPulseTool::initialize()
{
  // Get Cable Service
  m_cableSvc = svc<ICableSvc>(m_cableSvcName,true);

  // Get PMT simulation input data service
  m_simDataSvc = svc<ISimDataSvc>(m_simDataSvcName,true);

  // Random number service
  IRndmGenSvc *rgs = 0;
  if (service("RndmGenSvc",rgs,true).isFailure()) {
    fatal() << "Failed to get random service" << endreq;
    return StatusCode::FAILURE;        
  }
  
  StatusCode sc;
  sc = m_uni.initialize(rgs, Rndm::Flat(0,1));
  if (sc.isFailure()) {
    fatal() << "Failed to initialize uniform random numbers" << endreq;
    return StatusCode::FAILURE;
  }
  sc = m_gauss.initialize(rgs, Rndm::Gauss(0,1));
  if (sc.isFailure()) {
    fatal() << "Failed to initialize gaussian random numbers" << endreq;
    return StatusCode::FAILURE;
  }
  sc = m_exp.initialize(rgs, Rndm::Exponential(m_speExpDecay));
  if (sc.isFailure()) {
    fatal() << "Failed to initialize exponential random numbers" << endreq;
    return StatusCode::FAILURE;
  }
  if(m_enablePrePulse){
    info() << "Prepulse simulation enabled" << endreq;
    sc = m_randPrePulseTime.initialize(rgs, Rndm::DefinedPdf(m_prePulsePdf,0));
    if (sc.isFailure()) {
      fatal() << "Failed to initialize random numbers for timing distribution" << endreq;
      return StatusCode::FAILURE;
    }
  }
  else info() << "Prepulse simulation disabled" << endreq;
  
  if(m_enableAfterPulse){
    info() << "Afterpulse simulation enabled" << endreq;
    sc = m_randAfterPulseTime.initialize(rgs, Rndm::DefinedPdf(m_afterPulsePdf,0));
    if (sc.isFailure()) {
      fatal() << "Failed to initialize random numbers for timing distribution" << endreq;
      return StatusCode::FAILURE;
    }
  }
  else info() << "Afterpulse simulation disabled" << endreq;
  
  // Check user input for user-defined PDFs
  if(m_enablePrePulse){
    if (m_prePulsePdf.size()+1   != m_prePulseEdges.size()) { 
        // ||   m_afterPulsePdf.size()+1 != m_afterPulseEdges.size() ) {
      fatal() << "There should be 1 more number of bin edges than bins when defining pdf" << endreq;
      return StatusCode::FAILURE; 
    }
  }

  // DWL: Check user input for user-defined PMT hit counting time interval
  if (m_timeInterval <= 0) {
    fatal() << "PMT hit couting time interval should be great than zero !" << endreq;
    return StatusCode::FAILURE;
  }
  if(m_enableAfterPulse){
    if ((m_afterPulseAmpMode != "SinglePE") && (m_afterPulseAmpMode != "PDF")) {
      fatal() << "Not a recognized afterpulse amplitude distribution mode" << endreq;
      return StatusCode::FAILURE;
    }
    if (m_afterPulseAmpMode == "SinglePE") debug() << "Assume single p.e. amplitude for afterpulses" << endreq; 
    if (m_afterPulseAmpMode == "PDF" || m_darkPulseAmpMode == "PDF") {
      debug() << "Set up amplitude PDF for afterpulses" << endreq; 
      m_afterPulseAmpPdf.clear();
      m_afterPulseAmpEdges.clear();
      getAfterPulseAmpPdf(m_afterPulseAmpPdf);
      getAfterPulseAmpEdges(m_afterPulseAmpEdges);
    }  
  }
  return StatusCode::SUCCESS;
}

StatusCode EsPmtEffectPulseTool::finalize()
{
  return StatusCode::SUCCESS;
}

StatusCode EsPmtEffectPulseTool::generatePulses(DayaBay::SimHitCollection* hits,
                                                DayaBay::ElecPulseCollection* pulses)
{    
  debug() << "running generatePulses()" << endreq; 

  // Define simulation time window
  m_simTimeEarliest = pulses->header()->header()->earliest();
  m_simTimeLatest   = pulses->header()->header()->latest();
  TimeStamp headerTime = pulses->header()->header()->timeStamp();

  // Maintain a list of primary pulses for each PMT
  std::map<DayaBay::DetectorSensor, 
           std::vector<DayaBay::ElecPulse*> > pulsesByPmt;
  std::map<DayaBay::DetectorSensor, 
           std::vector<DayaBay::ElecPulse*> >::iterator pulsePmtIter, 
           pulsePmtEnd; 

  const DayaBay::SimHitCollection::hit_container& htcon = hits->collection();
  DayaBay::SimHitCollection::hit_container::const_iterator hcIter, 
    hcDone = htcon.end();
  DayaBay::ElecPulseCollection::PulseContainer& pulsecon = pulses->pulses();

  // Context, ServiceMode for this data
  Context context(pulses->detector().site(), SimFlag::kMC, 
                  pulses->header()->header()->timeStamp(),
                  pulses->detector().detectorId());
  int task = 0;
  ServiceMode svcMode(context, task);

  debug() << "Processing " << htcon.size() << " simulated hits from detector " 
          << hits->detector() << endreq; 


  for (hcIter=htcon.begin(); hcIter != hcDone; ++hcIter) {

    DayaBay::DetectorSensor sensDetId((*hcIter)->sensDetId());

    if (sensDetId.isAD()) {
        DayaBay::AdPmtSensor pmtid(sensDetId.fullPackedData());
        if (pmtid.bogus()) {
            warning() << "Got bogus AD PMT: " << pmtid << endreq;
        }
    }
    else {
        DayaBay::PoolPmtSensor pmtid(sensDetId.fullPackedData());
        if (pmtid.bogus()) {
            warning() << "Got bogus Pool PMT: " << pmtid << endreq;
        }
    }

    //temporarily ignore very long time hit to avoid program break
    //if((*hcIter)->hitTime() > 1e7) {
    if((*hcIter)->hitTime() > m_veryLongTime ) {
        warning() << "Skipping very long hit time of " << (*hcIter)->hitTime() 
                  << " in PMT " << sensDetId << endreq;
        continue;
    }

    verbose() << "Requesting pmtData for " << sensDetId << endreq;
    const DayaBay::PmtSimData* pmtData = 
      m_simDataSvc->pmtSimData(sensDetId, svcMode);
    if(!pmtData){
      warning() << "No Simulation input properties for PMT: " << sensDetId 
                << " Ignoring..." << endreq;
      continue;
    }
    // Ignore SimHit due to efficiency?
    if (m_uni() > pmtData->m_efficiency) continue;


    const DayaBay::ElecChannelId& elecChannelId = m_cableSvc->elecChannelId( 
                                                                     sensDetId,
                                                                     svcMode);
    if(elecChannelId.fullPackedData() == 0){
      warning() << "PMT: " << sensDetId 
                << " is not connected in cable map.  Ignoring..." << endreq;
      continue;
    }
    DayaBay::ElecPulse* simpulse = new DayaBay::ElecPmtPulse();

    // Primary vertex time (absolute, in seconds)
    TimeStamp vertexTime = (*hcIter)->hc()->header()->header()->timeStamp();
    // Need to subtract ElecHeader time (absolute, in seconds) to get 
    // pulse time wrt ElecHeader
    // Set pulse time, include transit time (nanoseconds)
    double transitTime = m_gauss() * pmtData->m_timeSpread 
      + pmtData->m_timeOffset;
    double headerOffset = double(vertexTime - headerTime)*Gaudi::Units::second;
    double pulseHitTime = headerOffset + (*hcIter)->hitTime() + transitTime;

    simpulse->setTime(pulseHitTime);

    simpulse->setAmplitude(PulseAmp((*hcIter)->weight(),pmtData->m_gain,
                                    pmtData->m_sigmaGain));  // Include SPE effects, relative gain
    simpulse->setAncestor((*hcIter));
    simpulse->setType(PulseType::kPmtHit);
    simpulse->setChannelId(elecChannelId);

    pulsecon.push_back(simpulse);

    // PMT pulse counting
    if(m_enableAfterPulse){
      if( m_enableNonlinearAfterpulse ){
        // Add pulse to the list for this PMT
        pulsesByPmt[sensDetId].push_back(simpulse);
      }else{
        // Immediatly add linear afterpulsing
        if (m_uni() < pmtData->m_afterPulseProb) 
          pulsecon.push_back(makeAfterPulse(simpulse, *pmtData));
      }
    }
    // Include pre-pulse?
    if (m_uni() < pmtData->m_prePulseProb && m_enablePrePulse){
      pulsecon.push_back(makePrePulse(simpulse, *pmtData));
    }
  } 

  // Add afterpulses in nonlinear mode
  if(m_enableAfterPulse){
    if( m_enableNonlinearAfterpulse){
      pulsePmtEnd = pulsesByPmt.end(); 
      for(pulsePmtIter = pulsesByPmt.begin();pulsePmtIter != pulsePmtEnd;pulsePmtIter++){
        const DayaBay::DetectorSensor& sensDetId = pulsePmtIter->first;
        std::vector<DayaBay::ElecPulse*>& pulseList = pulsePmtIter->second;
        // Get PMT properties
        const DayaBay::PmtSimData* pmtData = 
          m_simDataSvc->pmtSimData(sensDetId, svcMode);
        if(!pmtData){
          error() << "Nonlinear: No Simulation input properties for PMT: "
                  << sensDetId << endreq;
          return StatusCode::FAILURE;
        }
  
        // Prepare a handle for memory buffer
        int* pulseCount = 0;
  
        // Process hits to make afterpulses
        std::vector<DayaBay::ElecPulse*>::iterator pulseIter, 
          pulseEnd = pulseList.end();
        if(pulseList.size() > m_linearAfterpulseThreshold){
          // Add Nonlinear afterpulsing
          // Step 1: assemble pulse counts in a specified time window
          int numTimeSlots = double(m_simTimeLatest - headerTime)
            *Gaudi::Units::second/m_timeInterval + 1;
          // allocate memory only if needed
          if(!pulseCount) pulseCount = new int[numTimeSlots];
          // Clear memory buffer
          for(int i=0; i<numTimeSlots; i++) pulseCount[i] = 0;
          // Loop once to create pulse counts
          for(pulseIter=pulseList.begin(); pulseIter != pulseEnd; pulseIter++){
            int hitTimeSlot = int((*pulseIter)->time()/m_timeInterval);
            int pulseNo = int((*pulseIter)->amplitude()+0.5);
            pulseCount[hitTimeSlot] += pulseNo;
          }
          // Loop again to create afterpulses
          for(pulseIter=pulseList.begin(); pulseIter != pulseEnd; pulseIter++){
            int hitTimeSlot = int((*pulseIter)->time()/m_timeInterval);
            double nonlinearProb = pmtData->m_afterPulseProb/0.0164
              *NumAfterPulse(pulseCount[hitTimeSlot]) / pulseCount[hitTimeSlot];
            if (m_uni() < nonlinearProb) 
              pulsecon.push_back(makeAfterPulse(*pulseIter, *pmtData));
          }
        }else{
          // Add Linear afterpulsing
          for(pulseIter=pulseList.begin(); pulseIter != pulseEnd; pulseIter++){
            if (m_uni() < pmtData->m_afterPulseProb) 
              pulsecon.push_back(makeAfterPulse(*pulseIter, *pmtData));
          }
        }
        // Free memory buffer if it was used
        if(pulseCount) delete [] pulseCount;
      }
    }
  }
  // Process dark pulses
  // Get list of all detector sensors
  debug() << "Adding dark hits" << endreq; 
  
  const std::vector<DayaBay::DetectorSensor>& detSensors = 
    m_cableSvc->sensors( svcMode );

  for (size_t idet = 0; idet < detSensors.size(); ++idet) {
      DayaBay::DetectorSensor detSens = detSensors[idet];      
      if (detSens.isAD() ) {
          DayaBay::AdPmtSensor pmtid(detSens.fullPackedData());
          if (pmtid.bogus()) {
              warning() << "bogus AD PMT id: " << pmtid << endreq;
          }
      }
      else {
          DayaBay::PoolPmtSensor pmtid(detSens.fullPackedData());
          if (pmtid.bogus()) {
              warning() << "bogus Pool PMT id: " << pmtid << endreq;
          }
      }
  }


  for (unsigned int idet=0; idet < detSensors.size(); ++idet) {
    // Skip this detector sensor if not in detector of ElecPulseCollection/SimHitCollection
    DayaBay::DetectorSensor detSens = detSensors[idet];

    if (detSens.detName() != pulses->detector().detName()) continue;

    const DayaBay::ElecChannelId& channelId = m_cableSvc->elecChannelId( 
                                                                 detSens,
                                                                 svcMode);

    const DayaBay::PmtSimData* pmtData = 
      m_simDataSvc->pmtSimData(detSens, svcMode);
    if(!pmtData){
        error() << "No Simulation input properties for PMT #" << idet
                << ": id = " << detSens
                << " svcMode = " << svcMode
                << endreq;
        return StatusCode::FAILURE;
    }
    // Generate Poisson-distributed number around mean number of dark hits in simulation time window
    int Ndark = PoissonRand(pmtData->m_darkRate
                            *double(m_simTimeLatest-m_simTimeEarliest)
                            *Gaudi::Units::second);
    verbose() << "Adding " << Ndark << " dark hits on pmt " 
            << detSens << " with dark rate " << pmtData->m_darkRate
            << " and dt " 
            << double(m_simTimeLatest-m_simTimeEarliest)*Gaudi::Units::second 
            << endreq;
    for (int dummy = 0; dummy < Ndark; ++dummy) {
      pulsecon.push_back(makeDarkPulse(channelId, *pmtData));
    }
  }
  
  debug() << "Adding " << pulsecon.size() << " pulses to detector " 
          << pulses->detector() << endreq; 
  return StatusCode::SUCCESS;
}

DayaBay::ElecPulse* EsPmtEffectPulseTool::makePrePulse(
                                          DayaBay::ElecPulse* simpulse,
                                          const DayaBay::PmtSimData& pmtData) {
  DayaBay::ElecPulse* prepulse = new DayaBay::ElecPmtPulse();

  // Time offset from main pulse based on time PDF of pre-pulses
  double current_rand = m_randPrePulseTime();
  double prePulseTime = ConvertPdfRand(current_rand,m_prePulsePdf,m_prePulseEdges);

  prepulse->setTime(simpulse->time() + prePulseTime);
  prepulse->setAmplitude(PulseAmp(simpulse->amplitude(), pmtData.m_gain,
                                  pmtData.m_sigmaGain));
  prepulse->setAncestor(simpulse->ancestor());
  prepulse->setType(PulseType::kPrePulse);
  prepulse->setChannelId(simpulse->channelId());
  
  return prepulse;
}

DayaBay::ElecPulse* EsPmtEffectPulseTool::makeAfterPulse(
                                          DayaBay::ElecPulse* simpulse,
                                          const DayaBay::PmtSimData& pmtData) {
  DayaBay::ElecPulse* afterpulse = new DayaBay::ElecPmtPulse();

  // Time offset from main pulse based on time PDF of after-pulses
  // double current_rand = m_randAfterPulseTime();
  double current_rand = m_uni();
  double afterPulseTime = ConvertPdfRand01(current_rand,m_afterPulsePdf,m_afterPulseEdges);
  //  std::cout << "current_rand = " << current_rand  << ", time = " << afterPulseTime << std::endl;

  afterpulse->setTime(simpulse->time() + afterPulseTime); 
  if(m_afterPulseAmpMode == "SinglePE") 
    afterpulse->setAmplitude(PulseAmp(simpulse->amplitude(), pmtData.m_gain, pmtData.m_sigmaGain));
  if(m_afterPulseAmpMode == "PDF"){
    current_rand = m_uni();
    double amplitude = ConvertPdfRand01(current_rand, m_afterPulseAmpPdf, m_afterPulseAmpEdges);
    afterpulse->setAmplitude(PulseAmp(amplitude, pmtData.m_gain, pmtData.m_sigmaGain));
  }
  afterpulse->setAncestor(simpulse->ancestor()); 
  afterpulse->setType(PulseType::kAfterPulse);
  afterpulse->setChannelId(simpulse->channelId());
  
  return afterpulse;
}

DayaBay::ElecPulse* EsPmtEffectPulseTool::makeDarkPulse(
                                        const DayaBay::ElecChannelId& channelId,
                                        const DayaBay::PmtSimData& pmtData) {
  DayaBay::ElecPulse* darkpulse = new DayaBay::ElecPmtPulse();
  double darkPulseTime
    = m_uni()*double(m_simTimeLatest-m_simTimeEarliest)*Gaudi::Units::second;

  darkpulse->setTime(darkPulseTime);
  if(m_darkPulseAmpMode == "SinglePE") darkpulse->setAmplitude(PulseAmp(1.0, pmtData.m_gain, pmtData.m_sigmaGain));
  if(m_darkPulseAmpMode == "PDF"){
    double current_rand = m_uni();
    double amplitude = ConvertPdfRand01(current_rand, m_afterPulseAmpPdf, m_afterPulseAmpEdges);
    darkpulse->setAmplitude(PulseAmp(amplitude, pmtData.m_gain, pmtData.m_sigmaGain));
  }
  darkpulse->setAncestor(0);
  darkpulse->setType(PulseType::kDarkPulse);
  darkpulse->setChannelId(channelId);
  
  return darkpulse;
}

float EsPmtEffectPulseTool::PulseAmp(float weight, float gain, float sigmaGain){

  // Include relative gain
  float randW = m_uni();
  if (randW > m_expWeight ){
    return m_gauss() * sigmaGain + gain;
  }
  else {
    return (m_exp() + m_speCutoff) * gain;
  }
}


double EsPmtEffectPulseTool::ConvertPdfRand(double rand, std::vector<double> pdf, std::vector<double> edges) {
  // Defined PDF returns random number in [0,1] distributed according to user-defined histogram.
  // It assumes even bin sizes, so accomodate uneven bin sizes for generality.
  int current_bin;
  int Nbins = pdf.size();

  for (int bin=0;bin<Nbins;bin++) { // find which bin rand is in
    if (rand >= (double)(bin/Nbins) && rand < (double)((bin+1)/Nbins)) current_bin = bin;
    else current_bin = Nbins-1;     // else rand=1, so current_bin is last bin (Nbins-1)
  }
  return edges[current_bin] + (rand*Nbins - current_bin)*(edges[current_bin+1]-edges[current_bin]);
}

double EsPmtEffectPulseTool::ConvertPdfRand01(double rand, std::vector<double> pdf, std::vector<double> edges) {
  // Defined PDF returns random number in [0,1] distributed according to user-defined histogram.
  // It assumes even bin sizes, so accomodate uneven bin sizes for generality.
  int current_bin;
  int Nbins = pdf.size();

  for(int bin=0; bin<Nbins; bin++) {
    if(rand >= pdf[bin] && rand < pdf[bin+1]) {
      current_bin = bin;
      break;
    }
    else
      current_bin = Nbins-1;
  }

  return edges[current_bin] + (rand-pdf[current_bin])*(edges[current_bin+1]-edges[current_bin])
    /(pdf[current_bin+1]-pdf[current_bin]);
}

int EsPmtEffectPulseTool::PoissonRand(double mean) {
  // Using source code from ROOT's TRandom::Poisson
  // Note: ROOT uses different algorithms depending on the mean, but mean is small 
  //       for our purposes, so use algorithm for mean<25
 
  int n;
  if (mean <= 0) return 0;

  double expmean = exp(-mean);
  double pir = 1;
  n = -1;
  while(1) {
    n++;
    pir *= m_uni();
    if (pir <= expmean) break;
  }
  return n;
}


double EsPmtEffectPulseTool::NumAfterPulse(const int numPmtHit) {
  // Get number of after-pulses per PMT pulse based on the number of PMT hits

  double cnt;

  if (numPmtHit <= 400)
    cnt = 0.016 * numPmtHit;
  else if (numPmtHit <= 2500)
    cnt = -1.307853 + 0.02666278 * numPmtHit - 0.2001321e-4 * pow(numPmtHit, 2)
               +0.7330343e-8 * pow(numPmtHit, 3) - 0.102098e-11 * pow(numPmtHit, 4);
  else
    cnt = 15.;

  return cnt;
}

void EsPmtEffectPulseTool::getAfterPulseAmpEdges(std::vector<double>& edges){
  edges.push_back(0.5);edges.push_back(0.7);edges.push_back(0.9);edges.push_back(1.1);
  edges.push_back(1.3);edges.push_back(1.5);edges.push_back(1.7);edges.push_back(1.9);
  edges.push_back(2.1);edges.push_back(2.3);edges.push_back(2.5);edges.push_back(2.7);
  edges.push_back(3.05);edges.push_back(3.4);edges.push_back(3.75);edges.push_back(4.1);
  edges.push_back(4.6);edges.push_back(5.1);edges.push_back(5.6);edges.push_back(6.1);
  edges.push_back(7.1);edges.push_back(8.1);edges.push_back(9.1);edges.push_back(10.1);
  edges.push_back(11.1);edges.push_back(13.1);edges.push_back(15.1);edges.push_back(17.1);
  edges.push_back(19.1);edges.push_back(21.1);edges.push_back(23.1);edges.push_back(25.1);
  edges.push_back(27.1);edges.push_back(29.1);edges.push_back(31.1);edges.push_back(35.1);
  edges.push_back(39.1);edges.push_back(43.1);edges.push_back(47.1);edges.push_back(51.1);
  edges.push_back(55.1);edges.push_back(59.1);edges.push_back(63.1);edges.push_back(67.1);
  edges.push_back(71.1);edges.push_back(75.1);edges.push_back(79.1);edges.push_back(83.1);
  edges.push_back(87.1);edges.push_back(91.1);edges.push_back(95.1);edges.push_back(100);
}

void EsPmtEffectPulseTool::getAfterPulseAmpPdf(std::vector<double>& pdf){
  pdf.push_back(0);pdf.push_back(0.0219574);pdf.push_back(0.0931247);pdf.push_back(0.179757);
  pdf.push_back(0.264803);pdf.push_back(0.342568);pdf.push_back(0.411712);pdf.push_back(0.472498);
  pdf.push_back(0.525729);pdf.push_back(0.572335);pdf.push_back(0.613205);pdf.push_back(0.649137);
  pdf.push_back(0.702162);pdf.push_back(0.745131);pdf.push_back(0.780295);pdf.push_back(0.809341);
  pdf.push_back(0.842677);pdf.push_back(0.868769);pdf.push_back(0.889482);pdf.push_back(0.906131);
  pdf.push_back(0.930777);pdf.push_back(0.947671);pdf.push_back(0.95962);pdf.push_back(0.968294);
  pdf.push_back(0.974731);pdf.push_back(0.983341);pdf.push_back(0.988564);pdf.push_back(0.99189);
  pdf.push_back(0.994094);pdf.push_back(0.995601);pdf.push_back(0.996661);pdf.push_back(0.997423);
  pdf.push_back(0.997983);pdf.push_back(0.998401);pdf.push_back(0.998717);pdf.push_back(0.999151);
  pdf.push_back(0.999418);pdf.push_back(0.999591);pdf.push_back(0.999705);pdf.push_back(0.999783);
  pdf.push_back(0.999838);pdf.push_back(0.999876);pdf.push_back(0.999905);pdf.push_back(0.999926);
  pdf.push_back(0.999941);pdf.push_back(0.999953);pdf.push_back(0.999962);pdf.push_back(0.999969);
  pdf.push_back(0.999975);pdf.push_back(0.999979);pdf.push_back(0.999983);pdf.push_back(1.0);
}

---