Public Member Functions | Private Member Functions | Private Attributes | List of all members
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::LiveGeometryfLiveGeom
 
art::ServiceHandle< geo::GeometryfGeom
 

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().

25  {
27  this->reconfigure(pset);
28 
29  return;
30  }
void reconfigure(fhicl::ParameterSet const &pset)
trk::CosmicTrackUtilities::~CosmicTrackUtilities ( )
virtual

Definition at line 33 of file CosmicTrackUtilities.cxx.

34  {
35  }

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().

111  {
112  float azimuth = std::atan2(dcosX/magnitude, -dcosZ/magnitude) * 180./TMath::Pi();
113  if(azimuth < 0.) azimuth += 360.;
114 
115  return azimuth;
116  }
T atan2(T number)
Definition: d0nt_math.hpp:72
float trk::CosmicTrackUtilities::CosTheta ( float const &  dcosY,
float const &  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().

256  {
257 
258  planeCellMap.clear();
259 
260  // first make a plane cell map so that we can identify tri-cells. We will then fill the input map
261  // to return it back to the calling function
262  std::map<uint32_t, std::set<uint32_t> > tempMap;
263 
264  for(size_t h = 0; h < track.NCell(); ++h)
265  tempMap[track.Cell(h)->Plane()].insert(track.Cell(h)->Cell());
266 
267 
268  std::pair<uint32_t, uint32_t> planeCell;
269 
270  for(auto const& pacitr : tempMap){
271 
272  planeCell.first = pacitr.first;
273  // get the set for the current plane
274  // sets are presorted, go c++
275  auto cellSet = pacitr.second;
276 
277  // loop over the cells and look to see if there is
278  // a cell in this plane that has both its neighbors also
279  // hit on this track
280  for( auto citr : cellSet ){
281  planeCell.second = citr;
282  if( cellSet.find(citr - 1) != cellSet.end() &&
283  cellSet.find(citr + 1) != cellSet.end() )
284  planeCellMap[planeCell] = true;
285  else
286  planeCellMap[planeCell] = false;
287  }
288  }// end loop over the tempMap entries
289 
290  return;
291  }
Definition: event.h:19
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().

470  {
471  // loop over all the planes in the geometry and put their upper z boundary
472  // into a set
473  std::set<double> bounds;
474  bounds.insert(0.);
475  double xyz[3] = {0.};
476  fPlaneZDepth = 2.*fGeom->Plane(0)->Cell(0)->HalfD();
477 
478  for(size_t p = 0; p < fGeom->NPlanes(); ++p){
479  fGeom->Plane(p)->Cell(0)->GetCenter(xyz);
480  bounds.insert(xyz[2] + 0.5*fPlaneZDepth);
481  }
482 
483  fPlaneZBounds.swap(bounds);
484 
485  return;
486  }
double fPlaneZDepth
extent of one plane (cell) in the z direction
art::ServiceHandle< geo::Geometry > fGeom
::xsd::cxx::tree::bounds< char > bounds
Definition: Database.h:226
double HalfD() const
Definition: CellGeo.cxx:205
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
const PlaneGeo * Plane(unsigned int i) const
unsigned int NPlanes() const
std::set< double > fPlaneZBounds
boundaries of each plane in the z direction
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().

84  {
85  // not a lot of tracks have lengths greater than
86  // 2500 cm, so just use the correction for 2500 cm
87  // with longer tracks
88  float ln = (length < 2500.) ? length : 2500.;
89 
90  float correction = 1.;
91 
92  if(fPolCoeff.size() > 0){
93  for(size_t c = 0; c < fPolCoeff.size(); ++c)
94  correction -= fPolCoeff[c]*std::pow(ln, c);
95  }
96 
97  return 1./correction;
98  }
constexpr T pow(T x)
Definition: pow.h:75
std::vector< double > fPolCoeff
the coefficients in an 8th degree polynomial
length
Definition: demo0.py:21
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().

490  {
491  if(fPlaneZBounds.size() < 1) this->FindZBoundaries();
492 
493  zBoundMap zToPlaneBounds;
494 
495  double deltaZ = 0.;
496  double lowZ = 0.;
497  double highZ = 0.;
498 
499  auto litr = fPlaneZBounds.begin();
500  auto hitr = fPlaneZBounds.begin();
501 
502  TVector3 dir;
503  TVector3 lowBound;
504  TVector3 highBound;
505 
506  // for(auto tp : trajectory)
507  // LOG_VERBATIM("CosmicTrackUtilities") << "trajectory: " << tp.X() << " " << tp.Y() << " " << tp.Z();
508 
509  for(size_t tp = 0; tp < trajectory.size(); ++tp){
510 
511  // make sure the trajectory point is inside the detector
512  if( !this->PointInDetector(trajectory[tp]) ) continue;
513 
514  // find the intersection of the trajectoryectory with the boundary
515  // of the plane it represents
516  if(tp < trajectory.size() - 1) dir = trajectory[tp+1] - trajectory[tp];
517  else dir = trajectory[tp] - trajectory[tp-1];
518 
519  // find the plane boundaries for this point
520  litr = fPlaneZBounds.lower_bound(trajectory[tp].Z() - fPlaneZDepth);
521  hitr = fPlaneZBounds.upper_bound(trajectory[tp].Z());
522 
523  // in the muon catcher, the high side and low side iterators could be the same
524  if(hitr == litr){
525  if(litr != fPlaneZBounds.begin() ) --litr;
526  else ++hitr;
527  }
528  lowZ = *litr;
529  highZ = *hitr;
530 
531  // LOG_VERBATIM("CosmicTrackUtilities") << tp << " " << trajectory[tp].Z()
532  // << " low z: " << lowZ << " high z: " << highZ
533  // << "trajectory point: ("
534  // << trajectory[tp].X() << " " << trajectory[tp].Y() << " " << trajectory[tp].Z() << ")";
535 
536  // make sure we don't go past the detector boundaries
537  deltaZ = trajectory[tp].Z() - lowZ;
538  lowBound.SetXYZ(trajectory[tp].X() - deltaZ*dir.X()/dir.Z(),
539  trajectory[tp].Y() - deltaZ*dir.Y()/dir.Z(),
540  lowZ);
541 
542  // LOG_VERBATIM("CosmicTrackUtilities") << "is bounded by " << deltaZ << " " << dir.X()/dir.Z() << " " << dir.Y()/dir.Z()
543  // << "(" << lowBound.X() << " " << lowBound.Y() << " " << lowBound.Z() << ")";
544 
545  deltaZ = highZ - trajectory[tp].Z();
546  highBound.SetXYZ(trajectory[tp].X() + deltaZ*dir.X()/dir.Z(),
547  trajectory[tp].Y() + deltaZ*dir.Y()/dir.Z(),
548  highZ);
549 
550  // LOG_VERBATIM("CosmicTrackUtilities") << "and " << deltaZ << " ("
551  // << highBound.X() << " " << highBound.Y() << " " << highBound.Z() << ")";
552 
553  zToPlaneBounds[trajectory[tp].Z()] = std::make_pair(lowBound, highBound);
554 
555  }
556 
557  return zToPlaneBounds;
558  }
double fPlaneZDepth
extent of one plane (cell) in the z direction
std::map< double, zBounds > zBoundMap
bool PointInDetector(TVector3 const &point)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
Float_t Z
Definition: plot.C:38
TDirectory * dir
Definition: macro.C:5
Float_t X
Definition: plot.C:38
std::set< double > fPlaneZBounds
boundaries of each plane in the z direction
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().

320  {
321  double xyz[3] = {0.};
322  fGeom->Plane(pc.first)->Cell(pc.second)->GetCenter(xyz);
323  TVector3 cellCenter = TVector3(xyz);
324  TVector3 cellExtent = TVector3(fGeom->Plane(pc.first)->Cell(pc.second)->HalfW(),
325  fGeom->Plane(pc.first)->Cell(pc.second)->HalfL(),
326  fGeom->Plane(pc.first)->Cell(pc.second)->HalfD());
327  if(fGeom->Plane(pc.first)->View() == geo::kY)
328  cellExtent = TVector3(fGeom->Plane(pc.first)->Cell(pc.second)->HalfL(),
329  fGeom->Plane(pc.first)->Cell(pc.second)->HalfW(),
330  fGeom->Plane(pc.first)->Cell(pc.second)->HalfD());
331 
332  // find the trajectory points of the track that bound this
333  // hit in z and determine where the hit is expected to be
334  // based on the direction of the track between the trajectory
335  // points
336  auto bounds = this->ZBounds(zBounds, xyz[2]);
337 
338  TVector3 lowerPlaneBound = bounds.first;
339  TVector3 upperPlaneBound = bounds.second;
340 
341  // now determine where the trajectory intersects the cell boundaries
342  // evaluate the point where the line intersects the plane of each
343  // boundary and see if that point is within the cell
344  // see http://en.wikipedia.org/wiki/Line–plane_intersection
345 
346  float pathLength = 0.;
347  float d = 0.;
348  // all vectors are in order of -x, +x, -y, +y, -z, +z
349 
350  // determine a point on each plane bounding the cell
351  // define the edges of the cell using its center and extent values
352  auto pointsList = {TVector3(cellCenter.X() - cellExtent.X(), cellCenter.Y(), cellCenter.Z()),
353  TVector3(cellCenter.X() + cellExtent.X(), cellCenter.Y(), cellCenter.Z()),
354  TVector3(cellCenter.X(), cellCenter.Y() - cellExtent.Y(), cellCenter.Z()),
355  TVector3(cellCenter.X(), cellCenter.Y() + cellExtent.Y(), cellCenter.Z()),
356  TVector3(cellCenter.X(), cellCenter.Y(), cellCenter.Z() - cellExtent.Z()),
357  TVector3(cellCenter.X(), cellCenter.Y(), cellCenter.Z() + cellExtent.Z())};
358  std::vector<TVector3> pointsInPlanes(pointsList);
359 
360  auto normalList = { TVector3(-1., 0., 0.),
361  TVector3( 1., 0., 0.),
362  TVector3( 0., -1., 0.),
363  TVector3( 0., 1., 0.),
364  TVector3( 0., 0., -1.),
365  TVector3( 0., 0., 1.) };
366  std::vector<TVector3> normalVects(normalList);
367 
368  std::vector<TVector3> intersections;
369  TVector3 intersect;
370 
371  // vector describing the boundaries of this cell
372  auto boundList = {pointsInPlanes[0].X() - 1.e-4, pointsInPlanes[1].X() + 1.e-4,
373  pointsInPlanes[2].Y() - 1.e-4, pointsInPlanes[3].Y() + 1.e-4,
374  pointsInPlanes[4].Z() - 1.e-4, pointsInPlanes[5].Z() + 1.e-4};
375  std::vector<double> cellBounds(boundList);
376 
377  // LOG_VERBATIM("CosmicTrackUtilities") << "reco hit location: "
378  // << recoHitLoc.X() << " " << recoHitLoc.Y() << " " << recoHitLoc.Z()
379  // <<"\ncell center: "
380  // << cellCenter.X() << " " << cellCenter.Y() << " " << cellCenter.Z()
381  // << "\ncell boundaries: "
382  // << cellBounds[0] << " " << cellBounds[1] << " "
383  // << cellBounds[2] << " " << cellBounds[3] << " "
384  // << cellBounds[4] << " " << cellBounds[5];
385 
386  // check that the reco hit is between the bounds, return 0 if not
387  if(recoHitLoc.X() < cellBounds[0] || recoHitLoc.X() > cellBounds[1] ||
388  recoHitLoc.Y() < cellBounds[2] || recoHitLoc.Y() > cellBounds[3] ||
389  recoHitLoc.Z() < cellBounds[4] || recoHitLoc.Z() > cellBounds[5])
390  return 0.;
391 
392  // the recoHitLoc is the determined 3D position for the given hit
393  // the plane bounds are the 3D positions where the track crosses the
394  // plane containing this hit in the upstream and downstream locations
395 
396  // the vector from the upstream and downstream locations gives the direction
397  // from one side of the plane to the other. Grab the location of the hit
398  // and then walk along the direction vector until we step out of the cell
399  // be sure to check the boundary of the cell in both the z direction and
400  // in the width of the cell
401  TVector3 traj = (upperPlaneBound - lowerPlaneBound).Unit();
402 
403  // if(traj.Mag() == 0.)
404  // LOG_VERBATIM("CosmicTrackUtilities") << "upper bound: "
405  // << upperPlaneBound.X() << " "
406  // << upperPlaneBound.Y() << " "
407  // << upperPlaneBound.Z() << "\nlower bound: "
408  // << lowerPlaneBound.X() << " "
409  // << lowerPlaneBound.Y() << " "
410  // << lowerPlaneBound.Z() << "\ntrajectory: "
411  // << traj.X() << " " << traj.Y() << " " << traj.Z();
412 
413  // loop over the different planes bounding the cell and find the intersection points
414  // go in the reverse order and remove any surfaces where the intersection is outside
415  // of the bounds of the cell
416  for(size_t surf = 0; surf < normalVects.size(); ++surf){
417 
418  // check that you can do the calculation, if not, there is no
419  // intersectoin with that surface and we remove it.
420  if( std::abs(traj.Dot(normalVects[surf])) > 0. ){
421  d = (pointsInPlanes[surf] - lowerPlaneBound).Dot(normalVects[surf])/traj.Dot(normalVects[surf]);
422  intersect = d*traj + lowerPlaneBound;
423 
424  // LOG_VERBATIM("CosmicTrackUtilities") << "intersection: "
425  // << intersect.X() << " "
426  // << intersect.Y() << " "
427  // << intersect.Z() << " surf "
428  // << surf << " "
429  // << normalVects[surf].X() << " "
430  // << normalVects[surf].Y() << " "
431  // << normalVects[surf].Z() << " "
432  // << d << " "
433  // << d*traj.X() << " "
434  // << d*traj.Y() << " "
435  // << d*traj.Z();
436 
437  // check that the intersection is between the bounds of the cell
438  // if not, remove the intersection, allow the boundaries to be fuzzy by
439  // 0.1 cm to account for rounding errors
440  if(intersect.X() > cellBounds[0] && intersect.X() < cellBounds[1] &&
441  intersect.Y() > cellBounds[2] && intersect.Y() < cellBounds[3] &&
442  intersect.Z() > cellBounds[4] && intersect.Z() < cellBounds[5])
443  intersections.push_back(intersect);
444 
445  }// end if the denominator is non-zero
446  }// end loop over surfaces
447 
448  // for(auto i : intersections) LOG_VERBATIM("CosmicTrackUtilitites") << "intersection point: "
449  // << i.X() << " " << i.Y() << " " << i.Z();
450 
451  // better have just 2 intersections at this point
452  // possible to have only one, if the track appears to stop in the cell, then use the reco hit location as
453  // one end point
454  if(intersections.size() > 2){
455  throw cet::exception("CosmicTrackUtilities") << "too many cell surfaces intersected "
456  << "by trajectory: " << intersections.size();
457  }
458  else if(intersections.size() == 1)
459  pathLength = (intersections[0] - recoHitLoc).Mag();
460  else if(intersections.size() < 1)
461  pathLength = 0.;
462  else
463  pathLength = (intersections[0] - intersections[1]).Mag();
464 
465  return pathLength;
466  }
art::ServiceHandle< geo::Geometry > fGeom
double HalfL() const
Definition: CellGeo.cxx:198
::xsd::cxx::tree::bounds< char > bounds
Definition: Database.h:226
double HalfD() const
Definition: CellGeo.cxx:205
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
float Dot(const Proxy &v) const
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
double HalfW() const
Definition: CellGeo.cxx:191
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
const PlaneGeo * Plane(unsigned int i) const
float abs(float number)
Definition: d0nt_math.hpp:39
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
Float_t d
Definition: plot.C:236
zBounds ZBounds(zBoundMap const &bounds, double const &hitZ)
std::pair< TVector3, TVector3 > zBounds
TVector3 Unit() const
float Mag() const
surf
Definition: demo5.py:35
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().

120  {
121  bool inDetector = true;
122 
123  // make sure the trajectory point is inside the detector
125  if( std::min(fLiveGeom->DistToEdgeXY(point), fLiveGeom->DistToEdgeZ(point)) < 0.)
126  inDetector = false;
127  }
128  else{
129  if(fLiveGeom->InMuonCatcher(point)){
130  int wall = 0;
131  if( fLiveGeom->DistanceToEdgeInMC(point, wall) < 0 )
132  inDetector = false;
133  }
134  else if( std::min(fLiveGeom->DistToEdgeXY(point), fLiveGeom->DistToEdgeZ(point)) < 0.)
135  inDetector = false;
136  }
137 
138  return inDetector;
139  }
art::ServiceHandle< geo::Geometry > fGeom
art::ServiceHandle< geo::LiveGeometry > fLiveGeom
double DistToEdgeXY(TVector3 vertex)
double DistanceToEdgeInMC(TVector3 vertex, int &wall)
double DistToEdgeZ(TVector3 vertex)
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
Near Detector in the NuMI cavern.
T min(const caf::Proxy< T > &a, T b)
bool InMuonCatcher(TVector3 vertex)
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().

60  {
61  // get the momentum from range. The value of dE/dx is given
62  // in MeV/cm, without a density factor. The value was determined
63  // from data/MC comparisons such that at the longest track lengths,
64  // the fractional difference in data and MC was centered at zero
65 
66  float dist = 0.;
67  float rangeMom = 0.;
68 
69  std::vector<TVector3> const& traj = track.Trajectory();
70  for(size_t t = traj.size()-1; t > 0; --t){
71  dist = (traj[t] - traj[t-1]).Mag();
72  rangeMom += dist * fEffectivedEdx * 1.e-3; // 1.e-3 converts MeV to GeV
73  }
74 
75  // The exponential factor accounts for higher dE/dx for shorter, ie
76  // less energetic, muons. The values were determined fitting a
77  // histogram of Delta p/p vs track length
78  return rangeMom*this->LengthCorrection(track.TotalLength());
79 
80  }
float LengthCorrection(float const &length)
Definition: event.h:19
double dist
Definition: runWimpSim.h:113
float Mag() const
double fEffectivedEdx
effective dE/dx in MeV/cm
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().

39  {
40 
41  auto pol = { 0.584713,
42  -0.00369301,
43  1.17201e-05,
44  -2.3046e-08,
45  2.9472e-11,
46  -2.47581e-14,
47  1.35044e-17,
48  -4.59248e-21,
49  8.82865e-25,
50  -7.31569e-29 };
51 
52  fEffectivedEdx = pset.get< double >("EffectivedEdx", 2.122);
53  fPolCoeff = pset.get< std::vector<double> >("PolCoeff", pol );
54 
55  return;
56  }
std::vector< double > fPolCoeff
the coefficients in an 8th degree polynomial
double fEffectivedEdx
effective dE/dx in MeV/cm
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().

586  {
587  LOG_DEBUG("CosmicTrackUtilities") << "TrackResiduals";
588 
589  std::vector<float> residuals;
590  residuals.reserve(track.NCell());
591 
592  TVector3 cross;
593  TVector3 point;
594  TVector3 lowerBound;
595  TVector3 upperBound;
596  double xyz[3] = {0.};
597 
598  zBoundMap zToPlaneBounds = this->MakeZBoundaryMap(track.Trajectory());
599 
600  // Loop through both views and grab every CellHit in this track
601  for(size_t c = 0; c < track.NCell(); ++c){
602  const art::Ptr<rb::CellHit> chit = track.Cell(c);
603 
604  // LOG_VERBATIM("CosmicTrackAna") << chit->Plane() << " " << chit->Cell() << "cell hit view is " << view;
605 
606  const rb::RecoHit rhit = track.RecoHit(chit);
607  if(!rhit.IsCalibrated()){
608  //mf::LogError("CosmicTrackAna") << "Not calibrated?!";
609  continue;
610  }
611 
612  fGeom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
613  // find the trajectory points of the track that bound this
614  // hit in z and determine where the hit is expected to be
615  // based on the direction of the track between the trajectory
616  // points
617  auto bounds = this->ZBounds(zToPlaneBounds, xyz[2]);
618 
619  lowerBound = bounds.first;
620  upperBound = bounds.second;
621 
622  point = TVector3(rhit.X(), rhit.Y(), rhit.Z());
623 
624  cross = (point - lowerBound).Cross(point - upperBound);
625  residuals.push_back(cross.Mag()/(upperBound - lowerBound).Mag());
626 
627  }// end loop over cells in the track
628 
629  return residuals;
630  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
art::ServiceHandle< geo::Geometry > fGeom
::xsd::cxx::tree::bounds< char > bounds
Definition: Database.h:226
std::map< double, zBounds > zBoundMap
unsigned short Plane() const
Definition: CellHit.h:39
float Z() const
Definition: RecoHit.h:38
Definition: event.h:19
const PlaneGeo * Plane(unsigned int i) const
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
unsigned short Cell() const
Definition: CellHit.h:40
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
zBounds ZBounds(zBoundMap const &bounds, double const &hitZ)
zBoundMap MakeZBoundaryMap(std::vector< TVector3 > const &trajectory)
float X() const
Definition: RecoHit.h:36
float Mag() const
float Y() const
Definition: RecoHit.h:37
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().

298  {
299  float pathLength = -9999.;
300  auto planeDirMap = track.PlaneDirMap();
301  auto pdmitr = planeDirMap.find(plane);
302 
303  if( pdmitr != planeDirMap.end() ){
304  float dw = (view == geo::kX) ? std::abs(pdmitr->second.X()) : std::abs(pdmitr->second.Y());
305  pathLength = cellWidth/dw;
306  }
307 
308  return pathLength;
309  }
Vertical planes which measure X.
Definition: PlaneGeo.h:28
Definition: event.h:19
float abs(float number)
Definition: d0nt_math.hpp:39
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().

143  {
144  TVector3 start, end, prev;
145  bool setStart = false;
146  float length = 0.;
147  float momentum = 0.;
148 
149  // loop over the particle trajectory points and set the length and end points based on
150  // those points that are in the detector
151  for(auto const& traj : part.Trajectory() ){
152 
153  if( this->PointInDetector(traj.first.Vect()) ){
154  if(!setStart){
155  setStart = true;
156  start.SetXYZ(traj.first.X(), traj.first.Y(), traj.first.Z());
157  length += fLiveGeom->ProjectedDistanceToEdge(traj.first.Vect(), -traj.second.Vect().Unit());
158  momentum = std::abs(traj.second.P());
159  prev = start;
160  }
161  end.SetXYZ(traj.first.X(), traj.first.Y(), traj.first.Z());
162  length += (end - prev).Mag();
163  prev = end;
164  }
165  }
166 
167  return std::pair<std::pair<float, float>,
168  std::pair<TVector3, TVector3> >(std::pair<float, float>(length, momentum),
169  std::pair<TVector3, TVector3>(start,end));
170  }
bool PointInDetector(TVector3 const &point)
art::ServiceHandle< geo::LiveGeometry > fLiveGeom
float abs(float number)
Definition: d0nt_math.hpp:39
TString part[npart]
Definition: Style.C:32
length
Definition: demo0.py:21
void prev()
Definition: show_event.C:91
float Mag() const
double ProjectedDistanceToEdge(TVector3 vertex, TVector3 dir)
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().

175  {
176  LOG_DEBUG("CosmicTrackUtilities") << "TruePathLengthInCell";
177  std::map<geo::CellUniqueId, float> pathInCell;
178  std::vector< TVector3 > trajPoints;
179  std::vector< TVector3 > flsTrajectory;
180  TVector3 point;
181  float length = 0.;
182 
183  // loop over the particle trajectory points and get the points that are in the detector
184  // for(auto const& traj : part.Trajectory() ){
185 
186  // if( !this->PointInDetector(traj.first.Vect()) ) continue;
187 
188  // trajPoints.push_back(TVector3(traj.first.X(), traj.first.Y(), traj.first.Z()));
189  // }
190 
191  // the true trajectory likely does not have a point for every plane crossed in the detector.
192  // get the trajectory by using the average FLS hit position in each plane
193  std::map<std::pair<size_t, size_t>, TVector3> flsHitCenters;
194  std::map<size_t, std::pair<TVector3, double> > flsPlaneZPos;
195  auto fpzpItr = flsPlaneZPos.begin();
196  double cnt = 0.;
197  double xyz[3] = {0.};
198  double loc[3] = {0.};
199  for(auto fls : flsHits){
200 
201  cnt = 1.;
202  loc[0] = fls.GetXAverage();
203  loc[1] = fls.GetYAverage();
204  loc[2] = fls.GetZAverage();
205  fGeom->Plane(fls.GetPlaneID())->Cell(fls.GetCellID())->LocalToWorld(loc, xyz);
206  point.SetXYZ(xyz[0], xyz[1], xyz[2]);
207  flsHitCenters[std::make_pair(fls.GetPlaneID(), fls.GetCellID())] = point;
208 
209  // now make the map to get the average position in a plane mapped to the
210  // center z poisition of the plane
211  fGeom->Plane(fls.GetPlaneID())->Cell(fls.GetCellID())->GetCenter(xyz);
212  if(flsPlaneZPos.count(fls.GetPlaneID()) > 0){
213  fpzpItr = flsPlaneZPos.find(fls.GetPlaneID());
214  point += fpzpItr->second.first;
215  cnt += fpzpItr->second.second;
216  }
217 
218  flsPlaneZPos[fls.GetPlaneID()] = std::make_pair(point, cnt);
219 
220  } // end loop over fls hits
221 
222  // get the average 3D position at each plane
223  for(auto pzp : flsPlaneZPos){
224  point = pzp.second.first;
225  cnt = pzp.second.second;
226  point *= 1./cnt;
227  flsTrajectory.push_back(point);
228  }
229 
230  auto zBounds = MakeZBoundaryMap(flsTrajectory);
231 
232  // loop over the FLS hits to get the planes and cells for this particle
233  for(auto fhc : flsHitCenters){
234 
235  try{
236  length = this->PathLengthInCell(zBounds, fhc.second, fhc.first);
237  }
238  catch(cet::exception &e){
239  LOG_WARNING("CosmicTrackAna")
240  << "caught exception\n"
241  << e
242  << "\n when trying to calculate tracklength in cell. Set it to a std::numeric_limits<float>::min()";
244  }
245 
246  pathInCell[fGeom->Plane(fhc.first.first)->Cell(fhc.first.second)->Id()] = length;
247 
248  } // end loop over fls hits
249 
250  return pathInCell;
251  }
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
art::ServiceHandle< geo::Geometry > fGeom
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
const PlaneGeo * Plane(unsigned int i) const
const CellUniqueId & Id() const
Definition: CellGeo.h:55
length
Definition: demo0.py:21
std::pair< TVector3, TVector3 > zBounds
#define LOG_WARNING(category)
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
zBoundMap MakeZBoundaryMap(std::vector< TVector3 > const &trajectory)
Float_t e
Definition: plot.C:35
float PathLengthInCell(zBoundMap const &zBounds, TVector3 const &recoHitLoc, std::pair< uint32_t, uint32_t > const &pc)
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().

563  {
564  // for(auto itr = bounds.begin(); itr != bounds.end(); ++itr)
565  // LOG_VERBATIM("CosmicTrackUtilities") << "hit z: " << hitZ << " " << itr->first << " lower "
566  // << itr->second.first.X() << " "
567  // << itr->second.first.Y() << " "
568  // << itr->second.first.Z() << " upper "
569  // << itr->second.second.X() << " "
570  // << itr->second.second.Y() << " "
571  // << itr->second.second.Z();
572 
573  zBoundMap::const_iterator lbitr = bounds.lower_bound(hitZ - fPlaneZDepth);
574 
575  // if we did not find an iterator, then the z position passed is
576  // off the edge of the track, use the last entry in the map
577  if(lbitr == bounds.end() )
578  return zBounds(bounds.rbegin()->second.first, bounds.rbegin()->second.second);
579 
580  return zBounds(lbitr->second.first, lbitr->second.second);
581 
582  }
double fPlaneZDepth
extent of one plane (cell) in the z direction
::xsd::cxx::tree::bounds< char > bounds
Definition: Database.h:226
std::pair< TVector3, TVector3 > zBounds

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
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: