Supera/MCParticleList.cxx at develop · DeepLearnPhysics/Supera · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
#ifndef SUPERA_MCPARTICLELIST_CXX
#define SUPERA_MCPARTICLELIST_CXX

#include "MCParticleList.h"
namespace supera {

  void MCParticleList::Update(const std::vector<simb::MCParticle>& larmcp_v, int run, int event)
  {
    if (run == _run && event == _event) return;

    _run = run;
    _event = event;

    _trackid_v.resize(larmcp_v.size());
    _pdgcode_v.resize(larmcp_v.size());
    _parent_index_v.resize(larmcp_v.size());
    _parent_trackid_v.resize(larmcp_v.size());
    _parent_pdg_v.resize(larmcp_v.size());
    _ancestor_trackid_v.resize(larmcp_v.size());
    _ancestor_index_v.resize(larmcp_v.size());

    for (size_t i = 0; i < larmcp_v.size(); ++i)
      _trackid_v[i] = _pdgcode_v[i] = _parent_index_v[i] = _parent_trackid_v[i] = _parent_pdg_v[i] =
        _ancestor_trackid_v[i] = _ancestor_index_v[i] = -1;

    _trackid2index.resize(std::max(_trackid2index.size(), larmcp_v.size()));
    for (auto& v : _trackid2index)
      v = -1;

    for (size_t index = 0; index < larmcp_v.size(); ++index) {
      auto const& mcpart = larmcp_v[index];
      _trackid_v[index] = abs(mcpart.TrackId());
      _pdgcode_v[index] = abs(mcpart.PdgCode());
      _parent_trackid_v[index] = mcpart.Mother();
      if (mcpart.TrackId() >= ((int)(_trackid2index.size())))
        _trackid2index.resize(mcpart.TrackId() + 1, -1);
      _trackid2index[mcpart.TrackId()] = index;
    }
    // Parent/Ancestor info
    for (size_t index = 0; index < larmcp_v.size(); ++index) {
      auto const& mcpart = larmcp_v[index];
      int mother_id = mcpart.Mother();
      int mother_index = -1;
      if (mother_id == 0) mother_id = abs(mcpart.TrackId());
      if (mother_id < ((int)(_trackid2index.size()))) {
        mother_index = _trackid2index[mother_id];
        if (mother_index >= 0) {
          _parent_pdg_v[index] = larmcp_v[mother_index].PdgCode();
          _parent_index_v[index] = mother_index;
        }
      }

      int subject_track_id = abs(mcpart.TrackId());
      int parent_track_id = abs(mcpart.Mother());
      int ancestor_index = -1;
      int ancestor_track_id = -1;
      while (1) {
        if ((size_t)(parent_track_id) >= _trackid2index.size()) break;
        if (parent_track_id == 0 || parent_track_id == subject_track_id) {
          ancestor_index = _trackid2index[subject_track_id];
          ancestor_track_id = subject_track_id;
          break;
        }
        auto const& parent_index = _trackid2index[parent_track_id];
        if (parent_index < 0) break;
        auto const& parent = larmcp_v[parent_index];
        subject_track_id = abs(parent.TrackId());
        parent_track_id = abs(parent.Mother());
      }
      _ancestor_index_v[index] = ancestor_index;
      _ancestor_trackid_v[index] = ancestor_track_id;
    }
  }
}
#endif