#include <CpaCrossCosmicIdAlg.h>

Classes
struct	BeamTime

struct	Config

struct	Fiducial

Public Member Functions
	CpaCrossCosmicIdAlg (const Config &config)

	CpaCrossCosmicIdAlg (const fhicl::ParameterSet &pset)

	~CpaCrossCosmicIdAlg ()

void	reconfigure (const Config &config)

std::pair< double, bool >	T0FromCpaStitching (detinfo::DetectorPropertiesData const &detProp, recob::Track t1, std::vector< recob::Track > tracks)

bool	CpaCrossCosmicId (detinfo::DetectorPropertiesData const &detProp, recob::Track track, std::vector< recob::Track > tracks, art::FindManyP< recob::Hit > hitAssoc)

Private Attributes
double	fCpaStitchDistance

double	fCpaStitchAngle

double	fCpaXDifference

double	fMinX

double	fMinY

double	fMinZ

double	fMaxX

double	fMaxY

double	fMaxZ

double	fBeamTimeMin

double	fBeamTimeMax

TPCGeoAlg	fTpcGeo

Detailed Description

Definition at line 38 of file sbndcode/sbndcode/CosmicId/Algs/CpaCrossCosmicIdAlg.h.

Constructor & Destructor Documentation

sbnd::CpaCrossCosmicIdAlg::CpaCrossCosmicIdAlg ( const Config & config )

Definition at line 7 of file sbndcode/sbndcode/CosmicId/Algs/CpaCrossCosmicIdAlg.cc.

                                                             {
 
   this->reconfigure(config);
 
 }

sbnd::CpaCrossCosmicIdAlg::CpaCrossCosmicIdAlg ( const fhicl::ParameterSet & pset )

inline

Definition at line 102 of file sbndcode/sbndcode/CosmicId/Algs/CpaCrossCosmicIdAlg.h.

102 :

103 CpaCrossCosmicIdAlg(fhicl::Table<Config>(pset, {})()) {}

sbnd::CpaCrossCosmicIdAlg::CpaCrossCosmicIdAlg

CpaCrossCosmicIdAlg(const Config &config)

Definition: sbndcode/sbndcode/CosmicId/Algs/CpaCrossCosmicIdAlg.cc:7

sbnd::CpaCrossCosmicIdAlg::~CpaCrossCosmicIdAlg ( )

Definition at line 19 of file sbndcode/sbndcode/CosmicId/Algs/CpaCrossCosmicIdAlg.cc.

                                          {
 
 }

Member Function Documentation

bool sbnd::CpaCrossCosmicIdAlg::CpaCrossCosmicId	(	detinfo::DetectorPropertiesData const &	detProp,
		recob::Track	track,
		std::vector< recob::Track >	tracks,
		art::FindManyP< recob::Hit >	hitAssoc
	)

Definition at line 115 of file sbndcode/sbndcode/CosmicId/Algs/CpaCrossCosmicIdAlg.cc.

                                                                                                                            {
 
   // Sort tracks by tpc
   std::vector<recob::Track> tpcTracksTPC0;
   std::vector<recob::Track> tpcTracksTPC1;
   // Loop over the tpc tracks
   for(auto const& tpcTrack : tracks){
     // Work out where the associated wire hits were detected
     std::vector<art::Ptr<recob::Hit>> hits = hitAssoc.at(tpcTrack.ID());
     int tpc = fTpcGeo.DetectedInTPC(hits);
     double startX = tpcTrack.Start().X();
     double endX = tpcTrack.End().X();
     if(tpc == 0 && !(startX>0 || endX>0)) tpcTracksTPC0.push_back(tpcTrack);
     else if(tpc == 1 && !(startX<0 || endX<0)) tpcTracksTPC1.push_back(tpcTrack);
   }
 
   std::vector<art::Ptr<recob::Hit>> hits = hitAssoc.at(track.ID());
   int tpc = fTpcGeo.DetectedInTPC(hits);
 
   double stitchTime = -99999;
   bool stitchExit = false;
   // Try to match tracks from CPA crossers
   if(tpc == 0){
     std::pair<double, bool> stitchResults = T0FromCpaStitching(detProp, track, tpcTracksTPC1);
     stitchTime = stitchResults.first;
     stitchExit = stitchResults.second;
   }
   else if(tpc == 1){
     std::pair<double, bool> stitchResults = T0FromCpaStitching(detProp, track, tpcTracksTPC0);
     stitchTime = stitchResults.first;
     stitchExit = stitchResults.second;
   }
 
   // If tracks are stitched, get time and remove any outside of beam window
   if(stitchTime != -99999 && (stitchTime < fBeamTimeMin || stitchTime > fBeamTimeMax || stitchExit)) return true;
   
   return false;
 
 }

void sbnd::CpaCrossCosmicIdAlg::reconfigure ( const Config & config )

Definition at line 24 of file sbndcode/sbndcode/CosmicId/Algs/CpaCrossCosmicIdAlg.cc.

                                                          {
 
   fCpaStitchDistance = config.CpaStitchDistance(); 
   fCpaStitchAngle = config.CpaStitchAngle();
   fCpaXDifference = config.CpaXDifference();
   fMinX = config.FiducialCuts().MinX(); 
   fMinY = config.FiducialCuts().MinY(); 
   fMinZ = config.FiducialCuts().MinZ(); 
   fMaxX = config.FiducialCuts().MaxX(); 
   fMaxY = config.FiducialCuts().MaxY(); 
   fMaxZ = config.FiducialCuts().MaxZ();
   fBeamTimeMin = config.BeamTimeLimits().BeamTimeMin();
   fBeamTimeMax = config.BeamTimeLimits().BeamTimeMax();
 
   return;
 }

std::pair< double, bool > sbnd::CpaCrossCosmicIdAlg::T0FromCpaStitching	(	detinfo::DetectorPropertiesData const &	detProp,
		recob::Track	t1,
		std::vector< recob::Track >	tracks
	)

Definition at line 42 of file sbndcode/sbndcode/CosmicId/Algs/CpaCrossCosmicIdAlg.cc.

                                                                                                               {
   
   std::vector<std::pair<double, std::pair<double, bool>>> matchCandidates;
   double matchedTime = -99999;
   std::pair<double, bool> returnVal = std::make_pair(matchedTime, false);
 
   TVector3 trk1Front = t1.Vertex<TVector3>();
   TVector3 trk1Back = t1.End<TVector3>();
   double closestX1 = std::min(std::abs(trk1Front.X()), std::abs(trk1Back.X()));
 
   // Loop over all tracks in other TPC
   for(auto & track : tracks){
 
     TVector3 trk2Front = track.Vertex<TVector3>();
     TVector3 trk2Back = track.End<TVector3>();
     double closestX2 = std::min(std::abs(trk2Front.X()), std::abs(trk2Back.X()));
 
     // Try to match if their ends have similar x positions
     if(std::abs(closestX1-closestX2) > fCpaXDifference) continue;
 
     // Find which point is closest to CPA
     TVector3 t1Pos = trk1Front;
     TVector3 t1PosEnd = trk1Back;
     TVector3 t1Dir = t1.VertexDirection<TVector3>();
     if(std::abs(trk1Back.X()) == closestX1){ 
       t1Pos = trk1Back;
       t1PosEnd = trk1Front;
       t1Dir = t1.EndDirection<TVector3>();
     }
 
     TVector3 t2Pos = trk2Front;
     TVector3 t2PosEnd = trk2Back;
     TVector3 t2Dir = track.VertexDirection<TVector3>();
     if(std::abs(trk2Back.X()) == closestX2){ 
       t2Pos = trk2Back;
       t2PosEnd = trk2Front;
       t2Dir = track.EndDirection<TVector3>();
     }
 
     // Calculate the angle between the tracks
     double trkCos = std::abs(t1Dir.Dot(t2Dir));
     // Calculate the distance between the tracks at the middle of the CPA
     t1Pos[0] = 0.;
     t2Pos[0] = 0.;
     double dist = (t1Pos-t2Pos).Mag();
 
     // Does the track enter and exit the fiducial volume when merged?
     geo::Point_t mergeStart {t1PosEnd.X(), t1PosEnd.Y(), t1PosEnd.Z()};
     geo::Point_t mergeEnd {t2PosEnd.X(), t2PosEnd.Y(), t2PosEnd.Z()};
     bool exits = false;
     if(!fTpcGeo.InFiducial(mergeStart, fMinX, fMinY, fMinZ, fMaxX, fMaxY, fMaxZ) 
        && !fTpcGeo.InFiducial(mergeEnd, fMinX, fMinY, fMinZ, fMaxX, fMaxY, fMaxZ)) exits = true;
 
     // If the distance and angle are within the acceptable limits then record candidate
     if(dist < fCpaStitchDistance && trkCos > cos(TMath::Pi() * fCpaStitchAngle / 180.)){ 
       matchCandidates.push_back(std::make_pair(trkCos, std::make_pair(closestX1, exits)));
     }
   }
 
   // Choose the candidate with the smallest angle
   if(matchCandidates.size() > 0){
     std::sort(matchCandidates.begin(), matchCandidates.end(), [](auto& left, auto& right){
               return left.first < right.first;});
     double shiftX = matchCandidates[0].second.first;
     matchedTime = -((shiftX - fTpcGeo.CpaWidth())/detProp.DriftVelocity()); //subtract CPA width
     returnVal = std::make_pair(matchedTime, matchCandidates[0].second.second);
   }
 
   return returnVal;
 }