30{
31 if (!evt) { B2ERROR("HerwigFragHelper: null GenEvent passed"); return -1; }
32
33 const int nparticles = evt->particles_size();
34 if (nparticles <= 0) { B2ERROR("HepMC event has no particles"); return -1; }
35
36 const int event_offset = graph.
size();
37 int n_added = 0;
38 std::unordered_map<int, int> hash_index_map;
39
40
41
42
43
44 {
45 int mapped_count = 0;
46 auto it = evt->particles_begin();
47 for (int i = 0; i < nparticles; ++i, ++it) {
48 const int st = (*it)->status();
49 const int pdg = (*it)->pdg_id();
50 if (st == 4 || st < 1) continue;
51
52
53 if (pdg >= 81 && pdg <= 99) continue;
54 hash_index_map[(*it)->barcode()] = event_offset + mapped_count;
55 mapped_count++;
56 }
57 }
58
59 const double len_conv = HepMC::Units::conversion_factor(evt->length_unit(), HepMC::Units::CM);
60 const double mom_conv = HepMC::Units::conversion_factor(evt->momentum_unit(), HepMC::Units::GEV);
61
62
63 auto rp = evt->particles_begin();
64 for (int i = 0; i < nparticles; ++i, ++rp) {
65 auto* dv = (*rp)->end_vertex();
66 auto* pv = (*rp)->production_vertex();
67 const int status = (*rp)->status();
68 const int pdg_code = (*rp)->pdg_id();
69
70 if (status == 4 || status < 1) continue;
71 if (pdg_code >= 81 && pdg_code <= 99) continue;
72
74 MCParticleGraph::GraphParticle& p = graph[event_offset + n_added];
75
76 p.setPDG(pdg_code);
77
78 auto mom = (*rp)->momentum();
79 p.setMomentum(ROOT::Math::XYZVector(
83 p.setMass((*rp)->generated_mass() * mom_conv);
84
85 const double e2 = p.getMomentum().Mag2() + p.getMass() * p.getMass();
86 p.setEnergy(e2 > 0.0 ? std::sqrt(e2) : 0.0);
87
88 if (pv) {
89 auto pos = pv->position();
90 p.setProductionVertex(ROOT::Math::XYZVector(pos.x(), pos.y(), pos.z()) * len_conv *
Unit::cm);
92 p.setValidVertex(true);
93 }
94
96
97
98
99 const bool isFinal = !dv && status == 1;
102
103
104
105
106
107
108 if (pv) {
109 auto par_it = pv->particles_begin(HepMC::parents);
110 if (par_it != pv->particles_end(HepMC::parents)) {
111 auto found = hash_index_map.find((*par_it)->barcode());
112 if (found != hash_index_map.end()) {
113 const int pidx = found->second;
114 const int sidx = event_offset + n_added;
115 if (
static_cast<size_t>(pidx) < graph.
size() && pidx != sidx)
116 p.comesFrom(graph[pidx]);
117 } else {
118 bool linked = false;
119 for (int qidx : quarkIndices)
120 if (
static_cast<size_t>(qidx) < graph.
size() && graph[qidx].getPDG() == pdg_code)
121 { p.comesFrom(graph[qidx]); linked = true; break; }
122 if (!linked && !quarkIndices.empty())
123 p.comesFrom(graph[quarkIndices[0]]);
124 }
125 }
126 }
127
128 n_added++;
129 }
130
131
132 {
133 double fsE = 0, fsPx = 0, fsPy = 0, fsPz = 0;
134 for (auto pit = evt->particles_begin(); pit != evt->particles_end(); ++pit) {
135 if ((*pit)->status() != 1) continue;
136 const auto& m2 = (*pit)->momentum();
137 fsE += m2.t() * mom_conv; fsPx += m2.x() * mom_conv;
138 fsPy += m2.y() * mom_conv; fsPz += m2.z() * mom_conv;
139 }
140 const double fsM2 = fsE * fsE - fsPx * fsPx - fsPy * fsPy - fsPz * fsPz;
141 B2DEBUG(10, "Post-Herwig final-state sum: E=" << fsE
142 << " M_inv=" << (fsM2 > 0.0 ? std::sqrt(fsM2) : 0.0));
143 }
144
145 B2DEBUG(10, "Added " << n_added << " particles to graph");
146 return n_added;
147}
static const double speedOfLight
[cm/ns]
size_t size() const
Return the number of particles in the graph.
@ c_IsFSRPhoton
bit 7: Particle is from final state radiation
@ c_PrimaryParticle
bit 0: Particle is primary particle.
@ c_StableInGenerator
bit 1: Particle is stable, i.e., not decaying in the generator.
static const double cm
Standard units with the value = 1.
static const double GeV
Standard of [energy, momentum, mass].
GraphParticle & addParticle()
Add new particle to the graph.
1.13.2