All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Track3DKalmanSPS_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
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.
9 ////////////////////////////////////////////////////////////////////////
10 
11 // C++ includes
12 #include <cmath>
13 #include <vector>
14 #include <string>
15 #include <sstream>
16 #include <iterator> // std::distance()
17 #include <algorithm> // std::sort()
18 
19 // ROOT includes
20 #include "TVector3.h"
21 #include "TMath.h"
22 #include "TPrincipal.h"
23 #include "TDatabasePDG.h"
24 #include "TTree.h"
25 #include "TMatrixT.h"
26 
27 // Framework includes
28 #include "messagefacility/MessageLogger/MessageLogger.h"
29 #include "fhiclcpp/ParameterSet.h"
30 #include "fhiclcpp/exception.h"
31 
32 #include "art/Framework/Principal/Event.h"
33 #include "art/Framework/Principal/Handle.h"
34 #include "canvas/Persistency/Common/Ptr.h"
35 #include "canvas/Persistency/Common/PtrVector.h"
36 #include "art/Framework/Services/Registry/ServiceHandle.h"
37 #include "art_root_io/TFileService.h"
38 #include "art/Framework/Core/ModuleMacros.h"
39 #include "art/Framework/Core/EDProducer.h"
40 #include "canvas/Persistency/Common/FindManyP.h"
41 
42 // GENFIT includes
43 #include "larreco/Genfit/GFAbsTrackRep.h"
44 #include "larreco/Genfit/RKTrackRep.h"
45 #include "larreco/Genfit/GFConstField.h"
46 #include "larreco/Genfit/GFFieldManager.h"
47 #include "larreco/Genfit/PointHit.h"
48 #include "larreco/Genfit/GFTrack.h"
49 #include "larreco/Genfit/GFKalman.h"
50 
51 // LArSoft includes
52 #include "larcore/Geometry/Geometry.h"
53 #include "larcorealg/Geometry/WireGeo.h"
54 #include "larcorealg/Geometry/PlaneGeo.h"
55 #include "lardataobj/RecoBase/Cluster.h"
56 #include "lardataobj/RecoBase/Hit.h"
57 #include "lardataobj/RecoBase/Track.h"
58 #include "lardataobj/RecoBase/SpacePoint.h"
59 //\todo Reconstruction Producers should never include SimulationBase objects
60 #include "nusimdata/SimulationBase/MCTruth.h"
61 #include "lardata/Utilities/AssociationUtil.h"
62 
63 static bool sp_sort_3dz(const art::Ptr<recob::SpacePoint>& h1, const art::Ptr<recob::SpacePoint>& h2)
64 {
65  const double* xyz1 = h1->XYZ();
66  const double* xyz2 = h2->XYZ();
67  return xyz1[2] < xyz2[2];
68 }
69 static bool sp_sort_3dy(const art::Ptr<recob::SpacePoint>& h1, const art::Ptr<recob::SpacePoint>& h2)
70 {
71  const double* xyz1 = h1->XYZ();
72  const double* xyz2 = h2->XYZ();
73  return xyz1[1] < xyz2[1];
74 }
75 static bool sp_sort_3dx(const art::Ptr<recob::SpacePoint>& h1, const art::Ptr<recob::SpacePoint>& h2)
76 {
77  const double* xyz1 = h1->XYZ();
78  const double* xyz2 = h2->XYZ();
79  return xyz1[0] < xyz2[0];
80 }
81 
82 static bool sp_sort_nsppts(const art::PtrVector<recob::SpacePoint>& h1, const art::PtrVector<recob::SpacePoint>& h2)
83 {
84  const unsigned int s1 = h1.size();
85  const unsigned int s2 = h2.size();
86  return s1 > s2;
87 }
88 
89 
90 namespace trkf {
91 
92  class Track3DKalmanSPS : public art::EDProducer {
93  public:
94  explicit Track3DKalmanSPS(fhicl::ParameterSet const& pset);
95 
96  private:
97  void produce(art::Event& evt);
98  void beginJob();
99  void endJob();
100  double energyLossBetheBloch(const double& mass,
101  const double p
102  );
103 
104  void rotationCov(TMatrixT<Double_t> &cov, const TVector3 &u, const TVector3 &v);
105  std::vector <double> dQdxCalc(const art::FindManyP<recob::Hit> &h, const art::PtrVector<recob::SpacePoint> &s, const TVector3 &p, const TVector3 &d );
106 
107  std::string fClusterModuleLabel;// label for input collection
108  std::string fSpptModuleLabel;// label for input collection
109  std::string fGenieGenModuleLabel;// label for input MC single particle generator
110  std::string fG4ModuleLabel;// label for input MC single particle generator
111  std::string fSortDim; // direction in which to sort spacepoints
112 
113  // TFile *fileGENFIT;
114  TTree *tree;
115 
116  TMatrixT<Double_t> *stMCT;
117  TMatrixT<Double_t> *covMCT;
118  TMatrixT<Double_t> *stREC;
119  TMatrixT<Double_t> *covREC;
120  Float_t chi2;
121  Float_t chi2ndf;
122  int fcont;
123 
124  Float_t *fpRECL;
125  Float_t *fpREC;
126  Float_t *fpMCMom;
127  Float_t *fpMCPos;
128  Float_t *fState0;
129  Float_t *fCov0;
130  int nfail;
131  int ndf;
132  int nchi2rePass;
133  int ispptvec;
134  int nspptvec;
135  unsigned int evtt;
136  unsigned int nTrks;
137  unsigned int fptsNo;
138  Float_t *fshx;
139  Float_t *fshy;
140  Float_t *fshz;
141  Float_t *feshx;
142  Float_t *feshy;
143  Float_t *feshz;
144  Float_t *feshyz;
145  Float_t *fupdate;
146  Float_t *fchi2hit;
147  Float_t *fth;
148  Float_t *feth;
149  Float_t *fedudw;
150  Float_t *fedvdw;
151  Float_t *feu;
152  Float_t *fev;
153  Float_t *fsep;
154  Float_t *fdQdx;
155  unsigned int fDimSize; // if necessary will get this from pset in constructor.
156  Float_t *fPCmeans;
157  Float_t *fPCevals;
158  Float_t *fPCsigmas;
159  Float_t *fPC1;
160  Float_t *fPC2;
161  Float_t *fPC3;
162 
163  std::vector<double> fPosErr;
164  std::vector<double> fMomErr;
165  std::vector<double> fMomStart;
166  double fPerpLim;
167  bool fDoFit;
168  int fNumIt;
169  uint16_t fMinNumSppts;
170  double fErrScaleS;
171  double fErrScaleM;
172  int fDecimate;
173  double fMaxUpdate;
174  int fDecimateU;
175  double fDistanceU;
176  double fMaxUpdateU;
177  double fMomLow;
178  double fMomHigh;
179  int fPdg;
180  double fChi2Thresh;
181  int fMaxPass;
182 
183  genf::GFAbsTrackRep *repMC;
184  genf::GFAbsTrackRep *rep;
185 
186  }; // class Track3DKalmanSPS
187 
188 }
189 
190 namespace trkf {
191 
192 //-------------------------------------------------
193  Track3DKalmanSPS::Track3DKalmanSPS(fhicl::ParameterSet const& pset)
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  }
256 
257 
258 //-------------------------------------------------
259 // stolen, mostly, from GFMaterialEffects.
260  double Track3DKalmanSPS::energyLossBetheBloch(const double& mass,
261  const double p=1.5
262  )
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  }
291 
292  void Track3DKalmanSPS::rotationCov(TMatrixT<Double_t> &cov, const TVector3 &u, const TVector3 &v)
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  }
324 
325  std::vector<double> Track3DKalmanSPS::dQdxCalc(const art::FindManyP<recob::Hit> &h, const art::PtrVector<recob::SpacePoint> &s, const TVector3 &dir, const TVector3 &loc )
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  }
394 
395 //-------------------------------------------------
396  void Track3DKalmanSPS::beginJob()
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  }
504 
505 //-------------------------------------------------
506  void Track3DKalmanSPS::endJob()
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  }
550 
551 
552 //------------------------------------------------------------------------------------//
553 void Track3DKalmanSPS::produce(art::Event& evt)
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 }
1264 
1265  DEFINE_ART_MODULE(Track3DKalmanSPS)
1266 
1267 } // end namespace
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::GFTrack::getHitPlaneUxUyUz
std::vector< TVector3 > getHitPlaneUxUyUz()
Definition: GFTrack.h:369
genf::GFKalman::setMaxUpdate
void setMaxUpdate(Double_t f)
Definition: GFKalman.h:117
genf::GFTrack::getHitUpdate
std::vector< TMatrixT< Double_t > > getHitUpdate()
Definition: GFTrack.h:363
recob::Hit::SignalType
geo::SigType_t SignalType() const
Signal type for the plane of the hit.
Definition: Hit.h:231
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
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
genf::GFTrack::getHitState
std::vector< TMatrixT< Double_t > > getHitState()
Definition: GFTrack.h:365
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
recob::Hit::WireID
geo::WireID WireID() const
Definition: Hit.h:233
Hit.h
Declaration of signal hit object.
p
pdgs p
Definition: selectors.fcl:22
genf::GFKalman::setMomHigh
void setMomHigh(Double_t f)
Definition: GFKalman.h:116
recob::Track::Flags_t
TrackTrajectory::Flags_t Flags_t
Definition: bj/lardataobj/RecoBase/Track.h:68
recob::Hit::Integral
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:224
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
genf::GFTrack::getHitPlaneXYZ
std::vector< TVector3 > getHitPlaneXYZ()
Definition: GFTrack.h:368
genf::GFTrack::getHitCov
std::vector< TMatrixT< Double_t > > getHitCov()
Definition: GFTrack.h:366
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
trkf::Track3DKalmanSPS::produce
void produce(art::Event &evt)
Definition: Track3DKalmanSPS_module.cc:553
genf::GFAbsTrackRep
Base Class for genfit track representations. Defines interface for track parameterizations.
Definition: GFAbsTrackRep.h:83
h
while getopts h
Definition: getStatistics.sh:13
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
trkf::Track3DKalmanSPS::fMomStart
std::vector< double > fMomStart
Definition: Track3DKalmanSPS_module.cc:165
icarus::distance
double distance(geo::Point_t const &point, CathodeDesc_t const &cathode)
Returns the distance of a point from the cathode.
Definition: CathodeCrossingUtils.cxx:35
genf::GFTrack::getHitPlaneV
std::vector< TVector3 > getHitPlaneV()
Definition: GFTrack.h:371
genf::GFAbsTrackRep::getMom
virtual TVector3 getMom(const GFDetPlane &pl)=0
genf::GFTrack::getHitChi2
std::vector< Double_t > getHitChi2()
Definition: GFTrack.h:364
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
geo::SigType_t
enum geo::_plane_sigtype SigType_t
genf::GFKalman::setBlowUpFactor
void setBlowUpFactor(double f)
Set the blowup factor (see blowUpCovs() )
Definition: GFKalman.h:114
detsim::MF_LOG_WARNING
MF_LOG_WARNING("SimWire")<< "SimWire is an example module that works for the "<< "MicroBooNE detector. Each experiment should implement "<< "its own version of this module to simulate electronics "<< "response."
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
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
genf::GFTrack::setPDG
void setPDG(int pdgt)
Definition: GFTrack.h:361
genf::GFTrack::getFailedHits
int getFailedHits(int repId=-1)
return the number of failed Hits in track fit repId == -1 will use cardinal rep
Definition: GFTrack.h:149
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
genf::GFFieldManager::getInstance
static GFFieldManager * getInstance()
Definition: GFFieldManager.h:70
trkf::Track3DKalmanSPS::covMCT
TMatrixT< Double_t > * covMCT
Definition: Track3DKalmanSPS_module.cc:117
recob::tracking::convertCollToPoint
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
Definition: TrackingTypes.h:55
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
Cluster.h
Declaration of cluster object.
genf::GFAbsTrackRep::getCov
const TMatrixT< Double_t > & getCov() const
Definition: GFAbsTrackRep.h:197
genf::GFTrack::getHitPlaneU
std::vector< TVector3 > getHitPlaneU()
Definition: GFTrack.h:370
tca::dEdx
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
Definition: PFPUtils.cxx:2687
Track.h
Provides recob::Track data product.
genf::ROOTobjectToString
std::string ROOTobjectToString(const ROOTOBJ &obj)
Shortcut to write one ROOT object into a string.
Definition: GFException.h:116
gen_crt_frags.w
tuple w
Definition: gen_crt_frags.py:1416
make_calib_df.exp
def exp
Definition: make_calib_df.py:43
dropbox.dir
tuple dir
Definition: dropbox.py:28
WireGeo.h
Encapsulate the geometry of a wire.
genf::GFTrack::getHitCov7x7
std::vector< TMatrixT< Double_t > > getHitCov7x7()
Definition: GFTrack.h:367
one
then echo Cowardly refusing to create a new FHiCL file with the same name as the original one('${SourceName}')." >&2 exit 1 fi echo "'$
Definition: makeNoOverburdenVersionOf.sh:65
util::CreateAssn
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.
Definition: lardata/lardata/Utilities/AssociationUtil.h:586
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
PlaneGeo.h
Encapsulate the construction of a single detector plane.
larg4::dist
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Definition: OpFastScintillation.hh:563
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
s
then echo File list $list not found else cat $list while read file do echo $file sed s
Definition: file_to_url.sh:60
lib.dedx.mass
float mass
Definition: dedx.py:47
e
do i e
Definition: updateVersionList.sh:68
genf::GFFieldManager::init
void init(GFAbsBField *b)
set the magntic field here. Magnetic field classes must be derived from GFAbsBField ...
Definition: GFFieldManager.h:66
angle
finds tracks best matching by angle
Definition: trackfinderalgorithms.fcl:64
genf::GFKalman::setErrorScaleSTh
void setErrorScaleSTh(Double_t f)
Definition: GFKalman.h:118
recob::Hit
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
genf::GFAbsTrackRep::getChiSqu
double getChiSqu() const
Definition: GFAbsTrackRep.h:243
genf::GFTrack::addHit
void addHit(GFAbsRecoHit *theHit)
deprecated!
Definition: GFTrack.h:300
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
k
pdgs k
Definition: selectors.fcl:22
trkf::Track3DKalmanSPS::Track3DKalmanSPS
Track3DKalmanSPS(fhicl::ParameterSet const &pset)
Definition: Track3DKalmanSPS_module.cc:193
close
print OUTPUT<< EOF;< setup name="Default"version="1.0">< worldref="volWorld"/></setup ></gdml > EOF close(OUTPUT)
Geometry.h
art framework interface to geometry description
recob::Track
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a &quot;fitted&quot; track:
Definition: bj/lardataobj/RecoBase/Track.h:49
genf::GFKalman::setErrorScaleMTh
void setErrorScaleMTh(Double_t f)
Definition: GFKalman.h:119
genf::GFTrack::getHitMeasuredCov
std::vector< TMatrixT< Double_t > > getHitMeasuredCov()
Definition: GFTrack.h:362
geo::kCollection
Signal from collection planes.
Definition: geo_types.h:146
Generated by   doxygen 1.8.5