ReggeGribovPartonMCInterface/src/sibyll_epos.f

0001 c 15.02.2004 Link routines between SIBYLL2.1 and EPOS.
0002 c author T. Pierog
0003
0004
0005 c-----------------------------------------------------------------------
0006       subroutine IniSibyll
0007 c-----------------------------------------------------------------------
0008 c Primary initialization for Sibyll
0009 c-----------------------------------------------------------------------
0010       include 'epos.inc'
0011       COMMON /S_DEBUG/ Ncall, Ndebug
0012       COMMON /S_CSYDEC/ CBR(102), KDEC(612), LBARP(49), IDB(49)
0013
0014       call utpri('inisib',ish,ishini,6)
0015       write(ifmt,'(a,i6)')'initialize Sibyll ...'
0016
0017       Ndebug=0
0018       if(ish.ge.3)Ndebug=ish-2
0019
0020 C... SIBYLL initialization
0021       CALL SIBYLL_INI
0022
0023 C...Cross sections for nucleus-nucleus and hadron nucleus
0024       CALL NUC_NUC_INI
0025
0026 C...define all particles as unstable
0027       do i=1,49
0028         IDB(i) = abs(IDB(i))   ! >0 means unstable
0029       enddo
0030
0031
0032       egymin=10.1
0033       egymax=1e7
0034       irescl=0
0035
0036
0037       call utprix('inisib',ish,ishini,6)
0038       end
0039
0040 c-----------------------------------------------------------------------
0041       subroutine IniEvtSib
0042 c-----------------------------------------------------------------------
0043 c Initialization for each type of event (for given proj, targ and egy)
0044 c-----------------------------------------------------------------------
0045       include 'epos.inc'
0046       COMMON /S_CSYDEC/ CBR(102), KDEC(612), LBARP(49), IDB(49)
0047       common/geom/rmproj,rmtarg,bmax,bkmx
0048
0049
0050       if(matarg.gt.18.or.maproj.gt.64)
0051      &  call utstop('Mass too big for Sibyll (Mtrg<18, Mprj<64) !&')
0052       id=idtrafo('nxs','sib',idproj)
0053       ida=abs(id)
0054       if(ida.lt.6.or.ida.gt.14)
0055      &  call utstop('projectile no allowed in Sibyll !&')
0056       if(idtarg.ne.0.and.idtarg.ne.1120)
0057      &  call utstop('target no allowed in Sibyll !&')
0058       if(bminim.gt.0.or.bmaxim.lt.1000)
0059      &  write(ifmt,*)'Only min bias event in Sibyll ... no b !'
0060       call iclass(idproj,iclpro)
0061       bminim=0.
0062       bmaxim=10000.
0063       bmax=10.+maproj+matarg
0064       sibincs=fsibcrse(engy,maproj,matarg)
0065       if(engy.lt.egymin)sibincs=0.          !below egymin, no interaction
0066       call xsigma
0067
0068
0069 c Epos syntax to allow (or not) particle decay in Sibyll
0070 c (part taken from epos-dky: hdecas)
0071
0072       if(idecay.eq.1.and.ctaumin.le.0.)then
0073
0074       if(ndecay.eq.1.or.mod(ndecay/10,10).eq.1)then          !Kshort, Klong
0075         IDB(11) = -abs(IDB(11))   ! <0 means stable
0076         IDB(12) = -abs(IDB(12))
0077       endif
0078       if(ndecay.eq.1.or.mod(ndecay/100,10).eq.1)then         !Lambda
0079         IDB(39) = -abs(IDB(39))
0080       endif
0081       if(ndecay.eq.1.or.mod(ndecay/1000,10).eq.1)then        !sigma+
0082         IDB(34) = -abs(IDB(34))
0083       endif
0084       if(ndecay.eq.1.or.mod(ndecay/1000,10).eq.1)then        !sigma-
0085         IDB(36) = -abs(IDB(36))
0086       endif
0087       if(ndecay.eq.1.or.mod(ndecay/10000,10).eq.1)then       !Xi+/-
0088         IDB(38) = -abs(IDB(38))
0089       endif
0090       if(ndecay.eq.1.or.mod(ndecay/10000,10).eq.1)then       !Xi0
0091         IDB(37) = -abs(IDB(37))
0092       endif
0093       if(ndecay.eq.1.or.mod(ndecay/100000 ,10).eq.1)then     !omega
0094         IDB(49) = -abs(IDB(49))
0095       endif
0096       if(ndecay.eq.1.or.mod(ndecay/1000000,10).eq.1)then     !pi0
0097         IDB(6) = -abs(IDB(6))
0098       endif
0099
0100       if(nrnody.gt.0)then                      !all other particle
0101         do nod=1,nrnody
0102           idd=abs(idtrafo('nxs','sib',nody(nod)))
0103           if(idd.lt.50)IDB(idd) = -abs(IDB(idd))
0104         enddo
0105       endif
0106
0107       else
0108
0109 C...define all particles as stable
0110       do i=1,49
0111         IDB(i) = -abs(IDB(i))   ! <0 means stable
0112       enddo
0113
0114       endif
0115
0116       if(ish.ge.2)write(ifch,*)
0117      &  'Sibyll used with (E,proj,maproj,matarg)',engy,id,maproj
0118      &  ,matarg
0119
0120       return
0121       end
0122
0123 c-----------------------------------------------------------------------
0124       subroutine emssib(iret)
0125 c-----------------------------------------------------------------------
0126 c  call Sibyll to simulate interaction
0127 c-----------------------------------------------------------------------
0128       include 'epos.inc'
0129       common/geom/rmproj,rmtarg,bmax,bkmx
0130 C  SIBYLL
0131       COMMON /S_PLIST/ P(8000,5), LLIST(8000), NP
0132       COMMON /S_PLNUC/ PA(5,40000), LLA(40000), NPA
0133       PARAMETER (NW_max = 20)
0134       PARAMETER (NS_max = 20, NH_max = 50)
0135       PARAMETER (NJ_max = (NS_max+NH_max)*NW_max)
0136       COMMON /S_CHIST/ X1J(NJ_max),X2J(NJ_max),
0137      &    X1JSUM(NW_max),X2JSUM(NW_max),PTJET(NJ_max),PHIJET(NJ_max),
0138      &    NNPJET(NJ_max),NNPSTR(2*NW_max),NNSOF(NW_max),NNJET(NW_max),
0139      &    JDIF(NW_max),NW,NJET,NSOF
0140
0141       iret=0
0142       b1=bminim
0143       b2=min(bmax,bmaxim)
0144       a=pi*(b2**2-b1**2)
0145
0146       if(a.gt.0..and.rangen().gt.sibincs/10./a)goto 1001   !no interaction
0147       if(ish.ge.3)call alist('Determine Sibyll Production&',0,0)
0148
0149       nptl=0
0150       NP=0
0151       NPA=0
0152
0153       nevt=1
0154       kolevt=-1
0155       koievt=-1
0156       kohevt=-1
0157       npjevt=maproj
0158       ntgevt=matarg
0159       pmxevt=pnll
0160       egyevt=engy
0161       bimevt=0.
0162       bimp=0.
0163       phievt=0.
0164       phi=0.
0165       anintine=anintine+1.
0166       ns=0                      !number of projectile spectators
0167       npps=0
0168       npns=0
0169
0170       call conre
0171       call conwr
0172
0173       itrg=matarg
0174       if(idtargin.eq.0)itrg=0
0175       if(maproj.eq.1)then             !hadronic projectile
0176         L0=idtrafo('nxs','sib',idproj)
0177         CALL SIBYLL (L0, itrg, engy)
0178         CALL DECSIB
0179         if(ish.ge.5)write(ifch,'(a,i5)')
0180      $         ' number of particles from Sibyll :',NP
0181 c save interaction type
0182         if(JDIF(1).eq.0)then
0183           typevt=1                !ND
0184         elseif(JDIF(1).eq.3)then
0185           typevt=2                !DD
0186         else
0187           typevt=4                !SD
0188         endif
0189         do k=1,NP
0190
0191 c LLIST is the code of final particle, P - its 4-momentum and mass.
0192           ic=LLIST(k)
0193
0194           if(ish.ge.7)write(ifch,'(a,i5,a,i5,2a,4(e10.4,1x))')
0195      $       ' Sibyll particle ',k,' id :',ic,' before conversion'
0196      $     , ' momentum :',(P(k,i),i=1,5)
0197
0198           nptl=nptl+1
0199           if(nptl.gt.mxptl)call utstop('Sibyll: mxptl too small&')
0200
0201           if(abs(ic).ge.10000)then
0202             ic=ic-sign(10000,ic)
0203             istptl(nptl)=1
0204           else
0205             istptl(nptl)=0
0206           endif
0207
0208           id=idtrafo('sib','nxs',ic)
0209           if(ish.ge.7)write(ifch,'(a,i5,a,i5,a)')
0210      $         ' epos particle ',nptl,' id :',id,' after conversion'
0211
0212
0213           pptl(1,nptl)=P(k,1)   !P_x
0214           pptl(2,nptl)=P(k,2)   !P_y
0215           pptl(3,nptl)=P(k,3)   !P_z
0216           pptl(4,nptl)=abs(P(k,4))   !E
0217           pptl(5,nptl)=P(k,5)   !mass
0218           ityptl(nptl)=0
0219           iorptl(nptl)=1
0220           jorptl(nptl)=maproj+matarg
0221           ifrptl(1,nptl)=0
0222           ifrptl(2,nptl)=0
0223           xorptl(1,nptl)=0.
0224           xorptl(2,nptl)=0.
0225           xorptl(3,nptl)=0.
0226           xorptl(4,nptl)=0.
0227           tivptl(1,nptl)=0.
0228           tivptl(2,nptl)=0.
0229           idptl(nptl)=id
0230
0231
0232           if(ish.ge.5)write(ifch,'(a,i5,a,i5,a,4(e10.4,1x),f6.3)')
0233      $         ' particle from Sibyll ',nptl,' id :',idptl(nptl)
0234      $         , ' momentum :',(pptl(i,nptl),i=1,5)
0235
0236
0237         enddo
0238         if(NW.lt.matarg)then
0239           ntgevt=NW
0240           do is=maproj+1+NW,maproj+matarg !make the ns last target nucleon final
0241             iorptl(is)=0
0242             istptl(is)=0
0243           enddo
0244         endif
0245       else                          !for nucleus projectile
0246         nbar=0
0247         IAP = maproj
0248         CALL SIBNUC (IAP, itrg, engy)
0249         if(ish.ge.5)write(ifch,'(a,i5)')
0250      $         ' number of particles from Sibyll :',NPA
0251         do 100 k=1,NPA
0252
0253 c LLIST is the code of final particle, P - its 4-momentum and mass.
0254           ic=LLA(k)
0255
0256           if(ish.ge.7)write(ifch,'(a,i5,a,i5,2a,4(e10.4,1x))')
0257      $       ' Sibyll particle ',k,' id :',ic,' before conversion'
0258      $     , ' momentum :',(PA(i,k),i=1,5)
0259
0260
0261           nNuc=0
0262           if(ic.ge.1001) then                !count spectators
0263             nNuc=ic-1000
0264             if(infragm.le.1
0265      &         .or.PA(4,k).lt.0.5*egymin)then   !nuclear interaction only above min energy, otherwise : fragmentation
0266               ns=ns+nNuc
0267               goto 100
0268             elseif(ic.eq.1001)then
0269               if(drangen(dummy).lt.0.45d0) then
0270                 ic = 13
0271                 npps=npps+1
0272               else
0273                 ic = 14
0274                 npns=npns+1
0275                endif
0276               nNuc=0
0277             else
0278               ptm=sqrt(PA(1,k)*PA(1,k)+PA(2,k)*PA(2,k)+PA(5,k)*PA(5,k))
0279               PA(4,k)=PA(4,k)*float(nNuc)            !energy by nucleon
0280               PA(3,k)=sign(sqrt((PA(4,k)+ptm)*(PA(4,k)-ptm)),PA(3,k))
0281               npps=npps+nNuc/2
0282               npns=npns+nNuc-nNuc/2
0283             endif
0284           endif
0285           nptl=nptl+1
0286           if(nptl.gt.mxptl)call utstop('Sibyll: mxptl too small&')
0287           id=idtrafo('sib','nxs',ic)
0288           if(ish.ge.7)write(ifch,'(a,i5,a,i10,a)')
0289      $         ' epos particle ',nptl,' id :',id,' after conversion'
0290
0291
0292           nbar=nbar+nNuc
0293           if(abs(id).gt.1000.and.nNuc.eq.0)nbar=nbar+sign(1,id)
0294           pptl(1,nptl)=PA(1,k)  !P_x
0295           pptl(2,nptl)=PA(2,k)  !P_y
0296           pptl(3,nptl)=PA(3,k)  !P_z
0297           pptl(4,nptl)=PA(4,k)  !E
0298           pptl(5,nptl)=PA(5,k)  !mass
0299           istptl(nptl)=0
0300           ityptl(nptl)=0
0301           iorptl(nptl)=1
0302           jorptl(nptl)=maproj+matarg
0303           ifrptl(1,nptl)=0
0304           ifrptl(2,nptl)=0
0305           xorptl(1,nptl)=0.
0306           xorptl(2,nptl)=0.
0307           xorptl(3,nptl)=0.
0308           xorptl(4,nptl)=0.
0309           tivptl(1,nptl)=0.
0310           tivptl(2,nptl)=0.
0311           idptl(nptl)=id
0312
0313
0314           if(ish.ge.5)write(ifch,'(a,i5,a,i10,a,4(e10.4,1x),f6.3)')
0315      $         ' particle from Sibyll ',nptl,' id :',idptl(nptl)
0316      $         , ' momentum :',(pptl(i,nptl),i=1,5)
0317
0318  100    continue
0319
0320         ntw=nbar-(maproj-ns)
0321         nsf=0
0322         if(ntw.lt.matarg)then
0323           ntgevt=ntw
0324           do is=maproj+1+ntw,maproj+matarg !make the nts last target nucleon final (not wounded)
0325             iorptl(is)=0
0326             istptl(is)=0
0327           enddo
0328         else
0329           nsf=maproj+matarg-nbar
0330         endif
0331           if(ish.ge.5)write(ifch,'(a,i3,a,i3,a,i2,a)')
0332      $         ' target spectators :',matarg-ntw
0333      $        ,' projectile spectators (ns) :',nsf,' (',ns,')'
0334         if((infragm.le.1.or.ns.gt.0).and.nsf.le.maproj)then
0335           if(infragm.eq.2)ns=ns-nsf
0336           if(infragm.eq.1.and.ns.gt.0)then
0337 c  remaining nucleus is one fragment
0338             nptl=nptl+1
0339             istptl(nptl)=0
0340             pptl(1,nptl)=0.d0
0341             pptl(2,nptl)=0.d0
0342             pptl(4,nptl)=0.d0
0343             inucl=0
0344             do is=1,ns
0345               inucl=inucl+1
0346               pptl(4,nptl)=pptl(4,nptl)+pptl(4,is)
0347             enddo
0348             iprot= int(dble(inucl) / 2.15d0 + 0.7d0)
0349             idnucl=1000000000+iprot*10000+inucl*10 !code for nuclei
0350             npps=npps+iprot
0351             npns=npns+inucl-iprot
0352             call idmass(idnucl,am)
0353             pptl(5,nptl)=am  !mass
0354             ptot=(pptl(4,nptl)+am)*(pptl(4,nptl)-am)
0355             pptl(3,nptl)=sqrt(ptot)
0356             ityptl(nptl)=0
0357             istptl(nptl)=0
0358             iorptl(nptl)=1
0359             jorptl(nptl)=maproj
0360             ifrptl(1,nptl)=0
0361             ifrptl(2,nptl)=0
0362             xorptl(1,nptl)=xorptl(1,1)
0363             xorptl(2,nptl)=xorptl(2,1)
0364             xorptl(3,nptl)=xorptl(3,1)
0365             xorptl(4,nptl)=xorptl(4,1)
0366             tivptl(1,nptl)=tivptl(1,1)
0367             tivptl(2,nptl)=tivptl(2,1)
0368             idptl(nptl)=idnucl
0369           else
0370             do is=maproj-ns+1,maproj         !make the nsf first projectile nucleon final (not wounded)
0371               if(infragm.eq.2)istptl(is)=0
0372               if(idptl(is).eq.1120)npps=npps+1
0373               if(idptl(is).eq.1220)npns=npns+1
0374             enddo
0375           endif
0376         endif
0377 c number of participants
0378         if(laproj.gt.1)then
0379         npjevt=maproj-npps-npns
0380         npppar=max(0,laproj-npps)
0381         npnpar=npjevt-npppar
0382 c set participant projectile as non spectators
0383         do i=1,maproj
0384           if(idptl(i).eq.1120)then
0385             if(npppar.gt.0)then
0386               npppar=npppar-1
0387             else                !restore spectators
0388               iorptl(i)=0
0389               if(infragm.eq.0)istptl(i)=0
0390             endif
0391           endif
0392           if(idptl(i).eq.1220)then
0393             if(npnpar.gt.0)then
0394               npnpar=npnpar-1
0395             else                !restore spectators
0396               iorptl(i)=0
0397               if(infragm.eq.0)istptl(i)=0
0398             endif
0399           endif
0400         enddo
0401       endif
0402
0403       endif
0404
0405
0406 1000  return
0407
0408 1001  iret=-1
0409       goto 1000
0410
0411       end
0412
0413
0414 c------------------------------------------------------------------------------
0415       function fsibcrse(egy,mapro,matar)
0416 c------------------------------------------------------------------------------
0417 c hadron-proton particle production cross section with Sibyll.
0418 c egy - center of mass energy
0419 c------------------------------------------------------------------------------
0420       include 'epos.inc'
0421       dimension SDIF(3)
0422
0423       if(iclpro.eq.1)then
0424         L=2
0425       elseif(iclpro.eq.2)then
0426         L=1
0427       else
0428         L=3
0429       endif
0430       call SIB_SIGMA_HP(L,egy,ST,SEL,SINEL,SDIF,SL,RHO)
0431       if(matar.gt.1)then
0432 C  calculate hadron-A(matar) cross section
0433         CALL GLAUBER(matar,ST,SL,RHO,SIGTA,SIGELAdum,SIGQEA)
0434         fsibcrse=SIGTA-SIGQEA
0435       else
0436         fsibcrse=SINEL
0437       endif
0438
0439       if(mapro.gt.1)fsibcrse=ainfin !???????? temporary
0440
0441       return
0442       end
0443
0444 c------------------------------------------------------------------------------
0445       double precision function sibcrse(ek,mapro,matar,id)
0446 c------------------------------------------------------------------------------
0447 c inelastic cross section of Sibyll
0448 c if id=0, target = air
0449 c ek - kinetic energy in GeV
0450 c maproj - projec mass number     (1<maproj<64)
0451 c matarg - projec mass number
0452 c id - proj id (sibyll code)
0453 c------------------------------------------------------------------------------
0454       include 'epos.inc'
0455       double precision egy
0456       COMMON /CLENNN/ SSIGNUC(60), ALNUC(60)
0457
0458       sibcrse=0.d0
0459       call idmass(1120,amt1)
0460       call idmass(1220,amt2)
0461       amtar=0.5*(amt1+amt2)
0462       if(mapro.eq.1)call idmass(idproj,ampro)
0463       if(matar.eq.1)call idmass(idtarg,amtar)
0464       egy=dble(ek+ampro)
0465       if(id.eq.0)then
0466         if(maproj.eq.1)then
0467           sqs=sqrt( 2*real(egy)*amtar+amtar**2+ampro**2 )
0468           if(iclpro.eq.1)then
0469             L=2
0470           elseif(iclpro.eq.2)then
0471             L=1
0472           else
0473             L=3
0474           endif
0475           call SIB_SIGMA_HAIR (L,sqs,sibcr)
0476           sibcrse=dble(sibcr)
0477         else
0478           EO=real(egy)*1.e-3         !e0 in TeV
0479           CALL  SIGNUC_INI(mapro,E0) !  fills SSIGNUC and ALNUC
0480           sibcrse  = dble(SSIGNUC(mapro))
0481         endif
0482       else
0483         sqs=sqrt( 2*real(egy)*amtar+amtar**2+ampro**2 )
0484         sibcrse=dble(fsibcrse(sqs,mapro,matar))
0485       endif
0486
0487       return
0488       end
0489
0490 c--------------------------------------------------------------------
0491       function S_RNDM(idum)
0492 c--------------------------------------------------------------------
0493 c random number generator
0494 c--------------------------------------------------------------------
0495       include 'epos.inc'
0496
0497       S_RNDM=rangen()
0498       if(irandm.eq.1)write(ifch,*)'S_RNDM()= ',S_RNDM,idum
0499
0500       return
0501       end
0502
0503
0504