diff --git a/Control/CalypsoExample/SimHitExample/CMakeLists.txt b/Control/CalypsoExample/SimHitExample/CMakeLists.txt index 907c494c5bff07dfeb120ecc66165048fe2e693f..77084ea39ec8a8acdf4dae161c8abe0f6d175a94 100644 --- a/Control/CalypsoExample/SimHitExample/CMakeLists.txt +++ b/Control/CalypsoExample/SimHitExample/CMakeLists.txt @@ -1,7 +1,7 @@ atlas_subdir( SimHitExample ) atlas_add_component( SimHitExample - src/SimHitAlg.cxx + src/*.cxx src/components/SimHitExample_entries.cxx LINK_LIBRARIES AthenaBaseComps GeneratorObjects TrackerSimEvent ScintSimEvent FaserCaloSimEvent ) diff --git a/Control/CalypsoExample/SimHitExample/share/CaloSimHitExample_jobOptions.py b/Control/CalypsoExample/SimHitExample/share/CaloSimHitExample_jobOptions.py new file mode 100644 index 0000000000000000000000000000000000000000..e8d3b60077cfdd77ddfb146d2a066b4ae22bfa7f --- /dev/null +++ b/Control/CalypsoExample/SimHitExample/share/CaloSimHitExample_jobOptions.py @@ -0,0 +1,16 @@ +from AthenaCommon.GlobalFlags import globalflags + +globalflags.InputFormat.set_Value_and_Lock('pool') + +import AthenaPoolCnvSvc.ReadAthenaPool + +svcMgr.EventSelector.InputCollections = ["my.HITS.pool.root"] + +alg = CfgMgr.CaloSimHitAlg() +athAlgSeq += alg + +theApp.EvtMax=-1 +alg.McEventCollection = "TruthEvent" + +svcMgr += CfgMgr.THistSvc() +svcMgr.THistSvc.Output += ["HIST DATAFILE='myHistoFile.root' OPT='RECREATE'"] diff --git a/Control/CalypsoExample/SimHitExample/src/CaloSimHitAlg.cxx b/Control/CalypsoExample/SimHitExample/src/CaloSimHitAlg.cxx new file mode 100644 index 0000000000000000000000000000000000000000..dd131caa279def71f8054a3751cc593c9368417d --- /dev/null +++ b/Control/CalypsoExample/SimHitExample/src/CaloSimHitAlg.cxx @@ -0,0 +1,145 @@ +#include "CaloSimHitAlg.h" + + + +CaloSimHitAlg::CaloSimHitAlg(const std::string& name, ISvcLocator* pSvcLocator) +: AthHistogramAlgorithm(name, pSvcLocator) { m_eloss = nullptr; } + +CaloSimHitAlg::~CaloSimHitAlg() { } + +StatusCode CaloSimHitAlg::initialize() + +{ + + // initialize a histogram + // letter at end of TH1 indicated variable type (D double, F float etc) + // first string is root object name, second is histogram title + + m_eloss = new TH1D("eLoss", "Calo Hit Energy Loss fraction", 1000, 0, 0.0001); + m_elossTL = new TH1D("eLossTL", "Calo Hit Energy Loss Fraction Top Left Module", 1000, 0, 0.0001); + m_elossTR = new TH1D("eLossTR", "Calo Hit Energy Loss Fraction Top Right Module", 1000, 0, 0.0001); + m_elossBR = new TH1D("eLossBR", "Calo Hit Energy Loss Fraction Bottom Right Module", 1000, 0, 0.0001); + m_elossBL = new TH1D("eLossBL", "Calo Hit Energy Loss Fraction Bottom Left Module", 1000, 0, 0.0001); + + m_elossTot = new TH1D("eLossTot", "Total Energy Fraction Deposited in Calorimeter", 1000, 0, 1); + m_elossTLTot = new TH1D("eLossTLTot", "Total Energy Fraction Deposited in Top Left Module", 1000, 0, 1); + m_elossTRTot = new TH1D("eLossTRTot", "Total Energy Fraction Deposited in Top Right Module", 1000, 0, 1); + m_elossBRTot = new TH1D("eLossBRTot", "Total Energy Fraction Deposited in Bottom Right Module", 1000, 0, 1); + m_elossBLTot = new TH1D("eLossBLTot", "Total Energy Fraction Deposited in Bottom Left Module", 1000, 0, 1); + + m_module = new TH2D("module", "Calo Hit Module", 2, -0.5, 1.5, 2, -0.5, 1.5 ); + m_modulePos = new TH2D("modulePos", "Calo Hit Module Position", 13, -130, 130, 13, -130, 130 ); + + + ATH_CHECK(histSvc()->regHist("/HIST/eloss", m_eloss)); + ATH_CHECK(histSvc()->regHist("/HIST/elossTL", m_elossTL)); + ATH_CHECK(histSvc()->regHist("/HIST/elossTR", m_elossTR)); + ATH_CHECK(histSvc()->regHist("/HIST/elossBR", m_elossBR)); + ATH_CHECK(histSvc()->regHist("/HIST/elossBL", m_elossBL)); + + ATH_CHECK(histSvc()->regHist("/HIST/elossTot", m_elossTot)); + ATH_CHECK(histSvc()->regHist("/HIST/elossTLTot", m_elossTLTot)); + ATH_CHECK(histSvc()->regHist("/HIST/elossTRTot", m_elossTRTot)); + ATH_CHECK(histSvc()->regHist("/HIST/elossBRTot", m_elossBRTot)); + ATH_CHECK(histSvc()->regHist("/HIST/elossBLTot", m_elossBLTot)); + + ATH_CHECK(histSvc()->regHist("/HIST/modules", m_module)); + ATH_CHECK(histSvc()->regHist("/HIST/modulePos", m_modulePos)); + + // initialize data handle keys + ATH_CHECK( m_mcEventKey.initialize() ); + ATH_CHECK( m_faserCaloHitKey.initialize() ); + ATH_MSG_INFO( "Using GenEvent collection with key " << m_mcEventKey.key()); + ATH_MSG_INFO( "Using Faser CaloHit collection with key " << m_faserCaloHitKey.key()); + + return StatusCode::SUCCESS; +} + +StatusCode CaloSimHitAlg::execute() +{ + // Handles created from handle keys behave like pointers to the corresponding container + SG::ReadHandle<McEventCollection> h_mcEvents(m_mcEventKey); + ATH_MSG_INFO("Read McEventContainer with " << h_mcEvents->size() << " events"); + if (h_mcEvents->size() == 0) return StatusCode::FAILURE; + + SG::ReadHandle<CaloHitCollection> h_caloHits(m_faserCaloHitKey); + ATH_MSG_INFO("Read CaloHitCollection with " << h_caloHits->size() << " hits"); + + // Since we have no pile-up, there should always be a single GenEvent in the container + const HepMC::GenEvent* ev = (*h_mcEvents)[0]; + if (ev == nullptr) + { + ATH_MSG_FATAL("GenEvent pointer is null"); + return StatusCode::FAILURE; + } + ATH_MSG_INFO("Event contains " << ev->particles_size() << " truth particles" ); + + double ePrimary = 0; + if (ev->particles_size() > 0) + { + HepMC::GenEvent::particle_const_iterator p = ev->particles_begin(); + ePrimary = (*p)->momentum().e(); + } + + // The hit container might be empty because particles missed the modules + if (h_caloHits->size() == 0) return StatusCode::SUCCESS; + + // Loop over all hits; print and fill histogram + + // TL = top left calorimeter module + // TR = top right calorimeter module + // BL = bottom left calorimeter module + // BR = bottom right calorimeter module + + float e(0);; + float tl(0); + float tr(0); + float br(0); + float bl(0); + + for (const CaloHit& hit : *h_caloHits) { + //hit.print(); + float ehit = hit.energyLoss()/ePrimary; + m_eloss->Fill(ehit); + e += ehit; + + m_module->Fill(hit.getModule(), hit.getRow()); + m_modulePos->Fill(hit.getModule()==0 ? -121.9/2 + hit.localStartPosition().x() : 121.9/2 + hit.localStartPosition().x(), hit.getRow()==0 ? -121.9/2 + hit.localStartPosition().y() : 121.9/2 + hit.localStartPosition().y()); + + + if (hit.getModule() == 0 && hit.getRow() == 1) // 0 1 TL + { + m_elossTL->Fill(ehit); + tl += ehit; + } + if (hit.getModule() == 1 && hit.getRow() == 1) // 1 1 TR + { + m_elossTR->Fill(ehit); + tr += ehit; + } + if (hit.getModule() == 1 && hit.getRow() == 0) // 1 0 BR + { + m_elossBR->Fill(ehit); + br += ehit; + } + if (hit.getModule() == 0 && hit.getRow() == 0) // 0 0 BL + { + m_elossBL->Fill(ehit); + bl += ehit; + } + + } + + m_elossTot->Fill(e); + m_elossTLTot->Fill(tl); + m_elossTRTot->Fill(tr); + m_elossBRTot->Fill(br); + m_elossBLTot->Fill(bl); + + return StatusCode::SUCCESS; +} + +StatusCode CaloSimHitAlg::finalize() +{ + return StatusCode::SUCCESS; +} diff --git a/Control/CalypsoExample/SimHitExample/src/CaloSimHitAlg.h b/Control/CalypsoExample/SimHitExample/src/CaloSimHitAlg.h new file mode 100644 index 0000000000000000000000000000000000000000..6f34f71f14cb04df34e1c17b5256027d4f0d2a9b --- /dev/null +++ b/Control/CalypsoExample/SimHitExample/src/CaloSimHitAlg.h @@ -0,0 +1,45 @@ +#include "AthenaBaseComps/AthHistogramAlgorithm.h" +#include "GeneratorObjects/McEventCollection.h" +#include "FaserCaloSimEvent/CaloHitCollection.h" +#include <TH1.h> + +/* CaloSimHit reading example - Carl Gwilliam + Lottie Cavanagh, Liverpool */ + +class CaloSimHitAlg : public AthHistogramAlgorithm +{ + public: + CaloSimHitAlg(const std::string& name, ISvcLocator* pSvcLocator); + + virtual ~CaloSimHitAlg(); + + StatusCode initialize(); + StatusCode execute(); + StatusCode finalize(); + + private: + TH1* m_eloss; + TH2* m_module; + TH2* m_modulePos; + TH1* m_elossTot; + + // add new histograms + // TL = top left calorimeter module + // TR = top right calorimeter module + // BL = bottom left calorimeter module + // BR = bottom right calorimeter module + + TH1* m_elossTL; // TL + TH1* m_elossTR; // TR + TH1* m_elossBR; // BR + TH1* m_elossBL; // BL + TH1* m_elossTLTot; + TH1* m_elossTRTot; + TH1* m_elossBRTot; + TH1* m_elossBLTot; + + // Read handle keys for data containers + // Any other event data can be accessed identically + // Note the key names ("GEN_EVENT" or "SCT_Hits") are Gaudi properties and can be configured at run-time + SG::ReadHandleKey<McEventCollection> m_mcEventKey { this, "McEventCollection", "GEN_EVENT" }; + SG::ReadHandleKey<CaloHitCollection> m_faserCaloHitKey { this, "CaloHitCollection", "EcalHits" }; +}; diff --git a/Control/CalypsoExample/SimHitExample/src/components/SimHitExample_entries.cxx b/Control/CalypsoExample/SimHitExample/src/components/SimHitExample_entries.cxx index caba1d271a3b63d25e11c153cc741d8d4e106d18..382f43a0252c7308486f67b7b2fe2a42509b2bad 100644 --- a/Control/CalypsoExample/SimHitExample/src/components/SimHitExample_entries.cxx +++ b/Control/CalypsoExample/SimHitExample/src/components/SimHitExample_entries.cxx @@ -1,3 +1,5 @@ #include "../SimHitAlg.h" +#include "../CaloSimHitAlg.h" -DECLARE_COMPONENT( SimHitAlg ) \ No newline at end of file +DECLARE_COMPONENT( SimHitAlg ) +DECLARE_COMPONENT( CaloSimHitAlg )