GCC Code Coverage Report

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