LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
MCBTAlg.cxx
Go to the documentation of this file.
1 #include "MCBTAlg.h"
2 
11 
12 #include <string>
13 
14 namespace btutil {
15 
16  MCBTAlg::MCBTAlg(const std::vector<unsigned int>& g4_trackid_v,
17  const std::vector<sim::SimChannel>& simch_v)
18  {
19  Reset(g4_trackid_v, simch_v);
20  }
21 
22  void MCBTAlg::Reset(const std::vector<unsigned int>& g4_trackid_v,
23  const std::vector<sim::SimChannel>& simch_v)
24  {
25  _num_parts = 0;
26  _sum_mcq.clear();
27  _trkid_to_index.clear();
28  _event_info.clear();
29  //
30  for (auto const& id : g4_trackid_v)
31  Register(id);
32  _num_parts++;
33  ProcessSimChannel(simch_v);
34  }
35 
36  void MCBTAlg::Reset(const std::vector<std::vector<unsigned int>>& g4_trackid_v,
37  const std::vector<sim::SimChannel>& simch_v)
38  {
39  _num_parts = 0;
40  _sum_mcq.clear();
41  _trkid_to_index.clear();
42  _event_info.clear();
43  //
44  for (auto const& id : g4_trackid_v)
45  Register(id);
46  _num_parts++;
47  ProcessSimChannel(simch_v);
48  }
49 
50  void MCBTAlg::ProcessSimChannel(const std::vector<sim::SimChannel>& simch_v)
51  {
52  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
53  _sum_mcq.resize(wireReadoutGeom.Nplanes(), std::vector<double>(_num_parts, 0));
54 
55  for (auto const& sch : simch_v) {
56 
57  auto const ch = sch.Channel();
58  if (_event_info.size() <= ch) _event_info.resize(ch + 1);
59 
60  auto& ch_info = _event_info[ch];
61 
62  size_t plane = wireReadoutGeom.ChannelToWire(ch)[0].Plane;
63 
64  for (auto const& time_ide : sch.TDCIDEMap()) {
65 
66  auto const& time = time_ide.first;
67  auto const& ide_v = time_ide.second;
68 
69  auto& edep_info = ch_info[time];
70 
71  if (!edep_info.size()) edep_info.resize(_num_parts, 0);
72 
73  for (auto const& ide : ide_v) {
74 
75  size_t index = kINVALID_INDEX;
76  if (ide.trackID < (int)(_trkid_to_index.size())) { index = _trkid_to_index[ide.trackID]; }
77  if (_num_parts <= index) {
78  (*edep_info.rbegin()) += ide.numElectrons;
79  (*(_sum_mcq[plane]).rbegin()) += ide.numElectrons;
80  }
81  else {
82  edep_info[index] += ide.numElectrons;
83  _sum_mcq[plane][index] += ide.numElectrons;
84  }
85  }
86  }
87  }
88  }
89 
90  const std::vector<double>& MCBTAlg::MCQSum(const size_t plane_id) const
91  {
92  if (plane_id > _sum_mcq.size())
93  throw MCBTException(Form("Invalid plane requested: %zu", plane_id));
94  return _sum_mcq[plane_id];
95  }
96 
97  std::vector<double> MCBTAlg::MCQ(detinfo::DetectorClocksData const& clockData,
98  const WireRange_t& hit) const
99  {
100  std::vector<double> res(_num_parts, 0);
101 
102  if (_event_info.size() <= hit.ch) return res;
103 
104  auto const& ch_info = _event_info[hit.ch];
105 
106  auto itlow = ch_info.lower_bound((unsigned int)(clockData.TPCTick2TDC(hit.start)));
107  auto itup = ch_info.upper_bound((unsigned int)(clockData.TPCTick2TDC(hit.end)) + 1);
108 
109  while (itlow != ch_info.end() && itlow != itup) {
110 
111  auto const& edep_info = (*itlow).second;
112 
113  for (size_t part_index = 0; part_index < _num_parts; ++part_index)
114 
115  res[part_index] += edep_info[part_index];
116 
117  ++itlow;
118  }
119  return res;
120  }
121 
122  std::vector<double> MCBTAlg::MCQFrac(detinfo::DetectorClocksData const& clockData,
123  const WireRange_t& hit) const
124  {
125  auto res = MCQ(clockData, hit);
126  if (!res.size()) return res;
127 
128  double sum = 0;
129  for (auto const& v : res)
130  sum += v;
131  for (size_t i = 0; i < (res.size() - 1); ++i)
132  res[i] /= (sum - (*res.rbegin()));
133  (*res.rbegin()) /= sum;
134  return res;
135  }
136 
137  std::vector<double> MCBTAlg::MCQ(detinfo::DetectorClocksData const& clockData,
138  const std::vector<WireRange_t>& hit_v) const
139  {
140  std::vector<double> res(_num_parts, 0);
141  for (auto const& h : hit_v) {
142  auto tmp_res = MCQ(clockData, h);
143  for (size_t i = 0; i < res.size(); ++i)
144  res[i] += tmp_res[i];
145  }
146  return res;
147  }
148 
149  std::vector<double> MCBTAlg::MCQFrac(detinfo::DetectorClocksData const& clockData,
150  const std::vector<WireRange_t>& hit_v) const
151  {
152  auto res = MCQ(clockData, hit_v);
153  if (!res.size()) return res;
154 
155  double sum = 0;
156  for (auto const& v : res)
157  sum += v;
158  for (size_t i = 0; i < (res.size() - 1); ++i)
159  res[i] /= (sum - (*res.rbegin()));
160  (*res.rbegin()) /= sum;
161  return res;
162  }
163 
164  size_t MCBTAlg::Index(const unsigned int g4_track_id) const
165  {
166  if (g4_track_id >= _trkid_to_index.size()) return kINVALID_INDEX;
167  return _trkid_to_index[g4_track_id];
168  }
169 
170  void MCBTAlg::Register(const unsigned int& g4_track_id)
171  {
172  if (_trkid_to_index.size() <= g4_track_id)
173  _trkid_to_index.resize(g4_track_id + 1, kINVALID_INDEX);
174 
175  if (_trkid_to_index[g4_track_id] == kINVALID_INDEX) {
176  _trkid_to_index[g4_track_id] = _num_parts;
177  ++_num_parts;
178  }
179  }
180 
181  void MCBTAlg::Register(const std::vector<unsigned int>& track_id_v)
182  {
183  unsigned int max_id = 0;
184  for (auto const& id : track_id_v)
185  if (max_id < id) max_id = id;
186  if (_trkid_to_index.size() <= max_id) _trkid_to_index.resize(max_id + 1, kINVALID_INDEX);
187 
188  for (auto const& id : track_id_v) {
189 
190  if (_trkid_to_index[id] == kINVALID_INDEX)
191 
193 
194  else
195 
196  throw MCBTException(Form("Doubly used TrackID: %d", id));
197  }
198  ++_num_parts;
199  }
200 
201 }
unsigned int ch
Definition: MCBTAlg.h:34
const std::vector< double > & MCQSum(const size_t plane_id) const
Definition: MCBTAlg.cxx:90
std::vector< std::vector< double > > _sum_mcq
Definition: MCBTAlg.h:126
void Reset(const std::vector< unsigned int > &g4_trackid_v, const std::vector< sim::SimChannel > &simch_v)
Definition: MCBTAlg.cxx:22
Class def header for a class MCBTAlg.
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
std::vector<::btutil::ch_info_t > _event_info
Definition: MCBTAlg.h:124
pure virtual base interface for detector clocks
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
size_t Index(const unsigned int g4_track_id) const
Definition: MCBTAlg.cxx:164
std::vector< double > MCQ(detinfo::DetectorClocksData const &clockData, const WireRange_t &hit) const
Definition: MCBTAlg.cxx:97
static const size_t kINVALID_INDEX
Signifies invalid MCX index number.
void ProcessSimChannel(const std::vector< sim::SimChannel > &simch_v)
Definition: MCBTAlg.cxx:50
std::vector< size_t > _trkid_to_index
Definition: MCBTAlg.h:125
Definition of data types for geometry description.
Detector simulation of raw signals on wires.
Contains all timing reference information for the detector.
std::vector< double > MCQFrac(detinfo::DetectorClocksData const &clockData, const WireRange_t &hit) const
Definition: MCBTAlg.cxx:122
object containing MC truth information necessary for making RawDigits and doing back tracking ...
size_t _num_parts
Definition: MCBTAlg.h:127
Double_t sum
Definition: plot.C:31
void Register(const unsigned int &g4_track_id)
Definition: MCBTAlg.cxx:170
double TPCTick2TDC(double const tick) const
Class def header for exception classes in MCComp package.