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.14 subroutine h_fpp_tracking_ajp_alt(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 logical abort 15 character() err 16 integer4 dcset 17 18 character14 here 19 parameter (here= 'h_fpp_tracking') 20 c this routine is simply going to find the best combination of 21 c one hit per plane, with "best" defined by chi2: 22 puckett 1.1.2.14 23 logical sufficient_hits 24 25 integer i,j,k,l,m,n 26 integer itrack 27 28 integer4 bestcombo(h_fpp_n_dcinset,h_fpp_n_dclayers) 29 real4 HMStrack(4) 30 real4 FPPtrack(4),newFPPtrack(4) 31 real4 DCslopes(3),DCcoords(3), FPcoords(3), FPslopes(3) 32 33 real4 PI 34 parameter(PI=3.141592653) 35 36 real4 theta,mintheta,phi,sclose,minsclose,zclose 37 real4 criterion,mincriterion 38 real4 chi2,minchi2 39 integer4 conetest,bestref,jtrack,iterations,nhitsrequired 40 integer4 nhitsintrack,npoints 41 integer4 ichamber,ilayer,icluster 42 43 puckett 1.1.2.14 integer4 ntracksnew 44 c integer4 nhitsrequired_firsttrack 45 46 real4 SimpleTrack(6), FullTrack(7), temptrack(5) 47 48 logical ontrack(h_fpp_n_dcinset,h_fpp_n_dclayers) ! keep track of whether a cluster ends up on a track 49 integer4 BestClusters(H_FPP_N_DCINSET,H_FPP_N_DCLAYERS) 50 integer4 hitcluster(h_fpp_n_dcinset,h_fpp_n_dclayers) 51 52 logical firsttry,track_good 53 54 ABORT= .FALSE. 55 err= ' ' 56 57 c write(,) 'new event' 58 59 c initialize ntracks to zero: 60 hfpp_n_tracks(dcset) = 0 61 c determine whether enough hits for tracking: 62 call h_fpp_tracking_freehitcount(dcset,sufficient_hits) 63 c find the individual combination of one hit per plane which gives the smallest chi2 of drift tracking: 64 puckett 1.1.2.14 do while (sufficient_hits .and. (HFPP_N_tracks(DCset).lt.H_FPP_MAX_TRACKS)) 65 66 ntracksnew = 0 67 68 nhitsrequired = h_fpp_n_planes + 1 69 c iterations = 0 70 do while(nhitsrequired.ge.hfpp_minsethits) 71 iterations = 0 72 nhitsintrack = 0 73 call h_fpp_tracking_nexthitcombo(dcset,nhitsrequired,nhitsintrack,hitcluster) 74 75 firsttry = .true. 76 77 if(hfpp_n_tracks(dcset).gt.0.and. 78 $ nhitsrequired.lt. 79 $ hfpp_track_nlayers(dcset,hfpp_n_tracks(dcset))) exit 80 81 do while(nhitsintrack.gt.0) 82 83 c if(iterations.gt.hfpp_maxcombos+1) then 84 c write(,) 'WARNING: max FPP hit combos reached.' 85 puckett 1.1.2.14 c write(,) 'consider raising the limit or using ' 86 c write(,) 'tighter raw timing cuts' 87 c exit 88 c endif 89 90 iterations = iterations + 1 91 92 call h_fpp_fit_simple(dcset,hitcluster,npoints,temptrack,abort,err) 93 if (ABORT) then 94 call g_add_path(here,err) 95 return 96 endif 97 98 if(temptrack(5).ge.0.0.and.temptrack(5).le.hfpp_aok_chi2) 99 $ then ! passed simple chi2 cut 100 c call drift tracking: 101 do j=1,5 102 simpletrack(j) = temptrack(j) 103 enddo 104 simpletrack(6) = npoints 105 106 puckett 1.1.2.14 call h_fpp_tracking_drifttrack_ajp(dcset,simpletrack,hitcluster, 107 $ ontrack,track_good,fulltrack,nhitsrequired,1,abort,err) 108 109 if (ABORT) then 110 call g_add_path(here,err) 111 return 112 endif 113 114 if(track_good) then ! track passes chi2 cut: 115 if(firsttry.or.fulltrack(5).lt.minchi2) then 116 firsttry = .false. 117 minchi2 = fulltrack(5) 118 do ichamber=1,h_fpp_n_dcinset 119 do ilayer=1,h_fpp_n_dclayers 120 bestclusters(ichamber,ilayer) = hitcluster(ichamber,ilayer) 121 enddo 122 enddo 123 endif 124 endif 125 endif 126 c get the next free hit combo: 127 puckett 1.1.2.14 call h_fpp_tracking_nexthitcombo(dcset,nhitsrequired,nhitsintrack,hitcluster) 128 enddo ! end loop while hitsintrack>0 129 c write(,) 'finished loop over hit combos' 130 c at this point bestclusters contains the combination of one hit per plane with the best chi2: 131 if(.not. firsttry) then ! at least one track was found: add the best track to the track array: 132 c fit the simple and full tracks using the best clusters: 133 call h_fpp_fit_simple(dcset,bestclusters,npoints,temptrack,abort,err) 134 do j=1,5 135 simpletrack(j) = temptrack(j) 136 enddo 137 simpletrack(6) = npoints 138 call h_fpp_tracking_drifttrack_ajp(dcset,simpletrack,bestclusters, 139 $ ontrack,track_good,fulltrack,nhitsrequired,1,abort,err)
140 puckett 1.1.2.15
141 puckett 1.1.2.14 itrack = hfpp_n_tracks(dcset)+1 142 if(itrack.le.h_fpp_max_tracks) then 143 144 c if(itrack.eq.1) nhitsrequired_firsttrack = int(fulltrack(7)) 145 146 ntracksnew = ntracksnew + 1 147 148 hfpp_n_tracks(dcset) = itrack 149 hfpp_track_nlayers(dcset,itrack) = int(fulltrack(7)) 150 151 do j=1,4 152 hfpp_track_fine(dcset,itrack,j) = fulltrack(j) 153 enddo 154 155 c Mark the hits in the best candidate track as used! 156 do ichamber = 1,h_fpp_n_dcinset 157 do ilayer = 1,h_fpp_n_dclayers 158 icluster = bestclusters(ichamber,ilayer) 159 if(icluster.gt.0.and. 160 $ ontrack(ichamber,ilayer)) then 161 hfpp_clusterintrack(dcset,ichamber,ilayer,icluster) = itrack 162 puckett 1.1.2.14 hfpp_trackcluster(dcset,ichamber,ilayer,itrack) = icluster 163 endif 164 enddo
165 puckett 1.1.2.15 enddo
166 puckett 1.1.2.14 c 167 hfpp_track_chi2(dcset,itrack) = fulltrack(5) 168 hfpp_track_nhits(dcset,itrack) = int(fulltrack(6)) 169 170 do j=1,6 171 hfpp_track_rough(dcset,itrack,j) = simpletrack(j) 172 enddo 173 174 hfpp_track_uniq(dcset,itrack) = .true. 175 c convert from local to fp coords: 176 dccoords(1) = fulltrack(1) 177 dccoords(2) = fulltrack(3) 178 dccoords(3) = 1. 179 call h_fpp_dc2fp(dcset,.true.,dccoords,fpcoords) 180 181 fpptrack(3) = fpcoords(1) ! x' in fp coords 182 fpptrack(4) = fpcoords(2) ! y' in fp coords 183 184 * now transform coordinates from local FPP system to HMS fp coords: 185 186 dccoords(1) = fulltrack(2) ! x 187 puckett 1.1.2.14 dccoords(2) = fulltrack(4) ! y 188 dccoords(3) = 0.0 189 call h_fpp_dc2fp(dcset,.false.,dccoords,fpcoords) 190 191 fpptrack(1) = fpcoords(1) - fpcoords(3) * fpptrack(3) 192 fpptrack(2) = fpcoords(2) - fpcoords(3) * fpptrack(4) 193 194 call h_fpp_align(dcset,fpptrack,newfpptrack) 195 196 hfpp_track_x(dcset,itrack) = newfpptrack(1) 197 hfpp_track_y(dcset,itrack) = newfpptrack(2) 198 hfpp_track_dx(dcset,itrack) = newfpptrack(3) 199 hfpp_track_dy(dcset,itrack) = newfpptrack(4) 200 201 c calculate relative angles: 202 call h_fpp_relative_angles(hsxp_fp,hsyp_fp,hfpp_track_dx(dcset,itrack), 203 $ hfpp_track_dy(dcset,itrack),theta,phi) 204 hfpp_track_theta(dcset,itrack) = theta 205 hfpp_track_phi(dcset,itrack) = phi 206 207 * calculate closest approach parameters with respect to HMS, REGARDLESS of which FPP! 208 puckett 1.1.2.14 209 hmstrack(1) = hsxp_fp 210 hmstrack(2) = hsx_fp 211 hmstrack(3) = hsyp_fp 212 hmstrack(4) = hsy_fp 213 214 fpptrack(1) = hfpp_track_dx(dcset,itrack) 215 fpptrack(2) = hfpp_track_x(dcset,itrack) 216 fpptrack(3) = hfpp_track_dy(dcset,itrack) 217 fpptrack(4) = hfpp_track_y(dcset,itrack) 218 call h_fpp_closest(hmstrack,fpptrack,sclose,zclose) 219 hfpp_track_sclose(dcset,itrack) = sclose 220 hfpp_track_zclose(dcset,itrack) = zclose 221 222 conetest = 1 223 call h_fpp_conetest(hmstrack,dcset,zclose,theta,conetest) 224 hfpp_track_conetest(dcset,itrack) = conetest 225 226 if(dcset.eq.2.and.hfpp_n_tracks(1).gt.0) then 227 jtrack = hfpp_best_track(1) 228 229 puckett 1.1.2.14 call h_fpp_relative_angles(hfpp_track_dx(1,jtrack), 230 $ hfpp_track_dy(1,jtrack), 231 $ hfpp_track_dx(dcset,itrack), 232 $ hfpp_track_dy(dcset,itrack), 233 $ theta,phi) 234 hfpp_track_theta(dcset+1,itrack) = theta 235 hfpp_track_phi(dcset+1,itrack) = phi 236 237 hmstrack(1) = hfpp_track_dx(1,jtrack) 238 hmstrack(2) = hfpp_track_x(1,jtrack) 239 hmstrack(3) = hfpp_track_dy(1,jtrack) 240 hmstrack(4) = hfpp_track_y(1,jtrack) 241 242 fpptrack(1) = hfpp_track_dx(dcset,itrack) 243 fpptrack(2) = hfpp_track_x(dcset,itrack) 244 fpptrack(3) = hfpp_track_dy(dcset,itrack) 245 fpptrack(4) = hfpp_track_y(dcset,itrack) 246 247 call h_fpp_closest(hmstrack,fpptrack,sclose,zclose) 248 249 hfpp_track_sclose(dcset+1,itrack) = sclose 250 puckett 1.1.2.14 hfpp_track_zclose(dcset+1,itrack) = zclose 251 252 conetest = 1 253 254 call h_fpp_conetest(hmstrack,dcset,zclose,theta,conetest) 255 256 hfpp_track_conetest(dcset+1,itrack) = conetest 257 258 hfpp2_best_reference(itrack) = jtrack 259 260 if(hfpp_track_theta(dcset,itrack).lt. 261 $ hfpp_track_theta(dcset+1,itrack)) hfpp2_best_reference(itrack) = 0 262 endif 263 endif 264 265 exit ! break out of the loop and update the "free hit" count: 266 endif 267 268 c if(ntracksnew.eq.0 .and. hfpp_n_tracks(dcset).le.0 ) then 269 c nhitsrequired = nhitsrequired - 1 270 c endif 271 puckett 1.1.2.14 nhitsrequired = nhitsrequired - 1 272 273 enddo ! end loop while enough hits to do tracking 274 c write(,) 'finished loop while nhitsrequired' 275 c once the "nhits required" 276 c if enough free hits remain for a track, repeat: 277 c if no new tracks were found on this iteration, we are done: 278 if(ntracksnew.eq.0) exit 279 280 c write(,) 'ntracks=',hfpp_n_tracks(dcset) 281 282 call h_fpp_tracking_freehitcount(dcset,sufficient_hits) 283 enddo ! end loop while sufficient hits and tracks < max 284 285 c write(,) 'finished loop while suff. hits and tracks < max' 286 287 return 288 end 289
290 puckett 1.1.2.1 subroutine h_fpp_tracking_ajp(dcset,abort,err) 291 292 implicit none 293 save 294 295 include 'gen_detectorids.par' 296 include 'gen_decode_common.cmn' 297 INCLUDE 'hms_data_structures.cmn' 298 INCLUDE 'hms_geometry.cmn' 299 INCLUDE 'hms_fpp_event.cmn' 300 INCLUDE 'hms_fpp_params.cmn' 301 INCLUDE 'hms_id_histid.cmn' 302 303 character14 here 304 parameter (here= 'h_fpp_tracking') 305 306 integer4 DCset ! set of FPP DCs we are working on 307 308 logical ABORT 309 character() err 310 311 puckett 1.1.2.1 logical*4 sufficient_hits, track_good, any_track, any_good, any_great 312 313 c define arrays for candidate tracks: 314
315 puckett 1.1.2.4 integer*4 nhitsrequired 316
317 puckett 1.1.2.1 integer4 ngoodtracks 318 integer4 ncandidates,i,j,k,ichamber,ilayer 319 integer*4 nplanestemp,icluster,bestcandidate,best6,best5
320 puckett 1.1.2.9 integer*4 itrack,track,hit,ihit,wire
321 puckett 1.1.2.1 322 integer4 goodtracks(h_fpp_max_candidates) 323 real4 simpletracks(h_fpp_max_candidates,6) 324 real4 fulltracks(h_fpp_max_candidates,6) 325 real4 thetatracks(h_fpp_max_candidates) ! store track polar angle relative to HMS
326 puckett 1.1.2.3 real4 sclosetracks(h_fpp_max_candidates) ! store track distance of closest approach relative to HMS or FPP1 327 real4 zclosetracks(h_fpp_max_candidates) ! store track z of closest approach relative to HMS or FPP1 328 real4 phitracks(h_fpp_max_candidates) 329 integer4 conetesttracks(h_fpp_max_candidates) ! store cone test of tracks relative to HMS or FPP1 330 integer*4 bestreftracks(h_fpp_max_candidates)
331 puckett 1.1.2.6 logical greattrack(h_fpp_max_candidates)
332 puckett 1.1.2.1 real*4 chi2, minchi2 333
334 puckett 1.1.2.3 c logical ambiguity(h_fpp_max_candidates) 335 logical firsttheta,firstsclose
336 puckett 1.1.2.1 logical firsttry,any6,anytheta ! at this stage of the tracking, apply a maximum theta cut on all candidate tracks 337 logical any6theta ! any tracks with six planes AND passing theta cut?
338 puckett 1.1.2.3 logical anyzclose,any6zclose ! any tracks with zclose passing prune tests?
339 puckett 1.1.2.6 logical anyconetest,any6conetest
340 puckett 1.1.2.1 341 integer4 hitcombos(h_fpp_max_candidates,h_fpp_n_dcinset,h_fpp_n_dclayers) 342 343 logical hitsintrack(h_fpp_max_candidates,h_fpp_n_dcinset,h_fpp_n_dclayers) 344 logical ontrack(h_fpp_n_dcinset,h_fpp_n_dclayers) ! keep track of whether a cluster ends up on a track 345 346 integer4 candidate_nplanes(h_fpp_max_candidates) 347 integer*4 candidate_nhits(h_fpp_max_candidates) 348
349 puckett 1.1.2.8 c integer4 candidate_wires(h_fpp_max_candidates,h_fpp_max_fitpoints) 350 c real4 candidate_drifts(h_fpp_max_candidates,h_fpp_max_fitpoints) 351
352 puckett 1.1.2.1 real4 SimpleTrack(6), FullTrack(7) 353 354 integer4 BestClusters(H_FPP_N_DCINSET,H_FPP_N_DCLAYERS) 355 356 real*4 HMStrack(4)
357 puckett 1.1.2.4 real*4 FPPtrack(4),newFPPtrack(4)
358 puckett 1.1.2.3 real*4 DCslopes(3),DCcoords(3), FPcoords(3), FPslopes(3)
359 puckett 1.1.2.1 360 real*4 PI 361 parameter(PI=3.141592653) 362
363 puckett 1.1.2.3 real4 theta,mintheta,phi,sclose,minsclose,zclose 364 real4 criterion,mincriterion
365 puckett 1.1.2.1 integer*4 conetest,bestref,jtrack
366 puckett 1.1.2.3 integer*4 iteration 367
368 puckett 1.1.2.5 logical zgood
369 puckett 1.1.2.3 logical goodhms 370 371 real*4 prob
372 puckett 1.1.2.1 373 ABORT= .FALSE. 374 err= ' ' 375 376 any_track = .false. 377 any_good = .false. 378 any_great = .false. 379 380 hfpp_n_tracks(dcset) = 0 381 382 call h_fpp_tracking_freehitcount(dcset,sufficient_hits) 383
384 puckett 1.1.2.4 iteration = 0
385 puckett 1.1.2.3 386 goodhms = abs(hsxp_tar).le..08.and.abs(hsyp_tar).le..04.and.abs(hsy_tar)<5..and. 387 $ abs(hsdelta)<10. 388
389 puckett 1.1.2.1 do while (sufficient_hits .and. (HFPP_N_tracks(DCset).lt.H_FPP_MAX_TRACKS))
390 puckett 1.1.2.3 391 c$$$ if(goodhms) then 392 c$$$ write(,) 'FPP tracking: chamber = ',dcset,' iteration = ',iteration 393 c$$$ endif 394 iteration = iteration + 1
395 puckett 1.1.2.1 396 ncandidates = 0 397 any6 = .false. 398 anytheta = .false. 399 any6theta = .false.
400 puckett 1.1.2.3 anyzclose = .false. 401 any6zclose = .false.
402 puckett 1.1.2.6 anyconetest = .false. 403 any6conetest = .false.
404 puckett 1.1.2.3
405 puckett 1.1.2.6 any_great = .false. 406
407 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 408 c until we exceed the maximum number of candidate tracks: 409 410 c initialize local tracking arrays: 411 412 do i=1,h_fpp_max_candidates 413 candidate_nplanes(i) = 0 414 candidate_nhits(i) = 0 415 goodtracks(i) = 0 416 do ichamber=1,h_fpp_n_dcinset 417 do ilayer=1,h_fpp_n_dclayers 418 hitcombos(i,ichamber,ilayer) = 0 419 hitsintrack(i,ichamber,ilayer) = .false. 420 enddo 421 enddo 422 enddo 423 424 c write(,) 'Start of tracking, FPP,ntracks=',dcset,hfpp_n_tracks(dcset) 425
426 puckett 1.1.2.4 nhitsrequired = 0 427 428 c write(,) 'before simple tracking, nhits required=',nhitsrequired 429 430 call h_fpp_tracking_simple_ajp(dcset,hitcombos,simpletracks,ncandidates,nhitsrequired,abort,err)
431 puckett 1.1.2.1 if (ABORT) then 432 call g_add_path(here,err) 433 return 434 endif
435 puckett 1.1.2.3
436 puckett 1.1.2.4 c$$$ write(,) 'after simple tracking, nhits required,ncandidate=', 437 c$$$ $ nhitsrequired,ncandidates 438 439 if(iteration.eq.1) 440 $ hfpp_n_simple(dcset,1) = ncandidates
441 puckett 1.1.2.1
442 puckett 1.1.2.3 c$$$ if(goodhms) then 443 c$$$ write(,) 'number of candidates after simple tracking=',ncandidates 444 c$$$ endif
445 puckett 1.1.2.1 c$$$ if (.true.) then 446 c$$$ if (int(SimpleTrack(6)).le.0) then 447 c$$$ if(hidFPP_trkrough(DCset,6).gt.0) call hf1(hidFPP_trkrough(DCset,6),0.,1.) !Nraw 448 c$$$ else 449 c$$$ if(hidFPP_trkrough(DCset,1).gt.0) call hf1(hidFPP_trkrough(DCset,1),SimpleTrack(1),1.) !mx 450 c$$$ if(hidFPP_trkrough(DCset,2).gt.0) call hf1(hidFPP_trkrough(DCset,2),SimpleTrack(2),1.) !bx 451 c$$$ if(hidFPP_trkrough(DCset,3).gt.0) call hf1(hidFPP_trkrough(DCset,3),SimpleTrack(3),1.) !my 452 c$$$ if(hidFPP_trkrough(DCset,4).gt.0) call hf1(hidFPP_trkrough(DCset,4),SimpleTrack(4),1.) !by 453 c$$$ if(hidFPP_trkrough(DCset,5).gt.0) call hf1(hidFPP_trkrough(DCset,5),SimpleTrack(5),1.) !chi2 454 c$$$ if(hidFPP_trkrough(DCset,6).gt.0) call hf1(hidFPP_trkrough(DCset,6),SimpleTrack(6),1.) !Nraw 455 c$$$ endif 456 c$$$ endif 457 458 * * quit trying to make more tracks if we are out of hits 459 * * or if we couldnt make a good one now 460 c if (int(SimpleTrack(6)).le.0) exit 461 c if (int(SimpleTrack(5)).eq.H_FPP_BAD_CHI2) exit 462 463 c any_track = .true. 464 465 if(ncandidates.eq.0) exit 466 puckett 1.1.2.1 467 ngoodtracks = 0 468
469 puckett 1.1.2.14 c write(,) 'NEW EVENT' 470
471 puckett 1.1.2.4 134 do i=1,ncandidates
472 puckett 1.1.2.14 c write(,) 'icandidate,simple chi2=',i,simpletracks(i,5) 473 c write(,) 'simple track (dx,x,dy,y)=(',(simpletracks(i,j),j=1,4),')'
474 puckett 1.1.2.1 FullTrack(1) = H_FPP_BAD_COORD ! mx 475 FullTrack(2) = H_FPP_BAD_COORD ! bx 476 FullTrack(3) = H_FPP_BAD_COORD ! my 477 FullTrack(4) = H_FPP_BAD_COORD ! by 478 FullTrack(5) = H_FPP_BAD_CHI2 479 FullTrack(6) = 0. ! number of hits 480 FullTrack(7) = 0. ! number of planes 481 c write(,) 'candidate ',i,' simple chi2=',simpletracks(i,5) 482 483 do j=1,6 484 simpletrack(j) = simpletracks(i,j) 485 c write(,) 'i,simpletrack(i)=',j,simpletrack(j) 486 enddo 487 488 do ichamber=1,h_fpp_n_dcinset 489 do ilayer=1,h_fpp_n_dclayers 490 bestclusters(ichamber,ilayer) = hitcombos(i,ichamber,ilayer) 491 c write(,) 'chbr,pln,clstr=',ichamber,ilayer,bestclusters(ichamber,ilayer) 492 enddo 493 enddo 494 495 puckett 1.1.2.1 call h_fpp_tracking_drifttrack_ajp(dcset,simpletrack,bestclusters,
496 puckett 1.1.2.4 $ ontrack,track_good,fulltrack,nhitsrequired,1,abort,err)
497 puckett 1.1.2.1 498 499 c write(,) 'track_good = ',track_good 500 if (ABORT) then 501 call g_add_path(here,err) 502 return 503 endif 504
505 puckett 1.1.2.3 if(track_good) then ! tracking with drift resulted in a reasonable track: 506 507 * calculate everything so we can check its quality/value against quantities other than just chi2: 508 * transform from FPP to focal plane coordinate system: 509 * slopes first: 510 dcslopes(1) = fulltrack(1) 511 dcslopes(2) = fulltrack(3) 512 dcslopes(3) = 1.0 513 514 call h_fpp_dc2fp(dcset,.true.,dcslopes,fpslopes)
515 puckett 1.1.2.1
516 puckett 1.1.2.3 * Next, transform coordinates from FPP to focal plane system: 517 dccoords(1) = fulltrack(2) 518 dccoords(2) = fulltrack(4) 519 dccoords(3) = 0.0 520 521 call h_fpp_dc2fp(dcset,.false.,dccoords,fpcoords) 522 523 fpcoords(1) = fpcoords(1) - fpslopes(1)fpcoords(3) 524 fpcoords(2) = fpcoords(2) - fpslopes(2)fpcoords(3)
525 puckett 1.1.2.4 526 fpptrack(1) = fpcoords(1) 527 fpptrack(2) = fpcoords(2) 528 fpptrack(3) = fpslopes(1) 529 fpptrack(4) = fpslopes(2) 530 531 call h_fpp_align(dcset,fpptrack,newfpptrack)
532 puckett 1.1.2.3 * initialize bestref to zero: 533 bestref = 0 534 * initialize "hmstrack" to the actual hms track: 535 hmstrack(1) = hsxp_fp 536 hmstrack(2) = hsx_fp 537 hmstrack(3) = hsyp_fp 538 hmstrack(4) = hsy_fp 539 * initialize "fpptrack" to the current track:
540 puckett 1.1.2.4 c$$$ fpptrack(1) = fpslopes(1) 541 c$$$ fpptrack(2) = fpcoords(1) 542 c$$$ fpptrack(3) = fpslopes(2) 543 c$$$ fpptrack(4) = fpcoords(2) 544 fpptrack(1) = newfpptrack(3) ! dx/dz 545 fpptrack(2) = newfpptrack(1) ! x 546 fpptrack(3) = newfpptrack(4) ! dy/dz 547 fpptrack(4) = newfpptrack(2) ! y 548
549 puckett 1.1.2.3 * calculate closest approach parameters relative to hms track: 550 call h_fpp_closest(hmstrack,fpptrack,sclose,zclose) 551 call h_fpp_relative_angles(hmstrack(1),hmstrack(3),fpptrack(1),fpptrack(3),theta,phi)
552 puckett 1.1.2.4 * in FPP2, we use the "best" track from FPP1 as a reference track:
553 puckett 1.1.2.6 mintheta = theta 554 bestref = 0
555 puckett 1.1.2.4 if(dcset.eq.2.and.hfpp_n_tracks(1).gt.0) then
556 puckett 1.1.2.13 c do jtrack = 1, hfpp_n_tracks(1) 557 jtrack = hfpp_best_track(1) 558 hmstrack(1) = hfpp_track_dx(1,jtrack) 559 hmstrack(2) = hfpp_track_x(1,jtrack) 560 hmstrack(3) = hfpp_track_dy(1,jtrack) 561 hmstrack(4) = hfpp_track_y(1,jtrack) 562 call h_fpp_closest(hmstrack,fpptrack,sclose,zclose) 563 call h_fpp_relative_angles(hmstrack(1),hmstrack(3),fpptrack(1),fpptrack(3),theta,phi) 564 if(theta.lt.mintheta) then 565 bestref = jtrack 566 bestreftracks(i) = bestref 567 mintheta = theta 568 endif 569 c enddo
570 puckett 1.1.2.6 c now, calculate everything one more time using bestref: 571 if(bestref.gt.0) then 572 hmstrack(1) = hfpp_track_dx(1,bestref) 573 hmstrack(2) = hfpp_track_x(1,bestref) 574 hmstrack(3) = hfpp_track_dy(1,bestref) 575 hmstrack(4) = hfpp_track_y(1,bestref)
576 puckett 1.1.2.13 else ! revert to HMS track:
577 puckett 1.1.2.6 hmstrack(1) = hsxp_fp 578 hmstrack(2) = hsx_fp 579 hmstrack(3) = hsyp_fp 580 hmstrack(4) = hsy_fp 581 endif
582 puckett 1.1.2.13
583 puckett 1.1.2.3 call h_fpp_closest(hmstrack,fpptrack,sclose,zclose)
584 puckett 1.1.2.4 call h_fpp_relative_angles(hmstrack(1),hmstrack(3),fpptrack(1),fpptrack(3),theta,phi)
585 puckett 1.1.2.13
586 puckett 1.1.2.3 endif 587
588 puckett 1.1.2.6 bestreftracks(i) = bestref
589 puckett 1.1.2.13
590 puckett 1.1.2.3 conetest = 1
591 puckett 1.1.2.1
592 puckett 1.1.2.3 call h_fpp_conetest(hmstrack,dcset,zclose,theta,conetest) 593
594 puckett 1.1.2.1 thetatracks(i) = theta
595 puckett 1.1.2.3 phitracks(i) = phi
596 puckett 1.1.2.1
597 puckett 1.1.2.3 sclosetracks(i) = sclose 598 zclosetracks(i) = zclose 599 600 conetesttracks(i) = conetest 601 602 if(theta.le.hfpp_prune_thetamax(dcset)*PI/180.0) 603 $ anytheta = .true. 604
605 puckett 1.1.2.1 do j=1,6 606 fulltracks(i,j) = fulltrack(j) 607 enddo 608 609 nplanestemp = int(fulltrack(7)) 610 611 candidate_nplanes(i) = nplanestemp 612 candidate_nhits(i) = int(fulltrack(6)) 613 614 ngoodtracks = ngoodtracks + 1 615 goodtracks(ngoodtracks) = i 616
617 puckett 1.1.2.6 greattrack(i) = ( conetest.eq.1.and. 618 $ sclose.le.hfpp_prune_sclose(dcset) 619 $ .and.zclose.ge.hfpp_prune_zclose(2(dcset-1)+1)- 620 $ hfpp_prune_zslop(dcset)/tan(theta).and.zclose.le. 621 $ hfpp_prune_zclose(2(dcset-1)+2) + 622 $ hfpp_prune_zslop(dcset)/tan(theta) .and. 623 $ nplanestemp.eq.6 )
624 puckett 1.1.2.10 if(dcset.eq.2.and.bestref.gt.0) then 625 greattrack(i) = ( conetest.eq.1.and. 626 $ sclose.le.hfpp_prune_sclose(1) 627 $ .and.zclose.ge.hfpp_prune_zclose(2(dcset-1)+1)- 628 $ hfpp_prune_zslop(dcset)/tan(theta).and.zclose.le. 629 $ hfpp_prune_zclose(2(dcset-1)+2) + 630 $ hfpp_prune_zslop(dcset)/tan(theta) .and. 631 $ nplanestemp.eq.6 ) 632 endif 633
634 puckett 1.1.2.6 635 if(greattrack(i)) any_great = .true. 636
637 puckett 1.1.2.1 if(nplanestemp.eq.h_fpp_n_dcinset*h_fpp_n_dclayers) 638 $ any6 = .true. 639
640 puckett 1.1.2.6 if(conetest.eq.1) then 641 anyconetest = .true. 642 if(nplanestemp.eq.h_fpp_n_dcinset*h_fpp_n_dclayers) 643 $ then 644 any6conetest = .true. 645 endif 646 endif 647
648 puckett 1.1.2.5 if(zclose.ge.hfpp_prune_zclose(2(dcset-1)+1) - 649 $ hfpp_prune_zslop(dcset)/tan(theta).and. 650 $ zclose.le.hfpp_prune_zclose(2(dcset-1)+2) + 651 $ hfpp_prune_zslop(dcset)/tan(theta)) then
652 puckett 1.1.2.3 anyzclose = .true. 653 if(nplanestemp.eq.h_fpp_n_dcinset*h_fpp_n_dclayers) 654 $ then 655 any6zclose = .true. 656 endif 657 endif 658
659 puckett 1.1.2.1 if(theta.le.hfpp_prune_thetamax(dcset)PI/180.0 660 $ .and.nplanestemp.eq.h_fpp_n_dcinseth_fpp_n_dclayers) 661 $ any6theta = .true. 662 663 do ichamber=1,h_fpp_n_dcinset 664 do ilayer = 1,h_fpp_n_dclayers 665 if(ontrack(ichamber,ilayer)) 666 $ hitsintrack(i,ichamber,ilayer) = .true. 667 enddo 668 enddo
669 puckett 1.1.2.4
670 puckett 1.1.2.1 endif 671 enddo
672 puckett 1.1.2.6 673 c$$$ if(ngoodtracks.gt.hfpp_n_simple(dcset,2)) then
674 puckett 1.1.2.4 c$$$ hfpp_n_simple(dcset,2) = ngoodtracks 675 c$$$ endif 676 677 if(iteration.eq.1) then 678 hfpp_n_simple(dcset,2) = ngoodtracks 679 endif
680 puckett 1.1.2.6
681 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 682 c on the track, then don't require six planes:
683 puckett 1.1.2.6
684 puckett 1.1.2.3 if(hfpp_prune_nplanes.eq.0.and.hfppfixleftright.gt.0) then 685 any6 = .false. 686 endif 687 688 c$$$ if(goodhms) then 689 c$$$ write(,) 690 c$$$ $ 'number of good tracks after tracking with drift=',ngoodtracks 691 c$$$ endif
692 puckett 1.1.2.4 693 c write(,) 'ngood drift tracks = ',ngoodtracks 694 695 if(ngoodtracks.eq.0) then 696 if(nhitsrequired.gt.hfpp_minsethits) then 697 nhitsrequired = nhitsrequired - 1 698 goto 134 699 else 700 exit 701 endif 702 endif
703 puckett 1.1.2.1 704 bestcandidate = 0 705 firsttry = .true.
706 puckett 1.1.2.3 firsttheta = .true. 707 firstsclose = .true. 708 * get minimum chi2, sclose, and theta, always checking for existence of six-hit tracks first:
709 puckett 1.1.2.4 do i=1,ngoodtracks 710 track=goodtracks(i) 711 chi2 = fulltracks(track,5) 712 sclose = sclosetracks(track) 713 theta = thetatracks(track) 714 if(candidate_nplanes(track).eq.6.or..not.any6) then 715 if(firsttry) then 716 firsttry = .false. 717 mintheta = theta 718 minchi2 = chi2 719 minsclose = sclose 720 else 721 if(theta.lt.mintheta) mintheta = theta 722 if(chi2.lt.minchi2) minchi2 = chi2 723 if(sclose.lt.minsclose) minsclose = sclose 724 endif 725 endif 726 enddo
727 puckett 1.1.2.3 firsttry = .true. 728 do i=1,ngoodtracks ! here is where we try to pick the best track: 729
730 puckett 1.1.2.1 track = goodtracks(i) 731 732 chi2 = fulltracks(track,5)
733 puckett 1.1.2.3 sclose = sclosetracks(track) 734 theta = thetatracks(track) 735 zclose = zclosetracks(track) 736 conetest = conetesttracks(track) 737 bestref = bestreftracks(track) 738 phi = phitracks(track) 739 740 c if(ngoodtracks.gt.1) then 741 c if(.not.any6.or.candidate_nplanes(track).eq.6) then 742 c$$$ if(goodhms) then 743 c$$$ write(,) 'track ',track,' theta=',theta180.0/PI, 744 c$$$ $ ' chi2=',chi2, 745 c$$$ $ ' sclose=',sclose,' zclose=',zclose,' phi=',phi180.0/PI, 746 c$$$ $ ' nplanes=',candidate_nplanes(track), 747 c$$$ $ ' nhits=',candidate_nhits(track), 748 c$$$ $ ' conetest=',conetest 749 c$$$ if(dcset.eq.2) write(,) 'bestref=',bestref 750 c$$$ endif 751 c endif 752 753 c try chi2 + sclose^2/sigmasclose^2 754 puckett 1.1.2.3 c criterion = chi2 + (sclosehfpp_sclose_weight(dcset))*2
755 puckett 1.1.2.1
756 puckett 1.1.2.3 c criterion = chi2 + minchi2(hfpp_sclose_weight(dcset)(sclose/minsclose)*2 + 757 c $ hfpp_theta_weight(dcset)prob(dcset,mintheta)/prob(dcset,theta)) 758 759 c if(dcset.eq.2.and.bestref.gt.0) then 760 c criterion = chi2 + minchi2(hfpp_sclose_weight(1)(sclose/minsclose)*2 + 761 c $ hfpp_theta_weight(1)prob(1,mintheta)/prob(1,theta)) 762 c endif 763
764 puckett 1.1.2.5 zgood = zclose.ge.hfpp_prune_zclose(2(dcset-1)+1) - 765 $ hfpp_prune_zslop(dcset)/tan(theta) .and. 766 $ zclose.le.hfpp_prune_zclose(2(dcset-1)+2) + 767 $ hfpp_prune_zslop(dcset)/tan(theta) 768
769 puckett 1.1.2.3 c try picking smallest theta HERE instead:
770 puckett 1.1.2.4 if(hselectfpptrackprune.eq.1) then 771 criterion = theta
772 puckett 1.1.2.11 any_great = .false.
773 puckett 1.1.2.4 else if(hselectfpptrackprune.eq.2) then 774 criterion = chi2 + minchi2 * sclose**2 775 else if(hselectfpptrackprune.eq.3) then 776 criterion = chi2 777 else if(hselectfpptrackprune.eq.4) then 778 criterion = sclose 779 else 780 criterion = chi2 781 endif
782 puckett 1.1.2.3
783 puckett 1.1.2.11 any_great = .false. 784
785 puckett 1.1.2.6 c if(candidate_nplanes(track).eq.h_fpp_n_dcinset* 786 c $ h_fpp_n_dclayers.or..not.any6) then 787 if( greattrack(track).or. .not. any_great) then 788 if(candidate_nplanes(track).eq.6.or..not.any6) then
789 puckett 1.1.2.5 if(firsttry.or.criterion.lt.mincriterion) then 790 mincriterion = criterion 791 firsttry = .false. 792 bestcandidate = track 793 endif
794 puckett 1.1.2.1 endif 795 endif 796 enddo 797
798 puckett 1.1.2.3 c$$$ if(goodhms) then 799 c$$$ write(,) 'chosen track=',bestcandidate 800 c$$$ endif 801 if(hidFPP_trkrough(DCset,1).gt.0) call hf1(hidFPP_trkrough(DCset,1),SimpleTracks(bestcandidate,1),1.) !mx 802 if(hidFPP_trkrough(DCset,2).gt.0) call hf1(hidFPP_trkrough(DCset,2),SimpleTracks(bestcandidate,2),1.) !bx 803 if(hidFPP_trkrough(DCset,3).gt.0) call hf1(hidFPP_trkrough(DCset,3),SimpleTracks(bestcandidate,3),1.) !my 804 if(hidFPP_trkrough(DCset,4).gt.0) call hf1(hidFPP_trkrough(DCset,4),SimpleTracks(bestcandidate,4),1.) !by 805 if(hidFPP_trkrough(DCset,5).gt.0) call hf1(hidFPP_trkrough(DCset,5),SimpleTracks(bestcandidate,5),1.) !chi2 806 if(hidFPP_trkrough(DCset,6).gt.0) call hf1(hidFPP_trkrough(DCset,6),SimpleTracks(bestcandidate,6),1.) !Nraw 807 808 if(hidFPP_roughchi2vsnhit(dcset).gt.0) call hf2(hidFPP_roughchi2vsnhit(dcset),simpletracks(bestcandidate,6), 809 $ simpletracks(bestcandidate,5),1.) 810
811 puckett 1.1.2.1 do i=1,ncandidates 812 do ichamber=1,h_fpp_n_dcinset 813 do ilayer=1,h_fpp_n_dclayers 814 icluster = hitcombos(i,ichamber,ilayer) 815 if(icluster.gt.0) then ! mark all candidate hits as unused 816 hfpp_clusterintrack(dcset,ichamber,ilayer,icluster) = 0
817 puckett 1.1.2.9 c reset drift time and drift distance for all wire hits in the cluster: 818 c If they are not used in a track, do not store their drift time and distance: 819 do ihit=1,hfpp_nhitsincluster(dcset,ichamber,ilayer,icluster) 820 hit = hfpp_clusters(dcset,ichamber,ilayer,icluster,ihit) 821 wire = hfpp_raw_wire(hit) 822 hfpp_drift_time(dcset,ichamber,ilayer,wire) = h_fpp_bad_time 823 hfpp_drift_dist(dcset,ichamber,ilayer,wire) = h_fpp_bad_drift 824 enddo
825 puckett 1.1.2.1 endif 826 enddo 827 enddo 828 enddo
829 puckett 1.1.2.3
830 puckett 1.1.2.1 c write(,) 'best candidate track, chi2=',bestcandidate,fulltracks(bestcandidate,5)
831 puckett 1.1.2.3 c write(,) 'nplanes,nhits=',candidate_nplanes(bestcandidate),
832 puckett 1.1.2.1 c $ candidate_nhits(bestcandidate) 833
834 puckett 1.1.2.3 c$$$ do i=1,4 835 c$$$ write(,) 'i,besttrack(i)=',i,fulltracks(bestcandidate,i) 836 c$$$ enddo
837 puckett 1.1.2.1 838 if(bestcandidate.eq.0) exit
839 puckett 1.1.2.8 840 c fit the track one more time to get the correct drift distance, as it may have changed: 841 c also, improve drift time/dist calculation using the full track: 842 do j=1,6
843 puckett 1.1.2.9 simpletrack(j) = simpletracks(bestcandidate,j)
844 puckett 1.1.2.8 c write(,) 'i,simpletrack(i)=',j,simpletrack(j) 845 enddo 846 847 do ichamber=1,h_fpp_n_dcinset 848 do ilayer=1,h_fpp_n_dclayers 849 bestclusters(ichamber,ilayer) = hitcombos(bestcandidate,ichamber,ilayer) 850 c write(,) 'chbr,pln,clstr=',ichamber,ilayer,bestclusters(ichamber,ilayer) 851 enddo 852 enddo 853
854 puckett 1.1.2.9 c write(,) 'before refitting, track=',(fulltracks(bestcandidate,j),j=1,6) 855
856 puckett 1.1.2.8 call h_fpp_tracking_drifttrack_ajp(dcset,simpletrack,bestclusters, 857 $ ontrack,track_good,fulltrack,nhitsrequired,1,abort,err) 858 859 if(.not. track_good) exit 860 861 do j=1,6 862 fulltracks(bestcandidate,j) = fulltrack(j) 863 enddo
864 puckett 1.1.2.9 865 c write(,) 'after refitting, track=',(fulltracks(bestcandidate,j),j=1,6)
866 puckett 1.1.2.1 867 c Now that we have found the best candidate track, add it to the good track array: 868 869 itrack = hfpp_n_tracks(dcset) + 1 870 871 if(itrack.le.h_fpp_max_tracks) then
872 puckett 1.1.2.3
873 puckett 1.1.2.1 hfpp_n_tracks(dcset) = itrack 874 875 hfpp_track_nlayers(dcset,itrack) = candidate_nplanes(bestcandidate)
876 puckett 1.1.2.8 c hfpp_track_nlayers(dcset,itrack) = 877 878 do j=1,4 879 hfpp_track_fine(dcset,itrack,j) = fulltracks(bestcandidate,j) 880 enddo 881
882 puckett 1.1.2.1 c Mark the hits in the best candidate track as used! 883 do ichamber = 1,h_fpp_n_dcinset 884 do ilayer = 1,h_fpp_n_dclayers 885 icluster = hitcombos(bestcandidate,ichamber,ilayer) 886 if(icluster.gt.0.and. 887 $ hitsintrack(bestcandidate,ichamber,ilayer) ) then 888 hfpp_clusterintrack(dcset,ichamber,ilayer,icluster) = itrack
889 puckett 1.1.2.14 hfpp_trackcluster(dcset,ichamber,ilayer,itrack) = icluster
890 puckett 1.1.2.1 endif 891 enddo 892 enddo 893 894 hfpp_track_chi2(dcset,itrack) = fulltracks(bestcandidate,5) 895 hfpp_track_nhits(dcset,itrack) = int(fulltracks(bestcandidate,6)) 896 897 do j=1,6 898 hfpp_track_rough(dcset,itrack,j) = simpletracks(bestcandidate,j) 899 enddo 900 901 hfpp_track_uniq(dcset,itrack) = .true. 902 * transform slopes from local FPP coordinates to HMS fp coordinates: 903 dccoords(1) = fulltracks(bestcandidate,1) ! dx/dz 904 dccoords(2) = fulltracks(bestcandidate,3) ! dy/dz 905 dccoords(3) = 1.0 ! dz/dz 906 907 call h_fpp_dc2fp(dcset,.true.,dccoords,fpcoords) 908
909 puckett 1.1.2.4 c$$$ hfpp_track_dx(dcset,itrack) = fpcoords(1) 910 c$$$ hfpp_track_dy(dcset,itrack) = fpcoords(2) 911 912 fpptrack(3) = fpcoords(1) ! x' in fp coords 913 fpptrack(4) = fpcoords(2) ! y' in fp coords
914 puckett 1.1.2.1 915 * now transform coordinates from local FPP system to HMS fp coords: 916 917 dccoords(1) = fulltracks(bestcandidate,2) ! x 918 dccoords(2) = fulltracks(bestcandidate,4) ! y 919 dccoords(3) = 0.0 920 921 call h_fpp_dc2fp(dcset,.false.,dccoords,fpcoords) 922
923 puckett 1.1.2.4 c$$$ hfpp_track_x(dcset,itrack) = fpcoords(1) - fpcoords(3)* 924 c$$$ $ hfpp_track_dx(dcset,itrack) 925 c$$$ hfpp_track_y(dcset,itrack) = fpcoords(2) - fpcoords(3)* 926 c$$$ $ hfpp_track_dy(dcset,itrack) 927 928 fpptrack(1) = fpcoords(1) - fpcoords(3) * fpptrack(3) ! x in fp coords at z=0 929 fpptrack(2) = fpcoords(2) - fpcoords(3) * fpptrack(4) ! y in fp coords at z=0 930 931 call h_fpp_align(dcset,fpptrack,newfpptrack) 932 933 hfpp_track_x(dcset,itrack) = newfpptrack(1) 934 hfpp_track_y(dcset,itrack) = newfpptrack(2) 935 hfpp_track_dx(dcset,itrack) = newfpptrack(3) 936 hfpp_track_dy(dcset,itrack) = newfpptrack(4)
937 puckett 1.1.2.1 938 * calculate relative scattering angles theta and phi: 939 940 call h_fpp_relative_angles(hsxp_fp,hsyp_fp,hfpp_track_dx(dcset,itrack), 941 $ hfpp_track_dy(dcset,itrack),theta,phi) 942 943 hfpp_track_theta(dcset,itrack) = theta 944 hfpp_track_phi(dcset,itrack) = phi 945
946 puckett 1.1.2.4 * calculate closest approach parameters with respect to HMS, REGARDLESS of which FPP!
947 puckett 1.1.2.1 948 hmstrack(1) = hsxp_fp 949 hmstrack(2) = hsx_fp 950 hmstrack(3) = hsyp_fp 951 hmstrack(4) = hsy_fp 952 953 fpptrack(1) = hfpp_track_dx(dcset,itrack) 954 fpptrack(2) = hfpp_track_x(dcset,itrack) 955 fpptrack(3) = hfpp_track_dy(dcset,itrack) 956 fpptrack(4) = hfpp_track_y(dcset,itrack) 957 958 call h_fpp_closest(hmstrack,fpptrack,sclose,zclose) 959 960 hfpp_track_sclose(dcset,itrack) = sclose 961 hfpp_track_zclose(dcset,itrack) = zclose 962 963 conetest = 1 964 965 call h_fpp_conetest(hmstrack,dcset,zclose,theta,conetest) 966 967 hfpp_track_conetest(dcset,itrack) = conetest 968 puckett 1.1.2.1 969 if(dcset.eq.2) then
970 puckett 1.1.2.13 971 c firsttry=.true. 972 c bestref = 0 973 974 c mintheta = hfpp_track_theta(dcset,itrack) 975 c do jtrack=1,hfpp_n_tracks(1)
976 puckett 1.1.2.6 c how this works: always figure out the "best reference" track for this FPP2 track in FPP1: 977 c after that, make a subjective judgement, based on theta, of whether or not to compare this track 978 c to FPP1 or the HMS:
979 puckett 1.1.2.4 if(hfpp_n_tracks(1).gt.0) then
980 puckett 1.1.2.13 c do jtrack = 1, hfpp_n_tracks(1) 981 jtrack = hfpp_best_track(1) 982 983 c mintheta = hfpp_track_theta(dcset,itrack) 984 985 call h_fpp_relative_angles(hfpp_track_dx(1,jtrack), 986 $ hfpp_track_dy(1,jtrack),
987 puckett 1.1.2.1 $ hfpp_track_dx(dcset,itrack),
988 puckett 1.1.2.13 $ hfpp_track_dy(dcset,itrack), 989 $ theta,phi) 990 c$$$ if(theta.lt.mintheta) then 991 c$$$ bestref = jtrack 992 c$$$ mintheta = theta 993 c$$$ endif 994 c enddo 995 c endif 996 997 c if(bestref.gt.0) then 998 c call h_fpp_relative_angles(hfpp_track_dx(1,bestref), 999 c $ hfpp_track_dy(1,bestref), 1000 c $ hfpp_track_dx(dcset,itrack), 1001 c $ hfpp_track_dy(dcset,itrack),theta,phi)
1002 puckett 1.1.2.1 hfpp_track_theta(dcset+1,itrack) = theta 1003 hfpp_track_phi(dcset+1,itrack) = phi 1004
1005 puckett 1.1.2.13 hmstrack(1) = hfpp_track_dx(1,jtrack) 1006 hmstrack(2) = hfpp_track_x(1,jtrack) 1007 hmstrack(3) = hfpp_track_dy(1,jtrack) 1008 hmstrack(4) = hfpp_track_y(1,jtrack)
1009 puckett 1.1.2.1 1010 fpptrack(1) = hfpp_track_dx(dcset,itrack) 1011 fpptrack(2) = hfpp_track_x(dcset,itrack) 1012 fpptrack(3) = hfpp_track_dy(dcset,itrack) 1013 fpptrack(4) = hfpp_track_y(dcset,itrack) 1014 1015 call h_fpp_closest(hmstrack,fpptrack,sclose,zclose) 1016 1017 hfpp_track_sclose(dcset+1,itrack) = sclose 1018 hfpp_track_zclose(dcset+1,itrack) = zclose 1019 1020 conetest = 1 1021 1022 call h_fpp_conetest(hmstrack,dcset,zclose,theta,conetest) 1023 1024 hfpp_track_conetest(dcset+1,itrack) = conetest
1025 puckett 1.1.2.13 c endif 1026 hfpp2_best_reference(itrack) = jtrack
1027 puckett 1.1.2.6
1028 puckett 1.1.2.13 if(hfpp_track_theta(dcset,itrack).lt. 1029 $ hfpp_track_theta(dcset+1,itrack)) hfpp2_best_reference(itrack) = 0 1030 1031 endif
1032 puckett 1.1.2.4
1033 puckett 1.1.2.13 endif
1034 puckett 1.1.2.3 if(hfppfixleftright.gt.0) 1035 $ call h_fpp_check_leftright(dcset,itrack) 1036 1037 endif ! end filling track common block variables
1038 puckett 1.1.2.1 c NOW we need to call freehitcount and start over with the remaining free hits: 1039 1040 call h_fpp_tracking_freehitcount(dcset,sufficient_hits) 1041 1042 enddo 1043 1044 return 1045 end 1046 1047 subroutine h_fpp_tracking_simple_ajp(dcset,hitcombos,simpletracks,
1048 puckett 1.1.2.4 $ ncandidates,nhitsrequired,abort,err)
1049 puckett 1.1.2.1 1050 implicit none 1051 save 1052 1053 INCLUDE 'hms_data_structures.cmn' 1054 INCLUDE 'hms_geometry.cmn' 1055 include 'gen_detectorids.par' 1056 include 'gen_decode_common.cmn' 1057 INCLUDE 'hms_fpp_event.cmn' 1058 INCLUDE 'hms_fpp_params.cmn' 1059 1060 character21 here 1061 parameter (here= 'h_fpp_tracking_simple') 1062 1063 integer4 dcset 1064 integer4 hitcombos(h_fpp_max_candidates,h_fpp_n_dcinset,h_fpp_n_dclayers) 1065 real4 simpletracks(h_fpp_max_candidates,6) 1066
1067 puckett 1.1.2.14 integer4 temphitcombo(h_fpp_n_dcinset,h_fpp_n_dclayers) 1068 integer4 tempsimpletrack(6) 1069
1070 puckett 1.1.2.1 integer4 HitCluster(H_FPP_N_DCINSET,H_FPP_N_DCLAYERS) 1071 1072 integer4 ncandidates,npoints
1073 puckett 1.1.2.14 integer*4 nhitsrequired,iterations,ichamber,ilayer,i,j,k,l,m,n
1074 puckett 1.1.2.1 integer*4 nhitsintrack
1075 puckett 1.1.2.3 integer*4 n3hit,n2hit ! tabulate number of three-hit and two hit clusters.
1076 puckett 1.1.2.1 1077 logical abort 1078 character() err 1079 1080 real*4 temptrack(5) ! not including hit count 1081 1082 abort = .false. 1083 err= ' ' 1084 1085 do i=1,h_fpp_max_candidates 1086 do ichamber=1,h_fpp_n_dcinset 1087 do ilayer=1,h_fpp_n_dclayers 1088 hitcombos(i,ichamber,ilayer) = 0 1089 enddo 1090 enddo 1091 do j=1,5 1092 simpletracks(i,j) = h_fpp_bad_coord 1093 enddo 1094 simpletracks(i,6) = 0.0 1095 enddo 1096 1097 puckett 1.1.2.1 nhitsrequired = h_fpp_n_planes + 1 1098 iterations = 0 1099
1100 puckett 1.1.2.14 c do while(nhitsrequired .ge. hfpp_minsethits.and.ncandidates.lt. 1101 c $ h_fpp_max_candidates) 1102 do while(nhitsrequired.ge.hfpp_minsethits)
1103 puckett 1.1.2.1 1104 nhitsintrack = 0 1105 * call next combo with large nhitsrequired to get the first useful combo as well as the max 1106 * number of hits available 1107 ncandidates = 0 1108
1109 puckett 1.1.2.3 c$$$ do ichamber=1,h_fpp_n_dcinset 1110 c$$$ do ilayer=1,h_fpp_n_dclayers 1111 c$$$ hitcluster(ichamber,ilayer) = 0 1112 c$$$ enddo 1113 c$$$ enddo 1114
1115 puckett 1.1.2.1 call h_fpp_tracking_nexthitcombo(dcset,nhitsrequired,nhitsintrack,hitcluster) 1116 1117 * keep comparing permutations until all possibilities are tried: 1118
1119 puckett 1.1.2.14 c do while(nhitsintrack.gt.0.and.ncandidates.lt.h_fpp_max_candidates) 1120 do while(nhitsintrack.gt.0)
1121 puckett 1.1.2.1 1122 iterations = iterations + 1 1123
1124 puckett 1.1.2.3 if(iterations.gt.hfpp_maxcombos) then 1125 write(,) 'WARNING: max FPP hit combos reached.' 1126 write(,) 'consider raising the limit or using ' 1127 write(,) 'tighter raw timing cuts' 1128 exit 1129 endif
1130 puckett 1.1.2.1 1131 call h_fpp_fit_simple(dcset,hitcluster,npoints,temptrack,abort,err) 1132 1133 if (ABORT) then 1134 call g_add_path(here,err) 1135 return 1136 endif 1137
1138 puckett 1.1.2.3 n2hit = 0 1139 n3hit = 0 1140 1141 do ichamber=1,h_fpp_n_dcinset 1142 do ilayer=1,h_fpp_n_dclayers 1143 if(hitcluster(ichamber,ilayer).gt.0) then 1144 if(hfpp_nhitsincluster(dcset,ichamber,ilayer, 1145 $ hitcluster(ichamber,ilayer)).eq.2) n2hit = n2hit + 1 1146 if(hfpp_nhitsincluster(dcset,ichamber,ilayer,hitcluster(ichamber,ilayer)) 1147 $ .eq.3) n3hit = n3hit + 1 1148 endif 1149 enddo 1150 enddo 1151 * in contrast to Frank's algorithm, here we keep track of any hit combos which pass the "reasonable chi2" 1152 * criterion, and we don't choose a track until we look at the drift based tracking. 1153
1154 puckett 1.1.2.7 c$$$ if(temptrack(5).ge.0.0.and. 1155 c$$$ $ temptrack(5).le.hfpp_aok_chi2+ 1156 c$$$ $ (12.float(n2hit)+48.float(n3hit))/float(nhitsrequired-4)/ 1157 c$$$ $ float(nhitsrequired-4 + n2hit + 2*n3hit) 1158 c$$$ $ ) then ! add a new candidate track to the array: 1159 if(temptrack(5).ge.0.0.and.temptrack(5).le.hfpp_aok_chi2) 1160 $ then
1161 puckett 1.1.2.1 ncandidates = ncandidates + 1
1162 puckett 1.1.2.14 1163 do k=1,min(h_fpp_max_candidates,ncandidates) 1164 * make this array sorted by simple track chi2, so if there is overflow, 1165 * we will take the h_fpp_max_candidates combos with the smallest simple chi2: 1166 if(temptrack(5).lt.simpletracks(k,5).or. 1167 $ k.eq.min(h_fpp_max_candidates,ncandidates)) then 1168 c this hit combo either has a smaller chi2 than the kth combo or is the first 1169 c combo found or has a larger chis than all previous hit combos: 1170 c put this combo at position k, move kth combo to k+1, etc... 1171 c first, shift all combos from k to ncandidates-1 up by 1: 1172 do m=min(h_fpp_max_candidates-1,ncandidates-1),k,-1 1173 c start with last (mth) combo, work our way down to k: 1174 c copy mth combo into (m+1)th combo: 1175 do n=1,6 1176 simpletracks(m+1,n) = simpletracks(m,n) 1177 enddo 1178 do ichamber=1,h_fpp_n_dcinset 1179 do ilayer=1,h_fpp_n_dclayers 1180 hitcombos(m+1,ichamber,ilayer) = hitcombos(m,ichamber,ilayer) 1181 enddo 1182 enddo 1183 puckett 1.1.2.14 enddo 1184 c next, insert the new combo at position k: 1185 do ichamber=1,h_fpp_n_dcinset 1186 do ilayer=1,h_fpp_n_dclayers 1187 hitcombos(k,ichamber,ilayer) = hitcluster(ichamber,ilayer) 1188 enddo 1189 enddo 1190 do n=1,5 1191 simpletracks(k,n) = temptrack(n) 1192 enddo 1193 simpletracks(k,6) = npoints 1194 c break out of the loop when this condition is met! 1195 exit 1196 endif
1197 puckett 1.1.2.1 enddo
1198 puckett 1.1.2.14 1199 c$$$ do ichamber=1,h_fpp_n_dcinset 1200 c$$$ do ilayer=1,h_fpp_n_dclayers 1201 c$$$ hitcombos(ncandidates,ichamber,ilayer) = 1202 c$$$ $ hitcluster(ichamber,ilayer) 1203 c$$$ enddo 1204 c$$$ enddo 1205 c$$$ 1206 c$$$ do j=1,5 1207 c$$$ simpletracks(ncandidates,j) = temptrack(j) 1208 c$$$ enddo 1209 c$$$ simpletracks(ncandidates,6) = float(npoints)
1210 puckett 1.1.2.1 endif 1211 1212 call h_fpp_tracking_nexthitcombo(dcset,nhitsrequired,nhitsintrack,hitcluster) 1213 1214 enddo
1215 puckett 1.1.2.14 1216 ncandidates = min(h_fpp_max_candidates,ncandidates)
1217 puckett 1.1.2.1 1218 if(iterations.gt.hfpp_maxcombos) then 1219 ncandidates = 0 1220 exit 1221 endif 1222 1223 if(ncandidates.gt.0) exit 1224 1225 nhitsrequired = nhitsrequired - 1 1226 1227 enddo 1228 1229 return 1230 end
1231 puckett 1.1.2.3 1232 real4 function prob(chbr,angle) 1233 1234 integer4 chbr 1235 real4 angle 1236 c approximate cumulative probability that a track scatters in fpp with 1237 c an angle theta < angle 1238 c prob = min(1.0e-9,1.0 - .1662angle*(-0.2202)exp(-4.031angle0.981)) 1239 1240 if(chbr.eq.1) prob = exp(-2.946angle) 1241 else prob = exp(-3.485angle) 1242 1243 return 1244 end 1245 1246 subroutine h_fpp_check_leftright(ifpp,itrack) 1247 1248 implicit none 1249 save 1250 1251 include 'gen_detectorids.par' 1252 puckett 1.1.2.3 include 'gen_decode_common.cmn' 1253 include 'gen_event_info.cmn' 1254 INCLUDE 'hms_data_structures.cmn' 1255 INCLUDE 'hms_geometry.cmn' 1256 INCLUDE 'hms_fpp_event.cmn' 1257 INCLUDE 'hms_fpp_params.cmn' 1258 INCLUDE 'hms_id_histid.cmn' 1259 1260 c fortunately, this is made easier by the fact that we know which clusters are on a track and 1261 c we know which hits within a cluster are on a track from their drift distance variable, so 1262 c we can revisit the left-right combo here: 1263 1264 integer4 nambig,ambigcombos(64) 1265 integer4 sign,oldsign,j 1266 integer4 npoints,ncombos,point,ipoint,icombo,combo 1267 integer4 ifpp,itrack,track,jtrack 1268 integer4 ichamber,ilayer,iwire,wire 1269 integer4 icluster,cluster,ihit,hit 1270 integer4 oldbestcombo,newbestcombo 1271 real4 combochi2(64) 1272 real4 combotheta(64) 1273 puckett 1.1.2.3 real4 combophi(64) 1274 real4 combosclose(64) 1275 real4 combozclose(64) 1276 integer4 comboconetest(64) 1277 logical ambig(64),anyambig,first 1278 1279 real4 chi2,theta,phi,sclose,zclose 1280 real4 oldchi2,oldtheta,oldphi,oldsclose,oldzclose 1281 1282 integer4 oldconetest,conetest 1283 1284 real4 plusminustest(h_fpp_max_fitpoints) 1285 real4 plusminusbest(h_fpp_max_fitpoints) 1286 1287 integer4 chambers(h_fpp_max_fitpoints) 1288 integer4 layers(h_fpp_max_fitpoints) 1289 integer4 wires(h_fpp_max_fitpoints) 1290 1291 real4 points(h_fpp_max_fitpoints,2) ! w and z 1292 real4 drifts(h_fpp_max_fitpoints) ! drift distance 1293 real4 projects(h_fpp_max_fitpoints,2) ! Px and Py coefficients 1294 puckett 1.1.2.3 real4 sigma2s(h_fpp_max_fitpoints) ! sigma^2 for chi2 calc. 1295 real4 hitpos(h_fpp_max_fitpoints,2) ! arrays for track fitting routine 1296 1297 real4 dccoords(3),fpcoords(3)
1298 puckett 1.1.2.4 real*4 reftrack(4),fpptrack(4),hmstrack(4),newfpptrack(4)
1299 puckett 1.1.2.3 1300 real4 testtrack(5) 1301 1302 real4 minsclose,mintheta 1303 real4 minstest,minthetatest,maxambig,maxdtheta,maxthetadiff 1304 real4 wtrack,xtrack,ytrack,xptrack,yptrack,Px,Py,z 1305 real4 zmin,zmax 1306 1307 external jbit ! cernlib bit routine 1308 external jibset 1309 integer4 jbit 1310 integer4 jibset ! Declare to help f2c 1311 1312 hfpp_track_nambig(ifpp,itrack) = 0 1313 1314 if(itrack.le.0.or.itrack.gt.hfpp_n_tracks(ifpp).or. 1315 $ hfpp_track_nhits(ifpp,itrack).ne. 1316 $ hfpp_track_nlayers(ifpp,itrack)) return 1317 1318 initialize points on the track: 1319 1320 puckett 1.1.2.3 c write(,) 'event = ',gen_event_id_number 1321 1322 npoints = 0 1323 1324 oldbestcombo = 0 1325 1326 xptrack = hfpp_track_fine(ifpp,itrack,1) 1327 xtrack = hfpp_track_fine(ifpp,itrack,2) 1328 yptrack = hfpp_track_fine(ifpp,itrack,3) 1329 ytrack = hfpp_track_fine(ifpp,itrack,4) 1330 1331 do ichamber=1,h_fpp_n_dcinset 1332 do ilayer=1,h_fpp_n_dclayers 1333 icluster = hfpp_trackcluster(ifpp,ichamber,ilayer,itrack) 1334 1335 if(icluster.gt.0) then 1336 do ihit=1,hfpp_nhitsincluster(ifpp,ichamber,ilayer,icluster) 1337 hit = hfpp_clusters(ifpp,ichamber,ilayer,icluster,ihit) 1338 1339 wire = hfpp_raw_wire(hit) 1340 1341 puckett 1.1.2.3 if(hfpp_drift_dist(ifpp,ichamber,ilayer,wire).ne.h_fpp_bad_drift) 1342 $ then ! this hit is on the track: 1343 npoints = npoints + 1 1344 if(npoints.le.h_fpp_max_fitpoints) then 1345 points(npoints,1) = hfpp_layeroffset(ifpp,ichamber,ilayer) 1346 $ + hfpp_spacing(ifpp,ichamber,ilayer)float(wire) 1347 points(npoints,2) = hfpp_layerz(ifpp,ichamber,ilayer) 1348 sigma2s(npoints) = hfpp_resolution(ifpp,ichamber,ilayer) 1349 projects(npoints,1) = hfpp_direction(ifpp,ichamber,ilayer,1) 1350 projects(npoints,2) = hfpp_direction(ifpp,ichamber,ilayer,2) 1351 drifts(npoints) = abs(hfpp_drift_dist(ifpp,ichamber,ilayer,wire)) 1352 1353 chambers(npoints) = ichamber 1354 layers(npoints) = ilayer 1355 wires(npoints) = wire 1356 1357 Px = projects(npoints,1) 1358 Py = projects(npoints,2) 1359 1360 z = points(npoints,2) 1361 1362 puckett 1.1.2.3 wtrack = Px(xtrack + xptrackz) + Py(ytrack + yptrackz) 1363 1364 if(abs(wtrack - points(npoints,1) - drifts(npoints)) 1365 $ .le.abs(wtrack - points(npoints,1) + drifts(npoints)) ) 1366 $ then ! track crosses on + side of wire-->assume hit originally had a + sign! 1367 oldbestcombo = oldbestcombo + 2*(npoints-1) 1368 endif 1369 1370 endif 1371 endif 1372 enddo 1373 endif 1374 enddo 1375 enddo 1376 1377 if(npoints.ne.hfpp_track_nhits(ifpp,itrack)) return
1378 puckett 1.1.2.11 if(npoints.lt.6) return
1379 puckett 1.1.2.3 1380 oldchi2 = hfpp_track_chi2(ifpp,itrack) 1381 oldtheta = hfpp_track_theta(ifpp,itrack) 1382 oldphi = hfpp_track_phi(ifpp,itrack) 1383 oldsclose = hfpp_track_sclose(ifpp,itrack) 1384 oldzclose = hfpp_track_zclose(ifpp,itrack) 1385 oldconetest = hfpp_track_conetest(ifpp,itrack) 1386 1387 reftrack(1) = hsxp_fp 1388 reftrack(2) = hsx_fp 1389 reftrack(3) = hsyp_fp 1390 reftrack(4) = hsy_fp 1391 1392 if(ifpp.eq.2.and.hfpp2_best_reference(itrack).gt.0) then 1393 track = hfpp2_best_reference(itrack) 1394 1395 oldtheta = hfpp_track_theta(ifpp+1,itrack) 1396 oldphi = hfpp_track_phi(ifpp+1,itrack) 1397 oldsclose = hfpp_track_sclose(ifpp+1,itrack) 1398 oldzclose = hfpp_track_zclose(ifpp+1,itrack) 1399 oldconetest = hfpp_track_conetest(ifpp+1,itrack) 1400 puckett 1.1.2.3 1401 reftrack(1) = hfpp_track_dx(1,track) 1402 reftrack(2) = hfpp_track_x(1,track) 1403 reftrack(3) = hfpp_track_dy(1,track) 1404 reftrack(4) = hfpp_track_y(1,track) 1405 1406 endif 1407 1408 c now we have baseline against which to compare other left-right combos. 1409
1410 puckett 1.1.2.12 zmin = hfpp_prune_zclose(2*(ifpp-1)+1)
1411 puckett 1.1.2.3 zmax = hfpp_prune_zclose(2(ifpp-1)+2) 1412 1413 ncombos = 2npoints 1414 1415 first = .true. 1416 1417 newbestcombo = -1 1418 1419 maxthetadiff = 0.0 1420 1421 minsclose = oldsclose 1422 mintheta = oldtheta 1423 1424 nambig = 0 1425 1426 do icombo=0,ncombos-1 1427 1428 ambig(icombo+1) = .false. 1429 1430 if(icombo.ne.oldbestcombo) then 1431 do ipoint=1,npoints 1432 puckett 1.1.2.3 if(jbit(icombo,ipoint).eq.1) then 1433 plusminustest(ipoint) = 1.0 1434 else 1435 plusminustest(ipoint) = -1.0 1436 endif 1437 1438 hitpos(ipoint,1) = points(ipoint,1) + 1439 $ drifts(ipoint)plusminustest(ipoint) 1440 hitpos(ipoint,2) = points(ipoint,2) 1441 1442 enddo 1443 1444 call h_fpp_fit3d_ajp(npoints,hitpos,sigma2s,projects,testtrack) 1445 1446 chi2 = testtrack(5) 1447 1448 combochi2(icombo+1) = chi2 1449
1450 puckett 1.1.2.13 c$$$ if(chi2.le.hfpp_min_chi2.and.chi2-oldchi2.le.hfpp_ambig_chi2cut(1) 1451 c$$$ $ .and.chi2/oldchi2.le.hfpp_ambig_chi2cut(2)) then ! "ambiguity" 1452 1453 if(chi2.le.hfpp_min_chi2.and.chi2-oldchi2.le.hfpp_ambig_chi2cut(1) ) 1454 $ then
1455 puckett 1.1.2.3 1456 c first check whether any hits on this combo which disagree with the best chi2 combo have 1457 c "large" drift distances, i.e., big enough that changing the sign makes a non-trivial difference in track 1458 c reconstruction: 1459 1460 anyambig = .false. 1461 1462 maxambig = 0.0 1463 1464 do ipoint=1,npoints 1465 sign = jbit(icombo,ipoint) 1466 oldsign = jbit(oldbestcombo,ipoint) 1467 if(sign.ne.oldsign.and. 1468 $ abs(drifts(ipoint)).ge.1.25sqrt(sigma2s(ipoint))) 1469 $ then 1470 anyambig = .true. 1471 1472 if(2.0drifts(ipoint).gt.maxambig) then 1473 maxambig = 2.0drifts(ipoint) 1474 endif 1475 1476 puckett 1.1.2.3 endif 1477 enddo 1478 1479 c maximum angular displacement of new track from old track is approximately equal to atan(maxambig/dz) 1480 1481 maxdtheta = atan(maxambig/abs(points(npoints,2)-points(1,2)) ) 1482 c$$$ write(,*) 1483 c$$$ $ 'max. possible angular displacement this track=',maxdtheta 1484 1485 if(maxdtheta.gt.maxthetadiff) maxthetadiff = maxdtheta 1486 1487 if(anyambig) then 1488 ambig(icombo+1) = .true. 1489 c this combo is "ambiguous": similar chi2 to the best combo and differs in sign from the best 1490 c combo in at least one "large" drift distance hit. Calculate theta,phi,sclose,and zclose: 1491 1492 c first, transform from dc to fp coords: 1493 1494 nambig = nambig + 1 1495 ambigcombos(nambig) = icombo 1496 1497 puckett 1.1.2.3 dccoords(1) = testtrack(1) ! dx/dz 1498 dccoords(2) = testtrack(3) ! dy/dz 1499 dccoords(3) = 1.0 ! dz/dz 1500 1501 call h_fpp_dc2fp(ifpp,.true.,dccoords,fpcoords) 1502
1503 puckett 1.1.2.4 fpptrack(3) = fpcoords(1) ! dx/dz 1504 fpptrack(4) = fpcoords(2) ! dy/dz
1505 puckett 1.1.2.3 1506 dccoords(1) = testtrack(2) ! x 1507 dccoords(2) = testtrack(4) ! y 1508 dccoords(3) = 0.0 ! z 1509 1510 call h_fpp_dc2fp(ifpp,.false.,dccoords,fpcoords) 1511
1512 puckett 1.1.2.4 fpptrack(1) = fpcoords(1) - fpcoords(3)fpptrack(3) ! x - x'z 1513 fpptrack(2) = fpcoords(2) - fpcoords(3)fpptrack(4) ! y - y'z 1514 1515 call h_fpp_align(ifpp,fpptrack,newfpptrack) 1516 1517 fpptrack(1) = newfpptrack(3) 1518 fpptrack(2) = newfpptrack(1) 1519 fpptrack(3) = newfpptrack(4) 1520 fpptrack(4) = newfpptrack(2)
1521 puckett 1.1.2.3 1522 call h_fpp_relative_angles(reftrack(1),reftrack(3),fpptrack(1),fpptrack(3),theta,phi) 1523 call h_fpp_closest(reftrack,fpptrack,sclose,zclose) 1524 conetest = 1 1525 call h_fpp_conetest(reftrack,ifpp,zclose,theta,conetest) 1526 1527 if(theta.lt.mintheta) mintheta = theta 1528 if(sclose.lt.minsclose) minsclose = sclose 1529 1530 combotheta(icombo+1) = theta 1531 combophi(icombo+1) = phi 1532 combosclose(icombo+1) = sclose 1533 combozclose(icombo+1) = zclose 1534 comboconetest(icombo+1) = conetest 1535 1536 endif 1537 endif 1538 endif 1539 enddo 1540
1541 puckett 1.1.2.13 if(nambig.gt.0.and.(mintheta.lt.oldtheta.or.minsclose.lt.oldsclose) )
1542 puckett 1.1.2.3 $ then 1543 1544 c hfpp_track_nambig(ifpp,itrack) = nambig 1545 1546 c reset minsclose and oldsclose: 1547 if(mintheta.le.hfpp_ambig_smallthetatest) then ! track consistent with theta=0: 1548 c use smallest theta regardless of chi2,sclose,zclose: 1549 minthetatest = mintheta 1550 mintheta = oldtheta 1551 do icombo=1,nambig 1552 combo = ambigcombos(icombo)+1 1553 if(combotheta(combo).lt.mintheta) then 1554 newbestcombo = combo-1 1555 mintheta = combotheta(combo) 1556 endif 1557 enddo 1558 else ! track not consistent with theta=0: use sclose instead: 1559 minstest = hfpp_ambig_sclosetest 1560 do icombo=1,nambig 1561 combo = ambigcombos(icombo)+1
1562 puckett 1.1.2.12 1563 zmin = hfpp_prune_zclose(2(ifpp-1)+1) - hfpp_prune_zslop(ifpp)/tan(combotheta(combo)) 1564 zmax = hfpp_prune_zclose(2(ifpp-1)+2) + hfpp_prune_zslop(ifpp)/tan(combotheta(combo)) 1565
1566 puckett 1.1.2.3 if( (combosclose(combo)/oldsclose)2.lt.minstest .and. 1567 $ zmin.lt.combozclose(combo).and. 1568 $ zmax.gt.combozclose(combo) ) then 1569 minstest = (combosclose(combo)/oldsclose)2 1570 newbestcombo = combo-1 1571 endif 1572 enddo 1573 endif 1574 endif 1575 1576 if(newbestcombo.ge.0) then 1577 1578 hfpp_track_nambig(ifpp,itrack) = nambig 1579 c refit the track using the "new best combo" 1580 1581 c$$$ write(,) 'found new best combo FPP,itrack=',ifpp,itrack 1582 c$$$ write(,) 'nhits=',npoints 1583 c$$$ write(,) 'old combo, nambig = ',oldbestcombo,nambig 1584 c$$$ write(,) 'old chi2,sclose,zclose,theta,phi,conetest=',oldchi2,oldsclose, 1585 c$$$ $ oldzclose,oldtheta,oldphi,oldconetest 1586 1587 puckett 1.1.2.3 do ipoint=1,npoints 1588 if(jbit(newbestcombo,ipoint).eq.1) then 1589 plusminusbest(ipoint) = 1.0 1590 else 1591 plusminusbest(ipoint) = -1.0 1592 endif 1593 1594 ichamber = chambers(ipoint) 1595 ilayer = layers(ipoint) 1596 iwire = wires(ipoint) 1597 1598 hfpp_drift_dist(ifpp,ichamber,ilayer,iwire) = 1599 $ plusminusbest(ipoint)drifts(ipoint) 1600 1601 hitpos(ipoint,1) = points(ipoint,1) + 1602 $ plusminusbest(ipoint)drifts(ipoint) 1603 hitpos(ipoint,2) = points(ipoint,2) 1604 enddo 1605 1606 call h_fpp_fit3d_ajp(npoints,hitpos,sigma2s,projects,testtrack) 1607 c copy the chi2 of the new track to common block variables: 1608 puckett 1.1.2.3 hfpp_track_chi2(ifpp,itrack) = testtrack(5) 1609 1610 do j=1,4 1611 hfpp_track_fine(ifpp,itrack,j) = testtrack(j) 1612 enddo 1613 1614 c re-calculate everything that was modified by this routine: 1615 1616 dccoords(1) = testtrack(1) ! dx/dz 1617 dccoords(2) = testtrack(3) ! dy/dz 1618 dccoords(3) = 1.0 ! dz/dz 1619 1620 call h_fpp_dc2fp(ifpp,.true.,dccoords,fpcoords) 1621
1622 puckett 1.1.2.4 fpptrack(3) = fpcoords(1) 1623 fpptrack(4) = fpcoords(2)
1624 puckett 1.1.2.3 1625 dccoords(1) = testtrack(2) ! x 1626 dccoords(2) = testtrack(4) ! y 1627 dccoords(3) = 0.0 ! z 1628 1629 call h_fpp_dc2fp(ifpp,.false.,dccoords,fpcoords) 1630
1631 puckett 1.1.2.4 fpptrack(1) = fpcoords(1) - fpcoords(3)fpptrack(3) ! x - x'z 1632 fpptrack(2) = fpcoords(2) - fpcoords(3)fpptrack(4) ! y - y'z 1633 1634 call h_fpp_align(ifpp,fpptrack,newfpptrack) 1635 1636 fpptrack(1) = newfpptrack(3) 1637 fpptrack(2) = newfpptrack(1) 1638 fpptrack(3) = newfpptrack(4) 1639 fpptrack(4) = newfpptrack(2)
1640 puckett 1.1.2.3 1641 hfpp_track_dx(ifpp,itrack) = fpptrack(1) 1642 hfpp_track_x(ifpp,itrack) = fpptrack(2) 1643 hfpp_track_dy(ifpp,itrack) = fpptrack(3) 1644 hfpp_track_y(ifpp,itrack) = fpptrack(4) 1645 1646 hmstrack(1) = hsxp_fp 1647 hmstrack(2) = hsx_fp 1648 hmstrack(3) = hsyp_fp 1649 hmstrack(4) = hsy_fp 1650 1651 call h_fpp_relative_angles(hmstrack(1),hmstrack(3),fpptrack(1), 1652 $ fpptrack(3),theta,phi) 1653 1654 hfpp_track_theta(ifpp,itrack) = theta 1655 hfpp_track_phi(ifpp,itrack) = phi 1656 1657 call h_fpp_closest(hmstrack,fpptrack,sclose,zclose) 1658 1659 hfpp_track_sclose(ifpp,itrack) = sclose 1660 hfpp_track_zclose(ifpp,itrack) = zclose 1661 puckett 1.1.2.3 1662 conetest = 1 1663 call h_fpp_conetest(hmstrack,ifpp,zclose,theta,conetest) 1664 1665 hfpp_track_conetest(ifpp,itrack) = conetest 1666 1667 if(ifpp.eq.2.and.hfpp2_best_reference(itrack).gt.0) then 1668 c reftrack has already been initialized to the best track chosen from 1669 c FPP1:
1670 puckett 1.1.2.4 reftrack(1) = hfpp_track_dx(1,hfpp2_best_reference(itrack)) 1671 reftrack(2) = hfpp_track_x(1,hfpp2_best_reference(itrack)) 1672 reftrack(3) = hfpp_track_dy(1,hfpp2_best_reference(itrack)) 1673 reftrack(4) = hfpp_track_y(1,hfpp2_best_reference(itrack)) 1674
1675 puckett 1.1.2.3 call h_fpp_relative_angles(reftrack(1),reftrack(3),fpptrack(1), 1676 $ fpptrack(3),theta,phi) 1677 hfpp_track_theta(ifpp+1,itrack) = theta 1678 hfpp_track_phi(ifpp+1,itrack) = phi 1679 1680 call h_fpp_closest(reftrack,fpptrack,sclose,zclose)
1681 puckett 1.1.2.1
1682 puckett 1.1.2.3 hfpp_track_sclose(ifpp+1,itrack) = sclose 1683 hfpp_track_zclose(ifpp+1,itrack) = zclose 1684 1685 conetest = 1 1686 call h_fpp_conetest(reftrack,ifpp,zclose,theta,conetest) 1687 1688 hfpp_track_conetest(ifpp+1,itrack) = conetest 1689 1690 endif 1691 c$$$ write(,) 'new combo=',newbestcombo 1692 c$$$ write(,) 'new chi2,sclose,zclose,theta,phi,conetest=',testtrack(5),sclose,zclose,theta,phi,conetest 1693 1694 endif 1695 1696 return 1697 end

Analyzer/Replay: Mark Jones, Documents: Stephen Wood