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

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

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