NDRecoFxs.cxx

Go to the documentation of this file.

#include "NDReco/NDRecoFxs.h"

#include "Calibrator/Calibrator.h"
#include "Geometry/Geometry.h"
#include "Geometry/LiveGeometry.h"
#include "GeometryObjects/Geo.h"
#include "GeometryObjects/CellGeo.h"
#include "GeometryObjects/PlaneGeo.h"
#include "RecoBase/RecoHit.h"

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

#include "TMath.h"
#include "TVector3.h"


namespace ndreco
{

  //-------------------------------------------------------
  NDRecoFxs::NDRecoFxs() //constructor
  {
  }

  //-------------------------------------------------------
  NDRecoFxs::~NDRecoFxs() //destructor
  {
  }

  //-------------------------------------------------------

  float NDRecoFxs::getTrackScatt(const art::Ptr<rb::Track> track, TVector3 vert, float dist) {
    double tscattmax=0, tscatttot=0, tscatt;
    TVector3 traj1, traj2, dir, olddir;

    for(unsigned int p=0;p<track->NTrajectoryPoints()-1;p++) {
      traj1 = track->TrajectoryPoint(p);
      traj2 = track->TrajectoryPoint(p+1);
      dir = (traj2-traj1).Unit();
      if((dir-olddir).Mag()<1e-7) continue;
      if(p>0&&dir.Mag()<=1.0&&olddir.Mag()<=1.0&&(traj1-vert).Mag()>dist&&(traj2-vert).Mag()>dist) {
        tscatt = std::acos(dir.Dot(olddir)) * 180 / TMath::Pi();
        if(tscatt > tscattmax) tscattmax = tscatt;
        tscatttot += tscatt;
      }
      olddir = dir;
    }
    return tscatttot;  // calculated both maximum and total scatter, total seems more useful, ignore max for now
  }

  //-------------------------------------------------------
  // crude dE/dx: dE is GeV of hits, dx is a track length.  Complicated by excluding region within dist cm of vertex
  float NDRecoFxs::getTrackDedx(const art::Ptr<rb::Track> track, TVector3 vert, float dist) {
    art::ServiceHandle<geo::Geometry> geom;
    float gev = 0, length = 0;
    TVector3 traj1, traj2;

    art::PtrVector< rb::CellHit > const& trackHits = track->AllCells();

    rb::Cluster survivingCluster; // make cluster of hits further than dist away from vertex

    for (unsigned int trackHit = 0; trackHit < trackHits.size(); ++trackHit){ //fill cluster
      const geo::PlaneGeo* p = geom->Plane(trackHits[trackHit]->Plane());
      double xyz[3] = {0.,0.,0.}, junk = 0;
      p->Cell(trackHits[trackHit]->Cell())->GetCenter(xyz, junk);
      float dist2d = 0;
      if(trackHits[trackHit]->View()==geo::kX) dist2d = sqrt( pow(xyz[0]-vert.X(),2) + pow(xyz[2]-vert.Z(),2));
      else if(trackHits[trackHit]->View()==geo::kY) dist2d = sqrt( pow(xyz[1]-vert.Y(),2) + pow(xyz[2]-vert.Z(),2));
      if(dist2d>dist) survivingCluster.Add(trackHits[trackHit]);
    }

    if(survivingCluster.NCell()>0) gev = survivingCluster.TotalGeV(); // numerator is GeV of all these hits

    // dx isn't just old track length - dist.  To do right would check all places along all trajectory points.
    // we'll just worry about when one of the trajectory points are too close to vertex.  If a track segment went near vertex and it
    // didn't trigger formation of a new segment (ie a scatter) then it was probably unaffected by vertex activity, so ignore case where
    // line segment passes near vertex but both trajectory points are outside
    for(unsigned int p=0;p<track->NTrajectoryPoints()-1;p++) { // to get dx, look through trajectory points and see which are near vertex
      traj1 = track->TrajectoryPoint(p);
      traj2 = track->TrajectoryPoint(p+1);

      if((traj1-vert).Mag() > dist && (traj2-vert).Mag() > dist) length += (traj1-traj2).Mag(); // both end points of segment not near vertex
      else if((traj1-vert).Mag() < dist && (traj2-vert).Mag() < dist) continue; // both end points too close to vertex, ignore whole segment
      else {  // tricky case of one end point too close and other far away, want to add part of track length that is > dist from vertex
        TVector3 start, dir; // start at point near vertex, increment away from vertex until we roughly hit the edge of sphere of radius dist
        if((traj1-vert).Mag() < dist) {
          start = traj1;
          dir = (traj2-traj1).Unit();
        }
        else {
          start = traj2;
          dir = (traj1-traj2).Unit();
        }
        TVector3 point = start;
        bool keepgoing = true;
        while(keepgoing) { // go down line segment until we hit where we should stop
          point = point + dir;
          if( (point - vert).Mag() > (dist - 0.5)) { // this is where we'll stop
            length += ((traj1-traj2).Mag() - (point-start).Mag()); // add the length left, which is the part farther than dist from vertex
            keepgoing=false;
          }
          if((point-start).Mag() > (traj1-traj2).Mag() + 1) {
            keepgoing=false;
            std::cout << "SOMETHING WENT WRONG IN PION RECO DEDX CALC, TOO FAR DISTANCE" << std::endl;
          }
        } // end while loop
      } // end else
    } // end loop over trajectory points
    if(length <0) {
      std::cout << "SOMETHING WRONG IN PION RECO DEDX CALC, LENGTH = " << length << std::endl;
      return 0;
    }
    else if(length == 0) return 0;
    return gev/length;
  }

  //-----------------------------------------------------------------------------------
  float NDRecoFxs::getTrackActivity(const art::Ptr<rb::Track> track, TVector3 vert, float dist, art::PtrVector< rb::CellHit > const& sliceHits, float activitydist) {
    art::ServiceHandle<geo::Geometry> geom;
    TVector3 trackend = track->Stop();
    //    if((track->Start-vert).Mag() > 40 && (track->Stop()-vert).Mag() < 40) trackend = track->Start(); // doesn't do much
    //    in future may want to check start and end positions vs start/end positions of other tracks in slice instead

    art::PtrVector< rb::CellHit > const& trackHits = track->AllCells();
    rb::Cluster activityCluster;
    for (unsigned int sliceHit = 0; sliceHit < sliceHits.size(); ++sliceHit){
      bool match = 0;
      for (unsigned int trackHit = 0; trackHit < trackHits.size(); ++trackHit){
        match = (sliceHits[sliceHit] == trackHits[trackHit]);
        if (match) break;
      }//End of loop over hits in track

      const geo::PlaneGeo* p = geom->Plane(sliceHits[sliceHit]->Plane());
      double xyz[3] = {0.,0.,0.}, junk = 0;
      p->Cell(sliceHits[sliceHit]->Cell())->GetCenter(xyz, junk);
      float disttovert = 0, disttoend = 0;
      if(sliceHits[sliceHit]->View()==geo::kX) {
        disttovert = sqrt( pow(xyz[0]-vert.X(),2) + pow(xyz[2]-vert.Z(),2));
        disttoend = sqrt( pow(xyz[0]-trackend.X(),2) + pow(xyz[2]-trackend.Z(),2));
      }
      else if(sliceHits[sliceHit]->View()==geo::kY) {
        disttovert = sqrt( pow(xyz[1]-vert.Y(),2) + pow(xyz[2]-vert.Z(),2));
        disttoend = sqrt( pow(xyz[1]-trackend.Y(),2) + pow(xyz[2]-trackend.Z(),2));
      }
      if (!match && disttovert>dist && disttoend<activitydist) activityCluster.Add(sliceHits[sliceHit]);
    }//End of loop over hits in the slice
    if(activityCluster.NCell()>0) return activityCluster.TotalGeV();
    else return 0;
  }

  //--------------------------------------------------------------------------------------

  float NDRecoFxs::getTrackProximity(const std::vector< art::Ptr<rb::Track> > sliceTracks, int trknum, TVector3 vert, float dist) {

    art::Ptr<rb::Track> first = sliceTracks[trknum];

    TVector3 prox, min = {9999, 9999, 9999}, start1, start2, end1, end2;

    for(uint i=0; i<sliceTracks.size()-1; i++) {

      if((int)i==trknum) continue;
      
      art::Ptr<rb::Track> second = sliceTracks[i];

      start1 = first->Start();
      start2 = second->Start();

      end1 = first->Stop();
      end2 = second->Stop();

      // do all four cases
      // start with start to start

      prox = start1-start2;
      if(prox.Mag() < min.Mag() && (start1-vert).Mag() > dist && (start2-vert).Mag() > dist) min = prox;

      // start to end
      prox = start1-end2;
      if(prox.Mag() < min.Mag() && (start1-vert).Mag() > dist && (end2-vert).Mag() > dist)  min = prox;

      // end to start
      prox = end1-start2;
      if(prox.Mag() < min.Mag() && (end1-vert).Mag() > dist && (start2-vert).Mag() > dist)  min = prox;

      // now end to end
      prox = end1-end2;
      if(prox.Mag() < min.Mag() && (end1-vert).Mag() > dist && (end2-vert).Mag() > dist) min = prox;
    }
    if(min.Mag()>17000) return -1; //no conditions were ever met; min is still {9999,9999,9999)
    return min.Mag();
  }

  //--------------------------------------------------------------------------------------
  //dedx for prongs is simpler, only one line segment
  float NDRecoFxs::getProngDedx(const art::Ptr<rb::Prong> prong, TVector3 vert, float dist) {
    art::ServiceHandle<geo::Geometry> geom;
    float gev = 0, length = 0;
    TVector3 traj1, traj2;

    art::PtrVector< rb::CellHit > const& prongHits = prong->AllCells();

    rb::Cluster survivingCluster; // make cluster of hits further than dist away from vertex

    for (unsigned int prongHit = 0; prongHit < prongHits.size(); ++prongHit){ //fill cluster
      const geo::PlaneGeo* p = geom->Plane(prongHits[prongHit]->Plane());
      double xyz[3] = {0.,0.,0.}, junk = 0;
      p->Cell(prongHits[prongHit]->Cell())->GetCenter(xyz, junk);
      TVector3 cellcenter(xyz);
      float dist2d = 0;
      if(prongHits[prongHit]->View()==geo::kX) dist2d = sqrt( pow(xyz[0]-vert.X(),2) + pow(xyz[2]-vert.Z(),2));
      else if(prongHits[prongHit]->View()==geo::kY) dist2d = sqrt( pow(xyz[1]-vert.Y(),2) + pow(xyz[2]-vert.Z(),2));
      if(dist2d>dist) survivingCluster.Add(prongHits[prongHit]);
    }

    if(survivingCluster.NCell()>0) gev = survivingCluster.TotalGeV(); // numerator is GeV of all these hits

    // like for tracks but only one line segment
    traj1 = prong->Start();
    traj2 = prong->Start() + prong->Dir()*prong->TotalLength();

    if((traj1-vert).Mag() > dist && (traj2-vert).Mag() > dist) length += (traj1-traj2).Mag(); // both end points of segment not near vertex
    else if((traj1-vert).Mag() < dist && (traj2-vert).Mag() < dist) length = 0; // both end points too close to vertex, ignore whole segment
    else {  // tricky case of one end point too close and other far away, want to add part of track length that is > dist from vertex
      TVector3 start, dir; // start at point near vertex, increment away from vertex until we roughly hit the edge of sphere of radius dist
      if((traj1-vert).Mag() < dist) {
        start = traj1;
        dir = (traj2-traj1).Unit();
      }
      else {
        start = traj2;
        dir = (traj1-traj2).Unit();
      } 
      TVector3 point = start;
      bool keepgoing = true;
      while(keepgoing) { // go down line segment until we hit where we should stop
        point = point + dir;
        if( (point - vert).Mag() > (dist - 0.5)) { // this is where we'll stop
          length += ((traj1-traj2).Mag() - (point-start).Mag()); // add the length left, which is the part farther than dist from vertex
          keepgoing=false;
        }
        if((point-start).Mag() > (traj1-traj2).Mag() + 1) {
          keepgoing=false;
          std::cout << "SOMETHING WENT WRONG IN PION RECO DEDX CALC, TOO FAR DISTANCE" << std::endl;
        }       
      } // end while loop
    } // end else
    if(length <0) {
      std::cout << "SOMETHING WRONG IN PION RECO DEDX CALC, LENGTH = " << length << std::endl;
      return 0;
    }
    else if(length == 0) return 0;
    return gev/length;
  }

  //-----------------------------------------------------------------------------------
  float NDRecoFxs::getProngActivity(const art::Ptr<rb::Prong> prong, TVector3 vert, float dist, art::PtrVector< rb::CellHit > const& sliceHits, float activitydist) {
    art::ServiceHandle<geo::Geometry> geom;
    TVector3 prongend = prong->Start() + prong->Dir()*prong->TotalLength();

    art::PtrVector< rb::CellHit > const& prongHits = prong->AllCells();
    rb::Cluster activityCluster;
    for (unsigned int sliceHit = 0; sliceHit < sliceHits.size(); ++sliceHit){
      bool match = 0;
      for (unsigned int prongHit = 0; prongHit < prongHits.size(); ++prongHit){
        match = (sliceHits[sliceHit] == prongHits[prongHit]);
        if (match) break;
      }//End of loop over hits in prong

      const geo::PlaneGeo* p = geom->Plane(sliceHits[sliceHit]->Plane());
      double xyz[3] = {0.,0.,0.}, junk = 0;
      p->Cell(sliceHits[sliceHit]->Cell())->GetCenter(xyz, junk);
      float disttovert = 0, disttoend = 0;
      if(sliceHits[sliceHit]->View()==geo::kX) {
        disttovert = sqrt( pow(xyz[0]-vert.X(),2) + pow(xyz[2]-vert.Z(),2));
        disttoend = sqrt( pow(xyz[0]-prongend.X(),2) + pow(xyz[2]-prongend.Z(),2));
      }
      else if(sliceHits[sliceHit]->View()==geo::kY) {
        disttovert = sqrt( pow(xyz[1]-vert.Y(),2) + pow(xyz[2]-vert.Z(),2));
        disttoend = sqrt( pow(xyz[1]-prongend.Y(),2) + pow(xyz[2]-prongend.Z(),2));
      }
      if (!match && disttovert>dist && disttoend<activitydist) activityCluster.Add(sliceHits[sliceHit]);
    }//End of loop over hits in the slice
    if(activityCluster.NCell()>0) return activityCluster.TotalGeV();
    else return 0;
  }


  //--------------------------------------------------------------------------------------
  float NDRecoFxs::getProngProximity(const std::vector< art::Ptr<rb::Prong> > sliceProngs, int pngnum, TVector3 vert, float dist) {

    art::Ptr<rb::Prong> first = sliceProngs[pngnum];

    TVector3 prox, min = {9999, 9999, 9999}, start1, start2, end1, end2;

    for(uint i=0; i<sliceProngs.size()-1; i++) {

      if((int)i==pngnum) continue;
      
      art::Ptr<rb::Prong> second = sliceProngs[i];

      start1 = first->Start();
      start2 = second->Start();

      end1 = first->Start() + first->Dir()*first->TotalLength();
      end2 = second->Start() + first->Dir()*first->TotalLength();

      // do all four cases
      // start with start to start

      prox = start1-start2;
      if(prox.Mag() < min.Mag() && (start1-vert).Mag() > dist && (start2-vert).Mag() > dist) min = prox;

      // start to end
      prox = start1-end2;
      if(prox.Mag() < min.Mag() && (start1-vert).Mag() > dist && (end2-vert).Mag() > dist)  min = prox;

      // end to start
      prox = end1-start2;
      if(prox.Mag() < min.Mag() && (end1-vert).Mag() > dist && (start2-vert).Mag() > dist)  min = prox;

      // now end to end
      prox = end1-end2;
      if(prox.Mag() < min.Mag() && (end1-vert).Mag() > dist && (end2-vert).Mag() > dist) min = prox;
    }
    if(min.Mag()>17000) return -1; //no conditions were ever met; min is still {9999,9999,9999)
    return min.Mag();
  }


} //end namespace