Analyzer/HTRACKING/h_fpp_tracking_ajp.f - annotate

Return to h_fpp_tracking_ajp.f CVS log

Up to [HallC] / Analyzer / HTRACKING

1 puckett 1.1.2.1 subroutine h_fpp_tracking_ajp(dcset,abort,err) 2 3 implicit none 4 save 5 6 include 'gen_detectorids.par' 7 include 'gen_decode_common.cmn' 8 INCLUDE 'hms_data_structures.cmn' 9 INCLUDE 'hms_geometry.cmn' 10 INCLUDE 'hms_fpp_event.cmn' 11 INCLUDE 'hms_fpp_params.cmn' 12 INCLUDE 'hms_id_histid.cmn' 13 14 character14 here 15 parameter (here= 'h_fpp_tracking') 16 17 integer4 DCset ! set of FPP DCs we are working on 18 19 logical ABORT 20 character() err 21 22 puckett 1.1.2.1 logical*4 sufficient_hits, track_good, any_track, any_good, any_great 23 24 c define arrays for candidate tracks: 25
26 puckett 1.1.2.4 integer*4 nhitsrequired 27
28 puckett 1.1.2.1 integer4 ngoodtracks 29 integer4 ncandidates,i,j,k,ichamber,ilayer 30 integer4 nplanestemp,icluster,bestcandidate,best6,best5 31 integer4 itrack,track 32 33 integer4 goodtracks(h_fpp_max_candidates) 34 real4 simpletracks(h_fpp_max_candidates,6) 35 real4 fulltracks(h_fpp_max_candidates,6) 36 real4 thetatracks(h_fpp_max_candidates) ! store track polar angle relative to HMS
37 puckett 1.1.2.3 real4 sclosetracks(h_fpp_max_candidates) ! store track distance of closest approach relative to HMS or FPP1 38 real4 zclosetracks(h_fpp_max_candidates) ! store track z of closest approach relative to HMS or FPP1 39 real4 phitracks(h_fpp_max_candidates) 40 integer4 conetesttracks(h_fpp_max_candidates) ! store cone test of tracks relative to HMS or FPP1 41 integer*4 bestreftracks(h_fpp_max_candidates)
42 puckett 1.1.2.6 logical greattrack(h_fpp_max_candidates)
43 puckett 1.1.2.1 real*4 chi2, minchi2 44
45 puckett 1.1.2.3 c logical ambiguity(h_fpp_max_candidates) 46 logical firsttheta,firstsclose
47 puckett 1.1.2.1 logical firsttry,any6,anytheta ! at this stage of the tracking, apply a maximum theta cut on all candidate tracks 48 logical any6theta ! any tracks with six planes AND passing theta cut?
49 puckett 1.1.2.3 logical anyzclose,any6zclose ! any tracks with zclose passing prune tests?
50 puckett 1.1.2.6 logical anyconetest,any6conetest
51 puckett 1.1.2.1 52 integer4 hitcombos(h_fpp_max_candidates,h_fpp_n_dcinset,h_fpp_n_dclayers) 53 54 logical hitsintrack(h_fpp_max_candidates,h_fpp_n_dcinset,h_fpp_n_dclayers) 55 logical ontrack(h_fpp_n_dcinset,h_fpp_n_dclayers) ! keep track of whether a cluster ends up on a track 56 57 integer4 candidate_nplanes(h_fpp_max_candidates) 58 integer4 candidate_nhits(h_fpp_max_candidates) 59 60 real4 SimpleTrack(6), FullTrack(7) 61 62 integer4 BestClusters(H_FPP_N_DCINSET,H_FPP_N_DCLAYERS) 63 64 real4 HMStrack(4)
65 puckett 1.1.2.4 real*4 FPPtrack(4),newFPPtrack(4)
66 puckett 1.1.2.3 real*4 DCslopes(3),DCcoords(3), FPcoords(3), FPslopes(3)
67 puckett 1.1.2.1 68 real*4 PI 69 parameter(PI=3.141592653) 70
71 puckett 1.1.2.3 real4 theta,mintheta,phi,sclose,minsclose,zclose 72 real4 criterion,mincriterion
73 puckett 1.1.2.1 integer*4 conetest,bestref,jtrack
74 puckett 1.1.2.3 integer*4 iteration 75
76 puckett 1.1.2.5 logical zgood
77 puckett 1.1.2.3 logical goodhms 78 79 real*4 prob
80 puckett 1.1.2.1 81 ABORT= .FALSE. 82 err= ' ' 83 84 any_track = .false. 85 any_good = .false. 86 any_great = .false. 87 88 hfpp_n_tracks(dcset) = 0 89 90 call h_fpp_tracking_freehitcount(dcset,sufficient_hits) 91
92 puckett 1.1.2.4 iteration = 0
93 puckett 1.1.2.3 94 goodhms = abs(hsxp_tar).le..08.and.abs(hsyp_tar).le..04.and.abs(hsy_tar)<5..and. 95 $ abs(hsdelta)<10. 96
97 puckett 1.1.2.1 do while (sufficient_hits .and. (HFPP_N_tracks(DCset).lt.H_FPP_MAX_TRACKS))
98 puckett 1.1.2.3 99 c$$$ if(goodhms) then 100 c$$$ write(,) 'FPP tracking: chamber = ',dcset,' iteration = ',iteration 101 c$$$ endif 102 iteration = iteration + 1
103 puckett 1.1.2.1 104 ncandidates = 0 105 any6 = .false. 106 anytheta = .false. 107 any6theta = .false.
108 puckett 1.1.2.3 anyzclose = .false. 109 any6zclose = .false.
110 puckett 1.1.2.6 anyconetest = .false. 111 any6conetest = .false.
112 puckett 1.1.2.3
113 puckett 1.1.2.6 any_great = .false. 114
115 puckett 1.1.2.1 c the kind of loop we want here goes until we either run out of hit combos for a simple track, or 116 c until we exceed the maximum number of candidate tracks: 117 118 c initialize local tracking arrays: 119 120 do i=1,h_fpp_max_candidates 121 candidate_nplanes(i) = 0 122 candidate_nhits(i) = 0 123 goodtracks(i) = 0 124 do ichamber=1,h_fpp_n_dcinset 125 do ilayer=1,h_fpp_n_dclayers 126 hitcombos(i,ichamber,ilayer) = 0 127 hitsintrack(i,ichamber,ilayer) = .false. 128 enddo 129 enddo 130 enddo 131 132 c write(,) 'Start of tracking, FPP,ntracks=',dcset,hfpp_n_tracks(dcset) 133
134 puckett 1.1.2.4 nhitsrequired = 0 135 136 c write(,) 'before simple tracking, nhits required=',nhitsrequired 137 138 call h_fpp_tracking_simple_ajp(dcset,hitcombos,simpletracks,ncandidates,nhitsrequired,abort,err)
139 puckett 1.1.2.1 if (ABORT) then 140 call g_add_path(here,err) 141 return 142 endif
143 puckett 1.1.2.3
144 puckett 1.1.2.4 c$$$ write(,) 'after simple tracking, nhits required,ncandidate=', 145 c$$$ $ nhitsrequired,ncandidates 146 147 if(iteration.eq.1) 148 $ hfpp_n_simple(dcset,1) = ncandidates
149 puckett 1.1.2.1
150 puckett 1.1.2.3 c$$$ if(goodhms) then 151 c$$$ write(,) 'number of candidates after simple tracking=',ncandidates 152 c$$$ endif
153 puckett 1.1.2.1 c$$$ if (.true.) then 154 c$$$ if (int(SimpleTrack(6)).le.0) then 155 c$$$ if(hidFPP_trkrough(DCset,6).gt.0) call hf1(hidFPP_trkrough(DCset,6),0.,1.) !Nraw 156 c$$$ else 157 c$$$ if(hidFPP_trkrough(DCset,1).gt.0) call hf1(hidFPP_trkrough(DCset,1),SimpleTrack(1),1.) !mx 158 c$$$ if(hidFPP_trkrough(DCset,2).gt.0) call hf1(hidFPP_trkrough(DCset,2),SimpleTrack(2),1.) !bx 159 c$$$ if(hidFPP_trkrough(DCset,3).gt.0) call hf1(hidFPP_trkrough(DCset,3),SimpleTrack(3),1.) !my 160 c$$$ if(hidFPP_trkrough(DCset,4).gt.0) call hf1(hidFPP_trkrough(DCset,4),SimpleTrack(4),1.) !by 161 c$$$ if(hidFPP_trkrough(DCset,5).gt.0) call hf1(hidFPP_trkrough(DCset,5),SimpleTrack(5),1.) !chi2 162 c$$$ if(hidFPP_trkrough(DCset,6).gt.0) call hf1(hidFPP_trkrough(DCset,6),SimpleTrack(6),1.) !Nraw 163 c$$$ endif 164 c$$$ endif 165 166 * * quit trying to make more tracks if we are out of hits 167 * * or if we couldnt make a good one now 168 c if (int(SimpleTrack(6)).le.0) exit 169 c if (int(SimpleTrack(5)).eq.H_FPP_BAD_CHI2) exit 170 171 c any_track = .true. 172 173 if(ncandidates.eq.0) exit 174 puckett 1.1.2.1 175 ngoodtracks = 0 176
177 puckett 1.1.2.4 134 do i=1,ncandidates
178 puckett 1.1.2.1 FullTrack(1) = H_FPP_BAD_COORD ! mx 179 FullTrack(2) = H_FPP_BAD_COORD ! bx 180 FullTrack(3) = H_FPP_BAD_COORD ! my 181 FullTrack(4) = H_FPP_BAD_COORD ! by 182 FullTrack(5) = H_FPP_BAD_CHI2 183 FullTrack(6) = 0. ! number of hits 184 FullTrack(7) = 0. ! number of planes 185 c write(,) 'candidate ',i,' simple chi2=',simpletracks(i,5) 186 187 do j=1,6 188 simpletrack(j) = simpletracks(i,j) 189 c write(,) 'i,simpletrack(i)=',j,simpletrack(j) 190 enddo 191 192 do ichamber=1,h_fpp_n_dcinset 193 do ilayer=1,h_fpp_n_dclayers 194 bestclusters(ichamber,ilayer) = hitcombos(i,ichamber,ilayer) 195 c write(,) 'chbr,pln,clstr=',ichamber,ilayer,bestclusters(ichamber,ilayer) 196 enddo 197 enddo 198 199 puckett 1.1.2.1 call h_fpp_tracking_drifttrack_ajp(dcset,simpletrack,bestclusters,
200 puckett 1.1.2.4 $ ontrack,track_good,fulltrack,nhitsrequired,1,abort,err)
201 puckett 1.1.2.1 202 203 c write(,) 'track_good = ',track_good 204 if (ABORT) then 205 call g_add_path(here,err) 206 return 207 endif 208
209 puckett 1.1.2.3 if(track_good) then ! tracking with drift resulted in a reasonable track: 210 211 * calculate everything so we can check its quality/value against quantities other than just chi2: 212 * transform from FPP to focal plane coordinate system: 213 * slopes first: 214 dcslopes(1) = fulltrack(1) 215 dcslopes(2) = fulltrack(3) 216 dcslopes(3) = 1.0 217 218 call h_fpp_dc2fp(dcset,.true.,dcslopes,fpslopes)
219 puckett 1.1.2.1
220 puckett 1.1.2.3 * Next, transform coordinates from FPP to focal plane system: 221 dccoords(1) = fulltrack(2) 222 dccoords(2) = fulltrack(4) 223 dccoords(3) = 0.0 224 225 call h_fpp_dc2fp(dcset,.false.,dccoords,fpcoords) 226 227 fpcoords(1) = fpcoords(1) - fpslopes(1)fpcoords(3) 228 fpcoords(2) = fpcoords(2) - fpslopes(2)fpcoords(3)
229 puckett 1.1.2.4 230 fpptrack(1) = fpcoords(1) 231 fpptrack(2) = fpcoords(2) 232 fpptrack(3) = fpslopes(1) 233 fpptrack(4) = fpslopes(2) 234 235 call h_fpp_align(dcset,fpptrack,newfpptrack)
236 puckett 1.1.2.3 * initialize bestref to zero: 237 bestref = 0 238 * initialize "hmstrack" to the actual hms track: 239 hmstrack(1) = hsxp_fp 240 hmstrack(2) = hsx_fp 241 hmstrack(3) = hsyp_fp 242 hmstrack(4) = hsy_fp 243 * initialize "fpptrack" to the current track:
244 puckett 1.1.2.4 c$$$ fpptrack(1) = fpslopes(1) 245 c$$$ fpptrack(2) = fpcoords(1) 246 c$$$ fpptrack(3) = fpslopes(2) 247 c$$$ fpptrack(4) = fpcoords(2) 248 fpptrack(1) = newfpptrack(3) ! dx/dz 249 fpptrack(2) = newfpptrack(1) ! x 250 fpptrack(3) = newfpptrack(4) ! dy/dz 251 fpptrack(4) = newfpptrack(2) ! y 252
253 puckett 1.1.2.3 * calculate closest approach parameters relative to hms track: 254 call h_fpp_closest(hmstrack,fpptrack,sclose,zclose) 255 call h_fpp_relative_angles(hmstrack(1),hmstrack(3),fpptrack(1),fpptrack(3),theta,phi)
256 puckett 1.1.2.4 * in FPP2, we use the "best" track from FPP1 as a reference track:
257 puckett 1.1.2.6 mintheta = theta 258 bestref = 0
259 puckett 1.1.2.4 if(dcset.eq.2.and.hfpp_n_tracks(1).gt.0) then
260 puckett 1.1.2.6 do jtrack = 1, hfpp_n_tracks(1) 261 hmstrack(1) = hfpp_track_dx(1,jtrack) 262 hmstrack(2) = hfpp_track_x(1,jtrack) 263 hmstrack(3) = hfpp_track_dy(1,jtrack) 264 hmstrack(4) = hfpp_track_y(1,jtrack) 265 call h_fpp_closest(hmstrack,fpptrack,sclose,zclose) 266 call h_fpp_relative_angles(hmstrack(1),hmstrack(3),fpptrack(1),fpptrack(3),theta,phi) 267 if(theta.lt.mintheta) then 268 bestref = jtrack 269 bestreftracks(i) = bestref 270 mintheta = theta 271 endif 272 enddo 273 c now, calculate everything one more time using bestref: 274 if(bestref.gt.0) then 275 hmstrack(1) = hfpp_track_dx(1,bestref) 276 hmstrack(2) = hfpp_track_x(1,bestref) 277 hmstrack(3) = hfpp_track_dy(1,bestref) 278 hmstrack(4) = hfpp_track_y(1,bestref) 279 else ! revert to HMS track: 280 hmstrack(1) = hsxp_fp 281 puckett 1.1.2.6 hmstrack(2) = hsx_fp 282 hmstrack(3) = hsyp_fp 283 hmstrack(4) = hsy_fp 284 endif 285
286 puckett 1.1.2.3 call h_fpp_closest(hmstrack,fpptrack,sclose,zclose)
287 puckett 1.1.2.4 call h_fpp_relative_angles(hmstrack(1),hmstrack(3),fpptrack(1),fpptrack(3),theta,phi)
288 puckett 1.1.2.6
289 puckett 1.1.2.3 endif 290
291 puckett 1.1.2.6 bestreftracks(i) = bestref 292
293 puckett 1.1.2.3 conetest = 1
294 puckett 1.1.2.1
295 puckett 1.1.2.3 call h_fpp_conetest(hmstrack,dcset,zclose,theta,conetest) 296
297 puckett 1.1.2.1 thetatracks(i) = theta
298 puckett 1.1.2.3 phitracks(i) = phi
299 puckett 1.1.2.1
300 puckett 1.1.2.3 sclosetracks(i) = sclose 301 zclosetracks(i) = zclose 302 303 conetesttracks(i) = conetest 304 305 if(theta.le.hfpp_prune_thetamax(dcset)*PI/180.0) 306 $ anytheta = .true. 307
308 puckett 1.1.2.1 do j=1,6 309 fulltracks(i,j) = fulltrack(j) 310 enddo 311 312 nplanestemp = int(fulltrack(7)) 313 314 candidate_nplanes(i) = nplanestemp 315 candidate_nhits(i) = int(fulltrack(6)) 316 317 ngoodtracks = ngoodtracks + 1 318 goodtracks(ngoodtracks) = i 319
320 puckett 1.1.2.6 greattrack(i) = ( conetest.eq.1.and. 321 $ sclose.le.hfpp_prune_sclose(dcset) 322 $ .and.zclose.ge.hfpp_prune_zclose(2(dcset-1)+1)- 323 $ hfpp_prune_zslop(dcset)/tan(theta).and.zclose.le. 324 $ hfpp_prune_zclose(2(dcset-1)+2) + 325 $ hfpp_prune_zslop(dcset)/tan(theta) .and. 326 $ nplanestemp.eq.6 ) 327 328 if(greattrack(i)) any_great = .true. 329
330 puckett 1.1.2.1 if(nplanestemp.eq.h_fpp_n_dcinset*h_fpp_n_dclayers) 331 $ any6 = .true. 332
333 puckett 1.1.2.6 if(conetest.eq.1) then 334 anyconetest = .true. 335 if(nplanestemp.eq.h_fpp_n_dcinset*h_fpp_n_dclayers) 336 $ then 337 any6conetest = .true. 338 endif 339 endif 340
341 puckett 1.1.2.5 if(zclose.ge.hfpp_prune_zclose(2(dcset-1)+1) - 342 $ hfpp_prune_zslop(dcset)/tan(theta).and. 343 $ zclose.le.hfpp_prune_zclose(2(dcset-1)+2) + 344 $ hfpp_prune_zslop(dcset)/tan(theta)) then
345 puckett 1.1.2.3 anyzclose = .true. 346 if(nplanestemp.eq.h_fpp_n_dcinset*h_fpp_n_dclayers) 347 $ then 348 any6zclose = .true. 349 endif 350 endif 351
352 puckett 1.1.2.1 if(theta.le.hfpp_prune_thetamax(dcset)PI/180.0 353 $ .and.nplanestemp.eq.h_fpp_n_dcinseth_fpp_n_dclayers) 354 $ any6theta = .true. 355 356 do ichamber=1,h_fpp_n_dcinset 357 do ilayer = 1,h_fpp_n_dclayers 358 if(ontrack(ichamber,ilayer)) 359 $ hitsintrack(i,ichamber,ilayer) = .true. 360 enddo 361 enddo
362 puckett 1.1.2.4
363 puckett 1.1.2.1 endif 364 enddo
365 puckett 1.1.2.6 366 c$$$ if(ngoodtracks.gt.hfpp_n_simple(dcset,2)) then
367 puckett 1.1.2.4 c$$$ hfpp_n_simple(dcset,2) = ngoodtracks 368 c$$$ endif 369 370 if(iteration.eq.1) then 371 hfpp_n_simple(dcset,2) = ngoodtracks 372 endif
373 puckett 1.1.2.6
374 puckett 1.1.2.3 c if the left-right fixing routine is turned on and we don't prune on the number of planes 375 c on the track, then don't require six planes:
376 puckett 1.1.2.6
377 puckett 1.1.2.3 if(hfpp_prune_nplanes.eq.0.and.hfppfixleftright.gt.0) then 378 any6 = .false. 379 endif 380 381 c$$$ if(goodhms) then 382 c$$$ write(,) 383 c$$$ $ 'number of good tracks after tracking with drift=',ngoodtracks 384 c$$$ endif
385 puckett 1.1.2.4 386 c write(,) 'ngood drift tracks = ',ngoodtracks 387 388 if(ngoodtracks.eq.0) then 389 if(nhitsrequired.gt.hfpp_minsethits) then 390 nhitsrequired = nhitsrequired - 1 391 goto 134 392 else 393 exit 394 endif 395 endif
396 puckett 1.1.2.1 397 bestcandidate = 0 398 firsttry = .true.
399 puckett 1.1.2.3 firsttheta = .true. 400 firstsclose = .true. 401 * get minimum chi2, sclose, and theta, always checking for existence of six-hit tracks first:
402 puckett 1.1.2.4 do i=1,ngoodtracks 403 track=goodtracks(i) 404 chi2 = fulltracks(track,5) 405 sclose = sclosetracks(track) 406 theta = thetatracks(track) 407 if(candidate_nplanes(track).eq.6.or..not.any6) then 408 if(firsttry) then 409 firsttry = .false. 410 mintheta = theta 411 minchi2 = chi2 412 minsclose = sclose 413 else 414 if(theta.lt.mintheta) mintheta = theta 415 if(chi2.lt.minchi2) minchi2 = chi2 416 if(sclose.lt.minsclose) minsclose = sclose 417 endif 418 endif 419 enddo
420 puckett 1.1.2.3 firsttry = .true. 421 do i=1,ngoodtracks ! here is where we try to pick the best track: 422
423 puckett 1.1.2.1 track = goodtracks(i) 424 425 chi2 = fulltracks(track,5)
426 puckett 1.1.2.3 sclose = sclosetracks(track) 427 theta = thetatracks(track) 428 zclose = zclosetracks(track) 429 conetest = conetesttracks(track) 430 bestref = bestreftracks(track) 431 phi = phitracks(track) 432 433 c if(ngoodtracks.gt.1) then 434 c if(.not.any6.or.candidate_nplanes(track).eq.6) then 435 c$$$ if(goodhms) then 436 c$$$ write(,) 'track ',track,' theta=',theta180.0/PI, 437 c$$$ $ ' chi2=',chi2, 438 c$$$ $ ' sclose=',sclose,' zclose=',zclose,' phi=',phi180.0/PI, 439 c$$$ $ ' nplanes=',candidate_nplanes(track), 440 c$$$ $ ' nhits=',candidate_nhits(track), 441 c$$$ $ ' conetest=',conetest 442 c$$$ if(dcset.eq.2) write(,) 'bestref=',bestref 443 c$$$ endif 444 c endif 445 446 c try chi2 + sclose^2/sigmasclose^2 447 puckett 1.1.2.3 c criterion = chi2 + (sclosehfpp_sclose_weight(dcset))*2
448 puckett 1.1.2.1
449 puckett 1.1.2.3 c criterion = chi2 + minchi2(hfpp_sclose_weight(dcset)(sclose/minsclose)*2 + 450 c $ hfpp_theta_weight(dcset)prob(dcset,mintheta)/prob(dcset,theta)) 451 452 c if(dcset.eq.2.and.bestref.gt.0) then 453 c criterion = chi2 + minchi2(hfpp_sclose_weight(1)(sclose/minsclose)*2 + 454 c $ hfpp_theta_weight(1)prob(1,mintheta)/prob(1,theta)) 455 c endif 456
457 puckett 1.1.2.5 zgood = zclose.ge.hfpp_prune_zclose(2(dcset-1)+1) - 458 $ hfpp_prune_zslop(dcset)/tan(theta) .and. 459 $ zclose.le.hfpp_prune_zclose(2(dcset-1)+2) + 460 $ hfpp_prune_zslop(dcset)/tan(theta) 461
462 puckett 1.1.2.3 c try picking smallest theta HERE instead:
463 puckett 1.1.2.4 if(hselectfpptrackprune.eq.1) then 464 criterion = theta 465 else if(hselectfpptrackprune.eq.2) then 466 criterion = chi2 + minchi2 * sclose**2 467 else if(hselectfpptrackprune.eq.3) then 468 criterion = chi2 469 else if(hselectfpptrackprune.eq.4) then 470 criterion = sclose 471 else 472 criterion = chi2 473 endif
474 puckett 1.1.2.3
475 puckett 1.1.2.6 c if(candidate_nplanes(track).eq.h_fpp_n_dcinset* 476 c $ h_fpp_n_dclayers.or..not.any6) then 477 if( greattrack(track).or. .not. any_great) then 478 if(candidate_nplanes(track).eq.6.or..not.any6) then
479 puckett 1.1.2.5 if(firsttry.or.criterion.lt.mincriterion) then 480 mincriterion = criterion 481 firsttry = .false. 482 bestcandidate = track 483 endif
484 puckett 1.1.2.1 endif 485 endif 486 enddo 487
488 puckett 1.1.2.3 c$$$ if(goodhms) then 489 c$$$ write(,) 'chosen track=',bestcandidate 490 c$$$ endif 491 if(hidFPP_trkrough(DCset,1).gt.0) call hf1(hidFPP_trkrough(DCset,1),SimpleTracks(bestcandidate,1),1.) !mx 492 if(hidFPP_trkrough(DCset,2).gt.0) call hf1(hidFPP_trkrough(DCset,2),SimpleTracks(bestcandidate,2),1.) !bx 493 if(hidFPP_trkrough(DCset,3).gt.0) call hf1(hidFPP_trkrough(DCset,3),SimpleTracks(bestcandidate,3),1.) !my 494 if(hidFPP_trkrough(DCset,4).gt.0) call hf1(hidFPP_trkrough(DCset,4),SimpleTracks(bestcandidate,4),1.) !by 495 if(hidFPP_trkrough(DCset,5).gt.0) call hf1(hidFPP_trkrough(DCset,5),SimpleTracks(bestcandidate,5),1.) !chi2 496 if(hidFPP_trkrough(DCset,6).gt.0) call hf1(hidFPP_trkrough(DCset,6),SimpleTracks(bestcandidate,6),1.) !Nraw 497 498 if(hidFPP_roughchi2vsnhit(dcset).gt.0) call hf2(hidFPP_roughchi2vsnhit(dcset),simpletracks(bestcandidate,6), 499 $ simpletracks(bestcandidate,5),1.) 500
501 puckett 1.1.2.1 do i=1,ncandidates 502 do ichamber=1,h_fpp_n_dcinset 503 do ilayer=1,h_fpp_n_dclayers 504 icluster = hitcombos(i,ichamber,ilayer) 505 if(icluster.gt.0) then ! mark all candidate hits as unused 506 hfpp_clusterintrack(dcset,ichamber,ilayer,icluster) = 0 507 endif 508 enddo 509 enddo 510 enddo
511 puckett 1.1.2.3
512 puckett 1.1.2.1 c write(,) 'best candidate track, chi2=',bestcandidate,fulltracks(bestcandidate,5)
513 puckett 1.1.2.3 c write(,) 'nplanes,nhits=',candidate_nplanes(bestcandidate),
514 puckett 1.1.2.1 c $ candidate_nhits(bestcandidate) 515
516 puckett 1.1.2.3 c$$$ do i=1,4 517 c$$$ write(,) 'i,besttrack(i)=',i,fulltracks(bestcandidate,i) 518 c$$$ enddo
519 puckett 1.1.2.1 520 if(bestcandidate.eq.0) exit 521 522 c Now that we have found the best candidate track, add it to the good track array: 523 524 itrack = hfpp_n_tracks(dcset) + 1 525 526 if(itrack.le.h_fpp_max_tracks) then
527 puckett 1.1.2.3
528 puckett 1.1.2.1 hfpp_n_tracks(dcset) = itrack 529 530 hfpp_track_nlayers(dcset,itrack) = candidate_nplanes(bestcandidate) 531 c Mark the hits in the best candidate track as used! 532 do ichamber = 1,h_fpp_n_dcinset 533 do ilayer = 1,h_fpp_n_dclayers 534 icluster = hitcombos(bestcandidate,ichamber,ilayer) 535 if(icluster.gt.0.and. 536 $ hitsintrack(bestcandidate,ichamber,ilayer) ) then 537 hfpp_clusterintrack(dcset,ichamber,ilayer,icluster) = itrack 538 endif 539 hfpp_trackcluster(dcset,ichamber,ilayer,itrack) = icluster 540 enddo 541 enddo 542 543 do j=1,4 544 hfpp_track_fine(dcset,itrack,j) = fulltracks(bestcandidate,j) 545 enddo 546 547 hfpp_track_chi2(dcset,itrack) = fulltracks(bestcandidate,5) 548 hfpp_track_nhits(dcset,itrack) = int(fulltracks(bestcandidate,6)) 549 puckett 1.1.2.1 550 do j=1,6 551 hfpp_track_rough(dcset,itrack,j) = simpletracks(bestcandidate,j) 552 enddo 553 554 hfpp_track_uniq(dcset,itrack) = .true. 555 * transform slopes from local FPP coordinates to HMS fp coordinates: 556 dccoords(1) = fulltracks(bestcandidate,1) ! dx/dz 557 dccoords(2) = fulltracks(bestcandidate,3) ! dy/dz 558 dccoords(3) = 1.0 ! dz/dz 559 560 call h_fpp_dc2fp(dcset,.true.,dccoords,fpcoords) 561
562 puckett 1.1.2.4 c$$$ hfpp_track_dx(dcset,itrack) = fpcoords(1) 563 c$$$ hfpp_track_dy(dcset,itrack) = fpcoords(2) 564 565 fpptrack(3) = fpcoords(1) ! x' in fp coords 566 fpptrack(4) = fpcoords(2) ! y' in fp coords
567 puckett 1.1.2.1 568 * now transform coordinates from local FPP system to HMS fp coords: 569 570 dccoords(1) = fulltracks(bestcandidate,2) ! x 571 dccoords(2) = fulltracks(bestcandidate,4) ! y 572 dccoords(3) = 0.0 573 574 call h_fpp_dc2fp(dcset,.false.,dccoords,fpcoords) 575
576 puckett 1.1.2.4 c$$$ hfpp_track_x(dcset,itrack) = fpcoords(1) - fpcoords(3)* 577 c$$$ $ hfpp_track_dx(dcset,itrack) 578 c$$$ hfpp_track_y(dcset,itrack) = fpcoords(2) - fpcoords(3)* 579 c$$$ $ hfpp_track_dy(dcset,itrack) 580 581 fpptrack(1) = fpcoords(1) - fpcoords(3) * fpptrack(3) ! x in fp coords at z=0 582 fpptrack(2) = fpcoords(2) - fpcoords(3) * fpptrack(4) ! y in fp coords at z=0 583 584 call h_fpp_align(dcset,fpptrack,newfpptrack) 585 586 hfpp_track_x(dcset,itrack) = newfpptrack(1) 587 hfpp_track_y(dcset,itrack) = newfpptrack(2) 588 hfpp_track_dx(dcset,itrack) = newfpptrack(3) 589 hfpp_track_dy(dcset,itrack) = newfpptrack(4)
590 puckett 1.1.2.1 591 * calculate relative scattering angles theta and phi: 592 593 call h_fpp_relative_angles(hsxp_fp,hsyp_fp,hfpp_track_dx(dcset,itrack), 594 $ hfpp_track_dy(dcset,itrack),theta,phi) 595 596 hfpp_track_theta(dcset,itrack) = theta 597 hfpp_track_phi(dcset,itrack) = phi 598
599 puckett 1.1.2.4 * calculate closest approach parameters with respect to HMS, REGARDLESS of which FPP!
600 puckett 1.1.2.1 601 hmstrack(1) = hsxp_fp 602 hmstrack(2) = hsx_fp 603 hmstrack(3) = hsyp_fp 604 hmstrack(4) = hsy_fp 605 606 fpptrack(1) = hfpp_track_dx(dcset,itrack) 607 fpptrack(2) = hfpp_track_x(dcset,itrack) 608 fpptrack(3) = hfpp_track_dy(dcset,itrack) 609 fpptrack(4) = hfpp_track_y(dcset,itrack) 610 611 call h_fpp_closest(hmstrack,fpptrack,sclose,zclose) 612 613 hfpp_track_sclose(dcset,itrack) = sclose 614 hfpp_track_zclose(dcset,itrack) = zclose 615 616 conetest = 1 617 618 call h_fpp_conetest(hmstrack,dcset,zclose,theta,conetest) 619 620 hfpp_track_conetest(dcset,itrack) = conetest 621 puckett 1.1.2.1 622 if(dcset.eq.2) then 623
624 puckett 1.1.2.4 c firsttry=.true.
625 puckett 1.1.2.1 bestref = 0
626 puckett 1.1.2.4
627 puckett 1.1.2.6 c mintheta = hfpp_track_theta(dcset,itrack)
628 puckett 1.1.2.4 c do jtrack=1,hfpp_n_tracks(1)
629 puckett 1.1.2.6 c how this works: always figure out the "best reference" track for this FPP2 track in FPP1: 630 c after that, make a subjective judgement, based on theta, of whether or not to compare this track 631 c to FPP1 or the HMS:
632 puckett 1.1.2.4 if(hfpp_n_tracks(1).gt.0) then
633 puckett 1.1.2.6 do jtrack = 1, hfpp_n_tracks(1) 634 call h_fpp_relative_angles(hfpp_track_dx(1,jtrack), 635 $ hfpp_track_dy(1,jtrack), 636 $ hfpp_track_dx(dcset,itrack), 637 $ hfpp_track_dy(dcset,itrack), 638 $ theta,phi) 639 if(jtrack.eq.1.or.theta.lt.mintheta) then 640 bestref = jtrack 641 mintheta = theta 642 endif 643 enddo
644 puckett 1.1.2.4 endif
645 puckett 1.1.2.1 646 if(bestref.gt.0) then 647 call h_fpp_relative_angles(hfpp_track_dx(1,bestref), 648 $ hfpp_track_dy(1,bestref), 649 $ hfpp_track_dx(dcset,itrack), 650 $ hfpp_track_dy(dcset,itrack),theta,phi) 651 hfpp_track_theta(dcset+1,itrack) = theta 652 hfpp_track_phi(dcset+1,itrack) = phi 653 654 hmstrack(1) = hfpp_track_dx(1,bestref) 655 hmstrack(2) = hfpp_track_x(1,bestref) 656 hmstrack(3) = hfpp_track_dy(1,bestref) 657 hmstrack(4) = hfpp_track_y(1,bestref) 658 659 fpptrack(1) = hfpp_track_dx(dcset,itrack) 660 fpptrack(2) = hfpp_track_x(dcset,itrack) 661 fpptrack(3) = hfpp_track_dy(dcset,itrack) 662 fpptrack(4) = hfpp_track_y(dcset,itrack) 663 664 call h_fpp_closest(hmstrack,fpptrack,sclose,zclose) 665 666 puckett 1.1.2.1 hfpp_track_sclose(dcset+1,itrack) = sclose 667 hfpp_track_zclose(dcset+1,itrack) = zclose 668 669 conetest = 1 670 671 call h_fpp_conetest(hmstrack,dcset,zclose,theta,conetest) 672 673 hfpp_track_conetest(dcset+1,itrack) = conetest 674 endif
675 puckett 1.1.2.6 676 hfpp2_best_reference(itrack) = bestref 677 678 if(bestref.gt.0.and.hfpp_track_theta(dcset,itrack).lt. 679 $ hfpp_track_theta(dcset+1,itrack)) hfpp2_best_reference(itrack) = 0 680
681 puckett 1.1.2.1 endif
682 puckett 1.1.2.4
683 puckett 1.1.2.3 if(hfppfixleftright.gt.0) 684 $ call h_fpp_check_leftright(dcset,itrack) 685 686 endif ! end filling track common block variables
687 puckett 1.1.2.1 c NOW we need to call freehitcount and start over with the remaining free hits: 688 689 call h_fpp_tracking_freehitcount(dcset,sufficient_hits) 690 691 enddo 692 693 return 694 end 695 696 subroutine h_fpp_tracking_simple_ajp(dcset,hitcombos,simpletracks,
697 puckett 1.1.2.4 $ ncandidates,nhitsrequired,abort,err)
698 puckett 1.1.2.1 699 implicit none 700 save 701 702 INCLUDE 'hms_data_structures.cmn' 703 INCLUDE 'hms_geometry.cmn' 704 include 'gen_detectorids.par' 705 include 'gen_decode_common.cmn' 706 INCLUDE 'hms_fpp_event.cmn' 707 INCLUDE 'hms_fpp_params.cmn' 708 709 character21 here 710 parameter (here= 'h_fpp_tracking_simple') 711 712 integer4 dcset 713 integer4 hitcombos(h_fpp_max_candidates,h_fpp_n_dcinset,h_fpp_n_dclayers) 714 real4 simpletracks(h_fpp_max_candidates,6) 715 716 integer4 HitCluster(H_FPP_N_DCINSET,H_FPP_N_DCLAYERS) 717 718 integer4 ncandidates,npoints 719 puckett 1.1.2.1 integer4 nhitsrequired,iterations,ichamber,ilayer,i,j 720 integer4 nhitsintrack
721 puckett 1.1.2.3 integer*4 n3hit,n2hit ! tabulate number of three-hit and two hit clusters.
722 puckett 1.1.2.1 723 logical abort 724 character() err 725 726 real4 temptrack(5) ! not including hit count 727 728 abort = .false. 729 err= ' ' 730 731 do i=1,h_fpp_max_candidates 732 do ichamber=1,h_fpp_n_dcinset 733 do ilayer=1,h_fpp_n_dclayers 734 hitcombos(i,ichamber,ilayer) = 0 735 enddo 736 enddo 737 do j=1,5 738 simpletracks(i,j) = h_fpp_bad_coord 739 enddo 740 simpletracks(i,6) = 0.0 741 enddo 742 743 puckett 1.1.2.1 nhitsrequired = h_fpp_n_planes + 1 744 iterations = 0 745 746 do while(nhitsrequired .ge. hfpp_minsethits.and.ncandidates.lt. 747 $ h_fpp_max_candidates) 748 749 nhitsintrack = 0 750 call next combo with large nhitsrequired to get the first useful combo as well as the max 751 * number of hits available 752 ncandidates = 0 753
754 puckett 1.1.2.3 c$$$ do ichamber=1,h_fpp_n_dcinset 755 c$$$ do ilayer=1,h_fpp_n_dclayers 756 c$$$ hitcluster(ichamber,ilayer) = 0 757 c$$$ enddo 758 c$$$ enddo 759
760 puckett 1.1.2.1 call h_fpp_tracking_nexthitcombo(dcset,nhitsrequired,nhitsintrack,hitcluster) 761 762 * keep comparing permutations until all possibilities are tried: 763 764 do while(nhitsintrack.gt.0.and.ncandidates.lt.h_fpp_max_candidates) 765 766 iterations = iterations + 1 767
768 puckett 1.1.2.3 if(iterations.gt.hfpp_maxcombos) then 769 write(,) 'WARNING: max FPP hit combos reached.' 770 write(,) 'consider raising the limit or using ' 771 write(,) 'tighter raw timing cuts' 772 exit 773 endif
774 puckett 1.1.2.1 775 call h_fpp_fit_simple(dcset,hitcluster,npoints,temptrack,abort,err) 776 777 if (ABORT) then 778 call g_add_path(here,err) 779 return 780 endif 781
782 puckett 1.1.2.3 n2hit = 0 783 n3hit = 0 784 785 do ichamber=1,h_fpp_n_dcinset 786 do ilayer=1,h_fpp_n_dclayers 787 if(hitcluster(ichamber,ilayer).gt.0) then 788 if(hfpp_nhitsincluster(dcset,ichamber,ilayer, 789 $ hitcluster(ichamber,ilayer)).eq.2) n2hit = n2hit + 1 790 if(hfpp_nhitsincluster(dcset,ichamber,ilayer,hitcluster(ichamber,ilayer)) 791 $ .eq.3) n3hit = n3hit + 1 792 endif 793 enddo 794 enddo 795 * in contrast to Frank's algorithm, here we keep track of any hit combos which pass the "reasonable chi2" 796 * criterion, and we don't choose a track until we look at the drift based tracking. 797
798 puckett 1.1.2.1 if(temptrack(5).ge.0.0.and.
799 puckett 1.1.2.3 $ temptrack(5).le.hfpp_aok_chi2+ 800 $ (12.float(n2hit)+48.float(n3hit))/float(nhitsrequired-4)/ 801 $ float(nhitsrequired-4 + n2hit + 2*n3hit) 802 $ ) then ! add a new candidate track to the array:
803 puckett 1.1.2.1 ncandidates = ncandidates + 1 804 do ichamber=1,h_fpp_n_dcinset 805 do ilayer=1,h_fpp_n_dclayers 806 hitcombos(ncandidates,ichamber,ilayer) = 807 $ hitcluster(ichamber,ilayer) 808 enddo 809 enddo 810 811 do j=1,5 812 simpletracks(ncandidates,j) = temptrack(j) 813 enddo 814 simpletracks(ncandidates,6) = float(npoints) 815 endif 816 817 call h_fpp_tracking_nexthitcombo(dcset,nhitsrequired,nhitsintrack,hitcluster) 818 819 enddo 820 821 if(iterations.gt.hfpp_maxcombos) then 822 ncandidates = 0 823 exit 824 puckett 1.1.2.1 endif 825 826 if(ncandidates.gt.0) exit 827 828 nhitsrequired = nhitsrequired - 1 829 830 enddo 831 832 return 833 end
834 puckett 1.1.2.3 835 real4 function prob(chbr,angle) 836 837 integer4 chbr 838 real4 angle 839 c approximate cumulative probability that a track scatters in fpp with 840 c an angle theta < angle 841 c prob = min(1.0e-9,1.0 - .1662angle*(-0.2202)exp(-4.031angle0.981)) 842 843 if(chbr.eq.1) prob = exp(-2.946angle) 844 else prob = exp(-3.485angle) 845 846 return 847 end 848 849 subroutine h_fpp_check_leftright(ifpp,itrack) 850 851 implicit none 852 save 853 854 include 'gen_detectorids.par' 855 puckett 1.1.2.3 include 'gen_decode_common.cmn' 856 include 'gen_event_info.cmn' 857 INCLUDE 'hms_data_structures.cmn' 858 INCLUDE 'hms_geometry.cmn' 859 INCLUDE 'hms_fpp_event.cmn' 860 INCLUDE 'hms_fpp_params.cmn' 861 INCLUDE 'hms_id_histid.cmn' 862 863 c fortunately, this is made easier by the fact that we know which clusters are on a track and 864 c we know which hits within a cluster are on a track from their drift distance variable, so 865 c we can revisit the left-right combo here: 866 867 integer4 nambig,ambigcombos(64) 868 integer4 sign,oldsign,j 869 integer4 npoints,ncombos,point,ipoint,icombo,combo 870 integer4 ifpp,itrack,track,jtrack 871 integer4 ichamber,ilayer,iwire,wire 872 integer4 icluster,cluster,ihit,hit 873 integer4 oldbestcombo,newbestcombo 874 real4 combochi2(64) 875 real4 combotheta(64) 876 puckett 1.1.2.3 real4 combophi(64) 877 real4 combosclose(64) 878 real4 combozclose(64) 879 integer4 comboconetest(64) 880 logical ambig(64),anyambig,first 881 882 real4 chi2,theta,phi,sclose,zclose 883 real4 oldchi2,oldtheta,oldphi,oldsclose,oldzclose 884 885 integer4 oldconetest,conetest 886 887 real4 plusminustest(h_fpp_max_fitpoints) 888 real4 plusminusbest(h_fpp_max_fitpoints) 889 890 integer4 chambers(h_fpp_max_fitpoints) 891 integer4 layers(h_fpp_max_fitpoints) 892 integer4 wires(h_fpp_max_fitpoints) 893 894 real4 points(h_fpp_max_fitpoints,2) ! w and z 895 real4 drifts(h_fpp_max_fitpoints) ! drift distance 896 real4 projects(h_fpp_max_fitpoints,2) ! Px and Py coefficients 897 puckett 1.1.2.3 real4 sigma2s(h_fpp_max_fitpoints) ! sigma^2 for chi2 calc. 898 real4 hitpos(h_fpp_max_fitpoints,2) ! arrays for track fitting routine 899 900 real4 dccoords(3),fpcoords(3)
901 puckett 1.1.2.4 real*4 reftrack(4),fpptrack(4),hmstrack(4),newfpptrack(4)
902 puckett 1.1.2.3 903 real4 testtrack(5) 904 905 real4 minsclose,mintheta 906 real4 minstest,minthetatest,maxambig,maxdtheta,maxthetadiff 907 real4 wtrack,xtrack,ytrack,xptrack,yptrack,Px,Py,z 908 real4 zmin,zmax 909 910 external jbit ! cernlib bit routine 911 external jibset 912 integer4 jbit 913 integer4 jibset ! Declare to help f2c 914 915 hfpp_track_nambig(ifpp,itrack) = 0 916 917 if(itrack.le.0.or.itrack.gt.hfpp_n_tracks(ifpp).or. 918 $ hfpp_track_nhits(ifpp,itrack).ne. 919 $ hfpp_track_nlayers(ifpp,itrack)) return 920 921 initialize points on the track: 922 923 puckett 1.1.2.3 c write(,) 'event = ',gen_event_id_number 924 925 npoints = 0 926 927 oldbestcombo = 0 928 929 xptrack = hfpp_track_fine(ifpp,itrack,1) 930 xtrack = hfpp_track_fine(ifpp,itrack,2) 931 yptrack = hfpp_track_fine(ifpp,itrack,3) 932 ytrack = hfpp_track_fine(ifpp,itrack,4) 933 934 do ichamber=1,h_fpp_n_dcinset 935 do ilayer=1,h_fpp_n_dclayers 936 icluster = hfpp_trackcluster(ifpp,ichamber,ilayer,itrack) 937 938 if(icluster.gt.0) then 939 do ihit=1,hfpp_nhitsincluster(ifpp,ichamber,ilayer,icluster) 940 hit = hfpp_clusters(ifpp,ichamber,ilayer,icluster,ihit) 941 942 wire = hfpp_raw_wire(hit) 943 944 puckett 1.1.2.3 if(hfpp_drift_dist(ifpp,ichamber,ilayer,wire).ne.h_fpp_bad_drift) 945 $ then ! this hit is on the track: 946 npoints = npoints + 1 947 if(npoints.le.h_fpp_max_fitpoints) then 948 points(npoints,1) = hfpp_layeroffset(ifpp,ichamber,ilayer) 949 $ + hfpp_spacing(ifpp,ichamber,ilayer)float(wire) 950 points(npoints,2) = hfpp_layerz(ifpp,ichamber,ilayer) 951 sigma2s(npoints) = hfpp_resolution(ifpp,ichamber,ilayer) 952 projects(npoints,1) = hfpp_direction(ifpp,ichamber,ilayer,1) 953 projects(npoints,2) = hfpp_direction(ifpp,ichamber,ilayer,2) 954 drifts(npoints) = abs(hfpp_drift_dist(ifpp,ichamber,ilayer,wire)) 955 956 chambers(npoints) = ichamber 957 layers(npoints) = ilayer 958 wires(npoints) = wire 959 960 Px = projects(npoints,1) 961 Py = projects(npoints,2) 962 963 z = points(npoints,2) 964 965 puckett 1.1.2.3 wtrack = Px(xtrack + xptrackz) + Py(ytrack + yptrackz) 966 967 if(abs(wtrack - points(npoints,1) - drifts(npoints)) 968 $ .le.abs(wtrack - points(npoints,1) + drifts(npoints)) ) 969 $ then ! track crosses on + side of wire-->assume hit originally had a + sign! 970 oldbestcombo = oldbestcombo + 2(npoints-1) 971 endif 972 973 endif 974 endif 975 enddo 976 endif 977 enddo 978 enddo 979 980 if(npoints.ne.hfpp_track_nhits(ifpp,itrack)) return 981 982 oldchi2 = hfpp_track_chi2(ifpp,itrack) 983 oldtheta = hfpp_track_theta(ifpp,itrack) 984 oldphi = hfpp_track_phi(ifpp,itrack) 985 oldsclose = hfpp_track_sclose(ifpp,itrack) 986 puckett 1.1.2.3 oldzclose = hfpp_track_zclose(ifpp,itrack) 987 oldconetest = hfpp_track_conetest(ifpp,itrack) 988 989 reftrack(1) = hsxp_fp 990 reftrack(2) = hsx_fp 991 reftrack(3) = hsyp_fp 992 reftrack(4) = hsy_fp 993 994 if(ifpp.eq.2.and.hfpp2_best_reference(itrack).gt.0) then 995 track = hfpp2_best_reference(itrack) 996 997 oldtheta = hfpp_track_theta(ifpp+1,itrack) 998 oldphi = hfpp_track_phi(ifpp+1,itrack) 999 oldsclose = hfpp_track_sclose(ifpp+1,itrack) 1000 oldzclose = hfpp_track_zclose(ifpp+1,itrack) 1001 oldconetest = hfpp_track_conetest(ifpp+1,itrack) 1002 1003 reftrack(1) = hfpp_track_dx(1,track) 1004 reftrack(2) = hfpp_track_x(1,track) 1005 reftrack(3) = hfpp_track_dy(1,track) 1006 reftrack(4) = hfpp_track_y(1,track) 1007 puckett 1.1.2.3 1008 endif 1009 1010 c now we have baseline against which to compare other left-right combos. 1011 1012 zmin = hfpp_prune_zclose(2(ifpp-1)+1) 1013 zmax = hfpp_prune_zclose(2(ifpp-1)+2) 1014 1015 ncombos = 2npoints 1016 1017 first = .true. 1018 1019 newbestcombo = -1 1020 1021 maxthetadiff = 0.0 1022 1023 minsclose = oldsclose 1024 mintheta = oldtheta 1025 1026 nambig = 0 1027 1028 puckett 1.1.2.3 do icombo=0,ncombos-1 1029 1030 ambig(icombo+1) = .false. 1031 1032 if(icombo.ne.oldbestcombo) then 1033 do ipoint=1,npoints 1034 if(jbit(icombo,ipoint).eq.1) then 1035 plusminustest(ipoint) = 1.0 1036 else 1037 plusminustest(ipoint) = -1.0 1038 endif 1039 1040 hitpos(ipoint,1) = points(ipoint,1) + 1041 $ drifts(ipoint)plusminustest(ipoint) 1042 hitpos(ipoint,2) = points(ipoint,2) 1043 1044 enddo 1045 1046 call h_fpp_fit3d_ajp(npoints,hitpos,sigma2s,projects,testtrack) 1047 1048 chi2 = testtrack(5) 1049 puckett 1.1.2.3 1050 combochi2(icombo+1) = chi2 1051 1052 if(chi2.le.hfpp_min_chi2.and.chi2-oldchi2.le.hfpp_ambig_chi2cut(1) 1053 $ .and.chi2/oldchi2.le.hfpp_ambig_chi2cut(2)) then ! "ambiguity" 1054 1055 c first check whether any hits on this combo which disagree with the best chi2 combo have 1056 c "large" drift distances, i.e., big enough that changing the sign makes a non-trivial difference in track 1057 c reconstruction: 1058 1059 anyambig = .false. 1060 1061 maxambig = 0.0 1062 1063 do ipoint=1,npoints 1064 sign = jbit(icombo,ipoint) 1065 oldsign = jbit(oldbestcombo,ipoint) 1066 if(sign.ne.oldsign.and. 1067 $ abs(drifts(ipoint)).ge.1.25sqrt(sigma2s(ipoint))) 1068 $ then 1069 anyambig = .true. 1070 puckett 1.1.2.3 1071 if(2.0drifts(ipoint).gt.maxambig) then 1072 maxambig = 2.0drifts(ipoint) 1073 endif 1074 1075 endif 1076 enddo 1077 1078 c maximum angular displacement of new track from old track is approximately equal to atan(maxambig/dz) 1079 1080 maxdtheta = atan(maxambig/abs(points(npoints,2)-points(1,2)) ) 1081 c$$$ write(,*) 1082 c$$$ $ 'max. possible angular displacement this track=',maxdtheta 1083 1084 if(maxdtheta.gt.maxthetadiff) maxthetadiff = maxdtheta 1085 1086 if(anyambig) then 1087 ambig(icombo+1) = .true. 1088 c this combo is "ambiguous": similar chi2 to the best combo and differs in sign from the best 1089 c combo in at least one "large" drift distance hit. Calculate theta,phi,sclose,and zclose: 1090 1091 puckett 1.1.2.3 c first, transform from dc to fp coords: 1092 1093 nambig = nambig + 1 1094 ambigcombos(nambig) = icombo 1095 1096 dccoords(1) = testtrack(1) ! dx/dz 1097 dccoords(2) = testtrack(3) ! dy/dz 1098 dccoords(3) = 1.0 ! dz/dz 1099 1100 call h_fpp_dc2fp(ifpp,.true.,dccoords,fpcoords) 1101
1102 puckett 1.1.2.4 fpptrack(3) = fpcoords(1) ! dx/dz 1103 fpptrack(4) = fpcoords(2) ! dy/dz
1104 puckett 1.1.2.3 1105 dccoords(1) = testtrack(2) ! x 1106 dccoords(2) = testtrack(4) ! y 1107 dccoords(3) = 0.0 ! z 1108 1109 call h_fpp_dc2fp(ifpp,.false.,dccoords,fpcoords) 1110
1111 puckett 1.1.2.4 fpptrack(1) = fpcoords(1) - fpcoords(3)fpptrack(3) ! x - x'z 1112 fpptrack(2) = fpcoords(2) - fpcoords(3)fpptrack(4) ! y - y'z 1113 1114 call h_fpp_align(ifpp,fpptrack,newfpptrack) 1115 1116 fpptrack(1) = newfpptrack(3) 1117 fpptrack(2) = newfpptrack(1) 1118 fpptrack(3) = newfpptrack(4) 1119 fpptrack(4) = newfpptrack(2)
1120 puckett 1.1.2.3 1121 call h_fpp_relative_angles(reftrack(1),reftrack(3),fpptrack(1),fpptrack(3),theta,phi) 1122 call h_fpp_closest(reftrack,fpptrack,sclose,zclose) 1123 conetest = 1 1124 call h_fpp_conetest(reftrack,ifpp,zclose,theta,conetest) 1125 1126 if(theta.lt.mintheta) mintheta = theta 1127 if(sclose.lt.minsclose) minsclose = sclose 1128 1129 combotheta(icombo+1) = theta 1130 combophi(icombo+1) = phi 1131 combosclose(icombo+1) = sclose 1132 combozclose(icombo+1) = zclose 1133 comboconetest(icombo+1) = conetest 1134 1135 endif 1136 endif 1137 endif 1138 enddo 1139 1140 if(nambig.gt.0.and.(mintheta.lt.oldtheta.or.minsclose.lt.oldsclose)) 1141 puckett 1.1.2.3 $ then 1142 1143 c hfpp_track_nambig(ifpp,itrack) = nambig 1144 1145 c reset minsclose and oldsclose: 1146 if(mintheta.le.hfpp_ambig_smallthetatest) then ! track consistent with theta=0: 1147 c use smallest theta regardless of chi2,sclose,zclose: 1148 minthetatest = mintheta 1149 mintheta = oldtheta 1150 do icombo=1,nambig 1151 combo = ambigcombos(icombo)+1 1152 if(combotheta(combo).lt.mintheta) then 1153 newbestcombo = combo-1 1154 mintheta = combotheta(combo) 1155 endif 1156 enddo 1157 else ! track not consistent with theta=0: use sclose instead: 1158 minstest = hfpp_ambig_sclosetest 1159 do icombo=1,nambig 1160 combo = ambigcombos(icombo)+1 1161 if( (combosclose(combo)/oldsclose)2.lt.minstest .and. 1162 puckett 1.1.2.3 $ zmin.lt.combozclose(combo).and. 1163 $ zmax.gt.combozclose(combo) ) then 1164 minstest = (combosclose(combo)/oldsclose)2 1165 newbestcombo = combo-1 1166 endif 1167 enddo 1168 endif 1169 endif 1170 1171 if(newbestcombo.ge.0) then 1172 1173 hfpp_track_nambig(ifpp,itrack) = nambig 1174 c refit the track using the "new best combo" 1175 1176 c$$$ write(,) 'found new best combo FPP,itrack=',ifpp,itrack 1177 c$$$ write(,) 'nhits=',npoints 1178 c$$$ write(,) 'old combo, nambig = ',oldbestcombo,nambig 1179 c$$$ write(,) 'old chi2,sclose,zclose,theta,phi,conetest=',oldchi2,oldsclose, 1180 c$$$ $ oldzclose,oldtheta,oldphi,oldconetest 1181 1182 do ipoint=1,npoints 1183 puckett 1.1.2.3 if(jbit(newbestcombo,ipoint).eq.1) then 1184 plusminusbest(ipoint) = 1.0 1185 else 1186 plusminusbest(ipoint) = -1.0 1187 endif 1188 1189 ichamber = chambers(ipoint) 1190 ilayer = layers(ipoint) 1191 iwire = wires(ipoint) 1192 1193 hfpp_drift_dist(ifpp,ichamber,ilayer,iwire) = 1194 $ plusminusbest(ipoint)drifts(ipoint) 1195 1196 hitpos(ipoint,1) = points(ipoint,1) + 1197 $ plusminusbest(ipoint)drifts(ipoint) 1198 hitpos(ipoint,2) = points(ipoint,2) 1199 enddo 1200 1201 call h_fpp_fit3d_ajp(npoints,hitpos,sigma2s,projects,testtrack) 1202 c copy the chi2 of the new track to common block variables: 1203 hfpp_track_chi2(ifpp,itrack) = testtrack(5) 1204 puckett 1.1.2.3 1205 do j=1,4 1206 hfpp_track_fine(ifpp,itrack,j) = testtrack(j) 1207 enddo 1208 1209 c re-calculate everything that was modified by this routine: 1210 1211 dccoords(1) = testtrack(1) ! dx/dz 1212 dccoords(2) = testtrack(3) ! dy/dz 1213 dccoords(3) = 1.0 ! dz/dz 1214 1215 call h_fpp_dc2fp(ifpp,.true.,dccoords,fpcoords) 1216
1217 puckett 1.1.2.4 fpptrack(3) = fpcoords(1) 1218 fpptrack(4) = fpcoords(2)
1219 puckett 1.1.2.3 1220 dccoords(1) = testtrack(2) ! x 1221 dccoords(2) = testtrack(4) ! y 1222 dccoords(3) = 0.0 ! z 1223 1224 call h_fpp_dc2fp(ifpp,.false.,dccoords,fpcoords) 1225
1226 puckett 1.1.2.4 fpptrack(1) = fpcoords(1) - fpcoords(3)fpptrack(3) ! x - x'z 1227 fpptrack(2) = fpcoords(2) - fpcoords(3)fpptrack(4) ! y - y'z 1228 1229 call h_fpp_align(ifpp,fpptrack,newfpptrack) 1230 1231 fpptrack(1) = newfpptrack(3) 1232 fpptrack(2) = newfpptrack(1) 1233 fpptrack(3) = newfpptrack(4) 1234 fpptrack(4) = newfpptrack(2)
1235 puckett 1.1.2.3 1236 hfpp_track_dx(ifpp,itrack) = fpptrack(1) 1237 hfpp_track_x(ifpp,itrack) = fpptrack(2) 1238 hfpp_track_dy(ifpp,itrack) = fpptrack(3) 1239 hfpp_track_y(ifpp,itrack) = fpptrack(4) 1240 1241 hmstrack(1) = hsxp_fp 1242 hmstrack(2) = hsx_fp 1243 hmstrack(3) = hsyp_fp 1244 hmstrack(4) = hsy_fp 1245 1246 call h_fpp_relative_angles(hmstrack(1),hmstrack(3),fpptrack(1), 1247 $ fpptrack(3),theta,phi) 1248 1249 hfpp_track_theta(ifpp,itrack) = theta 1250 hfpp_track_phi(ifpp,itrack) = phi 1251 1252 call h_fpp_closest(hmstrack,fpptrack,sclose,zclose) 1253 1254 hfpp_track_sclose(ifpp,itrack) = sclose 1255 hfpp_track_zclose(ifpp,itrack) = zclose 1256 puckett 1.1.2.3 1257 conetest = 1 1258 call h_fpp_conetest(hmstrack,ifpp,zclose,theta,conetest) 1259 1260 hfpp_track_conetest(ifpp,itrack) = conetest 1261 1262 if(ifpp.eq.2.and.hfpp2_best_reference(itrack).gt.0) then 1263 c reftrack has already been initialized to the best track chosen from 1264 c FPP1:
1265 puckett 1.1.2.4 reftrack(1) = hfpp_track_dx(1,hfpp2_best_reference(itrack)) 1266 reftrack(2) = hfpp_track_x(1,hfpp2_best_reference(itrack)) 1267 reftrack(3) = hfpp_track_dy(1,hfpp2_best_reference(itrack)) 1268 reftrack(4) = hfpp_track_y(1,hfpp2_best_reference(itrack)) 1269
1270 puckett 1.1.2.3 call h_fpp_relative_angles(reftrack(1),reftrack(3),fpptrack(1), 1271 $ fpptrack(3),theta,phi) 1272 hfpp_track_theta(ifpp+1,itrack) = theta 1273 hfpp_track_phi(ifpp+1,itrack) = phi 1274 1275 call h_fpp_closest(reftrack,fpptrack,sclose,zclose)
1276 puckett 1.1.2.1
1277 puckett 1.1.2.3 hfpp_track_sclose(ifpp+1,itrack) = sclose 1278 hfpp_track_zclose(ifpp+1,itrack) = zclose 1279 1280 conetest = 1 1281 call h_fpp_conetest(reftrack,ifpp,zclose,theta,conetest) 1282 1283 hfpp_track_conetest(ifpp+1,itrack) = conetest 1284 1285 endif 1286 c$$$ write(,) 'new combo=',newbestcombo 1287 c$$$ write(,) 'new chi2,sclose,zclose,theta,phi,conetest=',testtrack(5),sclose,zclose,theta,phi,conetest 1288 1289 endif 1290 1291 return 1292 end

Analyzer/Replay: Mark Jones, Documents: Stephen Wood