GCC Code Coverage Report

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