Public Types | Public Member Functions | Static Public Member Functions | Protected Member Functions | Private Member Functions | Private Attributes | List of all members
novaddt::TrackFit Class Reference
Inheritance diagram for novaddt::TrackFit:
art::EDProducer art::ProducerBase art::Consumer art::EngineCreator art::ProductRegistryHelper

Public Types

using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
template<typename UserConfig , typename KeysToIgnore = void>
using Table = ProducerBase::Table< UserConfig, KeysToIgnore >
 

Public Member Functions

 TrackFit (fhicl::ParameterSet const &p)
 
virtual ~TrackFit ()
 
void endJob () override
 
virtual void produce (art::Event &e)
 
template<typename PROD , BranchType B = InEvent>
ProductID getProductID (std::string const &instanceName={}) const
 
template<typename PROD , BranchType B>
ProductID getProductID (ModuleDescription const &moduleDescription, std::string const &instanceName) const
 
bool modifiesEvent () const
 
template<typename T , BranchType = InEvent>
ProductToken< Tconsumes (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< Tconsumes (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TconsumesView (InputTag const &it)
 
template<typename T , BranchType = InEvent>
ProductToken< TmayConsume (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< TmayConsume (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TmayConsumeView (InputTag const &it)
 
base_engine_tcreateEngine (seed_t seed)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make, label_t const &engine_label)
 
seed_t get_seed_value (fhicl::ParameterSet const &pset, char const key[]="seed", seed_t const implicit_seed=-1)
 

Static Public Member Functions

static cet::exempt_ptr< Consumernon_module_context ()
 

Protected Member Functions

CurrentProcessingContext const * currentContext () const
 
void validateConsumedProduct (BranchType const bt, ProductInfo const &pi)
 
void prepareForJob (fhicl::ParameterSet const &pset)
 
void showMissingConsumes () const
 

Private Member Functions

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)
 
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)
 
virtual double DsqrToLine (double x0, double y0, double x1, double y1, double x2, double y2)
 
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)
 
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)
 
virtual double UtilLinFit (const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double &m, double &c)
 

Private Attributes

std::string fSliceLabel
 
std::string fSliceInstance
 
double fDHitGood
 
bool fSortOutputHits
 
std::vector< double > fR
 Range of weight function. More...
 
std::vector< double > fT
 Transition width of weight function. More...
 
std::string _assn_to_slice_instance
 
int n_events = 0
 
int n_slices_input = 0
 
int n_tracks_output = 0
 

Detailed Description

Definition at line 36 of file TrackFit_module.cc.

Member Typedef Documentation

using art::EDProducer::ModuleType = EDProducer
inherited

Definition at line 34 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>
using art::EDProducer::Table = ProducerBase::Table<UserConfig, KeysToIgnore>
inherited

Definition at line 43 of file EDProducer.h.

using art::EDProducer::WorkerType = WorkerT<EDProducer>
inherited

Definition at line 35 of file EDProducer.h.

Constructor & Destructor Documentation

novaddt::TrackFit::TrackFit ( fhicl::ParameterSet const &  p)
explicit

Definition at line 91 of file TrackFit_module.cc.

References _assn_to_slice_instance, fDHitGood, fR, and fT.

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 }
const char * p
Definition: xmltok.h:285
std::vector< double > fT
Transition width of weight function.
std::string fSliceLabel
std::string fSliceInstance
std::string _assn_to_slice_instance
std::vector< double > fR
Range of weight function.
enum BeamMode string
novaddt::TrackFit::~TrackFit ( )
virtual

Definition at line 119 of file TrackFit_module.cc.

120 {
121  // Clean up dynamic memory and other resources here.
122 }

Member Function Documentation

template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::consumes ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::consumes ( InputTag const &  it)
inherited

Definition at line 146 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

147 {
148  if (!moduleContext_)
149  return ProductToken<T>::invalid();
150 
151  consumables_[BT].emplace_back(ConsumableType::Product,
152  TypeID{typeid(T)},
153  it.label(),
154  it.instance(),
155  it.process());
156  return ProductToken<T>{it};
157 }
set< int >::iterator it
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::consumesMany ( )
inherited

Definition at line 161 of file Consumer.h.

References T.

162 {
163  if (!moduleContext_)
164  return;
165 
166  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
167 }
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::consumesView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::consumesView ( InputTag const &  it)
inherited

Definition at line 171 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

172 {
173  if (!moduleContext_)
174  return ViewToken<T>::invalid();
175 
176  consumables_[BT].emplace_back(ConsumableType::ViewElement,
177  TypeID{typeid(T)},
178  it.label(),
179  it.instance(),
180  it.process());
181  return ViewToken<T>{it};
182 }
set< int >::iterator it
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make 
)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make,
label_t const &  engine_label 
)
inherited
CurrentProcessingContext const* art::EDProducer::currentContext ( ) const
protectedinherited
double novaddt::TrackFit::DsqrToLine ( double  x0,
double  y0,
double  x1,
double  y1,
double  x2,
double  y2 
)
privatevirtual

Definition at line 436 of file TrackFit_module.cc.

References x1, and y1.

Referenced by FitView(), and LinFitMinDperp().

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 }
Float_t y1[n_points_granero]
Definition: compare.C:5
Float_t x1[n_points_granero]
Definition: compare.C:5
void novaddt::TrackFit::endJob ( )
overridevirtual

Reimplemented from art::EDProducer.

Definition at line 303 of file TrackFit_module.cc.

References om::cout, allTimeWatchdog::endl, n_events, n_slices_input, and n_tracks_output.

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 }
OStream cout
Definition: OStream.cxx:6
void novaddt::TrackFit::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 
)
privatevirtual

Definition at line 345 of file TrackFit_module.cc.

References d, DsqrToLine(), allTimeWatchdog::endl, stan::math::exp(), fR, fT, MECModelEnuComparisons::i, calib::j, LinFitMinDperp(), LOG_DEBUG, x1, and submit_syst::x2.

Referenced by produce().

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
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
Float_t x1[n_points_granero]
Definition: compare.C:5
std::vector< double > fT
Transition width of weight function.
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
virtual double DsqrToLine(double x0, double y0, double x1, double y1, double x2, double y2)
std::vector< double > fR
Range of weight function.
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)
seed_t art::EngineCreator::get_seed_value ( fhicl::ParameterSet const &  pset,
char const  key[] = "seed",
seed_t const  implicit_seed = -1 
)
inherited
template<typename PROD , BranchType B>
ProductID art::EDProducer::getProductID ( std::string const &  instanceName = {}) const
inlineinherited

Definition at line 123 of file EDProducer.h.

References art::EDProducer::moduleDescription_.

Referenced by skim::NueSkimmer::CopyMichelSlice(), and skim::NueSkimmer::CopyMichelTrack().

124  {
125  return ProducerBase::getProductID<PROD, B>(moduleDescription_,
126  instanceName);
127  }
ModuleDescription moduleDescription_
Definition: EDProducer.h:115
template<typename PROD , BranchType B>
ProductID art::ProducerBase::getProductID ( ModuleDescription const &  moduleDescription,
std::string const &  instanceName 
) const
inherited

Definition at line 56 of file ProducerBase.h.

References art::ModuleDescription::moduleLabel().

Referenced by art::ProducerBase::modifiesEvent().

58  {
59  auto const& pd =
60  get_ProductDescription<PROD>(B, md.moduleLabel(), instanceName);
61  return pd.productID();
62  }
double novaddt::TrackFit::LinFit ( const std::vector< double > &  x,
const std::vector< double > &  y,
const std::vector< double > &  w,
double *  x1,
double *  y1,
double *  x2,
double *  y2 
)
privatevirtual

Definition at line 595 of file TrackFit_module.cc.

References ana::assert(), b, om::cerr, chi2(), allTimeWatchdog::endl, MECModelEnuComparisons::i, m, UtilLinFit(), xhi, make_syst_table_plots::xlo, yhi, and ylo.

Referenced by LinFitMinDperp().

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 }
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
OStream cerr
Definition: OStream.cxx:7
double chi2()
virtual double UtilLinFit(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double &m, double &c)
const hit & b
Definition: hits.cxx:21
assert(nhit_max >=nhit_nbins)
Float_t w
Definition: plot.C:20
double novaddt::TrackFit::LinFitMinDperp ( const std::vector< double > &  x,
const std::vector< double > &  y,
const std::vector< double > &  w,
double *  x1,
double *  y1,
double *  x2,
double *  y2 
)
privatevirtual

Definition at line 451 of file TrackFit_module.cc.

References std::abs(), C, chi2(), DsqrToLine(), MECModelEnuComparisons::i, LinFit(), cet::pow(), std::sqrt(), xhi, make_syst_table_plots::xlo, yhi, and ylo.

Referenced by FitView().

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
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
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
double chi2()
float abs(float number)
Definition: d0nt_math.hpp:39
const double C
virtual double DsqrToLine(double x0, double y0, double x1, double y1, double x2, double y2)
Float_t w
Definition: plot.C:20
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)
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::mayConsume ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::mayConsume ( InputTag const &  it)
inherited

Definition at line 189 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

190 {
191  if (!moduleContext_)
192  return ProductToken<T>::invalid();
193 
194  consumables_[BT].emplace_back(ConsumableType::Product,
195  TypeID{typeid(T)},
196  it.label(),
197  it.instance(),
198  it.process());
199  return ProductToken<T>{it};
200 }
set< int >::iterator it
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::mayConsumeMany ( )
inherited

Definition at line 204 of file Consumer.h.

References T.

205 {
206  if (!moduleContext_)
207  return;
208 
209  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
210 }
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::mayConsumeView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::mayConsumeView ( InputTag const &  it)
inherited

Definition at line 214 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

215 {
216  if (!moduleContext_)
217  return ViewToken<T>::invalid();
218 
219  consumables_[BT].emplace_back(ConsumableType::ViewElement,
220  TypeID{typeid(T)},
221  it.label(),
222  it.instance(),
223  it.process());
224  return ViewToken<T>{it};
225 }
set< int >::iterator it
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
bool art::ProducerBase::modifiesEvent ( ) const
inlineinherited

Definition at line 40 of file ProducerBase.h.

References art::ProducerBase::getProductID(), and string.

41  {
42  return true;
43  }
static cet::exempt_ptr<Consumer> art::Consumer::non_module_context ( )
staticinherited
void art::Consumer::prepareForJob ( fhicl::ParameterSet const &  pset)
protectedinherited
void novaddt::TrackFit::produce ( art::Event e)
virtual

Implements art::EDProducer.

Definition at line 125 of file TrackFit_module.cc.

References _assn_to_slice_instance, std::abs(), hadd_many_files::counter, util::CreateAssn(), allTimeWatchdog::endl, art::EventID::event(), fDHitGood, FitView(), fSliceInstance, fSliceLabel, fSortOutputHits, art::DataViewImpl::getByLabel(), make_syst_table_plots::h, MECModelEnuComparisons::i, art::Event::id(), calib::j, LOG_DEBUG, n_events, n_slices_input, n_tracks_output, art::PtrVector< T >::push_back(), art::Event::put(), SeedWeights(), art::DataViewImpl::size(), string, POTSpillRate::view, daqchannelmap::X_VIEW, and daqchannelmap::Y_VIEW.

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
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
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.
std::vector< DAQHit > HitList
Definition: HitList.h:15
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)
std::string fSliceLabel
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
Identifier for the Y measuring view of the detector (side)
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
Identifier for the X measuring view of the detector (top)
const double j
Definition: BetheBloch.cxx:29
std::string fSliceInstance
EventNumber_t event() const
Definition: EventID.h:116
std::string _assn_to_slice_instance
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Definition: fwd.h:28
size_t size() const
EventID id() const
Definition: Event.h:56
enum BeamMode string
void novaddt::TrackFit::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 
)
privatevirtual

Definition at line 314 of file TrackFit_module.cc.

References std::abs(), d, MECModelEnuComparisons::i, calib::j, and std::sqrt().

Referenced by produce().

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
T sqrt(T number)
Definition: d0nt_math.hpp:156
float abs(float number)
Definition: d0nt_math.hpp:39
Float_t d
Definition: plot.C:236
const double j
Definition: BetheBloch.cxx:29
Float_t w
Definition: plot.C:20
void art::Consumer::showMissingConsumes ( ) const
protectedinherited

Referenced by art::RootOutput::endJob().

double novaddt::TrackFit::UtilLinFit ( const std::vector< double > &  x,
const std::vector< double > &  y,
const std::vector< double > &  w,
double &  m,
double &  c 
)
privatevirtual

Definition at line 652 of file TrackFit_module.cc.

References ana::assert(), chi2(), d, DEFINE_ART_MODULE(), and MECModelEnuComparisons::i.

Referenced by LinFit().

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 }
double chi2()
Float_t d
Definition: plot.C:236
assert(nhit_max >=nhit_nbins)
Float_t w
Definition: plot.C:20
void art::Consumer::validateConsumedProduct ( BranchType const  bt,
ProductInfo const &  pi 
)
protectedinherited

Member Data Documentation

std::string novaddt::TrackFit::_assn_to_slice_instance
private

Definition at line 84 of file TrackFit_module.cc.

Referenced by produce(), and TrackFit().

double novaddt::TrackFit::fDHitGood
private

Definition at line 80 of file TrackFit_module.cc.

Referenced by produce(), and TrackFit().

std::vector<double> novaddt::TrackFit::fR
private

Range of weight function.

Definition at line 82 of file TrackFit_module.cc.

Referenced by FitView(), and TrackFit().

std::string novaddt::TrackFit::fSliceInstance
private

Definition at line 79 of file TrackFit_module.cc.

Referenced by produce().

std::string novaddt::TrackFit::fSliceLabel
private

Definition at line 78 of file TrackFit_module.cc.

Referenced by produce().

bool novaddt::TrackFit::fSortOutputHits
private

Definition at line 81 of file TrackFit_module.cc.

Referenced by produce().

std::vector<double> novaddt::TrackFit::fT
private

Transition width of weight function.

Definition at line 83 of file TrackFit_module.cc.

Referenced by FitView(), and TrackFit().

int novaddt::TrackFit::n_events = 0
private

Definition at line 85 of file TrackFit_module.cc.

Referenced by endJob(), and produce().

int novaddt::TrackFit::n_slices_input = 0
private

Definition at line 86 of file TrackFit_module.cc.

Referenced by endJob(), and produce().

int novaddt::TrackFit::n_tracks_output = 0
private

Definition at line 87 of file TrackFit_module.cc.

Referenced by endJob(), and produce().


The documentation for this class was generated from the following file: