LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
TrackCheater_module.cc
Go to the documentation of this file.
1 //
3 // TrackCheater module
4 //
5 // brebel@fnal.gov
6 //
8 // Framework includes
10 
11 #include <string>
12 
13 // ROOT includes
14 #include "TVector3.h"
15 
16 // LArSoft includes
27 
28 // Framework includes
34 #include "fhiclcpp/ParameterSet.h"
36 
37 namespace trkf {
38  class TrackCheater : public art::EDProducer {
39  public:
40  explicit TrackCheater(fhicl::ParameterSet const& pset);
41 
42  private:
43  void produce(art::Event& evt) override;
44 
45  std::string const fCheatedClusterLabel;
46  std::string const fG4ModuleLabel;
48  };
49 }
50 
51 namespace trkf {
52 
53  //--------------------------------------------------------------------
55  : EDProducer{pset}
56  , fCheatedClusterLabel{pset.get<std::string>("CheatedClusterLabel", "cluster")}
57  , fG4ModuleLabel{pset.get<std::string>("G4ModuleLabel", "largeant")}
58  {
59  produces<std::vector<recob::Track>>();
60  produces<std::vector<recob::SpacePoint>>();
61  produces<art::Assns<recob::Track, recob::Cluster>>();
62  produces<art::Assns<recob::Track, recob::SpacePoint>>();
63  produces<art::Assns<recob::Track, recob::Hit>>();
64  }
65 
66  //--------------------------------------------------------------------
68  {
72 
73  // grab the clusters that have been reconstructed
75  evt.getByLabel(fCheatedClusterLabel, clustercol);
76 
78 
79  // make a vector of them - we aren't writing anything out to a file
80  // so no need for a art::PtrVector here
81  std::vector<art::Ptr<recob::Cluster>> clusters;
82  art::fill_ptr_vector(clusters, clustercol);
83 
84  // loop over the clusters and figure out which particle contributed to each one
85  std::vector<art::Ptr<recob::Cluster>>::iterator itr = clusters.begin();
86 
87  // make a map of vectors of art::Ptrs keyed by eveID values
88  std::map<int, std::vector<std::pair<size_t, art::Ptr<recob::Cluster>>>> eveClusterMap;
89 
90  // loop over all clusters and fill in the map
91  for (size_t c = 0; c < clusters.size(); ++c) {
92 
93  std::pair<size_t, art::Ptr<recob::Cluster>> idxPtr(c, clusters[c]);
94 
95  // in the ClusterCheater module we set the cluster ID to be
96  // the eve particle track ID*1000 + plane*100 + tpc*10 + cryostat number. The
97  // floor function on the cluster ID / 1000 will give us
98  // the eve track ID
99  int eveID = floor((*itr)->ID() / 1000.);
100 
101  eveClusterMap[eveID].push_back(idxPtr);
102  ++itr;
103  } // end loop over clusters
104 
105  // loop over the map and make prongs
106  std::unique_ptr<std::vector<recob::Track>> trackcol(new std::vector<recob::Track>);
107  std::unique_ptr<std::vector<recob::SpacePoint>> spcol(new std::vector<recob::SpacePoint>);
108  std::unique_ptr<art::Assns<recob::Track, recob::SpacePoint>> tspassn(
110  std::unique_ptr<art::Assns<recob::Track, recob::Cluster>> tcassn(
112  std::unique_ptr<art::Assns<recob::Track, recob::Hit>> thassn(
114 
115  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
116 
117  for (auto const& clusterMapItr : eveClusterMap) {
118 
119  // separate out the hits for each particle into the different views
120  std::vector<std::pair<size_t, art::Ptr<recob::Cluster>>> const& eveClusters =
121  clusterMapItr.second;
122 
123  simb::MCParticle* part = pi_serv->ParticleList()[clusterMapItr.first];
124 
125  // is this prong electro-magnetic in nature or
126  // hadronic/muonic? EM --> shower, everything else is a track
127  if (abs(part->PdgCode()) != 11 && abs(part->PdgCode()) != 22 && abs(part->PdgCode()) != 111) {
128 
129  // vectors to hold the positions and directions of the track
130  std::vector<TVector3> points;
131  std::vector<TVector3> moms;
132 
133  mf::LogInfo("TrackCheater")
134  << "G4 id " << clusterMapItr.first << " is a track with pdg code " << part->PdgCode();
135 
137  std::vector<size_t> idxs;
138  for (auto const& idxPtr : eveClusters) {
139  ptrvs.push_back(idxPtr.second);
140  idxs.push_back(idxPtr.first);
141  }
142 
143  // grab the hits associated with the collection plane
144  std::vector<art::Ptr<recob::Hit>> hits;
145  for (size_t p = 0; p < ptrvs.size(); ++p) {
146  std::vector<art::Ptr<recob::Hit>> chits = fmh.at(idxs[p]);
147  if (!chits.size()) continue;
148  if (chits[0]->SignalType() != geo::kCollection) continue;
149  hits.insert(hits.end(), chits.begin(), chits.end());
150  }
151 
152  // need at least 2 hits to make a track
153  if (hits.size() < 2) continue;
154 
155  // loop over the hits to get the positions and directions
156  size_t spStart = spcol->size();
157  for (size_t t = 0; t < hits.size(); ++t) {
158  std::vector<double> xyz = bt_serv->HitToXYZ(clockData, hits[t]);
159  TVector3 point(xyz[0], xyz[1], xyz[2]);
160  points.push_back(point);
161 
162  std::vector<double> xyz1;
163  double dx = 0.;
164  double sign = 1.;
165 
166  if (t < hits.size() - 1) { xyz1 = bt_serv->HitToXYZ(clockData, hits[t + 1]); }
167  else {
168  xyz1 = bt_serv->HitToXYZ(clockData, hits[t - 1]);
169  sign = -1.;
170  }
171 
172  // dx is always positive
173  dx = std::sqrt(std::pow(xyz1[0] - xyz[0], 2) + std::pow(xyz1[1] - xyz[1], 2) +
174  std::pow(xyz1[2] - xyz[2], 2));
175 
176  // figure out momentum
177  double mom = 0;
178  double drmin = std::numeric_limits<double>::max();
179  for (unsigned int itp = 0; itp < part->NumberTrajectoryPoints(); itp++) {
180  TVector3 p(part->Vx(itp), part->Vy(itp), part->Vz(itp));
181  double dr = (p - point).Mag();
182  if (dr < drmin) {
183  mom = part->P(itp);
184  drmin = dr;
185  }
186  }
187 
188  // direction is always from the previous point along track to
189  // the next point, that is why sign is there
190  moms.push_back(TVector3(mom * sign * (xyz1[0] - xyz[0]) / dx,
191  mom * sign * (xyz1[1] - xyz[1]) / dx,
192  mom * sign * (xyz1[2] - xyz[2]) / dx));
193 
194  /*************************************************************/
195  /* WARNING */
196  /*************************************************************/
197  /* The dQdx information in recob::Track has been deprecated */
198  /* since 2016 and in 11/2018 the recob::Track interface was */
199  /* changed and DQdxAtPoint and NumberdQdx were removed. */
200  /* Therefore the code below is now commented out */
201  /* (note that it was most likely not functional anyways). */
202  /* For any issue please contact: larsoft-team@fnal.gov */
203  /*************************************************************/
204  /*
205  dQdx[0].push_back(charge/dx);
206  dQdx[1].push_back(charge/dx);
207  dQdx[2].push_back(charge/dx);
208  */
209  /*************************************************************/
210 
211  // make the space point and set its ID and XYZ
212  double xyzerr[6] = {1.e-3};
213  double chisqr = 0.9;
214  recob::SpacePoint sp(&xyz[0], xyzerr, chisqr, clusterMapItr.first * 10000 + t);
215  spcol->push_back(sp);
216  }
217 
218  size_t spEnd = spcol->size();
219 
220  // add a track to the collection. Make the track
221  // ID be the same as the track ID for the eve particle
222  trackcol->push_back(
225  recob::Track::Flags_t(points.size()),
226  true),
227  0,
228  -1.,
229  0,
232  clusterMapItr.first));
233 
234  // associate the track with its clusters
235  util::CreateAssn(evt, *trackcol, ptrvs, *tcassn);
236 
237  // assume the input tracks were previously associated with hits
238  hits.clear();
239  for (size_t p = 0; p < ptrvs.size(); ++p) {
240  hits = fmh.at(idxs[p]);
241  util::CreateAssn(evt, *trackcol, hits, *thassn);
242  }
243 
244  // associate the track to the space points
245  util::CreateAssn(evt, *trackcol, *spcol, *tspassn, spStart, spEnd);
246 
247  mf::LogInfo("TrackCheater") << "adding track: \n" << trackcol->back() << "\nto collection.";
248 
249  } //end if this is a track
250 
251  } // end loop over the map
252 
253  evt.put(std::move(trackcol));
254  evt.put(std::move(spcol));
255  evt.put(std::move(tcassn));
256  evt.put(std::move(thassn));
257  evt.put(std::move(tspassn));
258 
259  return;
260 
261  } // end produce
262 
263 } // end namespace
264 
265 namespace trkf {
266 
268 
269 }
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:219
int PdgCode() const
Definition: MCParticle.h:213
std::string const fCheatedClusterLabel
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
std::vector< Vector_t > convertCollToVector(std::vector< Vector > const &coll)
Definition: TrackingTypes.h:78
Declaration of signal hit object.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
TrackTrajectory::Flags_t Flags_t
Definition: Track.h:67
constexpr auto abs(T v)
Returns the absolute value of the argument.
TrackCheater(fhicl::ParameterSet const &pset)
ROOT::Math::SMatrix< Double32_t, 5, 5, ROOT::Math::MatRepSym< Double32_t, 5 >> SMatrixSym55
std::vector< double > HitToXYZ(detinfo::DetectorClocksData const &clockData, const recob::Hit &hit) const
Particle class.
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
TString part[npart]
Definition: Style.C:32
void hits()
Definition: readHits.C:15
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
A trajectory in space reconstructed from hits.
double P(const int i=0) const
Definition: MCParticle.h:235
Provides recob::Track data product.
std::vector< Point_t > convertCollToPoint(std::vector< Point > const &coll)
Definition: TrackingTypes.h:68
Declaration of cluster object.
size_type size() const
Definition: PtrVector.h:302
const sim::ParticleList & ParticleList() const
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
int sign(double val)
Definition: UtilFunc.cxx:104
double Vx(const int i=0) const
Definition: MCParticle.h:222
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Utility object to perform functions of association.
double Vz(const int i=0) const
Definition: MCParticle.h:224
TCEvent evt
Definition: DataStructs.cxx:8
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:306
void produce(art::Event &evt) override
Namespace collecting geometry-related classes utilities.
std::string const fG4ModuleLabel
label for module running G4 and making particles, etc
double Vy(const int i=0) const
Definition: MCParticle.h:223
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:49
Signal from collection planes.
Definition: geo_types.h:152