ATLAS Offline Software
TruthNtupleTool.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // TruthNtupleTool.cxx
7 // Source file for class TruthNtupleTool
9 // (c) ATLAS Detector software
11 // Sebastian.Fleischmann@cern.ch
13 
14 
15 //Gaudi
16 #include "GaudiKernel/ITHistSvc.h"
17 
18 // Trk
19 #include "TrkValTools/TruthNtupleTool.h"
20 #include "TrkParameters/TrackParameters.h"
21 // Truth
22 #include "TrkTruthData/TrackTruth.h"
23 
24 #include "AtlasHepMC/GenParticle.h"
25 #include "AtlasHepMC/GenVertex.h"
26 
27 // root
28 #include "TTree.h"
29 #include <TH1D.h>
30 #include <TGraphAsymmErrors.h>
31 
32 
33 #include <cmath>
34 
35 
36 // constructor
37 Trk::TruthNtupleTool::TruthNtupleTool(
38  const std::string& t,
39  const std::string& n,
40  const IInterface* p )
41  :
42  AthAlgTool(t,n,p),
43  m_fillJets(false),
44  m_nt(nullptr),
45  m_numberOfTreeEntries(0),
46  m_runNumber{},
47  m_eventNumber{},
48  m_TrackLinkIndex{},
49  m_mc_prob{},
50  m_classifications{},
51  m_mc_d0{},
52  m_mc_z0{},
53  m_mc_phi0{},
54  m_mc_theta{},
55  m_mc_qOverP{},
56  m_mc_qOverPt{},
57  m_mc_eta{},
58 
59  m_mc_particleID{},
60  m_mc_barcode{},
61  m_mc_energy{},
62  m_mc_jetLinkIndex{},
63  m_mc_prodR{},
64  m_mc_prodz{}
65 {
66  declareInterface<ITruthNtupleTool>(this);
67 
68  m_etaBins.push_back(0.0);
69  m_etaBins.push_back(0.8);
70  m_etaBins.push_back(1.6);
71  m_etaBins.push_back(2.5);
72  // Declare the properties
73  declareProperty("NtupleFileName", m_ntupleFileName = "TRKVAL/Validation", "Ntuple file handle");
74  declareProperty("NtupleTreeName", m_ntupleTreeName = "Truth", "Name of the ntuple tree");
75 }
76 
77 // destructor
78 Trk::TruthNtupleTool::~TruthNtupleTool() = default;
79 
80 
81 
82 
83 // initialize
84 StatusCode Trk::TruthNtupleTool::initialize() {
85 
86  ATH_CHECK( m_evt.initialize() );
87 
88  if (m_etaBins.size()<2) {
89  ATH_MSG_ERROR ("Vector of eta bins too small");
90  return StatusCode::FAILURE;
91  }
92 
93  // ---------------------------
94  // retrive pointer to THistSvc
95  ITHistSvc *tHistSvc;
96  StatusCode sc = service("THistSvc", tHistSvc);
97  if (sc.isFailure()) {
98  ATH_MSG_ERROR ("Unable to retrieve pointer to THistSvc");
99  return sc;
100  }
101  // ---------------------------
102  // create tree and register it to THistSvc
103  m_nt = new TTree(TString(m_ntupleTreeName), "Track Validation Truth");
104  // NB: we must not delete the tree, this is done by THistSvc
105  std::string fullNtupleName = "/"+m_ntupleFileName+"/"+m_ntupleTreeName;
106  sc = tHistSvc->regTree(fullNtupleName, m_nt);
107  if (sc.isFailure()) {
108  ATH_MSG_ERROR ("Unable to register TTree : " << m_ntupleTreeName);
109  return sc;
110  }
111 
112  //-----------------
113  // add items
114  // event info:
115  m_nt->Branch( "RunNumber", &m_runNumber );
116  m_nt->Branch( "EventNumber", &m_eventNumber );
117  m_nt->Branch( "truth_d0", &m_mc_d0 );
118  m_nt->Branch( "truth_z0", &m_mc_z0 );
119  m_nt->Branch( "truth_phi0", &m_mc_phi0 );
120  m_nt->Branch( "truth_theta", &m_mc_theta );
121  m_nt->Branch( "truth_qOverP", &m_mc_qOverP );
122  m_nt->Branch( "truth_qOverPt", &m_mc_qOverPt );
123  m_nt->Branch( "truth_eta", &m_mc_eta );
124  m_nt->Branch( "truth_particleID", &m_mc_particleID );
125  m_nt->Branch( "truth_barcode", &m_mc_barcode );
126  m_nt->Branch( "truth_energy", &m_mc_energy );
127  m_nt->Branch( "truth_prod_R", &m_mc_prodR );
128  m_nt->Branch( "truth_prod_z", &m_mc_prodz );
129 
130  ATH_MSG_VERBOSE ("added items to ntuple");
131 
132  return StatusCode::SUCCESS;
133 
134 }
135 
136 StatusCode Trk::TruthNtupleTool::initBranches(const std::vector<const Trk::ITrackTruthClassifier*>& classifiers, bool fillJets, const std::vector<std::string> trackCollectionNames) {
137  m_trackTruthClassifiers = classifiers;
138  m_classifications.resize(classifiers.size(), 0);
139  m_recoTrackCounts.resize(classifiers.size());
140  m_truthTrackCounts.resize(classifiers.size());
141  m_fillJets = fillJets;
142 
143  if (m_fillJets)
144  m_nt->Branch("JetLinkIndex", &m_mc_jetLinkIndex );// link to MC jets
145 
146  double etaBins[100];
147  for (unsigned int bin = 0; bin < m_etaBins.size(); ++bin ) {
148  etaBins[bin] = m_etaBins[bin];
149  }
150  for (unsigned int classifierIndex = 0; classifierIndex < classifiers.size(); ++classifierIndex ) {
151  if (!(classifiers[classifierIndex])) {
152  ATH_MSG_ERROR ("Got NULL pointer to truth classifier tool!");
153  return StatusCode::FAILURE;
154  }
155  std::string classifierName = classifiers[classifierIndex]->nameOfClassifier();
156  // add branch for classifier
157  ATH_MSG_INFO ("using classifier " << classifierName);
158  std::string branchName = "Class_" + classifierName;
159  std::string branchDescription = "Class_" + classifierName + "/I";
160  m_nt->Branch(branchName.c_str(), &(m_classifications[classifierIndex]), branchDescription.c_str());
161  // init statistic counters
162  for (unsigned int clIndex = 0; clIndex < classifiers[classifierIndex]->numberOfClassifiers(); clIndex++) {
163  m_truthTrackCounts[classifierIndex].push_back( new TH1D(
164  ("true"+classifiers[classifierIndex]->classificationAsString(clIndex)).c_str(),
165  ("true"+classifiers[classifierIndex]->classificationAsString(clIndex)).c_str(),
166  m_etaBins.size()-1,
167  etaBins
168  ) );
169  // TODO: do statistics for each input track collection
170  m_recoTrackCounts[classifierIndex].push_back( new TH1D(
171  ("reco"+classifiers[classifierIndex]->classificationAsString(clIndex)).c_str(),
172  ("reco"+classifiers[classifierIndex]->classificationAsString(clIndex)).c_str(),
173  m_etaBins.size()-1,
174  etaBins
175  ) );
176  }
177  }
178 
179  // add branches with links to the reconstructed tracks:
180  if (!m_TrackLinkIndex.empty()) {
181  ATH_MSG_ERROR ("Track links seem to be already initialized! Do NOT call initBranches() twice!");
182  return StatusCode::FAILURE;
183  }
184  m_TrackLinkIndex.resize(trackCollectionNames.size(), nullptr);
185  m_mc_prob.resize(trackCollectionNames.size(), nullptr);
186  for (unsigned int trackColIndex = 0; trackColIndex < trackCollectionNames.size(); ++trackColIndex ) {
187  m_TrackLinkIndex[trackColIndex] = new std::vector<unsigned int>();
188  std::string branchName = "TrackLinkIndex_" + trackCollectionNames[trackColIndex];
189  m_nt->Branch(branchName.c_str(), &(m_TrackLinkIndex[trackColIndex])); // link to tracks
190  m_mc_prob[trackColIndex] = new std::vector<float>();
191  branchName = "truth_prob_" + trackCollectionNames[trackColIndex];
192  m_nt->Branch(branchName.c_str(), &(m_mc_prob[trackColIndex])); // matching probability
193  }
194  return StatusCode::SUCCESS;
195 }
196 
200 StatusCode Trk::TruthNtupleTool::finalize() {
201 
202  ATH_MSG_INFO ("start finalize() in " << name());
203 
204  // retrive pointer to THistSvc
205  ITHistSvc *tHistSvc;
206  StatusCode sc = service("THistSvc", tHistSvc);
207  if (sc.isFailure()) {
208  ATH_MSG_ERROR ("Unable to retrieve pointer to THistSvc");
209  return sc;
210  }
211 
212  ATH_MSG_INFO ( "Efficiencies:" );
213 
214  //Double_t* efficiencyValues;
215  for (unsigned int classifierIndex = 0; classifierIndex < m_trackTruthClassifiers.size(); ++classifierIndex ) {
216  for (unsigned int clIndex = 0; clIndex < m_trackTruthClassifiers[classifierIndex]->numberOfClassifiers(); clIndex++) {
217  // output:
218 
219  msg() << std::resetiosflags(std::ios::right) << std::setiosflags(std::ios::left) << std::setw(15) << m_trackTruthClassifiers[classifierIndex]->classificationAsString(clIndex);
220  msg() << std::resetiosflags(std::ios::left) << std::setiosflags(std::ios::right);
221  for (int etaBin = 0; etaBin < m_recoTrackCounts[classifierIndex][clIndex]->GetNbinsX(); etaBin++) {
222  double truth = m_truthTrackCounts[classifierIndex][clIndex]->GetBinContent(etaBin+1);
223  double reco = m_recoTrackCounts[classifierIndex][clIndex]->GetBinContent(etaBin+1);
224  msg() << std::resetiosflags(std::ios::fixed | std::ios::showpoint) << std::setprecision(8) << std::setw(8) << reco << "/" << std::setw(8) << truth;
225  msg() << std::setiosflags(std::ios::fixed | std::ios::showpoint | std::ios::right ) << std::setw(7)
226  << " (" << std::setprecision(4) << (truth>0. ? reco/truth : 0.) << ")";
227  }
228  // calc efficiency:
229  TGraphAsymmErrors* effPlot = new TGraphAsymmErrors(m_recoTrackCounts[classifierIndex][clIndex], m_truthTrackCounts[classifierIndex][clIndex], "");
230  // add graph to ROOT file
231  std::string graphName = "trackEff_" + m_trackTruthClassifiers[classifierIndex]->classificationAsString(clIndex);
232  effPlot->SetNameTitle(graphName.c_str(), graphName.c_str());
233  StatusCode sc = tHistSvc->regGraph("/"+m_ntupleFileName+"/"+graphName, effPlot);
234  if(sc.isFailure()){
235  ATH_MSG_ERROR ("ROOT Graph registration failed");
236  return sc;
237  }
238  msg() << "\n";
239  delete m_recoTrackCounts[classifierIndex][clIndex];
240  m_recoTrackCounts[classifierIndex][clIndex] = nullptr;
241  delete m_truthTrackCounts[classifierIndex][clIndex];
242  m_truthTrackCounts[classifierIndex][clIndex] = nullptr;
243  }
244  }
245  msg() << endmsg;
246 
247  // delete track links and probs
248  for (unsigned int trackColIndex = 0; trackColIndex < m_TrackLinkIndex.size(); ++trackColIndex ) {
249  delete m_TrackLinkIndex[trackColIndex];
250  delete m_mc_prob[trackColIndex];
251  }
252  return StatusCode::SUCCESS;
253 }
254 
258 StatusCode Trk::TruthNtupleTool::writeTruthData (
259  const std::vector< Trk::ValidationTrackTruthData >& truthData) {
260 
261  ATH_MSG_DEBUG ("in writeTruthData(...) with ValTrackTruthData size = "<<truthData.size());
262  // ---------------------------------------
263  // fill event data
264  SG::ReadHandle<xAOD::EventInfo> evt(m_evt);
265  if (!evt.isValid()) {
266  ATH_MSG_WARNING ("Could not retrieve event info");
267  return StatusCode::FAILURE;
268  }
269 
270  m_runNumber=evt->runNumber();
271  m_eventNumber=evt->eventNumber();
272 
273  HepMC::ConstGenParticlePtr genParticle{nullptr};
274  const Trk::TrackParameters* truePerigee = nullptr;
275  for (unsigned int index = 0; index < truthData.size(); index++) {
276  genParticle = truthData[index].genParticle;
277  truePerigee = truthData[index].truthPerigee;
278  if (genParticle==nullptr || truePerigee==nullptr) {
279  if (genParticle==nullptr) ATH_MSG_WARNING ("NULL pointer to gen particle at index "<<index<<
280  ", problem with truth selection logic?");
281  else ATH_MSG_DEBUG ("NULL pointer perigee from TruthToTrack. Index is "<<index);
282  m_mc_barcode = 0;
283  for (unsigned int trackColIndex = 0; trackColIndex < m_mc_prob.size(); ++trackColIndex ) {
284  (m_mc_prob[trackColIndex])->clear();
285  (m_TrackLinkIndex[trackColIndex])->clear();
286  }
287  //m_classifications->clear(); do not clear! vector needs to be the same
288  m_mc_particleID = 0;
289  m_mc_energy = 0.;
290 
291  m_mc_d0 = 0.;
292  m_mc_z0 = 0.;
293  m_mc_phi0 = 0.;
294  m_mc_theta = 0.;
295  m_mc_qOverP = 0.;
296  m_mc_qOverPt = 0.;
297  m_mc_eta = 0.;
298  m_mc_prodR = 0.;
299  m_mc_prodz = 0.;
300  } else {
301  m_mc_barcode = HepMC::barcode(genParticle);
302  m_mc_particleID = genParticle->pdg_id();
303  m_mc_energy = genParticle->momentum().e();
304 
305  m_mc_d0 = truePerigee->parameters()[Trk::d0];
306  m_mc_z0 = truePerigee->parameters()[Trk::z0];
307  m_mc_phi0 = truePerigee->parameters()[Trk::phi0];
308  m_mc_theta = truePerigee->parameters()[Trk::theta];
309  m_mc_qOverP = truePerigee->parameters()[Trk::qOverP];
310  m_mc_qOverPt = (std::sin( m_mc_theta ) != 0.) ? m_mc_qOverP / std::sin( m_mc_theta ) : 0.;
311  m_mc_eta = truePerigee->eta();
312  m_mc_prodR = fabs(genParticle->production_vertex()->position().perp());
313  m_mc_prodz = fabs(genParticle->production_vertex()->position().z());
314 
315  for (unsigned int trackColIndex = 0; trackColIndex < truthData[index].truthToTrackIndices.size(); ++trackColIndex ) {
316  (*(m_TrackLinkIndex[trackColIndex])) = truthData[index].truthToTrackIndices[trackColIndex];
317  (*(m_mc_prob[trackColIndex])) = truthData[index].truthToTrackMatchingProbabilities[trackColIndex];
318  }
319  // get classification
320  m_classifications = truthData[index].classifications;
321 
322  if (m_fillJets)
323  m_mc_jetLinkIndex = truthData[index].truthToJetIndex;
324 
325  // do statistics:
326  for (unsigned int classIndex = 0; classIndex < m_classifications.size(); ++classIndex ) {
327  m_truthTrackCounts[classIndex][m_classifications[classIndex]]->Fill(std::fabs(m_mc_eta));
328  if (!truthData[index].truthToTrackIndices[0].empty()) {
329  // TODO: do statistics for each input track collection
330  m_recoTrackCounts[classIndex][m_classifications[classIndex]]->Fill(std::fabs(m_mc_eta));
331  }
332 
333  }
334 
335  }
336  m_nt->Fill();
337  m_numberOfTreeEntries++;
338  }
339  ATH_MSG_DEBUG ("filled " << m_numberOfTreeEntries << " entries, returning...");
340  return StatusCode::SUCCESS;
341 }
342 
343 
344 unsigned int Trk::TruthNtupleTool::getNumberOfTreeRecords() const {
345  return m_numberOfTreeEntries;
346 }
python.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
TrackParameters.h
Trk::TruthNtupleTool::getNumberOfTreeRecords
virtual unsigned int getNumberOfTreeRecords() const
Definition: TruthNtupleTool.cxx:344
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:70
index
Definition: index.py:1
GenVertex.h
ConvertOldUJHistosToNewHistos.etaBins
list etaBins
Definition: ConvertOldUJHistosToNewHistos.py:145
TH1D
Definition: rootspy.cxx:342
bin
Definition: BinsDiffFromStripMedian.h:43
Trk::z0
@ z0
Definition: ParamDefs.h:70
LArG4FSStartPointFilter.evt
evt
Definition: LArG4FSStartPointFilter.py:42
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
Trk::TruthNtupleTool::writeTruthData
virtual StatusCode writeTruthData(const std::vector< Trk::ValidationTrackTruthData > &truthData)
write track data to ntuple
Definition: TruthNtupleTool.cxx:258
empty
bool empty(TH1 *h)
Definition: computils.cxx:294
GenParticle.h
Trk::TruthNtupleTool::initBranches
virtual StatusCode initBranches(const std::vector< const Trk::ITrackTruthClassifier * > &classifiers, bool, const std::vector< std::string > trackCollectionNames)
Definition: TruthNtupleTool.cxx:136
AthenaPoolTestRead.sc
sc
Definition: AthenaPoolTestRead.py:27
Trk::TruthNtupleTool::TruthNtupleTool
TruthNtupleTool(const std::string &, const std::string &, const IInterface *)
Definition: TruthNtupleTool.cxx:37
xAOD::etaBin
setSAddress setEtaMS setDirPhiMS setDirZMS setBarrelRadius setEndcapAlpha setEndcapRadius setInterceptInner setEtaMap etaBin
Definition: L2StandAloneMuon_v1.cxx:148
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
beamspotman.n
n
Definition: beamspotman.py:731
Trk::theta
@ theta
Definition: ParamDefs.h:72
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
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
HepMC::barcode
int barcode(const T *p)
Definition: Barcode.h:16
xAOD::JetTests::fillJets
void fillJets(JetContainer &jetCont, const std::vector< xAOD::JetFourMom_t > &jet4moms)
Fill input JetContainer with new jets which 4-momentum are given by jet4moms.
Definition: JetFactory.h:78
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Trk::ParametersBase
Definition: ParametersBase.h:55
HepMC::ConstGenParticlePtr
const GenParticle * ConstGenParticlePtr
Definition: GenParticle.h:38
TrackTruth.h
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:195
Trk::d0
@ d0
Definition: ParamDefs.h:69
plotBeamSpotVxVal.bin
int bin
Definition: plotBeamSpotVxVal.py:83
Trk::TruthNtupleTool::~TruthNtupleTool
~TruthNtupleTool()
DeMoScan.index
string index
Definition: DeMoScan.py:362
VKalVrtAthena::varHolder_detail::clear
void clear(T &var)
Definition: NtupleVars.h:48
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
Trk::qOverP
@ qOverP
perigee
Definition: ParamDefs.h:73
Trk::TruthNtupleTool::initialize
StatusCode initialize()
Definition: TruthNtupleTool.cxx:84
declareProperty
#define declareProperty(n, p, h)
Definition: BaseFakeBkgTool.cxx:15
Trk::ParametersBase::eta
double eta() const
Access method for pseudorapidity - from momentum.
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
AthAlgTool
Definition: AthAlgTool.h:26
Trk::TruthNtupleTool::finalize
StatusCode finalize()
finalize
Definition: TruthNtupleTool.cxx:200
TruthNtupleTool.h
python.AutoConfigFlags.msg
msg
Definition: AutoConfigFlags.py:7
Trk::phi0
@ phi0
Definition: ParamDefs.h:71
MuonSegmentReaderConfig.reco
reco
Definition: MuonSegmentReaderConfig.py:133
Generated on Tue Jul 2 2024 21:32:00 for ATLAS Offline Software by   doxygen 1.8.18