diff --git a/benchmarks/barrel_ecal/config.yml b/benchmarks/barrel_ecal/config.yml index 0612753a3121d41355ff4b8202d88cecfcdbc825..cfc9a88752d7ffa47aff4427243296128969d95a 100644 --- a/benchmarks/barrel_ecal/config.yml +++ b/benchmarks/barrel_ecal/config.yml @@ -1,8 +1,8 @@ -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 diff --git a/benchmarks/barrel_ecal/run_emcal_barrel_electrons.sh b/benchmarks/barrel_ecal/run_emcal_barrel_electrons.sh index b445a8053565f2213eccf0610d01aa6087265ff8..b70563d8273f975e5cc6ae448d440aa85720e02e 100755 --- a/benchmarks/barrel_ecal/run_emcal_barrel_electrons.sh +++ b/benchmarks/barrel_ecal/run_emcal_barrel_electrons.sh @@ -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" diff --git a/benchmarks/barrel_ecal/run_emcal_barrel_energy_scan.sh b/benchmarks/barrel_ecal/run_emcal_barrel_energy_scan.sh index 82172b286fe7768a6193e806a41f9f2de39be8fd..db154b6cbc9a233a3df6b5eb59334c85fde5fb53 100644 --- a/benchmarks/barrel_ecal/run_emcal_barrel_energy_scan.sh +++ b/benchmarks/barrel_ecal/run_emcal_barrel_energy_scan.sh @@ -1,10 +1,16 @@ #!/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" diff --git a/benchmarks/barrel_ecal/scripts/emcal_barrel_energy_scan_analysis.cxx b/benchmarks/barrel_ecal/scripts/emcal_barrel_energy_scan_analysis.cxx new file mode 100644 index 0000000000000000000000000000000000000000..065d121775bb2b18bd3d128f442e87f5729013b9 --- /dev/null +++ b/benchmarks/barrel_ecal/scripts/emcal_barrel_energy_scan_analysis.cxx @@ -0,0 +1,190 @@ +//////////////////////////////////////// +// 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"); +}