JMTrackMerge_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief Modified from TrackMerge
3 /// \author Jianming Bian
4 ////////////////////////////////////////////////////////////////////////
5 #include <cmath>
6 #include <fstream>
7 #include <iostream>
8 #include <string>
9 #include <vector>
10 
11 #include "TH1F.h"
12 #include "TTree.h"
13 #include "TVector3.h"
14 
24 #include "fhiclcpp/ParameterSet.h"
26 
27 #include "Geometry/Geometry.h"
29 #include "Calibrator/Calibrator.h"
30 #include "RecoBase/CellHit.h"
31 #include "RecoBase/Cluster.h"
32 #include "RecoBase/Track.h"
33 #include "Utilities/func/MathUtil.h"
34 
35 
36 namespace jmshower {
38  {
39  public:
40  JMTrackMerge(fhicl::ParameterSet const &pset);
41  ~JMTrackMerge();
42  void produce(art::Event &evt);
43  void reconfigure(const fhicl::ParameterSet& pset);
44  void beginJob();
45 
46 
47  private:
50  std::string fOutputFolder; ///< Output folder for Clusters, Prongs, Tracks
51  double fSigmaStart; ///<
52  double fSigmaStop; ///<
53  double fSigmaChi; ///<
54  double fMaxThresh; ///<
55  double fMinThresh; ///<
56  double fMaxScore; ///<
57  int fMinPlaneLength; ///<
58  bool fMCCheater; ///< True/false MCCheater SimHit/SimTrack input
59  TH1F* fHistoStart; ///<
60  TH1F* fHistoStop; ///<
61  TH1F* fHistoStartLog; ///< Easier to read in the reco validation.
62  TH1F* fHistoStopLog; ///< Easier to read in the reco validation.
63  TH1F* fHistoChi; ///<
64  TH1F* fHistoTrack; ///<
65  TH1F* fHistoChiZoom; ///<
66  TH1F* fHistoTrackScore; ///<
67  TH1F* fHistoTrackScoreLog; ///< Easier to read in the reco validation.
68  TH1F* fHistoPurity; ///<
70  TH1F* fHistoEfficiency; ///<
71 
72  TTree* fTrack;
78  double m_trackTrueNuP4[4];
79  double m_trackTrueNuVtx[3];
83 
84 
88 
93  double m_trackTrueVtx[3];
94  double m_trackTrueP4[4];
95  double m_trackTruePartHit[10];
97 
107  //double m_trackTrueRWElecPhotons;
108  //double m_trackTrueRWTotPhotons;
109  //double m_trackTrueTrkRWElecPhotons;
110  //double m_trackTrueTrkRWTotPhotons;
111 
119  double m_trackChisq;
125 
126  double m_trackGap1;
127  double m_trackGap2;
129  double m_trackOut1;
130  double m_trackOut2;
131  double m_trackPE1;
132  double m_trackPE2;
133  double m_trackPEM1;
134  double m_trackPEM2;
139 
140 
141  double m_trackDir[3];
142  double m_trackStart1[3];
143  double m_trackStop1[3];
144  double m_trackStart2[3];
145  double m_trackStop2[3];
146 
147 
151 
156 
157  double m_trackScore;
158  double m_trackCBScore;
164 
170 
177 
178 
179 
184  double m_trackBtTrueVtx[3];
185  double m_trackBtTrueP4[4];
188 
198  //double m_trackBtTrueRWElecPhotons;
199  //double m_trackBtTrueRWTotPhotons;
200  //double m_trackBtTrueTrkRWElecPhotons;
201  //double m_trackBtTrueTrkRWTotPhotons;
202 
213  double m_trackBtPE1;
214  double m_trackBtPE2;
215 
216  double m_trackBtDir[3];
217  double m_trackBtStart1[3];
218  double m_trackBtStop1[3];
219  double m_trackBtStart2[3];
220  double m_trackBtStop2[3];
221  };
222 }
223 
224 
225 namespace jmshower
226 {
227  //...............................................................
229  {
230  reconfigure(pset);
231  produces< std::vector<rb::Track> >();
232  }
233 
234  //......................................................................
236  {
237  pSliceDir = pset.get< std::string >("SliceDir" );
238  fTrackLabel = pset.get< std::string >("TrackLabel" );
239  fSigmaStart = pset.get< double >("SigmaStart" );
240  fSigmaStop = pset.get< double >("SigmaStop" );
241  fSigmaChi = pset.get< double >("SigmaChi" );
242  fMaxThresh = pset.get< double >("MaxThresh" );
243  fMinThresh = pset.get< double >("MinThresh" );
244  fMaxScore = pset.get< double >("MaxScore" );
245  fMinPlaneLength = pset.get< int >("MinPlaneLength");
246  fMCCheater = pset.get< bool >("MCCheater" );
247  }
248 
249  //......................................................................
251  {
252  }
253 
254 
255  //......................................................................
257  {
259  fHistoStart = tfs->make<TH1F>("fHistoStart","StartZ Diff Squared;StartZ Diff;events", 10000,0.,10000.);
260  fHistoStartLog = tfs->make<TH1F>("fHistoStartLog","StartZ Diff Squared;Log10[StartZDiff + 1];events", 50,0.,4.);
261  fHistoChiZoom = tfs->make<TH1F>("fHistoChiZoom","Zoom Chi2;Chi2 diff;events", 100,0.,5.);
262  fHistoStop = tfs->make<TH1F>("fHistoStop","StopZ Diff Squared;StartZ Diff;events", 10000,0.,10000.);
263  fHistoStopLog = tfs->make<TH1F>("fHistoStopLog","StopZ Diff Squared;Log10[StopZ Diff + 1];events", 50,0.,4.);
264  fHistoChi = tfs->make<TH1F>("fHistoChi","Chi2 Diff Squared;Chi2 Diff;events", 1000,0.,1000.);
265  fHistoTrack = tfs->make<TH1F>("fHistoTrack","Tracks;Tracks;events", 30,0.,30.);
266  fHistoTrackScore = tfs->make<TH1F>("fHistoTrackScore","Track Score;score;events", 1000,0.,1000.);
267  fHistoPurity = tfs->make<TH1F>("fHistoPurity",
268  "Purity of Clusters;(Hits of true particle/total hits in cluster);events",101,0.,1.01);
269  fHistoRemainingCells = tfs->make<TH1F>("fHistoRemainingCells",
270  "True hits not put into a 3D prong;Missing Cells;events",100,0.,200.);
271  fHistoEfficiency = tfs->make<TH1F>("fHistoEfficiency",
272  "Ratio of hits of one type in a prong to total hits of that type;Ratio;events",101,0.,1.01);
273  fHistoTrackScoreLog = tfs->make<TH1F>("fHistoTrackScoreLog","Track Score;Log10[score + 1];events", 50,0.,3.);
274 
275  fTrack = tfs->make<TTree>("fTrack","fTrack");
276  fTrack->Branch("trackTrueNuCCNC", &m_trackTrueNuCCNC);
277  fTrack->Branch("trackTrueNuMode", &m_trackTrueNuMode);
278  fTrack->Branch("trackTrueNuPdg", &m_trackTrueNuPdg);
279  fTrack->Branch("trackTrueNuEnergy", &m_trackTrueNuEnergy);
280  fTrack->Branch("trackTrueNuIsFid", &m_trackTrueNuIsFid);
281  fTrack->Branch("trackTrueNuP4[4]", m_trackTrueNuP4);
282  fTrack->Branch("trackTrueNuVtx[3]", m_trackTrueNuVtx);
283  fTrack->Branch("trackEvtTrueEnergy[20]", m_trackEvtTrueEnergy);
284  fTrack->Branch("trackEvtTruePdg[20]", m_trackEvtTruePdg);
285  fTrack->Branch("trackEvtTrueElecPhotons",&m_trackEvtTrueElecPhotons);
286 
287  fTrack->Branch("trackTrueDang", &m_trackTrueDang);
288  fTrack->Branch("trackTrueEDang", &m_trackTrueEDang);
289  fTrack->Branch("trackTruePdg", &m_trackTruePdg);
290  fTrack->Branch("trackTrueId", &m_trackTrueId);
291  fTrack->Branch("trackTrueVtx[3]", m_trackTrueVtx);
292  fTrack->Branch("trackTrueP4[4]", m_trackTrueP4);
293  fTrack->Branch("trackTruePartHit[10]", m_trackTruePartHit);
294  fTrack->Branch("trackTruePartEnergy[10]", m_trackTruePartEnergy);
295  fTrack->Branch("trackTrueDepEnergy", &m_trackTrueDepEnergy);
296  fTrack->Branch("trackTrueElecPhotons", &m_trackTrueElecPhotons);
297  fTrack->Branch("trackTrueTotPhotons", &m_trackTrueTotPhotons);
298  fTrack->Branch("trackTrueTrkElecPhotons", &m_trackTrueTrkElecPhotons);
299  fTrack->Branch("trackTrueTrkTotPhotons", &m_trackTrueTrkTotPhotons);
300  fTrack->Branch("trackTrueTrkElecPhotons1", &m_trackTrueTrkElecPhotons1);
301  fTrack->Branch("trackTrueTrkTotPhotons1", &m_trackTrueTrkTotPhotons1);
302  fTrack->Branch("trackTrueTrkElecPhotons2", &m_trackTrueTrkElecPhotons2);
303  fTrack->Branch("trackTrueTrkTotPhotons2", &m_trackTrueTrkTotPhotons2);
304  fTrack->Branch("trackTrueElecTag",&m_trackTrueElecTag);
305  fTrack->Branch("trackTrueElecTag1",&m_trackTrueElecTag1);
306  fTrack->Branch("trackTrueElecTag2",&m_trackTrueElecTag2);
307  fTrack->Branch("trackBestChisqTag", &m_trackBestChisqTag);
308  fTrack->Branch("trackBestPETag", &m_trackBestPETag);
309  fTrack->Branch("trackBestPlaneTag", &m_trackBestPlaneTag);
310  fTrack->Branch("trackBestOverlapTag", &m_trackBestOverlapTag);
311 
312  fTrack->Branch("trackEvtMinStart", &m_trackEvtMinStart);
313  fTrack->Branch("trackNHit1", &m_trackNHit1);
314  fTrack->Branch("trackNHit2",&m_trackNHit2);;
315  fTrack->Branch("trackNPlane1",&m_trackNPlane1);
316  fTrack->Branch("trackNPlane2",&m_trackNPlane2);
317  fTrack->Branch("trackChisq1",&m_trackChisq1);
318  fTrack->Branch("trackChisq2",&m_trackChisq2);
319  fTrack->Branch("trackChisq",&m_trackChisq);
320  fTrack->Branch("trackxzSlope",&m_trackxzSlope);
321  fTrack->Branch("trackyzSlope",&m_trackyzSlope);
322  fTrack->Branch("trackxzIntercept",&m_trackxzIntercept);
323  fTrack->Branch("trackyzIntercept",&m_trackyzIntercept);
324  fTrack->Branch("trackLength",&m_trackLength);
325 
326  fTrack->Branch("trackGap1",&m_trackGap1);
327  fTrack->Branch("trackGap2",&m_trackGap2);
328  fTrack->Branch("trackOverlap",&m_trackOverlap);
329  fTrack->Branch("trackOut1",&m_trackOut1);
330  fTrack->Branch("trackOut2",&m_trackOut2);
331  fTrack->Branch("trackPE1",&m_trackPE1);
332  fTrack->Branch("trackPE2",&m_trackPE2);
333  fTrack->Branch("trackPEM1",&m_trackPEM1);
334  fTrack->Branch("trackPEM2",&m_trackPEM2);
335  fTrack->Branch("trackPEOut1",&m_trackPEOut1);
336  fTrack->Branch("trackPEOut2",&m_trackPEOut2);
337  fTrack->Branch("trackPESub1",&m_trackPESub1);
338  fTrack->Branch("trackPESub2",&m_trackPESub2);
339 
340 
341  fTrack->Branch("trackDir[3]", m_trackDir);
342  fTrack->Branch("trackStart1[3]", m_trackStart1);
343  fTrack->Branch("trackStop1[3]", m_trackStop1);
344  fTrack->Branch("trackStart2[3]", m_trackStart2);
345  fTrack->Branch("trackStop2[3]", m_trackStop2);
346 
347  fTrack->Branch("trackScore", &m_trackScore);
348  fTrack->Branch("trackCBScore", &m_trackCBScore);
349  fTrack->Branch("trackPEMatch", &m_trackPEMatch);
350  fTrack->Branch("trackPEMatchM", &m_trackPEMatchM);
351  fTrack->Branch("trackPEMatchSub", &m_trackPEMatchSub);
352  fTrack->Branch("trackPlaneMatch", &m_trackPlaneMatch);
353  fTrack->Branch("trackOverlapMatch", &m_trackOverlapMatch);
354 
355  fTrack->Branch("trackBestScore1", &m_trackBestScore1);
356  fTrack->Branch("trackBestCBScore1", &m_trackBestCBScore1);
357  fTrack->Branch("trackBestPEMatch1", &m_trackBestPEMatch1);
358  fTrack->Branch("trackBestPlaneMatch1", &m_trackBestPlaneMatch1);
359  fTrack->Branch("trackBestOverlapMatch1", &m_trackBestOverlapMatch1);
360 
361  fTrack->Branch("trackBestScore2", &m_trackBestScore2);
362  fTrack->Branch("trackBestCBScore2", &m_trackBestCBScore2);
363  fTrack->Branch("trackBestPEMatch2", &m_trackBestPEMatch2);
364  fTrack->Branch("trackBestPlaneMatch2", &m_trackBestPlaneMatch2);
365  fTrack->Branch("trackBestOverlapMatch2", &m_trackBestOverlapMatch2);
366  fTrack->Branch("trackMergeTag", &m_trackMergeTag);
367 
368 
369 
370  fTrack->Branch("trackBtTrueDang", &m_trackBtTrueDang);
371  fTrack->Branch("trackBtTrueEDang", &m_trackBtTrueEDang);
372  fTrack->Branch("trackBtTruePdg", &m_trackBtTruePdg);
373  fTrack->Branch("trackBtTrueId", &m_trackBtTrueId);
374  fTrack->Branch("trackBtTrueVtx[3]", m_trackBtTrueVtx);
375  fTrack->Branch("trackBtTrueP4[4]", m_trackBtTrueP4);
376  fTrack->Branch("trackBtTruePartHit[10]", m_trackBtTruePartHit);
377  fTrack->Branch("trackBtTruePartEnergy[10]", m_trackBtTruePartEnergy);
378  fTrack->Branch("trackBtTrueDepEnergy", &m_trackBtTrueDepEnergy);
379  fTrack->Branch("trackBtTrueElecPhotons", &m_trackBtTrueElecPhotons);
380  fTrack->Branch("trackBtTrueTotPhotons", &m_trackBtTrueTotPhotons);
381  fTrack->Branch("trackBtTrueTrkElecPhotons", &m_trackBtTrueTrkElecPhotons);
382  fTrack->Branch("trackBtTrueTrkTotPhotons", &m_trackBtTrueTrkTotPhotons);
383  fTrack->Branch("trackBtTrueTrkElecPhotons1", &m_trackBtTrueTrkElecPhotons1);
384  fTrack->Branch("trackBtTrueTrkTotPhotons1", &m_trackBtTrueTrkTotPhotons1);
385  fTrack->Branch("trackBtTrueTrkElecPhotons2", &m_trackBtTrueTrkElecPhotons2);
386  fTrack->Branch("trackBtTrueTrkTotPhotons2", &m_trackBtTrueTrkTotPhotons2);
387  fTrack->Branch("trackBtTrueElecTag",&m_trackBtTrueElecTag);
388  fTrack->Branch("trackBtTrueElecTag1",&m_trackBtTrueElecTag1);
389  fTrack->Branch("trackBtTrueElecTag2",&m_trackBtTrueElecTag2);
390 
391  fTrack->Branch("trackBtNHit1", &m_trackBtNHit1);
392  fTrack->Branch("trackBtNHit2",&m_trackBtNHit2);;
393  fTrack->Branch("trackBtNPlane1",&m_trackBtNPlane1);
394  fTrack->Branch("trackBtNPlane2",&m_trackBtNPlane2);
395  fTrack->Branch("trackBtChisq1",&m_trackBtChisq1);
396  fTrack->Branch("trackBtChisq2",&m_trackBtChisq2);
397  fTrack->Branch("trackBtChisq",&m_trackBtChisq);
398  fTrack->Branch("trackBtGap1",&m_trackBtGap1);
399  fTrack->Branch("trackBtGap2",&m_trackBtGap2);
400  fTrack->Branch("trackBtOverlap",&m_trackBtOverlap);
401  fTrack->Branch("trackBtPE1",&m_trackBtPE1);
402  fTrack->Branch("trackBtPE2",&m_trackBtPE2);
403  fTrack->Branch("trackBtDir[3]", m_trackBtDir);
404  fTrack->Branch("trackBtStart1[3]", m_trackBtStart1);
405  fTrack->Branch("trackBtStop1[3]", m_trackBtStop1);
406  fTrack->Branch("trackBtStart2[3]", m_trackBtStart2);
407  fTrack->Branch("trackBtStop2[3]", m_trackBtStop2);
408  }
409 
410  //......................................................................
412  {
413 
415  evt.getByLabel(pSliceDir,slicecol);
416  art::PtrVector<rb::Cluster> slicelist;
417 
418  for(unsigned int i = 0; i < slicecol->size(); ++i){
419  art::Ptr<rb::Cluster>slice(slicecol,i);
420  slicelist.push_back(slice);
421  }
422 
424  evt.getByLabel(fTrackLabel,trackhandle);
425  std::unique_ptr< std::vector<rb::Track> > trackcol(new std::vector<rb::Track> );
426  std::vector<art::Ptr<rb::Track> > track;
427  std::vector<art::Ptr<rb::Track> > xzTrack;
428  std::vector<art::Ptr<rb::Track> > yzTrack;
429 
430  std::vector<rb::Cluster> vCluster;
431  unsigned int cc = 0;
435 
436 
437 
438  std::vector< int > trId;
439  std::vector< int > hitcount;
440  std::vector< std::vector< int > > trackeff;
441 
442  //get slices for timing information
443 
444  std::vector< std::pair< double, double> > tMap; //contains min and max time for each slice in event, first number is min, secont max
445  for (unsigned int i=0; i<slicecol->size(); ++i) {
446  const rb::Cluster& s = (*slicecol)[i];
447  if (s.IsNoise()) continue;
448  std::pair< double, double> p(s.MinTNS(), s.MaxTNS()); //add min and max time as a pair for a given slice
449  tMap.push_back(p);
450  }
451 
452  for(unsigned int i = 0; i < trackhandle->size(); i++){
453  art::Ptr<rb::Track> p(trackhandle, i);
454  track.push_back(p);
455  }
456  for(unsigned int i=0;i<track.size();i++){
457  if (track[i]->NXCell() > 0){
458  xzTrack.push_back(track[i]);
459  }
460  if (track[i]->NYCell() > 0){
461  yzTrack.push_back(track[i]);
462  }
463  }
464 
465  std::vector< std::vector<double> > starts;
466  starts.resize(xzTrack.size()+yzTrack.size());
467  std::vector< std::vector<double> > ends;
468  ends.resize(xzTrack.size()+yzTrack.size());
469  std::vector< std::vector<double> > dirs;
470  dirs.resize(xzTrack.size()+yzTrack.size());
471  std::vector<double> signal;
472  signal.resize(xzTrack.size()+yzTrack.size());
473 
474  fHistoTrack->Fill(track.size());
475 
476  std::vector<int> ntxtrkelecphotonscol;
477  std::vector<int> ntytrkelecphotonscol;
478  std::vector<int> ntxtrktotphotonscol;
479  std::vector<int> ntytrktotphotonscol;
480  ntxtrkelecphotonscol.clear();
481  ntytrkelecphotonscol.clear();
482  ntxtrktotphotonscol.clear();
483  ntytrktotphotonscol.clear();
484 
485 
486  std::vector<std::vector< art::Ptr<rb::CellHit> > > xzCells(xzTrack.size());
487  view = geo::kX;
488  for (unsigned int i=0;i<xzTrack.size();i++) {
489  int trackElecPhotons = 0;
490  int trackTotPhotons = 0;
491 
492  ntxtrkelecphotonscol.push_back(trackElecPhotons);
493  ntxtrktotphotonscol.push_back(trackTotPhotons);
494 
495 
496  for ( unsigned int nh=0;nh<xzTrack[i]->NCell();nh++) {
497  xzCells[i].push_back(xzTrack[i]->Cell(view,nh));
498  }
499  }
500  std::vector<std::vector< art::Ptr<rb::CellHit> > > yzCells(yzTrack.size());
501  view = geo::kY;
502  for (unsigned int i=0;i<yzTrack.size();i++) {
503  int trackElecPhotons = 0;
504  int trackTotPhotons = 0;
505 
506  ntytrkelecphotonscol.push_back(trackElecPhotons);
507  ntytrktotphotonscol.push_back(trackTotPhotons);
508 
509  for ( unsigned int nh=0;nh<yzTrack[i]->NCell();nh++) {
510  yzCells[i].push_back(yzTrack[i]->Cell(view,nh));
511  }
512  }
513 
514  double xzSlope[xzTrack.size()][yzTrack.size()];
515  double xzIntercept[xzTrack.size()][yzTrack.size()];
516  double yzSlope[xzTrack.size()][yzTrack.size()];
517  double yzIntercept[xzTrack.size()][yzTrack.size()];
518 
519  double xyz[3] = {0.,0.,0.};
520 
521  std::vector<art::Ptr<rb::CellHit> >::iterator itr;
522 
523  using util::sqr;
524 
525  double score[xzTrack.size()][yzTrack.size()];
526  double cbscore[xzTrack.size()][yzTrack.size()];
527  double pematch[xzTrack.size()][yzTrack.size()];
528  double planematch[xzTrack.size()][yzTrack.size()];
529  double overlapmatch[xzTrack.size()][yzTrack.size()];
530 
531  double minstart = 99999;
532  for (unsigned int i=0;i<xzTrack.size();i++){
533  for (unsigned int j=0;j<yzTrack.size();j++){
534 
535  double start1 = std::min(xzTrack[i]->Start().Z(), xzTrack[i]->Stop().Z());
536  double stop1 = std::max(xzTrack[i]->Start().Z(), xzTrack[i]->Stop().Z());
537  double start2 = std::min(yzTrack[j]->Start().Z(), yzTrack[j]->Stop().Z());
538  double stop2 = std::max(yzTrack[j]->Start().Z(), yzTrack[j]->Stop().Z());
539 
540  double gap1 = (start1-start2)/double(xzTrack[i]->ExtentPlane()+yzTrack[j]->ExtentPlane());
541  double gap2 = (stop1-stop2)/double(xzTrack[i]->ExtentPlane()+yzTrack[j]->ExtentPlane());
542 
543  if(fabs(gap1)>fabs(gap2)){
544  double gaptmp = gap1;
545  gap1 = gap2;
546  gap2 = gaptmp;
547  }
548 
549  art::PtrVector<rb::CellHit> mergecells;
550  double sig1 = 0.;
551  double sig2 = 0.;
552  double sig1match = 0.;
553  double sig2match = 0.;
554 
555  double wx = 0, wy=0;
556  if(yzTrack[j]->ExtentPlane()==1){
557  double ycellXYZ[3];
558  geom->Plane(yzTrack[j]->YCell(0)->Plane())->Cell(yzTrack[j]->YCell(0)->Cell())->GetCenter(ycellXYZ);
559  wx = ycellXYZ[1];
560  }
561 
562  if(xzTrack[i]->ExtentPlane()==1){
563  double xcellXYZ[3];
564  geom->Plane(xzTrack[i]->XCell(0)->Plane())->Cell(xzTrack[i]->XCell(0)->Cell())->GetCenter(xcellXYZ);
565  wy = xcellXYZ[0];
566  }
567 
568  itr = xzCells[i].begin();
569  std::vector<double> xzViewX;
570  std::vector<double> xzViewZ;
571  std::vector<double> xzViewW;
572  xzViewX.clear();
573  xzViewZ.clear();
574  xzViewW.clear();
575 
576  double xzStart[2] = {0.,0.};
577  double xzEnd[2] = {0.,0.};
578 
579  while (itr != xzCells[i].end()){
580  mergecells.push_back(*itr);
581  double witrk = 0;
582  if(yzTrack[j]->ExtentPlane()==1){witrk=wx;}
583  else{witrk = yzTrack[j]->W((*itr).get());}
584  const double pecor = cal->GetPECorr(*(*itr), witrk);
585  sig1 += pecor;
586  if((int)(*itr)->Plane()>=(int)yzTrack[j]->MinPlane()-1 && (int)(*itr)->Plane()<=(int)yzTrack[j]->MaxPlane()) sig1match +=pecor;
587 
588  geom->Plane((*itr)->Plane())->Cell((*itr)->Cell())->GetCenter(xyz);
589  xzViewX.push_back(xyz[0]);
590  xzViewZ.push_back(xyz[2]);
591  xzViewW.push_back(pecor);
592 
593  itr++;
594  }
595 
596  xzSlope[i][j] = 999999.;
597  xzIntercept[i][j] = xzEnd[1] - xzSlope[i][j]*xzEnd[0];
598  if(fabs(xzStart[0] - xzEnd[0]) > 1.e-2){
599  xzSlope[i][j] = (xzEnd[1] - xzStart[1])/(xzEnd[0] - xzStart[0]);
600  xzIntercept[i][j] = xzEnd[1] - xzSlope[i][j]*xzEnd[0];
601  }
602 
603 
604  itr = yzCells[j].begin();
605  std::vector<double> yzViewY;
606  std::vector<double> yzViewZ;
607  std::vector<double> yzViewW;
608  yzViewY.clear();
609  yzViewZ.clear();
610  yzViewW.clear();
611  double yzStart[2] = {0.,0.};
612  double yzEnd[2] = {0.,0.};
613 
614 
615  while (itr != yzCells[j].end()){
616  mergecells.push_back(*itr);
617  double witrk = 0;
618  if(xzTrack[i]->ExtentPlane()==1){witrk=wy;}
619  else{witrk = xzTrack[i]->W((*itr).get());}
620  const double pecor = cal->GetPECorr(*(*itr), witrk);
621  sig2 += pecor;
622  if((int)(*itr)->Plane()>=(int)xzTrack[i]->MinPlane()-1 && (int)(*itr)->Plane()<=(int)xzTrack[i]->MaxPlane()) sig2match +=pecor;
623 
624  geom->Plane((*itr)->Plane())->Cell((*itr)->Cell())->GetCenter(xyz);
625  yzViewY.push_back(xyz[1]);
626  yzViewZ.push_back(xyz[2]);
627  yzViewW.push_back(pecor);
628 
629  itr++;
630  }
631 
632  yzSlope[i][j] = 999999.;
633  yzIntercept[i][j] = yzEnd[1] - yzSlope[i][j]*yzEnd[0];
634  if(fabs(yzStart[0] - yzEnd[0]) > 1.e-2){
635  yzSlope[i][j] = (yzEnd[1] - yzStart[1])/(yzEnd[0] - yzStart[0]);
636  yzIntercept[i][j] = yzEnd[1] - yzSlope[i][j]*yzEnd[0];
637  }
638 
639 
640  if(sig1+sig2<0.00001)sig1=1;
641  if(sig1match+sig2match<0.00001)sig1match=1;
642 
643 
644  double trackOverlap = -9999;
645  if(start1<=start2&&stop1<=stop2){
646  trackOverlap = stop1 - start2;
647  }else if(start1<=start2&&stop1>stop2){
648  trackOverlap = stop2 - start2;
649  }else if(start1>start2&&stop1<=stop2){
650  trackOverlap = stop1 - start1;
651  }else if(start1>start2&&stop1>stop2){
652  trackOverlap = stop2 - start1;
653  }
654 
655 
656  score[i][j] = (gap1/1.73191e-01)*(gap1/1.73191e-01) + (gap2/1.01610e+0)*(gap2/1.01610e+0);
657  pematch[i][j] = fabs(sig1-sig2)/double(sig1+sig2);
658  cbscore[i][j] = log(score[i][j]+1)/1.57518e-01 + pematch[i][j]/3.45537e-02;
659  planematch[i][j] = fabs((double)xzTrack[i]->ExtentPlane()-(double)yzTrack[j]->ExtentPlane())/double((int)xzTrack[i]->ExtentPlane()+(int)yzTrack[j]->ExtentPlane());
660  overlapmatch[i][j] = trackOverlap/double(xzTrack[i]->ExtentPlane())+trackOverlap/double(yzTrack[j]->ExtentPlane());
661 
662  if(xzTrack[i]->IsNoise()==true&&yzTrack[j]->IsNoise()==true) score[i][j] = 999999.;
663 
664  if(start1<minstart){
665  minstart = start1;
666  }
667  if(start2<minstart){
668  minstart = start2;
669  }
670 
671  }
672  }
673 
674 
675  double xbestscore[xzTrack.size()];
676  double xbestcbscore[xzTrack.size()];
677  double xbestpematch[xzTrack.size()];
678  double xbestplanematch[xzTrack.size()];
679  double xbestoverlapmatch[xzTrack.size()];
680 
681  double ybestscore[yzTrack.size()];
682  double ybestcbscore[yzTrack.size()];
683  double ybestpematch[yzTrack.size()];
684  double ybestplanematch[yzTrack.size()];
685  double ybestoverlapmatch[yzTrack.size()];
686 
687  unsigned int jxbestscore[xzTrack.size()];
688  unsigned int jxbestcbscore[xzTrack.size()];
689  unsigned int jxbestpematch[xzTrack.size()];
690  unsigned int jxbestplanematch[xzTrack.size()];
691  unsigned int jxbestoverlapmatch[xzTrack.size()];
692 
693  unsigned int iybestscore[yzTrack.size()];
694  unsigned int iybestcbscore[yzTrack.size()];
695  unsigned int iybestpematch[yzTrack.size()];
696  unsigned int iybestplanematch[yzTrack.size()];
697  unsigned int iybestoverlapmatch[yzTrack.size()];
698 
699  bool loop = true;
700 
701  while (loop) {
702  for (unsigned int i=0;i<xzTrack.size();i++){
703  xbestscore[i] = 999999;
704  xbestcbscore[i] = 999999;
705  xbestpematch[i] = 999999;
706  xbestplanematch[i] = 999999;
707  xbestoverlapmatch[i] = -999999;
708  jxbestscore[i]=0;
709  jxbestcbscore[i]=0;
710  jxbestpematch[i]=0;
711  jxbestplanematch[i]=0;
712  jxbestoverlapmatch[i]=0;
713  }
714 
715  for (unsigned int j=0;j<yzTrack.size();j++){
716  ybestscore[j] = 999999;
717  ybestcbscore[j] = 999999;
718  ybestpematch[j] = 999999;
719  ybestplanematch[j] = 999999;
720  ybestoverlapmatch[j] = -999999;
721  iybestscore[j]=0;
722  iybestcbscore[j]=0;
723  iybestpematch[j]=0;
724  iybestplanematch[j]=0;
725  iybestoverlapmatch[j]=0;
726  }
727 
728  for (unsigned int i=0;i<xzTrack.size();i++){
729  for (unsigned int j=0;j<yzTrack.size();j++){
730  if (score[i][j] < xbestscore[i]){
731  xbestscore[i] = score[i][j];
732  jxbestscore[i] = j;
733  }
734  if (score[i][j] < ybestscore[j]){
735  ybestscore[j] = score[i][j];
736  iybestscore[j] = i;
737  }
738  if (cbscore[i][j] < xbestcbscore[i]){
739  xbestcbscore[i] = cbscore[i][j];
740  jxbestcbscore[i] = j;
741  }
742  if (cbscore[i][j] < ybestcbscore[j]){
743  ybestcbscore[j] = cbscore[i][j];
744  iybestcbscore[j] = i;
745  }
746 
747  if (pematch[i][j] < xbestpematch[i]){
748  xbestpematch[i] = pematch[i][j];
749  jxbestpematch[i] = j;
750  }
751  if (pematch[i][j] < ybestpematch[j]){
752  ybestpematch[j] = pematch[i][j];
753  iybestpematch[j] = i;
754  }
755 
756  if (planematch[i][j] < xbestplanematch[i]){
757  xbestplanematch[i] = planematch[i][j];
758  jxbestplanematch[i] = j;
759  }
760  if (planematch[i][j] < ybestplanematch[j]){
761  ybestplanematch[j] = planematch[i][j];
762  iybestplanematch[j] = i;
763  }
764 
765  if (overlapmatch[i][j] > xbestoverlapmatch[i]){
766  xbestoverlapmatch[i] = overlapmatch[i][j];
767  jxbestoverlapmatch[i] = j;
768  }
769  if (overlapmatch[i][j] > ybestoverlapmatch[j]){
770  ybestoverlapmatch[j] = overlapmatch[i][j];
771  iybestoverlapmatch[j] = i;
772  }
773  }
774  }
775 
776  bool matchtag = false;
777  for (unsigned int i=0;i<xzTrack.size();i++){
778  if(matchtag == true)break;
779  for (unsigned int j=0;j<yzTrack.size();j++){
780  if(matchtag == true)break;
781  if(i==iybestcbscore[j]||j==jxbestcbscore[i]){
782  if(score[i][j]<50000.&&pematch[i][j]<10.5&&overlapmatch[i][j]>-1000.1){
783 
784  double zStart = yzTrack[j]->Start().Z();
785  double zStop = yzTrack[j]->Stop().Z();
786  if (xzTrack[i]->Start().Z() < yzTrack[j]->Start().Z()) zStart = xzTrack[i]->Start().Z();
787  if (xzTrack[i]->Stop().Z() > yzTrack[j]->Stop().Z()) zStop = xzTrack[i]->Stop().Z();
788 
789  double xyz1[3] = {(xzSlope[i][j]*zStart +xzIntercept[i][j]),
790  (yzSlope[i][j]*zStart +yzIntercept[i][j]), zStart};
791  double xyz2[3] = {(xzSlope[i][j]*zStop +xzIntercept[i][j]),
792  (yzSlope[i][j]*zStop +yzIntercept[i][j]), zStop};
793 
794 
795 
796  double length = util::pythag(zStop-zStart,
797  xyz2[0] - xyz1[0],
798  xyz2[1] - xyz1[1]);
799 
800 
801  if(xzSlope[i][j]>900000||yzSlope[i][j]>900000||fabs(xyz1[0])>geom->DetHalfWidth()||fabs(xyz1[1])>geom->DetHalfHeight()||fabs(xyz2[0])>geom->DetHalfWidth()||fabs(xyz2[1])>geom->DetHalfHeight()||(length>600&&std::min(xzTrack[i]->ExtentPlane(),yzTrack[j]->ExtentPlane())<=20)){
802  xyz1[0] = xzTrack[i]->Start().X();
803  xyz1[1] = yzTrack[j]->Start().Y();
804  xyz1[2] = zStart;
805  xyz2[0] = xzTrack[i]->Stop().X();
806  xyz2[1] = yzTrack[j]->Stop().Y();
807  xyz2[2] = zStop;
808  }
809 
810  length = util::pythag(zStop-zStart,
811  xyz2[0] - xyz1[0],
812  xyz2[1] - xyz1[1]);
813 
814  double dcosx = (xyz2[0] - xyz1[0])/length;
815  double dcosy = (xyz2[1] - xyz1[1])/length;
816  double dcosz = (zStop-zStart)/length;
817 
818  TVector3 xydir(dcosx,dcosy,dcosz);
819 
820 
821  art::PtrVector<rb::CellHit> mergecells;
822 
823  double sig1 = 0.;
824  double sig2 = 0.;
825  double sig = 0.;
826  double sig1match = 0.;
827  double sig2match = 0.;
828 
829  double wx = 0, wy=0;
830  if(yzTrack[j]->ExtentPlane()==1){
831  double ycellXYZ[3];
832  geom->Plane(yzTrack[j]->YCell(0)->Plane())->Cell(yzTrack[j]->YCell(0)->Cell())->GetCenter(ycellXYZ);
833  wx = ycellXYZ[1];
834  }
835 
836  if(xzTrack[i]->ExtentPlane()==1){
837  double xcellXYZ[3];
838  geom->Plane(xzTrack[i]->XCell(0)->Plane())->Cell(xzTrack[i]->XCell(0)->Cell())->GetCenter(xcellXYZ);
839  wy = xcellXYZ[0];
840  }
841 
842  itr = xzCells[i].begin();
843  while (itr != xzCells[i].end()){
844  mergecells.push_back(*itr);
845  double witrk = 0;
846  if(yzTrack[j]->ExtentPlane()==1){witrk=wx;}
847  else{witrk = yzTrack[j]->W((*itr).get());}
848  const double pecor = cal->GetPECorr(*(*itr), witrk);
849  sig1 += pecor;
850  if((int)(*itr)->Plane()>=(int)yzTrack[j]->MinPlane()-1 && (int)(*itr)->Plane()<=(int)yzTrack[j]->MaxPlane()) sig1match +=pecor;
851 
852  itr++;
853  }
854  itr = yzCells[j].begin();
855  while (itr != yzCells[j].end()){
856  mergecells.push_back(*itr);
857  double witrk = 0;
858  if(xzTrack[i]->ExtentPlane()==1){witrk=wy;}
859  else{witrk = xzTrack[i]->W((*itr).get());}
860  const double pecor = cal->GetPECorr(*(*itr), witrk);
861  sig2 += pecor;
862  if((int)(*itr)->Plane()>=(int)xzTrack[i]->MinPlane()-1 && (int)(*itr)->Plane()<=(int)xzTrack[i]->MaxPlane()) sig2match +=pecor;
863  itr++;
864  }
865 
866  bool SliceMatch = false;
867  for( unsigned int itm = 0; itm < tMap.size(); ++itm){
868 
869  if ((xzTrack[i]->MinTNS() >= tMap[itm].first) && (xzTrack[i]->MaxTNS() <= tMap[itm].second)){
870  if ((yzTrack[j]->MinTNS() >= tMap[itm].first) && (yzTrack[j]->MaxTNS() <= tMap[itm].second)) SliceMatch = true;
871  break;
872  }
873  }
874  if (SliceMatch == false){
875  score[i][j] = 999999.;
876  cbscore[i][j] = 999999.;
877  pematch[i][j] = 999999.;
878  continue;
879  }
880 
881  rb::Cluster Clust(mergecells);
882 
883  vCluster.push_back(Clust);
884  starts[cc].push_back(xyz1[0]);
885  starts[cc].push_back(xyz1[1]);
886  starts[cc].push_back(xyz1[2]);
887  ends[cc].push_back(xyz2[0]);
888  ends[cc].push_back(xyz2[1]);
889  ends[cc].push_back(xyz2[2]);
890  dirs[cc].push_back(dcosx);
891  dirs[cc].push_back(dcosy);
892  dirs[cc].push_back(dcosz);
893  signal[cc] = sig;
894  cc++;
895 
896  for(unsigned int ii=0;ii<xzTrack.size();ii++){
897  score[ii][j] = 999999;
898  cbscore[ii][j] = 999999;
899  pematch[ii][j] = 999999;
900  }
901 
902  for(unsigned int jj=0;jj<yzTrack.size();jj++){
903  score[i][jj] = 999999;
904  cbscore[i][jj] = 999999;
905  pematch[i][jj] = 999999;
906  }
907 
908  matchtag = true;
909 
910  }
911  }
912  }
913  }
914  if(matchtag==false)loop=false;
915  }
916 
917  for (unsigned int i=0;i<xzTrack.size();i++){
918  xbestscore[i] = 999999;
919  xbestcbscore[i] = 999999;
920  xbestpematch[i] = 999999;
921  xbestplanematch[i] = 999999;
922  xbestoverlapmatch[i] = -999999;
923  jxbestscore[i]=0;
924  jxbestcbscore[i]=0;
925  jxbestpematch[i]=0;
926  jxbestplanematch[i]=0;
927  jxbestoverlapmatch[i]=0;
928  }
929 
930  for (unsigned int j=0;j<yzTrack.size();j++){
931  ybestscore[j] = 999999;
932  ybestcbscore[j] = 999999;
933  ybestpematch[j] = 999999;
934  ybestplanematch[j] = 999999;
935  ybestoverlapmatch[j] = -999999;
936  iybestscore[j]=0;
937  iybestcbscore[j]=0;
938  iybestpematch[j]=0;
939  iybestplanematch[j]=0;
940  iybestoverlapmatch[j]=0;
941  }
942 
943  for (unsigned int i=0;i<xzTrack.size();i++){
944  for (unsigned int j=0;j<yzTrack.size();j++){
945  if (score[i][j] < xbestscore[i]){
946  xbestscore[i] = score[i][j];
947  jxbestscore[i] = j;
948  }
949  if (score[i][j] < ybestscore[j]){
950  ybestscore[j] = score[i][j];
951  iybestscore[j] = i;
952  }
953  if (cbscore[i][j] < xbestcbscore[i]){
954  xbestcbscore[i] = cbscore[i][j];
955  jxbestcbscore[i] = j;
956  }
957  if (cbscore[i][j] < ybestcbscore[j]){
958  ybestcbscore[j] = cbscore[i][j];
959  iybestcbscore[j] = i;
960  }
961 
962  if (pematch[i][j] < xbestpematch[i]){
963  xbestpematch[i] = pematch[i][j];
964  jxbestpematch[i] = j;
965  }
966  if (pematch[i][j] < ybestpematch[j]){
967  ybestpematch[j] = pematch[i][j];
968  iybestpematch[j] = i;
969  }
970 
971  if (planematch[i][j] < xbestplanematch[i]){
972  xbestplanematch[i] = planematch[i][j];
973  jxbestplanematch[i] = j;
974  }
975  if (planematch[i][j] < ybestplanematch[j]){
976  ybestplanematch[j] = planematch[i][j];
977  iybestplanematch[j] = i;
978  }
979 
980  if (overlapmatch[i][j] > xbestoverlapmatch[i]){
981  xbestoverlapmatch[i] = overlapmatch[i][j];
982  jxbestoverlapmatch[i] = j;
983  }
984  if (overlapmatch[i][j] > ybestoverlapmatch[j]){
985  ybestoverlapmatch[j] = overlapmatch[i][j];
986  iybestoverlapmatch[j] = i;
987  }
988  }
989  }
990  bool dumpmatch = false;
991  if(dumpmatch==true){
992  for (unsigned int i=0;i<xzTrack.size();i++){
993  std::cout<<xbestscore[i]<<std::endl;
994  std::cout<<xbestcbscore[i]<<std::endl;
995  std::cout<<xbestpematch[i]<<std::endl;
996  std::cout<<xbestplanematch[i]<<std::endl;
997  std::cout<<xbestoverlapmatch[i]<<std::endl;
998  std::cout<<jxbestscore[i]<<std::endl;
999  std::cout<<jxbestcbscore[i]<<std::endl;
1000  std::cout<<jxbestpematch[i]<<std::endl;
1001  std::cout<<jxbestplanematch[i]<<std::endl;
1002  std::cout<<jxbestoverlapmatch[i]<<std::endl;
1003  }
1004 
1005 
1006  for (unsigned int j=0;j<yzTrack.size();j++){
1007  std::cout<<ybestscore[j]<<std::endl;
1008  std::cout<<ybestcbscore[j]<<std::endl;
1009  std::cout<<ybestpematch[j]<<std::endl;
1010  std::cout<<ybestplanematch[j]<<std::endl;
1011  std::cout<<ybestoverlapmatch[j]<<std::endl;
1012  std::cout<<iybestscore[j]<<std::endl;
1013  std::cout<<iybestcbscore[j]<<std::endl;
1014  std::cout<<iybestpematch[j]<<std::endl;
1015  std::cout<<iybestplanematch[j]<<std::endl;
1016  std::cout<<iybestoverlapmatch[j]<<std::endl;
1017  }
1018 
1019  }
1020 
1021 
1022  for (unsigned int iMin=0;iMin<xzTrack.size();iMin++){
1023  for (unsigned int jMin=0;jMin<yzTrack.size();jMin++){
1024  if(score[iMin][jMin]>900000)continue;
1025  if(cbscore[iMin][jMin]>900000)continue;
1026  if(pematch[iMin][jMin]>900000)continue;
1027 
1028  double zStart = yzTrack[jMin]->Start().Z();
1029  double zStop = yzTrack[jMin]->Stop().Z();
1030  if (xzTrack[iMin]->Start().Z() < yzTrack[jMin]->Start().Z()) zStart = xzTrack[iMin]->Start().Z();
1031  if (xzTrack[iMin]->Stop().Z() > yzTrack[jMin]->Stop().Z()) zStop = xzTrack[iMin]->Stop().Z();
1032 
1033  double xyz1[3] = {(xzSlope[iMin][jMin]*zStart +xzIntercept[iMin][jMin]),
1034  (yzSlope[iMin][jMin]*zStart +yzIntercept[iMin][jMin]), zStart};
1035  double xyz2[3] = {(xzSlope[iMin][jMin]*zStop +xzIntercept[iMin][jMin]),
1036  (yzSlope[iMin][jMin]*zStop +yzIntercept[iMin][jMin]), zStop};
1037 
1038  double length = util::pythag(zStop-zStart,
1039  xyz2[0] - xyz1[0],
1040  xyz2[1] - xyz1[1]);
1041 
1042 
1043  if(xzSlope[iMin][jMin]>900000||yzSlope[iMin][jMin]>900000||fabs(xyz1[0])>geom->DetHalfWidth()||fabs(xyz1[1])>geom->DetHalfHeight()||fabs(xyz2[0])>geom->DetHalfWidth()||fabs(xyz2[1])>geom->DetHalfHeight()||(length>600&&std::min(xzTrack[iMin]->ExtentPlane(),yzTrack[jMin]->ExtentPlane())<=20)){
1044  xyz1[0] = xzTrack[iMin]->Start().X();
1045  xyz1[1] = yzTrack[jMin]->Start().Y();
1046  xyz1[2] = zStart;
1047  xyz2[0] = xzTrack[iMin]->Stop().X();
1048  xyz2[1] = yzTrack[jMin]->Stop().Y();
1049  xyz2[2] = zStop;
1050  }
1051 
1052  length = util::pythag(zStop-zStart,
1053  xyz2[0] - xyz1[0],
1054  xyz2[1] - xyz1[1]);
1055 
1056 
1057 
1058  double dcosx = (xyz2[0] - xyz1[0])/length;
1059  double dcosy = (xyz2[1] - xyz1[1])/length;
1060  double dcosz = (zStop-zStart)/length;
1061 
1062 
1063  TVector3 xydir(dcosx,dcosy,dcosz);
1064 
1065 
1066  art::PtrVector<rb::CellHit> mergecells;
1067 
1068  double sig1 = 0.;
1069  double sig2 = 0.;
1070  double sig = 0.;
1071  double sig1match = 0.;
1072  double sig2match = 0.;
1073 
1074  double wx = 0, wy=0;
1075  if(yzTrack[jMin]->ExtentPlane()==1){
1076  double ycellXYZ[3];
1077  geom->Plane(yzTrack[jMin]->YCell(0)->Plane())->Cell(yzTrack[jMin]->YCell(0)->Cell())->GetCenter(ycellXYZ);
1078  wx = ycellXYZ[1];
1079  }
1080 
1081  if(xzTrack[iMin]->ExtentPlane()==1){
1082  double xcellXYZ[3];
1083  geom->Plane(xzTrack[iMin]->XCell(0)->Plane())->Cell(xzTrack[iMin]->XCell(0)->Cell())->GetCenter(xcellXYZ);
1084  wy = xcellXYZ[0];
1085  }
1086 
1087  itr = xzCells[iMin].begin();
1088  while (itr != xzCells[iMin].end()){
1089  mergecells.push_back(*itr);
1090  double witrk = 0;
1091  if(yzTrack[jMin]->ExtentPlane()==1){witrk=wx;}
1092  else{witrk = yzTrack[jMin]->W((*itr).get());}
1093 
1094  const double pecor = cal->GetPECorr(*(*itr), witrk);
1095 
1096  sig1 += pecor;
1097  if((int)(*itr)->Plane()>=(int)yzTrack[jMin]->MinPlane()-1 && (int)(*itr)->Plane()<=(int)yzTrack[jMin]->MaxPlane()) sig1match +=pecor;
1098 
1099  itr++;
1100  }
1101  itr = yzCells[jMin].begin();
1102  while (itr != yzCells[jMin].end()){
1103  mergecells.push_back(*itr);
1104  double witrk = 0;
1105  if(xzTrack[iMin]->ExtentPlane()==1){witrk=wy;}
1106  else{witrk = xzTrack[iMin]->W((*itr).get());}
1107  const double pecor = cal->GetPECorr(*(*itr), witrk);
1108  sig2 += pecor;
1109  if((int)(*itr)->Plane()>=(int)xzTrack[iMin]->MinPlane()-1 && (int)(*itr)->Plane()<=(int)xzTrack[iMin]->MaxPlane()) sig2match +=pecor;
1110  itr++;
1111  }
1112 
1113 
1114  int mtag = 0;
1115  if(pematch[iMin][jMin]<0.5&&score[iMin][jMin]<50.&&overlapmatch[iMin][jMin]>-1000.1&&(jMin==jxbestscore[iMin]||iMin==iybestscore[jMin]||jMin==jxbestpematch[iMin]||iMin==iybestpematch[jMin])) mtag = 1;
1116 
1117  if(mtag==0)continue;
1118 
1119  // check to make sure both track merge candidates are part of the same slice
1120  bool SliceMatch = false;
1121  for( unsigned int itm = 0; itm < tMap.size(); ++itm){
1122  //get slice window for xzTrack
1123  if ((xzTrack[iMin]->MinTNS() >= tMap[itm].first) && (xzTrack[iMin]->MaxTNS() <= tMap[itm].second)){
1124  //check to see if yzTrack is in same window, if so allow match
1125  if ((yzTrack[jMin]->MinTNS() >= tMap[itm].first) && (yzTrack[jMin]->MaxTNS() <= tMap[itm].second)) SliceMatch = true;
1126  break;
1127  }
1128  }
1129  //if score match is false, remove that combination from 3D track candidates and continue
1130  if (SliceMatch == false){
1131  score[iMin][jMin] = 999999.;
1132  cbscore[iMin][jMin] = 999999.;
1133  pematch[iMin][jMin] = 999999.;
1134  continue;
1135  }
1136 
1137 
1138  rb::Cluster Clust(mergecells);
1139 
1140  score[iMin][jMin] = 999999.;
1141  cbscore[iMin][jMin] = 999999.;
1142  pematch[iMin][jMin] = 999999.;
1143 
1144  vCluster.push_back(Clust);
1145  starts[cc].push_back(xyz1[0]);
1146  starts[cc].push_back(xyz1[1]);
1147  starts[cc].push_back(xyz1[2]);
1148  ends[cc].push_back(xyz2[0]);
1149  ends[cc].push_back(xyz2[1]);
1150  ends[cc].push_back(xyz2[2]);
1151  dirs[cc].push_back(dcosx);
1152  dirs[cc].push_back(dcosy);
1153  dirs[cc].push_back(dcosz);
1154  signal[cc] = sig;
1155  cc++;
1156  std::vector< int > trackcount;
1157 
1158  }
1159  }
1160 
1161 
1162  std::unique_ptr< std::vector<rb::Cluster> > clustcol(new std::vector<rb::Cluster> );
1163  for(unsigned int i = 0; i < vCluster.size(); i++)
1164  {
1165  clustcol->push_back(vCluster[i]);
1166  }
1167  // put the clusters into the event
1168  for( size_t i = 0; i < vCluster.size(); i++)
1169  {
1170  TVector3 start(starts[i][0], starts[i][1], starts[i][2]);
1171  TVector3 end(ends[i][0], ends[i][1], ends[i][2]);
1172  TVector3 dir(dirs[i][0], dirs[i][1], dirs[i][2]);
1173  rb::Track track(vCluster[i], start, dir);
1174  track.AppendTrajectoryPoint(end);
1175 
1176  trackcol->push_back(track);
1177 
1178  }
1179  evt.put(std::move(trackcol));
1180  return; // kFailed if you want to fail the event
1181 
1182  }
1183 
1184 }//end namespace jmshower
1185 
1186 namespace jmshower
1187 {
1189 }
T max(const caf::Proxy< T > &a, T b)
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
void produce(art::Event &evt)
Vertical planes which measure X.
Definition: PlaneGeo.h:28
A collection of associated CellHits.
Definition: Cluster.h:47
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
bool fMCCheater
True/false MCCheater SimHit/SimTrack input.
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
TODO.
Definition: FillPIDs.h:11
TH1F * fHistoStopLog
Easier to read in the reco validation.
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
Float_t Z
Definition: plot.C:38
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
void reconfigure(const fhicl::ParameterSet &pset)
const XML_Char * s
Definition: expat.h:262
TH1F * fHistoStartLog
Easier to read in the reco validation.
double MinTNS() const
Definition: Cluster.cxx:482
length
Definition: demo0.py:21
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
void AppendTrajectoryPoint(TVector3 pt)
Definition: Track.cxx:112
T get(std::string const &key) const
Definition: ParameterSet.h:231
const double j
Definition: BetheBloch.cxx:29
double DetHalfHeight() const
JMTrackMerge(fhicl::ParameterSet const &pset)
TH1F * fHistoTrackScoreLog
Easier to read in the reco validation.
double MaxTNS() const
Definition: Cluster.cxx:528
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
double DetHalfWidth() const
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
void cc()
Definition: test_ana.C:28
void geom(int which=0)
Definition: geom.C:163
double GetPECorr(rb::CellHit const &cellhit, double w)
T min(const caf::Proxy< T > &a, T b)
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
std::string fOutputFolder
Output folder for Clusters, Prongs, Tracks.
Float_t e
Definition: plot.C:35
Float_t track
Definition: plot.C:35
Encapsulate the geometry of one entire detector (near, far, ndos)