*********************************************************************** * subroutine filegen(nvtx,xvert,beta) * * generates kinematics by reading from input files * created on: 18/01/2000 by R. Schicker * * modified on 05/04/2005 by R. Holzmann * last modified on 10/12/2007 by R. Holzmann * * - optionally event plane is rotated randomly for each event (and each file!) * - beta calculated from beampar(3) is overridden *********************************************************************** implicit none #include "geant321/gconst.inc" #include "geant321/gcflag.inc" #include "geant321/gckine.inc" #include "geant321/gclist.inc" #include "geant321/gctrak.inc" #include "kinetups.inc" #include "user.inc" #include "geaevent.inc" #include "geatarget.inc" *----------------------------------------------------- * Read kinematics from unit(s) 24+ (ascii file(s)) * loop over all event input files *----------------------------------------------------- integer nvtx real xvert(3), beta(4), tof, tofsave real Mass, GeantMass, du3, du4, du5, du6 character*20 cdu1 integer idu2 real cphir, sphir, phir, px, py, etot, random real ubuf(5), p(3), pab, weight, dv(3), xv(3) integer i, ifile, itr, itrack, itype, npart, idsrc, idpar,ixpar integer icharge, iweight, iswitch real uval real genetype1, genetype2, genetype3 logical has_weight do ifile=1,NEvFiles if (iswit(6).gt.0) then ! sample event plane angle call grndm(random,1) phir = twopi*random ! twopi taken from gconst.inc cphir = cos(phir) sphir = sin(phir) else phir = 0. cphir = 1. sphir = 0. end if if (ifile.eq.1) phi_event = 57.2958*phir ! save event-plane angle 15 eventID=0 npart=0 read (23+ifile,*,end=900) eventID,npart,ebeam,bpar,iweight iswitch = iweight/10 iweight = iweight - iswitch*10 iswitch = iabs(iswitch) if (eventId.eq.0 .and. npart.eq.0) goto 900 ! end of file ebeam = ebeam*1000.0 ! convert to MeV (e.g. Pluto has GeV!) 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) *--------------------------------------------------------- * Check for weighted particles * *--------------------------------------------------------- if (iweight.eq.0) then has_weight = .false. else has_weight = .true. endif c genetype1 = float(ifile) ! bug genetype1 = 0. genetype2 = 0. genetype3 = -1. *----------------------------------------------------------------------- * Loop over particles in file ifile *----------------------------------------------------------------------- xv(1) = xvert(1) xv(2) = xvert(2) xv(3) = xvert(3) do itr=1,npart uval = -1. if (has_weight) then if(iweight.eq.1) then if(iswitch.eq.0) then read(23+ifile,*) Etot,p,itype,weight else read(23+ifile,*) Etot,p,itype,weight,uval endif endif if(iweight.eq.2) then if(iswitch.eq.0) then read(23+ifile,*) Etot,p,itype,idsrc,idpar,weight else read(23+ifile,*) Etot,p,itype,idsrc,idpar,weight, + uval endif genetype1 = float(idsrc) genetype2 = float(idpar) endif if(iweight.eq.-2) then if(iswitch.eq.0) then read(23+ifile,*) Etot,p,itype,idsrc,idpar,ixpar, + weight else read(23+ifile,*) Etot,p,itype,idsrc,idpar,ixpar, + weight,uval endif genetype1 = float(idsrc) ! source id genetype2 = float(idpar) ! parent id genetype3 = float(ixpar) ! parent index endif if(iweight.eq.3) then ! read with vertex if(iswitch.eq.0) then read(23+ifile,*) Etot,p,dv,itype,idsrc,idpar,weight else read(23+ifile,*) Etot,p,dv,itype,idsrc,idpar,weight, + uval endif genetype1 = float(idsrc) genetype2 = float(idpar) xv(1) = xvert(1) + 0.1*dv(1) ! adjust vertex (dv is in mm) xv(2) = xvert(2) + 0.1*dv(2) xv(3) = xvert(3) + 0.1*dv(3) call gsvert(xv,0,0,ubuf,0,nvtx) end if if(iweight.eq.-3) then ! read with vertex if(iswitch.eq.0) then read(23+ifile,*) + Etot,p,dv,itype,idsrc,idpar,ixpar,weight else read(23+ifile,*) + Etot,p,dv,itype,idsrc,idpar,ixpar,weight,uval endif genetype1 = float(idsrc) genetype2 = float(idpar) genetype3 = float(ixpar) xv(1) = xvert(1) + 0.1*dv(1) ! adjust vertex (dv is in mm) xv(2) = xvert(2) + 0.1*dv(2) xv(3) = xvert(3) + 0.1*dv(3) call gsvert(xv,0,0,ubuf,0,nvtx) end if if(iweight.eq.4) then ! read with vertex and tof if(iswitch.eq.0) then read(23+ifile,*) + Etot,p,tof,dv,itype,idsrc,idpar,weight else read(23+ifile,*) + Etot,p,tof,dv,itype,idsrc,idpar,weight,userval endif genetype1 = float(idsrc) genetype2 = float(idpar) xv(1) = xvert(1) + 0.1*dv(1) ! adjust vertex (dv is in mm) xv(2) = xvert(2) + 0.1*dv(2) xv(3) = xvert(3) + 0.1*dv(3) tofsave = tofg tofg = tof*1.e-9 ! set correct tof in /GCTRAK/ call gsvert(xv,0,0,ubuf,0,nvtx) tofg = tofsave end if if(iweight.eq.-4) then ! read with vertex and tof if(iswitch.eq.0) then read(23+ifile,*) + Etot,p,tof,dv,itype,idsrc,idpar,ixpar,weight else read(23+ifile,*) + Etot,p,tof,dv,itype,idsrc,idpar,ixpar,weight, + uval endif genetype1 = float(idsrc) genetype2 = float(idpar) genetype3 = float(ixpar) xv(1) = xvert(1) + 0.1*dv(1) ! adjust vertex (dv is in mm) xv(2) = xvert(2) + 0.1*dv(2) xv(3) = xvert(3) + 0.1*dv(3) tofsave = tofg tofg = tof*1.e-9 ! set correct tof in /GCTRAK/ call gsvert(xv,0,0,ubuf,0,nvtx) tofg = tofsave end if call gfpart(itype,cdu1,idu2,GeantMass,du3,du4,du5,du6) !get mass Mass = sqrt(Etot*Etot-p(1)*p(1)-p(2)*p(2)-p(3)*p(3)) if (abs(GeantMass-Mass).gt.0.005) then c particle is off-shell (this is an UrQMD bug!), flag it. genetype1 = -genetype1 genetype2 = -genetype2 if (idsrc.eq.0 .and. idpar.eq.0) genetype2 = -1. endif else if(iswitch.eq.0) then read(23+ifile,*) Etot, p, itype else read(23+ifile,*) Etot, p, itype, uval endif weight = 1.0 endif c if (iswit(5).gt.0 .and. idevt.lt.iswit(5)) goto 20 c ! start with event iswit(5) if (iswit(1).gt.0 .and. itype.ne.iswit(1)) goto 20 ! skip this particle, do not track it if (iswit(1).lt.0 .and. itype.eq.-iswit(1)) goto 20 if (Etot.lt.0.0) then Etot = -Etot icharge = -1 ! negative particle else icharge = 1 ! positive particle endif px = p(1) py = p(2) pab = sqrt(px*px + py*py + p(3)*p(3)) p(1) = px*cphir - py*sphir ! rotate event in x-y plane to make p(2) = px*sphir + py*cphir ! sure that impact vector changes ubuf(1) = genetype1 ubuf(2) = genetype2 ubuf(3) = genetype3 ubuf(4) = weight ubuf(5) = uval if (iswit(3).eq.1) then ! input is in CM, transform to lab first call cmtolab(beta,p,itype) end if call gskine(p,itype,nvtx,ubuf,5,itrack) if (itrack.gt.0) then ! save primary kinematics call copykine(itrack,itype,p,xv,0,0, + weight,uval,genetype1,genetype2,genetype3) endif 20 continue; enddo ! end loop over particles in event 900 if (eventId.eq.0 .and. npart.eq.0) then c rewind(unit=23+ifile) c write(6,'('' File `'',a,''` rewound'')') EvFileName(ifile) c goto 15 write(6,'('' End of File `'',a,''` reached!'')') EvFileName(ifile) NEVENT=0 endif 99 continue enddo ! end loop over event files return end