LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
Track3DKalmanSPS_module.cc
Go to the documentation of this file.
1 //
3 // \file Track3DKalmanSPS.cxx
4 //
5 // \author echurch@fnal.gov
6 //
7 // This algorithm is designed to reconstruct 3D tracks through
8 // GENFIT's Kalman filter.
10 
11 // C++ includes
12 #include <algorithm> // std::sort()
13 #include <cmath>
14 #include <iterator> // std::distance()
15 #include <sstream>
16 #include <string>
17 #include <vector>
18 
19 // ROOT includes
20 #include "TDatabasePDG.h"
21 #include "TMath.h"
22 #include "TMatrixT.h"
23 #include "TPrincipal.h"
24 #include "TTree.h"
25 #include "TVector3.h"
26 
27 // Framework includes
28 #include "fhiclcpp/ParameterSet.h"
29 #include "fhiclcpp/exception.h"
31 
37 #include "art_root_io/TFileService.h"
41 
42 // GENFIT includes
47 #include "larreco/Genfit/GFTrack.h"
50 
51 // LArSoft includes
59 //\todo Reconstruction Producers should never include SimulationBase objects
62 
65 {
66  const double* xyz1 = h1->XYZ();
67  const double* xyz2 = h2->XYZ();
68  return xyz1[2] < xyz2[2];
69 }
72 {
73  const double* xyz1 = h1->XYZ();
74  const double* xyz2 = h2->XYZ();
75  return xyz1[1] < xyz2[1];
76 }
79 {
80  const double* xyz1 = h1->XYZ();
81  const double* xyz2 = h2->XYZ();
82  return xyz1[0] < xyz2[0];
83 }
84 
87 {
88  const unsigned int s1 = h1.size();
89  const unsigned int s2 = h2.size();
90  return s1 > s2;
91 }
92 
93 namespace trkf {
94 
96  public:
97  explicit Track3DKalmanSPS(fhicl::ParameterSet const& pset);
98 
99  private:
100  void produce(art::Event& evt);
101  void beginJob();
102  void endJob();
103  double energyLossBetheBloch(const double& mass, const double p);
104 
105  void rotationCov(TMatrixT<Double_t>& cov, const TVector3& u, const TVector3& v);
106  std::vector<double> dQdxCalc(const art::FindManyP<recob::Hit>& h,
108  const TVector3& p,
109  const TVector3& d);
110 
111  std::string fClusterModuleLabel; // label for input collection
112  std::string fSpptModuleLabel; // label for input collection
113  std::string fGenieGenModuleLabel; // label for input MC single particle generator
114  std::string fG4ModuleLabel; // label for input MC single particle generator
115  std::string fSortDim; // direction in which to sort spacepoints
116 
117  // TFile *fileGENFIT;
118  TTree* tree;
119 
120  TMatrixT<Double_t>* stMCT;
121  TMatrixT<Double_t>* covMCT;
122  TMatrixT<Double_t>* stREC;
123  TMatrixT<Double_t>* covREC;
124  Float_t chi2;
125  Float_t chi2ndf;
126  int fcont;
127 
128  Float_t* fpRECL;
129  Float_t* fpREC;
130  Float_t* fpMCMom;
131  Float_t* fpMCPos;
132  Float_t* fState0;
133  Float_t* fCov0;
134  int nfail;
135  int ndf;
137  int ispptvec;
138  int nspptvec;
139  unsigned int evtt;
140  unsigned int nTrks;
141  unsigned int fptsNo;
142  Float_t* fshx;
143  Float_t* fshy;
144  Float_t* fshz;
145  Float_t* feshx;
146  Float_t* feshy;
147  Float_t* feshz;
148  Float_t* feshyz;
149  Float_t* fupdate;
150  Float_t* fchi2hit;
151  Float_t* fth;
152  Float_t* feth;
153  Float_t* fedudw;
154  Float_t* fedvdw;
155  Float_t* feu;
156  Float_t* fev;
157  Float_t* fsep;
158  Float_t* fdQdx;
159  unsigned int fDimSize; // if necessary will get this from pset in constructor.
160  Float_t* fPCmeans;
161  Float_t* fPCevals;
162  Float_t* fPCsigmas;
163  Float_t* fPC1;
164  Float_t* fPC2;
165  Float_t* fPC3;
166 
167  std::vector<double> fPosErr;
168  std::vector<double> fMomErr;
169  std::vector<double> fMomStart;
170  double fPerpLim;
171  bool fDoFit;
172  int fNumIt;
173  uint16_t fMinNumSppts;
174  double fErrScaleS;
175  double fErrScaleM;
177  double fMaxUpdate;
179  double fDistanceU;
180  double fMaxUpdateU;
181  double fMomLow;
182  double fMomHigh;
183  int fPdg;
184  double fChi2Thresh;
185  int fMaxPass;
186 
189 
190  }; // class Track3DKalmanSPS
191 
192 }
193 
194 namespace trkf {
195 
196  //-------------------------------------------------
198  : EDProducer{pset}
199  , fDoFit(true)
200  , fNumIt(5)
201  , fMinNumSppts(5)
202  , fErrScaleS(1.0)
203  , fErrScaleM(1.0)
204  , fDecimate(1)
205  , fMaxUpdate(0.10)
206  , fDecimateU(1)
207  , fDistanceU(1.)
208  , fMaxUpdateU(0.10)
209  , fMomLow(0.001)
210  , fMomHigh(100.)
211  , fPdg(-13)
212  , fChi2Thresh(12.0E12)
213  , fMaxPass(1)
214  {
215  fClusterModuleLabel = pset.get<std::string>("ClusterModuleLabel");
216  fSpptModuleLabel = pset.get<std::string>("SpptModuleLabel");
217  fGenieGenModuleLabel = pset.get<std::string>("GenieGenModuleLabel");
218  fG4ModuleLabel = pset.get<std::string>("G4ModuleLabel");
219  fPosErr = pset.get<std::vector<double>>("PosErr3"); // resolution. cm
220  fMomErr = pset.get<std::vector<double>>("MomErr3"); // GeV
221  fMomStart = pset.get<std::vector<double>>("MomStart3"); //
222  fPerpLim = pset.get<double>("PerpLimit", 1.e6); // PCA cut.
223  fDoFit = pset.get<bool>("DoFit", true); // Der.
224  fNumIt = pset.get<int>("NumIt", 5); // Number x2 passes per fit.
225  fMinNumSppts =
226  pset.get<int>("MinNumSppts", 5); // Min number of sppts in vector to bother fitting
227  fErrScaleS = pset.get<double>("ErrScaleSim", 1.0); // error scale.
228  fErrScaleM = pset.get<double>("ErrScaleMeas", 1.0); // error scale.
229  fDecimate = pset.get<int>("DecimateC", 40); // Sparsify data.
230  fMaxUpdate = pset.get<double>("MaxUpdateC", 0.1); // 0-out.
231  fDecimateU = pset.get<int>("DecimateU", 100); // Sparsify data.
232  fDistanceU =
233  pset.get<double>("DistanceU", 10.0); // Require this separation on uncontained 2nd pass.
234  fMaxUpdateU = pset.get<double>("MaxUpdateU", 0.02); // 0-out.
235  fMomLow = pset.get<double>("MomLow", 0.01); // Fit Range.
236  fMomHigh = pset.get<double>("MomHigh", 20.); // Fit Range.
237  fPdg = pset.get<int>("PdgCode", -13); // mu+ Hypothesis.
238  fChi2Thresh = pset.get<double>("Chi2HitThresh", 12.0E12); //For Re-pass.
239  fSortDim = pset.get<std::string>("SortDirection", "z"); // case sensitive
240  fMaxPass = pset.get<int>("MaxPass", 2); // mu+ Hypothesis.
241  bool fGenfPRINT;
242  if (pset.get_if_present("GenfPRINT", fGenfPRINT)) {
243  MF_LOG_WARNING("Track3DKalmanSPS_GenFit")
244  << "Parameter 'GenfPRINT' has been deprecated.\n"
245  "Please use the standard message facility to enable GenFit debug output.";
246  // A way to enable debug output is all of the following:
247  // - compile in debug mode (no optimization, no profiling)
248  // - if that makes everything too noisy, add to have everything else quiet
249  // services.message.debugModules: [ "Track3DKalmanSPS" ]
250  // - to print all the GenFit debug messages, set
251  // services.message.destinations.LogDebugFile.categories.Track3DKalmanSPS_GenFit.limit: -1
252  // (assuming there is a LogDebugFile destination already; for example
253  // see the settings in uboonecode/uboone/Utilities/services_microboone.fcl )
254  }
255 
256  produces<std::vector<recob::Track>>();
257  produces<art::Assns<recob::Track, recob::Cluster>>();
258  produces<art::Assns<recob::Track, recob::SpacePoint>>();
259  produces<art::Assns<recob::Track, recob::Hit>>();
260  }
261 
262  //-------------------------------------------------
263  // stolen, mostly, from GFMaterialEffects.
264  double Track3DKalmanSPS::energyLossBetheBloch(const double& mass, const double p = 1.5)
265  {
266  const double charge(1.0);
267  const double mEE(188.); // eV
268  const double matZ(18.);
269  const double matA(40.);
270  const double matDensity(1.4);
271  const double me(0.000511);
272 
273  double beta = p / std::sqrt(mass * mass + p * p);
274  double gammaSquare = 1. / (1.0 - beta * beta);
275  // 4pi.r_e^2.N.me = 0.307075, I think.
276  double dedx = 0.307075 * matDensity * matZ / matA / (beta * beta) * charge * charge;
277  double massRatio = me / mass;
278  // me=0.000511 here is in GeV. So mEE comes in here in eV.
279  double argument = gammaSquare * beta * beta * me * 1.E3 * 2. /
280  ((1.E-6 * mEE) * std::sqrt(1 + 2 * std::sqrt(gammaSquare) * massRatio +
281  massRatio * massRatio));
282 
283  if (mass == 0.0) return (0.0);
284  if (argument <= exp(beta * beta)) { dedx = 0.; }
285  else {
286  dedx *= (log(argument) - beta * beta); // Bethe-Bloch [MeV/cm]
287  dedx *= 1.E-3; // in GeV/cm, hence 1.e-3
288  if (dedx < 0.) dedx = 0.;
289  }
290  return dedx;
291  }
292 
293  void Track3DKalmanSPS::rotationCov(TMatrixT<Double_t>& cov, const TVector3& u, const TVector3& v)
294  {
295  TVector3 xhat(1.0, 0.0, 0.0);
296  TVector3 yhat(0.0, 1.0, 0.0);
297  TVector3 zhat(0.0, 0.0, 1.0);
298  TVector3 w(u.Cross(v));
299  TVector3 uprime(u);
300  TVector3 vprime(w.Cross(xhat)); // vprime now lies in yz plane
301  Double_t angle(v.Angle(vprime)); /* This is the angle through which v
302  must rotate. */
303  uprime.Rotate(angle, w); // u now is rotated the same amount
304  if (uprime * xhat < 0) {
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  }
324 
327  const TVector3& dir,
328  const TVector3& loc)
329  {
330  // For now just Collection plane.
331  // We should loop over all views, more generally.
335  // art::PtrVector<recob::SpacePoint>::const_iterator sppt = sstart;
337  std::vector<double> v;
338 
339  double mindist(100.0); // cm
340  auto spptminIt(sppt);
341  while (sppt != s.end()) {
342  if (((**sppt).XYZ() - loc).Mag() < mindist) {
343  double dist = ((**sppt).XYZ() - loc).Mag();
344 
345  // Jump out if we're as close as 0.1 mm away.
346  if (dist < mindist) {
347  mindist = dist;
348  spptminIt = sppt;
349  if (mindist < 0.01) break;
350  }
351  }
352  sppt++;
353  }
354  sstart = spptminIt; // for next time.
355  unsigned int ind(std::distance(s.begin(), spptminIt));
356 
357  std::vector<art::Ptr<recob::Hit>> hitlist = h.at(ind);
358 
359  double wirePitch = 0.;
360  double angleToVert = 0;
361  double charge = 0.;
362 
363  // FIXME: Currently the last iteration sets the charge, wirePitch, and angleToVert variables.
364  // Is this expected?
365  for (std::vector<art::Ptr<recob::Hit>>::const_iterator ihit = hitlist.begin();
366  ihit != hitlist.end();
367  ++ihit) {
368  const recob::Hit& hit1 = **ihit;
369  if (hit1.SignalType() != sig) continue;
370  geo::PlaneID const& hit1PlaneID = hit1.WireID().asPlaneID();
371  charge = hit1.Integral();
372  wirePitch = geom->WirePitch(hit1PlaneID);
373  angleToVert = geom->Plane(hit1PlaneID).Wire(0).ThetaZ(false) - 0.5 * TMath::Pi();
374  }
375 
376  double cosgamma =
377  TMath::Abs(TMath::Sin(angleToVert) * dir.Y() + TMath::Cos(angleToVert) * dir.Z());
378 
379  v.push_back(charge / wirePitch / cosgamma);
380 
381  return v;
382  }
383 
384  //-------------------------------------------------
386  {
387 
389 
390  stMCT = new TMatrixT<Double_t>(5, 1);
391  covMCT = new TMatrixT<Double_t>(5, 5);
392  stREC = new TMatrixT<Double_t>(5, 1);
393  covREC = new TMatrixT<Double_t>(5, 5);
394 
395  fpMCMom = new Float_t[4];
396  fpMCPos = new Float_t[4];
397  fpREC = new Float_t[4];
398  fpRECL = new Float_t[4];
399  fState0 = new Float_t[5];
400  fCov0 = new Float_t[25];
401  fDimSize = 20000; // if necessary will get this from pset in constructor.
402 
403  fshx = new Float_t[fDimSize];
404  fshy = new Float_t[fDimSize];
405  fshz = new Float_t[fDimSize];
406  feshx = new Float_t[fDimSize];
407  feshy = new Float_t[fDimSize];
408  feshz = new Float_t[fDimSize];
409  feshyz = new Float_t[fDimSize];
410  fupdate = new Float_t[fDimSize];
411  fchi2hit = new Float_t[fDimSize];
412  fth = new Float_t[fDimSize];
413  feth = new Float_t[fDimSize];
414  fedudw = new Float_t[fDimSize];
415  fedvdw = new Float_t[fDimSize];
416  feu = new Float_t[fDimSize];
417  fev = new Float_t[fDimSize];
418  fsep = new Float_t[fDimSize];
419  fdQdx = new Float_t[fDimSize];
420 
421  fPC1 = new Float_t[3];
422  fPC2 = new Float_t[3];
423  fPC3 = new Float_t[3];
424  fPCmeans = new Float_t[3];
425  fPCsigmas = new Float_t[3];
426  fPCevals = new Float_t[3];
427 
428  // //TFile fileGENFIT("GENFITout.root","RECREATE");
429 
430  tree = tfs->make<TTree>("GENFITttree", "GENFITttree");
431  //tree->Branch("stMCT",&stMCT,"stMCT[5]/F"); // "TMatrixT<Double_t>"
432 
433  //tree->Branch("stMCT","TMatrixD",&stMCT,64000,0);
434  //tree->Branch("covMCT",covMCT,"covMCT[25]/F");
435  tree->Branch("covMCT", "TMatrixD", &covMCT, 64000, 0);
436  tree->Branch("stREC", fState0, "stREC[5]/F");
437  //tree->Branch("stREC","TMatrixD",&stREC,64000,0);
438  tree->Branch("covREC", fCov0, "covREC[25]/F");
439  //tree->Branch("covREC","TMatrixD",&covREC,64000,0);
440 
441  tree->Branch("nchi2rePass", &nchi2rePass, "nchi2rePass/I");
442  tree->Branch("ispptvec", &ispptvec, "ispptvec/I");
443  tree->Branch("chi2", &chi2, "chi2/F");
444  tree->Branch("nfail", &nfail, "nfail/I");
445  tree->Branch("ndf", &ndf, "ndf/I");
446  tree->Branch("evtNo", &evtt, "evtNo/I");
447  tree->Branch("nspptvec", &nspptvec, "nspptvec/I");
448  tree->Branch("chi2ndf", &chi2ndf, "chi2ndf/F");
449 
450  tree->Branch("trkNo", &nTrks, "trkNo/I");
451  tree->Branch("ptsNo", &fptsNo, "ptsNo/I");
452  tree->Branch("cont", &fcont, "cont/I"); //O? Yes, O. Not 0, not L, ...
453  tree->Branch("shx", fshx, "shx[ptsNo]/F");
454  tree->Branch("shy", fshy, "shy[ptsNo]/F");
455  tree->Branch("shz", fshz, "shz[ptsNo]/F");
456  tree->Branch("sep", fsep, "sep[ptsNo]/F");
457  tree->Branch("dQdx", fdQdx, "dQdx[ptsNo]/F");
458  tree->Branch("eshx", feshx, "eshx[ptsNo]/F");
459  tree->Branch("eshy", feshy, "eshy[ptsNo]/F");
460  tree->Branch("eshz", feshz, "eshz[ptsNo]/F");
461  tree->Branch("eshyz", feshyz, "eshyz[ptsNo]/F");
462  tree->Branch("update", fupdate, "update[ptsNo]/F");
463  tree->Branch("chi2hit", fchi2hit, "chi2hit[ptsNo]/F");
464  tree->Branch("th", fth, "th[ptsNo]/F");
465  tree->Branch("eth", feth, "eth[ptsNo]/F");
466  tree->Branch("edudw", fedudw, "edudw[ptsNo]/F");
467  tree->Branch("edvdw", fedvdw, "edvdw[ptsNo]/F");
468  tree->Branch("eu", feu, "eu[ptsNo]/F");
469  tree->Branch("ev", fev, "ev[ptsNo]/F");
470 
471  tree->Branch("pcMeans", fPCmeans, "pcMeans[3]/F");
472  tree->Branch("pcSigmas", fPCsigmas, "pcSigmas[3]/F");
473  tree->Branch("pcEvals", fPCevals, "pcEvals[3]/F");
474  tree->Branch("pcEvec1", fPC1, "pcEvec1[3]/F");
475  tree->Branch("pcEvec2", fPC2, "pcEvec2[3]/F");
476  tree->Branch("pcEvec3", fPC3, "pcEvec3[3]/F");
477 
478  tree->Branch("pMCMom", fpMCMom, "pMCMom[4]/F");
479  tree->Branch("pMCPos", fpMCPos, "pMCPos[4]/F");
480  tree->Branch("pRECKalF", fpREC, "pRECKalF[4]/F");
481  tree->Branch("pRECKalL", fpRECL, "pRECKalL[4]/F");
482 
483  //TGeoManager* geomGENFIT = new TGeoManager("Geometry", "Geane geometry");
484  //TGeoManager::Import("config/genfitGeom.root");
485  // gROOT->Macro("config/Geane.C");
486  }
487 
488  //-------------------------------------------------
490  {
491  if (!rep) delete rep;
492  if (!repMC) delete repMC;
493 
494  /*
495  // not sure why I can't do these, but at least some cause seg faults.
496  delete[] stMCT;
497  delete[] covMCT;
498  delete[] stREC;
499  delete[] covREC;
500  */
501 
502  delete[] fpREC;
503  delete[] fpRECL;
504  delete[] fState0;
505  delete[] fCov0;
506 
507  delete[] fshx;
508  delete[] fshy;
509  delete[] fshz;
510  delete[] feshx;
511  delete[] feshy;
512  delete[] feshyz;
513  delete[] feshz;
514  delete[] fupdate;
515  delete[] fchi2hit;
516  delete[] fth;
517  delete[] feth;
518  delete[] fedudw;
519  delete[] fedvdw;
520  delete[] feu;
521  delete[] fev;
522  delete[] fsep;
523  delete[] fdQdx;
524 
525  delete[] fPCmeans;
526  delete[] fPCsigmas;
527  delete[] fPCevals;
528  delete[] fPC1;
529  delete[] fPC2;
530  delete[] fPC3;
531  }
532 
533  //------------------------------------------------------------------------------------//
535  {
536 
537  rep = 0;
538  repMC = 0;
539  // get services
541 
543  // Make a std::unique_ptr<> for the thing you want to put into the event
544  // because that handles the memory management for you
546  std::unique_ptr<std::vector<recob::Track>> tcol(new std::vector<recob::Track>);
547  std::unique_ptr<art::Assns<recob::Track, recob::SpacePoint>> tspassn(
549  std::unique_ptr<art::Assns<recob::Track, recob::Hit>> thassn(
551 
552  unsigned int tcnt = 0;
553 
554  // define TPC parameters
555  TString tpcName = geom->GetLArTPCVolumeName();
556 
557  // get input Hit object(s).
558  art::Handle<std::vector<recob::Cluster>> clusterListHandle;
559  evt.getByLabel(fClusterModuleLabel, clusterListHandle);
560 
562  evt.getByLabel(fSpptModuleLabel, spptListHandle);
563 
565 
569  if (!evt.isRealData()) {
570 
571  art::Handle<std::vector<simb::MCTruth>> mctruthListHandle;
572  evt.getByLabel(fGenieGenModuleLabel, mctruthListHandle);
573 
574  for (unsigned int ii = 0; ii < mctruthListHandle->size(); ++ii) {
575  art::Ptr<simb::MCTruth> mctparticle(mctruthListHandle, ii);
576  mclist.push_back(mctparticle);
577  }
578  }
579 
580  std::vector<art::PtrVector<recob::SpacePoint>> spptIn(spptListHandle->begin(),
581  spptListHandle->end());
582  // Get the spptvectors that are largest to be first, and smallest last.
583  std::sort(spptIn.begin(), spptIn.end(), sp_sort_nsppts);
584 
585  std::vector<art::PtrVector<recob::SpacePoint>>::const_iterator sppt = spptIn.begin();
586  auto spptB = sppt;
587 
588  TVector3 MCOrigin;
589  TVector3 MCMomentum;
590  // TVector3 posErr(.05,.05,.05); // resolution. 0.5mm
591  // TVector3 momErr(.1,.1,0.2); // GeV
592  TVector3 posErr(fPosErr[0], fPosErr[1], fPosErr[2]); // resolution. 0.5mm
593  TVector3 momErr(fMomErr[0], fMomErr[1], fMomErr[2]); // GeV
594  TVector3 momErrFit(fMomErr[0], fMomErr[1], fMomErr[2]); // GeV
595 
596  // This is strictly for MC
600  if (!evt.isRealData()) {
601  // Below breaks are stupid, I realize. But rather than keep all the MC
602  // particles I just take the first primary, e.g., muon and only keep its
603  // info in the branches of the Ttree. I could generalize later, ...
604  for (unsigned int ii = 0; ii < mclist.size(); ++ii) {
605  art::Ptr<simb::MCTruth> mc(mclist[ii]);
606  for (int jj = 0; jj < mc->NParticles(); ++jj) {
608  MCOrigin.SetXYZ(part.Vx(), part.Vy(), part.Vz()); // V for Vertex
609  MCMomentum.SetXYZ(part.Px(), part.Py(), part.Pz());
610  MF_LOG_DEBUG("Track3DKalmanSPS_GenFit")
611  << "FROM MC TRUTH, the particle's pdg code is: " << part.PdgCode()
612  << " with energy = " << part.E() << ", with energy = " << part.E()
613  << "\n vtx: " << genf::ROOTobjectToString(MCOrigin)
614  << "\n momentum: " << genf::ROOTobjectToString(MCMomentum)
615  << "\n (both in Global (not volTPC) coords)";
616 
617  repMC = new genf::RKTrackRep(MCOrigin, MCMomentum, posErr, momErr, part.PdgCode());
618  break;
619  }
620  break;
621  }
622  //for saving of MC truth
623  stMCT->ResizeTo(repMC->getState());
624  *stMCT = repMC->getState();
625  covMCT->ResizeTo(repMC->getCov());
626  *covMCT = repMC->getCov();
627  MF_LOG_DEBUG("Track3DKalmanSPS_GenFit")
628  << " repMC, covMC are ... \n"
630 
631  } // !isRealData
632  nTrks = 0;
633  TParticlePDG* part = TDatabasePDG::Instance()->GetParticle(fPdg);
634  Double_t mass = part->Mass();
635 
636  size_t cntp(0);
637  while (sppt != spptIn.end()) {
638 
639  const art::PtrVector<recob::SpacePoint>& spacepoints = *sppt;
640 
641  double fMaxUpdateHere(fMaxUpdateU);
642  int fDecimateHere(fDecimateU);
643  double fErrScaleSHere(fErrScaleS);
644  double fErrScaleMHere(fErrScaleM);
645  int rePass0(1);
646  unsigned int nTailPoints = 0; // 100;
647  if (spacepoints.size() < 5) {
648  sppt++;
649  rePass0 = 3;
650  continue;
651  } // for now...
652 
653  MF_LOG_DEBUG("Track3DKalmanSPS_GenFit")
654  << "\n\t found " << spacepoints.size()
655  << " 3D spacepoint(s) for this element of std::vector<art:PtrVector> spacepoints. \n";
656 
657  // Let's find track's principle components.
658  // We will sort along that direction, rather than z.
659  // Further, we will skip outliers away from main axis.
660 
661  TPrincipal* principal = new TPrincipal(3, "ND");
662 
663  // I need to shuffle these around, so use copy constructor
664  // to make non-const version spacepointss.
665  art::PtrVector<recob::SpacePoint> spacepointss(spacepoints);
666 
667  // What I need is a nearest neighbor sorting.
668  if (fSortDim.compare("y") && fSortDim.compare("x"))
669  std::sort(spacepointss.begin(), spacepointss.end(), sp_sort_3dz);
670  if (!fSortDim.compare("y")) std::sort(spacepointss.begin(), spacepointss.end(), sp_sort_3dy);
671  if (!fSortDim.compare("x")) std::sort(spacepointss.begin(), spacepointss.end(), sp_sort_3dx);
672 
673  for (unsigned int point = 0; point < spacepointss.size(); ++point) {
674  if (point < (spacepointss.size() - nTailPoints)) {
675  principal->AddRow(spacepointss[point]->XYZ());
676  }
677  }
678  principal->MakePrincipals();
679 
680  const TVectorD* evals = principal->GetEigenValues();
681  const TMatrixD* evecs = principal->GetEigenVectors();
682  const TVectorD* means = principal->GetMeanValues();
683  const TVectorD* sigmas = principal->GetSigmas();
684 
685  Double_t tmp[3], tmp2[3];
686  principal->X2P((Double_t*)(means->GetMatrixArray()), tmp);
687  principal->X2P((Double_t*)(sigmas->GetMatrixArray()), tmp2);
688  for (unsigned int ii = 0; ii < 3; ++ii) {
689  fPCmeans[ii] = (Float_t)(tmp[ii]);
690  fPCsigmas[ii] = (Float_t)(tmp2[ii]);
691  fPCevals[ii] = (Float_t)(evals->GetMatrixArray())[ii];
692  // This method requires apparently pulling all 9
693  // elements. Maybe 3 works.
694  // Certainly, w can't be a scalar, I discovered.
695  double w[9];
696  evecs->ExtractRow(ii, 0, w);
697  fPC1[ii] = w[0];
698  fPC2[ii] = w[1];
699  fPC3[ii] = w[2];
700  }
701  Double_t tmp3[3], tmp4[3], tmp5[3];
702  principal->X2P((Double_t*)fPC1, tmp3);
703  principal->X2P((Double_t*)fPC2, tmp4);
704  principal->X2P((Double_t*)fPC3, tmp5);
705 
706  // Use a mip approximation assuming straight lines
707  // and a small angle wrt beam.
708  fMomStart[0] = spacepointss[spacepointss.size() - 1]->XYZ()[0] - spacepointss[0]->XYZ()[0];
709  fMomStart[1] = spacepointss[spacepointss.size() - 1]->XYZ()[1] - spacepointss[0]->XYZ()[1];
710  fMomStart[2] = spacepointss[spacepointss.size() - 1]->XYZ()[2] - spacepointss[0]->XYZ()[2];
711  // This presumes a 0.8 GeV/c particle
712  double dEdx = energyLossBetheBloch(mass, 1.0);
713  // mom is really KE.
714  TVector3 mom(dEdx * fMomStart[0], dEdx * fMomStart[1], dEdx * fMomStart[2]);
715  double pmag2 = pow(mom.Mag() + mass, 2. - mass * mass);
716  mom.SetMag(std::sqrt(pmag2));
717  // Over-estimate by just enough for contained particles (5%).
718  mom.SetMag(1.0 * mom.Mag());
719  // If 1st/last point is close to edge of TPC, this track is
720  // uncontained.Give higher momentum starting value in
721  // that case.
722  bool uncontained(false);
723  double close(5.); // cm.
724  double epsMag(0.001); // cm.
725  double epsX(250.0); // cm.
726  double epsZ(0.001); // cm.
727 
728  if (spacepointss[spacepointss.size() - 1]->XYZ()[0] > (2. * geom->DetHalfWidth() - close) ||
729  spacepointss[spacepointss.size() - 1]->XYZ()[0] < close ||
730  spacepointss[0]->XYZ()[0] > (2. * geom->DetHalfWidth() - close) ||
731  spacepointss[0]->XYZ()[0] < close ||
732  spacepointss[spacepointss.size() - 1]->XYZ()[1] > (1. * geom->DetHalfHeight() - close) ||
733  (spacepointss[spacepointss.size() - 1]->XYZ()[1] < -1. * geom->DetHalfHeight() + close) ||
734  spacepointss[0]->XYZ()[1] > (1. * geom->DetHalfHeight() - close) ||
735  spacepointss[0]->XYZ()[1] < (-1. * geom->DetHalfHeight() + close) ||
736  spacepointss[spacepointss.size() - 1]->XYZ()[2] > (geom->DetLength() - close) ||
737  spacepointss[spacepointss.size() - 1]->XYZ()[2] < close ||
738  spacepointss[0]->XYZ()[2] > (geom->DetLength() - close) ||
739  spacepointss[0]->XYZ()[2] < close)
740  uncontained = true;
741  fErrScaleSHere = fErrScaleS;
742  fErrScaleMHere = fErrScaleM;
743 
744  if (uncontained) {
745  // Big enough to not run out of gas right at end of
746  // track and give large angular deviations which
747  // will kill the fit.
748  mom.SetMag(2.0 * mom.Mag());
749  MF_LOG_DEBUG("Track3DKalmanSPS_GenFit") << "Uncontained track ... ";
750  fDecimateHere = fDecimateU;
751  fMaxUpdateHere = fMaxUpdateU;
752  }
753  else {
754  MF_LOG_DEBUG("Track3DKalmanSPS_GenFit")
755  << "Contained track ... Run " << evt.run() << " Event " << evt.id().event();
756  // Don't decimate contained tracks as drastically,
757  // and omit only very large corrections ...
758  // which hurt only high momentum tracks.
759  fDecimateHere = fDecimate;
760  fMaxUpdateHere = fMaxUpdate;
761  }
762  fcont = (int)(!uncontained);
763 
764  // This seems like best place to jump back to for a re-pass.
765  unsigned short rePass = rePass0; // 1 by default;
766  unsigned short maxPass(fMaxPass);
767  unsigned short tcnt1(0);
768  while (rePass <= maxPass) {
769 
770  TVector3 momM(mom);
771  TVector3 momErrFit(momM[0] / 3.0, momM[1] / 3.0,
772  momM[2] / 3.0); // GeV
773 
775  genf::GFDetPlane planeG((TVector3)(spacepointss[0]->XYZ()), momM);
776 
777  // Initialize with 1st spacepoint location and ...
778  rep = new genf::RKTrackRep((TVector3)(spacepointss[0]->XYZ()),
779  momM,
780  posErr,
781  momErrFit,
782  fPdg); // mu+ hypothesis
783 
784  genf::GFTrack fitTrack(rep); //initialized with smeared rep
785  fitTrack.setPDG(fPdg);
786  // Gonna sort in z cuz I want to essentially transform here to volTPC coords.
787  // volTPC coords, cuz that's what the Geant3/Geane stepper wants, as that's its understanding
788  // from the Geant4 geometry, which it'll use. EC, 7-Jan-2011.
789  int ihit = 0;
790  fptsNo = 0;
791  std::vector<unsigned int> spptSurvivedIndex;
792  std::vector<unsigned int> spptSkippedIndex;
793  unsigned int ppoint(0);
794  for (unsigned int point = 0; point < spacepointss.size(); ++point) {
795  double sep;
796  // Calculate the distance in 2nd and 3rd PCs and
797  // reject spt if it's too far out. Remember, the
798  // sigmas are std::sqrt(eigenvals).
799  double tmp[3];
800  principal->X2P((Double_t*)(spacepointss[point]->XYZ()), tmp);
801  sep = std::sqrt(tmp[1] * tmp[1] / fPCevals[1] + tmp[2] * tmp[2] / fPCevals[2]);
802  if ((std::abs(sep) > fPerpLim) && (point < (spacepointss.size() - nTailPoints)) &&
803  rePass <= 1) {
804  spptSkippedIndex.push_back(point);
805  continue;
806  }
807  // If point is too close in Mag or Z or too far in X from last kept point drop it.
808  // I think this is largely redundant with PCA cut.
809  TVector3 one(spacepointss[point]->XYZ());
810  TVector3 two(spacepointss[ppoint]->XYZ());
811  if (rePass == 2 && uncontained) {
812  epsMag = fDistanceU; // cm
813  fNumIt = 4;
814  fErrScaleMHere = 0.1;
815  // Above allows us to pretend as though measurements
816  // are perfect, which we can ostensibly do now with
817  // clean set of sppts. This creates larger gains, bigger
818  // updates: bigger sensitivity to multiple scattering.
819  }
820  else if (rePass == 2 && !uncontained) {}
821  if (point > 0 &&
822  ((one - two).Mag() < epsMag || // too close
823  ((one - two).Mag() > 8.0 && rePass == 1) || // too far
824  std::abs(spacepointss[point]->XYZ()[2] - spacepointss[ppoint]->XYZ()[2]) < epsZ ||
825  std::abs(spacepointss[point]->XYZ()[0] - spacepointss[ppoint]->XYZ()[0]) > epsX)) {
826  spptSkippedIndex.push_back(point);
827  continue;
828  }
829 
830  if (
831  point % fDecimateHere &&
832  rePass <=
833  1) // Jump out of loop except on every fDecimate^th pt. fDecimate==1 never sees continue.
834  {
835  /* Replace continue with a counter that will be used to index into vector of GFKalman fits. */
836  continue;
837  }
838 
839  ppoint = point;
840  TVector3 spt3 = (TVector3)(spacepointss[point]->XYZ());
841  std::vector<double> err3;
842  err3.push_back(spacepointss[point]->ErrXYZ()[0]);
843  err3.push_back(spacepointss[point]->ErrXYZ()[2]);
844  err3.push_back(spacepointss[point]->ErrXYZ()[4]);
845  err3.push_back(spacepointss[point]->ErrXYZ()[5]); // lower triangle diags.
846  if (fptsNo < fDimSize) {
847  fshx[fptsNo] = spt3[0];
848  fshy[fptsNo] = spt3[1];
849  fshz[fptsNo] = spt3[2];
850  feshx[fptsNo] = err3[0];
851  feshy[fptsNo] = err3[1];
852  feshz[fptsNo] = err3[3];
853  feshyz[fptsNo] = err3[2];
854  fsep[fptsNo] = sep;
855 
856  if (fptsNo > 1) {
857  TVector3 pointer(fshx[fptsNo] - fshx[fptsNo - 1],
858  fshy[fptsNo] - fshy[fptsNo - 1],
859  fshz[fptsNo] - fshz[fptsNo - 1]);
860  TVector3 pointerPrev(fshx[fptsNo - 1] - fshx[fptsNo - 2],
861  fshy[fptsNo - 1] - fshy[fptsNo - 2],
862  fshz[fptsNo - 1] - fshz[fptsNo - 2]);
863  fth[fptsNo] = (pointer.Unit()).Angle(pointerPrev.Unit());
864  }
865  feth[fptsNo] = 0.0;
866  fedudw[fptsNo] = 0.0;
867  fedvdw[fptsNo] = 0.0;
868  feu[fptsNo] = 0.0;
869  fev[fptsNo] = 0.0;
870  fupdate[fptsNo] = 0.0;
871  }
872 
873  MF_LOG_DEBUG("Track3DKalmanSPS_GenFit")
874  << "ihit xyz..." << spt3[0] << "," << spt3[1] << "," << spt3[2];
875 
876  fitTrack.addHit(new genf::PointHit(spt3, err3),
877  1, //dummy detector id
878  ihit++);
879  spptSurvivedIndex.push_back(point);
880  fptsNo++;
881  } // end loop over spacepoints.
882 
883  if (fptsNo <=
884  fMinNumSppts) // Cuz 1st 2 in each direction don't count. Should have, say, 3 more.
885  {
886  MF_LOG_DEBUG("Track3DKalmanSPS_GenFit")
887  << "Bailing cuz only " << fptsNo << " spacepoints.";
888  rePass++;
889  continue;
890  }
891  MF_LOG_DEBUG("Track3DKalmanSPS_GenFit") << "Fitting on " << fptsNo << " spacepoints.";
892 
893  genf::GFKalman k;
894  k.setBlowUpFactor(5); // 500 out of box. EC, 6-Jan-2011.
895  k.setMomHigh(fMomHigh); // Don't fit above this many GeV.
896  k.setMomLow(fMomLow); // Don't fit below this many GeV.
897 
898  k.setInitialDirection(+1); // Instead of 1 out of box. EC, 6-Jan-2011.
900  k.setMaxUpdate(fMaxUpdateHere); // 0 out abs(update) bigger than this.
901  k.setErrorScaleSTh(fErrScaleSHere);
902  k.setErrorScaleMTh(fErrScaleMHere);
903 
904  bool skipFill = false;
905  std::vector<TMatrixT<double>> hitMeasCov;
906  std::vector<TMatrixT<double>> hitUpdate;
907  std::vector<TMatrixT<double>> hitCov;
908  std::vector<TMatrixT<double>> hitCov7x7;
909  std::vector<TMatrixT<double>> hitState;
910  std::vector<double> hitChi2;
911  std::vector<TVector3> hitPlaneXYZ;
912  std::vector<TVector3> hitPlaneUxUyUz;
913  std::vector<TVector3> hitPlaneU;
914  std::vector<TVector3> hitPlaneV;
915 
916  try {
917  if (fDoFit) k.processTrack(&fitTrack);
918  }
919  catch (cet::exception& e) {
920  MF_LOG_ERROR("Track3DKalmanSPS")
921  << "just caught a cet::exception: " << e.what()
922  << "\nExceptions won't be further handled; skip filling big chunks of the TTree.";
923  skipFill = true;
924  }
925 
926  if (rep->getStatusFlag() == 0) // 0 is successful completion
927  {
928  if (mf::isDebugEnabled()) {
929 
930  std::ostringstream dbgmsg;
931  dbgmsg << "Original plane:";
932  planeG.Print(dbgmsg);
933 
934  dbgmsg << "Current (fit) reference Plane:";
935  rep->getReferencePlane().Print(dbgmsg);
936 
937  dbgmsg << "Last reference Plane:";
938  rep->getLastPlane().Print(dbgmsg);
939 
940  if (planeG != rep->getReferencePlane())
941  dbgmsg << " => original hit plane (not surprisingly) not current reference Plane!";
942 
943  MF_LOG_DEBUG("Track3DKalmanSPS_GenFit") << dbgmsg.str();
944  }
945  if (!skipFill) {
946  hitMeasCov = fitTrack.getHitMeasuredCov();
947  hitUpdate = fitTrack.getHitUpdate();
948  hitCov = fitTrack.getHitCov();
949  hitCov7x7 = fitTrack.getHitCov7x7();
950  hitState = fitTrack.getHitState();
951  hitChi2 = fitTrack.getHitChi2();
952  hitPlaneXYZ = fitTrack.getHitPlaneXYZ();
953  hitPlaneUxUyUz = fitTrack.getHitPlaneUxUyUz();
954  hitPlaneU = fitTrack.getHitPlaneU();
955  hitPlaneV = fitTrack.getHitPlaneV();
956  unsigned int totHits = hitState.size();
957 
958  // Pick up info from last fwd Kalman pass.
959  unsigned int jhit = 0;
960  for (unsigned int ihit = totHits - 2 * totHits / (2 * fNumIt);
961  ihit < (totHits - totHits / (2 * fNumIt));
962  ihit++) // was ihit<ihit<(totHits-fptsNo)<7
963  {
964  feth[jhit] = (Float_t)(hitMeasCov.at(ihit)[0][0]); // eth
965  fedudw[jhit] = (Float_t)(hitMeasCov.at(ihit)[1][1]);
966  fedvdw[jhit] = (Float_t)(hitMeasCov.at(ihit)[2][2]);
967  feu[jhit] = (Float_t)(hitMeasCov.at(ihit)[3][3]);
968  fev[jhit] = (Float_t)(hitMeasCov.at(ihit)[4][4]);
969  fupdate[jhit] = (Float_t)(hitUpdate.at(ihit)[0][0]);
970  fchi2hit[jhit] = (Float_t)(hitChi2.at(ihit));
971  jhit++;
972  }
973 
974  stREC->ResizeTo(rep->getState());
975  *stREC = rep->getState();
976  covREC->ResizeTo(rep->getCov());
977  *covREC = rep->getCov();
978  double dum[5];
979  double dum2[5];
980  for (unsigned int ii = 0; ii < 5; ii++) {
981  stREC->ExtractRow(ii, 0, dum);
982  fState0[ii] = dum[0];
983  covREC->ExtractRow(ii, 0, dum2);
984  for (unsigned int jj = 0; jj < 5; jj++) {
985  fCov0[ii * 5 + jj] = dum2[jj];
986  }
987  }
988  MF_LOG_DEBUG("Track3DKalmanSPS_GenFit")
989  << " First State and Cov:" << genf::ROOTobjectToString(*stREC)
991  chi2 = (Float_t)(rep->getChiSqu());
992  ndf = rep->getNDF();
993  nfail = fitTrack.getFailedHits();
994  nchi2rePass = (int)rePass;
995  ispptvec = 1 + std::distance(spptB, sppt);
996  chi2ndf = (Float_t)(chi2 / ndf);
997 
998  nTrks++;
999  MF_LOG_DEBUG("Track3DKalmanSPS_GenFit")
1000  << "Track3DKalmanSPS about to do tree->Fill(). Chi2/ndf is " << chi2 / ndf << ".";
1001  fpMCMom[3] = MCMomentum.Mag();
1002  for (int ii = 0; ii < 3; ++ii) {
1003  fpMCMom[ii] = MCMomentum[ii];
1004  fpMCPos[ii] = MCOrigin[ii];
1005  fpREC[ii] = hitPlaneUxUyUz.at(totHits - 2 * totHits / (2 * fNumIt))[ii];
1006  fpRECL[ii] = hitPlaneUxUyUz.at(totHits - totHits / (2 * fNumIt) - 1)[ii];
1007  }
1008 
1009  evtt = (unsigned int)evt.id().event();
1010  nspptvec = (unsigned int)spptListHandle->size();
1011 
1012  cntp++;
1013  std::vector<std::vector<double>> dQdx;
1014  // Calculate LastFwdPass quantities.
1015  std::vector<TMatrixT<double>> hitCovLFP;
1016  std::vector<TVector3> hitPlaneXYZLFP;
1017  std::vector<TVector3> hitPlaneUxUyUzLFP;
1018  std::vector<TVector3> hitPlanePxPyPzLFP;
1019  std::vector<TVector3> hitPlaneULFP;
1020  std::vector<TVector3> hitPlaneVLFP;
1021  std::vector<double> pLFP;
1022  std::vector<TMatrixT<double>> c7x7LFP;
1023 
1024  art::FindManyP<recob::Hit> hitAssns(spacepointss, evt, fSpptModuleLabel);
1025  for (unsigned int ii = 0; ii < totHits / (2 * fNumIt); ii++) {
1026  pLFP.push_back(1. / hitState.at(totHits - 2 * totHits / (2 * fNumIt) + ii)[0][0]);
1027  // hitCov -> hitCov7x7 !! EC, 11-May-2012.
1028  c7x7LFP.push_back(hitCov7x7.at(totHits - 2 * totHits / (2 * fNumIt) + ii));
1029  hitCovLFP.push_back(hitCov.at(totHits - 2 * totHits / (2 * fNumIt) + ii));
1030  hitPlaneXYZLFP.push_back(hitPlaneXYZ.at(totHits - 2 * totHits / (2 * fNumIt) + ii));
1031  hitPlaneUxUyUzLFP.push_back(
1032  hitPlaneUxUyUz.at(totHits - 2 * totHits / (2 * fNumIt) + ii));
1033  hitPlanePxPyPzLFP.push_back(hitPlaneUxUyUzLFP.back() * pLFP.back());
1034  hitPlaneULFP.push_back(hitPlaneU.at(totHits - 2 * totHits / (2 * fNumIt) + ii));
1035  hitPlaneVLFP.push_back(hitPlaneV.at(totHits - 2 * totHits / (2 * fNumIt) + ii));
1036  // Transform cov appropriate for track rotated
1037  // about w, forcing
1038  // v to be in y-z plane and u pointing in
1039  // +-ive x direction, per TrackAna convention.
1040 
1041  rotationCov(hitCovLFP.back(), hitPlaneULFP.back(), hitPlaneVLFP.back());
1042  dQdx.push_back(
1043  dQdxCalc(hitAssns, spacepointss, hitPlaneUxUyUzLFP.back(), hitPlaneXYZLFP.back()));
1044  fdQdx[ii] = dQdx.back().back();
1045  }
1046  fpREC[3] = rep->getMom(rep->getReferencePlane()).Mag();
1047  fpRECL[3] = pLFP[1];
1048 
1049  tree->Fill();
1050 
1051  // Put newest track on stack for this set of sppts,
1052  // remove previous one.
1053  recob::tracking::SMatrixSym55 covVtx, covEnd;
1054  for (unsigned int i = 0; i < 5; i++) {
1055  for (unsigned int j = i; j < 5; j++) {
1056  covVtx(i, j) = hitCovLFP.front()(i, j);
1057  covEnd(i, j) = hitCovLFP.back()(i, j);
1058  }
1059  }
1060  recob::Track the3DTrack(
1062  recob::tracking::convertCollToVector(hitPlanePxPyPzLFP),
1063  recob::Track::Flags_t(hitPlaneXYZLFP.size()),
1064  true),
1065  0,
1066  -1.,
1067  0,
1068  covVtx,
1069  covEnd,
1070  tcnt++);
1071  if (rePass == 1) tcnt1++; // won't get here if Trackfit failed.
1072  if (rePass != 1 && tcnt1) tcol->pop_back();
1073  tcol->push_back(the3DTrack);
1074  util::CreateAssn(evt, *tcol, spacepointss, *tspassn);
1075  art::PtrVector<recob::Hit> hits; // = hitAssns;
1076  for (unsigned int ii = 0; ii < spacepointss.size(); ++ii) {
1077  hits.insert(hits.end(), hitAssns.at(ii).begin(), hitAssns.at(ii).end());
1078  }
1079  util::CreateAssn(evt, *tcol, hits, *thassn, tcol->size() - 1);
1080  } // end !skipFill
1081  } // getStatusFlag
1082 
1083  if (!rep) delete rep;
1084  rePass++;
1085  // need to first excise bad spacepoints.
1086  // Grab up large Chi2hits first.
1087  art::PtrVector<recob::SpacePoint> spacepointssExcise;
1088  for (unsigned int ind = 0; ind < spptSurvivedIndex.size(); ++ind) {
1089  // Stricter to chuck sppts from uncontained then contained trks.
1090  if ((uncontained && fchi2hit[ind] > fChi2Thresh) ||
1091  (!uncontained && fchi2hit[ind] > 1.e9) || fchi2hit[ind] < 0.0
1092  // =0 eliminates ruled-out large updates. Not obviously
1093  // helpful.
1094  // add a restriction on dQdx here ...
1095  ) {
1097  spacepointss.begin() + spptSurvivedIndex[ind];
1098  spacepointssExcise.push_back(*spptIt);
1099  }
1100  }
1101  // Now grab up those sppts which we skipped and don't want
1102  // to reconsider cuz they're too close to each other or
1103  // cuz they're too far in x, e.g.
1104  for (unsigned int ind = 0; ind < spptSkippedIndex.size(); ++ind) {
1106  spacepointss.begin() + spptSkippedIndex[ind];
1107  spacepointssExcise.push_back(*spptIt);
1108  }
1109  // Get rid of redundantly Excised sppts before proceeding.
1110  std::stable_sort(spacepointss.begin(), spacepointss.end());
1111  std::stable_sort(spacepointssExcise.begin(), spacepointssExcise.end());
1112  std::set_union(spacepointssExcise.begin(),
1113  spacepointssExcise.end(),
1114  spacepointssExcise.begin(),
1115  spacepointssExcise.end(),
1116  spacepointssExcise.begin());
1117  // Now excise. New spacepointss will be smaller for second pass.
1119  std::set_difference(spacepointss.begin(),
1120  spacepointss.end(),
1121  spacepointssExcise.begin(),
1122  spacepointssExcise.end(),
1123  spacepointss.begin());
1124  spacepointss.erase(diffSpptIt, spacepointss.end());
1125 
1126  // calculate new seed momentum, and errors as merited
1127  if (rePass == 2 /* && uncontained */) {
1128  if (fpREC[3] < fMomHigh && fpREC[3] > fMomLow) {
1129  double kick(0.9); //Try to get away with a smaller start
1130  // for contained tracks. While for uncontained tracks
1131  // let's start up at a higher momentum and come down.
1132  // Though, 2 (1) GeV/c tracks are too low (high), so
1133  // instead let's actually lower starting value on
1134  // this second pass. -- EC 7-Mar-2013
1135  if (uncontained) kick = 0.5;
1136  for (int ii = 0; ii < 3; ++ii) {
1137  mom[ii] = momM[ii] * kick;
1138  }
1139  }
1140  else if (uncontained) {
1141  double unstick(1.0);
1142  if (fpREC[3] >= fMomHigh) unstick = 0.3;
1143  for (int ii = 0; ii < 3; ++ii) {
1144  mom[ii] = momM[ii] * unstick;
1145  }
1146  }
1147  else
1148  for (int ii = 0; ii < 3; ++ii) {
1149  mom[ii] = 1.1 * momM[ii];
1150  }
1151  }
1152 
1153  } // end while rePass<=maxPass
1154 
1155  sppt++;
1156 
1157  } // end while on elements of std::vector of art::PtrVectors of SpPts.
1158 
1159  if (!repMC) delete repMC;
1160 
1161  evt.put(std::move(tcol));
1162  // and now the spacepoints
1163  evt.put(std::move(tspassn));
1164  // and the hits. Note that these are all the hits from all the spacepoints considered,
1165  // even though they're not all contributing to the tracks.
1166  evt.put(std::move(thassn));
1167  }
1168 
1170 
1171 } // end namespace
TMatrixT< Double_t > * stREC
unsigned int getNDF() const
std::vector< double > fMomErr
void setMomLow(Double_t f)
Definition: GFKalman.h:117
std::string GetLArTPCVolumeName(TPCID const &tpcid=tpc_zero) const
Return the name of specified LAr TPC volume.
std::vector< TVector3 > getHitPlaneUxUyUz()
Definition: GFTrack.h:373
void setMaxUpdate(Double_t f)
Definition: GFKalman.h:119
geo::SigType_t SignalType() const
Signal type for the plane of the hit.
Definition: Hit.h:272
TMatrixT< Double_t > * stMCT
static bool sp_sort_3dz(const art::Ptr< recob::SpacePoint > &h1, const art::Ptr< recob::SpacePoint > &h2)
typename data_t::iterator iterator
Definition: PtrVector.h:54
const GFDetPlane & getReferencePlane() const
GFDetPlane getLastPlane() const
void rotationCov(TMatrixT< Double_t > &cov, const TVector3 &u, const TVector3 &v)
WireGeo const & Wire(unsigned int iwire) const
Definition: PlaneGeo.cxx:465
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
Definition: TrackingTypes.h:78
iterator begin()
Definition: PtrVector.h:217
Declaration of signal hit object.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
Length_t DetHalfWidth(TPCID const &tpcid=tpc_zero) const
Returns the half width of the active volume of the specified TPC.
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
void setMomHigh(Double_t f)
Definition: GFKalman.h:118
iterator erase(iterator position)
Definition: PtrVector.h:504
TrackTrajectory::Flags_t Flags_t
Definition: Track.h:67
std::vector< Double_t > getHitChi2()
Definition: GFTrack.h:368
constexpr auto abs(T v)
Returns the absolute value of the argument.
#define MF_LOG_ERROR(category)
std::vector< TMatrixT< Double_t > > getHitCov7x7()
Definition: GFTrack.h:371
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 >> SMatrixSym55
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:244
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:280
Float_t tmp
Definition: plot.C:35
int NParticles() const
Definition: MCTruth.h:75
std::vector< double > fPosErr
void Print(std::ostream &out=std::cout) const
Definition: GFDetPlane.cxx:228
const TMatrixT< Double_t > & getState() const
bool isRealData() const
Definition: Event.cc:53
std::vector< TMatrixT< Double_t > > getHitState()
Definition: GFTrack.h:369
std::vector< TVector3 > getHitPlaneU()
Definition: GFTrack.h:374
static bool sp_sort_3dy(const art::Ptr< recob::SpacePoint > &h1, const art::Ptr< recob::SpacePoint > &h2)
void produce(art::Event &evt)
Length_t DetLength(TPCID const &tpcid=tpc_zero) const
Returns the length of the active volume of the specified TPC.
double ThetaZ() const
Returns angle of wire with respect to z axis in the Y-Z plane in radians.
Definition: WireGeo.h:248
Base Class for genfit track representations. Defines interface for track parameterizations.
Definition: GFAbsTrackRep.h:81
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
PlaneGeo const & Plane(PlaneID const &planeid) const
Returns the specified wire.
TString part[npart]
Definition: Style.C:32
void hits()
Definition: readHits.C:15
typename data_t::const_iterator const_iterator
Definition: PtrVector.h:55
std::vector< double > fMomStart
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
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.
enum geo::_plane_sigtype SigType_t
Enumerate the possible plane projections.
void setBlowUpFactor(double f)
Set the blowup factor (see blowUpCovs() )
Definition: GFKalman.h:116
static bool sp_sort_nsppts(const art::PtrVector< recob::SpacePoint > &h1, const art::PtrVector< recob::SpacePoint > &h2)
iterator end()
Definition: PtrVector.h:231
double energyLossBetheBloch(const double &mass, const double p)
void setNumIterations(Int_t i)
Set number of iterations for Kalman Filter.
Definition: GFKalman.h:91
void setPDG(int pdgt)
Definition: GFTrack.h:365
Float_t d
Definition: plot.C:235
int getFailedHits(int repId=-1)
return the number of failed Hits in track fit repId == -1 will use cardinal rep
Definition: GFTrack.h:147
TMatrixT< Double_t > * covREC
void processTrack(GFTrack *)
Performs fit on a GFTrack.
Definition: GFKalman.cxx:57
Provides recob::Track data product.
static GFFieldManager * getInstance()
TMatrixT< Double_t > * covMCT
bool isDebugEnabled()
const Double32_t * XYZ() const
Definition: SpacePoint.h:78
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
Definition: TrackingTypes.h:68
constexpr PlaneID const & asPlaneID() const
Conversion to PlaneID (for convenience of notation).
Definition: geo_types.h:520
Declaration of cluster object.
const TMatrixT< Double_t > & getCov() const
std::vector< TVector3 > getHitPlaneV()
Definition: GFTrack.h:375
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
Definition: PFPUtils.cxx:2675
std::string ROOTobjectToString(const ROOTOBJ &obj)
Shortcut to write one ROOT object into a string.
Definition: GFException.h:114
size_type size() const
Definition: PtrVector.h:302
TH1F * h2
Definition: plot.C:44
Encapsulate the geometry of a wire.
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:76
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
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)
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Utility object to perform functions of association.
Encapsulate the construction of a single detector plane.
in close()
TDirectory * dir
Definition: macro.C:5
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
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:112
#define MF_LOG_DEBUG(id)
void init(GFAbsBField *b)
set the magntic field here. Magnetic field classes must be derived from GFAbsBField ...
TH1F * h1
Definition: plot.C:41
EventNumber_t event() const
Definition: EventID.h:116
std::vector< TMatrixT< Double_t > > getHitUpdate()
Definition: GFTrack.h:367
Length_t DetHalfHeight(TPCID const &tpcid=tpc_zero) const
Returns the half height of the active volume of the specified TPC.
void setErrorScaleSTh(Double_t f)
Definition: GFKalman.h:120
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
double getChiSqu() const
void addHit(GFAbsRecoHit *theHit)
deprecated!
Definition: GFTrack.h:301
TCEvent evt
Definition: DataStructs.cxx:8
std::vector< TMatrixT< Double_t > > getHitCov()
Definition: GFTrack.h:370
Track3DKalmanSPS(fhicl::ParameterSet const &pset)
#define MF_LOG_WARNING(category)
std::vector< TMatrixT< Double_t > > getHitMeasuredCov()
Definition: GFTrack.h:366
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.cc:29
Float_t w
Definition: plot.C:20
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
EventID id() const
Definition: Event.cc:23
art framework interface to geometry description
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
void setErrorScaleMTh(Double_t f)
Definition: GFKalman.h:121
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
Signal from collection planes.
Definition: geo_types.h:152
std::vector< TVector3 > getHitPlaneXYZ()
Definition: GFTrack.h:372