GCC Code Coverage Report

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