ViewMatchAlg.cxx

Go to the documentation of this file.

///
/// \file    ViewMatchAlg.cxx
/// \brief   3d matching of prongs
/// \author  eniner@indiana.edu
/// \version $Id: ViewMatchAlg.cxx,v 1.6 2012-11-05 22:36:04 edniner Exp $
///
#include "TensorFlowEvaluator/NERDProng/ViewMatchAlg.h"

#include "Calibrator/Calibrator.h"
#include "GeometryObjects/PlaneGeo.h"
#include "GeometryObjects/CellGeo.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Prong.h"
#include "RecoBase/RecoHit.h"
#include "RecoBase/PID.h"
#include "Utilities/func/MathUtil.h"

#include "TVector3.h"

#include "art/Framework/Services/Registry/ServiceHandle.h"

#include <iostream>
#include <cmath>

namespace nerd{

  ViewMatchAlg::ViewMatchAlg(const ViewMatchParams& /*params*/)
  {
  }

  //......................................................................
  
  ViewMatchAlg::~ViewMatchAlg()
  {
  }

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

  /// Load in candidate 2D prongs to prepare for matching
  void ViewMatchAlg::LoadProngs(std::vector<rb::Prong> prongs, std::vector<std::vector<rb::PID>> scores)
  {
    fXProngs.clear();
    fXScores.clear();
    fYProngs.clear();
    fYScores.clear();
    fMatchedProngs.clear();
    fUnmatchedProngs.clear();
    fMatchedXScores.clear();
    fMatchedYScores.clear();
    fUnmatchedScores.clear();

    for (size_t i=0; i< prongs.size(); i++){
      if (prongs[i].View() == geo::kX){
        fXProngs.push_back(prongs[i]);
        fXScores.push_back(scores[i]);
      }
      if (prongs[i].View() == geo::kY){
        fYProngs.push_back(prongs[i]);
        fYScores.push_back(scores[i]);
      }
      if (prongs[i].View() == geo::kXorY)
        std::cout<<"A prong is already 3D ?!?!?!?!?"<<std::endl;
    }
  }
  
  //......................................................................
  
  /// Function to do the work and execute the matching
  void ViewMatchAlg::Matching()
  {
    art::ServiceHandle<geo::Geometry> fGeo;

    fMatchedProngs.clear();
    fMatchedXScores.clear();
    fMatchedYScores.clear();
    fUnmatchedProngs.clear();
    fUnmatchedScores.clear();

    //clear containers to hold summary debuggin information for matching ntuple
    xzid.clear();
    yzid.clear();
    xzs.clear();
    yzs.clear();
    xzpe.clear();
    yzpe.clear();
    ksscore.clear();
    ksshift.clear();
    
    //size matrix to hold matching score for all pairs of tracks
    double score[fXProngs.size()][fYProngs.size()];
     
    //loop over xz prongs to create matching pairs
    for (size_t i=0; i<fXProngs.size(); i++){
      int ixPlane = fXProngs[i].MinPlane();
      int fxPlane = fXProngs[i].MaxPlane();
      double dirX = fXProngs[i].Dir().Z();

      //if we have a  1 plane track, alter the direction for maximum tilt, defense against near vertical tracks
      if ((ixPlane == fxPlane)||(acos(fabs(fXProngs[i].Dir().X()))<0.01)) {
        if (dirX > 0.0){
          if (fXProngs[i].Dir().X() > 0.0){
            fXProngs[i].SetDir(fGeo->Plane(0)->Cell(0)->HalfW()+fXProngs[i].MaxX()-fXProngs[i].Start().X(),
                              fGeo->Plane(0)->Cell(0)->HalfW()+fXProngs[i].MaxZ()-fXProngs[i].Start().Z());
          }
          else if (fXProngs[i].Dir().X() < 0.0){
            fXProngs[i].SetDir(-fGeo->Plane(0)->Cell(0)->HalfW()+fXProngs[i].MinX()-fXProngs[i].Start().X(),
                              fGeo->Plane(0)->Cell(0)->HalfW()+fXProngs[i].MaxZ()-fXProngs[i].Start().Z());
          }
        }
        else if (dirX < 0.0){
          if (fXProngs[i].Dir().X() > 0.0){
            fXProngs[i].SetDir(fGeo->Plane(0)->Cell(0)->HalfW()+fXProngs[i].MaxX()-fXProngs[i].Start().X(),
                              -fGeo->Plane(0)->Cell(0)->HalfW()+fXProngs[i].MinZ()-fXProngs[i].Start().Z());
          }
          else if (fXProngs[i].Dir().X() < 0.0){
            fXProngs[i].SetDir(-fGeo->Plane(0)->Cell(0)->HalfW()+fXProngs[i].MinX()-fXProngs[i].Start().X(),
                              -fGeo->Plane(0)->Cell(0)->HalfW()+fXProngs[i].MinZ()-fXProngs[i].Start().Z());
          }
        }
      }

      //inner loop over yz prongs to make pairs
      for (size_t j=0; j<fYProngs.size(); j++){
        int iyPlane = fYProngs[j].MinPlane();
        int fyPlane = fYProngs[j].MaxPlane();
        double dirY = fYProngs[j].Dir().Z();

        //if we have a  1 plane track, alter the direction for maximum tilt, defense against near vertical tracks
        if (((iyPlane == fyPlane)||(acos(fabs(fYProngs[j].Dir().Y()))<0.01))&&(i==0)) {
          if (dirY > 0.0){
            if (fYProngs[j].Dir().Y() > 0.0){
              fYProngs[j].SetDir(fGeo->Plane(0)->Cell(0)->HalfW()+fYProngs[j].MaxY()-fYProngs[j].Start().Y(),
                                fGeo->Plane(0)->Cell(0)->HalfW()+fYProngs[j].MaxZ()-fYProngs[j].Start().Z());
            }
            else if (fYProngs[j].Dir().Y() < 0.0){
              fYProngs[j].SetDir(-fGeo->Plane(0)->Cell(0)->HalfW()+fYProngs[j].MinY()-fYProngs[j].Start().Y(),
                                fGeo->Plane(0)->Cell(0)->HalfW()+fYProngs[j].MaxZ()-fYProngs[j].Start().Z());
            }
          }
          else if (dirY < 0.0){
            if (fYProngs[j].Dir().Y() > 0.0){
              fYProngs[j].SetDir(fGeo->Plane(0)->Cell(0)->HalfW()+fYProngs[j].MaxY()-fYProngs[j].Start().Y(),
                                -fGeo->Plane(0)->Cell(0)->HalfW()+fYProngs[j].MinZ()-fYProngs[j].Start().Z());
            }
            else if (fYProngs[j].Dir().Y() < 0.0){
              fYProngs[j].SetDir(-fGeo->Plane(0)->Cell(0)->HalfW()+fYProngs[j].MinY()-fYProngs[j].Start().Y(),
                                -fGeo->Plane(0)->Cell(0)->HalfW()+fYProngs[j].MinZ()-fYProngs[j].Start().Z());
            }
          }
        }

        score[i][j] = -1.0;
        
        //make sure candidate prong in both views are either forwards are backwards
        if (dirX*dirY < 0.0 ) continue;

        //create 3d test prong to compute score
        TVector3 newDir;
        TVector3 vert;
        rb::Prong testprong = this->MakeTestProng(fXProngs[i], fYProngs[j], &newDir, &vert);

        //clamp vertex to detector volume to find which view we are in for shift
        double vx = vert.X();
        double vy = vert.Y();
        double vz = vert.Z();
        if (vert.X() < -fGeo->DetHalfWidth()) vx = -fGeo->DetHalfWidth()+15.0;
        if (vert.X() > fGeo->DetHalfWidth()) vx = fGeo->DetHalfWidth()-15.0;
        if (vert.Y() < -fGeo->DetHalfHeight()) vy = -fGeo->DetHalfHeight()+15.0;
        if (vert.Y() > fGeo->DetHalfHeight()) vy = fGeo->DetHalfHeight()-15.0;
        if (vert.Z() < 0.0) vz = 2.0;
        if (vert.Z() > fGeo->DetLength()) vz = fGeo->DetLength()-2.0;
        //get plane position of vertex, roam in cell if vertex was in the plastic
        geo::CellUniqueId cid = fGeo->CellId(vx, vy, vz, 1.5, 1.5, 6., 1.0);
        geo::View_t vertview = geo::kX;
        unsigned int vertplane = 999;
        if (cid != 0){
          vertplane = fGeo->getPlaneID(cid);
          vertview = (fGeo->Plane(vertplane))->View();
        };

        //calculate energy profiles in each view
        double xgev = 0.0;
        double ygev = 0.0;
        std::vector<std::pair<double, double> > xlist = this->CalcEnergyProfile(testprong, newDir, vert, &xgev, geo::kX);
        std::vector<std::pair<double, double> > ylist = this->CalcEnergyProfile(testprong, newDir, vert, &ygev, geo::kY);

        //if the vertex is in the x view, shift all points on that track the equivalent of one cell forward or backward,
        //based on the slope of the track.  
        double shift = 0.0;
        double xyz[3];
        double shiftz = vz;
        //shift only if the vertex plane was defined
        if ((vertview == geo::kX)&&(vertplane != 999)){
          //if track in positive direction shift one plane up
          if (dirX > 0.0){
            unsigned int nplane = fGeo->NextPlaneOtherView(vertplane,1);
            //make sure next plane is defined
            if (nplane != 999000){ 
              fGeo->Plane(nplane)->Cell(0)->GetCenter(xyz);
              shiftz = xyz[2];
            }
            double shiftx = (testprong.Dir().X()/testprong.Dir().Z())*(shiftz-vz) + vx;
            double shifty = (testprong.Dir().Y()/testprong.Dir().Z())*(shiftz-vz) + vy;
            shift = util::pythag(vx-shiftx,vy-shifty,vz-shiftz);
          }
          //shift one plane down for backwards track
          else{
            unsigned int nplane = fGeo->NextPlaneOtherView(vertplane,-1);
            if (nplane != 999000){
              fGeo->Plane(nplane)->Cell(0)->GetCenter(xyz);
              shiftz = xyz[2];
            }
            double shiftx = (testprong.Dir().X()/testprong.Dir().Z())*(shiftz-vz) + vx;
            double shifty = (testprong.Dir().Y()/testprong.Dir().Z())*(shiftz-vz) + vy;
            shift = -util::pythag(vx-shiftx,vy-shifty,vz-shiftz);
          }
        }
        //else if the plane was in the YZ view, shift xz view accordingly
        if (vertview == geo::kY){
           if (dirX > 0.0){
            unsigned int nplane = fGeo->NextPlaneOtherView(vertplane,1);
            if (nplane != 999000){
              fGeo->Plane(nplane)->Cell(0)->GetCenter(xyz);
              shiftz = xyz[2];
            }
            double shiftx = (testprong.Dir().X()/testprong.Dir().Z())*(shiftz-vz) + vx;
            double shifty = (testprong.Dir().Y()/testprong.Dir().Z())*(shiftz-vz) + vy;
            shift = -util::pythag(vx-shiftx,vy-shifty,vz-shiftz);
          }
          else{
            unsigned int nplane = fGeo->NextPlaneOtherView(vertplane,-1);
            if (nplane != 999000){
              fGeo->Plane(nplane)->Cell(0)->GetCenter(xyz);
              shiftz = xyz[2];
            }
            double shiftx = (testprong.Dir().X()/testprong.Dir().Z())*(shiftz-vz) + vx;
            double shifty = (testprong.Dir().Y()/testprong.Dir().Z())*(shiftz-vz) + vy;
            shift = util::pythag(vx-shiftx,vy-shifty,vz-shiftz);
          }
        }

        //if tracks from each view overlap by at least one plane, attempt a match
        if(iyPlane <= fxPlane && fyPlane >= ixPlane){
          
          //perform kuiper test on candidate
          double kScore = KuiperTest(xlist, ylist, shift, xgev, ygev);

          /////////////////////////
          //block to collect summary information to monitor view matching process
          //xzid.push_back(i);
          //yzid.push_back(j);
          //ksscore.push_back(kScore);
          //std::vector<double> xstemp;
          //std::vector<double> ystemp;
          //std::vector<double> xpetemp;
          //std::vector<double> ypetemp;
          //for (unsigned int ix=0; ix<xlist.size(); ++ix){
          //  xstemp.push_back(xlist[ix].first);
          //  xpetemp.push_back(xlist[ix].second);
          //}
          //for (unsigned int iy=0; iy<ylist.size(); ++iy){
          //  ystemp.push_back(ylist[iy].first);
          //  ypetemp.push_back(ylist[iy].second);
          //}
          //xzs.push_back(xstemp);
          //yzs.push_back(ystemp);
          //xzpe.push_back(xpetemp);
          //yzpe.push_back(ypetemp);
          //////////////////////////////

          //save score as long as the tracks overlap, or if the track is only 1 plane in one view, 
          //then in the other view that track can only occupy the plane on either side
          if((fYProngs[j].ExtentPlane()>1) && (fXProngs[i].ExtentPlane()>1)) score[i][j] = kScore;
          else if((fYProngs[j].ExtentPlane() == 1) && (ixPlane == (iyPlane-1)) && (fxPlane == (iyPlane+1))) score[i][j] = kScore;
          else if((fXProngs[i].ExtentPlane() == 1) && (iyPlane == (ixPlane-1)) && (fyPlane == (ixPlane+1))) score[i][j] = kScore;      
        }
        //case if both tracks are in a single plane
        else if((ixPlane == fxPlane) && (iyPlane == fyPlane) && (dirX*dirY > 0.0)){
          // if x and y prongs appear in only one plane, allow match only if it is adjacent plane
          if(ixPlane == (iyPlane-1) || ixPlane == (iyPlane+1)){
            //perform test
            double kScore = KuiperTest(xlist, ylist, shift, xgev, ygev);

            //////////////////////////////////////////
            //block to collect summary information to monitor view matching process
            //xzid.push_back(i);
            //yzid.push_back(j);
            //ksscore.push_back(kScore);
            //std::vector<double> xstemp;
            //std::vector<double> ystemp;
            //std::vector<double> xpetemp;
            //std::vector<double> ypetemp;
            //for (unsigned int ix=0; ix<xlist.size(); ++ix){
            //  xstemp.push_back(xlist[ix].first);
            //  xpetemp.push_back(xlist[ix].second);
            //}
            //for (unsigned int iy=0; iy<ylist.size(); ++iy){
            //  ystemp.push_back(ylist[iy].first);
            //  ypetemp.push_back(ylist[iy].second);
            //}
            //xzs.push_back(xstemp);
            //yzs.push_back(ystemp);
            //xzpe.push_back(xpetemp);
            //yzpe.push_back(ypetemp);
            ///////////////////////////////////////

            score[i][j] = kScore;      
          }

        }//end case for single plane tracks
      } // end yprong loop
    } //end xprong loop
      
    // find the best score of all the possible matches.
    // Now get the best matched prongs out of the group based on lowest score (not zero)
    bool doneMatching = false;
    std::vector<int> xMatched(fXProngs.size());
    for(unsigned int i = 0; i<xMatched.size();++i){
      xMatched[i] = 0;
    }
    std::vector<int> yMatched(fYProngs.size());
    for(unsigned int i = 0; i<yMatched.size();++i){
      yMatched[i] = 0;
    }
    while(!doneMatching){
      int xbest = -1;
      int ybest = -1;
      double bestscore = -1;
      // loop over all possible matched prongs to find the best one
      for(unsigned int ix = 0; ix<fXProngs.size();++ix){
        for(unsigned int iy = 0; iy<fYProngs.size();++iy){
          if((score[ix][iy] <= bestscore || bestscore == -1) && (score[ix][iy] >= 0.0)){
            xbest = ix;
            ybest = iy;
            bestscore = score[ix][iy];
          }
        }
      }

      // store the best matched prong and reset things for the next best prong
      if(bestscore >= 0){
        // make a 3d prong out of the x and y prongs
        rb::Prong newprong;

        // add all of the hits to the prong.
        for(unsigned int i = 0; i<fXProngs[xbest].NCell();++i){
          newprong.Add(fXProngs[xbest].Cell(i), fXProngs[xbest].Weight(i));
        }
        for(unsigned int i = 0; i<fYProngs[ybest].NCell();++i){
          newprong.Add(fYProngs[ybest].Cell(i), fYProngs[ybest].Weight(i));
        }
        fMatchedXScores.push_back(fXScores[xbest]);
        fMatchedYScores.push_back(fYScores[ybest]);

        newprong.SetID(fXProngs[xbest].ID());
        TVector3 S;
        S[0] = fXProngs[xbest].Start().X();
        S[1] = fYProngs[ybest].Start().Y();
        S[2] = fXProngs[xbest].Start().Z();
        newprong.SetStart(S);
        double tanThetax = tan(acos(fXProngs[xbest].Dir().X()));
        double tanThetay = tan(acos(fYProngs[ybest].Dir().Y()));
        double cosThetax = 1/(tanThetax*sqrt(1+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
        double cosThetay = 1/(tanThetay*sqrt(1+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
        double newz = sqrt(1-cosThetax*cosThetax-cosThetay*cosThetay);
        if ((fXProngs[xbest].Dir().Z() < 0.0) && (fYProngs[ybest].Dir().Z() < 0.0)) newz = -newz;
        TVector3 newDir(cosThetax,cosThetay,newz);      
        newprong.SetDir(newDir);

        fMatchedProngs.push_back(newprong);
        xMatched[xbest] = 1;
        yMatched[ybest] = 1;

        // set all the prongs that are matched to either this x or y prong to have a null match score.
        for(unsigned int i = 0; i<fYProngs.size();++i){
          score[xbest][i] = -1.0;
        }
        for(unsigned int i = 0; i<fXProngs.size();++i){
          score[i][ybest] = -1.0;
        }

      }
      else{ doneMatching = true; }
    }//end loop to make matches
    
    for (unsigned int i=0; i<xMatched.size(); i++){
      if (xMatched[i] == 0){
        fUnmatchedProngs.push_back(fXProngs[i]);
        fUnmatchedScores.push_back(fXScores[i]);
      }
    }
    for (unsigned int i=0; i<yMatched.size(); i++){
      if (yMatched[i] == 0){
        fUnmatchedProngs.push_back(fYProngs[i]);
        fUnmatchedScores.push_back(fYScores[i]);
      }
    }
    
    
  }

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

  /// Function to calculate the energy profile of the test prong for a given view
  std::vector<std::pair<double,double> > ViewMatchAlg::CalcEnergyProfile(rb::Prong testprong, TVector3 dir, TVector3 start, double *totgev, geo::View_t v)
  {
    art::ServiceHandle<geo::Geometry> fGeo;
    art::ServiceHandle<calib::Calibrator> cal;

    std::vector<std::pair<double,double> > list;

    unsigned int ncells = testprong.NCell(v);
    unsigned int ip, ic, n, m;
    double w; // Distance of closest approach from read-out end of cell
    double s; // Distance of closest approach from track start point
      
    //treat each cell as a 9 point space
    double dx[3] = {-0.67, 0.0, 0.67};
    double dy[3] = {-1.0,  0.0, 1.0};

    TVector3 pc; //point in cell
    TVector3 pt; //point on track

    *totgev = 0.0;

    //loop over hits to make profile
    for (unsigned int ih=0; ih<ncells; ++ih){
      const rb::CellHit* ch = testprong.Cell(v, ih).get();
      ip = ch->Plane();
      ic = ch->Cell();

      rb::RecoHit rhit(cal->MakeRecoHit(*ch, testprong.W(ch)));

      if (!rhit.IsCalibrated()) continue;
      //divide by 9.0 since we are breaking the cell into 9 pieces
      double gev = (rhit.GeV()*testprong.Weight(v, ih))/9.0;

      //break each cell into 9 points to smear out energy
      for (n=0; n<3; ++n) {
        for (m=0; m<3; ++m) {
          *totgev += gev;
          s = -1.0; //the geometry function will not calculate s if it has a value of 0.0
          fGeo->ClosestApproach(ip, ic, start, dir, dx[n], dy[m], &w, &s, &pc, &pt);
          list.push_back(std::make_pair(s, gev));
        }
      }
    }//end loop over cells

    //sort list from shortest to longest path
    std::sort(list.begin(), list.end());

    return list;
  }

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

  /// Function to combine a prong from each view into a test candidate
  rb::Prong   ViewMatchAlg::MakeTestProng(rb::Prong& xprong, rb::Prong& yprong, TVector3* newDir, TVector3* vert)
  {
    rb::Cluster clus;
    for (unsigned int k=0; k<xprong.NCell(); k++) clus.Add(xprong.Cell(k), xprong.Weight(k));
    for (unsigned int k=0; k<yprong.NCell(); k++) clus.Add(yprong.Cell(k), yprong.Weight(k));

    //make a direction for this track candidate in 3D
    double tanThetax = tan(acos(xprong.Dir().X()));
    double tanThetay = tan(acos(yprong.Dir().Y()));
    double cosThetax = 1.0/(tanThetax*sqrt(1.0+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
    double cosThetay = 1.0/(tanThetay*sqrt(1.0+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
    double newz = sqrt(1.0-cosThetax*cosThetax-cosThetay*cosThetay);
    if ((xprong.Dir().Z() < 0.0) && (yprong.Dir().Z() < 0.0)) newz = -newz;
    newDir->SetXYZ(cosThetax,cosThetay,newz);

    //construct vertex candidate.
    //At this point all input prongs were forced to start from the elastic arms vertex. 
    //So the z coordinate will be the same
    vert->SetXYZ(xprong.Start().X(),yprong.Start().Y(),xprong.Start().Z());

    rb::Prong testprong(clus, *vert, *newDir);

    return testprong;

  }

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

  /// Function to perform the kuiper test to compare the similarity of the energy profiles between views
  double ViewMatchAlg::KuiperTest(std::vector<std::pair<double,double>> xlist, 
                                 std::vector<std::pair<double,double>> ylist,  
                                 double shift, 
                                 double xgev, 
                                 double ygev)
  {

    //if the track in either view did not have any calibrated energy, do not attempt a match
    if ((xgev <= 0.0)||(ygev <= 0.0)){
      return -1.0;
    }

    //now perform ks test to determine maximum distance for these profiles

    double xsum = 0.0;//running sum of cumulative energy in xz view
    double ysum = 0.0;//running sum of cumulative energy in yz view
    double sx;//current point in xz track
    double sy;//current point in yz track
    double ksdistp;//positive difference between profiles
    double ksdistm;//negative difference between profiles
    double kuiper;//score metric
    double ksmaxp = -0.00001;//largest positive difference between profiles
    double ksmaxm = -0.00001;//largest negative difference between profiles

    unsigned int ix=0;
    unsigned int iy=0;

    double  kuiperBest = 2.0;
    double  shiftBest = -13.0;
    
    //loop to vary the shift around the estimate calculated earlier
    //calculate the kuiper statistic for each shift value, keep the result with the lowest score
    for (double shiftdist = shift-12.0; shiftdist<shift+12.01; shiftdist+=0.5){
      ksmaxp = -0.00001;
      ksmaxm = -0.00001;
      ix = 0;
      iy = 0;
      xsum = 0.0;
      ysum = 0.0;


      while ((ix < xlist.size()) && (iy < ylist.size())) {
        sx = xlist[ix].first + shiftdist;
        sy = ylist[iy].first;
        if (sx <= sy){
          xsum += xlist[ix].second;
          ix++;
        }
        else if (sy < sx){
          ysum += ylist[iy].second;
          iy++;
        }
      
        ksdistp = ysum/ygev - xsum/xgev;
        ksdistm = xsum/xgev - ysum/ygev;
    
        if (ksdistp > ksmaxp) ksmaxp = ksdistp;
        if (ksdistm > ksmaxm) ksmaxm = ksdistm;

      }

      if (ksmaxp < 0.0) ksmaxp = 0.0;
      if (ksmaxm < 0.0) ksmaxm = 0.0;
      kuiper = ksmaxp + ksmaxm;

      if (kuiper < kuiperBest){
        kuiperBest = kuiper;
        shiftBest = shiftdist;
      }
    }

    ksshift.push_back(shiftBest);
    return kuiperBest;
  }

  
}//end hough namespace

////////////////////////////////////////////////////////////////////////