c Multi-Poisson fit to SiPM charge spectra
c Richard Jones
c Feb. 21, 2007

      real function fitpeaks(x)
      real x
      common /pawpar/par(9)
      real par
      real q
      integer p,maxp
      parameter (maxp=12)
      integer s,maxs
      parameter (maxs=12)
      real root2pi
      real peak,sigma
      real expfact,expfacts
      real gamma
      external gamma
      root2pi=sqrt(2*3.1415926)
      expfact=exp(-par(2))
      q=(x-par(7))/par(6)
      fitpeaks=0
      do p=0,maxp
        peak=0
        do s=0,maxs
          expfacts=exp(-p*par(3))
          sigma=sqrt(par(4)**2+(p+s)*par(5)**2)
          peak=peak+expfact*(par(2)**p)/gamma(p+1.)
     +             *expfacts*((p*par(3))**s)/gamma(s+1.)
     +             *exp(-0.5*((q-p-s)/sigma)**2)/(root2pi*sigma)
        enddo
        fitpeaks=fitpeaks+peak
      enddo
      fitpeaks=fitpeaks*par(1)
      end

      real function fitpeaksf(x)
      real x
      common /pawpar/parf(9)
      real parf
      real fitpeaks
      external fitpeaks
      vector par(9)
      do i=1,9
        parf(i)=par(i)
        print *, i,par(i)
      enddo
      fitpeaksf=fitpeaks(x)
      end

      subroutine ppeaks
      print *, 'ppeaks loaded'
      end