LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
trkf::Track3DKalmanSPS Class Reference
Inheritance diagram for trkf::Track3DKalmanSPS:
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

 Track3DKalmanSPS (fhicl::ParameterSet const &pset)
 
virtual ~Track3DKalmanSPS ()
 
void produce (art::Event &evt)
 
void beginJob ()
 
void endJob ()
 
void reconfigure (fhicl::ParameterSet const &p)
 
double energyLossBetheBloch (const double &mass, const double p)
 
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< T > consumes (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< T > consumes (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< T > consumesView (InputTag const &it)
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< T > mayConsume (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< T > mayConsumeView (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

void rotationCov (TMatrixT< Double_t > &cov, const TVector3 &u, const TVector3 &v)
 
std::vector< double > dQdxCalc (const art::FindManyP< recob::Hit > &h, const art::PtrVector< recob::SpacePoint > &s, const TVector3 &p, const TVector3 &d)
 

Private Attributes

std::string fClusterModuleLabel
 
std::string fSpptModuleLabel
 
std::string fGenieGenModuleLabel
 
std::string fG4ModuleLabel
 
std::string fSortDim
 
TTree * tree
 
TMatrixT< Double_t > * stMCT
 
TMatrixT< Double_t > * covMCT
 
TMatrixT< Double_t > * stREC
 
TMatrixT< Double_t > * covREC
 
Float_t chi2
 
Float_t chi2ndf
 
int fcont
 
Float_t * fpRECL
 
Float_t * fpREC
 
Float_t * fpMCMom
 
Float_t * fpMCPos
 
Float_t * fState0
 
Float_t * fCov0
 
int nfail
 
int ndf
 
int nchi2rePass
 
int ispptvec
 
int nspptvec
 
unsigned int evtt
 
unsigned int nTrks
 
unsigned int fptsNo
 
Float_t * fshx
 
Float_t * fshy
 
Float_t * fshz
 
Float_t * feshx
 
Float_t * feshy
 
Float_t * feshz
 
Float_t * feshyz
 
Float_t * fupdate
 
Float_t * fchi2hit
 
Float_t * fth
 
Float_t * feth
 
Float_t * fedudw
 
Float_t * fedvdw
 
Float_t * feu
 
Float_t * fev
 
Float_t * fsep
 
Float_t * fdQdx
 
unsigned int fDimSize
 
Float_t * fPCmeans
 
Float_t * fPCevals
 
Float_t * fPCsigmas
 
Float_t * fPC1
 
Float_t * fPC2
 
Float_t * fPC3
 
std::vector< double > fPosErr
 
std::vector< double > fMomErr
 
std::vector< double > fMomStart
 
double fPerpLim
 
bool fDoFit
 
int fNumIt
 
uint16_t fMinNumSppts
 
double fErrScaleS
 
double fErrScaleM
 
int fDecimate
 
double fMaxUpdate
 
int fDecimateU
 
double fDistanceU
 
double fMaxUpdateU
 
double fMomLow
 
double fMomHigh
 
int fPdg
 
double fChi2Thresh
 
int fMaxPass
 
genf::GFAbsTrackReprepMC
 
genf::GFAbsTrackReprep
 

Detailed Description

Definition at line 98 of file Track3DKalmanSPS_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

trkf::Track3DKalmanSPS::Track3DKalmanSPS ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 207 of file Track3DKalmanSPS_module.cc.

References reconfigure().

208  : fDoFit(true)
209  , fNumIt(5)
210  , fMinNumSppts(5)
211  , fErrScaleS(1.0)
212  , fErrScaleM(1.0)
213  , fDecimate(1)
214  , fMaxUpdate(0.10)
215  , fDecimateU(1)
216  , fDistanceU(1.)
217  , fMaxUpdateU(0.10)
218  , fMomLow(0.001)
219  , fMomHigh(100.)
220  , fPdg(-13)
221  , fChi2Thresh(12.0E12)
222  , fMaxPass (1)
223  {
224 
225  this->reconfigure(pset);
226 
227  produces< std::vector<recob::Track> >();
228  produces<art::Assns<recob::Track, recob::Cluster> >();
229  produces<art::Assns<recob::Track, recob::SpacePoint> >();
230  produces<art::Assns<recob::Track, recob::Hit> >();
231 
232  }
void reconfigure(fhicl::ParameterSet const &p)
trkf::Track3DKalmanSPS::~Track3DKalmanSPS ( )
virtual

Definition at line 279 of file Track3DKalmanSPS_module.cc.

280  {
281  }

Member Function Documentation

void trkf::Track3DKalmanSPS::beginJob ( )
virtual

Reimplemented from art::EDProducer.

Definition at line 421 of file Track3DKalmanSPS_module.cc.

References chi2, chi2ndf, covMCT, covREC, evtt, fchi2hit, fcont, fCov0, fDimSize, fdQdx, fedudw, fedvdw, feshx, feshy, feshyz, feshz, feth, feu, fev, fPC1, fPC2, fPC3, fPCevals, fPCmeans, fPCsigmas, fpMCMom, fpMCPos, fpREC, fpRECL, fptsNo, fsep, fshx, fshy, fshz, fState0, fth, fupdate, ispptvec, art::TFileDirectory::make(), nchi2rePass, ndf, nfail, nspptvec, nTrks, stMCT, stREC, and tree.

422  {
423 
424 
426 
427  stMCT = new TMatrixT<Double_t>(5,1);
428  covMCT = new TMatrixT<Double_t>(5,5);
429  stREC = new TMatrixT<Double_t>(5,1);
430  covREC = new TMatrixT<Double_t>(5,5);
431 
432  fpMCMom = new Float_t[4];
433  fpMCPos = new Float_t[4];
434  fpREC = new Float_t[4];
435  fpRECL = new Float_t[4];
436  fState0 = new Float_t[5];
437  fCov0 = new Float_t[25];
438  fDimSize = 20000; // if necessary will get this from pset in constructor.
439 
440  fshx = new Float_t[fDimSize];
441  fshy = new Float_t[fDimSize];
442  fshz = new Float_t[fDimSize];
443  feshx = new Float_t[fDimSize];
444  feshy = new Float_t[fDimSize];
445  feshz = new Float_t[fDimSize];
446  feshyz = new Float_t[fDimSize];
447  fupdate = new Float_t[fDimSize];
448  fchi2hit = new Float_t[fDimSize];
449  fth = new Float_t[fDimSize];
450  feth = new Float_t[fDimSize];
451  fedudw = new Float_t[fDimSize];
452  fedvdw = new Float_t[fDimSize];
453  feu = new Float_t[fDimSize];
454  fev = new Float_t[fDimSize];
455  fsep = new Float_t[fDimSize];
456  fdQdx = new Float_t[fDimSize];
457 
458  fPC1 = new Float_t[3];
459  fPC2 = new Float_t[3];
460  fPC3 = new Float_t[3];
461  fPCmeans = new Float_t[3];
462  fPCsigmas = new Float_t[3];
463  fPCevals = new Float_t[3];
464 
465 // //TFile fileGENFIT("GENFITout.root","RECREATE");
466 
467 
468  tree = tfs->make<TTree>("GENFITttree","GENFITttree");
469  //tree->Branch("stMCT",&stMCT,"stMCT[5]/F"); // "TMatrixT<Double_t>"
470 
471  //tree->Branch("stMCT","TMatrixD",&stMCT,64000,0);
472  //tree->Branch("covMCT",covMCT,"covMCT[25]/F");
473  tree->Branch("covMCT","TMatrixD",&covMCT,64000,0);
474  tree->Branch("stREC",fState0,"stREC[5]/F");
475  //tree->Branch("stREC","TMatrixD",&stREC,64000,0);
476  tree->Branch("covREC",fCov0,"covREC[25]/F");
477  //tree->Branch("covREC","TMatrixD",&covREC,64000,0);
478 
479  tree->Branch("nchi2rePass",&nchi2rePass,"nchi2rePass/I");
480  tree->Branch("ispptvec",&ispptvec,"ispptvec/I");
481  tree->Branch("chi2",&chi2,"chi2/F");
482  tree->Branch("nfail",&nfail,"nfail/I");
483  tree->Branch("ndf",&ndf,"ndf/I");
484  tree->Branch("evtNo",&evtt,"evtNo/I");
485  tree->Branch("nspptvec",&nspptvec,"nspptvec/I");
486  tree->Branch("chi2ndf",&chi2ndf,"chi2ndf/F");
487 
488  tree->Branch("trkNo",&nTrks,"trkNo/I");
489  tree->Branch("ptsNo",&fptsNo,"ptsNo/I");
490  tree->Branch("cont",&fcont,"cont/I"); //O? Yes, O. Not 0, not L, ...
491  tree->Branch("shx",fshx,"shx[ptsNo]/F");
492  tree->Branch("shy",fshy,"shy[ptsNo]/F");
493  tree->Branch("shz",fshz,"shz[ptsNo]/F");
494  tree->Branch("sep",fsep,"sep[ptsNo]/F");
495  tree->Branch("dQdx",fdQdx,"dQdx[ptsNo]/F");
496  tree->Branch("eshx",feshx,"eshx[ptsNo]/F");
497  tree->Branch("eshy",feshy,"eshy[ptsNo]/F");
498  tree->Branch("eshz",feshz,"eshz[ptsNo]/F");
499  tree->Branch("eshyz",feshyz,"eshyz[ptsNo]/F");
500  tree->Branch("update",fupdate,"update[ptsNo]/F");
501  tree->Branch("chi2hit",fchi2hit,"chi2hit[ptsNo]/F");
502  tree->Branch("th",fth,"th[ptsNo]/F");
503  tree->Branch("eth",feth,"eth[ptsNo]/F");
504  tree->Branch("edudw",fedudw,"edudw[ptsNo]/F");
505  tree->Branch("edvdw",fedvdw,"edvdw[ptsNo]/F");
506  tree->Branch("eu",feu,"eu[ptsNo]/F");
507  tree->Branch("ev",fev,"ev[ptsNo]/F");
508 
509 
510  tree->Branch("pcMeans", fPCmeans,"pcMeans[3]/F");
511  tree->Branch("pcSigmas",fPCsigmas,"pcSigmas[3]/F");
512  tree->Branch("pcEvals", fPCevals,"pcEvals[3]/F");
513  tree->Branch("pcEvec1",fPC1,"pcEvec1[3]/F");
514  tree->Branch("pcEvec2",fPC2,"pcEvec2[3]/F");
515  tree->Branch("pcEvec3",fPC3,"pcEvec3[3]/F");
516 
517 
518  tree->Branch("pMCMom",fpMCMom,"pMCMom[4]/F");
519  tree->Branch("pMCPos",fpMCPos,"pMCPos[4]/F");
520  tree->Branch("pRECKalF",fpREC,"pRECKalF[4]/F");
521  tree->Branch("pRECKalL",fpRECL,"pRECKalL[4]/F");
522 
523 
524  //TGeoManager* geomGENFIT = new TGeoManager("Geometry", "Geane geometry");
525  //TGeoManager::Import("config/genfitGeom.root");
526  // gROOT->Macro("config/Geane.C");
527 
528  }
TMatrixT< Double_t > * stREC
TMatrixT< Double_t > * stMCT
TMatrixT< Double_t > * covREC
TMatrixT< Double_t > * covMCT
T * make(ARGS...args) const
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 147 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

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

Definition at line 162 of file Consumer.h.

163 {
164  if (!moduleContext_)
165  return;
166 
167  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
168 }
ConsumableProducts consumables_
Definition: Consumer.h:138
bool moduleContext_
Definition: Consumer.h:136
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 172 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

173 {
174  if (!moduleContext_)
175  return ViewToken<T>::invalid();
176 
177  consumables_[BT].emplace_back(ConsumableType::ViewElement,
178  TypeID{typeid(T)},
179  it.label(),
180  it.instance(),
181  it.process());
182  return ViewToken<T>{it};
183 }
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:138
bool moduleContext_
Definition: Consumer.h:136
EngineCreator::base_engine_t & EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make 
)
inherited

Definition at line 32 of file EngineCreator.cc.

References art::EngineCreator::rng().

34 {
35  return rng()->createEngine(
36  placeholder_schedule_id(), seed, kind_of_engine_to_make);
37 }
long seed
Definition: chem4.cc:68
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
EngineCreator::base_engine_t & EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make,
label_t const &  engine_label 
)
inherited

Definition at line 40 of file EngineCreator.cc.

References art::EngineCreator::rng().

43 {
44  return rng()->createEngine(
45  placeholder_schedule_id(), seed, kind_of_engine_to_make, engine_label);
46 }
long seed
Definition: chem4.cc:68
static art::ServiceHandle< art::RandomNumberGenerator > & rng()
CurrentProcessingContext const * art::EDProducer::currentContext ( ) const
protectedinherited

Definition at line 120 of file EDProducer.cc.

References art::EDProducer::current_context_.

121  {
122  return current_context_.get();
123  }
CPC_exempt_ptr current_context_
Definition: EDProducer.h:116
std::vector< double > trkf::Track3DKalmanSPS::dQdxCalc ( const art::FindManyP< recob::Hit > &  h,
const art::PtrVector< recob::SpacePoint > &  s,
const TVector3 &  p,
const TVector3 &  d 
)
private

Definition at line 350 of file Track3DKalmanSPS_module.cc.

References art::PtrVector< T >::begin(), art::PtrVector< T >::end(), recob::Hit::Integral(), geo::kCollection, geo::PlaneID::Plane, geo::GeometryCore::Plane(), recob::Hit::SignalType(), geo::WireGeo::ThetaZ(), lar::dump::vector(), geo::PlaneGeo::Wire(), recob::Hit::WireID(), and geo::GeometryCore::WirePitch().

Referenced by produce().

351  {
352  // For now just Collection plane.
353  // We should loop over all views, more generally.
357  // art::PtrVector<recob::SpacePoint>::const_iterator sppt = sstart;
359  std::vector <double> v;
360 
361 
362 
363  double mindist(100.0); // cm
364  auto spptminIt(sppt);
365  while (sppt != s.end())
366  {
367  if (((**sppt).XYZ() - loc).Mag() < mindist)
368  {
369  double dist = ((**sppt).XYZ() - loc).Mag();
370 
371  // Jump out if we're as close as 0.1 mm away.
372  if (dist<mindist)
373  {
374  mindist = dist;
375  spptminIt = sppt;
376  if (mindist < 0.01) break;
377  }
378  }
379  sppt++;
380  }
381  sstart = spptminIt; // for next time.
382  unsigned int ind(std::distance(s.begin(),spptminIt));
383 
384 
385 
386  std::vector< art::Ptr<recob::Hit> > hitlist = h.at(ind);
387 
388  double wirePitch = 0.;
389  double angleToVert = 0;
390  // unsigned int tpc1;
391  unsigned int plane1;
392  double charge = 0.;
393 
394  for(std::vector< art::Ptr<recob::Hit> >::const_iterator ihit = hitlist.begin();
395  ihit != hitlist.end(); ++ihit)
396  {
397  const recob::Hit& hit1 = **ihit;
398  // if (hit1.View() != view) continue;
399  if (hit1.SignalType() != sig) continue;
400  geo::WireID hit1WireID = hit1.WireID();
401  // tpc1 = hit1WireID.TPC;
402  plane1 = hit1WireID.Plane;
403  charge = hit1.Integral();
404  wirePitch = geom->WirePitch(plane1);
405  angleToVert = geom->Plane(plane1).Wire(0).ThetaZ(false) - 0.5*TMath::Pi();
406  }
407 
408  double cosgamma = TMath::Abs(TMath::Sin(angleToVert)*dir.Y() +
409  TMath::Cos(angleToVert)*dir.Z());
410  // if(cosgamma < 1.e-5)
411  // throw cet::exception("Track") << "cosgamma is basically 0, that can't be right\n";
412 
413  v.push_back(charge/wirePitch/cosgamma);
414  // std::cout << " Track3DKalmanSPS::dQdxCalc() : For loc.XYZ() hit is ... " << ind << " and v is " << v.back() << std::endl;
415 
416  return v;
417 
418  }
geo::SigType_t SignalType() const
Signal type for the plane of the hit.
Definition: Hit.h:232
PlaneGeo const & Plane(unsigned int const p, unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified wire.
WireGeo const & Wire(unsigned int iwire) const
Definition: PlaneGeo.cxx:506
geo::WireID WireID() const
Initial tdc tick for hit.
Definition: Hit.h:234
iterator begin()
Definition: PtrVector.h:223
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:225
double ThetaZ() const
Returns angle of wire with respect to z axis in the Y-Z plane in radians.
Definition: WireGeo.h:192
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
enum geo::_plane_sigtype SigType_t
Enumerate the possible plane projections.
iterator end()
Definition: PtrVector.h:237
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:258
data_t::const_iterator const_iterator
Definition: PtrVector.h:61
TDirectory * dir
Definition: macro.C:5
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:49
Signal from collection planes.
Definition: geo_types.h:93
void trkf::Track3DKalmanSPS::endJob ( )
virtual

Reimplemented from art::EDProducer.

Definition at line 531 of file Track3DKalmanSPS_module.cc.

References fchi2hit, fCov0, fdQdx, fedudw, fedvdw, feshx, feshy, feshyz, feshz, feth, feu, fev, fPC1, fPC2, fPC3, fPCevals, fPCmeans, fPCsigmas, fpREC, fpRECL, fsep, fshx, fshy, fshz, fState0, fth, fupdate, rep, and repMC.

532  {
533  if (!rep) delete rep;
534  if (!repMC) delete repMC;
535 
536  /*
537  // not sure why I can't do these, but at least some cause seg faults.
538  delete[] stMCT;
539  delete[] covMCT;
540  delete[] stREC;
541  delete[] covREC;
542  */
543 
544  delete[] fpREC;
545  delete[] fpRECL;
546  delete[] fState0;
547  delete[] fCov0;
548 
549  delete[] fshx;
550  delete[] fshy;
551  delete[] fshz;
552  delete[] feshx;
553  delete[] feshy;
554  delete[] feshyz;
555  delete[] feshz;
556  delete[] fupdate;
557  delete[] fchi2hit;
558  delete[] fth;
559  delete[] feth;
560  delete[] fedudw;
561  delete[] fedvdw;
562  delete[] feu;
563  delete[] fev;
564  delete[] fsep;
565  delete[] fdQdx;
566 
567  delete[] fPCmeans;
568  delete[] fPCsigmas;
569  delete[] fPCevals;
570  delete[] fPC1;
571  delete[] fPC2;
572  delete[] fPC3;
573 
574  }
double trkf::Track3DKalmanSPS::energyLossBetheBloch ( const double &  mass,
const double  p = 1.5 
)

Definition at line 285 of file Track3DKalmanSPS_module.cc.

References beta.

Referenced by produce().

288  {
289  const double charge(1.0);
290  const double mEE(188.); // eV
291  const double matZ(18.);
292  const double matA(40.);
293  const double matDensity(1.4);
294  const double me(0.000511);
295 
296  double beta = p/std::sqrt(mass*mass+p*p);
297  double gammaSquare = 1./(1.0 - beta*beta);
298  // 4pi.r_e^2.N.me = 0.307075, I think.
299  double dedx = 0.307075*matDensity*matZ/matA/(beta*beta)*charge*charge;
300  double massRatio = me/mass;
301  // me=0.000511 here is in GeV. So mEE comes in here in eV.
302  double argument = gammaSquare*beta*beta*me*1.E3*2./((1.E-6*mEE) * std::sqrt(1+2*std::sqrt(gammaSquare)*massRatio + massRatio*massRatio));
303 
304  if (mass==0.0) return(0.0);
305  if (argument <= exp(beta*beta))
306  {
307  dedx = 0.;
308  }
309  else{
310  dedx *= (log(argument)-beta*beta); // Bethe-Bloch [MeV/cm]
311  dedx *= 1.E-3; // in GeV/cm, hence 1.e-3
312  if (dedx<0.) dedx = 0.;
313  }
314  return dedx;
315  }
Double_t beta
EngineCreator::seed_t EngineCreator::get_seed_value ( fhicl::ParameterSet const &  pset,
char const  key[] = "seed",
seed_t const  implicit_seed = -1 
)
inherited

Definition at line 49 of file EngineCreator.cc.

References fhicl::ParameterSet::get().

Referenced by art::MixFilter< T >::initEngine_().

52 {
53  auto const& explicit_seeds = pset.get<std::vector<int>>(key, {});
54  return explicit_seeds.empty() ? implicit_seed : explicit_seeds.front();
55 }
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_.

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 B, and 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  }
Int_t B
Definition: plot.C:25
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 190 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

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

Definition at line 205 of file Consumer.h.

206 {
207  if (!moduleContext_)
208  return;
209 
210  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
211 }
ConsumableProducts consumables_
Definition: Consumer.h:138
bool moduleContext_
Definition: Consumer.h:136
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 215 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

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

Definition at line 40 of file ProducerBase.h.

References art::ProducerBase::getProductID().

41  {
42  return true;
43  }
void art::Consumer::prepareForJob ( fhicl::ParameterSet const &  pset)
protectedinherited

Definition at line 89 of file Consumer.cc.

References fhicl::ParameterSet::get_if_present().

Referenced by art::EDProducer::doBeginJob(), art::EDFilter::doBeginJob(), and art::EDAnalyzer::doBeginJob().

90 {
91  if (!moduleContext_)
92  return;
93 
94  pset.get_if_present("errorOnMissingConsumes", requireConsumes_);
95  for (auto& consumablesPerBranch : consumables_) {
96  cet::sort_all(consumablesPerBranch);
97  }
98 }
bool requireConsumes_
Definition: Consumer.h:137
ConsumableProducts consumables_
Definition: Consumer.h:138
bool moduleContext_
Definition: Consumer.h:136
void trkf::Track3DKalmanSPS::produce ( art::Event evt)
virtual
Todo:
Should never test whether the event is real data in reconstruction algorithms
Todo:
as that introduces potential data/MC differences that are very hard to track down
Todo:
Remove this test as soon as possible please
Todo:
Should never test whether the event is real data in reconstruction algorithms
Todo:
as that introduces potential data/MC differences that are very hard to track down
Todo:
Remove this test as soon as possible please

Implements art::EDProducer.

Definition at line 578 of file Track3DKalmanSPS_module.cc.

References genf::GFTrack::addHit(), art::PtrVector< T >::begin(), chi2, chi2ndf, close(), recob::tracking::convertCollToPoint(), recob::tracking::convertCollToVector(), covMCT, covREC, util::CreateAssn(), DEFINE_ART_MODULE, geo::GeometryCore::DetHalfHeight(), geo::GeometryCore::DetHalfWidth(), geo::GeometryCore::DetLength(), dQdxCalc(), e, art::PtrVector< T >::end(), energyLossBetheBloch(), art::PtrVector< T >::erase(), art::EventID::event(), evtt, fchi2hit, fChi2Thresh, fClusterModuleLabel, fcont, fCov0, fDecimate, fDecimateU, fDimSize, fDistanceU, fDoFit, fdQdx, fedudw, fedvdw, fErrScaleM, fErrScaleS, feshx, feshy, feshyz, feshz, feth, feu, fev, fGenieGenModuleLabel, fMaxPass, fMaxUpdate, fMaxUpdateU, fMinNumSppts, fMomErr, fMomHigh, fMomLow, fMomStart, fNumIt, fPC1, fPC2, fPC3, fPCevals, fPCmeans, fPCsigmas, fPdg, fPerpLim, fpMCMom, fpMCPos, fPosErr, fpREC, fpRECL, fptsNo, fsep, fshx, fshy, fshz, fSortDim, fSpptModuleLabel, fState0, fth, fupdate, art::DataViewImpl::getByLabel(), genf::GFAbsTrackRep::getChiSqu(), genf::GFAbsTrackRep::getCov(), genf::GFTrack::getFailedHits(), genf::GFTrack::getHitChi2(), genf::GFTrack::getHitCov(), genf::GFTrack::getHitCov7x7(), genf::GFTrack::getHitMeasuredCov(), genf::GFTrack::getHitPlaneU(), genf::GFTrack::getHitPlaneUxUyUz(), genf::GFTrack::getHitPlaneV(), genf::GFTrack::getHitPlaneXYZ(), genf::GFTrack::getHitState(), genf::GFTrack::getHitUpdate(), genf::GFFieldManager::getInstance(), geo::GeometryCore::GetLArTPCVolumeName(), genf::GFAbsTrackRep::getLastPlane(), genf::GFAbsTrackRep::getMom(), genf::GFAbsTrackRep::getNDF(), simb::MCTruth::GetParticle(), genf::GFAbsTrackRep::getReferencePlane(), genf::GFAbsTrackRep::getState(), genf::GFAbsTrackRep::getStatusFlag(), hits(), art::Event::id(), genf::GFFieldManager::init(), art::PtrVector< T >::insert(), mf::isDebugEnabled(), ispptvec, art::Event::isRealData(), LOG_DEBUG, LOG_ERROR, nchi2rePass, ndf, nfail, simb::MCTruth::NParticles(), nspptvec, nTrks, part, genf::GFDetPlane::Print(), genf::GFKalman::processTrack(), art::PtrVector< T >::push_back(), art::Event::put(), rep, repMC, genf::ROOTobjectToString(), rotationCov(), art::Event::run(), genf::GFKalman::setBlowUpFactor(), genf::GFKalman::setErrorScaleMTh(), genf::GFKalman::setErrorScaleSTh(), genf::GFKalman::setInitialDirection(), genf::GFKalman::setMaxUpdate(), genf::GFKalman::setMomHigh(), genf::GFKalman::setMomLow(), genf::GFKalman::setNumIterations(), genf::GFTrack::setPDG(), art::PtrVector< T >::size(), sp_sort_3dx(), sp_sort_3dy(), sp_sort_3dz(), sp_sort_nsppts(), stMCT, stREC, tmp, tree, and w.

579 {
580 
581  rep=0;
582  repMC=0;
583  // get services
585 
587  // Make a std::unique_ptr<> for the thing you want to put into the event
588  // because that handles the memory management for you
590  std::unique_ptr<std::vector<recob::Track> > tcol(new std::vector<recob::Track>);
591  std::unique_ptr< art::Assns<recob::Track, recob::SpacePoint> > tspassn(new art::Assns<recob::Track, recob::SpacePoint>);
592  std::unique_ptr< art::Assns<recob::Track, recob::Hit> > thassn(new art::Assns<recob::Track, recob::Hit>);
593 
594  unsigned int tcnt = 0;
595 
596  // define TPC parameters
597  TString tpcName = geom->GetLArTPCVolumeName();
598 
599 
600  // get input Hit object(s).
601  art::Handle< std::vector<recob::Cluster> > clusterListHandle;
602  evt.getByLabel(fClusterModuleLabel,clusterListHandle);
603 
605  evt.getByLabel(fSpptModuleLabel,spptListHandle);
606 
608 
612  if (!evt.isRealData()){
613 
614  // std::cout << "Track3DKalmanSPS: This is MC." << std::endl;
615  // std::cout<<"Run "<<evt.run()<<" Event "<<evt.id().event()<<std::endl;
616 
617  art::Handle< std::vector<simb::MCTruth> > mctruthListHandle;
618  evt.getByLabel(fGenieGenModuleLabel,mctruthListHandle);
619 
620  for (unsigned int ii = 0; ii < mctruthListHandle->size(); ++ii)
621  {
622  art::Ptr<simb::MCTruth> mctparticle(mctruthListHandle,ii);
623  mclist.push_back(mctparticle);
624  }
625 
626  }
627 
628 
629  // std::cout<<"Run "<<evt.run()<<" Event "<<evt.id().event()<<std::endl;
630 
631  // Put this back when Wes's reign of terror ends ...
632  // LOG_DEBUG("Track3DKalmanSPS") << "There are " << spptListHandle->size() << " Spacepoint PtrVectors (spacepoint clumps) in this event.";
633 
634  std::vector < art::PtrVector<recob::SpacePoint> > spptIn(spptListHandle->begin(),spptListHandle->end());
635  // Get the spptvectors that are largest to be first, and smallest last.
636  std::sort(spptIn.begin(), spptIn.end(), sp_sort_nsppts);
637 
638  std::vector < art::PtrVector<recob::SpacePoint> >::const_iterator sppt = spptIn.begin();
639  auto spptB = sppt;
640 
641 
642  TVector3 MCOrigin;
643  TVector3 MCMomentum;
644  // TVector3 posErr(.05,.05,.05); // resolution. 0.5mm
645  // TVector3 momErr(.1,.1,0.2); // GeV
646  TVector3 posErr(fPosErr[0],fPosErr[1],fPosErr[2]); // resolution. 0.5mm
647  TVector3 momErr(fMomErr[0],fMomErr[1],fMomErr[2]); // GeV
648  TVector3 momErrFit(fMomErr[0],fMomErr[1],fMomErr[2]); // GeV
649 
650  // This is strictly for MC
654  if (!evt.isRealData())
655  {
656  // Below breaks are stupid, I realize. But rather than keep all the MC
657  // particles I just take the first primary, e.g., muon and only keep its
658  // info in the branches of the Ttree. I could generalize later, ...
659  for( unsigned int ii = 0; ii < mclist.size(); ++ii )
660  {
661  //art::Ptr<const simb::MCTruth> mc(mctruthListHandle,i);
662  art::Ptr<simb::MCTruth> mc(mclist[ii]);
663  for(int jj = 0; jj < mc->NParticles(); ++jj)
664  {
665  simb::MCParticle part(mc->GetParticle(jj));
666  MCOrigin.SetXYZ(part.Vx(),part.Vy(),part.Vz()); // V for Vertex
667  MCMomentum.SetXYZ(part.Px(),part.Py(),part.Pz());
668  LOG_DEBUG("Track3DKalmanSPS_GenFit")
669  << "FROM MC TRUTH, the particle's pdg code is: "<<part.PdgCode()<< " with energy = "<<part.E() <<", with energy = "<<part.E()
670  << "\n vtx: " << genf::ROOTobjectToString(MCOrigin)
671  << "\n momentum: " << genf::ROOTobjectToString(MCMomentum)
672  << "\n (both in Global (not volTPC) coords)";
673 
674  repMC = new genf::RKTrackRep(MCOrigin,
675  MCMomentum,
676  posErr,
677  momErr,
678  part.PdgCode());
679  break;
680  }
681  break;
682  }
683  //for saving of MC truth
684  stMCT->ResizeTo(repMC->getState());
685  *stMCT = repMC->getState();
686  covMCT-> ResizeTo(repMC->getCov());
687  *covMCT = repMC->getCov();
688  LOG_DEBUG("Track3DKalmanSPS_GenFit") << " repMC, covMC are ... \n"
691 
692  } // !isRealData
693  nTrks = 0;
694  TParticlePDG * part = TDatabasePDG::Instance()->GetParticle(fPdg);
695  Double_t mass = part->Mass();
696 
697 
698  size_t cntp(0);
699  while (sppt!=spptIn.end())
700  {
701 
702  const art::PtrVector<recob::SpacePoint> & spacepoints = *sppt;
703 
704  double fMaxUpdateHere(fMaxUpdateU);
705  int fDecimateHere(fDecimateU);
706  double fErrScaleSHere(fErrScaleS);
707  double fErrScaleMHere(fErrScaleM);
708  int rePass0(1);
709  unsigned int nTailPoints = 0; // 100;
710  if (spacepoints.size()<5)
711  { sppt++; rePass0 = 3; continue;} // for now...
712 
713  LOG_DEBUG("Track3DKalmanSPS_GenFit")
714  <<"\n\t found "<<spacepoints.size()<<" 3D spacepoint(s) for this element of std::vector<art:PtrVector> spacepoints. \n";
715 
716  //const double resolution = posErr.Mag();
717  //
718 
719 
720  // Let's find track's principle components.
721  // We will sort along that direction, rather than z.
722  // Further, we will skip outliers away from main axis.
723 
724  TPrincipal* principal = new TPrincipal(3,"ND");
725 
726  // I need to shuffle these around, so use copy constructor
727  // to make non-const version spacepointss.
728  art::PtrVector<recob::SpacePoint> spacepointss(spacepoints);
729 
730  // What I need is a nearest neighbor sorting.
731  if (fSortDim.compare("y") && fSortDim.compare("x")) std::sort(spacepointss.begin(), spacepointss.end(), sp_sort_3dz);
732  if (!fSortDim.compare("y")) std::sort(spacepointss.begin(), spacepointss.end(), sp_sort_3dy);
733  if (!fSortDim.compare("x")) std::sort(spacepointss.begin(), spacepointss.end(), sp_sort_3dx);
734 
735  for (unsigned int point=0;point<spacepointss.size();++point)
736  {
737  // std::cout << "Spacepoint " << point << " added:" << spacepointss[point]->XYZ()[0]<< ", " << spacepointss[point]->XYZ()[1]<< ", " << spacepointss[point]->XYZ()[2]<< ". " << std::endl;
738  if (point<(spacepointss.size()-nTailPoints))
739  {
740  principal->AddRow(spacepointss[point]->XYZ());
741  }
742  }
743  principal->MakePrincipals();
744  /*
745  principal->Test();
746  principal->MakeHistograms();
747  principal->Print("MSEV");
748  */
749  const TVectorD* evals = principal->GetEigenValues();
750  const TMatrixD* evecs = principal->GetEigenVectors();
751  const TVectorD* means = principal->GetMeanValues();
752  const TVectorD* sigmas = principal->GetSigmas();
753  /*
754  std::vector<TVector3*> pcs;
755  Double_t* pc = new Double_t[3];
756  for (unsigned int point=0;point<spacepointss.size();++point)
757  {
758  principal->X2P((Double_t *)(spacepointss[point]->XYZ()),pc);
759  pcs.push_back((TVector3 *)pc); // !!!
760  }
761  delete [] pc;
762  */
763  Double_t tmp[3], tmp2[3];
764  principal->X2P((Double_t *)(means->GetMatrixArray()),tmp);
765  principal->X2P((Double_t *)(sigmas->GetMatrixArray()),tmp2);
766  for (unsigned int ii=0;ii<3;++ii)
767  {
768  fPCmeans[ii] = (Float_t )(tmp[ii]);
769  fPCsigmas[ii] = (Float_t )(tmp2[ii]);
770  fPCevals[ii] = (Float_t )(evals->GetMatrixArray())[ii];
771  // This method requires apparently pulling all 9
772  // elements. Maybe 3 works.
773  // Certainly, w can't be a scalar, I discovered.
774  double w[9];
775  evecs->ExtractRow(ii,0,w);
776  fPC1[ii] = w[0];
777  fPC2[ii] = w[1];
778  fPC3[ii] = w[2];
779  }
780  Double_t tmp3[3], tmp4[3], tmp5[3];
781  principal->X2P((Double_t *)fPC1,tmp3);
782  principal->X2P((Double_t *)fPC2,tmp4);
783  principal->X2P((Double_t *)fPC3,tmp5);
784 
785 
786  // Use a mip approximation assuming straight lines
787  // and a small angle wrt beam.
788  fMomStart[0] = spacepointss[spacepointss.size()-1]->XYZ()[0] - spacepointss[0]->XYZ()[0];
789  fMomStart[1] = spacepointss[spacepointss.size()-1]->XYZ()[1] - spacepointss[0]->XYZ()[1];
790  fMomStart[2] = spacepointss[spacepointss.size()-1]->XYZ()[2] - spacepointss[0]->XYZ()[2];
791  // This presumes a 0.8 GeV/c particle
792  double dEdx = energyLossBetheBloch(mass, 1.0);
793  // mom is really KE.
794  TVector3 mom(dEdx*fMomStart[0],dEdx*fMomStart[1],dEdx*fMomStart[2]);
795  double pmag2 = pow(mom.Mag()+mass, 2. - mass*mass);
796  mom.SetMag(std::sqrt(pmag2));
797  // Over-estimate by just enough for contained particles (5%).
798  mom.SetMag(1.0 * mom.Mag());
799  // My true 0.5 GeV/c muons need a yet bigger over-estimate.
800  //if (mom.Mag()<0.7) mom.SetMag(1.2*mom.Mag());
801  // if (mom.Mag()>2.0) mom.SetMag(10.0*mom.Mag());
802  // mom.SetMag(3*mom.Mag()); // EC, 15-Feb-2012. TEMPORARY!!!
803  // If 1st/last point is close to edge of TPC, this track is
804  // uncontained.Give higher momentum starting value in
805  // that case.
806  bool uncontained(false);
807  double close(5.); // cm.
808  double epsMag(0.001);// cm.
809  double epsX(250.0); // cm.
810  double epsZ(0.001); // cm.
811 
812  if (
813  spacepointss[spacepointss.size()-1]->XYZ()[0] > (2.*geom->DetHalfWidth(0,0)-close) || spacepointss[spacepointss.size()-1]->XYZ()[0] < close ||
814  spacepointss[0]->XYZ()[0] > (2.*geom->DetHalfWidth(0,0)-close) || spacepointss[0]->XYZ()[0] < close ||
815  spacepointss[spacepointss.size()-1]->XYZ()[1] > (1.*geom->DetHalfHeight(0,0)-close) || (spacepointss[spacepointss.size()-1]->XYZ()[1] < -1.*geom->DetHalfHeight(0,0)+close) ||
816  spacepointss[0]->XYZ()[1] > (1.*geom->DetHalfHeight(0,0)-close) || spacepointss[0]->XYZ()[1] < (-1.*geom->DetHalfHeight(0,0)+close) ||
817  spacepointss[spacepointss.size()-1]->XYZ()[2] > (geom->DetLength(0,0)-close) || spacepointss[spacepointss.size()-1]->XYZ()[2] < close ||
818  spacepointss[0]->XYZ()[2] > (geom->DetLength(0,0)-close) || spacepointss[0]->XYZ()[2] < close
819  )
820  uncontained = true;
821  fErrScaleSHere = fErrScaleS;
822  fErrScaleMHere = fErrScaleM;
823 
824  if (uncontained)
825  {
826  // Big enough to not run out of gas right at end of
827  // track and give large angular deviations which
828  // will kill the fit.
829  mom.SetMag(2.0 * mom.Mag());
830  LOG_DEBUG("Track3DKalmanSPS_GenFit")<<"Uncontained track ... ";
831  fDecimateHere = fDecimateU;
832  fMaxUpdateHere = fMaxUpdateU;
833  }
834  else
835  {
836  LOG_DEBUG("Track3DKalmanSPS_GenFit")<<"Contained track ... Run "<<evt.run()<<" Event "<<evt.id().event();
837  // Don't decimate contained tracks as drastically,
838  // and omit only very large corrections ...
839  // which hurt only high momentum tracks.
840  fDecimateHere = fDecimate;
841  fMaxUpdateHere = fMaxUpdate;
842  }
843  fcont = (int) (!uncontained);
844 
845  // This seems like best place to jump back to for a re-pass.
846  unsigned short rePass = rePass0; // 1 by default;
847  unsigned short maxPass(fMaxPass);
848  unsigned short tcnt1(0);
849  while (rePass<=maxPass)
850  {
851 
852  TVector3 momM(mom);
853  TVector3 momErrFit(momM[0]/3.0,
854  momM[1]/3.0,
855  momM[2]/3.0); // GeV
856 
858  genf::GFDetPlane planeG((TVector3)(spacepointss[0]->XYZ()),momM);
859 
860 
861  // std::cout<<"Track3DKalmanSPS about to do GAbsTrackRep."<<std::endl;
862  // Initialize with 1st spacepoint location and ...
863  rep = new genf::RKTrackRep(//posM-.5/momM.Mag()*momM,
864  (TVector3)(spacepointss[0]->XYZ()),
865  momM,
866  posErr,
867  momErrFit,
868  fPdg); // mu+ hypothesis
869  // std::cout<<"Track3DKalmanSPS: about to do GFTrack. repDim is " << rep->getDim() <<std::endl;
870 
871 
872  genf::GFTrack fitTrack(rep);//initialized with smeared rep
873  fitTrack.setPDG(fPdg);
874  // Gonna sort in z cuz I want to essentially transform here to volTPC coords.
875  // volTPC coords, cuz that's what the Geant3/Geane stepper wants, as that's its understanding
876  // from the Geant4 geometry, which it'll use. EC, 7-Jan-2011.
877  int ihit = 0;
878  fptsNo = 0;
879  std::vector <unsigned int> spptSurvivedIndex;
880  std::vector <unsigned int> spptSkippedIndex;
881  unsigned int ppoint(0);
882  for (unsigned int point=0;point<spacepointss.size();++point)
883  {
884  double sep;
885  // Calculate the distance in 2nd and 3rd PCs and
886  // reject spt if it's too far out. Remember, the
887  // sigmas are std::sqrt(eigenvals).
888  double tmp[3];
889  principal->X2P((Double_t *)(spacepointss[point]->XYZ()),tmp);
890  sep = std::sqrt(tmp[1]*tmp[1]/fPCevals[1]+tmp[2]*tmp[2]/fPCevals[2]);
891  if ((std::abs(sep) > fPerpLim) && (point<(spacepointss.size()-nTailPoints)) && rePass<=1)
892  {
893  // std::cout << "Spacepoint " << point << " DROPPED, cuz it's sufficiently far from the PCA major axis!!!:" << spacepointss[point]->XYZ()[0]<< ", " << spacepointss[point]->XYZ()[1]<< ", " << spacepointss[point]->XYZ()[2]<< ". " << std::endl;
894  spptSkippedIndex.push_back(point);
895  continue;
896  }
897  // If point is too close in Mag or Z or too far in X from last kept point drop it.
898  // I think this is largely redundant with PCA cut.
899  TVector3 one(spacepointss[point]->XYZ());
900  TVector3 two(spacepointss[ppoint]->XYZ());
901  if (rePass==2 && uncontained)
902  {
903  epsMag = fDistanceU; // cm
904  fNumIt = 4;
905  fErrScaleMHere = 0.1;
906  // Above allows us to pretend as though measurements
907  // are perfect, which we can ostensibly do now with
908  // clean set of sppts. This creates larger gains, bigger
909  // updates: bigger sensitivity to multiple scattering.
910 
911  // std::cout << "Spacepoint " << point << " ?DROPPED? magnitude and TV3 diff to ppoint is :" << (((TVector3)(spacepointss[point]->XYZ()-spacepointss[ppoint]->XYZ())).Mag()) << " and " << one[0] << ", " << one[1] << ", " << one[2] << two[0] << ", " << two[1] << ", " << two[2] << ". " << std::endl;
912  }
913  else if (rePass==2 && !uncontained)
914  {
915 
916  // fNumIt = 2;
917  // std::cout << "Spacepoint " << point << " ?DROPPED? magnitude and TV3 diff to ppoint is :" << (((TVector3)(spacepointss[point]->XYZ()-spacepointss[ppoint]->XYZ())).Mag()) << " and " << one[0] << ", " << one[1] << ", " << one[2] << two[0] << ", " << two[1] << ", " << two[2] << ". " << std::endl;
918  }
919  if (point>0 &&
920  (
921  (one-two).Mag()<epsMag || // too close
922  ((one-two).Mag()>8.0&&rePass==1) || // too far
923  std::abs(spacepointss[point]->XYZ()[2]-spacepointss[ppoint]->XYZ()[2])<epsZ ||
924  std::abs(spacepointss[point]->XYZ()[0]-spacepointss[ppoint]->XYZ()[0])>epsX
925  )
926  )
927  {
928  // std::cout << "Spacepoint " << point << " DROPPED, cuz it's too far in x or too close in magnitude or z to previous used spacepoint!!!:" << spacepointss[point]->XYZ()[0]<< ", " << spacepointss[point]->XYZ()[1]<< ", " << spacepointss[point]->XYZ()[2]<< ". " << std::endl;
929  // std::cout << "Prev used Spacepoint " << spacepointss[ppoint]->XYZ()[0]<< ", " << spacepointss[ppoint]->XYZ()[1]<< ", " << spacepointss[ppoint]->XYZ()[2]<< ". " << std::endl;
930  spptSkippedIndex.push_back(point);
931  continue;
932  }
933 
934  if (point%fDecimateHere && rePass<=1) // Jump out of loop except on every fDecimate^th pt. fDecimate==1 never sees continue.
935  {
936  /* Replace continue with a counter that will be used
937  to index into vector of GFKalman fits.
938  */
939  // spptSkippedIndex.push_back(point);
940  continue;
941  }
942 
943  ppoint=point;
944  TVector3 spt3 = (TVector3)(spacepointss[point]->XYZ());
945  std::vector <double> err3;
946  err3.push_back(spacepointss[point]->ErrXYZ()[0]);
947  err3.push_back(spacepointss[point]->ErrXYZ()[2]);
948  err3.push_back(spacepointss[point]->ErrXYZ()[4]);
949  err3.push_back(spacepointss[point]->ErrXYZ()[5]); // lower triangle diags.
950  if (fptsNo<fDimSize)
951  {
952  fshx[fptsNo] = spt3[0];
953  fshy[fptsNo] = spt3[1];
954  fshz[fptsNo] = spt3[2];
955  feshx[fptsNo] = err3[0];
956  feshy[fptsNo] = err3[1];
957  feshz[fptsNo] = err3[3];
958  feshyz[fptsNo] = err3[2];
959  fsep[fptsNo] = sep;
960 
961  if (fptsNo>1)
962  {
963  TVector3 pointer(fshx[fptsNo]-fshx[fptsNo-1],fshy[fptsNo]-fshy[fptsNo-1],fshz[fptsNo]-fshz[fptsNo-1]);
964  TVector3 pointerPrev(fshx[fptsNo-1]-fshx[fptsNo-2],fshy[fptsNo-1]-fshy[fptsNo-2],fshz[fptsNo-1]-fshz[fptsNo-2]);
965  fth[fptsNo] = (pointer.Unit()).Angle(pointerPrev.Unit());
966  }
967  feth[fptsNo] = 0.0;
968  fedudw[fptsNo] = 0.0;
969  fedvdw[fptsNo] = 0.0;
970  feu[fptsNo] = 0.0;
971  fev[fptsNo] = 0.0;
972  fupdate[fptsNo] = 0.0;
973  }
974 
975 
976  LOG_DEBUG("Track3DKalmanSPS_GenFit") << "ihit xyz..." << spt3[0]<<","<< spt3[1]<<","<< spt3[2];
977 
978  fitTrack.addHit(new genf::PointHit(spt3,err3),
979  1,//dummy detector id
980  ihit++
981  );
982  spptSurvivedIndex.push_back(point);
983  fptsNo++;
984  } // end loop over spacepoints.
985 
986  if (fptsNo<=fMinNumSppts) // Cuz 1st 2 in each direction don't count. Should have, say, 3 more.
987  {
988  LOG_DEBUG("Track3DKalmanSPS_GenFit") << "Bailing cuz only " << fptsNo << " spacepoints.";
989  rePass++;
990  continue;
991  }
992  LOG_DEBUG("Track3DKalmanSPS_GenFit") << "Fitting on " << fptsNo << " spacepoints.";
993  // std::cout<<"Track3DKalmanSPS about to do GFKalman."<<std::endl;
994  genf::GFKalman k;
995  k.setBlowUpFactor(5); // 500 out of box. EC, 6-Jan-2011.
996  k.setMomHigh(fMomHigh); // Don't fit above this many GeV.
997  k.setMomLow(fMomLow); // Don't fit below this many GeV.
998 
999  k.setInitialDirection(+1); // Instead of 1 out of box. EC, 6-Jan-2011.
1001  k.setMaxUpdate(fMaxUpdateHere); // 0 out abs(update) bigger than this.
1002  k.setErrorScaleSTh(fErrScaleSHere);
1003  k.setErrorScaleMTh(fErrScaleMHere);
1004 
1005  bool skipFill = false;
1006  // std::cout<<"Track3DKalmanSPS back from setNumIterations."<<std::endl;
1007  std::vector < TMatrixT<double> > hitMeasCov;
1008  std::vector < TMatrixT<double> > hitUpdate;
1009  std::vector < TMatrixT<double> > hitCov;
1010  std::vector < TMatrixT<double> > hitCov7x7;
1011  std::vector < TMatrixT<double> > hitState;
1012  std::vector < double > hitChi2;
1013  std::vector <TVector3> hitPlaneXYZ;
1014  std::vector <TVector3> hitPlaneUxUyUz;
1015  std::vector <TVector3> hitPlaneU;
1016  std::vector <TVector3> hitPlaneV;
1017 
1018  try{
1019  // std::cout<<"Track3DKalmanSPS about to processTrack."<<std::endl;
1020  if (fDoFit) k.processTrack(&fitTrack);
1021  //std::cout<<"Track3DKalmanSPS back from processTrack."<<std::endl;
1022  }
1023  //catch(GFException& e){
1024  catch(cet::exception &e){
1025  LOG_ERROR("Track3DKalmanSPS") << "just caught a cet::exception: " << e.what()
1026  << "\nExceptions won't be further handled; skip filling big chunks of the TTree.";
1027  skipFill = true;
1028  // exit(1);
1029  }
1030 
1031  if(rep->getStatusFlag()==0) // 0 is successful completion
1032  {
1033  if(mf::isDebugEnabled()) {
1034 
1035  std::ostringstream dbgmsg;
1036  dbgmsg << "Original plane:";
1037  planeG.Print(dbgmsg);
1038 
1039  dbgmsg << "Current (fit) reference Plane:";
1040  rep->getReferencePlane().Print(dbgmsg);
1041 
1042  dbgmsg << "Last reference Plane:";
1043  rep->getLastPlane().Print(dbgmsg);
1044 
1045  if (planeG != rep->getReferencePlane())
1046  dbgmsg <<" => original hit plane (not surprisingly) not current reference Plane!";
1047 
1048  LOG_DEBUG("Track3DKalmanSPS_GenFit") << dbgmsg.str();
1049  }
1050  if (!skipFill)
1051  {
1052  hitMeasCov = fitTrack.getHitMeasuredCov();
1053  hitUpdate = fitTrack.getHitUpdate();
1054  hitCov = fitTrack.getHitCov();
1055  hitCov7x7 = fitTrack.getHitCov7x7();
1056  hitState = fitTrack.getHitState();
1057  hitChi2 = fitTrack.getHitChi2();
1058  hitPlaneXYZ = fitTrack.getHitPlaneXYZ();
1059  hitPlaneUxUyUz = fitTrack.getHitPlaneUxUyUz();
1060  hitPlaneU = fitTrack.getHitPlaneU();
1061  hitPlaneV = fitTrack.getHitPlaneV();
1062  unsigned int totHits = hitState.size();
1063 
1064  // for (unsigned int ihit=0; ihit<fptsNo; ihit++)
1065  // Pick up info from last fwd Kalman pass.
1066  unsigned int jhit=0;
1067  for (unsigned int ihit=totHits-2*totHits/(2*fNumIt); ihit<(totHits-totHits/(2*fNumIt)); ihit++) // was ihit<ihit<(totHits-fptsNo)<7
1068  {
1069  feth[jhit] = (Float_t ) (hitMeasCov.at(ihit)[0][0]); // eth
1070  fedudw[jhit] = (Float_t ) (hitMeasCov.at(ihit)[1][1]);
1071  fedvdw[jhit] = (Float_t ) (hitMeasCov.at(ihit)[2][2]);
1072  feu[jhit] = (Float_t ) (hitMeasCov.at(ihit)[3][3]);
1073  fev[jhit] = (Float_t ) (hitMeasCov.at(ihit)[4][4]);
1074  fupdate[jhit] = (Float_t ) (hitUpdate.at(ihit)[0][0]);
1075  fchi2hit[jhit] = (Float_t ) (hitChi2.at(ihit));
1076  jhit++;
1077  }
1078 
1079  stREC->ResizeTo(rep->getState());
1080  *stREC = rep->getState();
1081  covREC->ResizeTo(rep->getCov());
1082  *covREC = rep->getCov();
1083  double dum[5];
1084  double dum2[5];
1085  for (unsigned int ii=0;ii<5;ii++)
1086  {
1087  stREC->ExtractRow(ii,0,dum);
1088  fState0[ii] = dum[0];
1089  covREC->ExtractRow(ii,0,dum2);
1090  for (unsigned int jj=0;jj<5;jj++)
1091  {
1092  fCov0[ii*5+jj] = dum2[jj];
1093  }
1094  }
1095  LOG_DEBUG("Track3DKalmanSPS_GenFit")
1096  << " First State and Cov:" << genf::ROOTobjectToString(*stREC)
1098  chi2 = (Float_t)(rep->getChiSqu());
1099  ndf = rep->getNDF();
1100  nfail = fitTrack.getFailedHits();
1101  nchi2rePass = (int)rePass;
1102  ispptvec=1+std::distance(spptB, sppt);
1103  chi2ndf = (Float_t)(chi2/ndf);
1104 
1105  nTrks++;
1106  LOG_DEBUG("Track3DKalmanSPS_GenFit") << "Track3DKalmanSPS about to do tree->Fill(). Chi2/ndf is " << chi2/ndf << ".";
1107  fpMCMom[3] = MCMomentum.Mag();
1108  for (int ii=0;ii<3;++ii)
1109  {
1110  fpMCMom[ii] = MCMomentum[ii];
1111  fpMCPos[ii] = MCOrigin[ii];
1112  fpREC[ii] = hitPlaneUxUyUz.at(totHits-2*totHits/(2*fNumIt))[ii];
1113  fpRECL[ii] = hitPlaneUxUyUz.at(totHits-totHits/(2*fNumIt)-1)[ii];
1114  }
1115 
1116  evtt = (unsigned int) evt.id().event();
1117  nspptvec = (unsigned int) spptListHandle->size();
1118 
1119  cntp++;
1120  std::vector < std::vector <double> > dQdx;
1121  // Calculate LastFwdPass quantities.
1122  std::vector < TMatrixT<double> > hitCovLFP;
1123  std::vector <TVector3> hitPlaneXYZLFP;
1124  std::vector <TVector3> hitPlaneUxUyUzLFP;
1125  std::vector <TVector3> hitPlanePxPyPzLFP;
1126  std::vector <TVector3> hitPlaneULFP;
1127  std::vector <TVector3> hitPlaneVLFP;
1128  std::vector <double> pLFP;
1129  std::vector < TMatrixT<double> > c7x7LFP;
1130 
1131  art::FindManyP<recob::Hit> hitAssns(spacepointss, evt, fSpptModuleLabel);
1132  for (unsigned int ii=0; ii<totHits/(2*fNumIt); ii++)
1133  {
1134  pLFP.push_back(1./hitState.at(totHits-2*totHits/(2*fNumIt)+ii)[0][0]);
1135  // hitCov -> hitCov7x7 !! EC, 11-May-2012.
1136  c7x7LFP.push_back(hitCov7x7.at(totHits-2*totHits/(2*fNumIt)+ii));
1137  hitCovLFP.push_back(hitCov.at(totHits-2*totHits/(2*fNumIt)+ii));
1138  hitPlaneXYZLFP.push_back(hitPlaneXYZ.at(totHits-2*totHits/(2*fNumIt)+ii));
1139  hitPlaneUxUyUzLFP.push_back(hitPlaneUxUyUz.at(totHits-2*totHits/(2*fNumIt)+ii));
1140  hitPlanePxPyPzLFP.push_back(hitPlaneUxUyUzLFP.back()*pLFP.back());
1141  hitPlaneULFP.push_back(hitPlaneU.at(totHits-2*totHits/(2*fNumIt)+ii));
1142  hitPlaneVLFP.push_back(hitPlaneV.at(totHits-2*totHits/(2*fNumIt)+ii));
1143  // Transform cov appropriate for track rotated
1144  // about w, forcing
1145  // v to be in y-z plane and u pointing in
1146  // +-ive x direction, per TrackAna convention.
1147 
1148  rotationCov(hitCovLFP.back(),
1149  hitPlaneULFP.back(),
1150  hitPlaneVLFP.back()
1151  );
1152  dQdx.push_back(dQdxCalc(hitAssns,
1153  spacepointss,
1154  hitPlaneUxUyUzLFP.back(),
1155  hitPlaneXYZLFP.back()
1156  )
1157  );
1158  fdQdx[ii] = dQdx.back().back();
1159 
1160  }
1161  fpREC[3] = rep->getMom(rep->getReferencePlane()).Mag();
1162  fpRECL[3] = pLFP[1];
1163 
1164  tree->Fill();
1165 
1166  // Put newest track on stack for this set of sppts,
1167  // remove previous one.
1168  recob::tracking::SMatrixSym55 covVtx, covEnd;
1169  for (unsigned int i=0; i<5; i++) {
1170  for (unsigned int j=i; j<5; j++) {
1171  covVtx(i,j) = hitCovLFP.front()(i,j);
1172  covEnd(i,j) = hitCovLFP.back()(i,j);
1173  }
1174  }
1176  recob::tracking::convertCollToVector(hitPlanePxPyPzLFP),
1177  recob::Track::Flags_t(hitPlaneXYZLFP.size()), true),
1178  0, -1., 0, covVtx, covEnd, tcnt++);
1179  if (rePass==1) tcnt1++; // won't get here if Trackfit failed.
1180  if (rePass!=1 && tcnt1) tcol->pop_back();
1181  tcol->push_back(the3DTrack);
1182  util::CreateAssn(*this, evt, *tcol, spacepointss, *tspassn);
1183  art::PtrVector<recob::Hit> hits;// = hitAssns;
1184  for (unsigned int ii=0; ii < spacepointss.size(); ++ii)
1185  {
1186  // for (unsigned int jj=0; jj < hitAssns.at(ii).size(); ++jj)
1187  // hits.push_back(hitAssns.at(ii).at(jj));
1188  hits.insert(hits.end(),hitAssns.at(ii).begin(),hitAssns.at(ii).end());
1189  }
1190  util::CreateAssn(*this, evt, *tcol, hits, *thassn, tcol->size()-1);
1191  } // end !skipFill
1192  } // getStatusFlag
1193 
1194 
1195  if (!rep) delete rep;
1196  rePass++;
1197  // need to first excise bad spacepoints.
1198  // Grab up large Chi2hits first.
1199  art::PtrVector<recob::SpacePoint> spacepointssExcise;
1200  for (unsigned int ind=0;ind<spptSurvivedIndex.size();++ind)
1201  {
1202  // Stricter to chuck sppts from uncontained then contained trks.
1203  if ((uncontained&&fchi2hit[ind] >fChi2Thresh) ||
1204  (!uncontained&&fchi2hit[ind]>1.e9) ||
1205  fchi2hit[ind]<0.0
1206  // =0 eliminates ruled-out large updates. Not obviously
1207  // helpful.
1208  // add a restriction on dQdx here ...
1209  )
1210  {
1211  art::PtrVector<recob::SpacePoint>::iterator spptIt = spacepointss.begin()+spptSurvivedIndex[ind];
1212  spacepointssExcise.push_back(*spptIt);
1213  }
1214  }
1215  // Now grab up those sppts which we skipped and don't want
1216  // to reconsider cuz they're too close to each other or
1217  // cuz they're too far in x, e.g.
1218  for (unsigned int ind=0;ind<spptSkippedIndex.size();++ind)
1219  {
1220  art::PtrVector<recob::SpacePoint>::iterator spptIt = spacepointss.begin()+spptSkippedIndex[ind];
1221  spacepointssExcise.push_back(*spptIt);
1222 
1223  }
1224  // Get rid of redundantly Excised sppts before proceeding.
1225  std::stable_sort(spacepointss.begin(),spacepointss.end());
1226  std::stable_sort(spacepointssExcise.begin(),spacepointssExcise.end());
1227  std::set_union(spacepointssExcise.begin(),spacepointssExcise.end(),
1228  spacepointssExcise.begin(),spacepointssExcise.end(),
1229  spacepointssExcise.begin()
1230  );
1231  // Now excise. New spacepointss will be smaller for second pass.
1233  std::set_difference(spacepointss.begin(),spacepointss.end(),
1234  spacepointssExcise.begin(),spacepointssExcise.end(),
1235  spacepointss.begin()
1236  );
1237  spacepointss.erase(diffSpptIt,spacepointss.end());
1238 
1239  // calculate new seed momentum, and errors as merited
1240  if (rePass==2/* && uncontained */)
1241  {
1242  if (fpREC[3]<fMomHigh && fpREC[3]>fMomLow)
1243  {
1244  double kick(0.9); //Try to get away with a smaller start
1245  // for contained tracks. While for uncontained tracks
1246  // let's start up at a higher momentum and come down.
1247  // Though, 2 (1) GeV/c tracks are too low (high), so
1248  // instead let's actually lower starting value on
1249  // this second pass. -- EC 7-Mar-2013
1250  if (uncontained) kick = 0.5;
1251  for (int ii=0;ii<3;++ii)
1252  {
1253  //mom[ii] = fpREC[ii]*fpREC[3]*kick;
1254  mom[ii] = momM[ii]*kick;
1255  }
1256  }
1257  else if (uncontained)
1258  {
1259  double unstick(1.0);
1260  if (fpREC[3]>=fMomHigh) unstick = 0.3;
1261  for (int ii=0;ii<3;++ii)
1262  {
1263  mom[ii] = momM[ii]*unstick;
1264  }
1265  }
1266  else
1267  for (int ii=0;ii<3;++ii)
1268  {
1269  mom[ii] = 1.1*momM[ii];
1270  }
1271 
1272  }
1273 
1274  } // end while rePass<=maxPass
1275 
1276  sppt++;
1277 
1278  } // end while on elements of std::vector of art::PtrVectors of SpPts.
1279 
1280  if (!repMC) delete repMC;
1281 
1282  evt.put(std::move(tcol));
1283  // and now the spacepoints
1284  evt.put(std::move(tspassn));
1285  // and the hits. Note that these are all the hits from all the spacepoints considered,
1286  // even though they're not all contributing to the tracks.
1287  evt.put(std::move(thassn));
1288 }
TMatrixT< Double_t > * stREC
unsigned int getNDF() const
std::vector< double > fMomErr
void setMomLow(Double_t f)
Definition: GFKalman.h:116
void setMaxUpdate(Double_t f)
Definition: GFKalman.h:118
TMatrixT< Double_t > * stMCT
static bool sp_sort_3dz(const art::Ptr< recob::SpacePoint > &h1, const art::Ptr< recob::SpacePoint > &h2)
geo::Length_t DetHalfWidth(geo::TPCID const &tpcid) const
Returns the half width of the active volume of the specified TPC.
const GFDetPlane & getReferencePlane() const
GFDetPlane getLastPlane() const
void rotationCov(TMatrixT< Double_t > &cov, const TVector3 &u, const TVector3 &v)
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 > > SMatrixSym55
Definition: TrackingTypes.h:85
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
Definition: TrackingTypes.h:58
iterator begin()
Definition: PtrVector.h:223
void setMomHigh(Double_t f)
Definition: GFKalman.h:117
TrackTrajectory::Flags_t Flags_t
Definition: Track.h:71
Float_t tmp
Definition: plot.C:37
std::vector< double > fPosErr
void Print(std::ostream &out=std::cout) const
Definition: GFDetPlane.cxx:242
const TMatrixT< Double_t > & getState() const
bool isRealData() const
Definition: Event.h:83
static bool sp_sort_3dy(const art::Ptr< recob::SpacePoint > &h1, const art::Ptr< recob::SpacePoint > &h2)
#define LOG_ERROR(category)
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
void hits()
Definition: readHits.C:15
geo::Length_t DetHalfHeight(geo::TPCID const &tpcid) const
Returns the half height of the active volume of the specified TPC.
std::vector< double > fMomStart
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
virtual TVector3 getMom(const GFDetPlane &pl)=0
static bool sp_sort_3dx(const art::Ptr< recob::SpacePoint > &h1, const art::Ptr< recob::SpacePoint > &h2)
A trajectory in space reconstructed from hits.
void setBlowUpFactor(double f)
Set the blowup factor (see blowUpCovs() )
Definition: GFKalman.h:115
static bool sp_sort_nsppts(const art::PtrVector< recob::SpacePoint > &h1, const art::PtrVector< recob::SpacePoint > &h2)
TString part[npart]
Definition: Style.C:32
iterator end()
Definition: PtrVector.h:237
double energyLossBetheBloch(const double &mass, const double p)
void setNumIterations(Int_t i)
Set number of iterations for Kalman Filter.
Definition: GFKalman.h:90
TMatrixT< Double_t > * covREC
void processTrack(GFTrack *)
Performs fit on a GFTrack.
Definition: GFKalman.cxx:48
geo::Length_t DetLength(geo::TPCID const &tpcid) const
Returns the length of the active volume of the specified TPC.
static GFFieldManager * getInstance()
TMatrixT< Double_t > * covMCT
bool CreateAssn(PRODUCER const &prod, art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t indx=UINT_MAX)
Creates a single one-to-one association.
bool isDebugEnabled()
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
Definition: TrackingTypes.h:55
const TMatrixT< Double_t > & getCov() const
data_t::iterator iterator
Definition: PtrVector.h:60
std::string ROOTobjectToString(const ROOTOBJ &obj)
Shortcut to write one ROOT object into a string.
Definition: GFException.h:118
size_type size() const
Definition: PtrVector.h:308
iterator insert(iterator position, Ptr< U > const &p)
std::vector< double > dQdxCalc(const art::FindManyP< recob::Hit > &h, const art::PtrVector< recob::SpacePoint > &s, const TVector3 &p, const TVector3 &d)
in close()
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void setInitialDirection(int d)
Sets the inital direction of the track fit (1 for inner to outer, or -1 for outer to inner)...
Definition: GFKalman.h:111
void init(GFAbsBField *b)
set the magntic field here. Magnetic field classes must be derived from GFAbsBField ...
#define LOG_DEBUG(id)
EventNumber_t event() const
Definition: EventID.h:117
void setErrorScaleSTh(Double_t f)
Definition: GFKalman.h:119
double getChiSqu() const
Float_t e
Definition: plot.C:34
RunNumber_t run() const
Definition: Event.h:77
Float_t w
Definition: plot.C:23
EventID id() const
Definition: Event.h:56
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
Definition: Track.h:52
void setErrorScaleMTh(Double_t f)
Definition: GFKalman.h:120
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::string GetLArTPCVolumeName(geo::TPCID const &tpcid) const
Return the name of specified LAr TPC volume.
void trkf::Track3DKalmanSPS::reconfigure ( fhicl::ParameterSet const &  p)

Definition at line 235 of file Track3DKalmanSPS_module.cc.

References fChi2Thresh, fClusterModuleLabel, fDecimate, fDecimateU, fDistanceU, fDoFit, fErrScaleM, fErrScaleS, fG4ModuleLabel, fGenieGenModuleLabel, fMaxPass, fMaxUpdate, fMaxUpdateU, fMinNumSppts, fMomErr, fMomHigh, fMomLow, fMomStart, fNumIt, fPdg, fPerpLim, fPosErr, fSortDim, fSpptModuleLabel, fhicl::ParameterSet::get(), fhicl::ParameterSet::get_if_present(), and LOG_WARNING.

Referenced by Track3DKalmanSPS().

236  {
237 
238  fClusterModuleLabel = pset.get< std::string >("ClusterModuleLabel");
239  fSpptModuleLabel = pset.get< std::string >("SpptModuleLabel");
240  fGenieGenModuleLabel = pset.get< std::string >("GenieGenModuleLabel");
241  fG4ModuleLabel = pset.get< std::string >("G4ModuleLabel");
242  fPosErr = pset.get< std::vector < double > >("PosErr3"); // resolution. cm
243  fMomErr = pset.get< std::vector < double > >("MomErr3"); // GeV
244  fMomStart = pset.get< std::vector < double > >("MomStart3"); //
245  fPerpLim = pset.get< double >("PerpLimit", 1.e6); // PCA cut.
246  fDoFit = pset.get< bool >("DoFit", true); // Der.
247  fNumIt = pset.get< int >("NumIt", 5); // Number x2 passes per fit.
248  fMinNumSppts = pset.get< int >("MinNumSppts", 5); // Min number of sppts in vector to bother fitting
249  fErrScaleS = pset.get< double >("ErrScaleSim", 1.0); // error scale.
250  fErrScaleM = pset.get< double >("ErrScaleMeas", 1.0); // error scale.
251  fDecimate = pset.get< int >("DecimateC", 40); // Sparsify data.
252  fMaxUpdate = pset.get< double >("MaxUpdateC", 0.1); // 0-out.
253  fDecimateU = pset.get< int >("DecimateU", 100);// Sparsify data.
254  fDistanceU = pset.get< double >("DistanceU", 10.0);// Require this separation on uncontained 2nd pass.
255  fMaxUpdateU = pset.get< double >("MaxUpdateU", 0.02); // 0-out.
256  fMomLow = pset.get< double >("MomLow", 0.01); // Fit Range.
257  fMomHigh = pset.get< double >("MomHigh", 20.); // Fit Range.
258  fPdg = pset.get< int >("PdgCode", -13); // mu+ Hypothesis.
259  fChi2Thresh = pset.get< double >("Chi2HitThresh", 12.0E12); //For Re-pass.
260  fSortDim = pset.get< std::string> ("SortDirection", "z"); // case sensitive
261  fMaxPass = pset.get< int >("MaxPass", 2); // mu+ Hypothesis.
262  bool fGenfPRINT;
263  if (pset.get_if_present("GenfPRINT", fGenfPRINT)) {
264  LOG_WARNING("Track3DKalmanSPS_GenFit")
265  << "Parameter 'GenfPRINT' has been deprecated.\n"
266  "Please use the standard message facility to enable GenFit debug output.";
267  // A way to enable debug output is all of the following:
268  // - compile in debug mode (no optimization, no profiling)
269  // - if that makes everything too noisy, add to have everything else quiet
270  // services.message.debugModules: [ "Track3DKalmanSPS" ]
271  // - to print all the GenFit debug messages, set
272  // services.message.destinations.LogDebugFile.categories.Track3DKalmanSPS_GenFit.limit: -1
273  // (assuming there is a LogDebugFile destination already; for example
274  // see the settings in uboonecode/uboone/Utilities/services_microboone.fcl )
275  }
276  }
std::vector< double > fMomErr
std::vector< double > fPosErr
std::vector< double > fMomStart
#define LOG_WARNING(category)
void trkf::Track3DKalmanSPS::rotationCov ( TMatrixT< Double_t > &  cov,
const TVector3 &  u,
const TVector3 &  v 
)
private

Definition at line 317 of file Track3DKalmanSPS_module.cc.

References s, and w.

Referenced by produce().

318  {
319  TVector3 xhat(1.0,0.0,0.0);
320  TVector3 yhat(0.0,1.0,0.0);
321  TVector3 zhat(0.0,0.0,1.0);
322  TVector3 w(u.Cross(v));
323  TVector3 uprime(u);
324  TVector3 vprime(w.Cross(xhat)); // vprime now lies in yz plane
325  Double_t angle(v.Angle(vprime));/* This is the angle through which v
326  must rotate. */
327  uprime.Rotate(angle,w);// u now is rotated the same amount
328  if (uprime*xhat<0)
329  {
330  uprime.Rotate(TMath::Pi(),w);
331  vprime.Rotate(TMath::Pi(),w);
332  angle+=TMath::Pi();
333  }
334  // Build the block-diagonal 5x5 matrix
335  double c = TMath::Cos(angle), s = TMath::Sin(angle);
336  TMatrixT<Double_t> rot(5,5);
337  rot[0][0] = 1.0;
338  rot[1][1] = c;
339  rot[1][2] = s;
340  rot[2][1] = -s;
341  rot[2][2] = c;
342  rot[3][3] = c;
343  rot[3][4] = s;
344  rot[4][3] = -s;
345  rot[4][4] = c;
346 
347  cov=rot*cov;
348  }
Float_t s
Definition: plot.C:23
Float_t w
Definition: plot.C:23
void art::Consumer::showMissingConsumes ( ) const
protectedinherited

Definition at line 125 of file Consumer.cc.

Referenced by art::EDProducer::doEndJob(), art::EDFilter::doEndJob(), art::EDAnalyzer::doEndJob(), and art::RootOutput::endJob().

126 {
127  if (!moduleContext_)
128  return;
129 
130  // If none of the branches have missing consumes statements, exit early.
131  if (std::all_of(cbegin(missingConsumes_),
132  cend(missingConsumes_),
133  [](auto const& perBranch) { return perBranch.empty(); }))
134  return;
135 
136  constexpr cet::HorizontalRule rule{60};
137  mf::LogPrint log{"MTdiagnostics"};
138  log << '\n'
139  << rule('=') << '\n'
140  << "The following consumes (or mayConsume) statements are missing from\n"
141  << module_context(moduleDescription_) << '\n'
142  << rule('-') << '\n';
143 
144  cet::for_all_with_index(
145  missingConsumes_, [&log](std::size_t const i, auto const& perBranch) {
146  for (auto const& pi : perBranch) {
147  log << " "
148  << assemble_consumes_statement(static_cast<BranchType>(i), pi)
149  << '\n';
150  }
151  });
152  log << rule('=');
153 }
cet::exempt_ptr< ModuleDescription const > moduleDescription_
Definition: Consumer.h:140
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
bool moduleContext_
Definition: Consumer.h:136
ConsumableProductSets missingConsumes_
Definition: Consumer.h:139
void art::Consumer::validateConsumedProduct ( BranchType const  bt,
ProductInfo const &  pi 
)
protectedinherited

Definition at line 101 of file Consumer.cc.

References art::errors::ProductRegistrationFailure.

103 {
104  // Early exits if consumes tracking has been disabled or if the
105  // consumed product is an allowed consumable.
106  if (!moduleContext_)
107  return;
108 
109  if (cet::binary_search_all(consumables_[bt], pi))
110  return;
111 
112  if (requireConsumes_) {
114  "Consumer: an error occurred during validation of a "
115  "retrieved product\n\n")
116  << "The following consumes (or mayConsume) statement is missing from\n"
117  << module_context(moduleDescription_) << ":\n\n"
118  << " " << assemble_consumes_statement(bt, pi) << "\n\n";
119  }
120 
121  missingConsumes_[bt].insert(pi);
122 }
cet::exempt_ptr< ModuleDescription const > moduleDescription_
Definition: Consumer.h:140
bool requireConsumes_
Definition: Consumer.h:137
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
constexpr T pi()
Returns the constant pi (up to 35 decimal digits of precision)
ConsumableProducts consumables_
Definition: Consumer.h:138
bool moduleContext_
Definition: Consumer.h:136
ConsumableProductSets missingConsumes_
Definition: Consumer.h:139

Member Data Documentation

Float_t trkf::Track3DKalmanSPS::chi2
private

Definition at line 131 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), and produce().

Float_t trkf::Track3DKalmanSPS::chi2ndf
private

Definition at line 132 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), and produce().

TMatrixT<Double_t>* trkf::Track3DKalmanSPS::covMCT
private

Definition at line 128 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), and produce().

TMatrixT<Double_t>* trkf::Track3DKalmanSPS::covREC
private

Definition at line 130 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), and produce().

unsigned int trkf::Track3DKalmanSPS::evtt
private

Definition at line 146 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::fchi2hit
private

Definition at line 157 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

double trkf::Track3DKalmanSPS::fChi2Thresh
private

Definition at line 191 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

std::string trkf::Track3DKalmanSPS::fClusterModuleLabel
private

Definition at line 118 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

int trkf::Track3DKalmanSPS::fcont
private

Definition at line 133 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::fCov0
private

Definition at line 140 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

int trkf::Track3DKalmanSPS::fDecimate
private

Definition at line 183 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

int trkf::Track3DKalmanSPS::fDecimateU
private

Definition at line 185 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

unsigned int trkf::Track3DKalmanSPS::fDimSize
private

Definition at line 166 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), and produce().

double trkf::Track3DKalmanSPS::fDistanceU
private

Definition at line 186 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

bool trkf::Track3DKalmanSPS::fDoFit
private

Definition at line 178 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

Float_t* trkf::Track3DKalmanSPS::fdQdx
private

Definition at line 165 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::fedudw
private

Definition at line 160 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::fedvdw
private

Definition at line 161 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

double trkf::Track3DKalmanSPS::fErrScaleM
private

Definition at line 182 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

double trkf::Track3DKalmanSPS::fErrScaleS
private

Definition at line 181 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

Float_t* trkf::Track3DKalmanSPS::feshx
private

Definition at line 152 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::feshy
private

Definition at line 153 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::feshyz
private

Definition at line 155 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::feshz
private

Definition at line 154 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::feth
private

Definition at line 159 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::feu
private

Definition at line 162 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::fev
private

Definition at line 163 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

std::string trkf::Track3DKalmanSPS::fG4ModuleLabel
private

Definition at line 121 of file Track3DKalmanSPS_module.cc.

Referenced by reconfigure().

std::string trkf::Track3DKalmanSPS::fGenieGenModuleLabel
private

Definition at line 120 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

int trkf::Track3DKalmanSPS::fMaxPass
private

Definition at line 192 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

double trkf::Track3DKalmanSPS::fMaxUpdate
private

Definition at line 184 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

double trkf::Track3DKalmanSPS::fMaxUpdateU
private

Definition at line 187 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

uint16_t trkf::Track3DKalmanSPS::fMinNumSppts
private

Definition at line 180 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

std::vector<double> trkf::Track3DKalmanSPS::fMomErr
private

Definition at line 175 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

double trkf::Track3DKalmanSPS::fMomHigh
private

Definition at line 189 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

double trkf::Track3DKalmanSPS::fMomLow
private

Definition at line 188 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

std::vector<double> trkf::Track3DKalmanSPS::fMomStart
private

Definition at line 176 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

int trkf::Track3DKalmanSPS::fNumIt
private

Definition at line 179 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

Float_t* trkf::Track3DKalmanSPS::fPC1
private

Definition at line 170 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::fPC2
private

Definition at line 171 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::fPC3
private

Definition at line 172 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::fPCevals
private

Definition at line 168 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::fPCmeans
private

Definition at line 167 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::fPCsigmas
private

Definition at line 169 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

int trkf::Track3DKalmanSPS::fPdg
private

Definition at line 190 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

double trkf::Track3DKalmanSPS::fPerpLim
private

Definition at line 177 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

Float_t* trkf::Track3DKalmanSPS::fpMCMom
private

Definition at line 137 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::fpMCPos
private

Definition at line 138 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), and produce().

std::vector<double> trkf::Track3DKalmanSPS::fPosErr
private

Definition at line 174 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

Float_t* trkf::Track3DKalmanSPS::fpREC
private

Definition at line 136 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::fpRECL
private

Definition at line 135 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

unsigned int trkf::Track3DKalmanSPS::fptsNo
private

Definition at line 148 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::fsep
private

Definition at line 164 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::fshx
private

Definition at line 149 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::fshy
private

Definition at line 150 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::fshz
private

Definition at line 151 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

std::string trkf::Track3DKalmanSPS::fSortDim
private

Definition at line 122 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

std::string trkf::Track3DKalmanSPS::fSpptModuleLabel
private

Definition at line 119 of file Track3DKalmanSPS_module.cc.

Referenced by produce(), and reconfigure().

Float_t* trkf::Track3DKalmanSPS::fState0
private

Definition at line 139 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::fth
private

Definition at line 158 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

Float_t* trkf::Track3DKalmanSPS::fupdate
private

Definition at line 156 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), endJob(), and produce().

int trkf::Track3DKalmanSPS::ispptvec
private

Definition at line 144 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), and produce().

int trkf::Track3DKalmanSPS::nchi2rePass
private

Definition at line 143 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), and produce().

int trkf::Track3DKalmanSPS::ndf
private

Definition at line 142 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), and produce().

int trkf::Track3DKalmanSPS::nfail
private

Definition at line 141 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), and produce().

int trkf::Track3DKalmanSPS::nspptvec
private

Definition at line 145 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), and produce().

unsigned int trkf::Track3DKalmanSPS::nTrks
private

Definition at line 147 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), and produce().

genf::GFAbsTrackRep* trkf::Track3DKalmanSPS::rep
private

Definition at line 195 of file Track3DKalmanSPS_module.cc.

Referenced by endJob(), and produce().

genf::GFAbsTrackRep* trkf::Track3DKalmanSPS::repMC
private

Definition at line 194 of file Track3DKalmanSPS_module.cc.

Referenced by endJob(), and produce().

TMatrixT<Double_t>* trkf::Track3DKalmanSPS::stMCT
private

Definition at line 127 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), and produce().

TMatrixT<Double_t>* trkf::Track3DKalmanSPS::stREC
private

Definition at line 129 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), and produce().

TTree* trkf::Track3DKalmanSPS::tree
private

Definition at line 125 of file Track3DKalmanSPS_module.cc.

Referenced by beginJob(), and produce().


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