1 gaskelld 1.1 	subroutine mc_shms_hut (m2,p,x_fp,dx_fp,y_fp,dy_fp,ms_flag,wcs_flag,
  2                   >		decay_flag,dflag,resmult,ok_hut,zinit,pathlen,cerflag)
  3              C----------------------------------------------------------------------
  4              C
  5              C Monte-Carlo of HMS detector hut.
  6              C	Note that we only pass on real*4 variables to the subroutine.
  7              C	This will make life easier for portability.
  8              C----------------------------------------------------------------------
  9              
 10              	implicit 	none
 11              
 12              	include 'struct_shms.inc'
 13              	include '../spectrometers.inc'
 14              	include 'hut.inc'
 15              
 16              C Math constants
 17              
 18              	real*8 pi,d_r,r_d,root
 19              
 20              	parameter (pi = 3.141592654)
 21              	parameter (d_r = pi/180.)
 22 gaskelld 1.1 	parameter (r_d = 180./pi)
 23              	parameter (root = 0.707106781)		!square root of 1/2
 24              
 25              C all parameters, later to take from .parm files
 26              C The arguments
 27              
 28              	real*8	p,m2			!momentum and mass of particle
 29              	real*8	x_fp,y_fp,dx_fp,dy_fp	!Focal plane values to return
 30              	logical	ms_flag, wcs_flag	!particle, m_scat, wc_smear
 31              	logical	ok_hut			!true if particle makes it
 32              	logical decay_flag,dflag
 33              	real*8 xcal,ycal		!Position of track at calorimeter.
 34              	real*8 pathlen
 35              	real*8 resmult			!DC resolution factor
 36              	real*8 zinit
 37              
 38              c external function
 39              	real*8 gauss1
 40              C Local declarations.
 41              
 42              	integer*4	i,iplane,jchamber,npl_off
 43 gaskelld 1.1 	integer*4	chan	/1/
 44              
 45              	real*8	radw,drift
 46              c need real*4 for cernlib routine lfut
 47              	real*4	xdc(12),ydc(12),zdc(12)		!positions at d.c. planes
 48              	real*4	x0,y0,dx,dy			!track fitting temporaries
 49              	real*4	badf				!temporaries
 50              
 51              	real*8  tmpran1,tmpran2
 52              	real*8   nsig_max
 53              	real*8 hdc_del_plane
 54              	parameter (nsig_max=99.0d0)
 55              c mkj
 56              	logical use_det_cut
 57              	parameter (use_det_cut=.true.)
 58              	logical cerflag
 59              
 60              C ================================ Executable Code =============================
 61              
 62              C Initialize some variables
 63              
 64 gaskelld 1.1 *	write (6,*) x_fp, y_fp, dx_fp, dy_fp
 65              	hdc_del_plane = hdc_thick + hdc_wire_thick + hdc_cath_thick
 66              
 67              C Initialize ok_hut to zero
 68              
 69              	ok_hut = .false.
 70              c
 71              	resmult= 1.0
 72              
 73              C Initialize the xdc and ydc arrays to zero
 74              
 75              	do i=1,12
 76              	  xdc(i) = 0.
 77              	  ydc(i) = 0.
 78              	enddo
 79              
 80              C Sanity check (and avoid unused variable on zinit)
 81              	if (abs(zinit).gt.1e-10) then
 82              	  write(6,*) 'zinit not equal to zero in call to shms/mc_shms_hut.f'
 83              	  write(6,*) 'code will ingore initial value: zinit=',zinit,'!!!'
 84              	endif
 85 gaskelld 1.1 
 86              C------------------------------------------------------------------------------C
 87              C                           Top of loop through hut                            C
 88              C------------------------------------------------------------------------------C
 89              
 90              C Go to the Cerenkov. (Drift back from focal plane). Cherenkov has 0.5 atm of Freon, 
 91              c and is located 4.0 meters behind the focal point.
 92              c The Cherenkov is 3m long, and has an Aluminum entrance window.
 93              c
 94              c only have cerenkov when p > 7500 MeV
 95              c
 96              	drift = hcer_1_zentrance
 97              c	call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
 98              c set decay flag false to avoid double counting
 99              	call project(xs,ys,drift,.false.,dflag,m2,p,pathlen)
100              
101              	radw = hcer_entr_thick/hcer_entr_radlen
102              	if(ms_flag .and. cerflag ) call musc(m2,p,radw,dydzs,dxdzs)
103              
104              
105              	drift = hcer_1_zmirror - hcer_1_zentrance
106 gaskelld 1.1 	radw = drift/hcer_1_radlen
107              c	call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
108              	call project(xs,ys,drift,.false.,dflag,m2,p,pathlen)
109              	if(ms_flag .and. cerflag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
110              
111              
112              	radw = hcer_mir_thick/hcer_mir_radlen
113              	if(ms_flag .and. cerflag) call musc(m2,p,radw,dydzs,dxdzs)
114              
115              	drift = hcer_1_zexit - hcer_1_zmirror
116              	radw = drift/hcer_1_radlen
117              c	call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
118              	call project(xs,ys,drift,.false.,dflag,m2,p,pathlen)
119              	if(ms_flag .and. cerflag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
120              
121              	radw = hcer_exit_thick/hcer_exit_radlen
122              	if(ms_flag .and. cerflag) call musc(m2,p,radw,dydzs,dxdzs)
123              	if(ms_flag .and. .not. cerflag) then
124              	  radw = hcer_entr_thick/hcer_entr_radlen
125                        call musc(m2,p,radw,dydzs,dxdzs)
126                      endif
127 gaskelld 1.1 
128              
129              
130              C Go to first drift chamber set
131              C For simplicity, perform air MS (probably negligeable) at before drift
132              C instead of 1/2 way through.
133              
134              *	write (6,*) 'made it to the first drift chamber'
135              	drift = (hdc_1_zpos - 0.5*hdc_nr_plan*hdc_del_plane) - hcer_1_zexit
136              	radw = drift/hair_radlen
137              	call project(xs,ys,drift,.false.,dflag,m2,p,pathlen)
138              	if (ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
139              
140              	jchamber = 1
141              	radw = hdc_entr_thick/hdc_entr_radlen
142              	if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
143              	npl_off = (jchamber-1)*hdc_nr_plan
144              	do iplane = 1,hdc_nr_plan
145              	  radw = hdc_cath_thick/hdc_cath_radlen
146              	  if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
147              	  drift = 0.5*hdc_thick
148 gaskelld 1.1 	  radw = drift/hdc_radlen
149              	  if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
150              	  drift = drift + hdc_cath_thick
151              	  call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
152              	  radw = hdc_wire_thick/hdc_wire_radlen
153              	  if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
154              	  if(wcs_flag) then
155              	    tmpran1 = gauss1(nsig_max)	!gaussian, truncated at 99 sigma.
156              	    tmpran2 = gauss1(nsig_max)
157              	  else
158              	    tmpran1 = 0.
159              	    tmpran2 = 0.
160              	  endif
161              	  xdc(npl_off+iplane) = xs + hdc_sigma(npl_off+iplane)*tmpran1*resmult
162              	  ydc(npl_off+iplane) = ys + hdc_sigma(npl_off+iplane)*tmpran2*resmult
163              	  if (iplane.eq.2 .or. iplane.eq.5) then
164              	    xdc(npl_off+iplane) = 0.   !y plane, no x information
165              	  else
166              c	     write(*,*) ' iplane = ',iplane
167              	    ydc(npl_off+iplane) = 0.   !x-like plane, no y info
168              	  endif
169 gaskelld 1.1 *	  write (6,*) 'i am still in the first drift chamber'
170              	  drift = 0.5*hdc_thick
171              	  radw = drift/hdc_radlen
172              	  if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
173              	  drift = drift + hdc_wire_thick
174              	  call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
175              	enddo
176              	radw = hdc_exit_thick/hdc_exit_radlen
177              	if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
178              	if (xs.gt.(hdc_1_bot-hdc_1x_offset) .or.
179                   >      xs.lt.(hdc_1_top-hdc_1x_offset) .or.
180                   >      ys.gt.(hdc_1_left-hdc_1y_offset) .or.
181                   >      ys.lt.(hdc_1_right-hdc_1y_offset) ) then
182              	  shmsSTOP_dc1 = shmsSTOP_dc1 + 1
183              c	  write(6,*) 'Lost in DC! delta,xp,yp',dpps,dxdzs,dydzs
184              	  if (use_det_cut) goto 500
185              	endif
186              	radw = hdc_cath_thick/hdc_cath_radlen
187              	if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
188              
189              C at last cathode foil of first drift chamber set, drift to next
190 gaskelld 1.1 
191              *	write (6,*) 'made it to the second drift chamber in the hut'
192              	drift = hdc_2_zpos - hdc_1_zpos - hdc_nr_plan*hdc_del_plane
193              C       Break this into 2 parts to properly account for decay.
194              C       We've already done decay up to the half-way point between the chambers.
195              	call project(xs,ys,drift/2.0,.false.,dflag,m2,p,pathlen)
196              	call project(xs,ys,drift/2.0,decay_flag,dflag,m2,p,pathlen)
197              	radw = drift/hair_radlen
198              	if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
199              
200              	jchamber = 2
201              	radw = hdc_entr_thick/hdc_entr_radlen
202              	if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
203              	npl_off = (jchamber-1)*hdc_nr_plan
204              	do iplane = 1,hdc_nr_plan
205              	  radw = hdc_cath_thick/hdc_cath_radlen
206              	  if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
207              	  drift = 0.5*hdc_thick
208              	  radw = drift/hdc_radlen
209              	  if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
210              	  drift = drift + hdc_cath_thick
211 gaskelld 1.1 	  call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
212              	  radw = hdc_wire_thick/hdc_wire_radlen
213              	  if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
214              	  if(wcs_flag) then
215              	    tmpran1 = gauss1(nsig_max)	!gaussian, truncated at 99 sigma.
216              	    tmpran2 = gauss1(nsig_max)
217              	  else
218              	    tmpran1 = 0.
219              	    tmpran2 = 0.
220              	  endif
221              	  xdc(npl_off+iplane) = xs + hdc_sigma(npl_off+iplane)*tmpran1*resmult
222              	  ydc(npl_off+iplane) = ys + hdc_sigma(npl_off+iplane)*tmpran2*resmult
223              	  if (iplane.eq.2 .or. iplane.eq.5) then
224              	    xdc(npl_off+iplane) = 0.   !y plane, no x information
225              	  else
226              	    ydc(npl_off+iplane) = 0.   !x-like plane, no y info
227              	  endif
228              	  drift = 0.5*hdc_thick
229              	  radw = drift/hdc_radlen
230              	  if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
231              	  drift = drift + hdc_wire_thick
232 gaskelld 1.1 	  call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
233              	enddo
234              	radw = hdc_exit_thick/hdc_exit_radlen
235              	if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
236              	if (xs.gt.(hdc_2_bot-hdc_2x_offset) .or.
237                   >      xs.lt.(hdc_2_top-hdc_2x_offset) .or.
238                   >      ys.gt.(hdc_2_left-hdc_2y_offset) .or.
239                   >      ys.lt.(hdc_2_right-hdc_2y_offset) ) then
240              	  shmsSTOP_dc2 = shmsSTOP_dc2 + 1
241              	  if (use_det_cut) goto 500
242              	endif
243              	radw = hdc_cath_thick/hdc_cath_radlen
244              	if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
245              
246              C at last cathode foil of second drift chamber set, drift to the 1st hodoscope
247              
248              	drift = hscin_1x_zpos - hdc_2_zpos - 0.5*hdc_nr_plan*hdc_del_plane
249              	radw = drift/hair_radlen
250              	call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
251              	if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
252              	if (ys.gt.(hscin_1x_left+hscin_1y_offset) .or.
253 gaskelld 1.1      >      ys.lt.(hscin_1x_right+hscin_1y_offset)) then
254              	  shmsSTOP_s1 = shmsSTOP_s1 + 1
255              	  if (use_det_cut) goto 500
256              	endif
257              	radw = hscin_1x_thick/hscin_radlen
258              	if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
259              	drift = hscin_1y_zpos - hscin_1x_zpos
260              	radw = drift/hair_radlen
261              	call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
262              	if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
263              	if (xs.gt.(hscin_1y_bot+hscin_1x_offset) .or.
264                   >      xs.lt.(hscin_1y_top+hscin_1x_offset)) then
265              	  shmsSTOP_s1 = shmsSTOP_s1 + 1
266              	  if (use_det_cut) goto 500
267              	endif
268               	radw = hscin_1y_thick/hscin_radlen
269              	if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
270              
271              C finished first hodoscope, drift to the second cherenkov
272              
273              *	write (6,*) 'made it to the second cherenkov in the hut'
274 gaskelld 1.1 	drift = hcer_2_zentrance - hscin_1y_zpos
275              	radw = drift/hair_radlen
276              	call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
277              	if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
278              
279              	radw = hcer_2_entr_thick/hcer_2_entr_radlen
280              	if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
281              
282              	drift = hcer_2_zmirror - hcer_2_zentrance
283              	radw = drift/hcer_2_radlen
284              	call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
285              	if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
286              
287              	radw = hcer_mir_thick/hcer_mir_radlen
288              
289              	if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
290              
291              	drift = hcer_2_zexit - hcer_2_zmirror
292              	radw = drift/hcer_2_radlen
293              	call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
294              	if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
295 gaskelld 1.1 
296              	radw = hcer_2_exit_thick/hcer_2_exit_radlen
297              	if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
298              
299              C drift to 2nd hodoscope
300              
301              	drift = hscin_2x_zpos - hcer_2_zexit
302              	radw = drift/hair_radlen
303              	call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
304              	if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
305              	if (ys.gt.(hscin_2x_left+hscin_2y_offset) .or.
306                   >      ys.lt.(hscin_2x_right+hscin_2y_offset)) then
307              	  shmsSTOP_s3 = shmsSTOP_s3 + 1
308              	  if (use_det_cut) goto 500
309              	endif
310              	radw = hscin_2x_thick/hscin_radlen
311              	if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
312              	drift = hscin_2y_zpos - hscin_2x_zpos
313              	radw = drift/hair_radlen
314              	call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
315              	if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
316 gaskelld 1.1 	if (xs.gt.(hscin_2y_bot+hscin_2x_offset) .or.
317                   >      xs.lt.(hscin_2y_top+hscin_2x_offset)) then
318              	  shmsSTOP_s2 = shmsSTOP_s2 + 1
319              	  if (use_det_cut) goto 500
320              	endif
321              	radw = hscin_2y_thick/hscin_radlen
322              	if(ms_flag) call musc(m2,p,radw,dydzs,dxdzs)
323              
324              C Don't need to drift to calorimeter unless it's required in your trigger.
325              C Note that even with the standard PID trigger, the calorimeter is NOT
326              C required, since the trigger needs either the cerenkov OR the calorimeter.
327              C There is a seperate fiducial cut needed if you require the calorimeter
328              C in you analysis.  That cut is applied AFTER fitting the track (see below).
329              
330              	drift = hcal_4ta_zpos - hscin_2y_zpos
331              	radw = drift/hair_radlen
332              	call project(xs,ys,drift,decay_flag,dflag,m2,p,pathlen)
333              	if(ms_flag) call musc_ext(m2,p,radw,drift,dydzs,dxdzs,ys,xs)
334              	if (ys.gt.hcal_left .or. ys.lt.hcal_right .or.
335                   >	   xs.gt.hcal_bottom .or. xs.lt.hcal_top) then
336              	  shmsSTOP_cal = shmsSTOP_cal + 1
337 gaskelld 1.1 	  if (use_det_cut) goto 500
338              	endif
339              
340              
341              C and fit track to give new focal plane values, use LFIT from GENLIB
342              
343              	do jchamber=1,hdc_nr_cham
344              	  npl_off = (jchamber-1)*hdc_nr_plan
345              	  do iplane=1,hdc_nr_plan
346              	    if (jchamber.eq.1) zdc(npl_off+iplane) = hdc_1_zpos +
347                   >          (iplane-0.5-0.5*hdc_nr_plan)*hdc_del_plane
348              	    if (jchamber.eq.2) zdc(npl_off+iplane) = hdc_2_zpos +
349                   >          (iplane-0.5-0.5*hdc_nr_plan)*hdc_del_plane
350              	  enddo
351              	enddo
352              	call lfit(zdc,xdc,12,0,dx,x0,badf)
353              	call lfit(zdc,ydc,12,0,dy,y0,badf)
354              
355              
356              	x_fp = x0
357              	y_fp = y0
358 gaskelld 1.1 	dx_fp = dx
359              	dy_fp = dy
360              
361              
362              C If you use a calorimeter cut in your analysis, the engine applied a
363              C a fiducial cut at the calorimeter.  This is based on the position of the
364              C TRACK at the calorimeter, not the real position of the event.  Go to the
365              C back of the calorimeter since engine uses a FID cut at the back.
366              C The standard fiducial cut is 5 cm from the edges of the block.
367              
368              *	write (6,*) 'at the shower counter in the hut'
369              	xcal = x_fp + dx_fp * hcal_4ta_zpos
370              	ycal = y_fp + dy_fp * hcal_4ta_zpos
371              	if (ycal.gt.(hcal_left-5.0) .or. ycal.lt.(hcal_right+5.0) .or.
372                   >	   xcal.gt.(hcal_bottom-5.0) .or. xcal.lt.(hcal_top+5.0)) then
373              	  shmsSTOP_cal = shmsSTOP_cal + 1
374              	  if (use_det_cut) goto 500
375              	endif
376              
377              	ok_hut = .true.
378              
379 gaskelld 1.1 C We are done with this event, whether GOOD or BAD.
380              
381              500	continue
382              
383              C ALL done!
384              	close (unit=chan)
385              
386              	return
387              	end