subroutine gukine c c c This is a subroutine to generate (x,y,theta,phi,delta-p) at c the focal plane for proton from the H(e,e'p) reaction. c It makes use of a simple kinematics routine, some knowledge c of the extended target acceptances, and polynomials for the c transport through the spectrometer from John Lerose. c implicit none integer ikine,itra,istak,ivert,ipart,itrtyp,napart,ipaold real pkine,amass,charge,tlife,vert,pvert common/gckine/ikine,pkine(10),itra,istak,ivert,ipart,itrtyp + ,napart(5),amass,charge,tlife,vert(3),pvert(4),ipaold integer nubuf parameter (nubuf=10) real ubuf(nubuf),rndm(3) integer ntbeam,nttarg,nvtx save nvtx real plab(3),ptot,etot real*8 ang,angr,phi,phir,psir real mp,ea,pa,remain real r2,r3 integer nt,ierri,i0,i1,i2,i3 integer i,j,k,ichoice,ichoice2,iremain real*8 rotmat,xyz(3),xyznew(3),termang integer ic,ii,jj integer*4 junk1,ikinsetting,isectrack real*8 xfp,yfp,tthfp,tphfp,pfp,junk2,junk3,thfp,phfp real*8 e0(10),eang(10),hang(10),targ_thick(10) real*8 xtgt,ytgt,thtgt,phtgt,ptgt,dptgt real yf,ff,thetaf,dptgt2 character*80 junkline c c The following variables are the SIMC variables c real*8 x,y,z !(cm) real*8 dpp !delta p/p (%) real*8 dxdz,dydz !X,Y slope in spectrometer real*8 x_fp,y_fp,dx_fp,dy_fp !Focal plane values to return real*8 p_spec,th_spec !spectrometer setting real*8 fry !vertical position@tgt (+y=down) real*8 pathlen logical ms_flag !mult. scattering flag logical wcs_flag !wire chamber smearing flag logical decay_flag !check for particle decay logical ok_spec !true if particle makes it real*8 resmult,m2 c include 'fpp_local.h' include 'geant_local.h' c c include 'parameter.h' c include 'espace_type.h' c include 'detector.h' c include 'transport.h' c include 'option.h' c common/kincom/rotmat(3,3) c nhu1=0 nhx1=0 nhv1=0 nhu2=0 nhx2=0 nhv2=0 nhu3=0 nhx3=0 nhv3=0 nhu4=0 nhx4=0 nhv4=0 c call grndm ( rndm , 3 ) c 111 format(a80) iremain=nevent/1000 remain=nevent/1000.0 c if(iremain.eq.remain) write(*,*)'nevent =',nevent write(*,*)'nevent =',nevent if(nevent.eq.0) then c write(*,*)'Getting kinematics setting ...' open(unit=1,file='geant_kinematics.dat',type='UNKNOWN') read(1,*)ikinsetting read(1,*)isectrack close(unit=1) if(isectrack.eq.1) then sectrack=.true. else sectrack=.false. endif open(unit=1,file='hdr_gep.dat',status='old') do i=1,10 read(1,*)e0(i),eang(i),hang(i),targ_thick(i) enddo einc=e0(ikinsetting) hrse_ang=eang(ikinsetting) hrsh_ang=-1.0*hang(ikinsetting) trg_thk=targ_thick(ikinsetting) close(unit=1) endif c c Now we have the incident electron energy and the angles of c the two spectrometers, as well as the target thickness. c We need to use a) kinematics routines and b) some knowledge c of the extended target acceptances to choose the x,y,theta,phi, c and dpmom at the target. c call kincalc(einc,hrse_ang,hrsh_ang,trg_thk,xtgt,ytgt, $ thtgt,phtgt,ptgt,dptgt) c p_spec=ptgt c if(nevent.le.2)write(*,*)'******** mom = ',p_spec,' ******' c call grndm(rndm,3) dptgt2=rndm(1)*0.200-0.100 xfp=rndm(2)*120.00-60.00 phfp=rndm(3)*80.0-40.0 c phfp=3.0 c phfp=rndm(3)*0.00001-0.000005 c call grndm(rndm,3) yfp=rndm(1)*60.00-30.00 thfp=rndm(2)*80.0-40.0 c thfp=3.0 c thfp=rndm(3)*0.00001-0.000005 c c write(6,*)'x,phi,y,theta,ok_spec =' c write(6,*)xfp,phfp,yfp,thfp,ok_spec c do ii=1,20 ntuple_array(ii)=0.0 enddo ntuple_array(11)=real(xfp) ntuple_array(12)=real(yfp) ntuple_array(13)=real(thfp) ntuple_array(14)=real(phfp) c pcentral=ptgt pfp=dptgt2 c c Finally, we convert this information over to a format that GEANT likes. c 1135 dpmom=pfp pfp=pcentral*(1.0+pfp) pmom=pfp ea=sqrt(pfp**2+938.2796**2) tinit=ea-938.2796 vert(1)=xfp vert(2)=yfp vert(3)=160.0 ntbeam = 0.0 nttarg = 0.0 ubuf(1) = 0.0 1000 continue c ipart = 14 ! geant pid (8=pi+,9,pi-, 14 =p) ptot=pfp/1000. call grndm ( rndm , 3 ) angr=thfp*3.14159265/180.0 phir=phfp*3.14159265/180.0 psir=datan(dtan(phir)*dcos(angr)) ang=thfp phi=phfp c c these are the actual parameters for the track. c xinit=vert(1) yinit=vert(2) thinit=angr phiinit=phir c sptransport.l.particle.fp_h.ph=dtan(angr) c sptransport.l.particle.fp_h.th=dtan(phir) c c sptransport.l.particle.fp_h.x=vert(1)/100.0 c c sptransport.l.particle.fp_h.y=vert(2)/100.0 c c c next we misalign the track to simulate misalignment of the c entire space frame with respect to the vdc's c c just include the translational offsets to start c c write(*,*)vert(1),vert(2),angr,phir,psir xofff=0.0 yofff=0.0 thofff=0.0 phofff=0.0 psofff=0.0 c c Define the inverse Euler rotation for (thofff,phiofff,psiofff) c c rotmat(1,1)=dcos(psofff)*dcos(thofff)+dsin(psofff)*dsin(thofff) $ *dsin(phofff) rotmat(1,2)=-dcos(phofff)*dsin(thofff) rotmat(1,3)=-dsin(psofff)*dcos(thofff)+dcos(psofff)*dsin(thofff) $ *dcos(phofff) rotmat(2,1)=dcos(psofff)*dsin(thofff)-dsin(psofff)*dcos(thofff) $ *dsin(phofff) rotmat(2,2)=dcos(phofff)*dcos(thofff) rotmat(2,3)=-dsin(psofff)*dsin(thofff)-dcos(psofff)*dcos(thofff) $ *dsin(phofff) rotmat(3,1)=dsin(psofff)*dcos(phofff) rotmat(3,2)=dsin(phofff) rotmat(3,3)=dcos(psofff)*dcos(phofff) c c vert(1)=vert(1)-xofff vert(2)=vert(2)-yofff c xyz(1)=dsin(psir) xyz(2)=dcos(psir)*dsin(angr) xyz(3)=dcos(psir)*dcos(angr) c do ic=1,3 xyznew(ic)=xyz(1)*rotmat(ic,1)+ & xyz(2)*rotmat(ic,2)+xyz(3)*rotmat(ic,3) enddo angr=datan(xyznew(2)/xyznew(3)) termang=xyznew(3)/dcos(angr) if(termang.gt.1.0) termang=1.0 if(termang.lt.-1.0) termang=-1.0 psir=dacos(termang)*abs(xyznew(1))/xyznew(1) phir=datan(dtan(psir)/dcos(angr)) c write(*,*)vert(1),vert(2),angr,phir,psir call gsvert ( vert,ntbeam,nttarg,ubuf,0,nvtx ) etot = ea - 938.2796 plab(1) = ptot*dsin(psir) plab(2) = ptot*dsin(angr)*dcos(psir) plab(3) = ptot*dcos(angr)*dcos(psir) call hfill ( 100, etot, 0., 1. ) call gskine ( plab,ipart,nvtx,ubuf,0,nt ) if ( nt.le.0 ) then write ( 6,* ) ' gukine: error defining track' write ( 6,* ) ' i=',i,' nt=',nt stop end if nevent=nevent+1 return end subroutine kincalc(e0,eang,hang,trg,xtgt,ytgt, $ thtgt,phtgt,ptgt,dptgt) implicit none real*8 e0,eang,hang,trg,xtgt,ytgt,thtgt,phtgt,ptgt,dptgt real*8 fg,gf integer i,j,k real*8 mt,mtg,mr,mpi,mn,mp,me,pi,mhe,alpha real*8 escat,pscat,pcentral,thetae,phie,thetap,phip real rndm(3) fg=3.14159265/180.0 gf=1.0/fg me = 0.511 mpi = 139.57 mp = 938.2796 mn = 939.5731 mhe = 2808.41 mt = mp mtg = mt/1.e3 alpha=1./137. 1432 call grndm ( rndm , 3 ) escat=mp/(1.0+mp/e0-cos(fg*eang)) pscat=sqrt(e0**2+escat**2- $ 2.0*e0*escat*cos(fg*eang)) pcentral=pscat c First, we assume that whichever are is the more backward will c determine the acceptance. We then randomly choose the theta c and phi for the arm that determines the acceptance within the c usual full HRS acceptance and then determine from kinematics c what the corresponding theta and phi are for the other arm. if (abs(eang).gt.abs(hang)) then c phie=-.025+rndm(1)*.050 c thetae=-.009+rndm(2)*.018 if(pcentral.ge.5000) then c phie=-.00+rndm(1)*.00 c thetae=-.00+rndm(2)*.00 phie=-.130+rndm(1)*.260 thetae=-.065+rndm(2)*.130 else if(pcentral.ge.3000.and.pcentral.lt.5000) then phie=-.067+rndm(1)*.135 thetae=-.034+rndm(2)*.067 else phie=-.087+rndm(1)*.174 thetae=-.044+rndm(2)*.087 endif c phie=-.00+rndm(1)*.00 c thetae=-.00+rndm(2)*.00 escat=mp/(1.0+mp/e0-cos(fg*eang+thetae)*cos(phie)) pscat=sqrt(e0**2+escat**2- $ 2.0*e0*escat*cos(fg*eang+thetae)*cos(phie)) phip=dasin(escat*dsin(phie)/pscat) thetap=dasin((escat*sin(fg*eang+thetae)*cos(phie))/ $ (pscat*cos(phip)))-fg*hang else c Hadron arm defining acceptance 1221 call grndm ( rndm, 3) phie=-.080+rndm(1)*.160 thetae=-.030+rndm(2)*.060 c phie=-.00+rndm(1)*.00 c thetae=-.00+rndm(2)*.00 escat=mp/(1.0+mp/e0-cos(fg*eang+thetae)*cos(phie)) pscat=sqrt(e0**2+escat**2- $ 2.0*e0*escat*cos(fg*eang+thetae)*cos(phie)) phip=-1.0*dasin(escat*sin(phie)/pscat) thetap=dasin((escat*sin(fg*eang+thetae)*cos(phie))/ $ (pscat*cos(phip)))-fg*hang if(abs(thetap).gt.0.030.or.abs(phip).gt.0.065) goto 1221 endif ptgt=pscat c write(*,*)escat,eang,pscat,hang dptgt=(pcentral-pscat)/pcentral c c Following statement to fill just the high and low dp bins c which have very low statistics. This is normally commented c out. c c if(abs(dptgt).le.0.030) goto 1432 c thtgt=thetap phtgt=phip xtgt=0.0 ytgt=0.0 return end