ATLAS Offline Software
HGTD_RDOAnalysis.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 "HGTD_RDOAnalysis.h"
6 #include "StoreGate/ReadHandle.h"
9 
10 #include "TruthUtils/AtlasPID.h"
12 
13 #include "TTree.h"
14 #include "TString.h"
15 
16 #include <algorithm>
17 #include <math.h>
18 #include <functional>
19 #include <iostream>
20 
21 
22 
24  ATH_MSG_DEBUG( "Initializing HGTD_RDOAnalysis" );
25 
26  // This will check that the properties were initialized
27  // properly by job configuration.
31 
32  // Grab HGTDID helper
35 
36  m_tree = new TTree(m_ntupleName.value().c_str(), "HGTD_RDOAnalysis");
37  ATH_CHECK(histSvc()->regTree(m_ntuplePath.value() + m_ntupleName.value(), m_tree));
38  m_tree->Branch("m_rdo_module_layer", &m_rdo_module_layer);
39  m_tree->Branch("m_rdo_module_x", &m_rdo_module_x);
40  m_tree->Branch("m_rdo_module_y", &m_rdo_module_y);
41  m_tree->Branch("m_rdo_module_z", &m_rdo_module_z);
42  m_tree->Branch("m_rdo_module_ID", &m_rdo_module_ID);
43  m_tree->Branch("m_rdo_hit_x", &m_rdo_hit_x);
44  m_tree->Branch("m_rdo_hit_y", &m_rdo_hit_y);
45  m_tree->Branch("m_rdo_hit_z", &m_rdo_hit_z);
46  m_tree->Branch("m_rdo_hit_toa", &m_rdo_hit_toa);
47  m_tree->Branch("m_rdo_hit_sdo_toa", &m_rdo_hit_sdo_toa);
48  m_tree->Branch("m_rdo_hit_sdo_truth_category", &m_rdo_hit_sdo_truth_category);
49  /*
50  // HISTOGRAMS
51 
53  m_h_rdoID = new TH1F("h_rdoID", "rdoID", 100, 0, 10e17);
54  m_h_rdoID->StatOverflows();
55  ATH_CHECK(histSvc()->regHist(m_histPath + m_h_rdoID->GetName(), m_h_rdoID));
56  */
57 
58  return StatusCode::SUCCESS;
59 }
60 
62  ATH_MSG_DEBUG(" In HGTD_RDOAnalysis::execute()" );
63  m_rdo_module_layer.clear();
64  m_rdo_module_x.clear();
65  m_rdo_module_y.clear();
66  m_rdo_module_z.clear();
67  m_rdo_module_ID.clear();
68  m_rdo_hit_x.clear();
69  m_rdo_hit_y.clear();
70  m_rdo_hit_z.clear();
71  m_rdo_hit_toa.clear();
72  m_rdo_hit_sdo_toa.clear();
74 
75  const EventContext& ctx{Gaudi::Hive::currentContext()};
76  // Raw HGTD Data
77  const HGTD_RDO_Container* p_RDO_cont{nullptr};
78  ATH_CHECK(SG::get(p_RDO_cont, m_inputKey, ctx));
79 
80  //Adding SimMap and McEvent here for added truthMatching checks
81  const InDetSimDataCollection* simDataMap{nullptr};
82  ATH_CHECK(SG::get(simDataMap, m_inputTruthKey, ctx));
83 
84  const McEventCollection* mcEventCollection{nullptr};
85  ATH_CHECK(SG::get(mcEventCollection, m_inputMcEventCollectionKey, ctx));
86 
87  bool doTruthMatching = true;
88  const HepMC::GenEvent* hardScatterEvent(nullptr);
89 
90  if (mcEventCollection->size()==0){
91  ATH_MSG_WARNING("Failed to retrieve a nonzero sized truth event collection, disabling truthMatching");
92  doTruthMatching = false;
93  }
94  if(doTruthMatching) hardScatterEvent = mcEventCollection->at(0);
95 
96  // loop over RDO container
97  for ( HGTD_RDO_Container::const_iterator rdoCont_itr = p_RDO_cont->begin(); rdoCont_itr != p_RDO_cont->end(); ++rdoCont_itr ) {
98 
99  const HGTD_RDO_Collection* p_RDO_coll(*rdoCont_itr);
100  const Identifier rdoIDColl((*rdoCont_itr)->identify());
101 
102  const InDetDD::HGTD_DetectorElement *elementColl = m_HGTD_Manager->getDetectorElement(rdoIDColl);
103  InDetDD::SiLocalPosition localPosColl = elementColl->rawLocalPositionOfCell(rdoIDColl);
104  Amg::Vector3D globalPosColl = elementColl->globalPosition(localPosColl);
105 
106 
107  for ( HGTD_RDO_Collection::const_iterator rdo_itr= p_RDO_coll->begin(); rdo_itr != p_RDO_coll->end(); ++rdo_itr ) {
108 
109  const Identifier rdoID((*rdo_itr)->identify());
110 
112  InDetDD::SiLocalPosition localPos_hit = rdo_element->rawLocalPositionOfCell(rdoID);
113  Amg::Vector3D globalPos_hit = rdo_element->globalPosition(localPos_hit);
114 
115  m_rdo_hit_x.push_back(globalPos_hit[Amg::x]);
116  m_rdo_hit_y.push_back(globalPos_hit[Amg::y]);
117  m_rdo_hit_z.push_back(globalPos_hit[Amg::z]);
118 
119  // Move global position storage for RDO module from outer to inner loop
120  m_rdo_module_x.push_back(globalPosColl[Amg::x]);
121  m_rdo_module_y.push_back(globalPosColl[Amg::y]);
122  m_rdo_module_z.push_back(globalPosColl[Amg::z]);
123 
124  m_rdo_module_layer.push_back(m_HGTD_ID->layer(rdoIDColl));
125  m_rdo_module_ID.push_back(rdoIDColl.get_compact());
126 
127  // Get the Time of Arrival which include the tof and the digitisation smearing of 25 ps
128  float rdo_toa = (*rdo_itr)->getTOA();
129  m_rdo_hit_toa.push_back(rdo_toa);
130 
131  // For truth matching studies we need to get the truth nature of each deposit.
132  // For each hit, there can be up to 7 deposit from different nature accessed through the SDO map.
133  if(doTruthMatching){
134  if(simDataMap) {
135  InDetSimDataCollection::const_iterator iter = (*simDataMap).find((*rdo_itr)->identify());
136  std::vector<SdoInfo> sdo_info;
137  if ( iter != (*simDataMap).end() ) {
138  const InDetSimData& sdo = iter->second;
139  const std::vector< InDetSimData::Deposit >& deposits = sdo.getdeposits();
140  for( std::vector< InDetSimData::Deposit >::const_iterator nextdeposit = deposits.begin() ; nextdeposit!=deposits.end(); ++nextdeposit) {
141  SdoInfo sdoi;
142  // The times of the deposit also include the smearing of digi of 25 ps
143  sdoi.time = nextdeposit->second;
144  const HepMcParticleLink& particleLink = nextdeposit->first;
145  if(particleLink.isValid()){
146  HepMC::ConstGenParticlePtr genPart(particleLink.cptr());
147  if(isHSGoodParticle(genPart,hardScatterEvent)) sdoi.truth = 1; // signal
148  else sdoi.truth = 2; // pile-up
149  }
150  sdo_info.push_back(sdoi);
151  }
152  } else {
153  // Secondaries are not saved in the SDO map, hence saving the TOA for them
154  SdoInfo sdoi;
155  sdoi.time = rdo_toa;
156  sdoi.truth = 3; // secondary
157  sdo_info.push_back(sdoi);
158  }
159 
160  // Only the particle that creates the first deposit is considered to be the correct particle
161  std::sort(sdo_info.begin(), sdo_info.end());
162  if (sdo_info.size() > 0) {
163  m_rdo_hit_sdo_toa.push_back(sdo_info.at(0).time);
164  m_rdo_hit_sdo_truth_category.push_back(sdo_info.at(0).truth);
165  }
166  }
167  else {
168  ATH_MSG_ERROR("No HGTD SDO map found!");
169  }
170  }
171  }
172  }
173 
174 
175  m_tree->Fill();
176 
177  return StatusCode::SUCCESS;
178 }
179 
180 bool HGTD_RDOAnalysis::isHSGoodParticle(HepMC::ConstGenParticlePtr particlePtr,const HepMC::GenEvent* hardScatterGenEvent, float min_pt_cut) {
181  bool decision = false;
182  if( MC::isGenStable(particlePtr) and
183  isCharged(particlePtr) and
184  particlePtr->momentum().perp() >= min_pt_cut and
185  particlePtr->momentum().eta() < 4. and
186  particlePtr->parent_event() == hardScatterGenEvent)
187  decision = true;
188 
189  return decision;
190 }
python.PyKernel.retrieve
def retrieve(aClass, aKey=None)
Definition: PyKernel.py:110
AthHistogramAlgorithm::histSvc
const ServiceHandle< ITHistSvc > & histSvc() const
The standard THistSvc (for writing histograms and TTrees and more to a root file) Returns (kind of) a...
Definition: AthHistogramAlgorithm.h:113
HGTD_RDOAnalysis::SdoInfo::truth
int truth
Definition: HGTD_RDOAnalysis.h:28
createLinkingScheme.iter
iter
Definition: createLinkingScheme.py:62
DataModel_detail::const_iterator
Const iterator class for DataVector/DataList.
Definition: DVLIterator.h:82
HGTD_RDOAnalysis::execute
virtual StatusCode execute() override final
Definition: HGTD_RDOAnalysis.cxx:61
HGTD_RDOAnalysis::m_rdo_hit_x
std::vector< float > m_rdo_hit_x
Definition: HGTD_RDOAnalysis.h:64
HGTD_RDOAnalysis::m_rdo_hit_y
std::vector< float > m_rdo_hit_y
Definition: HGTD_RDOAnalysis.h:65
InDetSimData::getdeposits
const std::vector< Deposit > & getdeposits() const
Definition: InDetSimData.h:74
Amg::y
@ y
Definition: GeoPrimitives.h:35
InDetDD::HGTD_DetectorElement
Definition: HGTD_DetectorElement.h:40
HGTD_RDOAnalysis::m_ntuplePath
Gaudi::Property< std::string > m_ntuplePath
Definition: HGTD_RDOAnalysis.h:53
Identifier::get_compact
value_type get_compact() const
Get the compact id.
InDetSimDataCollection
Definition: InDetSimDataCollection.h:25
HGTD_RDO_Container
Definition: HGTD_RDO_Container.h:18
HGTD_RDOAnalysis::m_inputTruthKey
SG::ReadHandleKey< InDetSimDataCollection > m_inputTruthKey
Definition: HGTD_RDOAnalysis.h:43
HGTD_RDOAnalysis::m_ntupleName
Gaudi::Property< std::string > m_ntupleName
Definition: HGTD_RDOAnalysis.h:54
HGTD_ID::layer
int layer(const Identifier &id) const
Definition: HGTD_ID.h:475
HGTD_RDO_Collection
Definition: HGTD_RDO_Collection.h:19
MC::isGenStable
bool isGenStable(const T &p)
Determine if the particle is stable at the generator (not det-sim) level,.
Definition: HepMCHelpers.h:54
AthCommonDataStore< AthCommonMsg< Algorithm > >::detStore
const ServiceHandle< StoreGateSvc > & detStore() const
The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.
Definition: AthCommonDataStore.h:95
InDetDD::SiLocalPosition
Definition: SiLocalPosition.h:31
HGTD_RDOAnalysis::m_HGTD_Manager
const HGTD_DetectorManager * m_HGTD_Manager
Definition: HGTD_RDOAnalysis.h:46
HGTD_RDOAnalysis::m_tree
TTree * m_tree
Definition: HGTD_RDOAnalysis.h:58
HGTD_RDOAnalysis::m_rdo_hit_toa
std::vector< float > m_rdo_hit_toa
Definition: HGTD_RDOAnalysis.h:67
HGTD_RDOAnalysis::m_rdo_module_z
std::vector< float > m_rdo_module_z
Definition: HGTD_RDOAnalysis.h:62
Amg::z
@ z
Definition: GeoPrimitives.h:36
HGTD_RDOAnalysis::SdoInfo
Definition: HGTD_RDOAnalysis.h:26
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
HGTD_DetectorManager::getDetectorElement
const InDetDD::HGTD_DetectorElement * getDetectorElement(const Identifier &id) const
access to individual elements : via Identifier
Definition: HGTD_DetectorManager.cxx:50
HGTD_RDOAnalysis::m_HGTD_Name
Gaudi::Property< std::string > m_HGTD_Name
Definition: HGTD_RDOAnalysis.h:47
SG::get
const T * get(const ReadCondHandleKey< T > &key, const EventContext &ctx)
Convenience function to retrieve an object given a ReadCondHandleKey.
Definition: ReadCondHandle.h:287
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
Amg::x
@ x
Definition: GeoPrimitives.h:34
HGTD_RDOAnalysis::m_rdo_hit_sdo_toa
std::vector< float > m_rdo_hit_sdo_toa
Definition: HGTD_RDOAnalysis.h:68
InDetSimData
Definition: InDetSimData.h:42
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
IdentifiableContainerMT::const_iterator
Definition: IdentifiableContainerMT.h:79
HGTD_RDOAnalysis::m_rdo_module_x
std::vector< float > m_rdo_module_x
Definition: HGTD_RDOAnalysis.h:60
SG::VarHandleKey::initialize
StatusCode initialize(bool used=true)
If this object is used as a property, then this should be called during the initialize phase.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleKey.cxx:103
McEventCollection
This defines the McEventCollection, which is really just an ObjectVector of McEvent objects.
Definition: McEventCollection.h:32
HGTD_RDOAnalysis::isHSGoodParticle
bool isHSGoodParticle(HepMC::ConstGenParticlePtr particlePtr, const HepMC::GenEvent *hardScatterEvent, float min_pt_cut=1000.)
Definition: HGTD_RDOAnalysis.cxx:180
HGTD_RDOAnalysis::m_rdo_module_layer
std::vector< int > m_rdo_module_layer
Definition: HGTD_RDOAnalysis.h:59
SiLocalPosition.h
HGTD_RDOAnalysis::SdoInfo::time
float time
Definition: HGTD_RDOAnalysis.h:27
HGTD_RDOAnalysis::m_inputMcEventCollectionKey
SG::ReadHandleKey< McEventCollection > m_inputMcEventCollectionKey
Definition: HGTD_RDOAnalysis.h:44
HepMC::ConstGenParticlePtr
const GenParticle * ConstGenParticlePtr
Definition: GenParticle.h:38
HGTD_RDOAnalysis::m_rdo_hit_z
std::vector< float > m_rdo_hit_z
Definition: HGTD_RDOAnalysis.h:66
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
SiDetectorElement.h
DataVector::end
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
InDetDD::SolidStateDetectorElementBase::rawLocalPositionOfCell
Amg::Vector2D rawLocalPositionOfCell(const SiCellId &cellId) const
Returns position (center) of cell.
Definition: SolidStateDetectorElementBase.cxx:230
isCharged
bool isCharged(const T &p)
Definition: AtlasPID.h:998
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
InDetDD::SolidStateDetectorElementBase::globalPosition
HepGeom::Point3D< double > globalPosition(const HepGeom::Point3D< double > &localPos) const
transform a reconstruction local position into a global position (inline):
HGTD_RDOAnalysis::m_rdo_module_ID
std::vector< unsigned long long > m_rdo_module_ID
Definition: HGTD_RDOAnalysis.h:63
HGTD_RDOAnalysis::initialize
virtual StatusCode initialize() override final
Definition: HGTD_RDOAnalysis.cxx:23
HGTD_RDOAnalysis::m_HGTDID_Name
Gaudi::Property< std::string > m_HGTDID_Name
Definition: HGTD_RDOAnalysis.h:49
HGTD_RDOAnalysis::m_rdo_module_y
std::vector< float > m_rdo_module_y
Definition: HGTD_RDOAnalysis.h:61
HGTD_RDOAnalysis.h
ReadHandle.h
Handle class for reading from StoreGate.
HGTD_RDOAnalysis::m_HGTD_ID
const HGTD_ID * m_HGTD_ID
Definition: HGTD_RDOAnalysis.h:48
HGTD_RDOAnalysis::m_inputKey
SG::ReadHandleKey< HGTD_RDO_Container > m_inputKey
Definition: HGTD_RDOAnalysis.h:42
HGTD_RDOAnalysis::m_rdo_hit_sdo_truth_category
std::vector< int > m_rdo_hit_sdo_truth_category
Definition: HGTD_RDOAnalysis.h:69
HepMCHelpers.h
DataVector::begin
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
AtlasPID.h
Identifier
Definition: IdentifierFieldParser.cxx:14