GCC Code Coverage Report

Directory: ./
File: src/iguana/algorithms/clas12/FiducialFilter/Algorithm.cc
Date: 2025-01-05 09:03:17
Exec Total Coverage
Lines: 143 159 89.9%
Functions: 9 9 100.0%
Branches: 125 180 69.4%
Line Branch Exec Source
1 #include "Algorithm.h"
2 #include "Pass1CutData.h"
3
4 namespace iguana::clas12 {
5
6 REGISTER_IGUANA_ALGORITHM(FiducialFilter);
7
8 2 void FiducialFilter::Start(hipo::banklist& banks)
9 {
10
11 // Read YAML config file
12 2 ParseYAMLConfig();
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_traj = GetBankIndex(banks, "REC::Traj");
15
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 b_config = GetBankIndex(banks, "RUN::config");
16
17
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 4 o_pass = GetCachedOption<int>("pass").value_or(1);
18
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if(o_pass!=1){
19 ✗ m_log->Warn("FiducialFilter only contains fiducial cuts for pass1...we will default to using those...");
20 }
21
22 2 }
23
24 2000 void FiducialFilter::Run(hipo::banklist& banks) const {
25 // get the banks
26
1/2
✓ Branch 0 taken 2000 times.
✗ Branch 1 not taken.
 2000 auto& particleBank = GetBank(banks, b_particle, "REC::Particle");
27
1/2
✓ Branch 0 taken 2000 times.
✗ Branch 1 not taken.
 2000 auto& trajBank = GetBank(banks, b_traj, "REC::Traj");
28
1/2
✓ Branch 0 taken 2000 times.
✗ Branch 1 not taken.
 2000 auto& configBank = GetBank(banks, b_config, "RUN::config");
29 2000 auto torus = configBank.getFloat("torus", 0);
30
31 // dump the bank
32
1/2
✓ Branch 0 taken 2000 times.
✗ Branch 1 not taken.
 2000 ShowBank(particleBank, Logger::Header("INPUT PARTICLES"));
33
34 // get a pindex'd map of the REC::Traj data
35 2000 auto traj_map = GetTrajMap(trajBank);
36
37 // filter the input bank for requested PDG code(s)
38
4/8
✓ Branch 0 taken 2000 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2000 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2000 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2000 times.
✗ Branch 7 not taken.
 4000 particleBank.getMutableRowList().filter([this, &traj_map, torus](hipo::bank& bank, int row) {
39 // Check if this particle has a REC::Traj component
40
2/2
✓ Branch 0 taken 4584 times.
✓ Branch 1 taken 4907 times.
 18982 if (auto it{traj_map.find(row)}; it != traj_map.end()) {
41 4584 auto traj_row = it->second;
42 // Reject particles with undefined sector
43
2/2
✓ Branch 0 taken 3274 times.
✓ Branch 1 taken 1310 times.
 4584 if(traj_row.sector==-1) return false;
44 3274 auto pid = bank.getInt("pid", row);
45 3274 return Filter(traj_row, torus, pid);
46 }
47 else return false; // particle not in REC::traj
48 });
49
50 // dump the modified bank
51
2/4
✓ Branch 0 taken 2000 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2000 times.
✗ Branch 3 not taken.
 4000 ShowBank(particleBank, Logger::Header("OUTPUT PARTICLES"));
52 2000 }
53
54 3274 bool FiducialFilter::Filter(FiducialFilter::traj_row_data const traj_row, float const torus, int const pid) const
55 {
56
57
2/4
✓ Branch 0 taken 3274 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 3274 times.
 3274 if(torus!=1&&torus!=-1){
58 ✗ m_log->Warn("torus={}...value must be either -1 or 1, otherwise fiducial cuts are not defined...filtering out all particles...",torus);
59 ✗ return false;
60 }
61
62
2/2
✓ Branch 0 taken 244 times.
✓ Branch 1 taken 3030 times.
 3274 if(pid==11) return DC_fiducial_cut_XY_pass1(traj_row, torus, pid);
63
4/4
✓ Branch 0 taken 636 times.
✓ Branch 1 taken 2394 times.
✓ Branch 2 taken 302 times.
✓ Branch 3 taken 334 times.
 3030 else if(pid==211 || pid==-211 || pid==2212){
64
1/2
✓ Branch 0 taken 2696 times.
✗ Branch 1 not taken.
 2696 if(torus<0) return DC_fiducial_cut_theta_phi_pass1(traj_row, torus, pid);
65 ✗ else if(torus>0) return DC_fiducial_cut_XY_pass1(traj_row, torus, pid);
66 else return false;
67 }
68 return true;
69 }
70
71
72
73 244 bool FiducialFilter::DC_fiducial_cut_XY_pass1(FiducialFilter::traj_row_data const traj_row, float const torus, int const pid) const
74 {
75
76
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 244 times.
 244 const auto minparams = ((torus<0) ? minparams_in_XY_pass1 : minparams_out_XY_pass1);
77 const auto maxparams = ((torus<0) ? maxparams_in_XY_pass1 : maxparams_out_XY_pass1);
78 double X=0;
79 double Y=0;
80
2/2
✓ Branch 0 taken 710 times.
✓ Branch 1 taken 228 times.
 938 for(int r = 0 ; r < 3; r++){
81 X=0;
82 Y=0;
83
3/3
✓ Branch 0 taken 244 times.
✓ Branch 1 taken 234 times.
✓ Branch 2 taken 232 times.
 710 switch(r){
84 244 case 0:
85 244 X = traj_row.x1;
86 244 Y = traj_row.y1;
87 244 break;
88 234 case 1:
89 234 X = traj_row.x2;
90 234 Y = traj_row.y2;
91 234 break;
92 232 case 2:
93 232 X = traj_row.x3;
94 232 Y = traj_row.y3;
95 232 break;
96 }
97
98 710 int sector = traj_row.sector;
99
100
2/2
✓ Branch 0 taken 162 times.
✓ Branch 1 taken 548 times.
 710 if(sector == 2)
101 {
102 162 const double X_new = X * std::cos(-60 * M_PI / 180) - Y * std::sin(-60 * M_PI / 180);
103 162 Y = X * std::sin(-60 * M_PI / 180) + Y * std::cos(-60 * M_PI / 180);
104 X = X_new;
105 }
106
107
2/2
✓ Branch 0 taken 140 times.
✓ Branch 1 taken 570 times.
 710 if(sector == 3)
108 {
109 140 const double X_new = X * std::cos(-120 * M_PI / 180) - Y * std::sin(-120 * M_PI / 180);
110 140 Y = X * std::sin(-120 * M_PI / 180) + Y * std::cos(-120 * M_PI / 180);
111 X = X_new;
112 }
113
114
2/2
✓ Branch 0 taken 62 times.
✓ Branch 1 taken 648 times.
 710 if(sector == 4)
115 {
116 62 const double X_new = X * std::cos(-180 * M_PI / 180) - Y * std::sin(-180 * M_PI / 180);
117 62 Y = X * std::sin(-180 * M_PI / 180) + Y * std::cos(-180 * M_PI / 180);
118 X = X_new;
119 }
120
121
2/2
✓ Branch 0 taken 102 times.
✓ Branch 1 taken 608 times.
 710 if(sector == 5)
122 {
123 102 const double X_new = X * std::cos(120 * M_PI / 180) - Y * std::sin(120 * M_PI / 180);
124 102 Y = X * std::sin(120 * M_PI / 180) + Y * std::cos(120 * M_PI / 180);
125 X = X_new;
126 }
127
128
2/2
✓ Branch 0 taken 104 times.
✓ Branch 1 taken 606 times.
 710 if(sector == 6)
129 {
130 104 const double X_new = X * std::cos(60 * M_PI / 180) - Y * std::sin(60 * M_PI / 180);
131 104 Y = X * std::sin(60 * M_PI / 180) + Y * std::cos(60 * M_PI / 180);
132 X = X_new;
133 }
134
135
136 int this_pid = 0;
137
138
1/7
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 710 times.
✗ Branch 6 not taken.
 710 switch (pid)
139 {
140 case 11:
141 this_pid = 0;
142 break;
143 ✗ case 2212:
144 this_pid = 1;
145 ✗ break;
146 ✗ case 211:
147 this_pid = 2;
148 ✗ break;
149 ✗ case -211:
150 this_pid = 3;
151 ✗ break;
152 ✗ case 321:
153 this_pid = 4;
154 ✗ break;
155 ✗ case -321:
156 this_pid = 5;
157 ✗ break;
158 default:
159 return false;
160 break;
161 }
162 710 double calc_min = minparams[this_pid][sector - 1][r][0] + minparams[this_pid][sector - 1][r][1] * X;
163 710 double calc_max = maxparams[this_pid][sector - 1][r][0] + maxparams[this_pid][sector - 1][r][1] * X;
164
2/4
✓ Branch 0 taken 710 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 710 times.
✗ Branch 3 not taken.
 710 if(std::isnan(calc_min)||std::isnan(calc_max)) return false;
165
4/4
✓ Branch 0 taken 702 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 694 times.
✓ Branch 3 taken 8 times.
 710 if((Y<calc_min) || (Y>calc_max)) { return false;}
166 }
167 return true;
168 }
169
170 2696 bool FiducialFilter::DC_fiducial_cut_theta_phi_pass1(FiducialFilter::traj_row_data const traj_row, float const torus, int const pid) const{
171
172
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2696 times.
 2696 const auto minparams = ((torus<0) ? minparams_in_theta_phi_pass1 : minparams_out_theta_phi_pass1);
173 const auto maxparams = ((torus<0) ? maxparams_in_theta_phi_pass1 : maxparams_out_theta_phi_pass1);
174 double theta_DCr = 5000;
175 double phi_DCr_raw = 5000;
176 double x=0;
177 double y=0;
178 double z=0;
179
2/2
✓ Branch 0 taken 8012 times.
✓ Branch 1 taken 2656 times.
 10668 for(int r = 0; r<3; r++){
180 x=0;y=0;z=0;
181
3/3
✓ Branch 0 taken 2696 times.
✓ Branch 1 taken 2658 times.
✓ Branch 2 taken 2658 times.
 8012 switch(r){
182 2696 case 0:
183 2696 x=traj_row.x1;
184 2696 y=traj_row.y1;
185 2696 z=traj_row.z1;
186 2696 break;
187 2658 case 1:
188 2658 x=traj_row.x2;
189 2658 y=traj_row.y2;
190 2658 z=traj_row.z2;
191 2658 break;
192 2658 case 2:
193 2658 x=traj_row.x3;
194 2658 y=traj_row.y3;
195 2658 z=traj_row.z3;
196 2658 break;
197 }
198 8012 int sector = traj_row.sector;
199
200 8012 theta_DCr = 180 / M_PI * acos(z / sqrt(pow(x,2) + pow(y,2) + pow(z,2)));
201 8012 phi_DCr_raw = 180 / M_PI * atan2(y / sqrt(pow(x,2) + pow(y,2) + pow(z,2)),
202 8012 x /sqrt(pow(x,2) + pow(y,2) + pow(z,2)));
203
204 double phi_DCr = 5000;
205
206
2/2
✓ Branch 0 taken 6676 times.
✓ Branch 1 taken 1336 times.
 8012 if (sector == 1) phi_DCr = phi_DCr_raw;
207
2/2
✓ Branch 0 taken 1284 times.
✓ Branch 1 taken 5392 times.
 6676 if (sector == 2) phi_DCr = phi_DCr_raw - 60;
208
2/2
✓ Branch 0 taken 1444 times.
✓ Branch 1 taken 5284 times.
 6728 if (sector == 3) phi_DCr = phi_DCr_raw - 120;
209
4/4
✓ Branch 0 taken 1244 times.
✓ Branch 1 taken 5324 times.
✓ Branch 2 taken 632 times.
✓ Branch 3 taken 612 times.
 6568 if (sector == 4 && phi_DCr_raw > 0) phi_DCr = phi_DCr_raw - 180;
210
2/2
✓ Branch 0 taken 612 times.
✓ Branch 1 taken 632 times.
 1244 if (sector == 4 && phi_DCr_raw < 0) phi_DCr = phi_DCr_raw + 180;
211
2/2
✓ Branch 0 taken 1362 times.
✓ Branch 1 taken 6650 times.
 8012 if (sector == 5) phi_DCr = phi_DCr_raw + 120;
212
2/2
✓ Branch 0 taken 1342 times.
✓ Branch 1 taken 5308 times.
 6650 if (sector == 6) phi_DCr = phi_DCr_raw + 60;
213
214 int this_pid = 0;
215
216
3/7
✗ Branch 0 not taken.
✓ Branch 1 taken 902 times.
✓ Branch 2 taken 5082 times.
✓ Branch 3 taken 2028 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 8012 switch (pid)
217 {
218 case 11: this_pid = 0; break;
219 902 case 2212: this_pid = 1; break;
220 5082 case 211: this_pid = 2; break;
221 2028 case -211: this_pid = 3; break;
222 ✗ case 321: this_pid = 4; break;
223 ✗ case -321: this_pid = 5; break;
224 default: return false; break;
225 }
226
227
228
229 8012 double calc_phi_min = minparams[this_pid][sector - 1][r][0] + minparams[this_pid][sector - 1][r][1] * std::log(theta_DCr)
230 8012 + minparams[this_pid][sector - 1][r][2] * theta_DCr + minparams[this_pid][sector - 1][r][3] * theta_DCr * theta_DCr;
231
232 8012 double calc_phi_max = maxparams[this_pid][sector - 1][r][0] + maxparams[this_pid][sector - 1][r][1] * std::log(theta_DCr)
233 8012 + maxparams[this_pid][sector - 1][r][2] * theta_DCr + maxparams[this_pid][sector - 1][r][3] * theta_DCr * theta_DCr;
234
235
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 8012 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 8012 times.
 8012 if(std::isnan(calc_phi_min)||std::isnan(calc_phi_max)) return false;
236
4/4
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 8000 times.
✓ Branch 2 taken 28 times.
✓ Branch 3 taken 7972 times.
 8012 if((phi_DCr < calc_phi_min) || (phi_DCr > calc_phi_max)) return false;
237 }
238 return true;
239 }
240
1/2
✓ Branch 0 taken 3000 times.
✗ Branch 1 not taken.
 3000 std::map<int, FiducialFilter::traj_row_data> FiducialFilter::GetTrajMap(hipo::bank const &bank)
241 {
242 std::map<int, FiducialFilter::traj_row_data> traj_map;
243
244
3/4
✓ Branch 0 taken 3000 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 112359 times.
✓ Branch 3 taken 3000 times.
 115359 for(auto const& row : bank.getRowList()){
245 112359 auto pindex = bank.getShort("pindex",row);
246 112359 auto x = bank.getFloat("x",row);
247 112359 auto y = bank.getFloat("y",row);
248 112359 auto z = bank.getFloat("z",row);
249 112359 auto layer = bank.getInt("layer",row);
250
251 // Ensure an entry exists in the map for the given pindex
252
2/2
✓ Branch 0 taken 8295 times.
✓ Branch 1 taken 104064 times.
 112359 if (traj_map.find(pindex) == traj_map.end()) {
253
1/2
✓ Branch 0 taken 8295 times.
✗ Branch 1 not taken.
 8295 traj_map[pindex] = FiducialFilter::traj_row_data();
254 }
255
256
4/4
✓ Branch 0 taken 10638 times.
✓ Branch 1 taken 5412 times.
✓ Branch 2 taken 5412 times.
✓ Branch 3 taken 90897 times.
 112359 switch(layer){
257 10638 case 6: // first DC
258
1/2
✓ Branch 0 taken 10638 times.
✗ Branch 1 not taken.
 10638 traj_map[pindex].x1 = x;
259
1/2
✓ Branch 0 taken 10638 times.
✗ Branch 1 not taken.
 10638 traj_map[pindex].y1 = y;
260
1/2
✓ Branch 0 taken 10638 times.
✗ Branch 1 not taken.
 10638 traj_map[pindex].z1 = z;
261 10638 break;
262 5412 case 18: // second DC
263
1/2
✓ Branch 0 taken 5412 times.
✗ Branch 1 not taken.
 5412 traj_map[pindex].x2 = x;
264
1/2
✓ Branch 0 taken 5412 times.
✗ Branch 1 not taken.
 5412 traj_map[pindex].y2 = y;
265
1/2
✓ Branch 0 taken 5412 times.
✗ Branch 1 not taken.
 5412 traj_map[pindex].z2 = z;
266 // Determine Sector from the center of the DC
267
5/10
✓ Branch 0 taken 5412 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 5412 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 5412 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 5412 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 5412 times.
✗ Branch 9 not taken.
 5412 traj_map[pindex].sector = determineSectorDC(traj_map[pindex].x2,traj_map[pindex].y2,traj_map[pindex].z2);
268 5412 break;
269 5412 case 36: // third DC
270
1/2
✓ Branch 0 taken 5412 times.
✗ Branch 1 not taken.
 5412 traj_map[pindex].x3 = x;
271
1/2
✓ Branch 0 taken 5412 times.
✗ Branch 1 not taken.
 5412 traj_map[pindex].y3 = y;
272
1/2
✓ Branch 0 taken 5412 times.
✗ Branch 1 not taken.
 5412 traj_map[pindex].z3 = z;
273 5412 break;
274 }
275 }
276
277
278
279 3000 return traj_map;
280 }
281
282 5412 int FiducialFilter::determineSectorDC(float x, float y, float z)
283 {
284 5412 float phi = 180 / M_PI * atan2(y / sqrt(pow(x,2) + pow(y,2) + pow(z,2)),
285 5412 x /sqrt(pow(x,2) + pow(y,2) + pow(z,2)));
286
4/4
✓ Branch 0 taken 3147 times.
✓ Branch 1 taken 2265 times.
✓ Branch 2 taken 2256 times.
✓ Branch 3 taken 891 times.
 5412 if(phi<30 && phi>=-30){return 1;}
287
4/4
✓ Branch 0 taken 3126 times.
✓ Branch 1 taken 1395 times.
✓ Branch 2 taken 2256 times.
✓ Branch 3 taken 870 times.
 4521 else if(phi<90 && phi>=30){return 2;}
288
4/4
✓ Branch 0 taken 3222 times.
✓ Branch 1 taken 429 times.
✓ Branch 2 taken 2256 times.
✓ Branch 3 taken 966 times.
 3651 else if(phi<150 && phi>=90){return 3;}
289
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 429 times.
✓ Branch 2 taken 1866 times.
✓ Branch 3 taken 390 times.
 2685 else if(phi>=150 || phi<-150){return 4;}
290
2/2
✓ Branch 0 taken 927 times.
✓ Branch 1 taken 939 times.
 1866 else if(phi<-90 && phi>=-150){return 5;}
291
1/2
✓ Branch 0 taken 927 times.
✗ Branch 1 not taken.
 927 else if(phi<-30 && phi>=-90){return 6;}
292
293 return -1;
294
295 }
296
297 1 void FiducialFilter::Stop()
298 {
299 1 }
300
301 }
302