ReggeGribovPartonMCInterface/src/qgsjet_epos.f

0001 c 25.04.2003 Link routines between QGSJet and epos.
0002 c author T. Pierog
0003
0004
0005 c-----------------------------------------------------------------------
0006       subroutine IniQGSJet
0007 c-----------------------------------------------------------------------
0008 c Primary initialization for QGSJet
0009 c-----------------------------------------------------------------------
0010       include 'epos.inc'
0011       COMMON /Q_DEBUG/  DEBUG
0012       COMMON /Q_AREA43/ MONIOU
0013       integer debug
0014       double precision BQGS,BMAXQGS,BMAXNEX,BMINNEX,XA(210,3),XB(210,3)
0015       COMMON /Q_QGSNEX1/ XA,XB,BQGS,BMAXQGS,BMAXNEX,BMINNEX  !ctp
0016
0017       call utpri('iniqgs',ish,ishini,6)
0018       write(ifmt,'(a,i6)')'initialize QGSJet ...'
0019
0020       debug=0
0021       if(ish.ge.3)debug=ish-2
0022       moniou=ifch
0023
0024 c common model parameters setting
0025       call psaset
0026 c particular model parameters setting
0027       call xxaset
0028 c common initialization procedure
0029       call qgspsaini
0030       BMAXNEX=dble(bmaxim)
0031       BMINNEX=dble(bminim)
0032       egymin=0.1
0033       egymax=egymax
0034       irescl=0
0035
0036
0037       call utprix('iniqgs',ish,ishini,6)
0038       end
0039
0040 c-----------------------------------------------------------------------
0041       subroutine IniEvtQGS
0042 c-----------------------------------------------------------------------
0043 c Initialization for each type of event (for given proj, targ and egy)
0044 c-----------------------------------------------------------------------
0045       include 'epos.inc'
0046       common/geom/rmproj,rmtarg,bmax,bkmx
0047 c QGSJet Common
0048       double precision XA(210,3),XB(210,3),BQGS,BMAXQGS,BMAXNEX,BMINNEX
0049      &  ,e0
0050       COMMON /Q_QGSNEX1/ XA,XB,BQGS,BMAXQGS,BMAXNEX,BMINNEX
0051
0052       if(matarg.gt.210.or.maproj.gt.210)
0053      &  call utstop('Nucleus too big for QGSJet (Mmax=210) !&')
0054       call iclass(idproj,iclpro)
0055       call iclass(idtarg,icltar)
0056       e0=dble(elab)
0057       icp=idtrafo('nxs','qgs',idproj)
0058       if(icp.eq.0)icp=1-2*int(rangen()+0.5)       !pi0=pi+ or p-
0059       call xxaini(e0,icp,maproj,matarg)
0060       bmax=BMAXQGS
0061       qgsincs=fqgscrse(ekin,maproj,matarg)
0062       if(engy.lt.egymin)qgsincs=0.          !below egymin, no interaction
0063       call xsigma
0064       bkmx=sqrt(sigine/10./pi)        !10= fm^2 -> mb
0065       if(ish.ge.2)write(ifch,*)
0066      &  'QGSJet used with (E,proj,maproj,matarg,bmax)',e0,icp,maproj
0067      &  ,matarg,bmax
0068
0069       return
0070       end
0071
0072 c-----------------------------------------------------------------------
0073       subroutine emsqgs(iret)
0074 c-----------------------------------------------------------------------
0075 c  call qgsjet to simulate interaction
0076 c-----------------------------------------------------------------------
0077       include 'epos.inc'
0078       include 'epos.incems'
0079       common/geom/rmproj,rmtarg,bmax,bkmx
0080       double precision XA(210,3),XB(210,3),BQGS,BMAXQGS,BMAXNEX,BMINNEX
0081       COMMON /Q_QGSNEX1/ XA,XB,BQGS,BMAXQGS,BMAXNEX,BMINNEX
0082       common/nucl3/phi,bimp
0083       common /q_area12/ nsp
0084       common /q_area14/ esp(4,95000),ich(95000)
0085       double precision esp
0086 c NSF - number of secondary fragments;
0087 c IAF(i) - mass of the i-th fragment
0088       COMMON /Q_AREA13/ NSF,IAF(210)
0089       COMMON /Q_AREA99/ NWT
0090
0091       iret=0
0092       b1=bminim
0093       b2=min(bmax,bmaxim)
0094       a=pi*(b2**2-b1**2)
0095
0096       if(a.gt.0..and.rangen().gt.qgsincs/10./a)goto 1001   !no interaction
0097       if(ish.ge.3)call alist('Determine QGSJet Production&',0,0)
0098
0099       nptl=0
0100       nsp=0
0101       call psconf
0102
0103       nevt=1
0104       kolevt=-1
0105       koievt=-1
0106       kohevt=-1
0107       npjevt=maproj
0108       ntgevt=matarg
0109       pmxevt=pnll
0110       egyevt=engy
0111       if(BQGS.ge.0.d0)then
0112         bimevt=real(BQGS)
0113         bimp=real(BQGS)
0114         phievt=2.*pi*rangen()
0115       else
0116         bimevt=0.
0117         bimp=0.
0118         phievt=0.
0119       endif
0120       anintine=anintine+1.
0121
0122       call conre
0123       call conwr
0124
0125 c keep the projectile spectators as fragments
0126       npns=0
0127       npps=0
0128       ns=0
0129       if(infragm.eq.2)then
0130
0131         if(NSF.gt.0)then
0132           do is=1,NSF           !count the number of spectators
0133           if(ish.ge.7)write(ifch,'(a,i5,a,i5)')
0134      $       ' Projecticle Fragment ',is,' Mass :',IAF(is)
0135             nptl=nptl+1
0136             istptl(nptl)=0
0137             if(IAF(is).eq.1)then
0138               id=idptl(is)
0139               pptl(3,nptl)=pptl(3,is)
0140               pptl(4,nptl)=pptl(4,is)
0141               pptl(5,nptl)=pptl(5,is)
0142               if(id.eq.1120)npps=npps+1
0143               if(id.eq.1220)npns=npns+1
0144             else
0145               if(IAF(is).eq.2)then
0146                 id=17
0147                 npps=npps+1
0148                 npns=npns+1
0149               elseif(IAF(is).eq.3)then
0150                 id=18
0151                 npps=npps+1
0152                 npns=npns+2
0153               elseif(IAF(is).eq.4)then
0154                 id=19
0155                 npps=npps+2
0156                 npns=npns+2
0157               else
0158                 inucl=IAF(is)
0159                 iprot= int(dble(inucl) / 2.15d0 + 0.7d0)
0160                 id=1000000000+iprot*10000+inucl*10 !code for nuclei
0161                 npps=npps+iprot
0162                 npns=npns+inucl-iprot
0163               endif
0164               call idmass(id,am)
0165               pptl(4,nptl)=dble(IAF(is))*pptl(4,is)      !Etot
0166               pptl(5,nptl)=am                             !mass
0167               pz2tmp=(pptl(4,nptl)+am)*(pptl(4,nptl)-am)
0168               if(pz2tmp.gt.0.d0)then
0169                 pptl(3,nptl)=sqrt(pz2tmp)                 !Pz
0170               else
0171                 write(*,*)'Warning in emsqgs !'
0172                 write(*,*)'energy of fragment too small :',IAF(is),am
0173      &                     ,pptl(4,nptl)
0174                 pptl(3,nptl)=pptl(4,nptl)
0175               endif
0176             endif
0177             pptl(1,nptl)=0.d0 !P_x
0178             pptl(2,nptl)=0.d0 !P_y
0179             ityptl(nptl)=0
0180             iorptl(nptl)=1
0181             jorptl(nptl)=maproj+matarg
0182             ifrptl(1,nptl)=0
0183             ifrptl(2,nptl)=0
0184             xorptl(1,nptl)=0.d0
0185             xorptl(2,nptl)=0.d0
0186             xorptl(3,nptl)=0.d0
0187             xorptl(4,nptl)=0.d0
0188             tivptl(1,nptl)=0.d0
0189             tivptl(2,nptl)=0.d0
0190             idptl(nptl)=id
0191             if(ish.ge.5)write(ifch,'(a,i5,a,i5,a,4(e10.4,1x),f6.3)')
0192      $       ' Fragment from qgsjet ',nptl,' id :',idptl(nptl)
0193      $  , ' momentum :',(pptl(k,nptl),k=1,5)
0194
0195           enddo
0196         endif
0197
0198 c make the projectile spectators as free nucleons
0199
0200       else
0201         ns=0
0202         if(NSF.gt.0)then
0203           do is=1,NSF           !count the number of spectators
0204             ns=ns+IAF(is)
0205           enddo
0206           if(infragm.eq.1)then
0207 c  remaining nucleus is one fragment
0208             nptl=nptl+1
0209             istptl(nptl)=0
0210             pptl(1,nptl)=0.d0
0211             pptl(2,nptl)=0.d0
0212             pptl(4,nptl)=0.d0
0213             inucl=0
0214             do is=1,ns
0215               inucl=inucl+1
0216               pptl(4,nptl)=pptl(4,nptl)+pptl(4,is)
0217             enddo
0218             iprot= int(dble(inucl) / 2.15d0 + 0.7d0)
0219             idnucl=1000000000+iprot*10000+inucl*10 !code for nuclei
0220             npps=npps+iprot
0221             npns=npns+inucl-iprot
0222             call idmass(idnucl,am)
0223             pptl(5,nptl)=am  !mass
0224             ptot=(pptl(4,nptl)+am)*(pptl(4,nptl)-am)
0225             pptl(3,nptl)=sqrt(ptot)
0226             ityptl(nptl)=0
0227             istptl(nptl)=0
0228             iorptl(nptl)=1
0229             jorptl(nptl)=maproj
0230             ifrptl(1,nptl)=0
0231             ifrptl(2,nptl)=0
0232             xorptl(1,nptl)=xorptl(1,1)
0233             xorptl(2,nptl)=xorptl(2,1)
0234             xorptl(3,nptl)=xorptl(3,1)
0235             xorptl(4,nptl)=xorptl(4,1)
0236             tivptl(1,nptl)=tivptl(1,1)
0237             tivptl(2,nptl)=tivptl(2,1)
0238             idptl(nptl)=idnucl
0239           else
0240             do is=maproj-ns+1,maproj          !make the ns last projectile nucleon final
0241               if(idptl(is).eq.1120)npps=npps+1
0242               if(idptl(is).eq.1220)npns=npns+1
0243             enddo
0244           endif
0245         endif
0246       endif
0247
0248 c number of participants
0249       if(laproj.gt.1)then
0250         npjevt=maproj-npps-npns
0251         npppar=max(0,laproj-npps)
0252         npnpar=npjevt-npppar
0253 c set participant projectile as non spectators
0254         do i=1,maproj
0255           if(idptl(i).eq.1120)then
0256             if(npppar.gt.0)then
0257               npppar=npppar-1
0258             else                !restore spectators
0259               iorptl(i)=0
0260               if(infragm.eq.0)istptl(i)=0
0261             endif
0262           endif
0263           if(idptl(i).eq.1220)then
0264             if(npnpar.gt.0)then
0265               npnpar=npnpar-1
0266             else                !restore spectators
0267               iorptl(i)=0
0268               if(infragm.eq.0)istptl(i)=0
0269             endif
0270           endif
0271         enddo
0272       endif
0273
0274 c restore target spectators
0275       ns=0
0276       if(NWT.lt.matarg)then
0277 c number of participants
0278         ntgevt=NWT
0279         do is=ntgevt+1,matarg         !make the last target nucleon final
0280           iorptl(maproj+is)=0
0281           istptl(maproj+is)=0
0282         enddo
0283       endif
0284
0285
0286       do is=1,nsp
0287
0288 c ich is the type of secondary hadron, esp - its transverse momentum,
0289 c and its energy
0290 c the following notations for the particles types are used: 0 - pi0, 1 -
0291 c pi+,
0292 c -1 - pi-, 2 - p, -2 - p, 3 - n, -3 - n, 4 - k+, -4 - k-, 5 - k0s, -5 -
0293 c k0l
0294           ic=ich(is)
0295           if(ish.ge.7)write(ifch,'(a,i5,a,i5,2a,4(e10.4,1x))')
0296      $       ' qgsjet particle ',is,' id :',ic,' before conversion'
0297      $     , ' momentum :',(esp(k,is),k=1,4)
0298
0299             nptl=nptl+1
0300             if(nptl.gt.mxptl)call utstop('qgsjet: mxptl too small&')
0301             id=idtrafo('qgs','nxs',ic)
0302             if(ish.ge.7)write(ifch,'(a,i5,a,i5,a)')
0303      $       ' epos particle ',nptl,' id :',id,' after conversion'
0304             call idmass(id,am)
0305
0306
0307             pptl(1,nptl)=real(esp(3,is)) !P_x
0308             pptl(2,nptl)=real(esp(4,is)) !P_y
0309             pptl(3,nptl)=real(esp(2,is)) !P_z
0310             pptl(4,nptl)=real(esp(1,is)) !E
0311             pptl(5,nptl)=am              !mass
0312             istptl(nptl)=0
0313             ityptl(nptl)=0
0314             iorptl(nptl)=1
0315             jorptl(nptl)=maproj+matarg
0316             ifrptl(1,nptl)=0
0317             ifrptl(2,nptl)=0
0318             xorptl(1,nptl)=0.
0319             xorptl(2,nptl)=0.
0320             xorptl(3,nptl)=0.
0321             xorptl(4,nptl)=0.
0322             tivptl(1,nptl)=0.
0323             tivptl(2,nptl)=0.
0324             idptl(nptl)=id
0325
0326             if(noebin.lt.0)pptl(3,nptl)=-pptl(3,nptl) !exchange projectile <-> target side in case of fake DIS
0327
0328             if(ish.ge.5)write(ifch,'(a,i5,a,i5,a,4(e10.4,1x),f6.3)')
0329      $       ' particle from qgsjet ',nptl,' id :',idptl(nptl)
0330      $  , ' momentum :',(pptl(k,nptl),k=1,5)
0331
0332
0333       enddo
0334
0335 1000  return
0336
0337 1001  iret=-1
0338       goto 1000
0339
0340       end
0341
0342
0343 c------------------------------------------------------------------------------
0344       function fqgscrse(ek,mapr,matg)
0345 c------------------------------------------------------------------------------
0346 c hadron-nucleus (hadron-proton) and nucl-nucl particle production cross section
0347 c with qgsjet.
0348 c ek - kinetic lab energy
0349 c maproj - projec mass number     (1<maproj<210)
0350 c matarg - target mass number     (1<matarg<210)
0351 c------------------------------------------------------------------------------
0352       dimension wk(3),wa(3),wb(3)
0353       include 'epos.inc'
0354       double precision gsect,qgsasect
0355       COMMON /Q_XSECT/  GSECT(10,5,4)
0356       COMMON /Q_AREA48/ QGSASECT(10,6,4)
0357
0358       fqgscrse=0.
0359       call idmass(1120,amt1)
0360       call idmass(1220,amt2)
0361       amtar=0.5*(amt1+amt2)
0362       if(matg.eq.1)amtar=amt1
0363       if(mapr.eq.1)then
0364         call idmass(idproj,ampro)
0365       else
0366         ampro=mapr*amtar
0367       endif
0368       egy=ek+ampro
0369       ye=max(1.,log10(egy))
0370       je=min(8,int(ye))
0371
0372       wk(2)=ye-je
0373       wk(3)=wk(2)*(wk(2)-1.)*.5
0374       wk(1)=1.-wk(2)+wk(3)
0375       wk(2)=wk(2)-2.*wk(3)
0376
0377       ya=real(matg)
0378       ya=log(ya)/1.38629+1.
0379       ja=min(int(ya),2)
0380       wa(2)=ya-ja
0381       wa(3)=wa(2)*(wa(2)-1.)*.5
0382       wa(1)=1.-wa(2)+wa(3)
0383       wa(2)=wa(2)-2.*wa(3)
0384
0385       if(mapr.eq.1)then
0386
0387         do i=1,3
0388           do m=1,3
0389          fqgscrse=fqgscrse+real(gsect(je+i-1,iclpro,ja+m-1))*wk(i)*wa(m)
0390           enddo
0391         enddo
0392
0393       else
0394
0395         yb=mapr
0396         yb=log(yb/2.)/.69315+1.
0397         jb=min(int(yb),4)
0398         wb(2)=yb-jb
0399         wb(3)=wb(2)*(wb(2)-1.)*.5
0400         wb(1)=1.-wb(2)+wb(3)
0401         wb(2)=wb(2)-2.*wb(3)
0402
0403         do i=1,3
0404           do m=1,3
0405             do n=1,3
0406               fqgscrse=fqgscrse+real(qgsasect(je+i-1,jb+n-1,ja+m-1)
0407      &                     *wk(i)*wa(m)*wb(n))
0408             enddo
0409           enddo
0410         enddo
0411
0412       endif
0413
0414       fqgscrse=exp(fqgscrse)
0415       return
0416       end
0417
0418 c------------------------------------------------------------------------------
0419       function qgscrse(egy,mapro,matar,id)
0420 c------------------------------------------------------------------------------
0421 c inelastic cross section of qgsjet
0422 c (id=0 corresponds to air)
0423 c egy - kintetic energy
0424 c maproj - projec mass number     (1<maproj<210)
0425 c matarg - target mass number     (1<matarg<210)
0426 c------------------------------------------------------------------------------
0427       include 'epos.inc'
0428
0429       qgscrse=0.
0430       if(id.eq.0)then
0431         do k=1,3
0432           mt=int(airanxs(k))
0433           qgscrse=qgscrse+airwnxs(k)*fqgscrse(egy,mapro,mt)
0434         enddo
0435       else
0436         qgscrse=fqgscrse(egy,mapro,matar)
0437       endif
0438
0439       return
0440       end
0441
0442 c--------------------------------------------------------------------
0443       double precision function psran(b10)
0444
0445
0446 c--------------------------------------------------------------------
0447 c  Random number generator
0448 c--------------------------------------------------------------------
0449       double precision b10,drangen
0450       include 'epos.inc'
0451       psran=drangen(b10)
0452       if(irandm.eq.1)write(ifch,*)'psran()= ',psran
0453
0454       return
0455       end
0456
0457
0458