ViewMatchAlg.cxx
Go to the documentation of this file.
1 ///
2 /// \file ViewMatchAlg.cxx
3 /// \brief 3d matching of prongs
4 /// \author eniner@indiana.edu
5 /// \version $Id: ViewMatchAlg.cxx,v 1.6 2012-11-05 22:36:04 edniner Exp $
6 ///
8 
12 #include "RecoBase/Cluster.h"
13 #include "RecoBase/Prong.h"
14 #include "RecoBase/RecoHit.h"
15 #include "Utilities/func/MathUtil.h"
16 
17 #include "TVector3.h"
18 
20 
21 #include <iostream>
22 #include <cmath>
23 
24 namespace fuzz{
25 
27  {
28  }
29 
30  //......................................................................
31 
33  {
34  }
35 
36  //......................................................................
37 
38  /// Load in candidate 2D prongs to prepare for matching
39  void ViewMatchAlg::LoadProngs(std::vector<rb::Prong> prongs)
40  {
41  fXProngs.clear();
42  fYProngs.clear();
43  fMatchedProngs.clear();
44  fUnmatchedProngs.clear();
45 
46  for (size_t i=0; i< prongs.size(); i++){
47  if (prongs[i].View() == geo::kX) fXProngs.push_back(prongs[i]);
48  if (prongs[i].View() == geo::kY) fYProngs.push_back(prongs[i]);
49  }
50 
51  }
52 
53  //......................................................................
54 
55  /// Function to do the work and execute the matching
57  {
59 
60  fMatchedProngs.clear();
61  fUnmatchedProngs.clear();
62 
63  //clear containers to hold summary debuggin information for matching ntuple
64  xzid.clear();
65  yzid.clear();
66  xzs.clear();
67  yzs.clear();
68  xzpe.clear();
69  yzpe.clear();
70  ksscore.clear();
71  ksshift.clear();
72 
73  //size matrix to hold matching score for all pairs of tracks
74  double score[fXProngs.size()][fYProngs.size()];
75 
76  //loop over xz prongs to create matching pairs
77  for (size_t i=0; i<fXProngs.size(); i++){
78  int ixPlane = fXProngs[i].MinPlane();
79  int fxPlane = fXProngs[i].MaxPlane();
80  double dirX = fXProngs[i].Dir().Z();
81 
82  //if we have a 1 plane track, alter the direction for maximum tilt, defense against near vertical tracks
83  if ((ixPlane == fxPlane)||(acos(fabs(fXProngs[i].Dir().X()))<0.01)) {
84  if (dirX > 0.0){
85  if (fXProngs[i].Dir().X() > 0.0){
86  fXProngs[i].SetDir(fGeo->Plane(0)->Cell(0)->HalfW()+fXProngs[i].MaxX()-fXProngs[i].Start().X(),
87  fGeo->Plane(0)->Cell(0)->HalfW()+fXProngs[i].MaxZ()-fXProngs[i].Start().Z());
88  }
89  else if (fXProngs[i].Dir().X() < 0.0){
90  fXProngs[i].SetDir(-fGeo->Plane(0)->Cell(0)->HalfW()+fXProngs[i].MinX()-fXProngs[i].Start().X(),
91  fGeo->Plane(0)->Cell(0)->HalfW()+fXProngs[i].MaxZ()-fXProngs[i].Start().Z());
92  }
93  }
94  else if (dirX < 0.0){
95  if (fXProngs[i].Dir().X() > 0.0){
96  fXProngs[i].SetDir(fGeo->Plane(0)->Cell(0)->HalfW()+fXProngs[i].MaxX()-fXProngs[i].Start().X(),
97  -fGeo->Plane(0)->Cell(0)->HalfW()+fXProngs[i].MinZ()-fXProngs[i].Start().Z());
98  }
99  else if (fXProngs[i].Dir().X() < 0.0){
100  fXProngs[i].SetDir(-fGeo->Plane(0)->Cell(0)->HalfW()+fXProngs[i].MinX()-fXProngs[i].Start().X(),
101  -fGeo->Plane(0)->Cell(0)->HalfW()+fXProngs[i].MinZ()-fXProngs[i].Start().Z());
102  }
103  }
104  }
105 
106  //inner loop over yz prongs to make pairs
107  for (size_t j=0; j<fYProngs.size(); j++){
108  int iyPlane = fYProngs[j].MinPlane();
109  int fyPlane = fYProngs[j].MaxPlane();
110  double dirY = fYProngs[j].Dir().Z();
111 
112  //if we have a 1 plane track, alter the direction for maximum tilt, defense against near vertical tracks
113  if (((iyPlane == fyPlane)||(acos(fabs(fYProngs[j].Dir().Y()))<0.01))&&(i==0)) {
114  if (dirY > 0.0){
115  if (fYProngs[j].Dir().Y() > 0.0){
116  fYProngs[j].SetDir(fGeo->Plane(0)->Cell(0)->HalfW()+fYProngs[j].MaxY()-fYProngs[j].Start().Y(),
117  fGeo->Plane(0)->Cell(0)->HalfW()+fYProngs[j].MaxZ()-fYProngs[j].Start().Z());
118  }
119  else if (fYProngs[j].Dir().Y() < 0.0){
120  fYProngs[j].SetDir(-fGeo->Plane(0)->Cell(0)->HalfW()+fYProngs[j].MinY()-fYProngs[j].Start().Y(),
121  fGeo->Plane(0)->Cell(0)->HalfW()+fYProngs[j].MaxZ()-fYProngs[j].Start().Z());
122  }
123  }
124  else if (dirY < 0.0){
125  if (fYProngs[j].Dir().Y() > 0.0){
126  fYProngs[j].SetDir(fGeo->Plane(0)->Cell(0)->HalfW()+fYProngs[j].MaxY()-fYProngs[j].Start().Y(),
127  -fGeo->Plane(0)->Cell(0)->HalfW()+fYProngs[j].MinZ()-fYProngs[j].Start().Z());
128  }
129  else if (fYProngs[j].Dir().Y() < 0.0){
130  fYProngs[j].SetDir(-fGeo->Plane(0)->Cell(0)->HalfW()+fYProngs[j].MinY()-fYProngs[j].Start().Y(),
131  -fGeo->Plane(0)->Cell(0)->HalfW()+fYProngs[j].MinZ()-fYProngs[j].Start().Z());
132  }
133  }
134  }
135 
136  score[i][j] = -1.0;
137 
138  //make sure candidate prong in both views are either forwards are backwards
139  if (dirX*dirY < 0.0 ) continue;
140 
141  //create 3d test prong to compute score
142  TVector3 newDir;
143  TVector3 vert;
144  rb::Prong testprong = this->MakeTestProng(fXProngs[i], fYProngs[j], &newDir, &vert);
145 
146  //clamp vertex to detector volume to find which view we are in for shift
147  double vx = vert.X();
148  double vy = vert.Y();
149  double vz = vert.Z();
150  if (vert.X() < -fGeo->DetHalfWidth()) vx = -fGeo->DetHalfWidth()+15.0;
151  if (vert.X() > fGeo->DetHalfWidth()) vx = fGeo->DetHalfWidth()-15.0;
152  if (vert.Y() < -fGeo->DetHalfHeight()) vy = -fGeo->DetHalfHeight()+15.0;
153  if (vert.Y() > fGeo->DetHalfHeight()) vy = fGeo->DetHalfHeight()-15.0;
154  if (vert.Z() < 0.0) vz = 2.0;
155  if (vert.Z() > fGeo->DetLength()) vz = fGeo->DetLength()-2.0;
156  //get plane position of vertex, roam in cell if vertex was in the plastic
157  geo::CellUniqueId cid = fGeo->CellId(vx, vy, vz, 1.5, 1.5, 6., 1.0);
158  geo::View_t vertview = geo::kX;
159  unsigned int vertplane = 999;
160  if (cid != 0){
161  vertplane = fGeo->getPlaneID(cid);
162  vertview = (fGeo->Plane(vertplane))->View();
163  };
164 
165  //calculate energy profiles in each view
166  double xgev = 0.0;
167  double ygev = 0.0;
168  std::vector<std::pair<double, double> > xlist = this->CalcEnergyProfile(testprong, newDir, vert, &xgev, geo::kX);
169  std::vector<std::pair<double, double> > ylist = this->CalcEnergyProfile(testprong, newDir, vert, &ygev, geo::kY);
170 
171  //if the vertex is in the x view, shift all points on that track the equivalent of one cell forward or backward,
172  //based on the slope of the track.
173  double shift = 0.0;
174  double xyz[3];
175  double shiftz = vz;
176  //shift only if the vertex plane was defined
177  if ((vertview == geo::kX)&&(vertplane != 999)){
178  //if track in positive direction shift one plane up
179  if (dirX > 0.0){
180  unsigned int nplane = fGeo->NextPlaneOtherView(vertplane,1);
181  //make sure next plane is defined
182  if (nplane != 999000){
183  fGeo->Plane(nplane)->Cell(0)->GetCenter(xyz);
184  shiftz = xyz[2];
185  }
186  double shiftx = (testprong.Dir().X()/testprong.Dir().Z())*(shiftz-vz) + vx;
187  double shifty = (testprong.Dir().Y()/testprong.Dir().Z())*(shiftz-vz) + vy;
188  shift = util::pythag(vx-shiftx,vy-shifty,vz-shiftz);
189  }
190  //shift one plane down for backwards track
191  else{
192  unsigned int nplane = fGeo->NextPlaneOtherView(vertplane,-1);
193  if (nplane != 999000){
194  fGeo->Plane(nplane)->Cell(0)->GetCenter(xyz);
195  shiftz = xyz[2];
196  }
197  double shiftx = (testprong.Dir().X()/testprong.Dir().Z())*(shiftz-vz) + vx;
198  double shifty = (testprong.Dir().Y()/testprong.Dir().Z())*(shiftz-vz) + vy;
199  shift = -util::pythag(vx-shiftx,vy-shifty,vz-shiftz);
200  }
201  }
202  //else if the plane was in the YZ view, shift xz view accordingly
203  if (vertview == geo::kY){
204  if (dirX > 0.0){
205  unsigned int nplane = fGeo->NextPlaneOtherView(vertplane,1);
206  if (nplane != 999000){
207  fGeo->Plane(nplane)->Cell(0)->GetCenter(xyz);
208  shiftz = xyz[2];
209  }
210  double shiftx = (testprong.Dir().X()/testprong.Dir().Z())*(shiftz-vz) + vx;
211  double shifty = (testprong.Dir().Y()/testprong.Dir().Z())*(shiftz-vz) + vy;
212  shift = -util::pythag(vx-shiftx,vy-shifty,vz-shiftz);
213  }
214  else{
215  unsigned int nplane = fGeo->NextPlaneOtherView(vertplane,-1);
216  if (nplane != 999000){
217  fGeo->Plane(nplane)->Cell(0)->GetCenter(xyz);
218  shiftz = xyz[2];
219  }
220  double shiftx = (testprong.Dir().X()/testprong.Dir().Z())*(shiftz-vz) + vx;
221  double shifty = (testprong.Dir().Y()/testprong.Dir().Z())*(shiftz-vz) + vy;
222  shift = util::pythag(vx-shiftx,vy-shifty,vz-shiftz);
223  }
224  }
225 
226  //if tracks from each view overlap by at least one plane, attempt a match
227  if(iyPlane <= fxPlane && fyPlane >= ixPlane){
228 
229  //perform kuiper test on candidate
230  double kScore = KuiperTest(xlist, ylist, shift, xgev, ygev);
231 
232  /////////////////////////
233  //block to collect summary information to monitor view matching process
234  //xzid.push_back(i);
235  //yzid.push_back(j);
236  //ksscore.push_back(kScore);
237  //std::vector<double> xstemp;
238  //std::vector<double> ystemp;
239  //std::vector<double> xpetemp;
240  //std::vector<double> ypetemp;
241  //for (unsigned int ix=0; ix<xlist.size(); ++ix){
242  // xstemp.push_back(xlist[ix].first);
243  // xpetemp.push_back(xlist[ix].second);
244  //}
245  //for (unsigned int iy=0; iy<ylist.size(); ++iy){
246  // ystemp.push_back(ylist[iy].first);
247  // ypetemp.push_back(ylist[iy].second);
248  //}
249  //xzs.push_back(xstemp);
250  //yzs.push_back(ystemp);
251  //xzpe.push_back(xpetemp);
252  //yzpe.push_back(ypetemp);
253  //////////////////////////////
254 
255  //save score as long as the tracks overlap, or if the track is only 1 plane in one view,
256  //then in the other view that track can only occupy the plane on either side
257  if((fYProngs[j].ExtentPlane()>1) && (fXProngs[i].ExtentPlane()>1)) score[i][j] = kScore;
258  else if((fYProngs[j].ExtentPlane() == 1) && (ixPlane == (iyPlane-1)) && (fxPlane == (iyPlane+1))) score[i][j] = kScore;
259  else if((fXProngs[i].ExtentPlane() == 1) && (iyPlane == (ixPlane-1)) && (fyPlane == (ixPlane+1))) score[i][j] = kScore;
260  }
261  //case if both tracks are in a single plane
262  else if((ixPlane == fxPlane) && (iyPlane == fyPlane) && (dirX*dirY > 0.0)){
263  // if x and y prongs appear in only one plane, allow match only if it is adjacent plane
264  if(ixPlane == (iyPlane-1) || ixPlane == (iyPlane+1)){
265  //perform test
266  double kScore = KuiperTest(xlist, ylist, shift, xgev, ygev);
267 
268  //////////////////////////////////////////
269  //block to collect summary information to monitor view matching process
270  //xzid.push_back(i);
271  //yzid.push_back(j);
272  //ksscore.push_back(kScore);
273  //std::vector<double> xstemp;
274  //std::vector<double> ystemp;
275  //std::vector<double> xpetemp;
276  //std::vector<double> ypetemp;
277  //for (unsigned int ix=0; ix<xlist.size(); ++ix){
278  // xstemp.push_back(xlist[ix].first);
279  // xpetemp.push_back(xlist[ix].second);
280  //}
281  //for (unsigned int iy=0; iy<ylist.size(); ++iy){
282  // ystemp.push_back(ylist[iy].first);
283  // ypetemp.push_back(ylist[iy].second);
284  //}
285  //xzs.push_back(xstemp);
286  //yzs.push_back(ystemp);
287  //xzpe.push_back(xpetemp);
288  //yzpe.push_back(ypetemp);
289  ///////////////////////////////////////
290 
291  score[i][j] = kScore;
292  }
293 
294  }//end case for single plane tracks
295  } // end yprong loop
296  } //end xprong loop
297 
298  // find the best score of all the possible matches.
299  // Now get the best matched prongs out of the group based on lowest score (not zero)
300  bool doneMatching = false;
301  std::vector<int> xMatched(fXProngs.size());
302  for(unsigned int i = 0; i<xMatched.size();++i){
303  xMatched[i] = 0;
304  }
305  std::vector<int> yMatched(fYProngs.size());
306  for(unsigned int i = 0; i<yMatched.size();++i){
307  yMatched[i] = 0;
308  }
309  while(!doneMatching){
310  int xbest = -1;
311  int ybest = -1;
312  double bestscore = -1;
313  // loop over all possible matched prongs to find the best one
314  for(unsigned int ix = 0; ix<fXProngs.size();++ix){
315  for(unsigned int iy = 0; iy<fYProngs.size();++iy){
316  if((score[ix][iy] <= bestscore || bestscore == -1) && (score[ix][iy] >= 0.0)){
317  xbest = ix;
318  ybest = iy;
319  bestscore = score[ix][iy];
320  }
321  }
322  }
323 
324  // store the best matched prong and reset things for the next best prong
325  if(bestscore >= 0){
326  // make a 3d prong out of the x and y prongs
327  rb::Prong newprong;
328 
329  // add all of the hits to the prong.
330  for(unsigned int i = 0; i<fXProngs[xbest].NCell();++i){
331  newprong.Add(fXProngs[xbest].Cell(i), fXProngs[xbest].Weight(i));
332  }
333  for(unsigned int i = 0; i<fYProngs[ybest].NCell();++i){
334  newprong.Add(fYProngs[ybest].Cell(i), fYProngs[ybest].Weight(i));
335  }
336 
337  newprong.SetID(fXProngs[xbest].ID());
338  TVector3 S;
339  S[0] = fXProngs[xbest].Start().X();
340  S[1] = fYProngs[ybest].Start().Y();
341  S[2] = fXProngs[xbest].Start().Z();
342  newprong.SetStart(S);
343  double tanThetax = tan(acos(fXProngs[xbest].Dir().X()));
344  double tanThetay = tan(acos(fYProngs[ybest].Dir().Y()));
345  double cosThetax = 1/(tanThetax*sqrt(1+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
346  double cosThetay = 1/(tanThetay*sqrt(1+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
347  double newz = sqrt(1-cosThetax*cosThetax-cosThetay*cosThetay);
348  if ((fXProngs[xbest].Dir().Z() < 0.0) && (fYProngs[ybest].Dir().Z() < 0.0)) newz = -newz;
349  TVector3 newDir(cosThetax,cosThetay,newz);
350  newprong.SetDir(newDir);
351 
352  fMatchedProngs.push_back(newprong);
353  xMatched[xbest] = 1;
354  yMatched[ybest] = 1;
355 
356  // set all the prongs that are matched to either this x or y prong to have a null match score.
357  for(unsigned int i = 0; i<fYProngs.size();++i){
358  score[xbest][i] = -1.0;
359  }
360  for(unsigned int i = 0; i<fXProngs.size();++i){
361  score[i][ybest] = -1.0;
362  }
363 
364  }
365  else{ doneMatching = true; }
366  }//end loop to make matches
367 
368  for (unsigned int i=0; i<xMatched.size(); i++){
369  if (xMatched[i] == 0) fUnmatchedProngs.push_back(fXProngs[i]);
370  }
371  for (unsigned int i=0; i<yMatched.size(); i++){
372  if (yMatched[i] == 0) fUnmatchedProngs.push_back(fYProngs[i]);
373  }
374 
375 
376  }
377 
378  //......................................................................
379 
380  /// Function to calculate the energy profile of the test prong for a given view
381  std::vector<std::pair<double,double> > ViewMatchAlg::CalcEnergyProfile(rb::Prong testprong, TVector3 dir, TVector3 start, double *totgev, geo::View_t v)
382  {
385 
386  std::vector<std::pair<double,double> > list;
387 
388  unsigned int ncells = testprong.NCell(v);
389  unsigned int ip, ic, n, m;
390  double w; // Distance of closest approach from read-out end of cell
391  double s; // Distance of closest approach from track start point
392 
393  //treat each cell as a 9 point space
394  double dx[3] = {-0.67, 0.0, 0.67};
395  double dy[3] = {-1.0, 0.0, 1.0};
396 
397  TVector3 pc; //point in cell
398  TVector3 pt; //point on track
399 
400  *totgev = 0.0;
401 
402  //loop over hits to make profile
403  for (unsigned int ih=0; ih<ncells; ++ih){
404  const rb::CellHit* ch = testprong.Cell(v, ih).get();
405  ip = ch->Plane();
406  ic = ch->Cell();
407 
408  rb::RecoHit rhit(cal->MakeRecoHit(*ch, testprong.W(ch)));
409 
410  if (!rhit.IsCalibrated()) continue;
411  //divide by 9.0 since we are breaking the cell into 9 pieces
412  double gev = (rhit.GeV()*testprong.Weight(v, ih))/9.0;
413 
414  //break each cell into 9 points to smear out energy
415  for (n=0; n<3; ++n) {
416  for (m=0; m<3; ++m) {
417  *totgev += gev;
418  s = -1.0; //the geometry function will not calculate s if it has a value of 0.0
419  fGeo->ClosestApproach(ip, ic, start, dir, dx[n], dy[m], &w, &s, &pc, &pt);
420  list.push_back(std::make_pair(s, gev));
421  }
422  }
423  }//end loop over cells
424 
425  //sort list from shortest to longest path
426  std::sort(list.begin(), list.end());
427 
428  return list;
429  }
430 
431  //......................................................................
432 
433  /// Function to combine a prong from each view into a test candidate
434  rb::Prong ViewMatchAlg::MakeTestProng(rb::Prong& xprong, rb::Prong& yprong, TVector3* newDir, TVector3* vert)
435  {
436  rb::Cluster clus;
437  for (unsigned int k=0; k<xprong.NCell(); k++) clus.Add(xprong.Cell(k), xprong.Weight(k));
438  for (unsigned int k=0; k<yprong.NCell(); k++) clus.Add(yprong.Cell(k), yprong.Weight(k));
439 
440  //make a direction for this track candidate in 3D
441  double tanThetax = tan(acos(xprong.Dir().X()));
442  double tanThetay = tan(acos(yprong.Dir().Y()));
443  double cosThetax = 1.0/(tanThetax*sqrt(1.0+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
444  double cosThetay = 1.0/(tanThetay*sqrt(1.0+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
445  double newz = sqrt(1.0-cosThetax*cosThetax-cosThetay*cosThetay);
446  if ((xprong.Dir().Z() < 0.0) && (yprong.Dir().Z() < 0.0)) newz = -newz;
447  newDir->SetXYZ(cosThetax,cosThetay,newz);
448 
449  //construct vertex candidate.
450  //At this point all input prongs were forced to start from the elastic arms vertex.
451  //So the z coordinate will be the same
452  vert->SetXYZ(xprong.Start().X(),yprong.Start().Y(),xprong.Start().Z());
453 
454  rb::Prong testprong(clus, *vert, *newDir);
455 
456  return testprong;
457 
458  }
459 
460  //.....................................................................................
461 
462  /// Function to perform the kuiper test to compare the similarity of the energy profiles between views
463  double ViewMatchAlg::KuiperTest(std::vector<std::pair<double,double>> xlist,
464  std::vector<std::pair<double,double>> ylist,
465  double shift,
466  double xgev,
467  double ygev)
468  {
469 
470  //if the track in either view did not have any calibrated energy, do not attempt a match
471  if ((xgev <= 0.0)||(ygev <= 0.0)){
472  return -1.0;
473  }
474 
475  //now perform ks test to determine maximum distance for these profiles
476 
477  double xsum = 0.0;//running sum of cumulative energy in xz view
478  double ysum = 0.0;//running sum of cumulative energy in yz view
479  double sx;//current point in xz track
480  double sy;//current point in yz track
481  double ksdistp;//positive difference between profiles
482  double ksdistm;//negative difference between profiles
483  double kuiper;//score metric
484  double ksmaxp = -0.00001;//largest positive difference between profiles
485  double ksmaxm = -0.00001;//largest negative difference between profiles
486 
487  unsigned int ix=0;
488  unsigned int iy=0;
489 
490  double kuiperBest = 2.0;
491  double shiftBest = -13.0;
492 
493  //loop to vary the shift around the estimate calculated earlier
494  //calculate the kuiper statistic for each shift value, keep the result with the lowest score
495  for (double shiftdist = shift-12.0; shiftdist<shift+12.01; shiftdist+=0.5){
496  ksmaxp = -0.00001;
497  ksmaxm = -0.00001;
498  ix = 0;
499  iy = 0;
500  xsum = 0.0;
501  ysum = 0.0;
502 
503 
504  while ((ix < xlist.size()) && (iy < ylist.size())) {
505  sx = xlist[ix].first + shiftdist;
506  sy = ylist[iy].first;
507  if (sx <= sy){
508  xsum += xlist[ix].second;
509  ix++;
510  }
511  else if (sy < sx){
512  ysum += ylist[iy].second;
513  iy++;
514  }
515 
516  ksdistp = ysum/ygev - xsum/xgev;
517  ksdistm = xsum/xgev - ysum/ygev;
518 
519  if (ksdistp > ksmaxp) ksmaxp = ksdistp;
520  if (ksdistm > ksmaxm) ksmaxm = ksdistm;
521 
522  }
523 
524  if (ksmaxp < 0.0) ksmaxp = 0.0;
525  if (ksmaxm < 0.0) ksmaxm = 0.0;
526  kuiper = ksmaxp + ksmaxm;
527 
528  if (kuiper < kuiperBest){
529  kuiperBest = kuiper;
530  shiftBest = shiftdist;
531  }
532  }
533 
534  ksshift.push_back(shiftBest);
535  return kuiperBest;
536  }
537 
538 
539 }//end hough namespace
540 
541 ////////////////////////////////////////////////////////////////////////
542 
543 
void LoadProngs(std::vector< rb::Prong > prongs)
Function to load the vector of 2d prongs to be matched.
std::vector< rb::Prong > fYProngs
Container for all the candidate YZ prongs.
Definition: ViewMatchAlg.h:74
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
std::vector< std::pair< double, double > > CalcEnergyProfile(rb::Prong testprong, TVector3 dir, TVector3 start, double *totpecorr, geo::View_t v)
Function to calculate the cumulative energy profile of a test prong.
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
#define S(x, n)
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
Definition: Prong.cxx:135
rb::Prong MakeTestProng(rb::Prong &xprong, rb::Prong &yprong, TVector3 *newdir, TVector3 *vert)
Function to combine a 2d prong from each view into a 3d test prong to compute a matching score for...
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Definition: CellHit.h:39
void Matching()
Function to perform the view matching.
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
T sqrt(T number)
Definition: d0nt_math.hpp:156
Vertical planes which measure X.
Definition: PlaneGeo.h:28
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
double HalfW() const
Definition: CellGeo.cxx:191
A collection of associated CellHits.
Definition: Cluster.h:47
double DetLength() const
TString ip
Definition: loadincs.C:5
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
ViewMatchAlg(const ViewMatchParams &params)
virtual TVector3 Start() const
Definition: Prong.h:73
T acos(T number)
Definition: d0nt_math.hpp:54
virtual void SetStart(TVector3 start)
Definition: Prong.cxx:75
std::vector< rb::Prong > fMatchedProngs
Container for all the 3D prongs produced after matching.
Definition: ViewMatchAlg.h:77
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
Float_t Y
Definition: plot.C:38
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
Float_t Z
Definition: plot.C:38
unsigned short Cell() const
Definition: CellHit.h:40
const XML_Char * s
Definition: expat.h:262
double dy[NP][NC]
double dx[NP][NC]
virtual void SetDir(TVector3 dir)
Definition: Prong.cxx:104
int getPlaneID(const CellUniqueId &id) const
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
const CellUniqueId CellId(const double &x, const double &y, const double &z, double dxds=0., double dyds=0., double dzds=1., double step=0.01) const
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
const double j
Definition: BetheBloch.cxx:29
Fuzzy k-Means prong-finding algorithm.
void SetID(int id)
Definition: Cluster.h:74
double DetHalfHeight() const
Definition: View.py:1
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
unsigned long long int CellUniqueId
Definition: CellUniqueId.h:15
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
std::vector< std::vector< double > > yzpe
Definition: ViewMatchAlg.h:91
float GeV() const
Definition: RecoHit.cxx:69
A Cluster with defined start position and direction.
Definition: Prong.h:19
double DetHalfWidth() const
T const * get() const
Definition: Ptr.h:321
std::vector< double > ksshift
Definition: ViewMatchAlg.h:87
TDirectory * dir
Definition: macro.C:5
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
double KuiperTest(std::vector< std::pair< double, double >> xlist, std::vector< std::pair< double, double >> ylist, double shift, double xpecorr, double ypecorr)
Function to perform the Kuiper test to compare the energy profiles of the track in each view...
double ClosestApproach(unsigned int ip, unsigned int ic, const TVector3 &vtx, const TVector3 &dir, double dx=0.0, double dy=0.0, double *w=0, double *s=0, TVector3 *point_in_cell=0, TVector3 *point_on_track=0) const
Find the distance of closest approach between a cell and a track.
std::vector< rb::Prong > fXProngs
Container for all the candidate XZ prongs.
Definition: ViewMatchAlg.h:71
int ncells
Definition: geom.C:124
std::vector< int > xzid
store information related to the view matching for a summary ntuple
Definition: ViewMatchAlg.h:84
std::vector< rb::Prong > fUnmatchedProngs
Container for all the unmatched 2D prongs left over.
Definition: ViewMatchAlg.h:80
T tan(T number)
Definition: d0nt_math.hpp:144
std::vector< double > ksscore
Definition: ViewMatchAlg.h:86
Float_t X
Definition: plot.C:38
Float_t w
Definition: plot.C:20
std::vector< std::vector< double > > xzpe
Definition: ViewMatchAlg.h:90
std::vector< std::vector< double > > yzs
Definition: ViewMatchAlg.h:89
std::vector< std::vector< double > > xzs
Definition: ViewMatchAlg.h:88
double Weight(unsigned int globalIdx) const
Weight assigned to the cell.
Definition: Cluster.cxx:209
std::vector< int > yzid
Definition: ViewMatchAlg.h:85
const unsigned int NextPlaneOtherView(unsigned int p, int d=+1) const