All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
icaruscode/icaruscode/CRT/CRTUtils/CRTT0MatchAlg.cc
Go to the documentation of this file.
1 #include "CRTT0MatchAlg.h"
2 #include "larcore/CoreUtils/ServiceUtil.h" // lar::providerFrom()
3 
4 namespace icarus{
5 
6 
7  CRTT0MatchAlg::CRTT0MatchAlg(const fhicl::ParameterSet& pset)
8  {
9  this->reconfigure(pset);
10  return;
11  }
12 
13  CRTT0MatchAlg::CRTT0MatchAlg() = default;
14 
15 
16  void CRTT0MatchAlg::reconfigure(const fhicl::ParameterSet& pset){
17 
18  fMinTrackLength = pset.get<double>("MinTrackLength", 20.0);
19  fTrackDirectionFrac = pset.get<double>("TrackDirectionFrac", 0.5);
20  fDistanceLimit = pset.get<double>("DistanceLimit", 100);
21  fTSMode = pset.get<int>("TSMode", 2);
22  fTimeCorrection = pset.get<double>("TimeCorrection", 0.);
23  fSCEposCorr = pset.get<bool>("SCEposCorr", true);
24  fDirMethod = pset.get<int>("DirMethod", 1);
25  fDCAuseBox = pset.get<bool>("DCAuseBox",false);
26  fDCAoverLength = pset.get<bool>("DCAoverLength", false);
27  fDoverLLimit = pset.get<double>("DoverLLimit", 1);
28  fPEcut = pset.get<double>("PEcut", 0.0);
29  fMaxUncert = pset.get<double>("MaxUncert", 1000.);
30  fTPCTrackLabel = pset.get<std::vector<art::InputTag> >("TPCTrackLabel", {""});
31  // fDistEndpointAVedge = pset.get<double>(.DistEndpointAVedge();
32 
33  fGeometryService = lar::providerFrom<geo::Geometry>();//GeometryService;
34  fSCE = lar::providerFrom<spacecharge::SpaceChargeService>();
35  //fSCE = SCE;
36  return;
37 
38  }
39 
40 
41  matchCand makeNULLmc (){
42  sbn::crt::CRTHit hit;
43  matchCand null;
44  null.thishit = hit;
45  null.t0 = -99999;
46  null.dca = -99999;
47  null.extrapLen = -99999;
48  return null;
49  }
50 
51 
52 
53  // Utility function that determines the possible t0 range of a track
54  std::pair<double, double> CRTT0MatchAlg::TrackT0Range(detinfo::DetectorPropertiesData const& detProp,
55  double startX, double endX, int driftDirection,
56  std::pair<double, double> xLimits){
57 
58  // If track is stitched return zeros
59  if(driftDirection == 0) return std::make_pair(0, 0);
60 
61  //std::pair<double, double> result; // unused
62  double Vd = driftDirection * detProp.DriftVelocity();
63 
64  //std::cout << " [ driftdirn, vd ] = [ " << driftDirection << " , " << Vd << " ]" << std::endl;
65 
66  // Shift the most postive end to the most positive limit
67  double maxX = std::max(startX, endX);
68  double maxLimit = std::max(xLimits.first, xLimits.second);
69  double maxShift = maxLimit - maxX;
70  // Shift the most negative end to the most negative limit
71  double minX = std::min(startX, endX);
72  double minLimit = std::min(xLimits.first, xLimits.second);
73  double minShift = minLimit - minX;
74  // Convert to time
75  double t0max = maxShift/Vd;
76  double t0min = minShift/Vd;
77 
78  /*
79  std::cout << "[ driftdirn, vd, startx , endx, xlimits, xlimite, maxx, minx, maxl, minl, maxs, mins, t0max, t0min ] = [ "
80  << driftDirection << " , " << Vd << " , " << startX << " ," << endX << " ," << xLimits.first << " ," << xLimits.second
81  << " ," << maxX <<" ," << minX << " ," <<maxLimit << " ," << minLimit << " ," <<maxShift << " ," <<minShift
82  << " ," << t0max << " ," << t0min << " ]"<< std::endl;
83  */
84  // if (t0min>2500) std::cout << " t0 min " << t0min << " t0max " << t0max << std::endl;
85  return std::make_pair(std::min(t0min, t0max), std::max(t0min, t0max));
86 
87 
88  } // CRTT0MatchAlg::TrackT0Range()
89 
90 
91  double CRTT0MatchAlg::DistOfClosestApproach(detinfo::DetectorPropertiesData const& detProp,
92  TVector3 trackPos, TVector3 trackDir,
93  sbn::crt::CRTHit crtHit, int driftDirection, double t0){
94 
95  //double minDist = 99999;
96 
97  // Convert the t0 into an x shift
98  double xshift = driftDirection* t0 * detProp.DriftVelocity();
99  trackPos[0] += xshift;
100 
101  if (fSCE->EnableCalSpatialSCE() && fSCEposCorr) {
102  geo::Point_t temppt = {trackPos.X(),trackPos.Y(),trackPos.Z()};
103  geo::TPCID tpcid = fGeometryService->PositionToTPCID(temppt);
104  geo::Vector_t fPosOffsets = fSCE->GetCalPosOffsets(temppt,tpcid.TPC);
105  trackPos[0] += fPosOffsets.X();
106  trackPos[1] += fPosOffsets.Y();
107  trackPos[2] += fPosOffsets.Z();
108  }
109 
110  TVector3 end = trackPos + trackDir;
111  /*
112  std::cout << "[trackPosx, y, z, trackDirx, y, z, endx, y, z ] = [ "
113  << trackPos.X() << " , " << trackPos.Y() << " , " << trackPos.Z() << " , "
114  << trackDir.X() << " , " << trackDir.Y() << " , " << trackDir.Z() << " , "
115  << end.X() << " , " << end.Y() << " , " << end.Z() << " ]" << std::endl;
116  */
117  //-------- ADDED BY ME----------
118  TVector3 pos (crtHit.x_pos, crtHit.y_pos, crtHit.z_pos);
119  //double denominator = trackDir.Mag();
120  //double numerator = (pos - trackPos).Cross(pos - end).Mag();
121  /*
122  std::cout << "[ crt x, y, z, startx, directionx, endx, num, denom, distance, altnative, altdist ] = [ "
123  << crtHit.x_pos << " , " << crtHit.y_pos << " , " << crtHit.z_pos << " , "
124  << trackPos.X() << " , " << trackDir.X() << " , " << end.X() << " , "
125  << numerator << " , " << denominator << " , " << numerator/denominator << " , "
126  << (pos - trackPos).Cross(trackDir).Mag() << " , " << ( pos - trackPos).Cross(trackDir).Mag()/denominator << " ] " << std::endl;
127  */
128  //-----------------------
129 
130  // calculate distance of closest approach (DCA)
131  // default is the distance to the point specified by the CRT hit (Simple DCA)
132  // useBox is the distance to the closest edge of the rectangle with the CRT hit at the center and the sides defined
133  // the position uncertainties on the CRT hits.
134  double thisdca;
135 
136  if (fDCAuseBox) thisdca = DistToCrtHit(crtHit, trackPos, end);
137  else thisdca = SimpleDCA(crtHit, trackPos, trackDir);
138  return thisdca;
139 
140  } // CRTT0MatchAlg::DistToOfClosestApproach()
141 
142 
143  std::pair<TVector3, TVector3> CRTT0MatchAlg::TrackDirectionAverage(recob::Track track, double frac)
144  {
145  // Calculate direction as an average over directions
146  size_t nTrackPoints = track.NumberTrajectoryPoints();
147  recob::TrackTrajectory trajectory = track.Trajectory();
148  std::vector<geo::Vector_t> validDirections;
149  for(size_t i = 0; i < nTrackPoints; i++){
150  if(trajectory.FlagsAtPoint(i)!=recob::TrajectoryPointFlags::InvalidHitIndex) continue;
151  validDirections.push_back(track.DirectionAtPoint(i));
152  }
153 
154  size_t nValidPoints = validDirections.size();
155  int endPoint = (int)floor(nValidPoints*frac);
156  double xTotStart = 0; double yTotStart = 0; double zTotStart = 0;
157  double xTotEnd = 0; double yTotEnd = 0; double zTotEnd = 0;
158  for(int i = 0; i < endPoint; i++){
159  geo::Vector_t dirStart = validDirections.at(i);
160  geo::Vector_t dirEnd = validDirections.at(nValidPoints - (i+1));
161  xTotStart += dirStart.X();
162  yTotStart += dirStart.Y();
163  zTotStart += dirStart.Z();
164  xTotEnd += dirEnd.X();
165  yTotEnd += dirEnd.Y();
166  zTotEnd += dirEnd.Z();
167  }
168  TVector3 startDir = {-xTotStart/endPoint, -yTotStart/endPoint, -zTotStart/endPoint};
169  TVector3 endDir = {xTotEnd/endPoint, yTotEnd/endPoint, zTotEnd/endPoint};
170 
171  return std::make_pair(startDir, endDir);
172 
173  } // CRTT0MatchAlg::TrackDirectionAverage()
174 
175 
176  std::pair<TVector3, TVector3> CRTT0MatchAlg::TrackDirection(detinfo::DetectorPropertiesData const& detProp,
177  recob::Track track, double frac,
178  double CRTtime, int driftDirection){
179 
180  size_t nTrackPoints = track.NPoints();
181  int midPt = (int)floor(nTrackPoints*frac);
182  geo::Point_t startP = track.Start();
183  geo::Point_t endP = track.End();
184  geo::Point_t midP = track.LocationAtPoint(midPt);
185 
186  double xshift = driftDirection * CRTtime * detProp.DriftVelocity();
187  TVector3 startPoint = {startP.X()+xshift,startP.Y(),startP.Z()};
188  TVector3 endPoint = {endP.X()+xshift,endP.Y(),endP.Z()};
189  TVector3 midPoint = {midP.X()+xshift,midP.Y(),midP.Z()};
190 
191  // std::cout <<"[ nTrackPoints, midPt, startP, endP, midP, xshift, CRTtime, startPoint, endPoint, midPoint ] = [ "
192  // << nTrackPoints << " , " << midPt << " , " << startP.X() << " , " << endP.X() << " , " << midP.X() << " , "
193  // << xshift << " , " << CRTtime<< " , " << startPoint.X() << " , " << endPoint.X() << " , " << midPoint.X() << " ]"<<std::endl;
194 
195  if (fSCE->EnableCalSpatialSCE() && fSCEposCorr) {
196 
197  // Apply the shift depending on which TPC the track is in
198  geo::Point_t fTrackPos = startP;
199  //std::cout <<" before set fTrackPos " << fTrackPos.X() <<std::endl;
200  fTrackPos.SetX(startPoint.X());
201 
202  geo::TPCID tpcid = fGeometryService->PositionToTPCID(fTrackPos);
203  geo::Vector_t fPosOffsets = fSCE->GetCalPosOffsets(geo::Point_t{fTrackPos.X(),fTrackPos.Y(),fTrackPos.Z()},tpcid.TPC);
204 
205  startPoint.SetX(fTrackPos.X() + fPosOffsets.X());
206  startPoint.SetY(fTrackPos.Y() + fPosOffsets.Y());
207  startPoint.SetZ(fTrackPos.Z() + fPosOffsets.Z());
208  // std::cout <<" [ after set fTrackPos, tpcid, offset x, offset y, offset z, startx, starty, startz ] = [ " << fTrackPos.X() << " , " << tpcid.TPC
209  // << " , " << fPosOffsets.X() << " , " << fPosOffsets.Y() << " , " << fPosOffsets.Z()
210  // << " , " << startPoint.X()<< " , " << startPoint.Y() << " , " << startPoint.Z() << " ]" <<std::endl;
211  fTrackPos = endP;
212  fTrackPos.SetX(endPoint.X());
213  tpcid = fGeometryService->PositionToTPCID(fTrackPos);
214  // fPosOffsets = fSCE->GetCalPosOffsets(fTrackPos,tpcid.TPC);
215  fPosOffsets = fSCE->GetCalPosOffsets(geo::Point_t{fTrackPos.X(),fTrackPos.Y(),fTrackPos.Z()},tpcid.TPC);
216  endPoint.SetX(fTrackPos.X() + fPosOffsets.X());
217  endPoint.SetY(fTrackPos.Y() + fPosOffsets.Y());
218  endPoint.SetZ(fTrackPos.Z() + fPosOffsets.Z());
219  //std::cout <<" [ after set end fTrackPos, tpcid, offset x, offset y, offset z, startx, starty, startz ] = [ " << fTrackPos.X() << " , " << tpcid.TPC
220  // << " , " << fPosOffsets.X() << " , " << fPosOffsets.Y() << " , " << fPosOffsets.Z()
221  // << " , " << endPoint.X()<< " , " << endPoint.Y() << " , " << endPoint.Z() << " ]" <<std::endl;
222 
223  fTrackPos = midP;
224  fTrackPos.SetX(midPoint.X());
225  tpcid = fGeometryService->PositionToTPCID(fTrackPos);
226  //fPosOffsets = fSCE->GetCalPosOffsets(fTrackPos,tpcid.TPC);
227  fPosOffsets = fSCE->GetCalPosOffsets(geo::Point_t{fTrackPos.X(),fTrackPos.Y(),fTrackPos.Z()},tpcid.TPC);
228  midPoint.SetX(fTrackPos.X() + fPosOffsets.X());
229  midPoint.SetY(fTrackPos.Y() + fPosOffsets.Y());
230  midPoint.SetZ(fTrackPos.Z() + fPosOffsets.Z());
231  //std::cout <<" [ after set mid fTrackPos, tpcid, offset x, offset y, offset z, startx, starty, startz ] = [ " << fTrackPos.X() << " , " << tpcid.TPC
232  // << " , " << fPosOffsets.X() << " , " << fPosOffsets.Y() << " , " << fPosOffsets.Z()
233  // << " , " << midPoint.X()<< " , " << midPoint.Y() << " , " << midPoint.Z() << " ]" <<std::endl;
234  }
235 
236  TVector3 startDir = {midPoint.X()-startPoint.X(),midPoint.Y()-startPoint.Y(),midPoint.Z()-startPoint.Z()};
237  float norm = startDir.Mag();
238  if (norm>0) startDir *=(1.0/norm);
239  /*
240  std::cout <<" [ startDirx, startDiry, startDirz, mag, xcap, ycap, zcap ] = [ "
241  << midPoint.X()-startPoint.X() << " , " << midPoint.Y()-startPoint.Y() << " , " << midPoint.Z()-startPoint.Z()
242  << " , " << norm << " , " << startDir.X()<< " , " << startDir.Y() << " , " << startDir.Z() << " ]" <<std::endl;
243  */
244  TVector3 endDir = {midPoint.X()-endPoint.X(),midPoint.Y()-endPoint.Y(),midPoint.Z()-endPoint.Z()};
245  norm = endDir.Mag();
246  if (norm>0) endDir *=(1.0/norm);
247  /*
248  std::cout <<" [ endDirx, endDiry, endDirz, mag, xcap, ycap, zcap ] = [ "
249  << midPoint.X()-endPoint.X() << " , " << midPoint.Y()-endPoint.Y() << " , " << midPoint.Z()-endPoint.Z()
250  << " , " << norm<< " , " << endDir.X()<< " , " << endDir.Y() << " , " << endDir.Z() << " ]" <<std::endl;
251  */
252  return std::make_pair(startDir, endDir);
253 
254  } // CRTT0MatchAlg::TrackDirection()
255 
256  std::pair<TVector3, TVector3> CRTT0MatchAlg::TrackDirectionAverageFromPoints(recob::Track track, double frac){
257 
258  // Calculate direction as an average over directions
259  size_t nTrackPoints = track.NumberTrajectoryPoints();
260  recob::TrackTrajectory trajectory = track.Trajectory();
261  std::vector<TVector3> validPoints;
262  for(size_t i = 0; i < nTrackPoints; i++){
263  if(trajectory.FlagsAtPoint(i) != recob::TrajectoryPointFlags::InvalidHitIndex) continue;
264  validPoints.push_back(track.LocationAtPoint<TVector3>(i));
265  }
266 
267  size_t nValidPoints = validPoints.size();
268  int endPoint = (int)floor(nValidPoints*frac);
269  TVector3 startDir = validPoints.at(0) - validPoints.at(endPoint-1);
270  TVector3 endDir = validPoints.at(nValidPoints - 1) - validPoints.at(nValidPoints - (endPoint));
271 
272  return std::make_pair(startDir.Unit(), endDir.Unit());
273 
274  } // CRTT0MatchAlg::TrackDirectionAverageFromPoints()
275 
276 
277  // Keeping ClosestCRTHit function for backward compatibility only
278  // *** use GetClosestCRTHit instead
279 
280  std::vector<std::pair<sbn::crt::CRTHit, double> > CRTT0MatchAlg::ClosestCRTHit(detinfo::DetectorPropertiesData const& detProp,
281  recob::Track tpcTrack, std::vector<sbn::crt::CRTHit> crtHits,
282  const art::Event& event, uint64_t trigger_timestamp) {
283  // matchCand newmc = makeNULLmc();
284  std::vector<std::pair<sbn::crt::CRTHit, double> > crthitpair;
285 
286  for(const auto& trackLabel : fTPCTrackLabel){
287  auto tpcTrackHandle = event.getValidHandle<std::vector<recob::Track>>(trackLabel);
288  if (!tpcTrackHandle.isValid()) continue;
289 
290  art::FindManyP<recob::Hit> findManyHits(tpcTrackHandle, event, trackLabel);
291  for (auto const& tpcTrack : (*tpcTrackHandle)){
292  std::vector<art::Ptr<recob::Hit>> hits = findManyHits.at(tpcTrack.ID());
293 
294  crthitpair.push_back(ClosestCRTHit(detProp, tpcTrack, hits, crtHits, trigger_timestamp));
295  // return ClosestCRTHit(detProp, tpcTrack, hits, crtHits);
296  }
297  }
298 
299  return crthitpair;
300  //for(const auto& crthit : crthitpair)
301  //return std::make_pair(crthit.first, crthit.second);
302 
303  //return std::make_pair( newmc.thishit, -9999);
304  }
305 
306 
307  std::pair<sbn::crt::CRTHit, double> CRTT0MatchAlg::ClosestCRTHit(detinfo::DetectorPropertiesData const& detProp,
308  recob::Track tpcTrack, std::vector<art::Ptr<recob::Hit>> hits,
309  std::vector<sbn::crt::CRTHit> crtHits, uint64_t trigger_timestamp) {
310 
311  auto start = tpcTrack.Vertex<TVector3>();
312  auto end = tpcTrack.End<TVector3>();
313  // Get the drift direction from the TPC
314  int driftDirection = TPCGeoUtil::DriftDirectionFromHits(fGeometryService, hits);
315  std::pair<double, double> xLimits = TPCGeoUtil::XLimitsFromHits(fGeometryService, hits);
316  // Get the allowed t0 range
317  std::pair<double, double> t0MinMax = TrackT0Range(detProp, start.X(), end.X(), driftDirection, xLimits);
318 
319  return ClosestCRTHit(detProp, tpcTrack, t0MinMax, crtHits, driftDirection, trigger_timestamp);
320  }
321 
322  std::pair<sbn::crt::CRTHit, double> CRTT0MatchAlg::ClosestCRTHit(detinfo::DetectorPropertiesData const& detProp,
323  recob::Track tpcTrack, std::pair<double, double> t0MinMax,
324  std::vector<sbn::crt::CRTHit> crtHits, int driftDirection, uint64_t trigger_timestamp) {
325 
326  matchCand bestmatch = GetClosestCRTHit(detProp, tpcTrack,t0MinMax,crtHits,driftDirection, trigger_timestamp);
327  return std::make_pair(bestmatch.thishit,bestmatch.dca);
328 
329  }
330 
331 
332  matchCand CRTT0MatchAlg::GetClosestCRTHit(detinfo::DetectorPropertiesData const& detProp,
333  recob::Track tpcTrack, std::vector<art::Ptr<recob::Hit>> hits,
334  std::vector<sbn::crt::CRTHit> crtHits, uint64_t trigger_timestamp) {
335 
336  auto start = tpcTrack.Vertex<TVector3>();
337  auto end = tpcTrack.End<TVector3>();
338 
339 
340 
341  // Get the drift direction from the TPC
342  int driftDirection = TPCGeoUtil::DriftDirectionFromHits(fGeometryService, hits);
343  //std::cout << "size of hit in a track: " << hits.size() << ", driftDirection: "<< driftDirection
344  // << " , tpc: "<< hits[0]->WireID().TPC << std::endl; //<< " , intpc: "<< icarus::TPCGeoUtil::DetectedInTPC(hits) << std::endl;
345  std::pair<double, double> xLimits = TPCGeoUtil::XLimitsFromHits(fGeometryService, hits);
346  // Get the allowed t0 range
347  std::pair<double, double> t0MinMax = TrackT0Range(detProp, start.X(), end.X(), driftDirection, xLimits);
348 
349  return GetClosestCRTHit(detProp, tpcTrack, t0MinMax, crtHits, driftDirection, trigger_timestamp);
350 
351  }
352 
353  std::vector<matchCand> CRTT0MatchAlg::GetClosestCRTHit(detinfo::DetectorPropertiesData const& detProp,
354  recob::Track tpcTrack, std::vector<sbn::crt::CRTHit> crtHits,
355  const art::Event& event, uint64_t trigger_timestamp) {
356  // matchCand nullmatch = makeNULLmc();
357  std::vector<matchCand> matchcanvec;
358  //std::vector<std::pair<sbn::crt::CRTHit, double> > matchedCan;
359  for(const auto& trackLabel : fTPCTrackLabel){
360  auto tpcTrackHandle = event.getValidHandle<std::vector<recob::Track>>(trackLabel);
361  if (!tpcTrackHandle.isValid()) continue;
362 
363  art::FindManyP<recob::Hit> findManyHits(tpcTrackHandle, event, trackLabel);
364  for (auto const& tpcTrack : (*tpcTrackHandle)){
365  std::vector<art::Ptr<recob::Hit>> hits = findManyHits.at(tpcTrack.ID());
366  matchcanvec.push_back(GetClosestCRTHit(detProp, tpcTrack, hits, crtHits, trigger_timestamp));
367  //return ClosestCRTHit(detProp, tpcTrack, hits, crtHits);
368  //matchCand closestHit = GetClosestCRTHit(detProp, tpcTrack, hits, crtHits);
369 
370  }
371  }
372  return matchcanvec;
373  //auto tpcTrackHandle = event.getValidHandle<std::vector<recob::Track>>(fTPCTrackLabel);
374  //art::FindManyP<recob::Hit> findManyHits(tpcTrackHandle, event, fTPCTrackLabel);
375  //std::vector<art::Ptr<recob::Hit>> hits = findManyHits.at(tpcTrack.ID());
376  //return GetClosestCRTHit(detProp, tpcTrack, hits, crtHits);
377  // for (const auto& match : matchedCan)
378  //return match;
379  //return nullmatch;
380  }
381 
382 
383  matchCand CRTT0MatchAlg::GetClosestCRTHit(detinfo::DetectorPropertiesData const& detProp,
384  recob::Track tpcTrack, std::pair<double, double> t0MinMax,
385  std::vector<sbn::crt::CRTHit> crtHits, int driftDirection, uint64_t& trigger_timestamp) {
386 
387  auto start = tpcTrack.Vertex<TVector3>();
388  auto end = tpcTrack.End<TVector3>();
389 
390  // ====================== Matching Algorithm ========================== //
391  // std::vector<std::pair<sbn::crt::CRTHit, double>> t0Candidates;
392  std::vector<matchCand> t0Candidates;
393 
394  // if (crtHits.size() == 0) continue;
395  // Loop over all the CRT hits
396  for(auto &crtHit : crtHits){
397  // Check if hit is within the allowed t0 range
398  double crtTime = -99999.; // units are us
399  if (fTSMode == 1) {
400  crtTime = ((double)(int)crtHit.ts1_ns) * 1e-3 + fTimeCorrection;
401  }
402  else {
403  //std::cout << "trigger_timestamp: "<< trigger_timestamp << " , t0 " << (uint64_t)crtHit.ts0_ns << std::endl;
404  crtTime = double(crtHit.ts0_ns - trigger_timestamp%1'000'000'000)/1e3;
405  //'
406  if(crtTime<-0.5e6) crtTime+=1e6;
407  else if(crtTime>0.5e6) crtTime-=1e6;
408  //'// std::cout << "(trigger - t0)/1e3: " << crtTime << std::endl;
409  //crtTime = -crtTime+1e6;
410  //std::cout << "-crtTime+1e6: " << crtTime << std::endl;
411  //crtTime = ((double)(int)crtHit.ts0_ns) * 1e-3 + fTimeCorrection;
412  }
413  // if (crtTime < 3000 && crtTime > -3000) std::cout << "crt hit times " << crtTime << std::endl;
414  // std::cout << "[ tpc t0 min , tpc t0 max ] = [ "<< t0MinMax.first << " , " << t0MinMax.second << " ]" << std::endl;
415  // If track is stitched then try all hits
416  if (!((crtTime >= t0MinMax.first - 10. && crtTime <= t0MinMax.second + 10.)
417  || t0MinMax.first == t0MinMax.second)) continue;
418 
419  //std::cout << "[ tpc t0 min , tpc t0 max, crttime ] = [ "<< t0MinMax.first << " , " << t0MinMax.second
420  // << " , " << crtTime << " ]" << std::endl;
421 
422  //std::cout << "passed ....................... " << std::endl;
423 
424  // cut on CRT hit PE value
425  if (crtHit.peshit<fPEcut) continue;
426  if (crtHit.x_err>fMaxUncert) continue;
427  if (crtHit.y_err>fMaxUncert) continue;
428  if (crtHit.z_err>fMaxUncert) continue;
429 
430  TVector3 crtPoint(crtHit.x_pos, crtHit.y_pos, crtHit.z_pos);
431 
432  //std::cout << "[ tpc t0 min , tpc t0 max, crttime, crtx, crty, crtz ] = [ "<< t0MinMax.first << " , " << t0MinMax.second
433  // << " , " << crtTime << " , " << crtHit.x_pos << " , " << crtHit.y_pos << " , " << crtHit.z_pos << " ]" << std::endl;
434  //Calculate Track direction
435  std::pair<TVector3, TVector3> startEndDir;
436  // dirmethod=2 is original algorithm, dirmethod=1 is simple algorithm for which SCE corrections are possible
437  if (fDirMethod==2) startEndDir = TrackDirectionAverage(tpcTrack, fTrackDirectionFrac);
438  else startEndDir = TrackDirection(detProp, tpcTrack, fTrackDirectionFrac, crtTime, driftDirection);
439  TVector3 startDir = startEndDir.first;
440  TVector3 endDir = startEndDir.second;
441 
442  // Calculate the distance between the crossing point and the CRT hit, SCE corrections are done inside but dropped
443  double startDist = DistOfClosestApproach(detProp, start, startDir, crtHit, driftDirection, crtTime);
444  double endDist = DistOfClosestApproach(detProp, end, endDir, crtHit, driftDirection, crtTime);
445 
446 
447  double xshift = driftDirection * crtTime * detProp.DriftVelocity();
448  auto thisstart = start;
449  thisstart.SetX(start.X()+xshift);
450  auto thisend = end;
451  thisend.SetX(end.X()+xshift);
452 
453  // repeat SCE correction for endpoints
454  if (fSCE->EnableCalSpatialSCE() && fSCEposCorr) {
455  geo::Point_t temppt = {thisstart.X(),thisstart.Y(),thisstart.Z()};
456  geo::TPCID tpcid = fGeometryService->PositionToTPCID(temppt);
457  geo::Vector_t fPosOffsets = fSCE->GetCalPosOffsets(temppt,tpcid.TPC);
458  thisstart[0] += fPosOffsets.X();
459  thisstart[1] += fPosOffsets.Y();
460  thisstart[2] += fPosOffsets.Z();
461  temppt.SetX(thisend.X());
462  temppt.SetY(thisend.Y());
463  temppt.SetZ(thisend.Z());
464  tpcid = fGeometryService->PositionToTPCID(temppt);
465  fPosOffsets = fSCE->GetCalPosOffsets(temppt,tpcid.TPC);
466  thisend[0] += fPosOffsets.X();
467  thisend[1] += fPosOffsets.Y();
468  thisend[2] += fPosOffsets.Z();
469  }
470 
471 
472  matchCand newmc = makeNULLmc();
473  if (startDist<fDistanceLimit || endDist<fDistanceLimit) {
474  double distS = (crtPoint-thisstart).Mag();
475  double distE = (crtPoint-thisend).Mag();
476  // std::cout << " distS " << distS << " distE " << distE << std::endl;
477  // std::cout << "startdis " << startDist << " endDist " << endDist << " dca " << std::endl;
478  // std::cout << " doL start " << startDist/distS << " doL end " << endDist/distE << std::endl;
479  /*
480  std::cout << " distS " << distS << " distE " << distE
481  << "startdis " << startDist << " endDist " << endDist << " dca " << std::endl;
482  */
483  if (distS < distE){
484  newmc.thishit = crtHit;
485  newmc.t0= crtTime;
486  newmc.dca = startDist;
487  newmc.extrapLen = distS;
488  t0Candidates.push_back(newmc);
489  // std::cout << " hello inside the distS < distE, found " << t0Candidates.size() << " candidates"<< std::endl;
490  }
491  else{
492  newmc.thishit = crtHit;
493  newmc.t0= crtTime;
494  newmc.dca = endDist;
495  newmc.extrapLen = distE;
496  t0Candidates.push_back(newmc);
497  //std::cout << " hello outside the distS < distE, found " << t0Candidates.size() << " candidates"<< std::endl;
498  }
499  }
500  }
501 
502 
503  //std::cout << " found " << t0Candidates.size() << " candidates" << std::endl;
504  matchCand bestmatch = makeNULLmc();
505  if(t0Candidates.size() > 0){
506  // Find candidate with shortest DCA or DCA/L value
507  bestmatch=t0Candidates[0];
508  double sin_angle = bestmatch.dca/bestmatch.extrapLen;
509  if (fDCAoverLength) { // Use dca/extrapLen to judge best
510  for(auto &thisCand : t0Candidates){
511  double this_sin_angle = thisCand.dca/thisCand.extrapLen;
512  if (bestmatch.dca<0 )bestmatch=thisCand;
513  else if (this_sin_angle<sin_angle && thisCand.dca>=0)bestmatch=thisCand;
514  }
515  }
516  else { // use Dca to judge best
517  for(auto &thisCand : t0Candidates){
518  //std::cout << "[bestmatch, thiscand] = [ " << bestmatch.dca << " , " << thisCand.dca << " ] " << std::endl;
519  if (bestmatch.dca<0 )bestmatch=thisCand;
520  else if (thisCand.dca<bestmatch.dca && thisCand.dca>=0)bestmatch=thisCand;
521  }
522  }
523  }
524 
525  //std::cout << "best match has dca of " << bestmatch.dca << std::endl;
526  return bestmatch;
527 
528  }
529 
530 
531  std::vector<double> CRTT0MatchAlg::T0FromCRTHits(detinfo::DetectorPropertiesData const& detProp,
532  recob::Track tpcTrack, std::vector<sbn::crt::CRTHit> crtHits,
533  const art::Event& event, uint64_t trigger_timestamp){
534  std::vector<double> ftime;
535  for(const auto& trackLabel : fTPCTrackLabel){
536  auto tpcTrackHandle = event.getValidHandle<std::vector<recob::Track>>(trackLabel);
537  if (!tpcTrackHandle.isValid()) continue;
538 
539  art::FindManyP<recob::Hit> findManyHits(tpcTrackHandle, event, trackLabel);
540  for (auto const& tpcTrack : (*tpcTrackHandle)){
541  std::vector<art::Ptr<recob::Hit>> hits = findManyHits.at(tpcTrack.ID());
542  ftime.push_back(T0FromCRTHits(detProp, tpcTrack, hits, crtHits, trigger_timestamp));
543  // return T0FromCRTHits(detProp, tpcTrack, hits, crtHits);
544  }
545  }
546  return ftime;
547  //auto tpcTrackHandle = event.getValidHandle<std::vector<recob::Track>>(fTPCTrackLabel);
548  //art::FindManyP<recob::Hit> findManyHits(tpcTrackHandle, event, fTPCTrackLabel);
549  //std::vector<art::Ptr<recob::Hit>> hits = findManyHits.at(tpcTrack.ID());
550  //return T0FromCRTHits(detProp, tpcTrack, hits, crtHits);
551 
552  // for(const auto& t0 : ftime) return t0;
553  //return -99999;
554  }
555 
556  double CRTT0MatchAlg::T0FromCRTHits(detinfo::DetectorPropertiesData const& detProp,
557  recob::Track tpcTrack, std::vector<art::Ptr<recob::Hit>> hits,
558  std::vector<sbn::crt::CRTHit> crtHits, uint64_t& trigger_timestamp) {
559 
560  if (tpcTrack.Length() < fMinTrackLength) return -99999;
561 
562  matchCand closestHit = GetClosestCRTHit(detProp, tpcTrack, hits, crtHits, trigger_timestamp);
563  if(closestHit.dca <0) return -99999;
564 
565  double crtTime;
566  if (fTSMode == 1) {
567  crtTime = ((double)(int)closestHit.thishit.ts1_ns) * 1e-3 + fTimeCorrection;
568  }
569  else {
570  crtTime = ((double)(int)closestHit.thishit.ts0_ns) * 1e-3 + fTimeCorrection;
571  }
572  if (closestHit.dca < fDistanceLimit && (closestHit.dca/closestHit.extrapLen) < fDoverLLimit) return crtTime;
573 
574  return -99999;
575 
576  }
577 
578  std::vector<std::pair<double, double> > CRTT0MatchAlg::T0AndDCAFromCRTHits(detinfo::DetectorPropertiesData const& detProp,
579  recob::Track tpcTrack, std::vector<sbn::crt::CRTHit> crtHits,
580  const art::Event& event, uint64_t trigger_timestamp){
581 
582  std::vector<std::pair<double, double> > ft0anddca;
583  for(const auto& trackLabel : fTPCTrackLabel){
584  auto tpcTrackHandle = event.getValidHandle<std::vector<recob::Track>>(trackLabel);
585  if (!tpcTrackHandle.isValid()) continue;
586 
587  art::FindManyP<recob::Hit> findManyHits(tpcTrackHandle, event, trackLabel);
588  for (auto const& tpcTrack : (*tpcTrackHandle)){
589  std::vector<art::Ptr<recob::Hit>> hits = findManyHits.at(tpcTrack.ID());
590  ft0anddca.push_back(T0AndDCAFromCRTHits(detProp, tpcTrack, hits, crtHits, trigger_timestamp));
591  // return T0AndDCAFromCRTHits(detProp, tpcTrack, hits, crtHits);
592  }
593  }
594  return ft0anddca;
595  // auto tpcTrackHandle = event.getValidHandle<std::vector<recob::Track>>(fTPCTrackLabel);
596  //art::FindManyP<recob::Hit> findManyHits(tpcTrackHandle, event, fTPCTrackLabel);
597  //std::vector<art::Ptr<recob::Hit>> hits = findManyHits.at(tpcTrack.ID());
598  //return T0AndDCAFromCRTHits(detProp, tpcTrack, hits, crtHits);
599  // for(const auto& t0anddca : ft0anddca) return std::make_pair(t0anddca.first, t0anddca.second);
600  //return std::make_pair(-9999., -9999.);
601  }
602 
603  std::pair<double, double> CRTT0MatchAlg::T0AndDCAFromCRTHits(detinfo::DetectorPropertiesData const& detProp,
604  recob::Track tpcTrack, std::vector<art::Ptr<recob::Hit>> hits,
605  std::vector<sbn::crt::CRTHit> crtHits, uint64_t& trigger_timestamp) {
606 
607  if (tpcTrack.Length() < fMinTrackLength) return std::make_pair(-9999., -9999.);
608 
609  matchCand closestHit = GetClosestCRTHit(detProp, tpcTrack, hits, crtHits, trigger_timestamp);
610 
611  if(closestHit.dca < 0 ) return std::make_pair(-9999., -9999.);
612  if (closestHit.dca < fDistanceLimit && (closestHit.dca/closestHit.extrapLen) < fDoverLLimit) return std::make_pair(closestHit.t0, closestHit.dca);
613 
614  return std::make_pair(-9999., -9999.);
615 
616 
617  }
618 
619  // Simple distance of closest approach between infinite track and centre of hit
620  double CRTT0MatchAlg::SimpleDCA(sbn::crt::CRTHit hit, TVector3 start, TVector3 direction){
621 
622  TVector3 pos (hit.x_pos, hit.y_pos, hit.z_pos);
623  TVector3 end = start + direction;
624  double denominator = direction.Mag();
625  double numerator = (pos - start).Cross(pos - end).Mag();
626  return numerator/denominator;
627 
628  }
629 
630  // Minimum distance from infinite track to CRT hit assuming that hit is a 2D square
631  double CRTT0MatchAlg::DistToCrtHit(sbn::crt::CRTHit hit, TVector3 start, TVector3 end){
632 
633  // Check if track goes inside hit
634  TVector3 min (hit.x_pos - hit.x_err, hit.y_pos - hit.y_err, hit.z_pos - hit.z_err);
635  TVector3 max (hit.x_pos + hit.x_err, hit.y_pos + hit.y_err, hit.z_pos + hit.z_err);
636  if(CubeIntersection(min, max, start, end).first.X() != -99999) return 0;
637 
638  // Calculate the closest distance to each edge of the CRT hit
639  // Assume min error is the fixed position of tagger
640  TVector3 vertex1 (hit.x_pos, hit.y_pos - hit.y_err, hit.z_pos - hit.z_err);
641  TVector3 vertex2 (hit.x_pos, hit.y_pos + hit.y_err, hit.z_pos - hit.z_err);
642  TVector3 vertex3 (hit.x_pos, hit.y_pos - hit.y_err, hit.z_pos + hit.z_err);
643  TVector3 vertex4 (hit.x_pos, hit.y_pos + hit.y_err, hit.z_pos + hit.z_err);
644  if(hit.y_err < hit.x_err && hit.y_err < hit.z_err){
645  vertex1.SetXYZ(hit.x_pos - hit.x_err, hit.y_pos, hit.z_pos - hit.z_err);
646  vertex2.SetXYZ(hit.x_pos + hit.x_err, hit.y_pos, hit.z_pos - hit.z_err);
647  vertex3.SetXYZ(hit.x_pos - hit.x_err, hit.y_pos, hit.z_pos + hit.z_err);
648  vertex4.SetXYZ(hit.x_pos + hit.x_err, hit.y_pos, hit.z_pos + hit.z_err);
649  }
650  if(hit.z_err < hit.x_err && hit.z_err < hit.y_err){
651  vertex1.SetXYZ(hit.x_pos - hit.x_err, hit.y_pos - hit.y_err, hit.z_pos);
652  vertex2.SetXYZ(hit.x_pos + hit.x_err, hit.y_pos - hit.y_err, hit.z_pos);
653  vertex3.SetXYZ(hit.x_pos - hit.x_err, hit.y_pos + hit.y_err, hit.z_pos);
654  vertex4.SetXYZ(hit.x_pos + hit.x_err, hit.y_pos + hit.y_err, hit.z_pos);
655  }
656 
657  double dist1 = LineSegmentDistance(vertex1, vertex2, start, end);
658  double dist2 = LineSegmentDistance(vertex1, vertex3, start, end);
659  double dist3 = LineSegmentDistance(vertex4, vertex2, start, end);
660  double dist4 = LineSegmentDistance(vertex4, vertex3, start, end);
661 
662  return std::min(std::min(dist1, dist2), std::min(dist3, dist4));
663 
664  }
665 
666 
667  // Distance between infinite line (2) and segment (1)
668  // http://geomalgorithms.com/a07-_distance.html
669  double CRTT0MatchAlg::LineSegmentDistance(TVector3 start1, TVector3 end1, TVector3 start2, TVector3 end2){
670 
671  double smallNum = 0.00001;
672 
673  // 1 is segment
674  TVector3 direction1 = end1 - start1;
675  // 2 is infinite line
676  TVector3 direction2 = end2 - start2;
677 
678  TVector3 u = direction1;
679  TVector3 v = direction2;
680  TVector3 w = start1 - start2;
681 
682  double a = u.Dot(u);
683  double b = u.Dot(v);
684  double c = v.Dot(v);
685  double d = u.Dot(w);
686  double e = v.Dot(w);
687  double D = a * c - b * b;
688  double sc, sN, sD = D; // sc = sN/sD
689  double tc, tN, tD = D; // sc = sN/sD
690 
691  // Compute the line parameters of the two closest points
692  if(D < smallNum){ // Lines are almost parallel
693  sN = 0.0;
694  sD = 1.0;
695  tN = e;
696  tD = c;
697  }
698  else{
699  sN = (b * e - c * d)/D;
700  tN = (a * e - b * d)/D;
701  if(sN < 0.){ // sc < 0, the s = 0 edge is visible
702  sN = 0.;
703  tN = e;
704  tD = c;
705  }
706  else if(sN > sD){ // sc > 1, the s = 1 edge is visible
707  sN = sD;
708  tN = e + b;
709  tD = c;
710  }
711  }
712 
713  sc = (std::abs(sN) < smallNum ? 0.0 : sN / sD);
714  tc = (std::abs(tN) < smallNum ? 0.0 : tN / tD);
715  // Get the difference of the two closest points
716  TVector3 dP = w + (sc * u) - (tc * v);
717 
718  return dP.Mag();
719 
720  }
721 
722  // Intersection between axis-aligned cube and infinite line
723  // (https://www.scratchapixel.com/lessons/3d-basic-rendering/minimal-ray-tracer-rendering-simple-shapes/ray-box-intersection)
724  std::pair<TVector3, TVector3> CRTT0MatchAlg::CubeIntersection(TVector3 min, TVector3 max, TVector3 start, TVector3 end){
725 
726  TVector3 dir = (end - start);
727  TVector3 invDir (1./dir.X(), 1./dir.Y(), 1/dir.Z());
728 
729  double tmin, tmax, tymin, tymax, tzmin, tzmax;
730 
731  TVector3 enter (-99999, -99999, -99999);
732  TVector3 exit (-99999, -99999, -99999);
733 
734  // Find the intersections with the X plane
735  if(invDir.X() >= 0){
736  tmin = (min.X() - start.X()) * invDir.X();
737  tmax = (max.X() - start.X()) * invDir.X();
738  }
739  else{
740  tmin = (max.X() - start.X()) * invDir.X();
741  tmax = (min.X() - start.X()) * invDir.X();
742  }
743 
744  // Find the intersections with the Y plane
745  if(invDir.Y() >= 0){
746  tymin = (min.Y() - start.Y()) * invDir.Y();
747  tymax = (max.Y() - start.Y()) * invDir.Y();
748  }
749  else{
750  tymin = (max.Y() - start.Y()) * invDir.Y();
751  tymax = (min.Y() - start.Y()) * invDir.Y();
752  }
753 
754  // Check that it actually intersects
755  if((tmin > tymax) || (tymin > tmax)) return std::make_pair(enter, exit);
756 
757  // Max of the min points is the actual intersection
758  if(tymin > tmin) tmin = tymin;
759 
760  // Min of the max points is the actual intersection
761  if(tymax < tmax) tmax = tymax;
762 
763  // Find the intersection with the Z plane
764  if(invDir.Z() >= 0){
765  tzmin = (min.Z() - start.Z()) * invDir.Z();
766  tzmax = (max.Z() - start.Z()) * invDir.Z();
767  }
768  else{
769  tzmin = (max.Z() - start.Z()) * invDir.Z();
770  tzmax = (min.Z() - start.Z()) * invDir.Z();
771  }
772 
773  // Check for intersection
774  if((tmin > tzmax) || (tzmin > tmax)) return std::make_pair(enter, exit);
775 
776  // Find final intersection points
777  if(tzmin > tmin) tmin = tzmin;
778 
779  // Find final intersection points
780  if(tzmax < tmax) tmax = tzmax;
781 
782  // Calculate the actual crossing points
783  double xmin = start.X() + tmin * dir.X();
784  double xmax = start.X() + tmax * dir.X();
785  double ymin = start.Y() + tmin * dir.Y();
786  double ymax = start.Y() + tmax * dir.Y();
787  double zmin = start.Z() + tmin * dir.Z();
788  double zmax = start.Z() + tmax * dir.Z();
789 
790  // Return pair of entry and exit points
791  enter.SetXYZ(xmin, ymin, zmin);
792  exit.SetXYZ(xmax, ymax, zmax);
793  return std::make_pair(enter, exit);
794 
795  }
796 
797 }
sbn::crt::CRTHit::z_err
float z_err
position uncertainty in z-direction (cm).
Definition: CRTHit.hh:43
geo::Vector_t
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space.
Definition: geo_vectors.h:164
sbn::crt::CRTHit::x_err
float x_err
position uncertainty in x-direction (cm).
Definition: CRTHit.hh:39
spacecharge::SpaceCharge::GetCalPosOffsets
virtual geo::Vector_t GetCalPosOffsets(geo::Point_t const &point, int const &TPCid) const =0
icarus::TPCGeoUtil::XLimitsFromHits
std::pair< double, double > XLimitsFromHits(const geo::GeometryCore *GeometryService, std::vector< art::Ptr< recob::Hit >> hits)
Definition: icaruscode/icaruscode/CRT/CRTUtils/TPCGeoUtil.cc:15
if
then if[["$THISISATEST"==1]]
Definition: neoSmazza.sh:95
ServiceUtil.h
Utilities related to art service access.
icarus::CRTT0MatchAlg::TrackDirectionAverage
std::pair< TVector3, TVector3 > TrackDirectionAverage(recob::Track track, double frac)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTT0MatchAlg.cc:143
sbn::crt::CRTHit::ts1_ns
double ts1_ns
Timestamp T1 ([signal time w.r.t. Trigger time]), in UTC absolute time scale in nanoseconds from the ...
Definition: CRTHit.hh:34
exit
exit
Definition: generate_icarus_geometry_splitwire.pl:430
recob::Track::Trajectory
const recob::TrackTrajectory & Trajectory() const
Access to the stored recob::TrackTrajectory.
Definition: bj/lardataobj/RecoBase/Track.h:98
sbn::crt::CRTHit::y_err
float y_err
position uncertainty in y-direction (cm).
Definition: CRTHit.hh:41
recob::Track::LocationAtPoint
Point_t const & LocationAtPoint(size_t i) const
Definition: bj/lardataobj/RecoBase/Track.h:126
recob::Track::NumberTrajectoryPoints
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: bj/lardataobj/RecoBase/Track.h:102
geo::GeometryCore::PositionToTPCID
geo::TPCID PositionToTPCID(geo::Point_t const &point) const
Returns the ID of the TPC at specified location.
Definition: GeometryCore.cxx:452
track
process_name use argoneut_mc_hitfinder track
Definition: standard_reco.fcl:36
icarus::CRTT0MatchAlg::reconfigure
void reconfigure(const fhicl::ParameterSet &pset)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTT0MatchAlg.cc:16
hit
process_name hit
Definition: cheaterreco.fcl:51
icarus::CRTT0MatchAlg::T0AndDCAFromCRTHits
std::vector< std::pair< double, double > > T0AndDCAFromCRTHits(detinfo::DetectorPropertiesData const &detProp, recob::Track tpcTrack, std::vector< sbn::crt::CRTHit > crtHits, const art::Event &event, uint64_t trigger_timestamp)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTT0MatchAlg.cc:578
icarus::TPCGeoUtil::DriftDirectionFromHits
int DriftDirectionFromHits(const geo::GeometryCore *GeometryService, std::vector< art::Ptr< recob::Hit >> hits)
Definition: icaruscode/icaruscode/CRT/CRTUtils/TPCGeoUtil.cc:28
a
process_name gaushit a
Definition: decode_signalprocess_icarus.fcl:37
xmin
process_name pandoraGausCryo1 vertexChargeCryo1 vertexStubCryo1 xmin
Definition: run_vertexcharge_icarus.fcl:109
sbn::crt::CRTHit::z_pos
float z_pos
position in z-direction (cm).
Definition: CRTHit.hh:42
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
make_calib_df.t0
int t0
Definition: make_calib_df.py:42
sbn::crt::CRTHit::ts0_ns
double ts0_ns
Timestamp T0 (from White Rabbit), in UTC absolute time scale in nanoseconds from the Epoch...
Definition: CRTHit.hh:32
icarus::CRTT0MatchAlg::ClosestCRTHit
std::pair< sbn::crt::CRTHit, double > ClosestCRTHit(detinfo::DetectorPropertiesData const &detProp, recob::Track tpcTrack, std::pair< double, double > t0MinMax, std::vector< sbn::crt::CRTHit > crtHits, int driftDirection, uint64_t trigger_timestamp)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTT0MatchAlg.cc:322
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
icarus::CRTT0MatchAlg::T0FromCRTHits
std::vector< double > T0FromCRTHits(detinfo::DetectorPropertiesData const &detProp, recob::Track tpcTrack, std::vector< sbn::crt::CRTHit > crtHits, const art::Event &event, uint64_t trigger_timestamp)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTT0MatchAlg.cc:531
recob::Track::Length
double Length(size_t p=0) const
Access to various track properties.
Definition: bj/lardataobj/RecoBase/Track.h:167
recob::TrackTrajectory::FlagsAtPoint
PointFlags_t const & FlagsAtPoint(size_t i) const
Returns the flags for the specified trajectory point.
Definition: TrackTrajectory.h:165
recob::Track::Start
Point_t const & Start() const
Access to track position at different points.
Definition: bj/lardataobj/RecoBase/Track.h:123
icarus::CRTT0MatchAlg::LineSegmentDistance
double LineSegmentDistance(TVector3 start1, TVector3 end1, TVector3 start2, TVector3 end2)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTT0MatchAlg.cc:669
recob::TrackTrajectory
A trajectory in space reconstructed from hits.
Definition: TrackTrajectory.h:62
icarus::CRTT0MatchAlg::TrackDirectionAverageFromPoints
std::pair< TVector3, TVector3 > TrackDirectionAverageFromPoints(recob::Track track, double frac)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTT0MatchAlg.cc:256
recob::Track::Vertex
Point_t const & Vertex() const
Definition: bj/lardataobj/RecoBase/Track.h:124
icarus::ns::util::end
auto end(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:585
detinfo::DetectorPropertiesData::DriftVelocity
double DriftVelocity(double efield=0., double temperature=0.) const
cm/us
Definition: DetectorPropertiesData.cc:23
geo::TPCID
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
icarus::CRTT0MatchAlg::DistToCrtHit
double DistToCrtHit(sbn::crt::CRTHit hit, TVector3 start, TVector3 end)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTT0MatchAlg.cc:631
recob::TrajectoryPointFlags::InvalidHitIndex
static constexpr HitIndex_t InvalidHitIndex
Value marking an invalid hit index.
Definition: TrajectoryPointFlags.h:336
gen_crt_frags.w
tuple w
Definition: gen_crt_frags.py:1416
geo::vect::norm
auto norm(Vector const &v)
Return norm of the specified vector.
Definition: geo_vectors_utils.h:1202
dropbox.dir
tuple dir
Definition: dropbox.py:28
sbn::crt::CRTHit::y_pos
float y_pos
position in y-direction (cm).
Definition: CRTHit.hh:40
icarus::CRTT0MatchAlg::TrackT0Range
std::pair< double, double > TrackT0Range(detinfo::DetectorPropertiesData const &detProp, double startX, double endX, int driftDirection, std::pair< double, double > xLimits)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTT0MatchAlg.cc:54
icarus::CRTT0MatchAlg::GetClosestCRTHit
matchCand GetClosestCRTHit(detinfo::DetectorPropertiesData const &detProp, recob::Track tpcTrack, std::pair< double, double > t0MinMax, std::vector< sbn::crt::CRTHit > crtHits, int driftDirection, uint64_t &trigger_timestamp)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTT0MatchAlg.cc:383
sbn::crt::CRTHit::x_pos
float x_pos
position in x-direction (cm).
Definition: CRTHit.hh:38
icarus::CRTT0MatchAlg::SimpleDCA
double SimpleDCA(sbn::crt::CRTHit hit, TVector3 start, TVector3 direction)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTT0MatchAlg.cc:620
e
do i e
Definition: updateVersionList.sh:68
recob::Track::End
Point_t const & End() const
Definition: bj/lardataobj/RecoBase/Track.h:125
icarus::CRTT0MatchAlg::DistOfClosestApproach
double DistOfClosestApproach(detinfo::DetectorPropertiesData const &detProp, TVector3 trackPos, TVector3 trackDir, sbn::crt::CRTHit crtHit, int driftDirection, double t0)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTT0MatchAlg.cc:91
recob::Track::DirectionAtPoint
Vector_t DirectionAtPoint(size_t i) const
Definition: bj/lardataobj/RecoBase/Track.h:134
icarus::CRTT0MatchAlg::TrackDirection
std::pair< TVector3, TVector3 > TrackDirection(detinfo::DetectorPropertiesData const &detProp, recob::Track track, double frac, double CRTtime, int driftDirection)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTT0MatchAlg.cc:176
geo::TPCID::TPC
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
spacecharge::SpaceCharge::EnableCalSpatialSCE
virtual bool EnableCalSpatialSCE() const =0
geo::Point_t
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:184
detsim::detProp
auto const detProp
Definition: SimWire_module.cc:157
recob::Track
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a &quot;fitted&quot; track:
Definition: bj/lardataobj/RecoBase/Track.h:49
icarus::CRTT0MatchAlg::CubeIntersection
std::pair< TVector3, TVector3 > CubeIntersection(TVector3 min, TVector3 max, TVector3 start, TVector3 end)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTT0MatchAlg.cc:724
Generated by   doxygen 1.8.5