(file) Return to mc_sos_hut.f.OLD CVS log (file) (dir) Up to [HallC] / Poltar / sos

  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

Analyzer/Replay: Mark Jones, Documents: Stephen Wood
Powered by
ViewCVS 0.9.2-cvsgraph-1.4.0