LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
MCRecoPart.cxx
Go to the documentation of this file.
1 //
3 // MCRecoPart source
4 //
6 
7 #ifndef MCRECOPART_CXX
8 #define MCRECOPART_CXX
9 
10 #include "MCRecoPart.h"
11 
12 namespace sim {
13 
14  //--------------------------------------------------------------------------------------------
16  //--------------------------------------------------------------------------------------------
17  {
18  this->clear();
19  _track_index.clear();
20  _pdg_list.clear();
21  for(auto const& id : pset.get<std::vector<int> >("SavePathPDGList"))
22 
23  _pdg_list.insert(id);
24 
26  // Build "Fiducial" Volume Definition:
27  //
28  // Iterate over all TPC's to get bounding box that covers volumes of each individual TPC in the detector
29  _x_min = std::min_element(geo->begin_TPC(), geo->end_TPC(), [](auto const &lhs, auto const &rhs){ return lhs.BoundingBox().MinX() < rhs.BoundingBox().MinX();})->MinX();
30  _y_min = std::min_element(geo->begin_TPC(), geo->end_TPC(), [](auto const &lhs, auto const &rhs){ return lhs.BoundingBox().MinY() < rhs.BoundingBox().MinY();})->MinY();
31  _z_min = std::min_element(geo->begin_TPC(), geo->end_TPC(), [](auto const &lhs, auto const &rhs){ return lhs.BoundingBox().MinZ() < rhs.BoundingBox().MinZ();})->MinZ();
32  _x_max = std::max_element(geo->begin_TPC(), geo->end_TPC(), [](auto const &lhs, auto const &rhs){ return lhs.BoundingBox().MaxX() < rhs.BoundingBox().MaxX();})->MaxX();
33  _y_max = std::max_element(geo->begin_TPC(), geo->end_TPC(), [](auto const &lhs, auto const &rhs){ return lhs.BoundingBox().MaxY() < rhs.BoundingBox().MaxY();})->MaxY();
34  _z_max = std::max_element(geo->begin_TPC(), geo->end_TPC(), [](auto const &lhs, auto const &rhs){ return lhs.BoundingBox().MaxZ() < rhs.BoundingBox().MaxZ();})->MaxZ();
35  }
36 
37  //--------------------------------------------------------------------------------------------
38  unsigned int MCRecoPart::MotherTrackID(const unsigned int part_index) const
39  //--------------------------------------------------------------------------------------------
40  {
41  if(this->size() <= part_index) return ::sim::kINVALID_UINT;
42 
43  unsigned int result = this->at(part_index)._mother;
44 
45  if(!result) return this->at(part_index)._track_id;
46 
47  if(TrackToParticleIndex(result) != ::sim::kINVALID_UINT) return result;
48 
49  //std::cout<< "\033[95mWarning:\033[00m Mother particle not in the particle list!"<<std::endl;
50  // Do brute search
51  unsigned int daughter_id = this->at(part_index)._track_id;
52 
53  for(auto const& part : *this) {
54 
55  if(part._daughters.find(daughter_id) != part._daughters.end())
56 
57  return part._track_id;
58 
59  }
60  return result;
61  }
62 
63  //--------------------------------------------------------------------------------------------
64  unsigned int MCRecoPart::AncestorTrackID(const unsigned int part_index)
65  //--------------------------------------------------------------------------------------------
66  {
67  if(part_index >= this->size()) return kINVALID_UINT;
68 
69  if((*this)[part_index]._ancestor != kINVALID_UINT) return (*this)[part_index]._ancestor;
70 
71  unsigned int result = MotherTrackID(part_index);
72 
73  if(result == this->at(part_index)._track_id) return result;
74 
75  if(!result) return this->at(part_index)._track_id;
76 
77  auto mother_index = TrackToParticleIndex(result);
78 
79  while(1) {
80 
81  if(mother_index != kINVALID_UINT) {
82 
83  auto const new_result = MotherTrackID(mother_index);
84 
85  if(new_result == this->at(mother_index)._track_id) break;
86 
87  result = new_result;
88 
89  }else{
90 
91  // Look for a particle that has a daughter = this mother
92  auto const old_result = result;
93  for(auto const& p : *this) {
94 
95  if(p._daughters.find(result) != p._daughters.end()) {
96  result = p._track_id;
97  break;
98  }
99  }
100  if(result == old_result)
101  break;
102  }
103 
104  mother_index = TrackToParticleIndex(result);
105 
106  }
107 
108  (*this)[part_index]._ancestor = result;
109  return result;
110  }
111 
112  //--------------------------------------------------------------------------------------------
113  bool MCRecoPart::InDetector(const double& x,
114  const double& y,
115  const double& z) const
116  //--------------------------------------------------------------------------------------------
117  {
118  return !( x > _x_max || x < _x_min ||
119  z > _z_max || z < _z_min ||
120  y > _y_max || y < _y_min );
121  }
122 
123  //--------------------------------------------------------------------------------------------
124  void MCRecoPart::AddParticles(const std::vector<simb::MCParticle>& mcp_v,
125  const std::vector<simb::Origin_t>& orig_v)
126  //--------------------------------------------------------------------------------------------
127  {
128  if(orig_v.size() != mcp_v.size()) throw cet::exception(__FUNCTION__) << "MCParticle and Origin_t vector size not same!";
129 
130  this->clear();
131  _track_index.clear();
132 
133  for(size_t i=0; i < mcp_v.size(); ++i) {
134 
135  auto const& mcp = mcp_v[i];
136 
137  //std::cout<<" Track ID : "<<mcp.TrackId()<<" ... Index : " <<this->size()<<std::endl;
138 
139  _track_index.insert(std::make_pair((size_t)(mcp.TrackId()),(size_t)(this->size())));
140 
141  this->push_back(MCMiniPart());
142 
143  auto& mini_mcp = (*this->rbegin());
144 
145  for(size_t i=0; i<(size_t)(mcp.NumberDaughters()); ++i)
146  mini_mcp._daughters.insert(mcp.Daughter(i));
147 
148  mini_mcp._track_id = mcp.TrackId();
149  mini_mcp._pdgcode = mcp.PdgCode();
150  mini_mcp._mother = mcp.Mother();
151  mini_mcp._process = mcp.Process();
152  mini_mcp._start_vtx = mcp.Position();
153  mini_mcp._start_mom = mcp.Momentum();
154  mini_mcp._end_vtx = mcp.EndPosition();
155  mini_mcp._end_mom = mcp.EndMomentum();
156  mini_mcp._origin = orig_v[i];
157 
158  // Change units to LArSoft (MeV, cm, us)
159  for(size_t i=0; i<4; ++i) {
160  mini_mcp._start_mom[i] *= 1.e3;
161  mini_mcp._end_mom[i] *= 1.e3;
162  }
163  /*
164  for(size_t i=0; i<3; ++i) {
165  mini_mcp._start_vtx[i] /= 10.;
166  mini_mcp._end_vtx[i] /= 10.;
167  }
168  mini_mcp.start_vtx[3] /= 1.e-3;
169  mini_mcp.end_vtx[3] /= 1.e-3;
170  */
171 
172  if(_pdg_list.find(mcp.PdgCode()) != _pdg_list.end()) {
173 
174  std::set<size_t> det_path_index;
175 
176  for(size_t i=0; i<mcp.NumberTrajectoryPoints(); ++i) {
177 
178  if(InDetector(mcp.Vx(i),mcp.Vy(i),mcp.Vz(i)))
179 
180  det_path_index.insert(i);
181 
182  }
183 
184  if(det_path_index.size()) {
185  if( (*det_path_index.begin()) )
186  det_path_index.insert( (*det_path_index.begin())-1 );
187  if( det_path_index.size()>1 ) {
188  if( ((*det_path_index.rbegin())+1) < mcp.NumberTrajectoryPoints() )
189  det_path_index.insert( (*det_path_index.rbegin())+1 );
190  }
191  mini_mcp._det_path.reserve(det_path_index.size());
192  for(auto const& index : det_path_index) {
193 
194  TLorentzVector vec(mcp.Momentum(index));
195  for(size_t i=0; i<4; ++i) vec[i] *= 1.e3;
196 
197  mini_mcp._det_path.push_back(std::make_pair(mcp.Position(index),vec));
198 
199  }
200  }
201  }
202  }
203  }
204 }
205 #endif
Float_t x
Definition: compare.C:6
double _z_max
z-max of volume box used to determine whether to save track information
Definition: MCRecoPart.h:123
Float_t y
Definition: compare.C:6
Double_t z
Definition: plot.C:279
unsigned int MotherTrackID(const unsigned int part_index) const
Definition: MCRecoPart.cxx:38
double _y_max
y-max of volume box used to determine whether to save track information
Definition: MCRecoPart.h:121
double _y_min
y-min of volume box used to determine whether to save track information
Definition: MCRecoPart.h:122
std::map< unsigned int, unsigned int > _track_index
Track ID => Index Map.
Definition: MCRecoPart.h:111
TPC_iterator begin_TPC() const
Returns an iterator pointing to the first TPC in the detector.
double _z_min
z-min of volume box used to determine whether to save track information
Definition: MCRecoPart.h:124
unsigned int AncestorTrackID(const unsigned int part_index)
Definition: MCRecoPart.cxx:64
MCRecoPart(fhicl::ParameterSet const &pset)
Default constructor with fhicl parameters.
Definition: MCRecoPart.cxx:15
void AddParticles(const std::vector< simb::MCParticle > &mcp_v, const std::vector< simb::Origin_t > &orig_v)
Definition: MCRecoPart.cxx:124
T get(std::string const &key) const
Definition: ParameterSet.h:231
TString part[npart]
Definition: Style.C:32
Monte Carlo Simulation.
double _x_max
x-max of volume box used to determine whether to save track information
Definition: MCRecoPart.h:119
double _x_min
x-min of volume box used to determine whether to save track information
Definition: MCRecoPart.h:120
std::set< int > _pdg_list
PDG code list for which particle&#39;s trajectory within the detector is saved.
Definition: MCRecoPart.h:115
const unsigned int kINVALID_UINT
Definition: MCLimits.h:14
Namespace collecting geometry-related classes utilities.
bool InDetector(const double &x, const double &y, const double &z) const
Definition: MCRecoPart.cxx:113
unsigned int TrackToParticleIndex(const unsigned int track_id) const
Definition: MCRecoPart.h:97
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
vec_iX clear()
TPC_iterator end_TPC() const
Returns an iterator pointing after the last TPC in the detector.