LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
MCTrackRecoAlg.cxx
Go to the documentation of this file.
1 //
3 // MCTrackRecoAlg source
4 //
5 // dEdx and dQdx Estimates Added by Joseph Zennamo (jaz8600@fnal.gov)
6 //
8 
9 #include "MCTrackRecoAlg.h"
10 #include <iostream>
11 
12 #include "cetlib_except/exception.h"
13 #include "fhiclcpp/ParameterSet.h" // for ParameterSet
14 
15 #include "larcoreobj/SimpleTypesAndConstants/geo_types.h" // for PlaneID
16 #include "lardataobj/MCBase/MCLimits.h" // for kINVALID_UINT
17 #include "lardataobj/MCBase/MCStep.h" // for MCStep
18 #include "lardataobj/MCBase/MCTrack.h" // for MCTrack
19 #include "larsim/MCSTReco/MCRecoEdep.h" // for MCEdep
20 #include "larsim/MCSTReco/MCRecoPart.h" // for MCMiniPart
21 
22 #include "TString.h"
23 
24 namespace sim {
25 
26  //##################################################################
28  //##################################################################
29  {
30  fDebugMode = pset.get<bool>("DebugMode");
31  }
32 
33  std::unique_ptr<std::vector<sim::MCTrack>> MCTrackRecoAlg::Reconstruct(MCRecoPart& part_v,
34  MCRecoEdep& edep_v)
35  {
36  auto result = std::make_unique<std::vector<sim::MCTrack>>();
37  auto& mctracks = *result;
38  auto pindex = details::createPlaneIndexMap();
39 
40  for (size_t i = 0; i < part_v.size(); ++i) {
41  auto const& mini_part = part_v[i];
42  if (part_v._pdg_list.find(mini_part._pdgcode) == part_v._pdg_list.end()) continue;
43 
44  ::sim::MCTrack mini_track;
45 
46  std::vector<double> dEdx;
47  std::vector<std::vector<double>> dQdx;
48  dQdx.resize(3);
49 
50  mini_track.Origin(mini_part._origin);
51  mini_track.PdgCode(mini_part._pdgcode);
52  mini_track.TrackID(mini_part._track_id);
53  mini_track.Process(mini_part._process);
54  mini_track.Start(MCStep(mini_part._start_vtx, mini_part._start_mom));
55  mini_track.End(MCStep(mini_part._end_vtx, mini_part._end_mom));
56 
57  unsigned int mother_track = part_v.MotherTrackID(i);
58  unsigned int ancestor_track = part_v.AncestorTrackID(i);
59 
60  if (mother_track == kINVALID_UINT || ancestor_track == kINVALID_UINT)
61 
62  throw cet::exception(__FUNCTION__) << "LOGIC ERROR: mother/ancestor track ID is invalid!";
63 
64  MCMiniPart mother_part;
65  MCMiniPart ancestor_part;
66 
67  unsigned int mother_index = part_v.TrackToParticleIndex(mother_track);
68  unsigned int ancestor_index = part_v.TrackToParticleIndex(ancestor_track);
69 
70  if (mother_index != kINVALID_UINT)
71  mother_part = part_v[mother_index];
72  else
73  mother_part._track_id = mother_track;
74 
75  if (ancestor_index != kINVALID_UINT)
76  ancestor_part = part_v[ancestor_index];
77  else
78  ancestor_part._track_id = ancestor_track;
79 
80  mini_track.MotherPdgCode(mother_part._pdgcode);
81  mini_track.MotherTrackID(mother_part._track_id);
82  mini_track.MotherProcess(mother_part._process);
83  mini_track.MotherStart(MCStep(mother_part._start_vtx, mother_part._start_mom));
84  mini_track.MotherEnd(MCStep(mother_part._end_vtx, mother_part._end_mom));
85 
86  mini_track.AncestorPdgCode(ancestor_part._pdgcode);
87  mini_track.AncestorTrackID(ancestor_part._track_id);
88  mini_track.AncestorProcess(ancestor_part._process);
89  mini_track.AncestorStart(MCStep(ancestor_part._start_vtx, ancestor_part._start_mom));
90  mini_track.AncestorEnd(MCStep(ancestor_part._end_vtx, ancestor_part._end_mom));
91 
92  // Fill trajectory points
93 
94  for (auto const& vtx_mom : mini_part._det_path) {
95  mini_track.push_back(MCStep(vtx_mom.first, vtx_mom.second));
96  }
97 
98  // No calorimetry for zero length tracks...
99  // JZ : I think we should remove zero length MCTracks because I do not see their utility
100  // JZ : Someone could make this a fcl parameter, I did not
101  if (mini_track.size() == 0) {
102  mctracks.push_back(mini_track);
103  continue;
104  }
105 
106  auto const& edep_index = edep_v.TrackToEdepIndex(mini_part._track_id);
107  if (edep_index < 0) continue;
108  auto const& edeps = edep_v.GetEdepArrayAt(edep_index);
109 
110  //int n = 0; // unused
111 
112  for (auto const& step_trk : mini_track) {
113 
114  if (int(&step_trk - &mini_track[0]) + 1 ==
115  int(mini_track.size())) { //annoying way to check if this is last step
116  continue;
117  }
118 
119  auto const& nxt_step_trk = mini_track.at(int(&step_trk - &mini_track[0]) + 1);
120 
121  //Defining the track step-by-step dEdx and dQdx
122 
123  //Find the distance between two adjacent MCSteps
124  double dist =
125  sqrt(pow(step_trk.X() - nxt_step_trk.X(), 2) + pow(step_trk.Y() - nxt_step_trk.Y(), 2) +
126  pow(step_trk.Z() - nxt_step_trk.Z(), 2));
127 
128  //Make a plane at the step pointed at the next step
129 
130  //Need to define a plane through the first MCStep with a normal along the momentum vector of the step
131  //The plane will be defined in the typical way:
132  // a*x + b*y + c*z + d = 0
133  // where, a = dir_x, b = dir_y, c = dir_z, d = - (a*x_0+b*y_0+c*z_0)
134  // then the *signed* distance of any point (x_1, y_1, z_1) from this plane is:
135  // D = (a*x_1 + b*y_1 + c*z_1 + d )/sqrt( pow(a,2) + pow(b,2) + pow(c,2))
136 
137  double a = 0, b = 0, c = 0, d = 0;
138  a = nxt_step_trk.X() - step_trk.X();
139  b = nxt_step_trk.Y() - step_trk.Y();
140  c = nxt_step_trk.Z() - step_trk.Z();
141  d = -1 * (a * step_trk.X() + b * step_trk.Y() + c * step_trk.Z());
142 
143  //Make a line connecting the two points and find the distance from that line
144  //
145  // [A dot B]^2 [A dot B]^2
146  // distance^2 = A^2 - 2* ____________ + ______________
147  // B^2 B^2
148  // Test point == x_0
149  // First Step == x_1
150  // Next step == x_2
151  // A = x_1 - x_0
152  // B = x_2 - x_1
153 
154  // 'B' definition
155  TVector3 B(nxt_step_trk.Position().X() - step_trk.Position().X(),
156  nxt_step_trk.Position().Y() - step_trk.Position().Y(),
157  nxt_step_trk.Position().Z() - step_trk.Position().Z());
158 
159  //Initialize the step-by-step dEdx and dQdx containers
160  double step_dedx = 0;
161  std::vector<double> step_dqdx;
162  step_dqdx.resize(3);
163 
164  //Iterate through all the energy deposition points
165  for (auto const& edep : edeps) {
166  // 'x_0' definition
167  TVector3 x_0(edep.pos._x, edep.pos._y, edep.pos._z);
168  // 'A' definition
169  TVector3 A(step_trk.Position().X() - x_0.X(),
170  step_trk.Position().Y() - x_0.Y(),
171  step_trk.Position().Z() - x_0.Z());
172 
173  // Distance from the line connecting x_1 and x_2
174  double LineDist = 0;
175 
176  if (B.Mag2() != 0) {
177  LineDist = sqrt(A.Mag2() - 2 * pow(A * B, 2) / B.Mag2() + pow(A * B, 2) / B.Mag2());
178  }
179  else {
180  LineDist = 0;
181  }
182 
183  //Planar Distance and Radial Line Distance Cuts
184  // Add in a voxel before and after to account for MCSteps
185  // the line distance allows for 1mm GEANT multiple columb scattering correction,
186  // small compared to average MCStep-to-MCStep distance
187  if ((a * edep.pos._x + b * edep.pos._y + c * edep.pos._z + d) /
188  sqrt(pow(a, 2) + pow(b, 2) + pow(c, 2)) <=
189  dist + 0.03 &&
190  (a * edep.pos._x + b * edep.pos._y + c * edep.pos._z + d) /
191  sqrt(pow(a, 2) + pow(b, 2) + pow(c, 2)) >=
192  0 - 0.03 &&
193  LineDist < 0.1) {
194 
195  //dEdx Calculation
196  int npid = 0;
197  double engy = 0;
198 
199  for (auto const& pid_energy : edep.deps) {
200  engy += pid_energy.energy;
201  npid++;
202  }
203 
204  if (npid != 0) { engy /= npid; }
205  else {
206  engy = 0;
207  }
208 
209  step_dedx += engy;
210  auto const pid = edep.pid;
211  auto q_i = pindex.find(pid);
212  if (q_i != pindex.end())
213  step_dqdx[pid.Plane] += (double)(edep.deps[pindex[pid]].charge);
214  }
215  }
216 
217  // Normalize to the 3D distance between the MCSteps
218 
219  //Disregard any energy deposition when 2 MCSteps are separated less than the voxel size
220  if (dist > 0.03) {
221  step_dedx /= dist;
222  step_dqdx[0] /= dist;
223  step_dqdx[1] /= dist;
224  step_dqdx[2] /= dist;
225  }
226  else {
227  step_dedx = 0;
228  step_dqdx[0] = 0;
229  step_dqdx[1] = 0;
230  step_dqdx[2] = 0;
231  }
232 
233  // Build the vector(s) to add to data product
234  dEdx.push_back(step_dedx);
235  dQdx[0].push_back(step_dqdx[0]);
236  dQdx[1].push_back(step_dqdx[1]);
237  dQdx[2].push_back(step_dqdx[2]);
238  }
239 
240  //Add calorimetry to the data product
241  mini_track.dEdx(dEdx);
242  mini_track.dQdx(dQdx);
243 
244  mctracks.push_back(mini_track);
245  }
246 
247  if (fDebugMode) {
248 
249  for (auto const& prof : mctracks) {
250 
251  std::cout
252 
253  << Form(" Track particle: PDG=%d : Track ID=%d Start @ (%g,%g,%g,%g) with Momentum "
254  "(%g,%g,%g,%g)",
255  prof.PdgCode(),
256  prof.TrackID(),
257  prof.Start().X(),
258  prof.Start().Y(),
259  prof.Start().Z(),
260  prof.Start().T(),
261  prof.Start().Px(),
262  prof.Start().Py(),
263  prof.Start().Pz(),
264  prof.Start().E())
265  << std::endl
266  << Form(" Mother particle: PDG=%d : Track ID=%d Start @ (%g,%g,%g,%g) with Momentum "
267  "(%g,%g,%g,%g)",
268  prof.MotherPdgCode(),
269  prof.MotherTrackID(),
270  prof.MotherStart().X(),
271  prof.MotherStart().Y(),
272  prof.MotherStart().Z(),
273  prof.MotherStart().T(),
274  prof.MotherStart().Px(),
275  prof.MotherStart().Py(),
276  prof.MotherStart().Pz(),
277  prof.MotherStart().E())
278  << std::endl
279  << Form(" Ancestor particle: PDG=%d : Track ID=%d Start @ (%g,%g,%g,%g) with Momentum "
280  "(%g,%g,%g,%g)",
281  prof.AncestorPdgCode(),
282  prof.AncestorTrackID(),
283  prof.AncestorStart().X(),
284  prof.AncestorStart().Y(),
285  prof.AncestorStart().Z(),
286  prof.AncestorStart().T(),
287  prof.AncestorStart().Px(),
288  prof.AncestorStart().Py(),
289  prof.AncestorStart().Pz(),
290  prof.AncestorStart().E())
291  << std::endl
292  << Form(" ... with %zu trajectory points!", prof.size()) << std::endl;
293 
294  if (prof.size()) {
295  std::cout << Form(" Start @ (%g,%g,%g,%g) with Momentum (%g,%g,%g,%g)",
296  prof[0].X(),
297  prof[0].Y(),
298  prof[0].Z(),
299  prof[0].T(),
300  prof[0].Px(),
301  prof[0].Py(),
302  prof[0].Pz(),
303  prof[0].E())
304  << std::endl
305  << Form(" End @ (%g,%g,%g,%g) with Momentum (%g,%g,%g,%g)",
306  (*prof.rbegin()).X(),
307  (*prof.rbegin()).Y(),
308  (*prof.rbegin()).Z(),
309  (*prof.rbegin()).T(),
310  (*prof.rbegin()).Px(),
311  (*prof.rbegin()).Py(),
312  (*prof.rbegin()).Pz(),
313  (*prof.rbegin()).E())
314  << std::endl;
315  }
316  }
317 
318  std::cout << std::endl << std::endl;
319  }
320  return result;
321  }
322 }
std::string _process
Definition: MCRecoPart.h:31
simb::Origin_t Origin() const
Definition: MCTrack.h:38
const std::string & AncestorProcess() const
Definition: MCTrack.h:58
Float_t Y
Definition: plot.C:37
const MCStep & MotherEnd() const
Definition: MCTrack.h:54
unsigned int AncestorTrackID() const
Definition: MCTrack.h:57
std::map< geo::PlaneID, size_t > createPlaneIndexMap()
Definition: MCRecoEdep.cxx:26
unsigned int MotherTrackID(const unsigned int part_index) const
Definition: MCRecoPart.cxx:50
const std::vector< sim::MCEdep > & GetEdepArrayAt(size_t edep_index) const
Returns a vector of MCEdep object at the given index.
Definition: MCRecoEdep.cxx:47
int AncestorPdgCode() const
Definition: MCTrack.h:56
TLorentzVector _start_vtx
Definition: MCRecoPart.h:35
const MCStep & End() const
Definition: MCTrack.h:43
Class def header for mcstep data container.
unsigned int MotherTrackID() const
Definition: MCTrack.h:51
TLorentzVector _start_mom
Definition: MCRecoPart.h:36
Float_t Z
Definition: plot.C:37
unsigned int AncestorTrackID(const unsigned int part_index)
Definition: MCRecoPart.cxx:76
int TrackToEdepIndex(unsigned int track_id) const
Converts a track ID to MCEdep array index. Returns -1 if no corresponding array found ...
Definition: MCRecoEdep.h:113
Float_t E
Definition: plot.C:20
TLorentzVector _end_mom
Definition: MCRecoPart.h:38
T get(std::string const &key) const
Definition: ParameterSet.h:314
Float_t d
Definition: plot.C:235
TLorentzVector _end_vtx
Definition: MCRecoPart.h:37
Class def header for mctrack data container.
const MCStep & AncestorStart() const
Definition: MCTrack.h:59
Double_t edep
Definition: macro.C:13
Monte Carlo Simulation.
Definition of data types for geometry description.
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
Definition: PFPUtils.cxx:2675
int PdgCode() const
Definition: MCTrack.h:39
MCTrackRecoAlg(fhicl::ParameterSet const &pset)
Default constructor with fhicl parameters.
const MCStep & MotherStart() const
Definition: MCTrack.h:53
int MotherPdgCode() const
Definition: MCTrack.h:50
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
std::set< int > _pdg_list
PDG code list for which particle&#39;s trajectory within the detector is saved.
Definition: MCRecoPart.h:127
const std::string & Process() const
Definition: MCTrack.h:41
const std::string & MotherProcess() const
Definition: MCTrack.h:52
const unsigned int kINVALID_UINT
Definition: MCLimits.h:14
const MCStep & Start() const
Definition: MCTrack.h:42
unsigned int _track_id
Definition: MCRecoPart.h:30
std::unique_ptr< std::vector< sim::MCTrack > > Reconstruct(MCRecoPart &part_v, MCRecoEdep &edep_v)
unsigned int TrackID() const
Definition: MCTrack.h:40
Float_t X
Definition: plot.C:37
unsigned int TrackToParticleIndex(const unsigned int track_id) const
Definition: MCRecoPart.h:112
const MCStep & AncestorEnd() const
Definition: MCTrack.h:60
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33