*
* KUMAC pour montrer le fit min Q2 d'un histogramme
*
his/cre/1d 100 'Histogramme a fitter' 20 0. 1.0
his/cre/1d 200 'Coefficient a0' 25  0.0 0.5
his/cre/1d 201 'Coefficient a1' 25 -1.0 0.5
his/cre/1d 202 'Coefficient a2' 30  0.0 1.5
his/cre/1d 300 'Q2   minimum  ' 25  0.0 50.
his/cre/1d 301 'chi2 minimum  ' 25  0.0 50.
opt stat
set stat 1111
opt fit
set fit  1111111
*                             Generer un echantillon de histogrammes
vec/cre  par(3) r
sigma s=0.015
sigma alph0=0.20 ; sigma alph1=-.40 ; sigma alph2=+.60
sigma xmax=alph0+alph1+alph2
vec/cre xx(20) r .025 .075 .125 .175 .225 .275 .325 .375 .425 .475_
                 .525 .575 .625 .675 .725 .775 .825 .875 .925 .975
vec/cre yy(20) r
*                             Nombre de histogrammes (attention: c'est lent!)
ne = 1000
*                             Chaque histogramme contient un echantillon
do i=1,[ne]
   his/op/reset 100
   sigma x=rndm(array(1000))
   sigma mu=(alph0+alph1*x+alph2*x**2)/(alph0+alph1+alph2)
   sigma r=rndm(array(1000))
   do j=1,1000
      if (r([j])<mu([j])) then
         call hfill(100,x([j]),0.,1.)
      endif
   enddo
*                             Pour chaque histo estimer les coefficients
   sigma print [i]
   his/fit 100 p2 'qn' 3 par
   his/get/cont 100 yy
   sigma rr=yy-(par(1)+par(2)*xx+par(3)*xx**2)
   sigma q2=vsum(rr**2/yy)
   sigma par=par*(alph0+alph1/2.+alph2/3.)/($hinfo(100,'entries')*.05)
   call hfill(200,par(1),0.,1.)
   call hfill(201,par(2),0.,1.)
   call hfill(202,par(3),0.,1.)
   call hfill(300,q2,0.,1.)
enddo
*                             Calculer distribution de chi2
sigma c=array(25,1.#49.)
sigma fc=[ne]*(c/2.)**7.5*exp(-c/2.)/gamma(8.5)
his/put/con 301 fc
*                             Dernier echantillon, fit et vraie fonction
zone 1 1
max 100 80.
his/fit 100 p2 'q0' 3 par
his/plo 100 'e'
fun/plo $hinfo(100,'entries')*.05*(.2-.4*x+.6*x**2)/.2 0. 1. 's'
wait
*                             Distributions des coefficients
zone 2 2
his/fit 200 g
his/fit 201 g
his/fit 202 g
*                             Distributions de chi2
his/plot 300
his/plot 301 'cs'