(file) Return to target.f CVS log (file) (dir) Up to [HallC] / simc_semi

  1 gaskelld 1.1 	subroutine trip_thru_target (narm, zpos, energy, theta, Eloss, radlen,
  2                   &                               mass, typeflag)
  3              
  4              	implicit none
  5              	include 'simulate.inc'
  6              
  7              	integer narm
  8              	integer typeflag   !1=generate eloss, 2=min, 3=max, 4=most probable
  9              	real*8 zpos, energy, mass, theta
 10              	real*8 Eloss, radlen
 11              	real*8 forward_path, side_path
 12              	real*8 s_target, s_Al, s_kevlar, s_air, s_mylar	! distances travelled
 13              	real*8 Eloss_target, Eloss_Al,Eloss_air		! energy losses
 14              	real*8 Eloss_kevlar,Eloss_mylar			! (temporary)
 15              	real*8 z_can,t,atmp,btmp,ctmp,costmp,th_can	!for the pudding-can target.
 16              	logical liquid
 17              
 18              	s_Al = 0.0
 19              	liquid = targ.Z.lt.2.4
 20              
 21              	if (abs(zpos) .gt. (targ.length/2.+1.d-5)) then
 22 gaskelld 1.1 	  write(6,*) 'call to trip_thru_target has |zpos| > targ.length/2.'
 23              	  write(6,*) 'could be numerical error, or could be error in target offset'
 24              	  write(6,*) 'zpos=',zpos,'  targ.length/2.=',targ.length/2.
 25              	endif
 26              
 27              ! Which particle are we interested in?
 28              
 29              	goto (10,20,30) narm
 30              
 31              ! The incoming electron
 32              
 33              10	continue
 34              	s_target = (targ.length/2. + zpos) / abs(cos(targ.angle))
 35              	if (liquid) then			!liquid target
 36              	  if (targ.can .eq. 1) then		!beer can (2.8 mil endcap)
 37              	    s_Al = s_Al + 0.0028*inch_cm
 38              	  else if (targ.can .eq. 2) then	!pudding can (5 mil Al, for now)
 39              	    s_Al = s_Al + 0.0050*inch_cm
 40              	  endif
 41              	endif
 42              
 43 gaskelld 1.1 ! ... compute distance in radiation lengths and energy loss
 44              	radlen = s_target/targ.X0_cm + s_Al/X0_cm_Al
 45              	call enerloss_new(s_target,targ.rho,targ.Z,targ.A,energy,mass,
 46                   &                  typeflag,Eloss_target)
 47              	call enerloss_new(s_Al,rho_Al,Z_Al,A_Al,energy,mass,typeflag,Eloss_Al)
 48              	Eloss = Eloss_target + Eloss_Al
 49              	return
 50              
 51              ! The scattered electron
 52              
 53              ! ... compute distances travelled assuming HMS (SOS)
 54              ! ...... 16.0 (8.0) mil of Al-foil for the target chamber exit foil
 55              ! ......~15 cm air between scattering chamber and spectrometer vacuum
 56              ! ...... 15.0 (5.0) mil of kevlar + 5.0 (3.0) mil mylar for the
 57              ! ......	spectrometer entrance foil
 58              ! ...... ASSUMES liquid targets are 2.65" wide, and have 5.0 mil Al side walls.
 59              
 60              ! ... compute distances travelled assuming HRS (E. Schulte thesis)
 61              ! ...... 13.0 mil of Al-foil for the target chamber exit foil
 62              ! ...... WILD GUESS:  ~15 cm air between scattering chamber and HRS vacuum
 63              ! ...... 10.0 mil of Kapton for spectrometer entrance (Use mylar, since 
 64 gaskelld 1.1 ! .....		X0=28.6cm for Kapton, X0=28.7cm for Mylar) 
 65              ! ...... ASSUMES liquid targets are 2.65" wide, and have 5.0 mil Al side walls.
 66              
 67              ! ... For SHMS, use SOS windows for now (just a space filler for now).
 68              
 69              20	continue
 70              	if (electron_arm.eq.1) then		!electron is in HMS
 71              	  s_Al = 0.016*inch_cm
 72              	  s_air = 15
 73               	  s_kevlar = 0.015*inch_cm
 74              	  s_mylar = 0.005*inch_cm
 75              	  forward_path = (targ.length/2.-zpos) / abs(cos(theta+targ.angle))
 76              	else if (electron_arm.eq.2) then			!SOS
 77              	  s_Al = 0.008*inch_cm
 78              	  s_air = 15
 79               	  s_kevlar = 0.005*inch_cm
 80              	  s_mylar = 0.003*inch_cm
 81              	  forward_path = (targ.length/2.-zpos) / abs(cos(theta-targ.angle))
 82              	else if (electron_arm.eq.3) then			!HRS-R
 83              	  s_Al = 0.013*inch_cm
 84              	  s_air = 15
 85 gaskelld 1.1  	  s_kevlar = 0.*inch_cm
 86              	  s_mylar = 0.010*inch_cm
 87              	  forward_path = (targ.length/2.-zpos) / abs(cos(theta-targ.angle))
 88              	else if (electron_arm.eq.4) then			!HRS-L
 89              	  s_Al = 0.013*inch_cm
 90              	  s_air = 15
 91               	  s_kevlar = 0.*inch_cm
 92              	  s_mylar = 0.010*inch_cm
 93              	  forward_path = (targ.length/2.-zpos) / abs(cos(theta-targ.angle))
 94              	else if (electron_arm.eq.5 .or. electron_arm.eq.6) then	!SHMS
 95              	  s_Al = 0.008*inch_cm
 96              	  s_air = 15
 97               	  s_kevlar = 0.005*inch_cm
 98              	  s_mylar = 0.003*inch_cm
 99              	  forward_path = (targ.length/2.-zpos) / abs(cos(theta-targ.angle))
100              	endif
101              	s_target = forward_path
102              
103              	if (liquid) then
104              	  if (targ.can .eq. 1) then		!beer can
105              	    side_path = 1.325*inch_cm / abs(sin(theta))
106 gaskelld 1.1 	    if (forward_path.lt.side_path) then
107              	      s_Al = s_Al + 0.005*inch_cm / abs(cos(theta))
108              	    else
109              	      s_target = side_path
110              	      s_Al = s_Al + 0.005*inch_cm / abs(sin(theta))
111              	    endif
112              	  else if (targ.can .eq. 2) then	!pudding can (5 mil Al, for now)
113              
114              ! this is ugly.  Solve for z position where particle intersects can.  The
115              ! pathlength is then (z_intersect - z_scatter)/cos(theta)
116              ! Angle from center to z_intersect is acos(z_intersect/R).  Therefore the
117              ! angle between the particle and wall is pi/2 - (theta - theta_intersect)
118              
119              	    t=tan(theta)**2
120              	    atmp=1+t
121              	    btmp=-2*zpos*t
122              	    ctmp=zpos**2*t-(targ.length/2.)**2
123              	    z_can=(-btmp+sqrt(btmp**2-4.*atmp*ctmp))/2./atmp
124              	    side_path = (z_can - zpos)/abs(cos(theta))
125                          s_target = side_path
126              	    costmp=z_can/(targ.length/2.)
127 gaskelld 1.1 	    if (abs(costmp).le.1) then
128              	      th_can=acos(z_can/(targ.length/2.))
129              	    else if (abs(costmp-1.).le.0.000001) then
130              	      th_can=0.   !extreme_trip_thru_target can give z/R SLIGHTLY>1.0
131              	    else
132              	      stop 'z_can > can radius in target.f !!!'
133              	    endif
134              	    s_Al = s_Al + 0.0050*inch_cm/abs(sin(target_pi/2 - (theta - th_can)))
135              	  endif
136              
137              	endif		
138              
139              ! ... compute distance in radiation lengths and energy loss
140              	radlen = s_target/targ.X0_cm + s_Al/X0_cm_Al + s_air/X0_cm_air +
141                   >		s_kevlar/X0_cm_kevlar + s_mylar/X0_cm_mylar
142              	call enerloss_new(s_target,targ.rho,targ.Z,targ.A,energy,mass,
143                   &                    typeflag,Eloss_target)
144              	call enerloss_new(s_Al,rho_Al,Z_Al,A_Al,energy,mass,typeflag,Eloss_Al)
145              	call enerloss_new(s_air,rho_air,Z_air,A_air,energy,mass,typeflag,
146                   &                    Eloss_air)
147              	call enerloss_new(s_kevlar,rho_kevlar,Z_kevlar,A_kevlar,energy,
148 gaskelld 1.1      &                    mass,typeflag,Eloss_kevlar)
149              	call enerloss_new(s_mylar,rho_mylar,Z_mylar,A_mylar,energy,mass,
150                   &                    typeflag,Eloss_mylar)
151              	Eloss = Eloss_target + Eloss_Al + Eloss_air + Eloss_kevlar + Eloss_mylar
152              
153              	return
154              
155              ! The scattered proton
156              
157              30	continue
158              	if (hadron_arm.eq.1) then		!proton in HMS
159              	  s_Al = 0.016*inch_cm
160              	  s_air = 15
161               	  s_kevlar = 0.015*inch_cm
162              	  s_mylar = 0.005*inch_cm
163              	  forward_path = (targ.length/2.-zpos) / abs(cos(theta+targ.angle))
164              	else if (hadron_arm.eq.2) then				!SOS
165              	  s_Al = 0.008*inch_cm
166              	  s_air = 15
167               	  s_kevlar = 0.005*inch_cm
168              	  s_mylar = 0.003*inch_cm
169 gaskelld 1.1 	  forward_path = (targ.length/2.-zpos) / abs(cos(theta-targ.angle))
170              	else if (hadron_arm.eq.3) then				!HRS-R
171              	  s_Al = 0.013*inch_cm
172              	  s_air = 15
173               	  s_kevlar = 0.*inch_cm
174              	  s_mylar = 0.010*inch_cm
175              	  forward_path = (targ.length/2.-zpos) / abs(cos(theta-targ.angle))
176              	else if (hadron_arm.eq.4) then				!HRS-L
177              
178              !	  if (abs(zpos/targ.length).lt.0.00001) write(6,*) 'SOMEONE LEFT THE SAFETY WINDOW IN SIMC!!!!!!!'
179              !	  s_Al = 0.125*inch_cm
180              
181              	  s_Al = 0.013*inch_cm
182              	  s_air = 15
183               	  s_kevlar = 0.*inch_cm
184              	  s_mylar = 0.010*inch_cm
185              	  forward_path = (targ.length/2.-zpos) / abs(cos(theta-targ.angle))
186              	else if (hadron_arm.eq.5 .or. hadron_arm.eq.6) then	!SHMS
187              	  s_Al = 0.008*inch_cm
188              	  s_air = 15
189               	  s_kevlar = 0.005*inch_cm
190 gaskelld 1.1 	  s_mylar = 0.003*inch_cm
191              	  forward_path = (targ.length/2.-zpos) / abs(cos(theta-targ.angle))
192              	endif
193              
194              	s_target = forward_path
195              	if (liquid) then
196              	  if (targ.can .eq. 1) then		!beer can
197              	    side_path = 1.325*inch_cm/ abs(sin(theta))
198              	    if (forward_path.lt.side_path) then
199              	      s_Al = s_Al + 0.005*inch_cm / abs(cos(theta))
200              	    else
201              	      s_target = side_path
202              	      s_Al = s_Al + 0.005*inch_cm / abs(sin(theta))
203              	    endif
204              	  else if (targ.can .eq. 2) then	!pudding can (5 mil Al, for now)
205              
206              ! this is ugly.  Solve for z position where particle intersects can.  The
207              ! pathlength is then (z_intersect - z_scatter)/cos(theta)
208              ! Angle from center to z_intersect is acos(z_intersect/R).  Therefore the
209              ! angle between the particle and wall is pi/2 - (theta - theta_intersect)
210              
211 gaskelld 1.1 	    t=tan(theta)**2
212              	    atmp=1+t
213              	    btmp=-2*zpos*t
214              	    ctmp=zpos**2*t-(targ.length/2.)**2
215              	    z_can=(-btmp+sqrt(btmp**2-4.*atmp*ctmp))/2./atmp
216              	    side_path = (z_can - zpos)/abs(cos(theta))
217                          s_target = side_path
218              	    costmp=z_can/(targ.length/2.)
219              	    if (abs(costmp).le.1) then
220              	      th_can=acos(z_can/(targ.length/2.))
221              	    else if (abs(costmp-1.).le.0.000001) then
222              	      th_can=0.   !extreme_trip_thru_target can give z/R SLIGHTLY>1.0
223              	    else
224              	      stop 'z_can > can radius in target.f !!!'
225              	    endif
226              	    s_Al = s_Al + 0.0050*inch_cm/abs(sin(target_pi/2 - (theta - th_can)))
227              	  endif
228              
229              	endif
230              
231              ! ... compute energy losses
232 gaskelld 1.1 
233              	radlen = s_target/targ.X0_cm + s_Al/X0_cm_Al + s_air/X0_cm_air +
234                   >		s_kevlar/X0_cm_kevlar + s_mylar/X0_cm_mylar
235              	call enerloss_new(s_target,targ.rho,targ.Z,targ.A,energy,mass,
236                   &                    typeflag,Eloss_target)
237              	call enerloss_new(s_Al,rho_Al,Z_Al,A_Al,energy,mass,typeflag,Eloss_Al)
238              	call enerloss_new(s_air,rho_air,Z_air,A_air,energy,mass,typeflag,
239                   &                    Eloss_air)
240              	call enerloss_new(s_kevlar,rho_kevlar,Z_kevlar,A_kevlar,energy,mass,
241                   &                    typeflag,Eloss_kevlar)
242              	call enerloss_new(s_mylar,rho_mylar,Z_mylar,A_mylar,energy,mass,
243                   &                    typeflag,Eloss_mylar)
244              
245              	Eloss=Eloss_target+Eloss_Al+Eloss_air+Eloss_kevlar+Eloss_mylar
246              
247              	return
248              	end
249              
250              !-----------------------------------------------------------------
251              
252              	subroutine extreme_trip_thru_target(ebeam, the, thp, pe, pp, z, m)
253 gaskelld 1.1 
254              	implicit none
255              	include 'target.inc'
256              	include 'constants.inc'
257              	structure /limits/
258              	    real*8  min, max
259              	end structure
260              
261              	record /limits/ the, thp, pe, pp, z, betap
262              	integer	i
263              	real*8 th_corner, th_corner_min, th_corner_max
264              	real*8 E1, E2, E3, E4, t1, t2, t3, t4
265              	real*8 zz, th1, th2, m
266              	real*8 ebeam, energymin, energymax
267              	logical	liquid
268              
269              	real*8 zero
270              	parameter (zero=0.0d0)	!double precision zero for subroutine calls
271              
272              !Given limiting values for the electron/proton angles, the z-position in the
273              !target, and beta for the proton, determine min and max losses in target (and
274 gaskelld 1.1 !corresponding amounts of material traversed) This procedure requires knowledge
275              !of the shape of the target, of course, that's why I've put it in here!
276              !Bothering to compute this at all, given that it's going to be used for slop
277              !limits, just proves what a glowing example of Obsessive Compulsive Disorder at
278              !work i truly am.
279              
280              	liquid = targ.Z.lt.2.4
281              
282              ! Incoming electron
283              	call trip_thru_target(1, z.max, ebeam, zero, targ.Eloss(1).max,
284                   &                        targ.teff(1).max, Me, 3)
285              	call trip_thru_target(1, z.min, ebeam, zero, targ.Eloss(1).min,
286                   &                        targ.teff(1).min, Me, 2)
287              
288              ! Scattered electron
289              C Above a few MeV, the energy loss increases as a function of electron
290              C energy, so for any energies we're dealing the max(min) energy loss should
291              C correspond to the max(min) spectrometer energy.  Note that this isn't
292              C the perfect range, but it's easier than reproducing the generated limits here
293              
294              	energymin=pe.min
295 gaskelld 1.1 	energymax=pe.max
296              ! ... solid target is infinitely wide
297              	if (.not.liquid) then
298              	  call trip_thru_target(2, z.min, energymax, the.max,
299                   &           targ.Eloss(2).max, targ.teff(2).max, Me, 3)
300              	  call trip_thru_target(2, z.max, energymin, the.min,
301                   &           targ.Eloss(2).min, targ.teff(2).min, Me, 2)
302              
303              ! ... liquid
304              	else if (targ.can .eq. 1) then		!beer can
305              	  if (z.max.ge.targ.length/2.) then
306              	    th_corner_max = target_pi/2.
307              	  else
308              	    th_corner_max = atan(1.25*inch_cm/(targ.length/2.-z.max))
309              	  endif
310              	  th_corner_min = atan(1.25*inch_cm/(targ.length/2.-z.min))
311              
312              ! ... max loss: do we have access to a corner shot?  try a hair on either
313              ! ... side, front or side walls could be thicker (too lazy to check!)
314              	  if (th_corner_min.le.the.max .and. th_corner_max.ge.the.min) then
315              	    th_corner = max(th_corner_min,the.min)
316 gaskelld 1.1 	    zz = targ.length/2. - 1.25*inch_cm/tan(th_corner)
317              	    th1 = th_corner-.0001
318              	    th2 = th_corner+.0001
319              	  else
320              	    zz = z.min
321              	    th1 = the.min
322              	    th2 = the.max
323              	  endif
324              	  call trip_thru_target(2, zz, energymax, th1, E1, t1, Me, 3)
325              	  call trip_thru_target(2, zz, energymax, th2, E2, t2, Me, 3)
326              	  targ.Eloss(2).max = max(E1,E2)
327              	  targ.teff(2).max = max(t1,t2)
328              
329              ! ........ min loss: try all possibilities (in min case, no way
330              ! ........ an intermediate z or th will do the trick).
331              	  targ.Eloss(2).min = 1.d10
332              	  do i = 0, 3
333              	    call trip_thru_target (2, z.min+int(i/2)*(z.max-z.min), energymin,
334                   >			the.min+mod(i,2)*(the.max-the.min), E1, t1, Me, 2)
335              	    if (E1 .lt. targ.Eloss(2).min) then
336              	      targ.Eloss(2).min = E1
337 gaskelld 1.1 	      targ.teff(2).min = t1
338              	    endif
339              	  enddo
340              
341              	else if (targ.can .eq. 2) then		!pudding can
342              
343              ! ... for pudding can, max loss occurs at lowest scattering angle, where
344              ! ... the outgoing particle leaves the can at z=0.  Therefore, take minimum
345              ! ... scattering angle, project back from z=0, x=radius to get z_init.
346              ! ... Minimum z_init is -radius.
347              	  zz = -(targ.length/2.)/tan(the.min)
348              	  zz = max (zz,(-targ.length/2.))
349              	  call trip_thru_target(2, zz, energymax, the.min, targ.Eloss(2).max,
350                   &                          targ.teff(2).max, Me, 3)
351              
352              ! ........ min loss: try all possibilities (in min case, no way
353              ! ........ an intermediate z or th will do the trick).  Actually, this
354              ! ........ may not be true for the pudding can target, but I'm leaving the
355              ! ........ code alone, for now.
356              	  targ.Eloss(2).min = 1.d10
357              	  do i = 0, 3
358 gaskelld 1.1 	    call trip_thru_target (2, z.min+int(i/2)*(z.max-z.min), energymin,
359                   >			the.min+mod(i,2)*(the.max-the.min), E1, t1, Me, 2)
360              	    if (E1 .lt. targ.Eloss(2).min) then
361              	      targ.Eloss(2).min = E1
362              	      targ.teff(2).min = t1
363              	    endif
364              	  enddo
365              	endif
366              
367              ! Scattered proton. As you can see I'm sufficiently lazy to make
368              ! the code work out whether high or low beta
369              ! maximizes or minimizes loss ... always lower beta --> higher
370              ! losses, it seems
371              
372              ! ... compute extreme energy values
373              	energymin = sqrt(pp.min**2 + m**2)
374              	energymax = sqrt(pp.max**2 + m**2)
375              ! ... compute extreme betap values
376              	betap.min = pp.min / sqrt(pp.min**2 + m**2)
377              	betap.max = pp.max / sqrt(pp.max**2 + m**2)
378              ! ... solid target is infinitely wide
379 gaskelld 1.1 	if (.not.liquid) then
380              	  call trip_thru_target (3, z.min, energymin, thp.max, E1, t1, m, 3)
381              	  call trip_thru_target (3, z.min, energymax, thp.max, E2, t2, m, 3)
382              	  targ.Eloss(3).max = max(E1,E2)
383              	  targ.teff(3).max = max(t1,t2)
384              
385              	  call trip_thru_target (3, z.max, energymin, thp.min, E1, t1, m, 2)
386              	  call trip_thru_target (3, z.max, energymax, thp.min, E2, t2, m, 2)
387              	  targ.Eloss(3).min = min(E1,E2)
388              	  targ.teff(3).min = min(t1,t2)
389              
390              ! ... liquid
391              	else if (targ.can .eq. 1) then	!beer can
392              	  if (z.max .ge. targ.length/2.) then
393              	    th_corner_max = target_pi/2.
394              	  else
395              	    th_corner_max = atan(1.25*inch_cm/(targ.length/2.-z.max))
396              	  endif
397              	  th_corner_min = atan(1.25*inch_cm/(targ.length/2.-z.min))
398              ! ........ max loss: do we have access to a corner shot?
399              ! ........ try a hair on either side, front or side walls could be
400 gaskelld 1.1 ! thicker (too lazy to check!)
401              	  if (th_corner_min.le.thp.max .and. th_corner_max.ge.thp.min) then
402              	    th_corner = max(th_corner_min,thp.min)
403              	    zz = targ.length/2. - 1.25*inch_cm/tan(th_corner)
404              	    th1 = th_corner-.0001
405              	    th2 = th_corner+.0001
406              	  else
407              	    zz = z.min
408              	    th1 = thp.min
409              	    th2 = thp.max
410              	  endif
411              	  call trip_thru_target (3, zz, energymin, th1, E1, t1, m, 3)
412              	  call trip_thru_target (3, zz, energymin, th2, E2, t2, m, 3)
413              	  call trip_thru_target (3, zz, energymax, th1, E3, t3, m, 3)
414              	  call trip_thru_target (3, zz, energymax, th2, E4, t4, m, 3)
415              	  targ.Eloss(3).max = max(E1,E2,E3,E4)
416              	  targ.teff(3).max = max(t1,t2,t3,t4)
417              
418              ! ........ min loss: try all possibilities (in min case, no way
419              ! ........ an intermediate z or th will do the trick)
420              	  targ.Eloss(3).min = 1.d10
421 gaskelld 1.1 	  do i = 0, 3
422              	    call trip_thru_target (3, z.min+int(i/2)*(z.max-z.min), energymin,
423                   >		thp.min+mod(i,2)*(thp.max-thp.min), E1, t1, m, 2)
424              	    if (E1 .lt. targ.Eloss(3).min) then
425              	      targ.Eloss(3).min = E1
426              	      targ.teff(3).min = t1
427              	      zz = z.min+int(i/2)*(z.max-z.min)
428              	      th1 = thp.min+mod(i,2)*(thp.max-thp.min)
429              	    endif
430              	  enddo
431              	  call trip_thru_target (3, zz, energymax, th1, E1, t1, m, 2)
432              	  targ.Eloss(3).min = min(targ.Eloss(3).min, E1)
433              	else if (targ.can .eq. 2) then	!pudding can
434              
435              ! ... for pudding can, max loss occurs at lowest scattering angle, where
436              ! ... the outgoing particle leaves the can at z=0.  Therefore, take minimum
437              ! ... scattering angle, project back from z=0, x=radius to get z_init.
438              ! ... Minimum z_init is -radius.
439              	  zz = -(targ.length/2.)/tan(the.min)
440              	  zz = max (zz,(-targ.length/2.))
441              
442 gaskelld 1.1 	  call trip_thru_target (3, zz, energymin, thp.min, E1, t1, m, 3)
443              	  call trip_thru_target (3, zz, energymax, thp.min, E2, t2, m, 3)
444              	  targ.Eloss(3).max = max(E1,E2)
445              	  targ.teff(3).max = max(t1,t2)
446              
447              ! ........ min loss: try all possibilities (in min case, no way
448              ! ........ an intermediate z or th will do the trick).  This may be
449              ! ........ wrong for the pudding can targ.  Check it out later.
450              	  targ.Eloss(3).min = 1.d10
451              	  do i = 0, 3
452              	    call trip_thru_target (3, z.min+int(i/2)*(z.max-z.min), energymin,
453                   >		thp.min+mod(i,2)*(thp.max-thp.min), E1, t1, m, 2)
454              	    if (E1 .lt. targ.Eloss(3).min) then
455              	      targ.Eloss(3).min = E1
456              	      targ.teff(3).min = t1
457              	      zz = z.min+int(i/2)*(z.max-z.min)
458              	      th1 = thp.min+mod(i,2)*(thp.max-thp.min)
459              	    endif
460              	  enddo
461              	  call trip_thru_target (3, zz, energymax, th1, E1, t1, m, 2)
462              	  targ.Eloss(3).min = min(targ.Eloss(3).min, E1)
463 gaskelld 1.1 
464              	endif
465              
466              *JRA*! Extreme multiple scattering
467              *JRA*
468              *JRA*! ........ don't consider multiple scattering of incoming electron
469              *JRA*	targ.musc_max(1) = 0.
470              
471              ! Extreme multiple scattering.  Use nominal beam energy rather than minimum
472              !  (should be close enough)
473              
474              	call extreme_target_musc(ebeam,1.d0,
475                   >		targ.teff(1).max,targ.musc_max(1),targ.musc_nsig_max)
476              	call extreme_target_musc(pe.min,1.d0,
477                   >		targ.teff(2).max,targ.musc_max(2),targ.musc_nsig_max)
478              	call extreme_target_musc(pp.min,betap.min,
479                   >		targ.teff(3).max,targ.musc_max(3),targ.musc_nsig_max)
480              
481              	return
482              	end
483              
484 gaskelld 1.1 !------------------------------------------------------------------
485              
486              	subroutine target_musc(p,beta,teff,dangles)
487              
488              	implicit none
489              
490              	real*8 Es, epsilon, nsig_max
491              	parameter (Es = 13.6)		!MeV
492              	parameter (epsilon = 0.088)
493              	parameter (nsig_max = 3.5)
494              
495              	real*8 p, beta, teff, dangles(2), dangle, r
496              	real*8 theta_sigma
497              	real*8 gauss1
498              
499              	if (p.lt.25.) write(6,*)
500                   >		'Momentum passed to target_musc should be in MeV, but p=',p
501              
502              ! Compute rms value for planar scattering angle distribution, cf. PDB
503              ! Note teff is thickness of material, in radiation lengths.
504              
505 gaskelld 1.1 	theta_sigma = Es/p/beta * sqrt(teff) * (1+epsilon*log10(teff))
506              
507              ! Compute scattering angles in perpendicular planes.
508              ! Generate two Gaussian numbers BELOW nsig_max.
509              
510              	dangles(1) = theta_sigma * gauss1(nsig_max)
511              	dangles(2) = theta_sigma * gauss1(nsig_max)
512              
513              	return
514              
515              ! Return info about extreme multiple scattering
516              
517              	entry extreme_target_musc(p, beta, teff, dangle, r)
518              
519              	theta_sigma = Es/p/beta * sqrt(teff) * (1+epsilon*log10(teff))
520              	dangle = theta_sigma * nsig_max
521              	r = nsig_max
522              	return
523              	end

Analyzer/Replay: Mark Jones, Documents: Stephen Wood
Powered by
ViewCVS 0.9.2-cvsgraph-1.4.0