BreakPoint_module.cc
Go to the documentation of this file.
1 ///
2 /// \file BreakPoint_module.cc
3 /// \brief Run the break-point fitter
4 /// \author messier@indiana.edu
5 /// \version $Id:$
6 ///
7 #include "TH2F.h"
8 #include "TNtuple.h"
9 //
14 //
15 #include "Geometry/Geometry.h"
17 #include "RecoBase/Vertex.h"
18 #include "RecoBase/Cluster.h"
19 #include "RecoBase/Prong.h"
20 #include "RecoBase/Track.h"
21 #include "RecoBase/FitSum.h"
22 //
23 #include "Utilities/AssociationUtil.h"
24 //
31 
32 namespace bpfit { class BreakPoint; }
33 
35 {
36 public:
38 
39  // FHICL Configuration data
40  bool fDebugHitList0; ///< Debug the detector-frame hit list?
41  bool fDebugHitList1; ///< Debug the track-frame hit list?
42  bool fDebugScatSurf; ///< Debug the scattering surfaces
43  bool fDebugFitList0; ///< Debug the fitted track locations (det. frame)
44  bool fDebugFitList1; ///< Debug the fitted track locations (track frame)
45  bool fLogResiduals; ///< Make ntuple of residuals?
46  bool fFillTrackNt; ///< Fill ntuple of tracking results?
47  bool fFillExceptNt; ///< Fill ntuple of exception cases?
48  bool fVerbose; ///< Run in verbose mode (print messages about all exceptions?)
49  bool fFitAsMuon; ///< Fit tracks under the assumption of a muon?
50  bool fFitAsPion; ///< Fit tracks under the assumption of a pion?
51  bool fFitAsProton; ///< Fit tracks under the assumption of a proton?
52  double fExtendUp; ///< Fractionally, how far to extend track into upstream dead material
53  double fExtendDn; ///< Fractionally, how far to extend track into downstream dead material
54  std::string fVertexLabel; ///< Where to find vertices in event
55  std::string fClusterLabel; ///< Where to find clusters in event
56  std::string fClusterInstance; ///< Which instance of clusters to find
57  unsigned int fNhitX; ///< Require prongs to have this many x hits
58  unsigned int fNhitY; ///< Require prongs to have this many y hits
59  unsigned int fNhit; ///< Require this many hits total
60  double fChi2Cutoff; ///< Drop hits with this chi2 or greater
61  double fdStol; ///< Scattering surface ds tolerance
62  double fdXtol; ///< Scattering surface dx tolerance
63  double fdTtol; ///< Scattering surface dt tolerance
64  double fSigAdj; ///< Multiplier to adjust sigmaXYZ for hits (see comment in the FillHitList function.) Set to 1.0 if no adjustment is desired.
65  unsigned int fMaxSliceHits; ///< Don't bother processing a slice if it has more than this many hits. This is to prevent computing problems.
66 
67 public:
68  // Member data
69  int fRun;
70  int fSubrun;
71  int fEvent;
72  unsigned int fSlice;
73  unsigned int fProng;
74  TNtuple* fHitList0Nt;
75  TNtuple* fHitList1Nt;
76  TNtuple* fScatSurfNt;
77  TNtuple* fFitList0Nt;
78  TNtuple* fFitList1Nt;
79  TNtuple* fResidNt;
80  TNtuple* fTrackNt;
81  TNtuple* fExceptNt;
82 
83  std::vector<int> fPartList; ///< List of particle types (by pdg code) to be used for track fits.
84 
85 public:
86  //......................................................................
87 
89  {
90  fDebugHitList0 = p.get<bool>("DebugHitList0");
91  fDebugHitList1 = p.get<bool>("DebugHitList1");
92  fDebugScatSurf = p.get<bool>("DebugScatSurf");
93  fDebugFitList0 = p.get<bool>("DebugFitList0");
94  fDebugFitList1 = p.get<bool>("DebugFitList1");
95  fLogResiduals = p.get<bool>("LogResiduals");
96  fFillTrackNt = p.get<bool>("FillTrackNt");
97  fFillExceptNt = p.get<bool>("FillExceptNt");
98  fVerbose = p.get<bool>("Verbose");
99 
100  fExtendUp = p.get<double>("ExtendUp");
101  fExtendDn = p.get<double>("ExtendDn");
102  fFitAsMuon = p.get<bool>("FitAsMuon");
103  fFitAsPion = p.get<bool>("FitAsPion");
104  fFitAsProton = p.get<bool>("FitAsProton");
105  fPartList.clear();
106  if(fFitAsMuon) fPartList.push_back(13);
107  if(fFitAsPion) fPartList.push_back(211);
108  if(fFitAsProton) fPartList.push_back(2212);
109 
110  fVertexLabel = p.get<std::string>("VertexLabel");
111  fClusterLabel = p.get<std::string>("ClusterLabel");
112  fClusterInstance = p.get<std::string>("ClusterInstance");
113 
114  fNhitX = p.get<unsigned int>("NhitX");
115  fNhitY = p.get<unsigned int>("NhitY");
116  fNhit = p.get<unsigned int>("Nhit");
117  fChi2Cutoff = p.get<double>("Chi2Cutoff");
118  fdStol = p.get<double>("dStol");
119  fdXtol = p.get<double>("dXtol");
120  fdTtol = p.get<double>("dTtol");
121  fSigAdj = p.get<double>("SigmaAdjust");
122  fMaxSliceHits = p.get<unsigned int>("MaxSliceHits");
123  }
124 
125  //......................................................................
126 
127  explicit BreakPoint(fhicl::ParameterSet const& p) :
128  fSliceToken(consumes<std::vector<rb::Cluster>>(p.get<std::string>("SliceLabel"))),
129  fRun(0),
130  fSubrun(0),
131  fEvent(0),
132  fProng(0),
133  fHitList0Nt(0),
134  fHitList1Nt(0),
135  fScatSurfNt(0),
136  fFitList0Nt(0),
137  fFitList1Nt(0),
138  fResidNt(0),
139  fTrackNt(0),
140  fExceptNt(0)
141  {
142  produces< std::vector<rb::Track> >();
143  produces< std::vector<rb::FitSum> >();
144  produces< art::Assns<rb::Track, rb::Prong> >();
145  produces< art::Assns<rb::Track, rb::Cluster> >();
146  produces< art::Assns<rb::FitSum, rb::Track> >();
147  this->reconfigure(p);
148  }
149 
150  //......................................................................
151 
152  virtual ~BreakPoint() { }
153 
154  //......................................................................
155 
156  void FindSlices(const art::Event& evt,
157  art::Handle< std::vector<rb::Cluster> >& h,
159  {
160  unsigned int i;
161  evt.getByToken(fSliceToken, h);
162  for (i=0; i<h->size(); ++i) {
164  }
165  }
166 
167  //......................................................................
168 
170  unsigned int islice,
172  std::vector< art::Ptr<rb::Vertex> >& vertexp,
173  std::vector< art::Ptr<rb::Prong> >& prongp)
174  {
175  if (prongh.isValid()) {
176  prongp = prongh.at(islice);
177  }
178  else {
179  //
180  // Hmm... No prongs associated with slice - try the vertex.
181  //
182  art::FindManyP<rb::Prong> prongh2(vertexp, evt, fClusterLabel);
183  if (prongh2.isValid()) {
184  if (prongh2.size()==1) {
185  prongp = prongh2.at(0);
186  }
187  else {
188  std::cerr << __FILE__ << ":" << __LINE__
189  << " Multiple vertices breaks our assumptions."
190  << std::endl;
191  }
192  }
193  }
194  }
195 
196  //......................................................................
197 
199  {
200  if (t->Is2D()) return false;
201 
202  unsigned int nx = t->NXCell();
203  unsigned int ny = t->NYCell();
204  if (nx<fNhitX || ny<fNhitY) return false;
205  if (nx+ny<fNhit) return false;
206 
207  return true;
208  }
209 
210  //......................................................................
211 
214  {
216  unsigned int itr, i;
217  unsigned int n = prongp->NCell();
218  //
219  // Build the 3D hits in two passes. In the first we don't have
220  // information about where transversly we are along a cell. In the
221  // second iteration, we do.
222  //
223  for (itr=0; itr<2; ++itr) {
224  for (i=0; i<n; ++i) {
225  const rb::CellHit* h = prongp->Cell(i).get();
226 
227  unsigned int ip = h->Plane();
228  unsigned int ic = h->Cell();
229 
230  const geo::PlaneGeo* p = geom->Plane(ip);
231  const geo::CellGeo* c = p-> Cell(ic);
232 
233  double xyz[3] = {0,0,0};;
234  double sigx=0, sigy=0, sigz=0, sigt=50;
235  double w = 1.0;
236  if (itr==0) {
237  //
238  // On first iteration we don't have information from other
239  // view. Use cell centers for position estimates.
240  //
241 
242  // NOTE: BPF appears to be giving too much weight to the hits
243  // when doing the track fit (the total chi2 distribution
244  // does not have a mean of 1.0.) I believe this to be
245  // related to the fact that the true particle location
246  // within a cell is not Gaussian distributed around the
247  // cell center. Introducing the fSigAdj factor adjusts
248  // this and brings the mean of the chi2 distribution down
249  // to 1.0 allowing for the fit scattering angles to be better.
250  c->GetCenter(xyz);
251  if (p->View()==geo::kX) sigx = c->HalfW()*fSigAdj;
252  if (p->View()==geo::kY) sigy = c->HalfW()*fSigAdj;
253  sigz = c->HalfD()*fSigAdj;
254  }
255  else {
256  //
257  // On second iteration, we have more information about where
258  // in the cell we are. Use the estimated location from the
259  // other view. Otherwise, most information is a copy of what
260  // is already there.
261  //
262  if (p->View()==geo::kX) c->GetCenter(xyz, hits[i].fY);
263  if (p->View()==geo::kY) c->GetCenter(xyz, hits[i].fX);
264  sigx = hits[i].fSigX;
265  sigy = hits[i].fSigY;
266  sigz = hits[i].fSigZ;
267  sigt = hits[i].fSigT;
268  w = hits[i].fW;
269  }
270 
271  if (p->View()==geo::kX) {
272  hits.SetHit(i, geo::kX, xyz[0], xyz[1], xyz[2], h->TNS(),
273  sigx, sigy, sigz, sigt, w);
274  }
275  else if (p->View()==geo::kY) {
276  hits.SetHit(i,
277  geo::kY, xyz[0], xyz[1], xyz[2], h->TNS(),
278  sigx, sigy, sigz, sigt, w);
279  }
280  else abort();
281  } // loop on hits
282  hits.CompleteOrtho();
283  } // iterations
284  //
285  // Keep all hits inside the detector
286  //
287  double xyz[3];
288  for (i=0; i<hits.size(); ++i) {
289  xyz[0] = hits[i].fX;
290  xyz[1] = hits[i].fY;
291  xyz[2] = hits[i].fZ;
292  this->ClampToDetector(xyz);
293  hits[i].fX = xyz[0];
294  hits[i].fY = xyz[1];
295  hits[i].fZ = xyz[2];
296  }
297 
298  hits.SortByZ();
299  }
300 
301  //......................................................................
302  ///
303  /// If there are less than 10 hits in one view with the same z then
304  /// take the average x or y depending on view and use that when
305  /// fitting the track.
306  ///
308  {
309 
310  // NOTE: To the reader of this code - I think there is a mistake below.
311  // I don't think the intention is to set the weight of any hits
312  // to 1e-9. I think that the weights should be left as 1.0. But
313  // we aren't currently using this function...
314 
315  for(unsigned int i=0;i<h.size();i++){
316 
317  int lastj = i;
318  int n_xhits = 0;
319  int n_yhits = 0;
320  unsigned int i_xhits[384], i_yhits[384]; // Max number of hits at the same z = # cells per plane
321  double sumx = 0.0;
322  double sumy = 0.0;
323  double sigx2 = 0.0;
324  double sigy2 = 0.0;
325 
326  for(unsigned int j=i+1;j<h.size();j++){
327  // Look for hits in the same view that have the
328  // same z coordinate.
329  if((fabs(h[i].fZ-h[j].fZ)<0.1)){
330  if(h[i].fView == h[j].fView){
331  lastj = j;
332  if(h[i].fView == geo::kX){
333  if(n_xhits==0) i_xhits[n_xhits] = i;
334  else i_xhits[n_xhits] = j;
335  n_xhits++;
336  }
337  else{
338  if(h[i].fView == geo::kY){
339  if(n_yhits==0) i_yhits[n_yhits] = i;
340  else i_yhits[n_yhits] = j;
341  n_yhits++;
342  }
343  }
344  }
345  }
346  else continue;
347  } // end j loop over hitlist3D
348 
349  if(n_xhits>0 && n_xhits<10){
350  // Averaging x
351  for(int xi=0; xi<n_xhits; xi++){
352  sumx += h[i_xhits[xi]].fX;
353  sigx2 += h[i_xhits[xi]].fSigX*h[i_xhits[xi]].fSigX;
354 
355  // Since we will be adding the averaged hit to
356  // the hitlist we need to reweight the same z hits
357  // to something very small.
358  if(xi!=0)h[i_xhits[xi]].fW = 1.0E-9;
359  }
360  double xave = sumx/n_xhits;
361  // Replace first hit in series that has the same z with
362  // new average x coordinate.
363  h.SetHit(i_xhits[0], h[i_xhits[0]].fView, xave, h[i_xhits[0]].fY, h[i_xhits[0]].fZ,
364  h[i_xhits[0]].fT, sqrt(sigx2), h[i_xhits[0]].fSigY,
365  h[i_xhits[0]].fSigZ, h[i_xhits[0]].fSigT, h[i_xhits[0]].fW);
366  } // end n_xhits requirement
367 
368  if(n_yhits>0 && n_yhits<10){
369  // Averaging y
370  for(int yi=0; yi<n_yhits; yi++){
371  sumy += h[i_yhits[yi]].fY;
372  sigy2 += h[i_yhits[yi]].fSigY*h[i_yhits[yi]].fSigY;
373  // Since we will be adding the averaged hit to
374  // the hitlist we need to reweight the same z hits
375  // to something very small.
376  if(yi!=0) h[i_yhits[yi]].fW = 1.0E-9;
377  }
378  double yave = sumy/n_yhits;
379  // Replace first hit in series that has the same z with
380  // new average y coordinate.
381  h.SetHit(i_yhits[0], h[i_yhits[0]].fView, h[i_yhits[0]].fX, yave, h[i_yhits[0]].fZ,
382  h[i_yhits[0]].fT, h[i_yhits[0]].fSigX, sqrt(sigy2),
383  h[i_yhits[0]].fSigZ, h[i_yhits[0]].fSigT, h[i_yhits[0]].fW);
384  } // end n_yhits requirement
385 
386  i = lastj;
387  } // end i loop over hitlist3D
388 
389  h.CompleteOrtho();
390 
391  h.SortByZ();
392 
393  } // end function definition
394 
395  //......................................................................
396  ///
397  /// Dump the hit list to an ntuple we can study offline for
398  /// debugging
399  ///
400  void DebugHitList(bpfit::HitList3D& h, int which)
401  {
402  TNtuple* nt = 0;
403 
404  if (which==0) {
405  if (fHitList0Nt==0) {
407  fHitList0Nt = tfs->make<TNtuple>("hits0",
408  "Det Frame Hit List Debugging",
409  "run:subrun:evt:trk:v:x:y:z:sx:sy:sz:slice");
410  }
411  nt = fHitList0Nt;
412  }
413  if (which==1) {
414  if (fHitList1Nt==0) {
416  fHitList1Nt = tfs->make<TNtuple>("hits1",
417  "Trk Frame Hit List Debugging",
418  "run:subrun:evt:trk:v:x:y:z:sx:sy:sz:slice");
419  }
420  nt = fHitList1Nt;
421  }
422  float x[12];
423  x[0] = fRun;
424  x[1] = fSubrun;
425  x[2] = fEvent;
426  x[3] = fProng;
427  x[11]= fSlice;
428  for (unsigned int i=0; i<h.size(); ++i) {
429  x[4] = h[i].fView;
430  x[5] = h[i].fX;
431  x[6] = h[i].fY;
432  x[7] = h[i].fZ;
433  x[8] = h[i].fSigX;
434  x[9] = h[i].fSigY;
435  x[10] = h[i].fSigZ;
436  nt->Fill(x);
437  }
438  }
439 
440  //......................................................................
441 
442  void FillTrackNt(const HitList3D& h,
443  const Lutz& lx,
444  const Lutz& ly,
445  const Path& path,
446  int pdg_id)
447  {
448  if (fTrackNt==0) {
450  fTrackNt = tfs->make<TNtuple>("trks",
451  "Track summary",
452  "run:subrun:evt:slice:trk:"
453  "nhitx:nhity:x0:y0:z0:x1:y1:z1:length:p:"
454  "theta:phi:"
455  "chi2xh:chi2xa:chi2yh:chi2ya:pdgID");
456  }
457  unsigned int i;
458  unsigned int nhitx = 0, nhity = 0;
459  for (i=0; i<h.size(); ++i) {
460  if (h[i].fView==geo::kX && h[i].fW > 1.0E-9) ++nhitx;
461  if (h[i].fView==geo::kY && h[i].fW > 1.0E-9) ++nhity;
462  }
463 
464  double p, len;
465  len = path.Length();
466  if (pdg_id == 13) path.MuonParams(0, 0, &p, 0, 0);
467  else if (pdg_id == 211) path.PionParams(0, 0, &p, 0, 0);
468  else if (pdg_id == 2212) path.ProtonParams(0, 0, &p, 0, 0);
469  // Can't get here unless we are one of the above particle types,
470  // but why leave traps...?
471  else path.MuonParams(0, 0, &p, 0, 0);
472  i = h.size()-1;
473  double dx = path.fXYZ[1].X() - path.fXYZ[0].X();
474  double dy = path.fXYZ[1].Y() - path.fXYZ[0].Y();
475  double dz = path.fXYZ[1].Z() - path.fXYZ[0].Z();
476  if (dz==0.0) dz = 1.0E-9;
477 
478  float x[22];
479  x[0] = fRun;
480  x[1] = fSubrun;
481  x[2] = fEvent;
482  x[3] = fSlice;
483  x[4] = fProng;
484  x[5] = nhitx;
485  x[6] = nhity;
486  x[7] = h[0].fX;
487  x[8] = h[0].fY;
488  x[9] = h[0].fZ;
489  x[10] = h[i].fX;
490  x[11] = h[i].fY;
491  x[12] = h[i].fZ;
492  x[13] = len;
493  x[14] = p;
494  x[15] = 180.0*acos(dz/sqrt(dx*dx+dy*dy+dz*dz))/M_PI;
495  x[16] = 180.0*atan2(dy,dx)/M_PI;
496  x[17] = lx.Chi2XI();
497  x[18] = lx.Chi2Beta();
498  x[19] = ly.Chi2XI();
499  x[20] = ly.Chi2Beta();
500  x[21] = pdg_id;
501  fTrackNt->Fill(x);
502  }
503 
504  //......................................................................
505 
507  {
508 
509  if (fExceptNt == 0) {
511  fExceptNt = tfs->make<TNtuple>("except",
512  "Exception summary",
513  "run:subrun:evt:slice:prong:flag:value");
514  }
515 
516  float x[7];
517  x[0] = fRun;
518  x[1] = fSubrun;
519  x[2] = fEvent;
520  x[3] = fSlice;
521  x[4] = fProng;
522  x[5] = Ex.fFlag;
523  x[6] = Ex.fValue;
524  fExceptNt->Fill(x);
525 
526  }
527 
528  //......................................................................
529 
531  const HitList3D& hits,
532  bpfit::TrackBasis& sys)
533  {
534  //
535  // Extract hits into format expected by TrackBasis
536  //
537  unsigned int i;
538  std::vector<double> x(hits.size());
539  std::vector<double> y(hits.size());
540  std::vector<double> z(hits.size());
541  std::vector<double> dx(hits.size());
542  std::vector<double> dy(hits.size());
543  std::vector<double> dz(hits.size());
544 
545  for (i=0; i<hits.size(); ++i) {
546  x[i] = hits[i].fX;
547  y[i] = hits[i].fY;
548  z[i] = hits[i].fZ;
549  dx[i] = hits[i].fSigX;
550  dy[i] = hits[i].fSigY;
551  dz[i] = hits[i].fSigZ;
552  }
553 
554  sys.MakeBasis(x,y,z,dx,dy,dz,vtx->GetX(),vtx->GetY(),vtx->GetZ());
555  }
556 
557  //......................................................................
558 
559  void ClampToDetector(double* p)
560  {
562  if (p[0]>geom->DetHalfWidth()) p[0] = geom->DetHalfWidth();
563  if (p[1]>geom->DetHalfHeight()) p[1] = geom->DetHalfHeight();
564  if (p[2]>geom->DetLength()) p[2] = geom->DetLength();
565 
566  if (p[0]<-geom->DetHalfWidth()) p[0] = -geom->DetHalfWidth();
567  if (p[1]<-geom->DetHalfHeight()) p[1] = -geom->DetHalfHeight();
568  if (p[2]<0) p[2] = 0;
569  }
570 
571  //......................................................................
572 
574  bpfit::TrackBasis& sys,
575  double* p0,
576  double* p1)
577  {
578  unsigned int i;
579  double zmin = 9e9;
580  double zmax = -9e9;
581  for (i=0; i<h.size(); ++i) {
582  if (h[i].fZ<zmin) zmin = h[i].fZ;
583  if (h[i].fZ>zmax) zmax = h[i].fZ;
584  }
585  p0[0] = 0.0; p0[1] = 0.0; p0[2] = zmin;
586  p1[0] = 0.0; p1[1] = 0.0; p1[2] = zmax;
587  double dx;
588  double dy;
589  double dz;
590  sys.TrkToDet(0, 0, zmin, 0, 0, 0,
591  &p0[0], &p0[1], &p0[2], &dx, &dy, &dz);
592  sys.TrkToDet(0, 0, zmax, 0, 0, 0,
593  &p1[0], &p1[1], &p1[2], &dx, &dy, &dz);
594  this->ClampToDetector(p0);
595  this->ClampToDetector(p1);
596  }
597 
598  //......................................................................
599 
601  {
602  if (fScatSurfNt==0) {
604  fScatSurfNt = tfs->make<TNtuple>("scat",
605  "Scattering Surfaces",
606  "run:subrun:evt:trk:s:sigalpha:pdgid:slice");
607  }
608  float x[8];
609  x[0] = fRun;
610  x[1] = fSubrun;
611  x[2] = fEvent;
612  x[3] = fProng;
613  x[7] = fSlice;
614  for (unsigned int i=0; i<s.N(); ++i) {
615  x[4] = s.S(i);
616  x[5] = s.SigAlpha(i);
617  x[6] = pdgid;
618  fScatSurfNt->Fill(x);
619  }
620  }
621 
622  //......................................................................
623 
624  void DebugFitList(unsigned int which,
625  const std::vector<double>& x,
626  const std::vector<double>& y,
627  const std::vector<double>& z)
628  {
629  TNtuple* nt = 0;
630  if (which==0) {
631  if (fFitList0Nt==0) {
633  fFitList0Nt = tfs->make<TNtuple>("fit0",
634  "fitted track - det. frame",
635  "run:subrun:evt:trk:x:y:z:slice");
636  }
637  nt = fFitList0Nt;
638  }
639  if (which==1) {
640  if (fFitList1Nt==0) {
642  fFitList1Nt = tfs->make<TNtuple>("fit1",
643  "fitted track - track frame",
644  "run:subrun:evt:trk:x:y:z:slice");
645  }
646  nt = fFitList1Nt;
647  }
648 
649  unsigned int i;
650  float v[8];
651  v[0] = fRun;
652  v[1] = fSubrun;
653  v[2] = fEvent;
654  v[3] = fProng;
655  v[7] = fSlice;
656  for (i=0; i<z.size(); ++i) {
657  v[4] = x[i];
658  v[5] = y[i];
659  v[6] = z[i];
660  nt->Fill(v);
661  }
662  }
663 
664  //......................................................................
665  //
666  // Create an ntuple so we can examine the residuals. Variables are:
667  // m - measured value. x/y/alpha in the two views
668  // sig - estimated 1-sigma error on measured value
669  // fit - fitted value of parameter
670  // w - which kind of measurement is it. Code follows:
671  // 0 = X measurement from x view
672  // 1 = X measurement interpolated
673  // 2 = Y measurement from the y view
674  // 3 = Y measurement interpolated
675  // 4 = angle "measurement" from the x view
676  // 5 = angle "measurement" from the y view
677  // zs - for hits, the z value of the hit in the track frame,
678  // for angles, the s value along the track
679  //
682  bpfit::Lutz& xfit,
683  bpfit::Lutz& yfit,
684  double* alphax,
685  double* alphay,
686  float pdgid)
687  {
688  if (fResidNt==0) {
690  fResidNt = tfs->make<TNtuple>("resid",
691  "resid",
692  "run:subrun:evt:trk:m:sig:fit:w:pdgid:zs:slice");
693  }
694  unsigned int i;
695  float v[11];
696  v[0] = fRun;
697  v[1] = fSubrun;
698  v[2] = fEvent;
699  v[3] = fProng;
700  v[10]= fSlice;
701  v[8] = pdgid;
702  for (i=0; i<hit.size(); ++i) {
703  v[4] = hit[i].fX;
704  v[5] = hit[i].fSigX;
705  v[6] = xfit.X(hit[i].fZ);
706  v[9] = hit[i].fZ;
707  if (hit[i].fView==geo::kX) v[7] = 0;
708  if (hit[i].fView==geo::kY) v[7] = 1;
709  fResidNt->Fill(v);
710 
711  v[4] = hit[i].fY;
712  v[5] = hit[i].fSigY;
713  v[6] = yfit.X(hit[i].fZ);
714  v[9] = hit[i].fZ;
715  if (hit[i].fView==geo::kY) v[7] = 2;
716  if (hit[i].fView==geo::kX) v[7] = 3;
717  fResidNt->Fill(v);
718  }
719  for (i=0; i<scat.N(); ++i) {
720  v[4] = 0;
721  v[5] = scat.SigAlpha(i);
722  v[6] = alphax[i];
723  v[7] = 4;
724  v[9] = scat.S(i);
725  fResidNt->Fill(v);
726 
727  v[4] = 0;
728  v[5] = scat.SigAlpha(i);
729  v[6] = alphay[i];
730  v[7] = 5;
731  v[9] = scat.S(i);
732  fResidNt->Fill(v);
733  }
734 
735  }
736 
737  //......................................................................
738 
740  art::Ptr<rb::Prong>& prongp,
741  std::vector<rb::Track>& tracks,
742  std::vector<rb::FitSum>& fitsums)
743  {
744  unsigned int i, j;
745  bpfit::TrackBasis sys;
746  double p0[3], p1[3];
747 
748  bpfit::HitList3D hits0(prongp->NCell());
749  bpfit::HitList3D hits1(prongp->NCell());
750  std::vector<double> s(hits1.size());
751 
752  // This first try/catch block is to pick out hit list errors. These
753  // will cause all particle track fit assumptions to fail.
754  try {
755  // First check that this prong has hits in more than one plane in
756  // each view. If not, we won't try to fit it.
758  std::set<unsigned int> Xplanes;
759  std::set<unsigned int> Yplanes;
760  for(unsigned int p = 0; p < prongp->NCell(); ++p) {
761  unsigned int plane = prongp->Cell(p).get()->Plane();
762  if(geom->Plane(plane)->View() == geo::kX) Xplanes.insert(plane);
763  if(geom->Plane(plane)->View() == geo::kY) Yplanes.insert(plane);
764  }
765  if(Xplanes.size() == 1 && Yplanes.size() == 1) {
766  throw BPException(__FILE__, __LINE__, kNPLANES, 2);
767  }
768  if(Xplanes.size() == 1) {
769  throw BPException(__FILE__, __LINE__, kNPLANES, 0);
770  }
771  if(Yplanes.size() == 1) {
772  throw BPException(__FILE__, __LINE__, kNPLANES, 1);
773  }
774 
775  //
776  // Track fitting starts with a hit list
777  //
778  this->FillHitList(prongp, hits0);
779 
780  /*
781  // Only average hits if either view has a small number of hits.
782  if(prongp->NXCell()<10 || prongp->NYCell()<10){
783  // Currently commented out the use of the averagexy function...
784  // I have verified that for most fits nothing is affected but there are still
785  // a few I need to dive deeper into.
786 
787  this->AverageXY(hits0);
788  }
789  */
790 
791  for(i=0; i<hits0.size(); ++i){
792  hits1.SetHit(i,hits0[i].fView, hits0[i].fX, hits0[i].fY, hits0[i].fZ, hits0[i].fT,
793  hits0[i].fSigX, hits0[i].fSigY, hits0[i].fSigZ, hits0[i].fSigT,
794  hits0[i].fW);
795  }
796 
797  if (fDebugHitList0) this->DebugHitList(hits0, 0);
798  //
799  // Then we want to swing around to a track-based coordinate system
800  // so that we can properly locate the scattering planes
801  //
802  this->MakeBasis(vertexp, hits0, sys);
803 
804  for (i=0; i<hits0.size(); ++i) {
805  sys.DetToTrk(hits0[i].fX, hits0[i].fY, hits0[i].fZ,
806  hits0[i].fSigX, hits0[i].fSigY, hits0[i].fSigZ,
807  &hits1[i].fX, &hits1[i].fY, &hits1[i].fZ,
808  &hits1[i].fSigX, &hits1[i].fSigY, &hits1[i].fSigZ);
809  }
810  if (fDebugHitList1) this->DebugHitList(hits1, 1);
811  hits1.SortByZ();
812 
813  //
814  // Get the extent of the track and z1 locations of measurement
815  // planes
816  //
817  this->FindEndPoints(hits1, sys, p0, p1);
818 
819  for (i=0; i<hits1.size(); ++i) {
820  s[i] = hits1[i].fZ;
821  }
822 
823  } // end of try block
824  catch (BPException& Ex) {
825  if(fVerbose) Ex.Print();
826  if(fFillExceptNt) this->FillExceptNt(Ex);
827  return;
828  } // end of catch block
829 
830  //
831  // Loop over the list of particle assumptions requested for this fit,
832  // and make one track for each of those assumptions.
833  //
834  for(unsigned int t = 0; t < fPartList.size(); ++t) {
835  // This second try/catch block is to pick out track fitting errors.
836  // This will cause one particle fit assumption to fail but not all.
837  try {
838 
839  // Make the track to be fit.
840  //
841  ///\ TODO: For now we initialize it to the prong that it is fit from.
842  // We should be putting in the hits from the prongs with the
843  // weights used to calculate the track.
844  rb::Track track(*prongp);
845  rb::FitSum fitsum;
846 
848  bpfit::ScatteringSurfaces scat_surf(fPartList[t], fdStol, fdXtol, fdTtol);
849  scat_surf.MakeSurfaces(*geom, p0, p1, s);
850 
851  if (fDebugScatSurf) this->DebugScatSurf(scat_surf, fPartList[t]);
852 
853  bpfit::Lutz lutzx(hits1.size(), scat_surf.N());
854  bpfit::Lutz lutzy(hits1.size(), scat_surf.N());
855  for (i=0; i<hits1.size(); ++i) {
856  lutzx.SetMeasurement(i, hits1[i].fZ, hits1[i].fX, hits1[i].fSigX);
857  lutzy.SetMeasurement(i, hits1[i].fZ, hits1[i].fY, hits1[i].fSigY);
858  }
859  for (j=0; j<scat_surf.N(); ++j) {
860  lutzx.SetScatteringPlane(j, scat_surf.S(j), scat_surf.SigAlpha(j));
861  lutzy.SetScatteringPlane(j, scat_surf.S(j), scat_surf.SigAlpha(j));
862  }
863  double alphax[scat_surf.N()];
864  double alphay[scat_surf.N()];
865  lutzx.Fit(0, 0, alphax, 0);
866  lutzy.Fit(0, 0, alphay, 0);
867  //
868  // Check the fit results and remove far-away outliers
869  //
870  double chix, chiy, chim, fac;
871  for (i=0; i<hits1.size(); ++i) {
872  chix = lutzx.Chi2XIi(i);
873  chiy = lutzy.Chi2XIi(i);
874  chim = (chix>chiy) ? chix : chiy;
875  fac = sqrt(chim);
876  if (chim>fChi2Cutoff) {
877  hits1[i].fW = 1.0E-9;
878  hits1[i].fSigX *= fac;
879  hits1[i].fSigY *= fac;
880  hits1[i].fSigZ *= fac;
881  }
882  }
883  //
884  // Need to recompute the projection from orthogonal projection as
885  // outliers would have biased the estimates
886  //
887  hits1.CompleteOrtho();
888  for (i=0; i<hits1.size(); ++i) {
889  lutzx.SetMeasurement(i, hits1[i].fZ, hits1[i].fX, hits1[i].fSigX);
890  lutzy.SetMeasurement(i, hits1[i].fZ, hits1[i].fY, hits1[i].fSigY);
891  }
892  lutzx.Fit(0, 0, alphax, 0);
893  lutzy.Fit(0, 0, alphay, 0);
894 
895  //
896  // Get the fitted track positions in the track frame. Take the
897  // start and end points from the hits and insert the track
898  // locations from the scattering surfaces.
899  //
900  std::vector<double> xfit(scat_surf.N()+2);
901  std::vector<double> yfit(scat_surf.N()+2);
902  std::vector<double> zfit(scat_surf.N()+2);
903  zfit[0] = hits1[0].fZ;
904  for (j=0; j<scat_surf.N(); ++j) {
905  zfit[1+j] = scat_surf.S(j);
906  }
907  zfit[1+j] = hits1[hits1.size()-1].fZ;
908  for (j=0; j<zfit.size(); ++j) {
909  xfit[j] = lutzx.X(zfit[j]);
910  yfit[j] = lutzy.X(zfit[j]);
911  }
912 
913  if (fDebugFitList1) this->DebugFitList(1, xfit, yfit, zfit);
914  if (fLogResiduals) {
915  this->LogResiduals(hits1, scat_surf, lutzx, lutzy, alphax, alphay, fPartList[t]);
916  }
917 
918  //
919  // Translate back to the detector frame and build a path
920  //
921  std::vector<double> xfit_det;
922  std::vector<double> yfit_det;
923  std::vector<double> zfit_det;
924  double xyz[3], dx, dy, dz;
925  bpfit::Path path(*geom, zfit.size(), 0.0);
926  for (j=0; j<zfit.size(); ++j) {
927  sys.TrkToDet(xfit[j], yfit[j], zfit[j], 0, 0, 0,
928  &xyz[0], &xyz[1], &xyz[2], &dx, &dy, &dz);
929  this->ClampToDetector(xyz);
930  path.SetPoint(j, xyz);
931  xfit_det.push_back(xyz[0]);
932  yfit_det.push_back(xyz[1]);
933  zfit_det.push_back(xyz[2]);
934  }
935  path.Scrub();
936  if (fExtendUp!=0) path.AdjustStartPoint(fExtendUp);
937  if (fExtendDn!=0) path.AdjustEndPoint(fExtendDn);
938 
939  // Debugging fit parameters in track frame.
940  if (fDebugFitList0) this->DebugFitList(1, xfit, yfit, zfit);
941  // Debugging fit parameters in detector frame.
942  // if (fDebugFitList0) this->DebugFitList(0, xfit_det, yfit_det, zfit_det);
943 
944  //
945  // Log the track parameters
946  //
947  if (fFillTrackNt) this->FillTrackNt(hits1, lutzx, lutzy, path, fPartList[t]);
948 
949  //
950  // Build the track before returning
951  //
952  track.SetDir(path.Dir(0.0));
953  track.SetStart(path.fXYZ[0]);
954  for (i=1; i<path.fXYZ.size(); ++i) {
955  track.AppendTrajectoryPoint(path.fXYZ[i]);
956  }
957 
958  // Calculate the degrees of freedom for each track.
959  unsigned int nhitx = 0, nhity = 0;
960  for (unsigned int i=0; i<hits1.size(); ++i) {
961  if (hits1[i].fView==geo::kX && hits1[i].fW > 1.0E-9) ++nhitx;
962  if (hits1[i].fView==geo::kY && hits1[i].fW > 1.0E-9) ++nhity;
963  }
964 
965  // Check that the number of hits used in each view isn't too small.
966  // The number of degrees of freedom for each track fit is N Hits
967  // that went into making the track minus 4, so if used 0 hits in
968  // either view or <= 4 total, we will throw this exception.
969  if (nhitx < 1 || nhity < 1 || (nhitx+nhity) <= 4) {
970  throw BPException(__FILE__, __LINE__, kBAD_NDOF, nhitx+nhity);
971  }
972 
973  double momentum = 0.0;
974  double mass = 0.0;
975 
976  if(fPartList[t] == 13) {
977  path.MuonParams(0, 0, &momentum, 0, 0);
978  mass = 105.658E-3;
979  }
980  if(fPartList[t] == 211) {
981  path.PionParams(0, 0, &momentum, 0, 0);
982  mass = 139.570E-3;
983  }
984  if(fPartList[t] == 2212) {
985  path.ProtonParams(0, 0, &momentum, 0, 0);
986  mass = 938.272E-3;
987  }
988 
989  // set the kinematic information and chi2 and ndof...
990  fitsum.SetKin(path.Dir(0.0), momentum, mass);
991  fitsum.SetPDG(fPartList[t]);
992  fitsum.SetChi2(lutzx.Chi2XI(),lutzy.Chi2XI(),
993  lutzx.Chi2Beta(),lutzy.Chi2Beta(),
994  nhitx, nhity, scat_surf.N(), scat_surf.N());
995 
996  tracks.push_back(track);
997  fitsums.push_back(fitsum);
998 
999  } // end of try block
1000  catch (BPException& Ex) {
1001  if(fVerbose) Ex.Print();
1002  if(fFillExceptNt) this->FillExceptNt(Ex);
1003  } // end of catch block
1004  } // end loop over fPartList
1005 
1006  return;
1007 
1008  }
1009 
1010  //......................................................................
1011 
1013  {
1014  //
1015  // Allocations for the track products and their associations
1016  //
1017  std::unique_ptr<std::vector<rb::Track> >
1018  trackcol(new std::vector<rb::Track>);
1019  std::unique_ptr<std::vector<rb::FitSum> >
1020  fitcol(new std::vector<rb::FitSum>);
1021  std::unique_ptr< art::Assns<rb::Track, rb::Prong> >
1022  trkPrngAssns(new art::Assns<rb::Track, rb::Prong>);
1023  std::unique_ptr< art::Assns<rb::Track, rb::Cluster> >
1024  trkSlcAssns(new art::Assns<rb::Track, rb::Cluster>);
1025  std::unique_ptr< art::Assns<rb::FitSum, rb::Track> >
1026  fitTrkAssns(new art::Assns<rb::FitSum, rb::Track>);
1027 
1028  fRun = evt.run();
1029  fSubrun = evt.subRun();
1030  fEvent = evt.id().event();
1031  fSlice = 0;
1032  //
1033  // Extract the slices, vertices, and tracks we need as inputs
1034  //
1037  this->FindSlices(evt, sliceh, slicep);
1038  art::FindManyP<rb::Vertex> vertexh(sliceh, evt, fVertexLabel);
1040  prongh(sliceh, evt, art::InputTag(fClusterLabel,fClusterInstance));
1041  if (!vertexh.isValid()) {
1042  return;
1043  }
1044 
1045  unsigned int islice;
1046  for (islice=0; islice<slicep.size(); ++islice) {
1047  fSlice = islice;
1048 
1049  //
1050  // Skip the noise slice - no physics there
1051  //
1052  if (slicep[islice]->IsNoise()) continue;
1053 
1054  //
1055  // To avoid problems with jobs running out of memory and/or taking too
1056  // long to process a single event, skip the slice if it has too many
1057  // hits since the reconstruction will be junk anyway.
1058  //
1059  if (slicep[islice]->NCell() > fMaxSliceHits) continue;
1060 
1061  //
1062  // Find vertex associated with this slice
1063  //
1064  std::vector< art::Ptr<rb::Vertex> > vertexp;
1065  vertexp = vertexh.at(islice);
1066  if (vertexp.size()!=1) continue;
1067 
1068  //
1069  // Find the tracks associated with this slice
1070  //
1071  std::vector< art::Ptr<rb::Prong> > prongp;
1072  this->FindProngs(evt, islice, prongh, vertexp, prongp);
1073  if (prongp.size()==0) continue;
1074 
1075  //
1076  // At this point we should have a list of tracks to fit. Go
1077  // ahead and process them
1078  //
1079  unsigned int nprong_good = 0; // How many tracks pass selection cuts
1080  for (fProng=0; fProng<prongp.size(); ++fProng) {
1081  //
1082  // Figure out if the track has enough information to be
1083  // fit. If not, skip to the next track
1084  //
1085  if (!this->GoodProng(prongp[fProng])) continue;
1086  ++nprong_good;
1087 
1088  //
1089  // Fit the tracks under the requested assumptions.
1090  //
1091  std::vector<rb::Track> tracks;
1092  std::vector<rb::FitSum> fitsums;
1093 
1094  this->FitTracks(vertexp[0], prongp[fProng], tracks, fitsums);
1095 
1096  // Add all returned tracks to the full list and make the appropriate Assns.
1097  assert(tracks.size() == fitsums.size());
1098  for(unsigned int t = 0; t < tracks.size(); ++t) {
1099  trackcol->push_back(tracks[t]);
1100  fitcol->push_back(fitsums[t]);
1101  util::CreateAssn(*this, evt, *trackcol, prongp[fProng], *trkPrngAssns);
1102  util::CreateAssn(*this, evt, *trackcol, slicep[islice] , *trkSlcAssns);
1103  util::CreateAssn(*this, evt, *fitTrkAssns, *fitcol, *trackcol);
1104  }
1105 
1106  } // Loop on prongs
1107  } // Loop on slices
1108 
1109  evt.put(std::move(trackcol));
1110  evt.put(std::move(fitcol));
1111  evt.put(std::move(trkPrngAssns));
1112  evt.put(std::move(trkSlcAssns));
1113  evt.put(std::move(fitTrkAssns));
1114 
1115  } // produce() method
1116 };
1117 
1118 //......................................................................
1119 
1121 
1122 ////////////////////////////////////////////////////////////////////////
const XML_Char int len
Definition: expat.h:262
float TNS() const
Definition: CellHit.h:46
std::string fClusterInstance
Which instance of clusters to find.
SubRunNumber_t subRun() const
Definition: Event.h:72
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 fSigAdj
Multiplier to adjust sigmaXYZ for hits (see comment in the FillHitList function.) Set to 1...
double HalfD() const
Definition: CellGeo.cxx:205
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
void MakeBasis(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &z, const std::vector< double > &dx, const std::vector< double > &dy, const std::vector< double > &dz, double vx, double vy, double vz)
Definition: TrackBasis.cxx:61
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
double fdStol
Scattering surface ds tolerance.
double Length() const
Return the total path length.
Definition: Path.cxx:104
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
TVector3 Dir(double s=0.0)
Return the direction of the track after it has travelled a distance s.
Definition: Path.cxx:120
void SetMeasurement(unsigned int i, double z, double x, double sigx)
Add measurements.
Definition: Lutz.cxx:38
unsigned int fMaxSliceHits
Don&#39;t bother processing a slice if it has more than this many hits. This is to prevent computing prob...
void FillTrackNt(const HitList3D &h, const Lutz &lx, const Lutz &ly, const Path &path, int pdg_id)
unsigned int fNhitY
Require prongs to have this many y hits.
double fdTtol
Scattering surface dt tolerance.
bool Is2D() const
Definition: Cluster.h:96
double GetX() const
Definition: Vertex.h:23
unsigned short Plane() const
Definition: CellHit.h:39
bool fFillTrackNt
Fill ntuple of tracking results?
Construct scattering surfaces for Lutz.
unsigned int fNhitX
Require prongs to have this many x hits.
double fExtendDn
Fractionally, how far to extend track into downstream dead material.
bool GoodProng(art::Ptr< rb::Prong > &t)
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
Vertical planes which measure X.
Definition: PlaneGeo.h:28
void AdjustEndPoint(double f)
Adjust the track end point to account for dead material beyond the end point.
Definition: Path.cxx:208
double HalfW() const
Definition: CellGeo.cxx:191
double DetLength() const
OStream cerr
Definition: OStream.cxx:7
"Break-point" track fitter
void ProtonParams(double s, double *T=0, double *p=0, double *beta=0, double *gamma=0) const
Definition: Path.cxx:462
TString ip
Definition: loadincs.C:5
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
double GetY() const
Definition: Vertex.h:24
const PlaneGeo * Plane(unsigned int i) const
std::string fClusterLabel
Where to find clusters in event.
DEFINE_ART_MODULE(TestTMapFile)
bool fFitAsPion
Fit tracks under the assumption of a pion?
const art::ProductToken< std::vector< rb::Cluster > > fSliceToken
T acos(T number)
Definition: d0nt_math.hpp:54
#define M_PI
Definition: SbMath.h:34
void DebugFitList(unsigned int which, const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &z)
double X(double z) const
After fit, best-fit x location of track at location z.
Definition: Lutz.cxx:207
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
double GetZ() const
Definition: Vertex.h:25
void FindProngs(art::Event &evt, unsigned int islice, art::FindManyP< rb::Prong > &prongh, std::vector< art::Ptr< rb::Vertex > > &vertexp, std::vector< art::Ptr< rb::Prong > > &prongp)
bool fDebugFitList0
Debug the fitted track locations (det. frame)
std::vector< TVector3 > fXYZ
x/y/z points along the path
Definition: Path.h:155
unsigned int N() const
The number of scattering planes constructed.
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
unsigned short Cell() const
Definition: CellHit.h:40
const XML_Char * s
Definition: expat.h:262
void MakeSurfaces(const geo::GeometryBase &geo, const double *p0, const double *p1, const std::vector< double > &s)
Build scattering surfaces.
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
void hits()
Definition: readHits.C:15
double dy[NP][NC]
void FindSlices(const art::Event &evt, art::Handle< std::vector< rb::Cluster > > &h, art::PtrVector< rb::Cluster > &s)
double dx[NP][NC]
void CompleteOrtho(Interp3D interp=0)
Definition: HitList3D.cxx:145
Break-point track fitter.
Definition: Lutz.h:20
void AdjustStartPoint(double f)
Definition: Path.cxx:304
void PionParams(double s, double *T=0, double *p=0, double *beta=0, double *gamma=0) const
Definition: Path.cxx:449
double dz[NP][NC]
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
void FillHitList(art::Ptr< rb::Prong > &prongp, bpfit::HitList3D &hits)
Track parameters (s, X0, KE, beta, ...) along a particular path in the detector.
Definition: Path.h:18
Float_t E
Definition: plot.C:20
void ClampToDetector(double *p)
bool fFitAsMuon
Fit tracks under the assumption of a muon?
virtual void SetDir(TVector3 dir)
Definition: Prong.cxx:104
void SetPDG(int pdg)
Definition: FitSum.h:38
void AppendTrajectoryPoint(TVector3 pt)
Definition: Track.cxx:112
T get(std::string const &key) const
Definition: ParameterSet.h:231
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
double Chi2XI() const
Definition: Lutz.cxx:277
Contain 3D hit information by extrapolating ortho view coordinate from surrounding points in the orth...
const double j
Definition: BetheBloch.cxx:29
std::vector< int > fPartList
List of particle types (by pdg code) to be used for track fits.
double fExtendUp
Fractionally, how far to extend track into upstream dead material.
const unsigned int kNPLANES
too few planes in one (or both) views
Definition: BPException.h:28
double fChi2Cutoff
Drop hits with this chi2 or greater.
void AverageXY(bpfit::HitList3D &h)
double DetHalfHeight() const
void DetToTrk(double x, double y, double z, double dx, double dy, double dz, double *x1, double *y1, double *z1, double *dx1, double *dy1, double *dz1)
Definition: TrackBasis.cxx:194
Vertex location in position and time.
Perform a "2 point" Hough transform on a collection of hits.
z
Definition: test.py:28
bool fDebugScatSurf
Debug the scattering surfaces.
void SetChi2(double chi2_hit1, double chi2_hit2, double chi2_angle1, double chi2_angle2, int Ndof_hit1, int Ndof_hit2, int Ndof_angle1, int Ndof_angle2)
Definition: FitSum.h:43
size_type size() const
Definition: PtrVector.h:308
Construct a set of scattering surfaces for Lutz.
void MuonParams(double s, double *T=0, double *p=0, double *beta=0, double *gamma=0) const
Return the estimated muon parameters at path length s.
Definition: Path.cxx:436
void FitTracks(art::Ptr< rb::Vertex > &vertexp, art::Ptr< rb::Prong > &prongp, std::vector< rb::Track > &tracks, std::vector< rb::FitSum > &fitsums)
void FillExceptNt(BPException Ex)
const unsigned int kBAD_NDOF
too few hits used to make the track
Definition: BPException.h:27
void MakeBasis(art::Ptr< rb::Vertex > &vtx, const HitList3D &hits, bpfit::TrackBasis &sys)
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
const std::string path
Definition: plot_BEN.C:43
EventNumber_t event() const
Definition: EventID.h:116
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
void DebugScatSurf(bpfit::ScatteringSurfaces &s, int pdgid)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
bool fLogResiduals
Make ntuple of residuals?
T * make(ARGS...args) const
unsigned int fNhit
Require this many hits total.
double DetHalfWidth() const
T const * get() const
Definition: Ptr.h:321
void produce(art::Event &evt)
void reconfigure(fhicl::ParameterSet const &p)
unsigned int fFlag
Exception code thrown.
Definition: BPException.h:41
Definition: event.h:1
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
Construct a track-based (x,y,z) coordinate system.
void geom(int which=0)
Definition: geom.C:163
void TrkToDet(double x, double y, double z, double dx, double dy, double dz, double *x1, double *y1, double *z1, double *dx1, double *dy1, double *dz1)
Definition: TrackBasis.cxx:218
double fValue
(optional) Value of bad variable that caused the exception.
Definition: BPException.h:42
bool fVerbose
Run in verbose mode (print messages about all exceptions?)
assert(nhit_max >=nhit_nbins)
BreakPoint(fhicl::ParameterSet const &p)
bool fFitAsProton
Fit tracks under the assumption of a proton?
void SetKin(TVector3 dir, double momentum, double mass)
Definition: FitSum.cxx:48
void SetPoint(unsigned int i, const double *xyz)
Add a space point to the path.
Definition: Path.cxx:60
std::string fVertexLabel
Where to find vertices in event.
unsigned int Scrub(double eps=0.01)
Remove points from the trajectory that have negligible distances between them.
Definition: Path.cxx:90
bool getByToken(ProductToken< PROD > const &token, Handle< PROD > &result) const
Definition: DataViewImpl.h:387
bool fDebugHitList1
Debug the track-frame hit list?
double Chi2Beta() const
Definition: Lutz.cxx:287
void LogResiduals(bpfit::HitList3D &hit, bpfit::ScatteringSurfaces &scat, bpfit::Lutz &xfit, bpfit::Lutz &yfit, double *alphax, double *alphay, float pdgid)
RunNumber_t run() const
Definition: Event.h:77
Encapsulate the cell geometry.
Definition: CellGeo.h:25
TNtuple * nt
Definition: drawXsec.C:2
Float_t track
Definition: plot.C:35
Essential 3D information for tracking.
Definition: HitList3D.h:58
Float_t w
Definition: plot.C:20
ProductToken< T > consumes(InputTag const &)
bool fDebugFitList1
Debug the fitted track locations (track frame)
bool fFillExceptNt
Fill ntuple of exception cases?
EventID id() const
Definition: Event.h:56
Encapsulate the geometry of one entire detector (near, far, ndos)
void FindEndPoints(HitList3D &h, bpfit::TrackBasis &sys, double *p0, double *p1)
bool fDebugHitList0
Debug the detector-frame hit list?
double SigAlpha(unsigned int i) const
The predicted RMS scattering angle at plane i.
void SetHit(unsigned int i, geo::View_t, double x, double y, double z, double t, double sig_x, double sig_y, double sig_z, double sig_t, double w=1.0)
Definition: HitList3D.cxx:29
double S(unsigned int i) const
The location of the ith scattering plane.
T atan2(T number)
Definition: d0nt_math.hpp:72
A container for kinematic information.
Definition: FitSum.h:22
double fdXtol
Scattering surface dx tolerance.
void DebugHitList(bpfit::HitList3D &h, int which)