1 gaskelld 1.1 subroutine mc_hms (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, fry, ok_spec, pathlen)
4
5 C+______________________________________________________________________________
6 !
7 ! Monte-Carlo of HMS spectrometer.
8 ! Note that we only pass on real*8 variables to the subroutine.
9 ! This will make life easier for portability.
10 !
11 ! Author: David Potterveld, March-1993
12 !
13 ! Modification History:
14 !
15 ! 11-Aug-1993 (D. Potterveld) Modified to use new transformation scheme in
16 ! which each transformation begins at the pivot.
17 !
18 ! 19-AUG-1993 (D. Potterveld) Modified to use COSY INFINITY transformations.
19 !
20 ! 15-SEP-1997 MODIFY stepping through spectrometer so that all drifts
21 ! use project.f (not transp.f), and so that project.f and
22 gaskelld 1.1 ! transp.f take care of decay. Decay distances all assume
23 ! the pathlength for the CENTRAL RAY.
24 C-______________________________________________________________________________
25
26 implicit none
27
28 include 'apertures_hms.inc'
29 include 'struct_hms.inc'
30
31 include '../constants.inc'
32 include '../spectrometers.inc'
33
34 *G.&I.N. STUFF - for checking dipole exit apertures and vacuum pipe to HMS hut.
35 real*8 x_offset_pipes/2.8/,y_offset_pipes/0.0/
36
37 C Spectrometer definitions - for double arm monte carlo compatability
38
39 integer*4 spectr
40 parameter (spectr = 1) !HMS is spec #1.
41
42 ! Collimator (octagon) dimensions and offsets.
43 gaskelld 1.1
44 real*8 h_entr,v_entr,h_exit,v_exit
45 real*8 x_off,y_off,z_off
46
47 ! Option for mock sieve slit. Just take particles and forces trajectory
48 ! to put the event in the nearest hole. Must choose the "No collimator"
49 ! values for the h(v)_entr and h(v)_exit values.
50 ! Note that this will mess up the physics distributions somewhat, but it
51 ! should still be pretty good for optics. Physics limits (e.g. elastic
52 ! peak at x<=1) will not be preserved.
53
54 logical use_sieve /.false./ !use a fake sieve slit.
55
56 ! No collimator - wide open
57 ! parameter (h_entr = 99.)
58 ! parameter (v_entr = 99.)
59 ! parameter (h_exit = 99.)
60 ! parameter (v_exit = 99.)
61
62 ! Old collimator for HMS-1 tune (called 'large' or 'pion' at the time).
63 ! parameter (h_entr = 3.536)
64 gaskelld 1.1 ! parameter (v_entr = 9.003)
65 ! parameter (h_exit = 3.708)
66 ! parameter (v_exit = 9.444)
67
68 ! New collimator for HMS-100 tune.
69 parameter (h_entr = 4.575)
70 parameter (v_entr = 11.646)
71 parameter (h_exit = 4.759)
72 parameter (v_exit = 12.114)
73
74 c parameter (x_off=+0.496) !+ve is slit DOWN - HMS1 tune
75 c parameter (x_off=+0.126) !HMS-100 (preliminary survey)
76 parameter (x_off=+0.000) !HMS-100 (zeroed before Fpi)
77
78 c parameter (y_off=-0.004) !+ve is slit LEFT (as seen from target)
79 parameter (y_off=+0.028) !HMS-100 (number from gaskell)
80
81 c parameter (z_off=+0.00) !1995 position
82 c parameter (z_off=+1.50) !1996 position
83 parameter (z_off=+40.17) !HMS100 tune (dg 5/27/98)
84
85 gaskelld 1.1 ! z-position of important apertures.
86 real*8 z_entr,z_exit
87 parameter (z_entr = 126.2d0 + z_off) !nominally 1.262 m
88 parameter (z_exit = z_entr + 6.3d0) !6.3 cm thick
89
90 real*8 z_dip1,z_dip2,z_dip3 !post dipole apertures
91 parameter (z_dip1 = 64.77d0) ! end of 26.65 inch pipe.
92 parameter (z_dip2 = z_dip1 + 297.18d0) !~117 inch pipe.
93 parameter (z_dip3 = z_dip2 + 115.57d0) ! 45.5 inch pipe.
94
95 ! Math constants
96
97 real*8 d_r,r_d
98 parameter (d_r = pi/180.)
99 parameter (r_d = 180./pi)
100
101 ! The arguments
102
103 real*8 x,y,z !(cm)
104 real*8 dpp !delta p/p (%)
105 real*8 dxdz,dydz !X,Y slope in spectrometer
106 gaskelld 1.1 real*8 x_fp,y_fp,dx_fp,dy_fp !Focal plane values to return
107 real*8 p_spec,th_spec !spectrometer setting
108 real*8 fry !vertical position@tgt (+y=down)
109 real*8 pathlen
110 logical ms_flag !mult. scattering flag
111 logical wcs_flag !wire chamber smearing flag
112 logical decay_flag !check for particle decay
113 logical ok_spec !true if particle makes it
114
115 ! Local declarations.
116
117 integer*4 chan /1/,n_classes
118
119 logical first_time_here/.true./
120
121 real*8 dpp_recon,dth_recon,dph_recon !reconstructed quantities
122 real*8 y_recon
123 real*8 p,m2 !More kinematic variables.
124 real*8 xt,yt,rt,tht !temporaries
125 real*8 resmult !DC resolution factor
126 real*8 zdrift
127 gaskelld 1.1
128 logical dflag !has particle decayed?
129 logical ok
130
131 real*8 grnd
132
133 ! Gaby's dipole shape stuff
134 logical hit_dipole
135 external hit_dipole
136
137 save !Remember it all!
138
139 C ================================ Executable Code =============================
140
141 ! Initialize ok_spec to .false., reset decay flag
142
143 ok_spec = .false.
144 dflag = .false. !particle has not decayed yet
145 hSTOP_trials = hSTOP_trials + 1
146 xt = th_spec !avoid 'unused variable' error for th_spec
147
148 gaskelld 1.1 ! Force particles to go through the sieve slit holes, for mock sieve option.
149
150 if (use_sieve) then
151 xt = x + z_entr * dxdz !project to collimator
152 yt = y + z_entr * dydz
153 xt = 2.540*nint(xt/2.54) !shift to nearest hole.
154 yt = 1.524*nint(yt/1.524)
155 rt = 0.254*sqrt(grnd()) !distance from center of hole(r=2.54mm)
156 tht= 2*pi*grnd() !angle of offset.
157 xt = xt + rt*cos(tht)
158 yt = yt + rt*sin(tht)
159 dxdz = (xt-x)/z_entr !force to correct angle.
160 dydz = (yt-y)/z_entr
161 endif
162
163 ! Save spectrometer coordinates.
164
165 xs = x
166 ys = y
167 zs = z
168 dxdzs = dxdz
169 gaskelld 1.1 dydzs = dydz
170
171 ! particle momentum
172
173 dpps = dpp
174 p = p_spec*(1.+dpps/100.)
175
176 ! Read in transport coefficients.
177
178 if (first_time_here) then
179 call transp_init(spectr,n_classes)
180 close (unit=chan)
181 if (n_classes.ne.12) stop 'MC_HMS, wrong number of transport classes'
182 first_time_here = .false.
183 endif
184
185
186 C------------------------------------------------------------------------------C
187 C Top of Monte-Carlo loop C
188 C------------------------------------------------------------------------------C
189
190 gaskelld 1.1 ! Begin transporting particle.
191
192 ! Do transformations, checking against apertures.
193 ! Check front of fixed slit.
194 zdrift = z_entr
195 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
196 if (abs(ys-y_off).gt.h_entr) then
197 hSTOP_slit_hor = hSTOP_slit_hor + 1
198 goto 500
199 endif
200 if (abs(xs-x_off).gt.v_entr) then
201 hSTOP_slit_vert = hSTOP_slit_vert + 1
202 goto 500
203 endif
204 if (abs(xs-x_off).gt. (-v_entr/h_entr*abs(ys-y_off)+3*v_entr/2)) then
205 hSTOP_slit_oct = hSTOP_slit_oct + 1
206 goto 500
207 endif
208
209 !Check back of fixed slit.
210
211 gaskelld 1.1 zdrift = z_exit-z_entr
212 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
213 if (abs(ys-y_off).gt.h_exit) then
214 hSTOP_slit_hor = hSTOP_slit_hor + 1
215 goto 500
216 endif
217 if (abs(xs-x_off).gt.v_exit) then
218 hSTOP_slit_vert = hSTOP_slit_vert + 1
219 goto 500
220 endif
221 if (abs(xs-x_off).gt. (-v_exit/h_exit*abs(ys-y_off)+3*v_exit/2)) then
222 hSTOP_slit_oct = hSTOP_slit_oct + 1
223 goto 500
224 endif
225
226 ! Go to Q1 IN mag bound. Drift rather than using COSY matrices
227
228 if (.not.adrift(spectr,1)) write(6,*) 'Transformation #1 is NOT a drift'
229 zdrift = driftdist(spectr,1) - z_exit
230 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
231 if ((xs*xs + ys*ys).gt.r_Q1*r_Q1) then
232 gaskelld 1.1 hSTOP_Q1_in = hSTOP_Q1_in + 1
233 goto 500
234 endif
235
236 ! Check aperture at 2/3 of Q1.
237
238 call transp(spectr,2,decay_flag,dflag,m2,p,125.233d0,pathlen)
239 if ((xs*xs + ys*ys).gt.r_Q1*r_Q1) then
240 hSTOP_Q1_mid = hSTOP_Q1_mid + 1
241 goto 500
242 endif
243
244 ! Go to Q1 OUT mag boundary.
245
246 call transp(spectr,3,decay_flag,dflag,m2,p,62.617d0,pathlen)
247 if ((xs*xs + ys*ys).gt.r_Q1*r_Q1) then
248 hSTOP_Q1_out = hSTOP_Q1_out + 1
249 goto 500
250 endif
251
252 ! Go to Q2 IN mag bound. Drift rather than using COSY matrices
253 gaskelld 1.1
254 if (.not.adrift(spectr,4)) write(6,*) 'Transformation #4 is NOT a drift'
255 zdrift = driftdist(spectr,4)
256 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
257 if ((xs*xs + ys*ys).gt.r_Q2*r_Q2) then
258 hSTOP_Q2_in = hSTOP_Q2_in + 1
259 goto 500
260 endif
261
262 ! Check aperture at 2/3 of Q2.
263
264 call transp(spectr,5,decay_flag,dflag,m2,p,143.90d0,pathlen)
265 if ((xs*xs + ys*ys).gt.r_Q2*r_Q2) then
266 hSTOP_Q2_mid = hSTOP_Q2_mid + 1
267 goto 500
268 endif
269
270 ! Go to Q2 OUT mag boundary.
271
272 call transp(spectr,6,decay_flag,dflag,m2,p,71.95d0,pathlen)
273 if ((xs*xs + ys*ys).gt.r_Q2*r_Q2) then
274 gaskelld 1.1 hSTOP_Q2_out = hSTOP_Q2_out + 1
275 goto 500
276 endif
277
278 ! Go to Q3 IN mag bound. Drift rather than using COSY matrices
279
280 if (.not.adrift(spectr,7)) write(6,*) 'Transformation #7 is NOT a drift'
281 zdrift = driftdist(spectr,7)
282 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
283 if ((xs*xs + ys*ys).gt.r_Q3*r_Q3) then
284 hSTOP_Q3_in = hSTOP_Q3_in + 1
285 goto 500
286 endif
287
288 ! Check aperture at 2/3 of Q3.
289
290 call transp(spectr,8,decay_flag,dflag,m2,p,143.8d0,pathlen)
291 if ((xs*xs + ys*ys).gt.r_Q3*r_Q3) then
292 hSTOP_Q3_mid = hSTOP_Q3_mid + 1
293 goto 500
294 endif
295 gaskelld 1.1
296 ! Go to Q3 OUT mag boundary.
297
298 call transp(spectr,9,decay_flag,dflag,m2,p,71.9d0,pathlen)
299 if ((xs*xs + ys*ys).gt.r_Q3*r_Q3) then
300 hSTOP_Q3_out = hSTOP_Q3_out + 1
301 goto 500
302 endif
303
304 ! Go to D1 IN magnetic boundary, Find intersection with rotated aperture plane.
305 ! Aperture has complicated form (evaluated by function hit_dipole).
306
307 if (.not.adrift(spectr,10)) write(6,*) 'Transformation #10 is NOT a drift'
308 zdrift = driftdist(spectr,10)
309 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
310 xt=xs
311 yt=ys
312 call rotate_haxis(-6.0d0,xt,yt)
313 if (hit_dipole(xt,yt)) then
314 hSTOP_D1_in = hSTOP_D1_in + 1
315 goto 500
316 gaskelld 1.1 endif
317
318 ! Go to D1 OUT magnetic boundary.
319 ! Find intersection with rotated aperture plane.
320
321 call transp(spectr,11,decay_flag,dflag,m2,p,526.053d0,pathlen)
322 xt=xs
323 yt=ys
324 call rotate_haxis(6.0d0,xt,yt)
325 if (hit_dipole(xt,yt)) then
326 hSTOP_D1_out = hSTOP_D1_out + 1
327 goto 500
328 endif
329
330 ! Check a number of apertures in the vacuum pipes following the
331 ! dipole. First the odd piece interfacing with the dipole itself
332
333 if ( (((xt-x_offset_pipes)**2+(yt-y_offset_pipes)**2).gt.30.48**2)
334 > .or. (abs((yt-y_offset_pipes)).gt.20.5232) ) then
335 hSTOP_D1_out = hSTOP_D1_out + 1
336 goto 500
337 gaskelld 1.1 endif
338
339 ! Check the exit of the 26.65 inch pipe
340
341 zdrift = z_dip1
342 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
343 if (((xs-x_offset_pipes)**2+(ys-y_offset_pipes)**2).gt.1145.518)then
344 hSTOP_D1_out = hSTOP_D1_out + 1
345 goto 500
346 endif
347
348 ! check exit of long (117 inch) pipe (entrance is bigger than previous pipe)
349 ! note: Rolf claims its 117.5 but the dravings say more like 116.x
350 ! .. so i put 117 even. Should be a 30.62 diameter pipe
351
352 zdrift = z_dip2 - z_dip1
353 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
354 if (((xs-x_offset_pipes)**2+(ys-y_offset_pipes)**2).gt.1512.2299)then
355 hSTOP_D1_out = hSTOP_D1_out + 1
356 goto 500
357 endif
358 gaskelld 1.1
359 ! lastly check the exit of the last piece of pipe. 45.5 inches long, 30.62 dia.
360
361 zdrift = z_dip3 - z_dip2
362 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
363 if (((xs-x_offset_pipes)**2+(ys-y_offset_pipes)**2).gt.2162.9383)then
364 hSTOP_D1_out = hSTOP_D1_out + 1
365 goto 500
366 endif
367
368 ! Note that we do NOT transport (project) to focal plane. We will do this
369 ! in mc_hms_hut.f so that it can take care of all of the decay, mult. scatt,
370 ! and apertures. Pass the current z position so that mc_hms_hut knows
371 ! where to start. Initial zposition for mc_hms_hut is -147.48 cm so that the
372 ! sum of four drift lengths between pipe and focal plane is 625.0 cm
373 ! (64.77+297.18+115.57+147.48=625)
374
375 ! If we get this far, the particle is in the hut.
376
377 hSTOP_hut = hSTOP_hut + 1
378
379 gaskelld 1.1 ! and track through the detector hut
380
381 if (.not.adrift(spectr,12)) write (6,*) 'Transformation #12 is NOT a drift'
382
383 zdrift = driftdist(spectr,12) - z_dip3 !distance left to go.
384 call mc_hms_hut(m2,p,x_fp,dx_fp,y_fp,dy_fp,ms_flag,wcs_flag,
385 > decay_flag,dflag,resmult,ok,-zdrift,pathlen)
386 if (.not.ok) goto 500
387
388 ! replace xs,ys,... with 'tracked' quantities.
389 xs=x_fp
390 ys=y_fp
391 dxdzs=dx_fp
392 dydzs=dy_fp
393
394 ! Reconstruct target quantities.
395 call mc_hms_recon(dpp_recon,dth_recon,dph_recon,y_recon,fry)
396
397 ! Fill output to return to main code
398 dpp = dpp_recon
399 dxdz = dph_recon
400 gaskelld 1.1 dydz = dth_recon
401 y = y_recon
402
403 ok_spec = .true.
404 hSTOP_successes = hSTOP_successes + 1
405
406 ! We are done with this event, whether GOOD or BAD.
407
408 500 continue
409
410 return
411 end
412
413 *-----------------------------------------------------------------------
414
415 logical function hit_dipole(x,y)
416 implicit none
417
418 * Original version made on 04/01/98 by G.&I.Niculescu
419 * to more accurately model the size/shape of the HMS
420 * dipole ...
421 gaskelld 1.1
422 include 'apertures_hms.inc'
423 real*8 x,y
424 real*8 x_local,y_local
425 logical check1,check2,check3,check4,check5,check6
426 logical miss_dipole
427
428 hit_dipole=.false.
429
430 * Let us observe first the obvious symmetry of the problem
431 * This helps reduce the checks to the first quadrant only...
432
433 x_local=abs(x)
434 y_local=abs(y)
435
436 * Now compare the current position and compare it with the different
437 * apertures..
438
439 check1=((x_local.le.x_d1).and.(y_local.le.y_d1))
440 check2=((x_local.le.x_d2).and.(y_local.le.y_d2))
441 check3=((x_local.le.x_d3).and.(y_local.le.y_d3))
442 gaskelld 1.1 check4=((x_local.le.x_d4).and.(y_local.le.y_d4))
443
444 * now, the fifth check is the rounded corner
445
446 check5=(((x_local-x_d5)**2+(y_local-y_d5)**2).le.r_d5**2)
447
448 * lastly the slanted piece
449
450 check6=((x_local.ge.x_d4).and.(x_local.le.x_d3).and.
451 > ((y_local-a_d6*x_local-b_d6).le.0.0))
452
453 * now, if we OR all the above we should get the inside of the can
454
455 miss_dipole=check1.or.check2.or.check3.or.check4.or.check5.or.check6
456
457 * for whatever reason mc_hms expects us to return the OUTSIDE of the can so...
458
459 hit_dipole = .not.miss_dipole
460
461 return
462 end
|
Return to
mc_hms.f
CVS log
Up to [HallC] / simc_gfortran / hms