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

Add example code to make plots from calorimeter hits

parent 5f3cce98
No related branches found
No related tags found
No related merge requests found
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