GCC Code Coverage Report

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