ExcessNoiseMaker.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief Code for simulating the spread in photo-electrons
3 /// according to the theory of APDs
4 /// \author aurisaam@ucmail.uc.edu
5 /// \date
6 ////////////////////////////////////////////////////////////////////////
8 
9 #include "NovaDAQConventions/DAQConventions.h"
10 
15 
16 #include "TMath.h"
17 #include <cmath>
18 #include <iostream>
19 
20 namespace rsim{
21 
23  CommonParameters(pset)
24  {
25  fK = (fAPDExcessNoiseFactor + 1.0/fGain - 2)/(fGain + 1.0/fGain - 2);
26  }
27 
28  double ExcessNoiseMaker::ExcessNoisePDF(int nPE_in, double nPE_out)
29  {
30  //number of electrons
31  int m = int(nPE_out*fGain);
32  if (nPE_in < 0 || m < 0) return 0.0;
33  if (nPE_in == 0)
34  {
35  if (m == 0) return 1.0;
36  else return 0.0;
37  }
38  int r = m - nPE_in;
39  if (r < 0) return 0.0;
40 
41  double logProb = log(nPE_in);
42  logProb += r*log( (1-fK)*(fGain-1)/fGain );
43  logProb += ((nPE_in+fK*r)/(1-fK))*log( (1+fK*(fGain-1))/fGain );
44  logProb += lgamma( (nPE_in+r)/(1-fK) );
45  logProb -= log( nPE_in+fK*r);
46  logProb -= lgamma( r + 1 );
47  logProb -= lgamma( (nPE_in+fK*r)/(1-fK) );
48 
49  return fGain*exp( logProb );
50  }
51 
53  {
54  if (nPE == 0)
55  {
56  return 0.0;
57  }
58 
59  std::map<int, std::vector<double> >::iterator it = fPDFs.find(nPE);
60  if (it == fPDFs.end()) it = MakeTemplate(nPE);
61 
62  const Double_t u(flat->fire());
63  const double sigma = std::sqrt(fAPDExcessNoiseFactor*nPE);
64  const double xlo = (nPE - 7*sigma <0) ? 0.0 : nPE - 7*sigma;
65  const double xhi = nPE + 7*sigma;
66  const double binw(1.0/fGain);
67  const int nbins((xhi-xlo)/binw);
68  //const Double_t binw(1.0/fGain);
69  //const Int_t nbins(10.0*nPE/binw);
70  const Int_t ibin = TMath::BinarySearch(nbins, &(it->second)[0], u);
71  Double_t x = xlo + ibin*binw; //low edge
72  if (u > it->second[ibin]) x += binw*(u - it->second[ibin])/(it->second[ibin+1] - it->second[ibin]);
73  return x;
74  }
75 
76  std::map<int, std::vector<double> >::iterator ExcessNoiseMaker::MakeTemplate(int nPE)
77  {
78  if (nPE == 0) {
79  std::cerr << "nPE should never be zero in GetSmearedPE!" << std::endl;
80  exit(1);
81  }
82 
83  //std::cout << "MakeTemplate(" << nPE << ")" << std::endl;
84 
85  const double sigma = std::sqrt(fAPDExcessNoiseFactor*nPE);
86  const double xlo = (nPE - 7*sigma < 0) ? 0.0 : nPE - 7*sigma;
87  const double xhi = nPE + 7*sigma;
88  const double binw(1.0/fGain);
89  const int nbins((xhi-xlo)/binw);
90 
91  std::vector<double> pdf(nbins);
92  for (int i = 0; i < nbins; ++i) {
93  double nPE_out = xlo + i*binw;
94  if (i == 0) pdf[i] = ExcessNoisePDF(nPE, nPE_out);
95  else pdf[i] = pdf[i-1] + ExcessNoisePDF(nPE,nPE_out);
96  }
97 
98  for (int i = 0; i < nbins; ++i) {
99  pdf[i] /= pdf[nbins-1];
100  }
101 
102  std::pair<std::map<int, std::vector<double> >::iterator, bool> result = fPDFs.insert( make_pair(nPE, pdf) );
103 
104  return result.first;
105 
106  }
107 
108 }
double fAPDExcessNoiseFactor
APD&#39;s "excess noise factor" (see .cxx for more)
double DrawSmearedPE(int nPE, CLHEP::RandFlat *flat)
set< int >::iterator it
fvar< T > lgamma(const fvar< T > &x)
Definition: lgamma.hpp:20
T sqrt(T number)
Definition: d0nt_math.hpp:156
double ExcessNoisePDF(int nPE_in, double nPE_out)
OStream cerr
Definition: OStream.cxx:7
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
Common configuration params for SimpleReadout, FPGAAlgorithms, NoiseMaker.
const int nbins
Definition: cellShifts.C:15
ExcessNoiseMaker(const fhicl::ParameterSet &pset)
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
double fGain
APD gain (electrons per photoelectron)
double sigma(TH1F *hist, double percentile)
Can be used as either a member holding configurations, or a mix-in.
std::map< int, std::vector< double > >::iterator MakeTemplate(int nPE)
exit(0)
std::map< int, std::vector< double > > fPDFs
TRandom3 r(0)