ExcessNoiseMaker.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \brief   Code for simulating the spread in photo-electrons 
///          according to the theory of APDs
/// \author  aurisaam@ucmail.uc.edu
/// \date
////////////////////////////////////////////////////////////////////////
#include "ReadoutSim/ExcessNoiseMaker.h"

#include "NovaDAQConventions/DAQConventions.h"

#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Services/Optional/RandomNumberGenerator.h"
#include "CLHEP/Random/RandGaussQ.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "TMath.h"
#include <cmath>
#include <iostream>

namespace rsim{

  ExcessNoiseMaker::ExcessNoiseMaker(const fhicl::ParameterSet& pset) :
    CommonParameters(pset)
  {
    fK = (fAPDExcessNoiseFactor + 1.0/fGain - 2)/(fGain + 1.0/fGain - 2);
  }
  
  double ExcessNoiseMaker::ExcessNoisePDF(int nPE_in, double nPE_out)
  {
    //number of electrons                                                                                
    int m = int(nPE_out*fGain);
    if (nPE_in < 0 || m < 0) return 0.0;
    if (nPE_in == 0)
      {
        if (m == 0) return 1.0;
        else return 0.0;
      }
    int r = m - nPE_in;
    if (r < 0) return 0.0;
    
    double logProb = log(nPE_in);
    logProb += r*log( (1-fK)*(fGain-1)/fGain );
    logProb += ((nPE_in+fK*r)/(1-fK))*log( (1+fK*(fGain-1))/fGain );
    logProb += lgamma( (nPE_in+r)/(1-fK) );
    logProb -= log( nPE_in+fK*r);
    logProb -= lgamma( r + 1 );
    logProb -= lgamma( (nPE_in+fK*r)/(1-fK) );
    
    return fGain*exp( logProb );
  }
  
  double ExcessNoiseMaker::DrawSmearedPE(int nPE, CLHEP::RandFlat* flat)
  {
    if (nPE == 0)
      {
        return 0.0;
      }
    
    std::map<int, std::vector<double> >::iterator it = fPDFs.find(nPE);
    if (it == fPDFs.end()) it = MakeTemplate(nPE);
    
    const Double_t u(flat->fire());
    const double sigma = std::sqrt(fAPDExcessNoiseFactor*nPE);
    const double xlo = (nPE - 7*sigma <0) ? 0.0 : nPE - 7*sigma;
    const double xhi = nPE + 7*sigma;
    const double binw(1.0/fGain);
    const int nbins((xhi-xlo)/binw);
    //const Double_t binw(1.0/fGain);
    //const Int_t nbins(10.0*nPE/binw);
    const Int_t ibin = TMath::BinarySearch(nbins, &(it->second)[0], u);
    Double_t x = xlo + ibin*binw; //low edge                                                                   
    if (u > it->second[ibin]) x += binw*(u - it->second[ibin])/(it->second[ibin+1] - it->second[ibin]);
    return x;
  }
  
  std::map<int, std::vector<double> >::iterator ExcessNoiseMaker::MakeTemplate(int nPE)
  {
    if (nPE == 0) {
      std::cerr << "nPE should never be zero in GetSmearedPE!" << std::endl;
      exit(1);
    }

    //std::cout << "MakeTemplate(" << nPE << ")" << std::endl;
    
    const double sigma = std::sqrt(fAPDExcessNoiseFactor*nPE);
    const double xlo = (nPE - 7*sigma < 0) ? 0.0 : nPE - 7*sigma;
    const double xhi = nPE + 7*sigma;
    const double binw(1.0/fGain);
    const int nbins((xhi-xlo)/binw);
    
    std::vector<double> pdf(nbins);
    for (int i = 0; i < nbins; ++i) {
      double nPE_out = xlo + i*binw;
      if (i == 0) pdf[i] = ExcessNoisePDF(nPE, nPE_out);
      else pdf[i] = pdf[i-1] + ExcessNoisePDF(nPE,nPE_out);
    }

    for (int i = 0; i < nbins; ++i) {
      pdf[i] /= pdf[nbins-1];
    }
    
    std::pair<std::map<int, std::vector<double> >::iterator, bool> result = fPDFs.insert( make_pair(nPE, pdf) );
    
    return result.first;
    
  }

}