Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
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;
}