ATLAS Offline Software
TruthParticleRetriever.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "AnalysisJiveXML/TruthParticleRetriever.h"
6 
7 #include <string>
8 
9 #include "CLHEP/Units/SystemOfUnits.h"
10 
11 #include "McParticleEvent/TruthParticleContainer.h"
12 
13 namespace JiveXML {
14 
15  //--------------------------------------------------------------------------
16 
17  TruthParticleRetriever::TruthParticleRetriever(const std::string& type, const std::string& name, const IInterface* parent):
18  AthAlgTool(type, name, parent),
19  m_typeName("CompositeParticle")
20  {
21 
22  declareInterface<IDataRetriever>(this);
23 
24  declareProperty("StoreGateKey", m_sgKey = "SpclMC","StoreGate key for TruthParticles");
25  declareProperty("TruthStatus", m_truthStatus = 2,"Truth status cut, currently not applied"); // currently unused
26  declareProperty("TruthMaximumPdgId", m_truthMaximumPdgId = 40,
27  "Maximum Truth PDG-ID, default:< 40");
28  declareProperty("TruthPtCut", m_truthPtCut = 10.,"pT cut on Truth, default 10 GeV");
29  declareProperty("SkimTruth", m_skimTruth = true,"Some criteria to skim most important Truth info ");
30  }
31 
32  //--------------------------------------------------------------------------
33 
34  StatusCode TruthParticleRetriever::retrieve(ToolHandle<IFormatTool> &FormatTool) {
35 
36  const TruthParticleContainer *truthCont = NULL;
37 
38  if ( !evtStore()->contains<TruthParticleContainer>(m_sgKey)){
39  if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "No TruthParticle found in SG at "
40  << m_sgKey << endmsg;
41  return StatusCode::SUCCESS;
42  }
43  if ( evtStore()->retrieve(truthCont,m_sgKey).isFailure() ) {
44  if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "TruthParticle retrieval from SG failed "
45  << m_sgKey << endmsg;
46  return StatusCode::SUCCESS;
47  }
48  int nTruth = truthCont->size();
49 
50  if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << " Retrieving TruthParticles with size " << nTruth << endmsg;
51 
52 //code from: PhysicsAnalysis/HiggsPhys/HiggsToFourLeptons/src/H4lElectronPerformanceTool.cxx
53 
54  DataVect pt; pt.reserve(nTruth);
55  DataVect phi; phi.reserve(nTruth);
56  DataVect eta; eta.reserve(nTruth);
57  DataVect typeEV; typeEV.reserve(nTruth);
58  DataVect label; label.reserve(nTruth);
59  DataVect typeLabelStr; typeLabelStr.reserve(nTruth);
60  DataVect pdgId; pdgId.reserve(nTruth);
61  DataVect dataType; dataType.reserve(nTruth);
62 
63  TruthParticleContainer::const_iterator mcpartItr = truthCont->begin();
64  TruthParticleContainer::const_iterator mcpartItrE = truthCont->end();
65 
66  std::string truthLabels;
67  std::string statusList="_";
68  std::string typeLabel="n_a";
69  int pdgId2;
70  int motherPdgId;
71  int childPdgId;
72  int countTruth = 1;
73  bool samePdgIdFlag = false;
74  bool initialProcessFlag = false;
75  bool protectedParticleFlag = false;
76 
77  for(; (mcpartItr != mcpartItrE) ; ++mcpartItr) { // loop over TruthParticle
78  samePdgIdFlag = false;
79  initialProcessFlag = false;
80  protectedParticleFlag = false;
81  statusList += std::to_string((*mcpartItr)->status()) + "_";
82 
83  if ( (*mcpartItr)->et()/CLHEP::GeV < m_truthPtCut ){ continue; }
84  if ( abs( (*mcpartItr)->pdgId()) > m_truthMaximumPdgId ){ continue; }
85 
86  pdgId2 = (*mcpartItr)->pdgId();
87 
88  // important particles are protected (can come from gamma etc)
89  if (( abs(pdgId2) == 24 ) || ( abs(pdgId2) == 5 ) ||
90  ( abs(pdgId2) == 6 ) || ( abs(pdgId2) == 23 ) ||
91  ( abs(pdgId2) == 36 ) || ( abs(pdgId2) == 37 ) ||
92  ( abs(pdgId2) == 25 ) ){
93  protectedParticleFlag = true;
94  }
95 
96  // particle has to originate from b,t,W,Z,H,A,H+- only
97  for ( unsigned int iMother = 0; iMother < (*mcpartItr)->nParents(); ++iMother ) {
98  bool motherHasPdgId = (*mcpartItr)->mother(iMother)->hasPdgId();
99  if ( motherHasPdgId ){
100  motherPdgId = (*mcpartItr)->mother(iMother)->pdgId();
101  if ( motherPdgId == pdgId2 ){ samePdgIdFlag = true; }
102  if (( abs(motherPdgId) == 24 ) || ( abs(motherPdgId) == 5 ) ||
103  ( abs(motherPdgId) == 6 ) || ( abs(motherPdgId) == 23 ) ||
104  ( abs(motherPdgId) == 36 ) || ( abs(motherPdgId) == 37 ) ||
105  ( abs(pdgId2) == 25 ) ){
106  initialProcessFlag = true;
107  }
108  }
109  }
110  // jump to next particle (ignore) if PdgId is same as mother
111  if (m_skimTruth){ // some skimming criteria for simpler display of Truth
112  if ( !initialProcessFlag && !protectedParticleFlag ){
113  // if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << " Particle rejected: not initial process " << endmsg;
114  continue;
115  }
116  }
117  if ( samePdgIdFlag ){
118  // if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << " Particle rejected: Same PdgId " << endmsg;
119  continue;
120  }
121 
122  pt.push_back( DataType((*mcpartItr)->et()/CLHEP::GeV ) );
123  phi.push_back( DataType((*mcpartItr)->phi() ) );
124  eta.push_back( DataType((*mcpartItr)->eta() ) );
125  dataType.push_back( DataType( (*mcpartItr)->dataType() ));
126 
127 // if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << " TruthParticle No:" << countTruth << " pt:"
128 // << (*mcpartItr)->et()/CLHEP::GeV << endmsg;
129 
130  pdgId.push_back( DataType( pdgId2 ));
131 
132  typeLabel = "n_a";
133  if( abs(pdgId2) == 11) typeLabel = "Electron";
134  if( abs(pdgId2) == 12) typeLabel = "NeutrinoElectron";
135  if( abs(pdgId2) == 13) typeLabel = "Muon";
136  if( abs(pdgId2) == 14) typeLabel = "NeutrinoMuon";
137  if( abs(pdgId2) == 15) typeLabel = "Tau";
138  if( abs(pdgId2) == 16) typeLabel = "NeutrinoTau";
139  if( pdgId2 == 6) typeLabel = "Top";
140  if( pdgId2 == -6) typeLabel = "AntiTop";
141  if( pdgId2 == 5) typeLabel = "Bottom";
142  if( pdgId2 == -5) typeLabel = "AntiBottom";
143  if( pdgId2 == 22) typeLabel = "Photon";
144  if( pdgId2 == 23) typeLabel = "Z0";
145  if( pdgId2 == 224) typeLabel = "Wplus";
146  if( pdgId2 == -24) typeLabel = "Wminus";
147  if( pdgId2 == 36) typeLabel = "A0";
148  if( pdgId2 == 25) typeLabel = "Higgs0";
149  if(( abs(pdgId2) >= 1) && ( abs(pdgId2) <= 4)) typeLabel = "LightQuark";
150 
151  typeEV.push_back( DataType( typeLabel ) );
152 
153  truthLabels = "No" + DataType( countTruth ).toString() +"_Pdg="
154  + DataType( (*mcpartItr)->pdgId() ).toString()
155  + "_stat=" + DataType( (*mcpartItr)->status() ).toString()
156  + "_nParents=" + DataType( (*mcpartItr)->nParents() ).toString()
157  + "_nDecay=" + DataType( (*mcpartItr)->nDecay() ).toString();
158 
159  for ( unsigned int iChild = 0; iChild < (*mcpartItr)->nDecay(); ++iChild ) {
160  bool childHasPdgId = (*mcpartItr)->child(iChild)->hasPdgId();
161  if ( childHasPdgId ){
162  childPdgId = (*mcpartItr)->child(iChild)->pdgId();
163  truthLabels += "_CNo=" + DataType(iChild+1).toString() + "_CPdgId=" +
164  DataType( childPdgId ).toString() + "_CET=" +
165  DataType( (*mcpartItr)->child(iChild)->et()/CLHEP::GeV ).toString();
166  if ( (*mcpartItr)->child(iChild)->et()/CLHEP::GeV <= m_truthPtCut ){ continue; }
167  if (childPdgId == -11){
168  truthLabels += "_IntoElectron";
169  typeLabel += "_IntoElectron";
170  }
171  if (childPdgId == 11){
172  truthLabels += "_FromPositron";
173  typeLabel += "_FromPositron";
174  }
175  if (childPdgId == 13){
176  truthLabels += "_IntoMuonMinus";
177  typeLabel += "_IntoMuonMinus";
178  }
179  if (childPdgId == -13){
180  truthLabels += "_IntoMuonPlus";
181  typeLabel += "_IntoMuonPlus";
182  }
183  if (childPdgId == 15){
184  truthLabels += "_IntoTauMinus";
185  typeLabel += "_IntoTauMinus";
186  }
187  if (childPdgId == -15){
188  truthLabels += "_IntoTauPlus";
189  typeLabel += "_IntoTauPlus";
190  }
191  if (childPdgId == 5){
192  truthLabels += "_IntoBottom";
193  typeLabel += "_IntoBottom";
194  }
195  if (childPdgId == -5){
196  truthLabels += "_IntoAntiBottom";
197  typeLabel += "_IntoAntiBottom";
198  }
199  if (childPdgId == 6){
200  truthLabels += "_IntoTop";
201  typeLabel += "_IntoTop";
202  }
203  if (childPdgId == -6){
204  truthLabels += "_IntoAntiTop";
205  typeLabel += "_IntoAntiTop";
206  }
207  if ( ( abs(childPdgId) >= 1) && ( abs(childPdgId) <=4) ) {
208  truthLabels += "_IntoLightQuark";
209  typeLabel += "_IntoLightQuark";
210  }
211  if (childPdgId == -24){
212  truthLabels += "_IntoWminus";
213  typeLabel += "_IntoWminus";
214  }
215  if (childPdgId == 24){
216  truthLabels += "_IntoWplus";
217  typeLabel += "_IntoWplus";
218  }
219  if (childPdgId == 25){
220  truthLabels += "_IntoHiggs0";
221  typeLabel += "_IntoHiggs0";
222  }
223  if (childPdgId == 23){
224  truthLabels += "_IntoZ0";
225  typeLabel += "_IntoZ0";
226  }
227  }
228  }
229  for ( unsigned int iMother = 0; iMother < (*mcpartItr)->nParents(); ++iMother ) {
230  bool motherHasPdgId = (*mcpartItr)->mother(iMother)->hasPdgId();
231  if ( motherHasPdgId ){
232  int motherPdgId = (*mcpartItr)->mother(iMother)->pdgId();
233  truthLabels += "_MNo=" + DataType(iMother+1).toString() + "_MPdgId=" +
234  DataType( motherPdgId ).toString();
235  if ( (*mcpartItr)->mother(iMother)->et()/CLHEP::GeV <= m_truthPtCut ){ continue; }
236  if (motherPdgId == -24){
237  truthLabels += "_FromWminus";
238  typeLabel += "_FromWminus";
239  }
240  if (motherPdgId == 24){
241  truthLabels += "_FromWplus";
242  typeLabel += "_FromWplus";
243  }
244  if (motherPdgId == 25){
245  truthLabels += "_FromHiggs0";
246  typeLabel += "_FromHiggs0";
247  }
248  if (motherPdgId == 23){
249  truthLabels += "_FromZ0";
250  typeLabel += "_FromZ0";
251  }
252  if (motherPdgId == 6){
253  truthLabels += "_FromTop";
254  typeLabel += "_FromTop";
255  }
256  if (motherPdgId == -6){
257  truthLabels += "_FromAntiTop";
258  typeLabel += "_FromAntiTop";
259  }
260  if (motherPdgId == 5){
261  truthLabels += "_FromBottom";
262  typeLabel += "_FromBottom";
263  }
264  if (motherPdgId == -5){
265  truthLabels += "_FromAntiBottom";
266  typeLabel += "_FromAntiBottom";
267  }
268  if (motherPdgId == 36){
269  truthLabels += "_FromA0";
270  typeLabel += "_FromA0";
271  }
272  if (motherPdgId == 37){
273  truthLabels += "_FromHiggsPlus";
274  typeLabel += "_FromHiggsPlus";
275  }
276  if (motherPdgId == -37){
277  truthLabels += "_FromHiggsMinus";
278  typeLabel += "_FromHiggsMinus";
279  }
280  } // end of hasMotherPdgId
281  } // end of MotherParticle loop
282  countTruth++;
283  if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "TruthParticles label: " << DataType( truthLabels ).toString() << endmsg;
284  label.push_back( DataType( truthLabels ).toString() );
285  typeLabelStr.push_back( DataType( typeLabel ).toString() );
286 
287  // log << MSG::DEBUG << "TruthParticles status: " << statusList << endmsg;
288  } // end of TruthParticle Loop
289  DataMap myDataMap;
290  myDataMap["pt"] = pt;
291  myDataMap["phi"] = phi;
292  myDataMap["eta"] = eta;
293  myDataMap["typeEV"] = typeEV;
294 // myDataMap["label"] = label;
295  myDataMap["label"] = typeLabelStr;
296  myDataMap["pdgId"] = pdgId;
297  myDataMap["dataType"] = dataType;
298 
299  if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << dataTypeName() << ": "<< phi.size() << endmsg;
300 
301  //forward data to formating tool
302  return FormatTool->AddToEvent(dataTypeName(), m_sgKey, &myDataMap);
303  }
304 }
JiveXML::TruthParticleRetriever::m_truthMaximumPdgId
int m_truthMaximumPdgId
Definition: TruthParticleRetriever.h:38
DataModel_detail::const_iterator
Const iterator class for DataVector/DataList.
Definition: DVLIterator.h:82
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:64
TruthParticleContainer.h
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:79
AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
PlotCalibFromCool.label
label
Definition: PlotCalibFromCool.py:78
JiveXML::DataVect
std::vector< DataType > DataVect
Defines a map with a key and a vector of DataType objects e.g.
Definition: DataType.h:58
DataType
OFFLINE_FRAGMENTS_NAMESPACE::PointerType DataType
Definition: RoIBResultByteStreamTool.cxx:25
test_pyathena.pt
pt
Definition: test_pyathena.py:11
JiveXML::TruthParticleRetriever::m_sgKey
std::string m_sgKey
properties:
Definition: TruthParticleRetriever.h:36
downloadSingle.dataType
string dataType
Definition: downloadSingle.py:18
JiveXML::TruthParticleRetriever::retrieve
virtual StatusCode retrieve(ToolHandle< IFormatTool > &FormatTool)
Retrieve all the data.
Definition: TruthParticleRetriever.cxx:34
AthCommonMsg< AlgTool >::msgLvl
bool msgLvl(const MSG::Level lvl) const
Definition: AthCommonMsg.h:30
plotBeamSpotCompare.statusList
statusList
Definition: plotBeamSpotCompare.py:200
JiveXML::DataMap
std::map< std::string, DataVect > DataMap
Definition: DataType.h:59
TruthParticleContainer
Definition: PhysicsAnalysis/TruthParticleID/McParticleEvent/McParticleEvent/TruthParticleContainer.h:42
PowhegPy8EG_H2a.pdgId
dictionary pdgId
Definition: PowhegPy8EG_H2a.py:128
AthCommonDataStore< AthCommonMsg< AlgTool > >::evtStore
ServiceHandle< StoreGateSvc > & evtStore()
The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.
Definition: AthCommonDataStore.h:85
Amg::toString
std::string toString(const Translation3D &translation, int precision=4)
GeoPrimitvesToStringConverter.
Definition: GeoPrimitivesToStringConverter.h:40
endmsg
#define endmsg
Definition: AnalysisConfig_Ntuple.cxx:63
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
test_pyathena.parent
parent
Definition: test_pyathena.py:15
JiveXML::TruthParticleRetriever::m_truthStatus
int m_truthStatus
Definition: TruthParticleRetriever.h:37
JiveXML::TruthParticleRetriever::dataTypeName
virtual std::string dataTypeName() const
Return the name of the data type.
Definition: TruthParticleRetriever.h:28
JiveXML
This header is shared inbetween the C-style server thread and the C++ Athena ServerSvc.
Definition: BadLArRetriever.cxx:21
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:192
ActsTrk::to_string
std::string to_string(const DetectorType &type)
Definition: GeometryDefs.h:34
JiveXML::TruthParticleRetriever::m_truthPtCut
float m_truthPtCut
Definition: TruthParticleRetriever.h:39
DataVector::end
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
DEBUG
#define DEBUG
Definition: page_access.h:11
AthCommonMsg< AlgTool >::msg
MsgStream & msg() const
Definition: AthCommonMsg.h:24
JiveXML::TruthParticleRetriever::m_skimTruth
bool m_skimTruth
Definition: TruthParticleRetriever.h:40
JiveXML::TruthParticleRetriever::TruthParticleRetriever
TruthParticleRetriever(const std::string &type, const std::string &name, const IInterface *parent)
Standard Constructor.
Definition: TruthParticleRetriever.cxx:17
AthAlgTool
Definition: AthAlgTool.h:26
GeV
#define GeV
Definition: CaloTransverseBalanceVecMon.cxx:30
DataVector::size
size_type size() const noexcept
Returns the number of elements in the collection.
DataVector::begin
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
TruthParticleRetriever.h
Generated on Sun Jun 30 2024 21:32:00 for ATLAS Offline Software by   doxygen 1.8.18