StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
tpxFCF_2D.cxx
1 #include <stdio.h>
2 #include <sys/types.h>
3 #include <string.h>
4 #include <stdlib.h>
5 #include <math.h>
6 
7 #include <rtsLog.h>
8 
9 #include <TPC/rowlen.h>
10 #include <DAQ_TPX/tpxCore.h>
11 #include <DAQ_TPX/tpxFCF_flags.h>
12 #include <DAQ_TPX/tpxGain.h>
13 
14 #include "tpxFCF_2D.h"
15 
16 #include <DAQ_READER/daq_dta_structs.h>
17 
18 #define CHECK_SANITY
19 
20 //#include <GL3/profiler.hh>
21 //PROFILER_DECLARE ;
22 #define PROFILER(x) (x)
23 
24 #define likely(x) __builtin_expect((x),1)
25 #define unlikely(x) __builtin_expect((x),0)
26 
27 
28 #define DTA(i,j) (dta + (dt_2)*(i) + (j))
29 #define DTA_T(i,j) (dta_t + (dt_2)*(i) + (j))
30 #define DTA_S(i,j) (dta_s + (dt_2)*(i) + (j))
31 #define DTA_ID(i,j) (dta_id + (dt_2)*(i) + (j))
32 
33 static inline u_int get10(u_int *l, u_int p)
34 {
35  u_int ret ;
36 
37  l -= 2*(p/4) ;
38 
39  switch(p&3) {
40  case 0 :
41  ret = (l[-1] & 0xFFC00) >> 10 ;
42  break ;
43  case 1 :
44  ret = (l[-1] & 0x3FF) ;
45  break ;
46  case 2 :
47  ret = (l[0] & 0xFFC00) >> 10 ;
48  break ;
49  case 3 :
50  default:
51  ret = l[0] & 0x3FF ;
52  break ;
53  }
54 
55  //printf("P %d, data 0x%X\n",p,ret) ;
56 
57  return ret ;
58 }
59 
60 u_int *tpxFCF_2D_scan_to_next(tpxFCF_2D *fcf, u_int *end, u_int *start, tpx_altro_struct *a)
61 {
62  u_int *h = end ;
63  int wc, lo, hi ;
64 
65  a->count = 0 ;
66  a->err = 0 ;
67  a->row = 0 ;
68  a->pad = 0 ;
69  a->id = 0 ;
70  a->ch = 0 ;
71 
72  lo = *h-- ;
73  hi = *h-- ;
74 
75  if((lo & 0xCFF00000) || (hi & 0xCFF00000)) {
76  a->err = 1 ;
77  return 0 ;
78  }
79 
80 
81  if((hi & 0xFFFC0) != 0xAAA80) { // typical error for random data...
82  a->err = 2 ;
83  return 0 ;
84  }
85 
86  if((lo & 0x0F000) != 0x0A000) {
87  a->err = 3 ;
88  return 0 ;
89  }
90 
91  // need for monitoring!
92  a->id = (lo & 0xFF0) >> 4 ;
93  a->ch = lo & 0xF ;
94 
95 
96  int rrow, ppad ;
97 
98  tpx_from_altro(fcf->rdo-1,a->id,a->ch,rrow,ppad) ;
99 
100  // need for monitoring!
101  a->row = rrow ;
102  a->pad = ppad ;
103 
104  wc = ((hi&0x3F)<<4) | ((lo&0xF0000)>>16) ;
105 
106 
107  if((wc > 529) || (wc < 0)) {
108  a->err = 4 ;
109  return 0 ;
110  }
111 
112 
113 
114  if(wc==0) return h ;
115  if(fcf->is_pad_valid(rrow,ppad)==0) return h ;
116 
117 
118  int l10 = wc ;
119 
120  while(l10 % 4) l10++ ;
121 
122  int p10 = 0 ;
123 
124  switch(wc&3) {
125  case 0 :
126  break ;
127  case 1 :
128  p10 += 3 ;
129  break ;
130  case 2 :
131  p10 += 2 ;
132  break ;
133  case 3 :
134  p10++ ;
135  break ;
136  }
137 
138  *fcf->data_raw_p++ = rrow ;
139  *fcf->data_raw_p++ = ppad ;
140 
141  short *remember_data_raw_p = fcf->data_raw_p ;
142 
143  int tb_prev = 512 ;
144 
145  while(p10 < l10) {
146  int tb_cou = get10(h,p10++) ;
147  int tb_last = get10(h,p10++) ;
148 
149  tb_cou -= 2 ;
150 
151  if((tb_cou > wc) || (tb_cou <= 0) || (tb_cou > 512)) {
152  a->err = 10 ;
153  return 0 ;
154  }
155 
156  if((tb_last >= 512) || (tb_last < 0) || (tb_last >= tb_prev)) {
157  a->err = 11 ;
158  return 0 ;
159  }
160 
161  a->count += tb_cou ;
162  if((a->count >= 512)) {
163  a->err = 12 ;
164  return 0 ;
165  }
166 
167  tb_prev = tb_last - tb_cou ;
168 
169  *fcf->data_raw_p++ = tb_cou ;
170  *fcf->data_raw_p++ = -1 ;
171  *fcf->data_raw_p++ = tb_last ;
172 
173 
174  for(;tb_last>tb_prev;tb_last--) {
175 
176  u_short adc = get10(h,p10++) ;
177  *fcf->data_raw_p++ = adc ;
178  }
179 
180 
181 
182  }
183 
184 
185  // all OK here
186  if(a->count) {
187  fcf->p_row_pad[rrow][ppad] = remember_data_raw_p ;
188  *fcf->data_raw_p++ = 0 ;
189  *fcf->data_raw_p = 0 ;
190  }
191 
192  l10 /= 2 ;
193 
194  h -= l10 ;
195 
196  return h ;
197 
198 }
199 
200 int tpxFCF_2D::do_print(int row)
201 {
202  printf("++++ Doing row %2d\n",row) ;
203 
204  for(int p=1;p<=tpx_rowlen[row];p++) {
205  short *b1 = p_row_pad[row][p] ;
206 
207  if(b1==0) {
208  printf("row %d, pad %d -- empty\n",row,p) ;
209  continue ;
210  }
211 
212  int seq_cou = 0 ;
213 
214  for(;;) {
215  short tb_cou = *b1++ ;
216 
217  if(tb_cou == 0) {
218  printf("row %d, pad %d -- done\n",row,p) ;
219  break ;
220  }
221 
222 
223  seq_cou++ ;
224 
225  short ix = *b1++ ;
226  short tb_hi = *b1++ ;
227 
228  printf("row %d, pad %d, tb_cou %d, ix %d, tb_hi %d\n",row,p,tb_cou,ix,tb_hi) ;
229 
230  for(int i=0;i<tb_cou;i++) {
231  printf(" adc %4d\n",*b1++) ;
232  }
233  }
234 
235  printf("Row %2d, pad %3d: seq_cou %d\n",row,p,seq_cou) ;
236  }
237 
238  return 0 ;
239 
240 }
241 
242 
243 int tpxFCF_2D::stage_2d(u_int *buff, int max_bytes)
244 {
245  u_int *locbuff = buff ;
246  u_int *startbuff = buff ;
247 
248 
249  int gprof0 = PROFILER(0) ;
250 
251  for(row=1;row<=45;row++) { // need row count here!!!
252 
253  int do_debug = 0 ;
254 
255 
256  // check to see if we have this row...
257  if(gain_storage[sector-1][row] ==0) continue ;
258 
259 
260 // LOG(TERR,"doing sector %d, rdo %d, row %d",sector,rdo,row) ;
261 
262  int gprof1 = PROFILER(gprof0) ;
263 
264 
265  struct {
266  u_short hi, lo ;
267  } late_merge[MAX_LATE_MERGE] ;
268 
269  int late_merge_cou = 0 ;
270 
271  u_int *row_cache = locbuff++ ; // reserve space
272  u_int *clust_count_cache = locbuff++ ;
273  u_int *outbuff = locbuff ;
274 
275  short blob_ix = 0 ;
276 
277  short b1_cou, t1_lo, t1_hi ;
278  short b2_cou, t2_lo, t2_hi ;
279 
280  short *ix1, *ix2 ;
281 
282 #ifdef DO_DBG1
283  printf("Doing row %2d\n",row) ;
284 
285  for(int p=1;p<=tpx_rowlen[row];p++) {
286  short *b1 = p_row_pad[row][p] ;
287 
288  if(b1==0) {
289  printf("row %d, pad %d -- empty\n",row,p) ;
290  continue ;
291  }
292 
293  int seq_cou = 0 ;
294 
295  for(;;) {
296  short tb_cou = *b1++ ;
297 
298  if(tb_cou == 0) {
299  #ifdef DO_DBG1
300  printf("row %d, pad %d -- done\n",row,p) ;
301  #endif
302  break ;
303  }
304 
305 
306  seq_cou++ ;
307 
308  #ifdef DO_DBG1
309  short ix = *b1++ ;
310  short tb_hi = *b1++ ;
311 
312  printf("row %d, pad %d, tb_cou %d, ix %d, tb_hi %d\n",row,p,tb_cou,ix,tb_hi) ;
313 
314  for(int i=0;i<tb_cou;i++) {
315  printf(" adc %4d\n",*b1++) ;
316  }
317  #else
318  b1 += tb_cou + 2 ;
319  #endif
320  }
321 
322  printf("Row %2d, pad %3d: seq_cou %d\n",row,p,seq_cou) ;
323  }
324 #endif
325 
326 
327  for(int p=1;p<tpx_rowlen[row];p++) { // the < is on purpose!
328  short *b1 = p_row_pad[row][p] ;
329 
330  if(b1==0) continue ; // pad missing?
331 
332 #ifdef CHECK_SANITY
333  int rr, aa, cc ;
334  if(rdo) {
335 
336  tpx_to_altro(row,p,rr,aa,cc) ;
337  if(rr != rdo) {
338  LOG(ERR,"Huh? rdo %d: rp %d:%d",rdo,row,p) ;
339  continue ;
340  }
341  }
342 
343 
344  if(get_static(row,p)->f & FCF_KILLED_PAD) {
345  LOG(ERR,"Killed pad in data: sec %d, rp %d:%d; count %d",sector,row,p,*b1) ;
346 
347  continue ;
348  }
349 #endif
350 
351 
352  for(;;) {
353 
354  b1_cou = *b1++ ; // count of adcs
355 
356  if(b1_cou == 0) break ; ;
357 
358  ix1 = b1++ ;
359  t1_hi = *b1++ ; // tb_hi ;
360  t1_lo = t1_hi - b1_cou + 1 ;
361 
362 
363  b1 += b1_cou ; // advance to next...
364 #ifdef DO_SIMULATION
365  b1 += b1_cou ; // to skip the track_ids as well
366 #endif
367 
368  short *b2 = p_row_pad[row][p+1] ;
369 
370  if(b2 == 0) { // no second pad?
371  goto sweep_unmerged ;
372  }
373 
374 #ifdef CHECK_SANITY
375  if(rdo) {
376  tpx_to_altro(row,p+1,rr,aa,cc) ;
377  if(rr != rdo) {
378  LOG(ERR,"Huh?") ;
379  goto sweep_unmerged ; // different RDO
380  }
381  }
382 
383  if(get_static(row,p+1)->f & FCF_KILLED_PAD) {
384  LOG(ERR,"Killed pad in data: sec %d, rp %d:%d; count %d",sector,row,p+1,*b2) ;
385  goto sweep_unmerged ;
386  }
387 #endif
388 
389  for(;;) {
390  b2_cou = *b2++ ;
391 
392  if(b2_cou == 0) break ;
393 
394  ix2 = b2++ ;
395  t2_hi = *b2++ ;
396  t2_lo = t2_hi - b2_cou + 1 ;
397 
398 
399 
400  b2 += b2_cou ;
401 
402 #ifdef DO_SIMULATION
403  b2 += b2_cou ; // to skip the track ids
404 #endif
405  int merge ;
406 
407  if(t1_lo > t2_hi) merge = 0 ;
408  else if(t2_lo > t1_hi) merge = 0 ;
409  else merge = 1 ;
410 
411  #ifdef DO_DBG1
412  printf("MMMM: %d: [%d,%d]%d vs. %d [%d,%d]%d %s\n",
413  p,t1_lo,t1_hi,*ix1,
414  p+1,t2_lo,t2_hi,*ix2,
415  merge?"-- merged":"") ;
416  #endif
417 
418 
419  if(merge) {
420  if(*ix1 >= 0) {
421  if(*ix2 >= 0) {
422  if(*ix1 != *ix2) {
423 
424  //printf("Late merge: r:p %d:%d, merges %d, ix1 %d, ix2 %d\n",row,p,late_merge_cou,*ix1,*ix2) ;
425  if(late_merge_cou >= MAX_LATE_MERGE) {
426  LOG(WARN,"Too many late merges %d/%d: rp %d:%d",late_merge_cou,MAX_LATE_MERGE,row,p) ;
427  do_debug = 1 ;
428 
429  #ifdef DO_DBG3
430  printf("***** WARN Too many late merges\n") ;
431  #endif
432  }
433  else {
434  late_merge[late_merge_cou].lo = *ix1 ;
435  late_merge[late_merge_cou].hi = *ix2 ;
436  late_merge_cou++ ;
437  }
438 
439  #if DO_DBG3
440  printf(" late merge pad %d: %d %d\n",p,*ix1,*ix2) ;
441  #endif
442 
443  }
444  }
445  else {
446  *ix2 = *ix1 ;
447  }
448  }
449  else {
450  if(*ix2 >= 0) {
451  *ix1 = *ix2 ;
452  }
453  else {
454  *ix1 = blob_ix ;
455  *ix2 = blob_ix ;
456  blobs[blob_ix].cou = 0 ;
457  blob_ix++ ;
458  }
459  }
460  #ifdef DO_DBG1
461  printf(" merge: %d: %d %d\n",p,*ix1,*ix2) ;
462  #endif
463  }
464 
465  }
466 
467  sweep_unmerged: ;
468 
469  // sweep unmerged sequences
470  if(*ix1 < 0) { // unassigned
471  *ix1 = blob_ix ;
472  blobs[blob_ix].cou = 0 ;
473  blob_ix++ ;
474 
475  #ifdef DO_DBG1
476  printf(" was unmerged: %d: %d\n",p,*ix1) ;
477  #endif
478  }
479 
480 
481  }
482 
483  }
484 
485  // do very last pad in padrow by hand...
486  // PROBLEM for row 8!
487  int p = tpx_rowlen[row] ;
488  short *b1 = p_row_pad[row][p] ;
489 
490  if(b1) {
491  int rr = rdo ;
492 #ifdef CHECK_SANITY
493  if(rdo) {
494  int aa, cc ;
495  tpx_to_altro(row,p,rr,aa,cc) ;
496  }
497 
498  if(rr != rdo) {
499  LOG(ERR,"Huh?") ;
500  } ;
501 #endif
502  if(rr==rdo) {
503  short b1_cou ;
504  short *ix1 ;
505 
506  for(;;) {
507 
508  b1_cou = *b1++ ; // count of adcs
509 
510  if(b1_cou == 0) break ; ;
511 
512  ix1 = b1++ ;
513 
514  if(*ix1 < 0) {
515  *ix1 = blob_ix ;
516  blobs[blob_ix].cou = 0 ;
517  blob_ix++ ;
518 
519  #ifdef DO_DBG1
520  printf(" last pad: %d: %d\n",p,*ix1) ;
521  #endif
522  }
523 
524  b1 += b1_cou + 1 ;
525 
526 #ifdef DO_SIMULATION
527  b1 += b1_cou ;
528 #endif
529  }
530  }
531 
532  }
533 
534  LOG(DBG,"Doing 1") ;
535  // end of row; asign blobs
536  int gprof2 = PROFILER(gprof1) ;
537 
538 
539  // first; wrap up the late merges
540 #define MAX_CHAIN_COU 32
541 #define MAX_CHAIN_IX 32
542  struct chain_t {
543  short ix[MAX_CHAIN_IX] ;
544  short cou ;
545  short leader ;
546  } chain[MAX_CHAIN_COU] ;
547 
548  int chain_cou = 0 ;
549 
550  for(int i=0;i<late_merge_cou;i++) {
551  int lo = late_merge[i].lo ;
552  int hi = late_merge[i].hi ;
553 
554  int got_lo = 0 ;
555  int got_hi = 0 ;
556 
557  for(int j=0;j<chain_cou;j++) {
558 
559  for(int k=0;k<chain[j].cou;k++) {
560  if(chain[j].ix[k] == lo) got_lo = j + 1 ;
561  if(chain[j].ix[k] == hi) got_hi = j + 1;
562  }
563 
564  if(got_hi || got_lo) break ;
565  }
566 
567  if(got_hi && got_lo) continue ;
568 
569  if(got_hi) {
570  chain[got_hi-1].ix[chain[got_hi-1].cou] = lo ;
571  chain[got_hi-1].cou++ ;
572 
573  if(chain[got_hi-1].cou >= MAX_CHAIN_IX) {
574  do_debug = 1 ;
575  LOG(WARN,"Too many hi") ;
576  goto start_chain ;
577  }
578 
579  }
580  else if(got_lo) {
581  chain[got_lo-1].ix[chain[got_lo-1].cou] = hi ;
582  chain[got_lo-1].cou++ ;
583 
584  if(chain[got_lo-1].cou >= MAX_CHAIN_IX) {
585  do_debug = 1 ;
586  LOG(WARN,"Too many lo") ;
587  goto start_chain ;
588  }
589 
590 
591  }
592  else {
593  if(chain_cou >= MAX_CHAIN_COU) {
594  LOG(WARN,"Too many chains") ;
595  goto start_chain ;
596  }
597  else {
598  chain[chain_cou].ix[0] = lo ;
599  chain[chain_cou].ix[1] = hi ;
600  chain[chain_cou].cou = 2 ;
601  chain_cou++ ;
602  }
603  }
604 
605  }
606 
607  start_chain: ;
608 
609  // find the chain leader: probably unnecessary, can use the 1st one...
610  for(int i=0;i<chain_cou;i++) {
611  short min = 0x7FFF ;
612  for(int j=0;j<chain[i].cou;j++) {
613  if(chain[i].ix[j] < min) {
614  min = chain[i].ix[j];
615  }
616  chain[i].leader = min ;
617  }
618  }
619 
620  #ifdef DO_DBG3
621  for(int i=0;i<chain_cou;i++) {
622  printf("chain %d: ",i) ;
623  for(int j=0;j<chain[i].cou;j++) {
624  printf("%2d ",chain[i].ix[j]) ;
625  }
626  printf("\n") ;
627  }
628  #endif
629 
630  if(blob_ix >= MAX_BLOB_COUNT) {
631  LOG(WARN,"Too many blobs %d/%d",blob_ix,MAX_BLOB_COUNT) ;
632  do_debug = 1 ;
633  blob_ix = MAX_BLOB_COUNT - 1 ; // cut it off by hand!
634  }
635 
636  LOG(DBG,"Doing 2") ;
637  int gprof3 = PROFILER(gprof2) ;
638 
639  for(int pad=1;pad<=tpx_rowlen[row];pad++) {
640  register short *b = p_row_pad[row][pad] ;
641 
642  if(b==0) continue ; // pad missing?
643 
644 #ifdef CHECK_SANITY
645  if(rdo) {
646  int rr, aa, cc ;
647  tpx_to_altro(row,pad,rr,aa,cc) ;
648  if(rr != rdo) {
649  LOG(ERR,"Huh?") ;
650  continue ;
651  }
652  }
653 #endif
654  short flags = get_static(row,pad)->f & 0x0FFF ;
655  flags &= ~FCF_ONEPAD ; // need to fix the default ala tpxFCF...
656  if(rdo==0) flags &= ~FCF_BROKEN_EDGE ;
657 
658  register short cou ;
659 
660  while((cou = *b)) {
661 // int prof4 = PROFILER(0) ;
662 
663  short *bs = b ;
664 
665  short ix = bs[1] ; // get the index out
666 
667  b += cou + 3 ; // advance pointer to next sequence
668 
669 #ifdef DO_SIMULATION
670  b += cou ;
671 #endif
672  // check for late merge via chains
673  for(int l=0;l<chain_cou;l++) {
674  for(int k=0;k<chain[l].cou;k++) {
675  if(ix == chain[l].ix[k]) {
676 
677  #ifdef DO_DBG3
678  printf(" pad %3d, tb_hi %3d: chain merge of %d into %d\n",pad,bs[2],ix,chain[l].leader) ;
679  #endif
680  ix = chain[l].leader ;
681  goto chain_done ;
682  }
683  }
684  }
685 
686  chain_done:;
687 
688 
689  if((ix < 0) || (ix > blob_ix)){
690  LOG(WARN,"Too many ix: r:p %d:%d: ix %d/%d",
691  row,pad,ix,blob_ix) ;
692  #ifdef DO_DBG
693  printf("***** WARN Too many ix: r:p %d:%d: ix %d/%d\n",
694  row,pad,ix,blob_ix) ;
695 
696  #endif
697  continue ;
698  }
699 
700 
701 
702  blob_t *bl = blobs + ix ;
703 
704  short seq_cou = bl->cou ;
705 
706  if(seq_cou >= MAX_SEQ_PER_BLOB) {
707  LOG(WARN,"Too many seqs: r:p %d:%d: ix %d: seq_cou %d/%d",row,pad,ix,seq_cou,MAX_SEQ_PER_BLOB) ;
708 
709  #ifdef DO_DBG
710  printf("***** WARN Too many seqs: r:p %d:%d: ix %d: %d/%d\n",
711  row,pad,ix,seq_cou,MAX_SEQ_PER_BLOB) ;
712 
713  #endif
714  continue ;
715  }
716 
717  if(seq_cou ==0) {
718  bl->flags = flags ;
719  }
720  else {
721  bl->flags |= flags ;
722  }
723 
724  bl->cou++ ;
725 
726 
727  blob_seq_t *seq = bl->seq + seq_cou ;
728 
729  seq->s = bs ;
730  seq->pad = pad ;
731 
732 
733 // PROFILER(prof4) ;
734  }
735 
736  }
737 
738 
739  // do stage 2: extract stuff out...
740  LOG(DBG,"Doing 3") ;
741  int gprof4 = PROFILER(gprof3) ;
742 
743 
744  // ************************ loop over all the blobs
745  for(int ix=0;ix<blob_ix;ix++) {
746  #ifdef DO_DBG3
747  printf("===> BLOB %2d: row %d, %d pads, flags 0x%X\n",ix,row,blobs[ix].cou,blobs[ix].flags) ;
748  #endif
749 
750  if(blobs[ix].cou == 0) continue ; // can only be from a late merge
751 
752  short flags = blobs[ix].flags ;
753 
754  short p1 = 1000 ;
755  short t1 = 1000 ;
756  short p2 = 0 ;
757  short t2 = 0 ;
758 
759  int prof1 = PROFILER(0) ;
760 
761  // get the extents while loooping over sequences/strips
762  for(int j=0;j<blobs[ix].cou;j++) {
763  short *s = blobs[ix].seq[j].s ;
764  short pad = blobs[ix].seq[j].pad ;
765 
766  if(pad < p1) p1 = pad ;
767  if(pad > p2) p2 = pad ;
768 
769  short tb_cou = *s++ ;
770 
771  s++ ; // skip ix
772 
773  int tb_hi = (int) *s++ ;
774 
775  int tb_lo = tb_hi - tb_cou + 1 ;
776 
777 #ifdef CHECK_SANITY
778  // cases exist when there's just 1 pixel on timebin 0 or last
779  if(tb_hi>=500 || (tb_lo > tb_hi)) {
780  LOG(ERR,"tb_cou %d, tb_hi %d, tb_lo %d",tb_cou,tb_hi,tb_lo) ;
781  }
782 #endif
783  if(tb_hi > t2) t2 = tb_hi ;
784  if(tb_lo < t1) t1 = tb_lo ;
785 
786  }
787 
788  // p1,p2 & t1,t2 are the real extents of the blob
789 
790  short dp = p2 - p1 + 1 ;
791  short dt = t2 - t1 + 1 ;
792  short dt_2 = dt + 2 ; // for indexing
793 
794 
795 #ifdef CHECK_SANITY
796  if((t1>500)||(t2>500)||(dp<=0)||(dt<=0)) {
797  LOG(ERR,"What dp %d, dt %d: %d:%d p, %d:%d t???",dp,dt,p1,p2,t1,t2) ;
798  }
799 #endif
800  // cut off blobs where dt is very small
801  if((dt <= 1) && do_cuts && !(flags & FCF_BROKEN_EDGE)) continue ;
802 
803  // ONEPAD cut but only if not before the trigger (aka tb=15, aka wire hits))
804  if(dp <= 1) {
805  if((t1>15) && do_cuts && !(flags & FCF_BROKEN_EDGE)) continue ;
806 
807  flags |= FCF_ONEPAD ;
808  }
809 
810 
811  int prof2 = PROFILER(prof1) ; // 0.1 us
812 
813 
814  // ************** copy over to local storage in a NxM matrix form
815 
816  // check total bytes
817  u_int tot_size = (dp+2)*dt_2*sizeof(short)*2 ;
818  if(tot_size > sizeof(dta)) {
819  LOG(WARN,"Cluster too big: sec %d:%d, row %d, dp %d, dt %d: %d/%d-- skipping",sector,rdo,row,dp,dt,tot_size,sizeof(dta)) ;
820  continue ;
821  }
822 
823  // clear local storage
824  memset(dta, 0, (dp+2)*dt_2*sizeof(short)*2) ; // clear both storages in one shot
825 
826  dta_s = dta + (dp+2)*dt_2 ; // for the filtered data
827 
828 #ifdef DO_SIMULATION
829  memset(dta_t,0,(dp+2)*dt_2*sizeof(short)) ;
830  memset(dta_id,0,(dp+2)*dt_2*sizeof(short)) ;
831 #endif
832 
833  int prof3 = PROFILER(prof2) ; // now 0.15us (was 0.47 us)
834 
835 
836 
837  for(int j=0;j<blobs[ix].cou;j++) {
838  short *s = blobs[ix].seq[j].s ;
839  short p0 = blobs[ix].seq[j].pad - p1 + 1; // start from 1
840 
841 
842  short tb_cou = *s++ ;
843 
844  s++ ; // skip ix
845 
846  int tb_hi = (int) *s++ ;
847 
848  int tb_lo = tb_hi - tb_cou + 1 ;
849 
850 
851  int tb = tb_hi - t1 + 1 ; // start from 1
852 
853  short *udd = DTA(p0,0) ;
854 
855 // this is the critical timing part and should not be present in real-time code!
856 #ifdef DO_SIMULATION
857  u_short *tdd = DTA_T(p0,0) ;
858 #endif
859  for(;tb_hi>=tb_lo;tb_hi--) {
860  short adc = *s++ ;
861 #ifdef DO_SIMULATION
862  u_short track_id = *s++ ;
863 
864  *(tdd + tb) = track_id ;
865 #endif
866 
867  *(udd + tb) = adc ;
868  tb-- ;
869 
870  }
871  }
872 
873  // dta starts at [1,1]
874 
875  int prof4 = PROFILER(prof3) ; // 0.17 us
876 
877 
878  // ******************* now do the filtering (average) to get rid of the noise for peak-finding...
879  for(int i=1;i<=(dp);i++) {
880  for(int j=1;j<=(dt);j++) {
881 
882  short sum = 0 ;
883 
884  int iy = j - 1 ;
885 
886 #if 1
887  for(int ii=-1;ii<=1;ii++) {
888  register short *udd = DTA(i+ii,iy) ;
889 
890  for(int jj=0;jj<3;jj++) {
891  sum += *udd++ ;
892  }
893  }
894 #else
895  register short *udd = DTA(i,iy) ;
896 
897  for(int jj=0;jj<3;jj++) {
898  sum += *udd++ ;
899  }
900 
901 #endif
902 
903  *DTA_S(i,j) = sum ;
904  }
905  }
906 
907  int prof5 = PROFILER(prof4) ; // 0.29 0.26 (0.36 us)
908 
909 
910  // ********************* find the count of peaks...
911  int peaks_cou = 0 ;
912 
913 #ifdef DO_DBG
914  int pix_cou = 0 ;
915 #endif
916  for(int i=1;i<=(dp);i++) {
917  for(int j=1;j<=(dt);j++) {
918  short adc ;
919  short adc_peak ;
920 
921  // ignore the peak location if the original ADC is smaller
922  adc_peak = *DTA(i,j) ;
923 
924 #ifdef DO_DBG
925  if(adc_peak) pix_cou++ ;
926 #endif
927 
928  if(unlikely(adc_peak < 1)) continue ;
929 
930  adc = *DTA_S(i,j) ;
931 
932  if(unlikely(adc < 1)) continue ; // this can cause the number of peaks to be 0...
933 
934  int iy = j - 1 ;
935 
936  for(int ii=-1;ii<=1;ii++) {
937  short *dd = DTA_S(i+ii,iy) ;
938 
939  for(int jj=0;jj<3;jj++) {
940  short s_adc = *dd++ ;
941 
942  if(likely(adc < s_adc)) goto skip_calc ;
943  if(unlikely(s_adc < 0)) goto skip_calc ;
944  }
945  }
946 
947  *DTA_S(i,j) = -100 ; // mark as used!
948 
949  peaks[peaks_cou].i = i ;
950  peaks[peaks_cou].j = j ;
951  peaks[peaks_cou].aux_flags = adc_peak ;
952  peaks_cou++ ;
953 
954  if(unlikely(peaks_cou >= MAX_PEAKS_PER_BLOB)) goto peaks_done ;
955 
956  j += 2 ;
957 
958  skip_calc:;
959 
960  }
961  }
962 
963  peaks_done:;
964 
965  int prof6 = PROFILER(prof5) ; // 0.49 0.53 us
966 
967  // It is possible that I didn't find any peaks in the blob:
968  // - peak finding can miss a wide peak
969  // - the peak(s) were beyond the cuts
970 
971  // if I see no peaks, just use the maximum as the starting point
972  if((peaks_cou==0)||(peaks_cou>=MAX_PEAKS_PER_BLOB)) {
973  //if((dp==1) && (dt==1)) ;
974 // LOG(WARN,"Peaks_cou = %d: dp %d, dt %d, flags 0x%X",peaks_cou,dp,dt,flags) ;
975 
976  #ifdef DO_DBG3
977  printf("Warning: peaks_cou %d vs %d\n",peaks_cou,MAX_PEAKS_PER_BLOB) ;
978  #endif
979  }
980 
981 
982 
983 
984 #ifdef DO_DBG
985 
986 if(peaks_cou <= 1) {
987  static u_int ccc ;
988 
989  int max_adc = 0 ;
990 
991  for(int i=1;i<=(dp);i++) {
992  for(int j=1;j<=(dt);j++) {
993  short adc ;
994 
995 
996  adc = *DTA(i,j) ;
997 
998  if(adc > max_adc) max_adc = adc ;
999  }
1000  }
1001 
1002  if(max_adc >= 10) {
1003  ccc++ ;
1004  for(int j=1;j<=(dt);j++) {
1005  int sum = 0 ;
1006  for(int i=1;i<=(dp);i++) {
1007  short adc ;
1008 
1009 
1010  adc = *DTA(i,j) ;
1011  sum += adc ;
1012 
1013  }
1014 
1015  printf("FFF %d %d %d %d %d %d %d\n",event,sector,row,p1,ccc,t1+j-1,sum) ;
1016  }
1017  }
1018 
1019 }
1020 
1021 
1022 if(1) {
1023 //if((peaks_cou<=1) || (t2==512)) {
1024  printf("+++++ sec %2d, row %2d: p [%d,%d], t [%d,%d]; peaks %d\n",sector,row,p1,p2,t1,t2,peaks_cou) ;
1025  printf("peaks: ") ;
1026  for(int i=0;i<peaks_cou;i++) {
1027  printf("[%d,%d] ",peaks[i].i,peaks[i].j) ;
1028  }
1029  printf("\n") ;
1030 
1031  printf(" ") ;
1032  for(int i=0;i<=(dp+1);i++) {
1033  printf("p%03d ",p1+i-1) ;
1034  }
1035  printf("\n") ;
1036 
1037  for(int j=(dt+1);j>=0;j--) {
1038  printf("t%03d ",t1+j-1) ;
1039 
1040  for(int i=0;i<=(dp+1);i++) {
1041  printf("%4d ",*DTA(i,j)) ;
1042  }
1043  printf("\n") ;
1044  }
1045 
1046 
1047  printf("##### p [%d,%d], t [%d,%d]; peaks %d\n",p1,p2,t1,t2,peaks_cou) ;
1048 
1049  printf(" ") ;
1050  for(int i=0;i<=(dp+1);i++) {
1051  printf("p%03d ",p1+i-1) ;
1052  }
1053  printf("\n") ;
1054 
1055  for(int j=(dt+1);j>=0;j--) {
1056  printf("t%03d ",t1+j-1) ;
1057 
1058  for(int i=0;i<=(dp+1);i++) {
1059  printf("%4d ",*DTA_S(i,j)) ;
1060  }
1061  printf("\n") ;
1062  }
1063 
1064 
1065  printf("\n") ;
1066 
1067 }
1068 #endif
1069 
1070 
1071  int prof7 = 0 ;
1072 
1073 
1074  blob_c.p1 = p1 ;
1075  blob_c.p2 = p2 ;
1076  blob_c.t1 = t1 ;
1077  blob_c.t2 = t2 ;
1078  blob_c.dp = dp ;
1079  blob_c.dt = dt ;
1080  blob_c.flags = flags ;
1081  blob_c.dt_2 = dt_2 ;
1082 
1083  peaks[0].p1 = p1 ;
1084  peaks[0].p2 = p2 ;
1085  peaks[0].t1 = t1 ;
1086  peaks[0].t2 = t2 ;
1087  peaks[0].flags = flags ;
1088 
1089  static u_int peak_counter ;
1090 
1091  if(peaks_cou <= 1 ) { // single peak, full mean!
1092 
1093  double f_charge = 0.0 ;
1094  double f_t_ave = 0.0 ;
1095  double f_p_ave = 0.0 ;
1096 
1097  int pix_cou = 0 ;
1098 
1099  peak_counter++ ;
1100 
1101  for(int i=1;i<=(dp);i++) { // data starts on 2
1102 
1103 
1104  int tb = t1 ;
1105 
1106  int i_charge = 0 ;
1107  int i_t_ave = 0 ;
1108 
1109  short *adc_p = DTA(i,1) ;
1110 
1111  for(int j=1;j<=(dt);j++) { // data starts on 2
1112 
1113 
1114  register int adc = *adc_p++ ;
1115 
1116  if(adc) pix_cou++ ;
1117 
1118 #ifdef DO_SIMULATION
1119  if(adc) *DTA_ID(i,j) = cluster_id ;
1120 #endif
1121  //if(adc && (row<=13)) printf("BLOB1 %d %d %d %d\n",peak_counter,i,j,adc) ;
1122 
1123  //if(adc==0) {
1124  // tb++ ;
1125  // continue ;
1126  //}
1127 
1128  #ifdef DO_DBG2
1129  printf("pad:tb %d:%d, i:j %d:%d = %d\n",pad,tb,i,j,adc) ;
1130  #endif
1131 
1132  i_charge += adc ;
1133  i_t_ave += (tb++) * adc ;
1134  }
1135 
1136 
1137  if(i_charge == 0) continue ;
1138 
1139 
1140  int pad = p1 + i - 1 ;
1141 
1142  double gain = get_static(row,pad)->g ;
1143  double t0 = get_static(row,pad)->t0 ;
1144 
1145 #ifdef CHECK_SANITY
1146  if(gain < 0.8) {
1147  LOG(ERR,"BAD Gain 0: sector %d, rdo %d: rp %d %d, flags 0x%X",sector,rdo,row,pad,flags) ;
1148  }
1149 #endif
1150 
1151  double corr_charge = i_charge * gain ;
1152 
1153  f_charge += corr_charge ;
1154  f_t_ave += i_t_ave * gain + t0 * corr_charge ;
1155  f_p_ave += pad * corr_charge ;
1156 
1157  }
1158 
1159 
1160  peaks[0].f_t_ave = f_t_ave ;
1161  peaks[0].f_p_ave = f_p_ave ;
1162  peaks[0].f_charge = f_charge ;
1163  peaks[0].pix_cou = pix_cou ;
1164 
1165 
1166 
1167  peaks_cou = 1 ; // force it!
1168 
1169 #ifdef DO_SIMULATION
1170  peaks[0].cluster_id = cluster_id++ ;
1171  do_track_id(peaks_cou) ;
1172 #endif
1173  outbuff += do_dump(0,outbuff) ;
1174 
1175  prof7 = PROFILER(prof6) ; // 0.19 us
1176 
1177  }
1178  else {
1179 
1180  do_peaks(peaks_cou) ;
1181 #ifdef DO_SIMULATION
1182  do_track_id(peaks_cou) ;
1183 #endif
1184 
1185  for(int p=0;p<peaks_cou;p++) {
1186  outbuff += do_dump(p,outbuff) ;
1187  }
1188 
1189  prof7 = PROFILER(prof6) ; // 3.3 us
1190 
1191  }
1192 
1193  PROFILER(prof7) ; // 0.12 us
1194  PROFILER(prof1) ; // 2.0 us total per blob
1195  }
1196 
1197 
1198 
1199  u_int cl_cou ;
1200  if(modes) {
1201  cl_cou = (outbuff - locbuff)/3 ;
1202  }
1203  else {
1204  cl_cou = (outbuff - locbuff)/2 ;
1205  }
1206 
1207 
1208  #ifdef DO_DBG
1209  //printf("Row %2d: %d hits\n",row,cl_cou) ;
1210  #endif
1211 
1212  if(cl_cou) {
1213  *row_cache = (FCF_2D_V_FY13 << 16) | row ;
1214  *clust_count_cache = cl_cou ;
1215  locbuff = outbuff ;
1216  }
1217  else { // roll back
1218  locbuff -= 2 ;
1219  }
1220 
1221 
1222  PROFILER(gprof4) ;
1223 
1224  PROFILER(gprof1) ; // start of row
1225 
1226  if(do_debug) do_print(row) ;
1227 
1228  } // end of all rows
1229 
1230  PROFILER(gprof0) ; // start of RDO
1231 
1232  return (u_int *)locbuff - (u_int *)startbuff ; // return number of intes stored
1233 }
1234 
1235 
1236 /*
1237  Used only in re-doing clusters later in offline or for simulated
1238  data.
1239 */
1240 int tpxFCF_2D::do_pad_2d(tpx_altro_struct *a, daq_sim_adc_tb *sim_adc)
1241 {
1242  int row = a->row ;
1243  int pad = a->pad ;
1244 
1245 
1246  if(is_pad_valid(row,pad)==0) {
1247  //LOG(WARN,"Sec %d, rp %d:%d -- not valid",sector,row,pad) ;
1248  return 0 ;
1249  }
1250 
1251  if(a->count == 0) return 0 ;
1252 
1253 
1254  // check size
1255  int shorts_so_far = data_raw_p - data_raw ;
1256  if(shorts_so_far >((data_raw_shorts*9)/10)) {
1257  LOG(ERR,"Too much data %d/%d",shorts_so_far,data_raw_shorts) ;
1258  return 0 ;
1259  }
1260 
1261 
1262  // NEED to sort in falling timebin!
1263 
1264  short *start = data_raw_p ;
1265 
1266  *data_raw_p++ = row ;
1267  *data_raw_p++ = pad ;
1268 
1269  p_row_pad[row][pad] = data_raw_p ; // remember pointer
1270 
1271  int tb_cou = 0 ;
1272  short *tb_cou_p = data_raw_p ;
1273  int tb_last = 600 ;
1274 
1275  for(int i=0;i<a->count;i++) {
1276  int tb = a->tb[i] ;
1277  int aa = a->adc[i];
1278 
1279 #ifdef DO_SIMULATION
1280  // every pixel has an associated track_id
1281  int track_id ;
1282  if(sim_adc && modes) track_id = sim_adc[i].track_id ;
1283  else track_id = 0xFFFF ;
1284 #endif
1285 
1286  if(tb >= tb_last) {
1287  LOG(ERR,"%d: %d %d",i,tb,aa) ;
1288  }
1289 
1290  if(tb != (tb_last-1)) { // new sequence
1291 
1292  *tb_cou_p = tb_cou ;
1293 
1294  tb_cou_p = data_raw_p++ ;
1295  *data_raw_p++ = -1 ;
1296  *data_raw_p++ = tb ;
1297 
1298  tb_cou = 1 ;
1299 
1300  *data_raw_p++ = aa ;
1301 #ifdef DO_SIMULATION
1302  *data_raw_p++ = track_id ;
1303 #endif
1304  }
1305  else {
1306 
1307  tb_cou++ ;
1308  *data_raw_p++ = aa ;
1309 #ifdef DO_SIMULATION
1310  *data_raw_p++ = track_id ;
1311 #endif
1312  }
1313 
1314  tb_last = tb ;
1315 
1316  }
1317 
1318 
1319  *tb_cou_p = tb_cou ;
1320 
1321 
1322  *data_raw_p++ = 0 ;
1323  *data_raw_p = 0 ;
1324 
1325  int s_cou = data_raw_p - start ;
1326 
1327  return s_cou ;
1328 
1329 }
1330 
1331 
1332 tpxFCF_2D::tpxFCF_2D()
1333 {
1334 #ifndef TPX_ONLINE
1335  LOG(WARN,"TPX_ONLINE undefined -- running in Offline mode.") ;
1336 #endif
1337 #ifdef DO_SIMULATION
1338  LOG(WARN,"DO_SIMULATION defined.") ;
1339 #endif
1340 #ifdef CHECK_SANITY
1341  LOG(WARN,"CHECK_SANITY defined.") ;
1342 #endif
1343 
1344  event = 0 ;
1345  data_raw_shorts = 2*MAX_PADS_PER_RDO*512 ; // we need this sizable because of simulated data with track_id!
1346  data_raw = (short *) valloc(2*data_raw_shorts) ; // pad_num X 512tb X 2bytes (short)
1347 }
1348 
1349 int tpxFCF_2D::do_dump(int ix, u_int *obuff)
1350 {
1351  int ret = 0 ; // how many ints did we use...
1352 
1353 
1354 
1355  double f_p_ave, f_t_ave ;
1356 
1357 
1358  short flags = peaks[ix].flags ;
1359  short p1 = peaks[ix].p1 ;
1360  short p2 = peaks[ix].p2 ;
1361  short t1 = peaks[ix].t1 ;
1362  short t2 = peaks[ix].t2 ;
1363 
1364  if(peaks[ix].f_charge) {
1365 
1366  f_p_ave = peaks[ix].f_p_ave / peaks[ix].f_charge ;
1367  f_t_ave = peaks[ix].f_t_ave / peaks[ix].f_charge ;
1368  }
1369  else {
1370  f_t_ave = 0 ;
1371  f_p_ave = 0 ;
1372  LOG(WARN,"no charge: ix %d: %d-%d, %d-%d: fla 0x%X: ij %d-%d, pixs %d, aux 0x%X",
1373  ix,
1374  p1,p2,t1,t2,
1375  flags,
1376  peaks[ix].i,peaks[ix].j,
1377  peaks[ix].pix_cou,
1378  peaks[ix].aux_flags) ;
1379  return 0 ;
1380  }
1381 
1382  u_int time_c = (u_int)(f_t_ave * 64.0 + 0.5) ;
1383  u_int pad_c = (u_int)(f_p_ave * 64.0 + 0.5) ;
1384  u_int cha = (u_int) (peaks[ix].f_charge + 0.5) ;
1385 
1386 
1387  if(flags & FCF_BROKEN_EDGE) goto keep ;
1388  if(do_cuts == 0) goto keep ;
1389 
1390  if(do_cuts != 2) {
1391  if(flags & (FCF_DEAD_EDGE | FCF_ROW_EDGE)) return 0 ;
1392  }
1393 
1394  if((t2-t1)<=2) return 0 ;
1395  if(cha < 10.0) return 0 ;
1396 
1397  keep:;
1398 
1399 #if 0
1400 if(flags & 2) ;
1401 else {
1402  if((f_t_ave > 26.0) && (f_t_ave < 415.0)) ;
1403  else {
1404  if(cha < 12) ;
1405  else printf("PROMPT %d %d %d %f %f %d\n",event,sector,row,f_p_ave,f_t_ave,cha) ;
1406  }
1407 }
1408 #endif
1409  u_int tmp_fl ;
1410 
1411  // get extents
1412  u_int tmp_p = pad_c / 64 ;
1413 
1414  p1 = tmp_p - p1 ;
1415  p2 = p2 - tmp_p ;
1416 
1417  if(p1 > 7) p1 = 7 ;
1418  if(p2 > 7) p2 = 7 ;
1419 
1420 
1421  if(p1 < 0) p1 = 0 ;
1422  if(p2 < 0) p2 = 0 ;
1423 
1424  tmp_fl = (p1 << 8) | (p2 << 11) ;
1425 
1426  tmp_p = time_c / 64 ;
1427 
1428  t1 = tmp_p - t1 ;
1429  t2 = t2 - tmp_p ;
1430 
1431 
1432  if(t1 > 15) t1 = 15 ;
1433  if(t2 > 15) t2 = 15 ;
1434 
1435  if(t1 < 0) t1 = 0 ;
1436  if(t2 < 0) t2 = 0 ;
1437 
1438  tmp_fl |= (t2 << 4) | t1 ;
1439 
1440  if(flags & FCF_ONEPAD) time_c |= 0x8000 ;
1441 
1442  if(flags & FCF_MERGED) pad_c |= 0x8000 ;
1443  if(flags & FCF_DEAD_EDGE) pad_c |= 0x4000 ;
1444 
1445  if(flags & FCF_ROW_EDGE) tmp_fl |= 0x8000 ;
1446  if(flags & FCF_BROKEN_EDGE) tmp_fl |= 0x4000 ;
1447 
1448  if(cha > 0x7FFF) cha = 0x8000 | (cha/1024) ;
1449 
1450  *obuff++ = (time_c << 16) | pad_c ;
1451  *obuff++ = (cha << 16) | tmp_fl ;
1452 
1453  ret = 2 ;
1454 
1455 #ifdef DO_SIMULATION
1456  if(modes) {
1457  int quality = peaks[ix].quality ;
1458  int track_id = peaks[ix].track_id ;
1459 
1460  *obuff++ = (quality << 16) | track_id ;
1461  ret++ ;
1462  }
1463 #else
1464  if(modes) {
1465  *obuff++ = 0 ;
1466  ret++ ;
1467  }
1468 #endif
1469 
1470  return ret ;
1471 }
1472 
1473 int tpxFCF_2D::do_peaks(int peaks_cou)
1474 {
1475  // ix==0 contains the extents of the blob...
1476  short p1 = blob_c.p1 ;
1477  short p2 = blob_c.p2 ;
1478  short t1 = blob_c.t1 ;
1479  short t2 = blob_c.t2 ;
1480 
1481  short dt, dt_2 ;
1482  short dp ;
1483 
1484 
1485 #ifdef CHECK_SANITY
1486  if((p1>p2)||(t1>t2)||(p1<1)||(t1<0)) {
1487  LOG(ERR,"What??? %d %d %d %d",p1,p2,t1,t2) ;
1488  }
1489 #endif
1490 
1491  dp = blob_c.dp ;
1492  dt = blob_c.dt ;
1493  dt_2 = blob_c.dt_2 ; // MUST be called dt_2 due to the DTA() define!!!
1494 
1495  int flags = 0 ;
1496  if(peaks_cou==0) peaks_cou = 1 ; // sanitize...
1497  if(peaks_cou > 1) flags = FCF_MERGED ;
1498 
1499 
1500 
1501  for(int p=0;p<peaks_cou;p++) {
1502  peaks[p].f_charge = 0.0 ;
1503  peaks[p].f_t_ave = 0.0 ;
1504  peaks[p].f_p_ave = 0.0 ;
1505  peaks[p].pix_cou = 0 ;
1506  peaks[p].t1 = 10000 ;
1507  peaks[p].p1 = 10000 ;
1508  peaks[p].t2 = 0 ;
1509  peaks[p].p2 = 0 ;
1510  peaks[p].flags = flags ;
1511 #ifdef DO_SIMULATION
1512  peaks[p].cluster_id++ ;
1513 #endif
1514  }
1515 
1516  for(int i=1;i<=(dp);i++) {
1517  int pad = p1 + i - 1 ;
1518 
1519  double gain = get_static(row,pad)->g ;
1520  double t0 = get_static(row,pad)->t0 ;
1521 
1522  short flags = get_static(row,pad)->f & 0xFFF ;
1523  if(rdo==0) flags &= ~FCF_BROKEN_EDGE ; // clear broken_edge if not running rdo-per-rdo
1524 
1525 #ifdef CHECK_SANITY
1526  if(gain < 0.8) {
1527  LOG(ERR,"What: lo gain %f %f, flags 0x%X: sec %d: rp %d:%d; i %d, dp %d",gain,t0,flags,sector,row,pad,i,dp) ;
1528  }
1529 #endif
1530 
1531  for(int j=1;j<=(dt);j++) {
1532 
1533  short adc = *DTA(i,j) ;
1534 
1535  if(adc==0) continue ;
1536 
1537  int min_dist = 10000000 ;
1538  int min_p = -1 ;
1539 
1540  // find the closest peak...
1541  for(int p=0;p<peaks_cou;p++) {
1542 
1543  int dx = i - peaks[p].i ;
1544  int dy = j - peaks[p].j ;
1545 
1546  int dist = dx*dx + dy*dy ;
1547 
1548  if(dist < min_dist) {
1549  min_dist = dist ;
1550  min_p = p ;
1551  }
1552  }
1553 
1554 
1555  if(min_p<0) {
1556  LOG(ERR,"What: %d: rp %d:%d",min_p,row,pad) ;
1557  continue ;
1558  }
1559 
1560 #ifdef DO_SIMULATION
1561  *DTA_ID(i,j) = peaks[min_p].cluster_id ;
1562 #endif
1563  double corr_charge = adc * gain ;
1564 
1565  peaks[min_p].f_charge += corr_charge ;
1566  peaks[min_p].f_t_ave += corr_charge*(t0+(double)(j+t1-1)) ;
1567  peaks[min_p].f_p_ave += corr_charge*(pad) ;
1568  peaks[min_p].pix_cou++ ;
1569 
1570 
1571  peaks[min_p].flags |= flags ;
1572 
1573  if(i < peaks[min_p].p1) peaks[min_p].p1 = i ;
1574  if(i > peaks[min_p].p2) peaks[min_p].p2 = i ;
1575  if(j < peaks[min_p].t1) peaks[min_p].t1 = j ;
1576  if(j > peaks[min_p].t2) peaks[min_p].t2 = j ;
1577  }
1578 
1579  }
1580 
1581  for(int p=0;p<peaks_cou;p++) {
1582  peaks[p].t1 += t1 - 1 ;
1583  peaks[p].p1 += p1 - 1 ;
1584  peaks[p].t2 += t1 - 1 ;
1585  peaks[p].p2 += p1 - 1 ;
1586 
1587  if((peaks[p].p2 - peaks[p].p1)==0) peaks[p].flags |= FCF_ONEPAD ;
1588  else peaks[p].flags &= ~FCF_ONEPAD ;
1589  }
1590 
1591  return 0 ;
1592 }
1593 
1594 
1595 void tpxFCF_2D::do_track_id(int peaks_cou)
1596 {
1597  short dt, dt_2 ;
1598  short dp ;
1599 
1600  dp = blob_c.dp ;
1601  dt = blob_c.dt ;
1602  dt_2 = blob_c.dt_2 ; // MUST be called dt_2 due to the DTA() define!!!
1603 
1604 
1605  for(int c=0;c<peaks_cou;c++) {
1606  // now to the track_id and Q
1607  const int MAX_TRACK_IDS = 10 ;
1608 
1609  struct {
1610  int cou ;
1611  u_int track_id ;
1612  } t_id[MAX_TRACK_IDS] ;
1613 
1614  memset(t_id,0,sizeof(t_id)) ;
1615 
1616 
1617  u_short track_id = 123;
1618  u_int blob_adc = 0 ;
1619 
1620  for(int i=1;i<=dp;i++) {
1621  for(int j=1;j<=dt;j++) {
1622  if(*DTA_ID(i,j) == peaks[c].cluster_id) {
1623  track_id = *DTA_T(i,j) ;
1624 
1625  //if(track_id != 0xFFFF) {
1626  // LOG(ERR,"Eh %d %d",track_id,*DTA(i,j)) ;
1627  //}
1628 
1629  blob_adc += *DTA(i,j) ;
1630 
1631  // did I see it before?
1632  int found = 0 ;
1633  int t_use = -1 ;
1634 
1635  for(int t=0;t<MAX_TRACK_IDS;t++) {
1636  if(t_id[t].track_id == track_id) {
1637  t_id[t].cou++ ;
1638  found = 1 ;
1639  break ;
1640  }
1641  else if(t_id[t].cou == 0) {
1642  t_use = t ;
1643  }
1644  }
1645 
1646  if(found) continue ;
1647  if(t_use < 0) continue ;
1648 
1649  // new track id...
1650  t_id[t_use].track_id = track_id ;
1651  t_id[t_use].cou++ ;
1652 
1653  }
1654  }}
1655 
1656  // get the track id with maximum counts
1657  int max_cou = 0 ;
1658  for(int t=0;t<MAX_TRACK_IDS;t++) {
1659  if(t_id[t].cou >= max_cou) {
1660  max_cou = t_id[t].cou ;
1661  track_id = t_id[t].track_id ;
1662  }
1663  }
1664 
1665  if(max_cou == 0) {
1666  LOG(ERR,"Wha?") ;
1667  }
1668 
1669  int track_adc = 0 ;
1670  for(int i=1;i<=dp;i++) {
1671  for(int j=1;j<=dt;j++) {
1672  if(*DTA_ID(i,j) == peaks[c].cluster_id) {
1673  if(*DTA_T(i,j) == track_id) {
1674  track_adc += *DTA(i,j) ;
1675  }
1676  }
1677  }
1678  }
1679 
1680 
1681  peaks[c].track_id = track_id ;
1682 
1683  if(blob_adc) {
1684  peaks[c].quality = (int)((double)track_adc/(double)blob_adc * 100.0 + 0.5) ;
1685  }
1686  else peaks[c].quality = 0 ;
1687 
1688 
1689  if(peaks[c].quality == 0) {
1690  LOG(ERR,"What??? %d %d %d",max_cou,track_adc,blob_adc) ;
1691  }
1692 
1693  }
1694 
1695 }
tpx_altro_struct
Definition: tpxCore.h:60
tpxFCF_2D
Definition: tpxFCF_2D.h:25
daq_sim_adc_tb
Definition: daq_dta_structs.h:26
Generated by   doxygen 1.8.5