GCC Code Coverage Report

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