1 jones 1.1 subroutine mc_hms_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 HMS detector hut.
7 C Note that we only pass on real*8 variables to the subroutine.
8 C This will make life easier for portability.
9 C
10 C The particle is stepped through the detector (using project), and
11 C multiple scattering is applied for each detector or air gap.
12 C If particle decay in enabled, then project.f also checks for
13 C decay of particle. The particle starts at z=zinit. This
14 C needs to be before the first mult. scattering point (the exit window)
15 C or the decay distance is negative, and things are BAD.
16 C
17 C----------------------------------------------------------------------
18
19 implicit none
20
21 include 'struct_hms.inc'
22 jones 1.1 include '../spectrometers.inc'
23
24 C Math constants
25
26 real*8 pi,d_r,r_d,root
27 parameter (pi = 3.141592654)
28 parameter (d_r = pi/180.)
29 parameter (r_d = 180./pi)
30 parameter (root = 0.707106781) !square root of 1/2
31
32 real*8 grnd,gauss1 !external functions
33
34 C----------------------------------------------------------------------
35 C HMS_MATERIALS
36 C CTP parameter file containing the materials of all the HMS 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
43 jones 1.1 C spectrometer exit window, 15 mil Kevlar, 5 mil Mylar (use 20 mil,X0=53.3)
44 real*8 hfoil_exit_radlen,hfoil_exit_thick
45 parameter (hfoil_exit_radlen = 53.3)
46 parameter (hfoil_exit_thick = 0.020*2.54)
47
48 C spectrometer air gaps
49 real*8 hair_radlen
50 parameter (hair_radlen = 30420.)
51
52 C chamber entrance foil, 1 mil of Mylar
53 real*8 hdc_entr_radlen,hdc_entr_thick
54 parameter (hdc_entr_radlen = 28.7)
55 parameter (hdc_entr_thick = 0.001*2.54)
56
57 C chamber gas, 50/50 ethane/argon
58 real*8 hdc_radlen,hdc_thick
59 parameter (hdc_radlen = 16700.0)
60 parameter (hdc_thick = 1.8)
61
62 C effective wire plane, 25 micron W+Be/Cu gives <t>=pi*(.0025/2)**2
63 real*8 hdc_wire_radlen,hdc_wire_thick
64 jones 1.1 parameter (hdc_wire_radlen = 0.35) !Assuming all Tungsten
65 parameter (hdc_wire_thick = 0.0000049)
66
67 C effective cathode plane, Be/Cu
68 real*8 hdc_cath_radlen,hdc_cath_thick
69 parameter (hdc_cath_radlen = 7.2) !'Ave' of Be/Cu
70 parameter (hdc_cath_thick = 0.000177)
71
72 C chamber exit foil, 1 mil of Mylar
73 real*8 hdc_exit_radlen,hdc_exit_thick
74 parameter (hdc_exit_radlen = 28.7)
75 parameter (hdc_exit_thick = 0.001*2.54)
76
77 C hodoscopes
78 real*8 hscin_radlen
79 parameter (hscin_radlen = 42.4)
80
81 C Cherenkov entrance foil, 40 mil of Al
82 real*8 hcer_entr_radlen,hcer_entr_thick
83 parameter (hcer_entr_radlen = 8.90)
84 parameter (hcer_entr_thick = 0.040*2.54)
85 jones 1.1
86 C Cerenkov gas.
87 real*8 hcer_radlen
88 C 0.5 atm of CO2 for e/pi
89 C parameter (hcer_radlen = 36620.0)
90 C 0.5 atm of Freon for better e/pi
91 parameter (hcer_radlen = 9620.0)
92 C 2 atm of Freon for pi/p
93 C parameter (hcer_radlen = 2405.0)
94
95 C Cherenkov mirror, 75 mu plus 2 cm of Rotacell +
96 real*8 hcer_mir_radlen,hcer_mir_thick
97 parameter (hcer_mir_radlen = 400.0)
98 parameter (hcer_mir_thick = 2.0)
99
100 C Cherenkov exit foil
101 real*8 hcer_exit_radlen,hcer_exit_thick
102 parameter (hcer_exit_radlen = 8.90)
103 parameter (hcer_exit_thick = 0.040*2.54)
104
105 C shower counter
106 jones 1.1 real*8 hcal_radlen
107 parameter (hcal_radlen = 2.50)
108
109 C Wire chamber resolutions (sigma)
110 real*8 hdc_sigma(1:12)/ 0.030,0.030,0.030,0.030,0.030,0.030,
111 > 0.030,0.030,0.030,0.030,0.030,0.030/
112
113 C Wire plane positions, construct hdc_zpos array using these parameters
114 integer*4 hdc_nr_cham,hdc_nr_plan
115 parameter (hdc_nr_cham = 2)
116 parameter (hdc_nr_plan = 6)
117
118 real*8 hdc_1_zpos,hdc_1_left,hdc_1_right,hdc_1_top,hdc_1_bot
119 real*8 hdc_1x_offset,hdc_1y_offset
120 real*8 hdc_2_zpos,hdc_2_left,hdc_2_right,hdc_2_top,hdc_2_bot
121 real*8 hdc_2x_offset,hdc_2y_offset
122 real*8 hdc_del_plane
123
124 parameter (hdc_1_zpos = -51.923)
125 parameter (hdc_2_zpos = 29.299)
126 parameter (hdc_del_plane = hdc_thick + hdc_wire_thick + hdc_cath_thick)
127 jones 1.1 parameter (hdc_1_left = 26.0)
128 parameter (hdc_1_right = -26.0)
129 parameter (hdc_1y_offset = 1.443) ! changed to match hdc.pos.nov95
130 parameter (hdc_1_top = -56.5)
131 parameter (hdc_1_bot = 56.5)
132 parameter (hdc_1x_offset = 1.670)
133 parameter (hdc_2_left = 26.0)
134 parameter (hdc_2_right = -26.0)
135 parameter (hdc_2y_offset = 2.753) !changed as above,on 03/03,DD
136 parameter (hdc_2_top = -56.5)
137 parameter (hdc_2_bot = 56.5)
138 parameter (hdc_2x_offset = 2.758)
139
140 C Scintillator positions and thicknesses
141 real*8 hscin_1x_zpos,hscin_1y_zpos
142 real*8 hscin_1x_thick,hscin_1y_thick
143 real*8 hscin_1x_left,hscin_1x_right,hscin_1x_offset
144 real*8 hscin_1y_top,hscin_1y_bot,hscin_1y_offset
145 real*8 hscin_2x_zpos,hscin_2y_zpos
146 real*8 hscin_2x_thick,hscin_2y_thick
147 real*8 hscin_2x_left,hscin_2x_right,hscin_2x_offset
148 jones 1.1 real*8 hscin_2y_top,hscin_2y_bot,hscin_2y_offset
149
150 parameter (hscin_1x_zpos = 77.830)
151 parameter (hscin_1y_zpos = 97.520)
152 parameter (hscin_2x_zpos = 298.820)
153 parameter (hscin_2y_zpos = 318.510)
154 parameter (hscin_1x_thick = 1.067) !1 cm thick, 6.7% overlap
155 parameter (hscin_1y_thick = 1.067)
156 parameter (hscin_2x_thick = 1.067)
157 parameter (hscin_2y_thick = 1.067)
158 parameter (hscin_1x_left = 37.75)
159 parameter (hscin_1x_right = -37.75)
160 parameter (hscin_1x_offset = -1.55)
161 parameter (hscin_1y_top = -60.25)
162 parameter (hscin_1y_bot = 60.25)
163 parameter (hscin_1y_offset = -0.37)
164 parameter (hscin_2x_left = 37.75)
165 parameter (hscin_2x_right = -37.75)
166 parameter (hscin_2x_offset = -0.63)
167 parameter (hscin_2y_top = -60.25)
168 parameter (hscin_2y_bot = 60.25)
169 jones 1.1 parameter (hscin_2y_offset = -1.46)
170
171 C Cherenkov position
172 real*8 hcer_zentrance,hcer_zmirror,hcer_zexit
173 parameter (hcer_zentrance = 110.000)
174 parameter (hcer_zmirror = 245.000)
175 parameter (hcer_zexit = 265.000)
176
177 C Calorimeter position
178 real*8 hcal_1pr_zpos,hcal_2ta_zpos,hcal_3ta_zpos,hcal_4ta_zpos
179 real*8 hcal_left,hcal_right,hcal_top,hcal_bottom
180 parameter (hcal_1pr_zpos = 338.69)
181 parameter (hcal_2ta_zpos = 349.69)
182 parameter (hcal_3ta_zpos = 360.69)
183 parameter (hcal_4ta_zpos = 371.69)
184 parameter (hcal_left = 35.00) !actual size of calorimeter
185 parameter (hcal_right = -35.00)
186 parameter (hcal_top = -69.66)
187 parameter (hcal_bottom = 60.34)
188
189 C The arguments
190 jones 1.1
191 real*8 p,m2 !momentum and mass of particle
192 real*8 x_fp,y_fp,dx_fp,dy_fp !Focal plane values to return
193 real*8 xcal,ycal !Position of track at calorimeter.
194 real*8 zinit !Initial z-position (Not at F.P.)
195 real*8 pathlen
196 logical ms_flag !mult. scattering flag.
197 logical wcs_flag !wire chamber smearing flag
198 logical decay_flag !check for decay
199 logical ok_hut !true if particle makes it
200
201 C Local declarations.
202
203 integer*4 i,iplane,jchamber,npl_off
204 integer*4 scintrig, scincount
205 parameter (scintrig = 3) !set trigger to 3/4
206
207 logical dflag !has particle decayed?
208
209 real*8 resmult,tmpran,tmplim
210 real*8 tmpran1,tmpran2 !temporary random numbers
211 jones 1.1 real*8 radw,drift
212
213 real*8 nsig_max
214 parameter(nsig_max=99.0d0) !max #/sigma for gaussian ran #s.
215
216 C These have to be real*4 for the CERNLIB lfit routine.
217 real*4 badf !temporaries
218 real*4 xfp4,yfp4,dxfp4,dyfp4 !real*4 versions of fp track.
219 real*4 xdc(12),ydc(12),zdc(12) !positions at d.c. planes
220
221 C ================================ Executable Code =============================
222
223 C Initialize some variables
224 C These come from Doug's examination of events with >6 hits per chamber.
225 C In some fraction of events (~10% ?) there were >6 hits per chamber.
226 C The position resolution for these event had approx. 3 times worse resolution.
227 C There was also some position (or delta) depdendence. The numbers below
228 C are based on examining one set of runs. If you're going to use these
229 C to try and reproduce tails in the resolution, you should probably
230 C see if the fraction with >6 hits, the resolution, and the delta depdendance
231 C are consistant with these numbers.
232 jones 1.1 C Note that if you play with the nominal DC resolutions, then you
233 C may want to reduce the resolution multiplier here.
234
235 ! These values are from Doug's analysis of resolutions for the Kaon run.
236 ! tmpran = grnd()
237 ! tmplim = 0.13+0.008666*abs(dpps) !fraction of events with >6 hits.
238 ! if (tmpran.lt.tmplim) then
239 ! resmult = 3.2 !resolution is 3.2x worse (both DCs)
240 ! else
241 ! resmult = 1.0 + 0.05333*dpps
242 ! endif
243 ! tmpran = grnd()
244
245 ! These are more conservative values (some mix of Kaon and NucPi values).
246 ! Best to check what you see in your experiment.
247 tmpran = grnd()
248 tmplim = 0.15 !fraction of events with >6 hits.
249 if (tmpran.lt.tmplim) then
250 resmult = 2.0
251 else
252 resmult = 1.0
253 jones 1.1 endif
254
255 C Initialize ok_hut
256
257 ok_hut = .false.
258
259 C Initialize the xdc and ydc arrays to zero
260
261 do i=1,12
262 xdc(i) = 0.
263 ydc(i) = 0.
264 enddo
265
266 C Initialize scincount to zero
267
268 scincount = 0
269
270 C------------------------------------------------------------------------------C
271 C Top of loop through hut C
272 C------------------------------------------------------------------------------C
273
274 jones 1.1 C Go to spectrometer exit foil, 25cm before DC1 (Drift forwards from zinit).
275 C As usual, neglect effect of nonzero dydzs and dxdzs on radw.
276
277 drift = (hdc_1_zpos-25.000)-zinit !Need to adjust number if known
278 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen) !drift and decay
279 radw = hfoil_exit_thick/hfoil_exit_radlen
280 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
281
282 C Go to first drift chamber set
283 C For simplicity, perform air mult. scattering (probably negligeable)
284 C before drifting instead of 1/2 way through.
285
286 drift = (hdc_1_zpos - 0.5*hdc_nr_plan*hdc_del_plane) -
287 > (hdc_1_zpos-25.000)
288 radw = drift/hair_radlen
289 if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
290 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen) !drift and decay
291
292 jchamber = 1
293 radw = hdc_entr_thick/hdc_entr_radlen
294 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
295 jones 1.1 npl_off = (jchamber-1)*hdc_nr_plan
296 do iplane = 1,hdc_nr_plan
297 radw = hdc_cath_thick/hdc_cath_radlen
298 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
299 drift = 0.5*hdc_thick
300 radw = drift/hdc_radlen
301 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
302 drift = drift + hdc_cath_thick
303 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen) !drift and decay
304 radw = hdc_wire_thick/hdc_wire_radlen
305 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
306 if(wcs_flag) then
307 tmpran1 = gauss1(nsig_max) !gaussian, truncated at 99 sigma.
308 tmpran2 = gauss1(nsig_max)
309 else
310 tmpran1 = 0.
311 tmpran2 = 0.
312 endif
313 xdc(npl_off+iplane)=xs+hdc_sigma(npl_off+iplane)*tmpran1*resmult
314 ydc(npl_off+iplane)=ys+hdc_sigma(npl_off+iplane)*tmpran2*resmult
315 if (iplane.eq.2 .or. iplane.eq.5) then
316 jones 1.1 xdc(npl_off+iplane) = 0. !y plane, no x information
317 else
318 ydc(npl_off+iplane) = 0. !x-like plane, no y info
319 endif
320 drift = 0.5*hdc_thick
321 radw = drift/hdc_radlen
322 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
323 drift = drift + hdc_wire_thick
324 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen) !drift and decay
325 enddo
326 radw = hdc_exit_thick/hdc_exit_radlen
327 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
328 if ( xs.gt.(hdc_1_bot-hdc_1x_offset) .or.
329 > xs.lt.(hdc_1_top-hdc_1x_offset) .or.
330 > ys.gt.(hdc_1_left-hdc_1y_offset) .or.
331 > ys.lt.(hdc_1_right-hdc_1y_offset) ) then
332 hSTOP.dc1 = hSTOP.dc1 + 1
333 goto 500
334 endif
335 radw = hdc_cath_thick/hdc_cath_radlen
336 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
337 jones 1.1
338 C at last cathode foil of first drift chamber set, drift to next
339
340 drift = (hdc_2_zpos - 0.5*hdc_nr_plan*hdc_del_plane) -
341 > (hdc_1_zpos + 0.5*hdc_nr_plan*hdc_del_plane)
342 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen) !drift and decay
343 radw = drift/hair_radlen
344 if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
345
346 jchamber = 2
347 radw = hdc_entr_thick/hdc_entr_radlen
348 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
349 npl_off = (jchamber-1)*hdc_nr_plan
350 do iplane = 1,hdc_nr_plan
351 radw = hdc_cath_thick/hdc_cath_radlen
352 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
353 drift = 0.5*hdc_thick
354 radw = drift/hdc_radlen
355 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
356 drift = drift + hdc_cath_thick
357 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen) !drift and decay
358 jones 1.1 radw = hdc_wire_thick/hdc_wire_radlen
359 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
360 if(wcs_flag) then
361 tmpran1 = gauss1(nsig_max) !gaussian, truncated at 99 sigma.
362 tmpran2 = gauss1(nsig_max)
363 else
364 tmpran1 = 0.
365 tmpran2 = 0.
366 endif
367 xdc(npl_off+iplane)=xs+hdc_sigma(npl_off+iplane)*tmpran1*resmult
368 ydc(npl_off+iplane)=ys+hdc_sigma(npl_off+iplane)*tmpran2*resmult
369 if (iplane.eq.2 .or. iplane.eq.5) then
370 xdc(npl_off+iplane) = 0. !y plane, no x information
371 else
372 ydc(npl_off+iplane) = 0. !x-like plane, no y info
373 endif
374
375 drift = 0.5*hdc_thick
376 radw = drift/hdc_radlen
377 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
378 drift = drift + hdc_wire_thick
379 jones 1.1 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen) !drift and decay
380 enddo
381 radw = hdc_exit_thick/hdc_exit_radlen
382 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
383 if ( xs.gt.(hdc_2_bot-hdc_2x_offset) .or.
384 > xs.lt.(hdc_2_top-hdc_2x_offset) .or.
385 > ys.gt.(hdc_2_left-hdc_2y_offset) .or.
386 > ys.lt.(hdc_2_right-hdc_2y_offset) ) then
387 hSTOP.dc2 = hSTOP.dc2 + 1
388 goto 500
389 endif
390 radw = hdc_cath_thick/hdc_cath_radlen
391 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
392
393 C fit track to give new focal plane values, use LFIT from GENLIB
394 do jchamber=1,hdc_nr_cham
395 npl_off = (jchamber-1)*hdc_nr_plan
396 do iplane=1,hdc_nr_plan
397 if (jchamber.eq.1) zdc(npl_off+iplane) = hdc_1_zpos +
398 > (iplane-0.5-0.5*hdc_nr_plan)*hdc_del_plane
399 if (jchamber.eq.2) zdc(npl_off+iplane) = hdc_2_zpos +
400 jones 1.1 > (iplane-0.5-0.5*hdc_nr_plan)*hdc_del_plane
401 enddo
402 enddo
403
404 call lfit(zdc,xdc,12,0,dxfp4,xfp4,badf)
405 call lfit(zdc,ydc,12,0,dyfp4,yfp4,badf)
406
407 x_fp = dble(xfp4)
408 y_fp = dble(yfp4)
409 dx_fp = dble(dxfp4)
410 dy_fp = dble(dyfp4)
411
412 C at last cathode foil of second drift chamber set, drift to hodoscopes
413
414 drift = hscin_1x_zpos - hdc_2_zpos - 0.5*hdc_nr_plan*hdc_del_plane
415 radw = drift/hair_radlen
416 if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
417 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen) !drift and decay
418 if (ys.lt.(hscin_1x_left+hscin_1y_offset) .and.
419 > ys.gt.(hscin_1x_right+hscin_1y_offset) .and.
420 > xs.lt.(hscin_1y_bot+hscin_1x_offset) .and.
421 jones 1.1 > xs.gt.(hscin_1y_top+hscin_1x_offset) ) scincount = scincount + 1
422 radw = hscin_1x_thick/hscin_radlen
423 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
424
425 drift = hscin_1y_zpos - hscin_1x_zpos
426 radw = drift/hair_radlen
427 if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
428 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen) !drift and decay
429 if (ys.lt.(hscin_1x_left+hscin_1y_offset) .and.
430 > ys.gt.(hscin_1x_right+hscin_1y_offset) .and.
431 > xs.lt.(hscin_1y_bot+hscin_1x_offset) .and.
432 > xs.gt.(hscin_1y_top+hscin_1x_offset) ) scincount = scincount + 1
433 radw = hscin_1y_thick/hscin_radlen
434 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
435
436 C finished first hodoscope, drift to cherenkov
437
438 drift = hcer_zentrance - hscin_1y_zpos
439 radw = drift/hair_radlen
440 if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
441 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen) !drift and decay
442 jones 1.1
443 radw = hcer_entr_thick/hcer_entr_radlen
444 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
445
446 drift = hcer_zmirror - hcer_zentrance
447 radw = drift/hcer_radlen
448 if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
449 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen) !drift and decay
450
451 radw = hcer_mir_thick/hcer_mir_radlen
452 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
453
454 drift = hcer_zexit - hcer_zmirror
455 radw = hcer_exit_thick/hcer_exit_radlen
456 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
457 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen) !drift and decay
458
459 radw = hcer_exit_thick/hcer_exit_radlen
460 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
461
462 C drift to second hodoscope
463 jones 1.1
464 drift = hscin_2x_zpos - hcer_zexit
465 radw = drift/hair_radlen
466 if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
467 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen) !drift and decay
468 if (ys.lt.(hscin_1x_left+hscin_1y_offset) .and.
469 > ys.gt.(hscin_1x_right+hscin_1y_offset) .and.
470 > xs.lt.(hscin_1y_bot+hscin_1x_offset) .and.
471 > xs.gt.(hscin_1y_top+hscin_1x_offset) ) scincount = scincount + 1
472 radw = hscin_2x_thick/hscin_radlen
473 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
474
475 drift = hscin_2y_zpos - hscin_2x_zpos
476 radw = drift/hair_radlen
477 if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
478 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen) !drift and decay
479 if (ys.lt.(hscin_1x_left+hscin_1y_offset) .and.
480 > ys.gt.(hscin_1x_right+hscin_1y_offset) .and.
481 > xs.lt.(hscin_1y_bot+hscin_1x_offset) .and.
482 > xs.gt.(hscin_1y_top+hscin_1x_offset) ) scincount = scincount + 1
483 radw = hscin_2y_thick/hscin_radlen
484 jones 1.1 if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
485
486 C Test on scintillator trigger (DJG 8/18/98)
487 if( scincount .lt. scintrig ) then
488 hSTOP.scin = hSTOP.scin + 1
489 goto 500
490 endif
491
492 C Drift to calorimeter
493 drift = hcal_4ta_zpos - hscin_2y_zpos
494 radw = drift/hair_radlen
495 if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
496 call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen) !drift and decay
497
498 C Note that even with the standard PID trigger, the calorimeter is NOT
499 C required, since the trigger needs either the cerenkov OR the calorimeter.
500 C If you require the calorimeter in you analysis, then you need to make
501 C sure that the calorimeter is hit here, AND apply a seperate post-tracking
502 C fiducial cut (whatever is required by the engine or in your analysis).
503 C The fiducial cut comes later in the code.
504 C
505 jones 1.1 C This cut simply requires that the event hit the calorimeter. It does
506 C not insure that the event is within the good region of the calorimeter.
507 C The seperate fiducial cut is required to insure that the entire shower
508 C energy is contained in the calorimeter. Or, if you like, you can require
509 C some distance between the track and the edge. 2-3 cm seems to be enough
510 C to get most or all of the energy in the calorimeter
511
512 !c if (ys.gt.(hcal_left-2.0) .or. ys.lt.(hcal_right+2.0) .or.
513 !c > xs.gt.(hcal_bottom-2.0) .or. xs.lt.(hcal_top+2.0)) then
514 ! if (ys.gt.hcal_left .or. ys.lt.hcal_right .or.
515 ! > xs.gt.hcal_bottom .or. xs.lt.hcal_top) then
516 ! hSTOP.cal = hSTOP.cal + 1
517 ! goto 500
518 ! endif
519
520 C If you use a calorimeter cut in your analysis, the engine applied a
521 C a fiducial cut at the calorimeter. This is based on the position of the
522 C TRACK at the calorimeter, not the real position of the event. Go to the
523 C back of the calorimeter since engine uses a FID cut at the back.
524 C The standard fiducial cut is 5 cm from the edges of the block.
525
526 jones 1.1 xcal = x_fp + dx_fp * hcal_4ta_zpos
527 ycal = y_fp + dy_fp * hcal_4ta_zpos
528 ! if (ycal.gt.(hcal_left-5.0) .or. ycal.lt.(hcal_right+5.0) .or.
529 ! > xcal.gt.(hcal_bottom-5.0) .or. xcal.lt.(hcal_top+5.0)) then
530
531 ok_hut = .true.
532
533 C We are done with this event, whether GOOD or BAD.
534
535 500 continue
536
537 return
538 end
|
Return to
mc_hms_hut.f
CVS log
Up to [HallC] / Poltar / hms