LiveGeometry_service.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file LiveGeometry_service.cc
3 /// \brief
4 /// \author bays@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
7 #include <algorithm>
8 #include <cassert>
9 #include <iostream>
10 
11 // Framework includes
15 
16 // NOvA includes
17 #include "DAQChannelMap/DAQChannelMap.h"
19 #include "Geometry/Geometry.h"
20 #include "Geometry/LiveGeometry.h"
22 #include "GeometryObjects/Geo.h"
25 
26 #include "TGeoManager.h"
27 
28 namespace geo
29 {
30  //......................................................................
32  fParams(params())
33  {
34  // tell the activity registry to check in after runs have been loaded
35  // into the event processing so that we can update configurations if necessary
36  // SetupLiveGeo must be called before anything will work; will be called automatically for all existing functions
37  // If making a new function, call by hand at beginning or make sure IsPointLive(Info) or ProjectedDistance is called first
38  // as they also call it
39  reg.sPostBeginSubRun.watch(this, &LiveGeometry::postBeginSubRun); // beware of race conditions with RH and BC.
41  fSetup = false;
42  }
43 
44  //......................................................................
46  {
47  if(fSetup)
48  mf::LogInfo("LiveGeometry") << "LiveGeometry found "
49  << fTotalDropouts << " dropped DCMs out of "
50  << fNEvents << " events.";
51  }
52 
53  //......................................................................
55  {
56  // if we change runs, redo relevant parts of Setup (check RH and bad channels for the run). But, wait until we are using LiveGeo (so never do it for first run)
58  if(fSetup) {
59  for(int i = 0; i<12; i++) { // initialize DCMs
60  for(int j = 0; j<14; j++) {
61  fDCMs[i][j] = 1;
62  }
63  }
64  // if(RunHist->GetDataType()!=1) { // don't check RH or Bad Channels for MC
65  GetRHInfo(); // get bad boxes based on RH (what is 'uninstrumented' or fDCMs=2)
66  if(fParams.CheckBadChannels()) GetBCInfo(); // get bad boxes based on BadChannels (What is 'bad' or fDCMs=0)
67  // }
68  }
69  }// end of preBeginRun
70 
71  //----------------------------------------------------------
73  {
75 
76  fWall = 0;
77  fEvFlag = false;
78 
79  }
80 
81  //...................................................................
82  bool LiveGeometry::IsPointLiveInfo(TVector3 vertex, int &info) { // are we in a good part of the detector, or a bad/dead one? return true = good
83  if(!fEvFlag) GetEvInfo(); // setup hasn't been called yet!
84 
85  info = 0; // outside detector
86  bool badxz=false, badyz = false, edgexz=false, edgeyz=false;
87 
88  if(vertex.X() > fEdgeXPlusXZ || vertex.X() < fEdgeXMinusXZ || // are we even inside the detector?
89  vertex.Y() > fEdgeTopYZ || vertex.Y() < fEdgeBottomYZ ||
90  vertex.Z() < std::min(fEdgeFrontXZ,fEdgeFrontYZ) || vertex.Z() > std::max(fEdgeBackXZ,fEdgeBackYZ) ){
91  LOG_DEBUG("LiveGeometry") << "outside of the detector";
92  return false;
93  }
94 
95 
96  for(unsigned int i = 0; i < fBadBoxXZ.size(); i++) { // check XZ view
97  if(fBadBoxXZ[i].Z1 <= vertex.Z() && fBadBoxXZ[i].Z2 >= vertex.Z() &&
98  fBadBoxXZ[i].XY1 <= vertex.X() && fBadBoxXZ[i].XY2 >= vertex.X()) {
99  if(fBadBoxXZ[i].edge == true) edgexz = true;
100  badxz = true;
101  }
102  }
103 
104  for(unsigned int i = 0; i < fBadBoxYZ.size(); i++) { // check YZ view
105  if(fBadBoxYZ[i].Z1 <= vertex.Z() && fBadBoxYZ[i].Z2 >= vertex.Z() &&
106  fBadBoxYZ[i].XY1 <= vertex.Y() && fBadBoxYZ[i].XY2 >= vertex.Y()) {
107  if(fBadBoxYZ[i].edge == true) edgeyz = true;
108  badyz = true;
109  }
110  }
111 
112  if(badxz&&badyz) { // both views bad
113  if(edgexz&&edgeyz) info = 5; // and both are edges
114  if(edgexz||edgeyz) info = 4; // only one is an edge
115  else info = 3; // neither are edges
116  return false;
117  }
118  else if(badxz&&edgexz) { // one view bad and edge
119  info = 2;
120  return false;
121  }
122  else if(badyz&&edgeyz) { // one view bad and edge
123  info = 2;
124  return false;
125  }
126  else if(badxz||badyz) { // one view bad, neither are edge
127  info = 1;
128  return false;
129  }
130 
131  return true;
132  }
133 
134  //...................................................................
136  int info;
137  bool live = this->IsPointLiveInfo(vertex, info);
138  if(!live)
139  LOG_DEBUG("LiveGeometry")
140  << "(" << vertex.X()
141  << ", " << vertex.Y()
142  << ", " << vertex.Z()
143  << ") is not live, info: "
144  << info;
145 
146  return live;
147  }
148 
149  //...................................................................
150  int LiveGeometry::ProjectedWall(TVector3 vertex, TVector3 dir) { // what wall would we exit assuming the downstream instrumented detector ends?
152  double tdist, mindist = 9999999;
153  int wall = 0, twall = 0;
155 
156  if(dir.Z() > 0) {
157  tdist = (geom->DetLength()*(float(fBackDownstream)/float(fNumDiblocks))-vertex.Z())/dir.Z(); // going +z project to back
158  twall = 6;
159  }
160  else if(dir.Z() < 0) {
161  tdist = (geom->DetLength()*(float(fFrontDownstream)/float(fNumDiblocks))-vertex.Z())/dir.Z(); // going -z project to front
162  twall = 5;
163  }
164  else { // dz is zero. This means something is wrong with the track
165  return 0;
166  }
167  if(tdist < mindist) {
168  mindist = tdist;
169  wall = twall;
170  }
171 
172  if(dir.Y() > 0) {
173  tdist = (fEdgeTopYZ-vertex.Y())/dir.Y(); // going +y project to top
174  twall = 4;
175  }
176  else if(dir.Y() < 0) {
177  tdist = (fEdgeBottomYZ-vertex.Y())/dir.Y(); // going -y project to bottom
178  twall = 3;
179  }
180  else tdist = 999999; // not moving in Y? May be 2D track... we can still proceed though
181  if(tdist < mindist) {
182  mindist = tdist;
183  wall = twall;
184  }
185 
186  if(dir.X() > 0) {
187  tdist = (fEdgeXPlusXZ-vertex.X())/dir.X(); // going +x project +x wall
188  twall = 2;
189  }
190  else if(dir.X() < 0) {
191  tdist = (fEdgeXMinusXZ-vertex.X())/dir.X(); // going -x project to -x wall
192  twall = 1;
193  }
194  else tdist = 999999; // not moving in X? may be 2D track... we can still proceed though
195  if(tdist < mindist) {
196  mindist = tdist;
197  wall = twall;
198  }
199 
200  return wall;
201 
202  }
203 
204  //...................................................................
205  double LiveGeometry::ProjectedDistanceToEdge(TVector3 vertex, TVector3 dir) {
206  double distToEdge, distToNextBad, distDead;
207  if(ProjectedDistance(vertex, dir, distToEdge, distToNextBad, distDead)) return(distToEdge);
208  else return -5;
209  }
210 
211  //...................................................................
213  double distToEdge, distToNextBad, distDead;
214  if(!IsPointLive(vertex)) return -5;
215  if(ProjectedDistance(vertex, dir, distToEdge, distToNextBad, distDead)) return(distToNextBad);
216  else return -5;
217  }
218 
219  //...................................................................
221  double distToEdge, distToNextBad, distDead;
222  if(ProjectedDistance(vertex, dir, distToEdge, distToNextBad, distDead)) return(distToEdge-distDead);
223  else return -5;
224  }
225 
226  //...................................................................
227  bool LiveGeometry::ProjectedDistance(TVector3 vertex, TVector3 dir, double &distToEdge, double &distToNextBad, double &distDead) {
228 
229  if(!fEvFlag) GetEvInfo(); // setup hasn't been called yet!
230 
231  float minedgedist = 999999.0, minbaddist=99999.0, tdist; // initial values for minima and temporary distance variable
232  double r[3], d[3];
233  vertex.GetXYZ(r);
234  dir.GetXYZ(d);
235  distDead = 0;
236  fWall = 0;
237 
238  if(IsPointLive(vertex)==false) { // want a negative number here that actually reflects distance to live edge; for simplicity use fully good region; if there are dropped DCMs this may be off
239  distToEdge = -999.;
240  distToNextBad = -999.;
241  distDead = 0;
242 
243  if(vertex.Z() < InstrumentedDetFront()) { // vertex in front of detector
244  if(dir.Z() >= 0) return false; // going wrong way!
245  tdist = (vertex.Z() - InstrumentedDetFront()) / dir.Z(); // work in positive numbers, make negative at end
246  if(tdist < minedgedist) minedgedist = tdist; // z check; want smallest positive number here
247  }
248 
249  if(vertex.Z() > InstrumentedDetBack()) { // vertex behind detector
250  if(dir.Z() <= 0) return false; // wrong direction
251  tdist = (vertex.Z() - InstrumentedDetBack()) / dir.Z();
252  if(tdist < minedgedist) minedgedist = tdist; // z check; want smallest positive number here
253  }
254 
255  if(vertex.Y() > fEdgeTopYZ) {
256  if(dir.Y() <= 0) return false;
257  tdist = -(fEdgeTopYZ-vertex.Y())/dir.Y(); // project to top
258  if(tdist < minedgedist) minedgedist = tdist; // y check
259  }
260  if(vertex.Y() < fEdgeBottomYZ) {
261  if(dir.Y() >= 0) return false;
262  tdist = -(fEdgeBottomYZ-vertex.Y())/dir.Y(); // project to bottom
263  if(tdist < minedgedist) minedgedist = tdist; // y check
264  }
265 
266  if(vertex.X() > fEdgeXPlusXZ) {
267  if(dir.X() <= 0) return false;
268  tdist = -(fEdgeXPlusXZ-vertex.X())/dir.X();
269  if(tdist < minedgedist) minedgedist = tdist; // x check
270  }
271 
272  if(vertex.X() < fEdgeXMinusXZ) {
273  if(dir.X() >= 0) return false;
274  tdist = -(fEdgeXMinusXZ-vertex.X())/dir.X();
275  if(tdist < minedgedist) minedgedist = tdist; // x check
276  }
277 
278  distToEdge = -minedgedist;
279  distToNextBad = -minedgedist;
280  distDead = 0;
281 
282  return false;
283  } // end vertex outside live detector
284 
285  // first check distances to full detector boundaries
286  if(dir.Z() > 0) tdist = (std::max(fEdgeBackYZ,fEdgeBackXZ)-vertex.Z())/dir.Z(); // going +z project to back
287  else if(dir.Z() < 0) tdist = (std::min(fEdgeFrontYZ,fEdgeFrontXZ)-vertex.Z())/dir.Z(); // going -z project to front
288  else { // dz is zero. This means something is wrong with the track
289  distToEdge = -9999.;
290  distToNextBad = -9999.;
291  distDead = -9999.;
292  return false;
293  }
294  if(tdist < minedgedist) minedgedist = tdist;
295 
296  if(dir.Y() > 0) tdist = (fEdgeTopYZ-vertex.Y())/dir.Y(); // going +y project to top
297  else if(dir.Y() < 0) tdist = (fEdgeBottomYZ-vertex.Y())/dir.Y(); // going -y project to bottom
298  else tdist = 999999; // not moving in Y? May be 2D track... we can still proceed though
299  if(tdist < minedgedist) minedgedist = tdist;
300 
301  if(dir.X() > 0) tdist = (fEdgeXPlusXZ-vertex.X())/dir.X(); // going +x project +x wall
302  else if(dir.X() < 0) tdist = (fEdgeXMinusXZ-vertex.X())/dir.X(); // going -x project to -x wall
303  else tdist = 999999; // not moving in X? may be 2D track... we can still proceed though
304  if(tdist < minedgedist) minedgedist = tdist;
305 
306  // now check distances to bad boxes
307 
308  double xyz[3];
309 
310  for(unsigned int i = 0; i < fBadBoxXZ.size(); i++) { // look first in XZ view
311  geo::ProjectToBoxEdgeFromOutside(r,d,2,fBadBoxXZ[i].Z1,xyz); // check front face
312  if(xyz[0] > fBadBoxXZ[i].XY1 && xyz[0] < fBadBoxXZ[i].XY2 && xyz[1] < fEdgeTopYZ && xyz[1] > fEdgeBottomYZ) {
313  tdist = (xyz[2] - vertex.Z())/dir.Z();
314  if(fBadBoxXZ[i].edge == false) {
315  if(tdist < minbaddist && tdist > 0) minbaddist = tdist; // distance to a non-edge bad region
316  }
317  else if(tdist < minedgedist && tdist > 0) minedgedist = tdist; // distance to an edge
318  }
319 
320  geo::ProjectToBoxEdgeFromOutside(r,d,2,fBadBoxXZ[i].Z2,xyz); // check back face
321  if(xyz[0] > fBadBoxXZ[i].XY1 && xyz[0] < fBadBoxXZ[i].XY2 && xyz[1] < fEdgeTopYZ && xyz[1] > fEdgeBottomYZ) {
322  tdist = (xyz[2] - vertex.Z())/dir.Z();
323  if(fBadBoxXZ[i].edge == false) {
324  if(tdist < minbaddist && tdist > 0) minbaddist = tdist; // distance to a non-edge bad region
325  }
326  else if(tdist < minedgedist && tdist > 0) minedgedist = tdist; // distance to an edge
327  }
328 
329  geo::ProjectToBoxEdgeFromOutside(r,d,0,fBadBoxXZ[i].XY1,xyz); // -x face
330  if(xyz[2] > fBadBoxXZ[i].Z1 && xyz[2] < fBadBoxXZ[i].Z2 && xyz[1] < fEdgeTopYZ && xyz[1] > fEdgeBottomYZ) {
331  tdist = (xyz[2] - vertex.Z())/dir.Z();
332  if(fBadBoxXZ[i].edge == false) {
333  if(tdist < minbaddist && tdist > 0) minbaddist = tdist; // distance to a non-edge bad region
334  }
335  else if(tdist < minedgedist && tdist > 0) minedgedist = tdist; // distance to an edge
336  }
337 
338  geo::ProjectToBoxEdgeFromOutside(r,d,0,fBadBoxXZ[i].XY2,xyz); // +x face
339  if(xyz[2] > fBadBoxXZ[i].Z1 && xyz[2] < fBadBoxXZ[i].Z2 && xyz[1] < fEdgeTopYZ && xyz[1] > fEdgeBottomYZ) {
340  tdist = (xyz[2] - vertex.Z())/dir.Z();
341  if(fBadBoxXZ[i].edge == false) {
342  if(tdist < minbaddist && tdist > 0) minbaddist = tdist; // distance to a non-edge bad region
343  }
344  else if(tdist < minedgedist && tdist > 0) minedgedist = tdist; // distance to an edge
345  }
346  }
347 
348  for(unsigned int i = 0; i < fBadBoxYZ.size(); i++) { // now look in YZ view
349  if(fBadBoxYZ[i].edge == false) continue;
350  geo::ProjectToBoxEdgeFromOutside(r,d,2,fBadBoxYZ[i].Z1,xyz); // check front face
351  if(xyz[0] > fEdgeXMinusXZ && xyz[0] < fEdgeXPlusXZ && xyz[1] > fBadBoxYZ[i].XY1 && xyz[1] < fBadBoxYZ[i].XY2) {
352  tdist = (xyz[2] - vertex.Z())/dir.Z();
353  if(fBadBoxYZ[i].edge == false) {
354  if(tdist < minbaddist && tdist > 0) minbaddist = tdist; // distance to a non-edge bad region
355  }
356  else if(tdist < minedgedist && tdist > 0) minedgedist = tdist; // distance to an edge
357  }
358 
359  geo::ProjectToBoxEdgeFromOutside(r,d,2,fBadBoxYZ[i].Z2,xyz); // check back face
360  if(xyz[0] > fEdgeXMinusXZ && xyz[0] < fEdgeXPlusXZ && xyz[1] > fBadBoxYZ[i].XY1 && xyz[1] < fBadBoxYZ[i].XY2) {
361  tdist = (xyz[2] - vertex.Z())/dir.Z();
362  if(fBadBoxYZ[i].edge == false) {
363  if(tdist < minbaddist && tdist > 0) minbaddist = tdist; // distance to a non-edge bad region
364  }
365  else if(tdist < minedgedist && tdist > 0) minedgedist = tdist; // distance to an edge
366  }
367 
368  geo::ProjectToBoxEdgeFromOutside(r,d,1,fBadBoxYZ[i].XY1,xyz); // bottom face
369  if(xyz[2] > fBadBoxYZ[i].Z1 && xyz[2] < fBadBoxYZ[i].Z2 && xyz[0] > fEdgeXMinusXZ && xyz[0] < fEdgeXPlusXZ) {
370  tdist = (xyz[2] - vertex.Z())/dir.Z();
371  if(fBadBoxYZ[i].edge == false) {
372  if(tdist < minbaddist && tdist > 0) minbaddist = tdist; // distance to a non-edge bad region
373  }
374  else if(tdist < minedgedist && tdist > 0) minedgedist = tdist; // distance to an edge
375  }
376 
377  geo::ProjectToBoxEdgeFromOutside(r,d,1,fBadBoxYZ[i].XY2,xyz); // top face
378  if(xyz[2] > fBadBoxYZ[i].Z1 && xyz[2] < fBadBoxYZ[i].Z2 && xyz[0] > fEdgeXMinusXZ && xyz[0] < fEdgeXPlusXZ) {
379  tdist = (xyz[2] - vertex.Z())/dir.Z();
380  if(fBadBoxYZ[i].edge == false) {
381  if(tdist < minbaddist && tdist > 0) minbaddist = tdist; // distance to a non-edge bad region
382  }
383  else if(tdist < minedgedist && tdist > 0) minedgedist = tdist; // distance to an edge
384  }
385  }
386 
387  distToEdge = minedgedist;
388  if(!IsPointLive(vertex)) { // are we starting in a live region?
389  distToNextBad = 0;
390  minbaddist = 0;
391  }
392  else distToNextBad = std::min(minbaddist,minedgedist);
393 
394  if(distToNextBad >= distToEdge) { // the closest bad region is an edge, and we didn't start in a bad region - we're done
395  distDead = 0;
396  return true;
397  }
398 
399  else { // we passed through non-edge bad regions, or started in a bad region - need to figure out how much, so we'll step through region - slower when we have to do this
400  bool liveYet = true; // have we encountered any live material yet
401  if(!IsPointLive(vertex)) liveYet = false; // starting in a bad region
402  int info; // for IsPointLiveInfo return
403  TVector3 position = vertex + minbaddist*dir;
404  double liveDist = 0;
405  double ss = 0.1; // step size (cm)
406  unsigned int maxstep = 100000; // keep it at 100 m
407  for(unsigned int istep = 0; istep < maxstep; istep++) { // should never get to (100 m), but prevent infinite loop if error
408  position = position + ss*dir; // step by step size
409  if(!IsPointLiveInfo(position,info)) {
410  if(liveYet) { // once we've seen live material, the next bad edge region means we're done
411  if(info==5||info==4||info==2||info==0) {
412  if(distToNextBad == 0) distToEdge = minbaddist+ss*istep; // if we starting in bad material, we'll have a messed up distToEdge calculation, since it will just be to the edge of the bad box we started in. We had to step through anyway, so count proper distance at the same time.
413  return true; // we've arrived at an edge, we're done
414  }
415  if(info==3) distDead+=ss; // dead in both views;
416  if(info==1) distDead+=ss/2.0; // dead in only one view, count for half
417  }
418  else { // we didn't start in a live region... need to get to one before we start looking for edge
419  if(info==5||info==4||info==3) distDead+=ss;
420  else distDead+=ss/2.0;
421  }
422  }
423  else { // we found live material
424  liveYet = true;
425  liveDist+=ss;
426  continue;
427  }
428  }
429  } // end determine distDead else
430  distToNextBad = -5.;
431  distToEdge = -5.;
432  return false; // how did we get here? Must have stepped all the way through. Probably means we started in a bad region, then at least one view was bad along the whole path. This is a failure mode, return false.
433 
434  }
435 
436  //...................................................................
439  double distToEdge, distToNextBad, distDead;
440  if(ProjectedDistance(vertex, dir, distToEdge, distToNextBad, distDead)) {
441  TVector3 r2 = vertex + distToEdge*dir;
442  std::vector< OfflineChan > xhits, yhits;
443  geom->CountCellsOnLineFast(vertex, r2, xhits, yhits);
444  return(xhits.size() + yhits.size());
445  }
446  else return -5;
447  }
448 
449  //...................................................................
452 
453  double distToEdge, distToNextBad, distDead;
454  if(!IsPointLive(vertex)) return -5;
455  if(ProjectedDistance(vertex, dir, distToEdge, distToNextBad, distDead)) {
456  TVector3 r2 = vertex + distToNextBad*dir;
457  std::vector< OfflineChan > xhits, yhits;
458  geom->CountCellsOnLineFast(vertex, r2, xhits, yhits);
459  return(xhits.size() + yhits.size());
460  }
461  else return -5;
462  }
463 
464  //...................................................................
467 
468  double distToEdge, distToNextBad, distDead;
469  if(ProjectedDistance(vertex, dir, distToEdge, distToNextBad, distDead)) {
470  TVector3 r2 = vertex + distToEdge*dir;
471  std::vector< OfflineChan > xhits, yhits;
472  geom->CountCellsOnLineFast(vertex, r2, xhits, yhits);
473  int nlivecells = (xhits.size() + yhits.size()) * (1- (distDead/distToEdge) ) + 1; // scale live cells by live distance, round up
474  return(nlivecells);
475  }
476  else return -5;
477  }
478 
479  //...................................................................
481 
482  return std::min(DistToBack(vertex),DistToFront(vertex));
483  }
484 
485  //...................................................................
487 
488  if(!IsPointLive(vertex)) return vertex.Z()-InstrumentedDetFront(); // already in a dead region
489  double fDistFront = vertex.Z()-std::min(fEdgeFrontXZ,fEdgeFrontYZ);
490  for(unsigned int i = 0; i < fBadBoxXZ.size(); i++) { // check XZ view bad regions
491  if(fBadBoxXZ[i].XY1 <= vertex.X() && fBadBoxXZ[i].XY2 >= vertex.X()) {
492  if(vertex.Z() >= fBadBoxXZ[i].Z2) fDistFront = std::min(fDistFront,vertex.Z()-fBadBoxXZ[i].Z2);
493  }
494  }
495  for(unsigned int i = 0; i < fBadBoxYZ.size(); i++) { // check YZ view bad regions
496  if(fBadBoxYZ[i].XY1 <= vertex.Y() && fBadBoxYZ[i].XY2 >= vertex.Y()) {
497  if(vertex.Z() >= fBadBoxYZ[i].Z2) fDistFront = std::min(fDistFront,vertex.Z()-fBadBoxYZ[i].Z2);
498  }
499  }
500 
501  return fDistFront;
502  }
503 
504  //...................................................................
505  double LiveGeometry::DistToBack(TVector3 vertex) {
506 
507  if(!IsPointLive(vertex)) return InstrumentedDetBack()-vertex.Z(); // already in a dead region
508  double fDistBack = 0;
509  fDistBack = std::max(fEdgeBackXZ,fEdgeBackYZ)-vertex.Z();
510 
511  for(unsigned int i = 0; i < fBadBoxXZ.size(); i++) { // check XZ view bad regions
512  if(fBadBoxXZ[i].XY1 <= vertex.X() && fBadBoxXZ[i].XY2 >= vertex.X()) {
513  if(vertex.Z() <= fBadBoxXZ[i].Z1) fDistBack = std::min(fDistBack,fBadBoxXZ[i].Z1-vertex.Z());
514  }
515  }
516 
517  for(unsigned int i = 0; i < fBadBoxYZ.size(); i++) { // check YZ view bad regions
518  if(fBadBoxYZ[i].XY1 <= vertex.Y() && fBadBoxYZ[i].XY2 >= vertex.Y()) {
519  if(vertex.Z() <= fBadBoxYZ[i].Z1) fDistBack = std::min(fDistBack,fBadBoxYZ[i].Z1-vertex.Z());
520  }
521  }
522  return fDistBack;
523  }
524 
525  //...................................................................
527 
528  return std::min(DistToTop(vertex),DistToBottom(vertex));
529  }
530 
531  //...................................................................
532  double LiveGeometry::DistToTop(TVector3 vertex) {
533 
534  double fDistTop = fEdgeTopYZ-vertex.Y();
535  if(!IsPointLive(vertex)) return fDistTop; // already in a dead region
536 
537  for(unsigned int i = 0; i < fBadBoxYZ.size(); i++) { // check YZ view bad regions
538  if(fBadBoxYZ[i].Z1 <= vertex.Z() && fBadBoxYZ[i].Z2 >= vertex.Z()) {
539  if(vertex.Y() <= fBadBoxYZ[i].XY1) fDistTop = std::min(fDistTop,fBadBoxYZ[i].XY1-vertex.Y());
540  }
541  }
542 
543  return fDistTop;
544  }
545 
546  //...................................................................
548 
549  double fDistBottom = vertex.Y()-fEdgeBottomYZ;
550  if(!IsPointLive(vertex)) return fDistBottom; // already in a dead region
551 
552  for(unsigned int i = 0; i < fBadBoxYZ.size(); i++) { // check YZ view bad regions
553  if(fBadBoxYZ[i].Z1 <= vertex.Z() && fBadBoxYZ[i].Z2 >= vertex.Z()) {
554  if(vertex.Y() >= fBadBoxYZ[i].XY1) fDistBottom = std::min(fDistBottom,vertex.Y()-fBadBoxYZ[i].XY1);
555  }
556  }
557 
558  return fDistBottom;
559  }
560 
561  //...................................................................
563 
564  return std::min(DistToEdgeXPlus(vertex),DistToEdgeXMinus(vertex));
565  }
566 
567  //...................................................................
569 
570  return std::min(DistToEdgeX(vertex),DistToEdgeY(vertex));
571  }
572 
573 
574  //...................................................................
576 
577  double fDistXPlus = fEdgeXPlusXZ-vertex.X(); // basic projected dists
578  if(!IsPointLive(vertex)) return fDistXPlus; // already in a dead region
579 
580  for(unsigned int i = 0; i < fBadBoxXZ.size(); i++) { // check XZ view bad regions
581  if(fBadBoxXZ[i].Z1 <= vertex.Z() && fBadBoxXZ[i].Z2 >= vertex.Z()) {
582  if(vertex.X() <= fBadBoxXZ[i].XY1) fDistXPlus = std::min(fDistXPlus,fBadBoxXZ[i].XY1-vertex.X());
583  }
584  }
585  return fDistXPlus;
586  }
587 
588  //...................................................................
590 
591  double fDistXMinus = vertex.X()-fEdgeXMinusXZ;
592  if(!IsPointLive(vertex)) return fDistXMinus; // already in a dead region
593 
594  for(unsigned int i = 0; i < fBadBoxXZ.size(); i++) { // check XZ view bad regions
595  if(fBadBoxXZ[i].Z1 <= vertex.Z() && fBadBoxXZ[i].Z2 >= vertex.Z()) {
596  if(vertex.X() >= fBadBoxXZ[i].XY1) fDistXMinus = std::min(fDistXMinus,vertex.X()-fBadBoxXZ[i].XY1);
597  }
598  }
599  return fDistXMinus;
600  }
601 
602  //...................................................................
604 
605  if(DistToEdgeXPlus(vertex) >= 0) return (int)(DistToEdgeXPlus(vertex)*float(fNumCells+1)/(fEdgeXPlusXZ-fEdgeXMinusXZ));
606  else return -999;
607  }
608 
609  //...................................................................
611 
612  if(DistToEdgeXMinus(vertex) >= 0) return (int)(DistToEdgeXMinus(vertex)*float(fNumCells+1)/(fEdgeXPlusXZ-fEdgeXMinusXZ));
613  else return -999;
614  }
615 
616  //...................................................................
618 
619  return std::min(CellsToEdgeXMinus(vertex),CellsToEdgeXPlus(vertex));
620  }
621 
622  //...................................................................
624 
625  if(DistToTop(vertex) >= 0) return (int)(DistToTop(vertex)*float(fNumCells+1)/(fEdgeTopYZ-fEdgeBottomYZ));
626  else return -999;
627  }
628 
629  //...................................................................
631 
632  if(DistToBottom(vertex) >= 0) return (int)(DistToBottom(vertex)*float(fNumCells+1)/(fEdgeTopYZ-fEdgeBottomYZ));
633  else return -999;
634  }
635 
636  //...................................................................
638 
639  return std::min(CellsToTop(vertex),CellsToBottom(vertex));
640  }
641 
642  //...................................................................
644 
645  return std::min(CellsToEdgeY(vertex),CellsToEdgeX(vertex));
646  }
647 
648  //...................................................................
650 
651  if(DistToFront(vertex) >= 0) return (int)(DistToFront(vertex)*float(fNumDiblocks*64)/(fEdgeBackXZ-fEdgeFrontXZ));
652  return -999;
653  }
654 
655  //...................................................................
657 
658  if(DistToBack(vertex) >= 0) return (int)(DistToBack(vertex)*float(fNumDiblocks*64)/(fEdgeBackXZ-fEdgeFrontXZ));
659  else return -999;
660  }
661 
662  //...................................................................
664 
665  return std::min(PlanesToFront(vertex),PlanesToBack(vertex));
666  }
667 
668  //...................................................................
669  double LiveGeometry::LiveDetectorVolume() { // in ktons
670 
672  if(!fEvFlag) GetEvInfo(); // setup hasn't been called yet!
673  return (double)(RunHist->NActiveChannels())*fMassKtons/float(fNumChannels);
674  }
675 
676  //...................................................................
678  {
679  if(!fEvFlag) GetEvInfo(); // setup hasn't been called yet!
680  return fBadBoxXZ.size();
681  }
682 
683  //...................................................................
685  {
686  if(!fEvFlag) GetEvInfo(); // setup hasn't been called yet!
687  return fBadBoxYZ.size();
688  }
689 
690  //...................................................................
691 
692  double LiveGeometry::GetBadBoxCorner(bool view, int coord, int i) {
693  if(!fEvFlag) GetEvInfo(); // setup hasn't been called yet!
694  if (view==false) { // XZ view
695  if(coord==1) return fBadBoxXZ[i].XY1;
696  else if(coord==2) return fBadBoxXZ[i].XY2;
697  else if(coord==3) return fBadBoxXZ[i].Z1;
698  else if(coord==4) return fBadBoxXZ[i].Z2;
699  else return -999;
700  }
701 
702  if (view==true) { // YZ view
703  if(coord==1) return fBadBoxYZ[i].XY1;
704  else if(coord==2) return fBadBoxYZ[i].XY2;
705  else if(coord==3) return fBadBoxYZ[i].Z1;
706  else if(coord==4) return fBadBoxYZ[i].Z2;
707  else return -999;
708  }
709  return -999; // should never get here
710  } // end GetBadBoxCorner
711 
712  // Private Functions
713 
714  //...................................................................
716  {
717 
720 
721  for(int i = 0; i<12; i++) { // initialize DCMs
722  for(int j = 0; j<14; j++) {
723  fDCMs[i][j] = 1;
724  }
725  }
726 
727  if(geom->DetId() == novadaq::cnv::kFARDET) {
728  fNumDiblocks = 14;
729  fNumDCMs = 12;
730  fNumCells = 383;
731  fNumChannels = 344064;
732  fMassKtons = 14.0;
733  if(fParams.Verbose()) std::cout << "LiveGeometry: Detector set as Far Detector" << std::endl;
734  }
735  else if (geom->DetId() == novadaq::cnv::kNEARDET) {
736  fNumDiblocks = 3; // ignore muon catcher for now
737  fNumDCMs = 4;
738  fNumCells = 95;
739  fNumChannels = 18432; // no muon catcher
740  fMassKtons = 0.3;
741  if(fParams.Verbose()) std::cout << "LiveGeometry: Detector set as Near Detector" << std::endl;
742  }
743  else if(geom->DetId() == novadaq::cnv::kTESTBEAM) {
744  fNumDiblocks = 1;
745  fNumDCMs = 3;
746  fNumCells = 63;
747  fNumChannels = 4032;
748  fMassKtons = 0.03; // this should be fixed once detector is filled
749  if(fParams.Verbose()) std::cout << "LiveGeometry: Detector set as Test Beam" << std::endl;
750  }
751 
752  // START GETTING CORRECT DETECTOR BOUNDARIES
753 
754  // Get real boundaries of 2-D detectors - XZ det and YZ det. Get from cell positions - cell positions may change, but cell locations relative to other cells will not
755  // Can't use fHalfWidth etc as these are 3D edges - horizontal planes may stick out in X compared to the edge of the vertical planes, and we don't want to
756  // have that ~12 cm of extra space with no cells in it when thinking about what happens in the XZ view, for instance.
757 
759 
760  if(fParams.Verbose()) {
761  std::cout << "LiveGeometry: Detector Edges: XZ front/back/-X/+X " << fEdgeFrontXZ << " " << fEdgeBackXZ << " " << fEdgeXMinusXZ << " " << fEdgeXPlusXZ << std::endl;
762  std::cout << "LiveGeometry: Detector Edges: YZ front/back/top/bottom " << fEdgeFrontYZ << " " << fEdgeBackYZ << " " << fEdgeTopYZ << " " << fEdgeBottomYZ << std::endl;
763  }
764 
765  // END GETTING CORRECT DETECTOR BOUNDARIES
766 
767  GetRHInfo(); // get bad boxes based on RH (what is 'uninstrumented' or fDCMs=2)
768  // if(fParams.CheckBadChannels()) GetBCInfo(); // get bad boxes based on BadChannels (What is 'bad' or fDCMs=0)
769  // now we know which DCMs are 'good' (1), 'bad' (0), or 'uninstrumented' (2). Now we need to make bad boxes in coordinate space.
770 
771  // FillBadBoxes();
772 
773  fSetup = true;
774 
775  }
776 
777  //...................................................................
779  {
780  // Check RunHistory to see if the DCM is instrumented. If not, set array to 2
783  int numinst[14][12] = {{0}}; // one entry for each DCM in full det; how many instrumented FEBs are in that DCM? start at zero
784  if(RunHist->NDiBlocks()==0) {
785  for(int i = 0; i<12; i++) { // initialize DCMs
786  for(int j = 0; j<14; j++) {
787  fDCMs[i][j] = 2;
788  }
789  }
790  return; // No RH info for this run! Ignore RH.
791  }
792  for(int idib = 0; idib < RunHist->NDiBlocks(); idib++) {// NOTE we loop over RH array, NOT real diblocks. These can differ; RH only keeps info with data in its index.
793  // so, if diblocks 2, 4, 5, and 6 had data but nothing else, then this will be an array from 0 to 4, with RH->NDiBlocks = 4.
794  // when we do something like take diblock 1 out of partition 1 these will not match up, so treat generally
795  nova::dbi::RunHistory::DiBlock dib = RunHist->GetDiBlock(idib, false); // false flag prevents loading in unnecessary pixel masks which is slooooooow
796  for(int idcm = 0; idcm < int(dib.dcm.size()); idcm++) { // loop over all DCMs in diblocks
797  nova::dbi::RunHistory::DCM dcm = dib.dcm[idcm];
798  if(fParams.Verbose()) std::cout << "LiveGeometry: Got dcm: " << idcm << " on diblock " << dib.num << " (which is number " << idib+1 << " of " << RunHist->NDiBlocks() << ") with " << dcm.nInstrumentedFEBs << " instrumented " << std::endl;
799  numinst[dib.num-1][idcm] = dcm.nInstrumentedFEBs;
800  }
801  }
802 
803  for(int idib = 0; idib < fNumDiblocks; idib++) { // now loop over REAL diblocks to check so things match up later
804  for(int idcm = 0; idcm < fNumDCMs; idcm++) {
805  if(geom->DetId() == novadaq::cnv::kNEARDET && (idcm == 0 || idcm == 3)) {if(numinst[idib][idcm] < 26) fDCMs[idcm][idib] = 2;}
806  else if(numinst[idib][idcm] < 56) fDCMs[idcm][idib] = 2;
807  }
808  }
809  }
810 
811  //...................................................................
813  {
814  // Now check only instrumented DCMs for bad channels. Use criteria:
815  // require 16 good pixels for a good FEB, require 56 good FEBs for a good DCM
816  // add criteria: no more than 200 total bad pixels out of 2048, about 10%)
817  // Mapping notes:
818  // DCMs 1-6 verticals, DCMs 7-12 horizontal
819  // plane 0-63 diblock 1, 64-127 diblock 2 ... 832-895 diblock 14
820  // even number planes (0,2,4...62) horizontals (YZ), odd planes (1,3,5..63) verticals (XZ)
821  // for vertical planes: cell 0-63 -> DCM 1, 64-127 -> DCM 2 ... 320-383 -> DCM 6
822  // for horizontal planes: cell 0-63 -> DCM 12, 64-127 -> DCM 11 ... 320-383 -> DCM 7
823  // for DCM 12: cells 0-31 is a FEB; plane 0 = FEB 32... plane 62 = FEB 63; cells 64-95 for DCM 11 etc
824  // for DCM 12: cells 32-63 is a FEB; plane 0 = FEB 0... plane 62 = FEB 31; cells 96-127 for DCM 11 etc
825  // for DCM 1: cells 0-31 is a FEB; plane 1 = FEB 0... plane 63 = FEB 31; cells 64-95 for DCM 2 etc
826  // for DCM 1: cells 32-63 is a FEB; plane 1 = FEB 32... plane 63 = FEB 63; cells 96-127 for DCM 2 etc
827  // bad channel map should show uninstrumented cells as bad, and we will track difference
830  fNBadChan = 0;
831  fNDropouts = 0;
832 
833  if(geom->DetId() == novadaq::cnv::kTESTBEAM){
834  std::cout<<"GetBCInfo() has incorrect mapping for TestBeam. Setting NBadChan and NDropouts to zero!"<<std::endl;
835  return;
836  }
837 
838  for(int i = 0; i<fNumDCMs; i++) {
839  for(int j = 0; j<fNumDiblocks; j++) {
840  if(fDCMs[i][j]==0) fDCMs[i][j]=1; // reset what is 'bad' every event
841  }
842  }
843 
844  for(int idib = 0; idib < fNumDiblocks; idib++) {
845 
846  int jplane, jcell1, jcell2; // actual plane and cell number
847  int nfeb1, nfeb2; // number of dead channels in FEB
848  int nbadDCM[12] = {0}; // number of bad FEBs in DCM
849  int nbadcount[12] = {0}; // number of bad pixels in DCM
850 
851  for(int idcm = 0; idcm < fNumDCMs/2; idcm++) { // verticals first (DCMs 1-6, indices 0-5)
852  if(fDCMs[idcm][idib]==1) { // check if already tagged as uninstrumented
853  for(int iplane = 0; iplane < 32; iplane++) { // each plane is 2 FEBs, for 64 FEBs per DCM
854  nfeb1 = 0; // reset these values for new FEBs
855  nfeb2 = 0;
856  jplane = (idib*64) + (2*iplane) + 1; // odd numbered planes
857  for(int icell = 0; icell < 32; icell++) { // this loops over the pixels of those 2 FEBs
858  jcell1 = icell + 64*idcm;
859  jcell2 = icell + 32 + 64*idcm;
860  if(BadList->IsBad(jplane,jcell1)) nfeb1++; // check if this pixel on FEB 1 is bad
861  if(BadList->IsBad(jplane,jcell2)) nfeb2++; // check if this pixel on FEB 2 is bad
862  } // end verticals pixels loop
863  if(nfeb1 > 6) nbadDCM[idcm]++; // if the FEB had too many bad channels, count it as a bad FEB
864  if(nfeb2 > 6) nbadDCM[idcm]++;
865  nbadcount[idcm]+=(nfeb1+nfeb2); // keep track of total bad pixels
866  } // end verticals planes loop
867  } // end if instrumented
868  } // end verticals dcm loop
869 
870  for(int idcm = fNumDCMs/2; idcm < fNumDCMs; idcm++) { // now horizontals (DCMs 7-12, indices 6-11)
871  if(fDCMs[idcm][idib]==1) { // check if already tagged as uninstrumented
872  for(int iplane = 0; iplane < 32; iplane++) { // each plane is 2 FEBs, for 64 FEBs per DCM
873  nfeb1 = 0; // reset these values for new FEBs
874  nfeb2 = 0;
875  jplane = (idib*64) + (2*iplane); // even numbered planes
876  for(int icell = 0; icell < 32; icell++) { // this loops over the pixels of those 2 FEBs
877  jcell1 = icell + 64*(11-idcm); // arranged opposite of verticals
878  jcell2 = icell + 32 + 64*(11-idcm);
879  if(BadList->IsBad(jplane,jcell1)) nfeb1++; // check if this pixel on FEB 1 is bad
880  if(BadList->IsBad(jplane,jcell2)) nfeb2++; // check if this pixel on FEB 2 is bad
881  } // end horizontals pixels loop
882  if(nfeb1 > 6) nbadDCM[idcm]++; // if the FEB had too many bad channels, count it as a bad FEB
883  if(nfeb2 > 6) nbadDCM[idcm]++;
884  nbadcount[idcm]+=(nfeb1+nfeb2); // keep track of total bad pixels
885  } // end horizontals planes loop
886  } // end if instrumented
887  } // end horizontals dcm loop
888 
889  for(int idcm = 0; idcm < fNumDCMs; idcm++) { // we looked at all pixels in this diblock, check results then go to next one
890  if(nbadDCM[idcm]>=63) fNDropouts++;
891  if(geom->DetId() == novadaq::cnv::kNEARDET) { // for ND need different criteria for DCMs that are only half instrumented on purpose (1 and 4 of every diblock)
892  if(idcm==0||idcm==3) {
893  if(nbadDCM[idcm] > 4 && fDCMs[idcm][idib]==1) {
894  fDCMs[idcm][idib] = 0; // half as strict criteria for half instrumented DCMs
895  fNBadChan++;
896  }
897  }
898  else if((nbadDCM[idcm] > 8) && fDCMs[idcm][idib]==1) {
899  fDCMs[idcm][idib] = 0; // full DCMs same as FD
900  fNBadChan++;
901  }
902  }
903  else if((nbadDCM[idcm] > 8) && fDCMs[idcm][idib]==1) {
904  fDCMs[idcm][idib] = 0; // set DCM to bad if it fails criteria and isn't already uninstrumented
905  fNBadChan++;
906  }
907  if(fParams.Verbose()) { // optional debugging information
908  std::cout << "LiveGeometry: fDCMs: dib " << idib+1 << " DCM " << idcm+1 << " value " << fDCMs[idcm][idib] << " # bc: " << nbadcount[idcm] << " # bad FEBs " << nbadDCM[idcm] << " count of bad DCMs: " << fNBadChan << std::endl;
909  }
910 
911  } // end checking DCMs loop
912 
913  } // end diblock loop
914  }
915 
916  //...................................................................
918  {
920  double xyz1[3], xyz2[3];
921 
922  if(!(geom->DetId() == novadaq::cnv::kFARDET ||
923  geom->DetId() == novadaq::cnv::kNEARDET ||
924  geom->DetId() == novadaq::cnv::kTESTBEAM )) {
925  std::cout << "UNKNOWN DETECTOR!!!" << std::endl;
926  fEdgeXPlusXZ = 0;
927  fEdgeXMinusXZ = 0;
928  fEdgeFrontXZ = 0;
929  fEdgeBackXZ = 0;
930  fEdgeTopYZ = 0;
931  fEdgeBottomYZ = 0;
932  fEdgeFrontYZ = 0;
933  fEdgeBackYZ = 0;
934  }
935 
936  int firstXZ = -1; // first XZ plane in diblock (vertical)
937  int lastXZ = -1; // last XZ plane in diblock
938  int firstYZ = -1; // first YZ plane in diblock (horizontal)
939  int lastYZ = -1; // last YZ plane in diblock
940  // TB starts and ends with vertical plane
941  if (geom->DetId() == novadaq::cnv::kTESTBEAM) {
942  firstXZ = 0;
943  lastXZ = 62;
944  firstYZ = 1;
945  lastYZ = 61; // no plane 63 in TB
946  }
947  // ND/FD start with a horizontal plane, end with vertical
948  else {
949  firstXZ = 1;
950  lastXZ = 63;
951  firstYZ = 0;
952  lastYZ = 62;
953  }
954 
955  // start with XZ view, which is easiest since z position does NOT change with cell number; this is odd plane numbers (except for TestBeam)
956  geom->Plane(firstXZ)->Cell(0)->GetCenter(xyz1);
957  fEdgeFrontXZ = xyz1[2]-3.5; // the half width of the cell in Z is ~3.5 cm
958  geom->Plane(lastXZ)->Cell(0)->GetCenter(xyz2);
959  fEdgeXMinusXZ = (xyz1[0]+xyz2[0])/2.0 - 2.3; // Plane 0 is 'higher' than plane 63; this repeats for each diblock exactly. We don't want to worry about this though, so average it. We pretend things are actual boxes with right angles.
960  geom->Plane(firstXZ)->Cell(fNumCells)->GetCenter(xyz1); // cell 383 is most +x, cell 0 is most -X
961  geom->Plane(lastXZ)->Cell(fNumCells)->GetCenter(xyz2);
962  fEdgeXPlusXZ = (xyz1[0]+xyz2[0])/2.0 + 2.3; // 2.28 is approximate half width of cell here; again, round up slightly
963  if (geom->DetId()== novadaq::cnv::kTESTBEAM) {
964  geom->Plane(62)->Cell(0)->GetCenter(xyz1); // last XZ plane for TB is 62
965  }
966  else {
967  geom->Plane((fNumDiblocks*64)-1)->Cell(0)->GetCenter(xyz1); // any cell in last plane will do for finding the back
968  }
969  fEdgeBackXZ = xyz1[2]+3.5;
970 
971  // now YZ view, which has Y AND Z changing with cell number for the same plane. The Z from cell 0 to cell 383 differs by about 2.5 cm.
972  // I'll average the Z values over a DCM to find the DCM edges. I never want a gap between a DCM edge and the detector edge, since that could unwittingly be considered 'good' even if the closest DCM is 'bad'. So, I'll use the DCMs that stick out the least to define the edge. This makes us 'lose' about 2 cm from parts of the edge at each end, but it's safest.
973 
974  geom->Plane(firstYZ)->Cell(0)->GetCenter(xyz1); // front bottom most cell
975  geom->Plane(firstYZ)->Cell(63)->GetCenter(xyz2); // front top cell of frontmost DCM
976  fEdgeFrontYZ = (xyz1[2]+xyz2[2])/2.0 - 3.5; // average the 2 to get the 'front'
977  geom->Plane(lastYZ)->Cell(0)->GetCenter(xyz2); // back most bottom cell of frontmost bottom DCM
978  fEdgeBottomYZ = (xyz1[1]+xyz2[1])/2.0 - 2.3; // reuse same point 1 for bottom calculation
979  geom->Plane(firstYZ)->Cell(fNumCells)->GetCenter(xyz1); // topmost frontmost cell, in topmost frontmost DCM
980  geom->Plane(lastYZ)->Cell(fNumCells)->GetCenter(xyz2); // topmost farther back cell in topmost frontmost DCM
981  fEdgeTopYZ = (xyz1[1]+xyz2[1])/2.0 + 2.3; // top edge
982  if (geom->DetId() == novadaq::cnv::kTESTBEAM) {
983  geom->Plane(61)->Cell(fNumCells-63)->GetCenter(xyz1); // frontmost topmost cell in backmost topmost DCM
984  geom->Plane(61)->Cell(fNumCells)->GetCenter(xyz2); // backmost topmost cell, in backmost topmost DCM
985  }
986  else {
987  geom->Plane((fNumDiblocks*64)-2)->Cell(fNumCells-63)->GetCenter(xyz1); // frontmost topmost cell in backmost topmost DCM
988  geom->Plane((fNumDiblocks*64)-2)->Cell(fNumCells)->GetCenter(xyz2); // backmost topmost cell, in backmost topmost DCM
989  }
990  fEdgeBackYZ = (xyz1[2]+xyz2[2])/2.0 + 3.5; // back edge
991 
992  // Muon catcher. TODO: Make these more precise
993 
994  // Based on the gdml
995  //fMCEdge = 199.65;
996  //fMCTop = 75.3;
997  //fMCFront = 1277.65;
998  //fMCBack = 1577.54;
999 
1000  // based on the event display
1001  fMCEdge = 190;
1002  fMCTop = 67.5;
1003  fMCFront = 1278.4;
1004  fMCBack = 1585;
1005  }
1006 
1007  //...................................................................
1009  {
1010  // First, cluster bad DCMs - if half of the detector is bad, don't want 72 bad boxes, just want one. Clustering is very dumb, but should get most of it.
1011 
1012  fBadBoxXZ.clear(); // get rid of old bad boxes - we're redoing them
1013  fBadBoxYZ.clear();
1014 
1015  int DCM[12][14] = {{0}}; // a third array that merges the information of the other two. This way fDCMs is never altered.
1016 
1017  for(int i = 0; i<fNumDCMs; i++) {
1018  for(int j = 0; j<fNumDiblocks; j++) {
1019  if(fDCMs[i][j]==0||fDCMs[i][j]==2) DCM[i][j]=0;
1020  else DCM[i][j]=1;
1021  }
1022  }
1023 
1025 
1026  if(geom->DetId() == novadaq::cnv::kTESTBEAM){
1027  std::cout<<"FillBadBoxes() not valid for TestBeam. Not filling any."<<std::endl;
1028  return;
1029  }
1030 
1031  double xyz1[3], xyz2[3];
1032  // do some clustering. Only cluster uninstrumented; DCMs that fail event check always left as single boxes. This helps distinguish
1033  bool xzview[fNumDiblocks+1],yzview[fNumDiblocks+1]; // look at each diblock/view for if all DCMs are uninstrumented. If they all are, this is true. Size 15 (ndib+1) just to prevent crash later.
1034  for(int i = 0; i<fNumDiblocks+1; i++) { // initialize to true
1035  xzview[i] = true;
1036  yzview[i] = true;
1037  }
1038 
1039  for(int idib = 0; idib < fNumDiblocks; idib++) {
1040  for(int idcm = 0; idcm < fNumDCMs/2; idcm++) {
1041  if(fDCMs[idcm][idib]!=2) xzview[idib] = false; // if even a single DCM is not bad, the whole thing is false
1042  }
1043  if(xzview[idib]==true) for(int idcm = 0; idcm < fNumDCMs/2; idcm++) {
1044  DCM[idcm][idib] = 1; // we're handling these bad boxes now, don't want to again later
1045  }
1046  for(int idcm = fNumDCMs/2; idcm < fNumDCMs; idcm++) {
1047  if(fDCMs[idcm][idib]!=2) yzview[idib] = false; // if even a single DCM is not bad, the whole thing is false
1048  }
1049  if(yzview[idib]==true) for(int idcm = fNumDCMs/2; idcm < fNumDCMs; idcm++) {
1050  DCM[idcm][idib] = 1; // we're handling these bad boxes now, don't want to again later
1051  }
1052  }
1053 
1054  LiveGeometry::BadBox tBox;
1055 
1056  int cstart=-5; // front of clusters
1057  for(int idib = 0; idib<fNumDiblocks+1; idib++) { // lets make the bad boxes, starting with XZ
1058 
1059  if(xzview[idib]==false||idib==fNumDiblocks) { // nope this diblock isn't all bad, or at end
1060  if(cstart==-5) continue; // there wasn't a cluster ongoing, we'll let the individual DCM bad boxes sort anything else out later
1061  else if(cstart>=0) { // there was a cluster that we have now just ended. Make the bad box!
1062 
1063  //verticals = XZ
1064  geom->Plane((cstart*64)+1)->Cell(fNumCells)->GetCenter(xyz1); // Z is static with cell number in this view; any cell at front face works
1065  tBox.Z1 = xyz1[2]-3.5; // set front (corner 1)
1066  tBox.XY2 = xyz1[0]+2.3; // set +X (corner 2); first plane in DCM, highest cell number has largest +X
1067  geom->Plane((idib*64)-1)->Cell(0)->GetCenter(xyz1); // back-most (of all detector) cell, -X most (of DCM, repeats with DCM) cell;
1068  tBox.Z2 = xyz1[2]+3.5; // set back (corner 2)
1069  tBox.XY1 = xyz1[0]-2.3; // set -X (corner 1)
1070  tBox.type = 0;
1071  tBox.edge = true;
1072  fBadBoxXZ.push_back(tBox);
1073  cstart=-5;
1074  continue;
1075  }
1076  }
1077 
1078  if(xzview[idib]==true) { // yes this diblock is completely uninstrumented in this view
1079  if(cstart==-5) { // no we haven't started any kind of cluster, so start it
1080  cstart = idib;
1081  continue;
1082  }
1083  else continue; // there is already a cluster ongoing, so just keep going
1084  }
1085  }
1086 
1087  cstart = -5;
1088  for(int idib = 0; idib<fNumDiblocks+1; idib++) { // lets make the bad boxes, now YZ
1089 
1090  if(yzview[idib]==false||idib==fNumDiblocks) { // nope this diblock isn't all bad
1091  if(cstart==-5) continue; // there wasn't a cluster ongoing, we'll let the individual DCM bad boxes sort anything else out later
1092  else if(cstart>=0) { // there was a cluster that we have now just ended. Make the bad box!
1093 
1094  LiveGeometry::BadBox tBox;
1095  //horizontals = YZ
1096  geom->Plane(cstart*64)->Cell(fNumCells)->GetCenter(xyz1); // front topmost most cell of uninstrumented region in YZ; most-front cell
1097  tBox.Z1 = xyz1[2]-3.5; // set front (corner 1)
1098  geom->Plane((idib*64)-2)->Cell(fNumCells)->GetCenter(xyz1);
1099  tBox.XY2 = xyz1[1]+2.3; // set top (corner 2)
1100  geom->Plane((idib*64)-2)->Cell(0)->GetCenter(xyz1); // back bottom-most cell in XZ; is furthest back cell
1101  tBox.Z2 = xyz1[2]+3.5; // set back (corner 2)
1102  geom->Plane((fNumDiblocks-1)*64)->Cell(0)->GetCenter(xyz1); // front bottom-most cell in back bottom-most DCM in XZ; is cell
1103  tBox.XY1 = xyz1[1]-2.3; // set bottom (corner 1)
1104  tBox.type = 0;
1105  tBox.edge = true;
1106  fBadBoxYZ.push_back(tBox);
1107  cstart=-5;
1108  continue;
1109  }
1110  }
1111 
1112  if(yzview[idib]==true) { // yes this diblock is completely uninstrumented in this view
1113  if(cstart==-5) { // no we haven't started any kind of cluster, so start it
1114  cstart = idib;
1115  continue;
1116  }
1117  else continue; // there is already a cluster ongoing, so just keep going
1118  }
1119 
1120  }
1121 
1122  // now get bad boxes for any remaining bad or uninstrumented regions, one DCM at a time:
1123 
1124  for(int idib = 0; idib < fNumDiblocks; idib++) { // verticals first = XZ
1125  for(int idcm = 0; idcm < fNumDCMs/2; idcm++) {
1126  if(DCM[idcm][idib]==0) {
1127  if(geom->DetId() == novadaq::cnv::kFARDET) {
1128  geom->Plane(1+64*idib)->Cell(idcm*64)->GetCenter(xyz1); // -x-most front
1129  geom->Plane(63+64*idib)->Cell(idcm*64)->GetCenter(xyz2); // -x-most back
1130  tBox.Z1 = xyz1[2]-3.5; // Z positions are easy in this view; front
1131  tBox.Z2 = xyz2[2]+3.5; // back
1132  tBox.XY1 = (xyz1[0]+xyz2[0])/2.0 - 2.3; // average DCM -x
1133  geom->Plane(1+64*idib)->Cell(63+idcm*64)->GetCenter(xyz1); // +x-most front
1134  geom->Plane(63+64*idib)->Cell(63+idcm*64)->GetCenter(xyz2); // +x-most back
1135  tBox.XY2 = (xyz1[0]+xyz2[0])/2.0 + 2.3; // average DCM +x
1136  if(fDCMs[idcm][idib] == 2) tBox.type = 0; // uninstrumented
1137  else tBox.type = 1; // bad
1138  bool connected = true; // are we connected to an edge?
1139  tBox.edge = false;
1140  if(idib==13 || idib == 0 || idcm == 0 || idcm == 5) tBox.edge = true; // are we immediately connected to the edge?
1141  //current check is just seeing if up, left, down, or right is an unbroken series of bad DCMs to the edge
1142  for(int i = idib+1; i < fNumDiblocks; i++) { // check for unbroken connection to the back
1143  if(DCM[idcm][i]==1) connected = false;
1144  }
1145  if(connected == true) tBox.edge = true;
1146  connected = true;
1147  for(int i = 0; i < idib; i++) { // check for unbroken connection to front
1148  if(DCM[idcm][i]==1) connected = false;
1149  }
1150  if(connected == true) tBox.edge = true;
1151  connected = true;
1152  for(int i = 0; i < idcm; i++) { // check for unbroken connection to bottom
1153  if(DCM[i][idib]==1) connected = false;
1154  }
1155  if(connected == true) tBox.edge = true;
1156  connected = true;
1157  for(int i = idcm+1; i < fNumDCMs/2; i++) { // check for unbroken connection to top
1158  if(DCM[i][idib]==1) connected = false;
1159  }
1160  if(connected == true) tBox.edge = true;
1161  }
1162 
1163  if(geom->DetId() == novadaq::cnv::kNEARDET) {
1164  if(idcm==0) {
1165  geom->Plane(1+64*idib)->Cell(0)->GetCenter(xyz1); // -x-most front
1166  geom->Plane(63+64*idib)->Cell(0)->GetCenter(xyz2); // -x-most back
1167  }
1168  else if(idcm==1) {
1169  geom->Plane(1+64*idib)->Cell(32)->GetCenter(xyz1); // -x-most front
1170  geom->Plane(63+64*idib)->Cell(32)->GetCenter(xyz2); // -x-most back
1171  }
1172  tBox.Z1 = xyz1[2]-3.5; // Z positions are easy in this view; front
1173  tBox.Z2 = xyz2[2]+3.5; // back
1174  tBox.XY1 = (xyz1[0]+xyz2[0])/2.0 - 2.3; // average DCM -x
1175  geom->Plane(1+64*idib)->Cell(idcm*64+31)->GetCenter(xyz1); // +x-most front
1176  geom->Plane(63+64*idib)->Cell(idcm*64+31)->GetCenter(xyz2); // +x-most back
1177  tBox.XY2 = (xyz1[0]+xyz2[0])/2.0 + 2.3; // average DCM +x
1178  if(fDCMs[idcm][idib] == 2) tBox.type = 0; // uninstrumented
1179  else tBox.type = 1; // bad
1180  tBox.edge = true; // every ND DCM is an edge
1181  }
1182  fBadBoxXZ.push_back(tBox);
1183  }
1184  }
1185  }
1186 
1187  for(int idib = 0; idib < fNumDiblocks; idib++) { // horizontals second = YZ
1188  for(int idcm = fNumDCMs/2; idcm < fNumDCMs; idcm++) {
1189  if(DCM[idcm][idib]==0) {
1190  if(geom->DetId() == novadaq::cnv::kFARDET) {
1191  geom->Plane(64*idib)->Cell((11-idcm)*64)->GetCenter(xyz1); // bottom front into 1
1192  geom->Plane(64*idib)->Cell(63+(11-idcm)*64)->GetCenter(xyz2); // top front into 2
1193  tBox.Z1 = (xyz1[2]+xyz2[2])/2.0 - 3.5; // average front
1194  geom->Plane(62+64*idib)->Cell((11-idcm)*64)->GetCenter(xyz2); // bottom back into 2
1195  tBox.XY1 = (xyz1[1]+xyz2[1])/2.0 - 2.3; // average DCM bottom
1196  geom->Plane(62+64*idib)->Cell(63+(11-idcm)*64)->GetCenter(xyz1); // top back into 1
1197  tBox.Z2 = (xyz1[2]+xyz2[2])/2.0 + 3.5; // average back
1198  geom->Plane(64*idib)->Cell(63+(11-idcm)*64)->GetCenter(xyz2); // top front into 2
1199  tBox.XY2 = (xyz1[1]+xyz2[1])/2.0 + 2.3; // average DCM top
1200  if(fDCMs[idcm][idib] == 2) tBox.type = 0; // uninstrumented
1201  else tBox.type = 1; // bad
1202  bool connected = true; // are we connected to an edge?
1203  tBox.edge = false;
1204  if(idib==13 || idib == 0 || idcm == fNumDCMs/2 || idcm == fNumDCMs) tBox.edge = true; // are we immediately connected to the edge?
1205  for(int i = idib+1; i < fNumDiblocks; i++) { // check for unbroken connection to the back
1206  if(DCM[idcm][i]==1) connected = false;
1207  }
1208  if(connected == true) tBox.edge = true;
1209  connected = true;
1210  for(int i = 0; i < idib; i++) { // check for unbroken connection to front
1211  if(DCM[idcm][i]==1) connected = false;
1212  }
1213  if(connected == true) tBox.edge = true;
1214  connected = true;
1215  for(int i = fNumDCMs/2; i < idcm; i++) { // check for unbroken connection to +x side
1216  if(DCM[i][idib]==1) connected = false;
1217  }
1218  if(connected == true) tBox.edge = true;
1219  connected = true;
1220  for(int i = idcm+1; i < fNumDCMs; i++) { // check for unbroken connection -x side
1221  if(DCM[i][idib]==1) connected = false;
1222  }
1223  if(connected == true) tBox.edge = true;
1224  }
1225 
1226  else if(geom->DetId() == novadaq::cnv::kNEARDET) {
1227  if(idcm==2) {
1228  geom->Plane(64*idib)->Cell(32)->GetCenter(xyz1); // bottom front into 1
1229  geom->Plane(64*idib)->Cell(95)->GetCenter(xyz2); // top front into 2
1230  tBox.Z1 = (xyz1[2]+xyz2[2])/2.0 - 3.5; // average front
1231  geom->Plane(62+64*idib)->Cell(32)->GetCenter(xyz2); // bottom back into 2
1232  tBox.XY1 = (xyz1[1]+xyz2[1])/2.0 - 2.3; // average DCM bottom
1233  geom->Plane(62+64*idib)->Cell(95)->GetCenter(xyz1); // top back into 1
1234  tBox.Z2 = (xyz1[2]+xyz2[2])/2.0 + 3.5; // average back
1235  geom->Plane(64*idib)->Cell(95)->GetCenter(xyz2); // top front into 2
1236  tBox.XY2 = (xyz1[1]+xyz2[1])/2.0 + 2.3; // average DCM top
1237  }
1238  if(idcm==3) {
1239  geom->Plane(64*idib)->Cell(0)->GetCenter(xyz1); // bottom front into 1
1240  geom->Plane(64*idib)->Cell(31)->GetCenter(xyz2); // top front into 2
1241  tBox.Z1 = (xyz1[2]+xyz2[2])/2.0 - 3.5; // average front
1242  geom->Plane(62+64*idib)->Cell(0)->GetCenter(xyz2); // bottom back into 2
1243  tBox.XY1 = (xyz1[1]+xyz2[1])/2.0 - 2.3; // average DCM bottom
1244  geom->Plane(62+64*idib)->Cell(31)->GetCenter(xyz1); // top back into 1
1245  tBox.Z2 = (xyz1[2]+xyz2[2])/2.0 + 3.5; // average back
1246  geom->Plane(64*idib)->Cell(31)->GetCenter(xyz2); // top front into 2
1247  tBox.XY2 = (xyz1[1]+xyz2[1])/2.0 + 2.3; // average DCM top
1248  }
1249 
1250  if(fDCMs[idcm][idib] == 2) tBox.type = 0; // uninstrumented
1251  else tBox.type = 1; // bad
1252  tBox.edge = true; // all ND DCMs are edges
1253  }
1254  fBadBoxYZ.push_back(tBox);
1255  }
1256  }
1257  }
1258 
1259  if(geom->DetId() == novadaq::cnv::kNEARDET) { // air gap atop muon catcher always bad; this is just cosmetic, as ND detector is considered only 3 DB and MC functions are separate
1260  tBox.Z1 = 1277.2;
1261  tBox.Z2 = 1585.5;
1262  tBox.XY1 = 70.0;
1263  tBox.XY2 = 199.7;
1264  tBox.edge = true;
1265  fBadBoxYZ.push_back(tBox);
1266  }
1267 
1268  // TODO: change algorithm for check if connected to edge to get L shapes and such
1269 
1270  if(fParams.Verbose()) {
1271  // optional debugging printouts:
1272  for(int i = 0; i < int(fBadBoxXZ.size()); i++) { // check XZ view bad regions // to see the bad boxes for debugging
1273  std::cout << "LiveGeometry: BadBoxXZ: X1 : " << fBadBoxXZ[i].XY1 << " X2: " << fBadBoxXZ[i].XY2 << " Z1: " << fBadBoxXZ[i].Z1 << " Z2: " << fBadBoxXZ[i].Z2 << std::endl;
1274  }
1275 
1276  for(int i = 0; i < int(fBadBoxYZ.size()); i++) { // check XZ view bad regions
1277  std::cout << "LiveGeometry: BadBoxYZ: Y1 : " << fBadBoxYZ[i].XY1 << " Y2: " << fBadBoxYZ[i].XY2 << " Z1: " << fBadBoxYZ[i].Z1 << " Z2: " << fBadBoxYZ[i].Z2 << std::endl;
1278  }
1279 
1280  }
1281  }
1282 
1283  //......................................................................
1285  {
1289 
1290  if(!fEvFlag) GetEvInfo(); // setup hasn't been called yet!
1291 
1292  int dibs[14] = {0}, dibsUD[14] = {0}, dibsRH[14] = {0};
1293 
1294  for(int i=0;i<fNumDiblocks;i++) dibsUD[i]=1;
1295 
1296  if(RunHist->GetDataType() == 2 && RunHist->NDiBlocks()>0) { // if it's data
1297  for(int idib = RunHist->NDiBlocks()-1; idib >= 0; idib--) { // start from back; this is downstream part
1298  nova::dbi::RunHistory::DiBlock dib = RunHist->GetDiBlock(idib, false);
1299  if (RunHist->IsDiBlockFullyInstrumented(dib.num)) dibsRH[dib.num-1] = 1;
1300  }
1301  }
1302  else for(int i=0;i<fNumDiblocks;i++) dibsRH[i]=1;
1303 
1304  for(int i=0;i<fNumDiblocks;i++) {
1305  if(dibsUD[i]==1&&dibsRH[i]==1) dibs[i]=1;
1306  }
1307 
1308  // look from front to back for upstream part
1309  for (int i = 0; i<fNumDiblocks; i++) {
1310  if(dibs[i] == 1 && fFrontUpstream == -5) fFrontUpstream = i; // found the front
1311  if(dibs[i] == 1 && fFrontUpstream != -5) continue; // still in good part, keep going
1312  if(dibs[i] == 0 && fFrontUpstream == -5) continue; // nothing yet, keep looking
1313  if(dibs[i] == 0 && fFrontUpstream != -5) { // found the back
1314  fBackUpstream = i;
1315  break;
1316  }
1317  }
1318  if(fBackUpstream == -5 && fFrontUpstream != -5) { // if the back is diblock 14
1320  }
1321  if(fBackUpstream == -5 && fFrontUpstream == -5) { // no part of the detector is good
1322  fBackUpstream = 0;
1323  fFrontUpstream = 0;
1324  }
1325 
1326  // look from back to front for upstream part
1327  for (int i = fNumDiblocks-1; i>=0; i--) {
1328  if(dibs[i] == 1 && fBackDownstream == -5) fBackDownstream = i+1; // found the back
1329  if(dibs[i] == 1 && fBackDownstream != -5) continue; // still in good part, keep going
1330  if(dibs[i] == 0 && fBackDownstream == -5) continue; // nothing yet, keep looking
1331  if(dibs[i] == 0 && fBackDownstream != -5) { // found the front
1332  fFrontDownstream = i+1;
1333  break;
1334  }
1335  }
1336  if(fFrontDownstream == -5 && fBackDownstream != -5) { // if the front is diblock 0
1337  fFrontDownstream = 0;
1338  }
1339  if(fFrontDownstream == -5 && fBackDownstream == -5) { // no part of the detector is good
1340  fBackDownstream = 0;
1341  fFrontDownstream = 0;
1342  }
1343 
1344  // End search
1345 
1346  if(fBackDownstream==-5 && fBackUpstream==-5) { // there were NO good blocks. Uh oh. Set them all to 0 so we know we looked (-5 is hasn't looked yet)
1347  fBackUpstream = 0;
1348  fFrontUpstream = 0;
1349  fBackDownstream = 0;
1350  fFrontDownstream = 0;
1351  }
1352  else if(fBackDownstream==fBackUpstream && fFrontDownstream==fFrontUpstream) { // this is the one contiguous detector case, both searches find same thing
1353  fBackUpstream = 0; // this is how we'll know there's only one contiguous block
1354  fFrontUpstream = 0;
1355  }
1356  else if(fBackUpstream < fFrontDownstream) {} // this is the two (or more, but then middle chunks ignored) contiguous detectors case
1357  else {
1358  std::cout << "Problem in instrumented detector logic; this should never happen. " << fBackUpstream << " " << fFrontUpstream << " " << fBackDownstream << " " << fFrontDownstream << std::endl;
1359  fBackUpstream = 0;
1360  fFrontUpstream = 0;
1361  fBackDownstream = 0;
1362  fFrontDownstream = 0;
1363  }
1364 
1365  return;
1366  }
1367 
1368  //......................................................................
1370  {
1375  }
1376 
1377  //......................................................................
1379  {
1384  }
1385 
1386  //......................................................................
1388  {
1393  }
1394 
1395  //......................................................................
1396  int LiveGeometry::InstrumentedDetEnds(double &frontDwnstrm, double &backDwnstrm, double &frontUpstrm, double &backUpstrm)
1397  {
1401 
1406 
1407  if(backDwnstrm>0&&backUpstrm>0) return 2; // two contiguous detectors
1408  else if(backDwnstrm>0&&backUpstrm==0) return 1; // one contiguous detector
1409  else return 0; // something is wrong, or no good diblocks
1410  }
1411 
1412  //......................................................................
1414  {
1416  if(geom->DetId() != novadaq::cnv::kNEARDET) {
1417  std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
1418  return false;
1419  }
1420  bool ret = false;
1421  if(vertex.X() > -fMCEdge && vertex.X() < fMCEdge &&
1422  vertex.Y() > -fMCEdge && vertex.Y() < fMCTop &&
1423  vertex.Z() > fMCFront && vertex.Z() < fMCBack) ret = true;
1424  return ret;
1425  }
1426 
1427  //......................................................................
1428  int LiveGeometry::DoWeEnterMC(TVector3 vertex, TVector3 dir)
1429  {
1431  if(geom->DetId() != novadaq::cnv::kNEARDET) {
1432  std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
1433  return 0;
1434  }
1435  if(dir.Z() < 0 && vertex.Z() < fMCFront) return 0; // this track couldn't possible reach the muon catcher
1436  if(dir.Z() > 0 && vertex.Z() > fMCFront) return 0; // this track couldn't possible enter the main ND
1437  double xyz[3], dxyz[3];
1438  xyz[0] = vertex.X();
1439  xyz[1] = vertex.Y();
1440  xyz[2] = vertex.Z();
1441  dxyz[0] = dir.X();
1442  dxyz[1] = dir.Y();
1443  dxyz[2] = dir.Z();
1445  if(InMuonCatcher(vertex)) return 1; // fully contained in muon catcher?
1446  else return 0; // nope, doesn't intersect the muon catcher
1447  }
1448  double ypos = YPositionAtTransition(vertex,dir);
1449  double xpos = XPositionAtTransition(vertex,dir);
1450  if(ypos > fMCTop && ypos < fEdgeTopYZ && xpos > fEdgeXMinusXZ && xpos < fEdgeXPlusXZ) return 2; // intersects but goes through the air gap
1451  else return 1; // seems like it properly enters the MC
1452  }
1453 
1454  //......................................................................
1456  {
1458  if(geom->DetId() != novadaq::cnv::kNEARDET) {
1459  std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
1460  return 0;
1461  }
1462 
1463  double distX, distY, distZ;
1464 
1465  if(InMuonCatcher(vertex)) { // if vertex is in the muon catcher
1466 
1467  if(dir.Z() > 0) distZ = (fMCBack-vertex.Z())/dir.Z(); // going +z project to back
1468  else if(dir.Z() < 0) distZ = (fMCFront-vertex.Z())/dir.Z(); // going -z project to front
1469  else { // dz is zero. This means something is wrong with the track
1470  return 0;
1471  }
1472 
1473  if(dir.Y() > 0) distY = (fMCTop-vertex.Y())/dir.Y(); // going +y project to top
1474  else if(dir.Y() < 0) distY = (-fMCEdge-vertex.Y())/dir.Y(); // going -y project to bottom
1475  else distY = 999999; // not moving in Y? May be 2D track... we can still proceed though
1476 
1477  if(dir.X() > 0) distX = (fMCEdge-vertex.X())/dir.X(); // going +x project +x wall
1478  else if(dir.X() < 0) distX = (-fMCEdge-vertex.X())/dir.X(); // going -x project to -x wall
1479  else distX = 999999; // not moving in X? may be 2D track... we can still proceed though
1480 
1481  return std::min(std::min(distX,distY),distZ);
1482  }
1483 
1484  else { // if vertex is outside the muon catcher (still has to be in the rest of the ND; won't work for any arbitrary vertex and direction outside of ND)
1485  if(vertex.Z() > fMCFront) return 0; // it's past the main ND but not in the muon catcher, so it is 'outside' the detector and we don't support it yet
1486  if(dir.Z() < 0) return 0; // this can't point into the muon catcher
1487  double Ytran = YPositionAtTransition(vertex, dir);
1488  double Xtran = XPositionAtTransition(vertex, dir);
1489  if(fabs(Xtran) < fMCEdge && Ytran < fMCTop && Ytran > -fMCEdge) { // this enters the muon catcher (but ignores air gap events)
1490  TVector3 vertex2(Xtran, Ytran, fMCFront+0.01);
1491  if(!InMuonCatcher(vertex2)) {
1492  std::cout << "Something weird with muon catcher distance projection in LiveGeometry... email this situation to bays@caltech.edu" << std::endl;
1493  return 0;
1494  }
1495  return ProjectedDistanceInMC(vertex2, dir);
1496  }
1497  }
1498 
1499  return 0; // should only reach this if not near detector or if weird error
1500  }
1501 
1502  //......................................................................
1504  {
1506  if(geom->DetId() != novadaq::cnv::kNEARDET) {
1507  std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
1508  return 0;
1509  }
1510 
1511  std::vector< OfflineChan > xhits, yhits;
1512 
1513  if(InMuonCatcher(vertex)) {
1514  TVector3 end = vertex + ProjectedDistanceInMC(vertex,dir)*dir;
1515  geom->CountCellsOnLineFast(vertex, end, xhits, yhits);
1516  return xhits.size() + yhits.size();
1517  }
1518 
1519  else {
1520  if(vertex.Z() > fMCFront) return 0; // it's past the main ND but not in the muon catcher, so it is 'outside' the detector and we don't support it yet
1521  if(dir.Z() < 0) return 0; // this can't point into the muon catcher
1522  double Ytran = YPositionAtTransition(vertex, dir);
1523  double Xtran = XPositionAtTransition(vertex, dir);
1524  if(fabs(Xtran) < fMCEdge && Ytran < fMCTop && Ytran > -fMCEdge) { // this enters the muon catcher (but ignores air gap events)
1525  TVector3 vertex2(Xtran, Ytran, fMCFront+0.01);
1526  if(!InMuonCatcher(vertex2)) {
1527  std::cout << "Something weird with muon catcher distance projection in LiveGeometry... email this situation to bays@caltech.edu" << std::endl;
1528  return 0;
1529  }
1530  return ProjectedCellsInMC(vertex2, dir);
1531  }
1532  }
1533 
1534  return 0;
1535  }
1536 
1537  //......................................................................
1538  int LiveGeometry::FullNDProjectedCells(TVector3 vertex, TVector3 dir) //combine main ND and muon catcher into one
1539  {
1541  if(geom->DetId() != novadaq::cnv::kNEARDET) {
1542  std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
1543  return 0;
1544  }
1545 
1546  double xyz[3], dxyz[3];
1547  xyz[0] = vertex.X();
1548  xyz[1] = vertex.Y();
1549  xyz[2] = vertex.Z();
1550  dxyz[0] = dir.X();
1551  dxyz[1] = dir.Y();
1552  dxyz[2] = dir.Z();
1553 
1554  if(IsPointLive(vertex)==false && InMuonCatcher(vertex)==false) return -5; // outside both
1555  if(IsPointLive(vertex)==true && DoWeEnterMC(vertex, dir)==false) return ProjectedLiveCellsToEdge(vertex,dir); // just main ND
1556  if(IsPointLive(vertex)==true && DoWeEnterMC(vertex, dir)==true) return (ProjectedLiveCellsToEdge(vertex,dir)+ProjectedCellsInMC(vertex,dir)); // main to MC
1558  return ProjectedCellsInMC(vertex,dir); }
1559  else { // MC to main
1560  TVector3 v2 = vertex + (ProjectedDistanceInMC(vertex,dir))*dir; //jump to boundary with main detector
1561  for(int i = 0; i<100; i++) {
1562  if(IsPointLive(v2)) return (ProjectedCellsInMC(vertex,dir) + ProjectedLiveCellsToEdge(v2,dir));
1563  v2 = v2 + 0.05*dir; //step along path until we're in main detector for sure
1564  }
1565  return ProjectedCellsInMC(vertex,dir); // never found main detector; something weird, just use MC part
1566  }
1567  return -111; // this should never happen
1568  }
1569 
1570  //......................................................................
1572  {
1574  if(geom->DetId() != novadaq::cnv::kNEARDET) {
1575  std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
1576  return 0;
1577  }
1578  if(dir.Z() == 0) return -9999;
1579  return vertex.X() + dir.X()*(fMCFront-vertex.Z())/dir.Z();
1580  }
1581 
1582  //......................................................................
1584  {
1586  if(geom->DetId() != novadaq::cnv::kNEARDET) {
1587  std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
1588  return 0;
1589  }
1590  if(dir.Z() == 0) return -9999;
1591  return vertex.Y() + dir.Y()*(fMCFront-vertex.Z())/dir.Z();
1592  }
1593 
1594  //......................................................................
1595  double LiveGeometry::DistanceToEdgeInMC(TVector3 vertex, int &wall)
1596  {
1598  if(geom->DetId() != novadaq::cnv::kNEARDET) {
1599  std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
1600  return 0;
1601  }
1602  wall = 0;
1603  if(!InMuonCatcher(vertex)) return 0; //NOT IN MUON CATCHER
1604  double mindist = 9999;
1605 
1606  if(mindist > (fMCEdge+vertex.X())) { // check -X direction
1607  mindist = fMCEdge+vertex.X();
1608  wall = 1;
1609  }
1610 
1611  if(mindist > (fMCEdge-vertex.X())) { // check +X direction
1612  mindist = fMCEdge-vertex.X();
1613  wall = 2;
1614  }
1615 
1616  if(mindist > (fMCEdge+vertex.Y())) { // check -Y direction
1617  mindist = fMCEdge+vertex.Y();
1618  wall = 3;
1619  }
1620 
1621  if(mindist > (fMCTop-vertex.Y())) { // check +Y direction
1622  mindist = fMCTop-vertex.Y();
1623  wall = 4;
1624  }
1625 
1626  if(mindist > (fMCBack-vertex.Z())) { // check +Z direction
1627  mindist = fMCBack-vertex.Z();
1628  wall = 6;
1629  }
1630 
1631  return mindist; // no need to check -Z because that is always live detector
1632  }
1633 
1634  //......................................................................
1635  int LiveGeometry::CellsToEdgeInMC(TVector3 vertex, int &wall)
1636  {
1638  if(geom->DetId() != novadaq::cnv::kNEARDET) {
1639  std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
1640  return 0;
1641  }
1642  std::vector< OfflineChan > xhits, yhits;
1643  TVector3 out = vertex;
1644  double dist = DistanceToEdgeInMC(vertex, wall);
1645  if(wall == 0) return 0;
1646  if(wall == 1) out.SetX(vertex.X()+dist);
1647  if(wall == 2) out.SetX(vertex.X()-dist);
1648  if(wall == 3) out.SetY(vertex.Y()+dist);
1649  if(wall == 4) out.SetY(vertex.Y()-dist);
1650  if(wall == 6) out.SetZ(vertex.Z()+dist);
1651  geom->CountCellsOnLineFast(vertex,out,xhits,yhits);
1652  return xhits.size()+yhits.size();
1653  }
1654 
1655  //......................................................................
1656  double LiveGeometry::ProjectedAirDist(TVector3 vertex, TVector3 dir)
1657  {
1659  if(geom->DetId() != novadaq::cnv::kNEARDET) {
1660  std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
1661  return 0;
1662  }
1663 
1664  if(DoWeEnterMC(vertex, dir) != 2) return 0;
1665 
1666  TVector3 posattrans; // position along line at transition between main ND and muon catcher
1667  posattrans.SetX(vertex.X() + dir.X()*(fMCFront-vertex.Z())/dir.Z());
1668  posattrans.SetY(vertex.Y() + dir.Y()*(fMCFront-vertex.Z())/dir.Z());
1669  posattrans.SetZ(fMCFront);
1670 
1671  TVector3 posatmctop; // position along line at the top of the muon catcher. Anything going through air and into muon catcher must go through top face of muon catcher.
1672  posatmctop.SetX(vertex.X() + dir.X()*(fMCTop-vertex.Y())/dir.Y());
1673  posatmctop.SetY(fMCTop);
1674  posatmctop.SetZ(vertex.Z() + dir.Z()*(fMCTop-vertex.Y())/dir.Y());
1675 
1676  if(dir.Z() > 0 && vertex.Z() > fMCFront) posattrans = vertex; // these two are for the odd case of a vertex actually in the air gap
1677  if(dir.Z() < 0 && vertex.Z() > fMCFront) posatmctop = vertex;
1678 
1679  return (posatmctop-posattrans).Mag(); // distance between these points is the distance through the air
1680  }
1681 
1682  //......................................................................
1683  double LiveGeometry::ProjectedSteelDist(TVector3 vertex, TVector3 dir)
1684  {
1686  if(geom->DetId() != novadaq::cnv::kNEARDET) {
1687  std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
1688  return 0;
1689  }
1690 
1691  if(DoWeEnterMC(vertex,dir) == 0) return 0;
1692 
1693  double p1[3] = {0}, p2[3], linedist; // start and end points of a line
1694 
1695  if(InMuonCatcher(vertex)) {
1696  p1[0] = vertex.X();
1697  p1[1] = vertex.Y();
1698  p1[2] = vertex.Z();
1699  }
1700 
1701  else {
1702  p1[0] = XPositionAtTransition(vertex,dir); // must have started in main ND; let's start counting steel at front of muon catcher
1703  p1[1] = YPositionAtTransition(vertex,dir);
1704  p1[2] = fMCFront;
1705  }
1706 
1707  linedist = std::sqrt((8*fMCEdge*fMCEdge) + ((fMCBack-fMCFront)*(fMCBack-fMCFront))) + 10; //longest possible distance through muon catcher (including possibility of air gap)
1708 
1709  p2[0] = p1[0] + linedist*dir.X(); // end of line is just line start point advanced along line direction the necessary distance to encompass the muon catcher
1710  p2[1] = p1[1] + linedist*dir.Y();
1711  p2[2] = p1[2] + linedist*dir.Z();
1712 
1713  std::vector< double > ds;
1714  std::vector < const TGeoMaterial * > mat;
1715 
1716  geom->MaterialsBetweenPoints(p1,p2,ds,mat);
1717 
1718  double steeldist = 0;
1719 
1720  for(unsigned int i = 0; i < mat.size(); i++) {
1721  if(strncmp("Steel",mat[i]->GetName(),5) == 0){
1722  steeldist+=ds[i];
1723  }
1724  }
1725 
1726  return steeldist;
1727 
1728  }
1729 
1730  //......................................................................
1732  {
1733  if(!fEvFlag) GetEvInfo(); // setup hasn't been called yet!
1734  return fNBadChan;
1735  }
1736 
1737  //......................................................................
1739  {
1740  if(!fEvFlag) GetEvInfo(); // setup hasn't been called yet!
1741  return fNDropouts;
1742  }
1743 
1744  //......................................................................
1746  {
1747  if(!fEvFlag) GetEvInfo(); // setup hasn't been called yet!
1748  return fTotalDropouts;
1749  }
1750 
1751  //......................................................................
1752  int LiveGeometry::DCMStatus(unsigned int dib, unsigned int dcm)
1753  {
1754  if(!fEvFlag) GetEvInfo(); // setup hasn't been called yet!
1755  if(fDCMs[dcm-1][dib-1]==2) return 0;
1756  if(fDCMs[dcm-1][dib-1]==0) return 2;
1757  if(fDCMs[dcm-1][dib-1]==1) return 1;
1758  return -1;
1759  }
1760 
1761  //......................................................................
1763  {
1764  if(!fSetup) SetupLiveGeo(); // setup hasn't been called yet!
1765  if(fParams.CheckBadChannels()) GetBCInfo(); // check bad channels every event in case of dropped DCMs
1766  FillBadBoxes(); // if we're checking each event, have to remake bad boxes every event. If we haven't set up yet, we aren't using yet, can avoid.
1767  fEvFlag = true;
1769  fNEvents++;
1770  }
1771 
1773  {
1774  if(!fEvFlag) GetEvInfo(); // setup hasn't been called yet!
1775 
1777  double mindist = 999999, tdist;
1778 
1779  if(geom->DetId() == novadaq::cnv::kNEARDET && InMuonCatcher(vertex)) {
1780  int wall;
1781  return DistanceToEdgeInMC(vertex, wall); // assuming the DCMs closest to the muon catcher are fine and the closest distance to bad in the muon catcher is just the detector edge
1782  }
1783 
1784  if(IsPointLive(vertex) == false) return 0; // already in a bad region
1785 
1786  // first get distances to detector edges; trivial
1787  tdist = vertex.Z() - std::min(fEdgeFrontXZ,fEdgeFrontYZ); // dist to front
1788  if(tdist < mindist) mindist = tdist;
1789 
1790  tdist = std::max(fEdgeBackXZ,fEdgeBackYZ) - vertex.Z(); // dist to back
1791  if(tdist < mindist) mindist = tdist;
1792 
1793  tdist = vertex.X() - fEdgeXMinusXZ;
1794  if(tdist < mindist) mindist = tdist;
1795 
1796  tdist = fEdgeXPlusXZ - vertex.X();
1797  if(tdist < mindist) mindist = tdist;
1798 
1799  tdist = vertex.Y() - fEdgeBottomYZ;
1800  if(tdist < mindist) mindist = tdist;
1801 
1802  tdist = fEdgeTopYZ = vertex.Y();
1803  if(tdist < mindist) mindist = tdist;
1804 
1805  // check bad boxes in XZ view
1806  for(unsigned int i = 0; i < fBadBoxXZ.size(); i++) {
1807 
1808  if(fBadBoxXZ[i].Z1 <= vertex.Z() && fBadBoxXZ[i].Z2 >= vertex.Z()) { // lined up in Z
1809  tdist = std::min(fabs(vertex.X()-fBadBoxXZ[i].XY1),fabs(vertex.X()-fBadBoxXZ[i].XY2));
1810  if(tdist < mindist) mindist = tdist;
1811  }
1812 
1813  else if(fBadBoxXZ[i].XY1 <= vertex.X() && fBadBoxXZ[i].XY2 >= vertex.X()) { // lined up in X
1814  tdist = std::min(fabs(vertex.Z()-fBadBoxXZ[i].Z1),fabs(vertex.Z()-fBadBoxXZ[i].Z2));
1815  if(tdist < mindist) mindist = tdist;
1816  }
1817 
1818  else { // not lined up in either x or z; closest will be one of the four corners
1819  float td1 = std::min((sqrt(pow(vertex.Z()-fBadBoxXZ[i].Z1,2)+pow(vertex.X()-fBadBoxXZ[i].XY1,2))),(sqrt(pow(vertex.Z()-fBadBoxXZ[i].Z1,2)+pow(vertex.X()-fBadBoxXZ[i].XY2,2))));
1820  float td2 = std::min((sqrt(pow(vertex.Z()-fBadBoxXZ[i].Z2,2)+pow(vertex.X()-fBadBoxXZ[i].XY1,2))),(sqrt(pow(vertex.Z()-fBadBoxXZ[i].Z2,2)+pow(vertex.X()-fBadBoxXZ[i].XY2,2))));
1821  tdist = std::min(td1,td2);
1822  if(tdist < mindist) mindist = tdist;
1823  }
1824  }
1825 
1826  // check bad boxes in YZ view
1827  for(unsigned int i = 0; i < fBadBoxYZ.size(); i++) {
1828 
1829  if(fBadBoxYZ[i].Z1 <= vertex.Z() && fBadBoxYZ[i].Z2 >= vertex.Z()) { // lined up in Z
1830  tdist = std::min(fabs(vertex.Y()-fBadBoxYZ[i].XY1),fabs(vertex.Y()-fBadBoxYZ[i].XY2));
1831  if(tdist < mindist) mindist = tdist;
1832  }
1833 
1834  else if(fBadBoxYZ[i].XY1 <= vertex.Y() && fBadBoxYZ[i].XY2 >= vertex.Y()) { // lined up in X
1835  tdist = std::min(fabs(vertex.Z()-fBadBoxYZ[i].Z1),fabs(vertex.Z()-fBadBoxYZ[i].Z2));
1836  if(tdist < mindist) mindist = tdist;
1837  }
1838 
1839  else { // not lined up in either x or z; closest will be one of the four corners
1840  float td1 = std::min((sqrt(pow(vertex.Z()-fBadBoxYZ[i].Z1,2)+pow(vertex.Y()-fBadBoxYZ[i].XY1,2))),(sqrt(pow(vertex.Z()-fBadBoxYZ[i].Z1,2)+pow(vertex.Y()-fBadBoxYZ[i].XY2,2))));
1841  float td2 = std::min((sqrt(pow(vertex.Z()-fBadBoxYZ[i].Z2,2)+pow(vertex.Y()-fBadBoxYZ[i].XY1,2))),(sqrt(pow(vertex.Z()-fBadBoxYZ[i].Z2,2)+pow(vertex.Y()-fBadBoxYZ[i].XY2,2))));
1842  tdist = std::min(td1,td2);
1843  if(tdist < mindist) mindist = tdist;
1844  }
1845  }
1846 
1847  return mindist;
1848 
1849  }
1850 
1852 } // end namespace geo
1853 ////////////////////////////////////////////////////////////////////////
int PlanesToEdge(TVector3 vertex)
T max(const caf::Proxy< T > &a, T b)
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
const XML_Char XML_Encoding * info
Definition: expat.h:530
void preBeginEvent(art::Event const &evt)
double GetBadBoxCorner(bool view, int coord, int i)
Get the corner of a bad box (hook for evend display)
int NDiBlocks()
gives number of active diblocks only, may be less than 14
Definition: RunHistory.h:394
Definition: event.h:34
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
int CellsToEdgeY(TVector3 vertex)
double fMCTop
muon catcher boundaries
Definition: LiveGeometry.h:136
bool IntersectsBox(const double xyz[], const double dxyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi)
Determine if a particle starting at xyz with direction dxyz will intersect a box defined by xlo...
Definition: Geo.cxx:97
unsigned int type
0 means uninstrumented, 1 means bad
Definition: LiveGeometry.h:165
bool IsDiBlockFullyInstrumented(int idb)
returns true if nInstrumentedFEBs in diblock (which counts FEBs that are instrumented, active, and unmasked only) is >= 700
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
LiveGeometry(const Parameters &params, art::ActivityRegistry &reg)
double ProjectedSteelDist(TVector3 vertex, TVector3 dir)
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
#define DEFINE_ART_SERVICE(svc)
Definition: ServiceMacros.h:93
int CellsToEdgeXMinus(TVector3 vertex)
double DistToEdgeXY(TVector3 vertex)
double DistToEdgeY(TVector3 vertex)
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
The instrumented geometry of one entire detector (near, far, ipnd)
Definition: LiveGeometry.h:21
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
DiBlock GetDiBlock(int i, bool loadAll=true)
get ith diblock is RH list (which only includes diblocks with activity), starting with i=0...
int NumBadBoxesX()
Size of the vector of bad boxes in the XZ view.
int CellsToTop(TVector3 vertex)
int CellsToBottom(TVector3 vertex)
int ProjectedCellsToNextBadRegion(TVector3 vertex, TVector3 dir)
Float_t ss
Definition: plot.C:24
double DistanceToEdgeInMC(TVector3 vertex, int &wall)
double DetLength() const
int fDCMs[12][14]
Definition: LiveGeometry.h:143
GlobalSignal< detail::SignalResponseType::LIFO, void(SubRun const &)> sPostBeginSubRun
bool fSetup
have we done setup yet this run
Definition: LiveGeometry.h:140
double DistToClosestBadRegion(TVector3 vertex)
double DistToEdgeX(TVector3 vertex)
double DistToFront(TVector3 vertex)
std::vector< BadBox > fBadBoxYZ
Definition: LiveGeometry.h:169
const PlaneGeo * Plane(unsigned int i) const
double DistToBack(TVector3 vertex)
double LiveDetectorVolume()
in ktons
int ProjectedCellsToEdge(TVector3 vertex, TVector3 dir)
void MaterialsBetweenPoints(const double *p1, const double *p2, std::vector< double > &ds, std::vector< const TGeoMaterial * > &mat) const
void GetBCInfo()
Get bad regions from BadChannels.
int ProjectedLiveCellsToEdge(TVector3 vertex, TVector3 dir)
recommended projection to use
float Z1
point 1 for minima, point 2 for maxima
Definition: LiveGeometry.h:161
void CountCellsOnLineFast(double x1, double y1, double z1, double x2, double y2, double z2, std::vector< OfflineChan > &Xhitsonline, std::vector< OfflineChan > &Yhitsonline)
Make list of cells in each view traversed by a line segment.
double dist
Definition: runWimpSim.h:113
void ProjectToBoxEdgeFromOutside(const double xyz[], const double dxyz[], int axis, double edge, double xyzout[])
Project from a position outside of a box to an edge of the box with coordinate value edge for the axi...
Definition: Geo.cxx:148
void postBeginSubRun(art::SubRun const &subrun)
bool ProjectedDistance(TVector3 vertex, TVector3 dir, double &distToEdge, double &distToNextBad, double &distDead)
double DistToTop(TVector3 vertex)
fhicl::Atom< bool > CheckBadChannels
Definition: LiveGeometry.h:27
double InstrumentedDetLength()
get instrumented detector length of downstream part
double ProjectedLiveDistanceToEdge(TVector3 vertex, TVector3 dir)
Far Detector at Ash River, MN.
std::string GetName(int i)
void SetInstrumentation()
calculate the detector end information
int ProjectedCellsInMC(TVector3 vertex, TVector3 dir)
double DistToBottom(TVector3 vertex)
double DistToEdgeZ(TVector3 vertex)
double fEdgeBottomYZ
proper 2D detector boundaries in YZ
Definition: LiveGeometry.h:135
int CellsToEdgeXY(TVector3 vertex)
double YPositionAtTransition(TVector3 vertex, TVector3 dir)
double InstrumentedDetFront()
get instrumented detector front of downstream part
std::vector< BadBox > fBadBoxXZ
Definition: LiveGeometry.h:169
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
void GetEvInfo()
Get event by event information.
int FullNDProjectedCells(TVector3 vertex, TVector3 dir)
Near Detector in the NuMI cavern.
int CellsToEdgeInMC(TVector3 vertex, int &wall)
Collect Geo headers and supply basic geometry functions.
double InstrumentedDetBack()
get instrumented detector back of downstream part
Float_t d
Definition: plot.C:236
bool fEvFlag
have we read in information for this event yet? (don&#39;t use preBeginEvent since race w/ BadChan) ...
Definition: LiveGeometry.h:141
int CellsToEdgeXPlus(TVector3 vertex)
void GetRHInfo()
Get uninstrumented regions from RunHistory.
const double j
Definition: BetheBloch.cxx:29
bool IsPointLiveInfo(TVector3 vertex, int &info)
Note the muon catcher is considered bad; use in combination with InMuonCatcher if needed...
void SetupLiveGeo()
Necessary set up before anything; get boudnaries, find and fill bad boxes.
void GetDetectorEdges()
Get proper detector boundaries for each view.
A very simple service to remember what detector we&#39;re working in.
OStream cout
Definition: OStream.cxx:6
Double_t Z2
Definition: plot.C:264
int CellsToEdgeX(TVector3 vertex)
Double_t Z1
Definition: plot.C:264
int ProjectedWall(TVector3 vertex, TVector3 dir)
if we assume the downstream InstrumentedDetEnds, what wall would this projected track exit from...
std::vector< DCM > dcm
Definition: RunHistory.h:311
fhicl::Atom< bool > Verbose
Definition: LiveGeometry.h:31
double ProjectedDistanceInMC(TVector3 vertex, TVector3 dir)
double DistToEdgeXMinus(TVector3 vertex)
TDirectory * dir
Definition: macro.C:5
double ProjectedDistanceToNextBadRegion(TVector3 vertex, TVector3 dir)
GlobalSignal< detail::SignalResponseType::FIFO, void(Event const &)> sPreProcessEvent
bool edge
are we connected to the detector edge?
Definition: LiveGeometry.h:166
bool IsPointLive(TVector3 vertex)
Note the muon catcher is considered bad; use in combination with InMuonCatcher if needed...
void geom(int which=0)
Definition: geom.C:163
double XPositionAtTransition(TVector3 vertex, TVector3 dir)
double ProjectedAirDist(TVector3 vertex, TVector3 dir)
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
int DCMStatus(unsigned int dib, unsigned int dcm)
returns: 0: RH says is bad; 1: good; 2: bad, either BC or DCM drop
double fEdgeXPlusXZ
proper 2D detector boundaries in XZ
Definition: LiveGeometry.h:134
T min(const caf::Proxy< T > &a, T b)
nova::dbi::DataType GetDataType()
Definition: RunHistory.h:452
bool InMuonCatcher(TVector3 vertex)
int nInstrumentedFEBs
number of FEBs that were instrumented and active in run and not masked
Definition: RunHistory.h:287
bool IsBad(int plane, int cell)
Helper for AttenCurve.
Definition: Path.h:10
int InstrumentedDetEnds(double &frontDwnstrm, double &backDwnstrm, double &frontUpstrm, double &backUpstrm)
give all detector ends information; most general
Simple object representing a (plane, cell) pair.
void FillBadBoxes()
Fill bad and uninstrumented regions of detector.
int NumBadBoxesY()
Size of the vector of bad boxes in the YZ view.
Encapsulate the geometry of one entire detector (near, far, ndos)
int PlanesToBack(TVector3 vertex)
double ProjectedDistanceToEdge(TVector3 vertex, TVector3 dir)
Eigen::MatrixXd mat
int DoWeEnterMC(TVector3 vertex, TVector3 dir)
double DistToEdgeXPlus(TVector3 vertex)
int PlanesToFront(TVector3 vertex)