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