GCC Code Coverage Report


Directory: ./
File: src/iguana/algorithms/physics/InclusiveKinematics/Algorithm.cc
Date: 2025-03-24 18:50:00
Exec Total Coverage
Lines: 129 142 90.8%
Functions: 7 7 100.0%
Branches: 118 244 48.4%

Line Branch Exec Source
1 #include "Algorithm.h"
2
3 // ROOT
4 #include <Math/Vector4D.h>
5
6 namespace iguana::physics {
7
8 REGISTER_IGUANA_ALGORITHM(InclusiveKinematics, "physics::InclusiveKinematics");
9
10 8 void InclusiveKinematics::Start(hipo::banklist& banks)
11 {
12
1/2
✓ Branch 0 (3→4) taken 8 times.
✗ Branch 1 (3→114) not taken.
8 b_particle = GetBankIndex(banks, "REC::Particle");
13
1/2
✓ Branch 0 (6→7) taken 8 times.
✗ Branch 1 (6→116) not taken.
8 b_config = GetBankIndex(banks, "RUN::config");
14
15 // create the output bank
16
1/2
✓ Branch 0 (9→10) taken 8 times.
✗ Branch 1 (9→118) not taken.
8 auto result_schema = CreateBank(banks, b_result, GetClassName());
17
1/2
✓ Branch 0 (11→12) taken 8 times.
✗ Branch 1 (11→184) not taken.
8 i_pindex = result_schema.getEntryOrder("pindex");
18
1/2
✓ Branch 0 (12→13) taken 8 times.
✗ Branch 1 (12→184) not taken.
8 i_Q2 = result_schema.getEntryOrder("Q2");
19
1/2
✓ Branch 0 (13→14) taken 8 times.
✗ Branch 1 (13→184) not taken.
8 i_x = result_schema.getEntryOrder("x");
20
1/2
✓ Branch 0 (14→15) taken 8 times.
✗ Branch 1 (14→184) not taken.
8 i_y = result_schema.getEntryOrder("y");
21
1/2
✓ Branch 0 (15→16) taken 8 times.
✗ Branch 1 (15→184) not taken.
8 i_W = result_schema.getEntryOrder("W");
22
1/2
✓ Branch 0 (16→17) taken 8 times.
✗ Branch 1 (16→184) not taken.
8 i_nu = result_schema.getEntryOrder("nu");
23
1/2
✓ Branch 0 (17→18) taken 8 times.
✗ Branch 1 (17→184) not taken.
8 i_qx = result_schema.getEntryOrder("qx");
24
1/2
✓ Branch 0 (18→19) taken 8 times.
✗ Branch 1 (18→184) not taken.
8 i_qy = result_schema.getEntryOrder("qy");
25
1/2
✓ Branch 0 (19→20) taken 8 times.
✗ Branch 1 (19→184) not taken.
8 i_qz = result_schema.getEntryOrder("qz");
26
1/2
✓ Branch 0 (20→21) taken 8 times.
✗ Branch 1 (20→184) not taken.
8 i_qE = result_schema.getEntryOrder("qE");
27
1/2
✓ Branch 0 (21→22) taken 8 times.
✗ Branch 1 (21→184) not taken.
8 i_beamPz = result_schema.getEntryOrder("beamPz");
28
2/4
✓ Branch 0 (22→23) taken 8 times.
✗ Branch 1 (22→184) not taken.
✓ Branch 2 (23→24) taken 8 times.
✗ Branch 3 (23→184) not taken.
8 i_targetM = result_schema.getEntryOrder("targetM");
29
30 // instantiate RCDB reader
31
3/6
✓ Branch 0 (23→24) taken 8 times.
✗ Branch 1 (23→184) not taken.
✓ Branch 2 (24→25) taken 8 times.
✗ Branch 3 (24→184) not taken.
✓ Branch 4 (27→28) taken 8 times.
✗ Branch 5 (27→120) not taken.
24 m_rcdb = std::make_unique<RCDBReader>("RCDB|" + GetName(), m_log->GetLevel());
32
33 // parse config file
34
1/2
✓ Branch 0 (32→33) taken 8 times.
✗ Branch 1 (32→184) not taken.
8 ParseYAMLConfig();
35
1/2
✓ Branch 0 (33→34) taken 8 times.
✗ Branch 1 (33→184) not taken.
8 o_runnum = ConcurrentParamFactory::Create<int>();
36
1/2
✓ Branch 0 (37→38) taken 8 times.
✗ Branch 1 (37→184) not taken.
8 o_target_PxPyPzM = ConcurrentParamFactory::Create<std::vector<double>>();
37
2/4
✓ Branch 0 (41→42) taken 8 times.
✗ Branch 1 (41→184) not taken.
✓ Branch 2 (45→46) taken 8 times.
✗ Branch 3 (45→133) not taken.
16 o_beam_PxPyPzM = ConcurrentParamFactory::Create<std::vector<double>>();
38
39 // get reconstruction method configuration
40
5/10
✓ Branch 0 (47→48) taken 8 times.
✗ Branch 1 (47→128) not taken.
✓ Branch 2 (48→49) taken 8 times.
✗ Branch 3 (48→126) not taken.
✓ Branch 4 (49→50) taken 8 times.
✗ Branch 5 (49→124) not taken.
✓ Branch 6 (53→54) taken 16 times.
✓ Branch 7 (53→56) taken 8 times.
✗ Branch 8 (129→130) not taken.
✗ Branch 9 (129→132) not taken.
32 auto method_reconstruction_str = GetOptionScalar<std::string>("reconstruction", {"method", "reconstruction"});
41
1/2
✓ Branch 0 (56→57) taken 8 times.
✗ Branch 1 (56→59) not taken.
8 if(method_reconstruction_str == "scattered_lepton") {
42
1/2
✓ Branch 0 (57→58) taken 8 times.
✗ Branch 1 (57→151) not taken.
8 o_method_reconstruction = method_reconstruction::scattered_lepton;
43 }
44 else {
45 m_log->Error("Unknown reconstruction method {:?}", method_reconstruction_str);
46 throw std::runtime_error("Start failed");
47 }
48
49 // get scattered lepton finder configuration
50
5/10
✓ Branch 0 (66→67) taken 8 times.
✗ Branch 1 (66→146) not taken.
✓ Branch 2 (67→68) taken 8 times.
✗ Branch 3 (67→144) not taken.
✓ Branch 4 (68→69) taken 8 times.
✗ Branch 5 (68→142) not taken.
✓ Branch 6 (72→73) taken 16 times.
✓ Branch 7 (72→75) taken 8 times.
✗ Branch 8 (147→148) not taken.
✗ Branch 9 (147→150) not taken.
32 auto method_lepton_finder_str = GetOptionScalar<std::string>("lepton_finder", {"method", "lepton_finder"});
51
1/2
✓ Branch 0 (75→76) taken 8 times.
✗ Branch 1 (75→78) not taken.
8 if(method_lepton_finder_str == "highest_energy_FD_trigger") {
52 8 o_method_lepton_finder = method_lepton_finder::highest_energy_FD_trigger;
53 }
54 else {
55 m_log->Error("Unknown lepton finder method {:?}", method_lepton_finder_str);
56 throw std::runtime_error("Start failed");
57 }
58
59 // get beam PDG and mass
60
1/2
✓ Branch 0 (76→77) taken 8 times.
✗ Branch 1 (76→169) not taken.
8 o_beam_pdg = 0;
61
5/10
✓ Branch 0 (85→86) taken 8 times.
✗ Branch 1 (85→164) not taken.
✓ Branch 2 (86→87) taken 8 times.
✗ Branch 3 (86→162) not taken.
✓ Branch 4 (87→88) taken 8 times.
✗ Branch 5 (87→160) not taken.
✓ Branch 6 (91→92) taken 16 times.
✓ Branch 7 (91→94) taken 8 times.
✗ Branch 8 (165→166) not taken.
✗ Branch 9 (165→168) not taken.
32 auto beam_particle = GetOptionScalar<std::string>("beam_particle", {"method", "beam_particle"});
62
1/2
✓ Branch 0 (99→95) taken 80 times.
✗ Branch 1 (99→100) not taken.
80 for(auto const& [pdg, name] : particle::name) {
63
2/2
✓ Branch 0 (95→96) taken 8 times.
✓ Branch 1 (95→98) taken 72 times.
80 if(name == beam_particle) {
64
1/2
✓ Branch 0 (96→97) taken 8 times.
✗ Branch 1 (96→178) not taken.
8 o_beam_pdg = pdg;
65 8 o_beam_mass = particle::mass.at(pdg);
66 8 break;
67 }
68 }
69
1/2
✗ Branch 0 (100→101) not taken.
✓ Branch 1 (100→108) taken 8 times.
8 if(o_beam_pdg == 0) {
70 m_log->Error("Unknown beam particle {:?}", beam_particle);
71 throw std::runtime_error("Start failed");
72 }
73
74 8 m_log->Warn("the kinematic calculations in this algorithm need to be cross checked; use this algorithm at your own risk!");
75
0/6
✗ Branch 0 (134→135) not taken.
✗ Branch 1 (134→137) not taken.
✗ Branch 2 (152→153) not taken.
✗ Branch 3 (152→155) not taken.
✗ Branch 4 (170→171) not taken.
✗ Branch 5 (170→173) not taken.
8 }
76
77 ///////////////////////////////////////////////////////////////////////////////
78
79 8000 void InclusiveKinematics::Run(hipo::banklist& banks) const
80 {
81
1/2
✓ Branch 0 (3→4) taken 8000 times.
✗ Branch 1 (3→44) not taken.
8000 auto& particle_bank = GetBank(banks, b_particle, "REC::Particle");
82
1/2
✓ Branch 0 (6→7) taken 8000 times.
✗ Branch 1 (6→46) not taken.
16000 auto& config_bank = GetBank(banks, b_config, "RUN::config");
83
1/2
✓ Branch 0 (9→10) taken 8000 times.
✗ Branch 1 (9→48) not taken.
16000 auto& result_bank = GetBank(banks, b_result, GetClassName());
84
1/2
✓ Branch 0 (12→13) taken 8000 times.
✗ Branch 1 (12→50) not taken.
8000 ShowBank(particle_bank, Logger::Header("INPUT PARTICLES"));
85
86 8000 auto key = PrepareEvent(config_bank.getInt("run",0));
87
88 8000 auto lepton_pindex = FindScatteredLepton(particle_bank, key);
89
2/2
✓ Branch 0 (17→18) taken 7064 times.
✓ Branch 1 (17→22) taken 936 times.
8000 if(lepton_pindex < 0) {
90
1/2
✓ Branch 0 (19→20) taken 7064 times.
✗ Branch 1 (19→52) not taken.
7064 ShowBank(result_bank, Logger::Header("CREATED BANK IS EMPTY"));
91 7064 return;
92 }
93
94 936 auto result_vars = ComputeFromLepton(
95 936 particle_bank.getFloat("px", lepton_pindex),
96 936 particle_bank.getFloat("py", lepton_pindex),
97 936 particle_bank.getFloat("pz", lepton_pindex),
98 key);
99 result_vars.pindex = lepton_pindex; // FIXME: should be done in `ComputeFromLepton`, but need a proper action function first...
100
101 936 result_bank.setRows(1);
102 936 result_bank.putShort(i_pindex, 0, static_cast<int16_t>(result_vars.pindex));
103 936 result_bank.putDouble(i_Q2, 0, result_vars.Q2);
104 936 result_bank.putDouble(i_x, 0, result_vars.x);
105 936 result_bank.putDouble(i_y, 0, result_vars.y);
106 936 result_bank.putDouble(i_W, 0, result_vars.W);
107 936 result_bank.putDouble(i_nu, 0, result_vars.nu);
108 936 result_bank.putDouble(i_qx, 0, result_vars.qx);
109 936 result_bank.putDouble(i_qy, 0, result_vars.qy);
110 936 result_bank.putDouble(i_qz, 0, result_vars.qz);
111 936 result_bank.putDouble(i_qE, 0, result_vars.qE);
112 936 result_bank.putDouble(i_beamPz, 0, result_vars.beamPz);
113 936 result_bank.putDouble(i_targetM, 0, result_vars.targetM);
114
115
1/2
✓ Branch 0 (40→41) taken 936 times.
✗ Branch 1 (40→54) not taken.
1872 ShowBank(result_bank, Logger::Header("CREATED BANK"));
116 }
117
118 ///////////////////////////////////////////////////////////////////////////////
119
120 8000 int InclusiveKinematics::FindScatteredLepton(hipo::bank const& particle_bank, concurrent_key_t const key) const
121 {
122 int lepton_row = -1;
123 double lepton_energy = 0;
124
125
1/2
✓ Branch 0 (2→3) taken 8000 times.
✗ Branch 1 (2→20) not taken.
8000 switch(o_method_lepton_finder) {
126 8000 case method_lepton_finder::highest_energy_FD_trigger: {
127 // find highest energy lepton
128
2/2
✓ Branch 0 (14→5) taken 57144 times.
✓ Branch 1 (14→15) taken 8000 times.
65144 for(auto const& row : particle_bank.getRowList()) {
129
2/2
✓ Branch 0 (6→7) taken 51936 times.
✓ Branch 1 (6→8) taken 5208 times.
57144 if(particle_bank.getInt("pid", row) == o_beam_pdg) {
130 5208 double px = particle_bank.getFloat("px", row);
131 5208 double py = particle_bank.getFloat("py", row);
132
2/2
✓ Branch 0 (11→7) taken 512 times.
✓ Branch 1 (11→12) taken 4696 times.
5208 double pz = particle_bank.getFloat("pz", row);
133 5208 double en = std::sqrt(std::pow(px, 2) + std::pow(py, 2) + std::pow(pz, 2) + std::pow(o_beam_mass, 2));
134
2/2
✓ Branch 0 (11→7) taken 512 times.
✓ Branch 1 (11→12) taken 4696 times.
5208 if(en > lepton_energy) {
135 4696 lepton_row = row;
136 lepton_energy = en;
137 }
138 }
139 }
140 // require it to be in the FD trigger
141
2/2
✓ Branch 0 (15→16) taken 4664 times.
✓ Branch 1 (15→21) taken 3336 times.
8000 if(lepton_row >= 0) {
142 4664 auto status = particle_bank.getShort("status", lepton_row);
143
2/2
✓ Branch 0 (17→18) taken 936 times.
✓ Branch 1 (17→20) taken 3728 times.
4664 if(status <= -3000 || status > -2000)
144 lepton_row = -1;
145 }
146 break;
147 }
148 }
149 if(lepton_row >= 0)
150 936 m_log->Debug("Found scattered lepton: row={}, energy={}", lepton_row, lepton_energy);
151 else
152 7064 m_log->Debug("Scattered lepton not found");
153 8000 return lepton_row;
154 }
155
156 ///////////////////////////////////////////////////////////////////////////////
157
158 8000 concurrent_key_t InclusiveKinematics::PrepareEvent(int const runnum, double const beam_energy) const
159 {
160
2/2
✓ Branch 0 (3→4) taken 4000 times.
✓ Branch 1 (3→8) taken 4000 times.
8000 m_log->Trace("calling PrepareEvent({})", runnum);
161
2/2
✓ Branch 0 (3→4) taken 4000 times.
✓ Branch 1 (3→8) taken 4000 times.
8000 if(o_runnum->NeedsHashing()) {
162 std::hash<int> hash_ftn;
163 auto hash_key = hash_ftn(runnum);
164
2/2
✓ Branch 0 (5→6) taken 4 times.
✓ Branch 1 (5→7) taken 3996 times.
4000 if(!o_runnum->HasKey(hash_key))
165 4 Reload(runnum, beam_energy, hash_key);
166 return hash_key;
167 } else {
168
3/4
✓ Branch 0 (8→9) taken 3996 times.
✓ Branch 1 (8→11) taken 4 times.
✗ Branch 2 (10→11) not taken.
✓ Branch 3 (10→12) taken 3996 times.
4000 if(o_runnum->IsEmpty() || o_runnum->Load(0) != runnum)
169 4 Reload(runnum, beam_energy, 0);
170 4000 return 0;
171 }
172 }
173
174 ///////////////////////////////////////////////////////////////////////////////
175
176 8 void InclusiveKinematics::Reload(int const runnum, double const user_beam_energy, concurrent_key_t key) const
177 {
178 8 std::lock_guard<std::mutex> const lock(m_mutex);
179
2/4
✓ Branch 0 (3→4) taken 8 times.
✗ Branch 1 (3→148) not taken.
✓ Branch 2 (4→5) taken 8 times.
✗ Branch 3 (4→148) not taken.
8 m_log->Trace("-> calling Reload({}, {}, {})", runnum, user_beam_energy, key);
180
1/2
✓ Branch 0 (4→5) taken 8 times.
✗ Branch 1 (4→148) not taken.
8 o_runnum->Save(runnum, key);
181
182 // parse config params
183
2/4
✓ Branch 0 (5→6) taken 8 times.
✗ Branch 1 (5→7) not taken.
✓ Branch 2 (6→7) taken 8 times.
✗ Branch 3 (6→148) not taken.
8 auto beam_energy = user_beam_energy < 0 ? m_rcdb->GetBeamEnergy(runnum) : user_beam_energy;
184
10/20
✓ Branch 0 (7→8) taken 8 times.
✗ Branch 1 (7→104) not taken.
✓ Branch 2 (8→9) taken 8 times.
✗ Branch 3 (8→104) not taken.
✓ Branch 4 (9→10) taken 8 times.
✗ Branch 5 (9→104) not taken.
✓ Branch 6 (10→11) taken 8 times.
✗ Branch 7 (10→102) not taken.
✓ Branch 8 (14→15) taken 8 times.
✗ Branch 9 (14→94) not taken.
✓ Branch 10 (15→16) taken 8 times.
✗ Branch 11 (15→92) not taken.
✓ Branch 12 (16→17) taken 8 times.
✗ Branch 13 (16→90) not taken.
✓ Branch 14 (20→21) taken 24 times.
✓ Branch 15 (20→23) taken 8 times.
✓ Branch 16 (26→27) taken 8 times.
✗ Branch 17 (26→123) not taken.
✗ Branch 18 (95→96) not taken.
✗ Branch 19 (95→98) not taken.
56 auto beam_direction = GetOptionVector<double>("beam_direction", {"initial_state", GetConfig()->InRange("runs", runnum), "beam_direction"});
185
8/18
✓ Branch 0 (27→28) taken 8 times.
✗ Branch 1 (27→123) not taken.
✓ Branch 2 (28→29) taken 8 times.
✗ Branch 3 (28→123) not taken.
✓ Branch 4 (29→30) taken 8 times.
✗ Branch 5 (29→121) not taken.
✓ Branch 6 (33→34) taken 8 times.
✗ Branch 7 (33→113) not taken.
✓ Branch 8 (34→35) taken 8 times.
✗ Branch 9 (34→111) not taken.
✓ Branch 10 (35→36) taken 8 times.
✗ Branch 11 (35→109) not taken.
✓ Branch 12 (39→40) taken 24 times.
✓ Branch 13 (39→42) taken 8 times.
✗ Branch 14 (114→115) not taken.
✗ Branch 15 (114→117) not taken.
✗ Branch 16 (144→145) not taken.
✗ Branch 17 (144→147) not taken.
56 auto target_particle = GetOptionScalar<std::string>("target_particle", {"initial_state", GetConfig()->InRange("runs", runnum), "target_particle"});
186
187 // get the target mass and momentum
188 double target_mass = -1;
189
1/2
✓ Branch 0 (50→46) taken 56 times.
✗ Branch 1 (50→51) not taken.
56 for(auto const& [pdg, name] : particle::name) {
190
2/2
✓ Branch 0 (46→47) taken 8 times.
✓ Branch 1 (46→49) taken 48 times.
56 if(name == target_particle) {
191 8 target_mass = particle::mass.at(pdg);
192 8 break;
193 }
194 }
195
1/2
✗ Branch 0 (51→52) not taken.
✓ Branch 1 (51→59) taken 8 times.
8 if(target_mass < 0) {
196 m_log->Error("Unknown target particle {:?}", target_particle);
197 throw std::runtime_error("Reload failed");
198 }
199 double target_px = 0.0;
200 double target_py = 0.0;
201 double target_pz = 0.0;
202
203 // get the beam momentum
204 double beam_px, beam_py, beam_pz;
205
1/2
✗ Branch 0 (59→60) not taken.
✓ Branch 1 (59→64) taken 8 times.
8 if(beam_direction.size() != 3) {
206 m_log->Error("Beam direction is not a 3-vector; assuming it is (0, 0, 1) instead");
207 beam_direction = {0.0, 0.0, 1.0};
208 }
209
1/2
✗ Branch 0 (66→67) not taken.
✓ Branch 1 (66→68) taken 8 times.
8 auto dir_mag = std::hypot(beam_direction[0], beam_direction[1], beam_direction[2]);
210 8 auto beam_p = std::sqrt(std::pow(beam_energy, 2) - std::pow(o_beam_mass, 2));
211
1/2
✓ Branch 0 (69→70) taken 8 times.
✗ Branch 1 (69→71) not taken.
8 if(dir_mag > 0) {
212 8 beam_px = beam_direction[0] * beam_p / dir_mag;
213 8 beam_py = beam_direction[1] * beam_p / dir_mag;
214
1/2
✓ Branch 0 (70→75) taken 8 times.
✗ Branch 1 (70→142) not taken.
8 beam_pz = beam_direction[2] * beam_p / dir_mag;
215 }
216 else {
217 m_log->Error("Beam direction magnitude is not > 0");
218 throw ::std::runtime_error("Reload failed");
219 }
220
221 // save the configuration
222
3/6
✓ Branch 0 (70→75) taken 8 times.
✗ Branch 1 (70→142) not taken.
✓ Branch 2 (75→76) taken 8 times.
✗ Branch 3 (75→134) not taken.
✓ Branch 4 (79→80) taken 8 times.
✗ Branch 5 (79→142) not taken.
16 o_beam_PxPyPzM->Save({beam_px, beam_py, beam_pz, o_beam_mass}, key);
223
2/4
✓ Branch 0 (79→80) taken 8 times.
✗ Branch 1 (79→142) not taken.
✓ Branch 2 (80→81) taken 8 times.
✗ Branch 3 (80→138) not taken.
16 o_target_PxPyPzM->Save({target_px, target_py, target_pz, target_mass}, key);
224
0/4
✗ Branch 0 (105→106) not taken.
✗ Branch 1 (105→108) not taken.
✗ Branch 2 (124→125) not taken.
✗ Branch 3 (124→127) not taken.
8 }
225
226 ///////////////////////////////////////////////////////////////////////////////
227
228 936 InclusiveKinematicsVars InclusiveKinematics::ComputeFromLepton(
229 vector_element_t const lepton_px,
230 vector_element_t const lepton_py,
231 vector_element_t const lepton_pz,
232 concurrent_key_t const key
233 ) const
234 {
235 InclusiveKinematicsVars result;
236
237 936 m_log->Trace("Reconstruct inclusive kinematics from lepton with p=({}, {}, {})", lepton_px, lepton_py, lepton_pz);
238
239 enum {px, py, pz, m};
240 936 auto beam = o_beam_PxPyPzM->Load(key);
241
1/2
✓ Branch 0 (4→5) taken 936 times.
✗ Branch 1 (4→48) not taken.
936 auto target = o_target_PxPyPzM->Load(key);
242
243 ROOT::Math::PxPyPzMVector vec_beam(beam[px], beam[py], beam[pz], beam[m]);
244 ROOT::Math::PxPyPzMVector vec_target(target[px], target[py], target[pz], target[m]);
245 ROOT::Math::PxPyPzMVector vec_lepton(lepton_px, lepton_py, lepton_pz, beam[m]);
246
247 auto vec_q = vec_beam - vec_lepton;
248 936 result.qx = vec_q.Px();
249 936 result.qy = vec_q.Py();
250 936 result.qz = vec_q.Pz();
251
1/2
✗ Branch 0 (24→25) not taken.
✓ Branch 1 (24→26) taken 936 times.
936 result.qE = vec_q.E();
252 936 result.Q2 = -1 * vec_q.M2();
253 936 result.x = result.Q2 / (2 * vec_q.Dot(vec_target));
254
1/2
✓ Branch 0 (30→31) taken 936 times.
✗ Branch 1 (30→44) not taken.
936 result.y = vec_target.Dot(vec_q) / vec_target.Dot(vec_beam);
255 result.W = (vec_beam + vec_target - vec_lepton).M();
256
1/2
✗ Branch 0 (37→38) not taken.
✓ Branch 1 (37→39) taken 936 times.
1872 result.nu = vec_target.Dot(vec_q) / target[m];
257 936 result.beamPz = beam[pz];
258 936 result.targetM = target[m];
259
260 936 return result;
261 }
262
263 ///////////////////////////////////////////////////////////////////////////////
264
265 4 void InclusiveKinematics::Stop()
266 {
267 4 }
268
269 }
270