GCC Code Coverage Report

Directory: ./
File: src/iguana/algorithms/physics/DihadronKinematics/Algorithm.cc
Date: 2025-11-25 17:57:04
Coverage Exec Excl Total
Lines: 93.8% 121 0 129
Functions: 100.0% 5 0 5
Branches: 54.1% 93 0 172
Line Branch Exec Source
1 #include "Algorithm.h"
2 #include "TypeDefs.h"
3
4 #include "iguana/algorithms/physics/Tools.h"
5 #include <Math/Boost.h>
6
7 namespace iguana::physics {
8
9 REGISTER_IGUANA_ALGORITHM(DihadronKinematics, "physics::DihadronKinematics");
10
11 2 void DihadronKinematics::Start(hipo::banklist& banks)
12 {
13 // parse config file
14 2 ParseYAMLConfig();
15
2/4
✓ Branch 4 → 5 taken 2 times.
✗ Branch 4 → 130 not taken.
✓ Branch 5 → 6 taken 2 times.
✗ Branch 5 → 124 not taken.
 4 o_particle_bank = GetOptionScalar<std::string>("particle_bank");
16
2/4
✓ Branch 19 → 20 taken 2 times.
✗ Branch 19 → 138 not taken.
✓ Branch 20 → 21 taken 2 times.
✗ Branch 20 → 132 not taken.
 4 o_hadron_a_pdgs = GetOptionSet<int>("hadron_a_list");
17
2/4
✓ Branch 30 → 31 taken 2 times.
✗ Branch 30 → 146 not taken.
✓ Branch 31 → 32 taken 2 times.
✗ Branch 31 → 140 not taken.
 4 o_hadron_b_pdgs = GetOptionSet<int>("hadron_b_list");
18
2/4
✓ Branch 41 → 42 taken 2 times.
✗ Branch 41 → 154 not taken.
✓ Branch 42 → 43 taken 2 times.
✗ Branch 42 → 148 not taken.
 4 o_phi_r_method = GetOptionScalar<std::string>("phi_r_method");
19
2/4
✓ Branch 56 → 57 taken 2 times.
✗ Branch 56 → 162 not taken.
✓ Branch 57 → 58 taken 2 times.
✗ Branch 57 → 156 not taken.
 4 o_theta_method = GetOptionScalar<std::string>("theta_method");
20
21 // check phiR method
22
1/2
✓ Branch 70 → 71 taken 2 times.
✗ Branch 70 → 72 not taken.
 2 if(o_phi_r_method == "RT_via_covariant_kT")
23 2 m_phi_r_method = e_RT_via_covariant_kT;
24 else
25 ✗ throw std::runtime_error(fmt::format("unknown phi_r_method: {:?}", o_phi_r_method));
26
27 // check theta method
28
1/2
✓ Branch 71 → 83 taken 2 times.
✗ Branch 71 → 86 not taken.
 2 if(o_theta_method == "hadron_a")
29 2 m_theta_method = e_hadron_a;
30 else
31 ✗ throw std::runtime_error(fmt::format("unknown theta_method: {:?}", o_theta_method));
32
33 // get bank indices
34 2 b_particle = GetBankIndex(banks, o_particle_bank);
35
2/4
✓ Branch 85 → 97 taken 2 times.
✗ Branch 85 → 180 not taken.
✗ Branch 97 → 98 not taken.
✓ Branch 97 → 100 taken 2 times.
 4 b_inc_kin = GetBankIndex(banks, "physics::InclusiveKinematics");
36
37 // create the output bank
38
1/2
✓ Branch 103 → 104 taken 2 times.
✗ Branch 103 → 186 not taken.
 2 auto result_schema = CreateBank(banks, b_result, GetClassName());
39
1/2
✓ Branch 109 → 110 taken 2 times.
✗ Branch 109 → 192 not taken.
 2 i_pindex_a = result_schema.getEntryOrder("pindex_a");
40
1/2
✓ Branch 110 → 111 taken 2 times.
✗ Branch 110 → 192 not taken.
 2 i_pindex_b = result_schema.getEntryOrder("pindex_b");
41
1/2
✓ Branch 111 → 112 taken 2 times.
✗ Branch 111 → 192 not taken.
 2 i_pdg_a = result_schema.getEntryOrder("pdg_a");
42
1/2
✓ Branch 112 → 113 taken 2 times.
✗ Branch 112 → 192 not taken.
 2 i_pdg_b = result_schema.getEntryOrder("pdg_b");
43
1/2
✓ Branch 113 → 114 taken 2 times.
✗ Branch 113 → 192 not taken.
 2 i_Mh = result_schema.getEntryOrder("Mh");
44
1/2
✓ Branch 114 → 115 taken 2 times.
✗ Branch 114 → 192 not taken.
 2 i_z = result_schema.getEntryOrder("z");
45
1/2
✓ Branch 115 → 116 taken 2 times.
✗ Branch 115 → 192 not taken.
 2 i_PhPerp = result_schema.getEntryOrder("PhPerp");
46
1/2
✓ Branch 116 → 117 taken 2 times.
✗ Branch 116 → 192 not taken.
 2 i_MX2 = result_schema.getEntryOrder("MX2");
47
1/2
✓ Branch 117 → 118 taken 2 times.
✗ Branch 117 → 192 not taken.
 2 i_xF = result_schema.getEntryOrder("xF");
48
1/2
✓ Branch 118 → 119 taken 2 times.
✗ Branch 118 → 192 not taken.
 2 i_yB = result_schema.getEntryOrder("yB");
49
1/2
✓ Branch 119 → 120 taken 2 times.
✗ Branch 119 → 192 not taken.
 2 i_phiH = result_schema.getEntryOrder("phiH");
50
1/2
✓ Branch 120 → 121 taken 2 times.
✗ Branch 120 → 192 not taken.
 2 i_phiR = result_schema.getEntryOrder("phiR");
51
1/2
✓ Branch 121 → 122 taken 2 times.
✗ Branch 121 → 192 not taken.
 2 i_theta = result_schema.getEntryOrder("theta");
52 2 }
53
54 ///////////////////////////////////////////////////////////////////////////////
55
56 286 bool DihadronKinematics::Run(hipo::banklist& banks) const
57 {
58
1/2
✓ Branch 7 → 8 taken 286 times.
✗ Branch 7 → 19 not taken.
 572 return Run(
59
1/2
✓ Branch 6 → 7 taken 286 times.
✗ Branch 6 → 19 not taken.
 286 GetBank(banks, b_particle, o_particle_bank),
60
3/8
✓ Branch 4 → 5 taken 286 times.
✗ Branch 4 → 25 not taken.
✓ Branch 5 → 6 taken 286 times.
✗ Branch 5 → 19 not taken.
✗ Branch 13 → 14 not taken.
✓ Branch 13 → 16 taken 286 times.
✗ Branch 25 → 26 not taken.
✗ Branch 25 → 28 not taken.
 572 GetBank(banks, b_inc_kin, "physics::InclusiveKinematics"),
61
1/2
✓ Branch 3 → 4 taken 286 times.
✗ Branch 3 → 25 not taken.
 572 GetBank(banks, b_result, GetClassName()));
62 }
63
64 286 bool DihadronKinematics::Run(
65 hipo::bank const& particle_bank,
66 hipo::bank const& inc_kin_bank,
67 hipo::bank& result_bank) const
68 {
69 286 result_bank.reset(); // IMPORTANT: always first `reset` the created bank(s)
70
1/2
✓ Branch 6 → 7 taken 286 times.
✗ Branch 6 → 136 not taken.
 572 ShowBank(particle_bank, Logger::Header("INPUT PARTICLES"));
71
72
2/4
✓ Branch 13 → 14 taken 286 times.
✗ Branch 13 → 16 not taken.
✗ Branch 15 → 16 not taken.
✓ Branch 15 → 18 taken 286 times.
 286 if(particle_bank.getRowList().empty() || inc_kin_bank.getRowList().empty()) {
73 ✗ m_log->Debug("skip this event, since not all required banks have entries");
74 ✗ return false;
75 }
76
77 // get beam and target momenta
78 // FIXME: makes some assumptions about the beam; this should be generalized...
79 ROOT::Math::PxPyPzMVector p_beam(
80 0.0,
81 0.0,
82 286 inc_kin_bank.getDouble("beamPz", 0),
83 286 particle::mass.at(particle::electron));
84 ROOT::Math::PxPyPzMVector p_target(
85 0.0,
86 0.0,
87 0.0,
88 286 inc_kin_bank.getDouble("targetM", 0));
89
90 // get virtual photon momentum
91 ROOT::Math::PxPyPzEVector p_q(
92 286 inc_kin_bank.getDouble("qx", 0),
93 286 inc_kin_bank.getDouble("qy", 0),
94 286 inc_kin_bank.getDouble("qz", 0),
95 286 inc_kin_bank.getDouble("qE", 0));
96
97 // get additional inclusive variables
98 286 auto x = inc_kin_bank.getDouble("x", 0);
99 286 auto W = inc_kin_bank.getDouble("W", 0);
100
101 // boosts
102 286 ROOT::Math::Boost boost__qp((p_q + p_target).BoostToCM()); // CoM frame of target and virtual photon
103 286 ROOT::Math::Boost boost__breit((p_q + 2 * x * p_target).BoostToCM()); // Breit frame
104 286 auto p_q__qp = boost__qp(p_q);
105 286 auto p_q__breit = boost__breit(p_q);
106
107 // build list of dihadron rows (pindices)
108 286 auto dih_rows = PairHadrons(particle_bank);
109
110 // loop over dihadrons
111
1/2
✓ Branch 42 → 43 taken 286 times.
✗ Branch 42 → 149 not taken.
 286 result_bank.setRows(dih_rows.size());
112 int dih_row = 0;
113
2/2
✓ Branch 121 → 44 taken 61 times.
✓ Branch 121 → 122 taken 286 times.
 347 for(auto const& [row_a, row_b] : dih_rows) {
114
115 // get hadron momenta
116 61 Hadron had_a{.row = row_a};
117 61 Hadron had_b{.row = row_b};
118
2/2
✓ Branch 53 → 45 taken 122 times.
✓ Branch 53 → 54 taken 61 times.
 183 for(auto& had : {&had_a, &had_b}) {
119 122 had->pdg = particle_bank.getInt("pid", had->row);
120 122 had->p = ROOT::Math::PxPyPzMVector(
121
1/2
✗ Branch 50 → 51 not taken.
✓ Branch 50 → 52 taken 122 times.
 122 particle_bank.getFloat("px", had->row),
122 122 particle_bank.getFloat("py", had->row),
123 122 particle_bank.getFloat("pz", had->row),
124
1/2
✓ Branch 46 → 47 taken 122 times.
✗ Branch 46 → 149 not taken.
 244 particle::mass.at(static_cast<particle::PDG>(had->pdg)));
125 }
126
127 // calculate dihadron momenta and boosts
128 auto p_Ph = had_a.p + had_b.p;
129
1/2
✓ Branch 56 → 57 taken 61 times.
✗ Branch 56 → 149 not taken.
 61 auto p_Ph__qp = boost__qp(p_Ph);
130
1/2
✓ Branch 57 → 58 taken 61 times.
✗ Branch 57 → 149 not taken.
 61 auto p_Ph__breit = boost__breit(p_Ph);
131 61 ROOT::Math::Boost boost__dih(p_Ph.BoostToCM()); // CoM frame of dihadron
132
133 // calculate z
134
1/2
✓ Branch 63 → 64 taken 61 times.
✗ Branch 63 → 149 not taken.
 61 double z = p_target.Dot(p_Ph) / p_target.Dot(p_q);
135
136 // calculate PhPerp
137
2/4
✓ Branch 63 → 64 taken 61 times.
✗ Branch 63 → 149 not taken.
✓ Branch 64 → 65 taken 61 times.
✗ Branch 64 → 66 not taken.
 61 auto opt_PhPerp = tools::RejectVector(p_Ph.Vect(), p_q.Vect());
138
1/2
✓ Branch 64 → 65 taken 61 times.
✗ Branch 64 → 66 not taken.
 61 double PhPerp = opt_PhPerp.has_value() ? opt_PhPerp.value().R() : tools::UNDEF;
139
140 // calculate Mh
141 double Mh = p_Ph.M();
142
143 // calculate MX2
144
1/2
✓ Branch 73 → 74 taken 61 times.
✗ Branch 73 → 149 not taken.
 61 double MX2 = (p_target + p_q - p_Ph).M2();
145
146 // calculate xF
147
1/2
✓ Branch 73 → 74 taken 61 times.
✗ Branch 73 → 149 not taken.
 61 double xF = 2 * p_Ph__qp.Vect().Dot(p_q__qp.Vect()) / (W * p_q__qp.Vect().R());
148
149 // calculate yB
150
2/4
✓ Branch 73 → 74 taken 61 times.
✗ Branch 73 → 149 not taken.
✓ Branch 76 → 77 taken 61 times.
✗ Branch 76 → 149 not taken.
 122 double yB = tools::ParticleRapidity(p_Ph__breit, p_q__breit.Vect()).value_or(tools::UNDEF);
151
152 // calculate phiH
153
1/2
✓ Branch 76 → 77 taken 61 times.
✗ Branch 76 → 149 not taken.
 61 double phiH = tools::PlaneAngle(
154 61 p_q.Vect(),
155 61 p_beam.Vect(),
156 61 p_q.Vect(),
157 61 p_Ph.Vect())
158 .value_or(tools::UNDEF);
159
160 // calculate PhiR
161 double phiR = tools::UNDEF;
162
1/2
✗ Branch 79 → 80 not taken.
✓ Branch 79 → 81 taken 61 times.
 61 switch(m_phi_r_method) {
163 61 case e_RT_via_covariant_kT: {
164
2/2
✓ Branch 94 → 82 taken 122 times.
✓ Branch 94 → 95 taken 61 times.
 183 for(auto& had : {&had_a, &had_b}) {
165
1/2
✓ Branch 84 → 85 taken 122 times.
✗ Branch 84 → 149 not taken.
 122 had->z = p_target.Dot(had->p) / p_target.Dot(p_q);
166
1/2
✓ Branch 84 → 85 taken 122 times.
✗ Branch 84 → 149 not taken.
 244 had->p_perp = tools::RejectVector(had->p.Vect(), p_q.Vect());
167 }
168
2/4
✗ Branch 95 → 80 not taken.
✓ Branch 95 → 96 taken 61 times.
✗ Branch 96 → 80 not taken.
✓ Branch 96 → 97 taken 61 times.
 61 if(had_a.p_perp.has_value() && had_b.p_perp.has_value()) {
169
1/2
✓ Branch 97 → 98 taken 61 times.
✗ Branch 97 → 149 not taken.
 61 auto RT = (had_b.z * had_a.p_perp.value() - had_a.z * had_b.p_perp.value()) / (had_a.z + had_b.z);
170
1/2
✓ Branch 97 → 98 taken 61 times.
✗ Branch 97 → 149 not taken.
 61 phiR = tools::PlaneAngle(
171 61 p_q.Vect(),
172 61 p_beam.Vect(),
173 61 p_q.Vect(),
174 RT)
175 .value_or(tools::UNDEF);
176 }
177 break;
178 }
179 }
180
181 // calculate theta
182 double theta = tools::UNDEF;
183
1/2
✓ Branch 101 → 102 taken 61 times.
✗ Branch 101 → 107 not taken.
 61 switch(m_theta_method) {
184
1/2
✓ Branch 102 → 103 taken 61 times.
✗ Branch 102 → 142 not taken.
 61 case e_hadron_a: {
185
1/2
✓ Branch 103 → 104 taken 61 times.
✗ Branch 103 → 142 not taken.
 61 theta = tools::VectorAngle(
186
1/2
✓ Branch 102 → 103 taken 61 times.
✗ Branch 102 → 142 not taken.
 61 boost__dih(had_a.p).Vect(),
187 ✗ p_Ph.Vect())
188 .value_or(tools::UNDEF);
189 61 break;
190 }
191 }
192
193 61 result_bank.putShort(i_pindex_a, dih_row, static_cast<int16_t>(had_a.row));
194 61 result_bank.putShort(i_pindex_b, dih_row, static_cast<int16_t>(had_b.row));
195 61 result_bank.putInt(i_pdg_a, dih_row, had_a.pdg);
196 61 result_bank.putInt(i_pdg_b, dih_row, had_b.pdg);
197 61 result_bank.putDouble(i_Mh, dih_row, Mh);
198 61 result_bank.putDouble(i_z, dih_row, z);
199 61 result_bank.putDouble(i_PhPerp, dih_row, PhPerp);
200 61 result_bank.putDouble(i_MX2, dih_row, MX2);
201 61 result_bank.putDouble(i_xF, dih_row, xF);
202 61 result_bank.putDouble(i_yB, dih_row, yB);
203 61 result_bank.putDouble(i_phiH, dih_row, phiH);
204 61 result_bank.putDouble(i_phiR, dih_row, phiR);
205 61 result_bank.putDouble(i_theta, dih_row, theta);
206
207 61 dih_row++;
208 }
209
210
4/8
✓ Branch 122 → 123 taken 286 times.
✗ Branch 122 → 149 not taken.
✓ Branch 125 → 126 taken 286 times.
✗ Branch 125 → 143 not taken.
✓ Branch 131 → 132 taken 38 times.
✓ Branch 131 → 134 taken 248 times.
✗ Branch 149 → 150 not taken.
✗ Branch 149 → 152 not taken.
 858 ShowBank(result_bank, Logger::Header("CREATED BANK"));
211
2/2
✓ Branch 131 → 132 taken 38 times.
✓ Branch 131 → 134 taken 248 times.
 286 return result_bank.getRows() > 0;
212 }
213
214 ///////////////////////////////////////////////////////////////////////////////
215
216
1/2
✓ Branch 2 → 3 taken 286 times.
✗ Branch 2 → 50 not taken.
 286 std::vector<std::pair<int, int>> DihadronKinematics::PairHadrons(hipo::bank const& particle_bank) const
217 {
218 std::vector<std::pair<int, int>> result;
219 // loop over particle bank rows, for hadron A
220
3/4
✓ Branch 2 → 3 taken 286 times.
✗ Branch 2 → 50 not taken.
✓ Branch 38 → 4 taken 1896 times.
✓ Branch 38 → 39 taken 286 times.
 2182 for(auto const& row_a : particle_bank.getRowList()) {
221 // check PDG is in the hadron-A list
222
2/2
✓ Branch 14 → 15 taken 142 times.
✓ Branch 14 → 37 taken 1754 times.
 3792 if(auto pdg_a{particle_bank.getInt("pid", row_a)}; o_hadron_a_pdgs.find(pdg_a) != o_hadron_a_pdgs.end()) {
223 // loop over particle bank rows, for hadron B
224
3/4
✓ Branch 15 → 16 taken 142 times.
✗ Branch 15 → 49 not taken.
✓ Branch 36 → 17 taken 1130 times.
✓ Branch 36 → 37 taken 142 times.
 1272 for(auto const& row_b : particle_bank.getRowList()) {
225 // don't pair a particle with itself
226
2/2
✓ Branch 17 → 18 taken 142 times.
✓ Branch 17 → 19 taken 988 times.
 1130 if(row_a == row_b)
227 142 continue;
228 // check PDG is in the hadron-B list
229
2/2
✓ Branch 29 → 30 taken 75 times.
✓ Branch 29 → 35 taken 913 times.
 1976 if(auto pdg_b{particle_bank.getInt("pid", row_b)}; o_hadron_b_pdgs.find(pdg_b) != o_hadron_b_pdgs.end()) {
230 // if the PDGs of hadrons A and B are the same, don't double count
231
4/4
✓ Branch 30 → 31 taken 28 times.
✓ Branch 30 → 33 taken 47 times.
✓ Branch 31 → 32 taken 14 times.
✓ Branch 31 → 33 taken 14 times.
 75 if(pdg_a == pdg_b && row_b < row_a)
232 continue;
233 // we have a unique dihadron, add it to the list
234 61 result.push_back({row_a, row_b});
235 }
236 }
237 }
238 }
239 // trace logging
240
2/4
✓ Branch 39 → 40 taken 286 times.
✗ Branch 39 → 50 not taken.
✗ Branch 40 → 41 not taken.
✓ Branch 40 → 48 taken 286 times.
 286 if(m_log->GetLevel() <= Logger::Level::trace) {
241 ✗ if(result.empty())
242 ✗ m_log->Trace("=> no dihadrons in this event");
243 else
244 ✗ m_log->Trace("=> number of dihadrons found: {}", result.size());
245 }
246 286 return result;
247 }
248
249 ///////////////////////////////////////////////////////////////////////////////
250
251 1 void DihadronKinematics::Stop()
252 {
253 1 }
254
255 }
256