PlaneOutline.cxx

Go to the documentation of this file.

#include "DataUtil/PlaneOutline.h"
#include <cmath>
#include <list>
#include <cassert>
#include <algorithm>
//#include <iostream>

#include "TPolyLine.h"
#include "TMath.h"
#include "TH2.h"
#include "TExec.h"
#include "TList.h"

#include "Conventions/Munits.h"

std::vector< std::pair<float, float> > PlaneOutline::fNearPartialV  = PlaneOutline::InitNearPartialV();
std::vector< std::pair<float, float> > PlaneOutline::fNearPartialU  = PlaneOutline::InitNearPartialU();
std::vector< std::pair<float, float> > PlaneOutline::fNearFullV     = PlaneOutline::InitNearFullV();
std::vector< std::pair<float, float> > PlaneOutline::fNearFullU     = PlaneOutline::InitNearFullU();
std::vector< std::pair<float, float> > PlaneOutline::fNearFullCoilV = PlaneOutline::InitNearFullCoilV();
std::vector< std::pair<float, float> > PlaneOutline::fNearFullCoilU = PlaneOutline::InitNearFullCoilU();
std::vector< std::pair<float, float> > PlaneOutline::fNearSteel     = PlaneOutline::InitNearSteel();

std::vector< std::pair<float, float> > PlaneOutline::fFarV          = PlaneOutline::InitFarV();
std::vector< std::pair<float, float> > PlaneOutline::fFarU          = PlaneOutline::InitFarU();
std::vector< std::pair<float, float> > PlaneOutline::fFarCoilV      = PlaneOutline::InitFarCoilV();
std::vector< std::pair<float, float> > PlaneOutline::fFarCoilU      = PlaneOutline::InitFarCoilU();
std::vector< std::pair<float, float> > PlaneOutline::fFarSteel      = PlaneOutline::InitFarSteel();

PlaneOutline::PlaneOutline() {}

std::vector< std::pair<float,float> > PlaneOutline::InitFarCoilU(){
  // coil hole of full u plane
  using Munits::mm;
  std::vector< std::pair<float,float> > v;
  // input in u,v coordinates
  v.push_back( std::pair<float,float>(41*mm, -300*mm) );
  v.push_back( std::pair<float,float>(243*mm, -92*mm) );
  v.push_back( std::pair<float,float>(243*mm, 92*mm) );
  v.push_back( std::pair<float,float>(41*mm, 300*mm) );
  v.push_back( std::pair<float,float>(-41*mm, 300*mm) );
  v.push_back( std::pair<float,float>(-243*mm, 92*mm) );
  v.push_back( std::pair<float,float>(-243*mm, -92*mm) );
  v.push_back( std::pair<float,float>(-41*mm, -300*mm) );
  v.push_back( std::pair<float,float>(41*mm, -300*mm) );
  
  for (unsigned int i=0; i<v.size(); i++){    
    float x=sqrt(1.0/2.0)*(v[i].first - v[i].second);
    float y=sqrt(1.0/2.0)*(v[i].first + v[i].second);
    v[i].first = x;
    v[i].second = y;
  }
  

  return v;
}

std::vector< std::pair<float,float> > PlaneOutline::InitFarCoilV(){
  // coil hole of full u plane
  using Munits::mm;
  std::vector< std::pair<float,float> > v;
  // input in u,v coordinates

  v.push_back( std::pair<float,float>(-300*mm, 41*mm) );
  v.push_back( std::pair<float,float>(-92*mm, 243*mm) );
  v.push_back( std::pair<float,float>(92*mm, 243*mm) );
  v.push_back( std::pair<float,float>(300*mm, 41*mm) );
  v.push_back( std::pair<float,float>(300*mm, -41*mm) );
  v.push_back( std::pair<float,float>(92*mm, -243*mm) );
  v.push_back( std::pair<float,float>(-92*mm, -243*mm) );
  v.push_back( std::pair<float,float>(-300*mm, -41*mm) );
  v.push_back( std::pair<float,float>(-300*mm, 41*mm) );
  
  for (unsigned int i=0; i<v.size(); i++){    
    float x=sqrt(1.0/2.0)*(v[i].first - v[i].second);
    float y=sqrt(1.0/2.0)*(v[i].first + v[i].second);
    v[i].first = x;
    v[i].second = y;
  }
  

  return v;
}


std::vector< std::pair<float,float> > PlaneOutline::InitNearFullCoilU(){
  // coil hole of full u plane
  using Munits::inch;
  std::vector< std::pair<float,float> > v;
  v.push_back( std::pair<float,float>(-12.22*inch, 12.56*inch) );
  v.push_back( std::pair<float,float>(2.89*inch, 15.10*inch) );
  v.push_back( std::pair<float,float>(14.93*inch, 3.06*inch) );
  v.push_back( std::pair<float,float>(12.22*inch, -12.22*inch) );
  v.push_back( std::pair<float,float>(-3.06*inch, -14.93*inch) );
  v.push_back( std::pair<float,float>(-14.93*inch, -3.22*inch) );
  v.push_back( std::pair<float,float>(-12.22*inch, 12.56*inch) );

  return v;
}


std::vector< std::pair<float,float> > PlaneOutline::InitNearFullCoilV(){
  // coil hole of full v plane
  using Munits::inch;
  std::vector< std::pair<float,float> > v;
  v.push_back( std::pair<float,float>(-3.13*inch, 14.60*inch) );
  v.push_back( std::pair<float,float>(12.16*inch, 12.34*inch) );
  v.push_back( std::pair<float,float>(14.25*inch, -3.30*inch) );
  v.push_back( std::pair<float,float>(1.74*inch, -14.77*inch) );
  v.push_back( std::pair<float,float>(-12.69*inch, -12.34*inch) );
  v.push_back( std::pair<float,float>(-14.77*inch, 2.78*inch) );
  v.push_back( std::pair<float,float>(-3.13*inch, 14.60*inch) );

  return v;
}


std::vector< std::pair<float,float> > PlaneOutline::InitNearPartialV(){
  // partial v plane
  using Munits::inch;
  std::vector< std::pair<float,float> > v;
  v.push_back( std::pair<float,float>(-1.96*inch, 12.72*inch) );
  v.push_back( std::pair<float,float>(58.00*inch, 73.15*inch) );
  v.push_back( std::pair<float,float>(103.75*inch, 27.89*inch) );
  v.push_back( std::pair<float,float>(108.15*inch, 32.04*inch) );
  v.push_back( std::pair<float,float>(108.15*inch, -32.07*inch) );
  v.push_back( std::pair<float,float>(62.90*inch, -77.58*inch) );
  v.push_back( std::pair<float,float>(-1.79*inch, -12.50*inch) );
  v.push_back( std::pair<float,float>(10.05*inch, 0.0*inch) );
  v.push_back( std::pair<float,float>(-1.96*inch, 12.72*inch) );

  return v;
}

std::vector< std::pair<float,float> > PlaneOutline::InitNearPartialU(){
  // partial u plane
  using Munits::inch;
  std::vector< std::pair<float,float> > v;
  v.push_back( std::pair<float,float>(-1.79*inch, 12.30*inch) );
  v.push_back( std::pair<float,float>(62.92*inch, 77.46*inch) );
  v.push_back( std::pair<float,float>(108.36*inch, 32.00*inch) );
  v.push_back( std::pair<float,float>(108.36*inch, -32.03*inch) );
  v.push_back( std::pair<float,float>(103.88*inch, -27.55*inch) );
  v.push_back( std::pair<float,float>(58.45*inch, -72.98*inch) );
  v.push_back( std::pair<float,float>(-1.79*inch, -12.54*inch) );
  v.push_back( std::pair<float,float>(10.09*inch, 0.0*inch) );
  v.push_back( std::pair<float,float>(-1.79*inch, 12.30*inch) );

  return v;
}

std::vector< std::pair<float,float> > PlaneOutline::InitNearFullV(){
  // full v plane
  using Munits::inch;
  std::vector< std::pair<float,float> > v;
  v.push_back( std::pair<float,float>(-67.52*inch, 4.93*inch) );
  v.push_back( std::pair<float,float>(4.17*inch, 76.10*inch) );
  v.push_back( std::pair<float,float>(72.18*inch, 72.18*inch) );
  v.push_back( std::pair<float,float>(113.11*inch, 31.52*inch) );
  v.push_back( std::pair<float,float>(113.11*inch, -31.80*inch) );
  v.push_back( std::pair<float,float>(71.66*inch, -73.77*inch) );
  v.push_back( std::pair<float,float>(9.63*inch, -72.45*inch) );
  v.push_back( std::pair<float,float>(-67.52*inch, 4.93*inch) );

  return v;
}

std::vector< std::pair<float,float> > PlaneOutline::InitNearFullU(){
  // full u plane
  using Munits::inch;
  std::vector< std::pair<float,float> > v;
  //  v.push_back( std::pair<float,float>(-68.22*inch, 5.36*inch) );
  v.push_back( std::pair<float,float>(-68.22*inch, -5.36*inch) );
  v.push_back( std::pair<float,float>(9.67*inch, 73.44*inch) );
  v.push_back( std::pair<float,float>(72.02*inch, 73.44*inch) );
  v.push_back( std::pair<float,float>(114.28*inch, 32.08*inch) );
  v.push_back( std::pair<float,float>(114.28*inch, -32.08*inch) );
  v.push_back( std::pair<float,float>(73.04*inch, -73.59*inch) );
  v.push_back( std::pair<float,float>(4.07*inch, -76.88*inch) );
  //  v.push_back( std::pair<float,float>(-68.22*inch, 5.36*inch) );
  v.push_back( std::pair<float,float>(-68.22*inch, -5.36*inch) );

  return v;
}

std::vector< std::pair<float,float> > PlaneOutline::InitNearSteel(){
  // near steel plane
  using Munits::inch;
  std::vector< std::pair<float,float> > v;

  // Steel pln (coords about the center, *not* MINOS coords)
  v.push_back( std::pair<float,float>(-121.43*inch,   0.47*inch) );
  v.push_back( std::pair<float,float>(-121.43*inch,   9.28*inch) );
  v.push_back( std::pair<float,float>( -69.80*inch,  60.91*inch) );
  v.push_back( std::pair<float,float>( -69.80*inch,  75.04*inch) );
  v.push_back( std::pair<float,float>(  69.80*inch,  75.04*inch) );
  v.push_back( std::pair<float,float>(  69.80*inch,  60.91*inch) );
  v.push_back( std::pair<float,float>( 121.43*inch,   9.28*inch) );
  v.push_back( std::pair<float,float>( 121.43*inch,   0.47*inch) );
  v.push_back( std::pair<float,float>( 120.25*inch,   0.47*inch) );
  v.push_back( std::pair<float,float>( 117.30*inch,   0.47*inch) );
  v.push_back( std::pair<float,float>( 110.16*inch,   3.42*inch) );
  v.push_back( std::pair<float,float>(  95.24*inch,  -1.48*inch) );
  v.push_back( std::pair<float,float>(  95.24*inch, -35.47*inch) );
  v.push_back( std::pair<float,float>(  69.80*inch, -60.91*inch) );
  v.push_back( std::pair<float,float>(  69.80*inch, -75.04*inch) );
  v.push_back( std::pair<float,float>( -69.80*inch, -75.04*inch) );
  v.push_back( std::pair<float,float>( -69.80*inch, -60.91*inch) );
  v.push_back( std::pair<float,float>( -95.24*inch, -35.47*inch) );
  v.push_back( std::pair<float,float>( -95.24*inch,  -1.48*inch) );
  v.push_back( std::pair<float,float>(-110.24*inch,   3.42*inch) );
  v.push_back( std::pair<float,float>(-115.14*inch,   3.42*inch) );
  v.push_back( std::pair<float,float>(-120.25*inch,   0.47*inch) );
  v.push_back( std::pair<float,float>(-121.43*inch,   0.47*inch) );

  // Correct to MINOS coords
  const double xOffset = +0.5578;
  for (unsigned int i = 0; i < v.size(); i++){
    v[i].first += xOffset;
  }

  return v;
}

std::vector< std::pair<float,float> > PlaneOutline::InitFarSteel(){
  // far steel plane
  using Munits::cm;
  std::vector< std::pair<float,float> > v;

  v.push_back( std::pair<float,float>(  143.7724*cm, -421.0036*cm) );
  v.push_back( std::pair<float,float>(  139.4268*cm, -421.0036*cm) );
  v.push_back( std::pair<float,float>(  116.8712*cm, -398.8619*cm) );
  v.push_back( std::pair<float,float>( -116.9622*cm, -398.8619*cm) );
  v.push_back( std::pair<float,float>( -139.1040*cm, -421.0036*cm) );
  v.push_back( std::pair<float,float>( -143.4496*cm, -421.0036*cm) );
  v.push_back( std::pair<float,float>( -421.5665*cm, -142.8867*cm) );
  v.push_back( std::pair<float,float>( -421.5665*cm, -137.7134*cm) );
  v.push_back( std::pair<float,float>( -399.8386*cm, -115.9855*cm) );
  v.push_back( std::pair<float,float>( -399.8386*cm,   84.7388*cm) );
  v.push_back( std::pair<float,float>( -406.6674*cm,   94.2577*cm) );
  v.push_back( std::pair<float,float>( -429.8438*cm,  101.9141*cm) );
  v.push_back( std::pair<float,float>( -439.9834*cm,  101.9141*cm) );
  v.push_back( std::pair<float,float>( -457.3658*cm,   94.2577*cm) );
  v.push_back( std::pair<float,float>( -457.3658*cm,   94.2577*cm) );
  v.push_back( std::pair<float,float>( -457.3658*cm,  109.3637*cm) );
  v.push_back( std::pair<float,float>( -143.4496*cm,  423.2799*cm) );
  v.push_back( std::pair<float,float>( -139.1040*cm,  423.2799*cm) );
  v.push_back( std::pair<float,float>( -116.9622*cm,  401.1382*cm) );
  v.push_back( std::pair<float,float>(  117.2850*cm,  401.1382*cm) );
  v.push_back( std::pair<float,float>(  139.4268*cm,  423.2799*cm) );
  v.push_back( std::pair<float,float>(  143.7724*cm,  423.2799*cm) );
  v.push_back( std::pair<float,float>(  457.6886*cm,  109.3637*cm) );
  v.push_back( std::pair<float,float>(  457.6886*cm,   94.2577*cm) );
  v.push_back( std::pair<float,float>(  454.7915*cm,   94.2577*cm) );
  v.push_back( std::pair<float,float>(  447.1350*cm,  101.5003*cm) );
  v.push_back( std::pair<float,float>(  429.1319*cm,  101.5003*cm) );
  v.push_back( std::pair<float,float>(  406.9902*cm,   94.0507*cm) );
  v.push_back( std::pair<float,float>(  400.7823*cm,   87.6358*cm) );
  v.push_back( std::pair<float,float>(  400.1614*cm, -115.9855*cm) );
  v.push_back( std::pair<float,float>(  421.8893*cm, -138.3342*cm) );
  v.push_back( std::pair<float,float>(  421.8893*cm, -142.8876*cm) );
  v.push_back( std::pair<float,float>(  143.7724*cm, -421.0036*cm) );

  return v;
} 

std::vector< std::pair<float,float> > PlaneOutline::InitFarU(){
  // full u plane
  const float d=4.0*Munits::meter;
  const float L=2.0*d*tan(TMath::Pi()/8.0);
  std::vector< std::pair<float,float> > v;
  v.push_back( std::pair<float,float>(-L/2.0, d) );
  v.push_back( std::pair<float,float>(L/2.0, d) );
  v.push_back( std::pair<float,float>(d, L/2.0) );
  v.push_back( std::pair<float,float>(d, -L/2.0) );
  v.push_back( std::pair<float,float>(L/2.0, -d) );
  v.push_back( std::pair<float,float>(-L/2.0, -d) );
  v.push_back( std::pair<float,float>(-d, -L/2.0) );
  v.push_back( std::pair<float,float>(-d, L/2.0) );
  v.push_back( std::pair<float,float>(-L/2.0, d) );

  return v;
}

std::vector< std::pair<float,float> > PlaneOutline::InitFarV(){
  // full u plane
  const float d=4.0*Munits::meter;
  const float L=2.0*d*tan(TMath::Pi()/8.0);
  std::vector< std::pair<float,float> > v;
  v.push_back( std::pair<float,float>(-L/2.0, d) );
  v.push_back( std::pair<float,float>(L/2.0, d) );
  v.push_back( std::pair<float,float>(d, L/2.0) );
  v.push_back( std::pair<float,float>(d, -L/2.0) );
  v.push_back( std::pair<float,float>(L/2.0, -d) );
  v.push_back( std::pair<float,float>(-L/2.0, -d) );
  v.push_back( std::pair<float,float>(-d, -L/2.0) );
  v.push_back( std::pair<float,float>(-d, L/2.0) );
  v.push_back( std::pair<float,float>(-L/2.0, d) );

  return v;
}



bool PlaneOutline::IsInside(float xp, float yp, 
                            const std::vector< std::pair<float,float> >& v) const
{
  // Is (xp,yp) inside the plane outline represented by v?


  // taken from TMath::IsInside and modified to use my data structure
  float xint;
  int inter = 0;
  int np=v.size(); 
   for (int i=0;i<np-1;i++) {
      if (v[i].second == v[i+1].second) continue;
      if (yp <= v[i].second && yp <= v[i+1].second) continue;
      if (v[i].second < yp && v[i+1].second < yp) continue;
      xint = v[i].first 
        + (yp-v[i].second)*(v[i+1].first-v[i].first)/(v[i+1].second-v[i].second);
      if (xp < xint) inter++;
   }
   if (inter%2) return true;
   return false;
}

bool PlaneOutline::IsInsideSteel(float x, float y, Detector::Detector_t det) const {
  // Is the (x,y) coord inside the outline of a steel plane in detector det (near/far)?

  bool result = false;
  if (det == Detector::kNear) result = IsInside(x,y,fNearSteel);
  if (det == Detector::kFar)  result = IsInside(x,y,fFarSteel);

  return result;
}

bool PlaneOutline::IsInside(float x, float y, PlaneView::PlaneView_t pv, 
                            PlaneCoverage::PlaneCoverage_t pc, bool bCoilOutside) const {
  // Is 'x,y' inside the outline of the plane type represented by 'pv' and 'pc'?

  bool result=false;
  if( pv==PlaneView::kU ){
    switch(pc){
    case PlaneCoverage::kNearPartial:
      result = IsInside(x,y,fNearPartialU);
      break;
    case PlaneCoverage::kNearFull:
      result = IsInside(x,y,fNearFullU);
      if(result && bCoilOutside){ result = !IsInside(x,y,fNearFullCoilU); } // coil hole region
      break;
    case PlaneCoverage::kComplete:
      result = IsInside(x,y,fFarU);
      if(result && bCoilOutside){ result = !IsInside(x,y,fFarCoilU); } // coil hole region      
      break;
    default:
      break;
    }
  }
  else if (pv==PlaneView::kV){
    switch(pc){
    case PlaneCoverage::kNearPartial:
      result = IsInside(x,y,fNearPartialV);
      break;
    case PlaneCoverage::kNearFull:
      result = IsInside(x,y,fNearFullV);
      if(result && bCoilOutside){ result = !IsInside(x,y,fNearFullCoilV); } // coil hole region
      break;
    case PlaneCoverage::kComplete:
      result = IsInside(x,y,fFarV);
      if(result && bCoilOutside){ result = !IsInside(x,y,fFarCoilV); } // coil hole region      
      break;
    default:
      break;      
    }
  }
  return result;
 
}


bool PlaneOutline::DistanceToEdge(float x, float y, PlaneView::PlaneView_t pv,
                                  PlaneCoverage::PlaneCoverage_t pc, 
                                  float& dist, float& xedge, float& yedge) const {
  // compute the distance from 'x,y' to the closest point on the edge 
  // of the plane denoted by the 'pc' and 'pv' arguments
  //
  // if a distance can be computed then set 'dist' to the minimum distance
  // and 'xedge', 'yedge' to the closest point on that edge and return true
  // otherwise return false

 
  bool result=false;
  if( pv==PlaneView::kU ){
    float dist2, x2,y2;
    switch(pc){
    case PlaneCoverage::kNearPartial:
      dist = DistanceToEdge(x,y,fNearPartialU,xedge,yedge);
      result=true;
      break;
    case PlaneCoverage::kNearFull:
      dist = DistanceToEdge(x,y,fNearFullU,xedge,yedge);
      dist2 = DistanceToEdge(x,y,fNearFullCoilU,x2,y2);
      if(dist2<dist) {yedge=y2; xedge=x2; dist=dist2;}
      result=true;
      break;
    case PlaneCoverage::kComplete:
      dist = DistanceToEdge(x,y,fFarU,xedge,yedge);
      dist2 = DistanceToEdge(x,y,fFarCoilU,x2,y2);
      if(dist2<dist) {yedge=y2; xedge=x2; dist=dist2;}
      result=true;
      break;
    default:
      break;
    }
  }
  else if (pv==PlaneView::kV){
    float dist2, x2, y2;
    switch(pc){
    case PlaneCoverage::kNearPartial:
      dist = DistanceToEdge(x,y,fNearPartialV,xedge,yedge);
      result=true;
      break;
    case PlaneCoverage::kNearFull:
      dist = DistanceToEdge(x,y,fNearFullV,xedge,yedge);
      dist2 = DistanceToEdge(x,y,fNearFullCoilV,x2,y2);
      if(dist2<dist) {yedge=y2; xedge=x2; dist=dist2;}
      result=true;
      break;
    case PlaneCoverage::kComplete:
      dist = DistanceToEdge(x,y,fFarV,xedge,yedge);
      dist2 = DistanceToEdge(x,y,fFarCoilV,x2,y2);
      if(dist2<dist) {yedge=y2; xedge=x2; dist=dist2;}
      result=true;
      break;
    default:
      break;      
    }
  }
  return result; 
}

bool PlaneOutline::DistanceToOuterEdge(float x, float y, PlaneView::PlaneView_t pv,
                                  PlaneCoverage::PlaneCoverage_t pc, 
                                  float& dist, float& xedge, float& yedge) const {
  // compute the distance from 'x,y' to the closest point on the edge 
  // of the plane denoted by the 'pc' and 'pv' arguments
  //
  // if a distance can be computed then set 'dist' to the minimum distance
  // and 'xedge', 'yedge' to the closest point on that edge and return true
  // otherwise return false

 
  bool result=false;
  if( pv==PlaneView::kU ){
    switch(pc){
    case PlaneCoverage::kNearPartial:
      dist = DistanceToEdge(x,y,fNearPartialU,xedge,yedge);
      result=true;
      break;
    case PlaneCoverage::kNearFull:
      dist = DistanceToEdge(x,y,fNearFullU,xedge,yedge);
      result=true;
      break;
    case PlaneCoverage::kComplete:
      dist = DistanceToEdge(x,y,fFarU,xedge,yedge);
      result=true;
      break;
    default:
      break;
    }
  }
  else if (pv==PlaneView::kV){
    switch(pc){
    case PlaneCoverage::kNearPartial:
      dist = DistanceToEdge(x,y,fNearPartialV,xedge,yedge);
      result=true;
      break;
    case PlaneCoverage::kNearFull:
      dist = DistanceToEdge(x,y,fNearFullV,xedge,yedge);
      result=true;
      break;
    case PlaneCoverage::kComplete:
      dist = DistanceToEdge(x,y,fFarV,xedge,yedge);
      result=true;
      break;
    default:
      break;      
    }
  }
  return result; 
}


bool PlaneOutline::DistanceToInnerEdge(float x, float y, PlaneView::PlaneView_t pv,
                                  PlaneCoverage::PlaneCoverage_t pc, 
                                  float& dist, float& xedge, float& yedge) const {
  // compute the distance from 'x,y' to the closest point on the edge 
  // of the plane denoted by the 'pc' and 'pv' arguments
  //
  // if a distance can be computed then set 'dist' to the minimum distance
  // and 'xedge', 'yedge' to the closest point on that edge and return true
  // otherwise return false

 
  bool result=false;
  if( pv==PlaneView::kU ){
    switch(pc){
    case PlaneCoverage::kNearFull:
      dist = DistanceToEdge(x,y,fNearFullCoilU,xedge,yedge);
      result=true;
      break;
    case PlaneCoverage::kComplete:
      dist = DistanceToEdge(x,y,fFarCoilU,xedge,yedge);
      result=true;
      break;
    default:
      break;
    }
  }
  else if (pv==PlaneView::kV){
    switch(pc){
    case PlaneCoverage::kNearFull:
      dist = DistanceToEdge(x,y,fNearFullCoilV,xedge,yedge);
      result=true;
      break;
    case PlaneCoverage::kComplete:
      dist = DistanceToEdge(x,y,fFarCoilV,xedge,yedge);
      result=true;
      break;
    default:
      break;      
    }
  }
  return result; 
}


float PlaneOutline::DistanceToEdge(float x, float y, 
                                   const std::vector< std::pair<float, float> >& v,
                                   float& xedge, float& yedge) const{

  // compute the distance from x,y to the closest point on the edge 
  // of the outline represented in v
  // works both inside and outside the 
  float save_dist = 0., save_x = 0., save_y = 0.;
  for (unsigned int i=1; i<v.size(); i++){
    float dist = ClosestPointToSegment(v[i-1].first, v[i-1].second,
                                    v[i].first, v[i].second,
                                    x,y, xedge, yedge);
    if(i==1){
      save_dist=dist;
      save_x=xedge;
      save_y=yedge;
    }
    else if(dist<save_dist){
      save_dist=dist;
      save_x=xedge;
      save_y=yedge;
    }
  }
  xedge=save_x;
  yedge=save_y;
  return save_dist;

}


float PlaneOutline::ClosestPointToLine(const float& x1, 
                                       const float& y1,
                                       const float& x2,
                                       const float& y2,
                                       const float& xpt,
                                       const float& ypt,
                                       float& x, float& y){
  
  // given an infinite line defined by points (x1,y1) (x2,y2) 
  // and a spatial point (xpt,ypt)
  // report the closest point on the line (x,y) to the spatial point
  // and compute the distance from (xpt,ypt) to (x,y)
  float u= (xpt-x1)*(x2-x1)+(ypt-y1)*(y2-y1);
  u/=sqrt( pow(x2-x1,2) + pow(y2-y1,2));
  x=x1+u*(x2-x1);
  y=y1+u*(y2-y1);
  float dist =sqrt(pow(xpt-x, 2)+pow(ypt-y,2));
  return dist;
}


float PlaneOutline::ClosestPointToSegment(const float& x1, 
                                          const float& y1,
                                          const float& x2,
                                          const float& y2,
                                          const float& xpt,
                                          const float& ypt,
                                          float& x, float& y){
  // given a line segment with end points (x1,y1) (x2,y2) 
  // and a spatial point (xpt,ypt)
  // report the closest point on the line (x,y) to the spatial point
  // and compute the distance from (xpt,ypt) to (x,y)
  //
  // if V(s) is a vector pointing to the point a distance s along the line
  // and W is the vector pointing to (xpt,ypt)
  // then (x,y) can be found by varying s to minimize |V-W|
  // i.e. d|V(s)-W|/ds = 0
  //
  

  float dist=0.0;

  float d = sqrt(pow(x2-x1,2)+pow(y2-y1,2) );
  if(d==0){
    x=x1;y=y1;
  }
  else {
    float nx = (x2-x1)/d;  float ny = (y2-y1)/d;
    float s = ( (xpt - x1)*nx + (ypt - y1)*ny )/(nx*nx + ny*ny);
    
    // check that s refers to a point inside the line segment
    // if not, set x,y to an end point
    if(s<0){
      x=x1;
      y=y1;
    }
    else if(s>d){
      x=x2;
      y=y2;
    }
    else{
      x=x1+s*nx;
      y=y1+s*ny;
    }

  }
  dist=sqrt( pow(x-xpt,2)+pow(y-ypt,2));
  return dist;
}

bool PlaneOutline::DepthInView(float x1, float y1, float x2, float y2, float& d,
                               std::vector< std::pair<float,float> >& v) const{
  d=0.0;
  // find all intersections
  // store their distance along segment defined by (x1,y1),(x2,y2)
  std::list<float> q;
  for (unsigned int i=0; i<(v.size()-1); i++){
    float s,l;
    if(WhereIntersect(x1,y1,x2,y2, v[i].first, v[i].second, 
                      v[i+1].first, v[i+1].second,s,l) ) {
      q.push_back(s);
    }
  }
  if(q.size()%2 !=0 ){
    //    std::cerr<<"["<<x1<<", "<<y1<<"]   ["<<x2<<", "<<y2<<"]"<<std::endl;
    for (unsigned int i=0; i<(v.size()-1); i++){
      //      std::cerr<<"("<<v[i].first<<", "<<v[i].second<<") --> ("
      //               <<v[i+1].first<<", "<<v[i+1].second<<")";
      float s,l;
      if(WhereIntersect(x1,y1,x2,y2, v[i].first, v[i].second, 
                        v[i+1].first, v[i+1].second,s,l) ) {
        //      std::cerr<<" : (s,l) = ("<<s<<", "<<l<<")";
      }
      //      std::cerr<<std::endl;

    }
    std::list<float>::iterator it = q.begin();
    while(it!=q.end()){
      //      std::cerr<<(*it)<<std::endl;
      it++;
    }    
  }
  assert(q.size()%2 == 0); // must be an even number of intersections
  //  std::sort(q.begin(),q.end());
  q.sort();
  std::list<float>::iterator it = q.begin();
  while(it!=q.end()){
    float a=*it;
    it++;
    float b=*it;
    assert( (b-a) > 0);
    d += (b-a);
    it++;
  }
  return true;
}

float PlaneOutline::DepthInView(float d, PlaneView::PlaneView_t myview,
                                PlaneView::PlaneView_t pv, PlaneCoverage::PlaneCoverage_t pc) const {

  const float is2 = 1.0/sqrt(2.0);
  float x1, y1, x2, y2; // will be used later
  float u1, v1, u2, v2; // for bookkeeping, defined here due to scoping
  switch(myview){
  case PlaneView::kX:
    y1=d; y2=d;
    x1=-8.2*Munits::meter; x2=8.2*Munits::meter;
    break;
  case PlaneView::kY:
    x1=d; x2=d;
    y1=-8.2*Munits::meter; y2=8.2*Munits::meter;
    break;
  case PlaneView::kU:
    v1=d; v2=d;
    u1=-8.2*Munits::meter; u2=8.2*Munits::meter;
    x1 = is2*(u1-v1);
    x2 = is2*(u2-v2);
    y1 = is2*(u1+v1);
    y2 = is2*(u2+v2);
    break;
  case PlaneView::kV:
    u1=d; u2=d;
    v1=-8.2*Munits::meter; v2=8.2*Munits::meter;
    x1 = is2*(u1-v1);
    x2 = is2*(u2-v2);
    y1 = is2*(u1+v1);
    y2 = is2*(u2+v2);
    break;
  default:
    return 0;
    break;
  }
  
  bool result=false;
  float dist=0.0;
  float dist2=0.0;
  if( pv==PlaneView::kU ){
    switch(pc){
    case PlaneCoverage::kNearPartial:
      if(DepthInView(x1,y1,x2,y2,dist,fNearPartialU)) result=true;
      break;
    case PlaneCoverage::kNearFull:
      if(DepthInView(x1,y1,x2,y2,dist,fNearFullU)) {
        if(DepthInView(x1,y1,x2,y2,dist2, fNearFullCoilU) ){
          dist-=dist2;
          result=true;
        }
      }
      break;
    case PlaneCoverage::kComplete:
      if(DepthInView(x1,y1,x2,y2,dist,fFarU)) {
        if(DepthInView(x1,y1,x2,y2,dist2, fFarCoilU) ){
          dist-=dist2;
          result=true;
        }
      }
      break;
    default:
      break;
    }
  }
  else if( pv==PlaneView::kV ){
    switch(pc){
    case PlaneCoverage::kNearPartial:
      if(DepthInView(x1,y1,x2,y2,dist,fNearPartialV)) result=true;
      break;
    case PlaneCoverage::kNearFull:
      if(DepthInView(x1,y1,x2,y2,dist,fNearFullV)) {
        if(DepthInView(x1,y1,x2,y2,dist2, fNearFullCoilV) ){
          dist-=dist2;
          result=true;
        }
      }
      break;
    case PlaneCoverage::kComplete:
      if(DepthInView(x1,y1,x2,y2,dist,fFarV)) {
        if(DepthInView(x1,y1,x2,y2,dist2, fFarCoilV) ){
          dist-=dist2;
          result=true;
        }
      }
      break;
    default:
      break;
    }
  }
  return dist; 
}


bool PlaneOutline::WhereIntersect(float x1, float y1, float x2, float y2,
                                  float x3, float y3, float x4, float y4,
                                  float& s, float& l){

  // x = ax + s*bx   =   cx + l*dx
  // y = ay + s*by   =   cy + l*dy

  // solve for s and l the distance along the two infinite lines, defined by
  // (x1,y1),(x2,y2) and (x3,y3),(x4,y4), where the lines intersect
  // 
  // if these two line *segments* intersect then return true

  float b=sqrt( pow(x2-x1,2) + pow(y2-y1,2) );
  float bx = (x2-x1)/b;  float by = (y2-y1)/b;

  float d=sqrt( pow(x4-x3,2) + pow(y4-y3,2) );
  float dx = (x4-x3)/d;  float dy = (y4-y3)/d;

  float ax = x1; float ay = y1;
  float cx = x3; float cy = y3;

  float den=dx*by-bx*dy;
  if(den!=0){
    s = (dx*(cy-ay) + dy*(ax-cx))/den;
    if(fabs(dx)>0) l = ((ax - cx) + s*bx)/dx;
    else  l = ((ay - cy) + s*by)/dy;
  }
  else return false;

  // is the intersection point within the line segments?
  bool result=false;
  float x = ax+bx*s;
  float y = ay+by*s;
  //  std::cout<<"(x,y)=("<<x<<", "<<y<<")"<<std::endl;

  /*
  if( ( (x>=x1 && x<=x2) || (x>=x2 && x<=x1) ) &&
      ( (y>=y1 && y<=y2) || (y>=y2 && y<=y1) ) &&
      ( (x>=x3 && x<=x4) || (x>=x4 && x<=x3) ) &&
      ( (y>=y3 && y<=y4) || (y>=y4 && y<=y3) ) ){

    result=true;
  }
  */

  bool t1 = ( (x>=x1 && x<=x2) || (x>=x2 && x<=x1) || (fabs(x-x2)<1e-6 && fabs(x-x1)<1e-6 && fabs(x1-x2)<1e-6));
  bool t2 = ( (y>=y1 && y<=y2) || (y>=y2 && y<=y1) || (fabs(y-y2)<1e-6 && fabs(y-y1)<1e-6 && fabs(y1-y2)<1e-6));
  bool t3 = ( (x>=x3 && x<=x4) || (x>=x4 && x<=x3) || (fabs(x-x4)<1e-6 && fabs(x-x3)<1e-6 && fabs(x3-x4)<1e-6));
  bool t4 = ( (y>=y3 && y<=y4) || (y>=y4 && y<=y3) || (fabs(y-y4)<1e-6 && fabs(y-y3)<1e-6 && fabs(y3-y4)<1e-6));
  

  /*
  std::cout<<"[";
  if(t1)std::cout<<"yes ";
  else std::cout<<"no ";
  if(t2)std::cout<<"yes ";
  else std::cout<<"no ";
  if(t3)std::cout<<"yes ";
  else std::cout<<"no ";
  if(t4)std::cout<<"yes ";
  else std::cout<<"no ";
  std::cout<<"]"<<std::endl;
  */
  if(t1 && t2 && t3 && t4) result=true;

  return result;
}


void PlaneOutline::GetSteelOutline(Detector::Detector_t det,
                                   TPolyLine*& outer       ){
  std::vector< std::pair<float,float> >* v = (det == Detector::kNear)? &fNearSteel : &fFarSteel;

  outer = new TPolyLine(v->size());
  for (unsigned int i=0; i<v->size(); i++){
    outer->SetPoint(i,(*v)[i].first,(*v)[i].second);
  }

  return;
}

void PlaneOutline::GetOutline(PlaneView::PlaneView_t pv,
                              PlaneCoverage::PlaneCoverage_t pc, 
                              TPolyLine*& outer, TPolyLine*& inner){
  
  std::vector< std::pair<float,float> >* v = &fNearFullU;
  std::vector< std::pair<float,float> >* vin = &fNearFullCoilU;
  bool result=true;
  bool doin=false;
  if( pv==PlaneView::kU ){
    switch(pc){
    case PlaneCoverage::kNearPartial:
      v=&fNearPartialU;
      break;
    case PlaneCoverage::kNearFull:
      //      std::cout<<"hello"<<std::endl;
      v=&fNearFullU;
      vin=&fNearFullCoilU; doin=true;
      break;
    case PlaneCoverage::kComplete:
      v=&fFarU;
      vin=&fFarCoilU; doin=true;
      break;
    default:
      result=false;
      break;
    }
  }
  else if (pv==PlaneView::kV){
    switch(pc){
    case PlaneCoverage::kNearPartial:
      v=&fNearPartialV;
      break;
    case PlaneCoverage::kNearFull:
      v=&fNearFullV;
      vin=&fNearFullCoilV; doin=true;
      break;
    case PlaneCoverage::kComplete:
      v=&fFarV;
      vin=&fFarCoilV; doin=true;
      break;
    default:
      result=false;
      break;    
    }
  }
  
  outer=0;
  inner=0;
  if(result){
    outer = new TPolyLine(v->size());
    for (unsigned int i=0; i<v->size(); i++){
      outer->SetPoint(i,(*v)[i].first,(*v)[i].second);
    }
  }
  if(doin){
    inner = new TPolyLine(vin->size());
    for (unsigned int i=0; i<vin->size(); i++){
      inner->SetPoint(i,(*vin)[i].first,(*vin)[i].second);
    }
  }

}


TH2* PlaneOutline::DepthInViewHist(PlaneCoverage::PlaneCoverage_t pc, 
                                   PlaneView::PlaneView_t pv){
  
  TH2F* h = new TH2F("h", ";distance (m); view direction;", 
                     410, -4.1*Munits::meter, 4.1*Munits::meter, 
                     9, -0.5, 8.5);
  h->GetYaxis()->SetBinLabel(1, ""); 
  h->GetYaxis()->SetBinLabel(2, "X"); 
  h->GetYaxis()->SetBinLabel(3, ""); 
  h->GetYaxis()->SetBinLabel(4, "Y"); 
  h->GetYaxis()->SetBinLabel(5, ""); 
  h->GetYaxis()->SetBinLabel(6, "U"); 
  h->GetYaxis()->SetBinLabel(7, ""); 
  h->GetYaxis()->SetBinLabel(8, "V"); 
  h->GetYaxis()->SetBinLabel(9, ""); 
  
  
 for(int i=1; i<=h->GetNbinsX(); i++){
   float s = h->GetXaxis()->GetBinCenter(i);
   float depth = DepthInView(s,PlaneView::kX,pv,pc);     
   h->Fill(s,1,depth);
   depth = DepthInView(s,PlaneView::kY,pv,pc); 
   h->Fill(s,3,depth);
   depth = DepthInView(s,PlaneView::kU,pv,pc);     
   h->Fill(s,5,depth);
   depth = DepthInView(s,PlaneView::kV,pv,pc);     
   h->Fill(s,7,depth);   
 }
 h->SetDirectory(0);
 return h;
}

TH2* PlaneOutline::GetNDPlanesHist(PlaneView::PlaneView_t view){
  // view here is the direction you are looking
  // so, to see V readout view=U
  if(!(view==PlaneView::kU || view==PlaneView::kV)) return 0;

  TH2* pu = DepthInViewHist(PlaneCoverage::kNearPartial,
                            PlaneView::kU);
  
  TH2* fu = DepthInViewHist(PlaneCoverage::kNearFull,
                            PlaneView::kU);
  
  TH2* pv = DepthInViewHist(PlaneCoverage::kNearPartial,
                            PlaneView::kV);
  
  TH2* fv = DepthInViewHist(PlaneCoverage::kNearFull,
                            PlaneView::kV);
  
  double offset=0.000;
  double lowlim=offset;
  int n=282;
  double pitch=0.0594;
  double highlim=n*pitch;

  
  TH2* h =0;
  if(view==PlaneView::kU) {
    h = new TH2F("po_vvz","Transverse vs Z view - V Planes",
                 n,lowlim,highlim,
                 fv->GetNbinsX(), 
                 fv->GetXaxis()->GetBinLowEdge(1),
                 fv->GetXaxis()->GetBinUpEdge(fv->GetNbinsX()));
  }
  else if(view==PlaneView::kV) {
    h = new TH2F("po_uvz","Transverse vs Z view - U Planes",
                 n,lowlim,highlim,
                 fu->GetNbinsX(), 
                 fu->GetXaxis()->GetBinLowEdge(1),
                 fu->GetXaxis()->GetBinUpEdge(fu->GetNbinsX()));
  }

  // kview: 6=U, 8=V
  int kview=6;
  int pstart=2;
  if(view==PlaneView::kV) { kview=8; pstart=1;}
  
  for(int plane=pstart; plane<n; plane+=2){
    TH1* htemp=0;
    if((plane-1)%5==0){
      // full
      htemp=fu;
      if(kview==8) htemp=fv;
    }
    else if(plane<=120){
      // partial
      htemp=pu;
      if(kview==8) htemp=pv;
    }
    if(htemp!=0){
      for(int j=1; j<=htemp->GetNbinsX(); j++){
        
        float bw=htemp->GetBinContent(j,kview);
        float tpos=htemp->GetXaxis()->GetBinCenter(j);
        int bin = h->GetYaxis()->FindBin(tpos);
        if(bw>0) bw=0.5;
        h->SetBinContent(plane+1,bin,bw);
      }
    }
  }
  
  int coil_bin_high=h->GetYaxis()->FindBin(0.5);
  int coil_bin_low=h->GetYaxis()->FindBin(-0.5);
  
  for(int plane=2; plane<n; plane+=1){
    //  for(int j=coil_bin_low; j<=coil_bin_high; j++){
    h->SetBinContent(plane+1,coil_bin_low,0.95);
    h->SetBinContent(plane+1,coil_bin_high,0.95);
    //  }
    
  }
  
  h->SetMaximum(1.0);
  h->SetMinimum(0.0);
  //h->Draw("col");

  TExec* ex = 0;
  if(view==6)
    ex = new TExec("ex_vvz","gStyle->SetPalette(8,0);");
  else
    ex = new TExec("ex_uvz","gStyle->SetPalette(8,0);");
  
  h->GetListOfFunctions()->Add(ex);
  
  delete pu;
  delete fu;
  delete pv;
  delete fv;
  
  h->SetDirectory(0);

  return h;
}

---