trk::CosmicTrackUtilities Class Reference

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N21-01-23/TrackFit/CosmicTrackUtilities.h"

Public Member Functions
	CosmicTrackUtilities ()
virtual	~CosmicTrackUtilities ()
void	reconfigure (fhicl::ParameterSet const &pset)
float	RangeMomentum (rb::Track const &track)
float	LengthCorrection (float const &length)
float	CosTheta (float const &dcosY, float const &magnitude)
float	Azimuth (float const &dcosX, float const &dcosZ, float const &magnitude)
std::pair< std::pair< float, float >, std::pair< TVector3, TVector3 > >	TrueLengthEndPoints (simb::MCParticle const &part)
std::map< geo::CellUniqueId, float >	TruePathLengthInCell (simb::MCParticle const &part, std::vector< sim::FLSHit > const &flsHits)
void	FindTriCells (rb::Track const &track, std::map< std::pair< uint32_t, uint32_t >, bool > &planeCellMap)
float	TriCellPathLength (rb::Track const &track, uint32_t const &plane, double const &cellWidth, geo::View_t const &view)
float	PathLengthInCell (zBoundMap const &zBounds, TVector3 const &recoHitLoc, std::pair< uint32_t, uint32_t > const &pc)
zBoundMap	MakeZBoundaryMap (std::vector< TVector3 > const &trajectory)
zBounds	ZBounds (zBoundMap const &bounds, double const &hitZ)
std::vector< float >	TrackResiduals (rb::Track const &track)
bool	PointInDetector (TVector3 const &point)

Private Member Functions
void	FindZBoundaries ()

Private Attributes
double	fEffectivedEdx
	effective dE/dx in MeV/cm More...
std::vector< double >	fPolCoeff
	the coefficients in an 8th degree polynomial More...
std::set< double >	fPlaneZBounds
	boundaries of each plane in the z direction More...
double	fPlaneZDepth
	extent of one plane (cell) in the z direction More...
art::ServiceHandle< geo::LiveGeometry >	fLiveGeom
art::ServiceHandle< geo::Geometry >	fGeom

Detailed Description

Definition at line 35 of file CosmicTrackUtilities.h.

Constructor & Destructor Documentation

trk::CosmicTrackUtilities::CosmicTrackUtilities ( )

explicit

Definition at line 24 of file CosmicTrackUtilities.cxx.

References reconfigure().

  {
    fhicl::ParameterSet pset;
    this->reconfigure(pset);

    return;
  }

trk::CosmicTrackUtilities::~CosmicTrackUtilities ( )

virtual

Definition at line 33 of file CosmicTrackUtilities.cxx.

  {
  }

Member Function Documentation

float trk::CosmicTrackUtilities::Azimuth	(	float const &	dcosX,
		float const &	dcosZ,
		float const &	magnitude
	)

Definition at line 108 of file CosmicTrackUtilities.cxx.

References std::atan2().

Referenced by upmuana::UpMuAnalysis::analyze(), trk::CosmicTrackAna::analyze(), trk::CosmicTrackAna::FillEventMCTruthHistograms(), trk::CosmicTrackAna::FillRecoInfo(), and trk::CosmicTrackAna::FillTrueInfo().

  {
    float azimuth  = std::atan2(dcosX/magnitude, -dcosZ/magnitude) * 180./TMath::Pi(); 
    if(azimuth < 0.) azimuth += 360.;

    return azimuth;
  }

float trk::CosmicTrackUtilities::CosTheta	(	float const &	dcosY,
		float const &	magnitude
	)

Definition at line 101 of file CosmicTrackUtilities.cxx.

Referenced by upmuana::UpMuAnalysis::analyze(), trk::CosmicTrackAna::FillEventMCTruthHistograms(), trk::CosmicTrackAna::FillRecoInfo(), and trk::CosmicTrackAna::FillTrueInfo().

  {
    return -dcosY/magnitude;
  }

void trk::CosmicTrackUtilities::FindTriCells	(	rb::Track const &	track,
		std::map< std::pair< uint32_t, uint32_t >, bool > &	planeCellMap
	)

Definition at line 254 of file CosmicTrackUtilities.cxx.

References rb::CellHit::Cell(), rb::Cluster::Cell(), make_syst_table_plots::h, rb::Cluster::NCell(), and rb::CellHit::Plane().

Referenced by trk::CosmicTrackAna::FillRecoInfo(), and trk::CosmicTrackAna::FillTrackHistograms().

  {

    planeCellMap.clear();

    // first make a plane cell map so that we can identify tri-cells.  We will then fill the input map
    // to return it back to the calling function
    std::map<uint32_t, std::set<uint32_t> > tempMap;

    for(size_t h = 0; h < track.NCell(); ++h)
      tempMap[track.Cell(h)->Plane()].insert(track.Cell(h)->Cell());


    std::pair<uint32_t, uint32_t> planeCell;

    for(auto const& pacitr : tempMap){
      
      planeCell.first = pacitr.first;
      // get the set for the current plane
      // sets are presorted, go c++
      auto cellSet = pacitr.second;
      
      // loop over the cells and look to see if there is 
      // a cell in this plane that has both its neighbors also
      // hit on this track
      for( auto citr : cellSet ){
        planeCell.second = citr;
        if( cellSet.find(citr - 1) != cellSet.end() &&
           cellSet.find(citr + 1) != cellSet.end() )
        planeCellMap[planeCell] = true;
        else
        planeCellMap[planeCell] = false;
      }
    }// end loop over the tempMap entries         

    return;
  }

void trk::CosmicTrackUtilities::FindZBoundaries ( )

private

Definition at line 469 of file CosmicTrackUtilities.cxx.

References geo::PlaneGeo::Cell(), fGeom, fPlaneZBounds, fPlaneZDepth, geo::CellGeo::GetCenter(), geo::CellGeo::HalfD(), geo::GeometryBase::NPlanes(), and geo::GeometryBase::Plane().

Referenced by MakeZBoundaryMap().

  {
    // loop over all the planes in the geometry and put their upper z boundary
    // into a set
    std::set<double> bounds;
    bounds.insert(0.);
    double xyz[3] = {0.};
    fPlaneZDepth  = 2.*fGeom->Plane(0)->Cell(0)->HalfD();

    for(size_t p = 0; p < fGeom->NPlanes(); ++p){
      fGeom->Plane(p)->Cell(0)->GetCenter(xyz);
      bounds.insert(xyz[2] + 0.5*fPlaneZDepth);
    }
    
    fPlaneZBounds.swap(bounds);
    
    return;
  }

float trk::CosmicTrackUtilities::LengthCorrection

(

float const &

length

)

Definition at line 83 of file CosmicTrackUtilities.cxx.

References plot_validation_datamc::c, fPolCoeff, and cet::pow().

Referenced by trk::CosmicTrackAna::FillRecoInfo(), trk::CosmicTrackAna::FillTrackMCTruthHistograms(), and RangeMomentum().

  {
    // not a lot of tracks have lengths greater than 
    // 2500 cm, so just use the correction for 2500 cm
    // with longer tracks
    float ln = (length < 2500.) ? length : 2500.;

    float correction = 1.;

    if(fPolCoeff.size() > 0){
      for(size_t c = 0; c < fPolCoeff.size(); ++c) 
        correction -= fPolCoeff[c]*std::pow(ln, c);
    }

    return 1./correction;
  }

zBoundMap trk::CosmicTrackUtilities::MakeZBoundaryMap

(

std::vector< TVector3 > const &

trajectory

)

Definition at line 489 of file CosmicTrackUtilities.cxx.

References dir, FindZBoundaries(), fPlaneZBounds, fPlaneZDepth, make_pair(), PointInDetector(), X, and Z.

Referenced by trk::CosmicTrackAna::FillTrackHistograms(), TrackResiduals(), and TruePathLengthInCell().

  {
    if(fPlaneZBounds.size() < 1) this->FindZBoundaries();

    zBoundMap zToPlaneBounds;

    double deltaZ = 0.;
    double lowZ   = 0.;
    double highZ  = 0.;

    auto litr = fPlaneZBounds.begin();
    auto hitr = fPlaneZBounds.begin();

    TVector3 dir;
    TVector3 lowBound;
    TVector3 highBound;

    // for(auto tp : trajectory) 
    //   LOG_VERBATIM("CosmicTrackUtilities") << "trajectory: " << tp.X() << " " << tp.Y() << " " << tp.Z();

    for(size_t tp = 0; tp < trajectory.size(); ++tp){
      
      // make sure the trajectory point is inside the detector
      if( !this->PointInDetector(trajectory[tp]) ) continue;

      // find the intersection of the trajectoryectory with the boundary
      // of the plane it represents
      if(tp < trajectory.size() - 1) dir = trajectory[tp+1] - trajectory[tp];
      else                           dir = trajectory[tp]   - trajectory[tp-1];

      // find the plane boundaries for this point
      litr = fPlaneZBounds.lower_bound(trajectory[tp].Z() - fPlaneZDepth);
      hitr = fPlaneZBounds.upper_bound(trajectory[tp].Z());
      
      // in the muon catcher, the high side and low side iterators could be the same
      if(hitr == litr){
        if(litr != fPlaneZBounds.begin() ) --litr;
        else ++hitr;
      }
      lowZ = *litr;
      highZ = *hitr;

      // LOG_VERBATIM("CosmicTrackUtilities") << tp << " " << trajectory[tp].Z() 
      //                                           << " low z: " << lowZ << " high z: " << highZ
      //                                           << "trajectory point: (" 
      //                                           << trajectory[tp].X() << " " << trajectory[tp].Y() << " " << trajectory[tp].Z() << ")"; 

      // make sure we don't go past the detector boundaries
      deltaZ = trajectory[tp].Z() - lowZ;
      lowBound.SetXYZ(trajectory[tp].X() - deltaZ*dir.X()/dir.Z(),
                      trajectory[tp].Y() - deltaZ*dir.Y()/dir.Z(),
                      lowZ);
      
      // LOG_VERBATIM("CosmicTrackUtilities") << "is bounded by " << deltaZ << " " << dir.X()/dir.Z() << " " << dir.Y()/dir.Z()
      //                                           << "(" << lowBound.X() << " " << lowBound.Y() << " " << lowBound.Z() << ")";

      deltaZ = highZ - trajectory[tp].Z();
      highBound.SetXYZ(trajectory[tp].X() + deltaZ*dir.X()/dir.Z(),
                       trajectory[tp].Y() + deltaZ*dir.Y()/dir.Z(),
                       highZ);

      // LOG_VERBATIM("CosmicTrackUtilities") << "and " << deltaZ << " (" 
      //                                           << highBound.X() << " " << highBound.Y() << " " << highBound.Z() << ")";

      zToPlaneBounds[trajectory[tp].Z()] = std::make_pair(lowBound, highBound);

    }

    return zToPlaneBounds;
  }

float trk::CosmicTrackUtilities::PathLengthInCell	(	zBoundMap const &	zBounds,
		TVector3 const &	recoHitLoc,
		std::pair< uint32_t, uint32_t > const &	pc
	)

Definition at line 317 of file CosmicTrackUtilities.cxx.

References std::abs(), geo::PlaneGeo::Cell(), d, Dot(), fGeom, geo::CellGeo::GetCenter(), geo::CellGeo::HalfD(), geo::CellGeo::HalfL(), geo::CellGeo::HalfW(), geo::kY, Mag(), geo::GeometryBase::Plane(), demo5::surf, Unit(), geo::PlaneGeo::View(), and ZBounds().

Referenced by trk::CosmicTrackAna::FillTrackHistograms(), and TruePathLengthInCell().

  {
    double xyz[3] = {0.};
    fGeom->Plane(pc.first)->Cell(pc.second)->GetCenter(xyz);
    TVector3 cellCenter = TVector3(xyz);
    TVector3 cellExtent = TVector3(fGeom->Plane(pc.first)->Cell(pc.second)->HalfW(),
                                   fGeom->Plane(pc.first)->Cell(pc.second)->HalfL(),
                                   fGeom->Plane(pc.first)->Cell(pc.second)->HalfD());
    if(fGeom->Plane(pc.first)->View() == geo::kY)
      cellExtent = TVector3(fGeom->Plane(pc.first)->Cell(pc.second)->HalfL(),
                            fGeom->Plane(pc.first)->Cell(pc.second)->HalfW(),
                            fGeom->Plane(pc.first)->Cell(pc.second)->HalfD());

    // find the trajectory points of the track that bound this
    // hit in z and determine where the hit is expected to be
    // based on the direction of the track between the trajectory 
    // points
    auto bounds = this->ZBounds(zBounds, xyz[2]);
    
    TVector3 lowerPlaneBound = bounds.first;
    TVector3 upperPlaneBound = bounds.second;
    
    // now determine where the trajectory intersects the cell boundaries
    // evaluate the point where the line intersects the plane of each 
    // boundary and see if that point is within the cell
    // see http://en.wikipedia.org/wiki/Line–plane_intersection

    float pathLength = 0.;
    float d = 0.;
    // all vectors are in order of -x, +x, -y, +y, -z, +z

    // determine a point on each plane bounding the cell
    // define the edges of the cell using its center and extent values
    auto pointsList = {TVector3(cellCenter.X() - cellExtent.X(), cellCenter.Y(),                  cellCenter.Z()),
                       TVector3(cellCenter.X() + cellExtent.X(), cellCenter.Y(),                  cellCenter.Z()),
                       TVector3(cellCenter.X(),                  cellCenter.Y() - cellExtent.Y(), cellCenter.Z()),
                       TVector3(cellCenter.X(),                  cellCenter.Y() + cellExtent.Y(), cellCenter.Z()),
                       TVector3(cellCenter.X(),                  cellCenter.Y(),                  cellCenter.Z() - cellExtent.Z()),
                       TVector3(cellCenter.X(),                  cellCenter.Y(),                  cellCenter.Z() + cellExtent.Z())};
    std::vector<TVector3> pointsInPlanes(pointsList);

    auto normalList = { TVector3(-1.,  0.,  0.),
                                          TVector3( 1.,  0.,  0.),
                                          TVector3( 0., -1.,  0.),
                                          TVector3( 0.,  1.,  0.),
                                          TVector3( 0.,  0., -1.),
                                          TVector3( 0.,  0.,  1.) };
    std::vector<TVector3> normalVects(normalList);

    std::vector<TVector3> intersections;
    TVector3 intersect;

    // vector describing the boundaries of this cell
    auto boundList = {pointsInPlanes[0].X() - 1.e-4, pointsInPlanes[1].X() + 1.e-4, 
                      pointsInPlanes[2].Y() - 1.e-4, pointsInPlanes[3].Y() + 1.e-4,
                      pointsInPlanes[4].Z() - 1.e-4, pointsInPlanes[5].Z() + 1.e-4};
    std::vector<double> cellBounds(boundList);

    // LOG_VERBATIM("CosmicTrackUtilities") << "reco hit location: " 
    //                                   << recoHitLoc.X() << " " << recoHitLoc.Y() << " " << recoHitLoc.Z()
    //                                   <<"\ncell center: "
    //                                   << cellCenter.X() << " " << cellCenter.Y() << " " << cellCenter.Z()
    //                                   << "\ncell boundaries: " 
    //                                   << cellBounds[0] << " " << cellBounds[1] << " " 
    //                                   << cellBounds[2] << " " << cellBounds[3] << " " 
    //                                   << cellBounds[4] << " " << cellBounds[5];

    // check that the reco hit is between the bounds, return 0 if not
    if(recoHitLoc.X() < cellBounds[0] || recoHitLoc.X() > cellBounds[1] ||
       recoHitLoc.Y() < cellBounds[2] || recoHitLoc.Y() > cellBounds[3] ||
       recoHitLoc.Z() < cellBounds[4] || recoHitLoc.Z() > cellBounds[5])
      return 0.;

    // the recoHitLoc is the determined 3D position for the given hit
    // the plane bounds are the 3D positions where the track crosses the 
    // plane containing this hit in the upstream and downstream locations
    
    // the vector from the upstream and downstream locations gives the direction
    // from one side of the plane to the other.  Grab the location of the hit
    // and then walk along the direction vector until we step out of the cell
    // be sure to check the boundary of the cell in both the z direction and 
    // in the width of the cell
    TVector3 traj = (upperPlaneBound - lowerPlaneBound).Unit();

    // if(traj.Mag() == 0.)
    //   LOG_VERBATIM("CosmicTrackUtilities") << "upper bound: " 
    //                                     << upperPlaneBound.X() << " " 
    //                                     << upperPlaneBound.Y() << " " 
    //                                     << upperPlaneBound.Z() << "\nlower bound: " 
    //                                     << lowerPlaneBound.X() << " " 
    //                                     << lowerPlaneBound.Y() << " " 
    //                                     << lowerPlaneBound.Z() << "\ntrajectory: " 
    //                                     << traj.X() << " " << traj.Y() << " " << traj.Z();

    // loop over the different planes bounding the cell and find the intersection points
    // go in the reverse order and remove any surfaces where the intersection is outside
    // of the bounds of the cell
    for(size_t surf = 0; surf < normalVects.size(); ++surf){

      // check that you can do the calculation, if not, there is no 
      // intersectoin with that surface and we remove it.
      if( std::abs(traj.Dot(normalVects[surf])) > 0. ){
        d = (pointsInPlanes[surf] - lowerPlaneBound).Dot(normalVects[surf])/traj.Dot(normalVects[surf]);
        intersect = d*traj + lowerPlaneBound;

        // LOG_VERBATIM("CosmicTrackUtilities") << "intersection: "
        //                                   << intersect.X() << " "
        //                                   << intersect.Y() << " "
        //                                   << intersect.Z() << " surf "
        //                                   << surf << " "
        //                                   << normalVects[surf].X() << " "
        //                                   << normalVects[surf].Y() << " "
        //                                   << normalVects[surf].Z() << " "
        //                                   << d << " "
        //                                   << d*traj.X() << " "
        //                                   << d*traj.Y() << " "
        //                                   << d*traj.Z();
        
        // check that the intersection is between the bounds of the cell
        // if not, remove the intersection, allow the boundaries to be fuzzy by
        // 0.1 cm to account for rounding errors
        if(intersect.X() > cellBounds[0] && intersect.X() < cellBounds[1] &&
           intersect.Y() > cellBounds[2] && intersect.Y() < cellBounds[3] &&
           intersect.Z() > cellBounds[4] && intersect.Z() < cellBounds[5])
        intersections.push_back(intersect);
        
      }// end if the denominator is non-zero      
    }// end loop over surfaces

    // for(auto i : intersections) LOG_VERBATIM("CosmicTrackUtilitites") << "intersection point: " 
    //                                                                << i.X() << " " << i.Y() << " " << i.Z();

    // better have just 2 intersections at this point
    // possible to have only one, if the track appears to stop in the cell, then use the reco hit location as 
    // one end point
    if(intersections.size() > 2){
      throw cet::exception("CosmicTrackUtilities") << "too many cell surfaces intersected "
                                                   << "by trajectory: " << intersections.size();
    }    
    else if(intersections.size() == 1)
      pathLength = (intersections[0] - recoHitLoc).Mag();
    else if(intersections.size() < 1)
      pathLength = 0.;
    else
      pathLength = (intersections[0] - intersections[1]).Mag();
      
    return pathLength;
  }

bool trk::CosmicTrackUtilities::PointInDetector

(

TVector3 const &

point

)

Definition at line 119 of file CosmicTrackUtilities.cxx.

References geo::GeometryBase::DetId(), geo::LiveGeometry::DistanceToEdgeInMC(), geo::LiveGeometry::DistToEdgeXY(), geo::LiveGeometry::DistToEdgeZ(), fGeom, fLiveGeom, geo::LiveGeometry::InMuonCatcher(), novadaq::cnv::kNEARDET, and std::min().

Referenced by MakeZBoundaryMap(), and TrueLengthEndPoints().

  {
    bool inDetector = true;

    // make sure the trajectory point is inside the detector
    if(fGeom->DetId() != novadaq::cnv::kNEARDET){
      if( std::min(fLiveGeom->DistToEdgeXY(point), fLiveGeom->DistToEdgeZ(point)) < 0.)
        inDetector = false;
    }
    else{
      if(fLiveGeom->InMuonCatcher(point)){
        int wall = 0;
        if( fLiveGeom->DistanceToEdgeInMC(point, wall) < 0 )
        inDetector = false;
      }
      else if( std::min(fLiveGeom->DistToEdgeXY(point), fLiveGeom->DistToEdgeZ(point)) < 0.)
        inDetector = false;
    }

    return inDetector;
  }

float trk::CosmicTrackUtilities::RangeMomentum

(

rb::Track const &

track

)

Definition at line 59 of file CosmicTrackUtilities.cxx.

References dist, fEffectivedEdx, LengthCorrection(), Mag(), confusionMatrixTree::t, rb::Track::TotalLength(), and rb::Track::Trajectory().

Referenced by trk::CosmicTrackAna::FillRecoInfo().

  {
    // get the momentum from range.  The value of dE/dx is given
    // in MeV/cm, without a density factor.  The value was determined
    // from data/MC comparisons such that at the longest track lengths,
    // the fractional difference in data and MC was centered at zero

    float dist     = 0.;
    float rangeMom = 0.;

    std::vector<TVector3> const& traj = track.Trajectory();
    for(size_t t = traj.size()-1; t > 0; --t){
      dist      = (traj[t] - traj[t-1]).Mag();
      rangeMom += dist * fEffectivedEdx * 1.e-3; // 1.e-3 converts MeV to GeV
    }
    
    // The exponential factor accounts for higher dE/dx for shorter, ie
    // less energetic, muons.  The values were determined fitting a 
    // histogram of Delta p/p vs track length
    return rangeMom*this->LengthCorrection(track.TotalLength()); 

  }

void trk::CosmicTrackUtilities::reconfigure

(

fhicl::ParameterSet const &

pset

)

Definition at line 38 of file CosmicTrackUtilities.cxx.

References fEffectivedEdx, fPolCoeff, and fhicl::ParameterSet::get().

Referenced by CosmicTrackUtilities(), and trk::CosmicTrackAna::reconfigure().

  {

    auto pol = { 0.584713,
      -0.00369301,
      1.17201e-05,
      -2.3046e-08,
      2.9472e-11,
      -2.47581e-14,
      1.35044e-17,
      -4.59248e-21,
      8.82865e-25,
      -7.31569e-29 };        

    fEffectivedEdx      = pset.get< double                    >("EffectivedEdx",            2.122);
    fPolCoeff                 = pset.get< std::vector<double> >("PolCoeff",                 pol  );

    return;
  }

std::vector< float > trk::CosmicTrackUtilities::TrackResiduals

(

rb::Track const &

track

)

Definition at line 585 of file CosmicTrackUtilities.cxx.

References plot_validation_datamc::c, rb::CellHit::Cell(), rb::Cluster::Cell(), fGeom, rb::RecoHit::IsCalibrated(), LOG_DEBUG, Mag(), MakeZBoundaryMap(), rb::Cluster::NCell(), rb::CellHit::Plane(), geo::GeometryBase::Plane(), rb::Cluster::RecoHit(), rb::Track::Trajectory(), rb::RecoHit::X(), rb::RecoHit::Y(), rb::RecoHit::Z(), and ZBounds().

Referenced by trk::CosmicTrackAna::FillRecoInfo(), and trk::CosmicTrackAna::FillTrackHistograms().

  {
    LOG_DEBUG("CosmicTrackUtilities") << "TrackResiduals";

    std::vector<float> residuals;
    residuals.reserve(track.NCell());

    TVector3 cross;
    TVector3 point;
    TVector3 lowerBound;
    TVector3 upperBound;
    double   xyz[3] = {0.};

    zBoundMap zToPlaneBounds = this->MakeZBoundaryMap(track.Trajectory());

    // Loop through both views and grab every CellHit in this track
    for(size_t c = 0; c < track.NCell(); ++c){  
      const art::Ptr<rb::CellHit> chit = track.Cell(c);
      
      // LOG_VERBATIM("CosmicTrackAna") << chit->Plane() << " " << chit->Cell() <<  "cell hit view is " << view;

      const rb::RecoHit rhit = track.RecoHit(chit);
      if(!rhit.IsCalibrated()){
          //mf::LogError("CosmicTrackAna") << "Not calibrated?!";
        continue;
      }

      fGeom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
      // find the trajectory points of the track that bound this
      // hit in z and determine where the hit is expected to be
      // based on the direction of the track between the trajectory 
      // points
      auto bounds = this->ZBounds(zToPlaneBounds, xyz[2]);

      lowerBound = bounds.first;
      upperBound = bounds.second;

      point = TVector3(rhit.X(), rhit.Y(), rhit.Z());
        
      cross  = (point - lowerBound).Cross(point - upperBound);
      residuals.push_back(cross.Mag()/(upperBound - lowerBound).Mag());

    }// end loop over cells in the track

    return residuals;
  }

float trk::CosmicTrackUtilities::TriCellPathLength	(	rb::Track const &	track,
		uint32_t const &	plane,
		double const &	cellWidth,
		geo::View_t const &	view
	)

Definition at line 294 of file CosmicTrackUtilities.cxx.

References std::abs(), geo::kX, and rb::Track::PlaneDirMap().

Referenced by trk::CosmicTrackAna::FillTrackHistograms().

  {
    float pathLength = -9999.;
    auto planeDirMap = track.PlaneDirMap();
    auto pdmitr      = planeDirMap.find(plane);

    if( pdmitr != planeDirMap.end() ){
      float dw   = (view == geo::kX) ? std::abs(pdmitr->second.X()) : std::abs(pdmitr->second.Y());
      pathLength = cellWidth/dw;
    }

    return pathLength;
  }

std::pair< std::pair< float, float >, std::pair< TVector3, TVector3 > > trk::CosmicTrackUtilities::TrueLengthEndPoints

(

simb::MCParticle const &

part

)

Definition at line 142 of file CosmicTrackUtilities.cxx.

References std::abs(), febshutoff_auto::end, fLiveGeom, demo0::length, Mag(), make_training::momentum, PointInDetector(), prev(), geo::LiveGeometry::ProjectedDistanceToEdge(), febshutoff_auto::start, and simb::MCParticle::Trajectory().

Referenced by trk::CosmicTrackAna::FillEventMCTruthHistograms(), and trk::CosmicTrackAna::FillTrueInfo().

  {
    TVector3 start, end, prev;
    bool     setStart = false;
    float    length   = 0.;
    float    momentum = 0.;

    // loop over the particle trajectory points and set the length and end points based on
    // those points that are in the detector
    for(auto const& traj : part.Trajectory() ){
       
      if( this->PointInDetector(traj.first.Vect()) ){
        if(!setStart){
          setStart = true;
          start.SetXYZ(traj.first.X(), traj.first.Y(), traj.first.Z());
          length += fLiveGeom->ProjectedDistanceToEdge(traj.first.Vect(), -traj.second.Vect().Unit());
          momentum = std::abs(traj.second.P());
          prev = start;
        }
        end.SetXYZ(traj.first.X(), traj.first.Y(), traj.first.Z());
        length += (end - prev).Mag();
        prev = end;
      }
    }

    return std::pair<std::pair<float, float>, 
                     std::pair<TVector3, TVector3> >(std::pair<float, float>(length, momentum),
                                                    std::pair<TVector3, TVector3>(start,end));
  } 

std::map< geo::CellUniqueId, float > trk::CosmicTrackUtilities::TruePathLengthInCell	(	simb::MCParticle const &	part,
		std::vector< sim::FLSHit > const &	flsHits
	)

Definition at line 173 of file CosmicTrackUtilities.cxx.

References geo::PlaneGeo::Cell(), e, fGeom, geo::CellGeo::Id(), demo0::length, PandAna.Demos.demo0::loc, LOG_DEBUG, LOG_WARNING, make_pair(), MakeZBoundaryMap(), min(), PathLengthInCell(), and geo::GeometryBase::Plane().

Referenced by trk::CosmicTrackAna::FillTrueInfo().

  {
    LOG_DEBUG("CosmicTrackUtilities") << "TruePathLengthInCell";
    std::map<geo::CellUniqueId, float> pathInCell;
    std::vector< TVector3 > trajPoints;
    std::vector< TVector3 > flsTrajectory;
    TVector3 point;
    float    length = 0.;

    // loop over the particle trajectory points and get the points that are in the detector
    // for(auto const& traj : part.Trajectory() ){    

    //   if( !this->PointInDetector(traj.first.Vect()) ) continue;
  
    //   trajPoints.push_back(TVector3(traj.first.X(), traj.first.Y(), traj.first.Z()));
    // }

    // the true trajectory likely does not have a point for every plane crossed in the detector.
    // get the trajectory by using the average FLS hit position in each plane
    std::map<std::pair<size_t, size_t>, TVector3>  flsHitCenters;
    std::map<size_t, std::pair<TVector3, double> > flsPlaneZPos;
    auto fpzpItr = flsPlaneZPos.begin();
    double cnt    = 0.;
    double xyz[3] = {0.};
    double loc[3] = {0.};
    for(auto fls : flsHits){

      cnt = 1.;
      loc[0] = fls.GetXAverage();
      loc[1] = fls.GetYAverage();
      loc[2] = fls.GetZAverage();
      fGeom->Plane(fls.GetPlaneID())->Cell(fls.GetCellID())->LocalToWorld(loc, xyz);
      point.SetXYZ(xyz[0], xyz[1], xyz[2]);
      flsHitCenters[std::make_pair(fls.GetPlaneID(), fls.GetCellID())] = point;

      // now make the map to get the average position in a plane mapped to the 
      // center z poisition of the plane
      fGeom->Plane(fls.GetPlaneID())->Cell(fls.GetCellID())->GetCenter(xyz);
      if(flsPlaneZPos.count(fls.GetPlaneID()) > 0){
        fpzpItr = flsPlaneZPos.find(fls.GetPlaneID());
        point  += fpzpItr->second.first;
        cnt    += fpzpItr->second.second;
      }

      flsPlaneZPos[fls.GetPlaneID()] = std::make_pair(point, cnt);

    } // end loop over fls hits

    // get the average 3D position at each plane
    for(auto pzp : flsPlaneZPos){
      point = pzp.second.first;
      cnt   = pzp.second.second;
      point *= 1./cnt;
      flsTrajectory.push_back(point);
    }

    auto zBounds = MakeZBoundaryMap(flsTrajectory);

    // loop over the FLS hits to get the planes and cells for this particle
    for(auto fhc : flsHitCenters){

      try{
        length = this->PathLengthInCell(zBounds, fhc.second, fhc.first);
      }
      catch(cet::exception &e){
        LOG_WARNING("CosmicTrackAna")
        << "caught exception\n"
        << e
        << "\n when trying to calculate tracklength in cell.  Set it to a std::numeric_limits<float>::min()";
        length = std::numeric_limits<float>::min();
      }

      pathInCell[fGeom->Plane(fhc.first.first)->Cell(fhc.first.second)->Id()] = length;

    } // end loop over fls hits

    return pathInCell;
  }

zBounds trk::CosmicTrackUtilities::ZBounds	(	zBoundMap const &	bounds,
		double const &	hitZ
	)

Definition at line 561 of file CosmicTrackUtilities.cxx.

References fPlaneZDepth.

Referenced by PathLengthInCell(), and TrackResiduals().

  {
    // for(auto itr = bounds.begin(); itr != bounds.end(); ++itr)
    //   LOG_VERBATIM("CosmicTrackUtilities") << "hit z: " << hitZ << " " << itr->first << " lower "
    //                                     << itr->second.first.X() << " " 
    //                                     << itr->second.first.Y() << " " 
    //                                     << itr->second.first.Z() << " upper " 
    //                                     << itr->second.second.X() << " " 
    //                                     << itr->second.second.Y() << " " 
    //                                     << itr->second.second.Z();

    zBoundMap::const_iterator lbitr = bounds.lower_bound(hitZ - fPlaneZDepth);

    // if we did not find an iterator, then the z position passed is
    // off the edge of the track, use the last entry in the map
    if(lbitr == bounds.end() )
      return zBounds(bounds.rbegin()->second.first, bounds.rbegin()->second.second);

    return zBounds(lbitr->second.first, lbitr->second.second);

  }

Member Data Documentation

double trk::CosmicTrackUtilities::fEffectivedEdx

private

effective dE/dx in MeV/cm

Definition at line 79 of file CosmicTrackUtilities.h.

Referenced by RangeMomentum(), and reconfigure().

art::ServiceHandle<geo::Geometry> trk::CosmicTrackUtilities::fGeom

private

Definition at line 85 of file CosmicTrackUtilities.h.

Referenced by FindZBoundaries(), PathLengthInCell(), PointInDetector(), TrackResiduals(), and TruePathLengthInCell().

art::ServiceHandle<geo::LiveGeometry> trk::CosmicTrackUtilities::fLiveGeom

private

Definition at line 84 of file CosmicTrackUtilities.h.

Referenced by PointInDetector(), and TrueLengthEndPoints().

std::set<double> trk::CosmicTrackUtilities::fPlaneZBounds

private

boundaries of each plane in the z direction

Definition at line 81 of file CosmicTrackUtilities.h.

Referenced by FindZBoundaries(), and MakeZBoundaryMap().

double trk::CosmicTrackUtilities::fPlaneZDepth

private

extent of one plane (cell) in the z direction

Definition at line 82 of file CosmicTrackUtilities.h.

Referenced by FindZBoundaries(), MakeZBoundaryMap(), and ZBounds().

std::vector<double> trk::CosmicTrackUtilities::fPolCoeff

private

the coefficients in an 8th degree polynomial

Definition at line 80 of file CosmicTrackUtilities.h.

Referenced by LengthCorrection(), and reconfigure().

---

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

• /cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N21-01-23/TrackFit/CosmicTrackUtilities.h
• /cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N21-01-23/TrackFit/CosmicTrackUtilities.cxx