1 gaskelld 1.1 subroutine mc_hrsr_hut (m2,p,x_fp,dx_fp,y_fp,dy_fp,ms_flag,wcs_flag,
2 > decay_flag,dflag,resmult,ok_hut,zinit,pathlen)
3
4 C----------------------------------------------------------------------
5 C
6 C Monte-Carlo of HRSR detector hut.
7 C
8 C The particle is stepped through the detector (using project), and
9 C multiple scattering is applied for each detector or air gap.
10 C If particle decay in enabled, then project.f also checks for
11 C decay of particle. The particle starts at z=zinit. This
12 C needs to be before the first mult. scattering point (the exit window)
13 C or the decay distance is negative, and things are BAD.
14 C
15 C----------------------------------------------------------------------
16
17 implicit none
18
19 include 'struct_hrsr.inc'
20 include '../spectrometers.inc'
21 include 'g_dump_all_events.inc'
22 gaskelld 1.1
23 C Math constants
24
25 real*8 pi,d_r,r_d
26
27 parameter (pi = 3.141592654)
28 parameter (d_r = pi/180.)
29 parameter (r_d = 180./pi)
30
31 real*8 gauss1 !external functions
32
33 C all parameters, later to take from .parm files
34 C----------------------------------------------------------------------
35 C HRSR_MATERIALS
36 C CTP parameter file containing the materials of all the HRSR detectors.
37 C For all materials AFTER the bend only have to do multiple scattering.
38 C radlen = 1 radiation length (in cm)
39 C thick = thickness in cm
40 C In case a "+" is added to the comment, the thickness is a guess only.
41 C----------------------------------------------------------------------
42 C spectrometer exit window, 0.1 mm of titanium (X0=3.56cm)
43 gaskelld 1.1 real*8 hfoil_exit_radlen,hfoil_exit_thick
44 parameter (hfoil_exit_radlen = 3.56)
45 parameter (hfoil_exit_thick = 0.01)
46
47 C spectrometer air gaps
48 real*8 hair_radlen
49 parameter (hair_radlen = 30420.)
50
51 C I DON'T KNOW WHERE THE VDC MATERIAL/THICKNESS NUMBERS COME FROM. DO YOU?
52
53 C chamber entrance foil, .18 mm of Kapton Foil (or polyimide film)
54 real*8 hdc_entr_radlen,hdc_entr_thick
55 parameter (hdc_entr_radlen = 34.4)
56 parameter (hdc_entr_thick = 0.00018*2.54)
57
58 C chamber gas, 35/65 ethane/argon
59 real*8 hdc_radlen,hdc_thick
60 parameter (hdc_radlen = 16700.0)
61 ! parameter (hdc_thick = 1.8)
62 parameter (hdc_thick = 1.5) !made up number, but want 1/2 thickness
63 ! < 5 cm to avoid negative drifts.
64 gaskelld 1.1
65 C effective wire plane, 25 micron W+Be/Cu gives <t>=pi*(.0025/2)**2
66 real*8 hdc_wire_radlen,hdc_wire_thick
67 parameter (hdc_wire_radlen = 0.35) !Assuming all Tungsten
68 parameter (hdc_wire_thick = 0.0000049)
69
70 C effective cathode plane, Be/Cu
71 real*8 hdc_cath_radlen,hdc_cath_thick
72 parameter (hdc_cath_radlen = 7.2) !'Ave' of Be/Cu
73 parameter (hdc_cath_thick = 0.000177)
74
75 C chamber exit foil, .18 mm of Kapton Foil (or polyimide film)
76 real*8 hdc_exit_radlen,hdc_exit_thick
77 parameter (hdc_exit_radlen = 34.4)
78 parameter (hdc_exit_thick = 0.00018*2.54)
79
80 C hodoscopes
81 real*8 hscin_radlen
82 parameter (hscin_radlen = 42.4)
83
84 C Cherenkov entrance foil, 40 mil of Tedlar Film
85 gaskelld 1.1 real*8 hcer_entr_radlen,hcer_entr_thick
86 parameter (hcer_entr_radlen = 8.90)
87 parameter (hcer_entr_thick = 0.040*2.54)
88
89 C Cherenkov, 0.5 atm of CO2 for better e/pi
90 real*8 hcer_radlen
91 parameter (hcer_radlen = 36620.0)
92
93 C Cherenkov, 2 atm of Freon for pi/p
94 C real*8 hcer_radlen
95 C parameter (hcer_radlen = 2405.0)
96
97 C Cherenkov mirror, 75 mu plus 2 cm of Rotacell +
98 real*8 hcer_mir_radlen,hcer_mir_thick
99 parameter (hcer_mir_radlen = 400.0)
100 parameter (hcer_mir_thick = 2.0)
101
102 C Cherenkov exit foil
103 real*8 hcer_exit_radlen,hcer_exit_thick
104 parameter (hcer_exit_radlen = 8.90)
105 parameter (hcer_exit_thick = 0.040*2.54)
106 gaskelld 1.1
107 C shower counter
108 real*8 hcal_radlen
109 parameter (hcal_radlen = 14.83)
110
111 C Wire chamber resolutions (sigma)
112
113 real*8 hdc_sigma(1:12)/ 0.0225,0.0225,0.0225,0.0225,0.0225,0.0225,
114 > 0.0225,0.0225,0.0225,0.0225,0.0225,0.0225/
115
116 C Wire plane positions, construct hdc_zpos array using these parameters
117
118 integer*4 hdc_nr_cham,hdc_nr_plan
119 parameter (hdc_nr_cham = 2)
120 parameter (hdc_nr_plan = 6)
121
122 real*8 hdc_1_zpos,hdc_1_left,hdc_1_right,hdc_1_top,hdc_1_bot
123 real*8 hdc_1x_offset,hdc_1y_offset
124 real*8 hdc_2_zpos,hdc_2_left,hdc_2_right,hdc_2_top,hdc_2_bot
125 real*8 hdc_2x_offset,hdc_2y_offset
126
127 gaskelld 1.1 real*8 hdc_del_plane
128
129 C Drift chamber 1 is the focal plane, so shift all zpos values by 25cm
130 parameter (hdc_1_zpos = -25.0 + 25.0)
131 parameter (hdc_2_zpos = 25.0 + 25.0)
132 parameter (hdc_del_plane = hdc_thick + hdc_wire_thick + hdc_cath_thick)
133 parameter (hdc_1_left = 14.4)
134 parameter (hdc_1_right = -14.4)
135 parameter (hdc_1y_offset = 0.000)
136 parameter (hdc_1_top = -105.6)
137 parameter (hdc_1_bot = 105.6)
138 parameter (hdc_1x_offset = 0.000)
139 parameter (hdc_2_left = 14.4)
140 parameter (hdc_2_right = -14.4)
141 parameter (hdc_2y_offset = 0.000)
142 parameter (hdc_2_top = -105.6)
143 parameter (hdc_2_bot = 105.6)
144 parameter (hdc_2x_offset = 0.000)
145
146 C Scintillator positions and thicknesses
147
148 gaskelld 1.1 real*8 hscin_1x_zpos
149 real*8 hscin_1x_thick
150 real*8 hscin_1x_left,hscin_1x_right,hscin_1x_offset
151 real*8 hscin_2x_zpos
152 real*8 hscin_2x_thick
153 real*8 hscin_2x_left,hscin_2x_right,hscin_2x_offset
154
155 parameter (hscin_1x_zpos = 95.0 + 25.0)
156 parameter (hscin_2x_zpos = 288.3 + 25.0)
157 parameter (hscin_1x_thick = 0.5*1.067) ! 1.067 for overlap
158 parameter (hscin_2x_thick = 0.5*1.067)
159 parameter (hscin_1x_left = 18.0)
160 parameter (hscin_1x_right = -18.0)
161 parameter (hscin_1x_offset = 0.00) !up-down offset to hscin_1x.
162 parameter (hscin_2x_left = 30.0)
163 parameter (hscin_2x_right = -30.0)
164 parameter (hscin_2x_offset = 0.00)
165
166
167 C Cherenkov position
168
169 gaskelld 1.1 real*8 hcer_zentrance,hcer_zmirror,hcer_zexit
170 parameter (hcer_zentrance = 137.0 + 25.0)
171 parameter (hcer_zmirror = 197.0 + 25.0)
172 parameter (hcer_zexit = 237.0 + 25.0)
173
174 C Calorimeter position
175
176 real*8 hcal_1pr_zpos,hcal_2ta_zpos,hcal_3ta_zpos,hcal_4ta_zpos
177 real*8 hcal_left,hcal_right,hcal_top,hcal_bottom
178 parameter (hcal_1pr_zpos = 302.3 + 25.0)
179 parameter (hcal_2ta_zpos = 337.3 + 25.0)
180 parameter (hcal_3ta_zpos = 372.3 + 25.0)
181 parameter (hcal_4ta_zpos = 407.3 + 25.0)
182 parameter (hcal_left = 7.5)
183 parameter (hcal_right = -7.5)
184 parameter (hcal_top = -7.5)
185 parameter (hcal_bottom = 7.5)
186
187 C The arguments
188
189 real*8 p,m2 !momentum and mass of particle
190 gaskelld 1.1 real*8 xt,yt !temporary variables.
191 real*8 x_fp,y_fp,dx_fp,dy_fp !Focal plane values to return
192 real*8 xcal,ycal !Position of track at calorimeter.
193 real*8 zinit !Initial z-position (Not at F.P.)
194 real*8 pathlen
195 logical ms_flag !mult. scattering flag.
196 logical wcs_flag !wire chamber smearing flag
197 logical decay_flag !check for decay
198 logical ok_hut !true if particle makes it
199
200 C Local declarations.
201
202 integer*4 i,iplane,jchamber,npl_off
203
204 logical dflag !has particle decayed?
205
206 real*8 resmult
207 real*8 tmpran1,tmpran2 !temporary random numbers
208 real*8 radw,drift
209
210 real*8 nsig_max
211 gaskelld 1.1 parameter(nsig_max=99.0d0) !max #/sigma for gaussian ran #s.
212
213 C These have to be real*4 for the CERNLIB lfit routine
214
215 real*4 badf !temporaries
216 real*4 xfp4,yfp4,dxfp4,dyfp4 !real*4 versions of fp track
217 real*4 xdc(12),ydc(12),zdc(12) !positions at d.c. planes
218
219 C ================================ Executable Code =============================
220
221 C Initialize ok_hut to zero
222
223 ok_hut = .false.
224
225 C Initialize the xdc and ydc arrays to zero
226
227 do i=1,12
228 xdc(i) = 0.
229 ydc(i) = 0.
230 enddo
231 resmult = 1.0
232 gaskelld 1.1
233 C------------------------------------------------------------------------------C
234 C Top of loop through hut C
235 C------------------------------------------------------------------------------C
236
237 C Scatter in spectrometer exit foil, which is located at zinit.
238 C As usual, neglect effect of nonzero dydzs and dxdzs on radw.
239
240 radw = hfoil_exit_thick/hfoil_exit_radlen
241 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
242
243 C Go to first drift chamber set
244 C For simplicity, perform air MS (probably negligeable) at before drift
245 C instead of 1/2 way through.
246
247 drift = (hdc_1_zpos - 0.5*hdc_nr_plan*hdc_del_plane) - zinit
248 radw = drift/hair_radlen
249 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
250 if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
251
252 jchamber = 1
253 gaskelld 1.1 radw = hdc_entr_thick/hdc_entr_radlen
254 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
255 npl_off = (jchamber-1)*hdc_nr_plan
256 do iplane = 1,hdc_nr_plan
257 radw = hdc_cath_thick/hdc_cath_radlen
258 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
259 drift = 0.5*hdc_thick
260 radw = drift/hdc_radlen
261 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
262 drift = drift + hdc_cath_thick
263 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
264 radw = hdc_wire_thick/hdc_wire_radlen
265 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
266 if(wcs_flag) then
267 tmpran1 = gauss1(nsig_max)
268 tmpran2 = gauss1(nsig_max)
269 else
270 tmpran1 = 0.
271 tmpran2 = 0.
272 endif
273 xdc(npl_off+iplane) = xs + hdc_sigma(npl_off+iplane)*tmpran1*resmult
274 gaskelld 1.1 ydc(npl_off+iplane) = ys + hdc_sigma(npl_off+iplane)*tmpran2*resmult
275 if (iplane.eq.1 .or. iplane.eq.3 .or. iplane.eq.5) then
276 xdc(npl_off+iplane) = 0. !y plane, no x information
277 else
278 ydc(npl_off+iplane) = 0. !x-like plane, no y info
279 endif
280
281 drift = 0.5*hdc_thick
282 radw = drift/hdc_radlen
283 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
284 drift = drift + hdc_wire_thick
285 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
286 enddo
287 radw = hdc_exit_thick/hdc_exit_radlen
288 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
289
290 !rotate 45 degrees to compare to VDCs. CHECK SIGN AND SIZE OF ROTATAION!!!
291 xt=xs
292 yt=ys
293 call rotate_haxis(45.0d0,xt,yt)
294
295 gaskelld 1.1 if (xt.gt.(hdc_1_bot-hdc_1x_offset) .or.
296 > xt.lt.(hdc_1_top-hdc_1x_offset) .or.
297 > yt.gt.(hdc_1_left-hdc_1y_offset) .or.
298 > yt.lt.(hdc_1_right-hdc_1y_offset) ) then
299 rSTOP_dc1 = rSTOP_dc1 + 1
300 stop_where=17.
301 x_stop=xs
302 y_stop=ys
303 goto 500
304 endif
305 radw = hdc_cath_thick/hdc_cath_radlen
306 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
307
308 C at last cathode foil of first drift chamber set, drift to next
309
310 drift = hdc_2_zpos - hdc_1_zpos - hdc_nr_plan*hdc_del_plane
311 radw = drift/hair_radlen
312 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
313 if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
314
315 jchamber = 2
316 gaskelld 1.1 radw = hdc_entr_thick/hdc_entr_radlen
317 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
318 npl_off = (jchamber-1)*hdc_nr_plan
319 do iplane = 1,hdc_nr_plan
320 radw = hdc_cath_thick/hdc_cath_radlen
321 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
322 drift = 0.5*hdc_thick
323 radw = drift/hdc_radlen
324 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
325 drift = drift + hdc_cath_thick
326 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
327 radw = hdc_wire_thick/hdc_wire_radlen
328 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
329 if(wcs_flag) then
330 tmpran1 = gauss1(nsig_max)
331 tmpran2 = gauss1(nsig_max)
332 else
333 tmpran1 = 0.
334 tmpran2 = 0.
335 endif
336 xdc(npl_off+iplane) = xs + hdc_sigma(npl_off+iplane)*tmpran1*resmult
337 gaskelld 1.1 ydc(npl_off+iplane) = ys + hdc_sigma(npl_off+iplane)*tmpran2*resmult
338 if (iplane.eq.1 .or. iplane.eq.3 .or. iplane.eq.5) then
339 xdc(npl_off+iplane) = 0. !y plane, no x information
340 else
341 ydc(npl_off+iplane) = 0. !x-like plane, no y info
342 endif
343
344 drift = 0.5*hdc_thick
345 radw = drift/hdc_radlen
346 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
347 drift = drift + hdc_wire_thick
348 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
349 enddo
350 radw = hdc_exit_thick/hdc_exit_radlen
351 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
352
353 !rotate 45 degrees to compare to VDCs. CHECK SIGN AND SIZE OF ROTATAION!!!
354 xt=xs
355 yt=ys
356 call rotate_haxis(45.0d0,xt,yt)
357
358 gaskelld 1.1 if (xt.gt.(hdc_2_bot-hdc_2x_offset) .or.
359 > xt.lt.(hdc_2_top-hdc_2x_offset) .or.
360 > yt.gt.(hdc_2_left-hdc_2y_offset) .or.
361 > yt.lt.(hdc_2_right-hdc_2y_offset) ) then
362 rSTOP_dc2 = rSTOP_dc2 + 1
363 stop_where=19.
364 x_stop=xs
365 y_stop=ys
366 goto 500
367 endif
368 radw = hdc_cath_thick/hdc_cath_radlen
369 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
370
371 C fit track to give new focal plane values, use LFIT from GENLIB
372
373 do jchamber=1,hdc_nr_cham
374 npl_off = (jchamber-1)*hdc_nr_plan
375 do iplane=1,hdc_nr_plan
376 if (jchamber.eq.1) zdc(npl_off+iplane) = hdc_1_zpos +
377 > (iplane-0.5-0.5*hdc_nr_plan)*hdc_del_plane
378 if (jchamber.eq.2) zdc(npl_off+iplane) = hdc_2_zpos +
379 gaskelld 1.1 > (iplane-0.5-0.5*hdc_nr_plan)*hdc_del_plane
380 enddo
381 enddo
382
383 call lfit(zdc,xdc,12,0,dxfp4,xfp4,badf)
384 call lfit(zdc,ydc,12,0,dyfp4,yfp4,badf)
385
386 x_fp = dble(xfp4)
387 y_fp = dble(yfp4)
388 dx_fp = dble(dxfp4)
389 dy_fp = dble(dyfp4)
390
391 C at last cathode foil of second drift chamber set, drift to hodoscopes
392
393 drift = hscin_1x_zpos - hdc_2_zpos - 0.5*hdc_nr_plan*hdc_del_plane
394 radw = drift/hair_radlen
395 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
396 if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
397 if (ys.gt.(hscin_1x_left) .or.
398 > ys.lt.(hscin_1x_right)) then
399 rSTOP_s1 = rSTOP_s1 + 1
400 gaskelld 1.1 stop_where=20.
401 x_stop=xs
402 y_stop=ys
403 goto 500
404 endif
405 radw = hscin_1x_thick/hscin_radlen
406 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
407
408 C finished first hodoscope, drift to cherenkov
409
410 drift = hcer_zentrance - hscin_1x_zpos
411 radw = drift/hair_radlen
412 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
413 if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
414
415 radw = hcer_entr_thick/hcer_entr_radlen
416 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
417
418 drift = hcer_zmirror - hcer_zentrance
419 radw = drift/hcer_radlen
420 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
421 gaskelld 1.1 if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
422
423 radw = hcer_mir_thick/hcer_mir_radlen
424 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
425
426 drift = hcer_zexit - hcer_zmirror
427 radw = drift/hcer_radlen
428 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
429 if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
430
431 radw = hcer_exit_thick/hcer_exit_radlen
432 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
433
434 C drift to second hodoscope
435
436 drift = hscin_2x_zpos - hcer_zexit
437 radw = drift/hair_radlen
438 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
439 if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
440 if (ys.gt.(hscin_2x_left) .or.
441 > ys.lt.(hscin_2x_right)) then
442 gaskelld 1.1 rSTOP_s2 = rSTOP_s2 + 1
443 stop_where=21.
444 x_stop=xs
445 y_stop=ys
446 goto 500
447 endif
448 radw = hscin_2x_thick/hscin_radlen
449 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
450
451 C Drift to back of calorimeter. This should be done even if we do not apply
452 C any calorimeter cuts, so that pathlen is calculated to the back of the spec.
453 drift = hcal_4ta_zpos - hscin_2x_zpos
454 radw = drift/hair_radlen
455 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
456 if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
457
458
459 C Don't need to drift to calorimeter unless it's required in your trigger.
460 C Note that even with the standard PID trigger, the calorimeter is NOT
461 C required, since the trigger needs either the cerenkov OR the calorimeter.
462 C There is a seperate fiducial cut needed if you require the calorimeter
463 gaskelld 1.1 C in you analysis. That cut is applied AFTER fitting the track (see below).
464 * if (ys.gt.hcal_left .or. ys.lt.hcal_right .or.
465 * > xs.gt.hcal_bottom .or. xs.lt.hcal_top) then
466 * rSTOP_cal = rSTOP_cal + 1
467 * stop_where=22.
468 * x_stop=xs
469 * y_stop=ys
470 * goto 500
471 * endif
472
473 C If you use a calorimeter cut in your analysis, the engine applied a
474 C a fiducial cut at the calorimeter. This is based on the position of the
475 C TRACK at the calorimeter, not the real position of the event. Go to the
476 C back of the calorimeter since engine uses a FID cut at the back.
477 C The standard fiducial cut is 5 cm from the edges of the block.
478
479 xcal = x_fp + dx_fp * hcal_4ta_zpos
480 ycal = y_fp + dy_fp * hcal_4ta_zpos
481 * if (ycal.gt.(hcal_left-5.0) .or. ycal.lt.(hcal_right+5.0) .or.
482 * > xcal.gt.(hcal_bottom-5.0) .or. xcal.lt.(hcal_top+5.0)) then
483 * rSTOP_cal = rSTOP_cal + 1
484 gaskelld 1.1 * stop_where=23.
485 * x_stop=xs
486 * y_stop=ys
487 * goto 500
488 * endif
489
490 ok_hut = .true.
491
492 C We are done with this event, whether GOOD or BAD.
493
494 500 continue
495
496 C ALL done!
497
498 return
499 end
|
Return to
mc_hrsr_hut.f
CVS log
Up to [HallC] / simc_gfortran / hrsr