Skip to content
Snippets Groups Projects
Commit 5c71ddcd authored by Sylvester Joosten's avatar Sylvester Joosten
Browse files

Merge branch 'get_event_file-update' into 'master' 

Modify get_file_event and event

See merge request jpsi007/simc-file-input!2
parents bbbea151 80e4066d
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
simc/event.f
+ 1
− 1
View file @ 5c71ddcd
......@@ -299,7 +299,7 @@ C modified 5/15/06 for poinct
endif !not (doing_hyd_elast)
if ( doing_2pi) then
ok_2pi = .true.
call get_file_event(spec%e%theta,spec%p%theta,
call get_file_event(electron_arm,spec%e%theta,spec%p%theta,
> vertex%e%xptar,vertex%e%yptar,vertex%e%p,vertex%e%E,
> vertex%p%xptar,vertex%p%yptar,vertex%p%p,vertex%p%E,
> main%gen_weight)
......
simc/get_file_event.f
+ 63
− 28
View file @ 5c71ddcd
subroutine get_file_event(th_spec_e,th_spec_p,
subroutine get_file_event(e_arm,th_spec_e,th_spec_p,
> dxdz,dydz,e_mom,e_E,dxdzp,dydzp,p_mom,p_E,weight)
c
c input variables:
c electron_arm : 1 means HMS is electron arm
c 5 means SHMS is the electron arm
c th_spec_e : central spec angle for electron (rad)
c th_spec_p : central spec angle for proton (rad)
c output variables:
......@@ -20,12 +22,14 @@ c
implicit none
include 'simulate.inc'
c
real*8 e_4v(4), e_4vr(4)
real*8 p_4v(4), p_4vr(4)
integer*4 e_arm,ii
real*8 SHMS_4v(4), e_4vr(4), e_4v(4)
real*8 HMS_4v(4), p_4vr(4), p_4v(4)
real*8 cth_hms, sth_hms, cth_shms, sth_shms
real*8 th_spec_e
real*8 th_spec_p
real*8 e_mom,e_E,dxdz,dydz,e_vz
real*8 sth_elec,cth_elec,sth_prot,cth_prot
real*8 e_mom,e_E,dxdz,dydz,e_vz,SHMS_vz,HMS_vz
real*8 p_mom,p_E,dxdzp,dydzp,p_vz
real*8 w, weight
character*80 multpifile
......@@ -47,40 +51,69 @@ c
c
c
end_of_2pi_file = .false.
read(51,*,end=999,err=999) e_4v,e_vz,p_4v,p_vz,w
888 continue
read(51,*,end=999,err=999) HMS_4v,HMS_vz,SHMS_4v,SHMS_vz,w
count = count + 1
c if (count .gt. 100) goto 999
c Rotate 4-vectors to their respective frames
c This is a rotation around the x-axis, which points downwards
c --> the HMS is at negative angles, and the SHMS at positive angles
if (e_arm .eq. 1) then
cth_elec = cos(th_spec_e)
sth_elec = sin(-th_spec_e)
cth_prot= cos(th_spec_p)
sth_prot = sin(th_spec_p)
e_vz = HMS_vz * 1.
p_vz = SHMS_vz * 1.
e_4v(1) = HMS_4v(1)
p_4v(1) = SHMS_4v(1)
e_4v(2) = HMS_4v(2)
p_4v(2) = SHMS_4v(2)
e_4v(3) = HMS_4v(3)
p_4v(3) = SHMS_4v(3)
e_4v(4) = HMS_4v(4)
p_4v(4) = SHMS_4v(4)
else
cth_prot = cos(th_spec_p)
sth_prot = sin(-th_spec_p)
cth_elec = cos(th_spec_e)
sth_elec = sin(th_spec_e)
e_vz = SHMS_vz * 1.
p_vz = HMS_vz * 1.
e_4v(1) = SHMS_4v(1)
p_4v(1) = HMS_4v(1)
e_4v(2) = SHMS_4v(2)
p_4v(2) = HMS_4v(2)
e_4v(3) = SHMS_4v(3)
p_4v(3) = HMS_4v(3)
e_4v(4) = SHMS_4v(4)
p_4v(4) = HMS_4v(4)
endif
if(debug(5)) then
write(*,*) ' '
write(*,*) ' NEW EVENT: ',count, weight
write(*,*) ' HMS: ',e_4v(1),e_4v(2),e_4v(3),e_4v(4)
write(*,*) ' SHMS: ',p_4v(1),p_4v(2),p_4v(3),p_4v(4)
write(*,*) ' e: ',e_4v(1),e_4v(2),e_4v(3),e_4v(4)
write(*,*) ' p: ',p_4v(1),p_4v(2),p_4v(3),p_4v(4)
write(*,*) ' vertex: ',e_vz,p_vz
endif !debug
c
c Rotate 4-vectors to their respective frames
c This is a rotation around the x-axis, which points downwards
c --> the HMS is at negative angles, and the SHMS at positive angles
cth_hms = cos(th_spec_e)
sth_hms = sin(-th_spec_e)
cth_shms = cos(th_spec_p)
sth_shms = sin(th_spec_p)
c Rotatation about the x-axis --> only y, and z change
cc Rotatation about the x-axis --> only y, and z change
e_4vr(1) = e_4v(1)
e_4vr(2) = cth_hms * e_4v(2) - sth_hms * e_4v(3)
e_4vr(3) = sth_hms * e_4v(2) + cth_hms * e_4v(3)
e_4vr(2) = cth_elec * e_4v(2) - sth_elec * e_4v(3)
e_4vr(3) = sth_elec * e_4v(2) + cth_elec * e_4v(3)
e_4vr(4) = e_4v(4)
p_4vr(1) = p_4v(1)
p_4vr(2) = cth_shms * p_4v(2) - sth_shms * p_4v(3)
p_4vr(3) = sth_shms * p_4v(2) + cth_shms * p_4v(3)
p_4vr(2) = cth_prot * p_4v(2) - sth_prot * p_4v(3)
p_4vr(3) = sth_prot * p_4v(2) + cth_prot * p_4v(3)
p_4vr(4) = p_4v(4)
if(debug(5)) then
write(*,*) ' '
write(*,*) ' SPEC:'
write(*,*) ' HMS: ',-th_spec_e/pi*180
write(*,*) ' SHMS: ',th_spec_p/pi*180
write(*,*) ' e: ',th_spec_e/pi*180
write(*,*) ' p: ',th_spec_p/pi*180
write(*,*) ' ROTATED: '
write(*,*) ' HMS: ',e_4vr(1),e_4vr(2),e_4vr(3),e_4vr(4)
write(*,*) ' SHMS: ',p_4vr(1),p_4vr(2),p_4vr(3),p_4vr(4)
write(*,*) ' e: ',e_4vr(1),e_4vr(2),e_4vr(3),e_4vr(4)
write(*,*) ' p: ',p_4vr(1),p_4vr(2),p_4vr(3),p_4vr(4)
endif !debug
c Calculate dxdz and dydz
e_mom = sqrt(e_4vr(1)**2 + e_4vr(2)**2 + e_4vr(3)**2)
......@@ -96,16 +129,18 @@ c Convert momenta and energies into MeV
p_E = p_4vr(4) * 1000.
c Weight and vertex positions are read in correctly.
c Vertez z is read in cm
e_vz = e_vz * 1.
p_vz = p_vz * 1.
weight = w
if(debug(5)) then
write(*,*) ' '
write(*,*) 'SIMC INPUT: '
write(*,*) ' HMS: ',dxdz,dydz,e_mom,e_E,e_vz
write(*,*) ' SHMS: ',dxdzp,dydzp,p_mom,p_E,p_vz
write(*,*) ' e: ',dxdz,dydz,e_mom,e_E,e_vz
write(*,*) ' p: ',dxdzp,dydzp,p_mom,p_E,p_vz
endif !debug
c
if ( abs(atan(e_4v(1)/e_4v(3))*57.3-th_spec_e)<3.
c .and. abs(atan(p_4v(1)/p_4v(3))*57.3-th_spec_p)<3. ) then
count_miss = count_miss+1
endif
return
c
999 write(*,*) ' reached end of file'
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment