ReggeGribovPartonMCInterface/src/qgsjetII_epos.f

0001 c 07.06.2006 Link routines between QGSJet-II-03 and EPOS.
0002 c author T. Pierog
0003
0004
0005 c-----------------------------------------------------------------------
0006       subroutine IniQGSJETII
0007 c-----------------------------------------------------------------------
0008 c Primary initialization for QGSJET-II
0009 c-----------------------------------------------------------------------
0010       parameter(iapmax=208)
0011       include 'epos.inc'
0012       integer debug
0013       common /qgdebug/ debug
0014       common /qgarr43/ moniou
0015       double precision bqgs2,bmaxqgs,bmaxnex,bminnex,xan,xbn
0016       common /qgsIInex1/xan(iapmax,3),xbn(iapmax,3)
0017      *,bqgs2,bmaxqgs,bmaxnex,bminnex
0018       CHARACTER DATDIR*(132)
0019
0020       call utpri('iniqgsII',ish,ishini,6)
0021       write(ifmt,'(a,i6)')'initialize QGSJET-II ...'
0022
0023       debug=0
0024       if(ish.ge.3)debug=ish-2
0025       moniou=ifch
0026
0027 c model parameter setting
0028       call qgset
0029 c common initialization procedure
0030       DATDIR="qgsjetII"
0031       call qgaini(DATDIR)
0032       BMAXNEX=dble(bmaxim)
0033       BMINNEX=dble(bminim)
0034       egymin=0.1
0035       egymax=egymax
0036       irescl=0
0037
0038
0039       call utprix('iniqgsII',ish,ishini,6)
0040       end
0041
0042 c-----------------------------------------------------------------------
0043       subroutine IniEvtQGSII
0044 c-----------------------------------------------------------------------
0045 c Initialization for each type of event (for given proj, targ and egy)
0046 c-----------------------------------------------------------------------
0047       parameter(iapmax=208)
0048       include 'epos.inc'
0049       common/geom/rmproj,rmtarg,bmax,bkmx
0050 c QGSJET-II Common
0051       double precision bqgs,bmaxqgs,bmaxnex,bminnex,xa,xb,e0
0052       common /qgsIInex1/xa(iapmax,3),xb(iapmax,3)
0053      *,bqgs,bmaxqgs,bmaxnex,bminnex
0054
0055
0056       if(matarg.gt.iapmax.or.maproj.gt.iapmax)
0057      &  call utstop('Nucleus too big for QGSJET-II (Mmax=208) !&')
0058       call iclass(idproj,iclpro)
0059       call iclass(idtarg,icltar)
0060       icp=idtrafo('nxs','qgs',idproj)
0061       if(icp.eq.0)icp=1-2*int(rangen()+0.5)       !pi0=pi+ or p-
0062       if(abs(icp).gt.5)
0063      &  call utstop('Projectile not allowed in QGSJET-II !&')
0064       e0=dble(elab)
0065       call qgini(e0,icp,maproj,matarg)
0066       call qgini(e0,icp,maproj,matarg)        !again to set bm properly
0067       bmax=BMAXQGS
0068       qgsIIincs=qgsIIcrse(ekin,maproj,matarg,idtarg)
0069       if(engy.lt.egymin)qgsIIincs=0.          !below egymin, no interaction
0070       call xsigma                             !change bm in qgfz
0071       bkmx=sqrt(sigine/10./pi)        !10= fm^2 -> mb
0072       if(ish.ge.2)write(ifch,*)
0073      &  'QGSJET-II used with (E,proj,maproj,matarg,bmax)',e0,icp
0074      &  ,maproj,matarg,bmax
0075
0076       return
0077       end
0078
0079 c-----------------------------------------------------------------------
0080       subroutine emsqgsII(iret)
0081 c-----------------------------------------------------------------------
0082 c  call qgsjet-II to simulate interaction
0083 c-----------------------------------------------------------------------
0084       parameter(iapmax=208,nptmax=95000)
0085       include 'epos.inc'
0086       include 'epos.incems'
0087       common/geom/rmproj,rmtarg,bmax,bkmx
0088       double precision bqgs,bmaxqgs,bmaxnex,bminnex,xa,xb,esp
0089       common /qgsIInex1/xa(iapmax,3),xb(iapmax,3)
0090      *,bqgs,bmaxqgs,bmaxnex,bminnex
0091       common/nucl3/phi,bimp
0092       common /qgarr12/ nsp
0093       common /qgarr14/ esp(4,nptmax),ich(nptmax)
0094 c nsf - number of secondary fragments;
0095 c iaf(i) - mass of the i-th fragment
0096       common /qgarr13/ nsf,iaf(iapmax)
0097       common /qgarr55/ nwt,nwp
0098
0099       iret=0
0100       b1=bminim
0101       b2=min(bmax,bmaxim)
0102       a=pi*(b2**2-b1**2)
0103
0104       if(a.gt.0..and.rangen().gt.qgsIIincs/10./a)goto 1001   !no interaction
0105       if(ish.ge.3)call alist('Determine QGSJET-II Production&',0,0)
0106
0107       nptl=0
0108       nsp=0
0109       call qgconf
0110
0111       nevt=1
0112       kolevt=-1
0113       koievt=-1
0114       kohevt=-1
0115       npjevt=maproj
0116       ntgevt=matarg
0117       pmxevt=pnll
0118       egyevt=engy
0119       if(BQGS.ge.0.d0)then
0120         bimevt=real(BQGS)
0121         bimp=real(BQGS)
0122         phi=2.*pi*rangen()
0123         phievt=phi
0124       else
0125         bimevt=0.
0126         bimp=0.
0127         phievt=0.
0128         phi=0.
0129       endif
0130       anintine=anintine+1.
0131
0132       call conre
0133       call conwr
0134       call conqgsII
0135
0136
0137
0138 c keep the projectile spectators as fragments
0139       npns=0
0140       npps=0
0141       ns=0
0142       if(infragm.eq.2)then
0143
0144         if(NSF.gt.0)then
0145           do is=1,NSF           !count the number of spectators
0146           if(ish.ge.7)write(ifch,'(a,i5,a,i5)')
0147      $       ' Projecticle Fragment ',is,' Mass :',IAF(is)
0148             nptl=nptl+1
0149             istptl(nptl)=0
0150             if(IAF(is).eq.1)then
0151               id=idptl(is)
0152               pptl(3,nptl)=pptl(3,is)
0153               pptl(4,nptl)=pptl(4,is)
0154               pptl(5,nptl)=pptl(5,is)
0155               if(id.eq.1120)npps=npps+1
0156               if(id.eq.1220)npns=npns+1
0157             else
0158               if(IAF(is).eq.2)then
0159                 id=17
0160                 npps=npps+1
0161                 npns=npns+1
0162               elseif(IAF(is).eq.3)then
0163                 id=18
0164                 npps=npps+1
0165                 npns=npns+2
0166               elseif(IAF(is).eq.4)then
0167                 id=19
0168                 npps=npps+2
0169                 npns=npns+2
0170               else
0171                 inucl=IAF(is)
0172                 iprot= int(dble(inucl) / 2.15d0 + 0.7d0)
0173                 id=1000000000+iprot*10000+inucl*10 !code for nuclei
0174                 npps=npps+iprot
0175                 npns=npns+inucl-iprot
0176               endif
0177               call idmass(id,am)
0178               pptl(4,nptl)=dble(IAF(is))*pptl(4,is)      !Etot
0179               pptl(5,nptl)=am                             !mass
0180               pz2tmp=(pptl(4,nptl)+am)*(pptl(4,nptl)-am)
0181               if(pz2tmp.gt.0.d0)then
0182                 pptl(3,nptl)=sqrt(pz2tmp)                 !Pz
0183               else
0184                 write(*,*)'Warning in emsqgsII !'
0185                 write(*,*)'energy of fragment too small :',IAF(is),am
0186      &                     ,pptl(4,nptl)
0187                 pptl(3,nptl)=pptl(4,nptl)
0188               endif
0189             endif
0190             pptl(1,nptl)=0.d0 !P_x
0191             pptl(2,nptl)=0.d0 !P_y
0192             ityptl(nptl)=0
0193             iorptl(nptl)=1
0194             jorptl(nptl)=maproj+matarg
0195             ifrptl(1,nptl)=0
0196             ifrptl(2,nptl)=0
0197             xorptl(1,nptl)=0.d0
0198             xorptl(2,nptl)=0.d0
0199             xorptl(3,nptl)=0.d0
0200             xorptl(4,nptl)=0.d0
0201             tivptl(1,nptl)=0.d0
0202             tivptl(2,nptl)=0.d0
0203             idptl(nptl)=id
0204             if(ish.ge.5)write(ifch,'(a,i5,a,i5,a,4(e10.4,1x),f6.3)')
0205      $       ' Fragment from qgsjetII ',nptl,' id :',idptl(nptl)
0206      $  , ' momentum :',(pptl(k,nptl),k=1,5)
0207
0208           enddo
0209         endif
0210
0211 c make the projectile spectators as free nucleons
0212
0213       else
0214         ns=0
0215         if(NSF.gt.0)then
0216           do is=1,NSF           !count the number of spectators
0217             ns=ns+IAF(is)
0218           enddo
0219           if(infragm.eq.1)then
0220 c  remaining nucleus is one fragment
0221             nptl=nptl+1
0222             istptl(nptl)=0
0223             pptl(1,nptl)=0.d0
0224             pptl(2,nptl)=0.d0
0225             pptl(4,nptl)=0.d0
0226             inucl=0
0227             do is=1,ns
0228               inucl=inucl+1
0229               pptl(4,nptl)=pptl(4,nptl)+pptl(4,is)
0230             enddo
0231             iprot= int(dble(inucl) / 2.15d0 + 0.7d0)
0232             idnucl=1000000000+iprot*10000+inucl*10 !code for nuclei
0233             npps=npps+iprot
0234             npns=npns+inucl-iprot
0235             call idmass(idnucl,am)
0236             pptl(5,nptl)=am  !mass
0237             ptot=(pptl(4,nptl)+am)*(pptl(4,nptl)-am)
0238             pptl(3,nptl)=sqrt(ptot)
0239             ityptl(nptl)=0
0240             istptl(nptl)=0
0241             iorptl(nptl)=1
0242             jorptl(nptl)=maproj
0243             ifrptl(1,nptl)=0
0244             ifrptl(2,nptl)=0
0245             xorptl(1,nptl)=xorptl(1,1)
0246             xorptl(2,nptl)=xorptl(2,1)
0247             xorptl(3,nptl)=xorptl(3,1)
0248             xorptl(4,nptl)=xorptl(4,1)
0249             tivptl(1,nptl)=tivptl(1,1)
0250             tivptl(2,nptl)=tivptl(2,1)
0251             idptl(nptl)=idnucl
0252           else
0253             do is=maproj-ns+1,maproj          !make the ns last projectile nucleon final
0254               if(idptl(is).eq.1120)npps=npps+1
0255               if(idptl(is).eq.1220)npns=npns+1
0256             enddo
0257           endif
0258         endif
0259       endif
0260
0261 c number of participants
0262       if(laproj.gt.1)then
0263         npjevt=maproj-npps-npns
0264         npppar=max(0,laproj-npps)
0265         npnpar=npjevt-npppar
0266 c set participant projectile as non spectators
0267         do i=1,maproj
0268           if(idptl(i).eq.1120)then
0269             if(npppar.gt.0)then
0270               npppar=npppar-1
0271             else                !restore spectators
0272               iorptl(i)=0
0273               if(infragm.eq.0)istptl(i)=0
0274             endif
0275           endif
0276           if(idptl(i).eq.1220)then
0277             if(npnpar.gt.0)then
0278               npnpar=npnpar-1
0279             else                !restore spectators
0280               iorptl(i)=0
0281               if(infragm.eq.0)istptl(i)=0
0282             endif
0283           endif
0284         enddo
0285       endif
0286
0287 c restore target spectators
0288       ns=0
0289       if(NWT.lt.matarg)then
0290 c number of participants
0291         ntgevt=NWT
0292         do is=ntgevt+1,matarg         !make the last target nucleon final
0293           iorptl(maproj+is)=0
0294           istptl(maproj+is)=0
0295         enddo
0296       endif
0297
0298
0299       do is=1,nsp
0300
0301 c ich is the type of secondary hadron, esp - its transverse momentum,
0302 c and its energy
0303 c the following notations for the particles types are used: 0 - pi0, 1 -
0304 c pi+,
0305 c -1 - pi-, 2 - p, -2 - p, 3 - n, -3 - n, 4 - k+, -4 - k-, 5 - k0s, -5 -
0306 c k0l
0307           ic=ich(is)
0308           if(ish.ge.7)write(ifch,'(a,i5,a,i5,2a,4(e10.4,1x))')
0309      $       ' qgsjet particle ',is,' id :',ic,' before conversion'
0310      $     , ' momentum :',(esp(k,is),k=1,4)
0311
0312             nptl=nptl+1
0313             if(nptl.gt.mxptl)call utstop('qgsjet: mxptl too small&')
0314             id=idtrafo('qgs','nxs',ic)
0315             if(ish.ge.7)write(ifch,'(a,i5,a,i5,a)')
0316      $       ' epos particle ',nptl,' id :',id,' after conversion'
0317             call idmass(id,am)
0318
0319
0320             pptl(1,nptl)=real(esp(3,is)) !P_x
0321             pptl(2,nptl)=real(esp(4,is)) !P_y
0322             pptl(3,nptl)=real(esp(2,is)) !P_z
0323             pptl(4,nptl)=real(esp(1,is)) !E
0324             pptl(5,nptl)=am              !mass
0325             istptl(nptl)=0
0326             ityptl(nptl)=0
0327             iorptl(nptl)=1
0328             jorptl(nptl)=maproj+matarg
0329             ifrptl(1,nptl)=0
0330             ifrptl(2,nptl)=0
0331             xorptl(1,nptl)=0.
0332             xorptl(2,nptl)=0.
0333             xorptl(3,nptl)=0.
0334             xorptl(4,nptl)=0.
0335             tivptl(1,nptl)=0.
0336             tivptl(2,nptl)=0.
0337             idptl(nptl)=id
0338
0339 c boost in CMS frame
0340             call utlob5(yhaha, pptl(1,nptl), pptl(2,nptl)
0341      .        , pptl(3,nptl), pptl(4,nptl), pptl(5,nptl))
0342
0343             if(ish.ge.5)write(ifch,'(a,i5,a,i5,a,4(e10.4,1x),f6.3)')
0344      $       ' particle from qgsjet ',nptl,' id :',idptl(nptl)
0345      $  , ' momentum :',(pptl(k,nptl),k=1,5)
0346
0347
0348       enddo
0349
0350 1000  return
0351
0352 1001  iret=-1
0353       goto 1000
0354
0355       end
0356
0357 c-----------------------------------------------------------------------
0358       subroutine conqgsII
0359 c-----------------------------------------------------------------------
0360 c  determines interaction configuration
0361 c-----------------------------------------------------------------------
0362
0363       include 'epos.inc'
0364       include 'epos.incems'
0365       common/geom/rmproj,rmtarg,bmax,bkmx
0366       common/nucl3/phi,bimp
0367       common /cncl/xproj(mamx),yproj(mamx),zproj(mamx)
0368      *            ,xtarg(mamx),ytarg(mamx),ztarg(mamx)
0369       common/cncl2/diffbpr(mamx,mamx),diffbtg(mamx,mamx)
0370       common/cncl3/iactpr(mamx),iacttg(mamx)
0371       common/cfacmss/facmss
0372
0373       call utpri('cqgsII ',ish,ishini,4)
0374
0375       iret=0
0376
0377 c     initialisations
0378 c     ---------------
0379
0380       vel=tanh(ypjtl-yhaha)+tanh(yhaha)
0381
0382 c     determine phi, bimp, coll, iproj, itarg, x/y/zproj, x/y/ztarg
0383 c     ---------------------------------------------------------------
0384
0385
0386       if(maproj.eq.1.and.matarg.eq.1)then
0387
0388       koll=1
0389       do n=1,4
0390         coord(n,1)=0.
0391       enddo
0392       bk(1)=bimp
0393       iproj(1)=1
0394       itarg(1)=1
0395       phi=phievt
0396       xproj(1)=bimp*cos(phi)
0397       yproj(1)=bimp*sin(phi)
0398       zproj(1)=0
0399       xtarg(1)=0
0400       ytarg(1)=0
0401       ztarg(1)=0
0402       lproj(1)=1
0403       ltarg(1)=1
0404       lproj3(1)=1
0405       ltarg3(1)=1
0406       kproj(1,1)=1
0407       ktarg(1,1)=1
0408
0409            else
0410
0411       bx=0
0412       by=0
0413       if(maproj.gt.0)then
0414         phi=phievt
0415         bimp=bimevt
0416         bx=cos(phi)*bimp
0417         by=sin(phi)*bimp
0418       endif
0419       if(maproj.eq.0)goto 1000
0420       koll=0
0421       do i=1,maproj
0422       lproj(i)=0
0423       enddo
0424       do j=1,matarg
0425       ltarg(j)=0
0426       enddo
0427       do 12 i=1,maproj
0428       do 11 j=1,matarg
0429       bij=sqrt((xproj(i)+bx-xtarg(j))**2+(yproj(i)+by-ytarg(j))**2)
0430       if(ish.ge.7)write(ifch,*)'i_p:',i,' i_t:',j,' b_ij:',bij
0431       if(bij.gt.bkmx)goto 11
0432
0433       if(koll.ge.kollmx)goto 1000
0434       koll=koll+1
0435
0436       bk(koll)=bij
0437       bkx(koll)=xproj(i)+bx-xtarg(j)
0438       bky(koll)=yproj(i)+by-ytarg(j)
0439       iproj(koll)=i
0440       itarg(koll)=j
0441       lproj(i)=lproj(i)+1
0442       ltarg(j)=ltarg(j)+1
0443       kproj(i,lproj(i))=koll
0444       ktarg(j,ltarg(j))=koll
0445
0446
0447 11    continue
0448 12    continue
0449
0450
0451       do k=1,koll
0452         do n=1,4
0453           coord(n,k)=0.
0454         enddo
0455       enddo
0456
0457           endif
0458
0459       if(ish.ge.3)write(ifch,*)'koll=',koll
0460       if(koll.eq.0)goto 1000
0461
0462
0463 c     determine coord
0464 c     ---------------
0465       do kl=1,koll
0466       i=iproj(kl)
0467       j=itarg(kl)
0468       dist=ztarg(j)-zproj(i)
0469       coord(1,kl)=(xproj(i)+xtarg(j))*0.5
0470       coord(2,kl)=(yproj(i)+ytarg(j))*0.5
0471       coord(3,kl)=(zproj(i)+ztarg(j))*0.5
0472       coord(4,kl)=dist/vel
0473       enddo
0474
0475 c     determine number of collisions according glauber model
0476 c     ------------
0477       nglevt=koll
0478
0479 c     exit
0480 c     ----
0481 1000  continue
0482       if(ish.ge.5)then
0483       write(ifch,*)'ia,x/y/zproj:'
0484       do mm=1,maproj
0485       write(ifch,*)mm,iactpr(mm),xproj(mm),yproj(mm),zproj(mm)
0486       enddo
0487       write(ifch,*)'ia,x/y/ztarg:'
0488       do mm=1,matarg
0489       write(ifch,*)mm,iacttg(mm),xtarg(mm),ytarg(mm),ztarg(mm)
0490       enddo
0491       write(ifch,*)'iret',iret
0492       endif
0493       call utprix('cqgsII ',ish,ishini,4)
0494       return
0495
0496       end
0497
0498
0499 c------------------------------------------------------------------------------
0500       function qgsIIcrse(egy,mapro,matar,id)
0501 c------------------------------------------------------------------------------
0502 c inelastic cross section of qgsjet-II
0503 c (id=0 corresponds to air)
0504 c egy - kinetic energy
0505 c maproj - projec mass number     (1<maproj<64)
0506 c matarg - target mass number     (1<matarg<64)
0507 c------------------------------------------------------------------------------
0508       include 'epos.inc'
0509       double precision egyl,qgsect
0510
0511       qgsIIcrse=0.
0512       call iclass(idproj,icpro)
0513       call idmass(1120,amt1)
0514       call idmass(1220,amt2)
0515       amtar=0.5*(amt1+amt2)
0516       if(matar.eq.1)amtar=amt1
0517       if(mapro.eq.1)then
0518         call idmass(idproj,ampro)
0519       else
0520         ampro=mapro*amtar
0521       endif
0522       egyl=dble(egy+ampro)
0523
0524       if(id.eq.0)then
0525         do k=1,3
0526           mt=int(airanxs(k))
0527           qgsIIcrse=qgsIIcrse+airwnxs(k)*qgsect(egyl,icpro,mapro,mt)
0528         enddo
0529       else
0530         qgsIIcrse=qgsect(egyl,icpro,mapro,matar)
0531       endif
0532
0533       return
0534       end
0535
0536 c--------------------------------------------------------------------
0537       double precision function qgran(b10)
0538 c--------------------------------------------------------------------
0539 c Random number generator
0540 c--------------------------------------------------------------------
0541       include 'epos.inc'
0542       double precision b10,drangen
0543       qgran=drangen(b10)
0544       if(irandm.eq.1)write(ifch,*)'qgran()= ',qgran
0545
0546       return
0547       end
0548