TrackFit_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: TrackFit
3 // Module Type: producer
4 // File: TrackFit_module.cc
5 //
6 // Generated at Tue Oct 30 14:07:52 2012 by Hongyue Duyang using artmod
7 // from art v1_01_01.
8 ////////////////////////////////////////////////////////////////////////
9 
17 #include "cetlib_except/exception.h"
18 
25 
26 #include "DAQChannelMap/DAQChannelMap.h"
27 
28 #include <cmath>
29 #include <string>
30 #include <vector>
31 
32 namespace novaddt {
33  class TrackFit;
34 }
35 
37 public:
38  explicit TrackFit(fhicl::ParameterSet const & p);
39  virtual ~TrackFit();
40  void endJob() override;
41  virtual void produce(art::Event & e);
42 
43 private:
44 
45  virtual void SeedWeights(const std::vector<double>& x,
46  const std::vector<double>& zx,
47  const std::vector<double>& y,
48  const std::vector<double>& zy,
49  std::vector<double>& w);
50  virtual void FitView(const std::vector<double>& x,
51  const std::vector<double>& zx,
52  std::vector<double>& wx,
53  double *start,
54  double *End,
55  const std::vector<double>& tdc,
56  double &eTDC,
57  double &lTDC
58  );
59  virtual double DsqrToLine(double x0, double y0,
60  double x1, double y1,
61  double x2, double y2);
62  virtual double LinFitMinDperp(const std::vector<double>& x,
63  const std::vector<double>& y,
64  const std::vector<double>& w,
65  double* x1, double* y1,
66  double* x2, double* y2);
67  virtual double LinFit(const std::vector<double>& x,
68  const std::vector<double>& y,
69  const std::vector<double>& w,
70  double* x1, double* y1,
71  double* x2, double* y2);
72  virtual double UtilLinFit(const std::vector<double>& x,
73  const std::vector<double>& y,
74  const std::vector<double>& w,
75  double& m, double& c);
76 
77  // Declare member data here.
80  double fDHitGood;
82  std::vector<double> fR; ///< Range of weight function
83  std::vector<double> fT; ///< Transition width of weight function
85  int n_events = 0;
86  int n_slices_input = 0;
87  int n_tracks_output = 0;
88 };
89 
90 //------------------------------------------------------------------------------
92  : fSliceLabel (p.get< std::string >("slice_label") )
93  , fSliceInstance (p.get< std::string >("slice_instance"))
94  , fDHitGood (p.get<double> ("DHitGood") )
95  , fSortOutputHits(p.get<bool> ("SortOutputHits"))
96  , fR(8)
97  , fT(8)
98  , _assn_to_slice_instance(p.get<std::string>("assn_to_slice_instance"))
99 {
100  fR[0] = 0.0; fT[0] = 0.0;
101  fR[1] = fDHitGood; fT[1] = 1000.0;//fR = 6.0 from 1
102  fR[2] = fDHitGood; fT[2] = 500.0;
103  fR[3] = fDHitGood; fT[3] = 250.0;
104  fR[4] = fDHitGood; fT[4] = 100.0;
105  fR[5] = fDHitGood; fT[5] = 50.0;
106  fR[6] = fDHitGood; fT[6] = 25.0;
107  fR[7] = fDHitGood; fT[7] = 10.0;
108 
109  // Call appropriate Produces<>() functions here.
110  produces< std::vector<novaddt::Track> >();
111  produces< std::vector<novaddt::HitList> >();
112 
113  // Associate Tracks with the Slice they came from
114  produces< art::Assns<novaddt::Track, novaddt::HitList> >();
115  produces< art::Assns<novaddt::Track, novaddt::HitList> >(_assn_to_slice_instance);
116 }
117 
118 //------------------------------------------------------------------------------
120 {
121  // Clean up dynamic memory and other resources here.
122 }
123 
124 //------------------------------------------------------------------------------
126 {
127  LOG_DEBUG("Tracking") << "event: " << e.id().event() << std::endl;
128  n_events++;
129  //get sliced hits
131  e.getByLabel(fSliceLabel, fSliceInstance, slices);//TimeSlice
132  //Now make it an art::PtrVector to use in the associations later
134  for(unsigned int i = 0; i < slices->size(); ++i){
135  art::Ptr<novaddt::HitList>slice(slices,i);
136  slicelist.push_back(slice);
137  }
138 
139  //a container of track object to be write into event
140  std::unique_ptr< std::vector<novaddt::Track> > trackcol (new std::vector<novaddt::Track> );
141  std::unique_ptr< std::vector<novaddt::HitList> > hitListcol(new std::vector<novaddt::HitList>);
142  std::unique_ptr< art::Assns<novaddt::Track, novaddt::HitList> > assn(new art::Assns<novaddt::Track, novaddt::HitList>);
143  std::unique_ptr< art::Assns<novaddt::Track, novaddt::HitList> > assn_b(new art::Assns<novaddt::Track, novaddt::HitList>);
144 
145  //loop over grouped hit list
146  for(size_t i = 0; i < slices->size(); ++i){
147  n_slices_input++;
148  art::Ptr<novaddt::HitList> slicedhits(slices,i);
149 
150  std::vector<double> xzViewX;
151  std::vector<double> xzViewZ;
152  std::vector<double> yzViewY;
153  std::vector<double> yzViewZ;
154  std::vector<double> xzViewW;
155  std::vector<double> yzViewW;
156 
157  std::vector<double> xzTDC;
158  std::vector<double> yzTDC;
159 
160  //earliest and latest TDC
161  double eTDC = 0;
162  double lTDC = slicedhits->at(0).TDC().val;
163 
164  LOG_DEBUG("Tracking") << "new slice" << std::endl;
165  unsigned int counter = 0;
166  //loop over hits in each slice
167  for(auto const& j : *slicedhits){
168  std::string view = "y";
169  if( j.View().val == daqchannelmap::X_VIEW ){
170  view = "x";
171  }
172  LOG_DEBUG("Tracking") << "\thit[" << counter
173  << "]: TDC: " << j.TDC().val
174  << ", ADC: " << j.ADC().val
175  << ", plane: " << j.Plane().val
176  << ", " << view
177  << "-cell: " << j.Cell().val
178  << std::endl;
179  counter++;
180  //DAQHit hit(slicedhits[j]);
181 
182  //initiallize eTDC and lTDC
183  eTDC = j.TDC().val;
184 
185  //get hit position in xz view and yz view
186  if( j.View().val == daqchannelmap::X_VIEW ){//xz view
187  xzViewX.push_back(j.Cell().val);
188  xzViewZ.push_back(j.Plane().val);
189  xzViewW.push_back(0);
190  xzTDC.push_back(j.TDC().val);
191  }
192  else if( j.View().val == daqchannelmap::Y_VIEW ){//yz view
193  yzViewY.push_back(j.Cell().val);
194  yzViewZ.push_back(j.Plane().val);
195  yzViewW.push_back(0);
196  yzTDC.push_back(j.TDC().val);
197  }
198  //LOG_DEBUG("Tracking") << "plane: " << j.Plane().val
199  //<< " cell: " << j.Cell().val
200  //<< " tdc: " << j.TDC().val
201  //<< std::endl;
202  }//end loop over hits in each slice
203 
204  //get seed weight
205  this->SeedWeights(xzViewX, xzViewZ, yzViewY, yzViewZ, xzViewW);
206  this->SeedWeights(yzViewY, yzViewZ, xzViewX, xzViewZ, yzViewW);
207 
208  //do a 2D fit, 0 is z (plane), 1 is x/y (cell)
209  double xzStart[2] = {0.};//[0]=z, [1]=x
210  double xzEnd[2] = {0.};
211  double yzStart[2] = {0.};//[0]=z, [1]=y
212  double yzEnd[2] = {0.};
213 
214  this->FitView(xzViewX, xzViewZ, xzViewW, xzStart, xzEnd, xzTDC, eTDC, lTDC);
215  this->FitView(yzViewY, yzViewZ, yzViewW, yzStart, yzEnd, yzTDC, eTDC, lTDC);
216 
217  // protect against bad fits
218  if(xzStart[0] == 0. &&
219  xzStart[1] == 0. &&
220  xzEnd[0] == 0. &&
221  xzEnd[1] == 0. &&
222  yzStart[0] == 0. &&
223  yzStart[1] == 0. &&
224  yzEnd[0] == 0. &&
225  yzEnd[1] == 0.) continue;
226 
227  double xzSlope = (xzStart[1] - xzEnd[1])/(xzStart[0] - xzEnd[0]);
228  double yzSlope = (yzStart[1] - yzEnd[1])/(yzStart[0] - yzEnd[0]);
229  double xzInt = xzStart[1] - xzSlope*xzStart[0];
230  double yzInt = yzStart[1] - yzSlope*yzStart[0];
231 
232  // go through and make an association between hits on track
233  // and the track
234  novaddt::HitList hl_xz;
235  novaddt::HitList hl_yz;
236 
237  for(auto const& h : *slicedhits){
238  if( h.View().val == daqchannelmap::X_VIEW ){//xz view
239  if(std::abs(h.Plane().val*xzSlope + xzInt - h.Cell().val) < fDHitGood){
240  hl_xz.push_back(h);
241  }
242  }
243  else if( h.View().val == daqchannelmap::Y_VIEW ){//yz view
244  if(std::abs(h.Plane().val*yzSlope + yzInt - h.Cell().val) < fDHitGood){
245  hl_yz.push_back(h);
246  }
247  }
248  }//end loop over hits in each slice
249 
250  // protect against bad fits
251  if((hl_xz.size() < 2) &&
252  (hl_yz.size() < 2)) continue;
253 
254  LOG_DEBUG("Tracking") << "--- Track found: x: " << xzStart[1]
255  << " - " << xzEnd[1]
256  << ", z: " << xzStart[0]
257  << " - " << xzEnd[0]
258  << ", n hits: " << hl_xz.size()
259  << std::endl;
260  LOG_DEBUG("Tracking") << "--- Track found: y: " << yzStart[1]
261  << " - " << yzEnd[1]
262  << ", z: " << yzStart[0]
263  << " - " << yzEnd[0]
264  << ", n hits: " << hl_yz.size()
265  << std::endl;
266 
267  if (fSortOutputHits){
268  sort(hl_xz.begin(), hl_xz.end(), CompareDAQHit<TDC>());
269  sort(hl_yz.begin(), hl_yz.end(), CompareDAQHit<TDC>());
270  }
271  //Implement the starting and end time of the track according to your own algorithm:
272  float startt =0;
273  float endt = 0;
274  n_tracks_output+=2;
276  xzStart[1], xzStart[0],startt,
277  xzEnd[1], xzEnd[0],endt
278  );
280  yzStart[1], yzStart[0],startt,
281  yzEnd[1], yzEnd[0],endt
282  );
283  //All tracks registered here are supposed to be c speed, but I made them temperarily all yz view as well. Will be changed.
284  trackcol->push_back(x_track);
285  hitListcol->push_back(hl_xz);
286  util::CreateAssn(*this, e, *trackcol, *hitListcol, *assn, hitListcol->size()-1, hitListcol->size());
287  util::CreateAssn(*this, e, *trackcol, slicelist[i], *assn_b);
288  trackcol->push_back(y_track); // push back a second time for the yz view
289  hitListcol->push_back(hl_yz);
290  util::CreateAssn(*this, e, *trackcol, *hitListcol, *assn, hitListcol->size()-1, hitListcol->size());
291  util::CreateAssn(*this, e, *trackcol, slicelist[i], *assn_b);
292 
293  }//end loop over groups
294 
295  //save the fitted track into the event
296  e.put(std::move(trackcol));
297  e.put(std::move(hitListcol));
298  e.put(std::move(assn));
299  e.put(std::move(assn_b), _assn_to_slice_instance);
300 
301 }//end produce
302 //---------------------------------------------------------------------
304 {
305  std::cout << "=== novaddt::TrackFit endJob" << std::endl;
306  std::cout << "\tNumber of events: " << n_events << std::endl;
307  std::cout << "\tNumber of slices input: " << n_slices_input << std::endl;
308  std::cout << "\tNumber of tracks output: " << n_tracks_output << std::endl;
309 }
310 
311 //---------------------------------------------------------------------
312 //................................................................
313 //Get seed weights. copied from offline RLFit.cxx
314 void novaddt::TrackFit::SeedWeights(const std::vector<double>& x,
315  const std::vector<double>& zx,
316  const std::vector<double>& y,
317  const std::vector<double>& zy,
318  std::vector<double>& w)
319 {
320  for (unsigned int i=0; i<zx.size(); ++i) {
321  double d;
322  w[i] = 1.0E-9;
323  //
324  // Find closest hit in same view
325  //
326  for (unsigned int j=0; j<zx.size(); ++j) {
327  if (i==j) continue;
328  //d = util::pythag(zx[i]-zx[j], x[i]-x[j]);
329  d = sqrt( (zx[i]-zx[j])*(zx[i]-zx[j]) + (x[i]-x[j])*(x[i]-x[j]) );
330  if (d<3.0) w[i] += 1.0;
331  }
332  //
333  // Find closest hit in orthogonal view
334  //
335  for (unsigned int j=0; j<zy.size(); ++j) {
336  d = std::abs(zx[i]-zy[j]);
337  if (d<3.0) w[i] += 0.1;
338  }
339  }
340 }//end SeedWeights
341 
342 
343 //...................................................................
344 //Modify RLFit to fit 2D track
345 void novaddt::TrackFit::FitView(const std::vector<double>& x,
346  const std::vector<double>& zx,
347  std::vector<double>& wx,
348  double *start,
349  double *end,
350  const std::vector<double>& tdc,
351  double &eTDC,
352  double &lTDC)
353 {
354  //
355  // Make a line fit to the hits. To improve preformance in a noisy
356  // environment, use a deterministic annealing algoritm
357  //
358  double x1 = 0;
359  double x2 = 0;
360  double zx1 = 0;
361  double zx2 = 0;
362 
363  double slope = 0;
364  double intercept = 0;
365 
366  //loop over fR
367  for (unsigned int i=0; i<fR.size(); ++i) {
368  static const double eps = 1.0E-6;
369 
370  // Calculate weights for this iteration. First iteration uses seed weights
371  if (i>0) {
372  for (unsigned int j=0; j<wx.size(); ++j) {
373  double d = this->DsqrToLine(zx[j], x[j], zx1, x1, zx2, x2);
374 
375  wx[j] = (1.0+exp(-fR[i]/fT[i]))/(1.0+exp((d-fR[i])/fT[i]));
376 
377  // Penalty for hits outside the detector in the other view
378  //double y = fMzy*fZx[j] + fBzy;
379  //if (y<fBoxYlo) fWx[j] *= exp(-(fBoxYlo-y)/fT[i]);
380  //if (y>fBoxYhi) fWx[j] *= exp(-(y-fBoxYhi)/fT[i]);
381 
382  if (wx[j]<eps) wx[j] = eps;
383  }
384  }
385 
386  // Perform the fit
387  if(zx.size() != x.size() ||
388  wx.size() != x.size() ||
389  x.size() < 2){
390  return;
391  }
392 
393  this->LinFitMinDperp(zx, x, wx, &zx1, &x1, &zx2, &x2);
394 
395  // Compute slopes and intercepts
396  double dzx = zx2-zx1;
397  if (dzx<=0.0) dzx = 1.0E-9;
398 
399  slope = (x2-x1)/dzx;
400  intercept = (x1*zx2-x2*zx1)/dzx;
401  }//end fR loop
402  LOG_DEBUG("Tracking") << "x1: " << x1
403  << ", x2: " << x2
404  << ", zx1: " << zx1
405  << ", zx2: " << zx2
406  << std::endl;
407  // Adjust z ranges to completely contain hits on the track.
408  //"On track" is defined as having significant weight
409 
410  start[0] = zx[0];
411  end[0] = zx[0]+0.01;
412 
413  for (size_t i = 0; i < wx.size(); ++i) {
414  if (wx[i]>0.001) {
415  if (zx[i] < start[0]) start[0] = zx[i];
416  if (zx[i] > end[0]) end[0] = zx[i];
417  if (tdc[i] < eTDC) eTDC = tdc[i];
418  if (tdc[i] > lTDC) lTDC = tdc[i];
419  }
420  }
421  // Use line fits to compute x and y locations at the start and end z positions
422  LOG_DEBUG("Tracking") << "start: " << start[0] << ", end: " << end[0] << ", slope: " << slope << ", intercept: " << intercept << std::endl;
423  start[1] = slope*start[0] + intercept;
424  end[1] = slope*end[0] + intercept;
425  // protect against lines exiting the detector,
426  // at the moment the upper cell number is hard coded for the FD
427  if (start[1]<=0.) start[1] = 1.0E-9;
428  if (end[1]<=0.) end[1] = 1.0E-9;
429  if (start[1]>383.) start[1] = 383.;
430  if (end[1]>383.) end[1] = 383.;
431 
432 }//end FitView
433 
434 //...................................................................
435 //
436 double novaddt::TrackFit::DsqrToLine(double x0, double y0,
437  double x1, double y1,
438  double x2, double y2)
439 {
440  // If we've been handed a point instead of a line, return distance
441  // to the point
442  if (x1==x2 && y1==y2) return (x0-x1)*(x0-x1) +(y0-y1)*(y0-y1);
443 
444  return
445  ((x2-x1)*(y1-y0)-(x1-x0)*(y2-y1))*((x2-x1)*(y1-y0)-(x1-x0)*(y2-y1)) /
446  ( (x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) );
447 }
448 
449 //....................................................................
450 //
451 double novaddt::TrackFit::LinFitMinDperp(const std::vector<double>& x,
452  const std::vector<double>& y,
453  const std::vector<double>& w,
454  double* x1, double* y1,
455  double* x2, double* y2)
456 {
457  // check array sizes before passing them to this method to avoid
458  // cases where we would be tempted to throw exceptions
459  // x, y, and w must all be the same size and have size > 2
460 
461  unsigned register int i;
462 
463  // Find the extremes of the inputs
464  double xlo = DBL_MAX;
465  double xhi = -DBL_MAX;
466  double ylo = DBL_MAX;
467  double yhi = -DBL_MAX;
468  for (i=0; i<w.size(); ++i) {
469  if(w[i] == 0) continue;
470  if (x[i]<xlo) xlo = x[i];
471  if (y[i]<ylo) ylo = y[i];
472  if (x[i]>xhi) xhi = x[i];
473  if (y[i]>yhi) yhi = y[i];
474  }
475 
476  // This algorithm really fails if the y values are identical
477  // return the low and hi values if it's too flat
478  if(std::abs(yhi-ylo) < 0.1){
479  *x1 = xlo;
480  *y1 = ylo;
481  *x2 = xhi;
482  *y2 = yhi;
483  return 0.; //its a flat line, chi^2 is 0 by definition
484  }
485 
486  // For stability, fit in a way that avoids infinite slopes,
487  // exchanging the x and y variables as needed.
488  if (std::abs(xlo - xhi) < 1.e-3) {
489  double chi2;
490  double xx1, yy1, xx2, yy2;
491  // Notice: (x,y) -> (y,x)
492  chi2 = LinFitMinDperp(y,x,w,&xx1, &yy1,&xx2, &yy2);
493 
494  *x1 = xlo; // Infinite slope means that your initial x positions should be the input
495  *y1 = xx1;
496  *x2 = xhi;
497  *y2 = xx2;
498  return chi2;
499  }
500 
501  // Accumulate the sums needed to compute the fit line
502  double Sw = 0.0;
503  double Swx = 0.0;
504  double Swy = 0.0;
505  double Swxy= 0.0;
506  double Swy2= 0.0;
507  double Swx2= 0.0;
508  for (i=0; i<w.size(); ++i) {
509  Sw += w[i];
510  Swx += w[i]*x[i];
511  Swy += w[i]*y[i];
512  Swxy += w[i]*x[i]*y[i];
513  Swy2 += w[i]*y[i]*y[i];
514  Swx2 += w[i]*x[i]*x[i];
515  }
516  if(Sw < 0.0)
517  throw cet::exception("TrackFit") << "Sw < 0";
518 
519  // Precalculate reused squares
520  //using util::sqr;
521  const double SwSq = Sw*Sw;//sqr(Sw);
522  const double SwxSq = Swx*Swx;//sqr(Swx);
523  const double SwySq = Swy*Swy;//sqr(Swy);
524 
525  // C and D are two handy temporary variables
526  double C =
527  +SwxSq*SwxSq // Swx^4
528  -2.0*SwxSq*Swx2*Sw
529  +2.0*SwxSq*SwySq
530  +2.0*SwxSq*Swy2*Sw
531  +Swx2*Swx2*SwSq
532  +2.0*Swx2*Sw*SwySq
533  -2.0*Swx2*SwSq*Swy2
534  +pow(Swy, 4) // Swy^4
535  -2.0*SwySq*Swy2*Sw
536  +Swy2*Swy2*SwSq
537  -8.0*Swx*Swy*Swxy*Sw
538  +4.0*Swxy*Swxy*SwSq;
539  if (C>0.0) C = sqrt(C);
540  else C = 0.0;
541  double D = -Swx*Swy+Swxy*Sw;
542 
543  // The first of the two solutions. This one is the best fit through
544  // the points
545  double M1;
546  double B1;
547  if (D==0.0) {
548  // D could be zero, this might happen in a segmented detector if
549  // the track fit passed through only a single readout plane and
550  // all hits share the same x coordinates. I've tried to avoid
551  // this by checking if it makes more sense to fit with x and y
552  // exchanged above. This is an extra precaution.
553  D = 1.0E-9;
554  }
555  M1 = -(-SwxSq + Swx2*Sw + SwySq - Swy2*Sw - C)/(2.0*D);
556  B1 = -(-Swy+M1*Swx)/Sw;
557 
558  // This second solution is perpendicular to the first and represents
559  // the worst fit
560  // M2 = -(-SwxSq + Swx2*Sw + SwySq - Swy2*Sw + C)/(2.0*D);
561  // B2 = -(-Swy+M2*Swx)/Sw;
562 
563  *x1 = xlo;
564  *y1 = M1*xlo + B1;
565  *x2 = xhi;
566  *y2 = M1*xhi + B1;
567 
568  // Compute the chi^2 function for return
569  double chi2 = 0.0;
570  for (i=0; i<w.size(); ++i) {
571  chi2 += w[i]*this->DsqrToLine(x[i],y[i],*x1,*y1,*x2, *y2);
572  }
573  // test if linfit is better
574 
575  double chi2lin, linx1, linx2, liny1, liny2;
576 
577  chi2lin = LinFit(x,y,w,&linx1,&liny1,&linx2,&liny2);
578 
579  if(chi2lin < chi2) { //something is wrong with the perp fit, use the lin fit
580  *x1 = linx1;
581  *y1 = liny1;
582  *x2 = linx2;
583  *y2 = liny2;
584 
585  chi2 = chi2lin;
586  }
587 
588 
589  return chi2;
590 }//end LinFitMinDperp
591 
592 
593 //...............................................................
594 //
595 double novaddt::TrackFit::LinFit(const std::vector<double>& x,
596  const std::vector<double>& y,
597  const std::vector<double>& w,
598  double* x1, double* y1,
599  double* x2, double* y2)
600 {
601  register unsigned int i;
602 
603  // Before going ahead, make sure we have sensible arrays
604  assert(x.size()==y.size());
605  assert(x.size()==w.size());
606  assert(x.size()>=2);
607 
608  // Find the ranges of the input variables
609  double xlo = x[0];
610  double xhi = x[0];
611  double ylo = y[0];
612  double yhi = y[0];
613 
614  for (i=1; i<w.size(); ++i) {
615  if (x[i]<xlo) xlo = x[i];
616  if (x[i]>xhi) xhi = x[i];
617  if (y[i]<ylo) ylo = y[i];
618  if (y[i]>yhi) yhi = y[i];
619  }
620  // Make sure the fit has some span
621  if (yhi-ylo==0.0 && xhi-xlo==0.0) {
622  std::cerr << __FILE__ << ":" << __LINE__
623  << " All points identical in line fit!" << std::endl;
624  *x1 = xlo;
625  *y1 = ylo;
626  *x2 = xhi;
627  *y2 = yhi;
628  return 0.0;
629  }
630 
631  //
632  // If the y extent is larger than the x extent, flip the variables
633  // and fit
634  //
635  if (yhi-ylo > xhi-xlo) {
636  return this->LinFit(y,x,w,y1,x1,y2,x2);
637  }
638 
639  double m, b;
640  const double chi2 = this->UtilLinFit(x, y, w, m, b);
641 
642  *x1 = xlo;
643  *x2 = xhi;
644  *y1 = m*(*x1)+b;
645  *y2 = m*(*x2)+b;
646 
647  return chi2;
648 }
649 
650 //..................................................................
651 //copied from offline Util::LinFit
652 double novaddt::TrackFit::UtilLinFit(const std::vector<double>& x,
653  const std::vector<double>& y,
654  const std::vector<double>& w,
655  double& m, double& c)
656 {
657  // Before going ahead, make sure we have sensible arrays
658  assert(x.size() == y.size());
659  assert(x.size() == w.size());
660  assert(x.size() >= 2);
661 
662  // Accumulate the sums for the fit
663  double Sw = 0;
664  double Swx = 0;
665  double Swy = 0;
666  double Swxy = 0;
667  double Swy2 = 0;
668  double Swx2 = 0;
669  for(unsigned int i = 0; i < w.size(); ++i) {
670  Sw += w[i];
671  Swx += w[i]*x[i];
672  Swy += w[i]*y[i];
673  Swx2 += w[i]*x[i]*x[i];
674  Swxy += w[i]*x[i]*y[i];
675  Swy2 += w[i]*y[i]*y[i];
676  }
677  const double d = Sw*Swx2 - Swx*Swx;
678  m = (Sw*Swxy - Swx*Swy)/d;
679  c = (Swy*Swx2 - Swx*Swxy)/d;
680 
681  const double chi2 =
682  Swy2 - 2.0*m*Swxy - 2.0*c*Swy + 2.0*m*c*Swx +
683  c*c*Sw + m*m*Swx2;
684 
685  return chi2;
686 }
687 
689 
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
virtual void produce(art::Event &e)
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.
Float_t y1[n_points_granero]
Definition: compare.C:5
virtual 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)
Float_t x1[n_points_granero]
Definition: compare.C:5
std::vector< DAQHit > HitList
Definition: HitList.h:15
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
OStream cerr
Definition: OStream.cxx:7
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
virtual void FitView(const std::vector< double > &x, const std::vector< double > &zx, std::vector< double > &wx, double *start, double *End, const std::vector< double > &tdc, double &eTDC, double &lTDC)
DEFINE_ART_MODULE(TestTMapFile)
std::vector< double > fT
Transition width of weight function.
std::string fSliceLabel
double chi2()
virtual void SeedWeights(const std::vector< double > &x, const std::vector< double > &zx, const std::vector< double > &y, const std::vector< double > &zy, std::vector< double > &w)
float abs(float number)
Definition: d0nt_math.hpp:39
virtual double UtilLinFit(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double &m, double &c)
Identifier for the Y measuring view of the detector (side)
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
const double C
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
Identifier for the X measuring view of the detector (top)
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
Float_t d
Definition: plot.C:236
const double j
Definition: BetheBloch.cxx:29
std::string fSliceInstance
void endJob() override
OStream cout
Definition: OStream.cxx:6
virtual double DsqrToLine(double x0, double y0, double x1, double y1, double x2, double y2)
EventNumber_t event() const
Definition: EventID.h:116
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::string _assn_to_slice_instance
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void End(void)
Definition: gXSecComp.cxx:210
std::vector< double > fR
Range of weight function.
const hit & b
Definition: hits.cxx:21
TrackFit(fhicl::ParameterSet const &p)
assert(nhit_max >=nhit_nbins)
Float_t e
Definition: plot.C:35
Float_t w
Definition: plot.C:20
Definition: fwd.h:28
size_t size() const
EventID id() const
Definition: Event.h:56
virtual double LinFitMinDperp(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double *x1, double *y1, double *x2, double *y2)