1 jones 1.1 subroutine mc_sos_hut (m2,p,x_fp,dx_fp,y_fp,dy_fp,ms_flag,wcs_flag,
2 > decay_flag,dflag,resmult,ok_hut,zinit)
3
4 C----------------------------------------------------------------------
5 C
6 C Monte-Carlo of SOS detector hut.
7 C Note that we only pass on real*4 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 '../ran.inc'
22 jones 1.1 include '../struct_sos.inc'
23 include '../track.inc'
24
25 C Math constants
26
27 real*4 pi,d_r,r_d,root
28 parameter (pi = 3.141592654)
29 parameter (d_r = pi/180.)
30 parameter (r_d = 180./pi)
31 parameter (root = 0.707106781) !square root of 1/2
32
33 C----------------------------------------------------------------------
34 C SOS_MATERIALS
35 C CTP parameter file containing the materials of all the SOS detectors.
36 C For all materials AFTER the bend only have to do multiple scattering.
37 C radlen = 1 radiation length (in cm)
38 C thick = thickness in cm
39 C In case a "+" is added to the comment, the thickness is a guess only.
40 C----------------------------------------------------------------------
41 C spectrometer exit window, 15 mil Kevlar, 5 mil Mylar (use 20 mil,X0=53.3)
42 real*4 sfoil_exit_radlen,sfoil_exit_thick
43 jones 1.1 real*4 sfoil_exit_zpos,sfoil_exit_ang
44 parameter (sfoil_exit_radlen = 53.3)
45 parameter (sfoil_exit_thick = 0.020*2.54)
46 parameter (sfoil_exit_zpos = -3.22)
47 parameter (sfoil_exit_ang = 0.*d_r) !50 deg. BEFORE vacuum extension,
48 ! 0 deg. AFTER installation.
49 C spectrometer air gaps
50 real*4 sair_radlen
51 parameter (sair_radlen = 30420.) !radlen, cm
52
53 C chamber gas
54 real*4 sdc_radlen,sdc_thick
55 parameter (sdc_radlen = 16700.0)
56 parameter (sdc_thick = 0.61775)
57
58 C effective wire plane ! 60 micron field wires+30 micron sense wires, Tungsten
59 real*4 sdc_wire_radlen,sdc_wire_thick
60 parameter (sdc_wire_radlen = 0.35)
61 parameter (sdc_wire_thick = 0.0000354)
62
63 C chamber cathode foils, 1/2 mil of Mylar
64 jones 1.1 real*4 sdc_cath_radlen,sdc_cath_thick
65 parameter (sdc_cath_radlen = 28.7)
66 parameter (sdc_cath_thick = 0.0005*2.54)
67
68 C hodoscopes
69 real*4 sscin_radlen
70 parameter (sscin_radlen = 42.4)
71
72 C Cherenkov entrance foil, 0.5 mm Al +
73 real*4 scer_entr_radlen,scer_entr_thick
74 parameter (scer_entr_radlen = 8.90)
75 parameter (scer_entr_thick = 0.050)
76
77 C Cherenkov, 1 atm Freon 12 +
78 real*4 scer_radlen
79 parameter (scer_radlen = 4810.0)
80
81 C Cherenkov mirror, 75 mu plus 2 cm Rotacell +
82 real*4 scer_mir_radlen,scer_mir_thick
83 parameter (scer_mir_radlen = 400.0)
84 parameter (scer_mir_thick = 2.0)
85 jones 1.1
86 C Cherenkov exit foil, 0.5 mm Al +
87 real*4 scer_exit_radlen,scer_exit_thick
88 parameter (scer_exit_radlen = 8.90)
89 parameter (scer_exit_thick = 0.050)
90
91 C Aerogel entrance foil
92 real*4 saer_entr_radlen,saer_entr_thick
93 parameter (saer_entr_radlen = 8.90)
94 parameter (saer_entr_thick = 0.0625)
95
96 C Aerogel +
97 real*4 saer_radlen,saer_thick
98 parameter (saer_radlen = 150.0)
99 parameter (saer_thick = 9.0)
100
101 C Aerogel exit foil
102 real*4 saer_exit_radlen,saer_exit_thick
103 parameter (saer_exit_radlen = 8.90)
104 parameter (saer_exit_thick = 0.0625)
105
106 jones 1.1 C shower counter
107 real*4 scal_radlen
108 parameter (scal_radlen = 2.50)
109
110 C Wire chamber resolutions (sigma)
111 real*4 sdc_sigma(1:12)/ 0.030,0.030,0.030,0.030,0.030,0.030,
112 > 0.030,0.030,0.030,0.030,0.030,0.030/
113
114 C Wire plane positions, construct sdc_zpos array using these parameters
115 integer*4 sdc_nr_cham,sdc_nr_plan
116 parameter (sdc_nr_cham = 2)
117 parameter (sdc_nr_plan = 6)
118
119 real*4 sdc_1_zpos,sdc_1_left,sdc_1_right,sdc_1_top,sdc_1_bot
120 real*4 sdc_1x_offset,sdc_1y_offset
121 real*4 sdc_2_zpos,sdc_2_left,sdc_2_right,sdc_2_top,sdc_2_bot
122 real*4 sdc_2x_offset,sdc_2y_offset
123 real*4 sdc_del_plane
124 parameter (sdc_1_zpos = 6.25)
125 parameter (sdc_2_zpos = 55.77)
126 parameter (sdc_del_plane = sdc_thick + sdc_wire_thick + sdc_cath_thick)
127 jones 1.1 parameter (sdc_1_left = 24.0)
128 parameter (sdc_1_right = -24.0)
129 c parameter (sdc_1y_offset = -1.769)
130 parameter (sdc_1y_offset = -1.822)
131 parameter (sdc_1_top = -32.0)
132 parameter (sdc_1_bot = 32.0)
133 c parameter (sdc_1x_offset = 8.0226)
134 parameter (sdc_1x_offset = 8.649)
135 parameter (sdc_2_left = 24.0)
136 parameter (sdc_2_right = -24.0)
137 c parameter (sdc_2y_offset = -1.794)
138 parameter (sdc_2y_offset = -1.976)
139 parameter (sdc_2_top = -32.0)
140 parameter (sdc_2_bot = 32.0)
141 c parameter (sdc_2x_offset = -2.231)
142 parameter (sdc_2x_offset = -1.532)
143
144 C Scintillator positions and thicknesses
145 real*4 sscin_1x_zpos,sscin_1y_zpos
146 real*4 sscin_1x_thick,sscin_1y_thick
147 real*4 sscin_1_left,sscin_1_right,sscin_1x_offset
148 jones 1.1 real*4 sscin_1_top,sscin_1_bot,sscin_1y_offset
149 real*4 sscin_2x_zpos,sscin_2y_zpos
150 real*4 sscin_2x_thick,sscin_2y_thick
151 real*4 sscin_2_left,sscin_2_right,sscin_2x_offset
152 real*4 sscin_2_top,sscin_2_bot,sscin_2y_offset
153 parameter (sscin_1y_zpos = 73.61)
154 parameter (sscin_1x_zpos = 97.11)
155 parameter (sscin_2y_zpos = 249.51)
156 parameter (sscin_2x_zpos = 290.81)
157 parameter (sscin_1x_thick = 1.040) !1.040 for overlap
158 parameter (sscin_1y_thick = 1.098)
159 parameter (sscin_2x_thick = 1.040)
160 parameter (sscin_2y_thick = 1.098)
161 parameter (sscin_1_left = 18.25)
162 parameter (sscin_1_right = -18.25)
163 parameter (sscin_1x_offset = 2.8)
164 parameter (sscin_1_top = -31.75)
165 parameter (sscin_1_bot = 31.75)
166 parameter (sscin_1y_offset = 2.25)
167 parameter (sscin_2_left = 18.25)
168 parameter (sscin_2_right = -18.25)
169 jones 1.1 parameter (sscin_2x_offset = 4.9)
170 parameter (sscin_2_top = -56.25)
171 parameter (sscin_2_bot = 56.25)
172 parameter (sscin_2y_offset = 2.9)
173
174 C Cherenkov position
175 real*4 scer_zentrance,scer_zmirror,scer_zexit
176 parameter (scer_zentrance = 130.000)
177 parameter (scer_zmirror = 155.000)
178 parameter (scer_zexit = 160.000)
179
180 C Calorimeter position
181 real*4 scal_1pr_zpos,scal_2ta_zpos,scal_3ta_zpos,scal_4ta_zpos
182 real*4 scal_left,scal_right,scal_top,scal_bottom
183 parameter (scal_1pr_zpos = 313.01)
184 parameter (scal_2ta_zpos = 324.01)
185 parameter (scal_3ta_zpos = 335.01)
186 parameter (scal_4ta_zpos = 346.01)
187 parameter (scal_left = 35.00)
188 parameter (scal_right = -35.00)
189 parameter (scal_top = -61.)
190 jones 1.1 parameter (scal_bottom = 49.)
191
192 C The arguments
193
194 real*4 p,m2 !momentum and mass of particle
195 real*4 x_fp,y_fp,dx_fp,dy_fp !Focal plane values to return
196 real*4 ms_flag, wcs_flag !particle, m_scat, wc_smear
197 real*4 ok_hut !true if particle makes it
198 real*4 xcal,ycal !Position of track at calorimeter.
199 real*4 zinit !Initial z-position (Not at F.P.)
200 logical decay_flag !check for decay
201
202 C Local declarations.
203
204 integer*4 i,iplane,jchamber,npl_off
205 integer*4 sscintrig,sscincount
206 parameter (sscintrig = 3) !set trigger to 3/4
207
208 logical*4 dflag !has particle decayed?
209
210 real*4 resmult,tmpran
211 jones 1.1 real*4 random(4) !space for 4 random numbers
212 real*4 xt,badf !temporaries
213 real*4 radw,drift
214 real*4 xdc(12),ydc(12),zdc(12) !positions at d.c. planes
215 real*4 x0,y0,dx,dy !track fitting temporaries
216
217 C ================================ Executable Code =============================
218
219 C Initialize some variables
220
221 resmult=1.
222 tmpran=ran(iran3)
223 if(tmpran.lt.0.15)resmult=3. !15% of events have >6 hits, 3x resolution
224
225 sscincount = 0
226
227 C Initialize ok_hut to zero
228
229 ok_hut = 0.
230
231 C Initialize the xdc and ydc arrays to zero
232 jones 1.1
233 do i=1,12
234 xdc(i) = 0.
235 ydc(i) = 0.
236 enddo
237
238 C------------------------------------------------------------------------------C
239 C Top of loop through hut C
240 C------------------------------------------------------------------------------C
241
242 C Go to spectrometer exit foil. (Drift forwards from zinit).
243 C Approximation is used in calculating xt = x at foil crossing.
244 C As usual, neglect effect of nonzero dydzs and dxdzs on radw.
245
246 xt = xs + dxdzs * sfoil_exit_zpos
247 drift = (sfoil_exit_zpos + xt * tan(sfoil_exit_ang)) - zinit
248 if (abs(drift).le.0.001) drift=0.001 !avoid drift<=0
249 call project(xs,ys,drift,decay_flag,dflag,m2,p) !drift and decay
250 radw = sfoil_exit_thick/sfoil_exit_radlen/cos(sfoil_exit_ang)
251 if(ms_flag.gt.0.1) call musc(m2,p,radw,dydzs,dxdzs)
252
253 jones 1.1 !Already decayed for the 3.22 cm, so it get double counted in the next
254 !step. Oh well. Fix this later.
255
256
257 C Go to first drift chamber set
258 C For simplicity, perform air MS (probably negligeable) before DCs
259 C instead of 1/2 way through.
260
261 drift = (sdc_1_zpos - 0.5*sdc_nr_plan*sdc_del_plane) -
262 > (sfoil_exit_zpos + xt * tan(sfoil_exit_ang))
263 radw = drift/sair_radlen
264 if(ms_flag.gt.0.1) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
265 call project(xs,ys,drift,decay_flag,dflag,m2,p) !drift and decay
266
267 jchamber = 1
268 npl_off = (jchamber-1)*sdc_nr_plan
269 do iplane = 1,sdc_nr_plan
270 radw = sdc_cath_thick/sdc_cath_radlen
271 if(ms_flag.gt.0.1) call musc(m2,p,radw,dydzs,dxdzs)
272 drift = 0.5*sdc_thick
273 radw = drift/sdc_radlen
274 jones 1.1 if(ms_flag.gt.0.1) call musc(m2,p,radw,dydzs,dxdzs)
275 drift = drift + sdc_cath_thick
276 call project(xs,ys,drift,decay_flag,dflag,m2,p) !drift and decay
277 radw = sdc_wire_thick/sdc_wire_radlen
278 if(ms_flag.gt.0.1) call musc(m2,p,radw,dydzs,dxdzs)
279 if(wcs_flag.gt.0.1) then
280 call gaussians(random(1),random(2))
281 else
282 random(1) = 0.
283 random(2) = 0.
284 endif
285 xdc(npl_off+iplane)=xs+sdc_sigma(npl_off+iplane)*random(1)*resmult
286 ydc(npl_off+iplane)=ys+sdc_sigma(npl_off+iplane)*random(2)*resmult
287 if (iplane.eq.1 .or. iplane.eq.3 .or. iplane.eq.5) then
288 xdc(npl_off+iplane) = 0. !assume 3 planes are x, 3 are y
289 else
290 ydc(npl_off+iplane) = 0.
291 endif
292 drift = 0.5*sdc_thick
293 radw = drift/sdc_radlen
294 if(ms_flag.gt.0.1) call musc(m2,p,radw,dydzs,dxdzs)
295 jones 1.1 drift = drift + sdc_wire_thick
296 call project(xs,ys,drift,decay_flag,dflag,m2,p) !drift and decay
297 enddo
298 if ( xs.gt.(sdc_1_bot-sdc_1x_offset) .or.
299 > xs.lt.(sdc_1_top-sdc_1x_offset) .or.
300 > ys.gt.(sdc_1_left-sdc_1y_offset) .or.
301 > ys.lt.(sdc_1_right-sdc_1y_offset) ) then
302 sSTOP.detec = sSTOP.detec + 1
303 goto 500
304 endif
305 radw = sdc_cath_thick/sdc_cath_radlen
306 if(ms_flag.gt.0.1) 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 = sdc_2_zpos - sdc_1_zpos - sdc_nr_plan*sdc_del_plane
311 radw = drift/sair_radlen
312 if(ms_flag.gt.0.1) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
313 call project(xs,ys,drift,decay_flag,dflag,m2,p) !drift and decay
314
315 jchamber = 2
316 jones 1.1 npl_off = (jchamber-1)*sdc_nr_plan
317 do iplane = 1,sdc_nr_plan
318 radw = sdc_cath_thick/sdc_cath_radlen
319 if(ms_flag.gt.0.1) call musc(m2,p,radw,dydzs,dxdzs)
320 drift = 0.5*sdc_thick
321 radw = drift/sdc_radlen
322 if(ms_flag.gt.0.1) call musc(m2,p,radw,dydzs,dxdzs)
323 drift = drift + sdc_cath_thick
324 call project(xs,ys,drift,decay_flag,dflag,m2,p) !drift and decay
325 radw = sdc_wire_thick/sdc_wire_radlen
326 if(ms_flag.gt.0.1) call musc(m2,p,radw,dydzs,dxdzs)
327 if(wcs_flag.gt.0.1) then
328 call gaussians(random(1),random(2))
329 else
330 random(1) = 0.
331 random(2) = 0.
332 endif
333 xdc(npl_off+iplane)=xs+sdc_sigma(npl_off+iplane)*random(1)*resmult
334 ydc(npl_off+iplane)=ys+sdc_sigma(npl_off+iplane)*random(2)*resmult
335 if (iplane.eq.1 .or. iplane.eq.3 .or. iplane.eq.5) then
336 xdc(npl_off+iplane) = 0. !assume 3 planes are x, 3 are y
337 jones 1.1 else
338 ydc(npl_off+iplane) = 0.
339 endif
340
341 drift = 0.5*sdc_thick
342 radw = drift/sdc_radlen
343 if(ms_flag.gt.0.1) call musc(m2,p,radw,dydzs,dxdzs)
344 drift = drift + sdc_wire_thick
345 call project(xs,ys,drift,decay_flag,dflag,m2,p) !drift and decay
346 enddo
347 if ( xs.gt.(sdc_2_bot-sdc_2x_offset) .or.
348 > xs.lt.(sdc_2_top-sdc_2x_offset) .or.
349 > ys.gt.(sdc_2_left-sdc_2y_offset) .or.
350 > ys.lt.(sdc_2_right-sdc_2y_offset) ) then
351 sSTOP.detec = sSTOP.detec + 1
352 goto 500
353 endif
354 radw = sdc_cath_thick/sdc_cath_radlen
355 if(ms_flag.gt.0.1) call musc(m2,p,radw,dydzs,dxdzs)
356
357 C at last cathode foil of second drift chamber set, drift to hodoscopes
358 jones 1.1
359 drift = sscin_1y_zpos - sdc_2_zpos - 0.5*sdc_nr_plan*sdc_del_plane
360 radw = drift/sair_radlen
361 if(ms_flag.gt.0.1) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
362 call project(xs,ys,drift,decay_flag,dflag,m2,p) !drift and decay
363 if (ys.lt.(sscin_1_left+sscin_1y_offset) .and.
364 > ys.gt.(sscin_1_right+sscin_1y_offset) .and.
365 > xs.lt.(sscin_1_bot+sscin_1x_offset) .and.
366 > xs.gt.(sscin_1_top+sscin_1x_offset)) sscincount = sscincount + 1
367 radw = sscin_1y_thick/sscin_radlen
368 if(ms_flag.gt.0.1) call musc(m2,p,radw,dydzs,dxdzs)
369
370 drift = sscin_1x_zpos - sscin_1y_zpos
371 radw = drift/sair_radlen
372 if(ms_flag.gt.0.1) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
373 call project(xs,ys,drift,decay_flag,dflag,m2,p) !drift and decay
374 if (ys.lt.(sscin_1_left+sscin_1y_offset) .and.
375 > ys.gt.(sscin_1_right+sscin_1y_offset) .and.
376 > xs.lt.(sscin_1_bot+sscin_1x_offset) .and.
377 > xs.gt.(sscin_1_top+sscin_1x_offset)) sscincount = sscincount + 1
378 radw = sscin_1x_thick/sscin_radlen
379 jones 1.1 if(ms_flag.gt.0.1) call musc(m2,p,radw,dydzs,dxdzs)
380
381 C finished first hodoscope, drift to cherenkov
382
383 drift = scer_zentrance - sscin_1x_zpos
384 radw = drift/sair_radlen
385 if(ms_flag.gt.0.1) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
386 call project(xs,ys,drift,decay_flag,dflag,m2,p) !drift and decay
387
388 radw = scer_entr_thick/scer_entr_radlen
389 if(ms_flag.gt.0.1) call musc(m2,p,radw,dydzs,dxdzs)
390
391 drift = scer_zmirror - scer_zentrance
392 radw = drift/scer_radlen
393 if(ms_flag.gt.0.1) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
394 call project(xs,ys,drift,decay_flag,dflag,m2,p) !drift and decay
395
396 radw = scer_mir_thick/scer_mir_radlen
397 if(ms_flag.gt.0.1) call musc(m2,p,radw,dydzs,dxdzs)
398
399 drift = scer_zexit - scer_zmirror
400 jones 1.1 radw = drift/scer_radlen
401 if(ms_flag.gt.0.1) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
402 call project(xs,ys,drift,decay_flag,dflag,m2,p) !drift and decay
403
404 radw = scer_exit_thick/scer_exit_radlen
405 if(ms_flag.gt.0.1) call musc(m2,p,radw,dydzs,dxdzs)
406
407 C drift to second hodoscope
408
409 drift = sscin_2y_zpos - scer_zexit
410 radw = drift/sair_radlen
411 if(ms_flag.gt.0.1) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
412 call project(xs,ys,drift,decay_flag,dflag,m2,p) !drift and decay
413 if (ys.lt.(sscin_2_left+sscin_2y_offset) .and.
414 > ys.gt.(sscin_2_right+sscin_2y_offset) .and.
415 > xs.lt.(sscin_2_bot+sscin_2x_offset) .and.
416 > xs.gt.(sscin_2_top+sscin_2x_offset)) sscincount = sscincount + 1
417 radw = sscin_2y_thick/sscin_radlen
418 if(ms_flag.gt.0.1) call musc(m2,p,radw,dydzs,dxdzs)
419
420 drift = sscin_2x_zpos - sscin_2y_zpos
421 jones 1.1 radw = drift/sair_radlen
422 if(ms_flag.gt.0.1) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
423 call project(xs,ys,drift,decay_flag,dflag,m2,p) !drift and decay
424 if (ys.lt.(sscin_2_left+sscin_2y_offset) .and.
425 > ys.gt.(sscin_2_right+sscin_2y_offset) .and.
426 > xs.lt.(sscin_2_bot+sscin_2x_offset) .and.
427 > xs.gt.(sscin_2_top+sscin_2x_offset)) sscincount = sscincount + 1
428 radw = sscin_2x_thick/sscin_radlen
429 if(ms_flag.gt.0.1) call musc(m2,p,radw,dydzs,dxdzs)
430
431 C Test on scintillator trigger (DJG 8/18/98)
432 if( sscincount .lt. sscintrig ) then
433 sSTOP.detec = sSTOP.detec + 1
434 goto 500
435 endif
436
437 C Drift to calorimeter
438 drift = scal_4ta_zpos - sscin_2x_zpos
439 radw = drift/sair_radlen
440 if(ms_flag.gt.0.1) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
441 call project(xs,ys,drift,decay_flag,dflag,m2,p) !drift and decay
442 jones 1.1
443 C Note that even with the standard PID trigger, the calorimeter is NOT
444 C required, since the trigger needs either the cerenkov OR the calorimeter.
445 C If you require the calorimeter in you analysis, then you need to make
446 C sure that the calorimeter is hit here, AND apply a seperate post-tracking
447 C fiducial cut (whatever is required by the engine or in your analysis).
448 C The fiducial cut comes later in the code.
449 C
450 C This cut simply requires that the event hit the calorimeter. It does
451 C not insure that the event is within the good region of the calorimeter.
452 C The seperate fiducial cut is required to insure that the entire shower
453 C energy is contained in the calorimeter. Or, if you like, you can require
454 C some distance between the track and the edge. 2-3 cm seems to be enough
455 C to get most or all of the energy in the calorimeter
456
457 !c if (ys.gt.(scal_left-2.0) .or. ys.lt.(scal_right+2.0) .or.
458 !c > xs.gt.(scal_bottom-2.0) .or. xs.lt.(scal_top+2.0)) then
459 ! if (ys.gt.scal_left .or. ys.lt.scal_right .or.
460 ! > xs.gt.scal_bottom .or. xs.lt.scal_top) then
461 ! sSTOP.detec = sSTOP.detec + 1
462 ! goto 500
463 jones 1.1 ! endif
464
465 C and fit track to give new focal plane values, use LFIT from GENLIB
466
467 do jchamber=1,sdc_nr_cham
468 npl_off = (jchamber-1)*sdc_nr_plan
469 do iplane=1,sdc_nr_plan
470 if (jchamber.eq.1) zdc(npl_off+iplane) = sdc_1_zpos +
471 > (iplane-0.5-0.5*sdc_nr_plan)*sdc_del_plane
472 if (jchamber.eq.2) zdc(npl_off+iplane) = sdc_2_zpos +
473 > (iplane-0.5-0.5*sdc_nr_plan)*sdc_del_plane
474 enddo
475 enddo
476
477 call lfit(zdc,xdc,12,0,dx,x0,badf)
478 call lfit(zdc,ydc,12,0,dy,y0,badf)
479
480 write(6,*) 'xdc=',xdc
481 write(6,*) 'ydc=',ydc
482 write(6,*) 'zdc=',zdc
483 write(6,*) 'x,y,dx,dy=',x0,y0,dx,dy
484 jones 1.1 if (x0.lt.-20) stop
485
486 x_fp = x0
487 y_fp = y0
488 dx_fp = dx
489 dy_fp = dy
490
491 C If you use a calorimeter cut in your analysis, the engine applied a
492 C a fiducial cut at the calorimeter. This is based on the position of the
493 C TRACK at the calorimeter, not the real position of the event. Go to the
494 C back of the calorimeter since engine uses a FID cut at the back.
495 C The standard fiducial cut is 5 cm from the edges of the block.
496
497 xcal = x_fp + dx_fp * scal_4ta_zpos
498 ycal = y_fp + dy_fp * scal_4ta_zpos
499 ! if (ycal.gt.(scal_left-5.0) .or. ycal.lt.(scal_right+5.0) .or.
500 ! > xcal.gt.(scal_bottom-5.0) .or. xcal.lt.(scal_top+5.0)) then
501 ! sSTOP.detec = sSTOP.detec + 1
502 ! goto 500
503 ! endif
504
505 jones 1.1 ok_hut = 1.
506
507 C We are done with this event, whether GOOD or BAD.
508
509 500 continue
510
511 return
512 end
|