StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
itpcFCF.cxx
1 #include <stdio.h>
2 #include <unistd.h>
3 #include <getopt.h>
4 #include <sys/types.h>
5 #include <stdlib.h>
6 #include <time.h>
7 #include <arpa/inet.h>
8 #include <sys/time.h>
9 #include <errno.h>
10 #include <assert.h>
11 
12 #include <rtsLog.h> // for my LOG() call
13 #include <rtsSystems.h>
14 
15 // this needs to be always included
16 //#include <DAQ_READER/daqReader.h>
17 #include <DAQ_READER/daq_dta.h>
18 
19 //#include <trgDataDefs.h>
20 //#include "trgConfNum.h"
21 
22 // only the detectors we will use need to be included
23 // for their structure definitions...
24 #include <DAQ_TPX/daq_tpx.h>
25 #include <DAQ_TPX/tpxFCF_flags.h>
26 #include <TPC/rowlen.h>
27 
28 #include <DAQ_ITPC/daq_itpc.h>
29 #include <DAQ_ITPC/itpcCore.h>
30 #include <DAQ_ITPC/itpcPed.h>
31 #include <DAQ_ITPC/itpcInterpreter.h>
32 #include <DAQ_ITPC/itpc_rowlen.h>
33 
34 #include <DAQ_ITPC/itpcFCF.h>
35 
36 #define VERSION 0x20180002
37 // 0x20180001 until Jun 6 -- had cuts in do_ch()
38 
39 
40 //#define DO_DBG1 1
41 
42 static double mark(void)
43 {
44  struct timeval tmval ;
45 
46  gettimeofday(&tmval,0) ;
47 
48  return ((double)tmval.tv_sec*1000000.0 + (double)tmval.tv_usec) ;
49 }
50 
51 static double delta(double v)
52 {
53  return mark() - v ;
54 }
55 
56 itpc_fcf_c::itpc_fcf_c()
57 {
58  want_data = 1 ; // ALWAYS!
59  my_id = 0 ; // normally starts from 1 so 0 is "un-initialized"
60 
61  version = VERSION ;
62  sector_id = 0 ; // 0 is ALL sectors
63  offline = 1 ; // new: offline is default
64 
65 
66  det_type = 1 ; // ITPC
67  y_is_timebin = 1 ; // normal
68  max_x = 120 ;
69  max_y = 512 ;
70  max_slice = 40 ;
71 
72  use_gain = 1 ;
73 
74  row_pad_store = 0 ;
75  words_per_cluster = 2 ;
76 
77  track_dta = 0 ;
78 
79  // just in case
80  run_start() ;
81 
82  f_stat.s1_found = 0 ;
83  blob_id = 0 ;
84 
85 } ;
86 
87 itpc_fcf_c::~itpc_fcf_c()
88 {
89 // LOG(TERR,"%s: destructor %d",__PRETTY_FUNCTION__,my_id) ;
90 
91  for(int i=0;i<=MAX_SEC;i++) {
92  if(sec_gains[i]) free(sec_gains[i]) ;
93  sec_gains[i] = 0 ;
94  }
95 
96  if(track_dta) {
97  free(track_dta) ;
98  track_dta = 0 ;
99  }
100 
101  if(row_pad_store) {
102  free(row_pad_store) ;
103  row_pad_store = 0 ;
104  }
105 } ;
106 
107 // generic "rowlen" function
108 int itpc_fcf_c::x_max(int slice, int y)
109 {
110  if(det_type == 0) { // TPX
111  if(y_is_timebin) { // normal
112  if(slice>45) return -1 ;
113  return tpc_rowlen[slice] ;
114  }
115  else { // y is row
116  if(y>45) return -1 ;
117  return 182 ; // tpc_rowlen[y] ;
118  }
119  }
120  else if(det_type == 1) {
121  if(y_is_timebin) {
122  if(slice>40) return -1 ;
123  return itpc_rowlen[slice] ;
124  }
125  else {
126  if(y>40) return -1 ;
127  return 120 ; //itpc_rowlen[slice] ;
128  }
129  }
130  else { // generic
131  return max_x ;
132  }
133 }
134 
135 int itpc_fcf_c::x_min(int slice, int y)
136 {
137  switch(det_type) {
138  case 0 :
139  case 1 :
140  if(y_is_timebin) return 1 ;
141  break ;
142  }
143 
144  return 1 ; // generic case
145 }
146 
147 int itpc_fcf_c::y_min()
148 {
149  return 0 ; // always!
150 }
151 
152 int itpc_fcf_c::y_max()
153 {
154  switch(det_type) {
155  case 0 :
156  if(y_is_timebin) return 512 ;
157  return 45 ;
158  case 1:
159  if(y_is_timebin) return 425 ;
160  return 40 ;
161  }
162 
163  return max_y ;
164 }
165 
166 
167 
168 // Go from raw channel (aka pad) to x position
169 // Needed for reverse clusterfinding because of the pizza shape
170 int itpc_fcf_c::pad_to_x(int pad, int row)
171 {
172  switch(det_type) {
173  case 0 :
174  if(y_is_timebin) return pad ;
175  return (182/2 - tpc_rowlen[row]/2 + pad) ;
176  case 1 :
177  if(y_is_timebin) return pad ;
178  return (120/2 - itpc_rowlen[row]/2 + pad) ;
179  }
180 
181  // all other detectors are assumed to be square
182  return pad ;
183 }
184 
185 
186 
187 //static member
188 //int itpc_fcf_c::rowlen[45+1] ;
189 
190 itpc_fcf_c::gain_rp_t *itpc_fcf_c::sec_gains[MAX_SEC+1] ;
191 
192 int itpc_fcf_c::get_bad(int sec1, int row1, int pad1)
193 {
194  gain_rp_t (*gain_p)[MAX_PHYS_PAD+1] ;
195 
196  if(sec_gains[sec1]==0) return 0 ; // good
197 
198  gain_p = (gain_rp_t (*)[MAX_PHYS_PAD+1]) sec_gains[sec1] ;
199 
200  if(gain_p[row1][pad1].gain==0.0) return 1 ;
201 
202  return 0 ;
203 }
204 
205 float itpc_fcf_c::get_gain(int sec1, int row1, int pad1)
206 {
207  if(row1==0 || pad1==0) return 0.0 ;
208 
209  gain_rp_t (*gain_p)[MAX_PHYS_PAD+1] ;
210 
211  if(sec_gains[sec1]==0) return 0.0 ; // bad
212 
213  gain_p = (gain_rp_t (*)[MAX_PHYS_PAD+1]) sec_gains[sec1] ;
214 
215  return gain_p[row1][pad1].gain ;
216 
217 }
218 
219 float itpc_fcf_c::get_t0(int sec1, int row1, int pad1)
220 {
221  if(row1==0 || pad1==0) return 0.0 ;
222 
223  gain_rp_t (*gain_p)[MAX_PHYS_PAD+1] ;
224 
225  if(sec_gains[sec1]==0) return 0.0 ; // bad
226 
227  gain_p = (gain_rp_t (*)[MAX_PHYS_PAD+1]) sec_gains[sec1] ;
228 
229  return gain_p[row1][pad1].t0 ;
230 
231 }
232 
233 void itpc_fcf_c::set_gain(int sec1, int row1, int pad1, float gain)
234 {
235 // if(row1==0 || pad1==0) return 0.0 ;
236 
237  gain_rp_t (*gain_p)[MAX_PHYS_PAD+1] ;
238 
239  if(sec_gains[sec1]==0) return ; // bad
240 
241  gain_p = (gain_rp_t (*)[MAX_PHYS_PAD+1]) sec_gains[sec1] ;
242 
243  gain_p[row1][pad1].gain = gain ;
244 }
245 
246 u_char itpc_fcf_c::get_flags(int sec1, int row1, int pad1)
247 {
248  if(row1==0 || pad1==0) return 255 ;
249 
250  gain_rp_t (*gain_p)[MAX_PHYS_PAD+1] ;
251 
252  if(sec_gains[sec1]==0) return 255 ; // bad
253 
254  gain_p = (gain_rp_t (*)[MAX_PHYS_PAD+1]) sec_gains[sec1] ;
255 
256  return gain_p[row1][pad1].flags ;
257 
258 }
259 
260 void itpc_fcf_c::set_flags(int sec1, int row1, int pad1, u_char flags)
261 {
262 
263  gain_rp_t (*gain_p)[MAX_PHYS_PAD+1] ;
264 
265  if(sec_gains[sec1]==0) return ; // bad
266 
267  gain_p = (gain_rp_t (*)[MAX_PHYS_PAD+1]) sec_gains[sec1] ;
268 
269  gain_p[row1][pad1].flags = flags ;
270 }
271 
272 void itpc_fcf_c::zap_fee(int sec1, int rdo1, int port1)
273 {
274  itpcInterpreter *itpc ;
275 
276  int padplane_id = itpc->itpc_fee_map[sec1-1][rdo1-1][port1-1] ;
277 
278  for(int c=0;c<64;c++) {
279  u_char flags ;
280  int row, pad ;
281 
282  itpc_ifee_to_rowpad(padplane_id,c,row,pad) ;
283 
284  int p1 = pad - 1 ;
285  int p2 = pad + 1 ;
286 
287  if(p1<1) p1 = 1 ;
288  if(p2<x_max(row,0)) p2 = x_max(row,0) ;
289 
290  set_gain(sec1,row,pad,0.0) ;
291 
292  flags = get_flags(sec1,row,pad) ;
293  set_flags(sec1,row,pad,flags|3) ;
294 
295  flags = get_flags(sec1,row,p1) ;
296  set_flags(sec1,row,p1,flags|2) ;
297 
298  flags = get_flags(sec1,row,p2) ;
299  set_flags(sec1,row,p2,flags|1) ;
300 
301  }
302 
303  return ;
304 }
305 
306 struct itpc_fcf_c::rp_t *itpc_fcf_c::get_row_pad(int row, int pad)
307 {
308  int max_pad_all = max_x + 1 ;
309 
310 // HACK
311  if(offline) s1_data_length = (1 + max_y) * 2 ; // need more for track_id
312 // if(offline) s1_data_length = 1 + max_y * 2 ; // need more for track_id
313  else s1_data_length = 1 + max_y ;
314 
315  if(row_pad_store==0) { // initialize on first use...
316  int max_row_all = max_slice + 1 ;
317  int chunks = max_row_all * max_pad_all ;
318 
319  row_pad_store = (u_short *) valloc(chunks *s1_data_length * sizeof(u_short)) ;
320  assert(row_pad_store) ;
321  memset(row_pad_store,0,chunks *s1_data_length * sizeof(u_short)) ; // superstition...
322  }
323 
324  u_short *p = row_pad_store + ((row*max_pad_all) + pad)*s1_data_length ;
325 
326  return (rp_t *)p ;
327 }
328 
329 
330 // This is only called if I'm in Offline and AFTER the usual init(sector,"")!!!
331 int itpc_fcf_c::init(daq_dta *gain)
332 {
333  if(gain==0) return -1 ;
334 
335  int bad_ch = 0 ;
336  int tot_ch = 0 ;
337 
338  while(gain->iterate()) {
339  int s = gain->sec ;
340  int row = gain->row ;
341 
342  if(sec_gains[s]==0) continue ; // not for me
343  if(row==0) continue ;
344 
345  gain_rp_t (*sr)[MAX_PHYS_PAD+1] = (gain_rp_t (*)[MAX_PHYS_PAD+1]) sec_gains[s] ;
346 
347  daq_det_gain *gp = (daq_det_gain *) gain->Void ;
348 
349  for(u_int p=0;p<gain->ncontent;p++) {
350  if(p==0) continue ;
351 // if(p>(u_int)rowlen[row]) continue ;
352  if(p>(u_int)x_max(row,0)) continue ;
353 
354  tot_ch++ ;
355 
356  sr[row][p].gain = gp[p].gain ;
357  sr[row][p].t0 = gp[p].t0 ;
358 
359 
360  if(gp[p].gain < 0.01) {
361  int p1 = p - 1;
362  int p2 = p + 1 ;
363 
364  if(p1<1) p1 = 1;
365  //if(p2>rowlen[row]) p2 = rowlen[row] ;
366  if(p2>x_max(row,0)) p2 = x_max(row,0) ;
367 
368  // mark this and neighbours as bad
369  sr[row][p].flags |= 3 ; // bad and edge
370  sr[row][p1].flags |= 2 ; // just edge
371  sr[row][p2].flags |= 2 ; // just edge
372 
373 
374 
375  bad_ch++ ;
376  }
377  }
378 
379 
380  }
381 
382  if(bad_ch) {
383  LOG(WARN,"%d/%d bad channels",bad_ch,tot_ch) ;
384  }
385 
386  return 0 ;
387 }
388 
389 // det_type and various other max_x,max_y,max_slice etc HAVE to be set before!!!!
390 int itpc_fcf_c::init(int sec, const char *fname)
391 {
392 
393  int s_min, s_max ;
394  gain_rp_t (*gain_p)[MAX_PHYS_PAD+1] ;
395 
396  if(sec>MAX_SEC) return -1 ;
397 
398  if(sec<=0) {
399  s_min = 1 ;
400  s_max = MAX_SEC ;
401  }
402  else {
403  s_min = s_max = sec ;
404  }
405 
406  const char *s_type, *s_orient ;
407 
408 
409  switch(det_type) {
410  case 0 :
411  s_type = "TPX" ;
412  if(y_is_timebin) {
413  use_gain = 1 ;
414 
415  max_x = 182 ;
416  max_y = 512 ;
417  max_slice = 45 ;
418 
419  }
420  else {
421  use_gain = 0 ;
422 
423  max_x = 182 ;
424  max_y = 45 ;
425  max_slice = 512;
426  }
427  break ;
428  case 1 :
429  s_type = "ITPC" ;
430  if(y_is_timebin) {
431  use_gain = 1 ;
432 
433  max_x = 120 ;
434  max_y = 512 ;
435  max_slice = 40 ;
436  }
437  else {
438  use_gain = 0 ;
439 
440  max_x = 120 ;
441  max_y = 40 ;
442  max_slice = 512 ;
443  }
444  break ;
445  default :
446  s_type = "OTHER" ;
447  use_gain = 0 ;
448  break ;
449  }
450 
451  if(y_is_timebin) s_orient="across-rows" ;
452  else s_orient = "across-timebins" ;
453 
454  LOG(INFO,"%s: det_type %s(%s): max_x(pad) %d, max_y(timebin) %d, max_slice(row) %d",__FUNCTION__,
455  s_type,s_orient,
456  max_x, max_y, max_slice) ;
457 
458  //alocation first
459  for(int s=s_min;s<=s_max;s++) {
460 
461  if(sec_gains[s]) continue ; // silently continue ;
462 
463  sec_gains[s] = (gain_rp_t *) malloc((MAX_PHYS_ROW+1)*(MAX_PHYS_PAD+1)*sizeof(gain_rp_t)) ;
464 
465  //defaults now
466  gain_p = (gain_rp_t (*)[MAX_PHYS_PAD+1]) sec_gains[s] ;
467 
468  for(int r=1;r<=MAX_PHYS_ROW;r++) {
469 
470  for(int p=1;p<=MAX_PHYS_PAD;p++) {
471 
472  gain_p[r][p].gain = 1.0 ;
473  gain_p[r][p].t0 = 0.0 ;
474 
475  //if(p==1 || p==rowlen[r]) gain_p[r][p].flags = 0x2 ; // edge
476 
477  if(p==1 || p==x_max(r,0)) gain_p[r][p].flags = 0x2 ; // edge
478  else gain_p[r][p].flags = 0 ;
479 
480  }
481  }
482 
483 
484  }
485 
486 
487 
488  // do we have a file?
489  if(fname==0) return 0 ;
490 
491 
492 
493  // yes, open it
494  FILE *f = fopen(fname,"r") ;
495  if(f==0) {
496  LOG(ERR,"%s: %s [%s]",__FUNCTION__,fname,strerror(errno)) ;
497  return -1 ;
498  }
499 
500  LOG(INFO,"Using gain file %s",fname) ;
501 
502  int ch_bad = 0 ;
503  int ch_all = 0 ;
504  int fee_bad = 0 ;
505 
506  // we now have an opened gain file
507  while(!feof(f)) {
508  int sec,rdo,port,ch,row,pad ;
509  float g, t ;
510 
511  char buff[128] ;
512 
513  if(fgets(buff,sizeof(buff),f)==0) continue ;
514 
515  if(buff[0]=='#') continue ;
516  if(strlen(buff)<1) continue ;
517 
518  int ret = sscanf(buff,"%d %d %d %d %d %d %f %f",&sec,&rdo,&port,&ch,&row,&pad,&g,&t) ;
519 
520  if(ret != 8) continue ;
521 
522  if(sec_gains[sec]==0) continue ; // not for me!
523 
524  if(ch<0) { // kill entire FEE!
525  zap_fee(sec,rdo,port) ;
526  fee_bad++ ;
527  LOG(TERR,"%d %d %d kill FEE implemented",sec,rdo,port) ;
528  continue ;
529  }
530 
531 
532  gain_rp_t (*sr)[MAX_PHYS_PAD+1] = (gain_rp_t (*)[MAX_PHYS_PAD+1]) sec_gains[sec] ;
533 
534  sr[row][pad].gain = g ;
535  sr[row][pad].t0 = t ;
536 
537  ch_all++ ;
538 
539  if(g<0.01) { // some really small number close to 0.0
540  int p1 = pad - 1;
541  int p2 = pad + 1 ;
542 
543  if(p1<1) p1 = 1;
544  //if(p2>rowlen[row]) p2 = rowlen[row] ;
545  if(p2>x_max(row,0)) p2 = x_max(row,0) ;
546 
547  // mark this and neighbours as bad
548  sr[row][pad].flags |= 3 ; // bad and edge
549  sr[row][p1].flags |= 2 ; // just edge
550  sr[row][p2].flags |= 2 ; // just edge
551 
552  ch_bad++ ;
553  }
554 
555  }
556 
557  fclose(f) ;
558 
559  if(ch_bad) LOG(WARN,"...with %d/%d bad channels, %d bad_fees",ch_bad,ch_all,fee_bad) ;
560  else LOG(TERR,"...with %d/%d no bad channels, %d bad_fees",ch_bad,ch_all,fee_bad) ;
561 
562  return 0 ;
563 }
564 
565 
566 
567 
568 
569 
570 // At start event
571 void itpc_fcf_c::event_start()
572 {
573  //for statistics, so not really necessary
574  f_stat.s1_found = 0 ;
575 
576  if(offline) {
577  for(int r=0;r<=max_slice;r++) {
578  for(int p=0;p<=max_x;p++) {
579  get_row_pad(r,p)->s1_len = 0 ;
580  }
581  }
582  }
583 
584 }
585 
586 
587 int itpc_fcf_c::fcf_decode(u_int *p_buff, daq_sim_cld_x *dc, u_int version)
588 {
589  daq_cld cld ;
590  int words_used ;
591 
592  words_used = fcf_decode(p_buff,&cld,version) ;
593 
594  memcpy(&dc->cld,&cld,sizeof(cld)) ;
595 
596  p_buff += words_used ; // point to track id
597 
598  dc->track_id = *p_buff & 0xFFFF ;
599  dc->quality = (*p_buff)>>16 ;
600 
601  p_buff++ ; // point to adc
602 
603  dc->max_adc = *p_buff & 0xFFFF ;
604  dc->pixels = (*p_buff)>>16 ;
605 
606  return words_used + 2 ;
607 
608 }
609 
610 int itpc_fcf_c::fcf_decode(u_int *p_buff, daq_cld *dc, u_int version)
611 {
612  double p, t ;
613  int p1,p2,t1,t2,cha,fla ;
614  u_int p_tmp, t_tmp ;
615 
616  // pad
617  p_tmp = *p_buff & 0xFFFF ;
618 
619  // time
620  t_tmp = *p_buff >> 16 ;
621 
622 
623  p = (double)(p_tmp & 0x3FFF) / 64.0 ;
624  t = (double)(t_tmp & 0x7FFF) / 64.0 ;
625 
626  fla = 0 ;
627  if(p_tmp & 0x8000) fla |= FCF_MERGED ;
628  if(p_tmp & 0x4000) fla |= FCF_DEAD_EDGE ;
629  if(t_tmp & 0x8000) fla |= FCF_ONEPAD ;
630 
631 
632  p_buff++ ; // advance to next word
633  cha = *p_buff >> 16 ;
634 
635  if(cha >= 0x8000) { // special case of very large charge...
636 //printf("1: Big cha: 0x%08X %d; %f %f\n",cha,cha,p,t) ;
637 
638  fla |= FCF_BIG_CHARGE ;
639  cha = (cha & 0x7FFF) * 1024 ;
640 // if(cha == 0) cha = 0x8000; // exaclty, but can't be I think...
641 
642 //printf("2: Big cha: 0x%08X %d\n",cha,cha) ;
643 
644  // quasi fix of the very large problem...
645  if(cha > 0xFFFF) cha = 0xFFFF ; // because the daq_cld structure has charge as a short... damn...
646 
647  }
648 
649  p_tmp = *p_buff & 0xFFFF ; // reuse p_tmp
650 
651  if(p_tmp & 0x8000) fla |= FCF_ROW_EDGE ;
652  if(p_tmp & 0x4000) fla |= FCF_BROKEN_EDGE ;
653 
654  t1 = p_tmp & 0xF ;
655  t2 = (p_tmp >> 4) & 0xF ;
656 
657 
658  p1 = (p_tmp >> 8) & 0x7 ;
659  p2 = (p_tmp >> 11) & 0x7 ;
660 
661 
662  t1 = (int)t - t1 ;
663  t2 = (int)t + t2 ;
664 
665 
666  p1 = (int)p - p1 ;
667  p2 = (int)p + p2 ;
668 
669  dc->t1 = t1 ;
670  dc->t2 = t2 ;
671  dc->p1 = p1 ;
672  dc->p2 = p2 ;
673  dc->charge = cha ; // this is a problem for BIG_CHARGE... it will strip the upper bits... unsolved.
674  dc->flags = fla ;
675  dc->pad = p ;
676  dc->tb = t ;
677 
678 
679  return 2 ; // 2 u_ints used
680 
681 }
682 
683 
684 // Called by the raw data unpacker: fee_id and fee_ch are _physical_ channels
685 // Main purpose is unpack the raw data into the canonical form ready for FCF.
686 // Canonical form is:
687 // seq_ix (set to 0 before cluster finder begins)
688 // timebin_count
689 // timebin_start
690 // adc...
691 // adc...
692 
693 int itpc_fcf_c::do_ch(int fee_id, int fee_ch, u_int *data, int words)
694 {
695  int row, pad ;
696  int seq_cou ;
697  int s_count ;
698  int t_stop_last ;
699  u_short *s1_data ;
700  int t_cou ;
701  int err = 0 ;
702 
703  seq_cou = 0;
704  int word32 = (words/3) + (words%3?1:0) ;
705 
706  itpc_ifee_to_rowpad(fee_id, fee_ch, row, pad) ;
707 
708  gain_rp_t (*gain_row_pad)[MAX_PHYS_PAD+1] = (gain_rp_t (*)[MAX_PHYS_PAD+1]) sec_gains[sector_id] ;
709 
710  if(row==0) goto err_ret ; // unphysical pad
711 
712 
713 // LOG(TERR,"RP %d:%d = %d %d",row,pad,gain_row_pad[row][pad].flags,word32) ;
714 
715  if(gain_row_pad[row][pad].flags & 1) {
716  //LOG(WARN,"RP %d:%d killed",row,pad) ;
717  goto err_ret ; // skip dead pads!
718  }
719 
720 
721 
722  if((word32*3)>=MAX_TB_EVER) {
723  LOG(ERR,"%d:#%d: FEE %d:%d, words %d",rdo,port,fee_id,fee_ch,words) ;
724  err |= 1 ;
725  goto err_ret ;
726  }
727 
728  t_cou = 0 ;
729  for(int i=0;i<word32;i++) {
730  u_int d = *data++ ;
731 
732  if(d & 0xC0000000) {
733  seq_cou = 0 ;
734  LOG(ERR,"%d:#%d: FEE %d:%d, words %d",rdo,port,fee_id,fee_ch,words) ;
735  err |= 2 ;
736  goto err_ret ;
737  }
738 
739  tb_buff[t_cou++] = (d>>20) & 0x3FF ;
740  tb_buff[t_cou++] = (d>>10) & 0x3FF ;
741  tb_buff[t_cou++] = d & 0x3FF ;
742  }
743 
744 
745 
746 // s1_data = row_pad[row][pad].s1_data ;
747  s1_data = get_row_pad(row,pad)->s1_data ;
748 
749  // I MUST put protection against broken data!!!
750  t_stop_last = -1 ;
751 
752  for(int i=0;i<words;) { // now timebins!
753  int t_cou = tb_buff[i++] ;
754  int t_start = tb_buff[i++] ;
755  int t_stop = t_start + t_cou - 1 ;
756 
757  // happens too often, get rid of ERR
758  if(t_start <= t_stop_last) {
759  LOG(NOTE,"%d:#%d: FEE %d:%d, words %d: %d <= %d",rdo,port,fee_id,fee_ch,words,t_start,t_stop_last) ;
760  seq_cou = 0 ;
761  err |= 4 ;
762  goto err_ret ;
763  }
764  if(t_stop > 511) {
765  LOG(NOTE,"%d:#%d: FEE %d:%d, words %d: %d > 511",rdo,port,fee_id,fee_ch,words,t_stop) ;
766  seq_cou = 0 ;
767  err |= 4 ;
768  goto err_ret ;
769  }
770 
771  t_stop_last = t_stop ;
772 
773  seq_cou++ ;
774 
775  *s1_data++ = 0 ; // index now 0!
776  *s1_data++ = t_cou ;
777  *s1_data++ = t_start ;
778 
779 
780 
781  for(int t=t_start;t<=t_stop;t++) {
782  // initial cuts, where I blow of data
783 // Test: was 0
784 #if 0
785  // cut timebin due to gating grid pickup
786  if(t>425) {
787  *s1_data = 0 ;
788  }
789  else if((t>=26)&&(t<=31)) *s1_data = 0 ;
790  else {
791  if(tb_buff[i]<=4) *s1_data = 0 ; // cut low ADC data
792  else *s1_data = tb_buff[i] ;
793  }
794 
795 
796  i++ ;
797  s1_data++ ;
798 #else
799  *s1_data++ = tb_buff[i++] ;
800 #endif
801 
802  }
803  }
804 
805  *s1_data++ = 0xFFFF ; // end sentinel
806 
807  // check for data overrun!
808 // s_count = s1_data - row_pad[row][pad].s1_data ;
809  s_count = s1_data - get_row_pad(row,pad)->s1_data ;
810  if(s_count >= s1_data_length) {
811  err |= 8 ;
812  LOG(ERR,"In trouble at RP %d:%d",row,pad) ;
813  }
814 
815  // for later optimization!
816  if(s_count >= (int)f_stat.max_s1_len) {
817  f_stat.max_s1_len = s_count ;
818  }
819 
820  f_stat.s1_found += seq_cou ;
821 
822 // if(seq_cou) LOG(TERR,"RP %d:%d = seq_cou %d",row,pad,seq_cou) ;
823 
824  err_ret:;
825 
826 // row_pad[row][pad].s1_len = seq_cou ;
827  get_row_pad(row,pad)->s1_len = seq_cou ; // sequence count!
828 
829  return err ;
830 }
831 
832 // Called by the raw data unpacker: fee_id and fee_ch are _physical_ channels
833 // Main purpose is unpack the raw data into the canonical form ready for FCF.
834 // Canonical form is:
835 // seq_ix (set to 0 before cluster finder begins)
836 // timebin_count
837 // timebin_start
838 // adc...
839 // adc...
840 
841 int itpc_fcf_c::do_ch_sim(int row, int pad, u_short *tb_buff, u_short *track_id)
842 {
843  int seq_cou ;
844  int s_count ;
845 
846  u_short *s1_data ;
847 
848  int t_start ;
849  int t_cou ;
850  u_short *p_start = 0 ;
851 
852  offline = 1 ; // juuuuust in case!
853  words_per_cluster = 4 ; // added stuff
854  if(track_dta==0) {
855  track_dta = (u_short *)malloc(MAX_BLOB_SIZE*2) ;
856  }
857 
858  seq_cou = 0;
859 
860 
861 
862  if(row==0) goto err_ret ; // unphysical pad
863 
864  if(use_gain) {
865  gain_rp_t (*gain_row_pad)[MAX_PHYS_PAD+1] = (gain_rp_t (*)[MAX_PHYS_PAD+1]) sec_gains[sector_id] ;
866  if(gain_row_pad[row][pad].flags & 1) {
867  goto err_ret ; // skip dead pads!
868  }
869  }
870 
871 // s1_data = row_pad[row][pad].s1_data ;
872  s1_data = get_row_pad(row,pad)->s1_data ;
873 
874  t_start = -1 ;
875  t_cou = 0 ;
876 
877  //condition data ala cuts
878  if(det_type==1 && y_is_timebin) { // only for ITPC in normal orientation!
879 
880 // also removed the cuts in the Offline version on 26-Jan-2019
881 #if 0
882  for(int i=0;i<512;i++) {
883  if(i>425) tb_buff[i] = 0 ;
884  else if((i>=26)&&(i<=31)) tb_buff[i]= 0 ;
885  else {
886  if(tb_buff[i]<=4) tb_buff[i] = 0 ;
887  }
888  }
889 #endif
890  }
891 
892  for(int i=0;i<MAX_TB_EVER;i++) {
893  if(tb_buff[i]) {
894  if(t_start<0) {
895  *s1_data++ = 0 ;
896  p_start = s1_data ;
897  *s1_data++ = 0 ; //t_cou
898  *s1_data++ = i ; //t_start ;
899 
900  t_start = i ;
901  seq_cou++ ;
902  }
903  //NEW in 11Jun2023
904  if(tb_buff[i]==0xFFFF) { // this is how I skip zeros
905  *s1_data++ = 0 ;
906  }
907  else {
908  *s1_data++ = tb_buff[i] ; // ADC data
909  }
910  *s1_data++ = track_id[i] ;
911  t_cou++ ;
912  }
913  else {
914  if(t_start >= 0) { //was started
915  *p_start = t_cou ;
916  t_start = -1 ;
917  t_cou = 0 ;
918  }
919  }
920  }
921 
922  if(t_start >= 0) {
923  *p_start = t_cou ;
924  t_start = -1 ;
925  }
926 
927 
928  *s1_data++ = 0xFFFF ; // end sentinel
929 
930  // check for data overrun!
931 // s_count = s1_data - row_pad[row][pad].s1_data ;
932  s_count = s1_data - get_row_pad(row,pad)->s1_data ;
933  if(s_count >= s1_data_length) {
934  LOG(ERR,"Data too big at RP %d:%d",row,pad) ;
935  }
936 
937 // LOG(TERR,"RP %d %d = %d",row,pad,s_count) ;
938 
939  // for later optimization!
940  if(s_count >= (int)f_stat.max_s1_len) {
941  f_stat.max_s1_len = s_count ;
942  }
943 
944  f_stat.s1_found += seq_cou ;
945 
946 // if(seq_cou) LOG(TERR,"RP %d:%d = seq_cou %d",row,pad,seq_cou) ;
947 
948  err_ret:;
949 
950 // row_pad[row][pad].s1_len = seq_cou ; // sequence count!
951  get_row_pad(row,pad)->s1_len = seq_cou ;
952 
953  return 0 ;
954 }
955 
956 
957 // Returns words(ints) of storage.
958 int itpc_fcf_c::do_fcf(void *v_store, int bytes)
959 {
960  out_store = (u_int *) v_store ;
961  max_out_bytes = bytes ;
962 
963  u_int *store_start = out_store ;
964 
965  f_stat.evt_cou++ ;
966 
967  blob_id = 0 ;
968 
969  for(int row=1;row<=max_slice;row++) {
970  double tmx ;
971  int found_ints ;
972 
973  u_int *row_store = out_store++ ; // leave space for row number
974  out_store++ ; // leave space for version
975  out_store++ ; // leave space for cluster count
976 
977 
978  tmx = mark() ;
979  do_blobs_stage1(row) ;
980  tmx = delta(tmx) ;
981 
982  f_stat.tm[1] += tmx ;
983 
984  //do_row_check(row) ;
985 
986  tmx = mark() ;
987  do_blobs_stage2(row) ;
988  tmx = delta(tmx) ;
989 
990  f_stat.tm[2] += tmx ;
991 
992  tmx = mark() ;
993  found_ints = do_blobs_stage3(row) ;
994  tmx = delta(tmx) ;
995 
996  f_stat.tm[3] += tmx ;
997 
998 
999  if(found_ints) {
1000 
1001  *row_store++ = (words_per_cluster<<16)|row ; // words-per-cluster | row
1002  *row_store++ = version ;
1003  *row_store++ = found_ints ; // in ints!
1004 
1005  out_store += found_ints ;
1006  }
1007  else {
1008  out_store = row_store ; // rewind
1009  }
1010 
1011  int ints_used = out_store - store_start ;
1012  //int ints_needed = (max_out_bytes/4)-1000 ;
1013  int ints_needed = (max_out_bytes/4)-400 ;
1014 
1015  if(ints_used && (ints_used>ints_needed)) {
1016  LOG(ERR,"not enough ints %d vs %d",out_store-store_start,max_out_bytes/4) ;
1017  break ;
1018  }
1019  }
1020 
1021 
1022  if(f_stat.s1_found > f_stat.max_s1_found) {
1023  f_stat.max_s1_found = f_stat.s1_found ;
1024  }
1025 
1026  return (out_store-store_start) ; // in ints
1027 }
1028 
1029 void itpc_fcf_c::run_start()
1030 {
1031  // all just for statistics and monitoring...
1032  LOG(NOTE,"%s: %d",__PRETTY_FUNCTION__,my_id) ;
1033 
1034  memset(&f_stat,0,sizeof(f_stat)) ;
1035 
1036  for(int r=0;r<=max_slice;r++) {
1037  for(int p=0;p<=max_x;p++) {
1038  get_row_pad(r,p)->s1_len = 0 ;
1039  }
1040  }
1041 
1042 
1043 }
1044 
1045 void itpc_fcf_c::run_stop()
1046 {
1047  // all just for statistics and moniroing
1048 
1049 // LOG(TERR,"%s: %d",__PRETTY_FUNCTION__,my_id) ;
1050 
1051  for(int i=0;i<10;i++) {
1052  f_stat.tm[i] /= f_stat.evt_cou ;
1053  }
1054 
1055  LOG(INFO,"itpcFCF: %d: events %d, times %f %f %f %f",my_id,
1056  f_stat.evt_cou,
1057  f_stat.tm[0],f_stat.tm[1],f_stat.tm[2],f_stat.tm[3]) ;
1058 
1059  LOG(INFO," times %f %f %f %f",
1060  f_stat.tm[4],f_stat.tm[5],f_stat.tm[6],f_stat.tm[7]) ;
1061 
1062 
1063  LOG(INFO," max s1_found %d, s1_len %d, blob_cou %d",f_stat.max_s1_found,f_stat.max_s1_len,f_stat.max_blob_cou) ;
1064  if(f_stat.toobigs) {
1065  LOG(WARN," toobigs %d",f_stat.toobigs) ;
1066  }
1067 
1068 }
1069 
1070 
1071 int itpc_fcf_c::do_blobs_stage1(int row)
1072 {
1073  int pads = 0 ;
1074  int late_merges = 0 ;
1075 
1076 // rp_t *rp = &(row_pad[row][0]) ;
1077 // rp_t *rp = get_row_pad(row,0) ;
1078 
1079 // int rl = rowlen[row] ;
1080  int rl = x_max(row,0) ;
1081 
1082  blob_cou = 1 ;
1083 
1084 // LOG(TERR,"stage1: row %d, x_max %d",row,rl) ;
1085 
1086  for(int p=1;p<=rl;p++) { // < is on purpose!!!
1087  int t1_cou, t1_lo, t1_hi ;
1088  int t2_cou, t2_lo, t2_hi ;
1089  u_short *ix1_p ;
1090  u_short *ix2_p ;
1091 
1092  rp_t *rp1 = get_row_pad(row,p) ;
1093 
1094  if(rp1->s1_len==0) continue ;
1095 
1096  //LOG(TERR,"Row %d: s1_len %d",row,rp[p].s1_len) ;
1097 
1098  u_short *d1 = rp1->s1_data ;
1099 
1100  for(int i=0;i<rp1->s1_len;i++) {
1101  u_short *d2 ;
1102 
1103  rp_t *rp2 = get_row_pad(row,p+1) ;
1104 
1105  ix1_p = d1++ ;
1106 
1107 
1108 
1109  t1_cou = *d1++ ; // actually tb_cou
1110  t1_lo = *d1++ ; // actually t_start ;
1111 
1112  t1_hi = t1_lo + t1_cou - 1 ; // now is really
1113 
1114  if(offline) {
1115  d1 += t1_cou * 2 ; // also track_id!
1116  }
1117  else {
1118  d1 += t1_cou ; // advance to next
1119  }
1120 
1121  if(p==rl) goto sweep_unmerged ; // since there is no pad after that
1122 
1123  d2 = rp2->s1_data ;
1124 
1125  for(int j=0;j<rp2->s1_len;j++) {
1126  ix2_p = d2++ ;
1127 
1128  t2_cou = *d2++ ;
1129  t2_lo = *d2++ ;
1130 
1131  t2_hi = t2_lo + t2_cou - 1 ;
1132 
1133  if(offline) d2 += t2_cou * 2 ;
1134  else d2 += t2_cou ;
1135 
1136  int merge = 0 ;
1137 
1138  if(t1_lo > t2_hi) merge = 0 ;
1139  else if(t2_lo > t1_hi) merge = 0 ;
1140  else merge = 1 ;
1141 
1142  if(merge==0) continue ;
1143 
1144 
1145  if(*ix1_p > 0) {
1146  if(*ix2_p > 0) {
1147  if(*ix1_p != *ix2_p) {
1148  // I have to merge 2 distinct blobs
1149  // I will have to decide later what to do
1150  //LOG(TERR,"Late merge %d %d",*ix1_p,*ix2_p) ;
1151  late_merges++ ;
1152 
1153  if(blob_ix[*ix1_p] < blob_ix[*ix2_p]) {
1154  blob_ix[*ix2_p] = blob_ix[*ix1_p] ;
1155 
1156  blob[blob_ix[*ix2_p]].merges++ ;
1157  }
1158  else {
1159  blob_ix[*ix1_p] = blob_ix[*ix2_p] ;
1160 
1161  blob[blob_ix[*ix1_p]].merges++ ;
1162  }
1163 
1164  }
1165  }
1166  else {
1167  *ix2_p = *ix1_p ;
1168  }
1169  }
1170  else {
1171  if(*ix2_p > 0) {
1172  *ix1_p = *ix2_p ;
1173  }
1174  else {
1175  // new blob!
1176  *ix1_p = blob_cou ;
1177  *ix2_p = blob_cou ;
1178  blob_ix[blob_cou] = blob_cou ;
1179  blob[blob_cou].merges = 0 ;
1180  blob_cou++ ;
1181  }
1182  }
1183 
1184  }
1185 
1186  sweep_unmerged: ;
1187 
1188  if(*ix1_p == 0) { // still unassigned!
1189  // new blob
1190  *ix1_p = blob_cou ;
1191  blob_ix[blob_cou] = blob_cou ;
1192  blob[blob_cou].merges = 0 ;
1193  blob_cou++ ;
1194  }
1195  }
1196 
1197 
1198  pads++ ;
1199  }
1200 
1201  if(blob_cou > f_stat.max_blob_cou) {
1202  f_stat.max_blob_cou = blob_cou ;
1203  }
1204 
1205 
1206 // LOG(TERR,"Row %d: got %d pads, %d blobs, %d late merges",row,pads,blob_cou,late_merges) ;
1207 
1208 
1209  return blob_cou ;
1210 }
1211 
1212 int itpc_fcf_c::do_blobs_stage2(int row)
1213 {
1214 // rp_t *rp = &(row_pad[row][0]) ;
1215 // rp_t *rp = get_row_pad(row,0) ;
1216 
1217 // int rl = rowlen[row] ;
1218  int rl = x_max(row,0) ;
1219 
1220  //LOG(TERR,"stage2: row %d",row) ;
1221 
1222  for(int i=0;i<blob_cou;i++) {
1223  blob[i].seq_cou = 0 ; // mark as unused
1224 
1225  // initialize extents
1226  blob[i].p1 = 0xFFFF ;
1227  blob[i].p2 = 0 ;
1228 
1229  blob[i].t1 = 0xFFFF ;
1230  blob[i].t2 = 0 ;
1231 
1232  blob[i].flags = 0 ;
1233  }
1234 
1235  for(int p=1;p<=rl;p++) {
1236  rp_t *rp = get_row_pad(row,p) ;
1237 
1238  if(rp->s1_len==0) {
1239  //LOG(WARN,"row:pad %d:%d is void",row,p) ;
1240  continue ;
1241  }
1242 
1243  u_short *d = rp->s1_data ;
1244 
1245  for(int i=0;i<rp->s1_len;i++) { // loop over sequnces
1246  int t_cou, t_start, t_stop ;
1247 
1248  int ix = *d++ ;
1249 
1250 
1251  int b_ix = blob_ix[ix] ;
1252 
1253  //LOG(TERR," using bix %d %d",ix,b_ix) ;
1254  if(b_ix != ix) {
1255  LOG(WARN,"Can't be: %d %d, RP %d:%d",b_ix,ix,row,p) ;
1256  }
1257 
1258  if(b_ix==0) {
1259  LOG(ERR,"Can't be: %d %d, RP %d:%d",b_ix,ix,row,p) ;
1260  }
1261 
1262  blob_t *bl = &(blob[b_ix]) ;
1263 
1264 
1265  t_cou = *d++ ;
1266  t_start = *d++ ;
1267  t_stop = t_start + t_cou - 1 ;
1268 
1269  if(offline) d += 2*t_cou ;
1270  else d += t_cou ;
1271 
1272 
1273  //LOG(TERR,"%d %d: %d %d",t_start,t_cou,bl->t1,bl->t2) ;
1274 
1275  if(p > bl->p2) bl->p2 = p ;
1276  if(p < bl->p1) bl->p1 = p ;
1277 
1278  if(t_stop > bl->t2) bl->t2 = t_stop ;
1279  if(t_start < bl->t1) bl->t1 = t_start ;
1280 
1281 
1282  bl->seq_cou++ ;
1283  }
1284  }
1285 
1286  //cleanup: cuts etc
1287  int blob_ok = 0 ;
1288  for(int i=0;i<blob_cou;i++) {
1289  if(blob[i].seq_cou == 0) continue ;
1290 
1291  int dp = blob[i].p2 - blob[i].p1 + 1 ;
1292  int dt = blob[i].t2 - blob[i].t1 + 1;
1293 
1294  if(dp<=0) {
1295  blob[i].seq_cou = 0 ;
1296  LOG(ERR,"%d: dp %d",rdo,dp) ;
1297  continue ;
1298  }
1299 
1300  if(dt<=0) {
1301  blob[i].seq_cou = 0 ;
1302  LOG(ERR,"%d: dt %d",rdo,dt) ;
1303  continue ;
1304  }
1305 
1306 
1307 // if(blob[i].flags & FCF_FLAGS_TRUNC) {
1308 // //LOG(WARN,"truncated") ;
1309 // blob[i].seq_cou = 0 ; // kill it
1310 // continue ;
1311 // }
1312 
1313  if(dp<=1) { // one pad
1314  //LOG(WARN,"%d: 1pad %d %d: %d",i,dp,dt,blob[i].seq_cou) ;
1315  //LOG(WARN,"%d %d %d %d",blob[i].p1,blob[i].p2,blob[i].t1,blob[i].t2) ;
1316  blob[i].seq_cou = 0 ; // kill it
1317  //blob[i].p1 *= 1000 ;
1318  continue ;
1319  }
1320 
1321 
1322 
1323  if(dt<=1 || (dt<=3 && y_is_timebin)) { // tb range < 3
1324  //LOG(WARN,"%d: 3tb %d %d %d",i,dp,dt,blob[i].seq_cou) ;
1325  //LOG(WARN,"%d %d %d %d",blob[i].p1,blob[i].p2,blob[i].t1,blob[i].t2) ;
1326  blob[i].seq_cou = 0 ; // kill it
1327  blob[i].p1 *= 2000 ;
1328  continue ;
1329  }
1330 
1331  u_int bytes_for_blob = (dp+2)*(dt+2)*2*2 ;
1332 
1333  if(bytes_for_blob > sizeof(smooth_dta)) {
1334 // if((u_int)(dt*dp) > sizeof(smooth_dta)/sizeof(smooth_dta[0])) { // too big!
1335  f_stat.toobigs++ ;
1336 
1337 // LOG(ERR,"row %d: %d: toobig %d X %d",row,i,dp,dt) ;
1338  //LOG(WARN,"%d %d %d %d",blob[i].p1,blob[i].p2,blob[i].t1,blob[i].t2) ;
1339  blob[i].seq_cou = 0 ;
1340  blob[i].p1 *= 3000 ;
1341  continue ;
1342  }
1343 
1344  blob_ok++ ;
1345  }
1346 
1347 
1348 
1349 
1350 #if 0
1351  LOG(TERR,"Blobs OK %d/%d in row %d",blob_ok,blob_cou,row) ;
1352 
1353  for(int i=0;i<blob_cou;i++) {
1354  //if(blob[i].seq_cou==0) continue ;
1355 
1356  LOG(TERR,"Blob %d: seq %d: flags 0x%X: pad %d:%d, tb %d:%d",i,blob[i].seq_cou,blob[i].flags,
1357  blob[i].p1,blob[i].p2,blob[i].t1,blob[i].t2) ;
1358  }
1359 
1360 #endif
1361 
1362  return blob_ok ;
1363 }
1364 
1365 // Returns the number of bytes
1366 int itpc_fcf_c::do_blobs_stage3(int row)
1367 {
1368  short *dta_s ;
1369 // rp_t *rp = &(row_pad[row][0]) ;
1370 // rp_t *rp = get_row_pad(row,0) ;
1371  double tm ;
1372 
1373  u_int *obuff = (u_int *)out_store ;
1374 
1375  //LOG(TERR,"stage3: row %d",row) ;
1376 
1377 
1378  gain_rp_t (*gain_row_pad)[MAX_PHYS_PAD+1] = (gain_rp_t (*)[MAX_PHYS_PAD+1]) sec_gains[sector_id] ;
1379 
1380 
1381  int clusters_cou = 0 ;
1382 
1383 // LOG(TERR,"stage3: row %d, blob_cou %d",row,blob_cou) ;
1384 
1385  for(int ix=0;ix<blob_cou;ix++) {
1386  if(blob[ix].seq_cou==0) continue ;
1387 
1388 
1389 
1390  tm = mark() ;
1391 
1392 
1393 
1394  blob[ix].tot_charge = 0 ;
1395  blob[ix].pixels = 0 ;
1396 
1397  // calculate the necessary size
1398  int dt = blob[ix].t2 - blob[ix].t1 + 1 ;
1399  int dp = blob[ix].p2 - blob[ix].p1 + 1 ;
1400 
1401 
1402  int dt_2 = dt + 2 ;
1403 
1404 
1405 #ifdef DO_DBG1
1406  printf("...bytes %lu vs %lu\n",dt_2*(dp+2)*sizeof(short)*2,sizeof(smooth_dta)) ;
1407 #endif
1408  memset(smooth_dta,0,dt_2*(dp+2)*sizeof(short)*2) ; // clear both storages which is why there's a 2
1409 
1410  if(offline) {
1411  memset(track_dta,0xFF,dt_2*(dp+2)*sizeof(short)*2) ; // set to 0xFF!
1412  }
1413 
1414  dta_s = smooth_dta + (dp+2)*dt_2 ; // for smoothed data
1415 
1416  int p1 = blob[ix].p1 ;
1417  int p2 = blob[ix].p2 ;
1418 
1419 #ifdef DO_DBG1
1420  printf("BLOB START: %d X %d\n",dp,dt) ;
1421 #endif
1422  // now let's rip though the data and bring it in our XY structure
1423  for(int p=p1;p<=p2;p++) {
1424  rp_t *rp = get_row_pad(row,p) ;
1425 
1426  if(rp->s1_len==0) continue ;
1427 
1428  u_short *d = rp->s1_data ;
1429 
1430  int pp = p - p1 + 1 ;
1431 
1432  short *adc_p = smooth_dta + dt_2 * pp ;
1433  u_short *track_p = track_dta + dt_2 * pp ;
1434 
1435  for(int j=0;j<rp->s1_len;j++) { // loop over sequnces
1436  int t_cou, t_start ;
1437 
1438  int ixx = *d++ ;
1439 
1440  int b_ix = blob_ix[ixx] ;
1441 
1442 
1443  t_cou = *d++ ;
1444  t_start = *d++ ;
1445 
1446 
1447  if(b_ix != ix) {
1448  if(offline) d += 2*t_cou ;
1449  else d += t_cou ;
1450 
1451  continue ;
1452  }
1453 
1454  blob_t *bl = &(blob[ix]) ;
1455 
1456  // asign ptr
1457  // d now points to ADC
1458  for(int t=0;t<t_cou;t++) {
1459  u_short adc = *d++ ;
1460 
1461  u_short tb = t_start++ ;
1462 
1463  int ttb = tb - bl->t1 + 1 ;
1464  adc_p[ttb] = adc ;
1465 
1466  if(offline) {
1467  track_p[ttb] = *d++ ;
1468  }
1469 
1470  bl->tot_charge += adc ;
1471  bl->pixels++ ;
1472  }
1473  }
1474  }
1475 
1476  if(blob[ix].pixels == 0) { // possible if I blew off pixels due to GG or ADC cuts
1477  //LOG(ERR,"Blob %d: No pixels???",ix) ;
1478  blob[ix].seq_cou = 0 ;
1479  continue ;
1480  }
1481 
1482  f_stat.tm[4] += delta(tm) ;
1483 
1484  // do 3x3 averaging
1485  for(int i=1;i<=dp;i++) {
1486  for(int j=1;j<=dt;j++) {
1487  short *adc_p = smooth_dta + dt_2 * i + j ;
1488 
1489  short sum = *adc_p ; // weight by 2*value!
1490 
1491  if(sum <= 10) { // a peak can't have less than that in ADC counts!
1492  sum = 0 ;
1493  goto do_store ;
1494  }
1495 
1496  // sum up around the real ADC
1497  for(int ii=-1;ii<=1;ii++) {
1498  for(int jj=-1;jj<=1;jj++) {
1499 
1500  int iii = ii*jj ;
1501  if(iii==1 || iii==-1) continue ;
1502 
1503  short *adc_p = smooth_dta + dt_2 * (i+ii) + (j+jj) ;
1504  sum += *adc_p ;
1505  }
1506  }
1507 
1508  do_store: ;
1509 
1510  // and store it in the smoothed ADC
1511  adc_p = dta_s + dt_2 * i + j ;
1512 
1513  *adc_p = sum ;
1514  }
1515  }
1516 
1517  f_stat.tm[5] += delta(tm) ;
1518 
1519  // do peak finding over the smoothed ADC
1520  int peaks_cou = 0 ;
1521 
1522  // I think this should not include the edges but shold go from [2,dp-1] and [2,dt-1]
1523  // I'll see later.
1524 
1525  for(int i=1;i<=dp;i++) {
1526  for(int j=1;j<=dt;j++) {
1527  short *adc_p = dta_s + dt_2 * i + j ; // smoothed
1528 
1529  short adc = *adc_p - 5 ; // the others around MUST be at least 3 ADCs greater!
1530 
1531  // chop off small peaks, below 1; this might cause peaks_cou to be 0 incorrectly!
1532  if(adc < 1) continue ;
1533 
1534  for(int ii=-1;ii<=1;ii++) {
1535  for(int jj=-1;jj<=1;jj++) {
1536  if(ii==0 && jj==0) continue ;
1537 
1538  short *adc_p = dta_s + dt_2 * (i+ii) + (j+jj) ;
1539  short s_adc = *adc_p ;
1540 
1541  if(adc < s_adc) goto skip_calc ;
1542  if(s_adc < 0) goto skip_calc ;
1543  }
1544  }
1545 
1546  #if 0
1547  // Mask out all around?
1548  for(int ii=-1;ii<=1;ii++) {
1549  for(int jj=-1;jj<=1;jj++) {
1550  short *adc_p = dta_s + dt_2 * (i+ii) + (j+jj) ;
1551 
1552  *adc_p = -1 ;
1553  }
1554  }
1555  #endif
1556 
1557  // Aha -- we have a peak!
1558 
1559  *adc_p = -adc ; // -100*adc ; // mark as used!
1560 
1561  peaks[peaks_cou].i = i ;
1562  peaks[peaks_cou].j = j ;
1563  peaks[peaks_cou].adc = adc ;
1564  peaks_cou++ ;
1565 
1566  if(peaks_cou >= MAX_PEAKS) {
1567  LOG(WARN,"Too many peaks %d/%d in row %d",peaks_cou,MAX_PEAKS,row) ;
1568  // At this point I could go over the found peaks and blow off the ones which are too low!
1569  // but later...
1570  goto peaks_done ;
1571  }
1572 
1573  j += 5 ; // skip some timebins (at least 2) as to not have them close together in time
1574 
1575  skip_calc: ;
1576  }
1577  }
1578 
1579  peaks_done: ;
1580 
1581  f_stat.tm[6] += delta(tm) ;
1582 
1583  //LOG(TERR,"Blob %d: peaks %d",ix,peaks_cou) ;
1584  blob_id++ ;
1585 
1586 
1587 
1588  if(peaks_cou <= 1) { // do usual averaging!
1589  double f_charge = 0.0 ;
1590  double f_t_ave = 0.0 ;
1591  double f_p_ave = 0.0 ;
1592  int flags = 0 ;
1593 
1594  int adc_max = 0 ;
1595  int pixels = 0 ;
1596  u_short track_id = 0xFFFF ;
1597  u_short quality = 0 ;
1598 
1599  if(offline) { // various track id things
1600  u_int i_charge = 0 ;
1601  u_int t_charge = 0 ;
1602 
1603  for(int i=1;i<=dp;i++) {
1604 
1605  short *adc_p = smooth_dta + dt_2 * i + 1;
1606  u_short *track_p = track_dta + dt_2 * i + 1 ;
1607 
1608 
1609  for(int j=1;j<=dt;j++) {
1610 
1611  int adc = *adc_p++ ;
1612  u_short t_id = *track_p++ ;
1613 
1614  // count pixels
1615  if(adc) pixels++ ;
1616  else continue ;
1617 
1618  // find the ADC max and asign its trackid
1619  if(adc > adc_max) {
1620  //also asign the trackid!
1621  track_id = t_id ;
1622  adc_max = adc ;
1623  }
1624 
1625  //sum up the pixels as well
1626  i_charge += adc ;
1627 
1628  }
1629  }
1630 
1631  //now sum up the only pixels which belong to the trackid
1632  for(int i=1;i<=dp;i++) {
1633 
1634  short *adc_p = smooth_dta + dt_2 * i + 1;
1635  u_short *track_p = track_dta + dt_2 * i + 1 ;
1636 
1637 
1638  for(int j=1;j<=dt;j++) {
1639 
1640  int adc = *adc_p++ ;
1641  u_short t_id = *track_p++ ;
1642 
1643  if(t_id == track_id) t_charge += adc ;
1644 
1645  }
1646  }
1647 
1648  if(i_charge) {
1649  quality = (u_short)(100.0*(double)t_charge/(double)i_charge+0.5) ;
1650  }
1651 
1652 
1653 
1654  }
1655 
1656  for(int i=1;i<=dp;i++) {
1657  int pad = p1 + i - 1 ;
1658 
1659  //LOG(TERR,"Gain RP %d:%d = %f %d",row,pad,gain_row_pad[row][pad].gain,gain_row_pad[row][pad].flags) ;
1660 
1661  if(use_gain) flags |= gain_row_pad[row][pad].flags ;
1662 
1663  short *adc_p = smooth_dta + dt_2 * i + 1;
1664 
1665  u_int i_charge = 0 ;
1666  u_int i_t_ave = 0 ;
1667 
1668  for(int j=1;j<=dt;j++) {
1669 
1670  int adc = *adc_p++ ;
1671 
1672  if(!adc) continue ;
1673 
1674  i_charge += adc ;
1675  i_t_ave += j * adc ;
1676  }
1677 
1678  if(i_charge==0) continue ;
1679 
1680 
1681 
1682  if(use_gain) {
1683  double corr_charge = (double) i_charge * gain_row_pad[row][pad].gain ;
1684 
1685  f_charge += corr_charge ;
1686  f_t_ave += i_t_ave * gain_row_pad[row][pad].gain + gain_row_pad[row][pad].t0 * corr_charge ;
1687  f_p_ave += i * corr_charge ;
1688  }
1689  else {
1690  f_charge += (double) i_charge ;
1691  f_p_ave += i * (double) i_charge ;
1692  f_t_ave += (double) i_t_ave ;
1693  }
1694  }
1695 
1696  if(f_charge<0.1) {
1697  goto done_peaks;
1698  }
1699 
1700  f_t_ave /= f_charge ;
1701  f_p_ave /= f_charge ;
1702 
1703  f_p_ave += p1 - 1 ;
1704  f_t_ave += blob[ix].t1 - 1 ;
1705 
1706  // and dump out to storage!
1707 
1708  // integerized values
1709  u_int time_c = (u_int)(f_t_ave * 64.0 + 0.5) ;
1710  u_int pad_c = (u_int)(f_p_ave * 64.0 + 0.5) ;
1711  u_int cha = (u_int)(f_charge + 0.5) ;
1712 
1713  // cant happen
1714  //if(cha==0) printf("WTF cha %f\n",f_charge*1000.0) ;
1715 
1716  //extents
1717  u_int tmp_fl ;
1718 
1719  int p_lo = (pad_c/64) - blob[ix].p1 ;
1720  int p_hi = blob[ix].p2 - (pad_c/64) ;
1721 
1722  if(p_lo < 0) p_lo = 0 ;
1723  if(p_hi < 0) p_hi = 0 ;
1724  if(p_lo > 7) p_lo = 7 ;
1725  if(p_hi > 7) p_hi = 7 ;
1726 
1727  tmp_fl = (p_lo<<8)|(p_hi<<11) ;
1728 
1729  int t_lo = (time_c/64) - blob[ix].t1 ;
1730  int t_hi = blob[ix].t2 - (time_c/64) ;
1731 
1732  if(t_lo < 0) t_lo = 0 ;
1733  if(t_hi < 0) t_hi = 0 ;
1734  if(t_lo > 15) t_lo = 15 ;
1735  if(t_hi > 15) t_hi = 15 ;
1736 
1737  tmp_fl |= (t_hi << 4) | t_lo ;
1738 
1739  if(cha > 0x7FFF) cha = 0x8000 | (cha/1024) ;
1740 
1741  if(flags & 3) tmp_fl |= 0x8000 ; // ROW_EDGE
1742  // NEW 11-Jan-2022
1743  else {
1744 
1745  *obuff++ = (time_c << 16) | pad_c ;
1746  *obuff++ = (cha << 16) | tmp_fl ;
1747 
1748  if(words_per_cluster>2) {
1749  *obuff++ = (quality<<16)|track_id ;
1750  }
1751  if(words_per_cluster>3) {
1752  *obuff++ = (pixels<<16)|adc_max ;
1753  }
1754 
1755 #ifdef DO_DBG1
1756 // LOG(TERR,"**** S %d: row %d: %f %f %f",clusters_cou,row,f_p_ave,f_t_ave,f_charge) ;
1757 #endif
1758  clusters_cou++ ;
1759  }
1760  }
1761  else { // multiple peak hauristics
1762  int ip1, ip2 ;
1763  int it1, it2 ;
1764 
1765  //LOG(TERR,"peaks_cou %d",peaks_cou) ;
1766 
1767  for(int pk=0;pk<peaks_cou;pk++) {
1768  double f_charge = 0.0 ;
1769  double f_t_ave = 0.0 ;
1770  double f_p_ave = 0.0 ;
1771  int flags = 0 ;
1772 
1773  int adc_max = 0 ;
1774  int pixels = 0 ;
1775  u_short track_id = 0xFFFF ;
1776  u_short quality = 0 ;
1777 
1778  ip1 = peaks[pk].i - 1 ;
1779  if(ip1 < 1) ip1 = 1 ;
1780  ip2 = peaks[pk].i + 1 ;
1781  if(ip2 > dp) ip2 = dp ;
1782 
1783  it1 = peaks[pk].j - 2 ;
1784  if(it1 < 1) it1 = 1 ;
1785 
1786  it2 = peaks[pk].j + 2 ;
1787  if(it2 > dt) it2 = dt ;
1788 
1789 #ifdef DO_DBG1
1790  printf("from %d:%d %d:%d\n",ip1,ip2,it1,it2) ;
1791 #endif
1792  for(int i=ip1;i<=ip2;i++) {
1793  int pad = p1 + i - 1 ;
1794 
1795  if(use_gain) {
1796  flags |= gain_row_pad[row][pad].flags ;
1797  }
1798 
1799  short *adc_p = smooth_dta + dt_2 * i + it1 ;
1800 
1801  u_int i_charge = 0 ;
1802  u_int i_t_ave = 0 ;
1803 
1804 
1805 
1806  for(int j=it1;j<=it2;j++) {
1807 
1808  int adc = *adc_p++ ;
1809 
1810 
1811 #ifdef DO_DBG1
1812  printf("%d %d = %d\n",i,j,adc) ;
1813 #endif
1814  if(!adc) continue ;
1815  else pixels++ ;
1816 
1817  if(adc>adc_max) adc_max = adc ;
1818 
1819  i_charge += adc ;
1820  i_t_ave += j * adc ;
1821  }
1822 
1823  if(i_charge==0) continue ;
1824 
1825  if(use_gain) {
1826  double corr_charge = (double) i_charge * gain_row_pad[row][pad].gain ;
1827 
1828  f_charge += corr_charge ;
1829  f_t_ave += (double) i_t_ave * gain_row_pad[row][pad].gain + gain_row_pad[row][pad].t0 * corr_charge ;
1830  f_p_ave += i * corr_charge ;
1831  }
1832  else {
1833  f_charge += (double) i_charge ;
1834  f_p_ave += i * (double)i_charge ;
1835  f_t_ave += (double) i_t_ave ;
1836  }
1837  }
1838 
1839 
1840  if(f_charge < 0.1) continue ;
1841 
1842  f_t_ave /= f_charge ;
1843  f_p_ave /= f_charge ;
1844 
1845  f_p_ave += p1 - 1 ;
1846  f_t_ave += blob[ix].t1 - 1 ;
1847 
1848  // and dump out to storage!
1849 
1850  // integerized values
1851  u_int time_c = (u_int)(f_t_ave * 64.0 + 0.5) ;
1852  u_int pad_c = (u_int)(f_p_ave * 64.0 + 0.5) ;
1853  u_int cha = (u_int)(f_charge + 0.5) ;
1854 
1855  //if(cha==0) {
1856  // printf("WTF Multi peak %f\n",f_charge) ;
1857  //}
1858 
1859  //extents
1860  u_int tmp_fl ;
1861 
1862  int p_lo = 1 ;
1863  int p_hi = 1 ;
1864 
1865  tmp_fl = (p_lo<<8)|(p_hi<<11) ;
1866 
1867  int t_lo = 2 ;
1868  int t_hi = 2 ;
1869 
1870  tmp_fl |= (t_hi << 4) | t_lo ;
1871 
1872  if(cha > 0x7FFF) cha = 0x8000 | (0xFFFF & (cha/1024)) ;
1873 
1874  // add flags here
1875  pad_c |= 0x8000 ; // merged flag
1876  if(flags & 3) tmp_fl |= 0x8000 ; // ROW_EDGE
1877  // NEW 11-Jan-2022
1878  else {
1879 
1880 
1881  *obuff++ = (time_c << 16) | pad_c ;
1882  *obuff++ = (cha << 16) | tmp_fl ;
1883 
1884  if(words_per_cluster > 2) {
1885  *obuff++ = (quality<<16)|track_id ;
1886  }
1887  if(words_per_cluster > 3) {
1888  *obuff++ = (pixels<<16)|adc_max ;
1889  }
1890 
1891 #ifdef DO_DBG1
1892 // LOG(TERR,"0x%X 0x%X 0x%X 0x%X - 0x%X 0x%X",pad_c,time_c,cha,tmp_fl, (time_c<16)|pad_c,(cha<<16)|tmp_fl) ;
1893 // LOG(TERR,"**** D %d: %f %f %f",clusters_cou,f_p_ave,f_t_ave,f_charge) ;
1894 #endif
1895  clusters_cou++ ;
1896  }
1897  }
1898  }
1899 
1900 
1901  done_peaks:;
1902  f_stat.tm[7] += delta(tm) ;
1903 
1904 #ifdef DO_DBG1
1905  printf("BLOB start: row %d, peaks %d: %d:%d, %d:%d\n",row,peaks_cou,blob[ix].p1,blob[ix].p2,blob[ix].t1,blob[ix].t2) ;
1906 
1907  for(int j=1;j<=dt;j++) {
1908 
1909  printf("TB %3d ",blob[ix].t1+j-1) ;
1910  for(int i=1;i<=dp;i++) {
1911  short *adc_s_p = dta_s + dt_2 * i + j ; // smoothed
1912 
1913  short *adc_p = smooth_dta + dt_2 * i + j ; // real
1914 
1915  if(*adc_s_p < 0) {
1916  printf("%s%4d(%4d)%s ",ANSI_RED,*adc_p,*adc_s_p,ANSI_RESET) ;
1917  }
1918  else {
1919  printf("%4d(%4d) ",*adc_p,*adc_s_p) ;
1920  }
1921  }
1922  printf("\n") ;
1923  }
1924  fflush(stdout) ;
1925 #endif
1926  }
1927 
1928 // LOG(TERR,"row %d: clusters %d, bytes %d",row,clusters_cou,(obuff-(u_int *)out_store)*4) ;
1929 
1930  return (obuff - out_store) ; // in ints!!!
1931 }
1932 
1933 int itpc_fcf_c::do_row_check(int row)
1934 {
1935 // rp_t *rp = &(row_pad[row][0]) ;
1936 // rp_t *rp = get_row_pad(row,0) ;
1937 
1938 // int rl = rowlen[row] ;
1939  int rl = x_max(row,0) ;
1940 
1941  for(int p=1;p<=rl;p++) { // < is on purpose!!!
1942  int t1_cou, t1_lo ;
1943  u_short *ix1_p ;
1944  rp_t *rp = get_row_pad(row,p) ;
1945 
1946  if(rp->s1_len==0) continue ;
1947 
1948  u_short *d1 = rp->s1_data ;
1949 
1950  for(int i=0;i<rp->s1_len;i++) {
1951  ix1_p = d1++ ;
1952 
1953  t1_cou = *d1++ ; // actually tb_cou
1954  t1_lo = *d1++ ; // actually t_start ;
1955 
1956  if(offline) d1 += 2*t1_cou ;
1957  else d1 += t1_cou ; // advance to next
1958 
1959  if(*ix1_p==0 || *ix1_p==0xFFFF) {
1960  LOG(ERR,"sequence unassigned %d:%d %d:%d = %u [%d]",p,rl,i,rp[p].s1_len,*ix1_p,t1_lo) ;
1961  }
1962  }
1963 
1964  }
1965 
1966 
1967  int max_ix = 0 ;
1968 
1969  for(int i=1;i<blob_cou;i++) {
1970  blob[i].cou = 0 ;
1971  }
1972 
1973  for(int i=1;i<blob_cou;i++) {
1974  blob[blob_ix[i]].cou++ ;
1975  }
1976 
1977  for(int i=1;i<blob_cou;i++) {
1978  LOG(DBG,"blob_ix[%d] = %d",i,blob_ix[i]) ;
1979 
1980  }
1981 
1982  for(int i=1;i<blob_cou;i++) {
1983  if(blob[i].cou) {
1984  LOG(DBG," Blob %d: %d",i,blob[i].cou) ;
1985  max_ix++ ;
1986  }
1987  }
1988 
1989  for(int i=1;i<max_ix;i++) {
1990  if(blob[i].merges) {
1991  LOG(DBG,"blob %d: merged %d",i,blob[i].merges) ;
1992  }
1993  }
1994 
1995 
1996  LOG(DBG,"Final blob count %d/%d",max_ix,blob_cou) ;
1997 
1998  return 0 ;
1999 
2000 }
2001 
daq_sim_cld_x
Definition: daq_dta_structs.h:58
daq_det_gain
Definition: daq_dta_structs.h:74
daq_dta
Definition: daq_dta.h:49
daq_cld
Definition: daq_dta_structs.h:37
Generated by   doxygen 1.8.5