ATLAS Offline Software
RecMomentumQualityValidation.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // RecMomentumQualityValidation.cxx, (c) ATLAS Detector software
8 
9 #include "TrkValAlgs/RecMomentumQualityValidation.h"
10 #include "TrkTrack/TrackCollection.h"
11 #include "TrkTruthData/TrackTruthCollection.h"
12 #include "TrkParameters/TrackParameters.h"
13 #include "AtlasHepMC/GenParticle.h"
14 #include "TrkToolInterfaces/ITruthToTrack.h"
15 #include "TrkToolInterfaces/ITrackSelectorTool.h"
16 #include "AtlasDetDescr/AtlasDetectorID.h"
17 #include "TrkEventUtils/IdentifierExtractor.h"
18 
19 
20 //================ Constructor =================================================
21 
22 Trk::RecMomentumQualityValidation::RecMomentumQualityValidation(const std::string& name, ISvcLocator* pSvcLocator):
23  AthAlgorithm(name,pSvcLocator),
24  m_inputTrackCollection("Tracks"),
25  m_trackTruthCollection("TrackTruthCollection"),
26  m_truthToTrack("Trk::TruthToTrack/TruthToTrack"),
27  m_trackSelector(""), //InDet::InDetTrackSelectorTool/InDetTrackSelectorTool),
28  m_idHelper(nullptr),
29  m_nHundred(4),m_nFifty(4),m_nTwenty(4),m_nTen(4),m_nFakeOrLost(4),
30  m_tHundred(4),m_tFifty(4),m_tTwenty(4),m_tTen(4),m_tFakeOrLost(4)
31 {
32 
33  // Get Parameter values from JobOptions file
34  declareProperty("InputTrackCollection", m_inputTrackCollection);
35  declareProperty("TrackTruthCollection", m_trackTruthCollection);
36  declareProperty("TruthToTrackTool", m_truthToTrack);
37  declareProperty("TrackSelectorTool", m_trackSelector);
38 
39 }
40 
41 //================ Destructor =================================================
42 
43 Trk::RecMomentumQualityValidation::~RecMomentumQualityValidation()
44 = default;
45 
46 
47 //================ Initialisation =================================================
48 
49 StatusCode Trk::RecMomentumQualityValidation::initialize()
50 {
51 
52  StatusCode sc = StatusCode::SUCCESS;
53  if (sc.isFailure()) { /* mute sc*/ }
54 
55  // Get the Track Selector Tool
56  if ( !m_trackSelector.empty() ) {
57  sc = m_trackSelector.retrieve();
58  if (sc.isFailure()) {
59  msg(MSG::FATAL) << "Could not retrieve "<< m_trackSelector <<" (to select the tracks which are to be counted) "<< endmsg;
60  msg(MSG::INFO) << "Set the ToolHandle to None if track selection is supposed to be disabled" << endmsg;
61  return sc;
62  }
63  }
64  sc = m_truthToTrack.retrieve();
65  if (sc.isFailure()) {
66  msg(MSG::FATAL) << "Could not retrieve "<< m_truthToTrack << endmsg;
67  return sc;
68  }
69 
70  if (detStore()->retrieve(m_idHelper, "AtlasID").isFailure()) {
71  ATH_MSG_ERROR ("Could not get AtlasDetectorID helper" );
72  return StatusCode::FAILURE;
73  }
74 
75  ATH_MSG_INFO( "initialize() successful in " << name() );
76  return StatusCode::SUCCESS;
77 
78 }
79 
80 //================ Finalisation =================================================
81 
82 StatusCode Trk::RecMomentumQualityValidation::finalize()
83 {
84  // Code entered here will be executed once at the end of the program run.
85 
86  if (msgLvl(MSG::INFO)) printTable();
87 
88  return StatusCode::SUCCESS;
89 }
90 
91 //================ Execution ====================================================
92 
93 StatusCode Trk::RecMomentumQualityValidation::execute()
94 {
95 
96  // Retrieving the Trackcollection specified via m_inputTrackCollection
97  StatusCode sc = StatusCode::SUCCESS;
98  const TrackCollection* trackCollection = nullptr;
99 
100  if (!m_inputTrackCollection.empty()) {
101  sc = evtStore()->retrieve(trackCollection,m_inputTrackCollection);
102  if (sc.isFailure())
103  ATH_MSG_ERROR( "TrackCollection "<<m_inputTrackCollection<<" not found!" );
104  else
105  ATH_MSG_VERBOSE("TrackCollection " << m_inputTrackCollection<<" found." );
106 
107  } else {
108  ATH_MSG_ERROR("No Track collection given!" );
109  return StatusCode::FAILURE;
110  }
111 
112  const TrackTruthCollection* trackTruthCollection = nullptr;
113  if (!m_trackTruthCollection.empty()) {
114  sc = evtStore()->retrieve(trackTruthCollection, m_trackTruthCollection);
115  if (sc.isFailure())
116  ATH_MSG_ERROR( "TruthCollection "<<m_trackTruthCollection<<" not found!" );
117  else
118  ATH_MSG_VERBOSE("TruthCollection " << m_trackTruthCollection<<" found." );
119  } else {
120  ATH_MSG_ERROR("No Truth collection given!" );
121  return StatusCode::FAILURE;
122  }
123 
124  // Looping over the tracks of retrieved trackcollection
125  TrackCollection::const_iterator trackIterator = trackCollection->begin();
126  for (;trackIterator!=trackCollection->end();++trackIterator) {
127 
128  if (!m_trackSelector.empty() && (*trackIterator)!=nullptr &&
129  m_trackSelector->decision(**trackIterator) ) {
130  const Trk::TrackParameters* generatedTrackPerigee(nullptr);
131 
132  // find matching truth particle
133  const TrackTruth* trackTruth = nullptr;
134  HepMC::ConstGenParticlePtr genParticle{nullptr};
135  TrackTruthCollection::const_iterator truthIterator
136  = trackTruthCollection->find( trackIterator - (*trackCollection).begin() );
137  if ( truthIterator == trackTruthCollection->end() ){
138  if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "No matching truth particle found for track" << endmsg;
139  } else {
140  trackTruth = &((*truthIterator).second);
141  if ( !(trackTruth->particleLink().isValid()) ) {
142  if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Link to generated particle information is not there - assuming a lost G4 particle ('fake fake')." << endmsg;
143  // genParticle = m_visibleParticleWithoutTruth; // with pdg_id 0
144  } else {
145 #ifdef HEPMC3
146  genParticle = trackTruth->particleLink().scptr();
147 #else
148  genParticle = trackTruth->particleLink().cptr();
149 #endif
150  if ( genParticle && genParticle->pdg_id() == 0 ) {
151  if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Associated Particle ID " << genParticle->pdg_id()
152  << " does not conform to PDG requirements... ignore it!" << endmsg;
153  genParticle = nullptr;
154  }
155  }
156  }
157  if (genParticle) {
158  if (msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) <<
159  "Associated Particle ID: " << genParticle->pdg_id() << endmsg;
160  // Perform extrapolation to generate perigee parameters
161  if ( genParticle->production_vertex() )
162  generatedTrackPerigee = m_truthToTrack->makePerigeeParameters( genParticle );
163  }
164 
165  // Getting the information if a track has TRT hits or not
166  bool track_has_trthits(false);
167  const Trk::TrackStates *currentTSOSList
168  = (*trackIterator)->trackStateOnSurfaces();
169  DataVector< const Trk::TrackStateOnSurface >::const_iterator itTSOS
170  = currentTSOSList->begin();
171  for (;itTSOS!=currentTSOSList->end();++itTSOS) {
172  const Trk::MeasurementBase* mb = (*itTSOS)->measurementOnTrack();
173  Identifier id = Trk::IdentifierExtractor::extract(mb);
174  if (id.is_valid() && m_idHelper->is_trt(id)) {
175  track_has_trthits=true;
176  break;
177  }
178  }
179  const TrackParameters* reconstructedPerigee = (*trackIterator)->perigeeParameters();
180 
181  if (generatedTrackPerigee != nullptr) {
182 
183  if (!reconstructedPerigee) return StatusCode::FAILURE;
184 
185  double this_eta = reconstructedPerigee->eta();
186 
187  double truth_over_recP = generatedTrackPerigee->parameters()[Trk::qOverP]
188  / reconstructedPerigee->parameters()[Trk::qOverP];
189  if ( (truth_over_recP > 0.9) && (truth_over_recP < 1.1)) {
190  monitorTrackFits(m_nTen,this_eta);
191  monitorTrackFits(m_nTwenty,this_eta);
192  monitorTrackFits(m_nFifty,this_eta);
193  monitorTrackFits(m_nHundred,this_eta);
194  if (track_has_trthits) {
195  monitorTrackFits(m_tTen, this_eta);
196  monitorTrackFits(m_tTwenty, this_eta);
197  monitorTrackFits(m_tFifty, this_eta);
198  monitorTrackFits(m_tHundred, this_eta);
199  }
200 
201  } else if ( (truth_over_recP > 0.8) && (truth_over_recP < 1.2) ) {
202  monitorTrackFits(m_nTwenty,this_eta);
203  monitorTrackFits(m_nFifty,this_eta);
204  monitorTrackFits(m_nHundred,this_eta);
205  if (track_has_trthits) {
206  monitorTrackFits(m_tTwenty, this_eta);
207  monitorTrackFits(m_tFifty, this_eta);
208  monitorTrackFits(m_tHundred, this_eta);
209  }
210  } else if ( (truth_over_recP > 0.5) && (truth_over_recP < 1.5) ) {
211  monitorTrackFits(m_nFifty,this_eta);
212  monitorTrackFits(m_nHundred,this_eta);
213  if (track_has_trthits) {
214  monitorTrackFits(m_tFifty, this_eta);
215  monitorTrackFits(m_tHundred, this_eta);
216  }
217  } else {
218  monitorTrackFits(m_nHundred,this_eta);
219  if (track_has_trthits) monitorTrackFits(m_tHundred, this_eta);
220  }
221 
222  delete generatedTrackPerigee;
223  } else {
224  double this_eta = reconstructedPerigee->eta();
225  monitorTrackFits(m_nFakeOrLost, this_eta);
226  if (track_has_trthits) monitorTrackFits(m_tFakeOrLost, this_eta);
227  }
228 
229  } // if it is from a certain track selector/pattern reco info
230 
231  } //loop trackCollection end
232 
233 
234  return StatusCode::SUCCESS;
235 }
236 
237 //============================================================================================
238 
239 void Trk::RecMomentumQualityValidation::monitorTrackFits(std::vector<unsigned int>& Ntracks,
240  const double& eta) {
241  Ntracks[Trk::RecMomentumQualityValidation::iAll] += 1;
242  if (std::abs(eta) < 0.80 ) ++Ntracks[Trk::RecMomentumQualityValidation::iBarrel];
243  else if (std::abs(eta) < 1.60) ++Ntracks[Trk::RecMomentumQualityValidation::iTransi];
244  else if (std::abs(eta) < 2.40) ++Ntracks[Trk::RecMomentumQualityValidation::iEndcap];
245 }
246 
247 
248 void Trk::RecMomentumQualityValidation::printTable() const {
249 
250  unsigned int iw = 2;
251  std::string d = " ";
252  for (unsigned int iref=100; iref<m_nHundred[iAll]; iref*=10, ++iw, d+=" ") {}
253  std::cout << "---------------------------------------------------------------------------------" << std::endl;
254  std::cout << " "<< name() << " results "
255  << (m_trackSelector.empty() ? " " : "(with track selection)") << std::endl;
256  std::cout << "---------------------------------------------------------------------------------" << std::endl;
257  std::cout << " q/p truth vicinity -- Any "<<d<<" 50% "<<d<<" 20% "<<d<<" 10%"<<d<<"noTruth " << std::endl;
258  std::cout << "---------------------------------------------------------------------------------" << std::endl;
259  std::cout << " total (Si+TRT) :" << std::setiosflags(std::ios::dec) << std::setw(iw+1)
260  << m_nHundred[iAll] << " ("<< std::setw(iw)<<m_tHundred[iAll]<<") "
261  << std::setiosflags(std::ios::dec) << std::setw(iw+1)
262  << m_nFifty[iAll] << " ("<< std::setw(iw)<<m_tFifty[iAll]<<") "
263  << std::setiosflags(std::ios::dec) << std::setw(iw+1)
264  << m_nTwenty[iAll] << " ("<< std::setw(iw)<<m_tTwenty[iAll]<<") "
265  << std::setiosflags(std::ios::dec) << std::setw(iw+1)
266  << m_nTen[iAll] << " ("<< std::setw(iw)<<m_tTen[iAll]<<") "
267  << std::setiosflags(std::ios::dec) << std::setw(iw+1)
268  << m_nFakeOrLost[iAll] << " ("<< std::setw(iw)<<m_tFakeOrLost[iAll]<<") "
269  << std::endl;
270  std::cout << " barrel (Si+TRT) :" << std::setiosflags(std::ios::dec) << std::setw(iw+1)
271  << m_nHundred[iBarrel] << " ("<< std::setw(iw)<<m_tHundred[iBarrel]<<") "
272  << std::setiosflags(std::ios::dec) << std::setw(iw+1)
273  << m_nFifty[iBarrel] << " ("<< std::setw(iw)<<m_tFifty[iBarrel]<<") "
274  << std::setiosflags(std::ios::dec) << std::setw(iw+1)
275  << m_nTwenty[iBarrel] <<" ("<< std::setw(iw)<<m_tTwenty[iBarrel]<<") "
276  << std::setiosflags(std::ios::dec) << std::setw(iw+1)
277  << m_nTen[iBarrel] << " ("<< std::setw(iw)<<m_tTen[iBarrel]<<") "
278  << std::setiosflags(std::ios::dec) << std::setw(iw+1)
279  << m_nFakeOrLost[iBarrel] << " ("<< std::setw(iw)<<m_tFakeOrLost[iBarrel]<<") "
280  << std::endl;
281  std::cout << " transition :" << std::setiosflags(std::ios::dec) << std::setw(iw+1)
282  << m_nHundred[iTransi] << " ("<< std::setw(iw)<<m_tHundred[iTransi]<<") "
283  << std::setiosflags(std::ios::dec) << std::setw(iw+1)
284  << m_nFifty[iTransi] << " ("<< std::setw(iw)<<m_tFifty[iTransi]<<") "
285  << std::setiosflags(std::ios::dec) << std::setw(iw+1)
286  << m_nTwenty[iTransi] <<" ("<< std::setw(iw)<<m_tTwenty[iTransi]<<") "
287  << std::setiosflags(std::ios::dec) << std::setw(iw+1)
288  << m_nTen[iTransi] << " ("<< std::setw(iw)<<m_tTen[iTransi]<<") "
289  << std::setiosflags(std::ios::dec) << std::setw(iw+1)
290  << m_nFakeOrLost[iTransi] << " ("<< std::setw(iw)<<m_tFakeOrLost[iTransi]<<") "
291  << std::endl;
292  std::cout << " endcap (Si+TRT) :" << std::setiosflags(std::ios::dec) << std::setw(iw+1)
293  << m_nHundred[iEndcap] << " ("<< std::setw(iw)<<m_tHundred[iEndcap]<<") "
294  << std::setiosflags(std::ios::dec) << std::setw(iw+1)
295  << m_nFifty[iEndcap] << " ("<< std::setw(iw)<<m_tFifty[iEndcap]<<") "
296  << std::setiosflags(std::ios::dec) << std::setw(iw+1)
297  << m_nTwenty[iEndcap] <<" ("<< std::setw(iw)<<m_tTwenty[iEndcap]<<") "
298  << std::setiosflags(std::ios::dec) << std::setw(iw+1)
299  << m_nTen[iEndcap] << " ("<< std::setw(iw)<<m_tTen[iEndcap]<<") "
300  << std::setiosflags(std::ios::dec) << std::setw(iw+1)
301  << m_nFakeOrLost[iEndcap] << " ("<< std::setw(iw)<<m_tFakeOrLost[iEndcap]<<") "
302  << std::endl;
303  std::cout << "---------------------------------------------------------------------------------" << std::endl << std::endl;
304 }
python.PyKernel.retrieve
def retrieve(aClass, aKey=None)
Definition: PyKernel.py:110
DataModel_detail::const_iterator
Const iterator class for DataVector/DataList.
Definition: DVLIterator.h:82
Trk::RecMomentumQualityValidation::iAll
@ iAll
Definition: RecMomentumQualityValidation.h:62
TrackParameters.h
Trk::RecMomentumQualityValidation::execute
StatusCode execute()
standard Athena-Algorithm method
Definition: RecMomentumQualityValidation.cxx:93
python.Constants.FATAL
int FATAL
Definition: Control/AthenaCommon/python/Constants.py:19
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:79
AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
hist_file_dump.d
d
Definition: hist_file_dump.py:137
Trk::RecMomentumQualityValidation::finalize
StatusCode finalize()
standard Athena-Algorithm method
Definition: RecMomentumQualityValidation.cxx:82
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
GenParticle.h
Trk::RecMomentumQualityValidation::monitorTrackFits
static void monitorTrackFits(std::vector< unsigned int > &, const double &)
Definition: RecMomentumQualityValidation.cxx:239
AthenaPoolTestRead.sc
sc
Definition: AthenaPoolTestRead.py:27
IdentifierExtractor.h
AtlasDetectorID.h
This class provides an interface to generate or decode an identifier for the upper levels of the dete...
Trk::RecMomentumQualityValidation::iBarrel
@ iBarrel
Definition: RecMomentumQualityValidation.h:62
TrackTruthCollection
Definition: TrackTruthCollection.h:21
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
RecMomentumQualityValidation.h
TrackTruthCollection.h
Identifier
Definition: DetectorDescription/Identifier/Identifier/Identifier.h:32
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
TrackCollection.h
Trk::RecMomentumQualityValidation::initialize
StatusCode initialize()
standard Athena-Algorithm method
Definition: RecMomentumQualityValidation.cxx:49
Trk::RecMomentumQualityValidation::printTable
void printTable() const
method: make the output table
Definition: RecMomentumQualityValidation.cxx:248
Trk::RecMomentumQualityValidation::m_inputTrackCollection
std::string m_inputTrackCollection
properties from JobOptions:
Definition: RecMomentumQualityValidation.h:55
Trk::ParametersBase
Definition: ParametersBase.h:55
DataVector< Trk::Track >
TrackTruth::particleLink
const HepMcParticleLink & particleLink() const
Definition: TrackTruth.h:26
AthAlgorithm
Definition: AthAlgorithm.h:47
Trk::RecMomentumQualityValidation::m_truthToTrack
ToolHandle< Trk::ITruthToTrack > m_truthToTrack
Tool handle to Trk::ITruthToTrack tool.
Definition: RecMomentumQualityValidation.h:57
Trk::MeasurementBase
Definition: MeasurementBase.h:58
python.PyKernel.detStore
detStore
Definition: PyKernel.py:41
HepMC::ConstGenParticlePtr
const GenParticle * ConstGenParticlePtr
Definition: GenParticle.h:38
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:195
Trk::RecMomentumQualityValidation::iEndcap
@ iEndcap
Definition: RecMomentumQualityValidation.h:62
Trk::RecMomentumQualityValidation::~RecMomentumQualityValidation
~RecMomentumQualityValidation()
Default Destructor.
DataVector::end
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
Trk::IdentifierExtractor::extract
static void extract(std::vector< Identifier > &ids, const std::vector< const MeasurementBase * > &measurements)
Definition: IdentifierExtractor.cxx:13
TrackTruth
MC particle associated with a reco track + the quality of match.
Definition: TrackTruth.h:14
ITruthToTrack.h
TRT_PAI_physicsConstants::mb
const double mb
1mb to cm2
Definition: TRT_PAI_physicsConstants.h:15
DEBUG
#define DEBUG
Definition: page_access.h:11
Trk::qOverP
@ qOverP
perigee
Definition: ParamDefs.h:73
Trk::RecMomentumQualityValidation::m_trackSelector
ToolHandle< Trk::ITrackSelectorTool > m_trackSelector
Tool handle to Trk::ITrackSelectorTool.
Definition: RecMomentumQualityValidation.h:58
pmontree.printTable
def printTable(compList, opt)
Definition: pmontree.py:325
Trk::ParametersBase::eta
double eta() const
Access method for pseudorapidity - from momentum.
python.Constants.VERBOSE
int VERBOSE
Definition: Control/AthenaCommon/python/Constants.py:14
Trk::RecMomentumQualityValidation::RecMomentumQualityValidation
RecMomentumQualityValidation(const std::string &name, ISvcLocator *pSvcLocator)
Standard Athena-Algorithm Constructor.
Definition: RecMomentumQualityValidation.cxx:22
ITrackSelectorTool.h
Trk::RecMomentumQualityValidation::iTransi
@ iTransi
Definition: RecMomentumQualityValidation.h:62
python.AutoConfigFlags.msg
msg
Definition: AutoConfigFlags.py:7
Trk::RecMomentumQualityValidation::m_trackTruthCollection
std::string m_trackTruthCollection
job option: the truth track collection name
Definition: RecMomentumQualityValidation.h:56
DataVector::begin
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
Generated on Thu Jul 4 2024 21:25:21 for ATLAS Offline Software by   doxygen 1.8.18