diff --git a/Simulation/Tools/RDOAnalysis/src/PixelRDOAnalysis.cxx b/Simulation/Tools/RDOAnalysis/src/PixelRDOAnalysis.cxx index 730da3a8ad9c428b65bf0187e564b4c816b868d4..7143e0ee41647d3f3b98450f2ddc3fc68ef80fe7 100644 --- a/Simulation/Tools/RDOAnalysis/src/PixelRDOAnalysis.cxx +++ b/Simulation/Tools/RDOAnalysis/src/PixelRDOAnalysis.cxx @@ -1,10 +1,12 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration */ #include "PixelRDOAnalysis.h" #include "StoreGate/ReadHandle.h" +#include "GeneratorObjects/McEventCollection.h" + #include "TTree.h" #include "TString.h" @@ -49,7 +51,6 @@ PixelRDOAnalysis::PixelRDOAnalysis(const std::string& name, ISvcLocator *pSvcLoc , m_barcode_vec(0) , m_eventIndex_vec(0) , m_charge_vec(0) - , m_h_rdoID(0) , m_h_rdoWord(0) , m_h_barrelEndcap(0) @@ -80,6 +81,7 @@ PixelRDOAnalysis::PixelRDOAnalysis(const std::string& name, ISvcLocator *pSvcLoc , m_h_ecBCID(0) , m_h_ecLVL1A(0) , m_h_ecLVL1ID(0) + , m_h_TruthMatchedRDOs(nullptr) , m_tree(0) , m_ntupleFileName("/ntuples/file1") @@ -316,6 +318,13 @@ StatusCode PixelRDOAnalysis::initialize() { m_h_charge->StatOverflows(); ATH_CHECK(m_thistSvc->regHist(m_path + m_h_charge->GetName(), m_h_charge)); + m_h_TruthMatchedRDOs = new TH1F("h_TruthMatchedPixelRDOs", "h_TruthMatchedPixelRDOs", 4, 1, 5); + TString truthMatchBinLables[4] = { "All RDOs", "Truth Matched", "HS Matched", "Unmatched" }; + for(unsigned int ibin = 1; ibin < 5; ibin++) { + m_h_TruthMatchedRDOs->GetXaxis()->SetBinLabel(ibin, truthMatchBinLables[ibin-1]); + } + ATH_CHECK(m_thistSvc->regHist(m_path + m_h_TruthMatchedRDOs->GetName(), m_h_TruthMatchedRDOs)); + return StatusCode::SUCCESS; } @@ -356,6 +365,18 @@ StatusCode PixelRDOAnalysis::execute() { // Raw Data SG::ReadHandle<PixelRDO_Container> p_pixelRDO_cont (m_inputKey); + //Adding SimMap and McEvent here for added truthMatching checks + SG::ReadHandle<InDetSimDataCollection> simDataMapPixel (m_inputTruthKey); + SG::ReadHandle<McEventCollection> mcEventCollection("TruthEvent"); + bool doTruthMatching = true; + const HepMC::GenEvent* hardScatterEvent(nullptr); + + if (mcEventCollection->size()==0){ + ATH_MSG_WARNING("Failed to retrieve a nonzero sized truth event collection, disabling truthMatching"); + doTruthMatching = false; + } + if(doTruthMatching) hardScatterEvent = mcEventCollection->at(0); + if(p_pixelRDO_cont.isValid()) { // loop over RDO container PixelRDO_Container::const_iterator rdoCont_itr(p_pixelRDO_cont->begin()); @@ -365,6 +386,30 @@ StatusCode PixelRDOAnalysis::execute() { PixelRDO_Collection::const_iterator rdo_itr(p_pixelRDO_coll->begin()); const PixelRDO_Collection::const_iterator rdo_end(p_pixelRDO_coll->end()); for ( ; rdo_itr != rdo_end; ++rdo_itr ) { + if(doTruthMatching){ + m_h_TruthMatchedRDOs->Fill(1.5); + bool findMatch = false; + if(simDataMapPixel.isValid()){ + InDetSimDataCollection::const_iterator iter = (*simDataMapPixel).find((*rdo_itr)->identify()); + + if ( iter != (*simDataMapPixel).end() ) { + const InDetSimData& sdo = iter->second; + const std::vector< InDetSimData::Deposit >& deposits = sdo.getdeposits(); + std::vector< InDetSimData::Deposit >::const_iterator nextdeposit = deposits.begin(); + std::vector< InDetSimData::Deposit >::const_iterator lastdeposit = deposits.end(); + for( ; nextdeposit!=lastdeposit; ++nextdeposit) { + const HepMcParticleLink& particleLink = nextdeposit->first; + if(particleLink.isValid() && !findMatch){ + const HepMC::GenParticle *genPart(particleLink.cptr()); + if(genPart->parent_event() == hardScatterEvent) m_h_TruthMatchedRDOs->Fill(3.5); + m_h_TruthMatchedRDOs->Fill(2.5); + findMatch = true; + } + } + } + } + if(!findMatch) m_h_TruthMatchedRDOs->Fill(4.5); + } const Identifier rdoID((*rdo_itr)->identify()); const int pixBrlEc(m_pixelID->barrel_ec(rdoID)); const unsigned int rdoWord((*rdo_itr)->getWord()); @@ -432,7 +477,6 @@ StatusCode PixelRDOAnalysis::execute() { } // Sim Data - SG::ReadHandle<InDetSimDataCollection> simDataMapPixel (m_inputTruthKey); if(simDataMapPixel.isValid()) { // loop over SDO container InDetSimDataCollection::const_iterator sdo_itr(simDataMapPixel->begin()); diff --git a/Simulation/Tools/RDOAnalysis/src/PixelRDOAnalysis.h b/Simulation/Tools/RDOAnalysis/src/PixelRDOAnalysis.h index c2371ca72506ad05162a96a57e0abf13c50bbab9..52ff0323b7ca6015009830901c151295b95db62d 100644 --- a/Simulation/Tools/RDOAnalysis/src/PixelRDOAnalysis.h +++ b/Simulation/Tools/RDOAnalysis/src/PixelRDOAnalysis.h @@ -123,6 +123,7 @@ private: TH1* m_h_barcode; TH1* m_h_eventIndex; TH1* m_h_charge; + TH1* m_h_TruthMatchedRDOs; TTree* m_tree; std::string m_ntupleFileName; diff --git a/Simulation/Tools/RDOAnalysis/src/SCT_RDOAnalysis.cxx b/Simulation/Tools/RDOAnalysis/src/SCT_RDOAnalysis.cxx index 217f16dfbc99c839d22fbd9bb5d194fa42422037..383484cdac8b4b5cbd4cbfe2040fd05f496424a5 100644 --- a/Simulation/Tools/RDOAnalysis/src/SCT_RDOAnalysis.cxx +++ b/Simulation/Tools/RDOAnalysis/src/SCT_RDOAnalysis.cxx @@ -5,6 +5,7 @@ #include "SCT_RDOAnalysis.h" #include "InDetIdentifier/SCT_ID.h" +#include "GeneratorObjects/McEventCollection.h" #include "StoreGate/ReadHandle.h" #include "TString.h" @@ -81,6 +82,7 @@ SCT_RDOAnalysis::SCT_RDOAnalysis(const std::string& name, ISvcLocator *pSvcLocat , m_h_eventIndex{nullptr} , m_h_charge{nullptr} , m_h_phi_v_eta_sdo{nullptr} + , m_h_TruthMatchedRDOs{nullptr} , m_tree{nullptr} , m_thistSvc("THistSvc", name) @@ -282,6 +284,15 @@ StatusCode SCT_RDOAnalysis::initialize() { m_h_phi_v_eta_sdo->StatOverflows(); ATH_CHECK(m_thistSvc->regHist(m_path.value() + m_h_phi_v_eta_sdo->GetName(), m_h_phi_v_eta_sdo)); + m_h_TruthMatchedRDOs = new TH1F("h_TruthMatchedSCTRDOs", "h_TruthMatchedSCTRDOs", 4, 1, 5); + TString truthMatchBinLables[4] = { "All RDOs", "Truth Matched", "HS Matched", "Unmatched" }; + for(unsigned int ibin = 1; ibin < 5; ibin++) { + m_h_TruthMatchedRDOs->GetXaxis()->SetBinLabel(ibin, truthMatchBinLables[ibin-1]); + } + ATH_CHECK(m_thistSvc->regHist(m_path + m_h_TruthMatchedRDOs->GetName(), m_h_TruthMatchedRDOs)); + + + return StatusCode::SUCCESS; } @@ -317,7 +328,19 @@ StatusCode SCT_RDOAnalysis::execute() { // RawData SG::ReadHandle<SCT_RDO_Container> p_SCT_RDO_cont (m_inputKey); - if (p_SCT_RDO_cont.isValid()) { + //Adding SimMap and McEvent here for added truthMatching checks + SG::ReadHandle<InDetSimDataCollection> simDataMapSCT (m_inputTruthKey); + SG::ReadHandle<McEventCollection> mcEventCollection("TruthEvent"); + + const HepMC::GenEvent* hardScatterEvent(nullptr); + bool doTruthMatching = true; + if (mcEventCollection->size()==0){ + ATH_MSG_WARNING("Failed to retrieve a nonzero sized truth event collection, disabling truthMatching"); + doTruthMatching = false; + } + if(doTruthMatching) hardScatterEvent = mcEventCollection->at(0); + + if(p_SCT_RDO_cont.isValid()) { // loop over RDO container SCT_RDO_Container::const_iterator rdoCont_itr(p_SCT_RDO_cont->begin()); const SCT_RDO_Container::const_iterator rdoCont_end(p_SCT_RDO_cont->end()); @@ -328,6 +351,30 @@ StatusCode SCT_RDOAnalysis::execute() { const SCT_RDO_Collection::const_iterator rdo_end(p_SCT_RDO_coll->end()); for ( ; rdo_itr != rdo_end; ++rdo_itr ) { + if(doTruthMatching){ + m_h_TruthMatchedRDOs->Fill(1.5); + bool findMatch = false; + if(simDataMapSCT.isValid()){ + InDetSimDataCollection::const_iterator iter = (*simDataMapSCT).find((*rdo_itr)->identify()); + + if ( iter != (*simDataMapSCT).end() ) { + const InDetSimData& sdo = iter->second; + const std::vector< InDetSimData::Deposit >& deposits = sdo.getdeposits(); + std::vector< InDetSimData::Deposit >::const_iterator nextdeposit = deposits.begin(); + std::vector< InDetSimData::Deposit >::const_iterator lastdeposit = deposits.end(); + for( ; nextdeposit!=lastdeposit; ++nextdeposit) { + const HepMcParticleLink& particleLink = nextdeposit->first; + if(particleLink.isValid() && !findMatch){ + const HepMC::GenParticle *genPart(particleLink.cptr()); + if(genPart->parent_event() == hardScatterEvent) m_h_TruthMatchedRDOs->Fill(3.5); + m_h_TruthMatchedRDOs->Fill(2.5); + findMatch = true; + } + } + } + } + if(!findMatch) m_h_TruthMatchedRDOs->Fill(4.5); + } const Identifier rdoID((*rdo_itr)->identify()); const unsigned int rdoWord((*rdo_itr)->getWord()); const int sctBrlEc(m_sctID->barrel_ec(rdoID)); @@ -383,7 +430,6 @@ StatusCode SCT_RDOAnalysis::execute() { } // SimData - SG::ReadHandle<InDetSimDataCollection> simDataMapSCT (m_inputTruthKey); if (simDataMapSCT.isValid()) { // loop over SDO container InDetSimDataCollection::const_iterator sdo_itr(simDataMapSCT->begin()); diff --git a/Simulation/Tools/RDOAnalysis/src/SCT_RDOAnalysis.h b/Simulation/Tools/RDOAnalysis/src/SCT_RDOAnalysis.h index 38f7864104db09c9a2d508dd7ed26daffb3ae8de..7b81fd25cd21e4d51e25c1721425c326b4436665 100644 --- a/Simulation/Tools/RDOAnalysis/src/SCT_RDOAnalysis.h +++ b/Simulation/Tools/RDOAnalysis/src/SCT_RDOAnalysis.h @@ -113,6 +113,8 @@ class SCT_RDOAnalysis : public AthAlgorithm { TH1* m_h_eventIndex; TH1* m_h_charge; TH2* m_h_phi_v_eta_sdo; + TH1* m_h_TruthMatchedRDOs; + TTree* m_tree; StringProperty m_ntupleFileName{this, "NtupleFileName", "/ntuples/file1"};