pepsiMaineRHIC_noradcorr.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pepsiMaineRHIC_noradcorr.f

C...As an example, consider a main program of the form
*======================================================================
        program pepsimainerhic

        implicit none
        include "leptou.inc"
        include "lujets.inc"
        include "ludat1.inc"
        include "mcRadCor.inc"
        include "phiout.inc"
        include "mc_set.inc"

        common/lypara/ipy(80),pypar(80),pyvar(80)
        INTEGER ipy
        REAL pypar, pyvar
        SAVE /lypara/

        integer nevgen, i, j, lastgenevent, genevent
        integer iev, nev, nprt
        integer idum1, idum2, initseed
        integer ltype, tracknr, t1, t2
        integer*4 today(3), now(3)
        real rval, alpha
        real mcfixedweight, mcextraweight
        double precision dilut, f1, f2, r, da1, da2
        double precision gendilut, genf1, genf2, genr, gena1, gena2
        DOUBLE PRECISION radgamenucl
        DOUBLE PRECISION pbeame, pbeta, pgamma, ebeame, epznucl
        DOUBLE PRECISION depol, gamma2, d, eta, eps, chi
        DOUBLE PRECISION gendepol, gend, geneta, geneps, genchi
        DOUBLE PRECISION dbrack, dxsec

        logical usepypar21
        real pypar21val

        CHARACTER param*20

c---------------------------------------------------------------------
c     ASCII output file
c ---------------------------------------------------------------------
        integer asciilun
        parameter(asciilun=29)
        CHARACTER*256 outputfilename
        CHARACTER*256 outname

        alpha = 1./137.
C Initialize event generation counter
        qedrad = 0
        nevgen = 0
        genevent = 0
        lastgenevent = 0.
        dilut= 1.
        r = 0.
        depol =0. 
        da1 = 0.
        da2 = 0.
        f1 = 0.
        f2 = 0.
        gamma2= 1.
        d = 0.
        eta = 0.
        eps = 0.
        chi = 0.
        gendilut= 1.
        genr = 0.
        gendepol =0.
        gena1 = 0.
        gena2 = 0.
        genf1 = 0.
        genf2 = 0.
        gend = 0.
        geneta = 0.
        geneps = 0.
        genchi = 0.
        radgamenucl = 0.
        mcfixedweight=1.
        mcextraweight=1.
        weight=1.
C Init all the parameters and read them in 
        ltype = 11
        pbeam = 100.
        ebeam = 4.

        usepypar21=.false.

C...Read parameters for LINIT call
C...Read output file name
       READ(*,*) outname
C...Read beam target particle type
       READ(*,*) ltype
C...Read beam and target particle energy
       READ(*,*) pbeam, ebeam
C...Read number of events to generate, and to print.
       READ(*,*) nev,nprt
C...Read min/max x of radgen lookup table
       READ(*,*) mcset_xmin, mcset_xmax
C...Read min/max y of generation range      
       READ(*,*) mcset_ymin, mcset_ymax
C...Read min/max Q2 of generation range      
       READ(*,*) mcset_q2min, mcset_q2max
C...Read information for cross section used in radgen
       READ(*,*) genset_fstruct, genset_r
C...Read switch for rad-corrections
       READ(*,*) qedrad
C...Read min/max W2 of generation range      
       READ(*,*) mcset_w2min,mcset_w2max
C...Read min/max Nu of generation range      
       READ(*,*) mcset_numin,mcset_numax
C...Read Beam and Target polarisation direction
       READ(*,*) mcset_pbvalue, mcset_ptvalue
C...Read Beam and Target polarisation 
       READ(*,*) mcset_pbeam, mcset_ptarget
C...Read all the pepsi/lepto parameters
  100  READ(*,*, end=200) param, j, rval
       if (param.eq.'LST') then
          lst(j)=rval
          write(*,*) 'LST(',j,')=',rval
       endif
       if (param.eq.'PARL') then
          parl(j)=rval
          write(*,*) 'PARL(',j,')=',rval
       endif
       if (param.eq.'PARJ') then
          parj(j)=rval
          write(*,*) 'PARJ(',j,')=',rval
       endif
       if (param.eq.'MSTJ') then
          mstj(j)=rval
          write(*,*) 'MSTJ(',j,')=',rval
       endif
       if (param.eq.'MSTU') then 
          mstu(j)=rval
          write(*,*) 'MSTU(',j,')=',rval
       endif
       if (param.eq.'PARU') then
          paru(j)=rval
          write(*,*) 'PARU(',j,')=',rval
       endif
       if (param.eq.'PYPAR') then
          pypar(j)=rval
          write(*,*) 'PYPAR(',j,')=',rval
         if (j.eq.21) then
            usepypar21=.true.
            pypar21val=rval
         endif
       endif
       goto 100
  200  write(*,*) '*********************************************'
       write(*,*) 'NOW all parameters are read by PEPSI'
       write(*,*) '*********************************************'

C ... PEPSI polarization
        if (mcset_pbvalue*mcset_ptvalue.lt.0.) then
          lst(40) = 1.
        elseif (mcset_pbvalue*mcset_ptvalue.gt.0.) then
          lst(40) = -1.
        else
          lst(40) = 0.
        endif

C ... polarization for helium 3 !
C     use lst(39)=1 for ideal helium 3 (no S prime and D state)
C     use lst(39)=3 or (02 should work too ?) for a real helium 3.
C       lst(39) = 0
        if (lst(40).ne.0.and.parl(1).eq.3.and.parl(2).eq.2) then
          lst(39) = 3
          lst(36) = 1
        endif

C ... xmin,xmax
        cut(1)=mcset_xmin
        cut(2)=mcset_xmax
C ... ymin,ymax
        cut(3)=mcset_ymin
        cut(4)=mcset_ymax
C ... Q2min, Q2max
        cut(5)=mcset_q2min
        cut(6)=mcset_q2max
C ... W2min, W2max
        cut(7)=mcset_w2min
C        cut(8)=mcSet_W2Max
C ... numin, numax
C        cut(9)=mcSet_NuMin
C        cut(10)=mcSet_NuMax

        mcset_tara=parl(1)
        mcset_tarz=parl(2)

C     Getting the date and time of the event generation
      call idate(today)   ! today(1)=day, (2)=month, (3)=year
      call itime(now)     ! now(1)=hour, (2)=minute, (3)=second
C      
C     Take date as the SEED for the random number generation
      initseed = today(1) + 10*today(2) + today(3) + now(1) + 5*now(3)
      write(6,*) 'SEED = ', initseed 
      call rndmq(idum1,idum2,initseed,' ')

      if (abs(ltype).eq.11) then
         masse=0.000510998910
      else
        write(*,*) "Only positron and electron beam supported"
      endif
      
      if ((parl(1).eq.1).and.(parl(2).eq.1)) then
         massp=0.938272013
      elseif ((parl(1).eq.1).and.(parl(2).eq.0)) then 
         massp=0.93956536
      elseif ((parl(1).eq.2).and.(parl(2).eq.1)) then
         massp=1.875612793
      elseif ((parl(1).eq.3).and.(parl(2).eq.2)) then
         massp=3.0160293  
      else
        write(*,*) "Currently only p,n,d,he-3 as target supported"
      endif

      pbeame=sqrt(pbeam**2+massp**2)
      pbeta=pbeam/pbeame
      pgamma=pbeame/massp
      ebeame=sqrt(ebeam**2+masse**2)
      ebeamenucl=pgamma*ebeame-pgamma*pbeta*(-ebeam)
      epznucl=-pgamma*pbeta*(ebeame)+pgamma*(-ebeam)
      write(*,*) ebeamenucl, ebeame, epznucl, -ebeam
      mcset_enebeam=sngl(ebeamenucl)

C Initialize PEPSI
       CALL timex(t1)
       call linit(0,ltype,(-ebeam),pbeam,lst(23))
       CALL timex(t2)

C PYPAR21 must be set after LINIT which sets it to PARP(1)
      if (usepypar21) then
         pypar(21)=pypar21val
      endif

c ---------------------------------------------------------------------
c     Open ascii output file
c ---------------------------------------------------------------------
       outputfilename=outname
       open(asciilun, file=outputfilename)
       write(*,*) 'the outputfile will be named: ', outname

C The total integrated xsec, in pbarns
       mcextraweight = parl(23)

! ============== Generate until happy (past cuts) =================

C   This is what we write in the ascii-file

        write(29,*)' PEPSI EVENT FILE '
        write(29,*)'============================================'
        write(29,30)
30      format('I, ievent, genevent, process, subprocess, nucleon,
     &  struckparton, partontrck, trueY, trueQ2, trueX, trueW2,
     &  trueNu, FixedWeight, weight, dxsec, Extraweight, 
     &  dilut, F1, F2, A1, A2, R, Depol, d, eta, eps, chi, 
     &  gendilut, genF1, genF2, genA1, genA2, genR, genDepol, gend,
     &  geneta, geneps, genchi, Sigcor, radgamEnucl, nrTracks')

        write(29,*)'============================================'
        write(29,*)' I  K(I,1)  K(I,2)  K(I,3)  K(I,4)  K(I,5)
     &  P(I,1)  P(I,2)  P(I,3)  P(I,4)  P(I,5)  V(I,1)  V(I,2)  V(I,3)'
        write(29,*)'============================================'

       DO 300 iev=1,nev

C Count generated events
10      nevgen = nevgen+1

        if (qedrad.eq.0) then
          If (iev.eq.1) then
              write(*,*) 'Generating full event with LEPTO/PEPSI.'
          endif
          DO i=1,mstu(4)
           DO j=1,5
              k(i,j)=0.
              p(i,j)=0.
              v(i,j)=0.
           enddo
          enddo
          lst(7) = 1
          call lepto
          if (lst(21).ne.0) then
              write(*,*) 'PEPSI failed, trying again, 
     &                    no rad-corrections'
              goto 10
          endif
          da2=0.
          da1=0.
          f1=0.
          f2=0.
          dilut=1.
          r=0.
          depol=0.
          gamma2=0.
          d=0.
          eta=0.
          eps=0.
          chi=0.
          genevent=nevgen-lastgenevent

        endif

        tracknr=n

C Get some info for the event
       call mkasym(dble(q2), dble(x), mcset_tara, mcset_tarz, da1, da2)
       call mkf2(dble(q2), dble(x), mcset_tara, mcset_tarz, f2, f1)
       call mkr(dble(q2), dble(x), r)
       call fdilut(dble(q2), dble(x), mcset_tara, dilut)

C... Calculate some kinematic variables
       gamma2=q2/u**2
       depol=(y*((1.+gamma2*y/2.)*(2.-y)-2*y**2*masse**2/q2))/
     &       (y**2*(1.-2.*masse**2/q2)*(1.+gamma2)+2.*(1.+r)*
     &       (1.-y-gamma2*(y**2)/4.))
       d=depol*((sqrt(1.-y-gamma2*(y**2)/4.)*(1+gamma2*y/2.)) /
     &   (1.-y/2.)*(1+gamma2*y/2.)-y**2*masse**2/q2)
       eta=(sqrt(gamma2)*(1.-y-gamma2*y**2/4.-y**2*masse**2/q2)) /
     &     ((1.+gamma2**2*y/2.)*(1-y/2.)-y**2*masse**2/q2)
       chi=(sqrt(gamma2)*(1.-y/2.-y**2.*masse**2/q2))/(1.+gamma2*y/2.)
       eps=1./(1.+1./2.*(1.-2.*masse**2/q2)*((y**2+gamma2*y**2)/
     &        (1.-y-gamma2*(y**2)/4.)))

        dbrack = (4.*paru(1)*alpha**2)/q2**2
c cross section in hbar^2c^2/GeV^2/dQ^2/dx
        dxsec = (1.-y-(massp*x*y/2./ebeame))*f2/x+y**2*f1
c cross sectin in mbarn/dQ^2/dx
        dxsec = dbrack * dxsec * 0.389379292

         write(29,32) 0, iev, genevent, lst(23), lst(24), lst(22),
     &   lst(25), lst(26), y, q2, x, w2, u, mcfixedweight, weight, 
     &   dxsec, mcextraweight, real(dilut), real(F1), 
     &   real(F2), real(dA1), real(dA2), real(R), real(Depol), real(d), 
     &   real(eta), real(eps), real(chi), real(gendilut), real(genF1),
     &   real(genF2), real(genA1), real(genA2), real(genR), 
     &   real(genDepol), real(gend), real(geneta), real(geneps), 
     &   real(genchi),  mcradcor_sigcor, radgamenucl, tracknr
 32      format((i4,1x,$),(i10,1x,$),6(i6,1x,$),
     &          33(e18.10,1x,$),i12,/)
         write(29,*)'============================================'

         DO i=1,tracknr
         write(29,34) i,k(i,1),k(i,2),k(i,3),k(i,4),k(i,5),
     &        p(i,1),p(i,2),p(i,3),p(i,4),p(i,5),
     &        v(i,1),v(i,2),v(i,3)
         ENDDO
 34      format(2(i6,1x,$),i10,1x,$,3(i6,1x,$),8(f15.6,1x,$),/)
         write(29,*)'=============== Event finished ==============='

         lastgenevent=nevgen

  300    CONTINUE

C...Print the Pythia cross section which is needed to get an absolut 
C   normalisation the number is in pbarns
       write(*,*)'==================================================='
       write(*,*)'Pepsi total cross section in pb from numerical', 
     &            ' integration', parl(23)
       write(*,*)'Pepsi total cross section in pb from MC simulation',
     &            parl(24)
       write(*,*)'Total Number of generated events', nev
       write(*,*)'Total Number of trials', nevgen
       write(*,*)'==================================================='
       close(29)
       END