* macro fitgp id sumevents chopt low high poly
* author Mark Lakata <lakata@att.net>
* Version 1.1
* 
* id         I 'Histogram ID'
* sumevents  I '1 = show integral of Gaussian, 0 don't show'
* chopt      C 'Chopt passed to h/fit'
* low        R 'low data range of fit region' = 0
* high       R 'high data range of fit region' = 0
* poly       I 'number of parameters in background polynomial'  =2
* 
* This macro simplies fitting a simple gaussian + polynomial background
* to a 1D histogram. The seed values for the fit are calculated automatically
* and the integral of the gaussian is calculated.  The macro requires the
* 'addstat' macro.  The default is a linear polynomial (a + bx).  To fit 
* with no polynomial, use poly=0 as an argument.
* 
* example:
* 
* exec fitgp 10
* 
* exec fitgp 10 low=120.0 high=150.0   ! fit only from 120MeV to 150MeV
* 
* 
* 
* 
 macro fitgp id=-1 sumevents=1 chopt='Q' low=0 high=0 poly=2

   if ([id] = -1) then
     message 'you need to give the histo id'
     message macro fitgp id=-1 sumevents=1 chopt=' ' low=0 high=0 poly=2
     message 'note: poly is the number of additional parameters, not DOF'
     exitm
   endif
   
   if [poly]>6 then
     mess Maximum of 6 polynomial parameters
     exitm
   endif
   if ($hinfo([id],'sum')=0) then
     message empty histogram
     h/pl [id]
     exitm
   endif
 
   rtp = 2.506628275
   binw  = ($hinfo([id],'xmax')-$hinfo([id],'xmin'))/$hinfo([id],'xbins')
   mean  = $hinfo([id],'mean')
   width = $hinfo([id],'rms')
   amp   = $hinfo([id],'sum')*[binw]/[rtp]/[width]
   
   fdim = [poly]+3
   
   if $vexist(p) then
     v/del p
   endif
   if $vexist(ep) then
     v/del ep
   endif
   if $vexist(pmin) then
     v/del pmin
   endif
   if $vexist(pmax) then
     v/del pmax
   endif
   if $vexist(step) then
     v/del step
   endif
   
   case ([poly]) in
     (0) v/create p(3) r [amp] [mean] [width]
     fname = 'g'
     bname = 0
     (1) v/create p(4) r [amp] [mean] [width] 0.0
     fname = 'g+p0'
     bname = '$rsigma(p(4))'
     (2) v/create p(5) r [amp] [mean] [width] 0.0 0.0
     fname = 'g+p1'
     bname = '$rsigma(p(4))+($rsigma(p(5))*x)'
     (3) v/create p(6) r [amp] [mean] [width] 0.0 0.0 0.0
     fname = 'g+p2'
     bname = '$sigma(p(4))+($sigma(p(5))*x)+($sigma(p(6))*x**2)'
     (4) v/create p(7) r [amp] [mean] [width] 0.0 0.0 0.0 0.0
     fname = 'g+p3'
     bname = '$sigma(p(4))+($sigma(p(5))*x)+($sigma(p(6))*x**2)+($sigma(p(7))*x**3)'
     (5) v/create p(8) r [amp] [mean] [width] 0.0 0.0 0.0 0.0 0.0
     fname = 'g+p4'
     bname = '$sigma(p(4))+($sigma(p(5))*x)+($sigma(p(6))*x**2)+($sigma(p(7))*x**3)+($sigma(p(8))*x**4)'
     (6) v/create p(9) r [amp] [mean] [width] 0.0 0.0 0.0 0.0 0.0 0.0
     fname = 'g+p5'
     bname = '$sigma(p(4))+($sigma(p(5))*x)+($sigma(p(6))*x**2)+($sigma(p(7))*x**3)+($sigma(p(8))*x**4)+($sigma(p(9))*x**5)'
   endcase
   
   v/create ep([fdim]) r
   v/create step([fdim]) r
   v/create pmin([fdim]) r
   v/create pmax([fdim]) r
   
   if [low]=[high] then
     hid=[id]
     low  = $hinfo([id],xmin)
     high = $hinfo([id],xmax)
   else
     hid=[id]($rsigma([low]):$rsigma([high]))
   endif
   
   h/fit id=[hid] func=[fname] np=[fdim] par=p errpar=ep step=step pmin=pmin pmax=pmax chopt=[chopt]
   
   if (bname.ne.'') then
     f/plot [bname] [low] [high] chopt=s
   endif
* mess bname = [bname]
     
   
   sum   = p(1)*[rtp]*p(3)/[binw]
   esum  = [sum]*ep(1)/p(1)
   
   message Events in gaussian = [sum] +- [esum]
   if [sumevents]=1 then
     do if=1,[fdim]
       fname  =$word('norm mean sigma a bx cx2 dx3 ex4 fx5',[if],1,' ')
       exec addstat 'P'//[if]//' '//[fname] '' relline=[if]
     enddo
* exec addstat 'P2 mean' '' relline=2
* exec addstat 'P3 sigma' '' relline=3
     exec addstat 'events' [sum] nfitpar=[fdim]
* exec addstat 'events' [sum]//'"a#'//[esum] npar=[fdim]
* igset txal 00
* igset chhe 0.40
* igset tang 0.0
* itx $hinfo([id],'xmin') $hinfo([id],'max') events=[sum] +- [esum]
   endif
   
 return