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