ATLAS Offline Software
MCTruthClassifierAthena.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #if !defined(XAOD_ANALYSIS) && !defined(GENERATIONBASE) // Can only be used in Athena
6 #include "MCTruthClassifier/MCTruthClassifier.h"
7 // xAOD EDM includes
8 #include "xAODCaloEvent/CaloCluster.h"
9 #include "xAODEgamma/EgammaxAODHelpers.h"
10 #include "xAODTracking/TrackParticle.h"
11 #include "xAODTruth/TruthParticle.h"
12 #include "xAODTruth/TruthParticleContainer.h"
13 // Athena only includes
14 #include "AthenaKernel/Units.h"
15 #include "AtlasHepMC/GenParticle.h"
16 #include "GeneratorObjects/xAODTruthParticleLink.h"
17 #include "RecoToolInterfaces/IParticleCaloExtensionTool.h"
18 #include "TrkEventPrimitives/PropDirection.h"
19 #include "TrkParametersIdentificationHelpers/TrackParametersIdHelper.h"
20 
21 #include <cmath>
22 
23 using namespace MCTruthPartClassifier;
24 
25 namespace {
26 std::unique_ptr<Trk::CurvilinearParameters> extractParamFromTruth(const xAOD::TruthParticle& particle) {
27  // get start parameters
28  const xAOD::TruthVertex* pvtx = particle.prodVtx();
29  if (pvtx == nullptr) return nullptr;
30  double charge = particle.charge();
31  Amg::Vector3D pos(pvtx->x(), pvtx->y(), pvtx->z());
32  Amg::Vector3D mom(particle.px(), particle.py(), particle.pz());
33  // Aproximate neutral particles as charged with infinite momentum
34  if (particle.isNeutral()) {
35  charge = 1.;
36  mom.normalize();
37  mom *= 1e10;
38  }
39  return std::make_unique<Trk::CurvilinearParameters>(pos, mom, charge);
40 }
41 }
42 
43 // Methods using directly the extrapolator usable only from Athena
44 const xAOD::TruthParticle* MCTruthClassifier::egammaClusMatch(const xAOD::CaloCluster* clus, bool isFwrdEle, MCTruthPartClassifier::Info* info) const
45 {
46  ATH_MSG_DEBUG("Executing egammaClusMatch ");
47  const xAOD::TruthParticle* theMatchPart = nullptr;
48  const EventContext& ctx = info ? info->eventContext : Gaudi::Hive::currentContext();
49  // retrieve collection and get a pointer
50  SG::ReadHandle<xAOD::TruthParticleContainer> truthParticleContainerReadHandle(m_truthParticleContainerKey, ctx);
51 
52  if (!truthParticleContainerReadHandle.isValid()) {
53  ATH_MSG_WARNING( " Invalid ReadHandle for xAOD::TruthParticleContainer with key: " << truthParticleContainerReadHandle.key());
54  return theMatchPart;
55  }
56 
57  SG::ReadCondHandle<CaloDetDescrManager> caloMgrHandle{ m_caloMgrKey, ctx };
58  if (!caloMgrHandle.isValid()) {
59  ATH_MSG_WARNING(" Invalid ReadCondHandle for CaloDetDescrManager with key: " << m_caloMgrKey.key());
60  return theMatchPart;
61  }
62 
63  const CaloDetDescrManager* caloDDMgr = *caloMgrHandle;
64  const xAOD::TruthParticle* theEgamma(nullptr);
65  const xAOD::TruthParticle* theLeadingPartInCone(nullptr);
66  const xAOD::TruthParticle* theBestPartOutCone(nullptr);
67  const xAOD::TruthParticle* theBestPartdR(nullptr);
68  double LeadingPhtPT(0);
69  double LeadingPartPT(0);
70  double LeadingPhtdR(999.);
71  double LeadingPartdR(999.);
72  double BestPartdR(999.);
73  double etaClus = clus->etaBE(2);
74  double phiClus = clus->phiBE(2);
75  if (etaClus < -900) {
76  etaClus = clus->eta();
77  }
78  if (phiClus < -900) {
79  phiClus = clus->phi();
80  }
81  std::vector<const xAOD::TruthParticle*> tps;
82  if (!m_truthInConeTool->particlesInCone(ctx, etaClus, phiClus, 0.5, tps)) {
83  ATH_MSG_WARNING("Truth Particle in Cone failed");
84  return theMatchPart;
85  }
86 
87  for (const auto* const thePart : tps) {
88  // loop over the stable particle
89  if (!MC::isStable(thePart)) continue;
90  // excluding G4 particle
91  if ((!isFwrdEle || (isFwrdEle && m_FwdElectronUseG4Sel)) && HepMC::is_simulation_particle(thePart)) continue;
92  long iParticlePDG = thePart->pdgId();
93  // excluding neutrino
94  if (std::abs(iParticlePDG) == 12 || std::abs(iParticlePDG) == 14 || std::abs(iParticlePDG) == 16) continue;
95  double pt = thePart->pt() / Athena::Units::GeV;
96  double q = thePart?thePart->charge():0.0;
97  // exclude charged particles with pT<1 GeV
98  if (q != 0 && pt < m_pTChargePartCut) continue;
99  if (q == 0 && pt < m_pTNeutralPartCut) continue;
100 
101  // eleptical cone for extrapolations m_partExtrConePhi X m_partExtrConeEta
102  if (!isFwrdEle && m_ROICone && std::hypot( detPhi(phiClus, thePart->phi())/m_partExtrConePhi, detEta(etaClus, thePart->eta())/m_partExtrConeEta) > 1.0) {
103  continue;
104  }
105  // Also check if the clus and true have different sign , i they need both to be <0 or >0
106  if (isFwrdEle && // It is forward and
107  (((etaClus < 0) - (thePart->eta() < 0) != 0)
108  // The truth eta has different sign wrt to the fwd electron
109  || (std::fabs(thePart->eta()) < m_FwdElectronTruthExtrEtaCut) // or the truth is less than 2.4 (default cut)
110  || (std::fabs(thePart->eta() - etaClus) > m_FwdElectronTruthExtrEtaWindowCut) // or if the delta Eta between el and truth is > 0.15
111  ) // then do no extrapolate this truth Particle for this fwd electron
112  ) {
113  continue;
114  }
115  double dR(-999.);
116  bool isNCone = false;
117  bool isExt = genPartToCalo(ctx, clus, thePart, isFwrdEle, dR, isNCone, *caloDDMgr);
118  if (!isExt) continue;
119  theMatchPart = MC::find_matching(truthParticleContainerReadHandle.ptr(), thePart);
120  if (info) {
121  info->egPartPtr.push_back(thePart);
122  info->egPartdR.push_back(dR);
123  info->egPartClas.push_back(particleTruthClassifier(theMatchPart, info));
124  }
125  // Gen particles Not forward
126  if (!isFwrdEle) {
127  // the leading photon or electron inside narrow eleptical cone
128  // m_phtClasConePhi X m_phtClasConeEta
129  if ((iParticlePDG == 22 || std::abs(iParticlePDG) == 11) && isNCone && pt > LeadingPhtPT) {
130  theEgamma = thePart;
131  LeadingPhtPT = pt;
132  LeadingPhtdR = dR;
133  }
134  // leading particle (excluding photon and electron) inside narrow eleptic
135  // cone m_phtClasConePhi X m_phtClasConeEta
136  if ((iParticlePDG != 22 && std::abs(iParticlePDG) != 11) && isNCone && pt > LeadingPartPT) {
137  theLeadingPartInCone = thePart;
138  LeadingPartPT = pt;
139  LeadingPartdR = dR;
140  };
141  // the best dR matched particle outside narrow eleptic cone cone
142  // m_phtClasConePhi X m_phtClasConeEta
143  if (!isNCone && dR < BestPartdR) {
144  theBestPartOutCone = thePart;
145  BestPartdR = dR;
146  };
147  } else {
148  if (dR < BestPartdR) {
149  theBestPartdR = thePart;
150  BestPartdR = dR;
151  };
152  }
153  }
154 
155  if (theEgamma != nullptr) {
156  theMatchPart = MC::find_matching(truthParticleContainerReadHandle.ptr(), theEgamma);
157  if (info) info->deltaRMatch = LeadingPhtdR;
158  } else if (theLeadingPartInCone != nullptr) {
159  theMatchPart = MC::find_matching(truthParticleContainerReadHandle.ptr(),theLeadingPartInCone);
160  if (info) info->deltaRMatch = LeadingPartdR;
161  } else if (theBestPartOutCone != nullptr) {
162  theMatchPart = MC::find_matching(truthParticleContainerReadHandle.ptr(),theBestPartOutCone);
163  if (info) info->deltaRMatch = BestPartdR;
164  } else if (isFwrdEle && theBestPartdR != nullptr) {
165  theMatchPart = MC::find_matching(truthParticleContainerReadHandle.ptr(),theBestPartdR );
166  if (info) info->deltaRMatch = BestPartdR;
167  } else {
168  theMatchPart = nullptr;
169  }
170  if (isFwrdEle || theMatchPart != nullptr || !m_inclG4part) return theMatchPart;
171 
172  // additional loop over G4 particles,
173  for (const auto* const thePart : tps) {
174  if (!MC::isStable(thePart)) continue;
175  if (!HepMC::is_simulation_particle(thePart)) continue;
176  long iParticlePDG = thePart->pdgId();
177  // exclude neutrino
178  if (std::abs(iParticlePDG) == 12 || std::abs(iParticlePDG) == 14 || std::abs(iParticlePDG) == 16) continue;
179  if (thePart->decayVtx() != nullptr) continue;
180  if (std::hypot( detPhi(phiClus, thePart->phi())/m_partExtrConePhi, detEta(etaClus, thePart->eta())/m_partExtrConeEta ) > 1.0) continue;
181 
182  double pt = thePart->pt() / Athena::Units::GeV;
183  double q = thePart->charge();
184  // exclude charged particles with pT<1 GeV
185  if (q != 0 && pt < m_pTChargePartCut) continue;
186  if (q == 0 && pt < m_pTNeutralPartCut) continue;
187 
188  double dR(-999.);
189  bool isNCone = false;
190  bool isExt = genPartToCalo(ctx, clus, thePart, isFwrdEle, dR, isNCone, *caloDDMgr);
191  if (!isExt) continue;
192 
193  theMatchPart = MC::find_matching(truthParticleContainerReadHandle.ptr(),thePart);
194  if (info) {
195  info->egPartPtr.push_back(thePart);
196  info->egPartdR.push_back(dR);
197  info->egPartClas.push_back(particleTruthClassifier(theMatchPart, info));
198  }
199  // the leading photon or electron inside narrow eleptical cone
200  // m_phtClasConePhi X m_phtClasConeEta
201  if ((iParticlePDG == 22 || std::abs(iParticlePDG) == 11) && isNCone && pt > LeadingPhtPT) {
202  theEgamma = thePart;
203  LeadingPhtPT = pt;
204  LeadingPhtdR = dR;
205  }
206  // leading particle (excluding photon or electron) inside narrow eleptic
207  // cone m_phtClasConePhi X m_phtClasConeEta
208  if ((iParticlePDG != 22 && std::abs(iParticlePDG) != 11) && isNCone && pt > LeadingPartPT) {
209  theLeadingPartInCone = thePart;
210  LeadingPartPT = pt;
211  LeadingPartdR = dR;
212  }
213  // the best dR matched particle outside narrow eleptic cone cone
214  // m_phtClasConePhi X m_phtClasConeEta
215  if (!isNCone && dR < BestPartdR) {
216  theBestPartOutCone = thePart;
217  BestPartdR = dR;
218  }
219  }
220 
221  if (theEgamma != nullptr) {
222  theMatchPart = MC::find_matching(truthParticleContainerReadHandle.ptr(),theEgamma);
223  if (info) info->deltaRMatch = LeadingPhtdR;
224  } else if (theLeadingPartInCone != nullptr) {
225  theMatchPart = MC::find_matching(truthParticleContainerReadHandle.ptr(),theLeadingPartInCone);
226  if (info) info->deltaRMatch = LeadingPartdR;
227  } else if (theBestPartOutCone != nullptr) {
228  theMatchPart = MC::find_matching(truthParticleContainerReadHandle.ptr(),theBestPartOutCone);
229  if (info) info->deltaRMatch = BestPartdR;
230  } else {
231  theMatchPart = nullptr;
232  }
233  ATH_MSG_DEBUG("succeeded egammaClusMatch ");
234  return theMatchPart;
235 }
236 
237 bool MCTruthClassifier::genPartToCalo(const EventContext& ctx,
238  const xAOD::CaloCluster* clus,
239  const xAOD::TruthParticle* thePart,
240  bool isFwrdEle,
241  double& dRmatch,
242  bool& isNarrowCone,
243  const CaloDetDescrManager& caloDDMgr) const
244 {
245  dRmatch = -999.;
246  isNarrowCone = false;
247  if (thePart == nullptr) return false;
248  double phiClus = clus->phiBE(2);
249  double etaClus = clus->etaBE(2);
250  if (etaClus < -900) {
251  etaClus = clus->eta();
252  }
253  if (phiClus < -900) {
254  phiClus = clus->phi();
255  }
256  //--FixMe
257  if (isFwrdEle || (etaClus == 0. && phiClus == 0.)) {
258  phiClus = clus->phi();
259  etaClus = clus->eta();
260  }
261  // define calo sample
262  CaloSampling::CaloSample sample = CaloSampling::EMB2;
263  if ((clus->inBarrel() && !clus->inEndcap()) ||
264  (clus->inBarrel() && clus->inEndcap() && clus->eSample(CaloSampling::EMB2) >= clus->eSample(CaloSampling::EME2))) {
265  // Barrel
266  sample = CaloSampling::EMB2;
267  } else if ((!clus->inBarrel() && clus->inEndcap() && !isFwrdEle) ||
268  (clus->inBarrel() && clus->inEndcap() && clus->eSample(CaloSampling::EME2) > clus->eSample(CaloSampling::EMB2))) {
269  // End-cap
270  sample = CaloSampling::EME2;
271  } else if (isFwrdEle && clus->inEndcap()) {
272  // FCAL
273  sample = CaloSampling::FCAL2;
274  } else {
275  return false;
276  }
277  std::unique_ptr<Trk::CurvilinearParameters> params = extractParamFromTruth(*thePart);
278  if (!params) return false;
279  // create extension to sample
280  std::vector<CaloSampling::CaloSample> samples = { sample };
281  auto extension = m_caloExtensionTool->layersCaloExtension(ctx, *params, samples, etaClus, caloDDMgr);
282  bool extensionOK = (!extension.empty());
283  if (!extensionOK) {
284  ATH_MSG_WARNING("extrapolation of Truth Particle with eta " << thePart->eta() << " , charge " << thePart->charge() << " , Pt " << thePart->pt() << " to calo failed");
285  return false;
286  }
287  double etaCalo = extension[0].second->position().eta();
288  double phiCalo = extension[0].second->position().phi();
289 
290  double dPhi = detPhi(phiCalo, phiClus);
291  double dEta = detEta(etaCalo, etaClus);
292  dRmatch = std::hypot(dPhi, dEta);
293 
294  if ((!isFwrdEle && dRmatch > m_phtdRtoTrCut) || (isFwrdEle && dRmatch > m_fwrdEledRtoTrCut)) return false;
295  if (!isFwrdEle && std::hypot( dPhi/m_phtClasConePhi, dEta/m_phtClasConeEta ) <= 1.0) isNarrowCone = true;
296  return true;
297 }
298 #endif
grepfile.info
info
Definition: grepfile.py:38
xAOD::CaloCluster_v1::phi
virtual double phi() const
The azimuthal angle ( ) of the particle.
Definition: CaloCluster_v1.cxx:256
Trk::ParticleSwitcher::particle
constexpr ParticleHypothesis particle[PARTICLEHYPOTHESES]
the array of masses
Definition: ParticleHypothesis.h:76
SG::ReadCondHandle
Definition: ReadCondHandle.h:44
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:70
TruthParticleContainer.h
test_pyathena.pt
pt
Definition: test_pyathena.py:11
PropDirection.h
xAOD::CaloCluster_v1::phiBE
float phiBE(const unsigned layer) const
Get the phi in one layer of the EM Calo.
Definition: CaloCluster_v1.cxx:680
MCTruthClassifier.h
MC::find_matching
T find_matching(C TruthTES, T bcin)
Function to find a particle in container.
Definition: HepMCHelpers.h:103
GenParticle.h
xAOD::TruthVertex_v1::y
float y() const
Vertex y displacement.
xAOD::CaloCluster_v1::etaBE
float etaBE(const unsigned layer) const
Get the eta in one layer of the EM Calo.
Definition: CaloCluster_v1.cxx:644
xAOD::CaloCluster_v1
Description of a calorimeter cluster.
Definition: CaloCluster_v1.h:59
xAOD::CaloCluster_v1::inEndcap
bool inEndcap() const
Returns true if at least one clustered cell in the endcap.
Definition: CaloCluster_v1.h:901
EgammaxAODHelpers.h
ParticleGun_EoverP_Config.mom
mom
Definition: ParticleGun_EoverP_Config.py:63
TrackParametersIdHelper.h
constants.EMB2
int EMB2
Definition: Calorimeter/CaloClusterCorrection/python/constants.py:54
FullCPAlgorithmsTest_eljob.sample
sample
Definition: FullCPAlgorithmsTest_eljob.py:100
HepMC::is_simulation_particle
bool is_simulation_particle(const T &p)
Method to establish if a particle (or barcode) was created during the simulation (TODO update to be s...
Definition: MagicNumbers.h:299
xAOD::CaloCluster_v1::eta
virtual double eta() const
The pseudorapidity ( ) of the particle.
Definition: CaloCluster_v1.cxx:251
CaloSampling::CaloSample
CaloSample
Definition: Calorimeter/CaloGeoHelpers/CaloGeoHelpers/CaloSampling.h:22
CaloCluster.h
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
01SubmitToGrid.samples
samples
Definition: 01SubmitToGrid.py:58
xAOD::TruthParticle_v1
Class describing a truth particle in the MC record.
Definition: TruthParticle_v1.h:41
TauGNNUtils::Variables::Track::dPhi
bool dPhi(const xAOD::TauJet &tau, const xAOD::TauTrack &track, double &out)
Definition: TauGNNUtils.cxx:530
xAOD::CaloCluster_v1::inBarrel
bool inBarrel() const
Returns true if at least one clustered cell in the barrel.
Definition: CaloCluster_v1.h:896
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
MCTruthClassifier::genPartToCalo
bool genPartToCalo(const EventContext &ctx, const xAOD::CaloCluster *clus, const xAOD::TruthParticle *thePart, bool isFwrdEle, double &dRmatch, bool &isNarrowCone, const CaloDetDescrManager &caloDDMgr) const
Definition: MCTruthClassifierAthena.cxx:237
xAOD::TruthVertex_v1
Class describing a truth vertex in the MC record.
Definition: TruthVertex_v1.h:41
charge
double charge(const T &p)
Definition: AtlasPID.h:494
xAOD::TruthParticle_v1::eta
virtual double eta() const override final
The pseudorapidity ( ) of the particle.
Definition: TruthParticle_v1.cxx:174
Units.h
Wrapper to avoid constant divisions when using units.
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
MCTruthClassifier::egammaClusMatch
const xAOD::TruthParticle * egammaClusMatch(const xAOD::CaloCluster *, bool, MCTruthPartClassifier::Info *info) const
Definition: MCTruthClassifierAthena.cxx:44
SG::VarHandleBase::key
virtual const std::string & key() const override final
Return the StoreGate ID for the referenced object.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleBase.cxx:64
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:18
SG::ReadHandle::ptr
const_pointer_type ptr()
Dereference the pointer.
xAOD::TruthVertex_v1::x
float x() const
Vertex x displacement.
TrackParticle.h
IParticleCaloExtensionTool.h
MC::isStable
bool isStable(const T &p)
Definition: HepMCHelpers.h:30
xAOD::CaloCluster_v1::eSample
float eSample(const CaloSample sampling) const
Definition: CaloCluster_v1.cxx:521
CaloDetDescrManager
This class provides the client interface for accessing the detector description information common to...
Definition: CaloDetDescrManager.h:473
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
MCTruthPartClassifier
Definition: TruthClassifiers.h:12
xAOD::TruthVertex_v1::z
float z() const
Vertex longitudinal distance along the beam line form the origin.
extractSporadic.q
list q
Definition: extractSporadic.py:98
xAOD::TruthParticle_v1::pt
virtual double pt() const override final
The transverse momentum ( ) of the particle.
Definition: TruthParticle_v1.cxx:166
CaloCell_ID_FCS::FCAL2
@ FCAL2
Definition: FastCaloSim_CaloCell_ID.h:42
PowhegControl_ttFCNC_NLO.params
params
Definition: PowhegControl_ttFCNC_NLO.py:226
TruthParticle.h
xAOD::TruthParticle_v1::pdgId
int pdgId() const
PDG ID code.
TauGNNUtils::Variables::Track::dEta
bool dEta(const xAOD::TauJet &tau, const xAOD::TauTrack &track, double &out)
Definition: TauGNNUtils.cxx:525
GeV
#define GeV
Definition: CaloTransverseBalanceVecMon.cxx:30
MCTruthPartClassifier::Info
Definition: IMCTruthClassifier.h:49
constants.EME2
int EME2
Definition: Calorimeter/CaloClusterCorrection/python/constants.py:56
xAOD::TruthParticle_v1::charge
double charge() const
Physical charge.
Generated on Sun Jun 30 2024 21:21:18 for ATLAS Offline Software by   doxygen 1.8.18