ReggeGribovPartonMCInterface/src/crmc-aaa.f

0001
0002 c 15.01.2009 Simplified Main program and random number generator for epos
0003
0004       subroutine crmc_set_f(iEvent,iSeed,pproj,ptarg,
0005      $     ipart,itarg,imodel,itab,ilheout,lheoutfile,param)
0006
0007 ***************************************************************
0008 *
0009 *  interface to epos subroutine
0010 *
0011 *   input: iEvent     - number of events to generate
0012 *          iseed      - random seed
0013 *          pproj      - beam momentum in GeV/c (detector frame)
0014 *          ptarg      - target momentum in GeV/c of nucleon! (detector frame)
0015 *          ipart      - primary particle type
0016 *          itarg      - target particle type
0017 *          imodel     - HE model switch
0018 *          itab       - force tables production or stop if missing
0019 *          ilheout    - output type
0020 *          lheoutfile - output file name
0021 *          param      - param file name
0022 *
0023 ***************************************************************
0024       implicit none
0025       include "epos.inc"
0026
0027 c     Input values
0028       integer iSeed,ipart,itarg,imodel,itab,ilheout,iEvent,iout
0029       double precision pproj, ptarg
0030       character*1000 param,lheoutfile,output
0031       common/lheoutput/iout,output
0032
0033       real   m1,m2
0034       double precision iecms,e1,e2
0035
0036       double precision ycm2det
0037       common/boostvars/ycm2det
0038
0039       common/producetab/ producetables              !used to link CRMC
0040       logical producetables                         !with EPOS and QII
0041
0042 c     Set parameters to default value
0043       call aaset(0)
0044
0045 c     Stop program if missing tables (after aaset)
0046       producetables=.false.
0047       if(itab.eq.1)producetables=.true.
0048
0049 c     Set common for crmc_init
0050       iout=ilheout
0051       output=lheoutfile
0052
0053 c     Calculations of energy of the center-of-mass in the detector frame
0054       call idmass(1120,m2)      !target mass = proton
0055       m1=m2                     !projectile mass
0056       if(abs(itarg).eq.120)call idmass(120,m2) !if pion as targ or proj
0057       if(abs(ipart).eq.120)call idmass(120,m1)
0058       e2=dsqrt(dble(m1)**2+ptarg**2)
0059       e1=dsqrt(dble(m2)**2+pproj**2)
0060       iecms=dsqrt((e1+e2)**2-(pproj+ptarg)**2)
0061 c     Later a rapidity boost back into the detector system will be performed
0062 c     ycm2det defines this rapidity
0063       ycm2det=0
0064 C       if (((e1+e2)-(pproj+ptarg)) .le. 0d0) then
0065 C          ycm2det=1d99
0066 C       elseif (((e1+e2)+(pproj+ptarg)) .le. 0d0) then
0067 C          ycm2det=-1d99
0068 C       else
0069 C          ycm2det=0.5d0*dlog(((e1+e2)+(pproj+ptarg))/
0070 C      +        ((e1+e2)-(pproj+ptarg)))
0071 C       endif
0072 C       if (pproj .le. 0d0) then
0073 C          ycm2det=-ycm2det
0074 C       endif
0075 c     Update some parameters value to run correctly
0076       call IniEpos(iEvent,iSeed,ipart,itarg,iecms,imodel)
0077
0078 c     The parameters can be changed optionnaly by reading a file
0079 c     (example.param) using the following subroutine call
0080       call EposInput(param)     !(it can be commented)
0081
0082 c     if you put what is in input.optns in example.param, you can even run
0083 c     exactly the same way (coded parameters are overwritten). Don't forget
0084 c     the command : "EndEposInput" at the end of example.param, otherwise it
0085 c     will not run.
0086       end
0087
0088       subroutine crmc_init_f()
0089 ***************************************************************
0090 *
0091 *  init models with values set in crmc_set_f
0092 *
0093 ***************************************************************
0094       implicit none
0095       integer iout
0096       character*1000 output
0097       common/lheoutput/iout,output
0098
0099 c     initialization for the given energy
0100       call ainit
0101
0102 c     Here the cross section sigineaa is defined
0103
0104 c     LHE type output done by EPOS
0105       if(iout.eq.1)call EposOutput(output)
0106
0107       end
0108
0109
0110       subroutine crmc_f(iout,ievent,noutpart,impactpar,outpart,outpx
0111      +                  ,outpy,outpz,oute,outm,outstat)
0112
0113 ***************************************************************
0114 *   output: iout      - output type
0115 *           ievent    - event number
0116 *           noutpart  - number of stable particles produced
0117 *           impactpar - impact parameter in fm
0118 *           outpart   - particle ID (pdg code)
0119 *           outpx     - particle momentum px
0120 *           outpy     - particle momentum py
0121 *           outpz     - particle momentum pz
0122 *           oute      - particle energy
0123 *           outm      - particle mass
0124 *           outstat   - particle status code
0125 *
0126 ***************************************************************
0127       implicit none
0128       include "epos.inc"
0129 c     Output quantities
0130       integer noutpart,iout,ievent
0131       double precision impactpar
0132       integer outpart(*)
0133       double precision outpx(*), outpy(*), outpz(*)
0134       double precision oute(*), outm(*)
0135       integer outstat(*)
0136
0137       double precision boostvec1,boostvec2,boostvec3,boostvec4,boostvec5
0138       double precision ycm2det
0139       common/boostvars/ycm2det !is this common block needed?
0140
0141       integer i!,k
0142
0143 c     Calculate an inelastic event
0144       call aepos(-1)
0145
0146 c     Fix final particles and some event parameters
0147       call afinal
0148
0149 c     Fill HEP common
0150       call hepmcstore
0151
0152 c     optional Statistic information (only with debug level ish=1)
0153       call astati
0154       if(ish.ge.1) call bfinal
0155
0156 c     Print out (same as the one defined in input.optns)
0157       if(nhep.gt.nmxhep)then
0158         print *,'Warning : produced number of particles is too high'
0159         print *,'          increase nmxhep : ',nhep,' > ',nmxhep
0160 c        stop
0161       endif
0162       noutpart=nhep
0163       impactpar=dble(bimevt)
0164 c     define vec to boost from cm. to cms frame
0165       boostvec1=0d0
0166       boostvec2=0d0
0167       boostvec3=dsinh(dble(ycm2det))
0168       boostvec4=dcosh(dble(ycm2det))
0169       boostvec5=1d0
0170 c      write(*,*)nevhep,nhep,boostvec3,boostvec4,ycm2det
0171       do i=1,nhep
0172 c     boost output to cms frame
0173          call utlob2(-1,boostvec1,boostvec2,boostvec3
0174      +        ,boostvec4,boostvec5
0175      +        ,vhep(1,i),vhep(2,i),vhep(3,i),vhep(4,i),-99)
0176          call utlob5dbl(-ycm2det
0177      +        ,phep(1,i), phep(2,i), phep(3,i), phep(4,i), phep(5,i))
0178          outpart(i)=idhep(i)
0179          outpx(i)=phep(1,i)
0180          outpy(i)=phep(2,i)
0181          outpz(i)=phep(3,i)
0182          oute(i)=phep(4,i)
0183          outm(i)=phep(5,i)
0184          outstat(i)=isthep(i)
0185 c      write(*,'(4x,i6,1x,4(e12.6,1x))')idhep(i),(vhep(k,i),k=1,4)
0186       enddo
0187
0188 c     Write lhe file
0189       if(iout.eq.1)call lhesave(ievent)
0190
0191       end
0192
0193
0194 c-----------------------------------------------------------------------
0195       subroutine IniEpos(iEvent,iSeed,ipart,itarg,iecms,iModel)
0196 c-----------------------------------------------------------------------
0197 c     Update some parameters and define path to tab files here can be set
0198 c     what particle is stable or not transfer number of event and debug
0199 c     level to main program (to avoid epos.inc in main)
0200 c-----------------------------------------------------------------------
0201       implicit none
0202       include "epos.inc"
0203
0204       double precision iecms
0205       integer iSeed,ipart,itarg,iModel,iadd,idtrafo,iEvent
0206       character*4 lhct
0207
0208       iframe=11                 !11 puts it always in nucleon nucleon reference
0209                                 !frame. This is ok because we use ecms
0210                                 !which is calculated in crmc_f.
0211
0212       model=max(1,iModel)              ! epos = 0,1 / qgsjet01 = 2 / gheisha = 3
0213                                 ! / pythia = 4 / hijing = 5 / sibyll 2.1
0214                                 ! = 6 / qgsjetII.04 = 7 / phojet = 8
0215                                 ! qgsjetII.03 = 11
0216       if(iModel.eq.0)then
0217         call LHCparameters      !LHC tune for EPOS
0218         isigma=1                !use analytic cross section for nuclear xs
0219         ionudi=1
0220         lhct=".lhc"
0221         iadd=4
0222       else
0223         isigma=2                !use pseudo-MC cross section for nuclear xs (slow)
0224         ionudi=3
0225         lhct=""
0226         iadd=0
0227       endif
0228
0229       nfnii=15                  ! epos tab file name lenght
0230       fnii="tabs/epos.initl"    ! epos tab file name
0231       nfnid=15
0232       fnid="tabs/epos.inidi"
0233       nfnie=15
0234       fnie="tabs/epos.iniev"
0235       nfnrj=15+iadd
0236       fnrj="tabs/epos.inirj"//lhct
0237       nfncs=15+iadd
0238       fncs="tabs/epos.inics"//lhct
0239
0240
0241       seedi=1.d0                !seed for random number generator: at start program
0242       seedj=iSeed               !seed for random number generator: for first event
0243       jwseed=0                  !print out seed in see file (1) or not
0244
0245
0246 c     Initialize decay of particles
0247       nrnody=0                  !number of particle types without decay
0248                                 !(if 0 (default) : all unstable
0249                                 !particles decay (at the end only
0250                                 !(anti)nucleons, (anti)electrons and
0251                                 !muons)
0252 c Particle code is given as
0253 c     id=+/-ijkl
0254 c
0255 c          mesons--
0256 c          i=0, j<=k, +/- is sign for j
0257 c          id=110 for pi0, id=220 for eta, etc.
0258 c
0259 c          baryons--
0260 c          i<=j<=k in general
0261 c          j<i<k for second state antisymmetric in (i,j), eg. l = 2130
0262 c
0263 c          other--
0264 c          id=1,...,6 for quarks
0265 c          id=9 for gluon
0266 c          id=10 for photon
0267 c          id=11,...,16 for leptons
0268 c          i=17 for deuteron
0269 c          i=18 for triton
0270 c          i=19 for alpha
0271 c          id=20 for ks, id=-20 for kl
0272 c
0273 c          i=21...26 for scalar quarks
0274 c          i=29 for gluino
0275 c
0276 c          i=30 for h-dibaryon
0277 c
0278 c          i=31...36 for scalar leptons
0279 c          i=39 for wino
0280 c          i=40 for zino
0281 c
0282 c          id=80 for w+
0283 c          id=81,...,83 for higgs mesons (h0, H0, A0, H+)
0284 c          id=84,...,87 for excited bosons (Z'0, Z''0, W'+)
0285 c          id=90 for z0
0286 c
0287 c          diquarks--
0288 c          id=+/-ij00, i<j for diquark composed of i,j.
0289 c
0290 c Examples : 2130 = lambda, 1330=xi0, 2330=xi-, 3331=omega
0291 c
0292 c Conversion from epos to pdg code can be done using
0293 c      id_pdg=idtrafo('nxs','pdg',id_epos)
0294
0295
0296 c$$$      nrnody=nrnody+1
0297 c$$$      nody(nrnody)=120     !pi+
0298 c$$$      nrnody=nrnody+1
0299 c$$$      nody(nrnody)=-120    !pi-
0300 c$$$      nrnody=nrnody+1
0301 c$$$      nody(nrnody)=130     !K+
0302 c$$$      nrnody=nrnody+1
0303 c$$$      nody(nrnody)=-130    !K-
0304 c$$$      nrnody=nrnody+1
0305 c$$$      nody(nrnody)=-20     !Kl
0306 c$$$      nrnody=nrnody+1
0307 c$$$      nody(nrnody)=-14     !mu+
0308 c$$$      nrnody=nrnody+1
0309 c$$$      nody(nrnody)=14      !mu-
0310
0311 c      nrnody=nrnody+1
0312 c      nody(nrnody)=idtrafo('pdg','nxs',3122)    !lambda using pdg code
0313
0314       iecho=0                   !"silent" reading mode
0315
0316 c     Debug
0317       ish=0                     !debug level
0318       ifch=6                    !debug output (screen)
0319 c     ifch=31    !debug output (file)
0320 c     fnch="epos.debug"
0321 c     nfnch=index(fnch,' ')-1
0322 c     open(ifcx,file=fnch(1:nfnch),status='unknown')
0323
0324       nevent=iEvent              !number of events
0325       modsho = 100000000                !printout every modsho events
0326
0327       ecms=sngl(iecms)          !center of mass energy in GeV/c2
0328 c     pnll=pproj                !beam momentum GeV/c
0329
0330       infragm=2                 !nuclear fragmentation (realistic)
0331
0332 c     Projecticle definitions
0333       if (abs(ipart) .eq. 1) then
0334 c     proton
0335          idprojin = sign(1120,ipart) !proton
0336          laproj =  sign(1,ipart)   !proj Z
0337          maproj = 1                !proj A
0338       elseif (ipart .eq. 12) then
0339 c     carbon
0340          idprojin=1120
0341          laproj = 6             !proj Z
0342          maproj = 12            !proj A
0343       elseif (ipart .eq. 208) then
0344 c     lead
0345          idprojin=1120
0346          laproj = 82            !proj Z
0347          maproj = 208           !proj A
0348       elseif (ipart .eq. 120) then
0349 c     pi+
0350          idprojin = ipart         !pi+
0351          laproj = -1            !proj Z
0352          maproj = 1             !proj A
0353       elseif (ipart .eq. -120) then
0354 c     pi-
0355          idprojin = ipart         !pi-
0356          laproj = -1            !proj Z
0357          maproj = 1             !proj A
0358 c nuclei
0359       elseif (ipart.gt.10000)then
0360          idprojin=1120
0361          maproj=mod(ipart,10000)/10           !proj A
0362          laproj=mod(ipart,10000000)/10000     !proj Z
0363 c PDG
0364       else
0365         idprojin=idtrafo('pdg','nxs',ipart)
0366         laproj = -1             !proj Z
0367         maproj = 1              !proj A
0368       endif
0369
0370       if(idprojin.eq.99)then
0371          print *,'Warning : projectile particle not known : ',ipart
0372          print *,'          id particle must be +/-120(pi+/-)'
0373          print *,'          1(proton) 12(carbon) 208(lead) or PDG'
0374          stop
0375       endif
0376
0377
0378 c     Target definitions : for nucleons, idtarg does not exist
0379 c     Mass number matarg as well as charge, latarg, must be defined
0380
0381 c     idtarg = 1120             !proton
0382       if ( itarg .eq. 1 ) then
0383 c     proton
0384          idtargin = sign(1120,itarg)
0385          latarg = sign(1,itarg) !targ Z
0386          matarg = 1             !targ A
0387       elseif ( itarg .eq. 12 ) then
0388 c     carbon
0389          idtargin=1120
0390          latarg = 6             !targ Z
0391          matarg = 12            !targ A
0392       elseif (ipart .eq. 120) then
0393 c     pi+
0394          idprojin = ipart         !pi+
0395          laproj = -1            !proj Z
0396          maproj = 1             !proj A
0397       elseif (ipart .eq. -120) then
0398 c     pi-
0399          idprojin = ipart         !pi-
0400          laproj = -1            !proj Z
0401          maproj = 1             !proj A
0402       elseif ( itarg .eq. 208 ) then
0403 c     lead
0404          idtargin=1120
0405          latarg = 82            !targ Z
0406          matarg = 208           !targ A
0407 c nuclei
0408       elseif (itarg.gt.10000)then
0409          idtargin=1120
0410          matarg=mod(itarg,10000)/10          !targ A
0411          latarg=mod(itarg,10000000)/10000    !targ Z
0412 c PDG
0413       elseif (abs(itarg).eq.2112.or.abs(itarg).eq.2212)then
0414         idtargin=idtrafo('pdg','nxs',itarg)
0415         latarg = -1             !targ Z
0416         matarg = 1              !targ A
0417       else
0418          print *,'Warning : target particle not known : ',itarg
0419          print *,'          id particle must be +/-120(pi+/-)'
0420          print *,'          1(proton) 12(carbon) 208(lead) or PDG'
0421          stop
0422
0423       endif
0424
0425       if ( model.eq.1 ) then   !model variable has no eposLHC
0426          istmax = 1             !final and mother particles (istmax=1 includes
0427                                 !mother particles)
0428       else
0429          istmax=0
0430       endif
0431
0432       end
0433
0434 c-----------------------------------------------------------------------
0435       subroutine lhesave(n)
0436 c-----------------------------------------------------------------------
0437 c     writes the results of a simulation into the file with unit ifdt
0438 c     contains a description of the stored variables.
0439 c     use Les Houches Event File as defined in hep-ph/0109068 for the
0440 c     common block and hep-ph/0609017 for the XML output.
0441 c     some code taken from example from Torbjrn Sjstrand
0442 c     in http://www.thep.lu.se/~torbjorn/lhef
0443 c-----------------------------------------------------------------------
0444       include 'epos.inc'
0445
0446       integer id
0447       real taugm
0448 C...User process event common block.
0449       INTEGER MAXNUP
0450       PARAMETER (MAXNUP=nmxhep)  !extend array for file production
0451 c      PARAMETER (MAXNUP=500)
0452       INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
0453       DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
0454       COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,IDUP(MAXNUP),
0455      &ISTUP(MAXNUP),MOTHUP(2,MAXNUP),ICOLUP(2,MAXNUP),PUP(5,MAXNUP),
0456      &VTIMUP(MAXNUP),SPINUP(MAXNUP)
0457       SAVE /HEPEUP/
0458
0459
0460 C...set event info and get number of particles.
0461       NUP=nhep             !number of particles
0462       IDPRUP=nint(typevt)  !type of event (ND,DD,SD)
0463       XWGTUP=1d0           !weight of event
0464       SCALUP=-1d0          !scale for PDF (not used)
0465       AQEDUP=-1d0          !alpha QED (not relevant)
0466       AQCDUP=-1d0          !alpha QCD (not relevant)
0467
0468 C...Copy event lines, omitting trailing blanks.
0469 C...Embed in <event> ... </event> block.
0470       write(ifdt,'(A)') '<event>'
0471       write(ifdt,*)NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP
0472       DO 220 i=1,nhep
0473
0474 c  store particle variables:
0475           IDUP(i)=idhep(i)
0476           if(isthep(i).eq.4)then
0477           ISTUP(i)=-9      !incoming particle
0478           else
0479           ISTUP(i)=isthep(i) !in LHEF:1=final, 2=decayed
0480           endif
0481           MOTHUP(1,i)=jmohep(1,i)
0482           MOTHUP(2,i)=jmohep(2,i)
0483           ICOLUP(1,i)=0        !color flow
0484           ICOLUP(2,i)=0        !color flow
0485           do J=1,5                !particle momentum (GeV/c)
0486             PUP(J,i)=phep(J,i)
0487           enddo
0488           id=idtrafo('pdg','nxs',idhep(i))
0489           call idtau(id,sngl(phep(4,i)),sngl(phep(5,i)),taugm)
0490           VTIMUP(i)=dble(taugm*(-alog(rangen())))*1d-12 !life time c*tau in mm
0491           if(VTIMUP(i).gt.dble(ainfin)
0492      &   .or.VTIMUP(i).ne.VTIMUP(i))VTIMUP(i)=ainfin
0493           SPINUP(i)=9           !polarization (not known)
0494           write(ifdt,*)IDUP(i),ISTUP(i),
0495      &      MOTHUP(1,i),MOTHUP(2,i),ICOLUP(1,i),ICOLUP(2,i),
0496      &      (PUP(J,i),J=1,5),VTIMUP(i),SPINUP(i)
0497   220 CONTINUE
0498
0499 c optional informations
0500       write(ifdt,*)'#geometry',bimevt,phievt
0501
0502       write(ifdt,'(A)') '</event>'
0503
0504       if(n.eq.nevent)then
0505 C...Successfully reached end of event loop: write closing tag
0506         write(ifdt,'(A)') '</LesHouchesEvents>'
0507         write(ifdt,'(A)') ' '
0508         close(ifdt)
0509       endif
0510
0511       return
0512       end
0513
0514
0515 c-----------------------------------------------------------------------
0516       subroutine EposOutput(iFile)
0517 c-----------------------------------------------------------------------
0518 c Use EPOS to create lhe file output
0519 c-----------------------------------------------------------------------
0520       include "epos.inc"
0521       character*1000 iFile
0522
0523       istore=4
0524       nfndt=index(iFile,'.lhe')+4   !'.lhe' extension added in bstora
0525       fndt(1:nfndt)=iFile(1:nfndt)
0526       kdtopen=0
0527
0528       call bstora
0529
0530 c      nopen=0
0531 c      ifop=52
0532 c      open(unit=ifop,file=TRIM(iParam),status='old')
0533 c      call aread
0534 c      close(ifop)
0535       end
0536
0537 c-----------------------------------------------------------------------
0538       subroutine EposInput(iParam)
0539 c-----------------------------------------------------------------------
0540 c     Read informations (new options or parameter change) in the file
0541 c     "epos.param". The unit "ifop" is used in aread. If not used, it will
0542 c     use the default value of all parameters.
0543 c-----------------------------------------------------------------------
0544       include "epos.inc"
0545       character*1000 iParam
0546
0547       nopen=0
0548       ifop=52
0549       open(unit=ifop,file=TRIM(iParam),status='old')
0550       call aread
0551       close(ifop)
0552       end
0553
0554
0555 c-----------------------------------------------------------------------
0556       subroutine utLob5dbl(yboost,x1,x2,x3,x4,x5)
0557 c-----------------------------------------------------------------------
0558 c     Same as utlob5 but in double precision
0559 c-----------------------------------------------------------------------
0560       implicit none
0561       double precision yboost,y,amt,x1,x2,x3,x4,x5
0562       amt=dsqrt(x5**2+x1**2+x2**2)
0563       y=dsign(1D0,x3)*dlog((x4+dabs(x3))/amt)
0564       y=y-yboost
0565       x4=amt*dcosh(y)
0566       x3=amt*dsinh(y)
0567       return
0568       end
0569
0570 c-----------------------------------------------------------------------
0571       subroutine crmc_xsection_f(xsigtot,xsigine,xsigela,xsigdd,xsigsd
0572      &                          ,xsloela,xsigtotaa,xsigineaa,xsigelaaa)
0573 c-----------------------------------------------------------------------
0574 c     Same as utlob5 but in double precision
0575 c-----------------------------------------------------------------------
0576       implicit none
0577       include 'epos.inc'
0578       double precision xsigtot,xsigine,xsigela,xsigdd,xsigsd
0579      &                ,xsloela,xsigtotaa,xsigineaa,xsigelaaa
0580
0581       xsigtot   = dble( sigtot   )
0582       xsigine   = dble( sigine   )
0583       xsigela   = dble( sigela   )
0584       xsigdd    = dble( sigdd    )
0585       xsigsd    = dble( sigsd    )
0586       xsloela   = dble( sloela   )
0587 c Nuclear cross section only if needed
0588       xsigtotaa = 0d0
0589       xsigineaa = 0d0
0590       xsigelaaa = 0d0
0591       if(maproj.gt.1.or.matarg.gt.1)then
0592         if(model.eq.1)then
0593           call crseaaEpos(sigtotaa,sigineaa,sigcutaa,sigelaaa)
0594         else
0595           call crseaaModel(sigtotaa,sigineaa,sigcutaa,sigelaaa)
0596         endif
0597         xsigtotaa = dble( sigtotaa )
0598         xsigineaa = dble( sigineaa )
0599         xsigelaaa = dble( sigelaaa )
0600       endif
0601
0602       return
0603       end
0604
0605 C c-----------------------------------------------------------------------
0606 C       function rangen()
0607 C c-----------------------------------------------------------------------
0608 C c     generates a random number
0609 C c-----------------------------------------------------------------------
0610 C       include 'epos.inc'
0611 C       double precision dranf
0612 C  1    rangen=sngl(dranf(dble(irandm)))
0613 C       if(rangen.le.0.)goto 1
0614 C       if(rangen.ge.1.)goto 1
0615 C       if(irandm.eq.1)write(ifch,*)'rangen()= ',rangen
0616
0617 C       return
0618 C       end
0619
0620 C c-----------------------------------------------------------------------
0621 C       double precision function drangen(dummy)
0622 C c-----------------------------------------------------------------------
0623 C c     generates a random number
0624 C c-----------------------------------------------------------------------
0625 C       include 'epos.inc'
0626 C       double precision dummy,dranf
0627 C       drangen=dranf(dummy)
0628 C       if(irandm.eq.1)write(ifch,*)'drangen()= ',drangen
0629
0630 C       return
0631 C       end
0632 C c-----------------------------------------------------------------------
0633 C       function cxrangen(dummy)
0634 C c-----------------------------------------------------------------------
0635 C c     generates a random number
0636 C c-----------------------------------------------------------------------
0637 C       include 'epos.inc'
0638 C       double precision dummy,dranf
0639 C       cxrangen=sngl(dranf(dummy))
0640 C       if(irandm.eq.1)write(ifch,*)'cxrangen()= ',cxrangen
0641
0642 C       return
0643 C       end
0644
0645
0646
0647 c Random number generator from CORSIKA *********************************
0648
0649
0650
0651
0652 C=======================================================================
0653
0654       DOUBLE PRECISION FUNCTION DRANF(dummy)
0655
0656 C-----------------------------------------------------------------------
0657 C  RAN(DOM  NUMBER) GEN(ERATOR) USED IN EPOS
0658 C  If calling this function within a DO-loop
0659 C  you should use an argument which prevents (dummy) to draw this function
0660 C  outside the loop by an optimizing compiler.
0661 C-----------------------------------------------------------------------
0662       implicit none
0663       common/eporansto2/irndmseq
0664       integer irndmseq
0665       double precision uni(1),dummy
0666 C-----------------------------------------------------------------------
0667
0668       call RMMARD( uni,1,irndmseq)
0669
0670       DRANF = UNI(1)
0671       UNI(1) = dummy        !to avoid warning
0672
0673       RETURN
0674       END
0675
0676
0677 c-----------------------------------------------------------------------
0678       subroutine ranfgt(seed)
0679 c-----------------------------------------------------------------------
0680 c Initialize seed in EPOS : read seed (output)
0681 c Since original output seed and EPOS seed are different,
0682 c define output seed as : seed=ISEED(3)*1E9+ISEED(2)
0683 c but only for printing. Important values stored in /eporansto/
0684 c Important : to be call before ranfst
0685 c-----------------------------------------------------------------------
0686       IMPLICIT NONE
0687       INTEGER          KSEQ
0688       PARAMETER        (KSEQ = 2)
0689       COMMON /CRRANMA3/CD,CINT,CM,TWOM24,TWOM48,MODCNS
0690       DOUBLE PRECISION CD,CINT,CM,TWOM24,TWOM48
0691       INTEGER          MODCNS
0692       COMMON /CRRANMA4/C,U,IJKL,I97,J97,NTOT,NTOT2,JSEQ
0693       DOUBLE PRECISION C(KSEQ),U(97,KSEQ)
0694       INTEGER          IJKL(KSEQ),I97(KSEQ),J97(KSEQ),
0695      *                 NTOT(KSEQ),NTOT2(KSEQ),JSEQ
0696       common/eporansto/diu0(100),iiseed(3)
0697       double precision    seed,diu0
0698       integer iiseed,i
0699
0700       iiseed(1)=IJKL(1)
0701       iiseed(2)=NTOT(1)
0702       iiseed(3)=NTOT2(1)
0703       seed=dble(iiseed(3))*dble(MODCNS)+dble(iiseed(2))
0704       diu0(1)=C(1)
0705       do i=2,98
0706         diu0(i)=U(i-1,1)
0707       enddo
0708       diu0(99)=dble(I97(1))
0709       diu0(100)=dble(J97(1))
0710       return
0711       end
0712
0713 c-----------------------------------------------------------------------
0714       subroutine ranfst(seed)
0715 c-----------------------------------------------------------------------
0716 c Initialize seed in EPOS :  restore seed (input)
0717 c Since original output seed and EPOS seed are different,
0718 c define output seed as : seed=ISEED(3)*1E9+ISEED(2)
0719 c but only for printing. Important values restored from /eporansto/
0720 c Important : to be call after ranfgt
0721 c-----------------------------------------------------------------------
0722       IMPLICIT NONE
0723       INTEGER          KSEQ
0724       PARAMETER        (KSEQ = 2)
0725       COMMON /CRRANMA3/CD,CINT,CM,TWOM24,TWOM48,MODCNS
0726       DOUBLE PRECISION CD,CINT,CM,TWOM24,TWOM48
0727       INTEGER          MODCNS
0728       COMMON /CRRANMA4/C,U,IJKL,I97,J97,NTOT,NTOT2,JSEQ
0729       DOUBLE PRECISION C(KSEQ),U(97,KSEQ)
0730       INTEGER          IJKL(KSEQ),I97(KSEQ),J97(KSEQ),
0731      *                 NTOT(KSEQ),NTOT2(KSEQ),JSEQ
0732       common/eporansto/diu0(100),iiseed(3)
0733       double precision    seedi,seed,diu0
0734       integer i,iiseed
0735
0736       seedi=seed
0737       IJKL(1)=iiseed(1)
0738       NTOT(1)=iiseed(2)
0739       NTOT2(1)=iiseed(3)
0740       C(1)=diu0(1)
0741       do i=2,98
0742         U(i-1,1)=diu0(i)
0743       enddo
0744       I97(1)=nint(diu0(99))
0745       J97(1)=nint(diu0(100))
0746       return
0747       end
0748
0749 c-----------------------------------------------------------------------
0750       subroutine ranflim(seed)
0751 c-----------------------------------------------------------------------
0752       double precision seed
0753       if(seed.gt.1d9)stop'seed larger than 1e9 not possible !'
0754       end
0755
0756 c-----------------------------------------------------------------------
0757       subroutine ranfcv(seed)
0758 c-----------------------------------------------------------------------
0759 c Convert input seed to EPOS random number seed
0760 c Since input seed and EPOS (from Corsika) seed are different,
0761 c define input seed as : seed=ISEED(3)*1E9+ISEED(2)
0762 c-----------------------------------------------------------------------
0763       IMPLICIT NONE
0764       COMMON /CRRANMA3/CD,CINT,CM,TWOM24,TWOM48,MODCNS
0765       DOUBLE PRECISION CD,CINT,CM,TWOM24,TWOM48
0766       INTEGER          MODCNS
0767       common/eporansto/diu0(100),iiseed(3)
0768       double precision    seed,diu0
0769       integer iiseed
0770
0771       iiseed(3)=nint(seed/dble(MODCNS))
0772       iiseed(2)=nint(mod(seed,dble(MODCNS)))
0773
0774       return
0775       end
0776
0777 c-----------------------------------------------------------------------
0778       subroutine ranfini(seed,iseq,iqq)
0779 c-----------------------------------------------------------------------
0780 c Initialize random number sequence iseq with seed
0781 c if iqq=-1, run first ini
0782 c    iqq=0 , set what sequence should be used
0783 c    iqq=1 , initialize sequence for initialization
0784 c    iqq=2 , initialize sequence for first event
0785 c-----------------------------------------------------------------------
0786       IMPLICIT NONE
0787       COMMON /CRRANMA3/CD,CINT,CM,TWOM24,TWOM48,MODCNS
0788       DOUBLE PRECISION CD,CINT,CM,TWOM24,TWOM48
0789       INTEGER          MODCNS
0790       common/eporansto/diu0(100),iiseed(3)
0791       double precision    seed,diu0
0792       integer iiseed
0793       common/eporansto2/irndmseq
0794       integer irndmseq
0795       integer iseed(3),iseq,iqq,iseqdum
0796
0797       if(iqq.eq.0)then
0798         irndmseq=iseq
0799       elseif(iqq.eq.-1)then
0800         iseqdum=0
0801         call RMMAQD(iseed,iseqdum,'R')   !first initialization
0802       elseif(iqq.eq.2)then
0803         irndmseq=iseq
0804         if(seed.ge.dble(MODCNS))then
0805            write(*,'(a,1p,e8.1)')'seedj larger than',dble(MODCNS)
0806            stop 'Forbidden !'
0807         endif
0808         iiseed(1)=nint(mod(seed,dble(MODCNS)))
0809 c iiseed(2) and iiseed(3) defined in aread
0810         call RMMAQD(iiseed,iseq,'S') !initialize random number generator
0811       elseif(iqq.eq.1)then        !dummy sequence for EPOS initialization
0812         irndmseq=iseq
0813         if(seed.ge.dble(MODCNS))then
0814            write(*,'(a,1p,e8.1)')'seedi larger than',dble(MODCNS)
0815            stop 'Forbidden !'
0816         endif
0817         iseed(1)=nint(mod(seed,dble(MODCNS)))
0818         iseed(2)=0
0819         iseed(3)=0
0820         call RMMAQD(iseed,iseq,'S') !initialize random number generator
0821       endif
0822       return
0823       end
0824
0825 C=======================================================================
0826
0827       SUBROUTINE RMMARD( RVEC,LENV,ISEQ )
0828
0829 C-----------------------------------------------------------------------
0830 C  C(ONE)X
0831 C  R(ANDO)M (NUMBER GENERATOR OF) MAR(SAGLIA TYPE) D(OUBLE PRECISION)
0832 C
0833 C  THESE ROUTINES (RMMARD,RMMAQD) ARE MODIFIED VERSIONS OF ROUTINES
0834 C  FROM THE CERN LIBRARIES. DESCRIPTION OF ALGORITHM SEE:
0835 C               http://consult.cern.ch/shortwrups/v113/top.html
0836 C  IT HAS BEEN CHECKED THAT RESULTS ARE BIT-IDENTICAL WITH CERN
0837 C  DOUBLE PRECISION RANDOM NUMBER GENERATOR RMM48, DESCRIBED IN
0838 C               http://consult.cern.ch/shortwrups/v116/top.html
0839 C  ARGUMENTS:
0840 C   RVEC   = DOUBLE PREC. VECTOR FIELD TO BE FILLED WITH RANDOM NUMBERS
0841 C   LENV   = LENGTH OF VECTOR (# OF RANDNUMBERS TO BE GENERATED)
0842 C   ISEQ   = # OF RANDOM SEQUENCE
0843 C
0844 C  VERSION OF D. HECK FOR DOUBLE PRECISION RANDOM NUMBERS.
0845 C  ADAPTATION  : T. PIEROG    IK  FZK KARLSRUHE FROM D. HECK VERSION
0846 C  DATE     : Feb  17, 2009
0847 C-----------------------------------------------------------------------
0848
0849       IMPLICIT NONE
0850       INTEGER          KSEQ
0851       PARAMETER        (KSEQ = 2)
0852       COMMON /CRRANMA3/CD,CINT,CM,TWOM24,TWOM48,MODCNS
0853       DOUBLE PRECISION CD,CINT,CM,TWOM24,TWOM48
0854       INTEGER          MODCNS
0855       COMMON /CRRANMA4/C,U,IJKL,I97,J97,NTOT,NTOT2,JSEQ
0856       DOUBLE PRECISION C(KSEQ),U(97,KSEQ),UNI
0857       INTEGER          IJKL(KSEQ),I97(KSEQ),J97(KSEQ),
0858      *                 NTOT(KSEQ),NTOT2(KSEQ),JSEQ
0859
0860       DOUBLE PRECISION RVEC(*)
0861       INTEGER          ISEQ,IVEC,LENV
0862       SAVE
0863
0864 C-----------------------------------------------------------------------
0865
0866       IF ( ISEQ .GT. 0  .AND.  ISEQ .LE. KSEQ ) JSEQ = ISEQ
0867
0868       DO   IVEC = 1, LENV
0869         UNI = U(I97(JSEQ),JSEQ) - U(J97(JSEQ),JSEQ)
0870         IF ( UNI .LT. 0.D0 ) UNI = UNI + 1.D0
0871         U(I97(JSEQ),JSEQ) = UNI
0872         I97(JSEQ)  = I97(JSEQ) - 1
0873         IF ( I97(JSEQ) .EQ. 0 ) I97(JSEQ) = 97
0874         J97(JSEQ)  = J97(JSEQ) - 1
0875         IF ( J97(JSEQ) .EQ. 0 ) J97(JSEQ) = 97
0876         C(JSEQ)    = C(JSEQ) - CD
0877         IF ( C(JSEQ) .LT. 0.D0 ) C(JSEQ)  = C(JSEQ) + CM
0878         UNI        = UNI - C(JSEQ)
0879         IF ( UNI .LT. 0.D0 ) UNI = UNI + 1.D0
0880 C  AN EXACT ZERO HERE IS VERY UNLIKELY, BUT LET'S BE SAFE.
0881         IF ( UNI .EQ. 0.D0 ) UNI = TWOM48
0882         RVEC(IVEC) = UNI
0883       ENDDO
0884
0885       NTOT(JSEQ) = NTOT(JSEQ) + LENV
0886       IF ( NTOT(JSEQ) .GE. MODCNS )  THEN
0887         NTOT2(JSEQ) = NTOT2(JSEQ) + 1
0888         NTOT(JSEQ)  = NTOT(JSEQ) - MODCNS
0889       ENDIF
0890
0891       RETURN
0892       END
0893
0894 C=======================================================================
0895
0896       SUBROUTINE RMMAQD( ISEED, ISEQ, CHOPT )
0897
0898 C-----------------------------------------------------------------------
0899 C  R(ANDO)M (NUMBER GENERATOR OF) MA(RSAGLIA TYPE INITIALIZATION) DOUBLE
0900 C
0901 C  SUBROUTINE FOR INITIALIZATION OF RMMARD
0902 C  THESE ROUTINE RMMAQD IS A MODIFIED VERSION OF ROUTINE RMMAQ FROM
0903 C  THE CERN LIBRARIES. DESCRIPTION OF ALGORITHM SEE:
0904 C               http://consult.cern.ch/shortwrups/v113/top.html
0905 C  FURTHER DETAILS SEE SUBR. RMMARD
0906 C  ARGUMENTS:
0907 C   ISEED  = SEED TO INITIALIZE A SEQUENCE (3 INTEGERS)
0908 C   ISEQ   = # OF RANDOM SEQUENCE
0909 C   CHOPT  = CHARACTER TO STEER INITIALIZE OPTIONS
0910 C
0911 C  CERN PROGLIB# V113    RMMAQ           .VERSION KERNFOR  1.0
0912 C  ORIG. 01/03/89 FCA + FJ
0913 C  ADAPTATION  : T. PIEROG    IK  FZK KARLSRUHE FROM D. HECK VERSION
0914 C  DATE     : Feb  17, 2009
0915 C-----------------------------------------------------------------------
0916
0917       IMPLICIT NONE
0918       INTEGER          KSEQ
0919       PARAMETER        (KSEQ = 2)
0920       COMMON /CRRANMA3/CD,CINT,CM,TWOM24,TWOM48,MODCNS
0921       DOUBLE PRECISION CD,CINT,CM,TWOM24,TWOM48
0922       INTEGER          MODCNS
0923       COMMON /CRRANMA4/C,U,IJKL,I97,J97,NTOT,NTOT2,JSEQ
0924       DOUBLE PRECISION C(KSEQ),U(97,KSEQ),UNI
0925       INTEGER          IJKL(KSEQ),I97(KSEQ),J97(KSEQ),
0926      *                 NTOT(KSEQ),NTOT2(KSEQ),JSEQ
0927
0928       DOUBLE PRECISION CC,S,T,UU(97)
0929       INTEGER          ISEED(3),I,IDUM,II,II97,IJ,IJ97,IORNDM,
0930      *                 ISEQ,J,JJ,K,KL,L,LOOP2,M,NITER
0931       CHARACTER        CHOPT*(*), CCHOPT*12
0932       LOGICAL          FIRST
0933       SAVE
0934       DATA             FIRST / .TRUE. /, IORNDM/11/, JSEQ/1/
0935
0936
0937 C-----------------------------------------------------------------------
0938
0939       IF ( FIRST ) THEN
0940         TWOM24 = 2.D0**(-24)
0941         TWOM48 = 2.D0**(-48)
0942         CD     = 7654321.D0*TWOM24
0943         CM     = 16777213.D0*TWOM24
0944         CINT   = 362436.D0*TWOM24
0945         MODCNS = 1000000000
0946         FIRST  = .FALSE.
0947         JSEQ   = 1
0948       ENDIF
0949       CCHOPT = CHOPT
0950       IF ( CCHOPT .EQ. ' ' ) THEN
0951         ISEED(1) = 54217137
0952         ISEED(2) = 0
0953         ISEED(3) = 0
0954         CCHOPT   = 'S'
0955         JSEQ     = 1
0956       ENDIF
0957
0958       IF     ( INDEX(CCHOPT,'S') .NE. 0 ) THEN
0959         IF ( ISEQ .GT. 0  .AND.  ISEQ .LE. KSEQ ) JSEQ = ISEQ
0960         IF ( INDEX(CCHOPT,'V') .NE. 0 ) THEN
0961           READ(IORNDM,'(3Z8)') IJKL(JSEQ),NTOT(JSEQ),NTOT2(JSEQ)
0962           READ(IORNDM,'(2Z8,Z16)') I97(JSEQ),J97(JSEQ),C(JSEQ)
0963           READ(IORNDM,'(24(4Z16,/),Z16)') U
0964           IJ = IJKL(JSEQ)/30082
0965           KL = IJKL(JSEQ) - 30082 * IJ
0966           I  = MOD(IJ/177, 177) + 2
0967           J  = MOD(IJ, 177)     + 2
0968           K  = MOD(KL/169, 178) + 1
0969           L  = MOD(KL, 169)
0970           CD =  7654321.D0 * TWOM24
0971           CM = 16777213.D0 * TWOM24
0972         ELSE
0973           IJKL(JSEQ)  = ISEED(1)
0974           NTOT(JSEQ)  = ISEED(2)
0975           NTOT2(JSEQ) = ISEED(3)
0976           IJ = IJKL(JSEQ) / 30082
0977           KL = IJKL(JSEQ) - 30082*IJ
0978           I  = MOD(IJ/177, 177) + 2
0979           J  = MOD(IJ, 177)     + 2
0980           K  = MOD(KL/169, 178) + 1
0981           L  = MOD(KL, 169)
0982           DO   II = 1, 97
0983             S = 0.D0
0984             T = 0.5D0
0985             DO   JJ = 1, 48
0986               M = MOD(MOD(I*J,179)*K, 179)
0987               I = J
0988               J = K
0989               K = M
0990               L = MOD(53*L+1, 169)
0991               IF ( MOD(L*M,64) .GE. 32 ) S = S + T
0992               T = 0.5D0 * T
0993             ENDDO
0994             UU(II) = S
0995           ENDDO
0996           CC    = CINT
0997           II97  = 97
0998           IJ97  = 33
0999 C  COMPLETE INITIALIZATION BY SKIPPING (NTOT2*MODCNS+NTOT) RANDOMNUMBERS
1000           NITER = MODCNS
1001           DO   LOOP2 = 1, NTOT2(JSEQ)+1
1002             IF ( LOOP2 .GT. NTOT2(JSEQ) ) NITER = NTOT(JSEQ)
1003             DO   IDUM = 1, NITER
1004               UNI = UU(II97) - UU(IJ97)
1005               IF ( UNI .LT. 0.D0 ) UNI = UNI + 1.D0
1006               UU(II97) = UNI
1007               II97     = II97 - 1
1008               IF ( II97 .EQ. 0 ) II97 = 97
1009               IJ97     = IJ97 - 1
1010               IF ( IJ97 .EQ. 0 ) IJ97 = 97
1011               CC       = CC - CD
1012               IF ( CC .LT. 0.D0 ) CC  = CC + CM
1013             ENDDO
1014           ENDDO
1015           I97(JSEQ) = II97
1016           J97(JSEQ) = IJ97
1017           C(JSEQ)   = CC
1018           DO   JJ = 1, 97
1019             U(JJ,JSEQ) = UU(JJ)
1020           ENDDO
1021         ENDIF
1022       ELSEIF ( INDEX(CCHOPT,'R') .NE. 0 ) THEN
1023         IF ( ISEQ .GT. 0 ) THEN
1024           JSEQ = ISEQ
1025         ELSE
1026           ISEQ = JSEQ
1027         ENDIF
1028         IF ( INDEX(CCHOPT,'V') .NE. 0 ) THEN
1029           WRITE(IORNDM,'(3Z8)') IJKL(JSEQ),NTOT(JSEQ),NTOT2(JSEQ)
1030           WRITE(IORNDM,'(2Z8,Z16)') I97(JSEQ),J97(JSEQ),C(JSEQ)
1031           WRITE(IORNDM,'(24(4Z16,/),Z16)') U
1032         ELSE
1033           ISEED(1) = IJKL(JSEQ)
1034           ISEED(2) = NTOT(JSEQ)
1035           ISEED(3) = NTOT2(JSEQ)
1036         ENDIF
1037       ENDIF
1038
1039       RETURN
1040       END