GCC Code Coverage Report

Directory: ./
File: src/iguana/algorithms/physics/SingleHadronKinematics/Algorithm.cc
Date: 2025-01-05 09:03:17
Exec Total Coverage
Lines: 85 86 98.8%
Functions: 3 3 100.0%
Branches: 52 94 55.3%
Line Branch Exec Source
1 #include "Algorithm.h"
2 #include "TypeDefs.h"
3 #include "iguana/algorithms/physics/Tools.h"
4
5 #include <Math/Boost.h>
6
7 namespace iguana::physics {
8
9 REGISTER_IGUANA_ALGORITHM(SingleHadronKinematics, "physics::SingleHadronKinematics");
10
11 2 void SingleHadronKinematics::Start(hipo::banklist& banks)
12 {
13
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 b_particle = GetBankIndex(banks, "REC::Particle");
14
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 b_inc_kin = GetBankIndex(banks, "physics::InclusiveKinematics");
15
16 // create the output bank
17 // FIXME: generalize the groupid and itemid
18 2 auto result_schema = CreateBank(
19 banks,
20 2 b_result,
21 2 GetClassName(),
22 {
23 "pindex/S",
24 "pdg/I",
25 "z/D",
26 "PhPerp/D",
27 "MX2/D",
28 "xF/D",
29 "yB/D",
30 "phiH/D",
31 "xi/D"
32 },
33 0xF000,
34
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 7);
35
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 i_pindex = result_schema.getEntryOrder("pindex");
36
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 i_pdg = result_schema.getEntryOrder("pdg");
37
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 i_z = result_schema.getEntryOrder("z");
38
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 i_PhPerp = result_schema.getEntryOrder("PhPerp");
39
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 i_MX2 = result_schema.getEntryOrder("MX2");
40
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 i_xF = result_schema.getEntryOrder("xF");
41
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 i_yB = result_schema.getEntryOrder("yB");
42
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 i_phiH = result_schema.getEntryOrder("phiH");
43
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 i_xi = result_schema.getEntryOrder("xi");
44
45 // parse config file
46
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 ParseYAMLConfig();
47
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_pdgs = GetOptionSet<int>("hadron_list");
48
49 2 }
50
51 ///////////////////////////////////////////////////////////////////////////////
52
53 2000 void SingleHadronKinematics::Run(hipo::banklist& banks) const
54 {
55
1/2
✓ Branch 0 taken 2000 times.
✗ Branch 1 not taken.
 2000 auto& particle_bank = GetBank(banks, b_particle, "REC::Particle");
56
1/2
✓ Branch 0 taken 2000 times.
✗ Branch 1 not taken.
 4000 auto& inc_kin_bank = GetBank(banks, b_inc_kin, "physics::InclusiveKinematics");
57
1/2
✓ Branch 0 taken 2000 times.
✗ Branch 1 not taken.
 4000 auto& result_bank = GetBank(banks, b_result, GetClassName());
58
1/2
✓ Branch 0 taken 2000 times.
✗ Branch 1 not taken.
 2000 ShowBank(particle_bank, Logger::Header("INPUT PARTICLES"));
59
60
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()) {
61 1766 m_log->Debug("skip this event, since not all required banks have entries");
62 1766 return;
63 }
64
65 // get beam and target momenta
66 // FIXME: makes some assumptions about the beam; this should be generalized...
67 ROOT::Math::PxPyPzMVector p_beam(
68 0.0,
69 0.0,
70 234 inc_kin_bank.getDouble("beamPz", 0),
71 234 particle::mass.at(particle::electron));
72 ROOT::Math::PxPyPzMVector p_target(
73 0.0,
74 0.0,
75 0.0,
76 234 inc_kin_bank.getDouble("targetM", 0));
77
78 // get virtual photon momentum
79 ROOT::Math::PxPyPzEVector p_q(
80 234 inc_kin_bank.getDouble("qx", 0),
81 234 inc_kin_bank.getDouble("qy", 0),
82 234 inc_kin_bank.getDouble("qz", 0),
83 234 inc_kin_bank.getDouble("qE", 0));
84
85 // get additional inclusive variables
86 234 auto x = inc_kin_bank.getDouble("x", 0);
87 234 auto W = inc_kin_bank.getDouble("W", 0);
88
89 // boosts
90 234 ROOT::Math::Boost boost__qp((p_q + p_target).BoostToCM()); // CoM frame of target and virtual photon
91 234 ROOT::Math::Boost boost__breit((p_q + 2 * x * p_target).BoostToCM()); // Breit frame
92 234 auto p_q__qp = boost__qp(p_q);
93 234 auto p_q__breit = boost__breit(p_q);
94
95 // banks' row lists
96 234 auto const& particle_bank_rowlist = particle_bank.getRowList();
97 hipo::bank::rowlist::list_t result_bank_rowlist{};
98
1/2
✓ Branch 0 taken 234 times.
✗ Branch 1 not taken.
 234 result_bank.setRows(particle_bank.getRows());
99
100 // loop over ALL rows of `particle_bank`
101 // - we will calculate kinematics for rows in `particle_bank_rowlist`, and zero out all the other rows
102 // - we want the `result_bank` to have the same number of rows as `particle_bank` and the same ordering,
103 // so that banks which reference `particle_bank` rows can be used to reference `result_bank` rows too
104
2/2
✓ Branch 0 taken 1430 times.
✓ Branch 1 taken 234 times.
 1664 for(int row = 0; row < particle_bank.getRows(); row++) {
105
106 // if the particle is in `o_hadron_pdgs` AND the row is in `particle_bank`'s filtered row list
107 1430 if(auto pdg{particle_bank.getInt("pid", row)};
108
3/4
✓ Branch 0 taken 200 times.
✓ Branch 1 taken 1230 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 200 times.
 1430 o_hadron_pdgs.find(pdg) != o_hadron_pdgs.end() &&
109
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 200 times.
 200 std::find(particle_bank_rowlist.begin(), particle_bank_rowlist.end(), row) != particle_bank_rowlist.end()) {
110
111 // hadron momentum
112 auto p_Ph = ROOT::Math::PxPyPzMVector(
113
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 200 times.
 200 particle_bank.getFloat("px", row),
114 200 particle_bank.getFloat("py", row),
115 200 particle_bank.getFloat("pz", row),
116
1/2
✓ Branch 0 taken 200 times.
✗ Branch 1 not taken.
 200 particle::mass.at(static_cast<particle::PDG>(pdg)));
117
1/2
✓ Branch 0 taken 200 times.
✗ Branch 1 not taken.
 200 auto p_Ph__qp = boost__qp(p_Ph);
118
1/2
✓ Branch 0 taken 200 times.
✗ Branch 1 not taken.
 200 auto p_Ph__breit = boost__breit(p_Ph);
119
120 // calculate z
121
1/2
✓ Branch 0 taken 200 times.
✗ Branch 1 not taken.
 200 double z = p_target.Dot(p_Ph) / p_target.Dot(p_q);
122
123 // calculate PhPerp
124
2/4
✓ Branch 0 taken 200 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 200 times.
✗ Branch 3 not taken.
 200 auto opt_PhPerp = tools::RejectVector(p_Ph.Vect(), p_q.Vect());
125
1/2
✓ Branch 0 taken 200 times.
✗ Branch 1 not taken.
 200 double PhPerp = opt_PhPerp.has_value() ? opt_PhPerp.value().R() : tools::UNDEF;
126
127 // calculate MX2
128
1/2
✓ Branch 0 taken 200 times.
✗ Branch 1 not taken.
 200 double MX2 = (p_target + p_q - p_Ph).M2();
129
130 // calculate xF
131
1/2
✓ Branch 0 taken 200 times.
✗ Branch 1 not taken.
 200 double xF = 2 * p_Ph__qp.Vect().Dot(p_q__qp.Vect()) / (W * p_q__qp.Vect().R());
132
133 // calculate yB
134
2/4
✓ Branch 0 taken 200 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 200 times.
✗ Branch 3 not taken.
 400 double yB = tools::ParticleRapidity(p_Ph__breit, p_q__breit.Vect()).value_or(tools::UNDEF);
135
136 // calculate phiH
137
1/2
✓ Branch 0 taken 200 times.
✗ Branch 1 not taken.
 200 double phiH = tools::PlaneAngle(
138 200 p_q.Vect(),
139 200 p_beam.Vect(),
140 200 p_q.Vect(),
141 ✗ p_Ph.Vect()).value_or(tools::UNDEF);
142
143 // calculate xi
144
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 199 times.
 200 double xi = p_q.Dot(p_Ph) / p_target.Dot(p_q);
145
146 // put this particle in `result_bank`'s row list
147 result_bank_rowlist.push_back(row);
148
149 // fill the bank
150 200 result_bank.putShort(i_pindex, row, static_cast<int16_t>(row));
151 200 result_bank.putInt(i_pdg, row, pdg);
152 200 result_bank.putDouble(i_z, row, z);
153 200 result_bank.putDouble(i_PhPerp, row, PhPerp);
154 200 result_bank.putDouble(i_MX2, row, MX2);
155 200 result_bank.putDouble(i_xF, row, xF);
156 200 result_bank.putDouble(i_yB, row, yB);
157 200 result_bank.putDouble(i_phiH, row, phiH);
158 200 result_bank.putDouble(i_xi, row, xi);
159 }
160 else {
161 // zero the row
162 1230 result_bank.putShort(i_pindex, row, static_cast<int16_t>(row));
163 1230 result_bank.putInt(i_pdg, row, pdg);
164 1230 result_bank.putDouble(i_z, row, 0);
165 1230 result_bank.putDouble(i_PhPerp, row, 0);
166 1230 result_bank.putDouble(i_MX2, row, 0);
167 1230 result_bank.putDouble(i_xF, row, 0);
168 1230 result_bank.putDouble(i_yB, row, 0);
169 1230 result_bank.putDouble(i_phiH, row, 0);
170 1230 result_bank.putDouble(i_xi, row, 0);
171 }
172 }
173
174 // apply the filtered rowlist to `result_bank`
175
2/4
✓ Branch 0 taken 234 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 234 times.
✗ Branch 3 not taken.
 234 result_bank.getMutableRowList().setList(result_bank_rowlist);
176
177
4/8
✓ Branch 0 taken 234 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 234 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 140 times.
✓ Branch 5 taken 94 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 468 ShowBank(result_bank, Logger::Header("CREATED BANK"));
178 }
179
180 ///////////////////////////////////////////////////////////////////////////////
181
182 1 void SingleHadronKinematics::Stop()
183 {
184 1 }
185
186 }
187