GCC Code Coverage Report

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