2 Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
4 #include "xAODTruth/TruthVertex.h"
7 template <bool IsDebug>
8 const xAOD::TruthParticle *ElasticDecayUtil<IsDebug>::getMother(const xAOD::TruthParticle &truth_particle, float max_energy_loss) const {
9 // from ElasticTruthTrajectoryBuilder.cxx
10 // Restrict to quasi-elastic processes (e.g. brems, delta-rays, pi->pi+Delta).
12 // Require not more than 2 outgoing particles. Note that
13 // delta-rays for primary==electron is a special case, because we have two
14 // outgoing particles with the same PDG id. The "correct" one
15 // is that with the higher energy (NOT pt).
17 // allow 1 outgoing to cover possible vertexes from interaction in detector material
18 // disallow hadronic pi->N*pi etc.
20 const xAOD::TruthParticle *last_parent = &truth_particle;
22 typename std::conditional<IsDebug, unsigned int,Empty>::type stat_n_parents {};
24 // follow decay chain upwards until PdgId changes, a particle has more than one
25 // parent, or a particle has more than one sibling.
27 const xAOD::TruthVertex* prod_vtx = last_parent->production_vertex();
29 if (prod_vtx->nIncomingParticles() != 1) break;
30 if (prod_vtx->nOutgoingParticles() > 2) break;
31 const xAOD::TruthParticle *grand_parent = prod_vtx->incomingParticle(0u);
32 if (!grand_parent) break;
33 if (grand_parent->pdgId() != last_parent->pdgId()) break;
34 if constexpr(IsDebug) {
35 unsigned int particle_type_i = kNoSibling;
36 std::lock_guard<std::mutex> lock(m_mutex);
37 if (prod_vtx->nOutgoingParticles()>1) {
38 int sibling_particle_type =0;
39 for (unsigned int part_out_i=0u; part_out_i<prod_vtx->nOutgoingParticles(); ++part_out_i) {
40 const xAOD::TruthParticle *a_sibling = prod_vtx->outgoingParticle(part_out_i);
42 if (a_sibling != last_parent && a_sibling) {
43 if (particle_type_i == kNoSibling) {
44 sibling_particle_type = a_sibling->absPdgId();
45 if (sibling_particle_type==11) {
46 particle_type_i = kElectron;
48 else if (sibling_particle_type==22) {
49 particle_type_i = kPhoton;
52 particle_type_i = kOther;
55 else if (sibling_particle_type != a_sibling->absPdgId() ) {
56 particle_type_i = kOther;
58 m_outgoingPdgIds.insert( a_sibling->pdgId() );
63 m_energyLossStat.at(particle_type_i).add( grand_parent->e() - last_parent->e() );
66 if (prod_vtx->nOutgoingParticles()>1) {
67 unsigned int part_out_i=0u;
68 for (; part_out_i<prod_vtx->nOutgoingParticles(); ++part_out_i) {
69 const xAOD::TruthParticle *a_sibling = prod_vtx->outgoingParticle(part_out_i);
70 if (a_sibling != last_parent) {
71 if (last_parent->pdgId() == a_sibling->pdgId() && (last_parent->pdgId() != 11 || last_parent->e() < a_sibling->e())) {
76 if (part_out_i < prod_vtx->nOutgoingParticles()) break;
78 if (last_parent->e() - grand_parent->e() > max_energy_loss) break;
80 if constexpr(IsDebug) {
81 m_outgoingMotherPdgIds.insert( grand_parent->pdgId() );
85 last_parent = grand_parent;
87 if constexpr(IsDebug) {
88 std::lock_guard<std::mutex> lock(m_mutex);
89 m_nParentsStat.add(stat_n_parents);
94 template <bool IsDebug>
95 inline void ElasticDecayUtil<IsDebug>::setEnergyLossBinning(const typename std::conditional<IsDebug,
97 EmptyProperty>::type &binning) {
98 if constexpr(IsDebug) {
99 if (binning.size()>=3) {
100 for (ActsUtils::StatHist &stat : m_energyLossStat ) {
101 stat.setBinning( static_cast<unsigned int>(binning[0]), binning[1],binning[2]);
107 template <bool IsDebug>
108 template <class T_OutStream>
109 inline void ElasticDecayUtil<IsDebug>::dumpStatistics(T_OutStream &out) const {
110 if constexpr(IsDebug) {
111 std::stringstream tmp_out;
112 tmp_out << m_nParentsStat << " number of parents" << std::endl
113 << m_nParentsStat.histogramToString() << std::endl;
114 std::array<std::string,kNParticleTypes> energyLossLabel{"no-sibling",
115 "electron as sibling",
119 unsigned int type_i=0;
120 for (const auto &stat : m_energyLossStat) {
122 tmp_out << stat << " energy loss, decay with " << energyLossLabel.at(type_i) << std::endl
123 << stat.histogramToString() << std::endl;
127 tmp_out << "Sibling PDG IDs : ";
128 for(const auto &elm : m_outgoingPdgIds) {
129 tmp_out << " " << elm;
131 tmp_out << std::endl;
132 tmp_out << "PDG IDs of followed particles : ";
133 for(const auto &elm : m_outgoingMotherPdgIds) {
134 tmp_out << " " << elm;
136 tmp_out << std::endl;
137 out << tmp_out.str();