LArSoft  v08_35_01
Liquid Argon Software toolkit - http://larsoft.org/
TrackProducerFromPFParticle_module.cc
Go to the documentation of this file.
13 #include "cetlib_except/exception.h"
14 //
15 #include <memory>
16 //
24 //
50 //
51 //
53 public:
55  // The compiler-generated destructor is fine for non-base
56  // classes without bare pointers or other resource use.
57  //
58  // Plugins should not be copied or assigned.
63 private:
64  // Required functions.
65  void produce(art::Event & e) override;
66  std::unique_ptr<trkmkr::TrackMaker> trackMaker_;
77 };
78 //
80  : EDProducer{p}
81  , trackMaker_{art::make_tool<trkmkr::TrackMaker>(p.get<fhicl::ParameterSet>("trackMaker"))}
82  , pfpInputTag{p.get<art::InputTag>("inputCollection")}
83  , doTrackFitHitInfo_{p.get<bool>("doTrackFitHitInfo")}
84  , doSpacePoints_{p.get<bool>("doSpacePoints")}
85  , spacePointsFromTrajP_{p.get<bool>("spacePointsFromTrajP")}
86  , trackFromPF_{p.get<bool>("trackFromPF")}
87  , showerFromPF_{p.get<bool>("showerFromPF")}
88  , seedFromPF_{p.get<bool>("seedFromPF")}
89 {
90  //
91  if (p.has_key("trackInputTag")) trkInputTag = p.get<art::InputTag>("trackInputTag");
92  else trkInputTag = pfpInputTag;
93  if (p.has_key("showerInputTag")) shwInputTag = p.get<art::InputTag>("showerInputTag");
94  else shwInputTag = pfpInputTag;
95  if (p.has_key("clusterInputTag")) clsInputTag = p.get<art::InputTag>("clusterInputTag");
96  else clsInputTag = pfpInputTag;
97  produces<std::vector<recob::Track> >();
98  produces<art::Assns<recob::Track, recob::Hit, recob::TrackHitMeta> >();
99  produces<art::Assns<recob::PFParticle, recob::Track> >();
100  if (doTrackFitHitInfo_) produces<std::vector<std::vector<recob::TrackFitHitInfo> > >();
101  if (doSpacePoints_) {
102  produces<std::vector<recob::SpacePoint> >();
103  produces<art::Assns<recob::Hit, recob::SpacePoint> >();
104  }
105 }
106 //
108 {
109  // Output collections
110  auto outputTracks = std::make_unique<std::vector<recob::Track> >();
111  auto outputHits = std::make_unique<art::Assns<recob::Track, recob::Hit, recob::TrackHitMeta> >();
112  auto outputPfpTAssn = std::make_unique<art::Assns<recob::PFParticle, recob::Track> >();
113  auto outputHitInfo = std::make_unique<std::vector<std::vector<recob::TrackFitHitInfo> > >();
114  auto outputSpacePoints = std::make_unique<std::vector<recob::SpacePoint> >();
115  auto outputHitSpacePointAssn = std::make_unique<art::Assns<recob::Hit, recob::SpacePoint> >();
116  //
117  // PtrMakers for Assns
118  art::PtrMaker<recob::Track> trackPtrMaker(e);
119  art::PtrMaker<recob::SpacePoint>* spacePointPtrMaker = nullptr;
120  if (doSpacePoints_) spacePointPtrMaker = new art::PtrMaker<recob::SpacePoint>(e);
121  //
122  // Input from event
123  art::ValidHandle<std::vector<recob::PFParticle> > inputPfps = e.getValidHandle<std::vector<recob::PFParticle> >(pfpInputTag);
124  std::unique_ptr<art::FindManyP<recob::Track> > assocTracks;
126  std::unique_ptr<art::FindManyP<recob::Shower> > assocShowers;
127  std::unique_ptr<art::FindManyP<recob::Seed> > assocSeeds;
128  if (trackFromPF_) {
129  assocTracks = std::unique_ptr<art::FindManyP<recob::Track> >(new art::FindManyP<recob::Track>(inputPfps, e, trkInputTag));
130  tkHitsAssn = *e.getValidHandle<art::Assns<recob::Track, recob::Hit> >(trkInputTag);
131  }
132  if (showerFromPF_) assocShowers = std::unique_ptr<art::FindManyP<recob::Shower> >(new art::FindManyP<recob::Shower>(inputPfps, e, shwInputTag));
133  if (seedFromPF_) assocSeeds = std::unique_ptr<art::FindManyP<recob::Seed> >(new art::FindManyP<recob::Seed>(inputPfps, e, pfpInputTag));
134  const auto& trackHitsGroups = util::associated_groups(tkHitsAssn);
135  //
136  std::unique_ptr<art::FindManyP<recob::Cluster> > assocClusters = std::unique_ptr<art::FindManyP<recob::Cluster> >(new art::FindManyP<recob::Cluster>(inputPfps, e, pfpInputTag));
137  auto const& clHitsAssn = *e.getValidHandle<art::Assns<recob::Cluster, recob::Hit> >(clsInputTag);
138  const auto& clusterHitsGroups = util::associated_groups(clHitsAssn);
139  //
140  // Initialize tool for this event
141  trackMaker_->initEvent(e);
142  //
143  // Loop over pfps to fit
144  for (unsigned int iPfp = 0; iPfp < inputPfps->size(); ++iPfp) {
145  //
146  const art::Ptr<recob::PFParticle> pfp(inputPfps, iPfp);
147  //
148  if (trackFromPF_) {
149  // Tracks associated to PFParticles
150  const std::vector<art::Ptr<recob::Track> >& tracks = assocTracks->at(iPfp);
151  // Loop over tracks to refit
152  for (art::Ptr<recob::Track> const& track: tracks) {
153  //
154  // Get track and its hits
155  std::vector<art::Ptr<recob::Hit> > inHits;
156  decltype(auto) hitsRange = util::groupByIndex(trackHitsGroups, track.key());
157  for (art::Ptr<recob::Hit> const& hit: hitsRange) inHits.push_back(hit);
158  //
159  // Declare output objects
160  recob::Track outTrack;
161  std::vector<art::Ptr<recob::Hit> > outHits;
162  trkmkr::OptionalOutputs optionals;
163  if (doTrackFitHitInfo_) optionals.initTrackFitInfos();
165  //
166  // Invoke tool to fit track and fill output objects
167  bool fitok = trackMaker_->makeTrack(track, inHits, outTrack, outHits, optionals);
168  if (!fitok) continue;
169  //
170  // Check that the requirement Nhits == Npoints is satisfied
171  // We also require the hits to the in the same order as the points (this cannot be enforced, can it?)
172  if (outTrack.NumberTrajectoryPoints()!=outHits.size()) {
173  throw cet::exception("TrackProducerFromPFParticle") << "Produced recob::Track required to have 1-1 correspondance between hits and points.\n";
174  }
175  //
176  // Fill output collections, including Assns
177  outputTracks->emplace_back(std::move(outTrack));
178  const art::Ptr<recob::Track> aptr = trackPtrMaker(outputTracks->size()-1);
179  outputPfpTAssn->addSingle(pfp, aptr);
180  unsigned int ip = 0;
181  for (auto const& trhit: outHits) {
182  recob::TrackHitMeta metadata(outputTracks->back().HasValidPoint(ip) ? ip : std::numeric_limits<int>::max(), -std::numeric_limits<double>::max());
183  outputHits->addSingle(aptr, trhit, metadata);
184  //
185  if (doSpacePoints_ && spacePointsFromTrajP_ && outputTracks->back().HasValidPoint(ip)) {
186  auto& tp = outputTracks->back().Trajectory().LocationAtPoint(ip);
187  const double fXYZ[3] = {tp.X(),tp.Y(),tp.Z()};
188  const double fErrXYZ[6] = {0};
189  recob::SpacePoint sp(fXYZ, fErrXYZ, -1.);
190  outputSpacePoints->emplace_back(std::move(sp));
191  const art::Ptr<recob::SpacePoint> apsp = (*spacePointPtrMaker)(outputSpacePoints->size()-1);
192  outputHitSpacePointAssn->addSingle(trhit, apsp);
193  }
194  ip++;
195  }
197  auto osp = optionals.spacePointHitPairs();
198  for (auto it = osp.begin(); it!=osp.end(); ++it ) {
199  outputSpacePoints->emplace_back(std::move(it->first));
200  const art::Ptr<recob::SpacePoint> apsp = (*spacePointPtrMaker)(outputSpacePoints->size()-1);
201  outputHitSpacePointAssn->addSingle(it->second,apsp);
202  }
203  }
204  if (doTrackFitHitInfo_) {
205  outputHitInfo->emplace_back(optionals.trackFitHitInfos());
206  }
207  }
208  }
209  //
210  if (showerFromPF_) {
211  //
212  // Showers associated to PFParticles
213  const std::vector<art::Ptr<recob::Shower> >& showers = assocShowers->at(iPfp);
214  // if there is more than one shower the logic below to get the hits does not work! this works, at least for uboone
215  if (showers.size()!=1) continue;
216  //
217  // Get hits for shower (through the chain pfp->clusters->hits)
218  std::vector<art::Ptr<recob::Hit> > inHits;
219  const std::vector<art::Ptr<recob::Cluster> > clustersRange = assocClusters->at(iPfp);
220  for (art::Ptr<recob::Cluster> const& cluster: clustersRange) {
221  // for hits we use groupByIndex since it preserves the order (and we can use it since each cluster must have associated hits)
222  decltype(auto) hitsRange = util::groupByIndex(clusterHitsGroups, cluster.key());
223  for (art::Ptr<recob::Hit> const& hit: hitsRange) inHits.push_back(hit);
224  }
225  // Loop over showers to refit (should be only one)
226  for (unsigned int iShower = 0; iShower < showers.size(); ++iShower) {
227  //
228  // Get the shower and convert/hack it into a trajectory so that the fit is initialized
229  art::Ptr<recob::Shower> shower = showers[iShower];
230  recob::tracking::Point_t pos(shower->ShowerStart().X(),shower->ShowerStart().Y(),shower->ShowerStart().Z());
231  recob::tracking::Vector_t dir(shower->Direction().X(),shower->Direction().Y(),shower->Direction().Z());
232  float mom = 1.;
233  if (shower->Energy().size()==3) mom = shower->Energy()[2]*0.001;
234  std::vector<recob::tracking::Point_t> p;
235  std::vector<recob::tracking::Vector_t> d;
236  for (unsigned int i=0; i<inHits.size(); ++i) {
237  p.push_back(pos);
238  d.push_back(mom*dir);
239  }
240  recob::TrackTrajectory traj(std::move(p), std::move(d), recob::TrackTrajectory::Flags_t(p.size()), false);
241  //
242  // Declare output objects
243  recob::Track outTrack;
244  std::vector<art::Ptr<recob::Hit> > outHits;
245  trkmkr::OptionalOutputs optionals;
246  if (doTrackFitHitInfo_) optionals.initTrackFitInfos();
248  //
249  // Invoke tool to fit track and fill output objects
250  bool fitok = trackMaker_->makeTrack(traj, iPfp, inHits, outTrack, outHits, optionals);
251  if (!fitok) continue;
252  //
253  // Check that the requirement Nhits == Npoints is satisfied
254  // We also require the hits to the in the same order as the points (this cannot be enforced, can it?)
255  if (outTrack.NumberTrajectoryPoints()!=outHits.size()) {
256  throw cet::exception("TrackProducerFromPFParticle") << "Produced recob::Track required to have 1-1 correspondance between hits and points.\n";
257  }
258  //
259  // Fill output collections, including Assns
260  outputTracks->emplace_back(std::move(outTrack));
261  const art::Ptr<recob::Track> aptr = trackPtrMaker(outputTracks->size()-1);
262  outputPfpTAssn->addSingle(pfp, aptr);
263  unsigned int ip = 0;
264  for (auto const& trhit: outHits) {
265  recob::TrackHitMeta metadata(outputTracks->back().HasValidPoint(ip) ? ip : std::numeric_limits<int>::max(), -std::numeric_limits<double>::max());
266  outputHits->addSingle(aptr, trhit, metadata);
267  //
268  if (doSpacePoints_ && spacePointsFromTrajP_ && outputTracks->back().HasValidPoint(ip)) {
269  auto& tp = outputTracks->back().Trajectory().LocationAtPoint(ip);
270  const double fXYZ[3] = {tp.X(),tp.Y(),tp.Z()};
271  const double fErrXYZ[6] = {0};
272  recob::SpacePoint sp(fXYZ, fErrXYZ, -1.);
273  outputSpacePoints->emplace_back(std::move(sp));
274  const art::Ptr<recob::SpacePoint> apsp = (*spacePointPtrMaker)(outputSpacePoints->size()-1);
275  outputHitSpacePointAssn->addSingle(trhit, apsp);
276  }
277  ip++;
278  }
280  auto osp = optionals.spacePointHitPairs();
281  for (auto it = osp.begin(); it!=osp.end(); ++it ) {
282  outputSpacePoints->emplace_back(std::move(it->first));
283  const art::Ptr<recob::SpacePoint> apsp = (*spacePointPtrMaker)(outputSpacePoints->size()-1);
284  outputHitSpacePointAssn->addSingle(it->second,apsp);
285  }
286  }
287  if (doTrackFitHitInfo_) {
288  outputHitInfo->emplace_back(optionals.trackFitHitInfos());
289  }
290  }
291  }
292  //
293  //
294  if (seedFromPF_) {
295  //
296  // Seeds associated to PFParticles
297  const std::vector<art::Ptr<recob::Seed> >& seeds = assocSeeds->at(iPfp);
298  // if there is more than one seed the logic below to get the hits does not work! this works, at least for uboone
299  if (seeds.size()!=1) continue;
300  //
301  // Get hits for pfp (through the chain pfp->clusters->hits)
302  std::vector<art::Ptr<recob::Hit> > inHits;
303  const std::vector<art::Ptr<recob::Cluster> > clustersRange = assocClusters->at(iPfp);
304  for (art::Ptr<recob::Cluster> const& cluster: clustersRange) {
305  // for hits we use groupByIndex since it preserves the order (and we can use it since each cluster must have associated hits)
306  decltype(auto) hitsRange = util::groupByIndex(clusterHitsGroups, cluster.key());
307  for (art::Ptr<recob::Hit> const& hit: hitsRange) inHits.push_back(hit);
308  }
309  if (inHits.size()<4) continue;
310  // Loop over seeds should be only one)
311  for (unsigned int iS = 0; iS < seeds.size(); ++iS) {
312  //
313  // Get the seed and convert/hack it into a trajectory so that the fit is initialized
314  art::Ptr<recob::Seed> seed = seeds[iS];
315  double p0[3], pe[3];
316  seed->GetPoint(p0,pe);
317  double d0[3], de[3];
318  seed->GetDirection(d0,de);
319  recob::tracking::Point_t pos(p0[0],p0[1],p0[2]);
320  recob::tracking::Vector_t dir(d0[0],d0[1],d0[2]);
321  std::vector<recob::tracking::Point_t> p;
322  std::vector<recob::tracking::Vector_t> d;
323  for (unsigned int i=0; i<inHits.size(); ++i) {
324  p.push_back(pos);
325  d.push_back(dir);
326  }
327  recob::TrackTrajectory traj(std::move(p), std::move(d), recob::TrackTrajectory::Flags_t(p.size()), false);
328  //
329  // Declare output objects
330  recob::Track outTrack;
331  std::vector<art::Ptr<recob::Hit> > outHits;
332  trkmkr::OptionalOutputs optionals;
333  if (doTrackFitHitInfo_) optionals.initTrackFitInfos();
335  //
336  // Invoke tool to fit track and fill output objects
337  bool fitok = trackMaker_->makeTrack(traj, iPfp, inHits, outTrack, outHits, optionals);
338  if (!fitok) continue;
339  //
340  // Check that the requirement Nhits == Npoints is satisfied
341  // We also require the hits to the in the same order as the points (this cannot be enforced, can it?)
342  if (outTrack.NumberTrajectoryPoints()!=outHits.size()) {
343  throw cet::exception("TrackProducerFromPFParticle") << "Produced recob::Track required to have 1-1 correspondance between hits and points.\n";
344  }
345  //
346  // Fill output collections, including Assns
347  outputTracks->emplace_back(std::move(outTrack));
348  const art::Ptr<recob::Track> aptr = trackPtrMaker(outputTracks->size()-1);
349  outputPfpTAssn->addSingle(pfp, aptr);
350  unsigned int ip = 0;
351  for (auto const& trhit: outHits) {
352  recob::TrackHitMeta metadata(outputTracks->back().HasValidPoint(ip) ? ip : std::numeric_limits<int>::max(), -std::numeric_limits<double>::max());
353  outputHits->addSingle(aptr, trhit, metadata);
354  //
355  if (doSpacePoints_ && spacePointsFromTrajP_ && outputTracks->back().HasValidPoint(ip)) {
356  auto& tp = outputTracks->back().Trajectory().LocationAtPoint(ip);
357  const double fXYZ[3] = {tp.X(),tp.Y(),tp.Z()};
358  const double fErrXYZ[6] = {0};
359  recob::SpacePoint sp(fXYZ, fErrXYZ, -1.);
360  outputSpacePoints->emplace_back(std::move(sp));
361  const art::Ptr<recob::SpacePoint> apsp = (*spacePointPtrMaker)(outputSpacePoints->size()-1);
362  outputHitSpacePointAssn->addSingle(trhit, apsp);
363  }
364  ip++;
365  }
367  auto osp = optionals.spacePointHitPairs();
368  for (auto it = osp.begin(); it!=osp.end(); ++it ) {
369  outputSpacePoints->emplace_back(std::move(it->first));
370  const art::Ptr<recob::SpacePoint> apsp = (*spacePointPtrMaker)(outputSpacePoints->size()-1);
371  outputHitSpacePointAssn->addSingle(it->second,apsp);
372  }
373  }
374  if (doTrackFitHitInfo_) {
375  outputHitInfo->emplace_back(optionals.trackFitHitInfos());
376  }
377  }
378  }
379  //
380  }
381  //
382  // Put collections in the event
383  e.put(std::move(outputTracks));
384  e.put(std::move(outputHits));
385  e.put(std::move(outputPfpTAssn));
386  if (doTrackFitHitInfo_) {
387  e.put(std::move(outputHitInfo));
388  }
389  if (doSpacePoints_) {
390  e.put(std::move(outputSpacePoints));
391  e.put(std::move(outputHitSpacePointAssn));
392  }
393  if (doSpacePoints_) delete spacePointPtrMaker;
394 }
395 //
const TVector3 & ShowerStart() const
Definition: Shower.h:192
void initTrackFitInfos()
initialize the output vector of TrackFitHitInfos
Definition: TrackMaker.h:102
EDProducer()=default
void GetPoint(double *Pt, double *Err) const
Definition: Seed.cxx:108
auto groupByIndex(Groups &&groups, std::size_t index) -> decltype(auto)
Returns the group within groups with the specified index.
const std::vector< double > & Energy() const
Definition: Shower.h:195
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.
Data related to recob::Hit associated with recob::Track.The purpose is to collect several variables t...
Definition: TrackHitMeta.h:43
Cluster finding and building.
TrackProducerFromPFParticle(fhicl::ParameterSet const &p)
auto associated_groups(A const &assns)
Helper functions to access associations in order.
std::unique_ptr< trkmkr::TrackMaker > trackMaker_
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< Coord_t >, ROOT::Math::GlobalCoordinateSystemTag > Vector_t
Type for representation of momenta in 3D space. See recob::tracking::Coord_t for more details on the ...
Definition: TrackingTypes.h:29
Int_t max
Definition: plot.C:27
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:56
long seed
Definition: chem4.cc:68
A trajectory in space reconstructed from hits.
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:482
Float_t d
Definition: plot.C:237
const TVector3 & Direction() const
Definition: Shower.h:189
std::vector< PointFlags_t > Flags_t
Type of point flag list.
Declaration of cluster object.
Detector simulation of raw signals on wires.
Produce a reco::Track collection, as a result of the fit of an existing recob::PFParticle collection...
std::vector< SpHitPair > spacePointHitPairs()
get the output vector of SpHitPair by releasing and moving
Definition: TrackMaker.h:119
std::vector< TrajPoint > seeds
Definition: DataStructs.cxx:13
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
Helper functions to access associations in order.
TrackProducerFromPFParticle & operator=(TrackProducerFromPFParticle const &)=delete
Float_t e
Definition: plot.C:34
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< Coord_t >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space. See recob::tracking::Coord_t for more detai...
Definition: TrackingTypes.h:26
Float_t track
Definition: plot.C:34
void GetDirection(double *Dir, double *Err) const
Definition: Seed.cxx:98
void initSpacePoints()
initialize the output vector of SpHitPair
Definition: TrackMaker.h:106
Struct holding optional TrackMaker outputs.
Definition: TrackMaker.h:73
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
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33