CosmicTrackUtilities.cxx
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////
2 /// \brief Cosmic track selection code
3 /// \author Brian Rebel - brebel@fnal.gov
4 /// \date
5 //////////////////////////////////////////////////////////////////////////
6 #include <vector>
7 #include <limits>
8 #include <initializer_list>
9 
10 #include "TVector3.h"
11 
13 
15 #include "RecoBase/Track.h"
16 #include "RecoBase/RecoHit.h"
18 #include "Simulation/FLSHit.h"
20 
21 namespace trk{
22 
23  //----------------------------------------------------------------------------
25  {
27  this->reconfigure(pset);
28 
29  return;
30  }
31 
32  //----------------------------------------------------------------------------
34  {
35  }
36 
37  //----------------------------------------------------------------------------
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  }
57 
58  //----------------------------------------------------------------------------
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  }
81 
82  //----------------------------------------------------------------------------
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  }
99 
100  //----------------------------------------------------------------------------
101  float CosmicTrackUtilities::CosTheta(float const& dcosY,
102  float const& magnitude)
103  {
104  return -dcosY/magnitude;
105  }
106 
107  //----------------------------------------------------------------------------
108  float CosmicTrackUtilities::Azimuth(float const& dcosX,
109  float const& dcosZ,
110  float const& magnitude)
111  {
112  float azimuth = std::atan2(dcosX/magnitude, -dcosZ/magnitude) * 180./TMath::Pi();
113  if(azimuth < 0.) azimuth += 360.;
114 
115  return azimuth;
116  }
117 
118  //----------------------------------------------------------------------------
119  bool CosmicTrackUtilities::PointInDetector(TVector3 const& point)
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  }
140 
141  //----------------------------------------------------------------------------
142  std::pair<std::pair<float, float>, std::pair<TVector3, TVector3> > CosmicTrackUtilities::TrueLengthEndPoints(simb::MCParticle const& part)
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  }
171 
172  //----------------------------------------------------------------------------
173  std::map<geo::CellUniqueId, float> CosmicTrackUtilities::TruePathLengthInCell(simb::MCParticle const& part,
174  std::vector<sim::FLSHit> const& flsHits)
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  }
252 
253  //----------------------------------------------------------------------------
255  std::map<std::pair<uint32_t, uint32_t>,bool> & planeCellMap)
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  }
292 
293  //----------------------------------------------------------------------------
295  uint32_t const& plane,
296  double const& cellWidth,
297  geo::View_t const& view)
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  }
310 
311 
312  //----------------------------------------------------------------------------
313  // the cellExtent values correspond to:
314  // X: the half width of the cell in X
315  // Y: the half width of the cell in Y
316  // Z: the half depth of the cell (ie extent in the detector z direction)
318  TVector3 const& recoHitLoc,
319  std::pair<uint32_t, uint32_t> const& pc)
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  }
467 
468  //----------------------------------------------------------------------------
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  }
487 
488  //----------------------------------------------------------------------------
489  zBoundMap CosmicTrackUtilities::MakeZBoundaryMap(std::vector<TVector3> const& trajectory)
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  }
559 
560  //----------------------------------------------------------------------------
562  double const& hitZ)
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  }
583 
584  //----------------------------------------------------------------------------
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  }
631 } // namespace trk
632 
double fPlaneZDepth
extent of one plane (cell) in the z direction
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
art::ServiceHandle< geo::Geometry > fGeom
float LengthCorrection(float const &length)
double HalfL() const
Definition: CellGeo.cxx:198
float TriCellPathLength(rb::Track const &track, uint32_t const &plane, double const &cellWidth, geo::View_t const &view)
::xsd::cxx::tree::bounds< char > bounds
Definition: Database.h:226
double HalfD() const
Definition: CellGeo.cxx:205
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
float Dot(const Proxy &v) const
std::map< double, zBounds > zBoundMap
bool PointInDetector(TVector3 const &point)
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
art::ServiceHandle< geo::LiveGeometry > fLiveGeom
const simb::MCTrajectory & Trajectory() const
Definition: MCParticle.h:252
unsigned short Plane() const
Definition: CellHit.h:39
double DistToEdgeXY(TVector3 vertex)
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
constexpr T pow(T x)
Definition: pow.h:75
Vertical planes which measure X.
Definition: PlaneGeo.h:28
void FindTriCells(rb::Track const &track, std::map< std::pair< uint32_t, uint32_t >, bool > &planeCellMap)
float Z() const
Definition: RecoHit.h:38
Definition: event.h:19
double HalfW() const
Definition: CellGeo.cxx:191
double DistanceToEdgeInMC(TVector3 vertex, int &wall)
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
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
Particle class.
float abs(float number)
Definition: d0nt_math.hpp:39
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
std::pair< std::pair< float, float >, std::pair< TVector3, TVector3 > > TrueLengthEndPoints(simb::MCParticle const &part)
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
double dist
Definition: runWimpSim.h:113
const CellUniqueId & Id() const
Definition: CellGeo.h:55
Float_t Z
Definition: plot.C:38
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
unsigned short Cell() const
Definition: CellHit.h:40
Track finder for cosmic rays.
std::map< geo::CellUniqueId, float > TruePathLengthInCell(simb::MCParticle const &part, std::vector< sim::FLSHit > const &flsHits)
TString part[npart]
Definition: Style.C:32
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
std::vector< double > fPolCoeff
the coefficients in an 8th degree polynomial
length
Definition: demo0.py:21
void prev()
Definition: show_event.C:91
double DistToEdgeZ(TVector3 vertex)
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
void reconfigure(fhicl::ParameterSet const &pset)
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
std::vector< float > TrackResiduals(rb::Track const &track)
Near Detector in the NuMI cavern.
float Azimuth(float const &dcosX, float const &dcosZ, float const &magnitude)
Float_t d
Definition: plot.C:236
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)
float CosTheta(float const &dcosY, float const &magnitude)
std::pair< TVector3, TVector3 > zBounds
TVector3 Unit() const
#define LOG_WARNING(category)
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
std::map< unsigned int, TVector3 > PlaneDirMap() const
map of the direction cosines at each plane (ie z)
Definition: Track.cxx:496
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
TDirectory * dir
Definition: macro.C:5
zBoundMap MakeZBoundaryMap(std::vector< TVector3 > const &trajectory)
float X() const
Definition: RecoHit.h:36
std::vector< TVector3 > const & Trajectory() const
return a constant reference to the track trajectory points
Definition: Track.h:86
float Mag() const
float Y() const
Definition: RecoHit.h:37
unsigned int NPlanes() const
T min(const caf::Proxy< T > &a, T b)
bool InMuonCatcher(TVector3 vertex)
Float_t e
Definition: plot.C:35
Float_t X
Definition: plot.C:38
double fEffectivedEdx
effective dE/dx in MeV/cm
std::set< double > fPlaneZBounds
boundaries of each plane in the z direction
float PathLengthInCell(zBoundMap const &zBounds, TVector3 const &recoHitLoc, std::pair< uint32_t, uint32_t > const &pc)
double ProjectedDistanceToEdge(TVector3 vertex, TVector3 dir)
float RangeMomentum(rb::Track const &track)
T atan2(T number)
Definition: d0nt_math.hpp:72
surf
Definition: demo5.py:35