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

Analyzer/Replay: Mark Jones, Documents: Stephen Wood