      PROGRAM READCWN
C--------------------------------------------------------------
C Reads columnwise ntuples.
C--------------------------------------------------------------
      IMPLICIT NONE
*------------------
* local variables
*------------------
      INTEGER LUN_OUT/80/
      INTEGER ICYCLE
      INTEGER IPASS,ISTA2
      INTEGER NPAR
*      
      LOGICAL L_EOD, L_EOF
      INTEGER ILINE,IERR,I,J,ISTA
      INTEGER      NHMEMOR
      REAL          HMEMOR
      PARAMETER   (NHMEMOR=1000000)
      COMMON /PAWC/ HMEMOR(NHMEMOR)
      
      INTEGER NFILES,IFILE
      INTEGER NOENT
      CHARACTER*80 NAMENTUPLE 
      INTEGER NTUPLID,LUNIT/1/
      PARAMETER (NTUPLID=5555)

*------------------
* functions used 
*------------------
      LOGICAL HEXIST
      REAL HI
*------------------
* variabili nelle ntuple
*------------------
      include "readcwn.inc"
*------------------
* altre variabili
*------------------

c---


C
C Initialise HBOOK
C
      CALL HLIMIT(NHMEMOR)
      CALL HCDIR('//PAWC',' ')
      CALL HBSET('BSIZE',4096,ISTA)

C
C Initialise event counting
C      
      CALL INIT

C
C Open the ntuple files
C      
      L_EOF=.FALSE.
      CALL HROPEN(LUNIT,'ALPHA',nomefile//'.hbook',' ',1024,ISTA)
      IF(ISTA.EQ.0) THEN
         CALL HRIN(NTUPLID,999999,0)
         CALL HBNAME(NTUPLID,' ',0,'$CLEAR')
         CALL HBNAME(NTUPLID,'TRANCH',CIRC,'$SET')
         CALL HBNAME(NTUPLID,'PTMIS',PTMT,'$SET')
         CALL HBNAME(NTUPLID,'LEPTONS',NLE_CWN,'$SET')
         CALL HBNAME(NTUPLID,'TOP',NTOP_CWN,'$SET')
         CALL HBNAME(NTUPLID,'GAMMAS',NGAM_CWN,'$SET')
         CALL HBNAME(NTUPLID,'JETS',NJG_CWN,'$SET')
      ELSE
         WRITE(*,*)'** Cannot open FILE ',nomefile//'.hbook'
         L_EOF = .TRUE.
      ENDIF
c inizializzare variabili a zero
      call init_inputs
      open(unit=19,file=nomefile//'.ascii0',status='UNKNOWN')
      open(unit=20,file=nomefile//'.ascii',status='UNKNOWN')
      open(unit=21,file=nomefile//'.ascii2',status='UNKNOWN')
c      open(unit=30,file=nomefile//'.btag',status='UNKNOWN')
*     
*     read the input file till its end
*     
      CALL HNOENT(NTUPLID,NOENT)
      WRITE(6,*) 'NUMBER OF ENTRIES',NOENT
      ILINE = 0
      DO WHILE (.NOT.L_EOF)
         ILINE = ILINE + 1
         CALL HCDIR('//ALPHA',' ')
         IF(HEXIST(NTUPLID)) THEN
            CALL HGNT(NTUPLID,ILINE,IERR) !<<<<***** read NTUPLE 
         ENDIF
         IF(IERR.NE.0)THEN  
            L_EOF = .TRUE.
            GOTO 999
         ENDIF
C     
C     Loop over particles and call analysis
C     
         if (mod(iline,5000).eq.0)write(6,*) 'Entry:',ILINE
         call count_lept
         call count_jet
         call corr
         call select_event
 999  ENDDO
      write(6,*)'end of file',ILINE-1
c     CALL HCDIR('//ALPHA',' ') 
      CALL HREND('ALPHA') 
      CLOSE(1)

c      call lfit(vect1,vect3,counter,1,a,b,var)
      call correlation(vect1,vect4,counter,r1) ! Correlazione Pmiss Plept
      call correlation(vect5,vect6,counter,r2) ! Correlazione eta_l eta_b


c Stampa dei risultati (su schermo e file):
      call STAMPA
      END

c----------------------------------------------------------

      subroutine INIT
      implicit none

      integer i,j,k
      include "readcwn.inc"
c---
*...Constants
      PI     = ACOS(-1.)
      TWOPI  = 2.* PI
      EPS    = 1.E-6
      EPSSQ  = 1.E-36
      BIGVAL = 1.E6

      write(6,*) "--------------------------------"
      write(6,*) "Kind of events analyzed: ", nomefile
      write(6,*) "--------------------------------"

      nltot=0
      nliso=0
      nltrue=0
      nltrueiso=0
      njtot=0
      njbtagged=0
      nbtot=0
      nctot=0
      nudstot=0
      nbbtagged=0
      ncbtagged=0
      nudsbtagged=0

      nevtot=0
      nevpresel=0
      nevsel=0
      do i=1,11
         ntrigger(i)=0
      enddo

      counter=0
      call vzero(vect1,nnmax)
      call vzero(vect2,nnmax)
      call vzero(vect3,nnmax)
      call vzero(vect4,nnmax)
      call vzero(vect5,nnmax)
      call vzero(vect6,nnmax)
      a=0.
      b=0.
      var=0.

      return
      end

c----------------------------------------------------------

      subroutine INIT_INPUTS
      implicit none

      integer i,j,k
      include "readcwn.inc"
c---
      circ=0.
      tren=0.
      trma=0.

      call vzero(ptmt,2)
      call vzero(ptm,2)
      call vzero(pznu,2)
      
      nle_cwn=0
      call vzero(pl_cwn,4*nlmax)
      do i=1,nlmax
         kl_cwn(i)=0
         mrl_cwn(i)=0
         lisol_cwn(i)=0
      enddo
      
      ntop_cwn=0
      call vzero(ptop_cwn,4*maxtop)

      ngam_cwn=0
      call vzero(pg_cwn,4*ngmax)
      do i=1,ngmax
         lisgam_cwn(i)=0
      enddo

      njg_cwn=0
      call vzero(pjg_cwn,4*nclus)
      call vzero(btag_cwn,nclus)
      do i=1,nclus
         mrkjet_cwn(i)=0
      enddo


      return
      end

c------------------------------------------

      subroutine STAMPA
      implicit none

      integer i,j,k
      integer iunit
      include "readcwn.inc"
      real fract, staterror
      character*1 mode
c---
      if (selectionmode.eq.1) mode='1'
      if (selectionmode.eq.2) mode='2'
      if (selectionmode.eq.2) mode='3'

      do i=1,2
         if (i.eq.1) then
            iunit=6
         elseif (i.eq.2) then
            iunit=10
            open(unit=10,file=nomefile//'.out'//mode,
     &           status='UNKNOWN')
            open(unit=30,file=nomefile//'.sn',
     &           status='UNKNOWN')
         endif
         write(iunit,*) "------------------------------"
         write(iunit,*) "Output for ",nomefile," events"
         write(iunit,*) "------------------------------"
         write(iunit,*) "* Leptons:"
         write(iunit,*) "total number: ", nltot, 
     &        " (W->l: ", nltrue, ")"
         write(iunit,*) "isolated leptons: ",nliso,
     &        " (W->l: ", nltrueiso, ")"
         write(iunit,*) "[", fract(nliso,nltot)," %  (",
     &        fract(nltrueiso,nltrue)," % ) ]"
         write(iunit,*) "* Jets:"
         write(iunit,*) "total number: ", njtot, 
     &        " (b: ", nbtot, ", c: ", nctot,", uds:", nudstot,")"
         write(iunit,*) "b-tagged (cut >", btagcut ,
     &        "): ",njbtagged,
     &        " (b: ", nbbtagged, ", c: ", ncbtagged,", uds:", 
     &        nudsbtagged,")"
         write(iunit,*) "[", fract(njbtagged,njtot)," %  (",
     &        fract(nbbtagged,nbtot)," %,", fract(ncbtagged,nctot),
     &        " %,", fract(nudsbtagged,nudstot),"% ) ]"
         write(iunit,*) "* Ptlept(isolated) vs Ptmiss:"
c         write(iunit,*) "slope = ",a
         write(iunit,*) "correlation coefficient = ",r1
         write(iunit,*) "* eta lept(isolated) vs best b-jet:"
         write(iunit,*) "correlation coefficient = ",r2
c         write(iunit,*) "b = ",b
c         write(iunit,*) "residual sum of squares (normalized): ",var
         write(iunit,*) "* Selection: "
         write(iunit,*) "method: 1l,2j(1b)"//
     &        ", with Ptjet > ", ptjetcut," ,btag > ", btagcut
         write(iunit,*) "total events: ", nevtot
         write(iunit,*) "triggered events: ", 
     &        (fract(ntrigger(j),nevtot),j=1,10), " %"
         write(iunit,*) "total triggered events: ", ntrigger(11)
         write(iunit,*) "[", fract(ntrigger(11),nevtot)," +-",
     &        staterror(ntrigger(11),nevtot)," % ]"
         write(iunit,*) "preselected events: ", nevpresel
         write(iunit,*) "[", fract(nevpresel,nevtot)," +-",
     &        staterror(nevpresel,nevtot)," % ]"
         if (i.eq.2) write(30,*) fract(nevpresel,nevtot)/100.
         if (i.eq.2) write(30,*) staterror(nevpresel,nevtot)/100.
         write(iunit,*) "selected events: ", nevsel
         write(iunit,*) "[", fract(nevsel,nevtot)," +-",
     &        staterror(nevsel,nevtot)," % ]"
* Write efficiency on file:
         if (i.eq.2) write(30,*) fract(nevsel,nevtot)/100.
         if (i.eq.2) write(30,*) staterror(nevsel,nevtot)/100.
c         write(iunit,*) " "
c         write(iunit,*) " "
c         write(iunit,*) " "
c         write(iunit,*) " "

      enddo

      return
      end

c----------------------------------------------------------

      subroutine COUNT_LEPT
      implicit none

      integer i,j,k
      include "readcwn.inc"
c---

      nltot=nltot+nle_cwn
      do i=1,nle_cwn
         if (lisol_cwn(i).eq.0) nliso=nliso+1
c         if (abs(kl_cwn(i)).eq.1.or.abs(kl_cwn(i)).eq.2) then ! true e,mu
         if (abs(mrl_cwn(i)).eq.24) then ! W->l
            nltrue=nltrue+1
            if (lisol_cwn(i).eq.0) nltrueiso=nltrueiso+1
         endif
      enddo

      return
      end

c----------------------------------------------------------

      subroutine COUNT_JET
      implicit none

      integer i,j,k
      include "readcwn.inc"
c---

c      njtot=njtot+njg_cwn
      do i=1,njg_cwn
         if (pt(1).gt.ptjetcut.and.eta(1).lt.2.4) then
            njtot=njtot+1
            if (btag_cwn(i).ge.btagcut) njbtagged=njbtagged+1
            if (mrkjet_cwn(i).eq.5) then
               nbtot=nbtot+1
               if (btag_cwn(i).ge.btagcut) nbbtagged=nbbtagged+1
            elseif (mrkjet_cwn(i).eq.4) then
               nctot=nctot+1
               if (btag_cwn(i).ge.btagcut) ncbtagged=ncbtagged+1
            else
               nudstot=nudstot+1
               if (btag_cwn(i).ge.btagcut) nudsbtagged=nudsbtagged+1
            endif
         endif
      enddo


      return
      end

c----------------------------------------------------------

      subroutine CORR
* correlation between Pmiss and Plept, eta_l and eta_b
      implicit none

      integer i,j,k
      include "readcwn.inc"

      real temp
c      real maxbtag
c      integer bestb
      real etafun
c---
      ptmiss=sqrt(ptm(1)*ptm(1)+ptm(2)*ptm(2))

      maxbtag=-10.
      bestb=-1
      do i=1,njg_cwn
         if (btag_cwn(i).ge.btagcut) then
            temp=temp+1
            if (btag_cwn(i).gt.maxbtag) then
               maxbtag=btag_cwn(i)
               bestb=i ! il miglior b-jet dell'evento
            endif
         endif
      enddo

      do i=1,nle_cwn
         if (lisol_cwn(i).eq.0) then
            counter=counter+1
            if (counter.eq.nnmax) then
               write(6,*) "WARNING!!!! counter=nnmax"
               goto 777
            endif
            vect1(counter)=ptmiss
c            vect2(counter)=pl_cwn(4,i)
            temp=pl_cwn(4,i)
            vect2(counter)=temp ! se non faccio cosi' fa segmentation fault! mah!!!
c            vect3(counter)=sqrt(pl_cwn(1,i)**2+pl_cwn(2,i)**2+
c     &           pl_cwn(3,i)**2)
            temp=sqrt(pl_cwn(1,i)**2+pl_cwn(2,i)**2+
     &           pl_cwn(3,i)**2)
            vect3(counter)=temp
            ptl=sqrt(pl_cwn(1,i)**2+pl_cwn(2,i)**2)
            vect4(counter)=ptl
            vect5(counter)=etafun(pl_cwn(1,i),pl_cwn(2,i),
     &           pl_cwn(3,i))
            vect6(counter)=etafun(pjg_cwn(1,bestb),pjg_cwn(2,bestb),
     &           pjg_cwn(3,bestb))
         endif
      enddo


 777  continue
      return
      end
c----------------------------------------------------------

      subroutine SELECT_EVENT
      implicit none

      integer i,j,k
      integer iunit
      integer temp
      include "readcwn.inc"
      real etafun,phifun,rdist
      real pawmlp
      logical highlumi
      logical trigger_ele
      logical trigger_diele
      logical trigger_mu
      logical trigger_dimu
      logical trigger_1jet
      logical trigger_3jet
      logical trigger_4jet
      logical trigger_jet_etmiss
      logical trigger_ele_jet
      logical trigger_mu_jet
      logical trigger_combo
      integer altro
c---  

      if ((nomefile.eq."ww_incl".or.nomefile.eq."zw_incl")
     &     .and.ntop_cwn.ne.0) then
* Questi sample inclusivi comprendono, nelle card, anche W*->tb
         write(6,*) 'ntop=',ntop_cwn
         goto 888
      endif

      nevtot=nevtot+1

* Trigger:
      highlumi=lumi
      if (trigger_ele(highlumi)) ntrigger(1)=ntrigger(1)+1
      if (trigger_diele(highlumi)) ntrigger(2)=ntrigger(2)+1
      if (trigger_mu(highlumi)) ntrigger(3)=ntrigger(3)+1
      if (trigger_dimu(highlumi)) ntrigger(4)=ntrigger(4)+1
      if (trigger_1jet(highlumi)) ntrigger(5)=ntrigger(5)+1
      if (trigger_3jet(highlumi)) ntrigger(6)=ntrigger(6)+1
      if (trigger_4jet(highlumi)) ntrigger(7)=ntrigger(7)+1
      if (trigger_jet_etmiss(highlumi)) ntrigger(8)=ntrigger(8)+1
      if (trigger_ele_jet(highlumi)) ntrigger(9)=ntrigger(9)+1
      if (trigger_mu_jet(highlumi)) ntrigger(10)=ntrigger(10)+1
      if (trigger_combo(highlumi)) then
         ntrigger(11)=ntrigger(11)+1
      else
         if (dotrigger) goto 888
      endif

* Cuts set:
      if (selectionmode.eq.1) then
         include "cuts_singletop.inc"
      elseif (selectionmode.eq.2) then
         include "cuts_schannel.inc"
      elseif (selectionmode.eq.3) then
         include "cuts_tchannel.inc"
      endif

* MC truth:
      if (preselswitch.eq.1) goto 1
      temp=0
      isol=-1
      do i=1,nle_cwn
         if (abs(kl_cwn(i)).eq.1.or.abs(kl_cwn(i)).eq.2
     &        .and.abs(mrl_cwn(i)).eq.24) then 
            temp=temp+1
            isol=i
         endif
      enddo
      if (temp.ne.1) goto 888   ! 1 W->l
      plept=sqrt(pl_cwn(1,isol)**2+pl_cwn(2,isol)**2+
     &     pl_cwn(3,isol)**2)
c      if (abs(kl_cwn(isol)).eq.1.and.plept.lt.ptelecut) goto 888
c      if (abs(kl_cwn(isol)).eq.2.and.plept.lt.ptmucut) goto 888
      temp=0
      ptjettot=0.
      jet1=-1
      jet2=-1
      jethighestpt=-1
      highestpt=0.
      do i=1,njg_cwn
         ptjet=sqrt(pjg_cwn(1,i)**2+pjg_cwn(2,i)**2)
         ptjettot=ptjettot+ptjet
         if (mrkjet_cwn(i).eq.5) then
            temp=temp+1
            if (jet1.eq.-1) then
               jet1=i
            else
               jet2=i
            endif
         else ! se non e' b, cerca il jet di Pt maggiore
            if (ptjet.gt.highestpt) then
               highestpt=ptjet
               jethighestpt=i
            endif
         endif
      enddo
      if (nomefile.eq."ewtop_t") then
         if (temp.ne.1) goto 888 ! 1 b-jet
         if (jethighestpt.eq.-1) goto 888 ! another jet
      else
         if (temp.ne.2) goto 888 ! 2 b-jets
      endif
      if (temp.eq.1) then
         bestb=jet1
         worstb=jethighestpt
      elseif (temp.eq.2) then
         if (btag_cwn(jet1).ge.btag_cwn(jet2)) then
            bestb=jet1
            worstb=jet2
         else
            bestb=jet2
            worstb=jet1
         endif
      endif
      pt(1)=sqrt(pjg_cwn(1,bestb)**2+pjg_cwn(2,bestb)**2)
      pt(2)=sqrt(pjg_cwn(1,worstb)**2+pjg_cwn(2,worstb)**2)
c$$$      if (pt(1).lt.ptjetcut_b) goto 888
 1    continue

* PRESELEZIONE DI TEST:
      if (preselswitch.eq.2) then

         temp=0
         isol=-1
         do i=1,nle_cwn
            if (lisol_cwn(i).eq.0) then 
               temp=temp+1
               isol=i
            endif
         enddo
         if (temp.ne.1) goto 888 ! 1 l isolato
         temp=0
         ptjettot=0.
         jet1=-1
         jet2=-1
         do i=1,njg_cwn
            ptjet=sqrt(pjg_cwn(1,i)**2+pjg_cwn(2,i)**2)
            ptjettot=ptjettot+ptjet
            if (ptjet.gt.ptjetcut) then
               temp=temp+1
               if (jet1.eq.-1) then
                  jet1=i
               else
                  jet2=i
               endif
            elseif (ptjet.le.ptjetcut.and.ptjet.gt.ptjetveto) then
               goto 888
            endif
         enddo
         if (temp.ne.2) goto 888 ! 2 high-pt jets         

         goto 2
      endif


* PRESELEZIONE: 1 l (isolato), 2 soli jets, almeno 1 b-jet
      if (preselswitch.eq.0) goto 2
      temp=0
      isol=-1
      do i=1,nle_cwn
         if (lisol_cwn(i).eq.0) then 
            temp=temp+1
            isol=i
         endif
      enddo
      if (temp.ne.1) goto 888   ! 1 l isolato
      plept=sqrt(pl_cwn(1,isol)**2+pl_cwn(2,isol)**2+
     &     pl_cwn(3,isol)**2)
      if (abs(kl_cwn(isol)).eq.1.and.plept.lt.ptelecut) goto 888
      if (abs(kl_cwn(isol)).eq.2.and.plept.lt.ptmucut) goto 888
      temp=0
      ptjettot=0.
      jet1=-1
      jet2=-1
      do i=1,njg_cwn
         ptjet=sqrt(pjg_cwn(1,i)**2+pjg_cwn(2,i)**2)
         ptjettot=ptjettot+ptjet
         if (ptjet.gt.ptjetcut) then
            temp=temp+1
            if (jet1.eq.-1) then
               jet1=i
            else
               jet2=i
            endif
         elseif (ptjet.le.ptjetcut.and.ptjet.gt.ptjetveto) then
            goto 888
         endif
      enddo
      if (temp.ne.2) goto 888   ! 2 high-pt jets         
      if (btag_cwn(jet1).ge.btag_cwn(jet2)) then
         bestb=jet1
         worstb=jet2
      else
         bestb=jet2
         worstb=jet1
      endif
      btag_before(1)=btag_cwn(bestb)
      btag_before(2)=btag_cwn(worstb)
      if (btag_cwn(bestb).lt.btagcut) goto 888 ! almeno 1 b-jet
c     write(30,*) btag_before(1),btag_before(2)
c     if (temp.eq.0) goto 888 ! almeno 1 b-jet
      pt(1)=sqrt(pjg_cwn(1,bestb)**2+pjg_cwn(2,bestb)**2)
      pt(2)=sqrt(pjg_cwn(1,worstb)**2+pjg_cwn(2,worstb)**2)
      if (pt(1).lt.ptjetcut_b) goto 888
 2    continue

* se e' arrivato fin qua vuol dire che l'evento e' stato preselezionato
      nevpresel=nevpresel+1

      if (preselswitch.eq.2) goto 888

* write on ascii:
      if (isol.eq.-1.or.bestb.eq.-1) goto 888
      pj(1)=sqrt(pjg_cwn(1,bestb)**2+pjg_cwn(2,bestb)**2+
     &     pjg_cwn(3,bestb)**2)
      pj(2)=sqrt(pjg_cwn(1,worstb)**2+pjg_cwn(2,worstb)**2+
     &     pjg_cwn(3,worstb)**2)
      q=-1.*float(kl_cwn(isol))/abs(float(kl_cwn(isol)))
      etal=etafun(pl_cwn(1,isol),pl_cwn(2,isol),pl_cwn(3,isol))
      phil=phifun(pl_cwn(1,isol),pl_cwn(2,isol))
      eta(1)=etafun(pjg_cwn(1,bestb),pjg_cwn(2,bestb),pjg_cwn(3,bestb))
      phi(1)=phifun(pjg_cwn(1,bestb),pjg_cwn(2,bestb))
      eta(2)=etafun(pjg_cwn(1,worstb),pjg_cwn(2,worstb),
     &     pjg_cwn(3,worstb))
      phi(2)=phifun(pjg_cwn(1,worstb),pjg_cwn(2,worstb))
      rlb=rdist(etal,eta(1),phil,phi(1))
      rlnonb=rdist(etal,eta(2),phil,phi(2))
      rbnonb=rdist(eta(1),eta(2),phi(1),phi(2))
      ptl=sqrt(pl_cwn(1,i)**2+pl_cwn(2,i)**2)
      ptmiss=sqrt(ptm(1)*ptm(1)+ptm(2)*ptm(2))
      ht=ptjettot+ptl+ptmiss
      btag_after(1)=btag_cwn(bestb)
      btag_after(2)=btag_cwn(worstb)
* Mreco computation:
      if (bestworstb.eq.1) ib=bestb
      if (bestworstb.eq.2) ib=worstb
      if (bestworstb.eq.1) ib2=worstb !temporary
      if (bestworstb.eq.2) ib2=bestb !temporary
      do i=1,4
         pj1(i)=pjg_cwn(i,ib)
c         pb(i)=pjg_cwn(i,ib)
         pj2(i)=pjg_cwn(i,ib2)
      enddo
      cosj1j2=vdotn(pj1,pj2,3)
      mj1j2=sqrt((pjg_cwn(4,ib)+pjg_cwn(4,ib2))**2-
     &     (pjg_cwn(1,ib)+pjg_cwn(1,ib2))**2-
     &     (pjg_cwn(2,ib)+pjg_cwn(2,ib2))**2-
     &     (pjg_cwn(3,ib)+pjg_cwn(3,ib2))**2)
      do i=1,2 ! loop over neutrino solutions
         plnu=pznu(i)
         enu=sqrt(ptmiss**2+plnu**2)
         mreco(i,1)=sqrt((pl_cwn(4,isol)+pjg_cwn(4,ib)+enu)**2-
     &        (pl_cwn(1,isol)+pjg_cwn(1,ib)+ptm(1))**2-
     &        (pl_cwn(2,isol)+pjg_cwn(2,ib)+ptm(2))**2-
     &        (pl_cwn(3,isol)+pjg_cwn(3,ib)+plnu)**2)
         mreco(i,2)=sqrt((pl_cwn(4,isol)+pjg_cwn(4,ib2)+enu)**2-
     &        (pl_cwn(1,isol)+pjg_cwn(1,ib2)+ptm(1))**2-
     &        (pl_cwn(2,isol)+pjg_cwn(2,ib2)+ptm(2))**2-
     &        (pl_cwn(3,isol)+pjg_cwn(3,ib2)+plnu)**2)
         if (toptruereco.eq.1) then
            ptop(1,i,1)=pl_cwn(1,isol)+pjg_cwn(1,ib)+ptm(1)
            ptop(1,i,2)=pl_cwn(1,isol)+pjg_cwn(1,ib2)+ptm(1)
            ptop(2,i,1)=pl_cwn(2,isol)+pjg_cwn(2,ib)+ptm(2)
            ptop(2,i,2)=pl_cwn(2,isol)+pjg_cwn(2,ib2)+ptm(2)
            ptop(3,i,1)=pl_cwn(3,isol)+pjg_cwn(3,ib)+pznu(i)
            ptop(3,i,2)=pl_cwn(3,isol)+pjg_cwn(3,ib2)+pznu(i)
            ptop(4,i,1)=pl_cwn(4,isol)+pjg_cwn(4,ib)+
     &           sqrt(ptmiss**2+pznu(i)**2)
            ptop(4,i,2)=pl_cwn(4,isol)+pjg_cwn(4,ib)+
     &           sqrt(ptmiss**2+pznu(i)**2)
            do k=1,4
               deltaptop(k,i,1)=(ptop(k,i,1)-ptop_cwn(k,1))/
     &              abs(ptop(k,i,1))
               deltaptop(k,i,2)=(ptop(k,i,2)-ptop_cwn(k,1))/
     &              abs(ptop(k,i,2))
            enddo
         elseif (toptruereco.eq.2.and.ntop_cwn.eq.1) then
            do k=1,4
               ptop(k,i,1)=ptop_cwn(k,1)
               ptop(k,i,2)=ptop_cwn(k,1)
            enddo
         elseif (toptruereco.eq.2.and.ntop_cwn.eq.2) then
            ran=ranf()
            if (ran.lt.0.5) iran=1
            if (ran.ge.0.5) iran=2
            do k=1,4
               ptop(k,i,1)=ptop_cwn(k,iran)
               ptop(k,i,2)=ptop_cwn(k,iran)
            enddo
         elseif (toptruereco.eq.2.and.ntop_cwn.eq.0) then
            do k=1,4
               ptop(k,i,1)=0
               ptop(k,i,2)=0
            enddo
         endif
         pw(1,i)=pl_cwn(1,isol)+ptm(1)
         pw(2,i)=pl_cwn(2,isol)+ptm(2)
         pw(3,i)=pl_cwn(3,isol)+pznu(i)
         pw(4,i)=pl_cwn(4,isol)+sqrt(ptmiss**2+pznu(i)**2)
         pnu(1,i)=ptm(1)
         pnu(2,i)=ptm(2)
         pnu(3,i)=pznu(i)
         pnu(4,i)=sqrt(ptmiss**2+pznu(i)**2)
      enddo

*
* choose the "best" top solution:
*
      if (mrecoswitch.eq.1) then  ! mreco(1) refers to the first neutrino solution (pz1>pz2)
         indexnu=1
      elseif (mrecoswitch.eq.2) then ! mreco(1) refers to the neutrino solution closest to the lepton 
         do i=1,2
            plnuv(i)=pznu(i)
            enuv(i)=sqrt(ptmiss**2+plnuv(i)**2)
            etanuv(i)=etafun(ptm(1),ptm(2),plnuv(i))
            phinuv(i)=phifun(ptm(1),ptm(2))
            deltar(i)=rdist(etal,etanuv(i),phil,phinuv(i))
         enddo
         if (deltar(1).gt.deltar(2)) then
            indexnu=1
         else
            indexnu=2
         endif
      elseif (mrecoswitch.eq.3) then ! mreco(1) refers to the neutrino solution that gives M(lnb) closest to M(top) 
         if (abs(mreco(1,1)-175.).lt.abs(mreco(2,1)-175.)) then
            indexnu=1
         else
            indexnu=2
         endif
      endif
      if (mrecoswitch2.eq.0) then
         indexb=1
         indexb2=2
       elseif (mrecoswitch2.eq.1) then
          do i=1,2
             if (abs(mreco(i,1)-175).lt.abs(mreco(i,2)-175)) then
                indexb=1
                indexb2=2
             else
                indexb=2
                indexb2=1
             endif
          enddo
      endif
      if (mrecoswitch3.eq.0) then
         continue
      elseif (mrecoswitch3.eq.1) then ! the best possible combination overall
         temp=99999.
         indexnu=1
         indexb=1
         indexb2=2
         do i=1,2
            do j=1,2
               if (abs(mreco(i,j)-175).lt.temp) then
                  temp=abs(mreco(i,j)-175)
                  indexnu=i
                  indexb=j
                  indexb2=altro(j)
               endif
            enddo
         enddo
      endif

* Invariant mass of the top-bottom system:
c      mtb=sqrt((pjg_cwn(4,ib2)+ptop(4))**2-
c     &     (pjg_cwn(1,ib2)+ptop(1))**2-
c     &     (pjg_cwn(2,ib2)+ptop(2))**2-
c     &     (pjg_cwn(3,ib2)+ptop(3))**2)

* Lorentz boost of the lepton in the top rest frame:
      do i=1,4
         pleptv(i)=pl_cwn(i,isol)
      enddo
      vbeam(1)=0.
      vbeam(2)=0.
      vbeam(3)=1.
      if (toptruereco.eq.1) then
         mtop=mreco(indexnu,indexb)
      elseif(toptruereco.eq.2) then
         mtop=175.
      endif
      do k=1,4
         ptopv(k)=ptop(k,indexnu,indexb)
      enddo
      CALL LORENF(mtop,ptopv,pleptv,plboost)
      CALL LORENF(mtop,ptopv,pj2,pj2boost)
      coslt=vdotn(plboost,ptopv,3)
      coslj2=vdotn(plboost,pj2boost,3)
      coslbeam=vdotn(plboost,vbeam,3)
      CALL LORENF(mtop,ptopv,pw,pwboost)
      coswt=vdotn(pwboost,ptopv,3)
      coswj2=vdotn(pwboost,pj2boost,3)
      coslw2=vdotn(plboost,pwboost,3) ! angular dist. of the lepton w.r.t. W in top rest frame
      CALL LORENF(mtop,pw,pleptv,plboost2)
      coslw=vdotn(plboost2,pw,3) ! angular dist. of the lepton w.r.t. W in W rest frame
* Lorentz boost of the b-jet (from t) in the top rest frame:
      CALL LORENF(mtop,ptopv,pj1,pbboost)
      cosbj2=vdotn(pbboost,pj2boost,3)
      costj2=vdotn(ptopv,pj2boost,3)
* Lorentz boost of the nu (from t) in the top rest frame:
      CALL LORENF(mtop,ptopv,pnu,pnuboost)
      cosnuj2=vdotn(pnuboost,pj2boost,3)
* Lorentz boost of the b-jet (from t) in the c.m. rest frame:
c      do i=1,4
c         pcm(i)=ptopv(i)+pj2(i)
c      enddo
c      CALL LORENF(mtop,pcm,pj1,pbboost)
c      CALL LORENF(mtop,pcm,pcm,pcmboost)
c      ebboost=pbboost(4)
c      ecm=pcmboost(4)

* Write in ascii file:
      if (preselswitch.eq.1) iunit=20
      if (preselswitch.eq.0) iunit=19
      if (preselswitch.eq.0.and.nevtot.gt.20000) goto 111
      write(iunit,*) q,etal,eta(1),eta(2),ptl,ptmiss,ht
     &     ,mreco(indexnu,indexb),mreco(indexnu,indexb2)
     &     ,rlb,rlnonb,rbnonb,pj(1),pj(2)
     &     ,btag_after(1),btag_after(2),coslt,coslj2
     &     ,deltaptop(1,indexnu,indexb),deltaptop(2,indexnu,indexb)
     &     ,deltaptop(3,indexnu,indexb),deltaptop(4,indexnu,indexb)
     &     ,cosbj2,cosnuj2,costj2,coslbeam
     &     ,cosj1j2,mj1j2
 111  continue

* Selection:
      do i=1,nle_cwn            ! controlla che non ci siano leptoni non isolati di pT eccessivo
         if (lisol_cwn(i).ne.0) then ! non deve contare di nuovo il leptone isolato
            temp=sqrt(pl_cwn(1,i)**2+pl_cwn(2,i)**2)
            if (temp.gt.ptlep_in_jet) goto 888
         endif
      enddo
      if (mreco(indexnu,indexb).lt.mreco_low.or.mreco(indexnu,indexb)
     &     .gt.mreco_high) goto 888
      if (ht.lt.ht_low.or.ht.gt.ht_high) goto 888
      if (ptmiss.lt.ptmiss_low.or.ptmiss.gt.ptmiss_high) goto 888
      if (pj(1).lt.pj1_low.or.pj(1).gt.pj1_high) goto 888
      if (pj(2).lt.pj2_low.or.pj(2).gt.pj2_high) goto 888
      if (btag_cwn(worstb).lt.btagcut2) goto 888 ! 2 b-jet (s-channel)
      if (btag_cwn(worstb).gt.antibtagcut) goto 888 ! 1 b-jet (t-channel)
      if (rbnonb.lt.rbnonb_low) goto 888
      if (rbnonb.gt.rbnonb_high) goto 888
      if (abs(eta(2)).gt.etanonb_high) goto 888
      if (mreco(indexnu,2)-mreco(indexnu,1)
     &     .gt.deltamreco_high) goto 888
c      if (pawmlp(pj(1),pj(2),mreco(1),ht,eta(2),rbnonb,ptmiss).lt.0.6) 
c     &     goto 888
      nevsel=nevsel+1

* Write in ascii file:
      write(21,*) q,etal,eta(1),eta(2),ptl,ptmiss,ht
     &     ,mreco(indexnu,indexb),mreco(indexnu,indexb2)
     &     ,rlb,rlnonb,rbnonb,pj(1),pj(2)
     &     ,btag_after(1),btag_after(2),coslt,coslj2
     &     ,deltaptop(1,indexnu,indexb),deltaptop(2,indexnu,indexb)
     &     ,deltaptop(3,indexnu,indexb),deltaptop(4,indexnu,indexb)
     &     ,cosbj2,cosnuj2,costj2,coslbeam
     &     ,cosj1j2,mj1j2

 888  continue
      return
      end

c----------------------------------------------------------
      integer function altro(j)
      implicit none
      integer j
c---

      if (j.eq.1) altro=2
      if (j.eq.2) altro=1

      return
      end
c----------------------------------------------------------

      real function FRACT(n,d)
      implicit none

      integer n,d
c---

      fract=100*float(n)/float(d)

      return
      end

c----------------------------------------------------------

      real function STATERROR(n,d)
      implicit none

      integer n,d
      real x
c---

      x=float(n)/float(d)
      staterror=100*sqrt(x*(1-x)/float(d))

      return
      end

c----------------------------------------------------------

      SUBROUTINE correlation(x,y,n,r)
* http://www.fundp.ac.be/~jpvigner/homepage/sources/recipes/PEARSN.F.html
      INTEGER n
      REAL prob,r,z,x(n),y(n),TINY
      PARAMETER (TINY=1.e-20)
CU    USES betai
      INTEGER j
      REAL ax,ay,df,sxx,sxy,syy,t,xt,yt,betai
      ax=0.
      ay=0.
      do 11 j=1,n
        ax=ax+x(j)
        ay=ay+y(j)
11    continue
      ax=ax/n
      ay=ay/n
      sxx=0.
      syy=0.
      sxy=0.
      do 12 j=1,n
        xt=x(j)-ax
        yt=y(j)-ay
        sxx=sxx+xt**2
        syy=syy+yt**2
        sxy=sxy+xt*yt
12    continue
      r=sxy/sqrt(sxx*syy)
      z=0.5*log(((1.+r)+TINY)/((1.-r)+TINY))
      df=n-2
      t=r*sqrt(df/(((1.-r)+TINY)*((1.+r)+TINY)))
c      prob=betai(0.5*df,0.5,df/(df+t**2))
C     prob=erfcc(abs(z*sqrt(n-1.))/1.4142136)
      return
      END

c----------------------------------------------------------

      REAL FUNCTION ETAFUN(PX,PY,PZ)

*...ETA  CALCULATIONS.

*KEEP,IMPNON.
      IMPLICIT NONE
*KEND.
      REAL PX, PY, PZ, COST

       COST = PZ / SQRT(PX**2 + PY**2 + PZ**2)
       IF(ABS(COST).LT.0.999999) THEN
          ETAFUN = 0.5 * ALOG((1.+ COST)/(1.- COST))
        ELSE
C          PRINT '('' ETAFUN : PX, PY, PX, COST = '',4E10.3)',
C      &   PZ, PY, PZ, COST
          ETAFUN=100.
       ENDIF

      END

*CMZ :          26/12/96  13.48.38  by  Salavat Abdullin
*-- Author :    Salavat Abdullin   09/12/94
*+DOC*******************************************************************
       REAL FUNCTION RDIST(ETA1,ETA2,PHI1,PHI2)
*
*  Description: distance in eta-phi space between two objects
*               defined by eta1,phi1, eta2,phi2
*
*  Calls: PHIDIS
*
*+DOC*******************************************************************

*KEEP,IMPNON.
      IMPLICIT NONE
*KEND.
      REAL ETA1, ETA2, PHI1, PHI2, PHIDIS

      RDIST = SQRT((ETA1 - ETA2)**2 + PHIDIS(PHI1,PHI2)**2)
      END
*CMZ :          01/05/2000  17.51.35  by  Salavat Abdullin
*-- Author :    N.Stepanov
************************************************************************

      REAL FUNCTION PHIDIS(PHI1,PHI2)

*...CALC. THE DISTANCE IN PHI

*KEEP,IMPNON.
      IMPLICIT NONE
*KEEP,CONSTA.
*...Useful constants.
      REAL            EPS, EPSSQ, BIGVAL, PI, TWOPI
      COMMON /CONSTA/ EPS, EPSSQ, BIGVAL, PI, TWOPI
*KEND.

      REAL PHI1, PHI2, D

      D = ABS(PHI1 - PHI2)
      PHIDIS = AMIN1(D,TWOPI-D)

      END
*CMZ :          01/05/2000  17.51.35  by  Salavat Abdullin
*-- Author :    N.Stepanov,       modified by S.Abdullin
************************************************************************

      REAL FUNCTION PHIFUN(PX,PY)

*...PHI  CALCULATIONS

*KEEP,IMPNON.
      IMPLICIT NONE
*KEEP,CONSTA.
*...Useful constants.
      REAL            EPS, EPSSQ, BIGVAL, PI, TWOPI
      COMMON /CONSTA/ EPS, EPSSQ, BIGVAL, PI, TWOPI
*KEND.

      REAL PX, PY

      IF(ABS(PX).LT.EPSSQ.OR.ABS(PY).LT.EPSSQ) THEN
         PRINT '(''PHIFUN : anomalous input : '',2(1PE13.5))',PX,PY
         IF(PY.GE.0.) THEN
            PHIFUN = 0.5 * PI
          ELSE
            PHIFUN = 1.5 * PI
         ENDIF
       ELSE
         PHIFUN = ATAN2(PY,PX)
      ENDIF
      IF(PHIFUN.LT.0.) PHIFUN = PHIFUN + TWOPI

      END


************************************************************************

      REAL FUNCTION pawmlp(pj1,pj2,mreco,ht,etanonb,rbnonb,ptmiss)
C
C    Function generated by the PAW - MLPfit interface
C
      real pj1,pj2,mreco,ht,etanonb,rbnonb,ptmiss
C
      OUT1 = TAN(PJ1-180)/5
      OUT2 = TAN(PJ2-220)/5.5
      OUT3 = TAN(MRECO-175)/5
      OUT4 = TAN(HT-350)/5.5
      OUT5 = TAN(ABS(ETANONB)-1.3)/4
      OUT6 = (RBNONB-2.6)/0.7
      OUT7 = TAN(PTMISS-50)/6
C
C     layer 2
      RIN1 = 1.170691e+00
     > +(-9.790955e-01) * OUT1
     > +(1.853928e-01) * OUT2
     > +(1.236500e+00) * OUT3
     > +(2.366655e+00) * OUT4
     > +(-1.115072e+00) * OUT5
     > +(-2.663857e+00) * OUT6
     > +(-2.407144e+00) * OUT7
      RIN2 = 6.996448e+00
     > +(6.969247e-01) * OUT1
     > +(-1.004530e-01) * OUT2
     > +(-2.708554e-01) * OUT3
     > +(-7.425382e-01) * OUT4
     > +(1.109086e+00) * OUT5
     > +(9.114064e+00) * OUT6
     > +(2.029438e+00) * OUT7
      RIN3 = -2.540457e+00
     > +(4.412740e-01) * OUT1
     > +(1.728730e+00) * OUT2
     > +(1.957380e+00) * OUT3
     > +(-4.889145e+00) * OUT4
     > +(-7.313916e-01) * OUT5
     > +(2.247548e+00) * OUT6
     > +(1.679114e+00) * OUT7
      RIN4 = 1.850075e+00
     > +(-1.541986e+00) * OUT1
     > +(-7.811122e-01) * OUT2
     > +(1.365098e+00) * OUT3
     > +(2.768455e+00) * OUT4
     > +(-1.427491e-01) * OUT5
     > +(-2.171442e-01) * OUT6
     > +(1.035885e+00) * OUT7
      RIN5 = 3.566954e+00
     > +(-3.092397e+00) * OUT1
     > +(8.288177e-01) * OUT2
     > +(3.813322e+00) * OUT3
     > +(2.972848e+00) * OUT4
     > +(2.175630e+00) * OUT5
     > +(-4.259762e-01) * OUT6
     > +(7.827602e+00) * OUT7
      RIN6 = 7.653780e-01
     > +(1.480447e-03) * OUT1
     > +(3.317100e-01) * OUT2
     > +(-3.701815e-01) * OUT3
     > +(3.732999e-01) * OUT4
     > +(5.913923e-01) * OUT5
     > +(7.445826e-01) * OUT6
     > +(9.727747e-01) * OUT7
      RIN7 = -4.753851e+00
     > +(-5.253307e-01) * OUT1
     > +(-7.268922e-02) * OUT2
     > +(9.336221e-01) * OUT3
     > +(1.887662e-01) * OUT4
     > +(2.519952e-01) * OUT5
     > +(3.908551e+00) * OUT6
     > +(8.255881e-01) * OUT7
C
      OUT1 = SIGMOID(RIN1)
      OUT2 = SIGMOID(RIN2)
      OUT3 = SIGMOID(RIN3)
      OUT4 = SIGMOID(RIN4)
      OUT5 = SIGMOID(RIN5)
      OUT6 = SIGMOID(RIN6)
      OUT7 = SIGMOID(RIN7)
C
C     layer 3
      RIN1 = -6.437652e-01
     > +(-1.779871e-02) * OUT1
     > +(-7.716693e-01) * OUT2
     > +(7.717387e-02) * OUT3
     > +(3.493250e-01) * OUT4
     > +(4.030627e-01) * OUT5
     > +(1.453769e-01) * OUT6
     > +(2.121023e-01) * OUT7
      RIN2 = -3.368318e+00
     > +(-1.123220e+00) * OUT1
     > +(-4.966046e+00) * OUT2
     > +(5.813270e-01) * OUT3
     > +(-5.070608e-01) * OUT4
     > +(6.732013e-01) * OUT5
     > +(-1.427610e+00) * OUT6
     > +(-1.310776e+00) * OUT7
      RIN3 = 4.504595e-01
     > +(7.384473e-01) * OUT1
     > +(-2.987122e-01) * OUT2
     > +(1.325858e+00) * OUT3
     > +(7.695643e-01) * OUT4
     > +(5.712001e-01) * OUT5
     > +(3.880785e-01) * OUT6
     > +(3.412979e-02) * OUT7
      RIN4 = 1.447447e+00
     > +(7.203953e-01) * OUT1
     > +(2.396933e+00) * OUT2
     > +(-8.622687e-01) * OUT3
     > +(-2.316570e-01) * OUT4
     > +(4.126376e-01) * OUT5
     > +(1.523980e+00) * OUT6
     > +(-1.886643e-01) * OUT7
      RIN5 = 3.751893e-01
     > +(4.025379e-01) * OUT1
     > +(1.078497e-01) * OUT2
     > +(9.354895e-01) * OUT3
     > +(1.012608e+00) * OUT4
     > +(3.347737e-01) * OUT5
     > +(-5.135098e-01) * OUT6
     > +(-1.324533e+00) * OUT7
      RIN6 = 9.799986e-02
     > +(1.309518e+00) * OUT1
     > +(-1.344192e+00) * OUT2
     > +(1.245869e+00) * OUT3
     > +(3.034470e-01) * OUT4
     > +(7.521171e-01) * OUT5
     > +(6.051920e-01) * OUT6
     > +(3.376608e-01) * OUT7
      RIN7 = 2.603099e+00
     > +(2.867437e-01) * OUT1
     > +(4.929474e+00) * OUT2
     > +(1.267490e+00) * OUT3
     > +(-2.068730e+00) * OUT4
     > +(-8.873354e-01) * OUT5
     > +(3.938049e-01) * OUT6
     > +(-3.663139e+00) * OUT7
C
      OUT1 = SIGMOID(RIN1)
      OUT2 = SIGMOID(RIN2)
      OUT3 = SIGMOID(RIN3)
      OUT4 = SIGMOID(RIN4)
      OUT5 = SIGMOID(RIN5)
      OUT6 = SIGMOID(RIN6)
      OUT7 = SIGMOID(RIN7)
C
C     layer 4
      pawmlp = 3.254260e-01
     > +(-3.595746e-01) * OUT1
     > +(-3.284619e+00) * OUT2
     > +(5.473597e-01) * OUT3
     > +(-6.840716e-01) * OUT4
     > +(5.710129e-01) * OUT5
     > +(-2.720685e-01) * OUT6
     > +(2.876209e-01) * OUT7
C
      END
      REAL FUNCTION SIGMOID(X)
      IF(X.GT.37.) THEN
         SIGMOID = 1.
      ELSE IF(X.LT.-37.) THEN
         SIGMOID = 0.
      ELSE
         SIGMOID = 1./(1.+EXP(-X))
      ENDIF
      END
