|
version 1.1, 2005/09/13 15:09:16
|
version 1.13, 2006/01/26 21:19:08
|
|
|
|
| subroutine gukine | subroutine gukine |
| c | c |
| c jpsullivan may 20, 1993: a simple example of gukine |
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 | c |
| implicit none | implicit none |
| integer ikine,itra,istak,ivert,ipart,itrtyp,napart,ipaold | integer ikine,itra,istak,ivert,ipart,itrtyp,napart,ipaold |
|
|
|
| save nvtx | save nvtx |
| real plab(3),ptot,etot | real plab(3),ptot,etot |
| real*8 ang,angr,phi,phir,psir | real*8 ang,angr,phi,phir,psir |
| real mp,ea,pa |
real mp,ea,pa,remain |
| real r2,r3 | real r2,r3 |
| integer nt,ierri,i0,i1,i2,i3 | integer nt,ierri,i0,i1,i2,i3 |
| integer i,j,k,ichoice,ichoice2 |
integer i,j,k,ichoice,ichoice2,iremain |
| real*8 rotmat,xyz(3),xyznew(3),termang | real*8 rotmat,xyz(3),xyznew(3),termang |
| integer ic,ii,jj | integer ic,ii,jj |
| integer*4 junk1 |
integer*4 junk1,ikinsetting,isectrack |
| real*8 xfp,yfp,tthfp,tphfp,pfp,junk2,junk3,thfp,phfp | real*8 xfp,yfp,tthfp,tphfp,pfp,junk2,junk3,thfp,phfp |
| real*8 e0(10),eang(10),hang(10),targ_thick(10) | real*8 e0(10),eang(10),hang(10),targ_thick(10) |
| real*8 escat,pscat |
real*8 xtgt,ytgt,thtgt,phtgt,ptgt,dptgt |
| integer*4 ikinsetting |
real yf,ff,thetaf |
| real*8 fg |
|
| real dptgt |
|
| character*80 junkline | character*80 junkline |
| c | 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 'fpp_local.h' |
| include 'geant_local.h' | include 'geant_local.h' |
| c | c |
|
|
|
| c include 'transport.h' | c include 'transport.h' |
| c include 'option.h' | c include 'option.h' |
| c | c |
| c common/kincom/rotmat(3,3) |
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 | c |
| 111 format(a80) | 111 format(a80) |
| |
iremain=nevent/1000 |
| |
remain=nevent/1000.0 |
| |
c if(iremain.eq.remain) write(*,*)'nevent =',nevent |
| write(*,*)'nevent =',nevent | write(*,*)'nevent =',nevent |
| mp=938.2785 |
|
| fg=3.14159265/180.0 |
|
| | |
| if(nevent.eq.0) then | if(nevent.eq.0) then |
| c write(*,*)'Which kinematics setting ?' |
c write(*,*)'Getting kinematics setting ...' |
| open(unit=1,file='geant_kinematics.dat',type='UNKNOWN') | open(unit=1,file='geant_kinematics.dat',type='UNKNOWN') |
| read(1,*)ikinsetting | read(1,*)ikinsetting |
| write(*,*)'Kinematics Setting: ',ikinsetting |
read(1,*)isectrack |
| close(unit=1) | close(unit=1) |
| |
if(isectrack.eq.1) then |
| |
sectrack=.true. |
| |
else |
| |
sectrack=.false. |
| |
endif |
| open(unit=1,file='hdr_gep.dat',status='old') | open(unit=1,file='hdr_gep.dat',status='old') |
| do i=1,10 | do i=1,10 |
| read(1,*)e0(i),eang(i),hang(i),targ_thick(i) | read(1,*)e0(i),eang(i),hang(i),targ_thick(i) |
| write(*,*)e0(i),eang(i),hang(i),targ_thick(i) |
|
| enddo | enddo |
| einc=e0(ikinsetting) | einc=e0(ikinsetting) |
| hrse_ang=eang(ikinsetting) | hrse_ang=eang(ikinsetting) |
| hrsh_ang=-1.0*hang(ikinsetting) | hrsh_ang=-1.0*hang(ikinsetting) |
| escat=mp/(1.0+mp/einc-cos(fg*hrse_ang)) |
trg_thk=targ_thick(ikinsetting) |
| pscat=sqrt(einc**2+escat**2- |
|
| $ 2.0*einc*escat*cos(fg*hrse_ang)) |
|
| pcentral=pscat |
|
| close(unit=1) | close(unit=1) |
| endif | endif |
| | |
| nevent=nevent+1 |
c |
| |
c Now we have the incident electron energy and the angles of |
| call grndm ( rndm , 3 ) |
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) |
| | |
| dptgt = -0.05+rndm(1)*0.100 |
|
| c | c |
| ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
c At this point, after using polynomials from Lerose, we should |
| |
c have x,y,theta,phi, and dpmom at the focal plane. |
| c | c |
| c make the vertex at z=0, (x,y) chosen randomly in a 20 cm x 20 cm square |
c write(*,*)'Back from kincalc ...' |
| |
c write(*,*)xtgt,ytgt,thtgt,phtgt,dptgt |
| |
c x(1)=tan(phtgt) |
| |
c x(2)=ytgt |
| |
c x(3)=tan(thtgt) |
| |
c x(4)=dptgt |
| |
y_target=ytgt*100.0 |
| |
c |
| |
c Here is where we call the SIMC routines to calculate the focal plane |
| |
c quantities from the target quantities. |
| |
c |
| |
c xfp=xfinal(x,4)*100.0 |
| |
c phfp=atan(thetaf(x,4))/3.14159265*180.0 |
| |
c yfp=yf(x,4)*100.0 |
| |
c thfp=atan(ff(x,4))/3.14159265*180.0 |
| |
p_spec=ptgt |
| |
th_spec=hrsh_ang |
| |
dpp=dptgt*100.00 |
| |
x=xtgt |
| |
y=ytgt |
| |
z=0.0 |
| |
dxdz=tan(phtgt*3.14159265/180.0) |
| |
dydz=tan(thtgt*3.14159265/180.0) |
| |
if(nevent.le.2)write(*,*)'******** mom = ',p_spec,' ******' |
| |
c write(*,*)x,y,z |
| |
c write(*,*)dxdz,dydz |
| |
c |
| |
c Call to SIMC routine |
| |
c |
| |
decay_flag=.false. |
| |
wcs_flag=.false. |
| |
ms_flag=.false. |
| |
m2=(938.2796**2) |
| |
call mc_hms (p_spec, th_spec, dpp, x, y, z, dxdz, dydz, |
| |
> x_fp, dx_fp, y_fp, dy_fp, m2, |
| |
> ms_flag, wcs_flag, decay_flag, resmult, fry, |
| |
> ok_spec, pathlen) |
| |
c write(*,*)'Momentum = ',p_spec,' Angle = ',th_spec |
| |
c write(*,*)x,y,z |
| |
c write(*,*)dxdz,dydz |
| |
c write(6,*)x_fp,dx_fp,y_fp,dy_fp,ok_spec |
| |
c write(*,*)m2,ms_flag,wcs_flag,decay_flag |
| |
c write(*,*)resmult,fry,ok_spec,pathlen |
| |
|
| |
xfp=x_fp |
| |
phfp=atan(dx_fp)/3.14159265*180.0 |
| |
yfp=y_fp*100.0 |
| |
thfp=atan(dy_fp)/3.14159265*180.0*100.0 |
| |
c write(6,*)'x,phi,y,theta,ok_spec =' |
| |
c write(6,*)xfp,phfp,yfp,thfp,ok_spec |
| c | c |
| c write ( 6,* ) ' gukine called',icount |
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 | c |
| vert(1)=-65.0+rndm(2)*110.0 |
pcentral=ptgt |
| vert(2)=-6.0+rndm(3)*12.0 |
|
| vert(3) = 0.0 |
|
| pfp=dptgt | pfp=dptgt |
| |
c |
| |
c |
| |
c write(*,*)'Focal plane quantities ...' |
| |
c write(*,*)xfp,yfp,thfp,phfp,pcentral,pfp |
| | |
| c write(*,*)xfp,phfp,yfp,thfp |
c xfp=xtgt ! in cm |
| |
c yfp=ytgt ! in cm |
| |
c thfp=thtgt/3.14159265*180.0 ! in degrees |
| |
c phfp=phtgt/3.14159265*180.0 ! in degrees |
| |
c pcentral=ptgt ! in MeV/c |
| |
c pfp=dptgt ! fraction of central momentum |
| |
c |
| |
c Finally, we convert this information over to a format that GEANT likes. |
| |
c |
| | |
| 1135 dpmom=pfp | 1135 dpmom=pfp |
| pfp=pcentral*(1.0+pfp) | pfp=pcentral*(1.0+pfp) |
| pmom=pfp | pmom=pfp |
| ea=sqrt(pfp**2+mp**2) |
ea=sqrt(pfp**2+938.2796**2) |
| tinit=ea-mp |
tinit=ea-938.2796 |
| |
vert(1)=xfp |
| |
vert(2)=yfp |
| |
vert(3)=0.0 |
| ntbeam = 0.0 | ntbeam = 0.0 |
| nttarg = 0.0 | nttarg = 0.0 |
| ubuf(1) = 0.0 | ubuf(1) = 0.0 |
| c write(6,*)'***********************************************' |
|
| c write ( 6,* ) ' guout: xyz =',vert(1),vert(2),vert(3) |
|
| c |
|
| cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|
| c |
|
| 1000 continue | 1000 continue |
| c | c |
| ipart = 14 ! geant pid (8=pi+,9,pi-, 14 =p) | ipart = 14 ! geant pid (8=pi+,9,pi-, 14 =p) |
| ptot=pfp/1000. | ptot=pfp/1000. |
| call grndm ( rndm , 3 ) | call grndm ( rndm , 3 ) |
| ang=-2.2+rndm(1)*4.4 |
angr=thfp*3.14159265/180.0 |
| phi=-7.0+rndm(2)*12.5 |
phir=phfp*3.14159265/180.0 |
| angr=ang*3.14159265/180.0 |
|
| phir=phi*3.14159265/180.0 |
|
| psir=datan(dtan(phir)*dcos(angr)) | psir=datan(dtan(phir)*dcos(angr)) |
| c ang=(-3.66+rndm(1)*7.32)/100. |
|
| c phi=(-9.98+rndm(2)*19.96)/100. |
ang=thfp |
| c ang=5.0 |
phi=phfp |
| c phi=5.0 |
|
| c | c |
| c these are the actual parameters for the track. | c these are the actual parameters for the track. |
| c | 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 ) | call gsvert ( vert,ntbeam,nttarg,ubuf,0,nvtx ) |
| etot = ea - mp |
etot = ea - 938.2796 |
| plab(1) = ptot*dsin(psir) | plab(1) = ptot*dsin(psir) |
| plab(2) = ptot*dsin(angr)*dcos(psir) | plab(2) = ptot*dsin(angr)*dcos(psir) |
| plab(3) = ptot*dcos(angr)*dcos(psir) | plab(3) = ptot*dcos(angr)*dcos(psir) |
| | |
| write(6,*)'initial momentum = ',ptot |
|
| write(6,*)'initial components are',plab(1),plab(2),plab(3) |
|
| c |
|
| c write ( 6,902 ) ipart |
|
| c902 format ( ' **********************************', |
|
| c x ' gukine: start new track: ipart=',i5 ) |
|
| c |
|
| c write ( 6,* ) ' filling etot histo' |
|
| call hfill ( 100, etot, 0., 1. ) | call hfill ( 100, etot, 0., 1. ) |
| c |
|
| c give the track to geant |
|
| c |
|
| c write ( 6,* ) ' gskine called' |
|
| call gskine ( plab,ipart,nvtx,ubuf,0,nt ) | call gskine ( plab,ipart,nvtx,ubuf,0,nt ) |
| c |
|
| if ( nt.le.0 ) then | if ( nt.le.0 ) then |
| write ( 6,* ) ' gukine: error defining track' | write ( 6,* ) ' gukine: error defining track' |
| write ( 6,* ) ' i=',i,' nt=',nt | write ( 6,* ) ' i=',i,' nt=',nt |
| stop | stop |
| end if | end if |
| c |
|
| c print the vertex and track info |
nevent=nevent+1 |
| c |
|
| c call gpvert(0) |
return |
| c call gpkine(0) |
|
| c |
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 | return |
| c |
|
| end | end |
Return to
gukine.f
CVS log
Up to [HallC] / geant_gep / src