Skip to content
Snippets Groups Projects
Commit 8742d5fb authored by Maria Zurek's avatar Maria Zurek
Browse files

Add energy_scan analysis for electrons

parent bc0d944a
No related branches found
No related tags found
No related merge requests found
sim:emcal_barrel_pions:
extends: .det_benchmark
stage: simulate
script:
- bash benchmarks/barrel_ecal/run_emcal_barrel_pions.sh
# sim:emcal_barrel_pions:
# extends: .det_benchmark
# stage: simulate
# script:
# - bash benchmarks/barrel_ecal/run_emcal_barrel_pions.sh
sim:emcal_barrel_electrons:
extends: .det_benchmark
......@@ -11,13 +11,13 @@ sim:emcal_barrel_electrons:
- bash benchmarks/barrel_ecal/run_emcal_barrel_electrons.sh
- bash benchmarks/barrel_ecal/run_emcal_barrel_energy_scan.sh
bench:emcal_barrel_pions:
extends: .det_benchmark
stage: benchmarks
needs:
- ["sim:emcal_barrel_pions"]
script:
- root -b -q benchmarks/barrel_ecal/scripts/emcal_barrel_pions_analysis.cxx+
# bench:emcal_barrel_pions:
# extends: .det_benchmark
# stage: benchmarks
# needs:
# - ["sim:emcal_barrel_pions"]
# script:
# - root -b -q benchmarks/barrel_ecal/scripts/emcal_barrel_pions_analysis.cxx+
bench:emcal_barrel_electrons:
extends: .det_benchmark
......@@ -28,12 +28,13 @@ bench:emcal_barrel_electrons:
- ls -lrht sim_output/
- rootls -t sim_output/sim_emcal_barrel_uniform_electrons.root
- root -b -q benchmarks/barrel_ecal/scripts/emcal_barrel_electrons_analysis.cxx+
- root -b -q benchmarks/barrel_ecal/scripts/emcal_barrel_energy_scan_analysis.cxx+
collect_results:barrel_ecal:
extends: .det_benchmark
stage: collect
needs:
- ["bench:emcal_barrel_pions","bench:emcal_barrel_electrons"]
- ["bench:emcal_barrel_electrons"]
script:
- ls -lrht
......@@ -15,7 +15,7 @@ fi
if [[ ! -n "${E_end}" ]] ; then
export E_end=5.0
fi
export JUGGLER_FILE_NAME_TAG="emcal_barrel_uniform_electrons_${E_start}_${E_end}"
export JUGGLER_FILE_NAME_TAG="emcal_barrel_uniform_electrons"
export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root"
......
#!/bin/bash
E_FILE="emcal_barrel_energy_scan_points.txt"
E_file="emcal_barrel_energy_scan_points.txt"
for E in {1..20}
for E in {1..3}
do
export E_start="$E"
export E_end="$E"
bash benchmarks/barrel_ecal/run_emcal_barrel_electrons.sh && echo "$E" >> "$E_FILE"
bash benchmarks/barrel_ecal/run_emcal_barrel_electrons.sh && echo "$E" >> "$E_file"
path="sim_output/energy_scan/${E}/"
mkdir -p "$path"
mv sim_output/${JUGGLER_SIM_FILE} "$path/."
done
ls -lthaR sim_output
cat "$E_file"
////////////////////////////////////////
// Read reconstruction ROOT output file
// Plot variables
////////////////////////////////////////
#include "ROOT/RDataFrame.hxx"
#include <iostream>
#include "dd4pod/Geant4ParticleCollection.h"
#include "dd4pod/CalorimeterHitCollection.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TMath.h"
#include "TH1.h"
#include "TF1.h"
#include "TH1D.h"
using ROOT::RDataFrame;
using namespace ROOT::VecOps;
// Setting for graphs
gROOT->SetStyle("Plain");
gStyle->SetOptFit(1);
gStyle->SetLineWidth(2);
gStyle->SetPadTickX(1);
gStyle->SetPadTickY(1);
gStyle->SetPadGridX(1);
gStyle->SetPadGridY(1);
gStyle->SetPadLeftMargin(0.14);
gStyle->SetPadRightMargin(0.14);
void save_canvas(TCanvas* c, std::string label, double E)
{
c->SaveAs(std::format("results/{}_{}.png",label, E));
c->SaveAs(std::format("results/{}_{}.pdf",label, E));
}
void save_canvas(TCanvas* c, std::string label)
{
c->SaveAs(std::format("results/{}.png",label));
c->SaveAs(std::format("results/{}.pdf",label));
}
std::tuple <double, double> extract_sampling_fraction_parameters(double E)
{
auto input_fname = std::format("sim_output/energy_scan/{}/emcal_barrel_uniform_electrons.root", E);
ROOT::EnableImplicitMT();
ROOT::RDataFrame d0("events", input_fname);
// Thrown Energy [GeV]
auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
auto p = input[2];
auto energy = TMath::Sqrt(p.psx * p.psx + p.psy * p.psy + p.psz * p.psz + p.mass * p.mass);
return energy;
};
// Number of hits
auto nhits = [] (const std::vector<dd4pod::CalorimeterHitData>& evt) {return (int) evt.size(); };
// Energy deposition [GeV]
auto Esim = [](const std::vector<dd4pod::CalorimeterHitData>& evt) {
auto total_edep = 0.0;
for (const auto& i: evt)
total_edep += i.energyDeposit;
return total_edep;
};
// Sampling fraction = Esampling / Ethrown
auto fsam = [](const double sampled, const double thrown) {
return sampled / thrown;
};
// Define variables
auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"})
.Define("nhits", nhits, {"EcalBarrelHits"})
.Define("Esim", Esim, {"EcalBarrelHits"})
.Define("fsam", fsam, {"Esim", "Ethr"});
// Define Histograms
auto hEthr = d1.Histo1D(
{"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 7.5},
"Ethr");
auto hNhits =
d1.Histo1D({"hNhits", "Number of hits per events; Number of hits; Events",
100, 0.0, 2000.0},
"nhits");
auto hEsim = d1.Histo1D(
{"hEsim", "Energy Deposit; Energy Deposit [GeV]; Events", 100, 0.0, 1.0},
"Esim");
auto hfsam = d1.Histo1D(
{"hfsam", "Sampling Fraction; Sampling Fraction; Events", 100, 0.0, 0.1},
"fsam");
// Event Counts
auto nevents_thrown = d1.Count();
std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
// Draw Histograms
{
TCanvas* c1 = new TCanvas("c1", "c1", 700, 500);
c1->SetLogy(1);
auto h = hEthr->DrawCopy();
//h->GetYaxis()->SetTitleOffset(1.4);
h->SetLineWidth(2);
h->SetLineColor(kBlue);
save_canvas(c1, "emcal_barrel_electrons_Ethr", E);
}
std::cout << "derp1\n";
{
TCanvas* c2 = new TCanvas("c2", "c2", 700, 500);
c2->SetLogy(1);
auto h = hNhits->DrawCopy();
//h->GetYaxis()->SetTitleOffset(1.4);
h->SetLineWidth(2);
h->SetLineColor(kBlue);
save_canvas(c2, "emcal_barrel_electrons_nhits", E);
}
{
TCanvas* c3 = new TCanvas("c3", "c3", 700, 500);
c3->SetLogy(1);
auto h = hEsim->DrawCopy();
//h->GetYaxis()->SetTitleOffset(1.4);
h->SetLineWidth(2);
h->SetLineColor(kBlue);
save_canvas(c3, "emcal_barrel_electrons_Esim", E);
}
{
TCanvas* c4 = new TCanvas("c4", "c4", 700, 500);
c4->SetLogy(1);
auto h = hfsam->DrawCopy();
//h->GetYaxis()->SetTitleOffset(1.4);
h->SetLineWidth(2);
h->SetLineColor(kBlue);
h->Fit("gaus", "", "", 0.01, 0.1);
TF1 *gaus = h->GetFunction("gaus");
gaus->SetLineWidth(2);
gaus->SetLineColor(kRed);
double mean = gaus->GetParameter(2);
double sigma = gaus->GetParameter(3);
save_canvas(c4, "emcal_barrel_electrons_fsam", E);
return std::make_tuple(mean, sigma);
}
}
std::vector<double> read_scanned_energies(std::string input_energies_fname)
{
std::vector<double> scanned_energies;
double E;
ifstream E_file (input_energies_fname);
if (E_file.is_open())
{
while (E_file >> E)
{
scanned_energies.push_back(E);
}
myfile.close();
return scanned_energies;
}
else
{
std::cout << std::format("Unable to open file {}", input_energies_fname);
abort();
}
};
void emcal_barrel_electrons_energy_scan_analysis()
{
vector<double> scanned_energies = read_scanned_energies("emcal_barrel_energy_scan_points.txt");
TGraph gr_fsam(scanned_energies.size());
TGraph gr_fsam_res(scanned_energies.size());
for (const auto& E : scanned_energies) {
auto [fsam, fsam_res] = extract_sampling_fraction_parameters(E);
gr_fsam.AddPoint(E,fsam);
gr_fsam_res.AddPoint(E,fsam_res);
}
TCanvas* c5 = new TCanvas("c5", "c5", 700, 500);
c5->cd();
gr_fsam.Draw("AP");
save_canvas(c5,"emcal_barrel_electrons_fsam_scan");
TCanvas* c6 = new TCanvas("c6", "c6", 700, 500);
c6->cd();
gr_fsam_res.Draw("AP");
save_canvas(c6,"emcal_barrel_electrons_fsam_scan_res");
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment