diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 943e1862ada5160e7e91e58e0538cb59ed6eb841..6c2ae4c18459ae9d3589ffb8a48c4710acba584a 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -132,7 +132,6 @@ include: - local: 'benchmarks/tracking_performances/config.yml' - local: 'benchmarks/tracking_performances_dis/config.yml' - local: 'benchmarks/barrel_ecal/config.yml' - - local: 'benchmarks/barrel_hcal/config.yml' - local: 'benchmarks/lfhcal/config.yml' - local: 'benchmarks/zdc/config.yml' - local: 'benchmarks/zdc_lyso/config.yml' @@ -158,7 +157,6 @@ deploy_results: - "collect_results:backgrounds" - "collect_results:backwards_ecal" - "collect_results:barrel_ecal" - - "collect_results:barrel_hcal" - "collect_results:calo_pid" - "collect_results:ecal_gaps" - "collect_results:lfhcal" diff --git a/benchmarks/barrel_hcal/config.yml b/benchmarks/barrel_hcal/config.yml deleted file mode 100644 index fa40a22fa6ff854bed70730fca2ccdc60fb90b99..0000000000000000000000000000000000000000 --- a/benchmarks/barrel_hcal/config.yml +++ /dev/null @@ -1,54 +0,0 @@ -sim:hcal_barrel: - extends: .det_benchmark - stage: simulate - parallel: - matrix: - - PARTICLE: ["piplus", "piminus", "kplus", "kminus", "kshort", "klong", "muon", "antimuon", "proton"] - script: - - bash benchmarks/barrel_hcal/run_hcal_barrel_particles.sh ${PARTICLE} - -sim:hcal_barrel:scan: - extends: .det_benchmark - stage: simulate - parallel: - matrix: - - PARTICLE: ["muon", "antimuon", "proton"] - E: ["0.25", "0.5", "1", "2", "3", "4", "7", "15", "20"] - script: - - bash benchmarks/barrel_hcal/run_hcal_barrel_energy_scan.sh ${PARTICLE} ${E} - -bench:hcal_barrel: - extends: .det_benchmark - stage: benchmarks - needs: - - ["sim:hcal_barrel"] - parallel: - matrix: - - PARTICLE: ["piplus", "piminus", "kplus", "kminus", "kshort", "klong", "muon", "antimuon", "proton"] - script: - - ls -lhtR sim_output/ - - rootls -t sim_output/sim_hcal_barrel_${PARTICLE}.edm4hep.root - - root -b -q 'benchmarks/barrel_hcal/scripts/hcal_barrel_particles_analysis.cxx+("'${PARTICLE}'")' - -bench:hcal_barrel:scan: - extends: .det_benchmark - stage: benchmarks - needs: - - ["sim:hcal_barrel:scan"] - parallel: - matrix: - - PARTICLE: ["muon", "antimuon", "proton"] - script: - - ls -lhtR sim_output/ - - sort -n sim_output/hcal_barrel_energy_scan_points_${PARTICLE}_*.txt > sim_output/hcal_barrel_energy_scan_points_${PARTICLE}.txt - - root -b -q 'benchmarks/barrel_hcal/scripts/hcal_barrel_energy_scan_analysis.cxx+("'${PARTICLE}'")' - -collect_results:barrel_hcal: - extends: .det_benchmark - stage: collect - needs: - - ["bench:hcal_barrel", "bench:hcal_barrel:scan"] - script: - - ls -lrht - - echo " FIX ME" - diff --git a/benchmarks/barrel_hcal/run_hcal_barrel_energy_scan.sh b/benchmarks/barrel_hcal/run_hcal_barrel_energy_scan.sh deleted file mode 100755 index 74282724280364300276177f4edc6b05355d554a..0000000000000000000000000000000000000000 --- a/benchmarks/barrel_hcal/run_hcal_barrel_energy_scan.sh +++ /dev/null @@ -1,31 +0,0 @@ -#!/bin/bash - -export PARTICLE=$1 -shift - -E_file="sim_output/hcal_barrel_energy_scan_points_${PARTICLE}_${CI_JOB_ID}.txt" - -if [ $# -gt 0 ] ; then - E_VALUES=("$@") -else - E_VALUES=(0.25 0.5 1 2 3 4 7 15 20) -fi - -for E in ${E_VALUES[@]} -do - export E_START="$E" - export E_END="$E" - export JUGGLER_FILE_NAME_TAG=hcal_barrel_${PARTICLE}_${E} - bash benchmarks/barrel_hcal/run_hcal_barrel_particles.sh "${PARTICLE}" && echo "$E" >> "$E_file" || exit 1 - path_rootfiles="sim_output/energy_scan/${E}/" - path_plots="results/energy_scan/${E}/" - mkdir -p "$path_rootfiles" - mkdir -p "$path_plots" - ls -lthaR sim_output/ - mv "sim_output/sim_${JUGGLER_FILE_NAME_TAG}.edm4hep.root" "$path_rootfiles/sim_hcal_barrel_${PARTICLE}.edm4hep.root" -done - -ls -lthaR sim_output -cat "$E_file" - -exit 0 diff --git a/benchmarks/barrel_hcal/run_hcal_barrel_particles.sh b/benchmarks/barrel_hcal/run_hcal_barrel_particles.sh deleted file mode 100755 index d9f459f85ef937f46a779749832563cd7a770f7e..0000000000000000000000000000000000000000 --- a/benchmarks/barrel_hcal/run_hcal_barrel_particles.sh +++ /dev/null @@ -1,66 +0,0 @@ -#!/bin/bash - -if [[ ! -n "${DETECTOR}" ]] ; then - export DETECTOR="topside" -fi - -if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then - export JUGGLER_N_EVENTS=100 -fi - -if [[ ! -n "${E_START}" ]] ; then - export E_START=5.0 -fi - -if [[ ! -n "${E_END}" ]] ; then - export E_END=5.0 -fi - -export PARTICLE=$1 -if [[ ! -n "${PARTICLE}" ]] ; then - export PARTICLE="electron" -fi - -if [[ ! -n "${JUGGLER_FILE_NAME_TAG}" ]] ; then - export JUGGLER_FILE_NAME_TAG="hcal_barrel_${PARTICLE}" -fi - -export JUGGLER_GEN_FILE="data/${JUGGLER_FILE_NAME_TAG}.hepmc" - -export JUGGLER_SIM_FILE="sim_output/sim_${JUGGLER_FILE_NAME_TAG}.edm4hep.root" -export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root" - -echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}" -echo "DETECTOR = ${DETECTOR}" - -# Generate the input events -root -b -q "benchmarks/barrel_hcal/scripts/hcal_barrel_particles_gen.cxx+(${JUGGLER_N_EVENTS}, ${E_START}, ${E_END}, \"${PARTICLE}\")" -if [[ "$?" -ne "0" ]] ; then - echo "ERROR running script: generating input events" - exit 1 -fi -# Plot the input events -root -b -q "benchmarks/barrel_hcal/scripts/hcal_barrel_particles_reader.cxx+(\"${PARTICLE}\")" -if [[ "$?" -ne "0" ]] ; then - echo "ERROR running script: plotting input events" - exit 1 -fi - -ls -ltRhL - -ddsim --runType batch \ - -v WARNING \ - --part.minimalKineticEnergy 0.5*GeV \ - --filter.tracker edep0 \ - --numberOfEvents ${JUGGLER_N_EVENTS} \ - --compactFile ${DETECTOR_PATH}/${DETECTOR_CONFIG}.xml \ - --inputFiles ${JUGGLER_GEN_FILE} \ - --outputFile ${JUGGLER_SIM_FILE} - -if [[ "$?" -ne "0" ]] ; then - echo "ERROR running npdet" - exit 1 -fi - -# Directory for plots -mkdir -p results diff --git a/benchmarks/barrel_hcal/scripts/hcal_barrel_common_functions.h b/benchmarks/barrel_hcal/scripts/hcal_barrel_common_functions.h deleted file mode 100644 index 3d3efac578ea88865e97b096cb469c468000c0df..0000000000000000000000000000000000000000 --- a/benchmarks/barrel_hcal/scripts/hcal_barrel_common_functions.h +++ /dev/null @@ -1,23 +0,0 @@ -////////////////////////// -// Common Particle Functions -// M. Scott 05/2021 -////////////////////////// - -std::tuple <int, double> extract_particle_parameters(std::string particle_name) { - if (particle_name == "electron") return std::make_tuple(11, 0.51099895e-3); - if (particle_name == "photon") return std::make_tuple(22, 0.0); - if (particle_name == "positron") return std::make_tuple(-11, 0.51099895e-3); - if (particle_name == "proton") return std::make_tuple(2212, 0.938272); - if (particle_name == "muon") return std::make_tuple(13, 0.1056583745); - if (particle_name == "antimuon") return std::make_tuple(-13, 0.1056583745); - if (particle_name == "pi0") return std::make_tuple(111, 0.1349768); - if (particle_name == "piplus") return std::make_tuple(211, 0.13957039); - if (particle_name == "piminus") return std::make_tuple(-211, 0.13957039); - if (particle_name == "kplus") return std::make_tuple(321, 0.493677); - if (particle_name == "kminus") return std::make_tuple(-321, 0.493677); - if (particle_name == "kshort") return std::make_tuple(310, 0.497648); - if (particle_name == "klong") return std::make_tuple(130, 0.497648); - - std::cout << "wrong particle name" << std::endl; - abort(); -} diff --git a/benchmarks/barrel_hcal/scripts/hcal_barrel_energy_scan_analysis.cxx b/benchmarks/barrel_hcal/scripts/hcal_barrel_energy_scan_analysis.cxx deleted file mode 100644 index 03c8af99d6130dcf57fa7f098167965b106dedec..0000000000000000000000000000000000000000 --- a/benchmarks/barrel_hcal/scripts/hcal_barrel_energy_scan_analysis.cxx +++ /dev/null @@ -1,239 +0,0 @@ -//////////////////////////////////////// -// Read reconstruction ROOT output file -// Plot variables -//////////////////////////////////////// - -#include "ROOT/RDataFrame.hxx" -#include <iostream> -#include <fstream> -#include <fmt/core.h> - -#include "edm4hep/MCParticleCollection.h" -#include "edm4hep/SimCalorimeterHitCollection.h" - -#include "TCanvas.h" -#include "TStyle.h" -#include "TMath.h" -#include "TH1.h" -#include "TF1.h" -#include "TH1D.h" -#include "TGraphErrors.h" - -using ROOT::RDataFrame; -using namespace ROOT::VecOps; - - -void save_canvas(TCanvas* c, std::string label) -{ - c->SaveAs(fmt::format("results/energy_scan/{}.png",label).c_str()); - c->SaveAs(fmt::format("results/energy_scan/{}.pdf",label).c_str()); -} - -void save_canvas(TCanvas* c, std::string label, double E) -{ - std::string label_with_E = fmt::format("{}/{}", E, label); - save_canvas(c, label_with_E); -} -void save_canvas(TCanvas* c, std::string label, std::string E_label) -{ - std::string label_with_E = fmt::format("{}/{}", E_label, label); - save_canvas(c, label_with_E); -} -void save_canvas(TCanvas* c, std::string var_label, std::string E_label, std::string particle_label) -{ - std::string label_with_E = fmt::format("{}/hcal_barrel_{}_{}", E_label, particle_label, var_label); - save_canvas(c, label_with_E); -} - -std::tuple <double, double, double, double> extract_sampling_fraction_parameters(std::string particle_label, std::string E_label) -{ - std::string input_fname = fmt::format("sim_output/energy_scan/{}/sim_hcal_barrel_{}.edm4hep.root", E_label, particle_label); - ROOT::EnableImplicitMT(); - ROOT::RDataFrame d0("events", input_fname); - - // Thrown Energy [GeV] - auto Ethr = [](std::vector<edm4hep::MCParticleData> const& input) { - auto p = input[2]; - auto energy = TMath::Sqrt(p.momentum.x * p.momentum.x + p.momentum.y * p.momentum.y + p.momentum.z * p.momentum.z + p.mass * p.mass); - return energy; - }; - - // Number of hits - auto nhits = [] (const std::vector<edm4hep::SimCalorimeterHitData>& evt) {return (int) evt.size(); }; - - // Energy deposition [GeV] - auto Esim = [](const std::vector<edm4hep::SimCalorimeterHitData>& evt) { - auto total_edep = 0.0; - for (const auto& i: evt) - total_edep += i.energy; - 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 assumptions - auto Ethr_max = 25.0; - auto fsam_estimate = 1.0; - if (d1.HasColumn("EcalBarrelScFiHits")) { - d1 = d1.Define("nhits", nhits, {"EcalBarrelImagingHits"}) - .Define("Esim", Esim, {"EcalBarrelImagingHits"}) - .Define("fsam", fsam, {"Esim", "Ethr"}); - fsam_estimate = 0.1; - } else { - d1 = d1.Define("nhits", nhits, {"EcalBarrelSciGlassHits"}) - .Define("Esim", Esim, {"EcalBarrelSciGlassHits"}) - .Define("fsam", fsam, {"Esim", "Ethr"}); - fsam_estimate = 1.0; - } - - // Define Histograms - auto hEthr = d1.Histo1D( - {"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, Ethr_max}, - "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", 500, 0.0, fsam_estimate * Ethr_max}, - "Esim"); - auto hfsam = d1.Histo1D( - {"hfsam", "Sampling Fraction; Sampling Fraction; Events", 200, 0.0, 2.0 * fsam_estimate}, - "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, "Ethr", E_label, particle_label); - } - - { - 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, "nhits", E_label, particle_label); - } - - { - 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); - double up_fit = h->GetMean() + 5*h->GetStdDev(); - double down_fit = h->GetMean() - 5*h->GetStdDev(); - h->GetXaxis()->SetRangeUser(0.,up_fit); - save_canvas(c3, "Esim", E_label, particle_label); - } - - { - 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); - double up_fit = h->GetMean() + 5*h->GetStdDev(); - double down_fit = h->GetMean() - 5*h->GetStdDev(); - h->Fit("gaus", "", "", down_fit, up_fit); - h->GetXaxis()->SetRangeUser(0.,up_fit); - TF1 *gaus = h->GetFunction("gaus"); - gaus->SetLineWidth(2); - gaus->SetLineColor(kRed); - double mean = gaus->GetParameter(1); - double sigma = gaus->GetParameter(2); - double mean_err = gaus->GetParError(1); - double sigma_err = gaus->GetParError(2); - save_canvas(c4, "fsam", E_label, particle_label); - return std::make_tuple(mean, sigma, mean_err, sigma_err); - } -} - -std::vector<std::string> read_scanned_energies(std::string input_energies_fname) -{ - std::vector<std::string> scanned_energies; - std::string E_label; - ifstream E_file (input_energies_fname); - if (E_file.is_open()) - { - while (E_file >> E_label) - { - scanned_energies.push_back(E_label); - } - E_file.close(); - return scanned_energies; - } - else - { - std::cout << "Unable to open file " << input_energies_fname << std::endl; - abort(); - } -} - -void hcal_barrel_energy_scan_analysis(std::string particle_label = "electron") -{ - // 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); - - auto scanned_energies = read_scanned_energies(fmt::format("sim_output/hcal_barrel_energy_scan_points_{}.txt", particle_label)); - - TGraphErrors gr_fsam(scanned_energies.size()-1); - TGraphErrors gr_fsam_res(scanned_energies.size()-1); - - for (const auto& E_label : scanned_energies) { - auto [fsam, fsam_res, fsam_err, fsam_res_err] = extract_sampling_fraction_parameters(particle_label, E_label); - auto E = std::stod(E_label); - - gr_fsam.SetPoint(gr_fsam.GetN(),E,100*fsam); - gr_fsam.SetPointError(gr_fsam.GetN()-1,0., 100*fsam_err); - gr_fsam_res.SetPoint(gr_fsam_res.GetN(),E,100.0*(fsam_res/fsam)); - auto fsam_res_rel_err = 100.0*(sqrt(pow((fsam_res_err/fsam),2)+pow((fsam_err*fsam_res)/(fsam*fsam),2))); - gr_fsam_res.SetPointError(gr_fsam_res.GetN()-1,0.,fsam_res_rel_err); - } - - TCanvas* c5 = new TCanvas("c5", "c5", 700, 500); - c5->cd(); - gr_fsam.SetTitle("Sampling Fraction Scan;True Energy [GeV];Sampling Fraction [%]"); - gr_fsam.SetMarkerStyle(20); - gr_fsam.Fit("pol0", "", "", 2., 20.); - gr_fsam.Draw("APE"); - save_canvas(c5, fmt::format("hcal_barrel_{}_fsam_scan", particle_label)); - - TCanvas* c6 = new TCanvas("c6", "c6", 700, 500); - c6->cd(); - TF1* func_res = new TF1("func_res", "[0]/sqrt(x) + [1]", 0.25, 20.); - func_res->SetLineWidth(2); - func_res->SetLineColor(kRed); - gr_fsam_res.SetTitle("Energy Resolution;True Energy [GeV];#Delta E/E [%]"); - gr_fsam_res.SetMarkerStyle(20); - gr_fsam_res.Fit(func_res,"R"); - gr_fsam_res.Draw("APE"); - save_canvas(c6,fmt::format("hcal_barrel_{}_fsam_scan_res", particle_label)); -} diff --git a/benchmarks/barrel_hcal/scripts/hcal_barrel_particles_analysis.cxx b/benchmarks/barrel_hcal/scripts/hcal_barrel_particles_analysis.cxx deleted file mode 100644 index 424802d33765b57ec223b745330f60a771ff2bc7..0000000000000000000000000000000000000000 --- a/benchmarks/barrel_hcal/scripts/hcal_barrel_particles_analysis.cxx +++ /dev/null @@ -1,149 +0,0 @@ -//////////////////////////////////////// -// Read reconstruction ROOT output file -// Plot variables -//////////////////////////////////////// - -#include "ROOT/RDataFrame.hxx" -#include <iostream> -#include <fstream> -#include <fmt/core.h> - -#include "edm4hep/MCParticleCollection.h" -#include "edm4hep/SimCalorimeterHitCollection.h" - -#include "TCanvas.h" -#include "TStyle.h" -#include "TMath.h" -#include "TH1.h" -#include "TF1.h" -#include "TH1D.h" - -#include "hcal_barrel_common_functions.h" - -using ROOT::RDataFrame; -using namespace ROOT::VecOps; - -void save_canvas(TCanvas* c, std::string label) -{ - c->SaveAs(fmt::format("results/{}.png",label).c_str()); - c->SaveAs(fmt::format("results/{}.pdf",label).c_str()); -} - -void save_canvas(TCanvas* c, std::string label, std::string particle_label) -{ - std::string label_with_E = fmt::format("hcal_barrel_{}_{}", particle_label, label); - save_canvas(c, label_with_E); -} - -void hcal_barrel_particles_analysis(std::string particle_name = "electron") -{ - // 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); - - ROOT::EnableImplicitMT(); - std::string input_fname = fmt::format("sim_output/sim_hcal_barrel_{}.edm4hep.root", particle_name); - ROOT::RDataFrame d0("events", input_fname); - - // Thrown Energy [GeV] - auto Ethr = [](std::vector<edm4hep::MCParticleData> const& input) { - auto p = input[2]; - auto energy = TMath::Sqrt(p.momentum.x * p.momentum.x + p.momentum.y * p.momentum.y + p.momentum.z * p.momentum.z + p.mass * p.mass); - return energy; - }; - - // Number of hits - auto nhits = [] (const std::vector<edm4hep::SimCalorimeterHitData>& evt) {return (int) evt.size(); }; - - // Energy deposition [GeV] - auto Esim = [](const std::vector<edm4hep::SimCalorimeterHitData>& evt) { - auto total_edep = 0.0; - for (const auto& i: evt) - total_edep += i.energy; - 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, {"HcalBarrelHits"}) - .Define("Esim", Esim, {"HcalBarrelHits"}) - .Define("fsam", fsam, {"Esim", "Ethr"}); - - // Define Histograms - auto hEthr = d1.Histo1D( - {"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 25.0}, - "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", 500, 0.0, 0.5}, - "Esim"); - auto hfsam = d1.Histo1D( - {"hfsam", "Sampling Fraction; Sampling Fraction; Events", 200, 0.0, 0.05}, - "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->SetLineWidth(2); - h->SetLineColor(kBlue); - save_canvas(c1,"Ethr",particle_name); - } - - { - TCanvas* c2 = new TCanvas("c2", "c2", 700, 500); - c2->SetLogy(1); - auto h = hNhits->DrawCopy(); - h->SetLineWidth(2); - h->SetLineColor(kBlue); - save_canvas(c2,"nhits",particle_name); - } - - { - TCanvas* c3 = new TCanvas("c3", "c3", 700, 500); - c3->SetLogy(1); - auto h = hEsim->DrawCopy(); - h->SetLineWidth(2); - h->SetLineColor(kBlue); - double up_fit = h->GetMean() + 5*h->GetStdDev(); - double down_fit = h->GetMean() - 5*h->GetStdDev(); - h->GetXaxis()->SetRangeUser(0.,up_fit); - save_canvas(c3,"Esim",particle_name); - } - - { - TCanvas* c4 = new TCanvas("c4", "c4", 700, 500); - c4->SetLogy(1); - auto h = hfsam->DrawCopy(); - h->SetLineWidth(2); - h->SetLineColor(kBlue); - double up_fit = h->GetMean() + 5*h->GetStdDev(); - double down_fit = h->GetMean() - 5*h->GetStdDev(); - h->Fit("gaus", "", "", down_fit, up_fit); - h->GetXaxis()->SetRangeUser(0.,up_fit); - TF1 *gaus = h->GetFunction("gaus"); - gaus->SetLineWidth(2); - gaus->SetLineColor(kRed); - save_canvas(c4,"fsam",particle_name); - } -} - diff --git a/benchmarks/barrel_hcal/scripts/hcal_barrel_particles_gen.cxx b/benchmarks/barrel_hcal/scripts/hcal_barrel_particles_gen.cxx deleted file mode 100644 index a35829db126657014b53a9656d077d80246e0359..0000000000000000000000000000000000000000 --- a/benchmarks/barrel_hcal/scripts/hcal_barrel_particles_gen.cxx +++ /dev/null @@ -1,90 +0,0 @@ -////////////////////////////////////////////////////////////// -// HCAL Barrel detector -// Generate particles with particle gun -// J.Kim 04/02/21 -// M.Zurek 05/05/21 -// W.Deconinck 05/27/21 -////////////////////////////////////////////////////////////// -#include "HepMC3/GenEvent.h" -#include "HepMC3/Print.h" -#include "HepMC3/ReaderAscii.h" -#include "HepMC3/WriterAscii.h" - -#include "TMath.h" -#include "TRandom.h" - -#include <cmath> -#include <iostream> -#include <fstream> -#include <math.h> -#include <random> -#include <cstdlib> -#include <fmt/core.h> - -#include "hcal_barrel_common_functions.h" - -using namespace HepMC3; - -void hcal_barrel_particles_gen(int n_events = 1e6, double e_start = 0.0, double e_end = 20.0, std::string particle_name = "muon") { - std::string out_fname; - auto env_fname = getenv("JUGGLER_GEN_FILE"); - if (env_fname != nullptr) - out_fname = env_fname; - else - out_fname = fmt::format("./data/hcal_barrel_{}.hepmc", particle_name); - WriterAscii hepmc_output(out_fname); - int events_parsed = 0; - GenEvent evt(Units::GEV, Units::MM); - - // Random number generator - TRandom* r1 = new TRandom(); - - // Constraining the solid angle, but larger than that subtended by the - // detector - // https://indico.bnl.gov/event/7449/contributions/35966/attachments/27177/41430/EIC-DWG-Calo-03192020.pdf - // See a figure on slide 26 - double cos_theta_min = std::cos(M_PI * (45.0 / 180.0)); - double cos_theta_max = std::cos(M_PI * (135.0 / 180.0)); - - for (events_parsed = 0; events_parsed < n_events; events_parsed++) { - // FourVector(px,py,pz,e,pdgid,status) - // type 4 is beam - // pdgid 11 - electron - // pdgid 2212 - proton - GenParticlePtr p1 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4); - GenParticlePtr p2 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4); - - // Define momentum - Double_t p = r1->Uniform(e_start, e_end); - Double_t phi = r1->Uniform(0.0, 2.0 * M_PI); - Double_t costheta = r1->Uniform(cos_theta_min, cos_theta_max); - Double_t theta = std::acos(costheta); - Double_t px = p * std::cos(phi) * std::sin(theta); - Double_t py = p * std::sin(phi) * std::sin(theta); - Double_t pz = p * std::cos(theta); - - // type 1 is final state - auto [id, mass] = extract_particle_parameters(particle_name); - GenParticlePtr p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (mass * mass))), id, 1); - - GenVertexPtr v1 = std::make_shared<GenVertex>(); - v1->add_particle_in(p1); - v1->add_particle_in(p2); - v1->add_particle_out(p3); - evt.add_vertex(v1); - - if (events_parsed == 0) { - std::cout << "First event: " << std::endl; - Print::listing(evt); - } - - hepmc_output.write_event(evt); - if (events_parsed % 10000 == 0) { - std::cout << "Event: " << events_parsed << std::endl; - } - evt.clear(); - } - hepmc_output.close(); - std::cout << "Events parsed and written: " << events_parsed << std::endl; -} - diff --git a/benchmarks/barrel_hcal/scripts/hcal_barrel_particles_reader.cxx b/benchmarks/barrel_hcal/scripts/hcal_barrel_particles_reader.cxx deleted file mode 100644 index f40f7c29c77c55165ecfca8350ee08924974370b..0000000000000000000000000000000000000000 --- a/benchmarks/barrel_hcal/scripts/hcal_barrel_particles_reader.cxx +++ /dev/null @@ -1,152 +0,0 @@ -////////////////////////// -// HCAL Barrel detector -// J.KIM 04/02/2021 -// M.Zurek 05/05/2021 -// W.Deconinck 05/28/21 -////////////////////////// -#include "HepMC3/GenEvent.h" -#include "HepMC3/Print.h" -#include "HepMC3/ReaderAscii.h" -#include "HepMC3/WriterAscii.h" - -#include "TROOT.h" -#include "TH1.h" -#include "TH2.h" -#include "TH3.h" -#include "TStyle.h" -#include "TCanvas.h" -#include "TMath.h" - -#include <iostream> -#include <fstream> -#include <cstdlib> -#include <fmt/core.h> - -#include "hcal_barrel_common_functions.h" - -using namespace HepMC3; - -void save_canvas(TCanvas* c, std::string label) -{ - c->SaveAs(fmt::format("results/{}.png",label).c_str()); - c->SaveAs(fmt::format("results/{}.pdf",label).c_str()); -} - -void save_canvas(TCanvas* c, std::string label, std::string particle_label) -{ - std::string label_with_E = fmt::format("input_hcal_barrel_{}_{}", particle_label, label); - save_canvas(c, label_with_E); -} - -void hcal_barrel_particles_reader(std::string particle_name = "electron") -{ - // Setting for graphs - gROOT->SetStyle("Plain"); - gStyle->SetOptFit(1); - gStyle->SetLineWidth(1); - gStyle->SetPadTickX(1); - gStyle->SetPadTickY(1); - gStyle->SetPadGridX(1); - gStyle->SetPadGridY(1); - gStyle->SetPadLeftMargin(0.14); - gStyle->SetPadRightMargin(0.17); - - std::string in_fname; - auto env_fname = getenv("JUGGLER_GEN_FILE"); - if (env_fname != nullptr) - in_fname = env_fname; - else - in_fname = fmt::format("./data/hcal_barrel_{}.hepmc", particle_name); - ReaderAscii hepmc_input(in_fname); - int events_parsed = 0; - GenEvent evt(Units::GEV, Units::MM); - - // Histograms - TH1F* h_energy = new TH1F(fmt::format("h_{}_energy",particle_name).c_str(), fmt::format("{} energy;E [GeV];Events",particle_name).c_str(), 100, -0.5, 30.5); - TH1F* h_eta = new TH1F(fmt::format("h_{}_eta",particle_name).c_str(), fmt::format("{} #eta;#eta;Events",particle_name).c_str(), 100, -10.0, 10.0); - TH1F* h_theta = new TH1F(fmt::format("h_{}_theta",particle_name).c_str(), fmt::format("{} #theta;#theta [degree];Events",particle_name).c_str(), 100, -0.5, 180.5); - TH1F* h_phi = new TH1F(fmt::format("h_{}_phi",particle_name).c_str(), fmt::format("{} #phi;#phi [degree];Events",particle_name).c_str(), 100, -180.5, 180.5); - TH2F* h_pzpt = new TH2F(fmt::format("h_{}_pzpt",particle_name).c_str(), fmt::format("{} pt vs pz;pt [GeV];pz [GeV]",particle_name).c_str(), 100, -0.5, 30.5, 100, -30.5, 30.5); - TH2F* h_pxpy = new TH2F(fmt::format("h_{}_pxpy",particle_name).c_str(), fmt::format("{} px vs py;px [GeV];py [GeV]",particle_name).c_str(), 100, -30.5, 30.5, 100, -30.5, 30.5); - TH3F* h_p = new TH3F(fmt::format("h_{}_p",particle_name).c_str(), fmt::format("{} p;px [GeV];py [GeV];pz [GeV]",particle_name).c_str(), 100, -30.5, 30.5, 100, -30.5, 30.5, 100, -30.5, 30.5); - - while (!hepmc_input.failed()) { - // Read event from input file - hepmc_input.read_event(evt); - // If reading failed - exit loop - if (hepmc_input.failed()) - break; - - auto [id, mass] = extract_particle_parameters(particle_name); - - for (const auto& v : evt.vertices()) { - for (const auto& p : v->particles_out()) { - if (p->pid() == id) { - h_energy->Fill(p->momentum().e()); - h_eta->Fill(p->momentum().eta()); - h_theta->Fill(p->momentum().theta() * TMath::RadToDeg()); - h_phi->Fill(p->momentum().phi() * TMath::RadToDeg()); - h_pzpt->Fill(TMath::Sqrt(p->momentum().px() * p->momentum().px() + p->momentum().py() * p->momentum().py()), p->momentum().pz()); - h_pxpy->Fill(p->momentum().px(), p->momentum().py()); - h_p->Fill(p->momentum().px(), p->momentum().py(), p->momentum().pz()); - } - } - } - evt.clear(); - events_parsed++; - } - std::cout << "Events parsed and written: " << events_parsed << std::endl; - - TCanvas* c = new TCanvas("c", "c", 500, 500); - h_energy->GetYaxis()->SetTitleOffset(1.8); - h_energy->SetLineWidth(2); - h_energy->SetLineColor(kBlue); - h_energy->DrawClone(); - save_canvas(c, "energy", particle_name); - - TCanvas* c1 = new TCanvas("c1", "c1", 500, 500); - h_eta->GetYaxis()->SetTitleOffset(1.9); - h_eta->SetLineWidth(2); - h_eta->SetLineColor(kBlue); - h_eta->DrawClone(); - save_canvas(c1, "eta", particle_name); - - TCanvas* c2 = new TCanvas("c2", "c2", 500, 500); - h_theta->GetYaxis()->SetTitleOffset(1.8); - h_theta->SetLineWidth(2); - h_theta->SetLineColor(kBlue); - h_theta->DrawClone(); - save_canvas(c2, "theta", particle_name); - - TCanvas* c3 = new TCanvas("c3", "c3", 500, 500); - h_phi->GetYaxis()->SetTitleOffset(1.8); - h_phi->SetLineWidth(2); - h_phi->GetYaxis()->SetRangeUser(0.0, h_phi->GetMaximum() + 100.0); - h_phi->SetLineColor(kBlue); - h_phi->DrawClone(); - save_canvas(c3, "phi", particle_name); - - TCanvas* c4 = new TCanvas("c4", "c4", 500, 500); - h_pzpt->GetYaxis()->SetTitleOffset(1.4); - h_pzpt->SetLineWidth(2); - h_pzpt->SetLineColor(kBlue); - h_pzpt->DrawClone("COLZ"); - save_canvas(c4, "pzpt", particle_name); - - TCanvas* c5 = new TCanvas("c5", "c5", 500, 500); - h_pxpy->GetYaxis()->SetTitleOffset(1.4); - h_pxpy->SetLineWidth(2); - h_pxpy->SetLineColor(kBlue); - h_pxpy->DrawClone("COLZ"); - save_canvas(c5, "pxpy", particle_name); - - TCanvas* c6 = new TCanvas("c6", "c6", 500, 500); - h_p->GetYaxis()->SetTitleOffset(1.8); - h_p->GetXaxis()->SetTitleOffset(1.6); - h_p->GetZaxis()->SetTitleOffset(1.6); - h_p->SetLineWidth(2); - h_p->SetLineColor(kBlue); - h_p->DrawClone(); - save_canvas(c6, "p", particle_name); -} -