c***************************************************************** subroutine thermalgen(ind,nvtx,xvert) c thermal source of particles in CM system defined by parameters c pmc2(1) ... pmc2(10) c c npart number of particles/event c ipart Geant particle code c T1 temperature of low-momentum component in MeV c T2 temperature of high-momentum component in MeV c frac fraction such that I = I1*frac + I2*(1-frac) (frac = 0-1) c beta blast velocity (v/c) c A2 polar distribution (cos**2 coefficient) c A4 polar distribution (cos**4 coefficient) c v1 directed flow c v2 elliptic flow (= squeeze-out) c c last modified on: 21/08/2000 by R.Holzmann c***************************************************************** implicit none #include "geant321/gclist.inc" #include "geant321/gckine.inc" #include "geant321/gcflag.inc" #include "hcons.inc" #include "kinetups.inc" #include "pmc.inc" integer ind,nvtx real xvert(3) real T1, T2, frac, blast, A2, A4, v1, v2, vmax real p(3), buf(4), beta(4) real theta, phi, Etot, ptot, restm, sint, cost, sinp, cosp integer npar, jpar, itrk, ip, i character*20 cdum1 real dum2, dum3, dum4, dum5, dum6 real rndm, hrndm1 real Ap, At, Apart, bmax, prob phi_event = 2.*pi*rndm(0) ! event plane angle if(pmc0(1).ge.0.0) then npar = ifix(pmc0(1)+0.01) ! fixed number of particles else call rnpssn(-pmc0(1),npar,i) ! Poisson with = |pcm0(1)| endif if(pmc0(1).eq.0.) then ! sample impact parameter, i.e. Apart Ap = pmc0(11) At = pmc0(12) prob = pmc0(13) bmax = 1.14*(Ap**0.333+At**0.333) bpar = bmax*sqrt(rndm(1)) ! sample impact parameter call rnpssn(prob*Apart(Ap,At,bpar),npar,i) ! sample nb of particles else bpar = 0.0 end if jpar = ifix(pmc0(2)+0.01) ! Geant particle code call gfpart(jpar,cdum1,dum2,restm,dum3,dum4,dum5,dum6) ! get mass (GeV) T1 = pmc0(3) T2 = pmc0(4) frac = pmc0(5) if(frac.gt.1.0) frac = 1.0 if(frac.lt.0.0) frac = 0.0 blast = pmc0(6) A2 = pmc0(7) A4 = pmc0(8) v1 = pmc0(9) v2 = pmc0(10) buf(1) = 400.0 ! generator type buf(2) = 0.0 buf(3) = -1. buf(4) = 1.0 ! weight beta(1) = 0.0 beta(2) = 0.0 beta(3) = -sqrt(Ebeam/(Ebeam+2*931.5)) ! beta(lab vs. cm) beta(4) = 1./sqrt(1.-beta(1)**2-beta(2)**2-beta(3)**2) ! gamma(cm) do ip=1,npar ! loop over particles now Etot = hrndm1(9900+ind) ! sample energy (MeV) theta = dtor*hrndm1(9910+ind) ! sample polar angle (degrees) sint = sin(theta) cost = cos(theta) vmax = 1. + cost*2.*v1 + sint*2.*v2 i = 0 1 phi = 2.*pi*rndm(i+100) !sample azimuthal angle (modulated w/ polar) i = i + 1 if(i.gt.100 .or. rndm(i)*vmax .gt. + 1.0 + cost*2.*v1*cos(phi) + sint*2.*v2*cos(2.*phi)) goto 1 sinp = sin(phi+phi_event) cosp = cos(phi+phi_event) Etot = 0.001*Etot ! go to GeV ptot = sqrt(Etot**2 - restm**2) ! total momentum (GeV/c) p(1) = ptot*sint*cosp p(2) = ptot*sint*sinp p(3) = ptot*cost call cmtolab(beta,p,jpar) ! transform to lab frame call gskine(p,jpar,nvtx,buf,3,itrk) ! pass particle to Geant if(itrk>0) + call copykine(itrk,jpar,p,xvert,0,0,1.,-1., + buf(1),buf(2),buf(3)) end do phi_event = rtod*phi_event ! event azimuth return end