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 character*14 here
15 parameter (here= 'h_fpp_tracking')
16
17 integer*4 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 integer*4 ngoodtracks
27 integer*4 ncandidates,i,j,k,ichamber,ilayer
28 integer*4 nplanestemp,icluster,bestcandidate,best6,best5
29 integer*4 itrack,track
30
31 integer*4 goodtracks(h_fpp_max_candidates)
32 real*4 simpletracks(h_fpp_max_candidates,6)
33 real*4 fulltracks(h_fpp_max_candidates,6)
34 real*4 thetatracks(h_fpp_max_candidates) ! store track polar angle relative to HMS
35 real*4 chi2, minchi2
36
37 logical firsttry,any6,anytheta ! at this stage of the tracking, apply a maximum theta cut on all candidate tracks
38 logical any6theta ! any tracks with six planes AND passing theta cut?
39
40 integer*4 hitcombos(h_fpp_max_candidates,h_fpp_n_dcinset,h_fpp_n_dclayers)
41
42 logical hitsintrack(h_fpp_max_candidates,h_fpp_n_dcinset,h_fpp_n_dclayers)
43 puckett 1.1.2.1 logical ontrack(h_fpp_n_dcinset,h_fpp_n_dclayers) ! keep track of whether a cluster ends up on a track
44
45 integer*4 candidate_nplanes(h_fpp_max_candidates)
46 integer*4 candidate_nhits(h_fpp_max_candidates)
47
48 real*4 SimpleTrack(6), FullTrack(7)
49
50 integer*4 BestClusters(H_FPP_N_DCINSET,H_FPP_N_DCLAYERS)
51
52 real*4 HMStrack(4)
53 real*4 FPPtrack(4)
54 real*4 DCcoords(3), FPcoords(3)
55
56 real*4 PI
57 parameter(PI=3.141592653)
58
59 real*4 theta,mintheta,phi,sclose,zclose
60 integer*4 conetest,bestref,jtrack
61
62 ABORT= .FALSE.
63 err= ' '
64 puckett 1.1.2.1
65 any_track = .false.
66 any_good = .false.
67 any_great = .false.
68
69 hfpp_n_tracks(dcset) = 0
70
71 call h_fpp_tracking_freehitcount(dcset,sufficient_hits)
72
73 do while (sufficient_hits .and. (HFPP_N_tracks(DCset).lt.H_FPP_MAX_TRACKS))
74
75 ncandidates = 0
76 any6 = .false.
77 anytheta = .false.
78 any6theta = .false.
79
80 c the kind of loop we want here goes until we either run out of hit combos for a simple track, or
81 c until we exceed the maximum number of candidate tracks:
82
83 c initialize local tracking arrays:
84
85 puckett 1.1.2.1 do i=1,h_fpp_max_candidates
86 candidate_nplanes(i) = 0
87 candidate_nhits(i) = 0
88 goodtracks(i) = 0
89 do ichamber=1,h_fpp_n_dcinset
90 do ilayer=1,h_fpp_n_dclayers
91 hitcombos(i,ichamber,ilayer) = 0
92 hitsintrack(i,ichamber,ilayer) = .false.
93 enddo
94 enddo
95 enddo
96
97 c write(*,*) 'Start of tracking, FPP,ntracks=',dcset,hfpp_n_tracks(dcset)
98
99 call h_fpp_tracking_simple_ajp(dcset,hitcombos,simpletracks,ncandidates,abort,err)
100 if (ABORT) then
101 call g_add_path(here,err)
102 return
103 endif
104
105 if(ncandidates.gt.hfpp_n_simple(dcset))
106 puckett 1.1.2.1 $ hfpp_n_simple(dcset) = ncandidates
107 c write(*,*) 'number of candidates after simple tracking=',ncandidates
108
109 c$$$ if (.true.) then
110 c$$$ if (int(SimpleTrack(6)).le.0) then
111 c$$$ if(hidFPP_trkrough(DCset,6).gt.0) call hf1(hidFPP_trkrough(DCset,6),0.,1.) !Nraw
112 c$$$ else
113 c$$$ if(hidFPP_trkrough(DCset,1).gt.0) call hf1(hidFPP_trkrough(DCset,1),SimpleTrack(1),1.) !mx
114 c$$$ if(hidFPP_trkrough(DCset,2).gt.0) call hf1(hidFPP_trkrough(DCset,2),SimpleTrack(2),1.) !bx
115 c$$$ if(hidFPP_trkrough(DCset,3).gt.0) call hf1(hidFPP_trkrough(DCset,3),SimpleTrack(3),1.) !my
116 c$$$ if(hidFPP_trkrough(DCset,4).gt.0) call hf1(hidFPP_trkrough(DCset,4),SimpleTrack(4),1.) !by
117 c$$$ if(hidFPP_trkrough(DCset,5).gt.0) call hf1(hidFPP_trkrough(DCset,5),SimpleTrack(5),1.) !chi2
118 c$$$ if(hidFPP_trkrough(DCset,6).gt.0) call hf1(hidFPP_trkrough(DCset,6),SimpleTrack(6),1.) !Nraw
119 c$$$ endif
120 c$$$ endif
121
122 * * quit trying to make more tracks if we are out of hits
123 * * or if we couldnt make a good one now
124 c if (int(SimpleTrack(6)).le.0) exit
125 c if (int(SimpleTrack(5)).eq.H_FPP_BAD_CHI2) exit
126
127 puckett 1.1.2.1 c any_track = .true.
128
129 if(ncandidates.eq.0) exit
130
131 ngoodtracks = 0
132
133 do i=1,ncandidates
134 FullTrack(1) = H_FPP_BAD_COORD ! mx
135 FullTrack(2) = H_FPP_BAD_COORD ! bx
136 FullTrack(3) = H_FPP_BAD_COORD ! my
137 FullTrack(4) = H_FPP_BAD_COORD ! by
138 FullTrack(5) = H_FPP_BAD_CHI2
139 FullTrack(6) = 0. ! number of hits
140 FullTrack(7) = 0. ! number of planes
141 c write(*,*) 'candidate ',i,' simple chi2=',simpletracks(i,5)
142
143 do j=1,6
144 simpletrack(j) = simpletracks(i,j)
145 c write(*,*) 'i,simpletrack(i)=',j,simpletrack(j)
146 enddo
147
148 puckett 1.1.2.1 do ichamber=1,h_fpp_n_dcinset
149 do ilayer=1,h_fpp_n_dclayers
150 bestclusters(ichamber,ilayer) = hitcombos(i,ichamber,ilayer)
151 c write(*,*) 'chbr,pln,clstr=',ichamber,ilayer,bestclusters(ichamber,ilayer)
152 enddo
153 enddo
154
155 call h_fpp_tracking_drifttrack_ajp(dcset,simpletrack,bestclusters,
156 $ ontrack,track_good,fulltrack,1,abort,err)
157
158
159 c write(*,*) 'track_good = ',track_good
160 if (ABORT) then
161 call g_add_path(here,err)
162 return
163 endif
164
165 if(track_good) then
166
167 dccoords(1) = fulltrack(1)
168 dccoords(2) = fulltrack(3)
169 puckett 1.1.2.1 dccoords(3) = 1.0
170
171 call h_fpp_dc2fp(dcset,.true.,dccoords,fpcoords)
172
173 call h_fpp_relative_angles(hsxp_fp,hsyp_fp,
174 $ fpcoords(1),fpcoords(2),theta,phi)
175
176 if(theta.le.hfpp_prune_thetamax(dcset)*PI/180.0)
177 $ anytheta = .true.
178
179 thetatracks(i) = theta
180
181 do j=1,6
182 fulltracks(i,j) = fulltrack(j)
183 enddo
184
185 nplanestemp = int(fulltrack(7))
186
187 candidate_nplanes(i) = nplanestemp
188 candidate_nhits(i) = int(fulltrack(6))
189
190 puckett 1.1.2.1 ngoodtracks = ngoodtracks + 1
191 goodtracks(ngoodtracks) = i
192
193 if(nplanestemp.eq.h_fpp_n_dcinset*h_fpp_n_dclayers)
194 $ any6 = .true.
195
196 if(theta.le.hfpp_prune_thetamax(dcset)*PI/180.0
197 $ .and.nplanestemp.eq.h_fpp_n_dcinset*h_fpp_n_dclayers)
198 $ any6theta = .true.
199
200 do ichamber=1,h_fpp_n_dcinset
201 do ilayer = 1,h_fpp_n_dclayers
202 if(ontrack(ichamber,ilayer))
203 $ hitsintrack(i,ichamber,ilayer) = .true.
204 enddo
205 enddo
206
207 endif
208 enddo
209
210 c write(*,*) 'number of good tracks after tracking with drift=',ngoodtracks
211 puckett 1.1.2.1
212 if(ngoodtracks.eq.0) exit
213
214 bestcandidate = 0
215 firsttry = .true.
216
217 do i=1,ngoodtracks
218 track = goodtracks(i)
219
220 chi2 = fulltracks(track,5)
221
222 c prioritize: tracks passing theta cut take precedence over six-plane tracks.
223 c we want tracks with six planes passing theta cut, but we will take a five-plane track
224 c passing theta over a six-plane track failing theta:
225 if(firsttry.or.chi2.lt.minchi2) then
226 if(any6theta) then ! there is at least one track passing six planes AND theta cut. Require both
227 if(candidate_nplanes(track).eq.h_fpp_n_dcinset*h_fpp_n_dclayers
228 $ .and.thetatracks(track).le.hfpp_prune_thetamax(dcset)*PI/180.0)
229 $ then
230 firsttry = .false.
231 minchi2 = chi2
232 puckett 1.1.2.1 bestcandidate = track
233 endif
234 else if(anytheta) then ! there is no track passing six planes AND theta, require theta cut
235 if(thetatracks(track).le.hfpp_prune_thetamax(dcset)*PI/180.0)
236 $ then
237 firsttry = .false.
238 minchi2 = chi2
239 bestcandidate = track
240 endif
241 else if(any6) then ! there is no track passing theta cut. Require six planes
242 if(candidate_nplanes(track).eq.h_fpp_n_dcinset*h_fpp_n_dclayers)
243 $ then
244 firsttry = .false.
245 minchi2 = chi2
246 bestcandidate = track
247 endif
248 else ! there is no track passing theta cut OR six planes. take best chi2
249 firsttry = .false.
250 minchi2 = chi2
251 bestcandidate = track
252 endif
253 puckett 1.1.2.1 endif
254 enddo
255
256 do i=1,ncandidates
257 do ichamber=1,h_fpp_n_dcinset
258 do ilayer=1,h_fpp_n_dclayers
259 icluster = hitcombos(i,ichamber,ilayer)
260 if(icluster.gt.0) then ! mark all candidate hits as unused
261 hfpp_clusterintrack(dcset,ichamber,ilayer,icluster) = 0
262 endif
263 enddo
264 enddo
265 enddo
266
267 c write(*,*) 'best candidate track, chi2=',bestcandidate,fulltracks(bestcandidate,5)
268 c write(*,*) 'nplanes,nhits=',candidate_nplanes(bestcandidate),
269 c $ candidate_nhits(bestcandidate)
270
271 c$$$ do i=1,4
272 c$$$ write(*,*) 'i,besttrack(i)=',i,fulltracks(bestcandidate,i)
273 c$$$ enddo
274 puckett 1.1.2.1
275 if(bestcandidate.eq.0) exit
276
277 c Now that we have found the best candidate track, add it to the good track array:
278
279 itrack = hfpp_n_tracks(dcset) + 1
280
281 if(itrack.le.h_fpp_max_tracks) then
282
283 hfpp_n_tracks(dcset) = itrack
284
285 hfpp_track_nlayers(dcset,itrack) = candidate_nplanes(bestcandidate)
286 c Mark the hits in the best candidate track as used!
287 do ichamber = 1,h_fpp_n_dcinset
288 do ilayer = 1,h_fpp_n_dclayers
289 icluster = hitcombos(bestcandidate,ichamber,ilayer)
290 if(icluster.gt.0.and.
291 $ hitsintrack(bestcandidate,ichamber,ilayer) ) then
292 hfpp_clusterintrack(dcset,ichamber,ilayer,icluster) = itrack
293 endif
294 hfpp_trackcluster(dcset,ichamber,ilayer,itrack) = icluster
295 puckett 1.1.2.1 enddo
296 enddo
297
298 do j=1,4
299 hfpp_track_fine(dcset,itrack,j) = fulltracks(bestcandidate,j)
300 enddo
301
302 hfpp_track_chi2(dcset,itrack) = fulltracks(bestcandidate,5)
303 hfpp_track_nhits(dcset,itrack) = int(fulltracks(bestcandidate,6))
304
305 do j=1,6
306 hfpp_track_rough(dcset,itrack,j) = simpletracks(bestcandidate,j)
307 enddo
308
309 hfpp_track_uniq(dcset,itrack) = .true.
310 * transform slopes from local FPP coordinates to HMS fp coordinates:
311 dccoords(1) = fulltracks(bestcandidate,1) ! dx/dz
312 dccoords(2) = fulltracks(bestcandidate,3) ! dy/dz
313 dccoords(3) = 1.0 ! dz/dz
314
315 call h_fpp_dc2fp(dcset,.true.,dccoords,fpcoords)
316 puckett 1.1.2.1
317 hfpp_track_dx(dcset,itrack) = fpcoords(1)
318 hfpp_track_dy(dcset,itrack) = fpcoords(2)
319
320 * now transform coordinates from local FPP system to HMS fp coords:
321
322 dccoords(1) = fulltracks(bestcandidate,2) ! x
323 dccoords(2) = fulltracks(bestcandidate,4) ! y
324 dccoords(3) = 0.0
325
326 call h_fpp_dc2fp(dcset,.false.,dccoords,fpcoords)
327
328 hfpp_track_x(dcset,itrack) = fpcoords(1) - fpcoords(3)*
329 $ hfpp_track_dx(dcset,itrack)
330 hfpp_track_y(dcset,itrack) = fpcoords(2) - fpcoords(3)*
331 $ hfpp_track_dy(dcset,itrack)
332
333 * calculate relative scattering angles theta and phi:
334
335 call h_fpp_relative_angles(hsxp_fp,hsyp_fp,hfpp_track_dx(dcset,itrack),
336 $ hfpp_track_dy(dcset,itrack),theta,phi)
337 puckett 1.1.2.1
338 hfpp_track_theta(dcset,itrack) = theta
339 hfpp_track_phi(dcset,itrack) = phi
340
341 * calculate closest approach parameters:
342
343 hmstrack(1) = hsxp_fp
344 hmstrack(2) = hsx_fp
345 hmstrack(3) = hsyp_fp
346 hmstrack(4) = hsy_fp
347
348 fpptrack(1) = hfpp_track_dx(dcset,itrack)
349 fpptrack(2) = hfpp_track_x(dcset,itrack)
350 fpptrack(3) = hfpp_track_dy(dcset,itrack)
351 fpptrack(4) = hfpp_track_y(dcset,itrack)
352
353 call h_fpp_closest(hmstrack,fpptrack,sclose,zclose)
354
355 hfpp_track_sclose(dcset,itrack) = sclose
356 hfpp_track_zclose(dcset,itrack) = zclose
357
358 puckett 1.1.2.1 conetest = 1
359
360 call h_fpp_conetest(hmstrack,dcset,zclose,theta,conetest)
361
362 hfpp_track_conetest(dcset,itrack) = conetest
363
364 if(dcset.eq.2) then
365
366 firsttry=.true.
367 bestref = 0
368
369 do jtrack=1,hfpp_n_tracks(1)
370 call h_fpp_relative_angles(hfpp_track_dx(1,jtrack),
371 $ hfpp_track_dy(1,jtrack),
372 $ hfpp_track_dx(dcset,itrack),
373 $ hfpp_track_dy(dcset,itrack),
374 $ theta,phi)
375 if(firsttry.or.theta.lt.mintheta) then
376 firsttry = .false.
377 mintheta = theta
378 bestref = jtrack
379 puckett 1.1.2.1 endif
380 enddo
381
382 hfpp2_best_reference(itrack) = bestref
383
384 if(bestref.gt.0) then
385 call h_fpp_relative_angles(hfpp_track_dx(1,bestref),
386 $ hfpp_track_dy(1,bestref),
387 $ hfpp_track_dx(dcset,itrack),
388 $ hfpp_track_dy(dcset,itrack),theta,phi)
389 hfpp_track_theta(dcset+1,itrack) = theta
390 hfpp_track_phi(dcset+1,itrack) = phi
391
392 hmstrack(1) = hfpp_track_dx(1,bestref)
393 hmstrack(2) = hfpp_track_x(1,bestref)
394 hmstrack(3) = hfpp_track_dy(1,bestref)
395 hmstrack(4) = hfpp_track_y(1,bestref)
396
397 fpptrack(1) = hfpp_track_dx(dcset,itrack)
398 fpptrack(2) = hfpp_track_x(dcset,itrack)
399 fpptrack(3) = hfpp_track_dy(dcset,itrack)
400 puckett 1.1.2.1 fpptrack(4) = hfpp_track_y(dcset,itrack)
401
402 call h_fpp_closest(hmstrack,fpptrack,sclose,zclose)
403
404 hfpp_track_sclose(dcset+1,itrack) = sclose
405 hfpp_track_zclose(dcset+1,itrack) = zclose
406
407 conetest = 1
408
409 call h_fpp_conetest(hmstrack,dcset,zclose,theta,conetest)
410
411 hfpp_track_conetest(dcset+1,itrack) = conetest
412 endif
413 endif
414 endif ! end filling track common block variables
415 c NOW we need to call freehitcount and start over with the remaining free hits:
416
417 call h_fpp_tracking_freehitcount(dcset,sufficient_hits)
418
419 enddo
420
421 puckett 1.1.2.1 return
422 end
423
424 subroutine h_fpp_tracking_simple_ajp(dcset,hitcombos,simpletracks,
425 $ ncandidates,abort,err)
426
427 implicit none
428 save
429
430 INCLUDE 'hms_data_structures.cmn'
431 INCLUDE 'hms_geometry.cmn'
432 include 'gen_detectorids.par'
433 include 'gen_decode_common.cmn'
434 INCLUDE 'hms_fpp_event.cmn'
435 INCLUDE 'hms_fpp_params.cmn'
436
437 character*21 here
438 parameter (here= 'h_fpp_tracking_simple')
439
440 integer*4 dcset
441 integer*4 hitcombos(h_fpp_max_candidates,h_fpp_n_dcinset,h_fpp_n_dclayers)
442 puckett 1.1.2.1 real*4 simpletracks(h_fpp_max_candidates,6)
443
444 integer*4 HitCluster(H_FPP_N_DCINSET,H_FPP_N_DCLAYERS)
445
446 integer*4 ncandidates,npoints
447 integer*4 nhitsrequired,iterations,ichamber,ilayer,i,j
448 integer*4 nhitsintrack
449
450 logical abort
451 character*(*) err
452
453 real*4 temptrack(5) ! not including hit count
454
455 abort = .false.
456 err= ' '
457
458 do i=1,h_fpp_max_candidates
459 do ichamber=1,h_fpp_n_dcinset
460 do ilayer=1,h_fpp_n_dclayers
461 hitcombos(i,ichamber,ilayer) = 0
462 enddo
463 puckett 1.1.2.1 enddo
464 do j=1,5
465 simpletracks(i,j) = h_fpp_bad_coord
466 enddo
467 simpletracks(i,6) = 0.0
468 enddo
469
470 nhitsrequired = h_fpp_n_planes + 1
471 iterations = 0
472
473 do while(nhitsrequired .ge. hfpp_minsethits.and.ncandidates.lt.
474 $ h_fpp_max_candidates)
475
476 nhitsintrack = 0
477 * call next combo with large nhitsrequired to get the first useful combo as well as the max
478 * number of hits available
479 ncandidates = 0
480
481 call h_fpp_tracking_nexthitcombo(dcset,nhitsrequired,nhitsintrack,hitcluster)
482
483 * keep comparing permutations until all possibilities are tried:
484 puckett 1.1.2.1
485 do while(nhitsintrack.gt.0.and.ncandidates.lt.h_fpp_max_candidates)
486
487 iterations = iterations + 1
488
489 if(iterations.gt.hfpp_maxcombos) exit
490
491 call h_fpp_fit_simple(dcset,hitcluster,npoints,temptrack,abort,err)
492
493 if (ABORT) then
494 call g_add_path(here,err)
495 return
496 endif
497 * in contrast to Frank's algorithm, here we keep track of any hit combos which pass the "reasonable chi2"
498 * criterion:
499
500 if(temptrack(5).ge.0.0.and.
501 $ temptrack(5).le.hfpp_aok_chi2) then ! add a new candidate track to the array:
502 ncandidates = ncandidates + 1
503 do ichamber=1,h_fpp_n_dcinset
504 do ilayer=1,h_fpp_n_dclayers
505 puckett 1.1.2.1 hitcombos(ncandidates,ichamber,ilayer) =
506 $ hitcluster(ichamber,ilayer)
507 enddo
508 enddo
509
510 do j=1,5
511 simpletracks(ncandidates,j) = temptrack(j)
512 enddo
513 simpletracks(ncandidates,6) = float(npoints)
514 endif
515
516 call h_fpp_tracking_nexthitcombo(dcset,nhitsrequired,nhitsintrack,hitcluster)
517
518 enddo
519
520 if(iterations.gt.hfpp_maxcombos) then
521 ncandidates = 0
522 exit
523 endif
524
525 if(ncandidates.gt.0) exit
526 puckett 1.1.2.1
527 nhitsrequired = nhitsrequired - 1
528
529 enddo
530
531 return
532 end
533
534
|
Return to
h_fpp_tracking_ajp.f
CVS log
Up to [HallC] / Analyzer / HTRACKING