52 auto result = std::make_unique<std::vector<sim::MCShower>>();
53 auto& mcshower = *result;
57 mcshower.reserve(shower_index_v.size());
58 std::vector<size_t> mcs_to_spart_v;
59 mcs_to_spart_v.reserve(shower_index_v.size());
61 bool daughter_stored =
false;
62 for (
size_t shower_index = 0; shower_index < shower_index_v.size(); ++shower_index) {
64 unsigned int shower_candidate = shower_index_v.at(shower_index);
65 auto const& shower_part = part_v.at(shower_candidate);
67 unsigned int mother_track = part_v.MotherTrackID(shower_candidate);
68 unsigned int ancestor_track = part_v.AncestorTrackID(shower_candidate);
72 throw cet::exception(__FUNCTION__) <<
"LOGIC ERROR: mother/ancestor track ID is invalid!";
74 MCMiniPart mother_part;
75 MCMiniPart ancestor_part;
77 unsigned int mother_index = part_v.TrackToParticleIndex(mother_track);
78 unsigned int ancestor_index = part_v.TrackToParticleIndex(ancestor_track);
81 mother_part = part_v[mother_index];
83 mother_part._track_id = mother_track;
86 ancestor_part = part_v[ancestor_index];
88 ancestor_part._track_id = ancestor_track;
90 double shower_g4_energy = shower_part._start_mom[3];
94 std::cout <<
"Found MCShower with mother energy: " << shower_g4_energy <<
" MeV";
98 if (
fDebugMode) std::cout <<
" ... below energy threshold: skipping!" << std::endl;
104 std::cout <<
" ... below # daughter particle count threshold: skipping!" << std::endl;
109 std::cout <<
" ... condition matched. Storing this MCShower..." << std::endl;
113 mcs_to_spart_v.push_back(shower_index);
117 std::cout <<
" Storage index " << mcshower.size() <<
" => Shower index " << shower_index
122 shower_prof.
Origin(shower_part._origin);
123 shower_prof.
PdgCode(shower_part._pdgcode);
124 shower_prof.
TrackID(shower_part._track_id);
125 shower_prof.
Process(shower_part._process);
135 shower_prof.
Start(MCStep(shower_part._start_vtx, shower_part._start_mom));
136 shower_prof.
End(MCStep(shower_part._end_vtx, shower_part._end_mom));
137 shower_prof.
MotherStart(MCStep(mother_part._start_vtx, mother_part._start_mom));
138 shower_prof.
MotherEnd(MCStep(mother_part._end_vtx, mother_part._end_mom));
139 shower_prof.
AncestorStart(MCStep(ancestor_part._start_vtx, ancestor_part._start_mom));
140 shower_prof.
AncestorEnd(MCStep(ancestor_part._end_vtx, ancestor_part._end_mom));
143 std::vector<unsigned int> daughter_track_id;
148 daughter_track_id.push_back(part_v.at(index)._track_id);
152 if (!daughter_stored && daughter_track_id.size() > 1) daughter_stored =
true;
154 mcshower.push_back(shower_prof);
158 std::cout <<
" Found " << mcshower.size()
159 <<
" MCShowers. Now computing DetProfile position..." << std::endl;
164 std::vector<TLorentzVector> mcs_daughter_vtx_v(
168 std::vector<TLorentzVector> mcs_daughter_mom_v(mcshower.size(), TLorentzVector());
170 std::vector<std::vector<double>> plane_charge_v(mcshower.size(), std::vector<double>(3, 0));
171 std::vector<std::vector<double>> plane_dqdx_v(mcshower.size(), std::vector<double>(3, 0));
174 std::vector<double> mcs_daughter_dedx_v(mcshower.size(), 0);
175 std::vector<double> mcs_daughter_dedxRAD_v(mcshower.size(), 0);
176 std::vector<TVector3> mcs_daughter_dir_v(mcshower.size(), TVector3());
178 for (
size_t mcs_index = 0; mcs_index < mcshower.size(); ++mcs_index) {
180 auto& mcs_daughter_vtx = mcs_daughter_vtx_v[mcs_index];
181 auto& mcs_daughter_mom = mcs_daughter_mom_v[mcs_index];
182 auto& mcs_daughter_dedx = mcs_daughter_dedx_v[mcs_index];
183 auto& mcs_daughter_dedxRAD = mcs_daughter_dedxRAD_v[mcs_index];
184 auto& mcs_daughter_dir = mcs_daughter_dir_v[mcs_index];
185 auto& plane_charge = plane_charge_v[mcs_index];
186 auto& plane_dqdx = plane_dqdx_v[mcs_index];
188 for (
auto const& daughter_trk_id : mcshower[mcs_index].DaughterTrackID()) {
190 auto const daughter_part_index = part_v.TrackToParticleIndex(daughter_trk_id);
192 auto const& daughter_part = part_v[daughter_part_index];
194 auto const daughter_edep_index = edep_v.TrackToEdepIndex(daughter_trk_id);
196 if (daughter_edep_index < 0)
continue;
198 auto const& daughter_edep = edep_v.GetEdepArrayAt(daughter_edep_index);
200 if (!(daughter_edep.size()))
continue;
204 for (
auto const&
edep : daughter_edep) {
206 double dist = sqrt(pow(
edep.pos._x - daughter_part._start_vtx[0], 2) +
207 pow(
edep.pos._y - daughter_part._start_vtx[1], 2) +
208 pow(
edep.pos._z - daughter_part._start_vtx[2], 2));
210 if (dist < min_dist) {
212 mcs_daughter_vtx[0] =
edep.pos._x;
213 mcs_daughter_vtx[1] =
edep.pos._y;
214 mcs_daughter_vtx[2] =
edep.pos._z;
215 mcs_daughter_vtx[3] = (dist / 100. / 2.998e8) * 1.e9 + daughter_part._start_vtx[3];
218 if (!daughter_stored) {
220 std::vector<double> shower_dir(3, 0);
221 shower_dir[0] = mcshower[mcs_index].Start().Px();
222 shower_dir[1] = mcshower[mcs_index].Start().Py();
223 shower_dir[2] = mcshower[mcs_index].Start().Pz();
224 double magnitude = 0;
225 for (
size_t i = 0; i < 3; ++i)
226 magnitude += pow(shower_dir[i], 2);
228 magnitude = sqrt(magnitude);
230 if (magnitude > 1.
e-10) {
234 for (
auto& v : shower_dir)
237 for (
auto const&
edep : daughter_edep) {
238 std::vector<double> shower_dep_dir(3, 0);
239 shower_dep_dir[0] =
edep.pos._x - mcshower[mcs_index].Start().X();
240 shower_dep_dir[1] =
edep.pos._y - mcshower[mcs_index].Start().Y();
241 shower_dep_dir[2] =
edep.pos._z - mcshower[mcs_index].Start().Z();
243 double dist = sqrt(pow(shower_dep_dir[0], 2) + pow(shower_dep_dir[1], 2) +
244 pow(shower_dep_dir[2], 2));
245 for (
auto& v : shower_dep_dir)
249 acos(shower_dep_dir[0] * shower_dir[0] + shower_dep_dir[1] * shower_dir[1] +
250 shower_dep_dir[2] * shower_dir[2]) /
253 if (dist < min_dist && angle < 10) {
256 mcs_daughter_vtx[0] =
edep.pos._x;
257 mcs_daughter_vtx[1] =
edep.pos._y;
258 mcs_daughter_vtx[2] =
edep.pos._z;
259 mcs_daughter_vtx[3] =
260 (dist / 100. / 2.998e8) * 1.e9 + mcshower[mcs_index].Start().T();
269 std::vector<double> mom(3, 0);
270 for (
auto const& daughter_trk_id : mcshower[mcs_index].DaughterTrackID()) {
277 auto const daughter_edep_index = edep_v.TrackToEdepIndex(daughter_trk_id);
279 if (daughter_edep_index < 0)
continue;
281 auto const& daughter_edep = edep_v.GetEdepArrayAt(daughter_edep_index);
283 if (!(daughter_edep.size()))
continue;
286 for (
auto const&
edep : daughter_edep) {
289 mom[0] =
edep.pos._x - mcs_daughter_vtx[0];
290 mom[1] =
edep.pos._y - mcs_daughter_vtx[1];
291 mom[2] =
edep.pos._z - mcs_daughter_vtx[2];
294 double magnitude = sqrt(pow(mom.at(0), 2) + pow(mom.at(1), 2) + pow(mom.at(2), 2));
298 for (
auto const& pid_energy :
edep.deps) {
300 energy += pid_energy.energy;
303 if (magnitude > 1.
e-10) {
304 mom.at(0) = mom.at(0) * energy / magnitude;
305 mom.at(1) = mom.at(1) * energy / magnitude;
306 mom.at(2) = mom.at(2) * energy / magnitude;
307 mcs_daughter_mom[0] += mom.at(0);
308 mcs_daughter_mom[1] += mom.at(1);
309 mcs_daughter_mom[2] += mom.at(2);
314 if (sqrt(pow(
edep.pos._x - mcs_daughter_vtx[0], 2) +
315 pow(
edep.pos._y - mcs_daughter_vtx[1], 2) +
316 pow(
edep.pos._z - mcs_daughter_vtx[2], 2)) < 2.4 &&
317 magnitude > 1.e-10) {
319 mcs_daughter_dir[0] += mom.at(0);
320 mcs_daughter_dir[1] += mom.at(1);
321 mcs_daughter_dir[2] += mom.at(2);
327 mcs_daughter_dedxRAD +=
E;
329 mcs_daughter_mom[3] +=
energy;
333 auto q_i = pindex.find(
pid);
334 if (q_i != pindex.end())
335 plane_charge[
pid.Plane] += (
double)(
edep.deps[pindex[
pid]].charge);
340 mcs_daughter_dedxRAD /= 2.4;
342 for (
auto const& daughter_trk_id : mcshower[mcs_index].DaughterTrackID()) {
349 auto const daughter_edep_index = edep_v.TrackToEdepIndex(daughter_trk_id);
351 if (daughter_edep_index < 0)
continue;
353 auto const& daughter_edep = edep_v.GetEdepArrayAt(daughter_edep_index);
355 if (!(daughter_edep.size()))
continue;
357 for (
auto const&
edep : daughter_edep) {
367 double p_mag = sqrt(pow(mcs_daughter_dir[0], 2) + pow(mcs_daughter_dir[1], 2) +
368 pow(mcs_daughter_dir[2], 2));
369 double a = 0, b = 0, c = 0,
d = 0;
370 if (p_mag > 1.
e-10) {
371 a = mcs_daughter_dir[0] / p_mag;
372 b = mcs_daughter_dir[1] / p_mag;
373 c = mcs_daughter_dir[2] / p_mag;
374 d = -1 * (a * mcs_daughter_vtx[0] + b * mcs_daughter_vtx[1] + c * mcs_daughter_vtx[2]);
377 mcs_daughter_dedx += 0;
381 if ((a *
edep.pos._x + b *
edep.pos._y + c *
edep.pos._z +
d) /
382 sqrt(pow(a, 2) + pow(b, 2) + pow(c, 2)) <
385 sqrt(pow(a, 2) + pow(b, 2) + pow(c, 2)) >
391 for (
auto const& pid_energy :
edep.deps) {
393 E += pid_energy.energy;
396 if (N > 0) { E /= N; }
401 mcs_daughter_dedx +=
E;
405 auto q_i = pindex.find(
pid);
406 if (q_i != pindex.end())
407 plane_dqdx[
pid.Plane] += (
double)(
edep.deps[pindex[
pid]].charge);
411 mcs_daughter_dedx /= 2.4;
412 plane_dqdx.at(0) /= 2.4;
413 plane_dqdx.at(1) /= 2.4;
414 plane_dqdx.at(2) /= 2.4;
419 std::cout <<
" Found " << mcshower.size() <<
" MCShowers. Now storing..." << std::endl;
422 for (
size_t mcs_index = 0; mcs_index < mcshower.size(); ++mcs_index) {
424 auto& daughter_vtx = mcs_daughter_vtx_v[mcs_index];
425 auto& daughter_mom = mcs_daughter_mom_v[mcs_index];
426 auto& daughter_dedx = mcs_daughter_dedx_v[mcs_index];
427 auto& daughter_dedxRAD = mcs_daughter_dedxRAD_v[mcs_index];
428 auto& daughter_dir = mcs_daughter_dir_v[mcs_index];
429 auto& plane_charge = plane_charge_v[mcs_index];
430 auto& plane_dqdx = plane_dqdx_v[mcs_index];
433 sqrt(pow(daughter_mom[0], 2) + pow(daughter_mom[1], 2) + pow(daughter_mom[2], 2));
435 if (daughter_mom[3] > 1.
e-10) {
436 daughter_mom[0] *= daughter_mom[3] / magnitude;
437 daughter_mom[1] *= daughter_mom[3] / magnitude;
438 daughter_mom[2] *= daughter_mom[3] / magnitude;
441 for (
size_t i = 0; i < 4; ++i)
444 mcshower.at(mcs_index).DetProfile(MCStep(daughter_vtx, daughter_mom));
445 mcshower.at(mcs_index).Charge(plane_charge);
446 mcshower.at(mcs_index).dQdx(plane_dqdx);
447 mcshower.at(mcs_index).dEdx(daughter_dedx);
448 mcshower.at(mcs_index).dEdxRAD(daughter_dedxRAD);
449 mcshower.at(mcs_index).StartDir(daughter_dir);
454 for (
auto const& prof : mcshower) {
458 << Form(
" Shower particle: PDG=%d : Track ID=%d Start @ (%g,%g,%g,%g) with Momentum " 471 << Form(
" Mother particle: PDG=%d : Track ID=%d Start @ (%g,%g,%g,%g) with Momentum " 473 prof.MotherPdgCode(),
474 prof.MotherTrackID(),
475 prof.MotherStart().X(),
476 prof.MotherStart().Y(),
477 prof.MotherStart().Z(),
478 prof.MotherStart().T(),
479 prof.MotherStart().Px(),
480 prof.MotherStart().Py(),
481 prof.MotherStart().Pz(),
482 prof.MotherStart().E())
484 << Form(
" Ancestor particle: PDG=%d : Track ID=%d Start @ (%g,%g,%g,%g) with Momentum " 486 prof.AncestorPdgCode(),
487 prof.AncestorTrackID(),
488 prof.AncestorStart().X(),
489 prof.AncestorStart().Y(),
490 prof.AncestorStart().Z(),
491 prof.AncestorStart().T(),
492 prof.AncestorStart().Px(),
493 prof.AncestorStart().Py(),
494 prof.AncestorStart().Pz(),
495 prof.AncestorStart().E())
497 << Form(
" ... with %zu daughters: Start @ (%g,%g,%g,%g) with Momentum (%g,%g,%g,%g)",
498 prof.DaughterTrackID().size(),
499 prof.DetProfile().X(),
500 prof.DetProfile().Y(),
501 prof.DetProfile().Z(),
502 prof.DetProfile().T(),
503 prof.DetProfile().Px(),
504 prof.DetProfile().Py(),
505 prof.DetProfile().Pz(),
506 prof.DetProfile().E())
508 <<
" Charge per plane: ";
509 size_t const nPlanes = prof.Charge().size();
510 for (
size_t i = 0; i < nPlanes; ++i) {
512 std::cout <<
" | Plane " << i << std::flush;
513 std::cout <<
" ... Q = " << prof.Charge(i) << std::flush;
515 std::cout << std::endl << std::endl;
const double kINVALID_DOUBLE
const MCStep & End() const
unsigned int TrackID() const
const std::vector< unsigned int > ShowerMothers() const
MCShowerRecoPart fPartAlg
std::map< geo::PlaneID, size_t > createPlaneIndexMap()
unsigned int fMinNumDaughters
const std::vector< unsigned int > & DaughterTrackID() const
simb::Origin_t Origin() const
int MotherPdgCode() const
const std::string & AncestorProcess() const
const MCStep & AncestorStart() const
const std::string & MotherProcess() const
const std::vector< unsigned int > & ShowerDaughters(const unsigned int shower_id) const
unsigned int AncestorTrackID() const
const MCStep & AncestorEnd() const
const MCStep & Start() const
const MCStep & MotherEnd() const
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
const unsigned int kINVALID_UINT
unsigned int MotherTrackID() const
const std::string & Process() const
const MCStep & MotherStart() const
int AncestorPdgCode() const
Namespace collecting geometry-related classes utilities.
cet::coded_exception< error, detail::translate > exception
void ConstructShower(const MCRecoPart &part_v)
Main function to read-in data and fill variables in this algorithm to reconstruct MC shower...