noe_module.cc
Go to the documentation of this file.
1 #include <signal.h>
2 
3 #include <vector>
4 
5 // For getting the event count when the file is opened
6 #include "TTree.h"
9 
12 
13 #include "RecoBase/CellHit.h"
14 
15 // For tracks
16 #include "Geometry/Geometry.h"
17 #include "RecoBase/Track.h"
18 #include "RecoBase/Vertex.h"
19 
20 #include "func/main.h"
21 #include "func/event.h"
22 
23 using std::vector;
24 
25 static uint64_t MAXHITS = 4*1024*1024*1024ULL/sizeof(hit);
26 
28 vector<noeevent> theevents;
29 
30 namespace noe {
31 class noe : public art::EDProducer {
32  public:
33  noe(fhicl::ParameterSet const& pset);
34  ~noe();
35  void produce(art::Event& evt);
36  void endJob();
37 
38  // Used to get the number of events in the file
39  void respondToOpenInputFile(art::FileBlock const &fb);
40 
41  // The art labels for tracks and vertices that we are going to display, or
42  // the empty string to display none.
46 };
47 
49 {
50  fCellHitLabel = pset.get< std::string >("cellhit_label");
51  fTrackLabel = pset.get< std::string >("track_label");
52  fVertexLabel= pset.get< std::string >("vertex_label");
53 }
54 
55 noe::~noe() { }
56 
57 void noe::respondToOpenInputFile(art::FileBlock const &fb)
58 {
59  // Get the number of events as soon as the file opens. This looks
60  // really fragile. It gets the number of entries in *some* tree, which
61  // at the moment that I'm testing this turns out to be the right one,
62  // but EDProducer::respondToOpenInputFile is totally undocumented as
63  // far as I can see.
64  //
65  // Anyway, if this is the wrong number, it just means that the status
66  // display will be wrong about what fraction of the file is loaded.
67  //
68  // If the job has more than one file, we don't know that until the
69  // second one triggers this function. This means the user is going to
70  // be disappointed when the percent-loaded is rewound. Oh well.
71  //
72  // If the user gave -n on the command line to limit the number of
73  // events, we don't pick that up either, so the status bar will just
74  // stop increasing when we stop reading without reaching 100%. Chris
75  // Backhouse says "you can introspect what the fcl configuration
76  // was for other modules in the path (see CAFMaker for an example
77  // of trying to dig out genie settings) so maybe you can get at the
78  // InputSource config (which is where that value goes)". So TODO.
79  auto const* rfb = dynamic_cast<art::RootFileBlock const*>(&fb);
80 
81  theevents.reserve(theevents.capacity() + rfb->tree()->GetEntries());
82 }
83 
84 void noe::endJob()
85 {
86  if(!theevents.empty()) realmain(true);
87 }
88 
89 // Inject a test event with all FD cells hit
90 __attribute__((unused)) static void add_test_fd_event()
91 {
92  noeevent ev;
93  ev.nevent = 0;
94  ev.nrun = 0;
95  ev.nsubrun = 0;
96 
97  for(unsigned int c = 0; c < 12 * 32; c++){
98  for(unsigned int p = 0; p < 32 * 28; p++){
99  hit thehit;
100  thehit.cell = c;
101  thehit.plane = p;
102  thehit.adc = (c*p)%124;
103  thehit.tdc = ((c*p)%432)*4;
104  ev.addhit(thehit);
105  }
106  }
107  theevents.push_back(ev);
108 }
109 
110 // Inject a test event with all ND cells hit
111 __attribute__((unused)) static void add_test_nd_event()
112 {
113  noeevent ev;
114  ev.nevent = 0;
115  ev.nrun = 0;
116  ev.nsubrun = 0;
117 
118  for(unsigned int c = 0; c < 3 * 32; c++){
119  for(unsigned int p = 0; p < 2 * (8 * 12 + 11); p++){
120  if(p >= 2 * 8 * 12 && c >= 2 * 32 && p%2 == 0) continue;
121  hit thehit;
122  thehit.cell = c;
123  thehit.plane = p;
124  thehit.adc = (c*p)%1234;
125  thehit.tdc = ((c*p)% 234)*4;
126  ev.addhit(thehit);
127  }
128  }
129  theevents.push_back(ev);
130 }
131 
132 static vector< vector<float> > xcell_z;
133 static vector< vector<float> > ycell_z;
134 static vector<float> xplane_z;
135 static vector<float> yplane_z;
136 static vector< vector<float> > xcell_x;
137 static vector< vector<float> > ycell_y;
138 
139 // Given a position in 3-space, return a plane and cell that is reasonably
140 // close in the given view assuming a regular detector geometry.
142  const double x, const double y, const double z, const geo::View_t view)
143 {
144  cppoint ans;
145 
146  // (1) First the plane
147  {
148  vector<float> & tplane_z = view == geo::kX?xplane_z:yplane_z;
149 
150  const vector<float>::iterator pi
151  = std::upper_bound(tplane_z.begin(), tplane_z.end(), z);
152 
153  const int add = (view == geo::kX);
154 
155  if(pi == tplane_z.end()){
156  // No plane had as big a 'z' as this, so use the last plane.
157  ans.plane = (tplane_z.size()-1) * 2 + add;
158  }
159  else if(pi == tplane_z.begin()){
160  // 'z' was smaller than the first plane, so that's the closest.
161  ans.plane = add;
162  }
163  else{
164  // Check which is closer, the first plane bigger than 'z', or the
165  // previous one.
166  const float thisone = fabs( *pi - z);
167  const float previous = fabs( *(pi - 1) - z);
168  if(thisone < previous) ans.plane = (pi - tplane_z.begin())*2 + add;
169  else ans.plane = (pi - tplane_z.begin())*2 + add - 2;
170  }
171  }
172 
173  // (2) Now find the cell
174  {
175  vector<float> & cells = (view == geo::kX? xcell_x:ycell_y)[ans.plane/2];
176 
177  const float t = (view == geo::kX? x: y);
178  const vector<float>::iterator ci
179  = std::upper_bound(cells.begin(), cells.end(), t);
180 
181  if(ci == cells.end()){
182  // No cell had as big a 't' as this, so use the last cell.
183  ans.cell = cells.size()-1;
184  }
185  else if(ci == cells.begin()){
186  // 't' was smaller than the first cell, so that's the closest.
187  ans.cell = 0;
188  }
189  else{
190  // Check which is closer
191  const float thisone = fabs( *ci - t);
192  const float previous = fabs( *(ci - 1) - t);
193  if(thisone < previous) ans.cell = (ci - cells.begin());
194  else ans.cell = (ci - cells.begin()) - 1;
195  }
196  }
197 
198  return ans;
199 }
200 
202 {
203  for(unsigned int pl = 0; pl < geo->NPlanes(); pl++){
204  const geo::PlaneGeo * const plane = geo->Plane(pl);
206  vector<float> cell_z, cell_t;
207  for(unsigned int ce = 0; ce < plane->Ncells(); ce++){
208  double cellcenter[3], dum[3];
209  geo->CellInfo(pl, ce, &view, cellcenter, dum);
210  cell_z.push_back(cellcenter[2]);
211  cell_t.push_back(geo::kX?cellcenter[0]:cellcenter[1]);
212  }
213 
214  (view == geo::kX?xcell_z:ycell_z).push_back(cell_z);
215  (view == geo::kX?xcell_x:ycell_y).push_back(cell_t);
216 
217  vector<float> & pz = (view == geo::kX?xplane_z:yplane_z);
218  pz.push_back(cell_z[cell_z.size()/2]);
219  }
220 }
221 
222 // Given a Cartesian position, tp, representing a track point, return the
223 // position in floating-point plane and cell number for both views where an
224 // integer means the cell center.
225 static std::pair<cppoint, cppoint> cart_to_cp(
226  art::ServiceHandle<geo::Geometry> & geo, const TVector3 &tp)
227 {
228  // With this lookup table (which isn't really a lookup table), finding
229  // track point uses ~15% of the time spent loading events. Probably
230  // could go faster, although that's not too bad.
231  {
232  static int first = true;
233  if(first) build_cell_lookup_table(geo);
234  first = false;
235  }
236 
237  // For each view, first find a plane and cell which is probably the closest
238  // one, or maybe one of the several closest. Then ask the geometry where that
239  // cell is and store the difference in the fractional part of the cppoint.
240  // This is right up to the difference between the mean plane and cell
241  // spacings and the actual spacing near the requested point. For purposes of
242  // the event display, it's fine.
243 
244  // Exact values are not very important
245  const double meanplanesep = 6.6681604;
246  const double meancellsep = 3.9674375;
247 
248  std::pair<cppoint, cppoint> answer;
249  answer.first = get_int_plane_and_cell(tp.X(), tp.Y(), tp.Z(), geo::kX);
250 
251  double cellcenter[3], dum[3];
252  geo::View_t dumv;
253 
254  geo->CellInfo(answer.first.plane, answer.first.cell, &dumv, cellcenter, dum);
255  answer.first.fcell = (tp.X() - cellcenter[0])/meancellsep;
256  answer.first.fplane = (tp.Z() - cellcenter[2])/meanplanesep;
257 
258  // Could optimize this by only checking the nearest two planes to the one
259  // found above, but I bet that geo::CellInfo is the hot spot, not
260  // std::upper_bound.
261  answer.second = get_int_plane_and_cell(tp.X(), tp.Y(), tp.Z(), geo::kY);
262 
263  geo->CellInfo(answer.second.plane, answer.second.cell, &dumv, cellcenter, dum);
264  answer.second.fcell = (tp.Y() - cellcenter[1])/meancellsep;
265  answer.second.fplane = (tp.Z() - cellcenter[2])/meanplanesep;
266 
267  return answer;
268 }
269 
270 void noe::produce(art::Event& evt)
271 {
272  signal(SIGINT, SIG_DFL); // just exit on Ctrl-C
273 
274  if(theevents_memory_limited) return;
275 
277 
278  if(!evt.getByLabel(fCellHitLabel, cellhits)){
279  fprintf(stderr, "NOE needs CellHits with label \"%s\" as configured in your "
280  "FCL, but event %d doesn't have those.\n", fCellHitLabel.c_str(), evt.event());
281  return;
282  }
283 
285  if(fTrackLabel != "" && !evt.getByLabel(fTrackLabel, tracks)){
286  fprintf(stderr,
287  "Warning: No tracks found with label \"%s\"\n", fTrackLabel.c_str());
288  fTrackLabel = "";
289  }
290 
292  if(fVertexLabel != "" && !evt.getByLabel(fVertexLabel, vertices)){
293  fprintf(stderr,
294  "Warning: No vertices found with label \"%s\"\n", fVertexLabel.c_str());
295  fVertexLabel = "";
296  }
297 
298  // Not needed for hits, just for reco. Aggressively don't load the
299  // Geometry if it isn't needed.
301  fVertexLabel == "" && fTrackLabel == ""? NULL:
303 
304 #if 0
305  if(theevents.empty()) add_test_nd_event();
306 #endif
307 
308  noeevent ev;
309  ev.nevent = evt.event();
310  ev.nrun = evt.run();
311  ev.nsubrun = evt.subRun();
312 
313  // When we're reading in an event, the GUI is unresponsive. This is
314  // a consequence of how we're working around art's design choices.
315  // But this loop is not the bottleneck. The delay is inside art, so
316  // there's no way to put hooks in the middle of it to keep the GUI
317  // responsive.
318  static uint64_t TOTHITS = 0;
319  for(unsigned int i = 0; i < cellhits->size(); i++){
320  const rb::CellHit & c = (*cellhits)[i];
321  hit thehit;
322  thehit.cell = c.Cell();
323  thehit.plane = c.Plane();
324  thehit.adc = c.ADC();
325  thehit.tdc = c.TDC();
326  thehit.tns = c.TNS();
327  thehit.good_tns = c.GoodTiming();
328  ev.addhit(thehit);
329  TOTHITS++;
330  }
331 
332  for(unsigned int i = 0; tracks.isValid() && i < tracks->size(); i++){
333  track thetrack;
334  thetrack.startx = 10*(*tracks)[i].Start().X();
335  thetrack.starty = 10*(*tracks)[i].Start().Y();
336  thetrack.startz = 10*(*tracks)[i].Start().Z();
337  thetrack.stopx = 10*(*tracks)[i].Stop ().X();
338  thetrack.stopy = 10*(*tracks)[i].Stop ().Y();
339  thetrack.stopz = 10*(*tracks)[i].Stop ().Z();
340  thetrack.tns = (*tracks)[i].MeanTNS();
341  thetrack.time = (*tracks)[i].MeanTNS()/1000*64.; // translate to TDC
342  for(unsigned int c = 0; c < (*tracks)[i].NCell(); c++){
343  hit thehit;
344  thehit.cell = (*tracks)[i].Cell(c)->Cell();
345  thehit.plane = (*tracks)[i].Cell(c)->Plane();
346  thetrack.hits.push_back(thehit);
347  }
348  for(unsigned int p = 0; p < (*tracks)[i].NTrajectoryPoints(); p++){
349  const TVector3 & tp = (*tracks)[i].TrajectoryPoint(p);
350  const std::pair<cppoint, cppoint> tps = cart_to_cp(*geo, tp);
351  thetrack.traj[geo::kX].push_back(tps.first);
352  thetrack.traj[geo::kY].push_back(tps.second);
353  }
354  ev.addtrack(thetrack);
355  }
356 
357  for(unsigned int i = 0; vertices.isValid() && i < vertices->size(); i++){
358  vertex thevertex;
359  const std::pair<cppoint, cppoint> cp = cart_to_cp(*geo, (*vertices)[i].GetXYZ());
360  thevertex.pos[0] = cp.first;
361  thevertex.pos[1] = cp.second;
362  thevertex.posx = 10*(*vertices)[i].GetX();
363  thevertex.posy = 10*(*vertices)[i].GetY();
364  thevertex.posz = 10*(*vertices)[i].GetZ();
365  thevertex.tns = (*vertices)[i].GetT();
366  thevertex.time = (*vertices)[i].GetT()/1000*64; // translate to TDC
367  ev.addvertex(thevertex);
368  }
369 
370  if(TOTHITS > MAXHITS){
371  fprintf(stderr, "Read over %lu hits. Will load no more events!\n", MAXHITS);
372  theevents.shrink_to_fit();
374  }
375  else{
376  theevents.push_back(ev);
377  }
378 
379  realmain(false);
380 }
381 
383 
384 }
static vector< vector< float > > xcell_z
Definition: noe_module.cc:132
float TNS() const
Definition: CellHit.h:46
static std::pair< cppoint, cppoint > cart_to_cp(art::ServiceHandle< geo::Geometry > &geo, const TVector3 &tp)
Definition: noe_module.cc:225
int16_t stopx
Definition: event.h:21
SubRunNumber_t subRun() const
Definition: Event.h:72
vector< noeevent > theevents
Definition: noe_module.cc:28
Definition: event.h:34
short posz
Definition: event.h:36
int32_t TDC() const
The time of the last baseline sample.
Definition: RawDigit.h:94
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
int32_t stopz
Definition: event.h:25
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Definition: CellHit.h:39
uint32_t nsubrun
Definition: event.h:45
uint16_t plane
Definition: event.h:2
const char * p
Definition: xmltok.h:285
static vector< vector< float > > ycell_y
Definition: noe_module.cc:137
void addhit(const hit &h)
Definition: event.h:83
int32_t tdc
Definition: event.h:5
uint16_t plane
Definition: event.h:11
Vertical planes which measure X.
Definition: PlaneGeo.h:28
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
std::string fVertexLabel
Definition: noe_module.cc:45
Definition: event.h:41
uint16_t cell
Definition: event.h:2
Definition: event.h:19
cppoint pos[2]
Definition: event.h:35
float tns
Definition: event.h:38
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
std::vector< cppoint > traj[2]
Definition: event.h:31
float tns
Definition: event.h:6
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
void CellInfo(unsigned int ip, unsigned int ic, View_t *view=0, double *pos=0, double *dpos=0) const
static vector< vector< float > > xcell_x
Definition: noe_module.cc:136
void addvertex(const vertex &v)
Definition: event.h:78
uint32_t nevent
Definition: event.h:45
bool isValid() const
Definition: Handle.h:189
short posy
Definition: event.h:36
unsigned short Cell() const
Definition: CellHit.h:40
uint32_t nrun
Definition: event.h:45
bool good_tns
Definition: event.h:4
int16_t startx
Definition: event.h:21
__attribute__((unused)) static void add_test_fd_event()
Definition: noe_module.cc:90
T get(std::string const &key) const
Definition: ParameterSet.h:231
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
int evt
base_types push_back(int_type())
int16_t adc
Definition: event.h:3
static void build_cell_lookup_table(art::ServiceHandle< geo::Geometry > &geo)
Definition: noe_module.cc:201
int32_t startz
Definition: event.h:25
static vector< float > xplane_z
Definition: noe_module.cc:134
static vector< float > yplane_z
Definition: noe_module.cc:135
EventNumber_t event() const
Definition: Event.h:67
static uint64_t MAXHITS
Definition: noe_module.cc:25
Vertex location in position and time.
z
Definition: test.py:28
std::vector< hit > hits
Definition: event.h:30
std::string fCellHitLabel
Definition: noe_module.cc:43
short posx
Definition: event.h:36
static vector< vector< float > > ycell_z
Definition: noe_module.cc:133
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
int32_t time
Definition: event.h:27
float tns
Definition: event.h:28
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Definition: structs.h:12
void realmain(const bool have_read_all)
Definition: main.cxx:962
static cppoint get_int_plane_and_cell(const double x, const double y, const double z, const geo::View_t view)
Definition: noe_module.cc:141
unsigned int NPlanes() const
bool GoodTiming() const
Definition: CellHit.h:48
RunNumber_t run() const
Definition: Event.h:77
Definition: event.h:9
Helper for AttenCurve.
Definition: Path.h:10
add("abs", expr_type(int_type()), expr_type(int_type()))
int16_t starty
Definition: event.h:21
void addtrack(const track &t)
Definition: event.h:73
uint16_t cell
Definition: event.h:11
Encapsulate the geometry of one entire detector (near, far, ndos)
int32_t time
Definition: event.h:37
bool theevents_memory_limited
Definition: noe_module.cc:27
std::string fTrackLabel
Definition: noe_module.cc:44
int16_t stopy
Definition: event.h:21
enum BeamMode string