#!/bin/csh
#BSUB -q 1nh
##BSUB -J jfp

if (${?WORKDIR} == 0)   then
     setenv WORKDIR $SCRATCH
endif

setenv ROOTSYS /afs/cern.ch/cms/external/lcg/external/root/5.17.04/slc3_ia32_gcc323/root
setenv LD_LIBRARY_PATH $ROOTSYS/lib
set path=($path $ROOTSYS/bin)

cp -r ~andrey/scratch1/pythia/lhapdf-5.3.0 .
cd lhapdf-5.3.0
make clean
./configure --prefix=$WORKDIR/lhapdf
make
make install
cd ../
ls -l

cat > pp14.f <<@EOF






      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)
      INTEGER PYK,PYCHGE,PYCOMP
      INTEGER Z
      REAL PT, NUMJ
       INTEGER njet2, nsets, IDproc
       REAL x1, x2, xf, Q2, cs
       REAL pxj2, pyj2, pzj2, pej2, pmass, etaj2, phij2, ptj2
       COMMON /JETM/ IDproc, njet2, x1, x2, xf, Q2,  pxj2(50),
     &pyj2(50), pzj2(50), etaj2(50), phij2(50), ptj2(50)

       COMMON /CROSS/ nsets, cs(3)

       EXTERNAL PYDATA
       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
       COMMON/PYDATR/MRPY(6),RRPY(100)
       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
      INTEGER NWPAWC
      PARAMETER (NWPAWC=300000)
      COMMON /PAWC/ PAW(NWPAWC)

*-------------------------------------------
      CALL HLIMIT(300000)
      CALL HROPEN(200,'CMS', 'newout.hbook', 'N', 1024, ISTAT)
      IDN = 101
      CALL HBNT (IDN,'',' ')
      CALL HBNAME(IDN,'JETM1',IDproc,'IDproc,njet2[0,50],
     &x1,x2,xf,Q2,pxj2(njet2),pyj2(njet2),pzj2(njet2),
     &etaj2(njet2),phij2(njet2),ptj2(njet2)')
      IDN = 102
      CALL HBNT (IDN,'',' ')
      CALL HBNAME(IDN,'JETM2',IDproc,'IDproc,njet2[0,50],
     &x1,x2,xf,Q2,pxj2(njet2),pyj2(njet2),pzj2(njet2),
     &etaj2(njet2),phij2(njet2),ptj2(njet2)')
      IDN = 103
      CALL HBNT (IDN,'',' ')
      CALL HBNAME(IDN,'JETM3',IDproc,'IDproc,njet2[0,50],
     &x1,x2,xf,Q2,pxj2(njet2),pyj2(njet2),pzj2(njet2),
     &etaj2(njet2),phij2(njet2),ptj2(njet2)')
      IDN = 111
      CALL HBNT (IDN,'',' ')
      CALL HBNAME(IDN,'CROSS1',nsets,'nsets[0,3],cs(nsets)')
      nsets = 3
*-----------------------------------------------------------------

      ECM = 14000
      CKIN(3) = 5.
C-----------------------------------------------------------------
      CALL SetPDFPath('./lhapdf/share/lhapdf/PDFsets')
      MSTP(52)=2

      MRPY(1) = $2          ! Generator value (D = 19780503)

      NEVT = $3                   ! Number of Events

      MSTU(51) = 114   ! number of eta bins
      MSTU(52) = 72    ! number of +phi bins
      PARU(51) = 5.    ! max pseudorapidity
      PARU(52) = 1.    ! min Et of initiator 
      PARU(53) = 10.   ! Et threshold for jet
      PARU(54) = 0.5   ! cone size
      MSTU(54) = 3     ! (px,py,pz,E,m) ; mass calculated

      MSTP(51) = 20550           ! MRST2006NNLO
      CALL PYINIT('CMS','p','p',ECM)
      DO 110 IEVT=1,NEVT
      CALL PYEVNT
      CALL PYEDIT(2)

      IDproc = MSTI(1)
      x1 = PARI(33)
      x2 = PARI(34)
      xf = PARI(35)
      Q2 = PARI(22)

      CALL PYCELL(NJET)
      njet2 = NJET
      IF(NJET.NE.0) THEN
       DO I_COUNT = N+1, N+NJET
       I_COUNT_1 = I_COUNT - N
       pxj2(I_COUNT_1) = P(I_COUNT,1)
       pyj2(I_COUNT_1) = P(I_COUNT,2)
       pzj2(I_COUNT_1) = P(I_COUNT,3)
       etaj2(I_COUNT_1) = PYP(I_COUNT,19)  
       phij2(I_COUNT_1) = PYP(I_COUNT,15)   
       ptj2(I_COUNT_1) = SQRT(P(I_COUNT,1)**2+P(I_COUNT,2)**2)     
       ENDDO
      ENDIF
      
      CALL HFNT(101)
110    CONTINUE
      cs(1) = PARI(1)
      CALL PYSTAT(1)

      MSTP(51) = 10050          ! CTEQ6M
      CALL PYINIT('CMS','p','p',ECM)       
      DO 210 IEVT=1,NEVT
      CALL PYEVNT
      CALL PYEDIT(2) 
      IDproc = MSTI(1)  
      x1 = PARI(33)
      x2 = PARI(34)
      xf = PARI(35)
      Q2 = PARI(22)
      CALL PYCELL(NJET)
      njet2 = NJET
      IF(NJET.NE.0) THEN
       DO I_COUNT = N+1, N+NJET
       I_COUNT_1 = I_COUNT - N
       pxj2(I_COUNT_1) = P(I_COUNT,1)
       pyj2(I_COUNT_1) = P(I_COUNT,2)
       pzj2(I_COUNT_1) = P(I_COUNT,3)
       etaj2(I_COUNT_1) = PYP(I_COUNT,19)
       phij2(I_COUNT_1) = PYP(I_COUNT,15)
       ptj2(I_COUNT_1) = SQRT(P(I_COUNT,1)**2+P(I_COUNT,2)**2)
       ENDDO
      ENDIF
      CALL HFNT(102)
210    CONTINUE
       cs(2) = PARI(1)
      CALL PYSTAT(1)

      MSTP(51) = 60300          ! ZEUS2005_ZJ
      CALL PYINIT('CMS','p','p',ECM)
      DO 310 IEVT=1,NEVT
      CALL PYEVNT
      CALL PYEDIT(2)
      IDproc = MSTI(1)
      x1 = PARI(33)
      x2 = PARI(34)
      xf = PARI(35)
      Q2 = PARI(22)
      CALL PYCELL(NJET)
      njet2 = NJET
      IF(NJET.NE.0) THEN
       DO I_COUNT = N+1, N+NJET
       I_COUNT_1 = I_COUNT - N
       pxj2(I_COUNT_1) = P(I_COUNT,1)
       pyj2(I_COUNT_1) = P(I_COUNT,2)
       pzj2(I_COUNT_1) = P(I_COUNT,3)
       etaj2(I_COUNT_1) = PYP(I_COUNT,19)
       phij2(I_COUNT_1) = PYP(I_COUNT,15)
       ptj2(I_COUNT_1) = SQRT(P(I_COUNT,1)**2+P(I_COUNT,2)**2)
       ENDDO
      ENDIF
      CALL HFNT(103)
310    CONTINUE
      cs(3) = PARI(1)
      CALL PYSTAT(1)

      CALL HFNT(111)

      CALL HROUT(0, ICYCLE, '')
      CALL HREND('CMS')
      END  

@EOF


echo ' ' 
echo Compliling pythia

g77 pp14.f -o pp14 ./lhapdf/lib/libLHAPDF.a `cernlib pythia6205 packlib mathlib`

./pp14

h2root newout.hbook nx.root
rfcp nx.root /castor/cern.ch/user/a/andrey/lhapdf/pp14TeV_lha_${1}.root

rm -r ./lha*