1 gaskelld 1.1 subroutine mc_hrsr (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 HRSR spectrometer. (USED TO BE 'HADRON ARM').
8 !
9 !
10 ! Author: David Meekins April 2000
11 ! based on code by David Potterveld
12 !
13 ! Modification History:
14 !
15 ! units for this file are percents, cm, and mrads.
16 !
17 C-______________________________________________________________________________
18
19 implicit none
20
21 include 'struct_hrsr.inc'
22 gaskelld 1.1 include 'apertures_hrsr.inc'
23 include 'g_dump_all_events.inc'
24
25 include '../constants.inc'
26 include '../spectrometers.inc'
27
28 C Spectrometer definitions - for double arm monte carlo compatability
29
30 integer*4 spectr
31 parameter (spectr = 3) !hrs-right is spec #3.
32
33 C Collimator (rectangle) dimensions and offsets.
34
35 real*8 h_entr,v_entr !horiz. and vert. 1/2 gap of fixed slit
36 real*8 h_exit,v_exit !horiz. and vert. 1/2 gap of fixed slit
37 real*8 y_off,x_off !horiz. and vert. position offset of slit
38 real*8 z_off !offset in distance from target to front of sli+
39
40 ! Option for mock sieve slit. Just take particles and forces trajectory
41 ! to put the event in the nearest hole. Must choose the "No collimator"
42 ! values for the h(v)_entr and h(v)_exit values.
43 gaskelld 1.1 ! Note that this will mess up the physics distributions somewhat, but it
44 ! should still be pretty good for optics. Physics limits (e.g. elastic
45 ! peak at x<=1) will not be preserved.
46
47 logical use_sieve /.false./ !use a fake sieve slit.
48
49 ! No collimator - wide open
50 ! parameter (h_entr = 99.)
51 ! parameter (v_entr = 99.)
52 ! parameter (h_exit = 99.)
53 ! parameter (v_exit = 99.)
54
55 ! Large collimator: (hallaweb.jlab.org/news/minutes/collimator-distance.html)
56 parameter (h_entr = 3.145)
57 parameter (v_entr = 6.090)
58 parameter (h_exit = 3.340) !0.1mm wider than 'electron arm'
59 parameter (v_exit = 6.485)
60
61 parameter (y_off = 0.0)
62 parameter (x_off = 0.0)
63 parameter (z_off = 0.0)
64 gaskelld 1.1
65 ! z-position of important apertures.
66 real*8 z_entr,z_exit
67 parameter (z_entr = 110.0d0 + z_off) !nominally 1.100 m
68 parameter (z_exit = z_entr + 8.0d0) !8.0 cm thick
69
70 C Math constants
71
72 real*8 d_r,r_d
73 parameter (d_r = pi/180.)
74 parameter (r_d = 180./pi)
75
76 C The arguments
77
78 real*8 x,y,z !(cm)
79 real*8 dpp !delta p/p (%)
80 real*8 dxdz,dydz !X,Y slope in spectrometer
81 real*8 x_fp,y_fp,dx_fp,dy_fp !Focal plane values to return
82 real*8 p_spec, th_spec !spectrometer setting
83 real*8 fry !vertical position@tgt (+y=down)
84 real*8 pathlen
85 gaskelld 1.1 logical ms_flag !mult. scattering flag
86 logical wcs_flag !wire chamber smearing flag
87 logical decay_flag !check for particle decay
88 logical ok_spec !true if particle makes it
89
90 C Local declarations.
91
92 integer*4 chan/1/,n_classes
93
94 logical*4 first_time_here/.true./
95
96 real*8 dpp_recon,dth_recon,dph_recon !reconstructed quantities
97 real*8 y_recon
98 real*8 p,m2 !More kinematic variables.
99 real*8 xt,yt,rt,tht !temporaries
100 real*8 resmult !DC resolution factor (unused)
101 real*8 zdrift,ztmp
102
103 logical dflag !has particle decayed?
104 logical ok
105
106 gaskelld 1.1 real*8 grnd
107
108 save !Remember it all!
109
110 C ================================ Executable Code =============================
111
112 ! Initialize ok_spec to .flase., restet decay flag
113
114 ok_spec = .false.
115 dflag = .false. !particle has not decayed yet
116 rSTOP_trials = rSTOP_trials + 1
117 xt = th_spec !avoid 'unused variable' error for th_spec
118
119 ! Force particles to go through the sieve slit holes, for mock sieve option.
120 ! Use z_exit, since sieve slit is ~8cm behind normal collimator.
121
122 if (use_sieve) then
123 xt = x + z_exit * dxdz !project to collimator
124 yt = y + z_exit * dydz
125 xt = 2.50*nint(xt/2.50) !shift to nearest hole.
126 yt = 1.25*nint(yt/1.25)
127 gaskelld 1.1 rt = 0.1*sqrt(grnd()) !distance from center of hole(r=1.0mm)
128 tht= 2*pi*grnd() !angle of offset.
129 dxdz = (xt-x)/z_exit !force to correct angle.
130 dydz = (yt-y)/z_exit
131 endif
132
133 ! Save spectrometer coordinates.
134
135 xs = x
136 ys = y
137 zs = z
138 dxdzs = dxdz
139 dydzs = dydz
140
141 ! particle momentum
142
143 dpps = dpp
144 p = p_spec*(1.+dpps/100.)
145
146 C Read in transport coefficients.
147
148 gaskelld 1.1 if (first_time_here) then
149 call transp_init(spectr,n_classes)
150 close (unit=chan)
151 if (n_classes.ne.12) stop 'MC_HRSR, wrong number of transport classes'
152 first_time_here = .false.
153 endif
154
155 ! Begin transporting particle.
156 ! Do transformations, checking against apertures.
157
158 ! Circular apertures before slitbox (only important for no collimator)
159 zdrift = 65.686
160 ztmp = zdrift
161 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen)
162 if (sqrt(xs*xs+ys*ys).gt.7.3787) then
163 rSTOP_slit_hor = rSTOP_slit_hor + 1
164 stop_where=20.
165 x_stop=xs
166 y_stop=ys
167 goto 500
168 endif
169 gaskelld 1.1
170 zdrift = 80.436 - ztmp
171 ztmp = 80.436
172 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen)
173 if (sqrt(xs*xs+ys*ys).gt.7.4092) then
174 rSTOP_slit_hor = rSTOP_slit_hor + 1
175 stop_where=21.
176 x_stop=xs
177 y_stop=ys
178 goto 500
179 endif
180
181
182 ! Check front of fixed slit.
183
184 zdrift = z_entr - ztmp
185 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen)
186 if (abs(ys-y_off).gt.h_entr) then
187 rSTOP_slit_hor = rSTOP_slit_hor + 1
188 stop_where=1.
189 x_stop=xs
190 gaskelld 1.1 y_stop=ys
191 goto 500
192 endif
193 if (abs(xs-x_off).gt.v_entr) then
194 rSTOP_slit_vert = rSTOP_slit_vert + 1
195 stop_where=2.
196 x_stop=xs
197 y_stop=ys
198 goto 500
199 endif
200
201 ! Check back of fixed slit.
202
203 zdrift = z_exit - z_entr
204 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen)
205 if (abs(ys-y_off).gt.h_exit) then
206 rSTOP_slit_hor = rSTOP_slit_hor + 1
207 stop_where=3.
208 x_stop=xs
209 y_stop=ys
210 goto 500
211 gaskelld 1.1 endif
212 if (abs(xs-x_off).gt.v_exit) then
213 rSTOP_slit_vert = rSTOP_slit_vert + 1
214 stop_where=4.
215 x_stop=xs
216 y_stop=ys
217 goto 500
218 endif
219
220 ! Aperture before Q1 (can only check this if next transformation is DRIFT).
221
222 ztmp = 135.064
223 zdrift = ztmp - z_exit
224 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
225 if (sqrt(xs*xs+ys*ys).gt.12.5222) then
226 rSTOP_Q1_in = rSTOP_Q1_in + 1
227 stop_where=22.
228 x_stop=xs
229 y_stop=ys
230 goto 500
231 endif
232 gaskelld 1.1
233 ! Go to Q1 IN mag bound.
234
235 if (.not.adrift(spectr,1)) write(6,*) 'Transformation #1 is NOT a drift'
236 zdrift = driftdist(spectr,1) - ztmp
237 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
238 if ((xs*xs + ys*ys).gt.r_Q1*r_Q1) then
239 rSTOP_Q1_in = rSTOP_Q1_in + 1
240 stop_where=5.
241 x_stop=xs
242 y_stop=ys
243 goto 500
244 endif
245
246 ! Check aperture at 2/3 of Q1.
247
248 call transp(spectr,2,decay_flag,dflag,m2,p,62.75333333d0,pathlen)
249 if ((xs*xs + ys*ys).gt.r_Q1*r_Q1) then
250 rSTOP_Q1_mid = rSTOP_Q1_mid + 1
251 stop_where=6.
252 x_stop=xs
253 gaskelld 1.1 y_stop=ys
254 goto 500
255 endif
256
257 ! Go to Q1 OUT mag boundary.
258
259 call transp(spectr,3,decay_flag,dflag,m2,p,31.37666667d0,pathlen)
260 if ((xs*xs + ys*ys).gt.r_Q1*r_Q1) then
261 rSTOP_Q1_out = rSTOP_Q1_out + 1
262 stop_where=7.
263 x_stop=xs
264 y_stop=ys
265 goto 500
266 endif
267
268 ! Apertures after Q1, before Q2 (can only check this if next trans. is DRIFT).
269
270 zdrift = 300.464 - 253.16 !Q1 exit is z=253.16
271 ztmp = zdrift !distance from Q1 exit
272 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
273 if (sqrt(xs*xs+ys*ys).gt.14.9225) then
274 gaskelld 1.1 rSTOP_Q1_out = rSTOP_Q1_out + 1
275 stop_where=23.
276 x_stop=xs
277 y_stop=ys
278 goto 500
279 endif
280
281 zdrift = 314.464 - 300.464
282 ztmp = ztmp + zdrift !distance from Q1 exit.
283 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
284 if (sqrt(xs*xs+ys*ys).gt.20.9550) then
285 rSTOP_Q2_in = rSTOP_Q2_in + 1
286 stop_where=24.
287 x_stop=xs
288 y_stop=ys
289 goto 500
290 endif
291
292 ! Go to Q2 IN mag bound.
293
294 if (.not.adrift(spectr,4)) write(6,*) 'Transformation #4 is NOT a drift'
295 gaskelld 1.1 zdrift = driftdist(spectr,4) - ztmp
296 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
297 if ((xs*xs + ys*ys).gt.r_Q2*r_Q2) then
298 rSTOP_Q2_in = rSTOP_Q2_in + 1
299 stop_where=8.
300 x_stop=xs
301 y_stop=ys
302 goto 500
303 endif
304
305 ! Check aperture at 2/3 of Q2.
306
307 call transp(spectr,5,decay_flag,dflag,m2,p,121.77333333d0,pathlen)
308 if ((xs*xs + ys*ys).gt.r_Q2*r_Q2) then
309 rSTOP_Q2_mid = rSTOP_Q2_mid + 1
310 stop_where=9.
311 x_stop=xs
312 y_stop=ys
313 goto 500
314 endif
315
316 gaskelld 1.1 ! Go to Q2 OUT mag boundary.
317
318 call transp(spectr,6,decay_flag,dflag,m2,p,60.88666667d0,pathlen)
319 if ((xs*xs + ys*ys).gt.r_Q2*r_Q2) then
320 rSTOP_Q2_out = rSTOP_Q2_out + 1
321 stop_where=10.
322 x_stop=xs
323 y_stop=ys
324 goto 500
325 endif
326
327 ! Apertures after Q2, before D1 (can only check this if next trans. is DRIFT).
328
329 zdrift = 609.664 - 553.020 !Q2 exit is z=553.02
330 ztmp = zdrift !distance from Q2 exit
331 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
332 if (sqrt(xs*xs+ys*ys).gt.30.0073) then
333 rSTOP_Q2_out = rSTOP_Q2_out + 1
334 stop_where=25.
335 x_stop=xs
336 y_stop=ys
337 gaskelld 1.1 goto 500
338 endif
339
340 zdrift = 641.800 - 609.664
341 ztmp = ztmp + zdrift !distance from Q2 exit.
342 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
343 if (sqrt(xs*xs+ys*ys).gt.30.0073) then
344 rSTOP_Q2_out = rSTOP_Q2_out + 1
345 stop_where=26.
346 x_stop=xs
347 y_stop=ys
348 goto 500
349 endif
350
351 zdrift = 819.489 - 641.800
352 ztmp = ztmp + zdrift !distance from Q2 exit.
353 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
354 if (abs(xs).gt.50.0 .or. abs(ys).gt.15.0) then
355 rSTOP_D1_in = rSTOP_D1_in + 1
356 stop_where=27.
357 x_stop=xs
358 gaskelld 1.1 y_stop=ys
359 goto 500
360 endif
361
362 ! Go to D1 IN magnetic boundary, Find intersection with rotated aperture plane.
363
364 if (.not.adrift(spectr,7)) write(6,*) 'Transformation #7 is NOT a drift'
365 zdrift = driftdist(spectr,7) - ztmp
366 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
367 xt=xs
368 yt=ys
369 call rotate_haxis(-30.0,xt,yt)
370 if (abs(xt-2.500).gt.52.5) then ! -50 < x < +55
371 rSTOP_D1_in = rSTOP_D1_in + 1
372 stop_where=11.
373 x_stop=xt
374 y_stop=yt
375 goto 500
376 endif
377 if ( (abs(yt)+0.01861*xt) .gt. 12.5 ) then !tan(1.066) ~ 0.01861
378 rSTOP_D1_in = rSTOP_D1_in +1
379 gaskelld 1.1 stop_where=12.
380 x_stop=xt
381 y_stop=yt
382 goto 500
383 endif
384
385 ! Go to D1 OUT magnetic boundary.
386 ! Find intersection with rotated aperture plane.
387
388 call transp(spectr,8,decay_flag,dflag,m2,p,659.73445725d0,pathlen)
389 xt=xs
390 yt=ys
391 call rotate_haxis(30.0,xt,yt)
392 if (abs(xt-2.500).gt.52.5) then ! -50 < x < +55
393 rSTOP_D1_out = rSTOP_D1_out + 1
394 stop_where=13.
395 x_stop=xt
396 y_stop=yt
397 goto 500
398 endif
399
400 gaskelld 1.1 if ( (abs(yt)+0.01861*xt) .gt. 12.5 ) then !tan(1.066) ~ 0.01861
401 rSTOP_D1_out = rSTOP_D1_out + 1
402 stop_where=14.
403 x_stop=xt
404 y_stop=yt
405 goto 500
406 endif
407
408
409 ! Apertures after D1, before Q3 (can only check this if next trans. is DRIFT).
410
411 zdrift = 1745.33546 - 1655.83446 !D1 exit is z=1655.83446
412 ztmp = zdrift !distance from D1 exit
413 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
414 if (sqrt(xs*xs+ys*ys).gt.30.3276) then
415 rSTOP_D1_out = rSTOP_D1_out + 1
416 stop_where=28.
417 x_stop=xs
418 y_stop=ys
419 goto 500
420 endif
421 gaskelld 1.1 if (abs(xs).gt.50.0 .or. abs(ys).gt.15.0) then
422 rSTOP_D1_out = rSTOP_D1_out + 1
423 stop_where=29.
424 x_stop=xs
425 y_stop=ys
426 goto 500
427 endif
428
429 zdrift = 1759.00946 - 1745.33546
430 ztmp = ztmp + zdrift !distance from D1 exit
431 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
432 if (sqrt(xs*xs+ys*ys).gt.30.3276) then
433 rSTOP_Q3_in = rSTOP_Q3_in + 1
434 stop_where=30.
435 x_stop=xs
436 y_stop=ys
437 goto 500
438 endif
439
440 ! Go to Q3 IN mag bound.
441
442 gaskelld 1.1 if (.not.adrift(spectr,9)) write(6,*) 'Transformation #9 is NOT a drift'
443 zdrift = driftdist(spectr,9) - ztmp
444 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
445 if ((xs*xs + ys*ys).gt.r_Q3*r_Q3) then
446 rSTOP_Q3_in = rSTOP_Q3_in + 1
447 stop_where=15.
448 x_stop=xs
449 y_stop=ys
450 goto 500
451 endif
452
453 ! Check aperture at 2/3 of Q3.
454
455 call transp(spectr,10,decay_flag,dflag,m2,p,121.7866667d0,pathlen)
456 if ((xs*xs + ys*ys).gt.r_Q3*r_Q3) then
457 rSTOP_Q3_mid = rSTOP_Q3_mid + 1
458 stop_where=16.
459 x_stop=xs
460 y_stop=ys
461 goto 500
462 endif
463 gaskelld 1.1
464 ! Go to Q3 OUT mag boundary.
465
466 call transp(spectr,11,decay_flag,dflag,m2,p,60.89333333d0,pathlen)
467 if ((xs*xs + ys*ys).gt.r_Q3*r_Q3) then
468 rSTOP_Q3_out = rSTOP_Q3_out + 1
469 stop_where=17.
470 x_stop=xs
471 y_stop=ys
472 goto 500
473 endif
474
475 ! Apertures after Q3 (can only check this if next trans. is DRIFT).
476
477 zdrift = 2080.38746 - 1997.76446 !Q3 exit is z=1997.76446
478 ztmp = zdrift !distance from Q3 exit
479 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
480 if (abs(xs).gt.35.56 .or. abs(ys).gt.17.145) then
481 rSTOP_Q3_out = rSTOP_Q3_out + 1
482 stop_where=31.
483 x_stop=xs
484 gaskelld 1.1 y_stop=ys
485 goto 500
486 endif
487
488 ! Vacuum window is 15.522cm before FP (which is at VDC1)
489
490 zdrift = 2327.47246 - 2080.38746
491 ztmp = ztmp + zdrift !distance from Q3 exit
492 call project(xs,ys,zdrift,decay_flag,dflag,m2,p,pathlen) !project and decay
493 if (abs(xs).gt.99.76635 .or. abs(ys).gt.17.145) then
494 rSTOP_Q3_out = rSTOP_Q3_out + 1
495 stop_where=32.
496 x_stop=xs
497 y_stop=ys
498 goto 500
499 endif
500
501 ! If we get this far, the particle is in the hut.
502
503 rSTOP_hut = rSTOP_hut + 1
504
505 gaskelld 1.1 if (.not.adrift(spectr,12)) write(6,*) 'Transformation #12 is NOT a drift'
506
507 zdrift = driftdist(spectr,12) - ztmp !distance left to go.
508 call mc_hrsr_hut(m2,p,x_fp,dx_fp,y_fp,dy_fp,ms_flag,wcs_flag,
509 > decay_flag,dflag,resmult,ok,-zdrift,pathlen)
510 if (.not.ok) goto 500
511
512 ! replace xs,ys,... with 'tracked' quantities.
513 xs=x_fp
514 ys=y_fp
515 dxdzs=dx_fp
516 dydzs=dy_fp
517
518 ! Apply offset to y_fp (detectors offset w.r.t optical axis).
519 ! In ESPACE, the offset is taken out for recon, but NOT for y_fp in ntuple,
520 ! so we do not apply it to ys (which goes to recon), but do shift it for y_fp.
521 ! IF THIS IS TRUE OFFSET, WE SHOULD SHIFT DETECTOR APERTURES - NEED TO CHECK!!!!
522 ! But in general the dectectors don't limit the acceptance, so we should be OK.
523
524 y_fp = y_fp - 0.48 !VDC center is at +4.8mm.
525
526 gaskelld 1.1 C Reconstruct target quantities.
527
528 call mc_hrsr_recon(dpp_recon,dth_recon,dph_recon,y_recon,fry)
529 if (.not.ok) then
530 write(6,*) 'mc_hrsr_recon returned ok=',ok
531 goto 500
532 endif
533
534 C Fill output to return to main code
535 dpp = dpp_recon
536 dxdz = dph_recon
537 dydz = dth_recon
538 y = y_recon
539
540 ok_spec = .true.
541 rSTOP_successes = rSTOP_successes + 1
542
543 C We are done with this event, whether GOOD or BAD.
544
545 500 continue
546
547 gaskelld 1.1 C ALL done!
548
549 return
550 end
|
Return to
mc_hrsr.f
CVS log
Up to [HallC] / simc_gfortran / hrsr