LiveGeometry_service.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \file    LiveGeometry_service.cc
/// \brief
/// \author  bays@caltech.edu
////////////////////////////////////////////////////////////////////////

#include <algorithm>
#include <cassert>
#include <iostream>

// Framework includes
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Principal/Event.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// NOvA includes
#include "DAQChannelMap/DAQChannelMap.h"
#include "ChannelInfo/BadChanList.h"
#include "Geometry/Geometry.h"
#include "Geometry/LiveGeometry.h"
#include "GeometryObjects/CellGeo.h"
#include "GeometryObjects/Geo.h"
#include "GeometryObjects/OfflineChan.h"
#include "RunHistory/service/RunHistoryService.h"

#include "TGeoManager.h"

namespace geo
{
  //......................................................................
  LiveGeometry::LiveGeometry(const Parameters& params, art::ActivityRegistry& reg) :
    fParams(params())
  {
    // tell the activity registry to check in after runs have been loaded
    // into the event processing so that we can update configurations if necessary
    // SetupLiveGeo must be called before anything will work; will be called automatically for all existing functions
    // If making a new function, call by hand at beginning or make sure IsPointLive(Info) or ProjectedDistance is called first
    // as they also call it
    reg.sPostBeginSubRun.watch(this, &LiveGeometry::postBeginSubRun);  // beware of race conditions with RH and BC.
    reg.sPreProcessEvent.watch(this, &LiveGeometry::preBeginEvent);
    fSetup = false;
  }

  //......................................................................
  LiveGeometry::~LiveGeometry()
  {
    if(fSetup)
      mf::LogInfo("LiveGeometry") << "LiveGeometry found "
                                  << fTotalDropouts << " dropped DCMs out of "
                                  << fNEvents << " events.";
  }

  //......................................................................
  void LiveGeometry::postBeginSubRun(art::SubRun const& subrun)
  {
    // 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)
    art::ServiceHandle<nova::dbi::RunHistoryService> RunHist;
    if(fSetup) {
      for(int i = 0; i<12; i++) { // initialize DCMs
        for(int j = 0; j<14; j++) {
          fDCMs[i][j] = 1;
        }
      }
      //      if(RunHist->GetDataType()!=1) { // don't check RH or Bad Channels for MC
        GetRHInfo(); // get bad boxes based on RH (what is 'uninstrumented' or fDCMs=2)
        if(fParams.CheckBadChannels()) GetBCInfo(); // get bad boxes based on BadChannels (What is 'bad' or fDCMs=0)
        //      }
    }
  }// end of preBeginRun

  //----------------------------------------------------------
  void LiveGeometry::preBeginEvent(art::Event const & evt)
  {
    art::ServiceHandle<geo::Geometry> geom;

    fWall = 0;
    fEvFlag = false;

  }

  //...................................................................
  bool LiveGeometry::IsPointLiveInfo(TVector3 vertex, int &info) {  // are we in a good part of the detector, or a bad/dead one?  return true = good
    if(!fEvFlag) GetEvInfo();  // setup hasn't been called yet!

    info = 0;  // outside detector
    bool badxz=false, badyz = false, edgexz=false, edgeyz=false;

    if(vertex.X() > fEdgeXPlusXZ || vertex.X() < fEdgeXMinusXZ ||  // are we even inside the detector?
       vertex.Y() > fEdgeTopYZ   || vertex.Y() < fEdgeBottomYZ ||
       vertex.Z() < std::min(fEdgeFrontXZ,fEdgeFrontYZ) || vertex.Z() > std::max(fEdgeBackXZ,fEdgeBackYZ) ){
      LOG_DEBUG("LiveGeometry") << "outside of the detector";
      return false;
    }


    for(unsigned int i = 0; i < fBadBoxXZ.size(); i++) {  // check XZ view
      if(fBadBoxXZ[i].Z1  <= vertex.Z() && fBadBoxXZ[i].Z2  >= vertex.Z() &&
         fBadBoxXZ[i].XY1 <= vertex.X() && fBadBoxXZ[i].XY2 >= vertex.X()) {
        if(fBadBoxXZ[i].edge == true) edgexz = true;
        badxz = true;
      }
    }

    for(unsigned int i = 0; i < fBadBoxYZ.size(); i++) {  // check YZ view
      if(fBadBoxYZ[i].Z1  <= vertex.Z() && fBadBoxYZ[i].Z2  >= vertex.Z() &&
         fBadBoxYZ[i].XY1 <= vertex.Y() && fBadBoxYZ[i].XY2 >= vertex.Y()) {
        if(fBadBoxYZ[i].edge == true) edgeyz = true;
        badyz = true;
      }
    }

    if(badxz&&badyz) {  // both views bad
      if(edgexz&&edgeyz) info = 5; // and both are edges
      if(edgexz||edgeyz) info = 4; // only one is an edge
      else info = 3; // neither are edges
      return false;
    }
    else if(badxz&&edgexz) { // one view bad and edge
      info = 2;
      return false;
    }
    else if(badyz&&edgeyz) { // one view bad and edge
      info = 2;
      return false;
    }
    else if(badxz||badyz) { // one view bad, neither are edge
      info = 1;
      return false;
    }

    return true;
  }

  //...................................................................
  bool LiveGeometry::IsPointLive(TVector3 vertex) {
    int info;
    bool live = this->IsPointLiveInfo(vertex, info);
    if(!live) 
      LOG_DEBUG("LiveGeometry") 
        << "(" << vertex.X() 
        << ", " << vertex.Y()
        << ", " << vertex.Z()
        << ") is not live, info: " 
        << info;

    return live;
  }

  //...................................................................
  int LiveGeometry::ProjectedWall(TVector3 vertex, TVector3 dir)  { // what wall would we exit assuming the downstream instrumented detector ends?
    if(fBackDownstream==-5) SetInstrumentation();
    double tdist, mindist = 9999999;
    int wall = 0, twall = 0;
    art::ServiceHandle<geo::Geometry> geom;

    if(dir.Z() > 0) {
      tdist = (geom->DetLength()*(float(fBackDownstream)/float(fNumDiblocks))-vertex.Z())/dir.Z();  // going +z project to back
      twall = 6;
    }
    else if(dir.Z() < 0) {
      tdist = (geom->DetLength()*(float(fFrontDownstream)/float(fNumDiblocks))-vertex.Z())/dir.Z(); // going -z project to front
      twall = 5;
    }
    else {  // dz is zero.  This means something is wrong with the track
      return 0;
    }
    if(tdist < mindist) {
      mindist = tdist;
      wall = twall;
    }

    if(dir.Y() > 0) {
      tdist = (fEdgeTopYZ-vertex.Y())/dir.Y();  // going +y project to top
      twall = 4;
    }
    else if(dir.Y() < 0) {
      tdist = (fEdgeBottomYZ-vertex.Y())/dir.Y(); // going -y project to bottom
      twall = 3;
    }
    else tdist = 999999; // not moving in Y?  May be 2D track... we can still proceed though
    if(tdist < mindist) {
      mindist = tdist;
      wall = twall;
    }

    if(dir.X() > 0) {
      tdist = (fEdgeXPlusXZ-vertex.X())/dir.X();  // going +x project +x wall
      twall = 2;
    }
    else if(dir.X() < 0) {
      tdist = (fEdgeXMinusXZ-vertex.X())/dir.X(); // going -x project to -x wall
      twall = 1;
    }
    else tdist = 999999; // not moving in X? may be 2D track... we can still proceed though
    if(tdist < mindist) {
      mindist = tdist;
      wall = twall;
    }

    return wall;

  }

  //...................................................................
  double LiveGeometry::ProjectedDistanceToEdge(TVector3 vertex, TVector3 dir)  {
    double distToEdge, distToNextBad, distDead;
    if(ProjectedDistance(vertex, dir, distToEdge, distToNextBad, distDead)) return(distToEdge);
    else return -5;
  }

  //...................................................................
  double LiveGeometry::ProjectedDistanceToNextBadRegion(TVector3 vertex, TVector3 dir)  {
    double distToEdge, distToNextBad, distDead;
    if(!IsPointLive(vertex)) return -5;
    if(ProjectedDistance(vertex, dir, distToEdge, distToNextBad, distDead)) return(distToNextBad);
    else return -5;
  }

  //...................................................................
  double LiveGeometry::ProjectedLiveDistanceToEdge(TVector3 vertex, TVector3 dir)  {
    double distToEdge, distToNextBad, distDead;
    if(ProjectedDistance(vertex, dir, distToEdge, distToNextBad, distDead)) return(distToEdge-distDead);
    else return -5;
  }

  //...................................................................
  bool LiveGeometry::ProjectedDistance(TVector3 vertex, TVector3 dir, double &distToEdge, double &distToNextBad, double &distDead)  {

    if(!fEvFlag) GetEvInfo();  // setup hasn't been called yet!

    float minedgedist = 999999.0, minbaddist=99999.0, tdist; // initial values for minima and temporary distance variable
    double r[3], d[3];
    vertex.GetXYZ(r);
    dir.GetXYZ(d);
    distDead = 0;
    fWall = 0;

    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
      distToEdge = -999.;
      distToNextBad = -999.;
      distDead = 0;

      if(vertex.Z() < InstrumentedDetFront()) {  // vertex in front of detector
        if(dir.Z() >= 0) return false; // going wrong way!
        tdist = (vertex.Z() - InstrumentedDetFront()) / dir.Z(); // work in positive numbers, make negative at end
        if(tdist < minedgedist) minedgedist = tdist;  // z check; want smallest positive number here
      }

      if(vertex.Z() > InstrumentedDetBack()) { // vertex behind detector
        if(dir.Z() <= 0) return false; // wrong direction
        tdist = (vertex.Z() - InstrumentedDetBack()) / dir.Z();
        if(tdist < minedgedist) minedgedist = tdist;  // z check; want smallest positive number here
      }

      if(vertex.Y() > fEdgeTopYZ) {
        if(dir.Y() <= 0) return false;
        tdist = -(fEdgeTopYZ-vertex.Y())/dir.Y();  // project to top
        if(tdist < minedgedist) minedgedist = tdist; // y check
      }
      if(vertex.Y() < fEdgeBottomYZ) {
        if(dir.Y() >= 0) return false;
        tdist = -(fEdgeBottomYZ-vertex.Y())/dir.Y(); // project to bottom
        if(tdist < minedgedist) minedgedist = tdist; // y check
      }

      if(vertex.X() > fEdgeXPlusXZ) {
        if(dir.X() <= 0) return false;
        tdist = -(fEdgeXPlusXZ-vertex.X())/dir.X();
        if(tdist < minedgedist) minedgedist = tdist; // x check
      }

      if(vertex.X() < fEdgeXMinusXZ) {
        if(dir.X() >= 0) return false;
        tdist = -(fEdgeXMinusXZ-vertex.X())/dir.X();
        if(tdist < minedgedist) minedgedist = tdist; // x check
      }

      distToEdge = -minedgedist;
      distToNextBad = -minedgedist;
      distDead = 0;

      return false;
    } // end vertex outside live detector

    // first check distances to full detector boundaries
    if(dir.Z() > 0) tdist = (std::max(fEdgeBackYZ,fEdgeBackXZ)-vertex.Z())/dir.Z();  // going +z project to back
    else if(dir.Z() < 0) tdist = (std::min(fEdgeFrontYZ,fEdgeFrontXZ)-vertex.Z())/dir.Z(); // going -z project to front
    else {  // dz is zero.  This means something is wrong with the track
      distToEdge = -9999.;
      distToNextBad = -9999.;
      distDead = -9999.;
      return false;
    }
    if(tdist < minedgedist) minedgedist = tdist;

    if(dir.Y() > 0) tdist = (fEdgeTopYZ-vertex.Y())/dir.Y();  // going +y project to top
    else if(dir.Y() < 0) tdist = (fEdgeBottomYZ-vertex.Y())/dir.Y(); // going -y project to bottom
    else tdist = 999999; // not moving in Y?  May be 2D track... we can still proceed though
    if(tdist < minedgedist) minedgedist = tdist;

    if(dir.X() > 0) tdist = (fEdgeXPlusXZ-vertex.X())/dir.X();  // going +x project +x wall
    else if(dir.X() < 0) tdist = (fEdgeXMinusXZ-vertex.X())/dir.X(); // going -x project to -x wall
    else tdist = 999999; // not moving in X? may be 2D track... we can still proceed though
    if(tdist < minedgedist) minedgedist = tdist;

    // now check distances to bad boxes

    double xyz[3];

    for(unsigned int i = 0; i < fBadBoxXZ.size(); i++) { // look first in XZ view
      geo::ProjectToBoxEdgeFromOutside(r,d,2,fBadBoxXZ[i].Z1,xyz);  // check front face
      if(xyz[0] > fBadBoxXZ[i].XY1 && xyz[0] < fBadBoxXZ[i].XY2 && xyz[1] < fEdgeTopYZ && xyz[1] > fEdgeBottomYZ) {
        tdist = (xyz[2] - vertex.Z())/dir.Z();
        if(fBadBoxXZ[i].edge == false) {
          if(tdist < minbaddist  && tdist > 0) minbaddist = tdist;  // distance to a non-edge bad region
        }
        else if(tdist < minedgedist  && tdist > 0) minedgedist = tdist;  // distance to an edge
      }

      geo::ProjectToBoxEdgeFromOutside(r,d,2,fBadBoxXZ[i].Z2,xyz);  // check back face
      if(xyz[0] > fBadBoxXZ[i].XY1 && xyz[0] < fBadBoxXZ[i].XY2 && xyz[1] < fEdgeTopYZ && xyz[1] > fEdgeBottomYZ) {
        tdist = (xyz[2] - vertex.Z())/dir.Z();
        if(fBadBoxXZ[i].edge == false) {
          if(tdist < minbaddist && tdist > 0) minbaddist = tdist;  // distance to a non-edge bad region
        }
        else if(tdist < minedgedist && tdist > 0) minedgedist = tdist;  // distance to an edge
      }

      geo::ProjectToBoxEdgeFromOutside(r,d,0,fBadBoxXZ[i].XY1,xyz);  // -x face
      if(xyz[2] > fBadBoxXZ[i].Z1 && xyz[2] < fBadBoxXZ[i].Z2 && xyz[1] < fEdgeTopYZ && xyz[1] > fEdgeBottomYZ) {
        tdist = (xyz[2] - vertex.Z())/dir.Z();
        if(fBadBoxXZ[i].edge == false) {
          if(tdist < minbaddist && tdist > 0) minbaddist = tdist;  // distance to a non-edge bad region
        }
        else if(tdist < minedgedist && tdist > 0) minedgedist = tdist;  // distance to an edge
      }

      geo::ProjectToBoxEdgeFromOutside(r,d,0,fBadBoxXZ[i].XY2,xyz);  // +x face
      if(xyz[2] > fBadBoxXZ[i].Z1 && xyz[2] < fBadBoxXZ[i].Z2 && xyz[1] < fEdgeTopYZ && xyz[1] > fEdgeBottomYZ) {
        tdist = (xyz[2] - vertex.Z())/dir.Z();
        if(fBadBoxXZ[i].edge == false) {
          if(tdist < minbaddist && tdist > 0) minbaddist = tdist;  // distance to a non-edge bad region
        }
        else if(tdist < minedgedist && tdist > 0) minedgedist = tdist;  // distance to an edge
      }
    }

    for(unsigned int i = 0; i < fBadBoxYZ.size(); i++) { // now look in YZ view
      if(fBadBoxYZ[i].edge == false) continue;
      geo::ProjectToBoxEdgeFromOutside(r,d,2,fBadBoxYZ[i].Z1,xyz);  // check front face
      if(xyz[0] > fEdgeXMinusXZ && xyz[0] < fEdgeXPlusXZ && xyz[1] > fBadBoxYZ[i].XY1 && xyz[1] < fBadBoxYZ[i].XY2) {
        tdist = (xyz[2] - vertex.Z())/dir.Z();
        if(fBadBoxYZ[i].edge == false) {
          if(tdist < minbaddist && tdist > 0) minbaddist = tdist;  // distance to a non-edge bad region
        }
        else if(tdist < minedgedist && tdist > 0) minedgedist = tdist;  // distance to an edge
      }

      geo::ProjectToBoxEdgeFromOutside(r,d,2,fBadBoxYZ[i].Z2,xyz);  // check back face
      if(xyz[0] > fEdgeXMinusXZ && xyz[0] < fEdgeXPlusXZ && xyz[1] > fBadBoxYZ[i].XY1 && xyz[1] < fBadBoxYZ[i].XY2) {
        tdist = (xyz[2] - vertex.Z())/dir.Z();
        if(fBadBoxYZ[i].edge == false) {
          if(tdist < minbaddist && tdist > 0) minbaddist = tdist;  // distance to a non-edge bad region
        }
        else if(tdist < minedgedist && tdist > 0) minedgedist = tdist;  // distance to an edge
      }

      geo::ProjectToBoxEdgeFromOutside(r,d,1,fBadBoxYZ[i].XY1,xyz);  // bottom face
      if(xyz[2] > fBadBoxYZ[i].Z1 && xyz[2] < fBadBoxYZ[i].Z2 && xyz[0] > fEdgeXMinusXZ && xyz[0] < fEdgeXPlusXZ) {
        tdist = (xyz[2] - vertex.Z())/dir.Z();
        if(fBadBoxYZ[i].edge == false) {
          if(tdist < minbaddist && tdist > 0) minbaddist = tdist;  // distance to a non-edge bad region
        }
        else if(tdist < minedgedist && tdist > 0) minedgedist = tdist;  // distance to an edge
      }

      geo::ProjectToBoxEdgeFromOutside(r,d,1,fBadBoxYZ[i].XY2,xyz);  // top face
      if(xyz[2] > fBadBoxYZ[i].Z1 && xyz[2] < fBadBoxYZ[i].Z2 && xyz[0] > fEdgeXMinusXZ && xyz[0] < fEdgeXPlusXZ) {
        tdist = (xyz[2] - vertex.Z())/dir.Z();
        if(fBadBoxYZ[i].edge == false) {
          if(tdist < minbaddist && tdist > 0) minbaddist = tdist;  // distance to a non-edge bad region
        }
        else if(tdist < minedgedist && tdist > 0) minedgedist = tdist;  // distance to an edge
      }
    }

    distToEdge = minedgedist;
    if(!IsPointLive(vertex)) {  // are we starting in a live region?
      distToNextBad = 0;
      minbaddist = 0;
    }
    else distToNextBad = std::min(minbaddist,minedgedist);

    if(distToNextBad >= distToEdge) {  // the closest bad region is an edge, and we didn't start in a bad region - we're done
      distDead = 0;
      return true;
    }

    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
      bool liveYet = true;  // have we encountered any live material yet
      if(!IsPointLive(vertex)) liveYet = false;  // starting in a bad region
      int info;  // for IsPointLiveInfo return
      TVector3 position = vertex + minbaddist*dir;
      double liveDist = 0;
      double ss = 0.1; // step size (cm)
      unsigned int maxstep = 100000; // keep it at 100 m
      for(unsigned int istep = 0; istep < maxstep; istep++) { // should never get to (100 m), but prevent infinite loop if error
        position = position + ss*dir;  // step by step size
        if(!IsPointLiveInfo(position,info)) {
          if(liveYet) {  // once we've seen live material, the next bad edge region means we're done
            if(info==5||info==4||info==2||info==0) {
              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.
              return true; // we've arrived at an edge, we're done
            }
            if(info==3) distDead+=ss; // dead in both views;
            if(info==1) distDead+=ss/2.0; // dead in only one view, count for half
          }
          else { // we didn't start in a live region... need to get to one before we start looking for edge
            if(info==5||info==4||info==3) distDead+=ss;
            else distDead+=ss/2.0;
          }
        }
        else {  // we found live material
          liveYet = true;
          liveDist+=ss;
          continue;
        }
      }
    } // end determine distDead else
    distToNextBad = -5.;
    distToEdge = -5.;
    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.

  }

  //...................................................................
  int LiveGeometry::ProjectedCellsToEdge(TVector3 vertex, TVector3 dir)  {
    art::ServiceHandle<geo::Geometry> geom;
    double distToEdge, distToNextBad, distDead;
    if(ProjectedDistance(vertex, dir, distToEdge, distToNextBad, distDead)) {
      TVector3 r2 = vertex + distToEdge*dir;
      std::vector< OfflineChan > xhits, yhits;
      geom->CountCellsOnLineFast(vertex, r2, xhits, yhits);
      return(xhits.size() + yhits.size());
    }
    else return -5;
  }

  //...................................................................
  int LiveGeometry::ProjectedCellsToNextBadRegion(TVector3 vertex, TVector3 dir)  {
    art::ServiceHandle<geo::Geometry> geom;

    double distToEdge, distToNextBad, distDead;
    if(!IsPointLive(vertex)) return -5;
    if(ProjectedDistance(vertex, dir, distToEdge, distToNextBad, distDead)) {
      TVector3 r2 = vertex + distToNextBad*dir;
      std::vector< OfflineChan > xhits, yhits;
      geom->CountCellsOnLineFast(vertex, r2, xhits, yhits);
      return(xhits.size() + yhits.size());
    }
    else return -5;
  }

  //...................................................................
  int LiveGeometry::ProjectedLiveCellsToEdge(TVector3 vertex, TVector3 dir)  {
    art::ServiceHandle<geo::Geometry> geom;

    double distToEdge, distToNextBad, distDead;
    if(ProjectedDistance(vertex, dir, distToEdge, distToNextBad, distDead)) {
      TVector3 r2 = vertex + distToEdge*dir;
      std::vector< OfflineChan > xhits, yhits;
      geom->CountCellsOnLineFast(vertex, r2, xhits, yhits);
      int nlivecells = (xhits.size() + yhits.size()) * (1- (distDead/distToEdge) ) + 1; // scale live cells by live distance, round up
      return(nlivecells);
    }
    else return -5;
  }

  //...................................................................
  double LiveGeometry::DistToEdgeZ(TVector3 vertex)  {

    return std::min(DistToBack(vertex),DistToFront(vertex));
  }

  //...................................................................
  double LiveGeometry::DistToFront(TVector3 vertex)  {

    if(!IsPointLive(vertex)) return vertex.Z()-InstrumentedDetFront();  // already in a dead region
    double fDistFront = vertex.Z()-std::min(fEdgeFrontXZ,fEdgeFrontYZ);
    for(unsigned int i = 0; i < fBadBoxXZ.size(); i++) {  // check XZ view bad regions
      if(fBadBoxXZ[i].XY1 <= vertex.X() && fBadBoxXZ[i].XY2 >= vertex.X()) {
        if(vertex.Z() >= fBadBoxXZ[i].Z2) fDistFront = std::min(fDistFront,vertex.Z()-fBadBoxXZ[i].Z2);
      }
    }
    for(unsigned int i = 0; i < fBadBoxYZ.size(); i++) {  // check YZ view bad regions
      if(fBadBoxYZ[i].XY1 <= vertex.Y() && fBadBoxYZ[i].XY2 >= vertex.Y()) {
        if(vertex.Z() >= fBadBoxYZ[i].Z2) fDistFront = std::min(fDistFront,vertex.Z()-fBadBoxYZ[i].Z2);
      }
    }

    return fDistFront;
  }

  //...................................................................
  double LiveGeometry::DistToBack(TVector3 vertex)  {

    if(!IsPointLive(vertex)) return InstrumentedDetBack()-vertex.Z();  // already in a dead region
    double fDistBack = 0;
    fDistBack = std::max(fEdgeBackXZ,fEdgeBackYZ)-vertex.Z();

    for(unsigned int i = 0; i < fBadBoxXZ.size(); i++) {  // check XZ view bad regions
      if(fBadBoxXZ[i].XY1 <= vertex.X() && fBadBoxXZ[i].XY2 >= vertex.X()) {
        if(vertex.Z() <= fBadBoxXZ[i].Z1) fDistBack = std::min(fDistBack,fBadBoxXZ[i].Z1-vertex.Z());
      }
    }

    for(unsigned int i = 0; i < fBadBoxYZ.size(); i++) {  // check YZ view bad regions
      if(fBadBoxYZ[i].XY1 <= vertex.Y() && fBadBoxYZ[i].XY2 >= vertex.Y()) {
        if(vertex.Z() <= fBadBoxYZ[i].Z1) fDistBack = std::min(fDistBack,fBadBoxYZ[i].Z1-vertex.Z());
      }
    }
    return fDistBack;
  }

  //...................................................................
  double LiveGeometry::DistToEdgeY(TVector3 vertex)  {

    return std::min(DistToTop(vertex),DistToBottom(vertex));
  }

  //...................................................................
  double LiveGeometry::DistToTop(TVector3 vertex)  {

    double fDistTop = fEdgeTopYZ-vertex.Y();
    if(!IsPointLive(vertex)) return fDistTop;  // already in a dead region

    for(unsigned int i = 0; i < fBadBoxYZ.size(); i++) {  // check YZ view bad regions
      if(fBadBoxYZ[i].Z1 <= vertex.Z() && fBadBoxYZ[i].Z2 >= vertex.Z()) {
        if(vertex.Y() <= fBadBoxYZ[i].XY1) fDistTop = std::min(fDistTop,fBadBoxYZ[i].XY1-vertex.Y());
      }
    }

    return fDistTop;
  }

  //...................................................................
  double LiveGeometry::DistToBottom(TVector3 vertex)  {

    double fDistBottom = vertex.Y()-fEdgeBottomYZ;
    if(!IsPointLive(vertex)) return fDistBottom;  // already in a dead region

    for(unsigned int i = 0; i < fBadBoxYZ.size(); i++) {  // check YZ view bad regions
      if(fBadBoxYZ[i].Z1 <= vertex.Z() && fBadBoxYZ[i].Z2 >= vertex.Z()) {
        if(vertex.Y() >= fBadBoxYZ[i].XY1) fDistBottom = std::min(fDistBottom,vertex.Y()-fBadBoxYZ[i].XY1);
      }
    }

    return fDistBottom;
  }

  //...................................................................
  double LiveGeometry::DistToEdgeX(TVector3 vertex)  {

    return std::min(DistToEdgeXPlus(vertex),DistToEdgeXMinus(vertex));
  }

  //...................................................................
  double LiveGeometry::DistToEdgeXY(TVector3 vertex)  {

    return std::min(DistToEdgeX(vertex),DistToEdgeY(vertex));
  }


  //...................................................................
  double LiveGeometry::DistToEdgeXPlus(TVector3 vertex)  {

    double fDistXPlus = fEdgeXPlusXZ-vertex.X(); // basic projected dists
    if(!IsPointLive(vertex)) return fDistXPlus;  // already in a dead region

    for(unsigned int i = 0; i < fBadBoxXZ.size(); i++) {  // check XZ view bad regions
      if(fBadBoxXZ[i].Z1 <= vertex.Z() && fBadBoxXZ[i].Z2 >= vertex.Z()) {
        if(vertex.X() <= fBadBoxXZ[i].XY1) fDistXPlus = std::min(fDistXPlus,fBadBoxXZ[i].XY1-vertex.X());
      }
    }
    return fDistXPlus;
  }

  //...................................................................
  double LiveGeometry::DistToEdgeXMinus(TVector3 vertex)  {

    double fDistXMinus = vertex.X()-fEdgeXMinusXZ;
    if(!IsPointLive(vertex)) return fDistXMinus;  // already in a dead region

    for(unsigned int i = 0; i < fBadBoxXZ.size(); i++) {  // check XZ view bad regions
      if(fBadBoxXZ[i].Z1 <= vertex.Z() && fBadBoxXZ[i].Z2 >= vertex.Z()) {
        if(vertex.X() >= fBadBoxXZ[i].XY1) fDistXMinus = std::min(fDistXMinus,vertex.X()-fBadBoxXZ[i].XY1);
      }
    }
    return fDistXMinus;
  }

  //...................................................................
  int LiveGeometry::CellsToEdgeXPlus(TVector3 vertex)  {

    if(DistToEdgeXPlus(vertex) >= 0) return (int)(DistToEdgeXPlus(vertex)*float(fNumCells+1)/(fEdgeXPlusXZ-fEdgeXMinusXZ));
    else return -999;
  }

  //...................................................................
  int LiveGeometry::CellsToEdgeXMinus(TVector3 vertex)  {

    if(DistToEdgeXMinus(vertex) >= 0) return (int)(DistToEdgeXMinus(vertex)*float(fNumCells+1)/(fEdgeXPlusXZ-fEdgeXMinusXZ));
    else return -999;
  }

  //...................................................................
  int LiveGeometry::CellsToEdgeX(TVector3 vertex)  {

    return std::min(CellsToEdgeXMinus(vertex),CellsToEdgeXPlus(vertex));
  }

  //...................................................................
  int LiveGeometry::CellsToTop(TVector3 vertex)  {

    if(DistToTop(vertex) >= 0) return (int)(DistToTop(vertex)*float(fNumCells+1)/(fEdgeTopYZ-fEdgeBottomYZ));
    else return -999;
  }

  //...................................................................
  int LiveGeometry::CellsToBottom(TVector3 vertex)  {

    if(DistToBottom(vertex) >= 0) return (int)(DistToBottom(vertex)*float(fNumCells+1)/(fEdgeTopYZ-fEdgeBottomYZ));
    else return -999;
  }

  //...................................................................
  int LiveGeometry::CellsToEdgeY(TVector3 vertex)  {

    return std::min(CellsToTop(vertex),CellsToBottom(vertex));
  }

  //...................................................................
  int LiveGeometry::CellsToEdgeXY(TVector3 vertex)  {

    return std::min(CellsToEdgeY(vertex),CellsToEdgeX(vertex));
  }

  //...................................................................
  int LiveGeometry::PlanesToFront(TVector3 vertex)  {

    if(DistToFront(vertex) >= 0) return (int)(DistToFront(vertex)*float(fNumDiblocks*64)/(fEdgeBackXZ-fEdgeFrontXZ));
    return -999;
  }

  //...................................................................
  int LiveGeometry::PlanesToBack(TVector3 vertex)  {

    if(DistToBack(vertex) >= 0) return (int)(DistToBack(vertex)*float(fNumDiblocks*64)/(fEdgeBackXZ-fEdgeFrontXZ));
    else return -999;
  }

  //...................................................................
  int LiveGeometry::PlanesToEdge(TVector3 vertex)  {

    return std::min(PlanesToFront(vertex),PlanesToBack(vertex));
  }

  //...................................................................
  double LiveGeometry::LiveDetectorVolume()  { // in ktons

    art::ServiceHandle<nova::dbi::RunHistoryService> RunHist;
    if(!fEvFlag) GetEvInfo();  // setup hasn't been called yet!
    return (double)(RunHist->NActiveChannels())*fMassKtons/float(fNumChannels);
  }

  //...................................................................
  int LiveGeometry::NumBadBoxesX()
  {
    if(!fEvFlag) GetEvInfo();  // setup hasn't been called yet!
    return fBadBoxXZ.size();
  }

  //...................................................................
  int LiveGeometry::NumBadBoxesY()
  {
    if(!fEvFlag) GetEvInfo();  // setup hasn't been called yet!
    return fBadBoxYZ.size();
  }

  //...................................................................

  double LiveGeometry::GetBadBoxCorner(bool view, int coord, int i) {
    if(!fEvFlag) GetEvInfo();  // setup hasn't been called yet!
    if (view==false) {  // XZ view
      if(coord==1) return fBadBoxXZ[i].XY1;
      else if(coord==2) return fBadBoxXZ[i].XY2;
      else if(coord==3) return fBadBoxXZ[i].Z1;
      else if(coord==4) return fBadBoxXZ[i].Z2;
      else return -999;
    }

    if (view==true) { // YZ view
      if(coord==1) return fBadBoxYZ[i].XY1;
      else if(coord==2) return fBadBoxYZ[i].XY2;
      else if(coord==3) return fBadBoxYZ[i].Z1;
      else if(coord==4) return fBadBoxYZ[i].Z2;
      else return -999;
    }
    return -999; // should never get here
  } // end GetBadBoxCorner

  // Private Functions

  //...................................................................
  void LiveGeometry::SetupLiveGeo()
  {

    art::ServiceHandle<nova::dbi::RunHistoryService> RunHist;
    art::ServiceHandle<geo::Geometry> geom;

    for(int i = 0; i<12; i++) { // initialize DCMs
      for(int j = 0; j<14; j++) {
        fDCMs[i][j] = 1;
      }
    }

    if(geom->DetId() == novadaq::cnv::kFARDET) {
      fNumDiblocks = 14;
      fNumDCMs = 12;
      fNumCells = 383;
      fNumChannels = 344064;
      fMassKtons = 14.0;
      if(fParams.Verbose()) std::cout << "LiveGeometry: Detector set as Far Detector" << std::endl;
    }
    else if (geom->DetId() == novadaq::cnv::kNEARDET) {
      fNumDiblocks = 3; // ignore muon catcher for now
      fNumDCMs = 4;
      fNumCells = 95;
      fNumChannels = 18432;  // no muon catcher
      fMassKtons = 0.3;
      if(fParams.Verbose()) std::cout << "LiveGeometry: Detector set as Near Detector" << std::endl;
    }
    else if(geom->DetId() == novadaq::cnv::kTESTBEAM) {
      fNumDiblocks = 1;
      fNumDCMs = 3;
      fNumCells = 63;
      fNumChannels = 4032;
      fMassKtons = 0.03;  // this should be fixed once detector is filled
      if(fParams.Verbose()) std::cout << "LiveGeometry: Detector set as Test Beam" << std::endl;
    }

    // START GETTING CORRECT DETECTOR BOUNDARIES

    // 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
    // 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
    // have that ~12 cm of extra space with no cells in it when thinking about what happens in the XZ view, for instance.

    GetDetectorEdges();

    if(fParams.Verbose()) {
      std::cout << "LiveGeometry: Detector Edges: XZ front/back/-X/+X " << fEdgeFrontXZ << " " << fEdgeBackXZ << " " << fEdgeXMinusXZ << " " << fEdgeXPlusXZ << std::endl;
      std::cout << "LiveGeometry: Detector Edges: YZ front/back/top/bottom " << fEdgeFrontYZ << " " << fEdgeBackYZ << " " << fEdgeTopYZ << " " << fEdgeBottomYZ << std::endl;
    }

    // END GETTING CORRECT DETECTOR BOUNDARIES

    GetRHInfo(); // get bad boxes based on RH (what is 'uninstrumented' or fDCMs=2)
    //    if(fParams.CheckBadChannels()) GetBCInfo(); // get bad boxes based on BadChannels (What is 'bad' or fDCMs=0)
    // now we know which DCMs are 'good' (1), 'bad' (0), or 'uninstrumented' (2).  Now we need to make bad boxes in coordinate space.

    //    FillBadBoxes();

    fSetup = true;

  }

  //...................................................................
  void LiveGeometry::GetRHInfo()
  {
    // Check RunHistory to see if the DCM is instrumented.  If not, set array to 2
    art::ServiceHandle<nova::dbi::RunHistoryService> RunHist;
    art::ServiceHandle<geo::Geometry> geom;
  int numinst[14][12] = {{0}};  // one entry for each DCM in full det; how many instrumented FEBs are in that DCM? start at zero
    if(RunHist->NDiBlocks()==0) {
      for(int i = 0; i<12; i++) { // initialize DCMs
        for(int j = 0; j<14; j++) {
          fDCMs[i][j] = 2;
        }
      }
      return; // No RH info for this run!  Ignore RH.
    }
    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.
      // 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.
      // when we do something like take diblock 1 out of partition 1 these will not match up, so treat generally
      nova::dbi::RunHistory::DiBlock dib = RunHist->GetDiBlock(idib, false);  // false flag prevents loading in unnecessary pixel masks which is slooooooow
      for(int idcm = 0; idcm < int(dib.dcm.size()); idcm++) {  // loop over all DCMs in diblocks
        nova::dbi::RunHistory::DCM dcm = dib.dcm[idcm];
        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;
        numinst[dib.num-1][idcm] = dcm.nInstrumentedFEBs;
      }
    }

    for(int idib = 0; idib < fNumDiblocks; idib++) { // now loop over REAL diblocks to check so things match up later
      for(int idcm = 0; idcm < fNumDCMs; idcm++) {
        if(geom->DetId() == novadaq::cnv::kNEARDET && (idcm == 0 || idcm == 3)) {if(numinst[idib][idcm] < 26) fDCMs[idcm][idib] = 2;}
        else if(numinst[idib][idcm] < 56) fDCMs[idcm][idib] = 2;
      }
    }
  }

  //...................................................................
  void LiveGeometry::GetBCInfo()
  {
    // Now check only instrumented DCMs for bad channels.  Use criteria:
    // require 16 good pixels for a good FEB, require 56 good FEBs for a good DCM
    // add criteria: no more than 200 total bad pixels out of 2048, about 10%)
    // Mapping notes:
    // DCMs 1-6 verticals, DCMs 7-12 horizontal
    // plane 0-63 diblock 1, 64-127 diblock 2 ... 832-895 diblock 14
    // even number planes (0,2,4...62) horizontals (YZ), odd planes (1,3,5..63) verticals (XZ)
    // for vertical planes: cell 0-63 -> DCM 1, 64-127 -> DCM 2 ... 320-383 -> DCM 6
    // for horizontal planes: cell 0-63 -> DCM 12, 64-127 -> DCM 11 ... 320-383 -> DCM 7
    // for DCM 12: cells 0-31 is a FEB; plane 0 = FEB 32... plane 62 = FEB 63; cells 64-95 for DCM 11 etc
    // for DCM 12: cells 32-63 is a FEB; plane 0 = FEB 0... plane 62 = FEB 31; cells 96-127 for DCM 11 etc
    // for DCM 1: cells 0-31 is a FEB; plane 1 = FEB 0... plane 63 = FEB 31; cells 64-95 for DCM 2 etc
    // for DCM 1: cells 32-63 is a FEB; plane 1 = FEB 32... plane 63 = FEB 63; cells 96-127 for DCM 2 etc
    // bad channel map should show uninstrumented cells as bad, and we will track difference
    art::ServiceHandle<geo::Geometry> geom;
    art::ServiceHandle<chaninfo::BadChanList> BadList;
    fNBadChan = 0;
    fNDropouts = 0;

    if(geom->DetId() == novadaq::cnv::kTESTBEAM){
      std::cout<<"GetBCInfo() has incorrect mapping for TestBeam. Setting NBadChan and NDropouts to zero!"<<std::endl;
      return;
    }

    for(int i = 0; i<fNumDCMs; i++) {
      for(int j = 0; j<fNumDiblocks; j++) {
        if(fDCMs[i][j]==0) fDCMs[i][j]=1;  // reset what is 'bad' every event
      }
    }

    for(int idib = 0; idib < fNumDiblocks; idib++) {

      int jplane, jcell1, jcell2;  // actual plane and cell number
      int nfeb1, nfeb2;  // number of dead channels in FEB
      int nbadDCM[12] = {0}; // number of bad FEBs in DCM
      int nbadcount[12] = {0}; // number of bad pixels in DCM

      for(int idcm = 0; idcm < fNumDCMs/2; idcm++) { // verticals first (DCMs 1-6, indices 0-5)
        if(fDCMs[idcm][idib]==1) { // check if already tagged as uninstrumented
          for(int iplane = 0; iplane < 32; iplane++) { // each plane is 2 FEBs, for 64 FEBs per DCM
            nfeb1 = 0;  // reset these values for new FEBs
            nfeb2 = 0;
            jplane = (idib*64) + (2*iplane) + 1; // odd numbered planes
            for(int icell = 0; icell < 32; icell++) { // this loops over the pixels of those 2 FEBs
              jcell1 = icell + 64*idcm;
              jcell2 = icell + 32 + 64*idcm;
              if(BadList->IsBad(jplane,jcell1)) nfeb1++; // check if this pixel on FEB 1 is bad
              if(BadList->IsBad(jplane,jcell2)) nfeb2++; // check if this pixel on FEB 2 is bad
            } // end verticals pixels loop
            if(nfeb1 > 6) nbadDCM[idcm]++;  // if the FEB had too many bad channels, count it as a bad FEB
            if(nfeb2 > 6) nbadDCM[idcm]++;
            nbadcount[idcm]+=(nfeb1+nfeb2); // keep track of total bad pixels
          } // end verticals planes loop
        } // end if instrumented
      } // end verticals dcm loop

      for(int idcm = fNumDCMs/2; idcm < fNumDCMs; idcm++) { // now horizontals (DCMs 7-12, indices 6-11)
        if(fDCMs[idcm][idib]==1) { // check if already tagged as uninstrumented
          for(int iplane = 0; iplane < 32; iplane++) {  // each plane is 2 FEBs, for 64 FEBs per DCM
            nfeb1 = 0;  // reset these values for new FEBs
            nfeb2 = 0;
            jplane = (idib*64) + (2*iplane); // even numbered planes
            for(int icell = 0; icell < 32; icell++) {  // this loops over the pixels of those 2 FEBs
              jcell1 = icell + 64*(11-idcm);   // arranged opposite of verticals
              jcell2 = icell + 32 + 64*(11-idcm);
              if(BadList->IsBad(jplane,jcell1)) nfeb1++; // check if this pixel on FEB 1 is bad
              if(BadList->IsBad(jplane,jcell2)) nfeb2++; // check if this pixel on FEB 2 is bad
            } // end horizontals pixels loop
            if(nfeb1 > 6) nbadDCM[idcm]++;  // if the FEB had too many bad channels, count it as a bad FEB
            if(nfeb2 > 6) nbadDCM[idcm]++;
            nbadcount[idcm]+=(nfeb1+nfeb2); // keep track of total bad pixels
          } // end horizontals planes loop
        } // end if instrumented
      } // end horizontals dcm loop

      for(int idcm = 0; idcm < fNumDCMs; idcm++) { // we looked at all pixels in this diblock, check results then go to next one
        if(nbadDCM[idcm]>=63) fNDropouts++;
        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)
          if(idcm==0||idcm==3) {
            if(nbadDCM[idcm] > 4 && fDCMs[idcm][idib]==1) {
              fDCMs[idcm][idib] = 0;  // half as strict criteria for half instrumented DCMs
              fNBadChan++;
            }
          }
          else if((nbadDCM[idcm] > 8) && fDCMs[idcm][idib]==1) {
            fDCMs[idcm][idib] = 0; // full DCMs same as FD
            fNBadChan++;
          }
        }
        else if((nbadDCM[idcm] > 8) && fDCMs[idcm][idib]==1) {
          fDCMs[idcm][idib] = 0; // set DCM to bad if it fails criteria and isn't already uninstrumented
          fNBadChan++;
        }
        if(fParams.Verbose()) { // optional debugging information
          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;
        }

      } // end checking DCMs loop

    }  // end diblock loop
   }

  //...................................................................
  void LiveGeometry::GetDetectorEdges()
  {
    art::ServiceHandle<geo::Geometry> geom;
    double xyz1[3], xyz2[3];

    if(!(geom->DetId() == novadaq::cnv::kFARDET ||
         geom->DetId() == novadaq::cnv::kNEARDET ||
         geom->DetId() == novadaq::cnv::kTESTBEAM )) {
      std::cout << "UNKNOWN DETECTOR!!!" << std::endl;
      fEdgeXPlusXZ = 0;
      fEdgeXMinusXZ = 0;
      fEdgeFrontXZ = 0;
      fEdgeBackXZ = 0;
      fEdgeTopYZ = 0;
      fEdgeBottomYZ = 0;
      fEdgeFrontYZ = 0;
      fEdgeBackYZ = 0;
    }

    int firstXZ = -1; // first XZ plane in diblock (vertical)
    int lastXZ  = -1; // last  XZ plane in diblock
    int firstYZ = -1; // first YZ plane in diblock (horizontal)
    int lastYZ  = -1; // last  YZ plane in diblock
    // TB starts and ends with vertical plane
    if (geom->DetId() == novadaq::cnv::kTESTBEAM) {
      firstXZ = 0;
      lastXZ  = 62;
      firstYZ = 1;
      lastYZ  = 61; // no plane 63 in TB
    }
    // ND/FD start with a horizontal plane, end with vertical
    else {
      firstXZ = 1;
      lastXZ  = 63;
      firstYZ = 0;
      lastYZ  = 62;
    }

    // start with XZ view, which is easiest since z position does NOT change with cell number; this is odd plane numbers (except for TestBeam)
    geom->Plane(firstXZ)->Cell(0)->GetCenter(xyz1);
    fEdgeFrontXZ = xyz1[2]-3.5;  // the half width of the cell in Z is ~3.5 cm
    geom->Plane(lastXZ)->Cell(0)->GetCenter(xyz2);
    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.
    geom->Plane(firstXZ)->Cell(fNumCells)->GetCenter(xyz1);  // cell 383 is most +x, cell 0 is most -X
    geom->Plane(lastXZ)->Cell(fNumCells)->GetCenter(xyz2);
    fEdgeXPlusXZ = (xyz1[0]+xyz2[0])/2.0 + 2.3;  // 2.28 is approximate half width of cell here; again, round up slightly
    if (geom->DetId()== novadaq::cnv::kTESTBEAM) {
      geom->Plane(62)->Cell(0)->GetCenter(xyz1);  // last XZ plane for TB is 62
    }
    else {
      geom->Plane((fNumDiblocks*64)-1)->Cell(0)->GetCenter(xyz1);  // any cell in last plane will do for finding the back
    }
    fEdgeBackXZ = xyz1[2]+3.5;

    // 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.
    // 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.

    geom->Plane(firstYZ)->Cell(0)->GetCenter(xyz1);  // front bottom most cell
    geom->Plane(firstYZ)->Cell(63)->GetCenter(xyz2); // front top cell of frontmost DCM
    fEdgeFrontYZ = (xyz1[2]+xyz2[2])/2.0 - 3.5;  // average the 2 to get the 'front'
    geom->Plane(lastYZ)->Cell(0)->GetCenter(xyz2); // back most bottom cell of frontmost bottom DCM
    fEdgeBottomYZ = (xyz1[1]+xyz2[1])/2.0 - 2.3; // reuse same point 1 for bottom calculation
    geom->Plane(firstYZ)->Cell(fNumCells)->GetCenter(xyz1); // topmost frontmost cell, in topmost frontmost DCM
    geom->Plane(lastYZ)->Cell(fNumCells)->GetCenter(xyz2); // topmost farther back cell in topmost frontmost DCM
    fEdgeTopYZ = (xyz1[1]+xyz2[1])/2.0 + 2.3;  // top edge
    if (geom->DetId() == novadaq::cnv::kTESTBEAM) {
      geom->Plane(61)->Cell(fNumCells-63)->GetCenter(xyz1);  // frontmost topmost cell in backmost topmost DCM
      geom->Plane(61)->Cell(fNumCells)->GetCenter(xyz2);  // backmost topmost cell, in backmost topmost DCM
    }
    else {
      geom->Plane((fNumDiblocks*64)-2)->Cell(fNumCells-63)->GetCenter(xyz1);  // frontmost topmost cell in backmost topmost DCM
      geom->Plane((fNumDiblocks*64)-2)->Cell(fNumCells)->GetCenter(xyz2);  // backmost topmost cell, in backmost topmost DCM
    }
    fEdgeBackYZ = (xyz1[2]+xyz2[2])/2.0 + 3.5; // back edge

    // Muon catcher.  TODO: Make these more precise

    // Based on the gdml
    //fMCEdge = 199.65;
    //fMCTop = 75.3;
    //fMCFront = 1277.65;
    //fMCBack = 1577.54;

    // based on the event display
    fMCEdge = 190;
    fMCTop = 67.5;
    fMCFront = 1278.4;
    fMCBack = 1585;
  }

  //...................................................................
  void LiveGeometry::FillBadBoxes()
  {
    // 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.

    fBadBoxXZ.clear();  // get rid of old bad boxes - we're redoing them
    fBadBoxYZ.clear();

    int DCM[12][14] = {{0}}; // a third array that merges the information of the other two.  This way fDCMs is never altered.

    for(int i = 0; i<fNumDCMs; i++) {
      for(int j = 0; j<fNumDiblocks; j++) {
        if(fDCMs[i][j]==0||fDCMs[i][j]==2) DCM[i][j]=0;
        else DCM[i][j]=1;
      }
    }

    art::ServiceHandle<geo::Geometry> geom;

    if(geom->DetId() == novadaq::cnv::kTESTBEAM){
      std::cout<<"FillBadBoxes() not valid for TestBeam. Not filling any."<<std::endl;
      return;
    }

    double xyz1[3], xyz2[3];
    // do some clustering.  Only cluster uninstrumented; DCMs that fail event check always left as single boxes.  This helps distinguish
    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.
    for(int i = 0; i<fNumDiblocks+1; i++) { // initialize to true
      xzview[i] = true;
      yzview[i] = true;
    }

    for(int idib = 0; idib < fNumDiblocks; idib++) {
      for(int idcm = 0; idcm < fNumDCMs/2; idcm++) {
        if(fDCMs[idcm][idib]!=2) xzview[idib] = false; // if even a single DCM is not bad, the whole thing is false
      }
      if(xzview[idib]==true) for(int idcm = 0; idcm < fNumDCMs/2; idcm++) {
        DCM[idcm][idib] = 1; // we're handling these bad boxes now, don't want to again later
      }
      for(int idcm = fNumDCMs/2; idcm < fNumDCMs; idcm++) {
        if(fDCMs[idcm][idib]!=2) yzview[idib] = false; // if even a single DCM is not bad, the whole thing is false
      }
      if(yzview[idib]==true) for(int idcm = fNumDCMs/2; idcm < fNumDCMs; idcm++) {
        DCM[idcm][idib] = 1; // we're handling these bad boxes now, don't want to again later
      }
    }

    LiveGeometry::BadBox tBox;

    int cstart=-5; // front of clusters
    for(int idib = 0; idib<fNumDiblocks+1; idib++) { // lets make the bad boxes, starting with XZ

      if(xzview[idib]==false||idib==fNumDiblocks) { // nope this diblock isn't all bad, or at end
        if(cstart==-5) continue; // there wasn't a cluster ongoing, we'll let the individual DCM bad boxes sort anything else out later
        else if(cstart>=0) { // there was a cluster that we have now just ended.  Make the bad box!

          //verticals = XZ
          geom->Plane((cstart*64)+1)->Cell(fNumCells)->GetCenter(xyz1);  // Z is static with cell number in this view; any cell at front face works
          tBox.Z1 = xyz1[2]-3.5;                         // set front (corner 1)
          tBox.XY2 = xyz1[0]+2.3;                      // set +X (corner 2); first plane in DCM, highest cell number has largest +X
          geom->Plane((idib*64)-1)->Cell(0)->GetCenter(xyz1);  // back-most (of all detector) cell, -X most (of DCM, repeats with DCM) cell;
          tBox.Z2 = xyz1[2]+3.5;                         // set back (corner 2)
          tBox.XY1 = xyz1[0]-2.3;                      // set -X (corner 1)
          tBox.type = 0;
          tBox.edge = true;
          fBadBoxXZ.push_back(tBox);
          cstart=-5;
          continue;
        }
      }

      if(xzview[idib]==true) { // yes this diblock is completely uninstrumented in this view
        if(cstart==-5) { // no we haven't started any kind of cluster, so start it
          cstart = idib;
          continue;
        }
        else continue; // there is already a cluster ongoing, so just keep going
      }
    }

    cstart = -5;
    for(int idib = 0; idib<fNumDiblocks+1; idib++) { // lets make the bad boxes, now YZ

      if(yzview[idib]==false||idib==fNumDiblocks) { // nope this diblock isn't all bad
        if(cstart==-5) continue; // there wasn't a cluster ongoing, we'll let the individual DCM bad boxes sort anything else out later
        else if(cstart>=0) { // there was a cluster that we have now just ended.  Make the bad box!

          LiveGeometry::BadBox tBox;
          //horizontals = YZ
          geom->Plane(cstart*64)->Cell(fNumCells)->GetCenter(xyz1);  // front topmost most cell of uninstrumented region in YZ; most-front cell
          tBox.Z1 = xyz1[2]-3.5;                         // set front (corner 1)
          geom->Plane((idib*64)-2)->Cell(fNumCells)->GetCenter(xyz1);
          tBox.XY2 = xyz1[1]+2.3;                       // set top (corner 2)
          geom->Plane((idib*64)-2)->Cell(0)->GetCenter(xyz1);  // back bottom-most cell in XZ; is furthest back cell
          tBox.Z2 = xyz1[2]+3.5;                         // set back (corner 2)
          geom->Plane((fNumDiblocks-1)*64)->Cell(0)->GetCenter(xyz1);  // front bottom-most cell in back bottom-most DCM in XZ; is  cell
          tBox.XY1 = xyz1[1]-2.3;                      // set bottom (corner 1)
          tBox.type = 0;
          tBox.edge = true;
          fBadBoxYZ.push_back(tBox);
          cstart=-5;
          continue;
        }
      }

      if(yzview[idib]==true) { // yes this diblock is completely uninstrumented in this view
        if(cstart==-5) { // no we haven't started any kind of cluster, so start it
          cstart = idib;
          continue;
        }
        else continue; // there is already a cluster ongoing, so just keep going
      }

    }

    // now get bad boxes for any remaining bad or uninstrumented regions, one DCM at a time:

    for(int idib = 0; idib < fNumDiblocks; idib++) { // verticals first = XZ
      for(int idcm = 0; idcm < fNumDCMs/2; idcm++) {
        if(DCM[idcm][idib]==0) {
          if(geom->DetId() == novadaq::cnv::kFARDET) {
            geom->Plane(1+64*idib)->Cell(idcm*64)->GetCenter(xyz1);  // -x-most front
            geom->Plane(63+64*idib)->Cell(idcm*64)->GetCenter(xyz2); // -x-most back
            tBox.Z1 = xyz1[2]-3.5;  // Z positions are easy in this view; front
            tBox.Z2 = xyz2[2]+3.5;  // back
            tBox.XY1 = (xyz1[0]+xyz2[0])/2.0 - 2.3; // average DCM -x
            geom->Plane(1+64*idib)->Cell(63+idcm*64)->GetCenter(xyz1);  // +x-most front
            geom->Plane(63+64*idib)->Cell(63+idcm*64)->GetCenter(xyz2); // +x-most back
            tBox.XY2 = (xyz1[0]+xyz2[0])/2.0 + 2.3; // average DCM +x
            if(fDCMs[idcm][idib] == 2) tBox.type = 0;  // uninstrumented
            else tBox.type = 1;  // bad
            bool connected = true;  // are we connected to an edge?
            tBox.edge = false;
            if(idib==13 || idib == 0 || idcm == 0 || idcm == 5) tBox.edge = true; // are we immediately connected to the edge?
            //current check is just seeing if up, left, down, or right is an unbroken series of bad DCMs to the edge
            for(int i = idib+1; i < fNumDiblocks; i++) { // check for unbroken connection to the back
              if(DCM[idcm][i]==1) connected = false;
            }
            if(connected == true) tBox.edge = true;
            connected = true;
            for(int i = 0; i < idib; i++) {  // check for unbroken connection to front
              if(DCM[idcm][i]==1) connected = false;
            }
            if(connected == true) tBox.edge = true;
            connected = true;
            for(int i = 0; i < idcm; i++) { // check for unbroken connection to bottom
              if(DCM[i][idib]==1) connected = false;
            }
            if(connected == true) tBox.edge = true;
            connected = true;
            for(int i = idcm+1; i < fNumDCMs/2; i++) { // check for unbroken connection to top
              if(DCM[i][idib]==1) connected = false;
            }
            if(connected == true) tBox.edge = true;
          }

          if(geom->DetId() == novadaq::cnv::kNEARDET) {
            if(idcm==0) {
              geom->Plane(1+64*idib)->Cell(0)->GetCenter(xyz1);  // -x-most front
              geom->Plane(63+64*idib)->Cell(0)->GetCenter(xyz2); // -x-most back
            }
            else if(idcm==1) {
              geom->Plane(1+64*idib)->Cell(32)->GetCenter(xyz1);  // -x-most front
              geom->Plane(63+64*idib)->Cell(32)->GetCenter(xyz2); // -x-most back
            }
            tBox.Z1 = xyz1[2]-3.5;  // Z positions are easy in this view; front
            tBox.Z2 = xyz2[2]+3.5;  // back
            tBox.XY1 = (xyz1[0]+xyz2[0])/2.0 - 2.3; // average DCM -x
            geom->Plane(1+64*idib)->Cell(idcm*64+31)->GetCenter(xyz1);  // +x-most front
            geom->Plane(63+64*idib)->Cell(idcm*64+31)->GetCenter(xyz2); // +x-most back
            tBox.XY2 = (xyz1[0]+xyz2[0])/2.0 + 2.3; // average DCM +x
            if(fDCMs[idcm][idib] == 2) tBox.type = 0;  // uninstrumented
            else tBox.type = 1;  // bad
            tBox.edge = true;  // every ND DCM is an edge
          }
          fBadBoxXZ.push_back(tBox);
        }
      }
    }

    for(int idib = 0; idib < fNumDiblocks; idib++) { // horizontals second = YZ
      for(int idcm = fNumDCMs/2; idcm < fNumDCMs; idcm++) {
        if(DCM[idcm][idib]==0) {
          if(geom->DetId() == novadaq::cnv::kFARDET) {
            geom->Plane(64*idib)->Cell((11-idcm)*64)->GetCenter(xyz1); // bottom front into 1
            geom->Plane(64*idib)->Cell(63+(11-idcm)*64)->GetCenter(xyz2); // top front into 2
            tBox.Z1 = (xyz1[2]+xyz2[2])/2.0 - 3.5;  // average front
            geom->Plane(62+64*idib)->Cell((11-idcm)*64)->GetCenter(xyz2);  // bottom back into 2
            tBox.XY1 = (xyz1[1]+xyz2[1])/2.0 - 2.3; // average DCM bottom
            geom->Plane(62+64*idib)->Cell(63+(11-idcm)*64)->GetCenter(xyz1);  // top back into 1
            tBox.Z2 = (xyz1[2]+xyz2[2])/2.0 + 3.5;  // average back
            geom->Plane(64*idib)->Cell(63+(11-idcm)*64)->GetCenter(xyz2); // top front into 2
            tBox.XY2 = (xyz1[1]+xyz2[1])/2.0 + 2.3; // average DCM top
            if(fDCMs[idcm][idib] == 2) tBox.type = 0;  // uninstrumented
            else tBox.type = 1;  // bad
            bool connected = true;  // are we connected to an edge?
            tBox.edge = false;
            if(idib==13 || idib == 0 || idcm == fNumDCMs/2 || idcm == fNumDCMs) tBox.edge = true; // are we immediately connected to the edge?
            for(int i = idib+1; i < fNumDiblocks; i++) { // check for unbroken connection to the back
              if(DCM[idcm][i]==1) connected = false;
            }
            if(connected == true) tBox.edge = true;
            connected = true;
            for(int i = 0; i < idib; i++) {  // check for unbroken connection to front
              if(DCM[idcm][i]==1) connected = false;
            }
            if(connected == true) tBox.edge = true;
            connected = true;
            for(int i = fNumDCMs/2; i < idcm; i++) { // check for unbroken connection to +x side
              if(DCM[i][idib]==1) connected = false;
            }
            if(connected == true) tBox.edge = true;
            connected = true;
            for(int i = idcm+1; i < fNumDCMs; i++) { // check for unbroken connection -x side
              if(DCM[i][idib]==1) connected = false;
            }
            if(connected == true) tBox.edge = true;
          }

          else if(geom->DetId() == novadaq::cnv::kNEARDET) {
            if(idcm==2) {
              geom->Plane(64*idib)->Cell(32)->GetCenter(xyz1); // bottom front into 1
              geom->Plane(64*idib)->Cell(95)->GetCenter(xyz2); // top front into 2
              tBox.Z1 = (xyz1[2]+xyz2[2])/2.0 - 3.5;  // average front
              geom->Plane(62+64*idib)->Cell(32)->GetCenter(xyz2);  // bottom back into 2
              tBox.XY1 = (xyz1[1]+xyz2[1])/2.0 - 2.3; // average DCM bottom
              geom->Plane(62+64*idib)->Cell(95)->GetCenter(xyz1);  // top back into 1
              tBox.Z2 = (xyz1[2]+xyz2[2])/2.0 + 3.5;  // average back
              geom->Plane(64*idib)->Cell(95)->GetCenter(xyz2); // top front into 2
              tBox.XY2 = (xyz1[1]+xyz2[1])/2.0 + 2.3; // average DCM top
            }
            if(idcm==3) {
              geom->Plane(64*idib)->Cell(0)->GetCenter(xyz1); // bottom front into 1
              geom->Plane(64*idib)->Cell(31)->GetCenter(xyz2); // top front into 2
              tBox.Z1 = (xyz1[2]+xyz2[2])/2.0 - 3.5;  // average front
              geom->Plane(62+64*idib)->Cell(0)->GetCenter(xyz2);  // bottom back into 2
              tBox.XY1 = (xyz1[1]+xyz2[1])/2.0 - 2.3; // average DCM bottom
              geom->Plane(62+64*idib)->Cell(31)->GetCenter(xyz1);  // top back into 1
              tBox.Z2 = (xyz1[2]+xyz2[2])/2.0 + 3.5;  // average back
              geom->Plane(64*idib)->Cell(31)->GetCenter(xyz2); // top front into 2
              tBox.XY2 = (xyz1[1]+xyz2[1])/2.0 + 2.3; // average DCM top
            }

            if(fDCMs[idcm][idib] == 2) tBox.type = 0;  // uninstrumented
            else tBox.type = 1;  // bad
            tBox.edge = true; // all ND DCMs are edges
          }
          fBadBoxYZ.push_back(tBox);
        }
      }
    }

    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
      tBox.Z1 = 1277.2;
      tBox.Z2 = 1585.5;
      tBox.XY1 = 70.0;
      tBox.XY2 = 199.7;
      tBox.edge = true;
      fBadBoxYZ.push_back(tBox);
    }

    // TODO: change algorithm for check if connected to edge to get L shapes and such

    if(fParams.Verbose()) {
    // optional debugging printouts:
      for(int i = 0; i < int(fBadBoxXZ.size()); i++) {  // check XZ view bad regions  // to see the bad boxes for debugging
        std::cout << "LiveGeometry: BadBoxXZ: X1 : " << fBadBoxXZ[i].XY1 << " X2: " << fBadBoxXZ[i].XY2 << " Z1: " << fBadBoxXZ[i].Z1 << " Z2: " << fBadBoxXZ[i].Z2 << std::endl;
      }

      for(int i = 0; i < int(fBadBoxYZ.size()); i++) {  // check XZ view bad regions
        std::cout << "LiveGeometry: BadBoxYZ: Y1 : " << fBadBoxYZ[i].XY1 << " Y2: " << fBadBoxYZ[i].XY2 << " Z1: " << fBadBoxYZ[i].Z1 << " Z2: " << fBadBoxYZ[i].Z2 << std::endl;
      }

    }
  }

  //......................................................................
  void LiveGeometry::SetInstrumentation()
  {
    art::ServiceHandle<geo::Geometry> geom;
    art::ServiceHandle<nova::dbi::RunHistoryService> RunHist;
    assert(geom->DetId() != novadaq::cnv::kUNKNOWN_DET);

    if(!fEvFlag) GetEvInfo();  // setup hasn't been called yet!

    int dibs[14] = {0}, dibsUD[14] = {0}, dibsRH[14] = {0};

    for(int i=0;i<fNumDiblocks;i++) dibsUD[i]=1;

    if(RunHist->GetDataType() == 2 && RunHist->NDiBlocks()>0) { // if it's data
      for(int idib = RunHist->NDiBlocks()-1; idib >= 0; idib--) { // start from back; this is downstream part
        nova::dbi::RunHistory::DiBlock dib = RunHist->GetDiBlock(idib, false);
        if (RunHist->IsDiBlockFullyInstrumented(dib.num)) dibsRH[dib.num-1] = 1;
      }
    }
    else for(int i=0;i<fNumDiblocks;i++) dibsRH[i]=1;

    for(int i=0;i<fNumDiblocks;i++) {
      if(dibsUD[i]==1&&dibsRH[i]==1) dibs[i]=1;
    }

    // look from front to back for upstream part
    for (int i = 0; i<fNumDiblocks; i++) {
      if(dibs[i] == 1 && fFrontUpstream == -5) fFrontUpstream = i; // found the front
      if(dibs[i] == 1 && fFrontUpstream != -5) continue; // still in good part, keep going
      if(dibs[i] == 0 && fFrontUpstream == -5) continue; // nothing yet, keep looking
      if(dibs[i] == 0 && fFrontUpstream != -5) {  // found the back
        fBackUpstream = i;
        break;
      }
    }
    if(fBackUpstream == -5 && fFrontUpstream != -5) { // if the back is diblock 14
      fBackUpstream = fNumDiblocks;
    }
    if(fBackUpstream == -5 && fFrontUpstream == -5) { // no part of the detector is good
      fBackUpstream = 0;
      fFrontUpstream = 0;
    }

    // look from back to front for upstream part
    for (int i = fNumDiblocks-1; i>=0; i--) {
      if(dibs[i] == 1 && fBackDownstream == -5) fBackDownstream = i+1; // found the back
      if(dibs[i] == 1 && fBackDownstream != -5) continue; // still in good part, keep going
      if(dibs[i] == 0 && fBackDownstream == -5) continue; // nothing yet, keep looking
      if(dibs[i] == 0 && fBackDownstream != -5) {  // found the front
        fFrontDownstream = i+1;
        break;
      }
    }
    if(fFrontDownstream == -5 && fBackDownstream != -5) { // if the front is diblock 0
      fFrontDownstream = 0;
    }
    if(fFrontDownstream == -5 && fBackDownstream == -5) { // no part of the detector is good
      fBackDownstream = 0;
      fFrontDownstream = 0;
    }

    // End search

    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)
      fBackUpstream = 0;
      fFrontUpstream = 0;
      fBackDownstream = 0;
      fFrontDownstream = 0;
    }
    else if(fBackDownstream==fBackUpstream && fFrontDownstream==fFrontUpstream) {  // this is the one contiguous detector case, both searches find same thing
      fBackUpstream = 0;  // this is how we'll know there's only one contiguous block
      fFrontUpstream = 0;
    }
    else if(fBackUpstream < fFrontDownstream) {}  // this is the two (or more, but then middle chunks ignored) contiguous detectors case
    else {
      std::cout << "Problem in instrumented detector logic; this should never happen. " << fBackUpstream << " " << fFrontUpstream << " " << fBackDownstream << " " << fFrontDownstream << std::endl;
      fBackUpstream = 0;
      fFrontUpstream = 0;
      fBackDownstream = 0;
      fFrontDownstream = 0;
    }

    return;
  }

  //......................................................................
  double LiveGeometry::InstrumentedDetFront()
  {
    art::ServiceHandle<geo::Geometry> geom;
    assert(geom->DetId() != novadaq::cnv::kUNKNOWN_DET);
    if(fBackDownstream==-5) SetInstrumentation();
    return std::max(fEdgeBackYZ,fEdgeBackXZ)*(float(fFrontDownstream)/float(fNumDiblocks));
  }

  //......................................................................
  double LiveGeometry::InstrumentedDetBack()
  {
    art::ServiceHandle<geo::Geometry> geom;
    assert(geom->DetId() != novadaq::cnv::kUNKNOWN_DET);
    if(fBackDownstream==-5) SetInstrumentation();
    return std::max(fEdgeBackYZ,fEdgeBackXZ)*(float(fBackDownstream)/float(fNumDiblocks));
  }

  //......................................................................
  double LiveGeometry::InstrumentedDetLength()
  {
    art::ServiceHandle<geo::Geometry> geom;
    assert(geom->DetId() != novadaq::cnv::kUNKNOWN_DET);
    if(fBackDownstream==-5) SetInstrumentation();
    return std::max(fEdgeBackYZ,fEdgeBackXZ)*(float(fBackDownstream-fFrontDownstream)/float(fNumDiblocks));
  }

  //......................................................................
  int LiveGeometry::InstrumentedDetEnds(double &frontDwnstrm, double &backDwnstrm, double &frontUpstrm, double &backUpstrm)
  {
    art::ServiceHandle<geo::Geometry> geom;
    assert(geom->DetId() != novadaq::cnv::kUNKNOWN_DET);
    if(fBackDownstream==-5) SetInstrumentation();

    frontDwnstrm=std::max(fEdgeBackYZ,fEdgeBackXZ)*(float(fFrontDownstream)/float(fNumDiblocks));
    backDwnstrm=std::max(fEdgeBackYZ,fEdgeBackXZ)*(float(fBackDownstream)/float(fNumDiblocks));
    frontUpstrm=std::max(fEdgeBackYZ,fEdgeBackXZ)*(float(fFrontUpstream)/float(fNumDiblocks));
    backUpstrm=std::max(fEdgeBackYZ,fEdgeBackXZ)*(float(fBackUpstream)/float(fNumDiblocks));

    if(backDwnstrm>0&&backUpstrm>0) return 2;  // two contiguous detectors
    else if(backDwnstrm>0&&backUpstrm==0) return 1; // one contiguous detector
    else return 0; // something is wrong, or no good diblocks
  }

  //......................................................................
  bool   LiveGeometry::InMuonCatcher(TVector3 vertex)
  {
    art::ServiceHandle<geo::Geometry> geom;
    if(geom->DetId() != novadaq::cnv::kNEARDET) {
      std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
      return false;
    }
    bool ret = false;
    if(vertex.X() > -fMCEdge && vertex.X() < fMCEdge &&
       vertex.Y() > -fMCEdge && vertex.Y() < fMCTop  &&
       vertex.Z() > fMCFront && vertex.Z() < fMCBack) ret = true;
    return ret;
  }

  //......................................................................
  int    LiveGeometry::DoWeEnterMC(TVector3 vertex, TVector3 dir)
  {
    art::ServiceHandle<geo::Geometry> geom;
    if(geom->DetId() != novadaq::cnv::kNEARDET) {
      std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
      return 0;
    }
    if(dir.Z() < 0 && vertex.Z() < fMCFront) return 0;  // this track couldn't possible reach the muon catcher
    if(dir.Z() > 0 && vertex.Z() > fMCFront) return 0;  // this track couldn't possible enter the main ND
    double xyz[3], dxyz[3];
    xyz[0] = vertex.X();
    xyz[1] = vertex.Y();
    xyz[2] = vertex.Z();
    dxyz[0] = dir.X();
    dxyz[1] = dir.Y();
    dxyz[2] = dir.Z();
    if(!geo::IntersectsBox(xyz,dxyz,-fMCEdge,fMCEdge,-fMCEdge,fMCTop,fMCFront,fMCBack)) {
      if(InMuonCatcher(vertex)) return 1;  // fully contained in muon catcher?
      else return 0; // nope, doesn't intersect the muon catcher
    }
    double ypos = YPositionAtTransition(vertex,dir);
    double xpos = XPositionAtTransition(vertex,dir);
    if(ypos > fMCTop && ypos < fEdgeTopYZ && xpos > fEdgeXMinusXZ && xpos < fEdgeXPlusXZ) return 2;  // intersects but goes through the air gap
    else return 1; // seems like it properly enters the MC
  }

  //......................................................................
  double LiveGeometry::ProjectedDistanceInMC(TVector3 vertex, TVector3 dir)
  {
    art::ServiceHandle<geo::Geometry> geom;
    if(geom->DetId() != novadaq::cnv::kNEARDET) {
      std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
      return 0;
    }

    double distX, distY, distZ;

    if(InMuonCatcher(vertex)) {  // if vertex is in the muon catcher

      if(dir.Z() > 0) distZ = (fMCBack-vertex.Z())/dir.Z();  // going +z project to back
      else if(dir.Z() < 0) distZ = (fMCFront-vertex.Z())/dir.Z(); // going -z project to front
      else {  // dz is zero.  This means something is wrong with the track
        return 0;
      }

      if(dir.Y() > 0) distY = (fMCTop-vertex.Y())/dir.Y();  // going +y project to top
      else if(dir.Y() < 0) distY = (-fMCEdge-vertex.Y())/dir.Y(); // going -y project to bottom
      else distY = 999999; // not moving in Y?  May be 2D track... we can still proceed though

      if(dir.X() > 0) distX = (fMCEdge-vertex.X())/dir.X();  // going +x project +x wall
      else if(dir.X() < 0) distX = (-fMCEdge-vertex.X())/dir.X(); // going -x project to -x wall
      else distX = 999999; // not moving in X? may be 2D track... we can still proceed though

      return std::min(std::min(distX,distY),distZ);
    }

    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)
      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
      if(dir.Z() < 0) return 0; // this can't point into the muon catcher
      double Ytran = YPositionAtTransition(vertex, dir);
      double Xtran = XPositionAtTransition(vertex, dir);
      if(fabs(Xtran) < fMCEdge && Ytran < fMCTop && Ytran > -fMCEdge) { // this enters the muon catcher (but ignores air gap events)
        TVector3 vertex2(Xtran, Ytran, fMCFront+0.01);
        if(!InMuonCatcher(vertex2)) {
          std::cout << "Something weird with muon catcher distance projection in LiveGeometry... email this situation to bays@caltech.edu" << std::endl;
          return 0;
        }
        return ProjectedDistanceInMC(vertex2, dir);
      }
    }

    return 0;  // should only reach this if not near detector or if weird error
  }

  //......................................................................
  int    LiveGeometry::ProjectedCellsInMC(TVector3 vertex, TVector3 dir)
  {
    art::ServiceHandle<geo::Geometry> geom;
    if(geom->DetId() != novadaq::cnv::kNEARDET) {
      std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
      return 0;
    }

    std::vector< OfflineChan > xhits, yhits;

    if(InMuonCatcher(vertex)) {
      TVector3 end = vertex + ProjectedDistanceInMC(vertex,dir)*dir;
      geom->CountCellsOnLineFast(vertex, end, xhits, yhits);
      return xhits.size() + yhits.size();
    }

    else {
      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
      if(dir.Z() < 0) return 0; // this can't point into the muon catcher
      double Ytran = YPositionAtTransition(vertex, dir);
      double Xtran = XPositionAtTransition(vertex, dir);
      if(fabs(Xtran) < fMCEdge && Ytran < fMCTop && Ytran > -fMCEdge) { // this enters the muon catcher (but ignores air gap events)
        TVector3 vertex2(Xtran, Ytran, fMCFront+0.01);
        if(!InMuonCatcher(vertex2)) {
          std::cout << "Something weird with muon catcher distance projection in LiveGeometry... email this situation to bays@caltech.edu" << std::endl;
          return 0;
        }
        return ProjectedCellsInMC(vertex2, dir);
      }
    }

    return 0;
  }

  //......................................................................
  int LiveGeometry::FullNDProjectedCells(TVector3 vertex, TVector3 dir) //combine main ND and muon catcher into one
  {
    art::ServiceHandle<geo::Geometry> geom;
    if(geom->DetId() != novadaq::cnv::kNEARDET) {
      std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
      return 0;
    }

    double xyz[3], dxyz[3];
    xyz[0] = vertex.X();
    xyz[1] = vertex.Y();
    xyz[2] = vertex.Z();
    dxyz[0] = dir.X();
    dxyz[1] = dir.Y();
    dxyz[2] = dir.Z();

    if(IsPointLive(vertex)==false && InMuonCatcher(vertex)==false) return -5;  // outside both
    if(IsPointLive(vertex)==true && DoWeEnterMC(vertex, dir)==false) return ProjectedLiveCellsToEdge(vertex,dir);  // just main ND
    if(IsPointLive(vertex)==true && DoWeEnterMC(vertex, dir)==true) return (ProjectedLiveCellsToEdge(vertex,dir)+ProjectedCellsInMC(vertex,dir)); // main to MC
    if(InMuonCatcher(vertex)==true && (dir.Z()>0 || !(geo::IntersectsBox(xyz,dxyz,fEdgeXMinusXZ,fEdgeXPlusXZ,fEdgeBottomYZ,fEdgeTopYZ,std::min(fEdgeFrontXZ,fEdgeFrontYZ),std::max(fEdgeBackXZ,fEdgeBackYZ))))) { // MC only
      return ProjectedCellsInMC(vertex,dir); }
    else { // MC to main
      TVector3 v2 = vertex + (ProjectedDistanceInMC(vertex,dir))*dir; //jump to boundary with main detector
      for(int i = 0; i<100; i++) {
        if(IsPointLive(v2)) return (ProjectedCellsInMC(vertex,dir) + ProjectedLiveCellsToEdge(v2,dir));
        v2 = v2 + 0.05*dir; //step along path until we're in main detector for sure
      }
      return ProjectedCellsInMC(vertex,dir); // never found main detector; something weird, just use MC part
    }
    return -111; // this should never happen
  }

  //......................................................................
  double LiveGeometry::XPositionAtTransition(TVector3 vertex, TVector3 dir)
  {
    art::ServiceHandle<geo::Geometry> geom;
    if(geom->DetId() != novadaq::cnv::kNEARDET) {
      std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
      return 0;
    }
    if(dir.Z() == 0) return -9999;
    return vertex.X() + dir.X()*(fMCFront-vertex.Z())/dir.Z();
  }

  //......................................................................
  double LiveGeometry::YPositionAtTransition(TVector3 vertex, TVector3 dir)
  {
    art::ServiceHandle<geo::Geometry> geom;
    if(geom->DetId() != novadaq::cnv::kNEARDET) {
      std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
      return 0;
    }
    if(dir.Z() == 0) return -9999;
    return vertex.Y() + dir.Y()*(fMCFront-vertex.Z())/dir.Z();
  }

  //......................................................................
  double LiveGeometry::DistanceToEdgeInMC(TVector3 vertex, int &wall)
  {
    art::ServiceHandle<geo::Geometry> geom;
    if(geom->DetId() != novadaq::cnv::kNEARDET) {
      std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
      return 0;
    }
    wall = 0;
    if(!InMuonCatcher(vertex)) return 0; //NOT IN MUON CATCHER
    double mindist = 9999;

    if(mindist > (fMCEdge+vertex.X())) { // check -X direction
      mindist = fMCEdge+vertex.X();
      wall = 1;
    }

    if(mindist > (fMCEdge-vertex.X())) { // check +X direction
      mindist = fMCEdge-vertex.X();
      wall = 2;
    }

    if(mindist > (fMCEdge+vertex.Y())) { // check -Y direction
      mindist = fMCEdge+vertex.Y();
      wall = 3;
    }

    if(mindist > (fMCTop-vertex.Y())) { // check +Y direction
      mindist = fMCTop-vertex.Y();
      wall = 4;
    }

    if(mindist > (fMCBack-vertex.Z())) { // check +Z direction
      mindist = fMCBack-vertex.Z();
      wall = 6;
    }

    return mindist; // no need to check -Z because that is always live detector
  }

  //......................................................................
  int LiveGeometry::CellsToEdgeInMC(TVector3 vertex, int &wall)
  {
    art::ServiceHandle<geo::Geometry> geom;
    if(geom->DetId() != novadaq::cnv::kNEARDET) {
      std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
      return 0;
    }
    std::vector< OfflineChan > xhits, yhits;
    TVector3 out = vertex;
    double dist = DistanceToEdgeInMC(vertex, wall);
    if(wall == 0) return 0;
    if(wall == 1) out.SetX(vertex.X()+dist);
    if(wall == 2) out.SetX(vertex.X()-dist);
    if(wall == 3) out.SetY(vertex.Y()+dist);
    if(wall == 4) out.SetY(vertex.Y()-dist);
    if(wall == 6) out.SetZ(vertex.Z()+dist);
    geom->CountCellsOnLineFast(vertex,out,xhits,yhits);
    return xhits.size()+yhits.size();
  }

  //......................................................................
  double LiveGeometry::ProjectedAirDist(TVector3 vertex, TVector3 dir)
  {
    art::ServiceHandle<geo::Geometry> geom;
    if(geom->DetId() != novadaq::cnv::kNEARDET) {
      std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
      return 0;
    }

    if(DoWeEnterMC(vertex, dir) != 2) return 0;

    TVector3 posattrans;  // position along line at transition between main ND and muon catcher
    posattrans.SetX(vertex.X() + dir.X()*(fMCFront-vertex.Z())/dir.Z());
    posattrans.SetY(vertex.Y() + dir.Y()*(fMCFront-vertex.Z())/dir.Z());
    posattrans.SetZ(fMCFront);

    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.
    posatmctop.SetX(vertex.X() + dir.X()*(fMCTop-vertex.Y())/dir.Y());
    posatmctop.SetY(fMCTop);
    posatmctop.SetZ(vertex.Z() + dir.Z()*(fMCTop-vertex.Y())/dir.Y());

    if(dir.Z() > 0 && vertex.Z() > fMCFront) posattrans = vertex; // these two are for the odd case of a vertex actually in the air gap
    if(dir.Z() < 0 && vertex.Z() > fMCFront) posatmctop = vertex;

    return (posatmctop-posattrans).Mag(); // distance between these points is the distance through the air
  }

  //......................................................................
  double LiveGeometry::ProjectedSteelDist(TVector3 vertex, TVector3 dir)
  {
    art::ServiceHandle<geo::Geometry> geom;
    if(geom->DetId() != novadaq::cnv::kNEARDET) {
      std::cout << "You are trying to use a near detector only function, but this isn't a near detector file!" << std::endl;
      return 0;
    }

    if(DoWeEnterMC(vertex,dir) == 0) return 0;

    double p1[3] = {0}, p2[3], linedist;  // start and end points of a line

    if(InMuonCatcher(vertex)) {
      p1[0] = vertex.X();
      p1[1] = vertex.Y();
      p1[2] = vertex.Z();
    }

    else {
      p1[0] = XPositionAtTransition(vertex,dir);  // must have started in main ND; let's start counting steel at front of muon catcher
      p1[1] = YPositionAtTransition(vertex,dir);
      p1[2] = fMCFront;
    }

    linedist = std::sqrt((8*fMCEdge*fMCEdge) + ((fMCBack-fMCFront)*(fMCBack-fMCFront))) + 10; //longest possible distance through muon catcher (including possibility of air gap)

    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
    p2[1] = p1[1] + linedist*dir.Y();
    p2[2] = p1[2] + linedist*dir.Z();

    std::vector< double > ds;
    std::vector < const TGeoMaterial * > mat;

    geom->MaterialsBetweenPoints(p1,p2,ds,mat);

    double steeldist = 0;

    for(unsigned int i = 0; i < mat.size(); i++) {
      if(strncmp("Steel",mat[i]->GetName(),5) == 0){
        steeldist+=ds[i];
      }
    }

    return steeldist;

  }

  //......................................................................
  int LiveGeometry::NBadDCMBC()
  {
    if(!fEvFlag) GetEvInfo();  // setup hasn't been called yet!
    return fNBadChan;
  }

  //......................................................................
  int LiveGeometry::NDropouts()
  {
    if(!fEvFlag) GetEvInfo();  // setup hasn't been called yet!
    return fNDropouts;
  }

  //......................................................................
  int LiveGeometry::TotalDropouts()
  {
    if(!fEvFlag) GetEvInfo();  // setup hasn't been called yet!
    return fTotalDropouts;
  }

  //......................................................................
  int LiveGeometry::DCMStatus(unsigned int dib, unsigned int dcm)
  {
    if(!fEvFlag) GetEvInfo();  // setup hasn't been called yet!
    if(fDCMs[dcm-1][dib-1]==2) return 0;
    if(fDCMs[dcm-1][dib-1]==0) return 2;
    if(fDCMs[dcm-1][dib-1]==1) return 1;
    return -1;
  }

  //......................................................................
  void LiveGeometry::GetEvInfo()
  {
    if(!fSetup) SetupLiveGeo();  // setup hasn't been called yet!
    if(fParams.CheckBadChannels()) GetBCInfo();    // check bad channels every event in case of dropped DCMs
    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.
    fEvFlag = true;
    fTotalDropouts+=fNDropouts;
    fNEvents++;
  }

  double LiveGeometry::DistToClosestBadRegion(TVector3 vertex)
  {
    if(!fEvFlag) GetEvInfo();  // setup hasn't been called yet!

    art::ServiceHandle<geo::Geometry> geom;
    double mindist = 999999, tdist;

    if(geom->DetId() == novadaq::cnv::kNEARDET && InMuonCatcher(vertex)) {
      int wall;
      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
    }

    if(IsPointLive(vertex) == false) return 0; // already in a bad region

    // first get distances to detector edges; trivial
    tdist = vertex.Z() - std::min(fEdgeFrontXZ,fEdgeFrontYZ); // dist to front
    if(tdist < mindist) mindist = tdist;

    tdist = std::max(fEdgeBackXZ,fEdgeBackYZ) - vertex.Z(); // dist to back
    if(tdist < mindist) mindist = tdist;

    tdist = vertex.X() - fEdgeXMinusXZ;
    if(tdist < mindist) mindist = tdist;

    tdist = fEdgeXPlusXZ - vertex.X();
    if(tdist < mindist) mindist = tdist;

    tdist = vertex.Y() - fEdgeBottomYZ;
    if(tdist < mindist) mindist = tdist;

    tdist = fEdgeTopYZ = vertex.Y();
    if(tdist < mindist) mindist = tdist;

    // check bad boxes in XZ view
    for(unsigned int i = 0; i < fBadBoxXZ.size(); i++) {

      if(fBadBoxXZ[i].Z1  <= vertex.Z() && fBadBoxXZ[i].Z2  >= vertex.Z()) { // lined up in Z
        tdist = std::min(fabs(vertex.X()-fBadBoxXZ[i].XY1),fabs(vertex.X()-fBadBoxXZ[i].XY2));
        if(tdist < mindist) mindist = tdist;
      }

      else if(fBadBoxXZ[i].XY1  <= vertex.X() && fBadBoxXZ[i].XY2  >= vertex.X()) { // lined up in X
        tdist = std::min(fabs(vertex.Z()-fBadBoxXZ[i].Z1),fabs(vertex.Z()-fBadBoxXZ[i].Z2));
        if(tdist < mindist) mindist = tdist;
      }

      else { // not lined up in either x or z; closest will be one of the four corners
        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))));
        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))));
        tdist = std::min(td1,td2);
        if(tdist < mindist) mindist = tdist;
      }
    }

    // check bad boxes in YZ view
    for(unsigned int i = 0; i < fBadBoxYZ.size(); i++) {

      if(fBadBoxYZ[i].Z1  <= vertex.Z() && fBadBoxYZ[i].Z2  >= vertex.Z()) { // lined up in Z
        tdist = std::min(fabs(vertex.Y()-fBadBoxYZ[i].XY1),fabs(vertex.Y()-fBadBoxYZ[i].XY2));
        if(tdist < mindist) mindist = tdist;
      }

      else if(fBadBoxYZ[i].XY1  <= vertex.Y() && fBadBoxYZ[i].XY2  >= vertex.Y()) { // lined up in X
        tdist = std::min(fabs(vertex.Z()-fBadBoxYZ[i].Z1),fabs(vertex.Z()-fBadBoxYZ[i].Z2));
        if(tdist < mindist) mindist = tdist;
      }

      else { // not lined up in either x or z; closest will be one of the four corners
        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))));
        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))));
        tdist = std::min(td1,td2);
        if(tdist < mindist) mindist = tdist;
      }
    }

    return mindist;

  }

  DEFINE_ART_SERVICE(LiveGeometry)
} // end namespace geo
////////////////////////////////////////////////////////////////////////