Public Member Functions | Private Member Functions | Private Attributes | List of all members
trkf::Track3DKalmanSPS Class Reference
Inheritance diagram for trkf::Track3DKalmanSPS:
art::EDProducer art::detail::Producer art::detail::LegacyModule art::Modifier art::ModuleBase art::ProductRegistryHelper

Public Member Functions

 Track3DKalmanSPS (fhicl::ParameterSet const &pset)
 
- Public Member Functions inherited from art::EDProducer
 EDProducer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 EDProducer (Table< Config > const &config)
 
std::string workerType () const
 
- Public Member Functions inherited from art::detail::Producer
virtual ~Producer () noexcept
 
 Producer (fhicl::ParameterSet const &)
 
 Producer (Producer const &)=delete
 
 Producer (Producer &&)=delete
 
Producer & operator= (Producer const &)=delete
 
Producer & operator= (Producer &&)=delete
 
void doBeginJob (SharedResources const &resources)
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
- Public Member Functions inherited from art::Modifier
 ~Modifier () noexcept
 
 Modifier ()
 
 Modifier (Modifier const &)=delete
 
 Modifier (Modifier &&)=delete
 
Modifier & operator= (Modifier const &)=delete
 
Modifier & operator= (Modifier &&)=delete
 
- Public Member Functions inherited from art::ModuleBase
virtual ~ModuleBase () noexcept
 
 ModuleBase ()
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 

Private Member Functions

void produce (art::Event &evt)
 
void beginJob ()
 
void endJob ()
 
double energyLossBetheBloch (const double &mass, const double p)
 
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::GFAbsTrackRep * repMC
 
genf::GFAbsTrackRep * rep
 

Additional Inherited Members

- Public Types inherited from art::EDProducer
using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
- Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using Table = Modifier::Table< UserConfig, KeysToIgnore >
 
- Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >
 
- Static Public Member Functions inherited from art::EDProducer
static void commitEvent (EventPrincipal &ep, Event &e)
 
- Protected Member Functions inherited from art::ModuleBase
ConsumesCollector & consumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 

Detailed Description

Definition at line 92 of file Track3DKalmanSPS_module.cc.

Constructor & Destructor Documentation

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

Definition at line 193 of file Track3DKalmanSPS_module.cc.

194  : EDProducer{pset}
195  , fDoFit(true)
196  , fNumIt(5)
197  , fMinNumSppts(5)
198  , fErrScaleS(1.0)
199  , fErrScaleM(1.0)
200  , fDecimate(1)
201  , fMaxUpdate(0.10)
202  , fDecimateU(1)
203  , fDistanceU(1.)
204  , fMaxUpdateU(0.10)
205  , fMomLow(0.001)
206  , fMomHigh(100.)
207  , fPdg(-13)
208  , fChi2Thresh(12.0E12)
209  , fMaxPass (1)
210  {
211  fClusterModuleLabel = pset.get< std::string >("ClusterModuleLabel");
212  fSpptModuleLabel = pset.get< std::string >("SpptModuleLabel");
213  fGenieGenModuleLabel = pset.get< std::string >("GenieGenModuleLabel");
214  fG4ModuleLabel = pset.get< std::string >("G4ModuleLabel");
215  fPosErr = pset.get< std::vector < double > >("PosErr3"); // resolution. cm
216  fMomErr = pset.get< std::vector < double > >("MomErr3"); // GeV
217  fMomStart = pset.get< std::vector < double > >("MomStart3"); //
218  fPerpLim = pset.get< double >("PerpLimit", 1.e6); // PCA cut.
219  fDoFit = pset.get< bool >("DoFit", true); // Der.
220  fNumIt = pset.get< int >("NumIt", 5); // Number x2 passes per fit.
221  fMinNumSppts = pset.get< int >("MinNumSppts", 5); // Min number of sppts in vector to bother fitting
222  fErrScaleS = pset.get< double >("ErrScaleSim", 1.0); // error scale.
223  fErrScaleM = pset.get< double >("ErrScaleMeas", 1.0); // error scale.
224  fDecimate = pset.get< int >("DecimateC", 40); // Sparsify data.
225  fMaxUpdate = pset.get< double >("MaxUpdateC", 0.1); // 0-out.
226  fDecimateU = pset.get< int >("DecimateU", 100);// Sparsify data.
227  fDistanceU = pset.get< double >("DistanceU", 10.0);// Require this separation on uncontained 2nd pass.
228  fMaxUpdateU = pset.get< double >("MaxUpdateU", 0.02); // 0-out.
229  fMomLow = pset.get< double >("MomLow", 0.01); // Fit Range.
230  fMomHigh = pset.get< double >("MomHigh", 20.); // Fit Range.
231  fPdg = pset.get< int >("PdgCode", -13); // mu+ Hypothesis.
232  fChi2Thresh = pset.get< double >("Chi2HitThresh", 12.0E12); //For Re-pass.
233  fSortDim = pset.get< std::string> ("SortDirection", "z"); // case sensitive
234  fMaxPass = pset.get< int >("MaxPass", 2); // mu+ Hypothesis.
235  bool fGenfPRINT;
236  if (pset.get_if_present("GenfPRINT", fGenfPRINT)) {
237  MF_LOG_WARNING("Track3DKalmanSPS_GenFit")
238  << "Parameter 'GenfPRINT' has been deprecated.\n"
239  "Please use the standard message facility to enable GenFit debug output.";
240  // A way to enable debug output is all of the following:
241  // - compile in debug mode (no optimization, no profiling)
242  // - if that makes everything too noisy, add to have everything else quiet
243  // services.message.debugModules: [ "Track3DKalmanSPS" ]
244  // - to print all the GenFit debug messages, set
245  // services.message.destinations.LogDebugFile.categories.Track3DKalmanSPS_GenFit.limit: -1
246  // (assuming there is a LogDebugFile destination already; for example
247  // see the settings in uboonecode/uboone/Utilities/services_microboone.fcl )
248  }
249 
250  produces< std::vector<recob::Track> >();
251  produces<art::Assns<recob::Track, recob::Cluster> >();
252  produces<art::Assns<recob::Track, recob::SpacePoint> >();
253  produces<art::Assns<recob::Track, recob::Hit> >();
254 
255  }
trkf::Track3DKalmanSPS::fMomErr
std::vector< double > fMomErr
Definition: Track3DKalmanSPS_module.cc:164
string
std::string string
Definition: nybbler.cc:12
art::EDProducer::EDProducer
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
trkf::Track3DKalmanSPS::fPosErr
std::vector< double > fPosErr
Definition: Track3DKalmanSPS_module.cc:163
trkf::Track3DKalmanSPS::fMomStart
std::vector< double > fMomStart
Definition: Track3DKalmanSPS_module.cc:165
MF_LOG_WARNING
#define MF_LOG_WARNING(category)
Definition: MessageLogger.h:227

Member Function Documentation

void trkf::Track3DKalmanSPS::beginJob ( )
privatevirtual

Reimplemented from art::EDProducer.

Definition at line 396 of file Track3DKalmanSPS_module.cc.

397  {
398 
399 
400  art::ServiceHandle<art::TFileService const> tfs;
401 
402  stMCT = new TMatrixT<Double_t>(5,1);
403  covMCT = new TMatrixT<Double_t>(5,5);
404  stREC = new TMatrixT<Double_t>(5,1);
405  covREC = new TMatrixT<Double_t>(5,5);
406 
407  fpMCMom = new Float_t[4];
408  fpMCPos = new Float_t[4];
409  fpREC = new Float_t[4];
410  fpRECL = new Float_t[4];
411  fState0 = new Float_t[5];
412  fCov0 = new Float_t[25];
413  fDimSize = 20000; // if necessary will get this from pset in constructor.
414 
415  fshx = new Float_t[fDimSize];
416  fshy = new Float_t[fDimSize];
417  fshz = new Float_t[fDimSize];
418  feshx = new Float_t[fDimSize];
419  feshy = new Float_t[fDimSize];
420  feshz = new Float_t[fDimSize];
421  feshyz = new Float_t[fDimSize];
422  fupdate = new Float_t[fDimSize];
423  fchi2hit = new Float_t[fDimSize];
424  fth = new Float_t[fDimSize];
425  feth = new Float_t[fDimSize];
426  fedudw = new Float_t[fDimSize];
427  fedvdw = new Float_t[fDimSize];
428  feu = new Float_t[fDimSize];
429  fev = new Float_t[fDimSize];
430  fsep = new Float_t[fDimSize];
431  fdQdx = new Float_t[fDimSize];
432 
433  fPC1 = new Float_t[3];
434  fPC2 = new Float_t[3];
435  fPC3 = new Float_t[3];
436  fPCmeans = new Float_t[3];
437  fPCsigmas = new Float_t[3];
438  fPCevals = new Float_t[3];
439 
440 // //TFile fileGENFIT("GENFITout.root","RECREATE");
441 
442 
443  tree = tfs->make<TTree>("GENFITttree","GENFITttree");
444  //tree->Branch("stMCT",&stMCT,"stMCT[5]/F"); // "TMatrixT<Double_t>"
445 
446  //tree->Branch("stMCT","TMatrixD",&stMCT,64000,0);
447  //tree->Branch("covMCT",covMCT,"covMCT[25]/F");
448  tree->Branch("covMCT","TMatrixD",&covMCT,64000,0);
449  tree->Branch("stREC",fState0,"stREC[5]/F");
450  //tree->Branch("stREC","TMatrixD",&stREC,64000,0);
451  tree->Branch("covREC",fCov0,"covREC[25]/F");
452  //tree->Branch("covREC","TMatrixD",&covREC,64000,0);
453 
454  tree->Branch("nchi2rePass",&nchi2rePass,"nchi2rePass/I");
455  tree->Branch("ispptvec",&ispptvec,"ispptvec/I");
456  tree->Branch("chi2",&chi2,"chi2/F");
457  tree->Branch("nfail",&nfail,"nfail/I");
458  tree->Branch("ndf",&ndf,"ndf/I");
459  tree->Branch("evtNo",&evtt,"evtNo/I");
460  tree->Branch("nspptvec",&nspptvec,"nspptvec/I");
461  tree->Branch("chi2ndf",&chi2ndf,"chi2ndf/F");
462 
463  tree->Branch("trkNo",&nTrks,"trkNo/I");
464  tree->Branch("ptsNo",&fptsNo,"ptsNo/I");
465  tree->Branch("cont",&fcont,"cont/I"); //O? Yes, O. Not 0, not L, ...
466  tree->Branch("shx",fshx,"shx[ptsNo]/F");
467  tree->Branch("shy",fshy,"shy[ptsNo]/F");
468  tree->Branch("shz",fshz,"shz[ptsNo]/F");
469  tree->Branch("sep",fsep,"sep[ptsNo]/F");
470  tree->Branch("dQdx",fdQdx,"dQdx[ptsNo]/F");
471  tree->Branch("eshx",feshx,"eshx[ptsNo]/F");
472  tree->Branch("eshy",feshy,"eshy[ptsNo]/F");
473  tree->Branch("eshz",feshz,"eshz[ptsNo]/F");
474  tree->Branch("eshyz",feshyz,"eshyz[ptsNo]/F");
475  tree->Branch("update",fupdate,"update[ptsNo]/F");
476  tree->Branch("chi2hit",fchi2hit,"chi2hit[ptsNo]/F");
477  tree->Branch("th",fth,"th[ptsNo]/F");
478  tree->Branch("eth",feth,"eth[ptsNo]/F");
479  tree->Branch("edudw",fedudw,"edudw[ptsNo]/F");
480  tree->Branch("edvdw",fedvdw,"edvdw[ptsNo]/F");
481  tree->Branch("eu",feu,"eu[ptsNo]/F");
482  tree->Branch("ev",fev,"ev[ptsNo]/F");
483 
484 
485  tree->Branch("pcMeans", fPCmeans,"pcMeans[3]/F");
486  tree->Branch("pcSigmas",fPCsigmas,"pcSigmas[3]/F");
487  tree->Branch("pcEvals", fPCevals,"pcEvals[3]/F");
488  tree->Branch("pcEvec1",fPC1,"pcEvec1[3]/F");
489  tree->Branch("pcEvec2",fPC2,"pcEvec2[3]/F");
490  tree->Branch("pcEvec3",fPC3,"pcEvec3[3]/F");
491 
492 
493  tree->Branch("pMCMom",fpMCMom,"pMCMom[4]/F");
494  tree->Branch("pMCPos",fpMCPos,"pMCPos[4]/F");
495  tree->Branch("pRECKalF",fpREC,"pRECKalF[4]/F");
496  tree->Branch("pRECKalL",fpRECL,"pRECKalL[4]/F");
497 
498 
499  //TGeoManager* geomGENFIT = new TGeoManager("Geometry", "Geane geometry");
500  //TGeoManager::Import("config/genfitGeom.root");
501  // gROOT->Macro("config/Geane.C");
502 
503  }
trkf::Track3DKalmanSPS::stREC
TMatrixT< Double_t > * stREC
Definition: Track3DKalmanSPS_module.cc:118
trkf::Track3DKalmanSPS::stMCT
TMatrixT< Double_t > * stMCT
Definition: Track3DKalmanSPS_module.cc:116
trkf::Track3DKalmanSPS::covREC
TMatrixT< Double_t > * covREC
Definition: Track3DKalmanSPS_module.cc:119
trkf::Track3DKalmanSPS::covMCT
TMatrixT< Double_t > * covMCT
Definition: Track3DKalmanSPS_module.cc:117
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 325 of file Track3DKalmanSPS_module.cc.

326  {
327  // For now just Collection plane.
328  // We should loop over all views, more generally.
329  geo::SigType_t sig(geo::kCollection);
330  art::ServiceHandle<geo::Geometry const> geom;
331  static art::PtrVector<recob::SpacePoint>::const_iterator sstart(s.begin());
332  // art::PtrVector<recob::SpacePoint>::const_iterator sppt = sstart;
333  art::PtrVector<recob::SpacePoint>::const_iterator sppt = s.begin();
334  std::vector <double> v;
335 
336 
337 
338  double mindist(100.0); // cm
339  auto spptminIt(sppt);
340  while (sppt != s.end())
341  {
342  if (((**sppt).XYZ() - loc).Mag() < mindist)
343  {
344  double dist = ((**sppt).XYZ() - loc).Mag();
345 
346  // Jump out if we're as close as 0.1 mm away.
347  if (dist<mindist)
348  {
349  mindist = dist;
350  spptminIt = sppt;
351  if (mindist < 0.01) break;
352  }
353  }
354  sppt++;
355  }
356  sstart = spptminIt; // for next time.
357  unsigned int ind(std::distance(s.begin(),spptminIt));
358 
359 
360 
361  std::vector< art::Ptr<recob::Hit> > hitlist = h.at(ind);
362 
363  double wirePitch = 0.;
364  double angleToVert = 0;
365  // unsigned int tpc1;
366  unsigned int plane1;
367  double charge = 0.;
368 
369  for(std::vector< art::Ptr<recob::Hit> >::const_iterator ihit = hitlist.begin();
370  ihit != hitlist.end(); ++ihit)
371  {
372  const recob::Hit& hit1 = **ihit;
373  // if (hit1.View() != view) continue;
374  if (hit1.SignalType() != sig) continue;
375  geo::WireID hit1WireID = hit1.WireID();
376  // tpc1 = hit1WireID.TPC;
377  plane1 = hit1WireID.Plane;
378  charge = hit1.Integral();
379  wirePitch = geom->WirePitch(plane1);
380  angleToVert = geom->Plane(plane1).Wire(0).ThetaZ(false) - 0.5*TMath::Pi();
381  }
382 
383  double cosgamma = TMath::Abs(TMath::Sin(angleToVert)*dir.Y() +
384  TMath::Cos(angleToVert)*dir.Z());
385  // if(cosgamma < 1.e-5)
386  // throw cet::exception("Track") << "cosgamma is basically 0, that can't be right\n";
387 
388  v.push_back(charge/wirePitch/cosgamma);
389  // std::cout << " Track3DKalmanSPS::dQdxCalc() : For loc.XYZ() hit is ... " << ind << " and v is " << v.back() << std::endl;
390 
391  return v;
392 
393  }
recob::Hit::SignalType
geo::SigType_t SignalType() const
Signal type for the plane of the hit.
Definition: Hit.h:231
geo::GeometryCore::Plane
PlaneGeo const & Plane(unsigned int const p, unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified wire.
Definition: GeometryCore.h:2994
geo::PlaneGeo::Wire
WireGeo const & Wire(unsigned int iwire) const
Definition: PlaneGeo.cxx:506
recob::Hit::WireID
geo::WireID WireID() const
Definition: Hit.h:233
art::PtrVector::begin
iterator begin()
Definition: PtrVector.h:217
vector
struct vector vector
recob::Hit::Integral
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:224
dir
string dir
Definition: INukeNucleonCorr.cxx:61
geo::WireGeo::ThetaZ
double ThetaZ() const
Returns angle of wire with respect to z axis in the Y-Z plane in radians.
Definition: WireGeo.h:250
geo::GeometryCore::WirePitch
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
Definition: GeometryCore.cxx:824
reco_momentum_tuples.h
h
Definition: reco_momentum_tuples.py:66
art::PtrVector::const_iterator
typename data_t::const_iterator const_iterator
Definition: PtrVector.h:55
geo::SigType_t
enum geo::_plane_sigtype SigType_t
art::PtrVector::end
iterator end()
Definition: PtrVector.h:231
protoana::distance
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
Definition: truepionXsection_module.cc:217
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
larg4::dist
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Definition: OpFastScintillation.hh:563
recob::Hit
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
geo::kCollection
Signal from collection planes.
Definition: geo_types.h:146
void trkf::Track3DKalmanSPS::endJob ( )
privatevirtual

Reimplemented from art::EDProducer.

Definition at line 506 of file Track3DKalmanSPS_module.cc.

507  {
508  if (!rep) delete rep;
509  if (!repMC) delete repMC;
510 
511  /*
512  // not sure why I can't do these, but at least some cause seg faults.
513  delete[] stMCT;
514  delete[] covMCT;
515  delete[] stREC;
516  delete[] covREC;
517  */
518 
519  delete[] fpREC;
520  delete[] fpRECL;
521  delete[] fState0;
522  delete[] fCov0;
523 
524  delete[] fshx;
525  delete[] fshy;
526  delete[] fshz;
527  delete[] feshx;
528  delete[] feshy;
529  delete[] feshyz;
530  delete[] feshz;
531  delete[] fupdate;
532  delete[] fchi2hit;
533  delete[] fth;
534  delete[] feth;
535  delete[] fedudw;
536  delete[] fedvdw;
537  delete[] feu;
538  delete[] fev;
539  delete[] fsep;
540  delete[] fdQdx;
541 
542  delete[] fPCmeans;
543  delete[] fPCsigmas;
544  delete[] fPCevals;
545  delete[] fPC1;
546  delete[] fPC2;
547  delete[] fPC3;
548 
549  }
double trkf::Track3DKalmanSPS::energyLossBetheBloch ( const double &  mass,
const double  p = 1.5 
)
private

Definition at line 260 of file Track3DKalmanSPS_module.cc.

263  {
264  const double charge(1.0);
265  const double mEE(188.); // eV
266  const double matZ(18.);
267  const double matA(40.);
268  const double matDensity(1.4);
269  const double me(0.000511);
270 
271  double beta = p/std::sqrt(mass*mass+p*p);
272  double gammaSquare = 1./(1.0 - beta*beta);
273  // 4pi.r_e^2.N.me = 0.307075, I think.
274  double dedx = 0.307075*matDensity*matZ/matA/(beta*beta)*charge*charge;
275  double massRatio = me/mass;
276  // me=0.000511 here is in GeV. So mEE comes in here in eV.
277  double argument = gammaSquare*beta*beta*me*1.E3*2./((1.E-6*mEE) * std::sqrt(1+2*std::sqrt(gammaSquare)*massRatio + massRatio*massRatio));
278 
279  if (mass==0.0) return(0.0);
280  if (argument <= exp(beta*beta))
281  {
282  dedx = 0.;
283  }
284  else{
285  dedx *= (log(argument)-beta*beta); // Bethe-Bloch [MeV/cm]
286  dedx *= 1.E-3; // in GeV/cm, hence 1.e-3
287  if (dedx<0.) dedx = 0.;
288  }
289  return dedx;
290  }
pionana::beta
double beta(double KE, const simb::MCParticle *part)
Definition: TruthAnalyzer_module.cc:59
test.p
p
Definition: test.py:223
void trkf::Track3DKalmanSPS::produce ( art::Event &  evt)
privatevirtual
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 553 of file Track3DKalmanSPS_module.cc.

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

Definition at line 292 of file Track3DKalmanSPS_module.cc.

293  {
294  TVector3 xhat(1.0,0.0,0.0);
295  TVector3 yhat(0.0,1.0,0.0);
296  TVector3 zhat(0.0,0.0,1.0);
297  TVector3 w(u.Cross(v));
298  TVector3 uprime(u);
299  TVector3 vprime(w.Cross(xhat)); // vprime now lies in yz plane
300  Double_t angle(v.Angle(vprime));/* This is the angle through which v
301  must rotate. */
302  uprime.Rotate(angle,w);// u now is rotated the same amount
303  if (uprime*xhat<0)
304  {
305  uprime.Rotate(TMath::Pi(),w);
306  vprime.Rotate(TMath::Pi(),w);
307  angle+=TMath::Pi();
308  }
309  // Build the block-diagonal 5x5 matrix
310  double c = TMath::Cos(angle), s = TMath::Sin(angle);
311  TMatrixT<Double_t> rot(5,5);
312  rot[0][0] = 1.0;
313  rot[1][1] = c;
314  rot[1][2] = s;
315  rot[2][1] = -s;
316  rot[2][2] = c;
317  rot[3][3] = c;
318  rot[3][4] = s;
319  rot[4][3] = -s;
320  rot[4][4] = c;
321 
322  cov=rot*cov;
323  }
s
static QCString * s
Definition: config.cpp:1042

Member Data Documentation

Float_t trkf::Track3DKalmanSPS::chi2
private

Definition at line 120 of file Track3DKalmanSPS_module.cc.

Float_t trkf::Track3DKalmanSPS::chi2ndf
private

Definition at line 121 of file Track3DKalmanSPS_module.cc.

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

Definition at line 117 of file Track3DKalmanSPS_module.cc.

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

Definition at line 119 of file Track3DKalmanSPS_module.cc.

unsigned int trkf::Track3DKalmanSPS::evtt
private

Definition at line 135 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fchi2hit
private

Definition at line 146 of file Track3DKalmanSPS_module.cc.

double trkf::Track3DKalmanSPS::fChi2Thresh
private

Definition at line 180 of file Track3DKalmanSPS_module.cc.

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

Definition at line 107 of file Track3DKalmanSPS_module.cc.

int trkf::Track3DKalmanSPS::fcont
private

Definition at line 122 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fCov0
private

Definition at line 129 of file Track3DKalmanSPS_module.cc.

int trkf::Track3DKalmanSPS::fDecimate
private

Definition at line 172 of file Track3DKalmanSPS_module.cc.

int trkf::Track3DKalmanSPS::fDecimateU
private

Definition at line 174 of file Track3DKalmanSPS_module.cc.

unsigned int trkf::Track3DKalmanSPS::fDimSize
private

Definition at line 155 of file Track3DKalmanSPS_module.cc.

double trkf::Track3DKalmanSPS::fDistanceU
private

Definition at line 175 of file Track3DKalmanSPS_module.cc.

bool trkf::Track3DKalmanSPS::fDoFit
private

Definition at line 167 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fdQdx
private

Definition at line 154 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fedudw
private

Definition at line 149 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fedvdw
private

Definition at line 150 of file Track3DKalmanSPS_module.cc.

double trkf::Track3DKalmanSPS::fErrScaleM
private

Definition at line 171 of file Track3DKalmanSPS_module.cc.

double trkf::Track3DKalmanSPS::fErrScaleS
private

Definition at line 170 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::feshx
private

Definition at line 141 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::feshy
private

Definition at line 142 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::feshyz
private

Definition at line 144 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::feshz
private

Definition at line 143 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::feth
private

Definition at line 148 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::feu
private

Definition at line 151 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fev
private

Definition at line 152 of file Track3DKalmanSPS_module.cc.

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

Definition at line 110 of file Track3DKalmanSPS_module.cc.

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

Definition at line 109 of file Track3DKalmanSPS_module.cc.

int trkf::Track3DKalmanSPS::fMaxPass
private

Definition at line 181 of file Track3DKalmanSPS_module.cc.

double trkf::Track3DKalmanSPS::fMaxUpdate
private

Definition at line 173 of file Track3DKalmanSPS_module.cc.

double trkf::Track3DKalmanSPS::fMaxUpdateU
private

Definition at line 176 of file Track3DKalmanSPS_module.cc.

uint16_t trkf::Track3DKalmanSPS::fMinNumSppts
private

Definition at line 169 of file Track3DKalmanSPS_module.cc.

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

Definition at line 164 of file Track3DKalmanSPS_module.cc.

double trkf::Track3DKalmanSPS::fMomHigh
private

Definition at line 178 of file Track3DKalmanSPS_module.cc.

double trkf::Track3DKalmanSPS::fMomLow
private

Definition at line 177 of file Track3DKalmanSPS_module.cc.

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

Definition at line 165 of file Track3DKalmanSPS_module.cc.

int trkf::Track3DKalmanSPS::fNumIt
private

Definition at line 168 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fPC1
private

Definition at line 159 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fPC2
private

Definition at line 160 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fPC3
private

Definition at line 161 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fPCevals
private

Definition at line 157 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fPCmeans
private

Definition at line 156 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fPCsigmas
private

Definition at line 158 of file Track3DKalmanSPS_module.cc.

int trkf::Track3DKalmanSPS::fPdg
private

Definition at line 179 of file Track3DKalmanSPS_module.cc.

double trkf::Track3DKalmanSPS::fPerpLim
private

Definition at line 166 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fpMCMom
private

Definition at line 126 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fpMCPos
private

Definition at line 127 of file Track3DKalmanSPS_module.cc.

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

Definition at line 163 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fpREC
private

Definition at line 125 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fpRECL
private

Definition at line 124 of file Track3DKalmanSPS_module.cc.

unsigned int trkf::Track3DKalmanSPS::fptsNo
private

Definition at line 137 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fsep
private

Definition at line 153 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fshx
private

Definition at line 138 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fshy
private

Definition at line 139 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fshz
private

Definition at line 140 of file Track3DKalmanSPS_module.cc.

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

Definition at line 111 of file Track3DKalmanSPS_module.cc.

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

Definition at line 108 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fState0
private

Definition at line 128 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fth
private

Definition at line 147 of file Track3DKalmanSPS_module.cc.

Float_t* trkf::Track3DKalmanSPS::fupdate
private

Definition at line 145 of file Track3DKalmanSPS_module.cc.

int trkf::Track3DKalmanSPS::ispptvec
private

Definition at line 133 of file Track3DKalmanSPS_module.cc.

int trkf::Track3DKalmanSPS::nchi2rePass
private

Definition at line 132 of file Track3DKalmanSPS_module.cc.

int trkf::Track3DKalmanSPS::ndf
private

Definition at line 131 of file Track3DKalmanSPS_module.cc.

int trkf::Track3DKalmanSPS::nfail
private

Definition at line 130 of file Track3DKalmanSPS_module.cc.

int trkf::Track3DKalmanSPS::nspptvec
private

Definition at line 134 of file Track3DKalmanSPS_module.cc.

unsigned int trkf::Track3DKalmanSPS::nTrks
private

Definition at line 136 of file Track3DKalmanSPS_module.cc.

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

Definition at line 184 of file Track3DKalmanSPS_module.cc.

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

Definition at line 183 of file Track3DKalmanSPS_module.cc.

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

Definition at line 116 of file Track3DKalmanSPS_module.cc.

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

Definition at line 118 of file Track3DKalmanSPS_module.cc.

TTree* trkf::Track3DKalmanSPS::tree
private

Definition at line 114 of file Track3DKalmanSPS_module.cc.

The documentation for this class was generated from the following file:
Generated by   doxygen 1.8.11