GCC Code Coverage Report

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