cmtool::CFAlgoVolumeOverlap Class Reference

#include <CFAlgoVolumeOverlap.h>

Inheritance diagram for cmtool::CFAlgoVolumeOverlap:

Public Member Functions
	CFAlgoVolumeOverlap ()
	Default constructor. More...
virtual	~CFAlgoVolumeOverlap ()
	Default destructor. More...
virtual float	Float (const std::vector< const cluster::ClusterParamsAlg * > &clusters)
virtual void	Report ()
virtual void	Reset ()
	Function to reset the algorithm instance, called together with manager's Reset() More...
void	SetVerbose (bool on)
	Function to set verbose output. More...
void	SetDebug (bool on)
	Function to set debug output. More...
void	SetUseAllPlanes (bool on)
	Function to set if _UseAllPlanes is on/off. More...
Public Member Functions inherited from cmtool::CFloatAlgoBase
	CFloatAlgoBase ()
	Default constructor. More...
virtual	~CFloatAlgoBase ()
	Default destructor. More...
Public Member Functions inherited from cmtool::CMAlgoBase
	CMAlgoBase ()
	Default constructor. More...
virtual	~CMAlgoBase ()
	Default destructor. More...
virtual void	EventBegin (const std::vector< cluster::ClusterParamsAlg > &clusters)
virtual void	EventEnd ()
virtual void	IterationBegin (const std::vector< cluster::ClusterParamsAlg > &clusters)
virtual void	IterationEnd ()
void	SetAnaFile (TFile *fout)
	Setter function for an output plot TFile pointer. More...

Private Attributes
double	_w2cm
double	_t2cm
bool	_verbose
bool	_debug
bool	_UseAllPlanes

Additional Inherited Members
Protected Attributes inherited from cmtool::CMAlgoBase
TFile *	_fout
	TFile pointer to an output file. More...
bool	_verbose
	Boolean to choose verbose mode. Turned on if CMergeManager/CMatchManager's verbosity level is >= kPerMerging. More...

Detailed Description

User implementation for CFloatAlgoBase class doxygen documentation!

Definition at line 25 of file CFAlgoVolumeOverlap.h.

Constructor & Destructor Documentation

cmtool::CFAlgoVolumeOverlap::CFAlgoVolumeOverlap

(

)

Default constructor.

Definition at line 15 of file CFAlgoVolumeOverlap.cxx.

                                           : CFloatAlgoBase()
  //-------------------------------------------------------
  {
    ::util::GeometryUtilities geou;
    _w2cm = geou.WireToCm();
    _t2cm = geou.TimeToCm();
    SetVerbose(false);
    SetDebug(false);
    SetUseAllPlanes(false); // Any plane combination OK
  }

virtual cmtool::CFAlgoVolumeOverlap::~CFAlgoVolumeOverlap ( )

inlinevirtual

Default destructor.

Definition at line 33 of file CFAlgoVolumeOverlap.h.

{};

Member Function Documentation

float cmtool::CFAlgoVolumeOverlap::Float ( const std::vector< const cluster::ClusterParamsAlg * > &  clusters )

virtual

Core function: given a set of CPANs, return a float which indicates the compatibility the cluster combination.

Reimplemented from cmtool::CFloatAlgoBase.

Definition at line 33 of file CFAlgoVolumeOverlap.cxx.

  {

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

    // Code-block by Kazu starts
    // This ensures the algorithm works only if # clusters is > 2 (and not =2)
    // You may take out this block if you want to allow matching using clusters from only 2 planes.
    if (_UseAllPlanes) { if(clusters.size()==2) return -1; }
    // Code-block by Kazu ends

    //This algorithm now works for 3 planes: find 3Dstart point from first 2 planes and find
    //How well that agrees with 3rd plane's start point location.

    //So first, make sure clusters vector has more than 1 element.
    //If not, return -1. Algorithm does not make sense
    if ( clusters.size() == 1 )
      return -1;

    if (_verbose) { std::cout << "Number of clusters taken into account: " << clusters.size() << std::endl; }

    //Nomenclature:
    //Planes: U == 0; V == 1; Y == 2.

    //Get hits for all 3 clusters
    std::vector<std::vector<util::PxHit> > Hits(3, std::vector<util::PxHit>());

    for (size_t c=0; c < clusters.size(); c++)
      Hits.at( clusters.at(c)->Plane() ) = clusters.at(c)->GetHitVector();

    //std::vector<util::PxHit> hits0 = clusters.at(0)->GetHitVector();
    //std::vector<util::PxHit> hits1 = clusters.at(1)->GetHitVector();
    //std::vector<util::PxHit> hits2 = clusters.at(2)->GetHitVector();

    //Wire Range Vector. Entries 0,1,2 correspond to planes
    std::vector<int> StartWires;
    std::vector<int> EndWires;
    std::vector<float> StartTimes;
    std::vector<float> EndTimes;
    //loop over number of planes and pre-fill StartWires and EndWires
    for (int pl=0; pl < 3; pl++){
      StartWires.push_back(9999);
      EndWires.push_back(0);
      StartTimes.push_back(9999);
      EndTimes.push_back(0);
    }

    for (size_t h=0; h < Hits.size(); h++){
      for ( auto& hit: Hits.at(h) ){
        if ( Hits.at(h).size() == 0 ){
          std::cout << "Need to insert fake hit ranges...";
        }
        else{
          if ( int(hit.w / _w2cm) < StartWires.at(h) ) { StartWires.at(h) = int(hit.w / _w2cm); }
          if ( int(hit.w / _w2cm) > EndWires.at(h) )   { EndWires.at(h)   = int(hit.w / _w2cm); }
          if ( hit.t < StartTimes.at(h) ) { StartTimes.at(h) = hit.t; }
          if ( hit.t > EndTimes.at(h) )   { EndTimes.at(h)   = hit.t; }
        }
      }//for all hits in range
    }//for all hit-lists (i.e. for all clusters)

    //if one of the plane's wire-range is still [9999, 0] then replace it
    //with min/max wire number from that plane
    //similarly for time-range
    //This allows for easy 2-Plane matching: if we are ignoring one plane
    //completely by giving the wire-range for the whole plane the intersection
    //will effectively be the intersection of the other two planes
    for (size_t h=0; h < Hits.size(); h++){
      if ( StartWires.at(h) == 9999 ) { StartWires.at(h) = 1; }
      if ( EndWires.at(h) == 0 ) { EndWires.at(h) = geo->Nwires(h)-1; }
      if ( StartTimes.at(h) == 9999 ) { StartTimes.at(h) = 0; }
      if ( EndTimes.at(h) == 0 ) { EndTimes.at(h) = 9999; }
    }


    //Find overlap in time between three planes
    if ( (StartTimes.at(1) > EndTimes.at(0)) or (StartTimes.at(1) > EndTimes.at(2)) or
         (StartTimes.at(2) > EndTimes.at(0)) or (StartTimes.at(2) > EndTimes.at(1)) or
         (StartTimes.at(0) > EndTimes.at(1)) or (StartTimes.at(0) > EndTimes.at(2)) ){
      if (_verbose) { std::cout << "No Time-Overlap!" << std::endl << std::endl; }
      return -1;
    }
    //do two planes at a time...start with U and V
    float Tmin=9999, Tmax=0;
    ( StartTimes.at(1) > StartTimes.at(0) ) ? (Tmin = StartTimes.at(1)) : (Tmin = StartTimes.at(0));
    if (StartTimes.at(2) > Tmin) { Tmin = StartTimes.at(2); }
    ( EndTimes.at(1) < EndTimes.at(0) ) ? (Tmax = EndTimes.at(1)) : (Tmax = EndTimes.at(0));
    if (EndTimes.at(2) < Tmax) { Tmax = EndTimes.at(2); }
    float Trange = Tmax-Tmin;

    //Now get start & end points of all these wires
    Double_t xyzStart0WireStart[3] = {0., 0., 0.}; //xyz array info of start point for smallest wire number on plane 0
    Double_t xyzStart0WireEnd[3]   = {0., 0., 0.};
    Double_t xyzEnd0WireStart[3] = {0., 0., 0.};
    Double_t xyzEnd0WireEnd[3]   = {0., 0., 0.};
    Double_t xyzStart1WireStart[3] = {0., 0., 0.}; //xyz array info of start point for smallest wire number on plane 1
    Double_t xyzStart1WireEnd[3]   = {0., 0., 0.};
    Double_t xyzEnd1WireStart[3] = {0., 0., 0.};
    Double_t xyzEnd1WireEnd[3]   = {0., 0., 0.};
    Double_t xyzStart2WireStart[3] = {0., 0., 0.}; //xyz array info of start point for smallest wire number on plane 2
    Double_t xyzStart2WireEnd[3]   = {0., 0., 0.};
    Double_t xyzEnd2WireStart[3] = {0., 0., 0.};
    Double_t xyzEnd2WireEnd[3]   = {0., 0., 0.};

    if (_debug) {
      std::cout << "Wire Ranges:" << std::endl;
      std::cout << "U-Plane: [ " << StartWires.at(0) << ", " << EndWires.at(0) << "]" << std::endl;
      std::cout << "V-Plane: [ " << StartWires.at(1) << ", " << EndWires.at(1) << "]" << std::endl;
      std::cout << "Y-Plane: [ " << StartWires.at(2) << ", " << EndWires.at(2) << "]" << std::endl;
      std::cout << std::endl;
    }

    /*
    geo->WireEndPoints(0, StartWires.at(0), xyzStart0WireStart, xyzStart0WireEnd);
    geo->WireEndPoints(0, EndWires.at(0), xyzEnd0WireStart, xyzEnd0WireEnd);
    geo->WireEndPoints(1, StartWires.at(1), xyzStart1WireStart, xyzStart1WireEnd);
    geo->WireEndPoints(1, EndWires.at(1), xyzEnd1WireStart, xyzEnd1WireEnd);
    geo->WireEndPoints(2, StartWires.at(2), xyzStart2WireStart, xyzStart2WireEnd);
    geo->WireEndPoints(2, EndWires.at(2), xyzEnd2WireStart, xyzEnd2WireEnd);
    */
    geo->WireEndPoints(0, 0, 0, StartWires.at(0), xyzStart0WireStart, xyzStart0WireEnd);
    geo->WireEndPoints(0, 0, 0, EndWires.at(0), xyzEnd0WireStart, xyzEnd0WireEnd);
    geo->WireEndPoints(0, 0, 1, StartWires.at(1), xyzStart1WireStart, xyzStart1WireEnd);
    geo->WireEndPoints(0, 0, 1, EndWires.at(1), xyzEnd1WireStart, xyzEnd1WireEnd);
    geo->WireEndPoints(0, 0, 2, StartWires.at(2), xyzStart2WireStart, xyzStart2WireEnd);
    geo->WireEndPoints(0, 0, 2, EndWires.at(2), xyzEnd2WireStart, xyzEnd2WireEnd);

    //check that z-positions for plane2 are the same...if not error!
    //if ( xyzStart2WireStart[2] != xyzStart2WireEnd[2] ) { std::cout << "Y-wire does not have same Z start and end..." << std::endl; }
    double zMin = xyzStart2WireStart[2];
    //if ( xyzEnd2WireStart[2] != xyzEnd2WireEnd[2] ) { std::cout << "Y-wire does not have same Z start and end..." << std::endl; }
    double zMax = xyzEnd2WireStart[2];

    //Plane U == Plane 0
    //Plane V == Plane 1
    //Plane Y == Plane 2

    //Each line can be described by function y = s*z + b (x,y,z coordinates same as uBooNE coord.)
    //Slopes known: Pl0 is +60 clockwise from vertical. Pl1 is -60, counterclockwise from vertical. Looking at TPC with beam from left
    double slopeU = tan(30*3.14/180.);
    double slopeV = tan(-30*3.14/180.);

    //Find intercepts:
    double bUup = xyzStart0WireStart[1] - xyzStart0WireStart[2] * slopeU;
    double bUdn   = xyzEnd0WireStart[1] - xyzEnd0WireStart[2] * slopeU;
    double bVdn = xyzStart1WireStart[1] - xyzStart1WireStart[2] * slopeV;
    double bVup   = xyzEnd1WireStart[1] - xyzEnd1WireStart[2] * slopeV;
    //make sure we know which line is above which
    if ( bUup < bUdn ) { std::cout << "Uup and Udn are mixed up!" << std::endl; }
    if ( bVdn > bVup ) { std::cout << "Vup and Vdn are mixed up!" << std::endl; }
    //Pl2 lines are vertical...slope is infinite and no intercept.

    //Find intercepts of U wire-ranges with Y plane (easy since vertical)
    //For now assume wire-ranges go to infinity, worry about TPC constraints later

    //Plug in Y-plane zMin and zMax coordinates into equations for U/V wire lines
    double VdnZmin = slopeV * zMin + bVdn;
    double VdnZmax = slopeV * zMax + bVdn;
    double VupZmin = slopeV * zMin + bVup;
    double VupZmax = slopeV * zMax + bVup;
    double UdnZmin = slopeU * zMin + bUdn;
    double UdnZmax = slopeU * zMax + bUdn;
    double UupZmin = slopeU * zMin + bUup;
    double UupZmax = slopeU * zMax + bUup;

    if (_debug){
      std::cout << "Y-Plane and U-Plane points [Z,Y]:" << std::endl;
      std::cout << "\t\t[ " << zMin << ", " << UdnZmin << "]" << std::endl;
      std::cout << "\t\t[ " << zMin << ", " << UupZmin << "]" << std::endl;
      std::cout << "\t\t[ " << zMax << ", " << UupZmax << "]" << std::endl;
      std::cout << "\t\t[ " << zMax << ", " << UdnZmax << "]" << std::endl;
      std::cout << "Y-Plane and V-Plane points [Z,Y]:" << std::endl;
      std::cout << "\t\t[ " << zMin << ", " << VdnZmin << "]" << std::endl;
      std::cout << "\t\t[ " << zMin << ", " << VupZmin << "]" << std::endl;
      std::cout << "\t\t[ " << zMax << ", " << VupZmax << "]" << std::endl;
      std::cout << "\t\t[ " << zMax << ", " << VdnZmax << "]" << std::endl;
    }
    //We now have Two polygons:
    //One is the intersection of Y-Plane wires with U-Plane wires
    //The other the intersection of planes Y and V.
    //The intersection points of these two polygons is the
    //overall intersection Area of the 3 clusters on the Y-Z plane.

    //Go through all segment intersections. If one is found add to
    //list of points.
    //Create a list of points for polygon
    std::vector< std::pair<float,float> > WireIntersection;
    double zInt; // temporary holder for z-intersection point of oblique sides
    //Check: Vup and Uup, Vup and Uright, Vup and Uleft, Vup and Udn.
    //Intersection between Vup and Uup: if within zMin, zMax then ok!
    zInt = (bUup-bVup)/(slopeV-slopeU);
    if ( (zInt > zMin) and (zInt < zMax) )
      WireIntersection.push_back( std::make_pair( zInt, slopeV*zInt+bVup ) );
    //Intersection between Vup and Uright:
    if ( (VupZmax < UupZmax) and (VupZmax > UdnZmax) )
      WireIntersection.push_back( std::make_pair( zMax, VupZmax ) );
    //Intersection between Vup and Uleft:
    if ( (VupZmin < UupZmin) and ( VupZmin > UdnZmin) )
      WireIntersection.push_back( std::make_pair( zMin, VupZmin ) );
    //Intersection between Vup and Udn:
    zInt = (bUdn-bVup)/(slopeV-slopeU);
    if ( (zInt > zMin) and (zInt < zMax) )
      WireIntersection.push_back( std::make_pair( zInt, slopeV*zInt+bVup ) );

    //Check: Vdn and Uup, Uright, Uleft, Udn:
    //Intersection between Vdn and Uup:
    zInt = (bUup-bVdn)/(slopeV-slopeU);
    if ( (zInt > zMin) and (zInt < zMax) )
      WireIntersection.push_back( std::make_pair( zInt, slopeV*zInt+bVdn ) );
    //Intersection between Vdn and Uright:
    if ( (VdnZmax < UupZmax) and (VdnZmax > UdnZmax) )
      WireIntersection.push_back( std::make_pair( zMax, VdnZmax ) );
    //Intersection between Vdn and Uleft:
    if ( (VdnZmin < UupZmin) and ( VdnZmin > UdnZmin) )
      WireIntersection.push_back( std::make_pair( zMin, VdnZmin ) );
    //Intersection between Vdn and Udn:
    zInt = (bUdn-bVdn)/(slopeV-slopeU);
    if ( (zInt > zMin) and (zInt < zMax) )
      WireIntersection.push_back( std::make_pair( zInt, slopeV*zInt+bVdn ) );

    //Check: Vright and Uup, Udn:
    //Intersection between Vright and Uup:
    if ( (UupZmax < VupZmax) and ( UupZmax > VdnZmax) )
      WireIntersection.push_back( std::make_pair( zMax, UupZmax ) );
    //Intersection between Vright and Udn:
    if ( (UdnZmax < VupZmax) and ( UdnZmax > VdnZmax) )
      WireIntersection.push_back( std::make_pair( zMax, UdnZmax ) );

    //Check Vleft and Uup, Udn:
    //Intersection between Vleft and Uup:
    if ( (UupZmin < VupZmin) and ( UupZmin > VdnZmin) )
      WireIntersection.push_back( std::make_pair( zMin, UupZmin ) );
    //Intersection between Vleft and Udn:
    if ( (UdnZmin < VupZmin) and ( UdnZmin > VdnZmin) )
      WireIntersection.push_back( std::make_pair( zMin, UdnZmin ) );

    //If length is 0 then no intersection...return -1
    if ( WireIntersection.size() == 0 ){
      if (_verbose) { std::cout << "No intersection..." << std::endl << std::endl; }
      return -1;
    }

    //Now our polygon is complete...
    //need to disentangle in case points added in incorrect order
    //then calculate area
    Polygon2D WirePolygon(WireIntersection);
    //Variable to hold final output Area
    double PolyArea = -1;
    //Check order
    WirePolygon.UntanglePolygon();

    //Create a Polygon for the Y-Z TPC Plane
    std::vector< std::pair<float,float> > TPCCorners;
    TPCCorners.push_back( std::make_pair(0., -geo->DetHalfHeight()) );
    TPCCorners.push_back( std::make_pair(geo->DetLength(), -geo->DetHalfHeight()) );
    TPCCorners.push_back( std::make_pair(geo->DetLength(), geo->DetHalfHeight()) );
    TPCCorners.push_back( std::make_pair(0., geo->DetHalfHeight()) );
    Polygon2D TPCPolygon(TPCCorners);
    if (TPCPolygon.Contained(WirePolygon) ){
      if (_verbose) {
        std::cout << "Wire Overlap contained in TPC" << std::endl;
        std::cout << "Intersection Polygon Coordinates [Z,Y]: " << std::endl;
        for (unsigned int s=0; s < WirePolygon.Size(); s++)
          std::cout << "\t\t[ " << WirePolygon.Point(s).first << ", " << WirePolygon.Point(s).second << "]" << std::endl;
        std::cout << std::endl;
      }
      PolyArea = WirePolygon.Area();
    }
    else {
      if (_verbose) {
        std::cout << "Wire overlap not fully contained in TPC" << std::endl;
        std::cout << "Intersection Polygon Coordinates [Z,Y]: " << std::endl;
        for (unsigned int s=0; s < WirePolygon.Size(); s++)
          std::cout << "\t\t[ " << WirePolygon.Point(s).first << ", " << WirePolygon.Point(s).second << "]" << std::endl;
        std::cout << std::endl;
      }
      //product polygon should be the intersection of WirePolygon and TPCPolygon
      Polygon2D IntersectionPolygon(TPCPolygon, WirePolygon);
      PolyArea = IntersectionPolygon.Area();
    }

    //return polygon area -> larger = better!
    if (_verbose) { std::cout << "Intersection area: " << PolyArea << std::endl << std::endl; }
    return  Trange*PolyArea;
  }

void cmtool::CFAlgoVolumeOverlap::Report ( )

virtual

Optional function: called after each iterative approach if a manager class is run with verbosity level <= kPerIteration. Maybe useful for debugging.

Reimplemented from cmtool::CMAlgoBase.

Definition at line 321 of file CFAlgoVolumeOverlap.cxx.

  {

  }

void cmtool::CFAlgoVolumeOverlap::Reset ( void  )

virtual

Function to reset the algorithm instance, called together with manager's Reset()

Reimplemented from cmtool::CMAlgoBase.

Definition at line 27 of file CFAlgoVolumeOverlap.cxx.

  {
  }

void cmtool::CFAlgoVolumeOverlap::SetDebug ( bool  on )

inline

Function to set debug output.

Definition at line 60 of file CFAlgoVolumeOverlap.h.

{ _debug = on; }

void cmtool::CFAlgoVolumeOverlap::SetUseAllPlanes ( bool  on )

inline

Function to set if _UseAllPlanes is on/off.

Definition at line 63 of file CFAlgoVolumeOverlap.h.

{ _UseAllPlanes = on; }

void cmtool::CFAlgoVolumeOverlap::SetVerbose ( bool  on )

inlinevirtual

Function to set verbose output.

Reimplemented from cmtool::CMAlgoBase.

Definition at line 57 of file CFAlgoVolumeOverlap.h.

{ _verbose = on; }

Member Data Documentation

bool cmtool::CFAlgoVolumeOverlap::_debug

private

Definition at line 69 of file CFAlgoVolumeOverlap.h.

double cmtool::CFAlgoVolumeOverlap::_t2cm

private

Definition at line 67 of file CFAlgoVolumeOverlap.h.

bool cmtool::CFAlgoVolumeOverlap::_UseAllPlanes

private

Definition at line 70 of file CFAlgoVolumeOverlap.h.

bool cmtool::CFAlgoVolumeOverlap::_verbose

private

Definition at line 68 of file CFAlgoVolumeOverlap.h.

double cmtool::CFAlgoVolumeOverlap::_w2cm

private

Definition at line 67 of file CFAlgoVolumeOverlap.h.

---

The documentation for this class was generated from the following files:

• larreco/larreco/RecoAlg/CMTool/CMTAlgMatch/CFAlgoVolumeOverlap.h
• larreco/larreco/RecoAlg/CMTool/CMTAlgMatch/CFAlgoVolumeOverlap.cxx