Public Member Functions | Protected Attributes | List of all members
sim::MCTrackRecoAlg Class Reference

#include <MCTrackRecoAlg.h>

Public Member Functions

 MCTrackRecoAlg (fhicl::ParameterSet const &pset)
 Default constructor with fhicl parameters. More...
 
std::unique_ptr< std::vector< sim::MCTrack > > Reconstruct (MCRecoPart &part_v, MCRecoEdep &edep_v)
 

Protected Attributes

bool fDebugMode
 

Detailed Description

Definition at line 21 of file MCTrackRecoAlg.h.

Constructor & Destructor Documentation

sim::MCTrackRecoAlg::MCTrackRecoAlg ( fhicl::ParameterSet const &  pset)
explicit

Default constructor with fhicl parameters.

Definition at line 24 of file MCTrackRecoAlg.cxx.

26  {
27  fDebugMode = pset.get<bool>("DebugMode");
28  }

Member Function Documentation

std::unique_ptr< std::vector< sim::MCTrack > > sim::MCTrackRecoAlg::Reconstruct ( MCRecoPart part_v,
MCRecoEdep edep_v 
)

Definition at line 30 of file MCTrackRecoAlg.cxx.

32  {
33  auto result = std::make_unique<std::vector<sim::MCTrack>>();
34  auto& mctracks = *result;
35  auto pindex = details::createPlaneIndexMap();
36 
37  for(size_t i=0; i<part_v.size(); ++i) {
38  auto const& mini_part = part_v[i];
39  if( part_v._pdg_list.find(mini_part._pdgcode) == part_v._pdg_list.end() ) continue;
40 
41  ::sim::MCTrack mini_track;
42 
43  std::vector<double> dEdx;
44  std::vector<std::vector<double> > dQdx;
45  dQdx.resize(3);
46 
47  mini_track.Origin ( mini_part._origin );
48  mini_track.PdgCode ( mini_part._pdgcode );
49  mini_track.TrackID ( mini_part._track_id );
50  mini_track.Process ( mini_part._process );
51  mini_track.Start ( MCStep( mini_part._start_vtx, mini_part._start_mom ) );
52  mini_track.End ( MCStep( mini_part._end_vtx, mini_part._end_mom ) );
53 
54  unsigned int mother_track = part_v.MotherTrackID(i);
55  unsigned int ancestor_track = part_v.AncestorTrackID(i);
56 
57  if(mother_track == kINVALID_UINT || ancestor_track == kINVALID_UINT)
58 
59  throw cet::exception(__FUNCTION__) << "LOGIC ERROR: mother/ancestor track ID is invalid!";
60 
61  MCMiniPart mother_part;
62  MCMiniPart ancestor_part;
63 
64  unsigned int mother_index = part_v.TrackToParticleIndex(mother_track);
65  unsigned int ancestor_index = part_v.TrackToParticleIndex(ancestor_track);
66 
67  if(mother_index != kINVALID_UINT) mother_part = part_v[mother_index];
68  else mother_part._track_id = mother_track;
69 
70  if(ancestor_index != kINVALID_UINT) ancestor_part = part_v[ancestor_index];
71  else ancestor_part._track_id = ancestor_track;
72 
73  mini_track.MotherPdgCode ( mother_part._pdgcode );
74  mini_track.MotherTrackID ( mother_part._track_id );
75  mini_track.MotherProcess ( mother_part._process );
76  mini_track.MotherStart ( MCStep( mother_part._start_vtx, mother_part._start_mom ) );
77  mini_track.MotherEnd ( MCStep( mother_part._end_vtx, mother_part._end_mom ) );
78 
79  mini_track.AncestorPdgCode ( ancestor_part._pdgcode );
80  mini_track.AncestorTrackID ( ancestor_part._track_id );
81  mini_track.AncestorProcess ( ancestor_part._process );
82  mini_track.AncestorStart ( MCStep( ancestor_part._start_vtx, ancestor_part._start_mom ) );
83  mini_track.AncestorEnd ( MCStep( ancestor_part._end_vtx, ancestor_part._end_mom ) );
84 
85 
86  // Fill trajectory points
87 
88  for(auto const& vtx_mom : mini_part._det_path){
89  mini_track.push_back(MCStep(vtx_mom.first,vtx_mom.second));
90  }
91 
92  // No calorimetry for zero length tracks...
93  // JZ : I think we should remove zero length MCTracks because I do not see their utility
94  // JZ : Someone could make this a fcl parameter, I did not
95  if(mini_track.size() == 0){
96  mctracks.push_back(mini_track);
97  continue;
98  }
99 
100  auto const& edep_index = edep_v.TrackToEdepIndex(mini_part._track_id);
101  if(edep_index < 0 ) continue;
102  auto const& edeps = edep_v.GetEdepArrayAt(edep_index);
103 
104  //int n = 0; // unused
105 
106  for(auto const& step_trk : mini_track){
107 
108  if( int(&step_trk - &mini_track[0])+1 == int(mini_track.size()) ){ //annoying way to check if this is last step
109  continue;}
110 
111 
112  auto const& nxt_step_trk = mini_track.at(int(&step_trk - &mini_track[0])+1);
113 
114  //Defining the track step-by-step dEdx and dQdx
115 
116  //Find the distance between two adjacent MCSteps
117  double dist = sqrt(pow(step_trk.X() - nxt_step_trk.X(),2) +
118  pow(step_trk.Y() - nxt_step_trk.Y(),2) +
119  pow(step_trk.Z() - nxt_step_trk.Z(),2));
120 
121  //Make a plane at the step pointed at the next step
122 
123  //Need to define a plane through the first MCStep with a normal along the momentum vector of the step
124  //The plane will be defined in the typical way:
125  // a*x + b*y + c*z + d = 0
126  // where, a = dir_x, b = dir_y, c = dir_z, d = - (a*x_0+b*y_0+c*z_0)
127  // then the *signed* distance of any point (x_1, y_1, z_1) from this plane is:
128  // D = (a*x_1 + b*y_1 + c*z_1 + d )/sqrt( pow(a,2) + pow(b,2) + pow(c,2))
129 
130  double a = 0, b = 0, c = 0, d = 0;
131  a = nxt_step_trk.X() - step_trk.X();
132  b = nxt_step_trk.Y() - step_trk.Y();
133  c = nxt_step_trk.Z() - step_trk.Z();
134  d = -1*(a*step_trk.X() + b*step_trk.Y() + c*step_trk.Z());
135 
136  //Make a line connecting the two points and find the distance from that line
137  //
138  // [A dot B]^2 [A dot B]^2
139  // distance^2 = A^2 - 2* ____________ + ______________
140  // B^2 B^2
141  // Test point == x_0
142  // First Step == x_1
143  // Next step == x_2
144  // A = x_1 - x_0
145  // B = x_2 - x_1
146 
147  // 'B' definition
148  TVector3 B(nxt_step_trk.Position().X() - step_trk.Position().X(),
149  nxt_step_trk.Position().Y() - step_trk.Position().Y(),
150  nxt_step_trk.Position().Z() - step_trk.Position().Z());
151 
152 
153  //Initialize the step-by-step dEdx and dQdx containers
154  double step_dedx = 0;
155  std::vector<double> step_dqdx;
156  step_dqdx.resize(3);
157 
158  //Iterate through all the energy deposition points
159  for(auto const& edep : edeps){
160  // 'x_0' definition
161  TVector3 x_0(edep.pos._x, edep.pos._y, edep.pos._z);
162  // 'A' definition
163  TVector3 A(step_trk.Position().X() - x_0.X(),
164  step_trk.Position().Y() - x_0.Y(),
165  step_trk.Position().Z() - x_0.Z());
166 
167  // Distance from the line connecting x_1 and x_2
168  double LineDist = 0;
169 
170  if(B.Mag2() != 0){
171  LineDist = sqrt(A.Mag2() - 2*pow(A*B,2)/B.Mag2() + pow(A*B,2)/B.Mag2());
172  }
173  else{LineDist = 0;}
174 
175  //Planar Distance and Radial Line Distance Cuts
176  // Add in a voxel before and after to account for MCSteps
177  // the line distance allows for 1mm GEANT multiple columb scattering correction,
178  // small compared to average MCStep-to-MCStep distance
179  if( (a*edep.pos._x + b*edep.pos._y + c*edep.pos._z + d)/sqrt( pow(a,2) + pow(b,2) + pow(c,2)) <= dist + 0.03 &&
180  (a*edep.pos._x + b*edep.pos._y + c*edep.pos._z + d)/sqrt( pow(a,2) + pow(b,2) + pow(c,2)) >= 0 - 0.03 &&
181  LineDist < 0.1){
182 
183  //dEdx Calculation
184  int npid = 0;
185  double engy = 0;
186 
187  for(auto const& pid_energy : edep.deps){
188  engy += pid_energy.energy;
189  npid++;
190  }
191 
192  if(npid != 0){
193  engy /= npid;}
194  else{engy = 0;}
195 
196  step_dedx += engy;
197  auto const pid = edep.pid;
198  auto q_i = pindex.find(pid);
199  if(q_i != pindex.end())
200  step_dqdx[pid.Plane] += (double)(edep.deps[pindex[pid]].charge);
201  }
202  }
203 
204  // Normalize to the 3D distance between the MCSteps
205 
206  //Disregard any energy deposition when 2 MCSteps are separated less than the voxel size
207  if(dist > 0.03){
208  step_dedx /= dist;
209  step_dqdx[0] /= dist;
210  step_dqdx[1] /= dist;
211  step_dqdx[2] /= dist;
212  }
213  else{
214  step_dedx = 0;
215  step_dqdx[0] = 0;
216  step_dqdx[1] = 0;
217  step_dqdx[2] = 0;
218  }
219 
220  // Build the vector(s) to add to data product
221  dEdx.push_back(step_dedx);
222  dQdx[0].push_back(step_dqdx[0]);
223  dQdx[1].push_back(step_dqdx[1]);
224  dQdx[2].push_back(step_dqdx[2]);
225 
226 
227 
228  }
229 
230  //Add calorimetry to the data product
231  mini_track.dEdx(dEdx);
232  mini_track.dQdx(dQdx);
233 
234 
235  mctracks.push_back(mini_track);
236 
237 
238 
239  }
240 
241  if(fDebugMode) {
242 
243  for(auto const& prof : mctracks) {
244 
245  std::cout
246 
247  << Form(" Track particle: PDG=%d : Track ID=%d Start @ (%g,%g,%g,%g) with Momentum (%g,%g,%g,%g)",
248  prof.PdgCode(),prof.TrackID(),
249  prof.Start().X(),prof.Start().Y(),prof.Start().Z(),prof.Start().T(),
250  prof.Start().Px(),prof.Start().Py(),prof.Start().Pz(),prof.Start().E())
251  << std::endl
252  << Form(" Mother particle: PDG=%d : Track ID=%d Start @ (%g,%g,%g,%g) with Momentum (%g,%g,%g,%g)",
253  prof.MotherPdgCode(),prof.MotherTrackID(),
254  prof.MotherStart().X(),prof.MotherStart().Y(),prof.MotherStart().Z(),prof.MotherStart().T(),
255  prof.MotherStart().Px(),prof.MotherStart().Py(),prof.MotherStart().Pz(),prof.MotherStart().E())
256  << std::endl
257  << Form(" Ancestor particle: PDG=%d : Track ID=%d Start @ (%g,%g,%g,%g) with Momentum (%g,%g,%g,%g)",
258  prof.AncestorPdgCode(),prof.AncestorTrackID(),
259  prof.AncestorStart().X(),prof.AncestorStart().Y(),prof.AncestorStart().Z(),prof.AncestorStart().T(),
260  prof.AncestorStart().Px(),prof.AncestorStart().Py(),prof.AncestorStart().Pz(),prof.AncestorStart().E())
261  << std::endl
262  << Form(" ... with %zu trajectory points!",prof.size())
263  << std::endl;
264 
265  if(prof.size()) {
266  std::cout
267  << Form(" Start @ (%g,%g,%g,%g) with Momentum (%g,%g,%g,%g)",
268  prof[0].X(), prof[0].Y(), prof[0].Z(), prof[0].T(),
269  prof[0].Px(), prof[0].Py(), prof[0].Pz(), prof[0].E())
270  << std::endl
271  << Form(" End @ (%g,%g,%g,%g) with Momentum (%g,%g,%g,%g)",
272  (*prof.rbegin()).X(), (*prof.rbegin()).Y(), (*prof.rbegin()).Z(), (*prof.rbegin()).T(),
273  (*prof.rbegin()).Px(), (*prof.rbegin()).Py(), (*prof.rbegin()).Pz(), (*prof.rbegin()).E())
274  << std::endl;
275  }
276  }
277 
278  std::cout<<std::endl<<std::endl;
279  }
280  return result;
281  }
simb::Origin_t Origin() const
Definition: MCTrack.h:40
static QCString result
const std::string & AncestorProcess() const
Definition: MCTrack.h:57
constexpr T pow(T x)
Definition: pow.h:72
const MCStep & MotherEnd() const
Definition: MCTrack.h:53
unsigned int AncestorTrackID() const
Definition: MCTrack.h:56
std::map< geo::PlaneID, size_t > createPlaneIndexMap()
Definition: MCRecoEdep.cxx:22
int AncestorPdgCode() const
Definition: MCTrack.h:55
const MCStep & End() const
Definition: MCTrack.h:45
unsigned int MotherTrackID() const
Definition: MCTrack.h:50
double dEdx(float dqdx, float Efield)
Definition: doAna.cpp:21
const double a
const MCStep & AncestorStart() const
Definition: MCTrack.h:58
int PdgCode() const
Definition: MCTrack.h:41
const MCStep & MotherStart() const
Definition: MCTrack.h:52
int MotherPdgCode() const
Definition: MCTrack.h:49
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
E
Definition: 018_def.c:13
const std::string & Process() const
Definition: MCTrack.h:43
const std::string & MotherProcess() const
Definition: MCTrack.h:51
const unsigned int kINVALID_UINT
Definition: MCLimits.h:14
#define A
Definition: memgrp.cpp:38
static bool * b
Definition: config.cpp:1043
const MCStep & Start() const
Definition: MCTrack.h:44
unsigned int TrackID() const
Definition: MCTrack.h:42
prof
Definition: tracks.py:485
const MCStep & AncestorEnd() const
Definition: MCTrack.h:59
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
QTextStream & endl(QTextStream &s)

Member Data Documentation

bool sim::MCTrackRecoAlg::fDebugMode
protected

Definition at line 30 of file MCTrackRecoAlg.h.


The documentation for this class was generated from the following files: