JMClusterMerge_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 ///
3 /// \brief Modified from ClusterMerge
4 ///
5 /// \author $Author: bianjm $
6 ////////////////////////////////////////////////////////////////////////
7 #include <cmath>
8 #include <iostream>
9 #include <list>
10 #include <string>
11 #include <vector>
12 
13 // ROOT includes
14 #include "TVector3.h"
15 #include "TH1D.h"
16 #include "TTree.h"
17 
18 // ART includes
29 #include "fhiclcpp/ParameterSet.h"
31 
32 // NOvA includes
33 #include "Calibrator/Calibrator.h"
34 #include "Geometry/Geometry.h"
35 #include "GeometryObjects/CellGeo.h"
36 #include "GeometryObjects/Geo.h"
38 #include "RecoBase/CellHit.h"
39 #include "RecoBase/Cluster.h"
40 #include "RecoBase/Track.h"
41 #include "Utilities/func/MathUtil.h"
42 
43 
44 
45 namespace jmshower {
47  {
48  public:
49  explicit JMClusterMerge(fhicl::ParameterSet const& pset);
50  ~JMClusterMerge();
51  void reconfigure(const fhicl::ParameterSet& pset);
52  void produce(art::Event& evt);
53  void beginJob();
54 
55 
56  private:
59  int fDebug; ///< Controls debug output 0=none, 99=scream
60  double fDHitGood; ///< Maximum distance from the line for a hit to be considered part of it
61  int fMinHit;
62  double fGoodness;
63  int fMaxPlane;
64  bool fMCCheater; ///< True/false MCCheater SimHit/SimTrack input
66  TH1F* fHistoGoodnessLog; ///< Easier to read in the reco validation.
67  TH1F* fHistoEnergy;
68  TH1F* fHistoPurity;
70 
71  TTree* fTrack;
76  double m_trackChisq1;
77  double m_trackChisq2;
78  double m_trackChisq;
80  double m_trackGap;
81  double m_trackDHit;
82 
83  //int m_trackTrueNuCCNC;
84  //int m_trackTrueNuMode;
85  //int m_trackTrueNuPdg;
86  //double m_trackTrueNuEnergy;
87  //int m_trackTrueNuIsFid;
88  //double m_trackTrueNuP4[4];
89  //double m_trackTrueNuVtx[3];
90  //int m_trackEvtTruePdg[20];
91  //double m_trackEvtTrueEnergy[20];
92  //double m_trackEvtTrueElecPhotons;
93 
94 
95  //double m_trackTrueDang;
96  //double m_trackTrueEDang;
97  //int m_trackTruePdg;
98  //int m_trackTrueId;
99  //double m_trackTrueVtx[3];
100  //double m_trackTrueP4[4];
101  //double m_trackTruePartHit[10];
102  //double m_trackTruePartEnergy[10];
103  //double m_trackTrueDepEnergy;
104  //double m_trackTrueElecPhotons;
105  //double m_trackTrueTotPhotons;
106  //double m_trackTrueTrkElecPhotons;
107  //double m_trackTrueTrkTotPhotons;
108  //double m_trackTrueRWElecPhotons;
109  //double m_trackTrueRWTotPhotons;
110  //double m_trackTrueTrkRWElecPhotons;
111  //double m_trackTrueTrkRWTotPhotons;
112 
113  };
114 }
115 
116 
117 namespace jmshower{
118 
120  {
121  return c1->NCell() > c2->NCell();
122  }
123 
125  {
126  reconfigure(pset);
127  produces< std::vector<rb::Track> >();
128 
129  }
130 
131 
133  {
134  }
135 
137  {
139 
140  fHistoGoodness = tfs->make<TH1F>("fHistoGoodness","Goodness of candidates for merging;Goodness;events", 1000,0.,1000.);
141  fHistoGoodnessLog = tfs->make<TH1F>("fHistoGoodnessLog","Goodness of candidates for merging;Log10[Goodness+1];events", 100,0.,3.);
142  fHistoEnergy = tfs->make<TH1F>("fHistoEnergy","Total Energy of Cells;PE;events", 100,10000.,20000.);
143  fHistoPurity = tfs->make<TH1F>("fHistoPurity","Purity of Clusters;(Hits of true particle/total hits in cluster);events",101,0.,1.01);
144  fHistoRemainingCells = tfs->make<TH1F>("fHistoRemainingCells","True hits not put into a cluster;Missing Cells;events",100,0.,100.);
145 
146  fTrack = tfs->make<TTree>("fTrack","fTrack");
147 
148  fTrack->Branch("trackNHit1", &m_trackNHit1);
149  fTrack->Branch("trackNHit2",&m_trackNHit2);;
150  fTrack->Branch("trackNPlane1",&m_trackNPlane1);
151  fTrack->Branch("trackNPlane2",&m_trackNPlane2);
152  fTrack->Branch("trackChisq1",&m_trackChisq1);
153  fTrack->Branch("trackChisq2",&m_trackChisq2);
154  fTrack->Branch("trackChisq",&m_trackChisq);
155  fTrack->Branch("trackDiffChisq",&m_trackDiffChisq);
156  fTrack->Branch("trackGap",&m_trackGap);
157  fTrack->Branch("trackDHit",&m_trackDHit);
158 
159  }
160 
161 
163  {
164  pSliceDir = pset.get< std::string >("SliceDir" );
165  fClusterLabel = pset.get< std::string >("ClusterLabel" );
166  fDebug = pset.get< int >("Debug" );
167  fDHitGood = pset.get< double >("DHitGood" );
168  fMinHit = pset.get< int >("MinHit" );
169  fGoodness = pset.get< double >("Goodness" );
170  fMaxPlane = pset.get< int >("MaxPlane" );
171  fMCCheater = pset.get< bool >("MCCheater");
172  }
173 
174 
175  //......................................................................
176 
178  {
180  evt.getByLabel(fClusterLabel, clusthandle);
181 
182  std::unique_ptr< std::vector<rb::Track> > trackcol(new std::vector<rb::Track> );
183 
184  std::vector< std::vector<double> > starts;
185  starts.resize(clusthandle->size());
186  std::vector< std::vector<double> > ends;
187  ends.resize(clusthandle->size());
188  std::vector< std::vector<double> > dirs;
189  dirs.resize(clusthandle->size());
190  std::vector<double> sig;
191  sig.resize(clusthandle->size());
192 
193  std::vector<bool> isnoise;
194  isnoise.resize(clusthandle->size());
195 
196  std::vector<rb::Cluster> vCluster;
197  std::vector< int > trId;
198 
199  unsigned int cc = 0;
200  double energy = 0;
201 
202  //get slices for timing information
204  evt.getByLabel(pSliceDir,slicecol);
205  std::vector< std::pair< double, double> > tMap; //contains min and max time for each slice in event, first number is min, secont max
206  for (unsigned int i=0; i<slicecol->size(); ++i) {
207  const rb::Cluster& s = (*slicecol)[i];
208  if (s.IsNoise()) continue;
209  std::pair< double, double> p(s.MinTNS(), s.MaxTNS()); //add min and max time as a pair for a given slice
210  tMap.push_back(p);
211  }
212 
213  for (unsigned int cv=0;cv<2;cv++) {
215  if (cv == 0) {
216  view = geo::kX;
217  }
218  else {
219  view = geo::kY;
220  }
221 
222  std::vector< art::Ptr<rb::Cluster> > cluster;
223 
224 
225 
226  for(unsigned int i = 0; i < clusthandle->size(); i++){
227  art::Ptr<rb::Cluster> clust(clusthandle, i);
228  if ((cv == 0) && (clust->NXCell() > 0))cluster.push_back(clust);
229  if ((cv == 1) && (clust->NYCell() > 0)) cluster.push_back(clust);
230  }
231 
232 
233  std::sort(cluster.begin(), cluster.end(), hitSort);
235  std::vector< std::vector< art::Ptr< rb::CellHit> > > cells(cluster.size());
236 
237  std::vector< std::vector<int> > motherids(cluster.size());
238  std::vector< std::vector<int> > sonids(cluster.size());
239  std::vector< unsigned int > samemotherson(cluster.size());
240 
241  for (unsigned int i=0;i<cluster.size();i++) {
242  samemotherson[i]=0;
243  for ( unsigned int nh=0;nh<cluster[i]->NCell();nh++) {
244  cells[i].push_back(cluster[i]->Cell(view,nh));
245  energy += cluster[i]->Cell(view,nh)->PE();
246  }
247  if (cells[i].size() > 1) rb::SortByPlane(cells[i]);
248  }
249 
250  for ( unsigned int nc=0;nc<cluster.size();nc++) {
251  if(cluster[nc]->NCell()==0)continue;
252  std::vector<art::Ptr< rb::CellHit> >::iterator itr = cells[nc].begin();
253  art::PtrVector<rb::CellHit> mergecells;
254  for(unsigned int i = 0; i < cluster[nc]->NCell(); i++){
255  mergecells.push_back(cluster[nc]->Cell(view,i));
256  }
257  std::vector<double> xzViewX;
258  std::vector<double> xzViewZ;
259  std::vector<double> yzViewY;
260  std::vector<double> yzViewZ;
261  std::vector<double> xzViewW; // W = weight
262  std::vector<double> yzViewW; // W = weight
263  double xyz[3]; // xyz[0]=x, xyz[1]=y, xyz[2]=z
264  int planeMax = -9999;
265  int planeMin = 9999;
266 
267  //initialize min and max time for a slice window
268  double minTNS = 0.;
269  double maxTNS = 0.;
270  //determine which slice the cluster candidate is a part of
271  for( unsigned int i = 0; i < tMap.size(); ++i){
272  if ((cluster[nc]->MinTNS() >= tMap[i].first) && (cluster[nc]->MaxTNS() <= tMap[i].second)){
273  //setting min and max time to the size of the slice the cluster belongs to
274  minTNS = tMap[i].first;
275  maxTNS = tMap[i].second;
276  break;
277  }
278  }
279 
280  while( itr != cells[nc].end() ){
281  double sig = 1.;
282  geom->Plane((*itr)->Plane())->Cell((*itr)->Cell())->GetCenter(xyz);
283  if((*itr)->View() == geo::kX){
284  if ((*itr)->Plane() > planeMax) planeMax = (*itr)->Plane();
285  if ((*itr)->Plane() < planeMin) planeMin = (*itr)->Plane();
286  xzViewX.push_back(xyz[0]);
287  xzViewZ.push_back(xyz[2]);
288  xzViewW.push_back(sig);
289  }
290  if((*itr)->View() == geo::kY){
291  if ((*itr)->Plane() > planeMax) planeMax = (*itr)->Plane();
292  if ((*itr)->Plane() < planeMin) planeMin = (*itr)->Plane();
293  yzViewY.push_back(xyz[1]);
294  yzViewZ.push_back(xyz[2]);
295  yzViewW.push_back(sig);
296  }
297  itr++;
298  }// end loop over hits in this track
299 
300  double Chi2 = 0.;
301 
302  double xzStart[2] = {0.,0.};
303  double xzEnd[2] = {0.,0.};
304  double slopeXZ = 0.;
305  double interceptXZ = 0.;
306  if (xzViewX.size() >1){
307  Chi2 = geo::LinFit(xzViewZ, xzViewX, xzViewW, &xzStart[0],
308  &xzStart[1], &xzEnd[0], &xzEnd[1]);
309  slopeXZ = 999999.;
310  interceptXZ = xzEnd[1] - slopeXZ*xzEnd[0];
311  if(fabs(xzStart[0] - xzEnd[0]) > 1.e-2){
312  slopeXZ = (xzEnd[1] - xzStart[1])/(xzEnd[0] - xzStart[0]);
313  interceptXZ = xzEnd[1] - slopeXZ*xzEnd[0];
314  }
315  }
316 
317  double yzStart[2] = {0.,0.};
318  double yzEnd[2] = {0.,0.};
319  double slopeYZ = 0.;
320  double interceptYZ = 0.;
321  if (yzViewY.size() >1){
322  Chi2 = geo::LinFit(yzViewZ, yzViewY, yzViewW, &yzStart[0],
323  &yzStart[1], &yzEnd[0], &yzEnd[1]);
324  slopeYZ = 999999.;
325  interceptYZ = yzEnd[1] - slopeYZ*yzEnd[0];
326  if(fabs(yzStart[0] - yzEnd[0]) > 1.e-2){
327  slopeYZ = (yzEnd[1] - yzStart[1])/(yzEnd[0] - yzStart[0]);
328  interceptYZ = yzEnd[1] - slopeYZ*yzEnd[0];
329  }
330  }
331 
332  ///////////////////////////////////////////////////////////////////
333  ///////////////////////////////////////////////////////////////////
334  unsigned int clustcount = cluster.size();
335 
336  for ( unsigned int numclust=0;numclust<cluster.size();numclust++) {
337  if(cluster[numclust]->NCell()==0)continue;
338  if (numclust <= nc) continue;
339 
340  //check to make sure merge candidate is part of the same time slice as cluster we are trying to add
341  if ((cluster[numclust]->MinTNS() < minTNS) || (cluster[numclust]->MaxTNS() > maxTNS)) continue;
342 
343  std::vector<double> xzViewXtest;
344  std::vector<double> xzViewZtest;
345  std::vector<double> yzViewYtest;
346  std::vector<double> yzViewZtest;
347  std::vector<double> xzViewWtest; // W = weight
348  std::vector<double> yzViewWtest; // W = weight
349  double xzStarttest[2] = {0.,0.};
350  double xzEndtest[2] = {0.,0.};
351  double yzStarttest[2] = {0.,0.};
352  double yzEndtest[2] = {0.,0.};
353  int tplaneMax = -9999;
354  int tplaneMin = 9999;
355  bool add = true;
356 
357  if (((int)cells[numclust].size() >= fMinHit) && ((int)cells[nc].size() >= fMinHit)){
358  itr = cells[numclust].begin();
359  while( itr != cells[numclust].end() ){
360  double sig = 1.;
361  geom->Plane((*itr)->Plane())->Cell((*itr)->Cell())->GetCenter(xyz);
362  if((*itr)->View() == geo::kX){
363  if ((*itr)->Plane() > tplaneMax) tplaneMax = (*itr)->Plane();
364  if ((*itr)->Plane() < tplaneMin) tplaneMin = (*itr)->Plane();
365  xzViewXtest.push_back(xyz[0]);
366  xzViewZtest.push_back(xyz[2]);
367  xzViewWtest.push_back(sig);
368  }
369  if((*itr)->View() == geo::kY){
370  if ((*itr)->Plane() > tplaneMax) tplaneMax = (*itr)->Plane();
371  if ((*itr)->Plane() < tplaneMin) tplaneMin = (*itr)->Plane();
372  yzViewYtest.push_back(xyz[1]);
373  yzViewZtest.push_back(xyz[2]);
374  yzViewWtest.push_back(sig);
375  }
376  itr++;
377  }// end loop over hits in this track
378 
379  double Chi2test = 0.;
380 
381  if (xzViewXtest.size() >1){
382  Chi2test = geo::LinFit(xzViewZtest, xzViewXtest, xzViewWtest, &xzStarttest[0], &xzStarttest[1], &xzEndtest[0], &xzEndtest[1]);
383  }
384 
385  if (yzViewYtest.size() >1){
386  Chi2test = geo::LinFit(yzViewZtest, yzViewYtest, yzViewWtest, &yzStarttest[0], &yzStarttest[1], &yzEndtest[0], &yzEndtest[1]);
387  }
388 
389  itr = cells[nc].begin();
390  while( itr != cells[nc].end() ){
391  double sig = 1.;
392  geom->Plane((*itr)->Plane())->Cell((*itr)->Cell())->GetCenter(xyz);
393  if((*itr)->View() == geo::kX){
394  xzViewXtest.push_back(xyz[0]);
395  xzViewZtest.push_back(xyz[2]);
396  xzViewWtest.push_back(sig);
397  }
398  if((*itr)->View() == geo::kY){
399  yzViewYtest.push_back(xyz[1]);
400  yzViewZtest.push_back(xyz[2]);
401  yzViewWtest.push_back(sig);
402  }
403  itr++;
404  }// end loop over hits in this track
405 
406  double Chi2merge = 0.;
407 
408  if (xzViewXtest.size() >1){
409  Chi2merge = geo::LinFit(xzViewZtest, xzViewXtest, xzViewWtest, &xzStarttest[0],
410  &xzStarttest[1], &xzEndtest[0], &xzEndtest[1]);
411  }
412 
413  if (yzViewYtest.size() >1){
414  Chi2merge = geo::LinFit(yzViewZtest, yzViewYtest, yzViewWtest, &yzStarttest[0],
415  &yzStarttest[1], &yzEndtest[0], &yzEndtest[1]);
416  }
417 
418  double diffChi2 = (Chi2merge - Chi2 - Chi2test)/(cells[numclust].size() +
419  cells[nc].size());
420 
421  if(cluster[nc]->NCell()+cluster[nc]->NCell()>=20){
422  if(diffChi2 > 7.) add = false;
423  }
424 
425  if(cluster[nc]->NCell()+cluster[nc]->NCell()<20){
426  if(diffChi2 > 30.) add = false;
427  }
428 
429  if (tplaneMax < planeMin) {
430  if (planeMin - tplaneMax > fMaxPlane) {add = false;}
431  }
432  if (tplaneMin > planeMax) {
433  if (tplaneMin - planeMax > fMaxPlane) {add = false;}
434  }
435 
436  for(unsigned int m1=0;m1<motherids[nc].size();m1++){
437  for(unsigned int m2=0;m2<motherids[numclust].size();m2++){
438  if(motherids[nc][m1]==motherids[numclust][m2]) add = false;
439  }
440  }
441 
442 
443  if (add) {
444 
445  motherids[numclust].push_back(nc);
446  sonids[nc].push_back(numclust);
447 
448  xzStart[0] = xzStarttest[0];
449  xzStart[1] = xzStarttest[1];
450  yzStart[0] = yzStarttest[0];
451  yzStart[1] = yzStarttest[1];
452  xzEnd[0] = xzEndtest[0];
453  xzEnd[1] = xzEndtest[1];
454  yzEnd[0] = yzEndtest[0];
455  yzEnd[1] = yzEndtest[1];
456 
457  for ( unsigned int nh=0;nh<cluster[numclust]->NCell();nh++) {
458  cells[nc].push_back(cluster[numclust]->Cell(view,nh));
459  mergecells.push_back(cluster[numclust]->Cell(view,nh));
460  }
461 
462  rb::SortByPlane(cells[nc]);
463  rb::SortByPlane(mergecells);
464 
465  Chi2 = Chi2merge;
466 
467 
468 
469  planeMax = -9999;
470  planeMin = 9999;
471  itr = cells[nc].begin();
472  while( itr != cells[nc].end() ){
473  if((*itr)->View() == geo::kX){
474  if ((*itr)->Plane() > planeMax) planeMax = (*itr)->Plane();
475  if ((*itr)->Plane() < planeMin) planeMin = (*itr)->Plane();
476  }
477  if((*itr)->View() == geo::kY){
478  if ((*itr)->Plane() > planeMax) planeMax = (*itr)->Plane();
479  if ((*itr)->Plane() < planeMin) planeMin = (*itr)->Plane();
480  }
481  itr++;
482  }// end loop over hits in this track
483 
484  slopeXZ = 999999.;
485  interceptXZ = xzEnd[1] - slopeXZ*xzEnd[0];
486  if(fabs(xzStart[0] - xzEnd[0]) > 1.e-2){
487  slopeXZ = (xzEnd[1] - xzStart[1])/(xzEnd[0] - xzStart[0]);
488  interceptXZ = xzEnd[1] - slopeXZ*xzEnd[0];
489  }
490 
491  slopeYZ = 999999.;
492  interceptYZ = yzEnd[1] - slopeYZ*yzEnd[0];
493  if(fabs(yzStart[0] - yzEnd[0]) > 1.e-2){
494  slopeYZ = (yzEnd[1] - yzStart[1])/(yzEnd[0] - yzStart[0]);
495  interceptYZ = yzEnd[1] - slopeYZ*yzEnd[0];
496  }
497  }
498  }
499 
500 
501  if (((int)cells[numclust].size() < fMinHit) && ((int)cells[nc].size() >= fMinHit)){
502 
503  double dHit = 0.; // distance of current hit from line
504  double slope = 0.; // slope for the hit view
505  double intercept = 0.; // intercept in the hit view
506  double point = 0.;
507  tplaneMax = -9999;
508  tplaneMin = 9999;
509  itr = cells[numclust].begin();
510  while( itr != cells[numclust].end() ){
511  geom->Plane((*itr)->Plane())->Cell((*itr)->Cell())->GetCenter(xyz);
512  if((*itr)->View() == geo::kX){
513  if ((*itr)->Plane() > tplaneMax) tplaneMax = (*itr)->Plane();
514  if ((*itr)->Plane() < tplaneMin) tplaneMin = (*itr)->Plane();
515  slope = slopeXZ;
516  intercept = interceptXZ;
517  point = xyz[0];
518 
519  }
520  else if((*itr)->View() == geo::kY){
521  if ((*itr)->Plane() > tplaneMax) tplaneMax = (*itr)->Plane();
522  if ((*itr)->Plane() < tplaneMin) tplaneMin = (*itr)->Plane();
523  slope = slopeYZ;
524  intercept = interceptYZ;
525  point = xyz[1];
526  }
527  dHit += fabs(-slope*xyz[2] + point - intercept)/sqrt(slope*slope + 1.);
528  itr++;
529  }
530  dHit = dHit / cells[numclust].size();
531 
532  if(fabs(cluster[numclust]->MinPlane()-cluster[nc]->MinPlane())<8){
533  if(dHit > 30.) add = false;
534  }
535  if(fabs(cluster[numclust]->MinPlane()-cluster[nc]->MinPlane())>=8){
536  if(dHit > 30.) add = false;
537  }
538 
539  if (tplaneMax < planeMin) {
540  if (planeMin - tplaneMax > fMaxPlane) {add = false;}
541  }
542  if (tplaneMin > planeMax) {
543  if (tplaneMin - planeMax > fMaxPlane) {add = false;}
544  }
545 
546  for(unsigned int m1=0;m1<motherids[nc].size();m1++){
547  for(unsigned int m2=0;m2<motherids[numclust].size();m2++){
548  if(motherids[nc][m1]==motherids[numclust][m2]) add = false;
549  }
550  }
551 
552 
553  if (add){
554 
555  motherids[numclust].push_back(nc);
556  sonids[nc].push_back(numclust);
557 
558  for ( unsigned int nh=0;nh<cluster[numclust]->NCell();nh++) {
559  cells[nc].push_back(cluster[numclust]->Cell(view,nh));
560  mergecells.push_back(cluster[numclust]->Cell(view,nh));
561  }
562 
563  rb::SortByPlane(cells[nc]);
564  rb::SortByPlane(mergecells);
565 
566  planeMax = -9999;
567  planeMin = 9999;
568  itr = cells[nc].begin();
569  while( itr != cells[nc].end() ){
570  double sig = 1.;
571  geom->Plane((*itr)->Plane())->Cell((*itr)->Cell())->GetCenter(xyz);
572 
573  if((*itr)->View() == geo::kX){
574  if ((*itr)->Plane() > planeMax) planeMax = (*itr)->Plane();
575  if ((*itr)->Plane() < planeMin) planeMin = (*itr)->Plane();
576  xzViewXtest.push_back(xyz[0]);
577  xzViewZtest.push_back(xyz[2]);
578  xzViewWtest.push_back(sig);
579  }
580  if((*itr)->View() == geo::kY){
581  if ((*itr)->Plane() > planeMax) planeMax = (*itr)->Plane();
582  if ((*itr)->Plane() < planeMin) planeMin = (*itr)->Plane();
583  yzViewYtest.push_back(xyz[1]);
584  yzViewZtest.push_back(xyz[2]);
585  yzViewWtest.push_back(sig);
586  }
587  itr++;
588  }// end loop over hits in this track
589 
590  if (xzViewXtest.size() >1){
591  Chi2 = geo::LinFit(xzViewZtest, xzViewXtest, xzViewWtest, &xzStarttest[0], &xzStarttest[1], &xzEndtest[0], &xzEndtest[1]);
592  }
593 
594  if (yzViewYtest.size() >1){
595  Chi2 = geo::LinFit(yzViewZtest, yzViewYtest, yzViewWtest, &yzStarttest[0], &yzStarttest[1], &yzEndtest[0], &yzEndtest[1]);
596  }
597 
598  xzStart[0] = xzStarttest[0];
599  xzStart[1] = xzStarttest[1];
600  yzStart[0] = yzStarttest[0];
601  yzStart[1] = yzStarttest[1];
602  xzEnd[0] = xzEndtest[0];
603  xzEnd[1] = xzEndtest[1];
604  yzEnd[0] = yzEndtest[0];
605  yzEnd[1] = yzEndtest[1];
606 
607  // reset the end points based on startz and endz
608  slopeXZ = 999999.;
609  interceptXZ = xzEnd[1] - slopeXZ*xzEnd[0];
610  if(fabs(xzStart[0] - xzEnd[0]) > 1.e-2){
611  slopeXZ = (xzEnd[1] - xzStart[1])/(xzEnd[0] - xzStart[0]);
612  interceptXZ = xzEnd[1] - slopeXZ*xzEnd[0];
613  }
614 
615  slopeYZ = 999999.;
616  interceptYZ = yzEnd[1] - slopeYZ*yzEnd[0];
617  if(fabs(yzStart[0] - yzEnd[0]) > 1.e-2){
618  slopeYZ = (yzEnd[1] - yzStart[1])/(yzEnd[0] - yzStart[0]);
619  interceptYZ = yzEnd[1] - slopeYZ*yzEnd[0];
620  }
621  }
622 
623  if ((numclust == (cluster.size()-1)) && (clustcount != cluster.size())){
624  numclust=0;
625  clustcount = cluster.size();
626  }
627 
628  }
629 
630 
631  if ((int)cells[numclust].size() < 3 && (int)cells[nc].size() < 3 && (int)cells[numclust].size()>0 && (int)cells[nc].size()>0 ) {
632  double dHit = 0.; // distance of current hit from line
633  double point = 0.;
634  tplaneMax = -9999;
635  tplaneMin = 9999;
636  dHit = 9999;
637  itr = cells[numclust].begin();
638  while( itr != cells[numclust].end() ){
639  geom->Plane((*itr)->Plane())->Cell((*itr)->Cell())->GetCenter(xyz);
640  if((*itr)->View() == geo::kX){
641  if ((*itr)->Plane() > tplaneMax) tplaneMax = (*itr)->Plane();
642  if ((*itr)->Plane() < tplaneMin) tplaneMin = (*itr)->Plane();
643  point = xyz[0];
644  }
645  else if((*itr)->View() == geo::kY){
646  if ((*itr)->Plane() > tplaneMax) tplaneMax = (*itr)->Plane();
647  if ((*itr)->Plane() < tplaneMin) tplaneMin = (*itr)->Plane();
648  point = xyz[1];
649  }
650  for(unsigned int ii = 0; ii<cells[nc].size(); ii++){
651  double point0 = 0;
652  double z0 = 0;
653  if((*itr)->View() == geo::kX){
654  point0 = xzViewX[ii];
655  z0 = xzViewZ[ii];
656  }else{
657  point0 = yzViewY[ii];
658  z0 = yzViewZ[ii];
659  }
660 
661  double Hitd = sqrt((point-point0)*(point-point0)+(xyz[2]-z0)*(xyz[2]-z0));
662  if(Hitd<dHit)dHit = Hitd;
663  }
664  itr++;
665  }
666 
667  int trackGap = 0;
668  if (tplaneMax <= planeMin) {
669  trackGap = planeMin - tplaneMax;
670  }
671  if (tplaneMin > planeMax) {
672  trackGap = tplaneMin - planeMax;
673  }
674 
675  if (trackGap ==0 && dHit > 15){add = false;}
676  if (trackGap !=0 && dHit > 20){add = false;}
677 
678  if (add){
679 
680  motherids[numclust].push_back(nc);
681  sonids[nc].push_back(numclust);
682 
683  for ( unsigned int nh=0;nh<cluster[numclust]->NCell();nh++) {
684  cells[nc].push_back(cluster[numclust]->Cell(view,nh));
685  mergecells.push_back(cluster[numclust]->Cell(view,nh));
686  }
687 
688  rb::SortByPlane(cells[nc]);
689  rb::SortByPlane(mergecells);
690 
691  planeMax = -9999;
692  planeMin = 9999;
693 
694 
695  itr = cells[nc].begin();
696  while( itr != cells[nc].end() ){
697  double sig = 1.;
698  geom->Plane((*itr)->Plane())->Cell((*itr)->Cell())->GetCenter(xyz);
699 
700  if((*itr)->View() == geo::kX){
701  if ((*itr)->Plane() > planeMax) planeMax = (*itr)->Plane();
702  if ((*itr)->Plane() < planeMin) planeMin = (*itr)->Plane();
703  xzViewXtest.push_back(xyz[0]);
704  xzViewZtest.push_back(xyz[2]);
705  xzViewWtest.push_back(sig);
706  }
707  if((*itr)->View() == geo::kY){
708  if ((*itr)->Plane() > planeMax) planeMax = (*itr)->Plane();
709  if ((*itr)->Plane() < planeMin) planeMin = (*itr)->Plane();
710  yzViewYtest.push_back(xyz[1]);
711  yzViewZtest.push_back(xyz[2]);
712  yzViewWtest.push_back(sig);
713  }
714  itr++;
715  }// end loop over hits in this track
716 
717  if (xzViewXtest.size() >1){
718  Chi2 = geo::LinFit(xzViewZtest, xzViewXtest, xzViewWtest, &xzStarttest[0], &xzStarttest[1], &xzEndtest[0], &xzEndtest[1]);
719  }
720 
721  if (yzViewYtest.size() >1){
722  Chi2 = geo::LinFit(yzViewZtest, yzViewYtest, yzViewWtest, &yzStarttest[0], &yzStarttest[1], &yzEndtest[0], &yzEndtest[1]);
723  }
724 
725  xzStart[0] = xzStarttest[0];
726  xzStart[1] = xzStarttest[1];
727  yzStart[0] = yzStarttest[0];
728  yzStart[1] = yzStarttest[1];
729  xzEnd[0] = xzEndtest[0];
730  xzEnd[1] = xzEndtest[1];
731  yzEnd[0] = yzEndtest[0];
732  yzEnd[1] = yzEndtest[1];
733 
734  // reset the end points based on startz and endz
735  slopeXZ = 999999.;
736  interceptXZ = xzEnd[1] - slopeXZ*xzEnd[0];
737  if(fabs(xzStart[0] - xzEnd[0]) > 1.e-2){
738  slopeXZ = (xzEnd[1] - xzStart[1])/(xzEnd[0] - xzStart[0]);
739  interceptXZ = xzEnd[1] - slopeXZ*xzEnd[0];
740  }
741 
742  slopeYZ = 999999.;
743  interceptYZ = yzEnd[1] - slopeYZ*yzEnd[0];
744  if(fabs(yzStart[0] - yzEnd[0]) > 1.e-2){
745  slopeYZ = (yzEnd[1] - yzStart[1])/(yzEnd[0] - yzStart[0]);
746  interceptYZ = yzEnd[1] - slopeYZ*yzEnd[0];
747  }
748  }
749 
750  if ((numclust == (cluster.size()-1)) && (clustcount != cluster.size())){
751  numclust=0;
752  clustcount = cluster.size();
753  }
754  }
755 
756  }//end numclust
757  ////////////////////////////////////////////////////////////////////////////////////////
758  ////////////////////////////////////////////////////////////////////////////////////////
759 
760 
761  // are there hits in each view? if not, reset the start and end z positions to
762  // be based on the view with hits
763 
764  if(mergecells.size()<3&&motherids[nc].size()>0)continue;
765  if(mergecells.size()<2)continue;
766 
767  bool sontag = false;
768 
769  if(motherids[nc].size()>0)continue;
770  double startz = 0.;
771  double endz = 0.;
772  if(xzViewZ.size() < 1){
773  startz = yzStart[0];
774  endz = yzEnd[0];
775  yzStart[1] = slopeYZ*startz + interceptYZ;
776  yzEnd[1] = slopeYZ*endz + interceptYZ;
777  xzStart[1] = geom->DetHalfWidth();
778  xzEnd[1] = geom->DetHalfWidth();
779  }
780  if(yzViewZ.size() < 1){
781  startz = xzStart[0];
782  endz = xzEnd[0];
783  xzStart[1] = slopeXZ*startz + interceptXZ;
784  xzEnd[1] = slopeXZ*endz + interceptXZ;
785  yzStart[1] = geom->DetHalfHeight();
786  yzEnd[1] = geom->DetHalfHeight();
787  }
788 
789 
790 
791  double xextremes[2] = {9999., -9999.}; // first is min, second max
792  double yextremes[2] = {9999., -9999.};
793  for(unsigned int i = 0; i < xzViewX.size(); ++i){
794  if(xzViewX[i] < xextremes[0]) xextremes[0] = xzViewX[i];
795  if(xzViewX[i] > xextremes[1]) xextremes[1] = xzViewX[i];
796  }
797 
798  for(unsigned int i = 0; i < yzViewY.size(); ++i){
799  if(yzViewY[i] < yextremes[0]) yextremes[0] = yzViewY[i];
800  if(yzViewY[i] > yextremes[1]) yextremes[1] = yzViewY[i];
801  }// end loop over hits in this track
802  // check that the end points are consistent with the extremes of hits used in each view
803  if(yzStart[1] > yzEnd[1] && yextremes[1] > yzStart[1]) yzStart[1] = yextremes[1];
804  if(yzStart[1] < yzEnd[1] && yextremes[0] < yzStart[1]) yzStart[1] = yextremes[0];
805  if(yzEnd[1] > yzStart[1] && yextremes[1] > yzEnd[1] ) yzEnd[1] = yextremes[1];
806  if(yzEnd[1] < yzStart[1] && yextremes[0] < yzEnd[1] ) yzEnd[1] = yextremes[0];
807  if(xzStart[1] > xzEnd[1] && xextremes[1] > xzStart[1]) xzStart[1] = xextremes[1];
808  if(xzStart[1] < xzEnd[1] && xextremes[0] < xzStart[1]) xzStart[1] = xextremes[0];
809  if(xzEnd[1] > xzStart[1] && xextremes[1] > xzEnd[1] ) xzEnd[1] = xextremes[1];
810  if(xzEnd[1] < xzStart[1] && xextremes[0] < xzEnd[1] ) xzEnd[1] = xextremes[0];
811 
812  //check that the end points are within the detector
813  if(xzStart[1] < -geom->DetHalfWidth()) xzStart[1] = -geom->DetHalfWidth();
814  if(xzStart[1] > geom->DetHalfWidth()) xzStart[1] = geom->DetHalfWidth();
815  if(xzEnd[1] < -geom->DetHalfWidth()) xzEnd[1] = -geom->DetHalfWidth();
816  if(xzEnd[1] > geom->DetHalfWidth()) xzEnd[1] = geom->DetHalfWidth();
817 
818  if(yzStart[1] < -geom->DetHalfHeight()) yzStart[1] = -geom->DetHalfHeight();
819  if(yzStart[1] > geom->DetHalfHeight()) yzStart[1] = geom->DetHalfHeight();
820  if(yzEnd[1] < -geom->DetHalfHeight()) yzEnd[1] = -geom->DetHalfHeight();
821  if(yzEnd[1] > geom->DetHalfHeight()) yzEnd[1] = geom->DetHalfHeight();
822 
823  double length = util::pythag(endz-startz, xzEnd[1]-xzStart[1], yzEnd[1]-yzStart[1]);
824  double dcosx = (xzEnd[1] - xzStart[1])/length;
825  double dcosy = (yzEnd[1] - yzStart[1])/length;
826  double dcosz = (endz - startz)/length;
827 
828  double evis = 0;
829  for(unsigned int i=0; i<mergecells.size(); i++){
830  evis += mergecells[i]->PE();
831  }
832  rb::Cluster Clust(mergecells);
833  vCluster.push_back(Clust);
834  starts[cc].push_back(xzStart[1]);
835  starts[cc].push_back(yzStart[1]);
836  starts[cc].push_back(startz);
837  ends[cc].push_back(xzEnd[1]);
838  ends[cc].push_back(yzEnd[1]);
839  ends[cc].push_back(endz);
840  dirs[cc].push_back(dcosx);
841  dirs[cc].push_back(dcosy);
842  dirs[cc].push_back(dcosz);
843  sig[cc] = evis;
844  isnoise[cc] = sontag;
845  cc++;
846  std::vector< int > trackcount;
847 
848  }
849  }
850 
851  std::unique_ptr< std::vector<rb::Cluster> > clustcol(new std::vector<rb::Cluster> );
852  for(unsigned int i = 0; i < vCluster.size(); i++)
853  {
854  TVector3 start(starts[i][0], starts[i][1], starts[i][2]);
855  TVector3 end(ends[i][0], ends[i][1], ends[i][2]);
856  TVector3 dir(dirs[i][0], dirs[i][1], dirs[i][2]);
857 
858  if(fDebug > 0){
859  std::cout << "3D end points are " << std::endl;
860  start.Print();
861  end.Print();
862  }
863 
864  rb::Track track(vCluster[i], start, dir);
865  track.AppendTrajectoryPoint(end);
866  track.SetNoise(isnoise[i]);
867  if(fDebug){
868  std::cout << "track start at " << track.Start().X() << " "
869  << track.Start().Y() << " " << track.Start().Z() << std::endl;
870  std::cout << "track end at " << track.Stop().X() << " "
871  << track.Stop().Y() << " " << track.Stop().Z() << std::endl;
872  }
873 
874  trackcol->push_back(track);
875  }
876  evt.put(std::move(trackcol));
877  return;
878  }//End of produce section
879 }//end namespace
880 
881 namespace jmshower{
883  }
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double LinFit(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double *x1, double *y1, double *x2, double *y2)
Find the best-fit line to a collection of points in 2-D by minimizing the squared vertical distance f...
Definition: Geo.cxx:373
void SetNoise(bool noise)
Declare the cluster to consist of noise hits or not.
Definition: Cluster.h:161
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
static bool hitSort(art::Ptr< rb::Cluster > c1, art::Ptr< rb::Cluster > c2)
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
JMClusterMerge(fhicl::ParameterSet const &pset)
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
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
TODO.
Definition: FillPIDs.h:11
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
TH1F * fHistoGoodnessLog
Easier to read in the reco validation.
int fDebug
Controls debug output 0=none, 99=scream.
c2
Definition: demo5.py:33
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
const XML_Char * s
Definition: expat.h:262
double MinTNS() const
Definition: Cluster.cxx:482
double fDHitGood
Maximum distance from the line for a hit to be considered part of it.
length
Definition: demo0.py:21
bool fMCCheater
True/false MCCheater SimHit/SimTrack input.
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
const std::string cv[Ncv]
void AppendTrajectoryPoint(TVector3 pt)
Definition: Track.cxx:112
T get(std::string const &key) const
Definition: ParameterSet.h:231
double energy
Definition: plottest35.C:25
Collect Geo headers and supply basic geometry functions.
double MaxTNS() const
Definition: Cluster.cxx:528
size_type size() const
Definition: PtrVector.h:308
OStream cout
Definition: OStream.cxx:6
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
static constexpr Double_t m2
Definition: Munits.h:145
void produce(art::Event &evt)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) 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
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
void geom(int which=0)
Definition: geom.C:163
c1
Definition: demo5.py:24
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
Float_t e
Definition: plot.C:35
Float_t track
Definition: plot.C:35
add("abs", expr_type(int_type()), expr_type(int_type()))
void reconfigure(const fhicl::ParameterSet &pset)
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