CalibUtil.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file CalibUtil.cxx
3 // \brief Functions of general use in calibration
4 // \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
7 #include <glob.h>
8 #include <cstring>
9 #include <sys/types.h>
10 #include <sys/stat.h>
11 
13 
14 // Framework includes
17 
18 // NOvASoft includes
21 #include "Geometry/Geometry.h"
22 #include "Geometry/LiveGeometry.h"
23 #include "GeometryObjects/Geo.h"
25 
26 namespace calib
27 {
28 
30  {
31  switch(det) {
32  case novadaq::cnv::kFARDET: return "fd";
33  case novadaq::cnv::kNEARDET: return "nd";
34  case novadaq::cnv::kTESTBEAM: return "tb";
35  default:
36  assert (!"The default case of getDetString switch was reached.");
37  return "NeverHere";
38  }
39  }
40 
41 
42 
43  /// Provide a simple utility function for our sets
44  struct ChannelSet: public std::set<geo::OfflineChan>
45  {
46  bool CellExists(int plane, int cell) const
47  {
48  return count(geo::OfflineChan(plane, cell));
49  }
50  };
51 
52  //......................................................................
53  /// Internal helper for \ref BestPathEstimates. Takes account of bad channels
54  bool HasXYAdjacents(int plane, int cell, const ChannelSet& hmap)
55  {
58 
59  // Look first below and then above this cell
60  for(int step = -1; step <= +1; step += 2){
61  int curCell = cell;
62  while(true){
63  curCell += step;
64 
65  // Ran out of detector
66  if(curCell < 0 ||
67  curCell >= int(geom->Plane(plane)->Ncells())) return false;
68 
69  // If this channel is bad we can keep on looking in this direction
70  if(bcl->IsBad(plane, curCell)) continue;
71 
72  // If there's a hit this side move on to the other side. If there's
73  // isn't then this cell fails the test.
74  if(hmap.CellExists(plane, curCell))
75  break;
76  else
77  return false;
78 
79  } // end while
80  } // end for step
81  return true;
82  }
83 
84  //......................................................................
85  // Internal helper for \ref BestPathEstimates. Needs to take account of bad
86  // channels.
87  bool HasZAdjacents(int plane, int cell,
88  const ChannelSet& hmap,
90  {
91  return (!hmap.CellExists(plane, cell-1) &&
92  !hmap.CellExists(plane, cell+1) &&
93  hmap.CellExists(geom->NextPlaneInView(plane, +1), cell) &&
94  hmap.CellExists(geom->NextPlaneInView(plane, -1), cell));
95  }
96 
97  //......................................................................
98  std::vector<double> BestPathEstimates(std::vector<geo::OfflineChan> const& chans,
99  std::map<unsigned int, TVector3> const& directions,
100  std::vector<EPathQuality> & quals,
101  std::vector<geo::OfflineChan> const& extras)
102  {
103  const size_t trkN = chans.size();
104  const size_t extrasN = extras.size();
105  const size_t hitMax = trkN+extrasN;
106 
107  std::vector<double> paths;
108  paths.resize(hitMax);
109 
110  quals.clear();
111  quals.resize(hitMax);
112 
114 
115  ChannelSet hmap;
116  for(size_t i = 0; i < trkN; ++i) hmap.insert(chans[i]);
117  hmap.insert(extras.begin(), extras.end());
118 
119 
120  // Look for ones that have perpendicular hits either side
121  for(size_t hitIdx = 0; hitIdx < hitMax; ++hitIdx){
122  int plane, cell;
124  if(hitIdx < trkN){
125  plane = chans[hitIdx].Plane();
126  cell = chans[hitIdx].Cell();
127  view = geom->Plane(plane)->View();
128  }
129  else{
130  plane = extras[hitIdx-trkN].Plane();
131  cell = extras[hitIdx-trkN].Cell();
132  view = geom->Plane(plane)->View();
133  }
134 
136 
137  // Mark as xy if it satisfies the criteria. Or if not then maybe as
138  // z. Otherwise falls back to the default kAverage.
139  if(HasXYAdjacents(plane, cell, hmap)){
140  qual = kXYAdjacents;
141  }
142  else{
143  if(HasZAdjacents(plane, cell, hmap, geom)){
144  qual = kZAdjacents;
145  }
146  }
147 
148  quals[hitIdx] = qual;
149 
150  // get the direction for this plane
151  TVector3 dr = calib::DirectionAtPlane(directions, plane);
152 
153  const double dx = std::abs(dr.X());
154  const double dy = std::abs(dr.Y());
155  const double dz = std::abs(dr.Z());
156  if(std::isnan(dx) || std::isnan(dy) || std::isnan(dz) ||
157  (view == geo::kX && dx == 0) ||
158  (view == geo::kY && dy == 0) ||
159  (qual == kZAdjacents && dz == 0)){
160  paths[hitIdx] = 0;
161  quals[hitIdx] = kBadTrack;
162  continue;
163  }
164 
165  const geo::CellGeo* cellgeo = geom->Plane(plane)->Cell(cell);
166  const double width = 2*cellgeo->HalfW();
167  const double depth = 2*cellgeo->HalfD();
168 
169  switch(qual){
170  case kXYAdjacents:
171  paths[hitIdx] = (view == geo::kX) ? width/dx : width/dy;
172  break;
173  case kZAdjacents:
174  paths[hitIdx] = depth/dz;
175  break;
176  case kAverage:
177  paths[hitIdx] = geo::AverageCellPathLength(view, dx, dy, dz);
178  break;
179  default:
180  assert(0 && "Shouldn't happen");
181  } // end switch
182  } // end for hitIdx
183 
184  return paths;
185  }
186 
187  //......................................................................
188  TVector3 const DirectionAtPlane(const std::map<unsigned int, TVector3>& directions,
189  unsigned int plane)
190  {
191  TVector3 dr;
192 
193  if(directions.count(plane) > 0){
194  dr = directions.find(plane)->second;
195  }
196  else{
197 
198  // the plane is not in the map, so now need to
199  // determine the direction of the track in z,
200  // and then use that result to determine whether
201  // to use the upper_bound or lower_bound methods
202  // remember that maps are ordered by the < operator
203  // of the key, in this case an unsigned int
204  // if the desired plane is outside the bounds of
205  // the map, just take the end point vector of the
206  // map in the direction of travel.
207 
208  auto ubitr = directions.upper_bound(plane);
209  auto lbitr = directions.lower_bound(plane);
210 
211  // the direction in z better not change or
212  // this track has problems
213  float dZ = directions.begin()->second.Z();
214  if( dZ > 0.){
215  if( ubitr != directions.end() )
216  dr = ubitr->second;
217  else
218  dr = directions.rbegin()->second;
219  }
220  if( dZ < 0. ){
221  if( lbitr != directions.end() )
222  dr = lbitr->second;
223  else
224  dr = directions.begin()->second;
225  }
226  } // end no key corresponding to plane in the map
227 
228  return dr;
229  }
230 
231  //......................................................................
232  void FindBelowThresholdCalibCandidates(const std::vector<geo::OfflineChan>& trkchans,
233  const std::map<unsigned int, TVector3>& directions,
234  std::vector<EPathQuality>& quals,
235  std::vector<double>& paths,
236  std::vector<geo::OfflineChan>& chans)
237  {
238  // In case user passed in junk
239  quals.clear();
240  paths.clear();
241  chans.clear();
242 
243  const int hitMax = trkchans.size();
244 
246 
247  ChannelSet hmap;
248  for(int n = 0; n < hitMax; ++n) hmap.insert(trkchans[n]);
249 
250  for(int hitIdx = 0; hitIdx < hitMax; ++hitIdx){
251  const int plane = trkchans[hitIdx].Plane();
252  const geo::View_t view = geom->Plane(plane)->View();
253 
254  auto const dr = calib::DirectionAtPlane(directions, plane);
255 
256  const double dx = std::abs(dr.X());
257  const double dy = std::abs(dr.Y());
258  const double dz = std::abs(dr.Z());
259 
260  if(std::isnan(dx) || std::isnan(dy) || std::isnan(dz)) continue;
261 
262  if((view == geo::kX && dx == 0) || (view == geo::kY && dy == 0)){
263  continue;
264  }
265 
266  // Need to search between each cell and the next hit cell in the same
267  // plane
268  int cell = trkchans[hitIdx].Cell();
269  while(true){
270  ++cell;
271 
272  // If it's hit, then it isn't thresholded. And that ends our search
273  // because we'll pick it up and start there next.
274  if(hmap.CellExists(plane, cell)) break;
275 
276  if(HasXYAdjacents(plane, cell, hmap)){
277  const geo::CellGeo* cellgeo = geom->Plane(plane)->Cell(cell);
278  const double width = 2*cellgeo->HalfW();
279 
280  quals.push_back(kXYAdjacents);
281  paths.push_back((view == geo::kX) ? width/dx : width/dy);
282  chans.push_back(geo::OfflineChan(plane, cell));
283  }
284  else{
285  // Definition of the XY criteria mean that nothing else in this range
286  // can work either
287  break;
288  }
289  } // end while
290  }
291  }
292 
293  //-----
294  void GetXYZD(geo::OfflineChan chan, double w, double *xyzd) {
295  // Return (via 'xyzd' argument) position of deposition in world
296  // coordinates and the distance to the readout, assuming deposition
297  // has cell-local-z position = w.
298 
300 
301  UInt_t plane = chan.Plane();
302  UInt_t cell = chan.Cell();
303  if ( plane >= geom->NPlanes() ) {
304  throw cet::exception("BAD_PLANE_NUMBER")
305  << "plane " << plane << " of " << geom->NPlanes() << " planes\n"
306  << __FILE__ << ":" << __LINE__ << "\n";
307  }
308  const geo::PlaneGeo* p = geom->Plane(plane);
309  if ( cell >= p->Ncells() ) {
310  throw cet::exception("BAD_CELL_NUMBER")
311  << "cell " << cell << " of " << p->Ncells()
312  << " in plane " << plane << "\n"
313  << __FILE__ << ":" << __LINE__ << "\n";
314  }
315  const geo::CellGeo* c = p->Cell(cell);
316 
317  // Should be replaced with alignment tables or whatnot. Okay for now.
318 
319  // w is not localZ, which makes this a lot messier than just
320  // c.GetCenter(xyzd,w) and c.DistToReadOut(w)
321  c->GetCenter(xyzd);
322  if (geom->Plane(plane)->View()==geo::kX)
323  xyzd[1] = w;
324  else
325  xyzd[0] = w;
326  Double_t localxyz[3];
327  c->WorldToLocal(xyzd,localxyz);
328  xyzd[3] = geom->DistToElectronics(localxyz[2], *c);
329  }
330 
331  //----------------------------------------------------------------------------
332  void FindZBoundaries(std::set<double> & planeZBounds)
333  {
334  // loop over all the planes in the geometry and put their upper z boundary
335  // into a set
336  std::set<double> bounds;
337  bounds.insert(0.);
338  double xyz[3] = {0.};
339 
341  double planeZDepth = 2.*geom->Plane(0)->Cell(0)->HalfD();
342 
343  for(size_t p = 0; p < geom->NPlanes(); ++p){
344  geom->Plane(p)->Cell(0)->GetCenter(xyz);
345  bounds.insert(xyz[2] + 0.5*planeZDepth);
346  }
347 
348  planeZBounds.swap(bounds);
349 
350  return;
351  }
352 
353  //----------------------------------------------------------------------------
354  zBoundMap MakeZBoundaryMap(std::set<double> const& planeZBounds,
355  std::vector<TVector3> const& trajectory)
356  {
357  zBoundMap zToPlaneBounds;
358 
359  std::vector<TVector3> localTraj(trajectory);
360 
361  // drop duplicate entries that are right next to each other to avoid zero
362  // vectors
363  auto itr1 = localTraj.begin();
364  auto itr2 = localTraj.begin();
365  while(itr1 != localTraj.end()){
366  // set the second iterator to be the next entry in the vector
367  itr2 = itr1;
368  ++itr2;
369  if(itr2 == localTraj.end()) break;
370 
371  // if the TVector3s are the same, erase itr2 and set
372  // itr1 to be the return value, ie the position after
373  // the original itr2. Otherwise increase itr1
374  if(std::abs((*itr1).X() - (*itr2).X()) < 0.1 &&
375  std::abs((*itr1).Y() - (*itr2).Y()) < 0.1 &&
376  std::abs((*itr1).Z() - (*itr2).Z()) < 0.1){
377  LOG_DEBUG("CalibUtil") << "erasing " << (*itr1).X() << " " << (*itr1).Y() << " " << (*itr1).Z();
378  itr1 = localTraj.erase(itr2);
379  }
380  else ++itr1;
381  }
382 
383  double deltaZ = 0.;
384  double lowZ = 0.;
385  double highZ = 0.;
386 
387  auto litr = planeZBounds.begin();
388  auto hitr = planeZBounds.begin();
389 
390  TVector3 dir;
391  TVector3 lowBound;
392  TVector3 highBound;
393 
396  double planeZDepth = 2.*geom->Plane(0)->Cell(0)->HalfD();
397 
398  for(size_t tp = 0; tp < localTraj.size(); ++tp){
399 
400  LOG_DEBUG("CalibUtils")
401  << tp << "/" << localTraj.size()
402  << " localTraj: " << localTraj[tp].X()
403  << " " << localTraj[tp].Y()
404  << " " << localTraj[tp].Z();
405 
406  // make sure the localTraj point is inside the detector
407  if( !liveGeo->IsPointLive(localTraj[tp]) ){
408  LOG_DEBUG("CalibUtils") << tp << " is not valid";
409  continue;
410  }
411 
412  // find the intersection of the localTrajectory with the boundary
413  // of the plane it represents
414  if(tp < localTraj.size() - 1) dir = localTraj[tp+1] - localTraj[tp];
415  else dir = localTraj[tp] - localTraj[tp-1];
416 
417  // find the plane boundaries for this point
418  litr = planeZBounds.lower_bound(localTraj[tp].Z() - planeZDepth);
419  hitr = planeZBounds.upper_bound(localTraj[tp].Z());
420 
421  // in the muon catcher, the high side and low side iterators could be the same
422  if(hitr == litr){
423  if(litr != planeZBounds.begin() ) --litr;
424  else ++hitr;
425  }
426  lowZ = *litr;
427  highZ = *hitr;
428 
429  LOG_DEBUG("CalibUtils")
430  << tp << " " << localTraj[tp].Z()
431  << " low z: " << lowZ << " high z: " << highZ
432  << " localTraj point: ("
433  << localTraj[tp].X() << " " << localTraj[tp].Y() << " " << localTraj[tp].Z() << ")";
434 
435  // make sure we don't go past the detector boundaries
436  deltaZ = localTraj[tp].Z() - lowZ;
437  lowBound.SetXYZ(localTraj[tp].X() - deltaZ*dir.X()/dir.Z(),
438  localTraj[tp].Y() - deltaZ*dir.Y()/dir.Z(),
439  lowZ);
440 
441  LOG_DEBUG("CalibUtils")
442  << "is bounded by " << deltaZ << " " << dir.X()/dir.Z() << " " << dir.Y()/dir.Z()
443  << " (" << lowBound.X() << " " << lowBound.Y() << " " << lowBound.Z() << ")";
444 
445  deltaZ = highZ - localTraj[tp].Z();
446  highBound.SetXYZ(localTraj[tp].X() + deltaZ*dir.X()/dir.Z(),
447  localTraj[tp].Y() + deltaZ*dir.Y()/dir.Z(),
448  highZ);
449 
450  LOG_DEBUG("CalibUtils")
451  << "and " << deltaZ << " ("
452  << highBound.X() << " " << highBound.Y() << " " << highBound.Z() << ")";
453 
454  zToPlaneBounds[localTraj[tp].Z()] = std::make_pair(lowBound, highBound);
455 
456  }
457 
458  return zToPlaneBounds;
459  }
460 
461  //----------------------------------------------------------------------------
463  double const& hitZ)
464  {
465  if(bounds.empty())
466  throw cet::exception("CalibUtil")
467  << "bounds map has no entries, cannot "
468  << "find bounding points for hit z location: "
469  << hitZ;
470 
471  // for(auto itr = bounds.begin(); itr != bounds.end(); ++itr)
472  // LOG_VERBATIM("CalibUtils") << "hit z: " << hitZ << " " << itr->first << " lower "
473  // << itr->second.first.X() << " "
474  // << itr->second.first.Y() << " "
475  // << itr->second.first.Z() << " upper "
476  // << itr->second.second.X() << " "
477  // << itr->second.second.Y() << " "
478  // << itr->second.second.Z();
480  double planeZDepth = 2.*geom->Plane(0)->Cell(0)->HalfD();
481 
482  zBoundMap::const_iterator lbitr = bounds.lower_bound(hitZ - planeZDepth);
483 
484  // if we did not find an iterator, then the z position passed is
485  // off the edge of the track, use the last entry in the map
486  if(lbitr == bounds.end() )
487  return zBounds(bounds.rbegin()->second.first, bounds.rbegin()->second.second);
488 
489  return zBounds(lbitr->second.first, lbitr->second.second);
490 
491  }
492 
493  //----------------------------------------------------------------------------
494  // the cellExtent values correspond to:
495  // X: the half width of the cell in X
496  // Y: the half width of the cell in Y
497  // Z: the half depth of the cell (ie extent in the detector z direction)
499  TVector3 const& recoHitLoc,
500  std::pair<uint32_t, uint32_t> const& pc)
501  {
503 
504  double xyz[3] = {0.};
505  geom->Plane(pc.first)->Cell(pc.second)->GetCenter(xyz);
506  TVector3 cellCenter = TVector3(xyz);
507  TVector3 cellExtent = TVector3(geom->Plane(pc.first)->Cell(pc.second)->HalfW(),
508  geom->Plane(pc.first)->Cell(pc.second)->HalfL(),
509  geom->Plane(pc.first)->Cell(pc.second)->HalfD());
510  if(geom->Plane(pc.first)->View() == geo::kY)
511  cellExtent = TVector3(geom->Plane(pc.first)->Cell(pc.second)->HalfL(),
512  geom->Plane(pc.first)->Cell(pc.second)->HalfW(),
513  geom->Plane(pc.first)->Cell(pc.second)->HalfD());
514 
515  // find the trajectory points of the track that bound this
516  // hit in z and determine where the hit is expected to be
517  // based on the direction of the track between the trajectory
518  // points
519  calib::zBounds boundaries = ZBounds(zBounds, xyz[2]);
520 
521 
522  TVector3 lowerPlaneBound = boundaries.first;
523  TVector3 upperPlaneBound = boundaries.second;
524 
525  // now determine where the trajectory intersects the cell boundaries
526  // evaluate the point where the line intersects the plane of each
527  // boundary and see if that point is within the cell
528  // see http://en.wikipedia.org/wiki/Line–plane_intersection
529 
530  float pathLength = 0.;
531  float d = 0.;
532  // all vectors are in order of -x, +x, -y, +y, -z, +z
533 
534  // determine a point on each plane bounding the cell
535  // define the edges of the cell using its center and extent values
536  auto pointsList = {TVector3(cellCenter.X() - cellExtent.X(), cellCenter.Y(), cellCenter.Z()),
537  TVector3(cellCenter.X() + cellExtent.X(), cellCenter.Y(), cellCenter.Z()),
538  TVector3(cellCenter.X(), cellCenter.Y() - cellExtent.Y(), cellCenter.Z()),
539  TVector3(cellCenter.X(), cellCenter.Y() + cellExtent.Y(), cellCenter.Z()),
540  TVector3(cellCenter.X(), cellCenter.Y(), cellCenter.Z() - cellExtent.Z()),
541  TVector3(cellCenter.X(), cellCenter.Y(), cellCenter.Z() + cellExtent.Z())};
542  std::vector<TVector3> pointsInPlanes(pointsList);
543 
544  auto normalList = {TVector3(-1., 0., 0.),
545  TVector3( 1., 0., 0.),
546  TVector3( 0., -1., 0.),
547  TVector3( 0., 1., 0.),
548  TVector3( 0., 0., -1.),
549  TVector3( 0., 0., 1.) };
550  std::vector<TVector3> normalVects(normalList);
551 
552  std::vector<TVector3> intersections;
553  TVector3 intersect;
554 
555  // vector describing the boundaries of this cell
556  auto boundList = {pointsInPlanes[0].X() - 1.e-4, pointsInPlanes[1].X() + 1.e-4,
557  pointsInPlanes[2].Y() - 1.e-4, pointsInPlanes[3].Y() + 1.e-4,
558  pointsInPlanes[4].Z() - 1.e-4, pointsInPlanes[5].Z() + 1.e-4};
559  std::vector<double> cellBounds(boundList);
560 
561  // check that the reco hit is between the bounds, return 0 if not
562  if(recoHitLoc.X() < cellBounds[0] || recoHitLoc.X() > cellBounds[1] ||
563  recoHitLoc.Y() < cellBounds[2] || recoHitLoc.Y() > cellBounds[3] ||
564  recoHitLoc.Z() < cellBounds[4] || recoHitLoc.Z() > cellBounds[5])
565  return 0.;
566 
567  // the recoHitLoc is the determined 3D position for the given hit
568  // the plane bounds are the 3D positions where the track crosses the
569  // plane containing this hit in the upstream and downstream locations
570 
571  // the vector from the upstream and downstream locations gives the direction
572  // from one side of the plane to the other. Grab the location of the hit
573  // and then walk along the direction vector until we step out of the cell
574  // be sure to check the boundary of the cell in both the z direction and
575  // in the width of the cell
576  TVector3 traj = (upperPlaneBound - lowerPlaneBound).Unit();
577 
578  // loop over the different planes bounding the cell and find the intersection points
579  // go in the reverse order and remove any surfaces where the intersection is outside
580  // of the bounds of the cell
581  for(size_t surf = 0; surf < normalVects.size(); ++surf){
582 
583  bool repeatIntersection = false;
584 
585  // check that you can do the calculation, if not, there is no
586  // intersectoin with that surface and we remove it.
587  if( std::abs(traj.Dot(normalVects[surf])) > 0. ){
588  d = (pointsInPlanes[surf] - lowerPlaneBound).Dot(normalVects[surf])/traj.Dot(normalVects[surf]);
589  intersect = d*traj + lowerPlaneBound;
590 
591  // LOG_VERBATIM("CalibUtils") << "intersection: "
592  // << intersect.X() << " "
593  // << intersect.Y() << " "
594  // << intersect.Z() << " surf "
595  // << surf << " "
596  // << normalVects[surf].X() << " "
597  // << normalVects[surf].Y() << " "
598  // << normalVects[surf].Z() << " "
599  // << d << " "
600  // << d*traj.X() << " "
601  // << d*traj.Y() << " "
602  // << d*traj.Z();
603 
604  // check that the intersection is between the bounds of the cell
605  // if not, remove the intersection, allow the boundaries to be fuzzy by
606  // 0.1 cm to account for rounding errors
607  if(intersect.X() > cellBounds[0] && intersect.X() < cellBounds[1] &&
608  intersect.Y() > cellBounds[2] && intersect.Y() < cellBounds[3] &&
609  intersect.Z() > cellBounds[4] && intersect.Z() < cellBounds[5]){
610 
611  // check that the intersection is not already in the vector
612  // that can happen if the intersection is in a corner of the cell, then
613  // two surfaces have the same intersection
614  for(auto p : intersections){
615  if(std::abs(p.X() - intersect.X()) < 0.01 &&
616  std::abs(p.Y() - intersect.Y()) < 0.01 &&
617  std::abs(p.Z() - intersect.Z()) < 0.01 ){
618  // LOG_VERBATIM("CalibUtils") << "intersection point repeats " << intersect.X()
619  // << " " << intersect.Y() << " " << intersect.Z();
620  repeatIntersection = true;
621  break;
622  }
623  } // end loop to look for a repeat intersection
624  if(!repeatIntersection) intersections.push_back(intersect);
625 
626  } // end if there is a legit intersection
627  }// end if the denominator is non-zero
628  }// end loop over surfaces
629 
630  // better have just 2 intersections at this point
631  // possible to have only one, if the track appears to stop in the cell, then use the reco hit location as
632  // one end point
633  if(intersections.size() > 2){
634  LOG_VERBATIM("CalibUtils")
635  << "reco hit location: "
636  << recoHitLoc.X() << " " << recoHitLoc.Y() << " " << recoHitLoc.Z()
637  <<"\ncell center: "
638  << cellCenter.X() << " " << cellCenter.Y() << " " << cellCenter.Z()
639  << "\ncell boundaries: "
640  << cellBounds[0] << " " << cellBounds[1] << " "
641  << cellBounds[2] << " " << cellBounds[3] << " "
642  << cellBounds[4] << " " << cellBounds[5];
643 
644 
645  LOG_VERBATIM("CalibUtils")
646  << "upper bound: "
647  << upperPlaneBound.X() << " "
648  << upperPlaneBound.Y() << " "
649  << upperPlaneBound.Z() << "\nlower bound: "
650  << lowerPlaneBound.X() << " "
651  << lowerPlaneBound.Y() << " "
652  << lowerPlaneBound.Z() << "\ntrajectory: "
653  << traj.X() << " " << traj.Y() << " " << traj.Z();
654 
655  for(auto i : intersections)
656  LOG_VERBATIM("CosmicTrackUtilitites")
657  << "intersection point: "
658  << i.X() << " " << i.Y() << " " << i.Z();
659 
660  throw cet::exception("CalibUtils")
661  << "too many cell surfaces intersected "
662  << "by trajectory: " << intersections.size();
663  }
664  else if(intersections.size() == 1)
665  pathLength = (intersections[0] - recoHitLoc).Mag();
666  else if(intersections.size() < 1)
667  pathLength = 0.;
668  else
669  pathLength = (intersections[0] - intersections[1]).Mag();
670 
671  return pathLength;
672  }
673 
674 
675 
676  std::vector<std::string> globWrapper(const std::string& pattern) {
677  using namespace std;
678 
679  // glob struct resides on the stack
680  glob_t glob_result;
681  memset(&glob_result, 0, sizeof(glob_result));
682 
683  // do the glob operation
684  int return_value = glob(pattern.c_str(), GLOB_TILDE, NULL, &glob_result);
685  if(return_value != 0) {
686  globfree(&glob_result);
687  stringstream ss;
688  ss << "glob(" << pattern.c_str() << ") failed with return_value " << return_value << endl;
689  throw std::runtime_error(ss.str());
690  }
691 
692  // collect all the filenames into a std::list<std::string>
693  vector<string> filenames;
694  for(size_t i = 0; i < glob_result.gl_pathc; ++i) {
695  filenames.push_back(string(glob_result.gl_pathv[i]));
696  }
697 
698  // cleanup
699  globfree(&glob_result);
700 
701  // done
702  return filenames;
703 }
704 
705 
706 std::string getCSVFilenameByParsingDirectory(int fCurrentRun, std::string prePattern, std::string postPattern) {
707  //This function just finds all the CSV files that match the pattern
708  // and thens earches through these files to find which one covers the fCurrentRun
709  std::string pattern=prePattern+std::string("*")+postPattern;
710  std::vector<std::string> answers=globWrapper(pattern);
711  int startRun=0; int endRun=0;
712 
713  //Check if there were any matches
714  if(answers.size()==0) {
715  std::cout << "We didn't find anything that matches the pattern for calibration CSV files\n";
716  std::cout << pattern << "\n";
717  assert(0 && "Shouldn't happen");
718  }
719 
720  //Loop through the matches
721  for(auto path : answers) {
722  if(path.size() > prePattern.size() + postPattern.size()) {
723  std::string epStr=path.substr(prePattern.size(),path.size()-(prePattern.size()+postPattern.size()));
724  // std::cout << epStr.c_str() << "\n";
725  if(epStr.c_str()[2]=='-') {
726  //Start run = 0
727  startRun=0;
728  }
729  else {
730  size_t hyphenPlace = epStr.find("-");
731  startRun=stoi( epStr.substr(2,hyphenPlace-2));
732  }
733  if(epStr.c_str()[epStr.size()-2]=='-') {
734  endRun=0;
735  }
736  else {
737  size_t hyphenPlace = epStr.find("-");
738  // std::cout << "Substring: " << epStr.substr(hyphenPlace+2) << "\n";
739  endRun=stoi(epStr.substr(hyphenPlace+2 ));
740  }
741 
742  }
743  // std::cout << startRun << "\t" << endRun << "\n";
744  if(fCurrentRun>=startRun && (endRun==0 || fCurrentRun<=endRun)) {
745  // std::cout << fCurrentRun << "\t" << path << "\n";
746  return path;
747  }
748 
749  }
750 
751  std::cout << "We didn't find anything that matches the pattern for calibration CSV files\n";
752  std::cout << "This shouldn't happen and might mean that the calibration CSV files are incomplete.\n";
753  std::cout << "Looking for run: " << fCurrentRun << " using pattern:\n";
754  std::cout << pattern << "\n";
755  assert(0 && "Shouldn't happen");
756  return std::string("");
757 }
758 
760  {
761  struct stat info;
762 
763  if(stat( dirSting.c_str(), &info ) != 0)
764  return false;
765  else if(info.st_mode & S_IFDIR)
766  return true;
767  else
768  return false;
769  }
770 
771 
772 } // end namespace
773 ////////////////////////////////////////////////////////////////////////
774 
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
const XML_Char XML_Encoding * info
Definition: expat.h:530
zBounds ZBounds(zBoundMap const &bounds, double const &hitZ)
Return the bounding points for a given z position along a track.
Definition: CalibUtil.cxx:462
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
bool doesDirectoryExist(std::string dirSting)
Return true if the path in dirString points to a directory.
Definition: CalibUtil.cxx:759
bool CellExists(int plane, int cell) const
Definition: CalibUtil.cxx:46
void WorldToLocal(const double *local, double *world) const
Definition: CellGeo.cxx:100
float Dot(const Proxy &v) const
zBoundMap MakeZBoundaryMap(std::set< double > const &planeZBounds, std::vector< TVector3 > const &trajectory)
Return a map of the z position of each trajectory point on a track to the bounding positions of where...
Definition: CalibUtil.cxx:354
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
Vertical planes which measure X.
Definition: PlaneGeo.h:28
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
std::vector< double > BestPathEstimates(std::vector< geo::OfflineChan > const &chans, std::map< unsigned int, TVector3 > const &directions, std::vector< EPathQuality > &quals, std::vector< geo::OfflineChan > const &extras)
Definition: CalibUtil.cxx:98
Float_t ss
Definition: plot.C:24
double HalfW() const
Definition: CellGeo.cxx:191
bool HasXYAdjacents(int plane, int cell, const ChannelSet &hmap)
Internal helper for BestPathEstimates. Takes account of bad channels.
Definition: CalibUtil.cxx:54
::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
std::pair< TVector3, TVector3 > zBounds
Definition: CalibUtil.h:106
float abs(float number)
Definition: d0nt_math.hpp:39
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
Float_t Z
Definition: plot.C:38
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
Far Detector at Ash River, MN.
CDPStorage service.
std::string getCSVFilenameByParsingDirectory(int fCurrentRun, std::string prePattern, std::string postPattern)
Definition: CalibUtil.cxx:706
double dy[NP][NC]
double dx[NP][NC]
double dz[NP][NC]
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
Near Detector in the NuMI cavern.
void GetXYZD(geo::OfflineChan chan, double w, double *xyzd)
Return position in world coordninates and distance to the readout.
Definition: CalibUtil.cxx:294
Collect Geo headers and supply basic geometry functions.
Float_t d
Definition: plot.C:236
EPathQuality
Methods used to estimate path length of a track through a cell.
Definition: CalibUtil.h:40
std::vector< std::string > globWrapper(const std::string &pattern)
Definition: CalibUtil.cxx:676
double AverageCellPathLength(geo::View_t view, double dx, double dy, double dz)
Mean path length of a ray with (unit) direction vector dx, dy, dz through a cell in view...
Definition: Geo.cxx:804
TVector3 Unit() const
const Cut qual
OStream cout
Definition: OStream.cxx:6
unsigned short Plane() const
Definition: OfflineChan.h:31
void FindBelowThresholdCalibCandidates(const std::vector< geo::OfflineChan > &trkchans, const std::map< unsigned int, TVector3 > &directions, std::vector< EPathQuality > &quals, std::vector< double > &paths, std::vector< geo::OfflineChan > &chans)
Find empty cells that would have contributed to the calibration if they&#39;d been hit. We assume they fell below threshold.
Definition: CalibUtil.cxx:232
const std::string path
Definition: plot_BEN.C:43
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
unsigned short Cell() const
Definition: OfflineChan.h:32
double DistToElectronics(double localz, const CellGeo &cell) const
get distance from local z position in cell to apd in cm, including pigtail
TDirectory * dir
Definition: macro.C:5
TVector3 const DirectionAtPlane(const std::map< unsigned int, TVector3 > &directions, unsigned int plane)
Return the direction at a plane using the map provided by rb::Track.
Definition: CalibUtil.cxx:188
const char * getDetString(novadaq::cnv::DetId det)
Utility function to get detector name as a string for file manipulations.
Definition: CalibUtil.cxx:29
const unsigned int NextPlaneInView(unsigned int p, int d=+1) const
A (plane, cell) pair.
Definition: OfflineChan.h:17
bool IsPointLive(TVector3 vertex)
Note the muon catcher is considered bad; use in combination with InMuonCatcher if needed...
void geom(int which=0)
Definition: geom.C:163
float Mag() const
assert(nhit_max >=nhit_nbins)
std::map< double, zBounds > zBoundMap
Definition: CalibUtil.h:107
#define LOG_VERBATIM(category)
unsigned int NPlanes() const
void FindZBoundaries(std::set< double > &planeZBounds)
Find the boundaries in the z direction of planes in the detector.
Definition: CalibUtil.cxx:332
Encapsulate the cell geometry.
Definition: CellGeo.h:25
Float_t X
Definition: plot.C:38
bool IsBad(int plane, int cell)
Float_t w
Definition: plot.C:20
Provide a simple utility function for our sets.
Definition: CalibUtil.cxx:44
Encapsulate the geometry of one entire detector (near, far, ndos)
bool HasZAdjacents(int plane, int cell, const ChannelSet &hmap, const art::ServiceHandle< geo::Geometry > geom)
Definition: CalibUtil.cxx:87
surf
Definition: demo5.py:35
float PathLengthInCell(zBoundMap const &zBounds, TVector3 const &recoHitLoc, std::pair< uint32_t, uint32_t > const &pc)
Return the path length of a track in the cell in question.
Definition: CalibUtil.cxx:498