      PROGRAM ESERC11
C----------------------------------------------------------------------
C-
C-   Purpose and Methods :
C-
C-   Inputs  :
C-   Outputs :
C-   Controls:
C-
C-   Created   7-MAR-2003 A. Porta
C-   Modified  6-APR-2003 L. Ramello
C-
C----------------------------------------------------------------------
      IMPLICIT NONE
C----------------------------------------------------------------------
      REAL tau,t,sigma,u(200),n(200),FUNC,x,fprob,erfc,m,fp,L(200),
     &   taum(200),d,L1,fp1,S,d1
      REAL Ls,sigmatau,tau1,tau2,taumedio
      INTEGER i,j,k,h,r,pippo,pluto,g,z
      INTEGER hmemor,icycle
      COMMON/pawc/hmemor(100000)
      CALL hlimit(100000)
CC      OPEN (unit=8,name='eserc11.hbk',form='unformatted',
CC     &  status='new',recl=1024,access='direct')
      CALL hbook1(1,'Lifetime - experiment 1',80,-4.,12.,0.)
      CALL hbook1(10,'Lifetime - experiment 2',80,-4.,12.,0.)
      CALL hbook1(2,'gaussian numbers - exp. 1',60,-3.,3.,0.)
      CALL hbook1(3,'gaussian numbers - exp. 2',60,-3.,3.,0.)
      CALL hbook1(20,'uniform numbers - exp. 1',50,0.,1.,0.)
      CALL hbook1(30,'uniform numbers - exp. 2',50,0.,1.,0.)
      CALL hbook1(4,'exponential * gaussian function',
     &            1600,-3.995,12.005,0.)
      CALL hbook1(5,'-ln(L) vs. tau - experiment 1',200,0.005,2.005,0.)
      CALL hbook1(6,'Pull',40,-4.,4.,0.)
C   inizializzazione dati
      tau=1.
      sigma=0.5
      S=0.
      d=0.
C   loop su 100 esperimenti 
      DO g=1,100
        CALL ranlux(u,200)  ! estraggo 200 numeri uniformi tra 0 e 1
        DO h=1,200
          CALL norran(n(h)) ! estraggo 200 numeri da gaussiana norm.
        ENDDO
        pippo=0
        pluto=0
        L(9)=1000000.
C   loop su tau da 0.1 a 2 per calcolo -log(L)
        DO k=10,200
         taum(k)=k/100.
         L1=0.
         L(k)=0.
C   loop su 200 misure di un esperimento
          DO i=1,200
            t=-tau*log(1-u(i)) ! estrazione da distr. esponenziale
            IF(g.EQ.1.AND.k.EQ.10)CALL hf1(20,u(i),1.) ! controllo numeri uniformi
            IF(g.EQ.2.AND.k.EQ.10)CALL hf1(30,u(i),1.) ! controllo numeri uniformi
            IF(g.EQ.1.AND.k.EQ.10)CALL hf1(2,n(i),1.) ! controllo numeri gaussiani
            IF(g.EQ.2.AND.k.EQ.10)CALL hf1(3,n(i),1.) ! controllo numeri gaussiani
            t=t+sigma*n(i)
            IF(g.eq.1.AND.taum(k).EQ.0.1)CALL hf1(1,t,1.) ! controllo tempi
            IF(g.eq.2.AND.taum(k).EQ.0.1)CALL hf1(10,t,1.) ! controllo tempi
C   calcolo probabilita'
            fp=(1/(2*taum(k)))*exp(((sigma**2)/(2*(taum(k))**2))-
     &        (t/taum(k)))*erfc((sigma/(1.414213562*taum(k)))-
     &        (t/(1.414213562*sigma)))
            IF(fp.LE.0.)THEN
              fp=0.00000000000001
            ENDIF
C   calcolo verosimiglianza 
            L(k)=-log(fp)+L(k)
          ENDDO ! fine loop su 200 misure
          IF(g.EQ.1)THEN
C            WRITE(3,*)L(k),'k',k
            CALL hf1(5,taum(k),L(k)) ! controllo curva -ln(L) vs. tau (esp. 1)
          ENDIF
C   trovo taumedio
          IF(L(k).GT.L(k-1).AND.pluto.NE.1.AND.taum(k).gt.0.5)THEN
            pluto=1
            taumedio=taum(k-1)
            Ls=L(k-1)+0.5    ! calcolo -Lmax+0.5 per avere sigma
          ENDIF
        ENDDO ! fine loop su tau 
C   loop per trovare le sigma
        DO r=10,200
          IF(Ls.GE.L(r+1).AND.Ls.LE.L(r).AND.PIPPO.NE.1)THEN
            tau1=taum(r)
            pippo=1
          ENDIF
          IF(Ls.LE.L(r+1).AND.Ls.GE.L(r).AND.PIPPO.EQ.1)THEN
            tau2=taum(r+1)
          ENDIF
        ENDDO ! fine loop per trovare le sigma
        sigmatau=(tau2-tau1)/2
C        WRITE(3,*)g,'tau',taumedio,'sigma',sigmatau,tau2,tau1
       d=(taumedio-tau)/sigmatau

C media pesata sugli errori
C        S=S+(1./sigmatau**2)
C        d=d+(taumedio/sigmatau**2)
C        d1=d/S
C        WRITE(3,*)'d',d1

        CALL hf1(6,d,1.) ! distribuzione del pull su tutti gli esperimenti
      ENDDO ! fine loop sugli esperimenti

C   verifica: distribuzione teorica di probabilita' per i tempi
      DO j=1,1600
        m=j
        x=(m-400.)/100.
        FUNC=(1/(sqrt(2*3.14)*sigma))*exp(-(x)**2/(2*(sigma)**2))*
     &   (1/tau)*EXP(-x/tau)
        fprob=(1/(2*tau))*exp(((sigma**2)/(2*tau**2))-(x/tau))
     &    *erfc((sigma/(1.414213562*tau))-(x/(1.414213562*sigma)))
        CALL hf1(4,x,fprob*40.)
      ENDDO

CC      call hrfile(8,'histo','n')
CC      call hrout(0,icycle,' ')
CC      call hrend('histo')
      call hrput(0,'eserc11.hbk','N')

      END
