program cocktail

   use genesis_m

   integer, parameter  :: nevents = 2000
   integer, parameter  :: nch = 150
   real, parameter     :: theta_lower = 0.1, theta_upper = 0.3
   logical, parameter  :: write_hb = .true., write_data = .false. 
   character(len=64)   :: hbfile = 'data'
   character(len=64)   :: outfile = 'none'
   
   integer             :: istat, icycle, idec, iele, hmemor
   integer             :: iquest, ievent, ipart, i
   integer             :: pairs = 0
   integer             :: clock_at_start, clock_at_end
   integer             :: clock_now, clock_rate
   integer, parameter  :: nwpawc = 4000000
   integer, parameter  :: ntags = 12
   integer, parameter  :: lunhb = 10, lunout = 20
   integer, parameter  :: id = 1
   real                :: tuple(ntags), the, thp
   character(len=8)    :: tags(ntags)
   type(tree), pointer :: ptree, pchild, pele, ppos
   real, external      :: atg
   common /pawc/ hmemor(nwpawc)
   common /quest/ iquest(100)

   !
   ! Init Hbook, open big Ntuple file
   !
   call hlimit(nwpawc)
   iquest(10) = 65000
   if (write_hb) then
      call hropen(lunhb,'HISTOS',hbfile,'NQ',1024,istat)
      if (istat.ne.0) then
         print *,'cocktail: error opening hbook file. stop.'
         stop
      endif
   endif

   !
   ! Book Ntuple
   !
   tags(1) = 'mass'   ! invariant mass of pair
   tags(2) = 'pe'     ! p of electron
   tags(3) = 'the'    ! theta of electron
   tags(4) = 'phe'    ! phi of electron
   tags(5) = 'pp'     ! p of positron
   tags(6) = 'thp'    ! theta of positron
   tags(7) = 'php'    ! phi of positron
   tags(8) = 'p'      ! p of parent particle
   tags(9) = 'th'     ! theta of parent particle
   tags(10) = 'ph'    ! phi of parent particle
   tags(11) = 'dec'   ! decay id
   tags(12) = 'wt'    ! weight factor for each decay
   call hbookn(id,'genesis',ntags,'HISTOS',1024,tags)
   print *,'booking histos...'
   call hbook1(10,'nch vs y',200,-1.,7.,0.)
   call hbook1(11,'pi0 vs y',200,-1.,7.,0.)
   call hbook1(12,'pi+ vs y',200,-1.,7.,0.)
   call hbook1(13,'pi- vs y',200,-1.,7.,0.)
   

   !
   ! Open data file
   !
   if (write_data) then
      open(lunout,file=outfile,form='unformatted',status='unknown',iostat=istat)
      if (istat.ne.0) then
	  print *,'cocktail: error opening output file. stop.'
	  stop
      endif
   endif

   !
   ! Init Genesis
   !
   call GeInit('d.genesisrc')

   !
   !  Print out
   !
   print *,'Run started: ',eurodate()
   print *,'----------------------------------------------------'
   print *,'Events................',nevents
   print *,'Multiplicity..........',nch
   print *,'Hbook file ...........',trim(hbfile)
   print *,'Output file ..........',trim(outfile)
   print *,'Acceptance (theta)....',theta_lower,theta_upper
   print *,'----------------------------------------------------'

   !
   ! Timing
   !
   call system_clock(count=clock_at_start)

   !
   ! Event loop
   !
   event_loop: do ievent = 1, nevents

      call GeTriggerEvent(nch,0.0)                     ! trigger one event

      if (mod(ievent,10) == 0) then                    ! print-out
         call system_clock(count=clock_now,count_rate=clock_rate)
         print '("event: ",i5,5x,"still ",i4," min")', &
                 ievent, &
	         nint((float(nevents)/ievent-1.)* &
		 float(clock_now-clock_at_start)/clock_rate/60.)
      endif

      particle_loop: do ipart = 1, lasttree          ! loop all particles

         ptree => treeList(ipart)               ! pointer to tree
	 
	 if (.not. ptree%decayed) then
	    if (ptree%parttype%charge /= 0 ) then
	       call hf1(10,sngl(etaof(ptree%p4%p)),ptree%mass_weight)
	       if (ptree%parttype%name == 'pion0') &
	          call hf1(11,sngl(etaof(ptree%p4%p)),ptree%mass_weight)
	       if (ptree%parttype%name == 'pion+') &
	          call hf1(12,sngl(etaof(ptree%p4%p)),ptree%mass_weight)
	       if (ptree%parttype%name == 'pion-') &
	          call hf1(13,sngl(etaof(ptree%p4%p)),ptree%mass_weight)
	    endif
	 endif

         !
         !  Loop children and check if e+e- pair is present
         !
         iele = 0
         child_loop: do i = 1, ptree%decaytype%nBody
            pchild => treeList(ptree%child(i))
            if (pchild%parttype%name == "electron") then
                pele => treeList(ptree%child(i))
                iele = iele + 1
            else if (pchild%parttype%name == "positron") then
                ppos => treeList(ptree%child(i))
                iele = iele + 1
            endif
         enddo child_loop
         if (iele < 2) cycle particle_loop

         !
         ! Assign decay id
         !
         idec = 0
         if (ptree%decaytype%name == 'eta->electron+positron+photon')   idec = 1
         if (ptree%decaytype%name == 'rho->electron+positron')          idec = 2
         if (ptree%decaytype%name == 'omega->pion0+electron+positron')  idec = 3
         if (ptree%decaytype%name == 'omega->electron+positron')        idec = 4
         if (ptree%decaytype%name == 'phi->electron+positron')          idec = 5
         if (ptree%decaytype%name == 'phi->eta+electron+positron')      idec = 6
         if (ptree%decaytype%name == 'pion0->photon+electron+positron') idec = 7
         if (ptree%decaytype%name == "eta'->electron+positron+photon")  idec = 8
         if (idec.eq.0) cycle particle_loop

         !
         ! Store decay and fill Ntuple if pair is within acceptance
         !
         the = thetaof(pele%p4%p)
         thp = thetaof(ppos%p4%p)
         if (the <= theta_lower .or. the >= theta_upper ) cycle particle_loop
         if (thp <= theta_lower .or. thp >= theta_upper ) cycle particle_loop
         tuple(1) = abs(pele%p4 + ppos%p4)                        ! invariant mass
         tuple(2) = abs(pele%p4%p)                                ! electron momentum
         tuple(3) = the                                           ! electron theta
         tuple(4) = atg(sngl(pele%p4%p%py),sngl(pele%p4%p%px))    ! electron phi
         tuple(5) = abs(ppos%p4%p)                                ! positron momentum
         tuple(6) = thp                                           ! positron theta
         tuple(7) = atg(sngl(ppos%p4%p%py),sngl(ppos%p4%p%px))    ! positron phi
         tuple(8) = abs(ptree%p4%p)                               ! parent p
         tuple(9) = thetaof(ptree%p4%p)                           ! parent theta
         tuple(10) = atg(sngl(ptree%p4%p%py),sngl(ptree%p4%p%px)) ! parent phi
         tuple(11) = idec                                         ! decay
         tuple(12) = pele%mass_weight                             ! weight

         if (write_hb) call hfn (id,tuple)
         if (write_data) write(lunout) tuple
	 pairs = pairs + 1

      enddo particle_loop

   enddo event_loop

   !
   !  Close files
   !
   if (write_hb) then
      call hrout (0,icycle,' ')
      call hrend ('HISTOS')
      close (lunhb)
   endif
   if (write_data) close (lunout)

   !
   ! Timing
   !
   call system_clock(count=clock_at_end,count_rate=clock_rate)

   !
   ! Say goodbye
   !
   print *,'Run ended: ',eurodate()
   print *,'----------------------------------------------------'
   print *,'Events processed......',ievent-1
   print *,'Pairs generated.......',pairs
   print *,'Pairs per event.......',pairs/float(ievent-1)
   print *,'CPU total (min).......', &
           float(clock_at_end-clock_at_start)/clock_rate/60.
   print *,'CPU event (s).........', &
           float(clock_at_end-clock_at_start)/clock_rate/ievent
   print *,'CPU pair (ms).........', &
           1000.*float(clock_at_end-clock_at_start)/clock_rate/pairs 
   print *,'----------------------------------------------------'
contains

   function eurodate()
      character(len=19) :: eurodate
      character(len=8)  :: cdate
      character(len=10) :: ctime

      call date_and_time(date=cdate,time=ctime)
      eurodate = cdate(7:8)//'/'//cdate(5:6)//'/'//cdate(1:4)//' '// &
                 ctime(1:2)//':'//ctime(3:4)//':'//ctime(5:6)
   end function eurodate

end program