RawUtil.cxx
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////
2 /// \brief Collection of basic utilities for working with raw digits
3 /// \author messier@indiana.edu
4 /// \date June 2011
5 //////////////////////////////////////////////////////////////////////////
6 #include "Utilities/RawUtil.h"
7 
9 
11 #include "Geometry/Geometry.h"
12 #include "RawData/RawDigit.h"
13 
14 #include <cmath>
15 #include <algorithm>
16 
17 //......................................................................
18 
21 {
22  return (d1->TDC()<d2->TDC());
23 }
24 
25 //......................................................................
26 
27 ///
28 /// Sort the list of digits to be in increasing time order, earliest
29 /// to lattest
30 ///
32 {
33  sort(d.begin(), d.end(), digi_sort);
34 }
35 
36 //......................................................................
37 
38 void util::CountXY(const std::vector< art::Ptr<rawdata::RawDigit> >& d,
39  unsigned int i1,
40  unsigned int i2,
41  unsigned int* nx,
42  unsigned int* ny)
43 {
46 
47  *nx = 0;
48  *ny = 0;
49  unsigned int i;
50  for (i=i1; i<=i2; ++i) {
51  unsigned int p = cmap->GetPlane(d[i].get());
52  geo::View_t v = geom->Plane(p)->View();
53  if (v==geo::kX) ++(*nx);
54  else if (v==geo::kY) ++(*ny);
55  }
56 }
57 
58 //......................................................................
59 
60 void util::EventBox(const std::vector< art::Ptr<rawdata::RawDigit> >& d,
61  unsigned int i1,
62  unsigned int i2,
63  unsigned int* plane1x,
64  unsigned int* plane2x,
65  unsigned int* cell1x,
66  unsigned int* cell2x,
67  unsigned int* plane1y,
68  unsigned int* plane2y,
69  unsigned int* cell1y,
70  unsigned int* cell2y)
71 {
74 
75  *plane1x = 999999;
76  *plane1y = 999999;
77  *cell1x = 999999;
78  *cell1y = 999999;
79  *plane2x = 0;
80  *cell2x = 0;
81  *plane2y = 0;
82  *cell2y = 0;
83 
84  unsigned int i;
85  unsigned int p, c;
86  geo::View_t v;
87  for (i=i1; i<=i2; ++i) {
88  p = cmap->GetPlane(d[i].get());
89  v = geom->Plane(p)->View();
90  c = cmap->GetCell(d[i].get());
91  if (v==geo::kX) {
92  if (p<(*plane1x)) (*plane1x) = p;
93  if (p>(*plane2x)) (*plane2x) = p;
94  if (c<(*cell1x)) (*cell1x) = c;
95  if (c>(*cell2x)) (*cell2x) = c;
96  }
97  else if (v==geo::kY) {
98  if (p<(*plane1y)) (*plane1y) = p;
99  if (p>(*plane2y)) (*plane2y) = p;
100  if (c<(*cell1y)) (*cell1y) = c;
101  if (c>(*cell2y)) (*cell2y) = c;
102  }
103  }
104 }
105 
106 //......................................................................
107 
108 ///
109 /// Find time windows with significant hit activity
110 ///
111 void util::TimeSlice(const std::vector< art::Ptr<rawdata::RawDigit> >& d,
112  unsigned int tdcwindow,
113  unsigned int nhit,
114  unsigned int nhitx,
115  unsigned int nhity,
116  std::vector<RawSlice>& slice)
117 {
118  //
119  // If the list can't meet our minimum hit requirements, we're done
120  // before we start. These checks also guarantee that the assumptions
121  // we make later regarding unsigned int arithmetic is correct
122  //
123  if (nhit<1 || nhit<(nhitx+nhity) || d.size()<nhit) return;
124 
125  unsigned int i, j;
126  std::vector<RawSlice> tmpslice;
127  for (i=0; i<d.size(); ++i) {
128  //
129  // Slide forward through the hits until we pass the limit of the
130  // time window
131  //
132  unsigned int dt = 0;
133  for (j=i+1; j<d.size(); ++j) {
134  dt = d[j]->TDC() - d[i]->TDC();
135  if (dt>=tdcwindow) {
136  //
137  // The last hit will have exceeded the time window so back up by
138  // one so that j is the index of the last hit to be included
139  //
140  --j;
141  //
142  // Do we have enough total hits?
143  //
144  if (j>=i && (j-i+1)>=nhit) {
145  //
146  // Do we have enough in the x and y views?
147  //
148  unsigned int nx=0, ny=0;
149  util::CountXY(d, i, j, &nx, &ny);
150  //
151  // Slice meets requirements. Put it in the list.
152  //
153  if (nx>=nhitx && ny>=nhity) {
154  RawSlice rs;
155  rs.first = i;
156  rs.second = j;
157  tmpslice.push_back(rs);
158  }
159  }
160  break; // Jump out of loop on j
161  } // if (dt>=tdcwindow)
162  } // for (j=...
163  //
164  // If we've bumped into the end of the digit list, we're done
165  //
166  if (j==d.size()) break;
167  }
168  //
169  // Merge overlaping slices
170  //
171  for (i=0; i<tmpslice.size(); ++i) {
172  if (tmpslice[i].first==0 && tmpslice[i].second==0) continue;
173  for (j=i+1; j<tmpslice.size(); ++j) {
174  if (tmpslice[j].first <= tmpslice[i].second) {
175  tmpslice[i].second = tmpslice[j].second;
176  tmpslice[j].first = 0;
177  tmpslice[j].second= 0;
178  }
179  }
180  }
181  //
182  // Push which ever slices survived to the output
183  //
184  for (unsigned int i=0; i<tmpslice.size(); ++i) {
185  if (!(tmpslice[i].first==0 && tmpslice[i].second==0)) {
186  slice.push_back(tmpslice[i]);
187  }
188  }
189 }
190 
191 //......................................................................
192 
193 unsigned int
195  int adc_sat,
196  int dt_tdc)
197 {
198  unsigned int i, j;
199 
200  //
201  // For efficiency, time sort the hits in advance
202  //
203  util::TimeSort(rd);
204 
205  //
206  // Flag all hits as good, look for ones that might be bad
207  //
208  std::vector<bool> isgood(rd.size());
209  for (i=0; i<rd.size(); ++i) isgood[i] = true;
210  for (i=0; (i+1)<rd.size(); ++i) {
211  //
212  // If the hit i is not in saturation, we can skip to the next
213  // hit
214  //
215  if (rd[i]->ADC()<adc_sat) continue;
216  //
217  // Get ID number for the FEB ID that has the hit in saturation
218  //
219  unsigned int chi = rd[i]->DaqChannel();
220  unsigned int diblocki = (chi&0x0FC00000)>>22;
221  unsigned int dcmi = (chi&0x003F0000)>>16;
222  unsigned int febi = (chi&0x0000FF00)>>8;
223  for (j=i+1; j<rd.size(); ++j) {
224  //
225  // Don't need to look beyond the time limit
226  //
227  if ((rd[j]->TDC()-rd[i]->TDC())>dt_tdc) break;
228  //
229  // Hit j is in time range of another hit which is in
230  // saturation. Check if hit j comes from the same FEB as hit
231  // i
232  //
233  unsigned int chj = rd[j]->DaqChannel();
234  unsigned int diblockj = (chj&0x0FC00000)>>22;
235  unsigned int dcmj = (chj&0x003F0000)>>16;
236  unsigned int febj = (chj&0x0000FF00)>>8;
237  bool sameFEB =
238  (diblocki==diblockj) && (dcmi==dcmj) && (febi==febj);
239  if (sameFEB) isgood[j] = false;
240  }
241  }
242  //
243  // Push all the good hits into the final list
244  //
245  std::vector< art::Ptr<rawdata::RawDigit> > good;
246  good.reserve(rd.size());
247  for (i=0; i<rd.size(); ++i) {
248  if (isgood[i]) good.push_back(rd[i]);
249  }
250  rd.swap(good);
251  return rd.size();
252 }
253 
254 ////////////////////////////////////////////////////////////////////////
int32_t TDC() const
The time of the last baseline sample.
Definition: RawDigit.h:94
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
static bool digi_sort(const art::Ptr< rawdata::RawDigit > d1, const art::Ptr< rawdata::RawDigit > d2)
Collection of basic utilities for working with raw digits.
Definition: RawUtil.cxx:19
const char * p
Definition: xmltok.h:285
Vertical planes which measure X.
Definition: PlaneGeo.h:28
nhit
Definition: demo1.py:25
const PlaneGeo * Plane(unsigned int i) const
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
std::pair< unsigned int, unsigned int > RawSlice
Definition: RawUtil.h:23
void CountXY(const std::vector< art::Ptr< rawdata::RawDigit > > &d, unsigned int i1, unsigned int i2, unsigned int *nx, unsigned int *ny)
Count the number of digits in each detector view.
Definition: RawUtil.cxx:38
Definition: Cand.cxx:23
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
void TimeSlice(const std::vector< art::Ptr< rawdata::RawDigit > > &d, unsigned int dt_tdc, unsigned int nhit, unsigned int nhitx, unsigned int nhity, std::vector< RawSlice > &slice)
Find windows in time that have significant activity in the detector.
Definition: RawUtil.cxx:111
Float_t d
Definition: plot.C:236
unsigned short GetPlane(const rawdata::RawDigit *dig)
Definition: CMap.cxx:285
const double j
Definition: BetheBloch.cxx:29
Functions for getting information about a collection of raw hits.
unsigned int FilterFEBFlash(std::vector< art::Ptr< rawdata::RawDigit > > &rd, int adc_sat=3400, int dt_tdc=1280)
Filter hits that are "FEB flash" candidates; that is, they occur within a specified time of an FEB hi...
Definition: RawUtil.cxx:194
void geom(int which=0)
Definition: geom.C:163
unsigned short GetCell(const rawdata::RawDigit *dig)
Definition: CMap.cxx:327
void EventBox(const std::vector< art::Ptr< rawdata::RawDigit > > &d, unsigned int i1, unsigned int i2, unsigned int *plane1x, unsigned int *plane2x, unsigned int *cell1x, unsigned int *cell2x, unsigned int *plane1y, unsigned int *plane2y, unsigned int *cell1y, unsigned int *cell2y)
Find boxes in plane/cell units that contain all the hits.
Definition: RawUtil.cxx:60
void TimeSort(std::vector< art::Ptr< rawdata::RawDigit > > &d)
Arrange the list of raw hits in time order (early to late)
Definition: RawUtil.cxx:31
Definition: fwd.h:28
Encapsulate the geometry of one entire detector (near, far, ndos)