*********************************************************************** * 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" #include "eventbuffer.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, ivertex real uval real genetype1, genetype2, genetype3 logical has_weight *----------------------------------------------------- * LOOP Repeat the pevious event nmaxloop1 times *----------------------------------------------------- if (nmaxloop1.gt.0) then eventnr1 = eventnr1 + 1 ! eventid from file will be off if (stop1.gt.0) then nloop1 = 0 ! repeat loop successful : read next stop1 = 0 ! reset the stop flag endif if (nloop1.gt.0) then ! do the loop for event repeat eventid = eventid1 npart = npart1 ebeam = ebeam1 bpar = bpar1 phi_event = phi_event1 do itr=1,npart1 xv(1) = vtx1(itr,1) xv(2) = vtx1(itr,2) xv(3) = vtx1(itr,3) p(1) = pmom1(itr,1) p(2) = pmom1(itr,2) p(3) = pmom1(itr,3) tof = tof1(itr) genetype1 = ubuf1(itr,1) genetype2 = ubuf1(itr,2) genetype3 = ubuf1(itr,3) weight = ubuf1(itr,4) uval = ubuf1(itr,5) ubuf (1) = ubuf1(itr,1) ubuf (2) = ubuf1(itr,2) ubuf (3) = ubuf1(itr,3) ubuf (4) = ubuf1(itr,4) ubuf (5) = ubuf1(itr,5) tofsave = tofg tofg = tof*1.e-9 ! set correct tof in /GCTRAK/ call gsvert(xv,0,0,ubuf,0,nvtx) tofg = tofsave call gskine(p,itype1(itr),nvtx,ubuf,5,itrack) if (itrack.gt.0) then ! save primary kinematics call copykine(itrack,itype1(itr),p,xv, + 0,0,weight,uval,genetype1,genetype2,genetype3) endif enddo nloop1 = nloop1 - 1 ! at nloop = 0 new event will be read if ( eventnr1.eq.NEVENT ) then ! Finish event reading NEVENT=0 endif goto 902 ! skip to return end if endif ! maxloop1 > 0 *----------------------------------------------------- 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 ivertex = iweight/100 iswitch = (iweight-ivertex*100)/10 iweight = iweight -ivertex*100 - iswitch*10 iswitch = iabs(iswitch) ivertex = iabs(ivertex) 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) *--------------------------------------------------------- * LOOP : make sure vertex is taken from file * if(nmaxloop1.gt.0) ivertex = 1 *--------------------------------------------------------- *--------------------------------------------------------- * 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) if(ivertex.eq.0) then 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) else xv(1) = 0.1*dv(1) ! adjust vertex (dv is in mm) xv(2) = 0.1*dv(2) xv(3) = 0.1*dv(3) endif 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) if(ivertex.eq.0) then 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) else xv(1) = 0.1*dv(1) ! adjust vertex (dv is in mm) xv(2) = 0.1*dv(2) xv(3) = 0.1*dv(3) endif 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,uval endif genetype1 = float(idsrc) genetype2 = float(idpar) if(ivertex.eq.0) then 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) else xv(1) = 0.1*dv(1) ! adjust vertex (dv is in mm) xv(2) = 0.1*dv(2) xv(3) = 0.1*dv(3) endif 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) if(ivertex.eq.0) then 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) else xv(1) = 0.1*dv(1) ! adjust vertex (dv is in mm) xv(2) = 0.1*dv(2) ! xv(3) = 0.1*dv(3) endif 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 *----------------------------------------------------- * LOOP Make a copy of current event *----------------------------------------------------- if (nmaxloop1.gt.0) then ! make a copy of event to repeat it nmaxloop1 eventid1 = eventid npart1 = npart ebeam1 = ebeam bpar1 = bpar phi_event1 = phi_event pmom1(itr,1) = p(1) pmom1(itr,2) = p(2) pmom1(itr,3) = p(3) vtx1(itr,1) = xv(1) vtx1(itr,2) = xv(2) vtx1(itr,3) = xv(3) tof1(itr) = tof itype1(itr) = itype ubuf1(itr,1) = ubuf(1) ubuf1(itr,2) = ubuf(2) ubuf1(itr,3) = ubuf(3) ubuf1(itr,4) = ubuf(4) ubuf1(itr,5) = ubuf(5) nloop1 = nmaxloop1 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 902 return end