173 def event(self):
174 """
175 Called at the beginning of each event.
176 """
177 b2.B2DEBUG(10, "---- Processing new event ----")
178
179
181
182
183 for candidate in candidate_list:
184
185 p_list = get_object_list(candidate.getFinalStateDaughters())
186
187
188 n_nodes = len(p_list)
189
190
191 masses = np.array([abs(p.getPDGCode()) for p in p_list])
192
193
194 graFEI_nFSP = n_nodes
195 graFEI_nPhotons_preFit = (masses == 22).sum()
196 graFEI_nCharged_preFit = graFEI_nFSP - graFEI_nPhotons_preFit
197 graFEI_nElectrons_preFit = (masses == 11).sum()
198 graFEI_nMuons_preFit = (masses == 13).sum()
199 graFEI_nPions_preFit = (masses == 211).sum()
200 graFEI_nKaons_preFit = (masses == 321).sum()
201 graFEI_nProtons_preFit = (masses == 2212).sum()
202 graFEI_nLeptons_preFit = graFEI_nElectrons_preFit + graFEI_nMuons_preFit
203 graFEI_nOthers_preFit = graFEI_nCharged_preFit - \
204 (graFEI_nLeptons_preFit + graFEI_nPions_preFit + graFEI_nKaons_preFit + graFEI_nProtons_preFit)
205
206 candidate.addExtraInfo("graFEI_nFSP", graFEI_nFSP)
207 candidate.addExtraInfo("graFEI_nCharged_preFit", graFEI_nCharged_preFit)
208 candidate.addExtraInfo("graFEI_nPhotons_preFit", graFEI_nPhotons_preFit)
209 candidate.addExtraInfo("graFEI_nElectrons_preFit", graFEI_nElectrons_preFit)
210 candidate.addExtraInfo("graFEI_nMuons_preFit", graFEI_nMuons_preFit)
211 candidate.addExtraInfo("graFEI_nPions_preFit", graFEI_nPions_preFit)
212 candidate.addExtraInfo("graFEI_nKaons_preFit", graFEI_nKaons_preFit)
213 candidate.addExtraInfo("graFEI_nProtons_preFit", graFEI_nProtons_preFit)
214 candidate.addExtraInfo("graFEI_nLeptons_preFit", graFEI_nLeptons_preFit)
215 candidate.addExtraInfo("graFEI_nOthers_preFit", graFEI_nOthers_preFit)
216
217
218 if n_nodes < 2:
219 b2.B2WARNING(
220 f"Skipping candidate with {n_nodes} reconstructed FSPs"
221 )
222
223 continue
224
225
226 x_nodes = np.empty((n_nodes, len(self.node_features)))
227 x_dis = np.empty((n_nodes, len(self.discarded_features)))
228
229
230 for p, particle in enumerate(p_list):
231 for f, feat in enumerate(self.node_features):
232 feat = feat[feat.find("feat_") + 5:]
233 x_nodes[p, f] = vm.evaluate(feat, particle)
234 for f, feat in enumerate(self.discarded_features):
235 feat = feat[feat.find("feat_") + 5:]
236 x_dis[p, f] = vm.evaluate(feat, particle)
237 b2.B2DEBUG(11, "Node features:\n", x_nodes)
238
239
240 x_edges = (compute_edge_features(self.edge_features, self.node_features + self.discarded_features,
241 np.concatenate([x_nodes, x_dis], axis=1)) if self.edge_features != [] else [])
242 edge_index = torch.tensor(list(itertools.permutations(range(n_nodes), 2)), dtype=torch.long)
243 b2.B2DEBUG(11, "Edge features:\n", x_edges)
244
245
246 x_global = (
247 np.array([[n_nodes]], dtype=float)
248 if self.glob_features != []
249 else []
250 )
251 b2.B2DEBUG(11, "Global features:\n", x_global)
252
253
254 torch_batch = torch.zeros(size=[n_nodes], dtype=torch.long)
255
256
257 np.nan_to_num(x_nodes, copy=False)
258 np.nan_to_num(x_edges, copy=False)
259 np.nan_to_num(x_global, copy=False)
260
261
262 if self.normalize is not None:
264 self.normalize,
265 self.node_features,
266 x_nodes,
267 self.edge_features,
268 x_edges,
269 self.glob_features,
270 x_global,
271 )
272
273
274 x = torch.tensor(x_nodes, dtype=torch.float).to(self.device)
275 edge_index = edge_index.t().contiguous().to(self.device)
276 edge_attr = torch.tensor(x_edges, dtype=torch.float).to(self.device)
277 u = torch.tensor(x_global, dtype=torch.float).to(self.device)
278 torch_batch = torch_batch.to(self.device)
279
280
281 batch = Batch(
282 x=x, edge_index=edge_index, edge_attr=edge_attr, u=u, batch=torch_batch
283 )
284
285
286 with torch.no_grad():
287 x_pred, e_pred, u_pred = self.model(batch)
288
289
290
291
292
293
294
295 edge_probs = torch.softmax(e_pred, dim=1)
296 edge_probability, predicted_LCA = edge_probs.max(dim=1)
297
298
299 mass_probs = torch.softmax(x_pred, dim=1)
300 mass_probability, predicted_masses = mass_probs.max(dim=1)
301 b2.B2DEBUG(10, "Predicted mass classes:\n", predicted_masses)
302 b2.B2DEBUG(11, "Mass class probabilities:\n", mass_probability)
303
304
305 graFEI_nPhotons_postFit = (predicted_masses == 6).sum()
306 graFEI_nCharged_postFit = graFEI_nFSP - graFEI_nPhotons_postFit
307 graFEI_nElectrons_postFit = (predicted_masses == 1).sum()
308 graFEI_nMuons_postFit = (predicted_masses == 2).sum()
309 graFEI_nPions_postFit = (predicted_masses == 3).sum()
310 graFEI_nKaons_postFit = (predicted_masses == 4).sum()
311 graFEI_nProtons_postFit = (predicted_masses == 5).sum()
312 graFEI_nLeptons_postFit = graFEI_nElectrons_postFit + graFEI_nMuons_postFit
313 graFEI_nOthers_postFit = (predicted_masses == 0).sum()
314
315
316 edge_probability_square = torch.sparse_coo_tensor(
317 edge_index, edge_probability
318 ).to_dense()
319 predicted_LCA_square = torch.sparse_coo_tensor(
320 edge_index, predicted_LCA, dtype=int
321 ).to_dense()
322 b2.B2DEBUG(10, "Predicted LCA:\n", predicted_LCA_square)
323 b2.B2DEBUG(11, "Edge class probabilities:\n", edge_probability_square)
324
325
326 edge_probability_unique = edge_probability_square[
327 edge_probability_square.tril(diagonal=-1) > 0
328 ]
329
330
331 predicted_matched = np.array(
332 [False if torch.all(i == 0) else True for i in predicted_LCA_square]
333 )
334 b2.B2DEBUG(10, "Predicted matched particles:\n", predicted_matched)
335
336 predicted_matched_noPhotons = predicted_matched[masses != 22]
337
338
339 graFEI_nPredictedUnmatched = (~predicted_matched).sum()
340 graFEI_nPredictedUnmatched_noPhotons = (
341 (~predicted_matched_noPhotons).sum()
342 if predicted_matched_noPhotons.size != 0
343 else 0
344 )
345
346
347 predicted_LCA_square_matched = predicted_LCA_square[predicted_matched]
348 predicted_LCA_square_matched = predicted_LCA_square_matched[:, predicted_matched]
349
350
351 predicted_masses_matched = predicted_masses[predicted_matched]
352
353
354 graFEI_validTree = 0
355 if not torch.all(predicted_LCA_square == 0):
356 try:
358 graFEI_validTree = 1
359 except InvalidLCAMatrix:
360 pass
361
362
363 graFEI_goodEvent = 0
364 if graFEI_validTree:
365
366 good_decay, root_level, sig_side_fsps = select_good_decay(predicted_LCA_square_matched,
367 predicted_masses_matched,
368 self.sig_side_lcas,
369 self.sig_side_masses)
370 graFEI_goodEvent = int((self.max_level == root_level) and good_decay)
371
372 if graFEI_goodEvent:
373
374 p_list_matched = list(np.array(p_list)[predicted_matched])
375 for i, particle in enumerate(p_list_matched):
376 if i in sig_side_fsps:
377 particle.addExtraInfo("graFEI_sigSide", 1)
378 else:
379 particle.addExtraInfo("graFEI_sigSide", 0)
380
381 b2.B2DEBUG(11, "This LCA describes a valid tree")
382 b2.B2DEBUG(
383 11,
384 "Predicted LCA on matched particles:\n",
385 predicted_LCA_square_matched,
386 )
387 b2.B2DEBUG(11, "Adjacency matrix:\n", adjacency)
388
389
390 for particle in p_list:
391 if not particle.hasExtraInfo("graFEI_sigSide"):
392 particle.addExtraInfo("graFEI_sigSide", -1)
393
394
395 graFEI_probEdgeProd = edge_probability_unique.prod().item()
396 graFEI_probEdgeMean = edge_probability_unique.mean().item()
397 graFEI_probEdgeGeom = torch.pow(edge_probability_unique.prod(), 1/n_nodes).item()
398
399
400 candidate.addExtraInfo("graFEI_probEdgeProd", graFEI_probEdgeProd)
401 candidate.addExtraInfo("graFEI_probEdgeMean", graFEI_probEdgeMean)
402 candidate.addExtraInfo("graFEI_probEdgeGeom", graFEI_probEdgeGeom)
403 candidate.addExtraInfo("graFEI_validTree", graFEI_validTree)
404 candidate.addExtraInfo("graFEI_goodEvent", graFEI_goodEvent)
405 candidate.addExtraInfo("graFEI_nPhotons_postFit", graFEI_nPhotons_postFit)
406 candidate.addExtraInfo("graFEI_nCharged_postFit", graFEI_nCharged_postFit)
407 candidate.addExtraInfo("graFEI_nElectrons_postFit", graFEI_nElectrons_postFit)
408 candidate.addExtraInfo("graFEI_nMuons_postFit", graFEI_nMuons_postFit)
409 candidate.addExtraInfo("graFEI_nPions_postFit", graFEI_nPions_postFit)
410 candidate.addExtraInfo("graFEI_nKaons_postFit", graFEI_nKaons_postFit)
411 candidate.addExtraInfo("graFEI_nProtons_postFit", graFEI_nProtons_postFit)
412 candidate.addExtraInfo("graFEI_nLeptons_postFit", graFEI_nLeptons_postFit)
413 candidate.addExtraInfo("graFEI_nOthers_postFit", graFEI_nOthers_postFit)
414 candidate.addExtraInfo("graFEI_nPredictedUnmatched", graFEI_nPredictedUnmatched)
415 candidate.addExtraInfo("graFEI_nPredictedUnmatched_noPhotons", graFEI_nPredictedUnmatched_noPhotons)
416
417
418 if self.storeTrueInfo:
419
420 parentID = np.array([vm.evaluate("ancestorBIndex", p) for p in p_list], dtype=int)
421 b2.B2DEBUG(10, "Ancestor true ID:\n", parentID)
422
423
424 p_indices = np.array(
425 [
426 p.getMCParticle().getArrayIndex() if parentID[i] >= 0 else -1
427 for (i, p) in enumerate(p_list)
428 ]
429 )
430
431 p_masses = masses_to_classes(
432 np.array(
433 [
434 p.getMCParticle().getPDG() if parentID[i] >= 0 else -1
435 for (i, p) in enumerate(p_list)
436 ]
437 )
438 )
439 b2.B2DEBUG(10, "True mass classes:\n", p_masses)
440
441 evt_primary = np.array(
442 [
443 p.getMCParticle().isPrimaryParticle()
444 if parentID[i] >= 0
445 else False
446 for (i, p) in enumerate(p_list)
447 ]
448 )
449 b2.B2DEBUG(10, "Is primary particle:\n", evt_primary)
450
451
452
453
454 B_indices = parentID[np.logical_and(evt_primary, predicted_matched)]
455 b2.B2DEBUG(
456 10, "Ancestor ID of predicted matched particles:\n", B_indices
457 )
458 B_indices = list(set(B_indices))
459
460
461 graFEI_truth_perfectLCA = 0
462 graFEI_truth_isSemileptonic = -1
463 graFEI_truth_nFSP = -1
464 graFEI_truth_perfectMasses = int((predicted_masses.numpy() == p_masses).all()
465 )
466 graFEI_truth_nPhotons = (p_masses == 6).sum()
467 graFEI_truth_nElectrons = (p_masses == 1).sum()
468 graFEI_truth_nMuons = (p_masses == 2).sum()
469 graFEI_truth_nPions = (p_masses == 3).sum()
470 graFEI_truth_nKaons = (p_masses == 4).sum()
471 graFEI_truth_nProtons = (p_masses == 5).sum()
472 graFEI_truth_nOthers = (p_masses == 0).sum()
473
474
476
477
478 if self.mc_particle == "Upsilon(4S):MC" and gen_list.getListSize() > 1:
479 b2.B2WARNING(
480 f"Found {gen_list.getListSize()} true Upsilon(4S) in the generated MC (??)")
481
482 if gen_list.getListSize() > 0:
483
484 for genP in gen_list.obj():
485 mcp = genP.getMCParticle()
486
487
488 if self.mc_particle != "Upsilon(4S):MC" and len(B_indices) != 1:
489 break
490
491
492 array_index = mcp.getArrayIndex()
493
494
495 if self.mc_particle != "Upsilon(4S):MC" and array_index != B_indices[0]:
496 continue
497
498
499 (
500 leaf_hist,
501 levels,
502 _,
503 _,
504 semilep_flag,
505 ) = write_hist(
506 particle=mcp,
507 leaf_hist={},
508 levels={},
509 hist=[],
510 pdg={},
511 leaf_pdg={},
512 semilep_flag=False,
513 )
514
515
516 if len(leaf_hist) < 2:
517 continue
518
519
520 true_LCA_square = np.zeros(
521 [len(leaf_hist), len(leaf_hist)], dtype=int
522 )
523
524
525 graFEI_truth_nFSP = len(leaf_hist)
526
527
528 for x, y in itertools.combinations(enumerate(leaf_hist), 2):
529 intersection = [
530 i for i in leaf_hist[x[1]] if i in leaf_hist[y[1]]
531 ]
532 true_LCA_square[x[0], y[0]] = levels[intersection[-1]]
533 true_LCA_square[y[0], x[0]] = levels[intersection[-1]]
534
535 x_leaves = p_indices
536 y_leaves = list(leaf_hist.keys())
537
538
539
540 locs = np.array(
541 [
542 np.where(y_leaves == i)[0].item()
543 if (i in y_leaves)
544 else 0
545 for i in x_leaves
546 ],
547 dtype=int,
548 )
549
550
551 true_LCA_square = true_LCA_square[locs, :][:, locs]
552
553
554
555 x_rows = np.array(
556 [
557 vm.evaluate("ancestorBIndex", p) == array_index
558 for p in p_list
559 ]
560 ) if self.mc_particle != "Upsilon(4S):MC" else evt_primary
561
562 primaries_from_right_cand = np.logical_and(evt_primary, x_rows)
563
564
565 true_LCA_square = np.where(
566 primaries_from_right_cand, true_LCA_square, 0
567 )
568
569 true_LCA_square = np.where(
570 primaries_from_right_cand[:, None], true_LCA_square, 0
571 )
572
573
574 true_LCA_square = torch.tensor(true_LCA_square, dtype=int)
575 b2.B2DEBUG(10, "True LCA:\n", true_LCA_square)
576
577
578 if (true_LCA_square == predicted_LCA_square).all():
579 graFEI_truth_perfectLCA = 1
580 b2.B2DEBUG(10, "LCA perfectly reconstructed!")
581
582
583 graFEI_truth_isSemileptonic = int(semilep_flag)
584
585
586 graFEI_truth_perfectEvent = int(graFEI_truth_perfectLCA and graFEI_truth_perfectMasses)
587
588
589 candidate.addExtraInfo("graFEI_truth_perfectLCA", graFEI_truth_perfectLCA)
590 candidate.addExtraInfo("graFEI_truth_perfectMasses", graFEI_truth_perfectMasses)
591 candidate.addExtraInfo("graFEI_truth_perfectEvent", graFEI_truth_perfectEvent)
592 candidate.addExtraInfo("graFEI_truth_isSemileptonic", graFEI_truth_isSemileptonic)
593 candidate.addExtraInfo("graFEI_truth_nFSP", graFEI_truth_nFSP)
594 candidate.addExtraInfo("graFEI_truth_nPhotons", graFEI_truth_nPhotons)
595 candidate.addExtraInfo("graFEI_truth_nElectrons", graFEI_truth_nElectrons)
596 candidate.addExtraInfo("graFEI_truth_nMuons", graFEI_truth_nMuons)
597 candidate.addExtraInfo("graFEI_truth_nPions", graFEI_truth_nPions)
598 candidate.addExtraInfo("graFEI_truth_nKaons", graFEI_truth_nKaons)
599 candidate.addExtraInfo("graFEI_truth_nProtons", graFEI_truth_nProtons)
600 candidate.addExtraInfo("graFEI_truth_nOthers", graFEI_truth_nOthers)
a (simplified) python wrapper for StoreObjPtr.