Analyzer/HTRACKING/h_physics.f - annotate

Return to h_physics.f CVS log

Up to [HallC] / Analyzer / HTRACKING

1 cdaq 1.2 SUBROUTINE H_PHYSICS(ABORT,err)
2 cdaq 1.1 -------------------------------------------------------- 3 - 4 - Purpose and Methods : Do final HMS physics analysis on HMS only part of 5 - event. 6 - 7 - to decoded information 8 - 9 - Required Input BANKS HMS_FOCAL_PLANE 10 - HMS_TARGET 11 - HMS_TRACK_TESTS 12 *-
13 cdaq 1.2 - Output BANKS HMS_PHYSICS_R4 14 - HMS_PHYSICS_I4
15 cdaq 1.1 - 16 - Output: ABORT - success or failure 17 - : err - reason for failure, if any 18 - 19 - Created 19-JAN-1994 D. F. Geesaman 20 - Dummy Shell routine
21 cdaq 1.2 * $Log: h_physics.f,v $
22 saw 1.18 * Revision 1.17 1996/04/30 12:46:06 saw 23 * (JRA) Add pathlength and rf calculations 24 *
25 saw 1.17 * Revision 1.16 1996/01/24 15:58:38 saw 26 * (JRA) Change cpbeam/cebeam to gpbeam/gebeam 27 *
28 saw 1.16 * Revision 1.15 1996/01/16 21:55:02 cdaq 29 * (JRA) Calculate q, W for electrons 30 *
31 cdaq 1.15 * Revision 1.14 1995/10/09 20:22:15 cdaq 32 * (JRA) Add call to h_dump_cal, change upper to lower case 33 *
34 cdaq 1.14 * Revision 1.13 1995/08/31 14:49:03 cdaq 35 * (JRA) Add projection to cerenkov mirror pos, fill hdc_sing_res array 36 *
37 cdaq 1.13 * Revision 1.12 1995/07/19 20:53:26 cdaq 38 * (SAW) Declare sind and tand for f2c compatibility 39 *
40 cdaq 1.12 * Revision 1.11 1995/05/22 19:39:15 cdaq 41 * (SAW) Split gen_data_data_structures into gen, hms, sos, and coin parts" 42 *
43 cdaq 1.11 * Revision 1.10 1995/05/11 17:15:07 cdaq 44 * (SAW) Add additional kinematics variables 45 *
46 cdaq 1.10 * Revision 1.9 1995/03/22 16:23:27 cdaq 47 * (SAW) Target track data is now slopes. 48 *
49 cdaq 1.9 * Revision 1.8 1995/02/23 13:37:31 cdaq 50 * (SAW) Reformat and cleanup 51 *
52 cdaq 1.8 * Revision 1.7 1995/02/10 18:44:47 cdaq 53 * (SAW) _tar values are now angles instead of slopes 54 *
55 cdaq 1.7 * Revision 1.6 1995/02/02 13:05:40 cdaq 56 * (SAW) Moved best track selection code into H_SELECT_BEST_TRACK (new) 57 *
58 cdaq 1.6 * Revision 1.5 1995/01/27 20:24:14 cdaq 59 * (JRA) Add some useful physics quantities 60 *
61 cdaq 1.5 * Revision 1.4 1995/01/18 16:29:26 cdaq 62 * (SAW) Correct some trig and check for negative arg in elastic kin calculation 63 *
64 cdaq 1.4 * Revision 1.3 1994/09/13 19:51:03 cdaq 65 * (JRA) Add HBETA_CHISQ 66 *
67 cdaq 1.3 * Revision 1.2 1994/06/14 03:49:49 cdaq 68 * (DFG) Calculate physics quantities 69 *
70 cdaq 1.2 * Revision 1.1 1994/02/19 06:16:08 cdaq 71 * Initial revision 72 *
73 cdaq 1.1 - 74 - 75 *--------------------------------------------------------
76 cdaq 1.2 IMPLICIT NONE 77 SAVE
78 cdaq 1.1 *
79 cdaq 1.10 character*9 here
80 cdaq 1.2 parameter (here= 'H_PHYSICS')
81 cdaq 1.1 *
82 cdaq 1.2 logical ABORT 83 character() err
84 cdaq 1.6 integer ierr
85 cdaq 1.1 *
86 cdaq 1.11 include 'gen_data_structures.cmn' 87 INCLUDE 'hms_data_structures.cmn'
88 cdaq 1.2 INCLUDE 'gen_routines.dec' 89 INCLUDE 'gen_constants.par' 90 INCLUDE 'gen_units.par' 91 INCLUDE 'hms_physics_sing.cmn'
92 cdaq 1.5 INCLUDE 'hms_calorimeter.cmn'
93 cdaq 1.8 INCLUDE 'hms_scin_parms.cmn'
94 cdaq 1.10 INCLUDE 'hms_tracking.cmn'
95 cdaq 1.13 INCLUDE 'hms_cer_parms.cmn' 96 INCLUDE 'hms_geometry.cmn'
97 cdaq 1.15 INCLUDE 'hms_id_histid.cmn' 98 INCLUDE 'hms_track_histid.cmn' 99 include 'gen_event_info.cmn'
100 saw 1.17 include 'hms_scin_tof.cmn'
101 cdaq 1.1 *
102 cdaq 1.2 * local variables
103 cdaq 1.15 integer4 i,ip,ihit 104 integer4 itrkfp
105 cdaq 1.6 real4 cosgamma,tandelphi,sinhphi,coshstheta,sinhstheta 106 real4 t1,ta,p3,t3,hminv2
107 cdaq 1.10 real*4 coshsthetaq
108 cdaq 1.12 real*4 sind,tand ! For f2c
109 cdaq 1.14 real*4 p_nonzero
110 cdaq 1.15 real4 xdist,ydist,dist(12),res(12) 111 real4 tmp,W2
112 saw 1.17 real*4 denom
113 cdaq 1.6 *
114 cdaq 1.1 -------------------------------------------------------- 115
116 saw 1.18 ierr=0
117 cdaq 1.15 hphi_lab=0.0
118 cdaq 1.14 if(hsnum_fptrack.gt.0) then ! Good track has been selected
119 cdaq 1.15 itrkfp=hsnum_fptrack
120 cdaq 1.14 hsp = hp_tar(hsnum_tartrack) 121 hsenergy = sqrt(hsphsp+hpartmasshpartmass)
122 cdaq 1.2 * Copy variables for ntuple so we can test on them
123 cdaq 1.14 hsdelta = hdelta_tar(hsnum_tartrack) 124 hsx_tar = hx_tar(hsnum_tartrack) 125 hsy_tar = hy_tar(hsnum_tartrack) 126 hsxp_tar = hxp_tar(hsnum_tartrack) ! This is an angle (radians) 127 hsyp_tar = hyp_tar(hsnum_tartrack) ! This is an angle (radians)
128 cdaq 1.15 hsbeta = hbeta(itrkfp) 129 hsbeta_chisq = hbeta_chisq(itrkfp) 130 hstime_at_fp = htime_at_fp(itrkfp)
131 cdaq 1.14
132 cdaq 1.15 hstrack_et = htrack_et(itrkfp) 133 hstrack_preshower_e = htrack_preshower_e(itrkfp)
134 cdaq 1.14 p_nonzero = max(.0001,hsp) !momentum (used to normalize calorim.) 135 hscal_suma = hcal_e1/p_nonzero !normalized cal. plane sums 136 hscal_sumb = hcal_e2/p_nonzero 137 hscal_sumc = hcal_e3/p_nonzero 138 hscal_sumd = hcal_e4/p_nonzero 139 hsprsum = hscal_suma 140 hsshsum = hcal_et/p_nonzero 141 hsprtrk = hstrack_preshower_e/p_nonzero 142 hsshtrk = hstrack_et/p_nonzero 143
144 cdaq 1.15 hsx_sp1 = hx_sp1(itrkfp) 145 hsy_sp1 = hy_sp1(itrkfp) 146 hsxp_sp1= hxp_sp1(itrkfp) 147 hsx_sp2 = hx_sp2(itrkfp) 148 hsy_sp2 = hy_sp2(itrkfp) 149 hsxp_sp2= hxp_sp2(itrkfp) 150 151 do ihit=1,hnum_scin_hit(itrkfp) 152 call hf1(hidscintimes,hscin_fptime(itrkfp,ihit),1.) 153 enddo 154 155 do ihit=1,hntrack_hits(itrkfp,1) 156 call hf1(hidcuttdc, 157 & float(hdc_tdc(hntrack_hits(itrkfp,ihit+1))),1.) 158 enddo 159 160 hsx_fp = hx_fp(itrkfp) 161 hsy_fp = hy_fp(itrkfp) 162 hsxp_fp = hxp_fp(itrkfp) ! This is a slope (dx/dz) 163 hsyp_fp = hyp_fp(itrkfp) ! This is a slope (dy/dz)
164 cdaq 1.13 hsx_dc1 = hsx_fp + hsxp_fp * hdc_1_zpos 165 hsy_dc1 = hsy_fp + hsyp_fp * hdc_1_zpos 166 hsx_dc2 = hsx_fp + hsxp_fp * hdc_2_zpos 167 hsy_dc2 = hsy_fp + hsyp_fp * hdc_2_zpos
168 cdaq 1.8 hsx_s1 = hsx_fp + hsxp_fp * hscin_1x_zpos 169 hsy_s1 = hsy_fp + hsyp_fp * hscin_1x_zpos
170 cdaq 1.13 hsx_cer = hsx_fp + hsxp_fp * hcer_mirror_zpos 171 hsy_cer = hsy_fp + hsyp_fp * hcer_mirror_zpos
172 cdaq 1.8 hsx_s2 = hsx_fp + hsxp_fp * hscin_2x_zpos 173 hsy_s2 = hsy_fp + hsyp_fp * hscin_2x_zpos 174 hsx_cal = hsx_fp + hsxp_fp * hcal_1pr_zpos 175 hsy_cal = hsy_fp + hsyp_fp * hcal_1pr_zpos
176 cdaq 1.5
177 saw 1.17 hsbeta_p = hsp/max(hsenergy,.00001)
178 saw 1.18 C old 'fit' value for pathlen correction 179 C hspathlength = -1.47e-2hsx_fp + 11.6hsxp_fp - 36hsxp_fp2 180 C new 'modeled' value. 181 hspathlength = 12.462hsxp_fp + 0.1138hsxp_fphsx_fp 182 & -0.0154hsx_fp - 72.292hsxp_fp*2 183 & -0.0000544hsx_fp*2 - 116.52hsyp_fp**2 184
185 saw 1.17 hspath_cor = hspathlength/hsbeta_p - 186 & hpathlength_central/speed_of_light*(1/max(.01,hsbeta_p) - 1) 187 188 hsrftime = hmisc_dec_data(8,1)/9.46 189 & - (hstime_at_fp-hstart_time_center) - hspath_cor
190 cdaq 1.8 do ip=1,4 191 hsscin_elem_hit(ip)=0 192 enddo
193 cdaq 1.15 do i=1,hnum_scin_hit(itrkfp) 194 ip=hscin_plane_num(hscin_hit(itrkfp,i))
195 cdaq 1.8 if (hsscin_elem_hit(ip).eq.0) then
196 cdaq 1.10 hsscin_elem_hit(ip)=hscin_counter_num(hscin_hit(
197 cdaq 1.15 $ itrkfp,i)) 198 hsdedx(ip)=hdedx(itrkfp,i)
199 cdaq 1.8 else ! more than 1 hit in plane 200 hsscin_elem_hit(ip)=18
201 cdaq 1.15 hsdedx(ip)=sqrt(hsdedx(ip)*hdedx(itrkfp,i))
202 cdaq 1.8 endif 203 enddo
204 cdaq 1.5
205 cdaq 1.15 hsnum_scin_hit = hnum_scin_hit(itrkfp) 206 hsnum_pmt_hit = hnum_pmt_hit(itrkfp) 207 208 hschi2perdeg = hchi2_fp(itrkfp) 209 $ /float(hnfree_fp(itrkfp)) 210 hsnfree_fp = hnfree_fp(itrkfp)
211 cdaq 1.6
212 cdaq 1.13 do ip = 1, hdc_num_planes
213 cdaq 1.15 hdc_sing_res(ip)=hdc_single_residual(itrkfp,ip) 214 hsdc_track_coord(ip)=hdc_track_coord(itrkfp,ip)
215 cdaq 1.13 enddo
216 cdaq 1.15 217 if (hntrack_hits(itrkfp,1).eq.12 .and. hschi2perdeg.le.4) then 218 xdist=hsx_dc1 219 ydist=hsy_dc1 220 do ip=1,12 221 if (hdc_readout_x(ip)) then 222 dist(ip) = ydist*hdc_readout_corr(ip)
223 saw 1.16 else !readout from top/bottom
224 cdaq 1.15 dist(ip) = xdist*hdc_readout_corr(ip) 225 endif 226 res(ip)=hdc_sing_res(ip)
227 saw 1.16 tmp = hdc_plane_wirecoord(itrkfp,ip) 228 $ -hdc_plane_wirecenter(itrkfp,ip) 229 if (tmp.eq.0) then !drift dist = 0
230 cdaq 1.15 res(ip)=abs(res(ip)) 231 else
232 saw 1.16 res(ip)=res(ip) * (abs(tmp)/tmp) !convert +/- res to near/far res
233 cdaq 1.15 endif 234 enddo 235 c write(37,'(12f7.2,12f8.3,12f8.5)') (hsdc_track_coord(ip),ip=1,12), 236 c & (dist(ip),ip=1,12),(res(ip),ip=1,12) 237 endif 238 239 hstheta =htheta_labTT/180. - hsyp_tar 240 HSP = hpcentral(1 + hsdelta/100.)
241 cdaq 1.14 sinhstheta = sin(hstheta)
242 saw 1.18 coshstheta = cos(hstheta)
243 cdaq 1.9 tandelphi = hsxp_tar / 244 & ( sinhthetas - coshthetas*hsyp_tar)
245 saw 1.16 hsphi = hphi_lab + atan(tandelphi) ! hphi_lab MUST BE MULTIPLE OF
246 cdaq 1.14 sinhphi = sin(hsphi) ! PI/2, OR ABOVE IS CRAP
247 saw 1.16 if(hpartmass .lt. 2*mass_electron) then ! Less than 1 MeV, HMS is elec 248 if(gtarg_z(gtarg_num).gt.0.)then 249 call total_eloss(1,.true.,gtarg_z(gtarg_num), 250 $ gtarg_a(gtarg_num),gtarg_thick(gtarg_num), 251 $ gtarg_dens(gtarg_num), 252 $ hstheta,gtarg_theta,1.0,hseloss)
253 saw 1.18 hsenergy=hsenergy- hseloss
254 saw 1.16 else 255 hseloss=0. 256 endif
257 saw 1.18 hqx = -HSPcos(HSxp_tar)sinhstheta
258 saw 1.16 hqy = -HSP*sin(Hsxp_tar)
259 saw 1.18 hqz = gpbeam - HSPcos(HSxp_tar)coshstheta
260 saw 1.16 hqabs= sqrt(hqx2+hqy2+hqz2) 261 W2 = gtarg_mass(gtarg_num)2 + 262 $ 2.gtarg_mass(gtarg_num)(gpbeam-hsp) - hqabs2 + 263 $ (gpbeam-hsp)2 264 if(W2.ge.0 ) then
265 saw 1.18 hinvmass = SQRT(W2)
266 saw 1.16 else
267 saw 1.18 hinvmass = 0.
268 saw 1.16 endif
269 cdaq 1.15 else
270 saw 1.16 if(gtarg_z(gtarg_num).gt.0.)then 271 call total_eloss(1,.false.,gtarg_z(gtarg_num), 272 $ gtarg_a(gtarg_num), 273 $ gtarg_thick(gtarg_num),gtarg_dens(gtarg_num), 274 $ hstheta,gtarg_theta,hsbeta,hseloss)
275 saw 1.18 hsenergy = hsenergy - hseloss
276 saw 1.16 else 277 hseloss=0. 278 endif
279 cdaq 1.15 endif
280 cdaq 1.2 * Calculate elastic scattering kinematics
281 saw 1.16 t1 = 2.hphysicsagpbeamcoshstheta 282 ta = 4gpbeam*2coshstheta2 - hphysicsb2
283 cdaq 1.8 ccc SAW 1/17/95. Add the stuff after the or.
284 cdaq 1.14 if(ta.eq.0.0 .or. ( hphysicab2 + hphysicsm3b * ta).lt.0.0) then
285 cdaq 1.8 p3=0. 286 else
287 cdaq 1.14 t3 = ta-hphysicsb*2 288 p3 = (T1 - sqrt( hphysicab2 + hphysicsm3b ta)) / ta
289 cdaq 1.8 endif
290 cdaq 1.2 * This is the difference in the momentum obtained by tracking 291 * and the momentum from elastic kinematics
292 cdaq 1.14 hselas_cor = hsp - P3 293 * invariant mass of the remaining particles
294 saw 1.16 hminv2 = ( (gebeam+gtarg_mass(gtarg_num)-hsenergy)*2 295 & - (gpbeam - hsp coshstheta)**2
296 cdaq 1.14 & - ( hsp * sinhstheta)**2 )
297 cdaq 1.8 if(hminv2.ge.0 ) then
298 cdaq 1.14 hsminv = sqrt(hminv2)
299 cdaq 1.8 else
300 cdaq 1.14 hsminv = 0.
301 cdaq 1.8 endif ! end test on positive arg of SQRT
302 cdaq 1.14 * hszbeam is the intersection of the beam ray with the spectrometer
303 cdaq 1.2 * as measured along the z axis.
304 cdaq 1.14 if( sinhstheta .eq. 0.) then 305 hszbeam = 0.
306 cdaq 1.8 else
307 saw 1.16 hszbeam = sinhphi * ( -hsy_tar + gbeam_y * coshstheta) /
308 cdaq 1.14 $ sinhstheta 309 endif ! end test on sinhstheta=0
310 cdaq 1.10 * 311 * More kinematics 312 * 313 if(hsbeta.gt.0) then
314 saw 1.18 hsmass2 = (1/hsbeta*2 - 1)hsp**2
315 cdaq 1.10 else
316 cdaq 1.14 hsmass2 = 1.0e10
317 cdaq 1.10 endif 318
319 saw 1.16 hst = (gebeam - hsenergy)*2 320 $ - (gpbeam - hspcoshstheta)*2 - (hspsinhstheta)2 321 hsu = (gtarg_mass(gtarg_num) - hsenergy)2 - hsp2 322 if(hseloss.eq.0.)then 323 hseloss = gebeam - hsenergy 324 endif 325 hsq3 = sqrt(gpbeam2 + hsp*2 - 2gpbeamhspcoshstheta) 326 if(gpbeam.ne.0.and.hsq3.ne.0.)then 327 coshsthetaq = (gpbeam*2 - gpbeamhsp*coshstheta)/gpbeam/hsq3 328 endif 329 if(coshsthetaq.le.1.and.coshsthetaq.ge.-1.)then 330 hsthetaq = acos(coshsthetaq) 331 endif
332 cdaq 1.14 hsphiq = hsphi + tt
333 cdaq 1.10 hsbigq2 = -hst 334 hsx = hsbigq2/(2mass_nucleonhseloss)
335 saw 1.16 hsy = hseloss/gebeam 336 hsw2 = gtarg_mass(gtarg_num)*2 + 337 $ 2gtarg_mass(gtarg_num)*hseloss - hsbigq2
338 cdaq 1.10 if(hsw2.ge.0.0) then 339 hsw = sqrt(hsw2) 340 else 341 hsw = 0.0 342 endif 343
344 saw 1.17 * Calculate photon energy in GeV (E89-012):
345 saw 1.18 denom = hphoto_mtarget - hsenergy + hsp*coshstheta
346 saw 1.17 if (abs(denom).le.1.e-10) then 347 hsegamma = -1000.0 348 else 349 hsegamma = ( hsenergy * hphoto_mtarget
350 saw 1.18 & - 0.5(hphoto_mtarget2 + hpartmass*2
351 saw 1.17 & - hphoto_mrecoil**2) ) / denom
352 saw 1.18 endif 353 354 * Photon energy (assuming D2 target, Proton or Deut. detected). 355 denom = 1.87561 - sqrt(hsp2 + 0.938272) + hspcoshstheta 356 hsegamma_p= ( sqrt(hsp2+0.938272) 1.87561 357 & - 0.5(1.875612 + 0.938272 - 0.939572) )/denom 358 359 denom = 1.87561 - sqrt(hsp2 + 1.875612) + hspcoshstheta 360 hsegamma_d= ( sqrt(hsp2+1.875612) * 1.87561 361 & - 0.5(1.875612 + 1.875612 - 0.13497*2) )/denom 362
363 saw 1.17 c------------------------------------------------------------------------ 364
365 cdaq 1.15 *
366 cdaq 1.10 * Write raw timing information for fitting. 367 if(hdebugdumptof.ne.0) call h_dump_tof
368 cdaq 1.14 if(hdebugdumpcal.ne.0) call h_dump_cal
369 cdaq 1.10 * 370 * Calculate physics statistics and wire chamber efficencies.
371 cdaq 1.8 call h_physics_stat(ABORT,err) 372 ABORT= ierr.ne.0 .or. ABORT 373 IF(ABORT) THEN 374 call G_add_path(here,err) 375 ENDIF
376 cdaq 1.2 endif ! end test on zero tracks
377 saw 1.16 return 378 end

Analyzer/Replay: Mark Jones, Documents: Stephen Wood