LArSoft  v08_35_01
Liquid Argon Software toolkit - http://larsoft.org/
KalmanFilterFinalTrackFitter_module.cc
Go to the documentation of this file.
1 
12 
13 #include "fhiclcpp/types/Atom.h"
14 #include "fhiclcpp/types/Table.h"
17 
27 
29 
35 
36 #include <memory>
37 
38 namespace trkf {
39 
41  public:
42 
43  struct Inputs {
44  using Name = fhicl::Name;
47  Name("inputPFParticleLabel"),
48  Comment("Label of recob::PFParticle Collection to be fit")
49  };
51  Name("inputTracksLabel"),
52  Comment("Label of recob::Track Collection to be fit")
53  };
55  Name("inputShowersLabel"),
56  Comment("Label of recob::Shower Collection (associated to PFParticles) to be fit")
57  };
59  Name("inputCaloLabel"),
60  Comment("Label of anab::Calorimetry Collection, matching inputTracksLabel, to be used for initial momentum estimate. Used only if momFromCalo is set to true.")
61  };
63  Name("inputMCLabel"),
64  Comment("Label of sim::MCTrack Collection to be used for initial momentum estimate. Used only if momFromMC is set to true.")
65  };
67  Name("inputPidLabel"),
68  Comment("Label of anab::ParticleID Collection, matching inputTracksLabel, to be used for particle Id.")
69  };
70  };
71 
72  struct Options {
73  using Name = fhicl::Name;
75  fhicl::Atom<bool> trackFromPF {
76  Name("trackFromPF"),
77  Comment("If true extract tracks from inputPFParticleLabel collection, if false from inputTracksLabel.")
78  };
79  fhicl::Atom<bool> showerFromPF {
80  Name("showerFromPF"),
81  Comment("If true extract showers from inputPFParticleLabel collection.")
82  };
83  fhicl::Atom<bool> pFromMSChi2 {
84  Name("momFromMSChi2"),
85  Comment("Flag used to get initial momentum estimate from trkf::TrackMomentumCalculator::GetMomentumMultiScatterChi2().")
86  };
87  fhicl::Atom<bool> pFromLength {
88  Name("momFromLength"),
89  Comment("Flag used to get initial momentum estimate from trkf::TrackMomentumCalculator::GetTrackMomentum().")
90  };
91  fhicl::Atom<bool> pFromCalo {
92  Name("momFromCalo"),
93  Comment("Flag used to get initial momentum estimate from inputCaloLabel collection.")
94  };
96  Name("momFromMC"),
97  Comment("Flag used to get initial momentum estimate from inputMCLabel collection.")
98  };
100  Name("momentumInGeV"),
101  Comment("Fixed momentum estimate value, to be used when momFromCalo, momFromMSChi2, momFromLength and momFromMC are all false, or if the estimate is not available.")
102  };
103  fhicl::Atom<bool> idFromPF {
104  Name("idFromPF"),
105  Comment("Flag used to get particle ID estimate from corresponding PFParticle. Needs trackFromPF=true.")
106  };
107  fhicl::Atom<bool> idFromCollection {
108  Name("idFromCollection"),
109  Comment("Flag used to get particle ID estimate from inputPidLabel collection.")
110  };
112  Name("pdgId"),
113  Comment("Default particle id hypothesis in case no valid id is provided either via PFParticle or in the ParticleId collection.")
114  };
115  fhicl::Atom<bool> dirFromVtxPF {
116  Name("dirFromVtxPF"),
117  Comment("Assume track direction from Vertex in PFParticle. Needs trackFromPF=true.")
118  };
119  fhicl::Atom<bool> dirFromMC {
120  Name("dirFromMC"),
121  Comment("Assume track direction from MC.")
122  };
123  fhicl::Atom<bool> dirFromVec {
124  Name("dirFromVec"),
125  Comment("Assume track direction from as the one giving positive dot product with vector specified by dirVec.")
126  };
128  Name("dirVec"),
129  Comment("Fhicl sequence defining the vector used when dirFromVec=true. It must have 3 elements.")
130  };
131  fhicl::Atom<bool> alwaysInvertDir {
132  Name("alwaysInvertDir"),
133  Comment("If true, fit all tracks from end to vertex assuming inverted direction.")
134  };
135  fhicl::Atom<bool> produceTrackFitHitInfo {
136  Name("produceTrackFitHitInfo"),
137  Comment("Option to produce (or not) the detailed TrackFitHitInfo.")
138  };
139  fhicl::Atom<bool> produceSpacePoints {
140  Name("produceSpacePoints"),
141  Comment("Option to produce (or not) the associated SpacePoints.")
142  };
143  fhicl::Atom<bool> keepInputTrajectoryPoints {
144  Name("keepInputTrajectoryPoints"),
145  Comment("Option to keep positions and directions from input trajectory/track. The fit will provide only covariance matrices, chi2, ndof, particle Id and absolute momentum. It may also modify the trajectory point flags. In order to avoid inconsistencies, it has to be used with the following fitter options all set to false: sortHitsByPlane, sortOutputHitsMinLength, skipNegProp.")
146  };
147  };
148 
149  struct Config {
150  using Name = fhicl::Name;
152  Name("inputs"),
153  };
155  Name("options")
156  };
158  Name("propagator")
159  };
161  Name("fitter")
162  };
163  };
165 
166  explicit KalmanFilterFinalTrackFitter(Parameters const & p);
167 
168  // Plugins should not be copied or assigned.
173 
174  private:
175  void produce(art::Event & e) override;
176 
182 
189 
190  std::unique_ptr<art::FindManyP<anab::Calorimetry> > trackCalo;
191  std::unique_ptr<art::FindManyP<anab::ParticleID> > trackId;
192  std::unique_ptr<art::FindManyP<recob::Track> > assocTracks;
193  std::unique_ptr<art::FindManyP<recob::Shower> > assocShowers;
194  std::unique_ptr<art::FindManyP<recob::Vertex> > assocVertices;
195 
196  double setMomValue(art::Ptr<recob::Track> ptrack, const std::unique_ptr<art::FindManyP<anab::Calorimetry> >& trackCalo, const double pMC, const int pId) const;
197  int setPId(const unsigned int iTrack, const std::unique_ptr<art::FindManyP<anab::ParticleID> >& trackId, const int pfPid = 0) const;
198  bool setDirFlip(const recob::Track& track, TVector3& mcdir, const std::vector<art::Ptr<recob::Vertex> >* vertices = 0) const;
199 
201  };
202 }
203 
205  : EDProducer{p}
206  , p_(p)
207  , prop{p_().propagator}
208  , kalmanFitter{&prop, p_().fitter}
209  , inputFromPF{p_().options().trackFromPF() || p_().options().showerFromPF()}
210 {
211 
212  if (inputFromPF) {
213  pfParticleInputTag = art::InputTag(p_().inputs().inputPFParticleLabel());
214  if (p_().options().showerFromPF()) showerInputTag = art::InputTag(p_().inputs().inputShowersLabel());
215  } else {
216  trackInputTag = art::InputTag(p_().inputs().inputTracksLabel());
217  if (p_().options().idFromCollection()) pidInputTag = art::InputTag(p_().inputs().inputPidLabel());
218  }
219  if (p_().options().pFromCalo()) caloInputTag = art::InputTag(p_().inputs().inputCaloLabel());
220  if (p_().options().pFromMC() || p_().options().dirFromMC()) simTrackInputTag = art::InputTag(p_().inputs().inputMCLabel());
221 
222  produces<std::vector<recob::Track> >();
223  produces<art::Assns<recob::Track, recob::Hit, recob::TrackHitMeta> >();
224  produces<art::Assns<recob::Track, recob::Hit> >();
225  if (inputFromPF) {
226  produces<art::Assns<recob::PFParticle, recob::Track> >();
227  }
228  if (p_().options().produceTrackFitHitInfo()) {
229  produces<std::vector<std::vector<recob::TrackFitHitInfo> > >();
230  }
231  if (p_().options().produceSpacePoints()) {
232  produces<std::vector<recob::SpacePoint> >();
233  produces<art::Assns<recob::Hit, recob::SpacePoint> >();
234  }
235 
236  //throw expections to avoid possible silent failures due to incompatible configuration options
237  if (p_().options().trackFromPF()==0 && p_().options().idFromPF())
238  throw cet::exception("KalmanFilterFinalTrackFitter") << "Incompatible configuration parameters: cannot use idFromPF=true with trackFromPF=false." << "\n";
239  if (p_().options().trackFromPF()==0 && p_().options().dirFromVtxPF())
240  throw cet::exception("KalmanFilterFinalTrackFitter") << "Incompatible configuration parameters: cannot use dirFromVtxPF=true with trackFromPF=false." << "\n";
241 
242  unsigned int nIds = 0;
243  if (p_().options().idFromPF()) nIds++;
244  if (p_().options().idFromCollection()) nIds++;
245  if (nIds>1) {
246  throw cet::exception("KalmanFilterFinalTrackFitter")
247  << "Incompatible configuration parameters: only at most one can be set to true among idFromPF and idFromCollection." << "\n";
248  }
249 
250  unsigned int nDirs = 0;
251  if (p_().options().dirFromVtxPF()) nDirs++;
252  if (p_().options().dirFromMC()) nDirs++;
253  if (p_().options().dirFromVec()) nDirs++;
254  if (p_().options().alwaysInvertDir()) nDirs++;
255  if (nDirs>1) {
256  throw cet::exception("KalmanFilterFinalTrackFitter")
257  << "Incompatible configuration parameters: only at most one can be set to true among dirFromVtxPF, dirFromMC, dirFromVec, and alwaysInvertDir." << "\n";
258  }
259 
260  unsigned int nPFroms = 0;
261  if (p_().options().pFromCalo()) nPFroms++;
262  if (p_().options().pFromMSChi2()) nPFroms++;
263  if (p_().options().pFromLength()) nPFroms++;
264  if (p_().options().pFromMC()) nPFroms++;
265  if (nPFroms>1) {
266  throw cet::exception("KalmanFilterFinalTrackFitter")
267  << "Incompatible configuration parameters: only at most one can be set to true among pFromCalo, pFromMSChi2, pFromLength, and pFromMC." << "\n";
268  }
269 
270  if (p_().options().keepInputTrajectoryPoints()) {
271  if (p_().fitter().sortHitsByPlane() || p_().fitter().sortOutputHitsMinLength() || p_().fitter().skipNegProp()) {
272  throw cet::exception("KalmanFilterTrajectoryFitter")
273  << "Incompatible configuration parameters: keepInputTrajectoryPoints needs the following fitter options all set to false: sortHitsByPlane, sortOutputHitsMinLength, skipNegProp." << "\n";
274  }
275  }
276 
277  if (p_().options().showerFromPF()) {
278  if (nPFroms>0 || nIds>0 || nDirs>0) {
279  throw cet::exception("KalmanFilterTrajectoryFitter")
280  << "Incompatible configuration parameters: showerFromPF currently does not support optional momentum values, particle hypotheses and directions." << "\n";
281  }
282  }
283 }
284 
286 {
287 
288  auto outputTracks = std::make_unique<std::vector<recob::Track> >();
289  auto outputHitsMeta = std::make_unique<art::Assns<recob::Track, recob::Hit, recob::TrackHitMeta> >();
290  auto outputHits = std::make_unique<art::Assns<recob::Track, recob::Hit> >();
291  auto outputHitInfo = std::make_unique<std::vector<std::vector<recob::TrackFitHitInfo> > >();
292 
293  auto const tid = e.getProductID<std::vector<recob::Track> >();
294  auto const tidgetter = e.productGetter(tid);
295 
296  auto outputSpacePoints = std::make_unique<std::vector<recob::SpacePoint> >();
297  auto outputHitSpacePointAssn = std::make_unique<art::Assns<recob::Hit, recob::SpacePoint> >();
298  auto const spid = e.getProductID<std::vector<recob::SpacePoint> >();
299  auto const spidgetter = e.productGetter(spid);
300 
301  //FIXME, eventually remove this (ok only for single particle MC)
302  double pMC = -1.;
303  TVector3 mcdir;
304  if (p_().options().pFromMC() || p_().options().dirFromMC()) {
305  art::ValidHandle<std::vector<sim::MCTrack> > simTracks = e.getValidHandle<std::vector<sim::MCTrack> >(simTrackInputTag);
306  for (unsigned int iMC = 0; iMC < simTracks->size(); ++iMC) {
307  const sim::MCTrack& mctrack = simTracks->at(iMC);
308  //fiducial cuts on MC tracks
309  if (mctrack.PdgCode()!=13) continue;
310  if (mctrack.Process()!="primary") continue;
311  pMC = mctrack.Start().Momentum().P()*0.001;
312  mcdir = TVector3(mctrack.Start().Momentum().X()*0.001/pMC,mctrack.Start().Momentum().Y()*0.001/pMC,mctrack.Start().Momentum().Z()*0.001/pMC);
313  break;
314  }
315  //std::cout << "mc momentum value = " << pval << " GeV" << std::endl;
316  }
317 
318  if (inputFromPF) {
319 
320  auto outputPFAssn = std::make_unique<art::Assns<recob::PFParticle, recob::Track> >();
321 
322  auto inputPFParticle = e.getValidHandle<std::vector<recob::PFParticle> >(pfParticleInputTag);
323  if (p_().options().trackFromPF()) assocTracks = std::make_unique<art::FindManyP<recob::Track>>(inputPFParticle, e, pfParticleInputTag);
324  if (p_().options().showerFromPF()) assocShowers = std::make_unique<art::FindManyP<recob::Shower>>(inputPFParticle, e, showerInputTag);
325  assocVertices = std::make_unique<art::FindManyP<recob::Vertex>>(inputPFParticle, e, pfParticleInputTag);
326 
327  for (unsigned int iPF = 0; iPF < inputPFParticle->size(); ++iPF) {
328 
329  if (p_().options().trackFromPF()) {
330  const std::vector<art::Ptr<recob::Track> >& tracks = assocTracks->at(iPF);
331  auto const& tkHitsAssn = *e.getValidHandle<art::Assns<recob::Track, recob::Hit> >(pfParticleInputTag);
332  const std::vector<art::Ptr<recob::Vertex> >& vertices = assocVertices->at(iPF);
333 
334  if (p_().options().pFromCalo()) {
335  trackCalo = std::make_unique<art::FindManyP<anab::Calorimetry>>(tracks, e, caloInputTag);
336  }
337 
338  for (unsigned int iTrack = 0; iTrack < tracks.size(); ++iTrack) {
339 
340  const recob::Track& track = *tracks[iTrack];
341  art::Ptr<recob::Track> ptrack = tracks[iTrack];
342  const int pId = setPId(iTrack, trackId, inputPFParticle->at(iPF).PdgCode());
343  const double mom = setMomValue(ptrack, trackCalo, pMC, pId);
344  const bool flipDir = setDirFlip(track, mcdir, &vertices);
345 
346  //this is not computationally optimal, but at least preserves the order unlike FindManyP
347  std::vector<art::Ptr<recob::Hit> > inHits;
348  for (auto it = tkHitsAssn.begin(); it!=tkHitsAssn.end(); ++it) {
349  if (it->first == ptrack) inHits.push_back(it->second);
350  else if (inHits.size()>0) break;
351  }
352 
353  recob::Track outTrack;
354  std::vector<art::Ptr<recob::Hit> > outHits;
355  trkmkr::OptionalOutputs optionals;
356  if (p_().options().produceTrackFitHitInfo()) optionals.initTrackFitInfos();
357  bool fitok = kalmanFitter.fitTrack(track.Trajectory(),track.ID(),
359  inHits, mom, pId, flipDir, outTrack, outHits, optionals);
360  if (!fitok) continue;
361 
362  if (p_().options().keepInputTrajectoryPoints()) {
363  restoreInputPoints(track.Trajectory().Trajectory(),inHits,outTrack,outHits);
364  }
365 
366  outputTracks->emplace_back(std::move(outTrack));
367  art::Ptr<recob::Track> aptr(tid, outputTracks->size()-1, tidgetter);
368  unsigned int ip = 0;
369  for (auto const& trhit: outHits) {
370  //the fitter produces collections with 1-1 match between hits and point
371  recob::TrackHitMeta metadata(ip,-1);
372  outputHitsMeta->addSingle(aptr, trhit, metadata);
373  outputHits->addSingle(aptr, trhit);
374  ip++;
375  }
376  outputPFAssn->addSingle(art::Ptr<recob::PFParticle>(inputPFParticle, iPF), aptr);
377  outputHitInfo->emplace_back(optionals.trackFitHitInfos());
378  }
379  }
380 
381  if (p_().options().showerFromPF()) {
382  art::Ptr<recob::PFParticle> pPF(inputPFParticle, iPF);
383  const std::vector<art::Ptr<recob::Shower> >& showers = assocShowers->at(iPF);
384  if (showers.size()==0) continue;
385  auto const& pfClustersAssn = *e.getValidHandle<art::Assns<recob::PFParticle, recob::Cluster> >(showerInputTag);
386  auto const& clHitsAssn = *e.getValidHandle<art::Assns<recob::Cluster, recob::Hit> >(showerInputTag);
387  std::vector<art::Ptr<recob::Hit> > inHits;
388  for (auto itpf = pfClustersAssn.begin(); itpf!=pfClustersAssn.end(); ++itpf) {
389  if (itpf->first == pPF) {
390  art::Ptr<recob::Cluster> clust = itpf->second;
391  for (auto it = clHitsAssn.begin(); it!=clHitsAssn.end(); ++it) {
392  if (it->first == clust) inHits.push_back(it->second);
393  }
394  } else if (inHits.size()>0) break;
395  }
396  //auto const& shHitsAssn = *e.getValidHandle<art::Assns<recob::Shower, recob::Hit> >(showerInputTag);
397  for (unsigned int iShower = 0; iShower < showers.size(); ++iShower) {
398  //
399  const recob::Shower& shower = *showers[iShower];
400  // art::Ptr<recob::Shower> pshower = showers[iShower];
401  // //this is not computationally optimal, but at least preserves the order unlike FindManyP
402  // std::vector<art::Ptr<recob::Hit> > inHits;
403  // for (auto it = shHitsAssn.begin(); it!=shHitsAssn.end(); ++it) {
404  // if (it->first == pshower) inHits.push_back(it->second);
405  // else if (inHits.size()>0) break;
406  // }
407 
408  recob::Track outTrack;
409  std::vector<art::Ptr<recob::Hit> > outHits;
410  trkmkr::OptionalOutputs optionals;
411  if (p_().options().produceTrackFitHitInfo()) optionals.initTrackFitInfos();
412  Point_t pos(shower.ShowerStart().X(),shower.ShowerStart().Y(),shower.ShowerStart().Z());
413  Vector_t dir(shower.Direction().X(),shower.Direction().Y(),shower.Direction().Z());
414  auto cov = SMatrixSym55();
415  auto pid = p_().options().pdgId();
416  auto mom = p_().options().pval();
417  bool fitok = kalmanFitter.fitTrack(pos, dir, cov, inHits, std::vector<recob::TrajectoryPointFlags>(),
418  shower.ID(), mom, pid,
419  outTrack, outHits, optionals);
420  if (!fitok) continue;
421 
422  outputTracks->emplace_back(std::move(outTrack));
423  art::Ptr<recob::Track> aptr(tid, outputTracks->size()-1, tidgetter);
424  unsigned int ip = 0;
425  for (auto const& trhit: outHits) {
426  // the fitter produces collections with 1-1 match between hits and point
427  recob::TrackHitMeta metadata(ip,-1);
428  outputHitsMeta->addSingle(aptr, trhit, metadata);
429  outputHits->addSingle(aptr, trhit);
430  if (p_().options().produceSpacePoints() && outputTracks->back().HasValidPoint(ip)) {
431  auto& tp = outputTracks->back().Trajectory().LocationAtPoint(ip);
432  double fXYZ[3] = {tp.X(),tp.Y(),tp.Z()};
433  double fErrXYZ[6] = {0};
434  recob::SpacePoint sp(fXYZ, fErrXYZ, -1.);
435  outputSpacePoints->emplace_back(std::move(sp));
436  art::Ptr<recob::SpacePoint> apsp(spid, outputSpacePoints->size()-1, spidgetter);
437  outputHitSpacePointAssn->addSingle(trhit, apsp);
438  }
439  ip++;
440  }
441  outputPFAssn->addSingle(art::Ptr<recob::PFParticle>(inputPFParticle, iPF), aptr);
442  outputHitInfo->emplace_back(optionals.trackFitHitInfos());
443  }
444  }
445 
446  }
447  e.put(std::move(outputTracks));
448  e.put(std::move(outputHitsMeta));
449  e.put(std::move(outputHits));
450  e.put(std::move(outputPFAssn));
451  if (p_().options().produceTrackFitHitInfo()) {
452  e.put(std::move(outputHitInfo));
453  }
454  if (p_().options().produceSpacePoints()) {
455  e.put(std::move(outputSpacePoints));
456  e.put(std::move(outputHitSpacePointAssn));
457  }
458  } else {
459 
460  art::ValidHandle<std::vector<recob::Track> > inputTracks = e.getValidHandle<std::vector<recob::Track> >(trackInputTag);
462 
463  if (p_().options().pFromCalo()) {
464  trackCalo = std::make_unique<art::FindManyP<anab::Calorimetry>>(inputTracks, e, caloInputTag);
465  }
466 
467  if (p_().options().idFromCollection()) {
468  trackId = std::make_unique<art::FindManyP<anab::ParticleID>>(inputTracks, e, pidInputTag);
469  }
470 
471  for (unsigned int iTrack = 0; iTrack < inputTracks->size(); ++iTrack) {
472 
473  const recob::Track& track = inputTracks->at(iTrack);
474  art::Ptr<recob::Track> ptrack(inputTracks, iTrack);
475  const int pId = setPId(iTrack, trackId);
476  const double mom = setMomValue(ptrack, trackCalo, pMC, pId);
477  const bool flipDir = setDirFlip(track, mcdir);
478 
479  //this is not computationally optimal, but at least preserves the order unlike FindManyP
480  std::vector<art::Ptr<recob::Hit> > inHits;
481  for (auto it = tkHitsAssn.begin(); it!=tkHitsAssn.end(); ++it) {
482  if (it->first == ptrack) inHits.push_back(it->second);
483  else if (inHits.size()>0) break;
484  }
485 
486  recob::Track outTrack;
487  std::vector<art::Ptr<recob::Hit> > outHits;
488  trkmkr::OptionalOutputs optionals;
489  if (p_().options().produceTrackFitHitInfo()) optionals.initTrackFitInfos();
490  bool fitok = kalmanFitter.fitTrack(track.Trajectory(),track.ID(),
492  inHits, mom, pId, flipDir, outTrack, outHits, optionals);
493  if (!fitok) continue;
494 
495  if (p_().options().keepInputTrajectoryPoints()) {
496  restoreInputPoints(track.Trajectory().Trajectory(),inHits,outTrack,outHits);
497  }
498 
499  outputTracks->emplace_back(std::move(outTrack));
500  art::Ptr<recob::Track> aptr(tid, outputTracks->size()-1, tidgetter);
501  unsigned int ip = 0;
502  for (auto const& trhit: outHits) {
503  //the fitter produces collections with 1-1 match between hits and point
504  recob::TrackHitMeta metadata(ip,-1);
505  outputHitsMeta->addSingle(aptr, trhit, metadata);
506  outputHits->addSingle(aptr, trhit);
507  if (p_().options().produceSpacePoints() && outputTracks->back().HasValidPoint(ip)) {
508  auto& tp = outputTracks->back().Trajectory().LocationAtPoint(ip);
509  double fXYZ[3] = {tp.X(),tp.Y(),tp.Z()};
510  double fErrXYZ[6] = {0};
511  recob::SpacePoint sp(fXYZ, fErrXYZ, -1.);
512  outputSpacePoints->emplace_back(std::move(sp));
513  art::Ptr<recob::SpacePoint> apsp(spid, outputSpacePoints->size()-1, spidgetter);
514  outputHitSpacePointAssn->addSingle(trhit, apsp);
515  }
516  ip++;
517  }
518  outputHitInfo->emplace_back(optionals.trackFitHitInfos());
519  }
520  e.put(std::move(outputTracks));
521  e.put(std::move(outputHitsMeta));
522  e.put(std::move(outputHits));
523  if (p_().options().produceTrackFitHitInfo()) {
524  e.put(std::move(outputHitInfo));
525  }
526  if (p_().options().produceSpacePoints()) {
527  e.put(std::move(outputSpacePoints));
528  e.put(std::move(outputHitSpacePointAssn));
529  }
530  }
531 }
532 
534  const auto np = outTrack.NumberTrajectoryPoints();
535  std::vector<Point_t> positions(np);
536  std::vector<Vector_t> momenta(np);
537  std::vector<recob::TrajectoryPointFlags> outFlags(np);
538  //
539  for (unsigned int p=0; p<np; ++p) {
540  auto flag = outTrack.FlagsAtPoint(p);
541  auto mom = outTrack.VertexMomentum();
542  auto op = flag.fromHit();
543  positions[op] = track.LocationAtPoint(op);
544  momenta[op] = mom*track.DirectionAtPoint(op);
545  auto mask = flag.mask();
547  outFlags[op] = recob::TrajectoryPointFlags(op,mask);
548  }
549  auto covs = outTrack.Covariances();
550  outTrack = recob::Track(recob::TrackTrajectory(std::move(positions),std::move(momenta),std::move(outFlags),true),
551  outTrack.ParticleId(),outTrack.Chi2(),outTrack.Ndof(),std::move(covs.first),std::move(covs.second),outTrack.ID());
552  //
553  outHits.clear();
554  for (auto h : inHits) outHits.push_back(h);
555 }
556 
557 double trkf::KalmanFilterFinalTrackFitter::setMomValue(art::Ptr<recob::Track> ptrack, const std::unique_ptr<art::FindManyP<anab::Calorimetry> >& trackCalo, const double pMC, const int pId) const {
558  double result = p_().options().pval();
559  if (p_().options().pFromMSChi2()) {
560  result = tmc.GetMomentumMultiScatterChi2(ptrack);
561  } else if (p_().options().pFromLength()) {
562  result = tmc.GetTrackMomentum(ptrack->Length(), pId);
563  } else if (p_().options().pFromCalo()) {
564  //take average energy from available views
565  const std::vector<art::Ptr<anab::Calorimetry> >& calo = trackCalo->at(ptrack.key());
566  double sumenergy = 0.;
567  int nviews = 0.;
568  for (auto caloit : calo) {
569  if (caloit->KineticEnergy()>0.) {
570  sumenergy+=caloit->KineticEnergy();
571  nviews+=1;
572  }
573  }
574  if (nviews!=0 && sumenergy!=0.) {
575  //protect against problematic cases
576  result = sumenergy/(nviews*1000.);
577  }
578  } else if (p_().options().pFromMC() && pMC>0.) {
579  result = pMC;
580  }
581  return result;
582 }
583 
584 int trkf::KalmanFilterFinalTrackFitter::setPId(const unsigned int iTrack, const std::unique_ptr<art::FindManyP<anab::ParticleID> >& trackId, const int pfPid) const {
585  int result = p_().options().pdgId();
586  if (p_().options().trackFromPF() && p_().options().idFromPF()) {
587  result = pfPid;
588  } else if (p_().options().idFromCollection()) {
589  //take the pdgId corresponding to the minimum chi2 (should we give preference to the majority? fixme)
590  double minChi2 = -1.;
591  for (auto idit : trackId->at(iTrack)) {
592  if ( idit->MinChi2()>0. && (minChi2<0. || idit->MinChi2()<minChi2) ) {
593  result = idit->Pdg();
594  minChi2 = idit->MinChi2();
595  }
596  }
597  }
598  return result;
599 }
600 
602  bool result = false;
603  if (p_().options().alwaysInvertDir()) {
604  return true;
605  } else if (p_().options().dirFromMC()) {
606  auto tdir = track.VertexDirection();
607  if ( (mcdir.X()*tdir.X() + mcdir.Y()*tdir.Y() + mcdir.Z()*tdir.Z())<0. ) result = true;
608  } else if (p_().options().dirFromVec()) {
609  std::array<float, 3> dir = p_().options().dirVec();
610  auto tdir = track.VertexDirection();
611  if ( (dir[0]*tdir.X() + dir[1]*tdir.Y() + dir[2]*tdir.Z())<0. ) result = true;
612  } else if (p_().options().trackFromPF() && p_().options().dirFromVtxPF() && vertices->size()>0) {
613  //if track end is closer to first vertex then track vertex, flip direction
614  double xyz[3];
615  (*vertices)[0]->XYZ(xyz);
616  auto& tv = track.Trajectory().Vertex();
617  auto& te = track.Trajectory().End();
618  if ( ((xyz[0]-te.X())*(xyz[0]-te.X()) + (xyz[1]-te.Y())*(xyz[1]-te.Y()) + (xyz[2]-te.Z())*(xyz[2]-te.Z())) >
619  ((xyz[0]-tv.X())*(xyz[0]-tv.X()) + (xyz[1]-tv.Y())*(xyz[1]-tv.Y()) + (xyz[2]-tv.Z())*(xyz[2]-tv.Z())) ) result = true;
620  }
621  return result;
622 }
623 
size_type size() const
Definition: Assns.h:447
Provides recob::Track data product.
const TVector3 & ShowerStart() const
Definition: Shower.h:192
double VertexMomentum() const
Definition: Track.h:142
Trajectory_t const & Trajectory() const
Returns the plain trajectory of this object.
Fit tracks using Kalman Filter fit+smooth.
void initTrackFitInfos()
initialize the output vector of TrackFitHitInfos
Definition: TrackMaker.h:102
ProductID getProductID(std::string const &instance_name="") const
Definition: DataViewImpl.h:349
double GetMomentumMultiScatterChi2(art::Ptr< recob::Track > const &trk)
static constexpr Flag_t NoPoint
The trajectory point is not defined.
int ParticleId() const
Access to various track properties.
Definition: Track.h:171
recob::tracking::Point_t Point_t
Definition: TrackState.h:18
recob::tracking::Vector_t Vector_t
Definition: TrackState.h:19
const recob::TrackTrajectory & Trajectory() const
Access to the stored recob::TrackTrajectory.
Definition: Track.h:98
Declaration of signal hit object.
EDProducer()=default
recob::tracking::SMatrixSym55 SMatrixSym55
Definition: TrackState.h:15
float Chi2() const
Access to various track properties.
Definition: Track.h:168
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:102
Class to keep data related to recob::Hit associated with recob::Track.
Vector_t VertexDirection() const
Access to track direction at different points.
Definition: Track.h:132
Class for propagation of a trkf::TrackState to a recob::tracking::Plane.
Data related to recob::Hit associated with recob::Track.The purpose is to collect several variables t...
Definition: TrackHitMeta.h:43
void restoreInputPoints(const recob::Trajectory &track, const std::vector< art::Ptr< recob::Hit > > &inHits, recob::Track &outTrack, std::vector< art::Ptr< recob::Hit > > &outHits) const
bool fitTrack(const recob::TrackTrajectory &traj, int tkID, const SMatrixSym55 &covVtx, const SMatrixSym55 &covEnd, const std::vector< art::Ptr< recob::Hit > > &hits, const double pval, const int pdgid, const bool flipDirection, recob::Track &outTrack, std::vector< art::Ptr< recob::Hit > > &outHits, trkmkr::OptionalOutputs &optionals) const
Fit track starting from TrackTrajectory.
std::pair< SMatrixSym55, SMatrixSym55 > Covariances() const
Access to position, momentum or covariance at the start and end of the track.
Definition: Track.h:162
std::unique_ptr< art::FindManyP< recob::Vertex > > assocVertices
const_iterator begin() const
Definition: Assns.h:454
Vector_t DirectionAtPoint(size_t i) const
Computes and returns the direction of the trajectory at a point.
Definition: Trajectory.cxx:117
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
bool setDirFlip(const recob::Track &track, TVector3 &mcdir, const std::vector< art::Ptr< recob::Vertex > > *vertices=0) const
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:167
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:56
A trajectory in space reconstructed from hits.
EDProductGetter const * productGetter(ProductID const pid) const
double setMomValue(art::Ptr< recob::Track > ptrack, const std::unique_ptr< art::FindManyP< anab::Calorimetry > > &trackCalo, const double pMC, const int pId) const
key_type key() const
Definition: Ptr.h:238
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:482
Point_t const & LocationAtPoint(size_t i) const
Returns the position at the specified trajectory point.
Definition: Trajectory.h:236
std::unique_ptr< art::FindManyP< anab::ParticleID > > trackId
const TVector3 & Direction() const
Definition: Shower.h:189
Declaration of cluster object.
Class def header for mctrack data container.
int Ndof() const
Access to various track properties.
Definition: Track.h:170
Point_t const & Vertex() const
Returns the position of the first valid point of the trajectory [cm].
int PdgCode() const
Definition: MCTrack.h:41
const TLorentzVector & Momentum() const
Definition: MCStep.h:38
const_iterator end() const
Definition: Assns.h:461
A trajectory in space reconstructed from hits.
Definition: Trajectory.h:67
int ID() const
Definition: Track.h:198
KalmanFilterFinalTrackFitter & operator=(KalmanFilterFinalTrackFitter const &)=delete
const SMatrixSym55 & EndCovarianceLocal5D() const
Accessors to track parameters and covariance matrices in Local5D and Global6D coordinates.
Definition: Track.h:207
std::vector< recob::TrackFitHitInfo > trackFitHitInfos()
get the output vector of TrackFitHitInfos by releasing and moving
Definition: TrackMaker.h:114
TDirectory * dir
Definition: macro.C:5
std::unique_ptr< art::FindManyP< anab::Calorimetry > > trackCalo
std::unique_ptr< art::FindManyP< recob::Shower > > assocShowers
Point_t const & End() const
Returns the position of the last valid point of the trajectory [cm].
const std::string & Process() const
Definition: MCTrack.h:43
int setPId(const unsigned int iTrack, const std::unique_ptr< art::FindManyP< anab::ParticleID > > &trackId, const int pfPid=0) const
PointFlags_t const & FlagsAtPoint(size_t i) const
Access to i-th TrajectoryPoint or its Flags.
Definition: Track.h:118
const MCStep & Start() const
Definition: MCTrack.h:44
TrackCollectionProxyElement< TrackCollProxy > Track
Proxy to an element of a proxy collection of recob::Track objects.
Definition: Track.h:1036
Float_t e
Definition: plot.C:34
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:692
Float_t track
Definition: plot.C:34
int ID() const
Definition: Shower.h:187
const SMatrixSym55 & VertexCovarianceLocal5D() const
Accessors to track parameters and covariance matrices in Local5D and Global6D coordinates.
Definition: Track.h:206
Struct holding optional TrackMaker outputs.
Definition: TrackMaker.h:73
Set of flags pertaining a point of the track.
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
Definition: Track.h:49
calorimetry
double GetTrackMomentum(double trkrange, int pdg) const
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::unique_ptr< art::FindManyP< recob::Track > > assocTracks