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