************************************************************************ subroutine refev implicit none cSEQ,GCONST,DRPARS,GUREF. #include "geant321/gconsp.inc" #include "gurefp.inc" #include "guref.inc" #include "gufil.inc" c c This generator is used for reference trajectories c 7-jan-1993 R.Schicker c update to z,rho parametrization of vertex c nref(1)=1 grid traj C nref(1)=3 evaluate find matrix for omega tracks C nref(1)=4 evaluate find matrix for conv tracks c c nref(2)=1 : vertex in x,y (continuous in x,y) c nref(2)=2 : vertex in z,rho (continuous in z,rho) c nref(2)=3 : vertex delz,rho (cont in rho, z on targ segm) c c 13-apr-1993 R.Schicker c HITS INPUTS c c PMCR(30) kine parameter c------------------------------------------------------------ c structure for n particles on the grid c PMCR(1) generator id c PMCR(2) part. type c c NREF(2)=1: c PMCR(3) magnitude of minimum vertex x-position c PMCR(4) step size in x-position c PMCR(5) number of gridpoints in x-position c PMCR(6) magnitude of minimum vertex y-position c PMCR(7) step size in y-position c PMCR(8) number of gridpoints in y-position c c NREF(2)=2: c PMCR(3) magnitude of minimum vertex z-position c PMCR(4) step size in z-position c PMCR(5) number of gridpoints in z-position c PMCR(6) magnitude of minimum vertex rho-position c PMCR(7) step size in rho-position c PMCR(8) number of gridpoints in rho-position c c PMCR(9) magnitude of minimum inverse momentum c PMCR(10) step size in inverse momentum c PMCR(11) number of gridpoints in inverse momentum c c PMCR(12) magnitude of minimum polar angle c PMCR(13) step size in polar angle c PMCR(14) number of gridpoints in polar angle c PMCR(15) magnitude of minimum azimuthal angle c PMCR(16) step size in azimuthal angle c PMCR(17) number of gridpoints in azimuthal angle c PMCR(18) number of trajectories per grid point c c------------------------------------------------------------- c logical variables c integer nubuf,nvtx,ntdumy,iz0 parameter(nubuf=3) real p1(3),vec4(4),vec5(5),vec6(6),pp(3),pe(3),vec2(2),VEC3(3) real p,pmin,delp,xmin,delx,ymin,dely,themin,delthe,phimin, 1 delphi,dx,dy,the,phi,pxl,pyl,pzl,x0(3),ubuf(nubuf),pt,pt2 real themax,pinv,weight,d1,dc,ds,phir,phi1,thel real rhomin,delrho,tphimin,tdelphi,rho,tphi,dm,dmr,zmin,delz integer icall,i0,ipmax,ixmax,iymax,ithemax,iphimax,intr,ik,ip,ix, 1 iy,ithe,iphi,ipold,k,idpart,itr,itype,i1,iz,irho,irhomax,izmax integer j1,j2,j3 logical has_weight save icall,ipold C if(icall.lt.0)then call hbook1(951,' e+ z ',200,-3.,3.,0.) call hbook1(952,' e+ rho ',200,-2.,2.,0.) call hbook1(953,' e+ mominv ',200,0.,25.,0.) call hbook1(954,' e+ the ',200,0.,90.,0.) call hbook1(955,' e+ phi ori ',200,-30.,30.,0.) call hbook1(956,' e+ mom ',200,0.,6.,0.) call hbook1(957,' e+ phiadd ',200,0.,360.,0.) call hbook1(958,' e+ phinew ',200,-200.,600.,0.) call hbook1(961,' e- z ',200,-3.,3.,0.) call hbook1(962,' e- rho ',200,-2.,2.,0.) call hbook1(963,' e- mominv ',200,0.,25.,0.) call hbook1(964,' e- the ',200,0.,90.,0.) call hbook1(965,' e- phi ori ',200,-30.,30.,0.) call hbook1(966,' e- mom ',200,0.,6.,0.) call hbook1(967,' e- phiadd ',200,0.,360.,0.) call hbook1(968,' e- phinew ',200,-200.,600.,0.) call hbook2(971,' e+ XY PAD1 ',100,-120.,120.,100,-150.,150.,0.) call hbook2(972,' e+ XY PAD2 ',100,-120.,120.,100,-150.,150.,0.) call hbook2(973,' e+ XY PAD3 ',100,-120.,120.,100,-150.,150.,0.) call hbook2(974,' e+ XY PAD4 ',100,-120.,120.,100,-150.,150.,0.) call hbook2(975,' e+ XY PAD5 ',100,-120.,120.,100,-150.,150.,0.) call hbook2(976,' e+ XY PAD6 ',100,-120.,120.,100,-150.,150.,0.) call hbook2(977,' e+ XY PAD7 ',100,-120.,120.,100,-150.,150.,0.) call hbook2(978,' e+ XY PAD8 ',100,-120.,120.,100,-150.,150.,0.) call hbook2(979,' e+ XY PAD9 ',100,-120.,120.,100,-150.,150.,0.) call hbook2(980,' e+ XY PAD10',100,-120.,120.,100,-150.,150.,0.) call hbook2(981,' e- XY PAD1 ',100,-120.,120.,100,-150.,150.,0.) call hbook2(982,' e- XY PAD2 ',100,-120.,120.,100,-150.,150.,0.) call hbook2(983,' e- XY PAD3 ',100,-120.,120.,100,-150.,150.,0.) call hbook2(984,' e- XY PAD4 ',100,-120.,120.,100,-150.,150.,0.) call hbook2(985,' e- XY PAD5 ',100,-120.,120.,100,-150.,150.,0.) call hbook2(986,' e- XY PAD6 ',100,-120.,120.,100,-150.,150.,0.) call hbook2(987,' e- XY PAD7 ',100,-120.,120.,100,-150.,150.,0.) call hbook2(988,' e- XY PAD8 ',100,-120.,120.,100,-150.,150.,0.) call hbook2(989,' e- XY PAD9 ',100,-120.,120.,100,-150.,150.,0.) call hbook2(990,' e- XY PAD10',100,-120.,120.,100,-150.,150.,0.) c IF(NREF(1).EQ.2)THEN c OPEN (UNIT=24,file='/u/hades/geant/evt/OMEGAS.EVT', c 1 STATUS='OLD') c END IF c IF(NREF(1).EQ.3)THEN c OPEN (UNIT=25,file='/d/hades2/pions/Pions25kw.evt', c 1 STATUS='OLD') c END IF end if dm=.000511 idpart=int(PMCR(2)) icall=icall+1 if(nref(2).eq.1)then xmin=PMCR(3) delx=PMCR(4) ixmax=int(PMCR(5)) ymin=PMCR(6) dely=PMCR(7) iymax=int(PMCR(8)) else if(nref(2).eq.2.or.nref(2).eq.3)then zmin=PMCR(3) delz=PMCR(4) izmax=int(PMCR(5)) rhomin=PMCR(6) delrho=PMCR(7) irhomax=int(PMCR(8)) end if pmin=PMCR(9) delp=PMCR(10) ipmax=int(PMCR(11)) themin=PMCR(12) delthe=PMCR(13) ithemax=int(PMCR(14)) themax=themin+float(ithemax-1)*delthe phimin=PMCR(15) delphi=PMCR(16) iphimax=int(PMCR(17)) intr=int(PMCR(18)) iref0(1)=int(pmcr(5)) iref0(2)=int(pmcr(8)) iref0(3)=int(pmcr(11)) iref0(4)=int(pmcr(14)) iref0(5)=int(pmcr(17)) nb(1)=iref0(5) do ik=2,4 nb(ik)=nb(ik-1)*iref0(6-ik) end do ir(5)=icall-1 do ik=1,4 ir(5-ik)=mod(ir(6-ik),nb(5-ik)) iref(ik)=(ir(6-ik)-ir(5-ik))/nb(5-ik)+1 end do iref(5)=ir(1)+1 if(iref(1).ne.ipold)then ipold=iref(1) write(6,101)(iref0(k),k=1,5),(iref(k),k=1,5) 101 format( ' ref0 101 ',5i3,' ref ',5i3) end if if(nref(2).eq.1)then ix=iref(1) iy=iref(2) else if(nref(2).eq.2)then iz=iref(1) irho=iref(2) end if ip=iref(3) ithe=iref(4) iphi=iref(5) if(nref(1).le.2)then c p=pmin+delp*float(ip-1) p=1./(pmin+delp*float(ip-1)) the=themin+delthe*float(iTHE-1) phi=PHImin+delPHI*float(iPHI-1) ubuf(1)=p ubuf(2)=the ubuf(3)=phi if(nref(2).eq.1)then x0(1)=xmin+delx*float(ix-1) x0(2)=ymin+dely*float(iy-1) x0(3)=0. else if(nref(2).eq.2)then rho=rhomin+delrho*float(irho-1) x0(1)=rho*cos(degrad*(phi+90.)) x0(2)=rho*sin(degrad*(phi+90.)) x0(3)=zmin+delz*float(iz-1) end if p1(1)=p*sin(degrad*the)*cos(degrad*phi) p1(2)=p*sin(degrad*the)*sin(degrad*phi) p1(3)=p*cos(degrad*the) thel=degrad*the c write(6,201)icall,(x0(j1),j1=1,3),(p1(j2),j2=1,3),the c 201 format(' refev=',i4,' vert ',3f6.2,' px ',3f6.3,f6.2) CALL GSVERT(X0,0,0,0.,0,NVTX) ! store do i0=1,intr CALL GSKINE(P1,idpart,nvtx,ubuf,nubuf,ntdumy) end do else if(nref(1).eq.3)then 200 continue read(24,*,end=910)intr if(mod(icall,1000).eq.0)write(6,800)icall 800 format(' event number ',i8) if(intr.eq.2)then call nran(vec4,4) phir=vec4(4)*360. do itr=1,2 read(24,*)p1,itype pt2=p1(1)**2+p1(2)**2 pt=sqrt(pt2) p=sqrt(pt2+p1(3)**2) the=raddeg*atan2(pt,p1(3)) * the=raddeg*asin(pt/p) phi=raddeg*atan2(p1(2),p1(1)) dc=cos(degrad*the) ds=sin(degrad*the) phi1=phi+phir if(phi1.gt.180.)phi1=phi1-360. pp(1)=p*ds*cos(degrad*phi1) pp(2)=p*ds*sin(degrad*phi1) pp(3)=p*dc if(nref(2).eq.1)then x0(1)=xmin+delx*float(ixmax-1)*vec4(1) x0(2)=ymin+dely*float(iymax-1)*vec4(2) x0(3)=0. else if(nref(2).eq.2)then rho=.5*delrho*float(irhomax-1)*sqrt(vec4(1)) if(vec4(2).lt.0.5)rho=-rho x0(1)=rho*cos(degrad*(phi+90.)) x0(2)=rho*sin(degrad*(phi+90.)) x0(3)=zmin+delz*float(izmax-1)*vec4(3) else if(nref(2).eq.3)then rho=.5*delrho*float(irhomax-1)*sqrt(vec4(1)) if(vec4(2).lt.0.5)rho=-rho x0(1)=rho*cos(degrad*(phi+90.)) x0(2)=rho*sin(degrad*(phi+90.)) iz=mod(icall,izmax) x0(3)=zmin+delz*float(iz) end if call gsvert(x0,0,0,0.,0,nvtx) d1=1./p if(itype.eq.2)then call hf1(951,x0(3),1.) call hf1(952,rho,1.) call hf1(953,d1,1.) call hf1(954,the,1.) call hf1(955,phi,1.) call hf1(956,p,1.) call hf1(957,phir,1.) call hf1(958,phi1,1.) else if(itype.eq.3)then call hf1(961,x0(3),1.) call hf1(962,rho,1.) call hf1(963,d1,1.) call hf1(964,the,1.) call hf1(965,phi,1.) call hf1(966,p,1.) call hf1(967,phir,1.) call hf1(968,phi1,1.) end if ubuf(1)=p ubuf(2)=the ubuf(3)=phi1 call gskine(pp,itype,nvtx,ubuf,nubuf,ntdumy) end do end if go to 900 910 rewind(unit=24) print *, ' rewinding unit 24 ' go to 200 else if(nref(1).eq.4)then 300 continue has_weight=.false. read(25,*,end=911)intr if(mod(icall,1000).eq.0)write(6,801)icall 801 format(' event number ',i8) if(intr.lt.0)then intr=-intr has_weight=.true. else has_weight=.false. end if if(intr.gt.0)then call nran(vec4,4) phir=vec4(4)*360. do itr=1,intr weight=0. if(has_weight)then read(25,*)p1,itype,weight else read(25,*)p1,itype end if pt2=p1(1)**2+p1(2)**2 pt=sqrt(pt2) p=sqrt(pt2+p1(3)**2) * the=raddeg*asin(pt/p) the=raddeg*atan2(pt,p1(3)) phi=raddeg*atan2(p1(2),p1(1)) dc=cos(degrad*the) ds=sin(degrad*the) phi1=phi+phir if(phi1.gt.180.)phi1=phi1-360. if(nref(2).eq.1)then x0(1)=xmin+delx*float(ixmax-1)*vec4(1) x0(2)=ymin+dely*float(iymax-1)*vec4(2) x0(3)=0. else if(nref(2).eq.2)then rho=.5*delrho*float(irhomax-1)*sqrt(vec4(1)) if(vec4(2).lt.0.5)rho=-rho x0(1)=rho*cos(degrad*(phi+90.)) x0(2)=rho*sin(degrad*(phi+90.)) x0(3)=zmin+delz*float(izmax-1)*vec4(3) else if(nref(2).eq.3)then rho=.5*delrho*float(irhomax-1)*sqrt(vec4(1)) if(vec4(2).lt.0.5)rho=-rho x0(1)=rho*cos(degrad*(phi+90.)) x0(2)=rho*sin(degrad*(phi+90.)) iz=mod(icall,izmax) x0(3)=zmin+delz*float(iz) end if call gsvert(x0,0,0,0.,0,nvtx) p1(1)=p*ds*cos(degrad*phi1) p1(2)=p*ds*sin(degrad*phi1) p1(3)=p*dc ubuf(1)=p ubuf(2)=the ubuf(3)=phi call gskine(p1,itype,nvtx,ubuf,nubuf,ntdumy) end do end if go to 900 911 rewind(unit=25) print *, ' rewinding unit 25 ' go to 300 end if 900 continue return end ************************************************************************