Skip to content
Snippets Groups Projects
Commit 1ee9dd99 authored by Carl Gwilliam's avatar Carl Gwilliam
Browse files

Merge branch 'master-22.0.29-ecal' into 'master'

Add example code to make plots from calorimeter hits

See merge request faser/calypso!122
parents 5f3cce98 9eb0d20b
No related branches found
No related tags found
1 merge request!122Add example code to make plots from calorimeter hits
Pipeline #2475052 passed
atlas_subdir( SimHitExample )
atlas_add_component( SimHitExample
src/SimHitAlg.cxx
src/*.cxx
src/components/SimHitExample_entries.cxx
LINK_LIBRARIES AthenaBaseComps GeneratorObjects TrackerSimEvent ScintSimEvent FaserCaloSimEvent
)
......
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'"]
#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;
}
#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" };
};
#include "../SimHitAlg.h"
#include "../CaloSimHitAlg.h"
DECLARE_COMPONENT( SimHitAlg )
\ No newline at end of file
DECLARE_COMPONENT( SimHitAlg )
DECLARE_COMPONENT( CaloSimHitAlg )
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment