GCC Code Coverage Report

Directory: ./
File: src/iguana/algorithms/physics/InclusiveKinematics/Algorithm.cc
Date: 2025-12-23 16:10:08
Coverage Exec Excl Total
Lines: 83.1% 147 0 177
Functions: 100.0% 8 0 8
Branches: 41.1% 132 0 321
Line Branch Exec Source
1 #include "Algorithm.h"
2
3 // ROOT
4 #include <Math/Vector3D.h>
5 #include <Math/Vector4D.h>
6 #include <TMath.h>
7
8 namespace iguana::physics {
9
10 REGISTER_IGUANA_ALGORITHM(InclusiveKinematics, "physics::InclusiveKinematics");
11
12 8 void InclusiveKinematics::Start(hipo::banklist& banks)
13 {
14 // parse config file
15 8 ParseYAMLConfig();
16
2/4
✓ Branch 4 → 5 taken 8 times.
✗ Branch 4 → 198 not taken.
✓ Branch 5 → 6 taken 8 times.
✗ Branch 5 → 192 not taken.
 16 o_particle_bank = GetOptionScalar<std::string>("particle_bank");
17 8 o_runnum = ConcurrentParamFactory::Create<int>();
18 8 o_target_PxPyPzM = ConcurrentParamFactory::Create<std::vector<double>>();
19 8 o_beam_PxPyPzM = ConcurrentParamFactory::Create<std::vector<double>>();
20
4/8
✓ Branch 31 → 32 taken 8 times.
✗ Branch 31 → 206 not taken.
✓ Branch 32 → 33 taken 8 times.
✗ Branch 32 → 200 not taken.
✗ Branch 33 → 34 not taken.
✓ Branch 33 → 36 taken 8 times.
✓ Branch 39 → 40 taken 8 times.
✗ Branch 39 → 221 not taken.
 16 o_theta_between_FD_and_FT = GetOptionScalar<double>("theta_between_FD_and_FT");
21
22 // get reconstruction method configuration
23
5/10
✓ Branch 41 → 42 taken 8 times.
✗ Branch 41 → 216 not taken.
✓ Branch 42 → 43 taken 8 times.
✗ Branch 42 → 214 not taken.
✓ Branch 43 → 44 taken 8 times.
✗ Branch 43 → 208 not taken.
✓ Branch 51 → 52 taken 16 times.
✓ Branch 51 → 54 taken 8 times.
✗ Branch 217 → 218 not taken.
✗ Branch 217 → 220 not taken.
 32 auto method_reconstruction_str = GetOptionScalar<std::string>("reconstruction", {"method", "reconstruction"});
24
1/2
✓ Branch 54 → 55 taken 8 times.
✗ Branch 54 → 57 not taken.
 8 if(method_reconstruction_str == "scattered_lepton") {
25
1/2
✓ Branch 55 → 56 taken 8 times.
✗ Branch 55 → 247 not taken.
 8 o_method_reconstruction = method_reconstruction::scattered_lepton;
26 }
27 else {
28 ✗ m_log->Error("Unknown reconstruction method {:?}", method_reconstruction_str);
29 ✗ throw std::runtime_error("Start failed");
30 }
31
32 // get scattered lepton finder configuration
33
5/12
✓ Branch 70 → 71 taken 8 times.
✗ Branch 70 → 242 not taken.
✓ Branch 71 → 72 taken 8 times.
✗ Branch 71 → 240 not taken.
✓ Branch 72 → 73 taken 8 times.
✗ Branch 72 → 234 not taken.
✓ Branch 80 → 81 taken 16 times.
✓ Branch 80 → 83 taken 8 times.
✗ Branch 243 → 244 not taken.
✗ Branch 243 → 246 not taken.
✗ Branch 320 → 321 not taken.
✗ Branch 320 → 323 not taken.
 32 auto method_lepton_finder_str = GetOptionScalar<std::string>("lepton_finder", {"method", "lepton_finder"});
34
1/2
✓ Branch 83 → 84 taken 8 times.
✗ Branch 83 → 85 not taken.
 8 if(method_lepton_finder_str == "highest_energy_FD_trigger")
35 8 o_method_lepton_finder = method_lepton_finder::highest_energy_FD_trigger;
36 ✗ else if(method_lepton_finder_str == "lund_beam_daughter")
37 ✗ o_method_lepton_finder = method_lepton_finder::lund_beam_daughter;
38 else {
39 ✗ m_log->Error("Unknown lepton finder method {:?}", method_lepton_finder_str);
40 ✗ throw std::runtime_error("Start failed");
41 }
42
43 // get beam PDG and mass
44
1/2
✓ Branch 100 → 101 taken 8 times.
✗ Branch 100 → 273 not taken.
 8 o_beam_pdg = 0;
45
5/12
✓ Branch 102 → 103 taken 8 times.
✗ Branch 102 → 268 not taken.
✓ Branch 103 → 104 taken 8 times.
✗ Branch 103 → 266 not taken.
✓ Branch 104 → 105 taken 8 times.
✗ Branch 104 → 260 not taken.
✓ Branch 112 → 113 taken 16 times.
✓ Branch 112 → 115 taken 8 times.
✗ Branch 269 → 270 not taken.
✗ Branch 269 → 272 not taken.
✗ Branch 314 → 315 not taken.
✗ Branch 314 → 317 not taken.
 32 auto beam_particle = GetOptionScalar<std::string>("beam_particle", {"method", "beam_particle"});
46
1/2
✓ Branch 120 → 116 taken 80 times.
✗ Branch 120 → 121 not taken.
 80 for(auto const& [pdg, name] : particle::name) {
47
2/2
✓ Branch 116 → 117 taken 8 times.
✓ Branch 116 → 119 taken 72 times.
 80 if(name == beam_particle) {
48
1/2
✓ Branch 117 → 118 taken 8 times.
✗ Branch 117 → 308 not taken.
 8 o_beam_pdg = pdg;
49 8 o_beam_mass = particle::mass.at(pdg);
50 8 break;
51 }
52 }
53
1/2
✗ Branch 121 → 122 not taken.
✓ Branch 121 → 135 taken 8 times.
 8 if(o_beam_pdg == 0) {
54 ✗ m_log->Error("Unknown beam particle {:?}", beam_particle);
55 ✗ throw std::runtime_error("Start failed");
56 }
57
58 // get bank indices
59
1/2
✓ Branch 135 → 136 taken 8 times.
✗ Branch 135 → 308 not taken.
 8 b_particle = GetBankIndex(banks, o_particle_bank);
60
4/8
✓ Branch 136 → 137 taken 8 times.
✗ Branch 136 → 308 not taken.
✓ Branch 137 → 138 taken 8 times.
✗ Branch 137 → 286 not taken.
✓ Branch 138 → 139 taken 8 times.
✗ Branch 138 → 141 not taken.
✓ Branch 143 → 144 taken 8 times.
✗ Branch 143 → 308 not taken.
 16 b_config = GetBankIndex(banks, "RUN::config");
61
62 // create the output bank
63
1/4
✓ Branch 144 → 145 taken 8 times.
✗ Branch 144 → 292 not taken.
✗ Branch 308 → 309 not taken.
✗ Branch 308 → 311 not taken.
 8 auto result_schema = CreateBank(banks, b_result, GetClassName());
64
1/2
✓ Branch 150 → 151 taken 8 times.
✗ Branch 150 → 306 not taken.
 8 i_pindex = result_schema.getEntryOrder("pindex");
65
1/2
✓ Branch 151 → 152 taken 8 times.
✗ Branch 151 → 306 not taken.
 8 i_Q2 = result_schema.getEntryOrder("Q2");
66
1/2
✓ Branch 152 → 153 taken 8 times.
✗ Branch 152 → 306 not taken.
 8 i_x = result_schema.getEntryOrder("x");
67
1/2
✓ Branch 153 → 154 taken 8 times.
✗ Branch 153 → 306 not taken.
 8 i_y = result_schema.getEntryOrder("y");
68
1/2
✓ Branch 154 → 155 taken 8 times.
✗ Branch 154 → 306 not taken.
 8 i_W = result_schema.getEntryOrder("W");
69
1/2
✓ Branch 155 → 156 taken 8 times.
✗ Branch 155 → 306 not taken.
 8 i_nu = result_schema.getEntryOrder("nu");
70
1/2
✓ Branch 156 → 157 taken 8 times.
✗ Branch 156 → 306 not taken.
 8 i_qx = result_schema.getEntryOrder("qx");
71
1/2
✓ Branch 157 → 158 taken 8 times.
✗ Branch 157 → 306 not taken.
 8 i_qy = result_schema.getEntryOrder("qy");
72
1/2
✓ Branch 158 → 159 taken 8 times.
✗ Branch 158 → 306 not taken.
 8 i_qz = result_schema.getEntryOrder("qz");
73
1/2
✓ Branch 159 → 160 taken 8 times.
✗ Branch 159 → 306 not taken.
 8 i_qE = result_schema.getEntryOrder("qE");
74
1/2
✓ Branch 160 → 161 taken 8 times.
✗ Branch 160 → 306 not taken.
 8 i_beamPz = result_schema.getEntryOrder("beamPz");
75
1/2
✓ Branch 161 → 162 taken 8 times.
✗ Branch 161 → 306 not taken.
 8 i_targetM = result_schema.getEntryOrder("targetM");
76
77 // instantiate RCDB reader `m_rcdb`
78
1/2
✓ Branch 162 → 163 taken 8 times.
✗ Branch 162 → 306 not taken.
 8 StartRCDBReader();
79
4/8
✓ Branch 163 → 164 taken 8 times.
✗ Branch 163 → 306 not taken.
✓ Branch 165 → 166 taken 8 times.
✗ Branch 165 → 304 not taken.
✓ Branch 166 → 167 taken 8 times.
✗ Branch 166 → 298 not taken.
✗ Branch 167 → 168 not taken.
✓ Branch 167 → 170 taken 8 times.
 24 o_override_beam_energy = GetOptionScalar<double>("override_beam_energy");
80
1/2
✗ Branch 173 → 174 not taken.
✓ Branch 173 → 175 taken 8 times.
 8 if(o_override_beam_energy > 0)
81 ✗ m_rcdb->SetBeamEnergyOverride(o_override_beam_energy);
82
0/6
✗ Branch 222 → 223 not taken.
✗ Branch 222 → 225 not taken.
✗ Branch 248 → 249 not taken.
✗ Branch 248 → 251 not taken.
✗ Branch 274 → 275 not taken.
✗ Branch 274 → 277 not taken.
 16 }
83
84 ///////////////////////////////////////////////////////////////////////////////
85
86 8000 bool InclusiveKinematics::Run(hipo::banklist& banks) const
87 {
88
1/2
✓ Branch 7 → 8 taken 8000 times.
✗ Branch 7 → 19 not taken.
 16000 return Run(
89
1/2
✓ Branch 6 → 7 taken 8000 times.
✗ Branch 6 → 19 not taken.
 8000 GetBank(banks, b_particle, o_particle_bank),
90
3/8
✓ Branch 4 → 5 taken 8000 times.
✗ Branch 4 → 25 not taken.
✓ Branch 5 → 6 taken 8000 times.
✗ Branch 5 → 19 not taken.
✗ Branch 13 → 14 not taken.
✓ Branch 13 → 16 taken 8000 times.
✗ Branch 25 → 26 not taken.
✗ Branch 25 → 28 not taken.
 16000 GetBank(banks, b_config, "RUN::config"),
91
1/2
✓ Branch 3 → 4 taken 8000 times.
✗ Branch 3 → 25 not taken.
 16000 GetBank(banks, b_result, GetClassName()));
92 }
93
94 8000 bool InclusiveKinematics::Run(
95 hipo::bank const& particle_bank,
96 hipo::bank const& config_bank,
97 hipo::bank& result_bank) const
98 {
99 8000 result_bank.reset(); // IMPORTANT: always first `reset` the created bank(s)
100
1/2
✓ Branch 6 → 7 taken 8000 times.
✗ Branch 6 → 54 not taken.
 16000 ShowBank(particle_bank, Logger::Header("INPUT PARTICLES"));
101
102 8000 auto key = PrepareEvent(config_bank.getInt("run", 0));
103
104 8000 auto const lepton_pindex = FindScatteredLepton(particle_bank, key);
105
2/2
✓ Branch 15 → 16 taken 6856 times.
✓ Branch 15 → 26 taken 1144 times.
 8000 if(!lepton_pindex.has_value()) {
106
1/2
✓ Branch 19 → 20 taken 6856 times.
✗ Branch 19 → 60 not taken.
 13712 ShowBank(result_bank, Logger::Header("CREATED BANK IS EMPTY"));
107 6856 return false;
108 }
109
110 1144 auto result_vars = ComputeFromLepton(
111 1144 particle_bank.getFloat("px", lepton_pindex.value()),
112 1144 particle_bank.getFloat("py", lepton_pindex.value()),
113 1144 particle_bank.getFloat("pz", lepton_pindex.value()),
114 key);
115 result_vars.pindex = lepton_pindex.value(); // FIXME: should be done in `ComputeFromLepton`, but need a proper action function first...
116
117 1144 result_bank.setRows(1);
118 1144 result_bank.putShort(i_pindex, 0, static_cast<int16_t>(result_vars.pindex));
119 1144 result_bank.putDouble(i_Q2, 0, result_vars.Q2);
120 1144 result_bank.putDouble(i_x, 0, result_vars.x);
121 1144 result_bank.putDouble(i_y, 0, result_vars.y);
122 1144 result_bank.putDouble(i_W, 0, result_vars.W);
123 1144 result_bank.putDouble(i_nu, 0, result_vars.nu);
124 1144 result_bank.putDouble(i_qx, 0, result_vars.qx);
125 1144 result_bank.putDouble(i_qy, 0, result_vars.qy);
126 1144 result_bank.putDouble(i_qz, 0, result_vars.qz);
127 1144 result_bank.putDouble(i_qE, 0, result_vars.qE);
128 1144 result_bank.putDouble(i_beamPz, 0, result_vars.beamPz);
129 1144 result_bank.putDouble(i_targetM, 0, result_vars.targetM);
130
131
1/2
✓ Branch 46 → 47 taken 1144 times.
✗ Branch 46 → 66 not taken.
 2288 ShowBank(result_bank, Logger::Header("CREATED BANK"));
132 1144 return true;
133 }
134
135 ///////////////////////////////////////////////////////////////////////////////
136
137 8000 std::optional<int> const InclusiveKinematics::FindScatteredLepton(hipo::bank const& particle_bank, concurrent_key_t const key) const
138 {
139 8000 std::optional<int> lepton_row = std::nullopt;
140 double lepton_energy = 0.0;
141
142
1/3
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 8000 times.
✗ Branch 2 → 47 not taken.
 8000 switch(o_method_lepton_finder) {
143 // ----------------------------------------------------------------------------------
144 // highest energy FD trigger lepton
145 // ----------------------------------------------------------------------------------
146 case method_lepton_finder::highest_energy_FD_trigger: {
147 // the `status` variable does not exist if we're looking at `MC::Particle`
148 8000 bool has_status = const_cast<hipo::bank&>(particle_bank).getSchema().exists("status");
149 // loop over ALL rows, not just filtered rows, since we don't want to accidentally pick the wrong lepton
150
2/2
✓ Branch 23 → 5 taken 55944 times.
✓ Branch 23 → 47 taken 8000 times.
 63944 for(int row = 0; row < particle_bank.getRows(); row++) {
151
2/2
✓ Branch 6 → 7 taken 50832 times.
✓ Branch 6 → 8 taken 5112 times.
 55944 if(particle_bank.getInt("pid", row) == o_beam_pdg) { // if beam PDG
152 // check if in FD: use `status` if we have it, otherwise rough theta cut
153 bool in_FD_trigger = false;
154
1/2
✓ Branch 8 → 9 taken 5112 times.
✗ Branch 8 → 11 not taken.
 5112 if(has_status) {
155 5112 auto status = particle_bank.getShort("status", row);
156 5112 in_FD_trigger = status > -3000 && status <= -2000; // trigger && in FD
157 }
158 else {
159 ROOT::Math::XYZVector p(
160 ✗ particle_bank.getFloat("px", row),
161 ✗ particle_bank.getFloat("py", row),
162 ✗ particle_bank.getFloat("pz", row));
163 ✗ in_FD_trigger = p.theta() * TMath::RadToDeg() > o_theta_between_FD_and_FT; // rough theta cut
164 }
165
2/2
✓ Branch 16 → 7 taken 3968 times.
✓ Branch 16 → 17 taken 1144 times.
 5112 if(in_FD_trigger) {
166 1144 m_log->Trace("row {} is in FD trigger", row);
167 1144 double en = std::sqrt(
168 1144 std::pow(particle_bank.getFloat("px", row), 2) +
169 1144 std::pow(particle_bank.getFloat("py", row), 2) +
170 1144 std::pow(particle_bank.getFloat("pz", row), 2) +
171 1144 std::pow(o_beam_mass, 2));
172
1/2
✗ Branch 21 → 7 not taken.
✓ Branch 21 → 22 taken 1144 times.
 1144 if(en > lepton_energy) { // select max-E
173 lepton_row = row;
174 lepton_energy = en;
175 }
176 }
177 }
178 }
179 break;
180 }
181 // ----------------------------------------------------------------------------------
182 // use MC::Lund to find the lepton that has a beam parent
183 // ----------------------------------------------------------------------------------
184 case method_lepton_finder::lund_beam_daughter: {
185 // find the beam lepton, assuming it has parent index == 0
186 // loop over ALL rows, in case the user filtered out beam particles
187 std::optional<int> beam_index = std::nullopt;
188 ✗ for(int row = 0; row < particle_bank.getRows(); row++) {
189 ✗ if(particle_bank.getInt("pid", row) == o_beam_pdg && particle_bank.getByte("parent", row) == 0) {
190 ✗ beam_index = particle_bank.getByte("index", row);
191 ✗ break;
192 // FIXME: should we check if there are more than 1?
193 }
194 }
195 // find the lepton with parent == beam lepton
196 // loop over ALL rows, not just filtered rows, since we don't want to accidentally pick the wrong lepton
197 ✗ if(beam_index.has_value()) {
198 ✗ for(int row = 0; row < particle_bank.getRows(); row++) {
199 ✗ if(particle_bank.getInt("pid", row) == o_beam_pdg && particle_bank.getByte("parent", row) == beam_index.value()) {
200 lepton_row = row;
201 break;
202 // FIXME: should we check if there are more than 1?
203 }
204 }
205 }
206 else
207 ✗ m_log->Debug("Failed to find beam lepton");
208 // complain if lepton not found
209 ✗ if(!lepton_row.has_value())
210 ✗ m_log->Debug("Failed to find scattered lepton");
211 break;
212 }
213 }
214
215 // make sure `lepton_row` was not filtered out
216
2/2
✓ Branch 47 → 48 taken 1144 times.
✓ Branch 47 → 60 taken 6856 times.
 8000 if(lepton_row.has_value()) {
217 1144 auto rowlist = particle_bank.getRowList();
218
1/2
✗ Branch 54 → 55 not taken.
✓ Branch 54 → 56 taken 1144 times.
 1144 if(std::find(rowlist.begin(), rowlist.end(), lepton_row.value()) == rowlist.end())
219 lepton_row = std::nullopt;
220 }
221
222 // return
223
2/2
✓ Branch 60 → 61 taken 1144 times.
✓ Branch 60 → 63 taken 6856 times.
 8000 if(lepton_row.has_value())
224 1144 m_log->Debug("Found scattered lepton: row={}", lepton_row.value());
225 else
226 6856 m_log->Debug("Scattered lepton not found");
227 8000 return lepton_row;
228 }
229
230 ///////////////////////////////////////////////////////////////////////////////
231
232 8000 concurrent_key_t InclusiveKinematics::PrepareEvent(int const runnum, double const beam_energy) const
233 {
234
2/2
✓ Branch 3 → 4 taken 4000 times.
✓ Branch 3 → 8 taken 4000 times.
 8000 m_log->Trace("calling PrepareEvent({})", runnum);
235
2/2
✓ Branch 3 → 4 taken 4000 times.
✓ Branch 3 → 8 taken 4000 times.
 8000 if(o_runnum->NeedsHashing()) {
236 std::hash<int> hash_ftn;
237 auto hash_key = hash_ftn(runnum);
238
2/2
✓ Branch 5 → 6 taken 4 times.
✓ Branch 5 → 7 taken 3996 times.
 4000 if(!o_runnum->HasKey(hash_key))
239 4 Reload(runnum, beam_energy, hash_key);
240 return hash_key;
241 }
242 else {
243
3/4
✓ Branch 8 → 9 taken 3996 times.
✓ Branch 8 → 11 taken 4 times.
✗ Branch 10 → 11 not taken.
✓ Branch 10 → 12 taken 3996 times.
 4000 if(o_runnum->IsEmpty() || o_runnum->Load(0) != runnum)
244 4 Reload(runnum, beam_energy, 0);
245 4000 return 0;
246 }
247 }
248
249 ///////////////////////////////////////////////////////////////////////////////
250
251 8 void InclusiveKinematics::Reload(int const runnum, double const user_beam_energy, concurrent_key_t key) const
252 {
253 8 std::lock_guard<std::mutex> const lock(m_mutex);
254
2/4
✓ Branch 3 → 4 taken 8 times.
✗ Branch 3 → 202 not taken.
✓ Branch 4 → 5 taken 8 times.
✗ Branch 4 → 202 not taken.
 8 m_log->Trace("-> calling Reload({}, {}, {})", runnum, user_beam_energy, key);
255
1/2
✓ Branch 4 → 5 taken 8 times.
✗ Branch 4 → 202 not taken.
 8 o_runnum->Save(runnum, key);
256
257 // parse config params
258
2/4
✓ Branch 5 → 6 taken 8 times.
✗ Branch 5 → 7 not taken.
✓ Branch 6 → 7 taken 8 times.
✗ Branch 6 → 202 not taken.
 8 auto beam_energy = user_beam_energy < 0 ? m_rcdb->GetBeamEnergy(runnum) : user_beam_energy;
259
10/20
✓ Branch 7 → 8 taken 8 times.
✗ Branch 7 → 142 not taken.
✓ Branch 8 → 9 taken 8 times.
✗ Branch 8 → 142 not taken.
✓ Branch 9 → 10 taken 8 times.
✗ Branch 9 → 142 not taken.
✓ Branch 10 → 11 taken 8 times.
✗ Branch 10 → 136 not taken.
✓ Branch 14 → 15 taken 8 times.
✗ Branch 14 → 128 not taken.
✓ Branch 15 → 16 taken 8 times.
✗ Branch 15 → 126 not taken.
✓ Branch 16 → 17 taken 8 times.
✗ Branch 16 → 120 not taken.
✓ Branch 24 → 25 taken 24 times.
✓ Branch 24 → 27 taken 8 times.
✓ Branch 34 → 35 taken 8 times.
✗ Branch 34 → 169 not taken.
✗ Branch 129 → 130 not taken.
✗ Branch 129 → 132 not taken.
 56 auto beam_direction = GetOptionVector<double>("beam_direction", {"initial_state", GetConfig()->InRange("runs", runnum), "beam_direction"});
260
8/18
✓ Branch 35 → 36 taken 8 times.
✗ Branch 35 → 169 not taken.
✓ Branch 36 → 37 taken 8 times.
✗ Branch 36 → 169 not taken.
✓ Branch 37 → 38 taken 8 times.
✗ Branch 37 → 163 not taken.
✓ Branch 41 → 42 taken 8 times.
✗ Branch 41 → 155 not taken.
✓ Branch 42 → 43 taken 8 times.
✗ Branch 42 → 153 not taken.
✓ Branch 43 → 44 taken 8 times.
✗ Branch 43 → 147 not taken.
✓ Branch 51 → 52 taken 24 times.
✓ Branch 51 → 54 taken 8 times.
✗ Branch 156 → 157 not taken.
✗ Branch 156 → 159 not taken.
✗ Branch 198 → 199 not taken.
✗ Branch 198 → 201 not taken.
 56 auto target_particle = GetOptionScalar<std::string>("target_particle", {"initial_state", GetConfig()->InRange("runs", runnum), "target_particle"});
261
262 // get the target mass and momentum
263 double target_mass = -1;
264
1/2
✓ Branch 66 → 62 taken 56 times.
✗ Branch 66 → 67 not taken.
 56 for(auto const& [pdg, name] : particle::name) {
265
2/2
✓ Branch 62 → 63 taken 8 times.
✓ Branch 62 → 65 taken 48 times.
 56 if(name == target_particle) {
266 8 target_mass = particle::mass.at(pdg);
267 8 break;
268 }
269 }
270
1/2
✗ Branch 67 → 68 not taken.
✓ Branch 67 → 81 taken 8 times.
 8 if(target_mass < 0) {
271 ✗ m_log->Error("Unknown target particle {:?}", target_particle);
272 ✗ throw std::runtime_error("Reload failed");
273 }
274 double target_px = 0.0;
275 double target_py = 0.0;
276 double target_pz = 0.0;
277
278 // get the beam momentum
279 double beam_px, beam_py, beam_pz;
280
1/2
✗ Branch 83 → 84 not taken.
✓ Branch 83 → 88 taken 8 times.
 8 if(beam_direction.size() != 3) {
281 ✗ m_log->Error("Beam direction is not a 3-vector; assuming it is (0, 0, 1) instead");
282 ✗ beam_direction = {0.0, 0.0, 1.0};
283 }
284 8 auto dir_mag = std::hypot(beam_direction[0], beam_direction[1], beam_direction[2]);
285 8 auto beam_p = std::sqrt(std::pow(beam_energy, 2) - std::pow(o_beam_mass, 2));
286
1/2
✓ Branch 93 → 94 taken 8 times.
✗ Branch 93 → 96 not taken.
 8 if(dir_mag > 0) {
287 8 beam_px = beam_direction[0] * beam_p / dir_mag;
288 8 beam_py = beam_direction[1] * beam_p / dir_mag;
289
1/2
✓ Branch 94 → 95 taken 8 times.
✗ Branch 94 → 192 not taken.
 8 beam_pz = beam_direction[2] * beam_p / dir_mag;
290 }
291 else {
292 ✗ m_log->Error("Beam direction magnitude is not > 0");
293 ✗ throw ::std::runtime_error("Reload failed");
294 }
295
296 // save the configuration
297
1/2
✓ Branch 95 → 100 taken 8 times.
✗ Branch 95 → 192 not taken.
 8 m_log->Trace("-> Reloaded beam: ({}, {}, {}, {})", beam_px, beam_py, beam_pz, o_beam_mass);
298
1/2
✓ Branch 100 → 101 taken 8 times.
✗ Branch 100 → 192 not taken.
 8 m_log->Trace("-> Reloaded target: ({}, {}, {}, {})", target_px, target_py, target_pz, target_mass);
299
3/6
✓ Branch 100 → 101 taken 8 times.
✗ Branch 100 → 192 not taken.
✓ Branch 101 → 102 taken 8 times.
✗ Branch 101 → 184 not taken.
✓ Branch 105 → 106 taken 8 times.
✗ Branch 105 → 192 not taken.
 16 o_beam_PxPyPzM->Save({beam_px, beam_py, beam_pz, o_beam_mass}, key);
300
3/8
✓ Branch 105 → 106 taken 8 times.
✗ Branch 105 → 192 not taken.
✓ Branch 106 → 107 taken 8 times.
✗ Branch 106 → 188 not taken.
✓ Branch 110 → 111 taken 8 times.
✗ Branch 110 → 113 not taken.
✗ Branch 192 → 193 not taken.
✗ Branch 192 → 195 not taken.
 16 o_target_PxPyPzM->Save({target_px, target_py, target_pz, target_mass}, key);
301
0/4
✗ Branch 143 → 144 not taken.
✗ Branch 143 → 146 not taken.
✗ Branch 170 → 171 not taken.
✗ Branch 170 → 173 not taken.
 8 }
302
303 ///////////////////////////////////////////////////////////////////////////////
304
305 1144 InclusiveKinematicsVars InclusiveKinematics::ComputeFromLepton(
306 vector_element_t const lepton_px,
307 vector_element_t const lepton_py,
308 vector_element_t const lepton_pz,
309 concurrent_key_t const key) const
310 {
311 1144 InclusiveKinematicsVars result;
312
313 1144 m_log->Trace("Reconstruct inclusive kinematics from lepton with p=({}, {}, {}), key={}", lepton_px, lepton_py, lepton_pz, key);
314
315 enum { px,
316 py,
317 pz,
318 m };
319 1144 auto beam = o_beam_PxPyPzM->Load(key);
320
1/2
✓ Branch 4 → 5 taken 1144 times.
✗ Branch 4 → 57 not taken.
 1144 auto target = o_target_PxPyPzM->Load(key);
321
322 ROOT::Math::PxPyPzMVector vec_beam(beam[px], beam[py], beam[pz], beam[m]);
323 ROOT::Math::PxPyPzMVector vec_target(target[px], target[py], target[pz], target[m]);
324 ROOT::Math::PxPyPzMVector vec_lepton(lepton_px, lepton_py, lepton_pz, beam[m]);
325
326 auto vec_q = vec_beam - vec_lepton;
327 1144 result.qx = vec_q.Px();
328 1144 result.qy = vec_q.Py();
329 1144 result.qz = vec_q.Pz();
330
1/2
✗ Branch 26 → 27 not taken.
✓ Branch 26 → 28 taken 1144 times.
 1144 result.qE = vec_q.E();
331 1144 result.Q2 = -1 * vec_q.M2();
332 1144 result.x = result.Q2 / (2 * vec_q.Dot(vec_target));
333
1/2
✓ Branch 32 → 33 taken 1144 times.
✗ Branch 32 → 53 not taken.
 1144 result.y = vec_target.Dot(vec_q) / vec_target.Dot(vec_beam);
334 ✗ result.W = (vec_beam + vec_target - vec_lepton).M();
335
1/2
✗ Branch 41 → 42 not taken.
✓ Branch 41 → 43 taken 1144 times.
 2288 result.nu = vec_target.Dot(vec_q) / target[m];
336 1144 result.beamPz = beam[pz];
337
1/2
✓ Branch 45 → 46 taken 1144 times.
✗ Branch 45 → 53 not taken.
 1144 result.targetM = target[m];
338
339
1/2
✓ Branch 46 → 47 taken 1144 times.
✗ Branch 46 → 49 not taken.
 1144 m_log->Trace("Result: Q2={} x={} W={}", result.Q2, result.x, result.W);
340
341 1144 return result;
342 }
343
344 ///////////////////////////////////////////////////////////////////////////////
345
346 4 void InclusiveKinematics::Stop()
347 {
348 4 }
349
350 }
351