From 2256039cffb7ace706a4ec9f06135aca96dd354b Mon Sep 17 00:00:00 2001
From: Lottie Cavanagh <ccavanagh@hep.ph.liv.ac.uk>
Date: Thu, 1 Jul 2021 14:00:46 +0100
Subject: [PATCH] Adding histograms and instructions for ECAL simulation
---
.../SimHitExample/CaloREADME.md | 47 +++++++
.../SimHitExample/src/CaloSimHitAlg.cxx | 118 +++++++++++++++---
.../SimHitExample/src/CaloSimHitAlg.h | 18 ++-
3 files changed, 165 insertions(+), 18 deletions(-)
create mode 100644 Control/CalypsoExample/SimHitExample/CaloREADME.md
diff --git a/Control/CalypsoExample/SimHitExample/CaloREADME.md b/Control/CalypsoExample/SimHitExample/CaloREADME.md
new file mode 100644
index 00000000..a93c9216
--- /dev/null
+++ b/Control/CalypsoExample/SimHitExample/CaloREADME.md
@@ -0,0 +1,47 @@
+
+# Running the ECAL simulation
+
+Use runEcal.py to specify particle type, energy, mass, z position in calo, theta and phi ranges etc
+runEcal.py location: calypso/Simulation/G4Faser/G4FaserAlg/test/runEcal.py
+ - ConfigFlags.Output.HITSFileName = ""
+specifies the output name of the HITS file
+ - pg.sampler.pid = -11
+particle ID (-11 = electron)
+ - pg.sampler.mom = PG.EThetaMPhiSampler(energy=100*GeV, theta=[0,pi/20], phi=[0,2*pi], mass=0.511)
+set particle energy, theta, phi and mass
+ - pg.sampler.pos = PG.PosSampler(x=[-5.0, 5.0], y=[-5.0, 5.0], z=2744.1, t=0.0)
+set position of beam in terms of x and y coordinates
+z = 2744.1 (mm) corresponds to calorimeter face
+x=60.9 and y=49.6 corresponds to centre of top left module
+ - sys.exit(int(acc.run(maxEvents=1000).isFailure()))
+set the number of events to simulate
+1000 event simulations at 1 TeV for electrons and photons need to be split
+change the random seed in this case
+ - pg.randomSeed = 123456
+
+Run this file in Calypso's run directory to genetate HITS file of simulation.
+
+# Plotting Histograms in CaloSimHitAlg
+
+Once the HITS file has been generated from runEcal.py:
+ - use CaloSimHitExample_jobOptions.py to loop over HITS and generate HistoFile
+ - the histogram plotting is specified in CaloSimHitAlg.h and CaloSimHitAlg.cxx
+
+CaloSimHitAlg.cxx plots variations of 6 histograms:
+ - eLoss - energy deposited per hit in the 4 ECAL modules (GeV)
+ - eLossBR - eLoss in bottom right ECAL module
+ - eLossTR - eLoss in top right ECAL module
+ - eLossBL - eLoss in bottom left ECAL module
+ - eLossTL - eLoss in top left ECAL module
+ - eLossTot - total energy deposited by particle across all 4 ECAL modules (GeV)
+ - eLossBRTot - eLossTot in bottom right ECAL module
+ - eLossTRTot - eLossTot in top right ECAL module
+ - eLossBLTot - eLossTot in bottom left ECAL module
+ - eLossTLTot - eLossTot in top left ECAL module
+ - module - cross section of hits in the calo, where (0,0), (1,0), (0,1) and (1,1) correspond to ECAL module centres
+ - modulePos - cross section of hits in the calo, showing their position in mm
+ - meanTime - events in the calo as a function of mean time (ns)
+ - weightedTime - energy weighted events in the calo as a function of mean time (ns)
+
+There are also fractional histograms for eLoss and eLossTot that give the deposited energy as a fraction of injected energy.
+Due to the histograms not depending on injected energy of simulated particle, may need to Rebin() when drawing eLoss and eLossTot histograms.
diff --git a/Control/CalypsoExample/SimHitExample/src/CaloSimHitAlg.cxx b/Control/CalypsoExample/SimHitExample/src/CaloSimHitAlg.cxx
index f99ce4b4..fb4049f7 100644
--- a/Control/CalypsoExample/SimHitExample/src/CaloSimHitAlg.cxx
+++ b/Control/CalypsoExample/SimHitExample/src/CaloSimHitAlg.cxx
@@ -14,22 +14,46 @@ 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;Energy Loss [GeV];Events", 1000000, 0, 100); // 100000 bins, can rebin if needed - 1000 GeV inj - max 100 per event
+ m_elossTL = new TH1D("eLossTL", "Calo Hit Energy Loss Top Left Module;Energy Loss [GeV];Events", 100000, 0, 100);
+ m_elossTR = new TH1D("eLossTR", "Calo Hit Energy Loss Top Right Module;Energy Loss [GeV];Events", 100000, 0, 100);
+ m_elossBR = new TH1D("eLossBR", "Calo Hit Energy Loss Bottom Right Module;Energy Loss [GeV];Events", 100000, 0, 100);
+ m_elossBL = new TH1D("eLossBL", "Calo Hit Energy Loss Bottom Left Module;Energy Loss [GeV];Events", 100000, 0, 100);
+// 1 million bins
+ m_elossTot = new TH1D("eLossTot", "Total Energy Deposited in Calorimeter;Deposited Energy [GeV];Events", 1000000, 0, 1000);
+ m_elossTLTot = new TH1D("eLossTLTot", "Total Energy Deposited in Top Left Module;Deposited Energy [GeV];Events", 1000000, 0, 1000);
+ m_elossTRTot = new TH1D("eLossTRTot", "Total Energy Deposited in Top Right Module;Deposited Energy [GeV];Events", 1000000, 0, 1000);
+ m_elossBRTot = new TH1D("eLossBRTot", "Total Energy Deposited in Bottom Right Module;Deposited Energy [GeV];Events", 1000000, 0, 1000);
+ m_elossBLTot = new TH1D("eLossBLTot", "Total Energy Deposited in Bottom Left Module;Deposited Energy [GeV];Events", 1000000, 0, 1000);
- m_eloss = new TH1D("eLoss", "Calo Hit Energy Loss Fraction;Fraction;Events", 1000, 0, 0.0001);
- m_elossTL = new TH1D("eLossTL", "Calo Hit Energy Loss Fraction Top Left Module;Fraction;Events", 1000, 0, 0.0001);
- m_elossTR = new TH1D("eLossTR", "Calo Hit Energy Loss Fraction Top Right Module;Fraction;Events", 1000, 0, 0.0001);
- m_elossBR = new TH1D("eLossBR", "Calo Hit Energy Loss Fraction Bottom Right Module;Fraction;Events", 1000, 0, 0.0001);
- m_elossBL = new TH1D("eLossBL", "Calo Hit Energy Loss Fraction Bottom Left Module;Fraction;Events", 1000, 0, 0.0001);
+ m_module = new TH2D("module;x;y", "Calo Hit Module", 2, -0.5, 1.5, 2, -0.5, 1.5 );
+ // m_modulePos = new TH2D("modulePos;x;y", "Calo Hit Module Position", 13, -130, 130, 13, -130, 130 );
+ // x4 original number of bins
+ m_modulePos = new TH2D("modulePos;x;y", "Calo Hit Module Position", 52, -130, 130, 52, -130, 130 );
- m_elossTot = new TH1D("eLossTot", "Total Energy Fraction Deposited in Calorimeter;Fraction;Events", 1000, 0, 1);
- m_elossTLTot = new TH1D("eLossTLTot", "Total Energy Fraction Deposited in Top Left Module;Fraction;Events", 1000, 0, 1);
- m_elossTRTot = new TH1D("eLossTRTot", "Total Energy Fraction Deposited in Top Right Module;Fraction;Events", 1000, 0, 1);
- m_elossBRTot = new TH1D("eLossBRTot", "Total Energy Fraction Deposited in Bottom Right Module;Fraction;Events", 1000, 0, 1);
- m_elossBLTot = new TH1D("eLossBLTot", "Total Energy Fraction Deposited in Bottom Left Module;Fraction;Events", 1000, 0, 1);
+ // as a fraction of injected energy (when known)
- m_module = new TH2D("module;x;y", "Calo Hit Module", 2, -0.5, 1.5, 2, -0.5, 1.5 );
- m_modulePos = new TH2D("modulePos;x;y", "Calo Hit Module Position", 13, -130, 130, 13, -130, 130 );
-
+ m_eloss_frac = new TH1D("eLoss_frac", "Calo Hit Energy Loss Fraction;Fraction;Events", 100000, 0, 0.1);
+ m_elossTL_frac = new TH1D("eLossTL_frac", "Calo Hit Energy Loss Fraction Top Left Module;Fraction;Events", 100000, 0, 0.01);
+ m_elossTR_frac = new TH1D("eLossTR_frac", "Calo Hit Energy Loss Fraction Top Right Module;Fraction;Events", 100000, 0, 0.01);
+ m_elossBR_frac = new TH1D("eLossBR_frac", "Calo Hit Energy Loss Fraction Bottom Right Module;Fraction;Events", 100000, 0, 0.01);
+ m_elossBL_frac = new TH1D("eLossBL_frac", "Calo Hit Energy Loss Fraction Bottom Left Module;Fraction;Events", 100000, 0, 0.01);
+
+ m_elossTot_frac = new TH1D("eLossTot_frac", "Total Energy Fraction Deposited in Calorimeter;Fraction;Events", 1000000, 0, 1);
+ m_elossTLTot_frac = new TH1D("eLossTLTot_frac", "Total Energy Fraction Deposited in Top Left Module;Fraction;Events", 1000000, 0, 1);
+ m_elossTRTot_frac = new TH1D("eLossTRTot_frac", "Total Energy Fraction Deposited in Top Right Module;Fraction;Events", 1000000, 0, 1);
+ m_elossBRTot_frac = new TH1D("eLossBRTot_frac", "Total Energy Fraction Deposited in Bottom Right Module;Fraction;Events", 1000000, 0, 1);
+ m_elossBLTot_frac = new TH1D("eLossBLTot_frac", "Total Energy Fraction Deposited in Bottom Left Module;Fraction;Events", 1000000, 0, 1);
+
+ // new histograms hits vs time 29/6/21
+
+ m_meanTime = new TH1D("meanTime", "Events vs Mean Time;Time [ns];Events",4000, 0, 2);
+ m_weightedTime = new TH1D("weightedTime", "Weighted Events vs Time;Time [ns];Events",4000, 0, 2); // weighted energy should be in GeV
+
+ ATH_CHECK(histSvc()->regHist("/HIST/meanTime", m_meanTime));
+ ATH_CHECK(histSvc()->regHist("/HIST/weightedTime", m_weightedTime));
ATH_CHECK(histSvc()->regHist("/HIST/eloss", m_eloss));
ATH_CHECK(histSvc()->regHist("/HIST/elossTL", m_elossTL));
@@ -42,10 +66,24 @@ StatusCode CaloSimHitAlg::initialize()
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));
+ // as a fraction of injected energy
+
+ ATH_CHECK(histSvc()->regHist("/HIST/eloss_frac", m_eloss_frac));
+ ATH_CHECK(histSvc()->regHist("/HIST/elossTL_frac", m_elossTL_frac));
+ ATH_CHECK(histSvc()->regHist("/HIST/elossTR_frac", m_elossTR_frac));
+ ATH_CHECK(histSvc()->regHist("/HIST/elossBR_frac", m_elossBR_frac));
+ ATH_CHECK(histSvc()->regHist("/HIST/elossBL_frac", m_elossBL_frac));
+
+ ATH_CHECK(histSvc()->regHist("/HIST/elossTot_frac", m_elossTot_frac));
+ ATH_CHECK(histSvc()->regHist("/HIST/elossTLTot_frac", m_elossTLTot_frac));
+ ATH_CHECK(histSvc()->regHist("/HIST/elossTRTot_frac", m_elossTRTot_frac));
+ ATH_CHECK(histSvc()->regHist("/HIST/elossBRTot_frac", m_elossBRTot_frac));
+ ATH_CHECK(histSvc()->regHist("/HIST/elossBLTot_frac", m_elossBLTot_frac));
+
// initialize data handle keys
ATH_CHECK( m_mcEventKey.initialize() );
ATH_CHECK( m_faserCaloHitKey.initialize() );
@@ -78,8 +116,18 @@ StatusCode CaloSimHitAlg::execute()
if (ev->particles_size() > 0)
{
HepMC::GenEvent::particle_const_iterator p = ev->particles_begin();
- ePrimary = (*p)->momentum().e();
+ //ePrimary = (*p)->momentum().e();
+
+ for ( ; p != ev->particles_end(); ++p) {
+ // Particle gun primaries have BCs starting at 10000, while others have 3 digit values
+ if ((*p)->barcode() > 11000) continue;
+ ePrimary += (*p)->momentum().e();
+ }
+ std::cout << "E Primary " << ePrimary << std::endl;
}
+
+ //ev->print();
+ //ePrimary = 1000000;
// The hit container might be empty because particles missed the modules
if (h_caloHits->size() == 0) return StatusCode::SUCCESS;
@@ -91,18 +139,38 @@ StatusCode CaloSimHitAlg::execute()
// BL = bottom left calorimeter module
// BR = bottom right calorimeter module
- float e(0);;
+ float e(0);
float tl(0);
float tr(0);
float br(0);
float bl(0);
+ float e_frac(0);
+ float tl_frac(0);
+ float tr_frac(0);
+ float br_frac(0);
+ float bl_frac(0);
+
+
+
for (const CaloHit& hit : *h_caloHits) {
//hit.print();
- float ehit = hit.energyLoss()/ePrimary;
+ float ehit_frac = hit.energyLoss()/ePrimary;
+ float ehit = hit.energyLoss()/1000; // to convert from MeV to GeV
+ m_eloss_frac->Fill(ehit_frac);
m_eloss->Fill(ehit);
+ e_frac += ehit_frac;
e += ehit;
+ // add extra time hist
+ float time = hit.meanTime();
+ m_meanTime->Fill(time);
+
+ // add weighted energy hist
+ m_weightedTime->Fill(time,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());
@@ -110,21 +178,29 @@ StatusCode CaloSimHitAlg::execute()
if (hit.getModule() == 0 && hit.getRow() == 1) // 0 1 TL
{
m_elossTL->Fill(ehit);
+ m_elossTL_frac->Fill(ehit_frac);
+ tl_frac += ehit_frac;
tl += ehit;
}
if (hit.getModule() == 1 && hit.getRow() == 1) // 1 1 TR
{
m_elossTR->Fill(ehit);
+ m_elossTR_frac->Fill(ehit_frac);
+ tr_frac += ehit_frac;
tr += ehit;
}
if (hit.getModule() == 1 && hit.getRow() == 0) // 1 0 BR
{
m_elossBR->Fill(ehit);
+ m_elossBR_frac->Fill(ehit_frac);
+ br_frac += ehit_frac;
br += ehit;
}
if (hit.getModule() == 0 && hit.getRow() == 0) // 0 0 BL
{
m_elossBL->Fill(ehit);
+ m_elossBL_frac->Fill(ehit_frac);
+ bl_frac += ehit_frac;
bl += ehit;
}
@@ -136,6 +212,14 @@ StatusCode CaloSimHitAlg::execute()
m_elossBRTot->Fill(br);
m_elossBLTot->Fill(bl);
+ // as a fraction of injected energy
+
+ m_elossTot_frac->Fill(e_frac);
+ m_elossTLTot_frac->Fill(tl_frac);
+ m_elossTRTot_frac->Fill(tr_frac);
+ m_elossBRTot_frac->Fill(br_frac);
+ m_elossBLTot_frac->Fill(bl_frac);
+
return StatusCode::SUCCESS;
}
diff --git a/Control/CalypsoExample/SimHitExample/src/CaloSimHitAlg.h b/Control/CalypsoExample/SimHitExample/src/CaloSimHitAlg.h
index 6f34f71f..2444115e 100644
--- a/Control/CalypsoExample/SimHitExample/src/CaloSimHitAlg.h
+++ b/Control/CalypsoExample/SimHitExample/src/CaloSimHitAlg.h
@@ -21,6 +21,9 @@ class CaloSimHitAlg : public AthHistogramAlgorithm
TH2* m_module;
TH2* m_modulePos;
TH1* m_elossTot;
+
+ TH1* m_meanTime;
+ TH1* m_weightedTime;
// add new histograms
// TL = top left calorimeter module
@@ -35,7 +38,20 @@ class CaloSimHitAlg : public AthHistogramAlgorithm
TH1* m_elossTLTot;
TH1* m_elossTRTot;
TH1* m_elossBRTot;
- TH1* m_elossBLTot;
+ TH1* m_elossBLTot;
+
+ // add in fractional histograms
+
+ TH1* m_eloss_frac;
+ TH1* m_elossTL_frac;
+ TH1* m_elossTR_frac;
+ TH1* m_elossBR_frac;
+ TH1* m_elossBL_frac;
+ TH1* m_elossTot_frac;
+ TH1* m_elossTLTot_frac;
+ TH1* m_elossTRTot_frac;
+ TH1* m_elossBRTot_frac;
+ TH1* m_elossBLTot_frac;
// Read handle keys for data containers
// Any other event data can be accessed identically
--
GitLab