1 jones 1.1 	subroutine mc_calo (p_spec, th_spec, dpp, x, y, z, dxdz, dydz,
  2                >		x_fp, dx_fp, y_fp, dy_fp, m2,
  3                >		ms_flag, wcs_flag, decay_flag, resmult, frx,fry, ok_spec, 
  4                >          pathlen,using_tgt_field,zi,ztemp,drift_to_cal)
  5           
  6           C+______________________________________________________________________________
  7           !
  8           ! Monte-Carlo of Ge/Gm calorimeter.
  9           !	Note that we only pass on real*8 variables to the subroutine.
 10           !	This will make life easier for portability.
 11           !
 12           ! Author: David Gaskell
 13           !
 14           ! Modification History:
 15           !
 16           !  30-Sep-2002  Just smear the input variables.	
 17           C-______________________________________________________________________________
 18           
 19           	implicit 	none
 20           
 21           	include '../constants.inc'
 22 jones 1.1         include '../spectrometers.inc'
 23           
 24           ! Spectrometer definitions
 25           
 26           	integer*4 spectr
 27           	parameter (spectr = 5)		!calo is spec #5.
 28           
 29           ! Collimator (octagon) dimensions.
 30           
 31           	real*8 h_entr,v_entr,h_exit,v_exit
 32           
 33           
 34           
 35           ! No collimator - wide open
 36           !	parameter (h_entr = 99.)
 37           !	parameter (v_entr = 99.)
 38           !	parameter (h_exit = 99.)
 39           !	parameter (v_exit = 99.)
 40           
 41           ! 'large' collimator.
 42           c	parameter (h_entr = 7.201)
 43 jones 1.1 c	parameter (v_entr = 4.696)
 44           c	parameter (h_exit = 7.567)
 45           c	parameter (v_exit = 4.935)
 46           
 47           ! No collimator, but use collimatro dimnsions to define calo
 48           	parameter (h_entr = 64.0)
 49           	parameter (v_entr = 109.0)
 50           	parameter (h_exit = 64.0)
 51           	parameter (v_exit = 109.0)
 52           
 53           ! z-position of important apertures.            
 54                   real*8 z_entr,z_exit
 55                   parameter (z_entr = 365.0d0)		!nominally 1.263 m
 56                   parameter (z_exit = z_entr + 6.3d0)	!6.3 cm thick
 57           
 58           ! Math constants
 59           
 60           	real*8 d_r,r_d
 61           	parameter (d_r = pi/180.)
 62           	parameter (r_d = 180./pi)
 63           
 64 jones 1.1 ! The arguments
 65           
 66           	real*8 x,y,z				!(cm)
 67           	real*8 dpp				!delta p/p (%)
 68           	real*8 dxdz,dydz			!X,Y slope in spectrometer
 69           	real*8 x_fp,y_fp,dx_fp,dy_fp		!Focal plane values to return
 70           	real*8 p_spec,th_spec			!spectrometer setting
 71           	real*8 fry,frx				!vertical position@tgt (+y=down)
 72           	real*8 pathlen
 73           	logical ms_flag				!mult. scattering flag
 74           	logical wcs_flag			!wire chamber smearing flag
 75           	logical decay_flag			!check for particle decay
 76           	logical ok_spec				!true if particle makes it
 77           	logical using_tgt_field
 78           
 79           ! Local declarations.
 80           
 81           	integer*4	chan	/1/,n_classes
 82           
 83           	logical	first_time_here/.true./
 84           
 85 jones 1.1 	real*8 dpp_recon,dth_recon,dph_recon	!reconstructed quantities
 86           	real*8 y_recon
 87           	real*8 ps,p,m2				!More kinematic variables.
 88           	real*8 xt,yt,rt,tht			!temporaries
 89           	real*8 resmult				!DC resolution factor
 90           	real*8 dummy
 91           	real*8 zdrift
 92           
 93           	logical dflag			!has particle decayed?
 94           	logical ok
 95           
 96           	real*8 ctheta,stheta                      !cos and sin of tgt angle
 97           
 98           	real*8 grnd
 99           
100           	real*8 tmpwidth				!used for aperture cut.
101           	real*8 gauss1
102           	real*8 ztemp,zi
103           	real*8 delta_z,delta_y
104           
105           	real*8 drift_to_cal
106 jones 1.1 
107           	save		!Remember it all!
108           
109           ! ================================ Executable Code =============================
110           	ctheta = cos(th_spec)	! GAW
111           	stheta = sin(th_spec)	! GAW
112           
113           
114           ! Initialize ok_spec to .false., reset decay flag
115           
116           	ok_spec = .false.
117           	dflag = .false.			!particle has not decayed yet
118           
119           ! Save spectrometer coordinates.
120           
121           	xs = x
122           	ys = y
123           	zs = z
124           	dxdzs = dxdz
125           	dydzs = dydz
126           
127 jones 1.1 ! particle momentum
128           
129           	dpps = dpp
130           	p = p_spec*(1.+dpps/100.)
131           
132           ! Begin transporting particle.
133           
134           ! Do transformations, checking against apertures.
135           ! Check front of fixed slit.  (here slit = calo)
136           
137           	  zdrift = drift_to_cal
138           	  call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen)
139           	  if (abs(ys).gt.h_entr) then
140           	    goto 500
141           	  endif
142           	  if (abs(xs).gt.v_entr) then
143           	    goto 500
144           	  endif
145           
146           ! replace xs,ys,... with 'tracked' quantities.
147           	  x_fp = xs + 0.4 * gauss1(99.0)
148 jones 1.1 	  y_fp = ys + 0.4 * gauss1(99.0)
149           
150           	  dx_fp = (x_fp-fry)/(drift_to_cal-zi*ctheta)
151           	  dy_fp = (y_fp-zi*stheta)/(drift_to_cal-zi*ctheta)
152           
153           	  xs=x_fp
154           	  ys=y_fp
155           	  dxdzs=dx_fp
156           	  dydzs=dy_fp
157           
158           	  ps = p*(1.0 + 0.05/sqrt(p/1000.0)*gauss1(99.0))
159           	  dpps = (ps/p_spec-1.0)*100.0
160           ! Fill output to return to main code
161           	  delta_y = -frx*ctheta - ztemp*stheta
162           	  delta_z = frx*stheta + ztemp*ctheta
163           	  call mc_calo_recon(dpp_recon,dth_recon,dph_recon,y_recon,fry,delta_y,delta_z,drift_to_cal)
164           
165           	  dpp = dpp_recon
166           	  dxdz = dph_recon
167           	  dydz = dth_recon
168           	  y = y_recon
169 jones 1.1 
170                     if (using_tgt_field) then
171           	     ok = .TRUE.
172           c	     write(6,*) 'calling track to target from mc_calo'
173           	     call track_to_tgt(dpp,y,dxdz,dydz,frx,-fry,
174           	1	  -p,sqrt(m2),ctheta,stheta,ztemp,-1,ok)
175                    endif
176           	  
177           	  ok_spec = .true.
178           
179           500	continue
180           
181           	return
182           	end