KalmanTrackMerge_module.cc
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////////
2 // \file KalmanTrackMerge_module.cc
3 // \brief A module for merging 2d tracks into 3d tracks
4 // \version $Id: KalmanTrackMerge_module.cc,v 1.8 2012-11-28 17:16:59 raddatz Exp $
5 // \author Nick Raddatz
6 ////////////////////////////////////////////////////////////////////////////
7 
8 #include <string>
9 #include <vector>
10 #include <cmath>
11 
12 // ROOT includes
13 #include "TMath.h"
14 #include "TVector3.h"
15 
16 
17 // NOvA includes
18 #include "Geometry/Geometry.h"
19 #include "GeometryObjects/Geo.h"
21 #include "RecoBase/Track.h"
22 #include "RecoBase/FilterList.h"
23 #include "Utilities/AssociationUtil.h"
24 #include "TrackFit/LocatedChan.h"
26 
27 
28 // Framework includes
32 
33 
34 //tracking algorithms
35 namespace trk {
37  public:
38  explicit KalmanTrackMerge(fhicl::ParameterSet const& pset);
39  virtual ~KalmanTrackMerge();
40  void produce(art::Event& evt);
41  void reconfigure(const fhicl::ParameterSet& p);
42  static bool sort_traj(const TVector3& a, const TVector3& b);
43 
44  private:
45  void MatchTracks(KalmanGeoHelper& kgeo,
46  std::vector<rb::Track> xtracks,
47  std::vector<rb::Track> ytracks,
48  std::vector<rb::Track>&matchedtracls,
49  std::vector<rb::Track>&unmatchedtracks);
50  bool CanJoinTracks(KalmanGeoHelper& kgeo, rb::Track tracka, rb::Track trackb);
52  double CalcMatchScore(KalmanGeoHelper& kgeo, rb::Track tracka, rb::Track trackb);
54  void ShiftInterpolationPoints(TVector3 endtraj, TVector3 nextraj,
55  unsigned int plane, rb::Track track,
56  TVector3 &lowztraj, TVector3 &highztraj);
57  std::string fTrackInput; ///< Input folder from track reco
58  std::string fSliceInput; ///< Input folder of slices
59  bool fBadChannels; ///< Account for bad channels in gaps ?
60  bool fObeyPreselection; ///< Check rb::IsFiltered?
61  art::ServiceHandle<geo::Geometry> fgeom; ///< Detector geometry
62  std::vector<double> avgXPos;
63  std::vector<std::vector<double> > avgXPosZRange;
64  std::vector<double> avgYPos;
65  std::vector<std::vector<double> > avgYPosZRange;
66  double avgX;
67  double avgY;
68 
69  };
70 }
71 
72 namespace trk{
73 
74  //......................................................................
76  {
77  reconfigure(pset);
78  produces< std::vector<rb::Track> >();
79  produces< art::Assns<rb::Track, rb::Cluster> >();
80  }
81 
82  //......................................................................
84  {
85  }
86 
87  //......................................................................
89  {
90  fTrackInput = pset.get< std::string >("TrackInput") ;
91  fSliceInput = pset.get< std::string >("SliceInput") ;
92  fBadChannels = pset.get< bool >("BadChannels") ;
93  fObeyPreselection = pset.get< bool >("ObeyPreselection");
94  }
95 
96 
97  //......................................................................
98  /// Reconstruction method for this module. Provides read-write access
99  /// to the event data.
100 
102  {
103  // Declare a container for Track objects to be stored in the art::event
104  std::unique_ptr< std::vector<rb::Track> > mergetrackcol(new std::vector<rb::Track>);
105  std::unique_ptr< art::Assns<rb::Track, rb::Cluster> > assns(new art::Assns<rb::Track, rb::Cluster>);
106 
108 
109  // define a vector holding slices
111  evt.getByLabel(fSliceInput,slicecol);
112 
113  art::PtrVector<rb::Cluster> slicelist;
114  for(unsigned int i = 0; i<slicecol->size();++i){
115  if(fObeyPreselection && rb::IsFiltered(evt,slicecol,i)){ continue; }
116  art::Ptr<rb::Cluster> clust(slicecol,i);
117  slicelist.push_back(clust);
118  }
119 
120  // loop over the slices and merge tracks
121  art::FindManyP<rb::Track> fmt(slicelist, evt, fTrackInput);
122 
123  for(unsigned int iSlice = 0; iSlice<slicelist.size(); ++iSlice){
124  // do nothing for noise slices
125  if(slicelist[iSlice]->IsNoise()){ continue; }
126  // Get all of the tracks in the slice
127  const std::vector< art::Ptr<rb::Track> > slTracks = fmt.at(iSlice);
128  // if no associated tracks do nothing
129  if(slTracks.empty()){ continue; }
130  avgXPosZRange.resize(0);
131  avgXPos.resize(0);
132  avgYPosZRange.resize(0);
133  avgYPos.resize(0);
134  avgX = slicelist[iSlice]->MeanX();
135  avgY = slicelist[iSlice]->MeanY();
136  std::vector<trk::LocatedChan> xcellchans;
137  double zxmin = 100000;
138  double zxmax = 0;
139  for(unsigned int nx = 0; nx < slicelist[iSlice]->NXCell(); ++nx){
140  art::Ptr<rb::CellHit> xhit = slicelist[iSlice]->XCell(nx);
141  LocatedChan xchan(xhit->Plane(),xhit->Cell());
142  xchan.SetHit(xhit);
143  double xyz[3];
144  const geo::CellGeo* c = fgeom->Plane(xchan.Plane())->Cell(xchan.Cell());
145  c->GetCenter(xyz);
146  xchan.SetCenter(xyz[0],xyz[1],xyz[2]);
147  xchan.SetHalfWidths(c->HalfW(),c->HalfL(),c->HalfD());
148  xcellchans.push_back(xchan);
149  if(xyz[2]>zxmax){ zxmax = xyz[2]; }
150  if(xyz[2]<zxmin){ zxmin = xyz[2]; }
151  }
152  // get a vector that holds the average x position for every 100 cm in z
153  trk::SortByPlane(xcellchans);
154  // how many break points in z
155  int nxzbreak = (int)floor((zxmax-zxmin)/100);
156  // decide where to break the z positions up
157  double dz = (zxmax-zxmin)/((double)nxzbreak);
158  std::vector<double> xzbreaks;
159  xzbreaks.push_back(zxmin);
160  for(int ibrk = 1; ibrk<=nxzbreak; ++ibrk){
161  xzbreaks.push_back(zxmin+((double)ibrk)*dz);
162  }
163  xzbreaks.push_back(zxmax);
164  for(unsigned int ibrk = 1; ibrk < xzbreaks.size(); ++ibrk){
165  double xtot = 0;
166  double ntot = 0;
167  for(unsigned int xchan = 0;xchan<xcellchans.size();++xchan){
168  double xyz[3];
169  xcellchans[xchan].GetCenter(xyz);
170  if(xyz[2] >= xzbreaks[ibrk-1] && xyz[2] <=xzbreaks[ibrk]){
171  xtot+=xyz[0];
172  ntot+=1.0;
173  }
174  }
175  std::vector<double> zrng;
176  zrng.push_back(xzbreaks[ibrk-1]);
177  zrng.push_back(xzbreaks[ibrk]);
178  avgXPosZRange.push_back(zrng);
179  if(ntot>0){ avgXPos.push_back(xtot/ntot); }
180  else{ avgXPos.push_back(avgX); }
181  }
182 
183 
184  // Get the y cell hits out of the slice
185  std::vector<trk::LocatedChan> ycellchans;
186  double zymin = 100000;
187  double zymax = 0;
188  for(unsigned int ny = 0; ny < slicelist[iSlice]->NYCell(); ++ny){
189  art::Ptr<rb::CellHit> yhit = slicelist[iSlice]->YCell(ny);
190  LocatedChan ychan(yhit->Plane(),yhit->Cell());
191  ychan.SetHit(yhit);
192  double xyz[3];
193  const geo::CellGeo* c = fgeom->Plane(ychan.Plane())->Cell(ychan.Cell());
194  c->GetCenter(xyz);
195  ychan.SetCenter(xyz[0],xyz[1],xyz[2]);
196  ychan.SetHalfWidths(c->HalfL(),c->HalfW(),c->HalfD());
197  ycellchans.push_back(ychan);
198  if(xyz[2]>zymax){ zymax = xyz[2]; }
199  if(xyz[2]<zymin){ zymin = xyz[2]; }
200  }
201 
202  // how many break points in yz
203  int nyzbreak = (int)floor((zymax-zymin)/100);
204  // decide where to break the z positions up
205  dz = (zymax-zymin)/((double)nyzbreak);
206  std::vector<double> yzbreaks;
207  yzbreaks.push_back(zymin);
208  for(int ibrk = 1; ibrk<=nyzbreak; ++ibrk){
209  yzbreaks.push_back(zymin+((double)ibrk)*dz);
210  }
211  yzbreaks.push_back(zymax);
212  for(unsigned int ibrk = 1; ibrk < yzbreaks.size(); ++ibrk){
213  double ytot = 0;
214  double ntot = 0;
215  for(unsigned int ychan = 0;ychan<ycellchans.size();++ychan){
216  double xyz[3];
217  ycellchans[ychan].GetCenter(xyz);
218  if(xyz[2] >= yzbreaks[ibrk-1] && xyz[2] <=yzbreaks[ibrk]){
219  ytot+=xyz[1];
220  ntot+=1.0;
221  }
222  }
223  std::vector<double> zrng;
224  zrng.push_back(yzbreaks[ibrk-1]);
225  zrng.push_back(yzbreaks[ibrk]);
226  avgYPosZRange.push_back(zrng);
227  if(ntot>0){ avgYPos.push_back(ytot/ntot); }
228  else{ avgYPos.push_back(avgY); }
229  }
230 
231  // seperate tracks by view
232  std::vector<rb::Track> xTracks;
233  std::vector<rb::Track> yTracks;
234  for(unsigned int j = 0;j<slTracks.size();++j){
235  if(slTracks[j]->View() == geo::kX){
236  xTracks.push_back(*slTracks[j]);
237  }
238  else if(slTracks[j]->View() == geo::kY){
239  yTracks.push_back(*slTracks[j]);
240  }
241  }
242 
243 
244  std::vector<rb::Track> matchedTracks;
245  std::vector<rb::Track> unmatchedTracks;
246  MatchTracks(kgeo, xTracks,yTracks,matchedTracks,unmatchedTracks);
247  for(unsigned int j = 0; j<matchedTracks.size();j++){
248  mergetrackcol->push_back(matchedTracks[j]);
249  util::CreateAssn(*this,evt,*mergetrackcol,slicelist[iSlice],*assns);
250  }
251  for(unsigned int j = 0; j<unmatchedTracks.size();j++){
252  mergetrackcol->push_back(unmatchedTracks[j]);
253  util::CreateAssn(*this,evt,*mergetrackcol,slicelist[iSlice],*assns);
254  }
255  }
256 
257 
258  evt.put(std::move(mergetrackcol));
259  evt.put(std::move(assns));
260 
261  }
262 
263  //......................................................................
264 
266  std::vector<rb::Track > xTracks, std::vector<rb::Track> yTracks,
267  std::vector<rb::Track> &matchedtracks, std::vector<rb::Track> &unmatchedtracks)
268  {
269 
270  std::vector<rb::Track> allXTracks;
271  std::vector<std::vector< int > > jointxtrkIDs;
272  // Do a first order joining of all tracks in each view
273  // First add all of the individual tracks
274  for(unsigned int ix = 0; ix<xTracks.size();++ix){
275  allXTracks.push_back(xTracks[ix]);
276  std::vector<int> xid;
277  xid.push_back(ix);
278  jointxtrkIDs.push_back(xid);
279  }
280  // now add all of the possible joint x tracks
281  unsigned int ixmaxold = 0;
282  unsigned int ixmax = allXTracks.size();
283  unsigned int ixmin = 0;
284  while(ixmaxold != ixmax && allXTracks.size()<100){
285  ixmaxold = ixmax;
286  for(unsigned int ix = 0; ix<ixmaxold; ++ix){
287  // get the list of tracks that contributed to this track
288  std::vector<int> ixids = jointxtrkIDs[ix];
289  std::vector<int> jntids = ixids;
290  for(unsigned int jx = ixmin; jx<ixmaxold; ++jx){
291  // get the list of tracks that contributed to this track
292  std::vector<int> jxids = jointxtrkIDs[jx];
293  // check if any of these ids overlap, if so don't try to join
294  bool idoverlap = false;
295  for(unsigned int i = 0; i<ixids.size(); ++i){
296  for(unsigned int j = 0; j<jxids.size(); ++j){
297  if(ixids[i] == jxids[j]){
298  idoverlap = true;
299  break;
300  }
301  }
302  if(idoverlap){ break; }
303  }
304  if(idoverlap){ continue; }
305  // check if this combination has been used yet
306  for(unsigned int j = 0; j<jxids.size(); ++j){
307  jntids.push_back(jxids[j]);
308  }
309  std::sort(jntids.begin(),jntids.end());
310  bool used = false;
311  for(unsigned int i = 0; i<jointxtrkIDs.size(); ++i){
312  if(jntids.size() != jointxtrkIDs[i].size()){ continue; }
313  bool same = true;
314  std::vector<int> jntidscheck = jointxtrkIDs[i];
315  std::sort(jntidscheck.begin(),jntidscheck.end());
316  for(unsigned int j = 0; j<jointxtrkIDs[i].size(); ++j){
317  if(jntids[j] != jntidscheck[j]){
318  same = false;
319  break;
320  }
321  }
322  if(same){
323  used = true;
324  break;
325  }
326  }
327  if(used){ continue; }
328  bool join = CanJoinTracks(kgeo,allXTracks[ix],allXTracks[jx]);
329  if(join){
330  rb::Track JointTrack = JoinTracks(kgeo,allXTracks[ix],allXTracks[jx]);
331  allXTracks.push_back(JointTrack);
332  std::vector<int> xjointid;
333  for(unsigned int i = 0; i<ixids.size(); ++i){
334  xjointid.push_back(ixids[i]);
335  }
336  for(unsigned int j = 0; j<jxids.size(); ++j){
337  xjointid.push_back(jxids[j]);
338  }
339  jointxtrkIDs.push_back(xjointid);
340  }
341  }
342  ixmax = allXTracks.size();
343  }
344  ixmax = allXTracks.size();
345  ixmin = ixmaxold;
346  }
347 
348  std::vector<rb::Track> allYTracks;
349  std::vector<std::vector< int > > jointytrkIDs;
350  // Do a first order joining of all tracks in each view
351  // First add all of the individual tracks
352  for(unsigned int iy = 0; iy<yTracks.size();++iy){
353  allYTracks.push_back(yTracks[iy]);
354  std::vector<int> yid;
355  yid.push_back(iy);
356  jointytrkIDs.push_back(yid);
357  }
358  // now add all of the possible joint y tracks
359  unsigned int iymaxold = 0;
360  unsigned int iymax = allYTracks.size();
361  unsigned int iymin = 0;
362  while(iymaxold != iymax && allYTracks.size()<100){
363  iymaxold = iymax;
364  for(unsigned int iy = 0; iy<iymaxold; ++iy){
365  // get the list of tracks that contributed to this track
366  std::vector<int> iyids = jointytrkIDs[iy];
367  std::vector<int> jntids = iyids;
368  for(unsigned int jy = iymin; jy<iymaxold; ++jy){
369  // get the list of tracks that contributed to this track
370  std::vector<int> jyids = jointytrkIDs[jy];
371  // check if any of these ids overlap, if so don't try to join
372  bool idoverlap = false;
373  for(unsigned int i = 0; i<iyids.size(); ++i){
374  for(unsigned int j = 0; j<jyids.size(); ++j){
375  if(iyids[i] == jyids[j]){
376  idoverlap = true;
377  break;
378  }
379  }
380  if(idoverlap){ break; }
381  }
382  if(idoverlap){ continue; }
383  // check if this combination has been used yet
384  for(unsigned int j = 0; j<jyids.size(); ++j){
385  jntids.push_back(jyids[j]);
386  }
387  std::sort(jntids.begin(),jntids.end());
388  bool used = false;
389  for(unsigned int i = 0; i<jointytrkIDs.size(); ++i){
390  if(jntids.size() != jointytrkIDs[i].size()){ continue; }
391  bool same = true;
392  std::vector<int> jntidscheck = jointytrkIDs[i];
393  std::sort(jntidscheck.begin(),jntidscheck.end());
394  for(unsigned int j = 0; j<jointytrkIDs[i].size(); ++j){
395  if(jntids[j] != jntidscheck[j]){
396  same = false;
397  break;
398  }
399  }
400  if(same){
401  used = true;
402  break;
403  }
404  }
405  if(used){ continue; }
406  bool join = CanJoinTracks(kgeo,allYTracks[iy],allYTracks[jy]);
407  if(join){
408  rb::Track JointTrack = JoinTracks(kgeo,allYTracks[iy],allYTracks[jy]);
409  allYTracks.push_back(JointTrack);
410  std::vector<int> yjointid;
411  for(unsigned int i = 0; i<iyids.size(); ++i){
412  yjointid.push_back(iyids[i]);
413  }
414  for(unsigned int j = 0; j<jyids.size(); ++j){
415  yjointid.push_back(jyids[j]);
416  }
417  jointytrkIDs.push_back(yjointid);
418  }
419  }
420  iymax = allYTracks.size();
421  }
422  iymax = allYTracks.size();
423  iymin = iymaxold - 1;
424  }
425 
426 
427  // create a vector that tells the score of the possible matched tracks, lowest score >0 is best.
428  double score[allXTracks.size()][allYTracks.size()];
429  for(unsigned int ix = 0; ix<allXTracks.size();++ix){
430  //find the begining and ending plane for each x track
431  for(unsigned int iy = 0; iy<allYTracks.size();++iy){
432  score[ix][iy] = CalcMatchScore(kgeo, allXTracks[ix],allYTracks[iy]);
433  }
434  }
435 
436  // find the best score of all the possible matches.
437  // Now get the best matched tracks out of the group based on lowest score (not zero)
438  bool doneMatching = false;
439  std::vector<int> xMatched(allXTracks.size());
440  for(unsigned int i = 0; i<xMatched.size();++i){
441  xMatched[i] = 0;
442  }
443  std::vector<int> yMatched(allYTracks.size());
444  for(unsigned int i = 0; i<yMatched.size();++i){
445  yMatched[i] = 0;
446  }
447  while(!doneMatching){
448  int xbest = -1;
449  int ybest = -1;
450  double bestscore = -1;
451  // loop over all possible matched tracks to find the best one
452  for(unsigned int ix = 0; ix<allXTracks.size();++ix){
453  for(unsigned int iy = 0; iy<allYTracks.size();++iy){
454  if((score[ix][iy] <= bestscore || bestscore == -1) && score[ix][iy] >= 0.0 && score[ix][iy] <= 5.0){
455  xbest = ix;
456  ybest = iy;
457  bestscore = score[ix][iy];
458  }
459  }
460  }
461 
462  // store the best matched track and reset things for the next best track
463  if(bestscore >= 0){
464  // make a 3d track out of the x and y tracks
465 
466  rb::Track newtrack = ViewMergeTracks(kgeo,allXTracks[xbest],allYTracks[ybest]);
467 
468  matchedtracks.push_back(newtrack);
469 
470  // loop over the ids of the x best track
471  for(unsigned int i = 0; i<jointxtrkIDs[xbest].size(); ++i){
472  int matchedtrkid = jointxtrkIDs[xbest][i];
473  // loop over all x tracks
474  for(unsigned int itrk = 0; itrk<allXTracks.size(); ++itrk){
475  // if this track has been matched, no need in checking it again
476  if(xMatched[itrk] == 0){
477  // loop over all ids that make up this track
478  bool usedtrk = false;
479  for(unsigned int itrkid = 0; itrkid<jointxtrkIDs[itrk].size(); ++itrkid){
480  if(jointxtrkIDs[itrk][itrkid] == matchedtrkid){
481  usedtrk = true;
482  break;
483  }
484  }
485  if(usedtrk){
486  // this track needs to be marked as used and all scores associated with it set to -1
487  xMatched[itrk] = 1;
488  for(unsigned int iytrk = 0; iytrk<allYTracks.size();++iytrk){
489  score[itrk][iytrk] = -1.0;
490  }
491  }
492  }
493  }
494  }
495 
496  // loop over the ids of the y best track
497  for(unsigned int i = 0; i<jointytrkIDs[ybest].size(); ++i){
498  int matchedtrkid = jointytrkIDs[ybest][i];
499  // loop over all x tracks
500  for(unsigned int itrk = 0; itrk<allYTracks.size(); ++itrk){
501  // if this track has been matched, no need in checking it again
502  if(yMatched[itrk] == 0){
503  // loop over all ids that make up this track
504  bool usedtrk = false;
505  for(unsigned int itrkid = 0; itrkid<jointytrkIDs[itrk].size(); ++itrkid){
506  if(jointytrkIDs[itrk][itrkid] == matchedtrkid){
507  usedtrk = true;
508  break;
509  }
510  }
511  if(usedtrk){
512  // this track needs to be marked as used and all scores associated with it set to -1
513  yMatched[itrk] = 1;
514  for(unsigned int ixtrk = 0; ixtrk<allXTracks.size();++ixtrk){
515  score[ixtrk][itrk] = -1.0;
516  }
517  }
518  }
519  }
520  }
521 
522  }
523  else{ doneMatching = true; }
524  }
525  // Now clean up any unmatched tracks
526  for(unsigned int i = 0; i<xTracks.size();++i){
527  if(xMatched[i] == 0){
528  // make a new track object that is the same as the current track object
529  rb::Cluster clust(geo::kX);
530  rb::Track track(clust,xTracks[i].Start().X(),xTracks[i].Start().Z(),
531  xTracks[i].Dir().X(),xTracks[i].Dir().Z());
532  for(unsigned int j = 1 ; j<xTracks[i].NTrajectoryPoints();++j){
533  track.AppendTrajectoryPoint(xTracks[i].TrajectoryPoint(j).X(),xTracks[i].TrajectoryPoint(j).Z());
534  }
535  for(unsigned int j = 0; j<xTracks[i].NCell();++j){
536  track.Add(xTracks[i].Cell(j));
537  }
538  track.SetID(xTracks[i].ID());
539  unmatchedtracks.push_back(track);
540  }
541  }
542  for(unsigned int i = 0; i<yTracks.size();++i){
543  if(yMatched[i] == 0){
544  // make a new track object that is the same as the current track object
545  rb::Cluster clust(geo::kY);
546  rb::Track track(clust,yTracks[i].Start().Y(),yTracks[i].Start().Z(),
547  yTracks[i].Dir().Y(),yTracks[i].Dir().Z());
548  for(unsigned int j = 1 ; j<yTracks[i].NTrajectoryPoints();++j){
549  track.AppendTrajectoryPoint(yTracks[i].TrajectoryPoint(j).Y(),yTracks[i].TrajectoryPoint(j).Z());
550  }
551  for(unsigned int j = 0; j<yTracks[i].NCell();++j){
552  track.Add(yTracks[i].Cell(j));
553  }
554  track.SetID(yTracks[i].ID());
555  unmatchedtracks.push_back(track);
556  }
557  }
558 
559  }
560 
561  ///////////////////////////////////////////////////////////////////////////
563  rb::Track trka, rb::Track trkb)
564  {
565 
566  // check if these tracks can be added together
567  // first check how many cells apart from closest ends
568  TVector3 start1 = trka.Start();
569  TVector3 stop1 = trka.Stop();
570  TVector3 start2 = trkb.Start();
571  TVector3 stop2 = trkb.Stop();
572  int view = 0;
573  if(trka.View() == geo::kY){ view = 1; }
574  // calculate the intersection point of these two tracks
575 
576  double stst = (start1-start2).Mag();
577  double stsp = (start1-stop2).Mag();
578  double spst = (stop1-start2).Mag();
579  double spsp = (stop1-stop2).Mag();
580  TVector3 joint1;
581  TVector3 joint2;
582  TVector3 end1;
583  TVector3 end2;
584 
585  if(trka.NTrajectoryPoints() < 2 || trkb.NTrajectoryPoints() < 2) return false;
586 
587  if(stst < stsp && stst < spst && stst< spsp){
588  joint1 = start1;
589  joint2 = start2;
590  end1 = trka.TrajectoryPoint(1);
591  end2 = trkb.TrajectoryPoint(1);
592  }
593  else if(stsp < stst && stsp < spst && stsp < spsp){
594  joint1 = start1;
595  joint2 = stop2;
596  end1 = trka.TrajectoryPoint(1);
597  end2 = trkb.TrajectoryPoint(trkb.NTrajectoryPoints()-2);
598  }
599  else if(spst< stst && spst< stsp && spst < spsp){
600  joint1 = stop1;
601  joint2 = start2;
602  end1 = trka.TrajectoryPoint(trka.NTrajectoryPoints()-2);
603  end2 = trkb.TrajectoryPoint(1);
604  }
605  else{
606  joint1 = stop1;
607  joint2 = stop2;
608  end1 = trka.TrajectoryPoint(trka.NTrajectoryPoints()-2);
609  end2 = trkb.TrajectoryPoint(trkb.NTrajectoryPoints()-2);
610  }
611  // Can't join if too large a gap between end points
612  double jointv1 = joint1.X();
613  double jointv2 = joint2.X();
614  double endv1 = end1.X();
615  double endv2 = end2.X();
616  if(view == 1){
617  jointv1 = joint1.Y();
618  jointv2 = joint2.Y();
619  endv1 = end1.Y();
620  endv2 = end2.Y();
621  }
622 
623  // first find the planes of the joints
624  int plane1;
625  int plane2;
626 
627  if(joint1.Z() == start1.Z()){plane1 = trka.MinPlane(); }
628  else{ plane1 = trka.MaxPlane(); }
629  if(joint2.Z() == start2.Z()){plane2 = trkb.MinPlane(); }
630  else{ plane2 = trkb.MaxPlane(); }
631 
632  // make simple cuts first
633  // 1) No back scatter in z
634  // no backscatter in z, need to be careful about the joining if hits are on same plane, might look backscattered when really not
635  // or might not look backscattered when they really are, treat this case seperate
636  if(plane1 == plane2){
637  // in this case use same joint z to determine backscatter
638  if( (joint1.Z()-end1.Z())*(end2.Z()-joint1.Z()) < 0 ){
639  return false;
640  }
641  }
642  else{
643  if( (joint1.Z()-end1.Z())*(joint2.Z()-joint1.Z()) < 0 ||
644  (joint2.Z()-joint1.Z())*(end2.Z()-joint2.Z()) < 0 ){
645  return false;
646  }
647  }
648 
649  if(trka.MinPlane() == trka.MaxPlane() && trka.MeanZ() < trkb.MeanZ()){ return false; }
650  else if(trkb.MinPlane() == trkb.MaxPlane() && trkb.MeanZ() < trka.MeanZ()){ return false; }
651 
652  // 2) No vertical scatter at low z (almost no backscatter)
653  // Also if one track is completely vertical only allow it to be linked if the vertical track is at high z
654  if(joint1.Z() == end1.Z() && joint1.Z() < end2.Z()){ return false; }
655  else if(joint2.Z() == end2.Z() && joint2.Z() < end1.Z()){ return false; }
656 
657  // Got rid of easy checks, now check gaps and gap probabilities
658  // First thing find all the cells on the track at the track ends to be joined
659  int celllow1 = 9999;
660  int cellhigh1 = 0;
661  for(unsigned int iCell = 0; iCell<trka.NCell(); ++iCell){
662  if(trka.Cell(iCell)->Plane() == plane1){
663  if(trka.Cell(iCell)->Cell()>cellhigh1){ cellhigh1 = trka.Cell(iCell)->Cell(); }
664  if(trka.Cell(iCell)->Cell()<celllow1){ celllow1 = trka.Cell(iCell)->Cell();}
665  }
666  }
667 
668  int celllow2 = 9999;
669  int cellhigh2 = 0;
670  for(unsigned int iCell = 0; iCell<trkb.NCell(); ++iCell){
671  if(trkb.Cell(iCell)->Plane() == plane2){
672  if(trkb.Cell(iCell)->Cell()>cellhigh2){ cellhigh2 = trkb.Cell(iCell)->Cell(); }
673  if(trkb.Cell(iCell)->Cell()<celllow2){ celllow2 = trkb.Cell(iCell)->Cell();}
674  }
675  }
676 
677  int cellDiff = 0;
678  // find the actual cell difference, if something strange happens where higher v has lower cell set cellDiff as 0
679  // if it overlaps the window in which the v occurs
680  if(jointv1 > jointv2){
681  // cell diff should be celllow1-cellhigh2
682  if(celllow1 <= cellhigh2){ cellDiff = celllow1- cellhigh2; }
683  else if(cellhigh2 >= celllow1 && cellhigh2 <= cellhigh1){ cellDiff = 0; }
684  else{ cellDiff = -(cellhigh2- celllow1); }
685  }
686  else{
687  // cell diff should be celllow1-cellhigh2
688  if(cellhigh1 <= celllow2){ cellDiff = celllow2-cellhigh1; }
689  else if(cellhigh1 >= celllow2 && cellhigh1 <= cellhigh2){ cellDiff = 0;}
690  else{ cellDiff = -(cellhigh1-celllow2); }
691  }
692  // #3 if the two tracks are to be joined at the same plane, make sure no gap exists between hits
693  if(plane1 == plane2 && abs(cellDiff) > 1){ return false; }
694 
695  // #4 Make sure there is not a gap in wrong direction to cause double scatters
696  TVector3 diff1 = joint1-end1;
697  TVector3 jointdiff = joint2-joint1;
698  TVector3 diff2 = end2-joint2;
699  double cos1 = diff1.Dot(jointdiff)/(diff1.Mag()*jointdiff.Mag());
700  double cos2 = diff2.Dot(jointdiff)/(diff2.Mag()*jointdiff.Mag());
701  double slope1 = diff1.X()/diff1.Z();
702  double slope2 = diff2.X()/diff2.Z();
703  double slopediff = jointdiff.X()/jointdiff.Z();
704  if(view == 1){
705  slope1 = diff1.Y()/diff1.Z();
706  slope2 = diff2.Y()/diff2.Z();
707  slopediff = jointdiff.Y()/jointdiff.Z();
708  }
709  // This is scattering cosines are correct unless it is backscattered , then it is off by 90 degrees
710  double rotAngle1 = TMath::ATan(slope1);
711  double rotAngle2 = TMath::ATan(slopediff);
712  double joint2zrot = TMath::Cos(rotAngle1)*(joint2.Z() - joint1.Z()) - TMath::Sin(rotAngle1)*(jointv2 - jointv1);
713  if(joint2zrot < 0){ cos1 = -cos1; }
714  double end2zrot = TMath::Cos(rotAngle2)*(end2.Z() - joint2.Z()) - TMath::Sin(rotAngle2)*(endv2 - jointv2);
715  if(end2zrot < 0){ cos2 = -cos2; }
716 
717 
718  // check that the intersection points happen within a box that would make a hard scatter a single scatter
719  double vint, zint;
720  geo::LineIntersection(endv1,end1.Z(),jointv1,joint1.Z(),endv2,end2.Z(),jointv2,joint2.Z(),vint,zint);
721  // if parallel lines don't join unless one cell apart (plastic), these would have been made one track by 2D tracking if possible
722  double dx = jointv1-endv1;
723  double dy = joint1.Z()-end1.Z();
724  double dX = jointv2-endv2;
725  double dY = joint2.Z()-end2.Z();
726  double denom = dX*dy-dY*dx;
727  if(denom == 0 && cellDiff>1){ return false; }
728  double xyzpl1[3];
729  fgeom->Plane(plane1)->Cell(celllow1)->GetCenter(xyzpl1);
730  double dz1 = fgeom->Plane(plane1)->Cell(celllow1)->HalfD();
731  double dv1 = fgeom->Plane(plane1)->Cell(celllow1)->HalfW();
732  double zlow = xyzpl1[2]-dz1;
733  double vlow = xyzpl1[view]-dv1;
734  double xyzpl2[3];
735  fgeom->Plane(plane2)->Cell(cellhigh2)->GetCenter(xyzpl2);
736  double dz2 = fgeom->Plane(plane2)->Cell(cellhigh2)->HalfD();
737  double dv2 = fgeom->Plane(plane2)->Cell(cellhigh2)->HalfW();
738  double zhigh = xyzpl2[2]+dz2;
739  double vhigh = xyzpl2[view]+dv2;
740 
741  if(joint2.Z() < joint1.Z()){
742  zlow = xyzpl2[2]-dz2;
743  zhigh = xyzpl1[2]+dz1;
744  }
745  if((slope1*slope2) > 0 && plane1 != plane2 && abs(cellDiff)>1){
746  if(jointv2 < jointv1){
747  fgeom->Plane(plane2)->Cell(celllow2)->GetCenter(xyzpl2);
748  fgeom->Plane(plane1)->Cell(cellhigh1)->GetCenter(xyzpl1);
749  vlow = xyzpl2[view]-dv2;
750  vhigh = xyzpl1[view]+dv1;
751  }
752  if(vint < vlow || vint > vhigh || zint < zlow || zint > zhigh){ return false; }
753  }
754  else if((slope1*slope2) <= 0 && plane1 != plane2 && abs(cellDiff)>1){
755  // if tracks have opposite slopes can only check that z falls between z1 and z2
756  // check that v location is in the detector
757  if(view == 0){
758  vlow = -fgeom->DetHalfWidth();
759  vhigh = fgeom->DetHalfWidth();
760  }
761  else{
762  vlow = -fgeom->DetHalfHeight();
763  vhigh = fgeom->DetHalfHeight();
764  }
765  if(zint < zlow || zint > zhigh || vint < vlow || vint > vhigh){ return false; }
766  }
767  else if(plane1 == plane2){
768  // give more cushion here to match the uncertainty in the case wher plane1 != plane2
769  if(zint < zlow-dz1 || zint > zhigh+dz1){ return false; }
770  }
771 
772  // now get the missed cells between the potential link
773  int numMissed = 0;
774  std::vector<geo::OfflineChan> hitsOnline;
775  std::vector<double> distance;
776  if(plane1 != plane2){
777  // need to use intersection point to get true number of cells missed
778  int numMissed2 = 0;
779  std::vector<geo::OfflineChan> hitsOnline2;
780  std::vector<double> distance2;
781  if(view == 0){
782  numMissed = kgeo.CountMissedCellsOnLine(geo::kX, joint1.X(), joint1.Z(), vint, zint, hitsOnline, distance);
783  numMissed2 = kgeo.CountMissedCellsOnLine(geo::kX, vint, zint, joint2.X(), joint2.Z(), hitsOnline2, distance2);
784  }
785  else{
786  numMissed = kgeo.CountMissedCellsOnLine(geo::kY, joint1.Y(), joint1.Z(), vint, zint, hitsOnline, distance);
787  numMissed2 = kgeo.CountMissedCellsOnLine(geo::kY, vint, zint, joint2.Y(), joint2.Z(), hitsOnline2, distance2);
788  }
789  // now add misses on the second bit of the scatter to the first if it doesn't exist already
790  for(int imiss2 = 0; imiss2 < numMissed2; ++imiss2){
791  ++numMissed;
792  hitsOnline.push_back(hitsOnline2[imiss2]);
793  distance.push_back(distance2[imiss2]);
794  }
795  // now loop over the hitsOnline to see if any fall into the cells at the end of either track
796  std::vector<int> missesToRemove;
797  for(unsigned int imiss = 0; imiss<hitsOnline.size();++imiss){
798  if(hitsOnline[imiss].Plane() == plane1 && hitsOnline[imiss].Cell() >= celllow1 && hitsOnline[imiss].Cell() <= cellhigh1){
799  missesToRemove.push_back(imiss);
800  }
801  else if(hitsOnline[imiss].Plane() == plane2 && hitsOnline[imiss].Cell() >= celllow2 && hitsOnline[imiss].Cell() <= cellhigh2){
802  missesToRemove.push_back(imiss);
803  }
804  }
805  // careful to not remove wrong ones by erasing from highest index in vector to lowest
806  reverse(missesToRemove.begin(),missesToRemove.end());
807  for(unsigned int imiss = 0; imiss < missesToRemove.size(); ++imiss){
808  --numMissed;
809  hitsOnline.erase(hitsOnline.begin()+missesToRemove[imiss]);
810  distance.erase(distance.begin()+missesToRemove[imiss]);
811  }
812  }
813  else{
814  // CountMissedCellsOnLine not quite right when both z values on the same plane
815  unsigned short lowestCell = cellhigh1;
816  unsigned short highestCell = celllow2;
817  if(jointv1 > jointv2){
818  lowestCell = cellhigh2;
819  highestCell = celllow1;
820  }
821  for(unsigned short imiss = lowestCell+1; imiss<highestCell;++imiss){
822  geo::OfflineChan miss(plane1,imiss);
823  hitsOnline.push_back(miss);
824  // get the width of cell to travel through
825  distance.push_back(2*fgeom->Plane(plane1)->Cell(imiss)->HalfW());
826  }
827  }
828  // got the missed cells, now get the probability of missing these cells
829  double probMiss = 1;
830  for(unsigned int i = 0; i < distance.size(); ++i){
831 
832  // find the approximate w of this cell
833  double w = 0;
834  double zave = 0.5*joint1.Z()+0.5*joint2.Z();
835  if(trka.View() == geo::kX){
836  // loop over all of the breakpoints and find what w to use
837  bool foundrng = false;
838  if(zave < avgYPosZRange[0][0]){
839  // use the first range
840  w = avgYPos[0];
841  foundrng = true;
842  }
843  else if(zave > avgYPosZRange[avgYPosZRange.size()-1][1]){
844  // use the last range
845  w = avgYPos[avgYPos.size()-1];
846  foundrng = true;
847  }
848  else{
849  // loop over all the ranges and find the one that encloses this z
850  for(unsigned int i = 0; i < avgYPosZRange.size();++i){
851  if(zave >= avgYPosZRange[i][0] && i <= avgYPosZRange[i][1]){
852  w = avgYPos[i];
853  foundrng = true;
854  break;
855  }
856  }
857  }
858  if(!foundrng){ w = avgYPos[0]; }
859  }
860  else if(trka.View() == geo::kY){
861  // loop over all of the breakpoints and find what w to use
862  bool foundrng = false;
863  if(zave < avgXPosZRange[0][0]){
864  // use the first range
865  w = avgXPos[0];
866  foundrng = true;
867  }
868  else if(zave > avgXPosZRange[avgXPosZRange.size()-1][1]){
869  // use the last range
870  w = avgXPos[avgXPos.size()-1];
871  foundrng = true;
872  }
873  else{
874  // loop over all the ranges and find the one that encloses this z
875  for(unsigned int i = 0; i < avgXPosZRange.size();++i){
876  if(zave >= avgXPosZRange[i][0] && i<= avgXPosZRange[i][1]){
877  w = avgXPos[i];
878  foundrng = true;
879  break;
880  }
881  }
882  }
883  if(!foundrng){ w = avgXPos[0]; }
884  }
885  probMiss*=(kgeo.CalcMissFrac(distance[i],w,trka.View()));
886  }
887 
888  double minProb = 0.00001;
889  double minScatProb = 0.1;
890  // #5 if prob miss is greater than the gap prob, then this has to be a hard scatter
891  // if it is a true hard scatter, be careful about joining, scale gap probability up by scattering angles
892  // #6 or if the prob miss is less than the gap prob, then the only way to get this should be by going through plastic
893  // if we are in plastic, then the cell difference between the two ends should be zero or one only
894  if(probMiss > minProb){
895  // remove any remaining backscatter
896  if((cos1+cos2)<=0){ return false; }
897  if(probMiss < 2.0*minProb/(cos1+cos2)){ return false; }
898  if(probMiss < minScatProb){ return false; }
899  }
900  else{
901  // check that slopes are relatively flat ie less than 10 degrees to be consistent with traveling through large sections of plastic
902  if(abs(cellDiff) > 1){ return false; }
903  double trkcos1 = TMath::Cos(TMath::ATan(slope1));
904  double trkcos2 = TMath::Cos(TMath::ATan(slope2));
905  double cosmax = trkcos1;
906  double cosmin = trkcos2;
907  if(trkcos2 > trkcos1){
908  cosmax = trkcos2;
909  cosmin = trkcos1;
910  }
911  if(cosmin < 0.98 || cosmax < 0.995){ return false; }
912  }
913 
914  // made it through all checks retrun that they can be joined
915  return true;
916  }
917 
918  ///////////////////////////////////////////////////////////////////////
919 
921  rb::Track trka, rb::Track trkb)
922  {
923  // check that these two tracks belong to the same view
924  assert(trka.View() == trkb.View());
925  assert(trka.View() != geo::kXorY);
926  rb::Cluster JointCluster(trka.View(),trka.ID());
927  // add hits from trka and trkb
928  JointCluster.Add(trka.AllCells());
929  JointCluster.Add(trkb.AllCells());
930 
931  // find which ends to join together
932  TVector3 start1 = trka.Start();
933  TVector3 stop1 = trka.Stop();
934  TVector3 start2 = trkb.Start();
935  TVector3 stop2 = trkb.Stop();
936  double stst = (start1-start2).Mag();
937  double stsp = (start1-stop2).Mag();
938  double spst = (stop1-start2).Mag();
939  double spsp = (stop1-stop2).Mag();
940  int ntraja = trka.NTrajectoryPoints();
941  int ntrajb = trkb.NTrajectoryPoints();
942  // find the start/stop planes
943  std::vector<double> zends;
944  zends.push_back(start1.Z());
945  zends.push_back(stop1.Z());
946  zends.push_back(start2.Z());
947  zends.push_back(stop2.Z());
948  std::vector<TVector3> trkends;
949  trkends.push_back(start1);
950  trkends.push_back(stop1);
951  trkends.push_back(start2);
952  trkends.push_back(stop2);
953  // Find the planes that correspond to the start and stop z
954  int trkaminpl = trka.MinPlane();
955  int trkamaxpl = trka.MaxPlane();
956  int trkbminpl = trkb.MinPlane();
957  int trkbmaxpl = trkb.MaxPlane();
958  int zplanes[4] = {trkaminpl, trkamaxpl, trkbminpl, trkbmaxpl};
959 
960  for(int iz = 0; iz < 4; ++iz){
961  if(iz < 2){
962  int plane = kgeo.MatchTrajectoryToPlaneInView(trkends[iz],trka.View(),trkaminpl,trkamaxpl);
963  zplanes[iz] = plane;
964  }
965  else{
966  int plane = kgeo.MatchTrajectoryToPlaneInView(trkends[iz],trkb.View(),trkbminpl,trkbmaxpl);
967  zplanes[iz] = plane;
968  }
969  }
970 
971  // always joint the ends that are closest together
972  std::vector<TVector3> traja(ntraja);
973  for(int i = 0; i<ntraja; ++i){
974  traja[i] = trka.TrajectoryPoint(i);
975  }
976  std::vector<TVector3> trajb(ntrajb);
977  for(int i = 0; i<ntrajb; ++i){
978  trajb[i] = trkb.TrajectoryPoint(i);
979  }
980  std::vector<TVector3> trajpts;
981  bool sameplane = false;
982  // get all of the hits on the end planes
983  double plhitspetot[4] = {0, 0, 0, 0};
984  for(unsigned int i = 0; i<trka.NCell(); ++i){
985  if(trka.Cell(i)->Plane() == zplanes[0]){ plhitspetot[0]+=trka.Cell(i)->PE(); }
986  if(trka.Cell(i)->Plane() == zplanes[1]){ plhitspetot[1]+=trka.Cell(i)->PE(); }
987  }
988  for(unsigned int i = 0; i<trkb.NCell(); ++i){
989  if(trkb.Cell(i)->Plane() == zplanes[2]){ plhitspetot[2]+=trkb.Cell(i)->PE(); }
990  if(trkb.Cell(i)->Plane() == zplanes[3]){ plhitspetot[3]+=trkb.Cell(i)->PE(); }
991  }
992 
993  if(stst < stsp && stst < spst && stst< spsp){
994  // link at start points
995  if(zplanes[0] == zplanes[2]){ sameplane = true; }
996  // check which end point is lower in z
997  if(trajb[ntrajb-1].Z() < traja[ntraja-1].Z()){
998  // trkb first then trka, but start with stop of b
999  std::reverse(trajb.begin(),trajb.end());
1000  for(int i = 0; i<ntrajb-1;++i){
1001  trajpts.push_back(trajb[i]);
1002  }
1003  if(!sameplane){
1004  trajpts.push_back(trajb[ntrajb-1]);
1005  trajpts.push_back(traja[0]);
1006  }
1007  else{
1008  TVector3 mergetrj = (plhitspetot[0]*traja[0]+plhitspetot[2]*trajb[ntrajb-1])*
1009  (1.0/(plhitspetot[0]+plhitspetot[2]));
1010  if(!std::isnan(mergetrj.X()) && !std::isnan(mergetrj.Y()) && !std::isnan(mergetrj.Z())){ trajpts.push_back(mergetrj); }
1011  else{ trajpts.push_back(traja[0]); }
1012  }
1013  for(int i = 1; i<ntraja;++i){
1014  trajpts.push_back(traja[i]);
1015  }
1016  }
1017  else{
1018  // trka first then trkb, but start with stop of a
1019  std::reverse(traja.begin(),traja.end());
1020  for(int i = 0; i<ntraja-1;++i){
1021  trajpts.push_back(traja[i]);
1022  }
1023  if(!sameplane){
1024  trajpts.push_back(traja[ntraja-1]);
1025  trajpts.push_back(trajb[0]);
1026  }
1027  else{
1028  TVector3 mergetrj = (plhitspetot[0]*traja[ntraja-1]+plhitspetot[2]*trajb[0])*
1029  (1.0/(plhitspetot[0]+plhitspetot[2]));
1030  if(!std::isnan(mergetrj.X()) && !std::isnan(mergetrj.Y()) && !std::isnan(mergetrj.Z())){ trajpts.push_back(mergetrj); }
1031  else{ trajpts.push_back(traja[ntraja-1]); }
1032  }
1033  for(int i = 1; i<ntrajb;++i){
1034  trajpts.push_back(trajb[i]);
1035  }
1036  }
1037  }
1038  else if(stsp < stst && stsp < spst && stsp < spsp){
1039  // link at start in a stop in b
1040  if(zplanes[0] == zplanes[3]){ sameplane = true; }
1041  // check which end point is lower
1042  if(trajb[0].Z() < traja[ntraja-1].Z()){
1043  // trkb first then trka, start with start of b
1044  for(int i = 0; i<ntrajb-1;++i){
1045  trajpts.push_back(trajb[i]);
1046  }
1047  if(!sameplane){
1048  trajpts.push_back(trajb[ntrajb-1]);
1049  trajpts.push_back(traja[0]);
1050  }
1051  else{
1052  TVector3 mergetrj = (plhitspetot[0]*traja[0]+plhitspetot[3]*trajb[ntrajb-1])*
1053  (1.0/(plhitspetot[0]+plhitspetot[3]));
1054  if(!std::isnan(mergetrj.X()) && !std::isnan(mergetrj.Y()) && !std::isnan(mergetrj.Z())){ trajpts.push_back(mergetrj); }
1055  else{ trajpts.push_back(traja[0]); }
1056  }
1057  for(int i = 1; i<ntraja;++i){
1058  trajpts.push_back(traja[i]);
1059  }
1060  }
1061  else{
1062  // trka first then trkb, start with stop of a
1063  std::reverse(traja.begin(),traja.end());
1064  for(int i = 0; i<ntraja-1;++i){
1065  trajpts.push_back(traja[i]);
1066  }
1067  // add from stop in b
1068  std::reverse(trajb.begin(),trajb.end());
1069  if(!sameplane){
1070  trajpts.push_back(traja[ntraja-1]);
1071  trajpts.push_back(trajb[0]);
1072  }
1073  else{
1074  TVector3 mergetrj = (plhitspetot[0]*traja[ntraja-1]+plhitspetot[3]*trajb[0])*
1075  (1.0/(plhitspetot[0]+plhitspetot[3]));
1076  if(!std::isnan(mergetrj.X()) && !std::isnan(mergetrj.Y()) && !std::isnan(mergetrj.Z())){ trajpts.push_back(mergetrj); }
1077  else{ trajpts.push_back(traja[ntraja-1]); }
1078  }
1079  for(int i = 1; i<ntrajb;++i){
1080  trajpts.push_back(trajb[i]);
1081  }
1082  }
1083  }
1084  else if(spst< stst && spst< stsp && spst < spsp){
1085  // link at stop in a start in b
1086  if(zplanes[1] == zplanes[2]){ sameplane = true; }
1087  // check which end point is lower
1088  if(trajb[ntrajb-1].Z() < traja[0].Z()){
1089  // trkb first then trka, start with stop of b
1090  std::reverse(trajb.begin(),trajb.end());
1091  for(int i = 0; i<ntrajb-1;++i){
1092  trajpts.push_back(trajb[i]);
1093  }
1094  // add from stop in a
1095  std::reverse(traja.begin(),traja.end());
1096  if(!sameplane){
1097  trajpts.push_back(trajb[ntrajb-1]);
1098  trajpts.push_back(traja[0]);
1099  }
1100  else{
1101  TVector3 mergetrj = (plhitspetot[1]*traja[0]+plhitspetot[2]*trajb[ntrajb-1])*
1102  (1.0/(plhitspetot[1]+plhitspetot[2]));
1103  if(!std::isnan(mergetrj.X()) && !std::isnan(mergetrj.Y()) && !std::isnan(mergetrj.Z())){ trajpts.push_back(mergetrj); }
1104  else{ trajpts.push_back(traja[0]); }
1105  }
1106  for(int i = 1; i<ntraja;++i){
1107  trajpts.push_back(traja[i]);
1108  }
1109  }
1110  else{
1111  // trka first then trkb, start with start of a
1112  for(int i = 0; i<ntraja-1;++i){
1113  trajpts.push_back(traja[i]);
1114  }
1115  if(!sameplane){
1116  trajpts.push_back(traja[ntraja-1]);
1117  trajpts.push_back(trajb[0]);
1118  }
1119  else{
1120  TVector3 mergetrj = (plhitspetot[1]*traja[ntraja-1]+plhitspetot[2]*trajb[0])*
1121  (1.0/(plhitspetot[1]+plhitspetot[2]));
1122  if(!std::isnan(mergetrj.X()) && !std::isnan(mergetrj.Y()) && !std::isnan(mergetrj.Z())){ trajpts.push_back(mergetrj); }
1123  else{ trajpts.push_back(traja[ntraja-1]); }
1124  }
1125  // add from start in b
1126  for(int i = 1; i<ntrajb;++i){
1127  trajpts.push_back(trajb[i]);
1128  }
1129  }
1130  }
1131  else{
1132  // link at stop points
1133  if(zplanes[1] == zplanes[3]){ sameplane = true; }
1134  // check which end point is lower
1135  if(trajb[0].Z() < traja[0].Z()){
1136  // trkb first then trka, start with start of b
1137  for(int i = 0; i<ntrajb-1;++i){
1138  trajpts.push_back(trajb[i]);
1139  }
1140  // add from stop in a
1141  std::reverse(traja.begin(),traja.end());
1142  if(!sameplane){
1143  trajpts.push_back(trajb[ntrajb-1]);
1144  trajpts.push_back(traja[0]);
1145  }
1146  else{
1147  TVector3 mergetrj = (plhitspetot[1]*traja[0]+plhitspetot[3]*trajb[ntrajb-1])*
1148  (1.0/(plhitspetot[1]+plhitspetot[3]));
1149  if(!std::isnan(mergetrj.X()) && !std::isnan(mergetrj.Y()) && !std::isnan(mergetrj.Z())){ trajpts.push_back(mergetrj); }
1150  else{ trajpts.push_back(traja[0]); }
1151  }
1152  for(int i = 1; i<ntraja;++i){
1153  trajpts.push_back(traja[i]);
1154  }
1155  }
1156  else{
1157  // trka first then trkb, start with start of a
1158  for(int i = 0; i<ntraja-1;++i){
1159  trajpts.push_back(traja[i]);
1160  }
1161  // add from stop in b
1162  std::reverse(trajb.begin(),trajb.end());
1163  if(!sameplane){
1164  trajpts.push_back(traja[ntraja-1]);
1165  trajpts.push_back(trajb[0]);
1166  }
1167  else{
1168  TVector3 mergetrj = (plhitspetot[1]*traja[ntraja-1]+plhitspetot[3]*trajb[0])*
1169  (1.0/(plhitspetot[1]+plhitspetot[3]));
1170  if(!std::isnan(mergetrj.X()) && !std::isnan(mergetrj.Y()) && !std::isnan(mergetrj.Z())){ trajpts.push_back(mergetrj); }
1171  else{ trajpts.push_back(traja[ntraja-1]); }
1172  }
1173  for(int i = 1; i<ntrajb;++i){
1174  trajpts.push_back(trajb[i]);
1175  }
1176  }
1177  }
1178  TVector3 dir = trajpts[1]-trajpts[0];
1179  rb::Track JointTrack(JointCluster,trajpts[0].X(),trajpts[0].Z(),dir.X(),dir.Z(),trka.ID());
1180  if(trka.View() == geo::kY){
1181  JointTrack.SetStart(trajpts[0].Y(),trajpts[0].Z());
1182  JointTrack.SetDir(dir.Y(),dir.Z());
1183  }
1184  for(unsigned int i = 1; i<trajpts.size();++i){
1185  if(trka.View() == geo::kX){ JointTrack.AppendTrajectoryPoint(trajpts[i].X(),trajpts[i].Z()); }
1186  else if(trka.View() == geo::kY){ JointTrack.AppendTrajectoryPoint(trajpts[i].Y(),trajpts[i].Z()); }
1187  }
1188  return JointTrack;
1189  }
1190 
1191  ///////////////////////////////////////////////////////////////////////
1193  {
1194  double score = -1.0;
1195  rb::Track xTrack;
1196  rb::Track yTrack;
1197  if(trka.View() == geo::kX && trkb.View() == geo::kY){
1198  xTrack = trka;
1199  yTrack = trkb;
1200  }
1201  else if(trka.View() == geo::kY && trkb.View() == geo::kX){
1202  xTrack = trkb;
1203  yTrack = trka;
1204  }
1205  else{
1206  // don't have one x track and one y track, can't merge these return meaningless value for score
1207  return score;
1208  }
1209 
1210  // check that trajectory points are good
1211  bool xgood = true;
1212  for(unsigned int i = 0; i<xTrack.NTrajectoryPoints();++i){
1213  if(std::isnan(xTrack.TrajectoryPoint(i).X()) || std::isnan(xTrack.TrajectoryPoint(i).Y()) || std::isnan(xTrack.TrajectoryPoint(i).Z())){
1214  xgood = false;
1215  }
1216  }
1217 
1218  bool ygood = true;
1219  for(unsigned int i = 0; i<yTrack.NTrajectoryPoints();++i){
1220  if(std::isnan(yTrack.TrajectoryPoint(i).X()) || std::isnan(yTrack.TrajectoryPoint(i).Y()) || std::isnan(yTrack.TrajectoryPoint(i).Z())){
1221  ygood = false;
1222  }
1223  }
1224  if(!xgood || !ygood){return score;}
1225 
1226  TVector3 ixtraj = xTrack.Start();
1227  TVector3 fxtraj = xTrack.Stop();
1228  TVector3 iytraj = yTrack.Start();
1229  TVector3 fytraj = yTrack.Stop();
1230 
1231  // make sure the iZ is before the fZ
1232  if(fxtraj.Z() < ixtraj.Z()){ std::swap(ixtraj,fxtraj); }
1233  if(fytraj.Z() < iytraj.Z()){ std::swap(iytraj,fytraj); }
1234 
1235  // get the corresponding planes to these z locations
1236  int ixPlane = kgeo.MatchTrajectoryToPlaneInView(ixtraj,xTrack.View(),xTrack.MinPlane(),xTrack.MaxPlane());
1237  int fxPlane = kgeo.MatchTrajectoryToPlaneInView(fxtraj,xTrack.View(),xTrack.MinPlane(),xTrack.MaxPlane());
1238  int iyPlane = kgeo.MatchTrajectoryToPlaneInView(iytraj,yTrack.View(),yTrack.MinPlane(),yTrack.MaxPlane());
1239  int fyPlane = kgeo.MatchTrajectoryToPlaneInView(fytraj,yTrack.View(),yTrack.MinPlane(),yTrack.MaxPlane());
1240 
1241  // apply basic z geometry cuts to ensure that the two tracks at least have some minimal overlap
1242  if(iyPlane <= fxPlane && fyPlane >= ixPlane){
1243  if(yTrack.ExtentPlane()>1 && xTrack.ExtentPlane()>1){
1244  // find the overlap between the tracks
1245  int minoverlapplane;
1246  int maxoverlapplane;
1247  if(ixPlane<iyPlane){
1248  minoverlapplane = iyPlane;
1249  if(fxPlane<fyPlane){ maxoverlapplane = fxPlane; }
1250  else{ maxoverlapplane = fyPlane; }
1251  }
1252  else{
1253  minoverlapplane = ixPlane;
1254  if(fxPlane<fyPlane){ maxoverlapplane = fxPlane; }
1255  else{ maxoverlapplane = fyPlane; }
1256  }
1257  // get all of the xtrack and ytrack hits
1258  std::vector<art::Ptr<rb::CellHit> > allhits;
1259  for(unsigned int ihit = 0; ihit<trka.NCell(); ++ihit){
1260  allhits.push_back(trka.Cell(ihit));
1261  }
1262  for(unsigned int ihit = 0; ihit<trkb.NCell(); ++ihit){
1263  allhits.push_back(trkb.Cell(ihit));
1264  }
1265  // sort the cellhits by plane
1266  rb::SortByPlane(allhits);
1267 
1268  double overlap = 0;
1269  for(int iPlane = minoverlapplane; iPlane<=maxoverlapplane; ++iPlane){
1270  // for this plane to be overlapped, must have a hit on previous plane and next plane, also need a hit this plane
1271  bool hitcurr = false;
1272  bool hitprev = false;
1273  bool hitnext = false;
1274  for(unsigned int ihit = 0; ihit < allhits.size();++ihit){
1275  if(!hitcurr && allhits[ihit]->Plane() == iPlane){ hitcurr = true; }
1276  if(!hitprev && allhits[ihit]->Plane() == iPlane-1){ hitprev = true; }
1277  if(!hitnext && allhits[ihit]->Plane() == iPlane+1){hitnext = true; }
1278  }
1279  if(hitcurr && hitprev && hitnext){ overlap+=1.0; }
1280 
1281  }
1282  int idiff = TMath::Abs(iyPlane-ixPlane);
1283  int fdiff = TMath::Abs(fyPlane-fxPlane);
1284 
1285  score = ((double)idiff+fdiff)/overlap;
1286  }
1287  else if(yTrack.ExtentPlane() == 1){
1288  // for a vertical y track, matching x track must only appear in the two or less surronding planes
1289  // check that the x track only has hits in the plane before and after the y track
1290  if(ixPlane == (iyPlane-1) && fxPlane == (iyPlane+1)){score = 1.0;}
1291  }
1292  else{
1293  // we are in the case of a vertical x track only, matching y track must only appear in the two surronding planes only
1294  // check that the y track only has hits in the plane before and after the y track
1295  if(iyPlane == (ixPlane-1) && fyPlane == (ixPlane+1)){score = 1.0;}
1296  }
1297  }
1298  else if(ixPlane == fxPlane && iyPlane == fyPlane){
1299  // if x and y tracks appear in only one plane, allow match only if it is adjacent plane
1300  if(ixPlane == (iyPlane-1) || ixPlane == (iyPlane+1)){score = 1.0;}
1301  }
1302  return score;
1303  }
1304 
1305  ///////////////////////////////////////////////////////////////////////////
1307  {
1308  // make a 3d track out of the x and y tracks
1309  rb::Track newtrack;
1310  // add all of the hits to the track.
1311  for(unsigned int i = 0; i<xTrack.NCell();++i){
1312  newtrack.Add(xTrack.Cell(i));
1313  }
1314  for(unsigned int i = 0; i<yTrack.NCell();++i){
1315  newtrack.Add(yTrack.Cell(i));
1316  }
1317  newtrack.SetID(xTrack.ID());
1318 
1319  std::vector<int> xPlanes;
1320  int minxPlane = xTrack.MinPlane();
1321  int maxxPlane = xTrack.MaxPlane();
1322 
1323  std::vector<int> yPlanes;
1324  int minyPlane = yTrack.MinPlane();
1325  int maxyPlane = yTrack.MaxPlane();
1326 
1327  if(xTrack.ExtentPlane() > 1 && yTrack.ExtentPlane() > 1){
1328  // Make a list of all planes that the trajectory points belong to
1329  std::vector<TVector3> xtrjpoints;
1330  double minxz = 99999;
1331  double maxxz = 0;
1332  for(unsigned int i = 0; i<xTrack.NTrajectoryPoints();++i){
1333  int xplane = kgeo.MatchTrajectoryToPlane(xTrack.TrajectoryPoint(i),minxPlane,maxxPlane);
1334  xPlanes.push_back(xplane);
1335  }
1336 
1337  for(unsigned int i =0; i<xPlanes.size(); ++i){
1338  if(i == 0 || i == xPlanes.size()-1){
1339  xtrjpoints.push_back(xTrack.TrajectoryPoint(i));
1340  if(xtrjpoints.back().Z() > maxxz){ maxxz = xtrjpoints.back().Z(); }
1341  if(xtrjpoints.back().Z() < minxz){ minxz = xtrjpoints.back().Z(); }
1342  }
1343  else if(!(xPlanes[i] == xPlanes[i-1] && xPlanes[i] == xPlanes[i+1])){
1344  xtrjpoints.push_back(xTrack.TrajectoryPoint(i));
1345  if(xtrjpoints.back().Z() > maxxz){ maxxz = xtrjpoints.back().Z(); }
1346  if(xtrjpoints.back().Z() < minxz){ minxz = xtrjpoints.back().Z(); }
1347  }
1348  }
1349 
1350  std::vector<TVector3> ytrjpoints;
1351  double minyz = 99999;
1352  double maxyz = 0;
1353  for(unsigned int i = 0; i<yTrack.NTrajectoryPoints();++i){
1354  int yplane = kgeo.MatchTrajectoryToPlane(yTrack.TrajectoryPoint(i),minyPlane,maxyPlane);
1355  yPlanes.push_back(yplane);
1356  }
1357 
1358  for(unsigned int i = 0; i<yPlanes.size(); ++i){
1359  if(i == 0 || i == yPlanes.size()-1){
1360  ytrjpoints.push_back(yTrack.TrajectoryPoint(i));
1361  if(ytrjpoints.back().Z() > maxyz){ maxyz = ytrjpoints.back().Z(); }
1362  if(ytrjpoints.back().Z() < minyz){ minyz = ytrjpoints.back().Z(); }
1363  }
1364  else if(!(yPlanes[i] == yPlanes[i-1] && yPlanes[i] == yPlanes[i+1])){
1365  ytrjpoints.push_back(yTrack.TrajectoryPoint(i));
1366  if(ytrjpoints.back().Z() > maxyz){ maxyz = ytrjpoints.back().Z(); }
1367  if(ytrjpoints.back().Z() < minyz){ minyz = ytrjpoints.back().Z(); }
1368  }
1369  }
1370 
1371  // check if these trajectory points will cause infinite lengths when interpolated
1372  // check low z side first
1373  if(minyz < minxz){
1374  if(xtrjpoints.size() >= 2){
1375  unsigned int xplane = xTrack.MinPlane();
1376  double deltaz = fgeom->Plane(xplane)->Cell(0)->HalfD();
1377  if(xtrjpoints[0].Z() <= xtrjpoints[xtrjpoints.size()-1].Z()){
1378  TVector3 endtraj = xtrjpoints[0];
1379  TVector3 nexttraj = xtrjpoints[1];
1380  if(fabs(endtraj.Z()-nexttraj.Z()) <= 2.0*deltaz){
1381  TVector3 lowztraj, highztraj;
1382  ShiftInterpolationPoints(endtraj,nexttraj,xplane,xTrack,lowztraj,highztraj);
1383  xtrjpoints[0] = lowztraj;
1384  xtrjpoints[1] = highztraj;
1385  }
1386  }
1387  else{
1388  TVector3 endtraj = xtrjpoints[xtrjpoints.size()-1];
1389  TVector3 nexttraj = xtrjpoints[xtrjpoints.size()-2];
1390  if(fabs(endtraj.Z()-nexttraj.Z()) <= 2.0*deltaz){
1391  TVector3 lowztraj, highztraj;
1392  ShiftInterpolationPoints(endtraj,nexttraj,xplane,xTrack,lowztraj,highztraj);
1393  xtrjpoints[xtrjpoints.size()-1] = lowztraj;
1394  xtrjpoints[xtrjpoints.size()-2] = highztraj;
1395  }
1396  }
1397  }
1398  }
1399  else{
1400  if(ytrjpoints.size() >= 2){
1401  unsigned int yplane = yTrack.MinPlane();
1402  double deltaz = fgeom->Plane(yplane)->Cell(0)->HalfD();
1403  if(ytrjpoints[0].Z() <= ytrjpoints[ytrjpoints.size()-1].Z()){
1404  TVector3 endtraj = ytrjpoints[0];
1405  TVector3 nexttraj = ytrjpoints[1];
1406  if(fabs(endtraj.Z()-nexttraj.Z()) <= 2.0*deltaz){
1407  TVector3 lowztraj, highztraj;
1408  ShiftInterpolationPoints(endtraj,nexttraj,yplane,yTrack,lowztraj,highztraj);
1409  ytrjpoints[0] = lowztraj;
1410  ytrjpoints[1] = highztraj;
1411  }
1412  }
1413  else{
1414  TVector3 endtraj = ytrjpoints[ytrjpoints.size()-1];
1415  TVector3 nexttraj = ytrjpoints[ytrjpoints.size()-2];
1416  if(fabs(endtraj.Z()-nexttraj.Z()) <= 2.0*deltaz){
1417  TVector3 lowztraj, highztraj;
1418  ShiftInterpolationPoints(endtraj,nexttraj,yplane,yTrack,lowztraj,highztraj);
1419  ytrjpoints[ytrjpoints.size()-1] = lowztraj;
1420  ytrjpoints[ytrjpoints.size()-2] = highztraj;
1421  }
1422  }
1423  }
1424  }
1425 
1426  // now check the high z side for interpolation problems
1427  if(maxyz > maxxz){
1428  if(xtrjpoints.size() >= 2){
1429  unsigned int xplane = xTrack.MaxPlane();
1430  double deltaz = fgeom->Plane(xplane)->Cell(0)->HalfD();
1431  if(xtrjpoints[0].Z() >= xtrjpoints[xtrjpoints.size()-1].Z()){
1432  TVector3 endtraj = xtrjpoints[0];
1433  TVector3 nexttraj = xtrjpoints[1];
1434  if(fabs(endtraj.Z()-nexttraj.Z()) <= 2.0*deltaz){
1435  TVector3 lowztraj, highztraj;
1436  ShiftInterpolationPoints(endtraj,nexttraj,xplane,xTrack,lowztraj,highztraj);
1437  xtrjpoints[0] = highztraj;
1438  xtrjpoints[1] = lowztraj;
1439  }
1440  }
1441  else{
1442  TVector3 endtraj = xtrjpoints[xtrjpoints.size()-1];
1443  TVector3 nexttraj = xtrjpoints[xtrjpoints.size()-2];
1444  if(fabs(endtraj.Z()-nexttraj.Z()) <= 2.0*deltaz){
1445  TVector3 lowztraj, highztraj;
1446  ShiftInterpolationPoints(endtraj,nexttraj,xplane,xTrack,lowztraj,highztraj);
1447  xtrjpoints[xtrjpoints.size()-1] = highztraj;
1448  xtrjpoints[xtrjpoints.size()-2] = lowztraj;
1449  }
1450  }
1451  }
1452  }
1453  else{
1454  if(ytrjpoints.size() >= 2){
1455  unsigned int yplane = yTrack.MaxPlane();
1456  double deltaz = fgeom->Plane(yplane)->Cell(0)->HalfD();
1457  if(ytrjpoints[0].Z() >= ytrjpoints[ytrjpoints.size()-1].Z()){
1458  TVector3 endtraj = ytrjpoints[0];
1459  TVector3 nexttraj = ytrjpoints[1];
1460  if(fabs(endtraj.Z()-nexttraj.Z()) <= 2.0*deltaz){
1461  TVector3 lowztraj, highztraj;
1462  ShiftInterpolationPoints(endtraj,nexttraj,yplane,yTrack,lowztraj,highztraj);
1463  ytrjpoints[0] = highztraj;
1464  ytrjpoints[1] = lowztraj;
1465  }
1466  }
1467  else{
1468  TVector3 endtraj = ytrjpoints[ytrjpoints.size()-1];
1469  TVector3 nexttraj = ytrjpoints[ytrjpoints.size()-2];
1470  if(fabs(endtraj.Z()-nexttraj.Z()) <= 2.0*deltaz){
1471  TVector3 lowztraj, highztraj;
1472  ShiftInterpolationPoints(endtraj,nexttraj,yplane,yTrack,lowztraj,highztraj);
1473  ytrjpoints[ytrjpoints.size()-1] = highztraj;
1474  ytrjpoints[ytrjpoints.size()-2] = lowztraj;
1475  }
1476  }
1477  }
1478  }
1479 
1480  // Replace track trajectory points with the new trajectory points avoid bad interpolation
1481  xTrack.ClearTrajectoryPoints();
1482  for(unsigned int i = 0; i<xtrjpoints.size(); ++i){
1483  if(i == 0){ xTrack.SetStart(xtrjpoints[0].X(),xtrjpoints[0].Z()); }
1484  else{ xTrack.AppendTrajectoryPoint(xtrjpoints[i].X(),xtrjpoints[i].Z()); }
1485  }
1486 
1487  yTrack.ClearTrajectoryPoints();
1488  for(unsigned int i = 0; i<ytrjpoints.size(); ++i){
1489  if(i == 0){ yTrack.SetStart(ytrjpoints[0].Y(),ytrjpoints[0].Z()); }
1490  else{ yTrack.AppendTrajectoryPoint(ytrjpoints[i].Y(),ytrjpoints[i].Z()); }
1491  }
1492 
1493  // also check if either of these two tracks backscatter in z
1494  bool backscatterx = false;
1495  std::vector<int> xbackscatind;
1496 
1497  for(unsigned int i = 0; i<xtrjpoints.size();++i){
1498  if(i > 0 && i < xtrjpoints.size()-1){
1499  double deltaz1 = xtrjpoints[i].Z()-xtrjpoints[i-1].Z();
1500  double deltaz2 = xtrjpoints[i+1].Z()-xtrjpoints[i].Z();
1501  if((deltaz1*deltaz2)<0){ xbackscatind.push_back(i); }
1502  }
1503  }
1504 
1505  if(xbackscatind.size()>0){ backscatterx = true; }
1506 
1507  bool backscattery = false;
1508  std::vector<int> ybackscatind;
1509  for(unsigned int i = 0; i<ytrjpoints.size();++i){
1510  if(i > 0 && i < ytrjpoints.size()-1){
1511  double deltaz1 = ytrjpoints[i].Z()-ytrjpoints[i-1].Z();
1512  double deltaz2 = ytrjpoints[i+1].Z()-ytrjpoints[i].Z();
1513  if((deltaz1*deltaz2)<0){ ybackscatind.push_back(i); }
1514  }
1515  }
1516 
1517  if(ybackscatind.size()>0){ backscattery = true; }
1518 
1519  std::vector<TVector3> trjpoints;
1520  if((!backscatterx && !backscattery) || (backscatterx && backscattery)){
1521  std::vector<TVector3> newtrjpoints;
1522  for(unsigned int i = 0; i<xtrjpoints.size();++i){
1523  double x, y;
1524  yTrack.InterpolateXY(xtrjpoints[i].Z(),x,y);
1525  xtrjpoints[i].SetY(y);
1526  newtrjpoints.push_back(xtrjpoints[i]);
1527  }
1528  for(unsigned int i = 0; i<ytrjpoints.size();++i){
1529  double x, y;
1530  xTrack.InterpolateXY(ytrjpoints[i].Z(),x,y);
1531  ytrjpoints[i].SetX(x);
1532  newtrjpoints.push_back(ytrjpoints[i]);
1533  }
1534  std::sort(newtrjpoints.begin(),newtrjpoints.end(),sort_traj);
1535  trjpoints = newtrjpoints;
1536  }
1537  if((backscatterx && !backscattery)){
1538  // no backscatter in y can safely use the interpolate function for the y position of x hits trajectories
1539  std::vector<TVector3> newtrjpoints;
1540  for(unsigned int i = 0; i<xtrjpoints.size();++i){
1541  double x, y;
1542  yTrack.InterpolateXY(xtrjpoints[i].Z(),x,y);
1543  xtrjpoints[i].SetY(y);
1544  newtrjpoints.push_back(xtrjpoints[i]);
1545  }
1546 
1547  // insert the start and end points of the x track into the list of backscatters, these serve as end points
1548  xbackscatind.insert(xbackscatind.begin(),0);
1549  xbackscatind.push_back(xtrjpoints.size()-1);
1550  // check how many backscatters are there
1551  int nback = xbackscatind.size();
1552 
1553  // Define a vector with one entry per y trajectory point. The entry is a vector holding
1554  // which xsection(s) the y trajectory point belongs to
1555  std::vector<std::vector<int> > xsectytrj(ytrjpoints.size());
1556  for(unsigned int iytraj = 0; iytraj<ytrjpoints.size(); ++iytraj){
1557  double zy = ytrjpoints[iytraj].Z();
1558  bool foundsect = false;
1559  // Is this point in the overlap of any of the x track sections?
1560  for(unsigned int ixbrk = 0; ixbrk<xbackscatind.size()-1; ++ixbrk){
1561  double z1 = xtrjpoints[xbackscatind[ixbrk]].Z();
1562  double z2 = xtrjpoints[xbackscatind[ixbrk+1]].Z();
1563  if( (z1 < zy && zy < z2) || (z1 > zy && zy > z2) ){
1564  xsectytrj[iytraj].push_back(ixbrk);
1565  foundsect = true;
1566  }
1567  }
1568  if(!foundsect){
1569  // can't find an overlapping x track section this is either higher or lower in z than any x hits, pick the closest section
1570  double closestdist = (xtrjpoints[xbackscatind[0]].Z() - zy);
1571  int closestxsect = 0;
1572  for(unsigned int ixbrk = 0; ixbrk<xbackscatind.size(); ++ixbrk){
1573  double dist = fabs(xtrjpoints[xbackscatind[ixbrk]].Z() - zy);
1574  if(dist<closestdist){
1575  closestxsect = ixbrk;
1576  closestdist = dist;
1577  }
1578  }
1579  xsectytrj[iytraj].push_back(closestxsect);
1580  }
1581  }
1582 
1583  std::vector<TVector3> orderedtrjpoints;
1584  // loop over the x sections in the order that we get them
1585  for(int isect = 0; isect<nback; ++isect){
1586  std::vector<TVector3> xsecttrjs;
1587  std::vector<TVector3> ysecttrjs;
1588  // get all the x trajectory points that belong to this section
1589  // this should just be all the x trajectory points between the backscatters
1590  if(isect != nback-1){
1591  for(int ix = xbackscatind[isect]; ix<=xbackscatind[isect+1]; ++ix){
1592  xsecttrjs.push_back(xtrjpoints[ix]);
1593  }
1594  }
1595  else{
1596  xsecttrjs.push_back(xtrjpoints[xbackscatind[isect]-1]);
1597  xsecttrjs.push_back(xtrjpoints[xbackscatind[isect]]);
1598  }
1599  // find all of the ytrajectory points that belong to this section
1600  for(unsigned int iy = 0; iy<xsectytrj.size(); ++iy){
1601  for(unsigned int iyx = 0; iyx<xsectytrj[iy].size(); ++iyx){
1602  if(xsectytrj[iy][iyx] == isect){ ysecttrjs.push_back(ytrjpoints[iy]); }
1603  }
1604  }
1605  // now we need to interpolate the x point of all the ytrajectory points
1606  for(unsigned int iy = 0; iy<ysecttrjs.size(); ++iy){
1607  // get the zvalue of this point
1608  double zy = ysecttrjs[iy].Z();
1609  bool foundngbr = false;
1610  // loop over the x trajectory points and find the closest points to this zy value
1611  for(unsigned int ix = 0; ix<xsecttrjs.size()-1; ++ix){
1612  double z1 = xsecttrjs[ix].Z();
1613  double z2 = xsecttrjs[ix+1].Z();
1614  if( (z1 < zy && zy < z2) || (z1 > zy && zy > z2) ){
1615  foundngbr = true;
1616  double xint = xsecttrjs[ix].X() + (zy- z1)*(xsecttrjs[ix+1].X()-xsecttrjs[ix].X())/(z2-z1);
1617  ysecttrjs[iy].SetX(xint);
1618  break;
1619  }
1620  }
1621  if(!foundngbr){
1622  // the z position is off the end of the track use two closest pts to interpolate
1623  if(fabs(xsecttrjs[0].Z() - zy) < fabs(xsecttrjs[xsecttrjs.size()-1].Z() - zy)){
1624  double x1 = xsecttrjs[0].X();
1625  double x2 = xsecttrjs[1].X();
1626  double z1 = xsecttrjs[0].Z();
1627  double z2 = xsecttrjs[1].Z();
1628  double xint = x1 + (zy-z1)*(x2-x1)/(z2-z1);
1629  ysecttrjs[iy].SetX(xint);
1630  }
1631  else{
1632  double x1 = xsecttrjs[xsecttrjs.size()-1].X();
1633  double x2 = xsecttrjs[xsecttrjs.size()-2].X();
1634  double z1 = xsecttrjs[xsecttrjs.size()-1].Z();
1635  double z2 = xsecttrjs[xsecttrjs.size()-2].Z();
1636  double xint = x1 + (zy-z1)*(x2-x1)/(z2-z1);
1637  ysecttrjs[iy].SetX(xint);
1638  }
1639  }
1640  }
1641 
1642  std::vector<TVector3> alltrjs;
1643  for(unsigned int i = 0; i<xsecttrjs.size();++i){
1644  if(isect < nback-1){
1645  if(i != xsecttrjs.size()-1){ alltrjs.push_back(xsecttrjs[i]); }
1646  }
1647  else{
1648  if(i != 0){ alltrjs.push_back(xsecttrjs[i]); }
1649  }
1650 
1651  }
1652  for(unsigned int i = 0; i<ysecttrjs.size();++i){
1653  alltrjs.push_back(ysecttrjs[i]);
1654  }
1655  // now order all of the trajectory points in this section
1656  double zDir = xsecttrjs[xsecttrjs.size()-1].Z() - xsecttrjs[0].Z();
1657  std::sort(alltrjs.begin(),alltrjs.end(),sort_traj);
1658  if(zDir < 0){ reverse(alltrjs.begin(),alltrjs.end()); }
1659  for(unsigned int i = 0; i<alltrjs.size(); ++i){
1660  orderedtrjpoints.push_back(alltrjs[i]);
1661  }
1662  }
1663  trjpoints = orderedtrjpoints;
1664  }
1665  if((!backscatterx && backscattery)){
1666  // no backscatter in x can safely use the interpolate function for the x position of y hits trajectories
1667  std::vector<TVector3> newtrjpoints;
1668  for(unsigned int i = 0; i<ytrjpoints.size();++i){
1669  double x, y;
1670  xTrack.InterpolateXY(ytrjpoints[i].Z(),x,y);
1671  ytrjpoints[i].SetX(x);
1672  newtrjpoints.push_back(ytrjpoints[i]);
1673  }
1674 
1675  // insert the start and end points of the y track into the list of backscatters, these serve as end points
1676  ybackscatind.insert(ybackscatind.begin(),0);
1677  ybackscatind.push_back(ytrjpoints.size()-1);
1678  // check how many backscatters are there
1679  int nback = ybackscatind.size();
1680 
1681  // Define a vector with one entry per x trajectory point. The entry is a vector holding
1682  // which ysection(s) the x trajectory point belongs to
1683  std::vector<std::vector<int> > ysectxtrj(xtrjpoints.size());
1684  for(unsigned int ixtraj = 0; ixtraj<xtrjpoints.size(); ++ixtraj){
1685  double zx = xtrjpoints[ixtraj].Z();
1686  bool foundsect = false;
1687  // Is this point in the overlap of any of the x track sections?
1688  for(unsigned int iybrk = 0; iybrk<ybackscatind.size()-1; ++iybrk){
1689  double z1 = ytrjpoints[ybackscatind[iybrk]].Z();
1690  double z2 = ytrjpoints[ybackscatind[iybrk+1]].Z();
1691  if( (z1 < zx && zx < z2) || (z1 > zx && zx > z2) ){
1692  ysectxtrj[ixtraj].push_back(iybrk);
1693  foundsect = true;
1694  }
1695  }
1696  if(!foundsect){
1697  // can't find an overlapping y track section this is either higher or lower in z than any y hits, pick the closest section
1698  double closestdist = (ytrjpoints[ybackscatind[0]].Z() - zx);
1699  int closestysect = 0;
1700  for(unsigned int iybrk = 0; iybrk<ybackscatind.size(); ++iybrk){
1701  double dist = fabs(ytrjpoints[ybackscatind[iybrk]].Z() - zx);
1702  if(dist<closestdist){
1703  closestysect = iybrk;
1704  closestdist = dist;
1705  }
1706  }
1707  ysectxtrj[ixtraj].push_back(closestysect);
1708  }
1709  }
1710 
1711  std::vector<TVector3> orderedtrjpoints;
1712  // loop over the y sections in order that we get them
1713  for(int isect = 0; isect<nback; ++isect){
1714  std::vector<TVector3> xsecttrjs;
1715  std::vector<TVector3> ysecttrjs;
1716  // get all the y trajectory points that belong to this section
1717  // this should just be all the y trajectory points between the backscatters
1718  if(isect != nback-1){
1719  for(int iy = ybackscatind[isect]; iy<=ybackscatind[isect+1]; ++iy){
1720  ysecttrjs.push_back(ytrjpoints[iy]);
1721  }
1722  }
1723  else{
1724  ysecttrjs.push_back(ytrjpoints[ybackscatind[isect]-1]);
1725  ysecttrjs.push_back(ytrjpoints[ybackscatind[isect]]);
1726  }
1727  // find all of the xtrajectory points that belong to this section
1728  for(unsigned int ix = 0; ix<ysectxtrj.size(); ++ix){
1729  for(unsigned int ixy = 0; ixy<ysectxtrj[ix].size(); ++ixy){
1730  if(ysectxtrj[ix][ixy] == isect){ xsecttrjs.push_back(xtrjpoints[ix]); }
1731  }
1732  }
1733  // now we need to interpolate the y point of all the xtrajectory points
1734  for(unsigned int ix = 0; ix<xsecttrjs.size(); ++ix){
1735  // get the zvalue of this point
1736  double zx = xsecttrjs[ix].Z();
1737  bool foundngbr = false;
1738  // loop over the y trajectory points and find the closest points to this zx value
1739  for(unsigned int iy = 0; iy<ysecttrjs.size()-1; ++iy){
1740  double z1 = ysecttrjs[iy].Z();
1741  double z2 = ysecttrjs[iy+1].Z();
1742  if( (z1 < zx && zx < z2) || (z1 > zx && zx > z2) ){
1743  foundngbr = true;
1744  double yint = ysecttrjs[iy].Y() + (zx- z1)*(ysecttrjs[iy+1].Y()-ysecttrjs[iy].Y())/(z2-z1);
1745  xsecttrjs[ix].SetY(yint);
1746  break;
1747  }
1748  }
1749  if(!foundngbr){
1750  // the z position is off the end of the track use two closest pts to interpolate
1751  if(fabs(ysecttrjs[0].Z() - zx) < fabs(ysecttrjs[ysecttrjs.size()-1].Z() - zx)){
1752  double y1 = ysecttrjs[0].Y();
1753  double y2 = ysecttrjs[1].Y();
1754  double z1 = ysecttrjs[0].Z();
1755  double z2 = ysecttrjs[1].Z();
1756  double yint = y1 + (zx-z1)*(y2-y1)/(z2-z1);
1757  xsecttrjs[ix].SetY(yint);
1758  }
1759  else{
1760  double y1 = ysecttrjs[ysecttrjs.size()-1].Y();
1761  double y2 = ysecttrjs[ysecttrjs.size()-2].Y();
1762  double z1 = ysecttrjs[ysecttrjs.size()-1].Z();
1763  double z2 = ysecttrjs[ysecttrjs.size()-2].Z();
1764  double yint = y1 + (zx-z1)*(y2-y1)/(z2-z1);
1765  xsecttrjs[ix].SetY(yint);
1766  }
1767  }
1768  }
1769 
1770  std::vector<TVector3> alltrjs;
1771  for(unsigned int i = 0; i<ysecttrjs.size();++i){
1772  if(isect < nback-1){
1773  if(i != ysecttrjs.size()-1){
1774  alltrjs.push_back(ysecttrjs[i]);
1775  }
1776  }
1777  else{
1778  if(i != 0){
1779  alltrjs.push_back(ysecttrjs[i]);
1780  }
1781  }
1782  }
1783  for(unsigned int i = 0; i<xsecttrjs.size();++i){
1784  alltrjs.push_back(xsecttrjs[i]);
1785  }
1786  // now order all of the trajectory points in this section
1787  double zDir = ysecttrjs[ysecttrjs.size()-1].Z() - ysecttrjs[0].Z();
1788  std::sort(alltrjs.begin(),alltrjs.end(),sort_traj);
1789  if(zDir < 0){ reverse(alltrjs.begin(),alltrjs.end()); }
1790  for(unsigned int i = 0; i<alltrjs.size(); ++i){
1791  orderedtrjpoints.push_back(alltrjs[i]);
1792  }
1793  }
1794  trjpoints = orderedtrjpoints;
1795  }
1796 
1797  newtrack.SetStart(trjpoints[0]);
1798  for(unsigned int i = 1; i<trjpoints.size();++i){
1799  newtrack.AppendTrajectoryPoint(trjpoints[i]);
1800  }
1801 
1802  double tanThetax = TMath::Tan(TMath::ACos(xTrack.Dir().X()));
1803  double tanThetay = TMath::Tan(TMath::ACos(yTrack.Dir().Y()));
1804  double cosThetax = 1/(tanThetax*sqrt(1+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
1805  double cosThetay = 1/(tanThetay*sqrt(1+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
1806  TVector3 newDir(cosThetax,cosThetay,sqrt(1-cosThetax*cosThetax-cosThetay*cosThetay));
1807  if(std::isnan(newDir.X()) || std::isnan(newDir.Y()) || std::isnan(newDir.Z())){
1808  if(trjpoints.size() >= 2){
1809  newDir.SetXYZ(trjpoints[1].X()-trjpoints[0].X(),trjpoints[1].Y()-trjpoints[0].Y(),trjpoints[1].Z()-trjpoints[0].Z());
1810  newDir.Unit();
1811  }
1812  }
1813  newtrack.SetDir(newDir);
1814  }
1815  else if(xTrack.ExtentPlane() == 1){
1816  // check if both tracks only exist in 1 plane
1817  if(yTrack.ExtentPlane() == 1){
1818  // both tracks are in one plane only, assume goes high to low cell
1819  double yHigh[3];
1820  double yLow[3];
1821  fgeom->Plane(yTrack.MinPlane())->Cell(yTrack.MinCell(geo::kY))->GetCenter(yLow);
1822  fgeom->Plane(yTrack.MaxPlane())->Cell(yTrack.MaxCell(geo::kY))->GetCenter(yHigh);
1823  double ydeltaz = fgeom->Plane(yTrack.MinPlane())->Cell(yTrack.MinCell(geo::kY))->HalfD();
1824  double yzlow = yLow[2]-ydeltaz;
1825  double ylow = yLow[1];
1826  double yzhigh = yHigh[2]+ydeltaz;
1827  double yhigh = yHigh[1];
1828  double xHigh[3];
1829  double xLow[3];
1830  fgeom->Plane(xTrack.MinPlane())->Cell(xTrack.MinCell(geo::kX))->GetCenter(xLow);
1831  fgeom->Plane(xTrack.MaxPlane())->Cell(xTrack.MaxCell(geo::kX))->GetCenter(xHigh);
1832  double xdeltaz = fgeom->Plane(xTrack.MinPlane())->Cell(xTrack.MinCell(geo::kX))->HalfD();
1833  double xzlow = xLow[2]-xdeltaz;
1834  double xlow = xLow[0];
1835  double xzhigh = xHigh[2]+xdeltaz;
1836  double xhigh = xHigh[0];
1837  if(yzlow<xzlow){
1838  TVector3 start(xlow+(xhigh-xlow)*(yzlow-xzlow)/(xzhigh-xzlow),ylow,yzlow);
1839  TVector3 stop(xhigh,yhigh+(yhigh-ylow)*(xzhigh-yzhigh)/(yzhigh-yzlow),xzhigh);
1840  newtrack.SetStart(start);
1841  newtrack.AppendTrajectoryPoint(stop);
1842  TVector3 diffs(stop.X()-start.X(),stop.Y()-start.Y(),stop.Z()-start.Z());
1843  newtrack.SetDir(diffs);
1844  }
1845  else{
1846  TVector3 start(xlow,(yhigh-ylow)*(xzlow-yzlow)/(yzhigh-yzlow)+ylow,xzlow);
1847  TVector3 stop(xhigh+(xhigh-xlow)*(yzhigh-xzhigh)/(xzhigh-xzlow),yhigh,yzhigh);
1848  newtrack.SetStart(start);
1849  newtrack.AppendTrajectoryPoint(stop);
1850  TVector3 diffs(stop.X()-start.X(),stop.Y()-start.Y(),stop.Z()-start.Z());
1851  newtrack.SetDir(diffs);
1852  }
1853  }
1854  else{
1855  // only the xtrack is in one plane
1856  // use the center of the top most and bottom most cells as the trajectory points (assume high to low going)
1857  double xHigh[3];
1858  double xLow[3];
1859  fgeom->Plane(xTrack.MinPlane())->Cell(xTrack.MinCell(geo::kX))->GetCenter(xLow);
1860  fgeom->Plane(xTrack.MaxPlane())->Cell(xTrack.MaxCell(geo::kX))->GetCenter(xHigh);
1861  double xdeltaz = fgeom->Plane(xTrack.MinPlane())->Cell(xTrack.MinCell(geo::kX))->HalfD();
1862  double xzlow = xLow[2]-xdeltaz;
1863  double xlow = xLow[0];
1864  double xzhigh = xHigh[2]+xdeltaz;
1865  double xhigh = xHigh[0];
1866  double xslope = (xhigh-xlow)/(xzhigh-xzlow);
1867  for(unsigned int ytraj = 0;ytraj<yTrack.NTrajectoryPoints();++ytraj){
1868  double y = yTrack.TrajectoryPoint(ytraj).Y();
1869  double z = yTrack.TrajectoryPoint(ytraj).Z();
1870  TVector3 start(xlow+xslope*(z-xzlow),y,z);
1871  if(ytraj == 0){newtrack.SetStart(start);}
1872  else{newtrack.AppendTrajectoryPoint(start);}
1873  }
1874  double tanThetax = TMath::Tan(TMath::ACos(xslope/sqrt(1+xslope*xslope)));
1875  double tanThetay = TMath::Tan(TMath::ACos(yTrack.Dir().Y()));
1876  double cosThetax = 1/(tanThetax*sqrt(1+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
1877  double cosThetay = 1/(tanThetay*sqrt(1+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
1878  TVector3 newDir(cosThetax,cosThetay,sqrt(1-cosThetax*cosThetax-cosThetay*cosThetay));
1879  newtrack.SetDir(newDir);
1880  }
1881  }
1882  else{
1883  // only the ytrack is in one plane
1884  // use the center of the top most and bottom most cells as the trajectory points (assume high to low going)
1885  double yHigh[3];
1886  double yLow[3];
1887  fgeom->Plane(yTrack.MinPlane())->Cell(yTrack.MinCell(geo::kY))->GetCenter(yLow);
1888  fgeom->Plane(yTrack.MaxPlane())->Cell(yTrack.MaxCell(geo::kY))->GetCenter(yHigh);
1889  double deltaz = fgeom->Plane(yTrack.MinPlane())->Cell(yTrack.MinCell(geo::kY))->HalfD();
1890  double yzlow = yLow[2]-deltaz;
1891  double ylow = yLow[1];
1892  double yzhigh = yHigh[2]+deltaz;
1893  double yhigh = yHigh[1];
1894  double yslope = (yhigh-ylow)/(yzhigh-yzlow);
1895  for(unsigned int xtraj = 0;xtraj<xTrack.NTrajectoryPoints();++xtraj){
1896  double x = xTrack.TrajectoryPoint(xtraj).X();
1897  double z = xTrack.TrajectoryPoint(xtraj).Z();
1898  TVector3 start(x,ylow+yslope*(z-yzlow),z);
1899  if(xtraj == 0){newtrack.SetStart(start);}
1900  else{ newtrack.AppendTrajectoryPoint(start);}
1901  }
1902  double tanThetay = TMath::Tan(TMath::ACos(yslope/sqrt(1+yslope*yslope)));
1903  double tanThetax = TMath::Tan(TMath::ACos(xTrack.Dir().X()));
1904  double cosThetax = 1/(tanThetax*sqrt(1+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
1905  double cosThetay = 1/(tanThetay*sqrt(1+1.0/(tanThetax*tanThetax)+1.0/(tanThetay*tanThetay)));
1906  TVector3 newDir(cosThetax,cosThetay,sqrt(1-cosThetax*cosThetax-cosThetay*cosThetay));
1907  newtrack.SetDir(newDir);
1908  }
1909  return newtrack;
1910  } // end of ViewMergeTracks function
1911 
1912  //////////////////////////////////////////////////////////
1913  bool KalmanTrackMerge::sort_traj(const TVector3& a, const TVector3& b)
1914  {
1915  return a.Z()<b.Z();
1916  }
1917  //////////////////////////////////////////////////////////
1918 
1919  void KalmanTrackMerge::ShiftInterpolationPoints(TVector3 endtraj, TVector3 nexttraj,
1920  unsigned int plane, rb::Track track,
1921  TVector3 &lowztraj, TVector3 &highztraj)
1922  {
1923 
1924  // need to adjust these trajectory points so interpolation does something sensible
1925  // get the min/max cells on this plane
1926  unsigned int mincell = 10000;
1927  unsigned int maxcell = 0;
1928  for(unsigned int icell = 0; icell < track.NCell(); ++icell){
1929  if(track.Cell(icell)->Plane() == plane){
1930  if(track.Cell(icell)->Cell() > maxcell){ maxcell = track.Cell(icell)->Cell(); }
1931  if(track.Cell(icell)->Cell() < mincell){ mincell = track.Cell(icell)->Cell(); }
1932  }
1933  }
1934  double xyzLow[3];
1935  double xyzHigh[3];
1936  fgeom->Plane(plane)->Cell(mincell)->GetCenter(xyzLow);
1937  fgeom->Plane(plane)->Cell(maxcell)->GetCenter(xyzHigh);
1938  double deltaz = fgeom->Plane(plane)->Cell(mincell)->HalfD();
1939  double zmid = 0.5*xyzLow[2]+0.5*xyzHigh[2];
1940  lowztraj.SetZ(zmid-deltaz);
1941  highztraj.SetZ(zmid+deltaz);
1942  if(track.View() == geo::kX){
1943  if((endtraj-track.TrajectoryPoint(0)).Mag() <= (endtraj-track.TrajectoryPoint(track.NTrajectoryPoints()-1)).Mag()){
1944  if(endtraj.X() <= nexttraj.X()){
1945  lowztraj.SetX(xyzLow[0]);
1946  highztraj.SetX(xyzHigh[0]);
1947  }
1948  else{
1949  lowztraj.SetX(xyzHigh[0]);
1950  highztraj.SetX(xyzLow[0]);
1951  }
1952  }
1953  else{
1954  if(endtraj.X() <= nexttraj.X()){
1955  highztraj.SetX(xyzLow[0]);
1956  lowztraj.SetX(xyzHigh[0]);
1957  }
1958  else{
1959  highztraj.SetX(xyzHigh[0]);
1960  lowztraj.SetX(xyzLow[0]);
1961  }
1962  }
1963  }
1964  else if(track.View() == geo::kY){
1965  if((endtraj-track.TrajectoryPoint(0)).Mag() <= (endtraj-track.TrajectoryPoint(track.NTrajectoryPoints()-1)).Mag()){
1966  if(endtraj.Y() <= nexttraj.Y()){
1967  lowztraj.SetY(xyzLow[1]);
1968  highztraj.SetY(xyzHigh[1]);
1969  }
1970  else{
1971  lowztraj.SetY(xyzHigh[1]);
1972  highztraj.SetY(xyzLow[1]);
1973  }
1974  }
1975  else{
1976  if(endtraj.Y() <= nexttraj.Y()){
1977  highztraj.SetY(xyzLow[1]);
1978  lowztraj.SetY(xyzHigh[1]);
1979  }
1980  else{
1981  highztraj.SetY(xyzHigh[1]);
1982  lowztraj.SetY(xyzLow[1]);
1983  }
1984  }
1985  }
1986 
1987  } // end of ShiftInterpolationPoints function
1988 
1989 
1990 }
1991 
1992 namespace trk{
1993 
1995 
1996 }
virtual geo::View_t View() const
kXorY for 3D clusters.
Definition: Cluster.h:99
size_t NTrajectoryPoints() const
Definition: Track.h:83
int CountMissedCellsOnLine(geo::View_t view, double vi, double zi, double vf, double zf, int maxCells)
Variation of geo::CountCellsOnlineFast function specific to KalmanTrack algorithm.
double HalfL() const
Definition: CellGeo.cxx:198
rb::Track ViewMergeTracks(KalmanGeoHelper &kgeo, rb::Track xTrack, rb::Track yTrack)
art::ServiceHandle< geo::Geometry > fgeom
Detector geometry.
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
virtual void SetStart(TVector3 start)
Definition: Track.cxx:168
double HalfD() const
Definition: CellGeo.cxx:205
bool LineIntersection(double x0, double y0, double x1, double y1, double X0, double Y0, double X1, double Y1, double &x, double &y)
Find the intersection between two line segments.
Definition: Geo.cxx:184
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
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
Float_t y1[n_points_granero]
Definition: compare.C:5
A collection of geometry functions used by the KalmanTrack tracking algorithm.
X or Y views.
Definition: PlaneGeo.h:30
unsigned short Plane() const
Definition: CellHit.h:39
Float_t x1[n_points_granero]
Definition: compare.C:5
void SetHit(art::Ptr< rb::CellHit > chit)
Definition: LocatedChan.cxx:17
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::string fSliceInput
Input folder of slices.
Vertical planes which measure X.
Definition: PlaneGeo.h:28
std::vector< std::vector< double > > avgYPosZRange
Definition: event.h:19
double HalfW() const
Definition: CellGeo.cxx:191
A collection of associated CellHits.
Definition: Cluster.h:47
double MeanZ(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.h:233
bool CanJoinTracks(KalmanGeoHelper &kgeo, rb::Track tracka, rb::Track trackb)
unsigned int MaxCell(geo::View_t view) const
Definition: Cluster.cxx:518
void abs(TH1 *hist)
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
const PlaneGeo * Plane(unsigned int i) const
Representation of a detector channel with positional information.
Definition: LocatedChan.h:15
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
unsigned distance(const T &t1, const T &t2)
static bool sort_traj(const TVector3 &a, const TVector3 &b)
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
Float_t Y
Definition: plot.C:38
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
double dist
Definition: runWimpSim.h:113
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
double CalcMatchScore(KalmanGeoHelper &kgeo, rb::Track tracka, rb::Track trackb)
Float_t Z
Definition: plot.C:38
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
unsigned short Cell() const
Definition: CellHit.h:40
Track finder for cosmic rays.
rb::Track JoinTracks(KalmanGeoHelper &kgeo, rb::Track tracka, rb::Track trackb)
double dy[NP][NC]
double dx[NP][NC]
KalmanTrackMerge(fhicl::ParameterSet const &pset)
double CalcMissFrac(double dist, double w, geo::View_t view)
Estimation of the mip missed hit fraction as a function of the dist through a cell and cell depth w...
double dz[NP][NC]
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
std::vector< std::vector< double > > avgXPosZRange
const double a
int MatchTrajectoryToPlaneInView(TVector3 traj, geo::View_t view, int minplane, int maxplane)
Variation of MatchTrajectoryToPlane that limits the matching plane to be in the input view...
virtual void SetDir(TVector3 dir)
Definition: Prong.cxx:104
void AppendTrajectoryPoint(TVector3 pt)
Definition: Track.cxx:112
T get(std::string const &key) const
Definition: ParameterSet.h:231
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
int evt
float PE() const
Definition: CellHit.h:42
void ClearTrajectoryPoints()
Forget about all trajectory points.
Definition: Track.cxx:359
Collect Geo headers and supply basic geometry functions.
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
Definition: FilterList.h:96
const double j
Definition: BetheBloch.cxx:29
void SetID(int id)
Definition: Cluster.h:74
double DetHalfHeight() const
Definition: View.py:1
z
Definition: test.py:28
size_type size() const
Definition: PtrVector.h:308
bool fBadChannels
Account for bad channels in gaps ?
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::vector< double > avgYPos
double DetHalfWidth() const
std::string fTrackInput
Input folder from track reco.
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
A (plane, cell) pair.
Definition: OfflineChan.h:17
void SortByPlane(std::vector< trk::LocatedChan > &c)
Definition: LocatedChan.cxx:67
void reconfigure(const fhicl::ParameterSet &p)
unsigned int MinCell(geo::View_t view) const
Definition: Cluster.cxx:472
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
float Mag() const
const hit & b
Definition: hits.cxx:21
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.h:250
void ShiftInterpolationPoints(TVector3 endtraj, TVector3 nextraj, unsigned int plane, rb::Track track, TVector3 &lowztraj, TVector3 &highztraj)
assert(nhit_max >=nhit_nbins)
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
const int ID() const
Definition: Cluster.h:75
virtual void InterpolateXY(double z, double &x, double &y) const
Definition: Track.cxx:325
void MatchTracks(KalmanGeoHelper &kgeo, std::vector< rb::Track > xtracks, std::vector< rb::Track > ytracks, std::vector< rb::Track > &matchedtracls, std::vector< rb::Track > &unmatchedtracks)
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
Encapsulate the cell geometry.
Definition: CellGeo.h:25
Float_t X
Definition: plot.C:38
Float_t track
Definition: plot.C:35
Float_t w
Definition: plot.C:20
bool fObeyPreselection
Check rb::IsFiltered?
Encapsulate the geometry of one entire detector (near, far, ndos)
void SortByPlane(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in plane order (low to high).
Definition: CellHit.cxx:124
std::vector< double > avgXPos
void produce(art::Event &evt)
int MatchTrajectoryToPlane(TVector3 traj, int minplane, int maxplane)
Finds the plane best matched to a trajectory point This method takes a trajectory point and finds the...