c hpd - a package for simulating the readout preamplifier electronics
c       for the hybrid photodiode (HPD)
c
c Richard Jones
c March 6, 2005
c
c Usage:  To generate one signal and 10 noise pulses do:
c            paw> call hpd.f77(100,10)
c         The output signal is saved in histogram 100 and histograms 101-110
c         will contain the noise traces.
c
c Model 1: fast readout with good timing and noise of 6 p.e.
c  This circuit couples the signal from the HPD (capacitance C0=200pF)
c  to the emitter of a fast bipolar transistor.  The bias current of the
c  transistor is set by a 500 Ohm resisitor and a 50 Ohm base network so
c  the effective discharge resistance is low (30mV / Icc) with Icc=1mA.
c  The collector current drives a R0=3KOhm gain resistor.  A fast op-amp
c  (specs taken from OPA657) reads the voltage across R0.  The feedback
c  circuit between the op-amp output and its negative input is a resistor
c  capacitor pair (R=5KOhm, C=1pF), with a reference resisitor Rt=100Ohms
c  between the negative pin and ground.  The output voltage for a 100 p.e.
c  pulse is about 350mV peak so no further amplification is needed before
c  digitization.  The following noise sources are included.
c   1. Johnson noise associated with the capacitance of the HPD
c   2. shot noise associated with the collector current of the transistor
c   3. input noise on the op-amp (as specified on the OPA657 data sheet)
c  The results give an integrated charge of 100 ± 6 p.e. and a leading-
c  edge timing resolution of 250ps r.m.s.  To use this simulation,
c  set the noise function to compute the sum of johnson_noise, shot_noise
c  and opamp_noise1, and set the hpd subroutine to call the filter1.  If
c  you change any component values, make sure that you change them in both
c  the noise_* and filter1 functions.
c
c Model 2: slow readout with poor timing and noise of 4 p.e.
c  This circuit couples the HPD (capacitance C0=200pF) directly to the
c  minus input pin of the op-amp.  The feedback of the op-amp is a RC
c  parallel circuit with R=50KOhms and C=2pF.  The positive input is
c  coupled directly to ground.  This is the configuration described in
c  the OPA657 data sheet, and turns out to be very slow because of the
c  large capacitance of the HPD.  The following noise sources are included.
c   1. Johnson noise associated with the capacitance of the HPD
c   2. input noise on the op-amp (as specified on the OPA657 data sheet)
c  The results give an integrated charge of 100 ± 6 p.e. and a leading-
c  edge timing resolution of 3ns r.m.s.  To use this simulation,
c  set the noise function to compute the sum of johnson_noise and
c  opamp_noise2, and set the hpd subroutine to call the filter2.  If
c  you change any component values, make sure that you change them in
c  both the noise_* and filter2 functions.

      subroutine hpd(id,n)
      integer id,n
      external pulse,pulse1,noise
      integer nbins
      real tstart,tend
c     parameter (nbins=1500,tstart=-50,tend=100)
c     call filter1(id,pulse,nbins,tstart,tend)
      parameter (nbins=1500,tstart=-100,tend=900)
      call filter2(id,pulse,nbins,tstart,tend)
      do i=1,n
c       call filter1(id+i,noise,nbins,tstart,tend)
        call filter2(id+i,noise,nbins,tstart,tend)
      enddo
      end

c  generic HPD pulse for 1 p.e. (assumes manufacturer's characteristics)
      real function pulse(t,dt)          ! Amperes as a function of ns
      real t,dt
      real pe,gain
      parameter (pe=1,gain=2700)
      real x
      x=t-10
      if (x.gt.0) then
         pulse=x**2*exp(-x/10-x**2/80)/125
      else
         pulse=0
      endif
      pulse=pulse*pe*gain*1.6e-19*1e9
      end

c  special shaped pulse coming from super-bright LED
c  [MRS Internet J. Nitride Semicond. Res. 5S1, W11.58 (2000)]
      real function pulse1(t,dt)          ! Amperes as a function of ns
      real t,dt
      real pe,gain
      parameter (pe=1,gain=2700)	  ! HPD characteristics
      real tau,beta
      parameter (tau=30,beta=0.45)	  ! LED decay characteristics
      real stime
      parameter (rtime=5)		  ! rising edge shaping time
      real x,b
      x=t-10
      if (x.gt.0) then
         b=1-(1-exp(-x/rtime))*(1-beta)
         pulse1=(exp(-(x/tau)**b)-exp(-x/rtime))/65
      else
         pulse1=0
      endif
      pulse1=pulse1*pe*gain*1.6e-19*1e9
      end

c  generic noise routine - customize to simulate either model 1 or 2
      real function noise(t,dt)		! Amperes as a function of ns
      real t,dt
      real johnson_noise,shot_noise,opamp_noise
      external johnson_noise,shot_noise,opamp_noise
      noise=johnson_noise(t,dt)
     +     +opamp_noise2(t,dt)
c     noise=johnson_noise(t,dt)
c    +     +shot_noise(t,dt)
c    +     +opamp_noise1(t,dt)
      end

c  noise current on the HPD from the equipartition theorem
      real function johnson_noise(t,dt)	! Amperes as a function of ns
      real t,dt
      real q,qsaved
      common /qsave1/qsaved
      data qsaved/0/
      real kT,C0
      parameter (kT=4.04e-21,C0=200e-12)
      real rangauss
      external rangauss
      if (qsaved.eq.0) then
        qsaved=rangauss()*sqrt(kT*C0)
      endif
      q=rangauss()*sqrt(kT*C0)
      johnson_noise=(q-qsaved)/(dt*1e-9)
      qsaved=q
      end

C  noise current on the collector from electron statistics
      real function shot_noise(t,dt)	! Amperes as a function of ns
      real t,dt
      real Icc
      parameter (Icc=1e-3)
      real rangauss
      external rangauss
      real ecount
      ecount=Icc*(dt*1e-9)/1.6e-19
      shot_noise=rangauss()*sqrt(ecount)*1.6e-19/(dt*1e-9)
      end

C  noise on the op-amp inputs from the manufacturer's specs (model 1)
      real function opamp_noise1(t,dt)	! Amperes as a function of ns
      real t,dt
      real Vin,R0
      parameter (Vin=4.8e-9,R0=5e3)
      real rangauss
      external rangauss
      opamp_noise1=rangauss()*Vin/sqrt(4*dt*1e-9)/R0
      end

C  noise on the op-amp inputs from the manufacturer's specs (model 2)
      real function opamp_noise2(t,dt)	! Amperes as a function of ns
      real t,dt
      real Vin,R,C,C0
      parameter (Vin=4.8e-9,R=50e3,C=2.0e-12,C0=200e-12)
      real q,qsaved
      common /qsave2/qsaved
      data qsaved/0/
      real rangauss
      external rangauss
      q=rangauss()*Vin/sqrt(4*dt*1e-9)*(C+C0)
      opamp_noise2=rangauss()*Vin/sqrt(4*dt*1e-9)/R
     +            +(q-qsaved)/(dt*1e-9)
      qsaved=q
      end

c simple two-pole model of the HPD+preamplifier response (model 1)
      subroutine filter1(id,f,nbins,xlo,xhi)
      integer id
      real f
      integer nbins
      real xlo,xhi
      real x,dx
      real sig(999999),amp
      real filtf(999999)
      real R,Rt,R0,R1,C,Gbw
      parameter (R0=3e+3,R=5e+3,C=1.0e-12,Rt=100,Gbw=2e9)
      real ceq(3),discrim
      real wpole(2),wfreq
      R1=1/(1/R+1/Rt)
      ceq(1)=Gbw/(R*C)*1e-18
      ceq(2)=(Gbw+1/(C*R1))*1e-9
      ceq(3)=1
      discrim=4*ceq(1)*ceq(3)-ceq(2)**2
      if (discrim.ge.0) then
        wpole(1)=ceq(2)/(2*ceq(3))
        wpole(2)=wpole(1)
        wfreq=sqrt(discrim)/(2*ceq(3))
      else
        wpole(1)=(ceq(2)-sqrt(-discrim))/(2*ceq(3))
        wpole(2)=(ceq(2)+sqrt(-discrim))/(2*ceq(3))
        wfreq=0
      endif
      print *, 'poles are ',wpole,wfreq
      dx=(xhi-xlo)/nbins
      do i=1,nbins
         if (wfreq.eq.0) then
           filtf(i)=-(Gbw*R0)/(R1*C)*1e-18
     +              *(exp(-wpole(1)*i*dx)-exp(-wpole(2)*i*dx))
     +              /(wpole(1)-wpole(2))
         else
           filtf(i)=+(Gbw*R0)/(C*R1)*1e-18
     +              *sin(wfreq*i*dx)*exp(-wpole(1)*i*dx)
     +              /wfreq
         endif
      enddo
      x=xlo-dx/2
      call hbook1(id,'filtered signal',nbins,xlo,xhi,0.)
      do i=2,nbins
         x=x+dx
         sig(i)=f(x,dx)
         amp=sig(1)*filtf(i-1)*dx
         do j=2,i-1
           amp=amp+(sig(j)*dx+R1*C*(sig(j)-sig(j-1)))*filtf(i-j)
         enddo
         call hfill(id,x,0.,amp)
      enddo
      end

c simple two-pole model of the HPD+preamplifier response (model 2)
      subroutine filter2(id,f,nbins,xlo,xhi)
      integer id
      real f
      integer nbins
      real xlo,xhi
      real x,dx
      real sig(999999),amp
      real filtf(999999)
      real R,C,R0,C0,Gbw
      parameter (R=50e+3,C=1.5e-12,R0=0,C0=200e-12,Gbw=1.6e9)
      real ceq(3),discrim
      real wpole(2),wfreq
      ceq(1)=Gbw*1e-18
      ceq(2)=(1+Gbw*(R*C+R0*C0))*1e-9
      ceq(3)=R*(C+C0)+R0*C0*(1+Gbw*R*C)
      discrim=4*ceq(1)*ceq(3)-ceq(2)**2
      if (discrim.ge.0) then
        wpole(1)=ceq(2)/(2*ceq(3))
        wpole(2)=wpole(1)
        wfreq=sqrt(discrim)/(2*ceq(3))
      else
        wpole(1)=(ceq(2)-sqrt(-discrim))/(2*ceq(3))
        wpole(2)=(ceq(2)+sqrt(-discrim))/(2*ceq(3))
        wfreq=0
      endif
      print *, 'poles are ',wpole,wfreq
      dx=(xhi-xlo)/nbins
      do i=1,nbins
         if (wfreq.eq.0) then
           filtf(i)=-Gbw*R/ceq(3)*1e-18
     +              *(exp(-wpole(1)*i*dx)-exp(-wpole(2)*i*dx))
     +              /(wpole(1)-wpole(2))
         else
           filtf(i)=+Gbw*R/ceq(3)*1e-18
     +              *sin(wfreq*i*dx)*exp(-wpole(1)*i*dx)
     +              /wfreq
         endif
      enddo
      x=xlo-dx/2
      call hbook1(id,'filtered signal',nbins,xlo,xhi,0.)
      do i=2,nbins
         x=x+dx
         sig(i)=f(x,dx)
         amp=sig(1)*filtf(i-1)*dx
         do j=2,i-1
           amp=amp+sig(j)*filtf(i-j)*dx
         enddo
         call hfill(id,x,0.,amp)
      enddo
      end

c  plugin replacement for filter functions, useful as a debugging check
      subroutine dintegrate(id,f,nbins,xlo,xhi)
      integer id
      real f
      integer nbins
      real xlo,xhi
      real R,C,CD,Gbw
      parameter (R=50e3,C=2.0e-12,CD=200e-12,Gbw=2e+9)
      real x,dx,Vout,Vout1,Vout2
      real Vin
      real c0,c1,c2
      c0=1/R
      c1=(C+1/(Gbw*R))*1e9
      c2=((C+CD)/Gbw)*1e18
      dx=(xhi-xlo)/nbins
      x=xlo-dx/2
      Vout=0
      Vout1=0
      Vout2=0
      call hbook1(id,'integrated signal',nbins,xlo,xhi,0.)
      do i=1,nbins
         x=x+dx
         Vout2=(f(x,dx)-c0*Vout-c1*Vout1)/c2
         Vout1=Vout1+Vout2*dx
         Vout=Vout+Vout1*dx
         Vin=-(Vout1/Gbw)*1e+9
         call hfill(id,x,0.,Vout)
      enddo
      end

c  Gaussian random number generator, based on CERNLIB rndm()
      real function rangauss()
      real r,phi,pi
      parameter (pi=3.1415926)
      external rndm
      r=sqrt(-2*log(rndm()))
      phi=2*pi*rndm()
      rangauss=r*cos(phi)
      end