diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index abb866c7..f6324f6d 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -124,74 +124,80 @@ common:setup: - runner_system_failure include: - - local: 'benchmarks/backgrounds/config.yml' - - local: 'benchmarks/backwards_ecal/config.yml' - - local: 'benchmarks/beamline/config.yml' - - local: 'benchmarks/calo_pid/config.yml' - - local: 'benchmarks/campaign/config.yml' - - local: 'benchmarks/ecal_gaps/config.yml' - - local: 'benchmarks/far_forward_dvcs/config.yml' - - local: 'benchmarks/lowq2_reconstruction/config.yml' - - local: 'benchmarks/tracking_detectors/config.yml' - - local: 'benchmarks/tracking_performances/config.yml' - - local: 'benchmarks/tracking_performances_dis/config.yml' - - local: 'benchmarks/barrel_ecal/config.yml' - - local: 'benchmarks/lfhcal/config.yml' - - local: 'benchmarks/zdc/config.yml' - - local: 'benchmarks/zdc_lyso/config.yml' - - local: 'benchmarks/zdc_neutron/config.yml' - - local: 'benchmarks/zdc_photon/config.yml' - - local: 'benchmarks/zdc_pi0/config.yml' - - local: 'benchmarks/material_scan/config.yml' - - local: 'benchmarks/pid/config.yml' - - local: 'benchmarks/rich/config.yml' - - local: 'benchmarks/b0_tracker/config.yml' - - local: 'benchmarks/insert_muon/config.yml' - - local: 'benchmarks/insert_tau/config.yml' - - local: 'benchmarks/zdc_sigma/config.yml' - - local: 'benchmarks/zdc_lambda/config.yml' - - local: 'benchmarks/insert_neutron/config.yml' - - local: 'benchmarks/femc_electron/config.yml' - - local: 'benchmarks/femc_photon/config.yml' - - local: 'benchmarks/femc_pi0/config.yml' - - local: 'benchmarks/nhcal_acceptance/config.yml' + # - local: 'benchmarks/backgrounds/config.yml' + # - local: 'benchmarks/backwards_ecal/config.yml' + # - local: 'benchmarks/beamline/config.yml' + # - local: 'benchmarks/calo_pid/config.yml' + # - local: 'benchmarks/campaign/config.yml' + # - local: 'benchmarks/ecal_gaps/config.yml' + # - local: 'benchmarks/far_forward_dvcs/config.yml' + # - local: 'benchmarks/lowq2_reconstruction/config.yml' + # - local: 'benchmarks/tracking_detectors/config.yml' + # - local: 'benchmarks/tracking_performances/config.yml' + # - local: 'benchmarks/tracking_performances_dis/config.yml' + # - local: 'benchmarks/barrel_ecal/config.yml' + # - local: 'benchmarks/lfhcal/config.yml' + # - local: 'benchmarks/zdc/config.yml' + # - local: 'benchmarks/zdc_lyso/config.yml' + # - local: 'benchmarks/zdc_neutron/config.yml' + # - local: 'benchmarks/zdc_photon/config.yml' + # - local: 'benchmarks/zdc_pi0/config.yml' + # - local: 'benchmarks/material_scan/config.yml' + # - local: 'benchmarks/pid/config.yml' + # - local: 'benchmarks/rich/config.yml' + # - local: 'benchmarks/b0_tracker/config.yml' + # - local: 'benchmarks/insert_muon/config.yml' + # - local: 'benchmarks/insert_tau/config.yml' + # - local: 'benchmarks/zdc_sigma/config.yml' + # - local: 'benchmarks/zdc_lambda/config.yml' + # - local: 'benchmarks/insert_neutron/config.yml' + # - local: 'benchmarks/femc_electron/config.yml' + # - local: 'benchmarks/femc_photon/config.yml' + # - local: 'benchmarks/femc_pi0/config.yml' + # - local: 'benchmarks/nhcal_acceptance/config.yml' - local: 'benchmarks/nhcal_basic_distribution/config.yml' + # - local: 'benchmarks/nhcal_sampling_fraction/config.yml' + - local: 'benchmarks/nhcal_dimuon_photoproduction/config.yml' + - local: 'benchmarks/nhcal_pion_rejection/config.yml' deploy_results: allow_failure: true stage: deploy needs: - - "collect_results:backgrounds" - - "collect_results:backwards_ecal" - - "collect_results:barrel_ecal" - - "collect_results:beamline" - - "collect_results:calo_pid" - - "collect_results:campaign" - - "collect_results:ecal_gaps" - - "collect_results:far_forward_dvcs" - - "collect_results:lfhcal" - - "collect_results:lowq2_reconstruction" - - "collect_results:material_scan" - - "collect_results:pid" - - "collect_results:rich" - - "collect_results:tracking_performance" - - "collect_results:tracking_performance_campaigns" - - "collect_results:zdc_sigma" - - "collect_results:zdc_lambda" - - "collect_results:insert_neutron" - - "collect_results:tracking_performances_dis" - - "collect_results:zdc" - - "collect_results:zdc_lyso" - - "collect_results:zdc_neutron" - - "collect_results:insert_muon" - - "collect_results:insert_tau" - - "collect_results:zdc_photon" - - "collect_results:zdc_pi0" - - "collect_results:femc_electron" - - "collect_results:femc_photon" - - "collect_results:femc_pi0" - - "collect_results:nhcal_acceptance" + # - "collect_results:backgrounds" + # - "collect_results:backwards_ecal" + # - "collect_results:barrel_ecal" + # - "collect_results:beamline" + # - "collect_results:calo_pid" + # - "collect_results:campaign" + # - "collect_results:ecal_gaps" + # - "collect_results:far_forward_dvcs" + # - "collect_results:lfhcal" + # - "collect_results:lowq2_reconstruction" + # - "collect_results:material_scan" + # - "collect_results:pid" + # - "collect_results:rich" + # - "collect_results:tracking_performance" + # - "collect_results:tracking_performance_campaigns" + # - "collect_results:zdc_sigma" + # - "collect_results:zdc_lambda" + # - "collect_results:insert_neutron" + # - "collect_results:tracking_performances_dis" + # - "collect_results:zdc" + # - "collect_results:zdc_lyso" + # - "collect_results:zdc_neutron" + # - "collect_results:insert_muon" + # - "collect_results:insert_tau" + # - "collect_results:zdc_photon" + # - "collect_results:zdc_pi0" + # - "collect_results:femc_electron" + # - "collect_results:femc_photon" + # - "collect_results:femc_pi0" + # - "collect_results:nhcal_acceptance" - "collect_results:nhcal_basic_distribution" + # - "collect_results:nhcal_sampling_fraction" + - "collect_results:nhcal_dimuon_photoproduction" + - "collect_results:nhcal_pion_rejection" script: - snakemake $SNAKEMAKE_FLAGS --cores 1 results/metadata.json - find results -print | sort | tee summary.txt diff --git a/Snakefile b/Snakefile index fbf6f5e0..a49a1146 100644 --- a/Snakefile +++ b/Snakefile @@ -43,33 +43,36 @@ def find_epic_libraries(): return libs -include: "benchmarks/backgrounds/Snakefile" -include: "benchmarks/backwards_ecal/Snakefile" -include: "benchmarks/barrel_ecal/Snakefile" -include: "benchmarks/beamline/Snakefile" -include: "benchmarks/calo_pid/Snakefile" -include: "benchmarks/campaign/Snakefile" -include: "benchmarks/ecal_gaps/Snakefile" -include: "benchmarks/far_forward_dvcs/Snakefile" -include: "benchmarks/lowq2_reconstruction/Snakefile" -include: "benchmarks/material_scan/Snakefile" -include: "benchmarks/tracking_performances/Snakefile" -include: "benchmarks/tracking_performances_dis/Snakefile" -include: "benchmarks/lfhcal/Snakefile" -include: "benchmarks/zdc_lyso/Snakefile" -include: "benchmarks/zdc_neutron/Snakefile" -include: "benchmarks/insert_muon/Snakefile" -include: "benchmarks/zdc_lambda/Snakefile" -include: "benchmarks/zdc_photon/Snakefile" -include: "benchmarks/zdc_pi0/Snakefile" -include: "benchmarks/zdc_sigma/Snakefile" -include: "benchmarks/insert_neutron/Snakefile" -include: "benchmarks/insert_tau/Snakefile" -include: "benchmarks/femc_electron/Snakefile" -include: "benchmarks/femc_photon/Snakefile" -include: "benchmarks/femc_pi0/Snakefile" -include: "benchmarks/nhcal_acceptance/Snakefile" +# include: "benchmarks/backgrounds/Snakefile" +# include: "benchmarks/backwards_ecal/Snakefile" +# include: "benchmarks/barrel_ecal/Snakefile" +# include: "benchmarks/beamline/Snakefile" +# include: "benchmarks/calo_pid/Snakefile" +# include: "benchmarks/campaign/Snakefile" +# include: "benchmarks/ecal_gaps/Snakefile" +# include: "benchmarks/far_forward_dvcs/Snakefile" +# include: "benchmarks/lowq2_reconstruction/Snakefile" +# include: "benchmarks/material_scan/Snakefile" +# include: "benchmarks/tracking_performances/Snakefile" +# include: "benchmarks/tracking_performances_dis/Snakefile" +# include: "benchmarks/lfhcal/Snakefile" +# include: "benchmarks/zdc_lyso/Snakefile" +# include: "benchmarks/zdc_neutron/Snakefile" +# include: "benchmarks/insert_muon/Snakefile" +# include: "benchmarks/zdc_lambda/Snakefile" +# include: "benchmarks/zdc_photon/Snakefile" +# include: "benchmarks/zdc_pi0/Snakefile" +# include: "benchmarks/zdc_sigma/Snakefile" +# include: "benchmarks/insert_neutron/Snakefile" +# include: "benchmarks/insert_tau/Snakefile" +# include: "benchmarks/femc_electron/Snakefile" +# include: "benchmarks/femc_photon/Snakefile" +# include: "benchmarks/femc_pi0/Snakefile" +# include: "benchmarks/nhcal_acceptance/Snakefile" include: "benchmarks/nhcal_basic_distribution/Snakefile" +# include: "benchmarks/nhcal_sampling_fraction/Snakefile" +include: "benchmarks/nhcal_dimuon_photoproduction/Snakefile" +include: "benchmarks/nhcal_pion_rejection/Snakefile" use_s3 = config["remote_provider"].lower() == "s3" use_xrootd = config["remote_provider"].lower() == "xrootd" diff --git a/benchmarks/nhcal_acceptance/scripts/acceptance_analysis.cxx b/benchmarks/nhcal_acceptance/scripts/acceptance_analysis.cxx index d8248894..7ef3ca16 100644 --- a/benchmarks/nhcal_acceptance/scripts/acceptance_analysis.cxx +++ b/benchmarks/nhcal_acceptance/scripts/acceptance_analysis.cxx @@ -32,10 +32,13 @@ int acceptance_analysis(TString filename, string outname_pdf, string outname_png int nPhiBins = 100; double etaMin = -5, etaMax = 0; - TH2D* hEtaPhiAll = new TH2D("hEtaPhiAll", "All #pi- (status==1); #eta[1]; #phi[rad]", + gStyle->SetTitleSize(0.05, "XYZ"); + gStyle->SetLabelSize(0.05, "XYZ"); + + TH2D* hEtaPhiAll = new TH2D("hEtaPhiAll", "All #pi- (status==1); #eta; #phi[rad]", nEtaBins, etaMin, etaMax, nPhiBins, -TMath::Pi(), TMath::Pi()); - TH2D* hEtaPhiDetected = new TH2D("hEtaPhiDetected", "#pi- detected in nHCal; #eta[1]; #phi[rad]", + TH2D* hEtaPhiDetected = new TH2D("hEtaPhiDetected", "#pi- detected in nHCal; #eta; #phi[rad]", nEtaBins, etaMin, etaMax, nPhiBins, -TMath::Pi(), TMath::Pi()); while (reader.Next()) diff --git a/benchmarks/nhcal_basic_distribution/Snakefile b/benchmarks/nhcal_basic_distribution/Snakefile index 9e672a71..ecca3d20 100644 --- a/benchmarks/nhcal_basic_distribution/Snakefile +++ b/benchmarks/nhcal_basic_distribution/Snakefile @@ -1,3 +1,8 @@ +def get_total_energy_neutron(kinetic_energy): + ke = float(kinetic_energy) + neutron_mass = 0.94 + return f"{ke + neutron_mass:.2f}" + rule nhcal_basic_distribution_simulate: input: warmup="warmup.edm4hep.root", @@ -9,9 +14,9 @@ rule nhcal_basic_distribution_simulate: SEED=lambda wildcards: "1" + wildcards.INDEX, DETECTOR_PATH=os.environ["DETECTOR_PATH"], DETECTOR_CONFIG=lambda wildcards: wildcards.DETECTOR_CONFIG, - ENERGY=lambda wildcards: wildcards.ENERGY, DD4HEP_HASH=get_spack_package_hash("dd4hep"), NPSIM_HASH=get_spack_package_hash("npsim"), + TOTAL_ENERGY=lambda wildcards: get_total_energy_neutron(wildcards.ENERGY), cache: True shell: """ @@ -25,15 +30,14 @@ exec npsim \ --gun.thetaMin 120*degree \ --gun.thetaMax 180*degree \ --gun.distribution uniform \ - --gun.energy "{params.ENERGY}*GeV" \ + --gun.energy "{params.TOTAL_ENERGY}*GeV" \ --outputFile {output} """ - rule nhcal_basic_distribution_combine: input: lambda wildcards: expand( - "sim_output/nhcal_basic_distribution/E{ENERGY:.1f}GeV/sim_{DETECTOR_CONFIG}.{INDEX:02d}.edm4hep.root", + "sim_output/nhcal_basic_distribution/E{ENERGY:.1f}GeV/sim_{DETECTOR_CONFIG}.{INDEX}.edm4hep.root", DETECTOR_CONFIG=wildcards.DETECTOR_CONFIG, ENERGY=float(wildcards.ENERGY), INDEX=range(int(wildcards.N)), @@ -42,20 +46,52 @@ rule nhcal_basic_distribution_combine: N=r"\d+", ENERGY=r"\d+(\.\d+)?" output: - temp("sim_output/nhcal_basic_distribution/sim_{DETECTOR_CONFIG}_E{ENERGY}GeV_combined_{N}files.edm4hep.root"), + temp("sim_output/nhcal_basic_distribution/sim_E{ENERGY}GeV_{DETECTOR_CONFIG}_{N}files.edm4hep.root"), shell: """ -hadd -f {output} {input} -""" + hadd -f {output} {input} + """ rule nhcal_basic_distribution_analysis: input: - combined="sim_output/nhcal_basic_distribution/sim_{DETECTOR_CONFIG}_E{ENERGY}GeV_combined_{N}files.edm4hep.root", + combined="sim_output/nhcal_basic_distribution/sim_E{ENERGY}GeV_{DETECTOR_CONFIG}_{N}files.edm4hep.root", script="benchmarks/nhcal_basic_distribution/scripts/basic_distribution_analysis.cxx", output: - pdf=f"results/nhcal_basic_distribution/analysis_{{DETECTOR_CONFIG}}_E{{ENERGY}}GeV_combined_{{N}}files.pdf", - png=f"results/nhcal_basic_distribution/analysis_{{DETECTOR_CONFIG}}_E{{ENERGY}}GeV_combined_{{N}}files.png", + png="results/nhcal_basic_distribution/analysis_E{ENERGY}GeV_{DETECTOR_CONFIG}_{N}files.png", + root="sim_output/nhcal_basic_distribution/energy_resolution_E{ENERGY}GeV_{DETECTOR_CONFIG}_{N}files.root", + wildcard_constraints: + ENERGY=r"\d+(\.\d+)?", + N=r"\d+", + params: + DETECTOR_PATH = os.environ["DETECTOR_PATH"], + DETECTOR_CONFIG = lambda wildcards: f"{wildcards.DETECTOR_CONFIG}.xml", shell: """ - root -l -b -q '{input.script}+("{input.combined}","{output.pdf}","{output.png}")' -""" + root -l -b -q '{input.script}+("{input.combined}","{output.png}","{output.root}","{params.DETECTOR_PATH}/{params.DETECTOR_CONFIG}")' + """ + +rule nhcal_basic_distribution_combine_energy_resolution: + input: + expand( + "sim_output/nhcal_basic_distribution/energy_resolution_E{ENERGY}GeV_{{DETECTOR_CONFIG}}_{N}files.root", + ENERGY=["0.5", "0.7", "1.0", "2.0", "5.0"], + N=["10"] + ) + output: + temp("sim_output/nhcal_basic_distribution/energy_resolution_combined_{DETECTOR_CONFIG}.root"), + shell: + """ + hadd -f {output} {input} + """ + +rule nhcal_basic_distribution_analysis_energy_resolution: + input: + root="sim_output/nhcal_basic_distribution/energy_resolution_combined_{DETECTOR_CONFIG}.root", + script="benchmarks/nhcal_basic_distribution/scripts/basic_distribution_energy_resolution.cxx", + output: + png="results/nhcal_basic_distribution/energy_resolution_analysis_{DETECTOR_CONFIG}.png", + shell: + """ + root -l -b -q '{input.script}("{input.root}","{output.png}")' + """ + diff --git a/benchmarks/nhcal_basic_distribution/config.yml b/benchmarks/nhcal_basic_distribution/config.yml index b3f7da69..4ad6c964 100644 --- a/benchmarks/nhcal_basic_distribution/config.yml +++ b/benchmarks/nhcal_basic_distribution/config.yml @@ -4,60 +4,69 @@ sim:nhcal_basic_distribution: timeout: 4 hours parallel: matrix: - - ENERGY: ["0.5GeV", "0.7GeV", "1.0GeV", "2.0GeV", "5.0GeV"] - INDEX_RANGE: ["0 4","5 9"] + - ENERGY: ["0.5", "0.7", "1.0", "2.0", "5.0"] + INDEX_RANGE: ["0 3","4 6", "7 9"] script: - - snakemake $SNAKEMAKE_FLAGS --cores $MAX_CORES_PER_JOB $(for INDEX in $(seq -f '%02.0f' $INDEX_RANGE); do echo sim_output/nhcal_basic_distribution/E${ENERGY}/sim_epic_backward_hcal_only.${INDEX}.edm4hep.root; done) + - snakemake $SNAKEMAKE_FLAGS --cores $MAX_CORES_PER_JOB $(for INDEX in $(seq $INDEX_RANGE); do echo sim_output/nhcal_basic_distribution/E${ENERGY}GeV/sim_epic_backward_hcal_only.${INDEX}.edm4hep.root; done) -sim:nhcal_basic_distribution_full: - extends: .det_benchmark - stage: simulate - timeout: 4 hours - parallel: - matrix: - - ENERGY: ["0.5GeV", "0.7GeV", "1.0GeV", "2.0GeV"] - INDEX_RANGE: ["0 4","5 9"] - - ENERGY: ["5.0GeV"] - INDEX_RANGE: ["0 3", "4 6", "7 9"] - script: - - snakemake $SNAKEMAKE_FLAGS --cores $MAX_CORES_PER_JOB $(for INDEX in $(seq -f '%02.0f' $INDEX_RANGE); do echo sim_output/nhcal_basic_distribution/E${ENERGY}/sim_epic_full.${INDEX}.edm4hep.root; done) +# sim:nhcal_basic_distribution_full: +# extends: .det_benchmark +# stage: simulate +# timeout: 4 hours +# parallel: +# matrix: +# - ENERGY: ["0.5", "0.7", "1.0", "2.0"] +# INDEX_RANGE: ["0 4","5 9"] +# - ENERGY: ["5.0"] +# INDEX_RANGE: ["0 2", "3 5", "6 7", "8 9"] + +# script: +# - snakemake $SNAKEMAKE_FLAGS --cores $MAX_CORES_PER_JOB $(for INDEX in $(seq $INDEX_RANGE); do echo sim_output/nhcal_basic_distribution/E${ENERGY}GeV/sim_epic_full.${INDEX}.edm4hep.root; done) bench:nhcal_basic_distribution_analysis: extends: .det_benchmark stage: benchmarks needs: - "sim:nhcal_basic_distribution" - parallel: - matrix: - - ENERGY: ["0.5GeV", "0.7GeV", "1.0GeV", "2.0GeV", "5.0GeV"] - script: - - snakemake $SNAKEMAKE_FLAGS --cores 1 results/nhcal_basic_distribution/analysis_epic_backward_hcal_only_E${ENERGY}_combined_10files.pdf - -bench:nhcal_basic_distribution_analysis_full: - extends: .det_benchmark - stage: benchmarks - needs: - - "sim:nhcal_basic_distribution_full" - parallel: - matrix: - - ENERGY: ["0.5GeV", "0.7GeV", "1.0GeV", "2.0GeV", "5.0GeV"] + variables: + ENERGIES: "0.5 0.7 1.0 2.0 5.0" script: - - snakemake $SNAKEMAKE_FLAGS --cores 1 results/nhcal_basic_distribution/analysis_epic_full_E${ENERGY}_combined_10files.pdf + - | + for ENERGY in $ENERGIES; do + snakemake $SNAKEMAKE_FLAGS --cores 1 \ + results/nhcal_basic_distribution/analysis_E${ENERGY}GeV_epic_backward_hcal_only_10files.png + done + - snakemake $SNAKEMAKE_FLAGS --cores 1 \ + results/nhcal_basic_distribution/energy_resolution_analysis_epic_backward_hcal_only.png +# bench:nhcal_basic_distribution_analysis_full: +# extends: .det_benchmark +# stage: benchmarks +# needs: +# - "sim:nhcal_basic_distribution_full" +# variables: +# ENERGIES: "0.5 0.7 1.0 2.0 5.0" +# script: +# - | +# for ENERGY in $ENERGIES; do +# snakemake $SNAKEMAKE_FLAGS --cores 1 \ +# results/nhcal_basic_distribution/analysis_E${ENERGY}GeV_epic_full_10files.png +# done +# - snakemake $SNAKEMAKE_FLAGS --cores 1 \ +# results/nhcal_basic_distribution/energy_resolution_analysis_epic_full.png collect_results:nhcal_basic_distribution: extends: .det_benchmark stage: collect needs: - "bench:nhcal_basic_distribution_analysis" - - "bench:nhcal_basic_distribution_analysis_full" + # - "bench:nhcal_basic_distribution_analysis_full" parallel: matrix: - - ENERGY: ["0.5GeV", "0.7GeV", "1.0GeV", "2.0GeV", "5.0GeV"] - DETECTOR_CONFIG: ["epic_backward_hcal_only", "epic_full"] + - DETECTOR_CONFIG: ["epic_backward_hcal_only"] script: - ls -lrht - mv results{,_save}/ - - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/nhcal_basic_distribution/analysis_${DETECTOR_CONFIG}_E${ENERGY}_combined_10files.pdf + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/nhcal_basic_distribution/energy_resolution_analysis_{DETECTOR_CONFIG}.png - mv results{_save,}/ diff --git a/benchmarks/nhcal_basic_distribution/scripts/HistogramsNeutronThresholds.h b/benchmarks/nhcal_basic_distribution/scripts/HistogramsNeutronThresholds.h new file mode 100644 index 00000000..2f4ce605 --- /dev/null +++ b/benchmarks/nhcal_basic_distribution/scripts/HistogramsNeutronThresholds.h @@ -0,0 +1,59 @@ +/* + * HistogramsSim.h + * + * Created on: 21 apr 2025 + * Author: Sam Corey + */ + +#ifndef HISTOGRAMSNEUTRONTHRESHOLDS_H_ +#define HISTOGRAMSNEUTRONTHRESHOLDS_H_ + +#include +#include + +#include +#include + +#include + +using namespace TMath; + +void CreateHistogamsNeutronThresholds(); + +// HcalEndcapNHit + +TH1D *h_nHCal_hit_contrib_energy; +TH1D *h_nHCal_hit_contrib_time; + +TH2D *h_nHCal_hit_contrib_energy_vs_time; +TH2D *h_nHCal_hit_contrib_energy_vs_telap; +TH2D *h_nHCal_hit_contrib_energy_vs_time_total; + + +void CreateHistogamsNeutronThresholds() +{ + + // HcalEndcapNHit + h_nHCal_hit_contrib_energy = new TH1D("h_nHCal_hit_contrib_time", "Backwards HCal Hit Contribution Time; t; counts", 50, 0, 1000); + h_nHCal_hit_contrib_time = new TH1D("h_nHCal_hit_contrib_energy", "Backwards HCal Hit Contribution Energy; Energy [GeV]; counts", 50, 0, 0.00000002); + + h_nHCal_hit_contrib_energy_vs_time = new TH2D("h_nHCal_hit_contrib_2D_E_vs_t", "Backwards HCal Hit Contribution Energy vs. time; Energy [GeV]; time [ns]; counts", 50, 0, 0.00000002, 50, 0, 1000); + h_nHCal_hit_contrib_energy_vs_telap = new TH2D("h_nHCal_hit_contrib_energy_vs_telap", "Backwards HCal Hit Contribution Energy vs. Time; Energy threshold [GeV]; t max; counts", 50, 0, 0.0005, 50, 0, 1100); + h_nHCal_hit_contrib_energy_vs_time_total = new TH2D("h_nHCal_hit_contrib_energy_vs_time_total", "Backwards HCal Hit Contribution Energy vs. Time; Energy threshold [GeV]; t max; counts", 50, 0, 0.0005, 50, 0, 1100); + +} + +void DeleteHistogamsNeutronThresholds() +{ + + delete h_nHCal_hit_contrib_energy; + delete h_nHCal_hit_contrib_time; + + delete h_nHCal_hit_contrib_energy_vs_time; + delete h_nHCal_hit_contrib_energy_vs_telap; + delete h_nHCal_hit_contrib_energy_vs_time_total; + +} + + +#endif /* HISTOGRAMSDEPTHCHECK_H_ */ diff --git a/benchmarks/nhcal_basic_distribution/scripts/NeutronThresholdUtil.h b/benchmarks/nhcal_basic_distribution/scripts/NeutronThresholdUtil.h new file mode 100644 index 00000000..e22277be --- /dev/null +++ b/benchmarks/nhcal_basic_distribution/scripts/NeutronThresholdUtil.h @@ -0,0 +1,93 @@ +// NeutronThresholdsUtil.h +#ifndef NEUTRON_THRESHOLDS_UTIL_H +#define NEUTRON_THRESHOLDS_UTIL_H + +#include +#include +#include + + +// Constants for the neutron threshold analysis +namespace NeutronThresholds { + const double E_MIP = 0.00075; + + // Time thresholds in ns + const std::vector T_MAX = {25, 50, 100, 200, 500}; + + // Energy thresholds based on MIP + const std::vector E_TH = {0.5 * E_MIP, 0.25 * E_MIP, 0.1 * E_MIP, 0.05 * E_MIP, 0}; + + // Initialize the hits passed counters + inline std::vector> createHitsPassedMatrix() { + return std::vector>(5, std::vector(5, 0)); + } + + // Process the contributions for a hit + template + void processContributions(const ContribCollection& contrib, std::vector>& hits_passed, std::vector>& hits_passed_telap, + TH1* h_contrib_time = nullptr, TH1* h_contrib_energy = nullptr, double E_TH_correction = 1.0, bool debug = false + ) { + std::vector E_sum(5, 0); + std::vector E_sum_telap(5, 0); + double t_min = 1100; + + for (unsigned c = 0; c < contrib.size(); ++c) { + if (contrib.at(c).getTime() < t_min) { t_min = contrib.at(c).getTime(); } + } + + if(debug) std::cout << contrib.at(0).getTime() << std::endl; + + for (unsigned c = 0; c < contrib.size(); ++c) { + if(debug) std::cout << "hit time = " << contrib.at(c).getTime() << std::endl; + + // double t_min_per_contrib = std::min(1100, contrib.at(c).getTime()); + + if (h_contrib_time) h_contrib_time->Fill(contrib.at(c).getTime()); + if (h_contrib_energy) h_contrib_energy->Fill(contrib.at(c).getEnergy()); + + for (int itm = 0; itm < T_MAX.size(); itm++) { + if (contrib.at(c).getTime() < T_MAX[itm]) { + E_sum[itm] += contrib.at(c).getEnergy(); + } + if ((contrib.at(c).getTime() - t_min) < T_MAX[itm]) { + E_sum_telap[itm] += contrib.at(c).getEnergy(); + } + } + } + + for (int iet = 0; iet < E_TH.size(); iet++) { + for (int ies = 0; ies < E_sum.size(); ies++) { + if (E_sum[ies] > E_TH[iet]*E_TH_correction) { + hits_passed[iet][ies] += 1; + } + if (E_sum_telap[ies] > E_TH[iet]*E_TH_correction) { + hits_passed_telap[iet][ies] += 1; + } + } + } + } + + // Fill histograms with the hits_passed data + void fillThresholdHistograms( + const std::vector>& hits_passed, + const std::vector>& hits_passed_telap, + TH2* h_energy_vs_time, + TH2* h_energy_vs_telap, + TH2* h_energy_vs_time_total, + double E_TH_correction = 1.0 + ) { + for (int iet = 0; iet < E_TH.size(); iet++) { + for (int itm = 0; itm < T_MAX.size(); itm++) { + if (hits_passed[iet][itm] > 0) { + h_energy_vs_time->Fill(E_TH[iet]*E_TH_correction, T_MAX[itm]); + } + if (hits_passed_telap[iet][itm] > 0) { + h_energy_vs_telap->Fill(E_TH[iet]*E_TH_correction, T_MAX[itm]); + } + h_energy_vs_time_total->Fill(E_TH[iet]*E_TH_correction, T_MAX[itm]); + } + } + } +} + +#endif // NEUTRON_THRESHOLDS_UTIL_H diff --git a/benchmarks/nhcal_basic_distribution/scripts/basic_distribution_analysis.cxx b/benchmarks/nhcal_basic_distribution/scripts/basic_distribution_analysis.cxx index 85717f86..3c40d989 100644 --- a/benchmarks/nhcal_basic_distribution/scripts/basic_distribution_analysis.cxx +++ b/benchmarks/nhcal_basic_distribution/scripts/basic_distribution_analysis.cxx @@ -14,6 +14,7 @@ #include #include #include +#include #include "TROOT.h" #include "TRandom.h" @@ -59,40 +60,149 @@ #include "edm4eic/CalorimeterHit.h" #include "edm4eic/CalorimeterHitObj.h" - #include -#include - -// #include "FileList.h" -// #include "EICutil.h" -// #include "BasicUtil.h" +#include -// #pragma link C++ class vector+; -// #pragma link C++ class vector+; -// #pragma link C++ class vector+; +#include "NeutronThresholdUtil.h" +#include "HistogramsNeutronThresholds.h" using namespace std; using namespace ROOT; using namespace TMath; -using namespace edm4hep; -int basic_distribution_analysis(const string &filename, string outname_pdf, string outname_png) +bool debug = false; +dd4hep::Detector* det = NULL; + +inline string addPrefixAfterSlash(const string& path, const string& prefix) { + const auto slash = path.find_last_of('/'); + if (slash == string::npos) + return prefix + path; + return path.substr(0, slash + 1) + prefix + path.substr(slash + 1); +} + +string changeExtension(const string& path, const string& new_ext) +{ + size_t pos = path.find_last_of('.'); + if (pos != string::npos) + return path.substr(0, pos) + new_ext; + return path + new_ext; +} + +inline TVector3 getPlasticDimensionsCM(dd4hep::Detector& det, + const dd4hep::DDSegmentation::BitFieldCoder* dec, + dd4hep::DDSegmentation::CellID cid, + int slice_idx, + int plastic_slice_value = 3) +{ + const double NaN = numeric_limits::quiet_NaN(); + TVector3 dims(NaN, NaN, NaN); + try { + if (!dec) throw runtime_error("decoder==nullptr"); + if (dec->get(cid, slice_idx) != plastic_slice_value) + throw runtime_error("cell is not plastic"); + + auto de = det.volumeManager().lookupDetElement(cid); + if (!de.isValid()) throw runtime_error("lookupDetElement failed"); + + auto pv = de.placement(); + if (!pv.isValid()) throw runtime_error("DetElement has no placement"); + + auto vol = pv.volume(); + if (!vol.isValid()) throw runtime_error("Invalid Volume"); + + auto* box = dynamic_cast(vol.solid().ptr()); + if (!box) throw runtime_error("Solid is not TGeoBBox"); + + dims.SetXYZ(2.0 * box->GetDX(), + 2.0 * box->GetDY(), + 2.0 * box->GetDZ()); + + //const double dz_cm = box->GetDZ(); + //const double thickness_cm = 2.0 * dz_cm; + //return thickness_cm; + return dims; + } catch (const exception& e) { + cerr << "[WARN] getPlasticThicknessMM: " << e.what() << " (cellID=" << cid << ")\n"; + return dims; + } +} + +int basic_distribution_analysis(const string &filename, string outname_png, string outname_root, TString compact_file) { podio::ROOTReader *reader = new podio::ROOTReader(); reader->openFile(filename); unsigned nEvents = reader->getEntries("events"); cout << "Number of events: " << nEvents << endl; - vector xPosition; - vector yPosition; - vector zPosition; - vector rPosition; - vector energy; - vector totalEventEnergy; - vector totalEventHits; - vector zLayerValues; - vector nHitsPerLayer; - vector nEnergyPerLayer; + det = &(dd4hep::Detector::getInstance()); + det->fromCompact(compact_file.Data()); + det->volumeManager(); + det->apply("DD4hepVolumeManager", 0, 0); + + dd4hep::rec::CellIDPositionConverter cellid_converter(*det); + auto idSpec = det->readout("HcalEndcapNHits").idSpec(); + auto decoder = idSpec.decoder(); + const int slice_index = decoder->index("slice"); + if (slice_index < 0) { + cerr << "ERROR: 'slice' field not found in cell ID spec!" << endl; + return 1; + } + + constexpr double NBINS = 50; + constexpr double MIN_X_MM = -3000, MAX_X_MM = 3000; + constexpr double MIN_Y_MM = -3000, MAX_Y_MM = 3000; + constexpr double MIN_Z_MM = -4500, MAX_Z_MM = -3600; + constexpr double MIN_R_MM = 0, MAX_R_MM = 3000; + constexpr double MIN_TOTAL_ENERGY_GEV = 0, MAX_TOTAL_ENERGY_GEV = 0.3; + constexpr double MIN_HITS = 0, MAX_HITS = 300; + constexpr double MIN_LAYER_HITS = 0, MAX_LAYER_HITS = 50; + constexpr double MIN_LAYER_ENERGY_GEV = 0, MAX_LAYER_ENERGY_GEV = 0.2; + constexpr double MIN_KINETIC_GEV = 0, MAX_KINETIC_GEV = 6; + + gStyle->SetTitleSize(0.03, "XYZ"); + gStyle->SetLabelSize(0.03, "XYZ"); + gStyle->SetPadLeftMargin(0.20); + gStyle->SetPadRightMargin(0.15); + gStyle->SetPadBottomMargin(0.15); + gStyle->SetPadTopMargin(0.10); + + string outname_pdf = changeExtension(outname_png, ".pdf"); + + TH1D *h_energyTotal = new TH1D("h_energyTotal", "Total Energy per Event; Energy per Event [GeV]; N_{events}", + NBINS, MIN_TOTAL_ENERGY_GEV, MAX_TOTAL_ENERGY_GEV); + + TH2D *h_layerEnergy = new TH2D("h_layerEnergy", "Energy in Layers; Z [mm]; Energy [GeV]", + NBINS, MIN_Z_MM, MAX_Z_MM, NBINS, MIN_LAYER_ENERGY_GEV, MAX_LAYER_ENERGY_GEV); + + TProfile *p_layerEnergy = new TProfile("p_layerEnergy", "Energy in Layers; Z [mm]; Mean Energy [GeV]", + NBINS, MIN_Z_MM, MAX_Z_MM); + + TH1D *h_hitCount = new TH1D("h_hitCount", "Hits per Event; Hits per Event; N_{hits}", + NBINS, MIN_HITS, MAX_HITS); + + TH2D *h_layerHits = new TH2D("h_layerHits", "Hits in Layers; Z [mm]; N_{hits}", + NBINS, MIN_Z_MM, MAX_Z_MM, NBINS, MIN_LAYER_HITS, MAX_LAYER_HITS); + + TProfile *p_layerHits = new TProfile("p_layerHits", "Hits in Layers; Z [mm]; Mean Hits", + NBINS, MIN_Z_MM, MAX_Z_MM); + + TH2D *h_XYPos = new TH2D("h_XYPos", "Hits position;X [mm];Y [mm]", + NBINS*2, MIN_X_MM, MAX_X_MM, NBINS*2, MIN_Y_MM, MAX_Y_MM); + + TH2D *h_ZRPos = new TH2D("h_ZRPos", "Hits position;Z [mm];R [mm]", + NBINS, MIN_Z_MM, MAX_Z_MM, NBINS, MIN_R_MM, MAX_R_MM); + + TH2D *h_XYEnergy = new TH2D("h_XYEnergy", "Hits energy;X [mm];Y [mm]", + NBINS*2, MIN_X_MM, MAX_X_MM, NBINS*2, MIN_Y_MM, MAX_Y_MM); + + TH2D *h_energyRes = new TH2D("h_energyRes", "Kinetic energy vs #sum hit energy;E_{kin} [GeV]; #sum E_{hit} [GeV]", + NBINS, MIN_KINETIC_GEV, MAX_KINETIC_GEV, + NBINS, MIN_TOTAL_ENERGY_GEV, MAX_TOTAL_ENERGY_GEV); + + TProfile *p_energyRes = new TProfile("p_energyRes", "Kinetic energy vs #sum hit energy", + NBINS, MIN_KINETIC_GEV, MAX_KINETIC_GEV, "s"); + + CreateHistogamsNeutronThresholds(); for (unsigned ev = 0; ev < nEvents; ev++) { @@ -106,127 +216,133 @@ int basic_distribution_analysis(const string &filename, string outname_pdf, stri podio::Frame frame(std::move(frameData)); auto& hits = frame.get("HcalEndcapNHits"); - if (!hits.isValid()) + auto& mcCol = frame.get("MCParticles"); + if (!hits.isValid() && !mcCol.isValid()) { - cerr << "HcalEndcapNHits collection is not valid!" << endl; + cerr << "HcalEndcapNHits or MCParticles collection is not valid!" << endl; + } + + double Ekin = mcCol[0].getEnergy() - mcCol[0].getMass(); + + double MCEta = -999; + if (mcCol.size() > 0) { + auto mcpart = mcCol[0]; + TVector3 mcMom(mcpart.getMomentum().x, mcpart.getMomentum().y, mcpart.getMomentum().z); + MCEta = mcMom.Eta(); } map> layerData; double totalEnergy = 0; int totalHits = 0; + if (MCEta < -1.5 && MCEta > -3.3) { + auto hits_passed = NeutronThresholds::createHitsPassedMatrix(); + auto hits_passed_telap = NeutronThresholds::createHitsPassedMatrix(); + double thickness_cm = 1.0; + + for (const auto& hit : hits) + { + thickness_cm = getPlasticDimensionsCM(*det, decoder, hit.getCellID(), slice_index, 3).z(); + auto contrib = hit.getContributions(); + NeutronThresholds::processContributions(contrib, hits_passed, hits_passed_telap, + h_nHCal_hit_contrib_time, h_nHCal_hit_contrib_energy, thickness_cm,debug); + } + + NeutronThresholds::fillThresholdHistograms(hits_passed, hits_passed_telap, + h_nHCal_hit_contrib_energy_vs_time, h_nHCal_hit_contrib_energy_vs_telap, + h_nHCal_hit_contrib_energy_vs_time_total, thickness_cm); + } + for (const auto& hit : hits) { totalHits++; - energy.push_back(hit.getEnergy()); totalEnergy += hit.getEnergy(); auto pos = hit.getPosition(); - xPosition.push_back(pos.x); - yPosition.push_back(pos.y); - zPosition.push_back(pos.z); - rPosition.push_back(sqrt(pos.x * pos.x + pos.y * pos.y)); + double r = sqrt(pos.x * pos.x + pos.y * pos.y); + + h_XYPos->Fill(pos.x, pos.y); + h_ZRPos->Fill(pos.z, r); + h_XYEnergy->Fill(pos.x, pos.y, hit.getEnergy()); double zBin = round(pos.z); layerData[zBin].first++; layerData[zBin].second += hit.getEnergy(); } - totalEventEnergy.push_back(totalEnergy); - totalEventHits.push_back(totalHits); + h_energyTotal->Fill(totalEnergy); + h_hitCount->Fill(totalHits); + h_energyRes->Fill(Ekin, totalEnergy); + p_energyRes->Fill(Ekin, totalEnergy); for (const auto& [zValue, stats] : layerData) { - zLayerValues.push_back(zValue); - nHitsPerLayer.push_back(stats.first); - nEnergyPerLayer.push_back(stats.second); + h_layerHits->Fill(zValue, stats.first); + p_layerHits->Fill(zValue, stats.first); + h_layerEnergy->Fill(zValue, stats.second); + p_layerEnergy->Fill(zValue, stats.second); } - - } - - auto minmax_xPosition = minmax_element(xPosition.begin(), xPosition.end()); - auto minmax_yPosition = minmax_element(yPosition.begin(), yPosition.end()); - auto minmax_zPosition = minmax_element(zPosition.begin(), zPosition.end()); - auto minmax_rPosition = minmax_element(rPosition.begin(), rPosition.end()); - auto minmax_totalEventEnergy = minmax_element(totalEventEnergy.begin(), totalEventEnergy.end()); - auto minmax_totalEventHits = minmax_element(totalEventHits.begin(), totalEventHits.end()); - auto minmax_energy = minmax_element(energy.begin(), energy.end()); - auto minmax_zLayerValues = minmax_element(zLayerValues.begin(), zLayerValues.end()); - auto minmax_nHitsPerLayer = minmax_element(nHitsPerLayer.begin(), nHitsPerLayer.end()); - auto minmax_nEnergyPerLayer = minmax_element(nEnergyPerLayer.begin(), nEnergyPerLayer.end()); - - int nbins = nEvents/1000; - - TH1D *h_energyTotal = new TH1D("h_energyTotal", "Total Energy per Event; Energy per Event", - nbins/2, *minmax_totalEventEnergy.first, *minmax_totalEventEnergy.second); - TH2D *h_layerEnergy = new TH2D("h_layerEnergy", "Energy in Layers (Z); Z; Energy", - nbins/5, *minmax_zLayerValues.first, *minmax_zLayerValues.second, - nbins/4, *minmax_nEnergyPerLayer.first, *minmax_nEnergyPerLayer.second); - TProfile *p_layerEnergy = new TProfile("p_layerEnergy", "Energy in Layers (Z); Z; Mean Energy", - nbins/5, *minmax_zLayerValues.first, *minmax_zLayerValues.second); - TH1D *h_hitCount = new TH1D("h_hitCount", "Number of Hits per Event; Hits per Event", - nbins/2, *minmax_totalEventHits.first, *minmax_totalEventHits.second); - TH2D *h_layerHits = new TH2D("h_layerHits", "Number of Hits in Layers (Z); Layer(Z); Hits", - nbins/5, *minmax_zLayerValues.first, *minmax_zLayerValues.second, - *minmax_nHitsPerLayer.second, *minmax_nHitsPerLayer.first, *minmax_nHitsPerLayer.second); - TProfile *p_layerHits = new TProfile("p_layerHits", "Number of Hits in Layers (Z); Z; Mean Hits", - nbins/5, *minmax_zLayerValues.first, *minmax_zLayerValues.second); - TH2D *h_XYPos = new TH2D("h_XYPos", "Hits position X,Y;X [mm];Y [mm]", - nbins, *minmax_xPosition.first, *minmax_xPosition.second, - nbins, *minmax_yPosition.first, *minmax_yPosition.second); - TH2D *h_ZRPos = new TH2D("h_ZRPos", "Hits position Z,R;Z [mm];R [mm]", - nbins/5, *minmax_zPosition.first, *minmax_zPosition.second, - nbins, *minmax_rPosition.first, *minmax_rPosition.second); - TH2D *h_XYEnergy = new TH2D("h_XYEnergy", "Hits energy X,Y;X [mm];Y [mm]", - nbins, *minmax_xPosition.first, *minmax_xPosition.second, - nbins, *minmax_yPosition.first, *minmax_yPosition.second); - - for(int i = 0; i < xPosition.size(); i++) - { - - h_XYPos->Fill(xPosition[i], yPosition[i]); - h_ZRPos->Fill(zPosition[i], rPosition[i]); - h_XYEnergy->Fill(xPosition[i], yPosition[i], energy[i]); - } - - for (int i = 0; i < zLayerValues.size(); i++) - { - h_layerHits->Fill(zLayerValues[i], nHitsPerLayer[i]); - p_layerHits->Fill(zLayerValues[i], nHitsPerLayer[i]); - - h_layerEnergy->Fill(zLayerValues[i], nEnergyPerLayer[i]); - p_layerEnergy->Fill(zLayerValues[i], nEnergyPerLayer[i]); } - - for(int i = 0; i < nEvents; i++) - { - h_energyTotal->Fill(totalEventEnergy[i]); - h_hitCount->Fill(totalEventHits[i]); - } - - gStyle->SetPadLeftMargin(0.13); - - TCanvas *canvas = new TCanvas("canvas", "canvas", 1600, 800); - canvas->Divide(4,2); - canvas->cd(1); + + TFile *outFile = new TFile(outname_root.c_str(), "RECREATE"); + h_energyRes->Write(); + p_energyRes->Write(); + outFile->Close(); + delete outFile; + + TCanvas *c_evLayers = new TCanvas("c_evLayers", "c_evLayers", 1600, 800); + c_evLayers->Divide(2,2); + c_evLayers->cd(1); h_energyTotal->Draw(); - canvas->cd(2); + c_evLayers->cd(2); h_layerEnergy->Draw("COLZ"); + p_layerEnergy->SetLineWidth(3); p_layerEnergy->SetLineColor(kRed); p_layerEnergy->SetMarkerColor(kRed); p_layerEnergy->Draw("SAME"); - canvas->cd(3); + c_evLayers->cd(3); h_hitCount->Draw(); - canvas->cd(4); + c_evLayers->cd(4); h_layerHits->Draw("COLZ"); + p_layerHits->SetLineWidth(3); p_layerHits->SetLineColor(kRed); p_layerHits->SetMarkerColor(kRed); p_layerHits->Draw("SAME"); - canvas->cd(5); + + c_evLayers->SaveAs(outname_png.c_str()); + c_evLayers->SaveAs(outname_pdf.c_str()); + + TCanvas *c_hit_posE = new TCanvas("c_hit_posE", "c_hit_posE", 1600, 800); + c_hit_posE->Divide(2,2); + c_hit_posE->cd(1); h_XYPos->Draw("COLZ"); - canvas->cd(6); + c_hit_posE->cd(2); h_ZRPos->Draw("COLZ"); - canvas->cd(7); + c_hit_posE->cd(3); h_XYEnergy->Draw("COLZ"); + c_hit_posE->cd(4); + + c_hit_posE->SaveAs(addPrefixAfterSlash(outname_png, "hit_posE_").c_str()); + c_hit_posE->SaveAs(addPrefixAfterSlash(outname_pdf, "hit_posE_").c_str()); + + // TCanvas *c_neutronThresholds = new TCanvas("c_neutronThresholds", "c_neutronThresholds", 1600, 800); + // c_neutronThresholds->Divide(3,2); + // c_neutronThresholds->cd(1); + // h_nHCal_hit_contrib_time->Draw(); + // c_neutronThresholds->cd(2); + // h_nHCal_hit_contrib_energy->Draw(); + // c_neutronThresholds->cd(3); + // h_nHCal_hit_contrib_energy_vs_time->Draw("COLZ"); + // c_neutronThresholds->cd(4); + // h_nHCal_hit_contrib_energy_vs_telap->Draw("COLZ"); + // c_neutronThresholds->cd(5); + // h_nHCal_hit_contrib_energy_vs_time_total->Draw("COLZ"); + + //c_neutronThresholds->SaveAs(addPrefixAfterSlash(outname_png, "neutronThresholds_").c_str()); + //c_neutronThresholds->SaveAs(addPrefixAfterSlash(outname_pdf, "neutronThresholds_").c_str()); + - canvas->SaveAs(outname_pdf.c_str()); - canvas->SaveAs(outname_png.c_str()); + delete reader; + delete c_evLayers; + delete c_hit_posE; + // delete c_neutronThresholds; + DeleteHistogamsNeutronThresholds(); return 0; -} +} \ No newline at end of file diff --git a/benchmarks/nhcal_basic_distribution/scripts/basic_distribution_energy_resolution.cxx b/benchmarks/nhcal_basic_distribution/scripts/basic_distribution_energy_resolution.cxx new file mode 100644 index 00000000..7fa52a10 --- /dev/null +++ b/benchmarks/nhcal_basic_distribution/scripts/basic_distribution_energy_resolution.cxx @@ -0,0 +1,162 @@ +#include +#include +#include +#include +#include + +#include "ROOT/RDataFrame.hxx" +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "TROOT.h" +#include "TRandom.h" +#include "TH3.h" + + +#include "DD4hep/Detector.h" +#include "DDRec/CellIDPositionConverter.h" + +#include +#include +#include "podio/ROOTReader.h" +//#include +#include "podio/CollectionIDTable.h" +#include "podio/ObjectID.h" + +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/MCParticleCollectionData.h" +#include "edm4hep/MCParticle.h" +#include "edm4hep/MCParticleData.h" + +#include "edm4hep/SimCalorimeterHitCollectionData.h" +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4hep/SimCalorimeterHitData.h" +#include "edm4hep/SimCalorimeterHit.h" + +#include "edm4hep/CalorimeterHit.h" +#include "edm4hep/CalorimeterHitCollectionData.h" +#include "edm4hep/CalorimeterHitCollection.h" +#include "edm4hep/CalorimeterHitData.h" +#include "edm4hep/CalorimeterHit.h" +#include "edm4hep/CalorimeterHitObj.h" + +#include "edm4eic/ClusterCollection.h" +#include "edm4eic/Cluster.h" +#include "edm4eic/ClusterData.h" + +#include "edm4eic/CalorimeterHit.h" +#include "edm4eic/CalorimeterHitCollectionData.h" +#include "edm4eic/CalorimeterHitCollection.h" +#include "edm4eic/CalorimeterHitData.h" +#include "edm4eic/CalorimeterHit.h" +#include "edm4eic/CalorimeterHitObj.h" + + +#include +#include +#include +#include + +using namespace std; +using namespace ROOT; +using namespace TMath; +using namespace edm4hep; + +string changeExtension(const string& path, const string& new_ext) +{ + size_t pos = path.find_last_of('.'); + if (pos != string::npos) + return path.substr(0, pos) + new_ext; + return path + new_ext; +} + +int basic_distribution_energy_resolution(const string &filename, string outname_png){ + + gStyle->SetTitleSize(0.045, "XYZ"); + gStyle->SetLabelSize(0.04, "XYZ"); + gStyle->SetPadLeftMargin(0.25); + gStyle->SetPadRightMargin(0.15); + gStyle->SetPadBottomMargin(0.15); + gStyle->SetPadTopMargin(0.10); + + string outname_pdf = changeExtension(outname_png, ".pdf"); + + TFile *file = TFile::Open(filename.c_str(), "READ"); + if (!file || file->IsZombie()) + { + cerr << "Cannot open file: " << filename << endl; + return 1; + } + + TH2D *h_energyRes = (TH2D*)file->Get("h_energyRes"); + TProfile *p_energyRes = (TProfile*)file->Get("p_energyRes"); + + if (!h_energyRes || !p_energyRes) + { + cerr << "Cannot retrieve histograms from file" << endl; + file->Close(); + return 1; + } + + TGraphErrors *g_resolution = new TGraphErrors(); + + for (int i = 1; i <= p_energyRes->GetNbinsX(); i++) { + double Ekin = p_energyRes->GetBinCenter(i); + double mean = p_energyRes->GetBinContent(i); + double rms = p_energyRes->GetBinError(i); + int entries = p_energyRes->GetBinEntries(i); + + if (entries < 10 || mean <= 0) continue; + + int n = g_resolution->GetN(); + g_resolution->SetPoint(n, Ekin, rms/mean); + g_resolution->SetPointError(n, 0, rms/mean * sqrt(1.0/entries)); + } + + TF1 *fit = new TF1("fit", "[0]*pow(x, -0.5) + [1] + [2]*pow(x, 1)", 0.1, 6); + fit->SetParameters(0.5, 0.05, 0.001); + fit->SetParNames("a (stochastic)", "b (constant)", "c (noise)"); + g_resolution->Fit(fit, "R"); + + TCanvas *c = new TCanvas("c", "Energy Resolution", 1600, 800); + c->Divide(2, 1); + c->cd(1); + + h_energyRes->Draw("COLZ"); + p_energyRes->SetLineWidth(3); + p_energyRes->SetLineColor(kRed); + p_energyRes->SetMarkerColor(kRed); + p_energyRes->Draw("SAME"); + + c->cd(2); + g_resolution->SetTitle("Energy Resolution;E_{kin} [GeV];#sigma/E"); + g_resolution->SetMarkerStyle(20); + g_resolution->SetLineWidth(2); + g_resolution->SetLineColor(kBlue); + g_resolution->Draw("APE"); + fit->Draw("SAME"); + + double par_a = fit->GetParameter(0); + double par_b = fit->GetParameter(1); + double par_c = fit->GetParameter(2); + + TLegend *leg = new TLegend(0.35, 0.65, 0.80, 0.80); + leg->SetTextSize(0.025); + leg->AddEntry(g_resolution, "Data", "p"); + leg->AddEntry(fit, Form("Fit: %.3f/#sqrt{E} + %.3f + %.3f#timesE", par_a, par_b, par_c), "l"); + leg->Draw(); + + c->SaveAs(outname_png.c_str()); + c->SaveAs(outname_pdf.c_str()); + + file->Close(); + + return 0; +} \ No newline at end of file diff --git a/benchmarks/nhcal_dimuon_photoproduction/Snakefile b/benchmarks/nhcal_dimuon_photoproduction/Snakefile new file mode 100644 index 00000000..733ab06b --- /dev/null +++ b/benchmarks/nhcal_dimuon_photoproduction/Snakefile @@ -0,0 +1,86 @@ +rule nhcal_dimuon_photoproduction_presimulate: + output: + "sim_input/DiMuon_ep_18x275GeV.{INDEX}.hepmc3", + params: + outdir=lambda wc, output: os.path.dirname(output[0]), + shell: + """ + pwd + echo "Output dir: {params.outdir}" + make -C benchmarks/nhcal_dimuon_photoproduction/scripts + benchmarks/nhcal_dimuon_photoproduction/scripts/run.sh {wildcards.INDEX} {params.outdir} + ls -la {params.outdir}/ + """ + +rule nhcal_dimuon_photoproduction_simulate: + input: + "sim_input/DiMuon_ep_18x275GeV.{INDEX}.hepmc3", + output: + "sim_output/nhcal_dimuon_photoproduction/sim_{DETECTOR_CONFIG}.{INDEX}.edm4hep.root", + params: + N_EVENTS=1000, + SEED=lambda wildcards: "1" + wildcards.INDEX, + DETECTOR_PATH=os.environ["DETECTOR_PATH"], + DETECTOR_CONFIG=lambda wildcards: wildcards.DETECTOR_CONFIG, + DD4HEP_HASH=get_spack_package_hash("dd4hep"), + NPSIM_HASH=get_spack_package_hash("npsim"), + wildcard_constraints: + INDEX = r"\d+", + cache: True + shell: + """ + exec npsim \ + --compactFile {params.DETECTOR_PATH}/{params.DETECTOR_CONFIG}.xml \ + --inputFile {input} \ + -v WARNING \ + --random.seed {params.SEED} \ + --numberOfEvents {params.N_EVENTS} \ + --outputFile {output} + """ + +rule nhcal_dimuon_photoproduction_recon: + input: + "sim_output/nhcal_dimuon_photoproduction/sim_{DETECTOR_CONFIG}.{INDEX}.edm4hep.root" + output: + "sim_output/nhcal_dimuon_photoproduction/sim_{DETECTOR_CONFIG}.{INDEX}.eicrecon.edm4hep.root" + params: + DETECTOR_CONFIG = lambda wildcards: wildcards.DETECTOR_CONFIG, + cache: True + shell: + """ + exec env DETECTOR_CONFIG={params.DETECTOR_CONFIG} \ + eicrecon {input} \ + -Ppodio:output_file={output} \ + + """ + +rule nhcal_dimuon_photoproduction_combine: + input: + lambda wildcards: expand( + "sim_output/nhcal_dimuon_photoproduction/sim_{DETECTOR_CONFIG}.{INDEX}.eicrecon.edm4hep.root", + DETECTOR_CONFIG=wildcards.DETECTOR_CONFIG, + INDEX=range(int(wildcards.N)), + ) + output: + "sim_output/nhcal_dimuon_photoproduction/sim_{DETECTOR_CONFIG}_{N}files.eicrecon.edm4hep.root", + wildcard_constraints: + N = r"\d+", + shell: + """ + hadd -f {output} {input} + """ + +rule nhcal_dimuon_photoproduction_analysis: + input: + combined="sim_output/nhcal_dimuon_photoproduction/sim_{DETECTOR_CONFIG}_{N}files.eicrecon.edm4hep.root", + script="benchmarks/nhcal_dimuon_photoproduction/scripts/dimuon_photoproduction_analysis.cxx", + output: + png="results/nhcal_dimuon_photoproduction/analysis_{DETECTOR_CONFIG}_{N}files.png", + pdf="results/nhcal_dimuon_photoproduction/analysis_{DETECTOR_CONFIG}_{N}files.pdf", + params: + DETECTOR_PATH = os.environ["DETECTOR_PATH"], + DETECTOR_CONFIG = lambda wildcards: f"{wildcards.DETECTOR_CONFIG}.xml", + shell: + """ + root -l -b -q '{input.script}+("{input.combined}","{output.pdf}","{output.png}","{params.DETECTOR_PATH}/{params.DETECTOR_CONFIG}")' + """ diff --git a/benchmarks/nhcal_dimuon_photoproduction/config.yml b/benchmarks/nhcal_dimuon_photoproduction/config.yml new file mode 100644 index 00000000..fb1805b7 --- /dev/null +++ b/benchmarks/nhcal_dimuon_photoproduction/config.yml @@ -0,0 +1,28 @@ +sim:nhcal_dimuon_photoproduction: + extends: .det_benchmark + stage: simulate + parallel: + matrix: + - INDEX_RANGE: ["0 1", "2 3", "4 5", "6 7", "8 9"] + script: + - snakemake $SNAKEMAKE_FLAGS --cores $MAX_CORES_PER_JOB $(for INDEX in $INDEX_RANGE; do echo sim_output/nhcal_dimuon_photoproduction/sim_epic_full.${INDEX}.edm4hep.root; done) + - snakemake $SNAKEMAKE_FLAGS --cores $MAX_CORES_PER_JOB $(for INDEX in $INDEX_RANGE; do echo sim_output/nhcal_dimuon_photoproduction/sim_epic_full.${INDEX}.eicrecon.edm4hep.root; done) + +bench:nhcal_dimuon_photoproduction_analysis: + extends: .det_benchmark + stage: benchmarks + needs: + - "sim:nhcal_dimuon_photoproduction" + script: + - snakemake $SNAKEMAKE_FLAGS --cores $MAX_CORES_PER_JOB results/nhcal_dimuon_photoproduction/analysis_epic_full_10files.pdf + +collect_results:nhcal_dimuon_photoproduction: + extends: .det_benchmark + stage: collect + needs: + - "bench:nhcal_dimuon_photoproduction_analysis" + script: + - ls -lrht + - mv results{,_save}/ + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/nhcal_dimuon_photoproduction/analysis_epic_full_10files.pdf + - mv results{_save,}/ diff --git a/benchmarks/nhcal_dimuon_photoproduction/scripts/HistogramsPythia.h b/benchmarks/nhcal_dimuon_photoproduction/scripts/HistogramsPythia.h new file mode 100755 index 00000000..e1c07ac4 --- /dev/null +++ b/benchmarks/nhcal_dimuon_photoproduction/scripts/HistogramsPythia.h @@ -0,0 +1,692 @@ +/* + * HistogramsPythia.h + * + * Created on: 30 sie 2024 + * Author: Khaless + */ + +#ifndef HISTOGRAMSPYTHIA_H_ +#define HISTOGRAMSPYTHIA_H_ + +#include "TH1.h" +#include "TH2.h" +#include "TH3.h" +#include "TMath.h" + +int CreateHistograms(); + + +void MakeLogBins(double *array, int nbins, double binLo, double binHi) +{ + //double exp = (nbinEdg-1)/nBinPerDecimal; + + // Calculate the logarithmic bin edges + double logMin = log10(binLo); + double logMax = log10(binHi); + double binWidth = (logMax - logMin) / nbins; + + for (int i = 0; i <= nbins; ++i) { + array[i] = pow(10, logMin + i * binWidth); + } +} + +// Event + +TH1F *h_Events; +TH1F *h_Events_types; +TH1F *h_Events_cuts; + +TH1F *h_Events_Diffractive; + +TH1F *h_Events_nPartonsOut; + +TH1F *h_XsecGen; +TH1F *h_XsecSel; + +TH2F *h_XsecGen_err; +TH2F *h_XsecSel_err; + +TH1F *h_Event_nPart_final; +TH1F *h_Event_nJets; +TH1F *h_Event_nJets_meas; +TH1F *h_Event_nJets_meas_no_nHCal; + +TH2F *h_Event_xQ2; +TH2F *h_Event_yQ2; +TH2F *h_Event_xy; + +TH2F *h_Event_nHCal_xQ2; +TH2F *h_Event_nHCal_yQ2; +TH2F *h_Event_nHCal_xy; + +TH1F *h_Event_nPion_p; +TH1F *h_Event_nPion_n; +TH1F *h_Event_nKaon_p; +TH1F *h_Event_nKaon_n; +TH1F *h_Event_nProton_p; +TH1F *h_Event_nProton_n; +TH1F *h_Event_nElectron_p; +TH1F *h_Event_nElectron_n; + +TH1F *h_Event_nNeutron; +TH1F *h_Event_nGamma; + +// special + + +TH1F *h_Event_Q2; +TH1F *h_Event_x; +TH1F *h_Event_y; + +TH1F *h_Event_jets_Q2; +TH1F *h_Event_jets_x; +TH1F *h_Event_jets_y; + +TH1F *h_Event_nHCal_0_Q2; +TH1F *h_Event_nHCal_0_x; +TH1F *h_Event_nHCal_0_y; + +TH1F *h_Event_nHCal_1_Q2; +TH1F *h_Event_nHCal_1_x; +TH1F *h_Event_nHCal_1_y; + +TH1F *h_Event_nHCal_2_Q2; +TH1F *h_Event_nHCal_2_x; +TH1F *h_Event_nHCal_2_y; + +TH1F *h_Event_AllHCal_Q2; +TH1F *h_Event_AllHCal_x; +TH1F *h_Event_AllHCal_y; + + +TH1F *h_Event_JetMeas_nHCal_0_Q2; +TH1F *h_Event_JetMeas_nHCal_0_x; +TH1F *h_Event_JetMeas_nHCal_0_y; + +TH1F *h_Event_JetMeas_nHCal_1_Q2; +TH1F *h_Event_JetMeas_nHCal_1_x; +TH1F *h_Event_JetMeas_nHCal_1_y; + +TH1F *h_Event_JetMeas_nHCal_2_Q2; +TH1F *h_Event_JetMeas_nHCal_2_x; +TH1F *h_Event_JetMeas_nHCal_2_y; + +TH1F *h_Event_JetMeas_AllHCal_Q2; +TH1F *h_Event_JetMeas_AllHCal_x; +TH1F *h_Event_JetMeas_AllHCal_y; + +TH1F *h_Event_JetMeas_no_nHCal_Q2; +TH1F *h_Event_JetMeas_no_nHCal_x; +TH1F *h_Event_JetMeas_no_nHCal_y; + +TH2F *h_Event_HCal_jets; +TH2F *h_Event_HCal_jets_meas; +TH2F *h_Event_HCal_jets_meas_no_nHCal; + +//full +TH2F *h_Event_HCal_jets_meas_full; + + +// Outgoing partons + +TH1F *h_Partons_status; +TH1F *h_Partons_types; +TH1F *h_Partons_types_anti; + +TH2F *h_Partons_eta; +TH2F *h_Partons_phi; +TH2F *h_Partons_p; +TH2F *h_Partons_pT; + + +TH2F *h_Parton_eta_pT; +TH2F *h_Parton_eta_p; +TH2F *h_Parton_eta_E; + +TH2F *h_Parton_x_eta; +TH2F *h_Parton_y_eta; + +TH2F *h_Parton_x_eta1; +TH2F *h_Parton_x_eta2; + +TH2F *h_Parton_y_eta1; +TH2F *h_Parton_y_eta2; + +// Particles + +TH1F *h_Particle_eta; + +TH2F *h_Particle_eta_p; +TH2F *h_Particle_eta_pT; +TH2F *h_Particle_eta_E; + + +// weighted +TH1F *h_Particle_eta_wE; + +// eta, momentum +TH2F *h_Particle_pion_p_eta_p; +TH2F *h_Particle_pion_n_eta_p; +TH2F *h_Particle_Kaon_p_eta_p; +TH2F *h_Particle_Kaon_n_eta_p; +TH2F *h_Particle_proton_p_eta_p; +TH2F *h_Particle_proton_n_eta_p; +TH2F *h_Particle_Electron_p_eta_p; +TH2F *h_Particle_Electron_n_eta_p; + +TH2F *h_Particle_Neutron_eta_p; +TH2F *h_Particle_Gamma_eta_p; + +// eta, transverse momentum pT +TH2F *h_Particle_pion_p_eta_pT; +TH2F *h_Particle_pion_n_eta_pT; +TH2F *h_Particle_Kaon_p_eta_pT; +TH2F *h_Particle_Kaon_n_eta_pT; +TH2F *h_Particle_proton_p_eta_pT; +TH2F *h_Particle_proton_n_eta_pT; +TH2F *h_Particle_Electron_p_eta_pT; +TH2F *h_Particle_Electron_n_eta_pT; + +TH2F *h_Particle_Neutron_eta_pT; +TH2F *h_Particle_Gamma_eta_pT; + +// eta, energy +TH2F *h_Particle_Pion_p_eta_E; +TH2F *h_Particle_Pion_n_eta_E; +TH2F *h_Particle_Kaon_p_eta_E; +TH2F *h_Particle_Kaon_n_eta_E; +TH2F *h_Particle_Proton_p_eta_E; +TH2F *h_Particle_Proton_n_eta_E; +TH2F *h_Particle_Electron_p_eta_E; +TH2F *h_Particle_Electron_n_eta_E; + +TH2F *h_Particle_Neutron_eta_E; +TH2F *h_Particle_Gamma_eta_E; + + +// Jets +TH1D *h_Jet_nPart; +TH1D *h_Jet_mass; +TH1D *h_Jet_charge; +TH1D *h_Jet_E; +TH1D *h_Jet_p; +TH1D *h_Jet_pT; +TH1D *h_Jet_eta; +TH1D *h_Jet_deta; + +TH2F *h_Jets_eta; +TH2F *h_Jets_phi; +TH2F *h_Jets_p; +TH2F *h_Jets_pT; +TH2F *h_Jets_E; + +// Jets measured +TH1D *h_Jet_meas_nPart; +TH1D *h_Jet_meas_mass; +TH1D *h_Jet_meas_charge; +TH1D *h_Jet_meas_E; +TH1D *h_Jet_meas_p; +TH1D *h_Jet_meas_pT; +TH1D *h_Jet_meas_eta; +TH1D *h_Jet_meas_deta; + +TH2F *h_Jets_meas_eta; +TH2F *h_Jets_meas_phi; +TH2F *h_Jets_meas_p; +TH2F *h_Jets_meas_pT; +TH2F *h_Jets_meas_E; + +TH1D *h_Jet_bHCal_part_eta; +TH1D *h_Jet_meas_bHCal_part_eta; + +TH2D *h_Jet_HCal_part_eta; +TH2D *h_Jet_meas_HCal_part_eta; + + +// Jets measured without nHCal +TH1D *h_Jet_meas_no_nHCal_nPart; +TH1D *h_Jet_meas_no_nHCal_mass; +TH1D *h_Jet_meas_no_nHCal_charge; +TH1D *h_Jet_meas_no_nHCal_E; +TH1D *h_Jet_meas_no_nHCal_p; +TH1D *h_Jet_meas_no_nHCal_pT; +TH1D *h_Jet_meas_no_nHCal_eta; +TH1D *h_Jet_meas_no_nHCal_deta; + +TH2F *h_Jets_meas_no_nHCal_eta; +TH2F *h_Jets_meas_no_nHCal_phi; +TH2F *h_Jets_meas_no_nHCal_p; +TH2F *h_Jets_meas_no_nHCal_pT; +TH2F *h_Jets_meas_no_nHCal_E; + +// measured jets vs. partons +TH2F *h_Jets_meas_Partons_eta; +TH2F *h_Jets_meas_Partons_phi; +TH2F *h_Jets_meas_Partons_E; + +TH2F *h_Jet_meas_Parton_eta1; +TH2F *h_Jet_meas_Parton_phi1; +TH2F *h_Jet_meas_Parton_E1; +TH2F *h_Jet_meas_Parton_eta2; +TH2F *h_Jet_meas_Parton_phi2; +TH2F *h_Jet_meas_Parton_E2; + + +// temp + +TH1F *hist_eta_energy_tmp; +TH1F *hist_eta_energy_denom_tmp; + + + + +int CreateHistograms() +{ + + + const int nbins_x = 200; + + const int nbinEdg_x = nbins_x+1; + + double *logBinsArray_x = new double[nbins_x]; + + MakeLogBins(logBinsArray_x, nbins_x, 10e-7, 1.0); + + + // Event + + h_Events = new TH1F("h_Events", "Number of events; events; counts", 10, 0.0, 10.0); + h_Events_types = new TH1F("h_Events_types", "Number of events; type; counts", 10, 0.0, 10.0); + + h_Events_types->GetXaxis()->SetBinLabel(1, "All events"); + h_Events_types->GetXaxis()->SetBinLabel(2, "1 parton in ePIC acc."); + h_Events_types->GetXaxis()->SetBinLabel(3, "2 partons in ePIC acc."); + h_Events_types->GetXaxis()->SetBinLabel(4, "1 jet in ePIC acc."); + h_Events_types->GetXaxis()->SetBinLabel(5, "2 jets in ePIC acc."); + h_Events_types->GetXaxis()->SetBinLabel(6, "2 or more jets in ePIC acc."); + + h_Events_cuts = new TH1F("h_Events_cuts", "Number of events; selection; counts", 10, 0.0, 10.0); + + h_Events_cuts->GetXaxis()->SetBinLabel(1, "Generated"); + h_Events_cuts->GetXaxis()->SetBinLabel(2, "Diffractive A or B"); + h_Events_cuts->GetXaxis()->SetBinLabel(3, "Hard diffractive A or B"); + h_Events_cuts->GetXaxis()->SetBinLabel(4, "2 partons in ePIC acc."); + + h_Events_Diffractive = new TH1F("h_Events_Diffractive", "Number of diffractive events; events; counts", 10, 0.0, 10.0); + + h_Events_Diffractive->GetXaxis()->SetBinLabel(1, "Diffractive A"); + h_Events_Diffractive->GetXaxis()->SetBinLabel(2, "Diffractive B"); + h_Events_Diffractive->GetXaxis()->SetBinLabel(3, "Hard diffractive A"); + h_Events_Diffractive->GetXaxis()->SetBinLabel(4, "Hard diffractive B"); + + h_XsecGen = new TH1F("h_XsecGen", "Generated event cross-section; #sigma [mb]; counts", 10000, 0.0, 0.001); + h_XsecSel = new TH1F("h_XsecSel", "Selected event cross-section; #sigma [mb]; counts", 10000, 0.0, 0.001); + + h_XsecGen_err = new TH2F("h_XsecGen_err", "Generated event cross-section vs. uncertainty; #sigma [mb]; #sigma_{#sigma} [mb]; counts", 10000, 0.0, 0.001, 1000, 0.0, 0.00001); + h_XsecSel_err = new TH2F("h_XsecSel_err", "Selected event cross-section vs. uncertainty; #sigma [mb]; #sigma_{#sigma} [mb]; counts", 10000, 0.0, 0.001, 1000, 0.0, 0.00001); + + + h_Events_nPartonsOut = new TH1F("h_Events_nPartonsOut", "Number of outgoing partons; N_{out} [1]; counts", 21, -0.5, 20.5); + + h_Event_nPart_final = new TH1F("h_Event_nPart_final", "Number of final MC particles; N_{MC} [1]; counts", 2001, -0.5, 2000.5); + h_Event_nJets = new TH1F("h_Event_nJets", "Number of jets; N_{jet} [1]; counts", 21, -0.5, 20.5); + h_Event_nJets_meas = new TH1F("h_Event_nJets_meas", "Number of measured jets; N_{jet} [1]; counts", 21, -0.5, 20.5); + h_Event_nJets_meas_no_nHCal = new TH1F("h_Event_nJets_meas_no_nHCal", "Number of measured jets wihout nHCal; N_{jet} [1]; counts", 21, -0.5, 20.5); + + + h_Event_xQ2 = new TH2F("h_Event_xQ2", "Event Q^{2} vs. x; x [1]; Q^{2} [GeV^{2}/c^{2}]; counts", nbins_x, logBinsArray_x, 500, 0.0, 5.0); + h_Event_yQ2 = new TH2F("h_Event_yQ2", "Event Q^{2} vs. inelasticity y; y [1]; Q^{2} [GeV^{2}/c^{2}]; counts", 1000, 0.0, 1.0, 500, 0.0, 5.0); + h_Event_xy = new TH2F("h_Event_xy", "Event inelasticity y vs. x; x [1]; y [1]; counts", nbins_x, logBinsArray_x, 1000, 0.0, 1.0); + + + h_Event_nHCal_xQ2 = new TH2F("h_Event_nHCal_xQ2", "Event with nHCal activity Q^{2} vs. x; x [1]; Q^{2} [GeV^{2}/c^{2}]; counts", nbins_x, logBinsArray_x, 500, 0.0, 5.0); + h_Event_nHCal_yQ2 = new TH2F("h_Event_nHCal_yQ2", "Event with nHCal activity Q^{2} vs. inelasticity y; y [1]; Q^{2} [GeV^{2}/c^{2}]; counts", 1000, 0.0, 1.0, 500, 0.0, 5.0); + h_Event_nHCal_xy = new TH2F("h_Event_nHCal_xy", "Event with nHCal activity inelasticity y vs. x; x [1]; y [1]; counts", nbins_x, logBinsArray_x, 1000, 0.0, 1.0); + + + h_Event_nPion_p = new TH1F("h_Event_nPion_p", "Number of MC particles #pi^{+}; N_{MC} [1]; counts", 2001, -0.5, 2000.5); + h_Event_nPion_n = new TH1F("h_Event_nPion_n", "Number of MC particles #pi^{-}; N_{MC} [1]; counts", 2001, -0.5, 2000.5); + h_Event_nKaon_p = new TH1F("h_Event_nKaon_p", "Number of MC particles K^{+}; N_{MC} [1]; counts", 2001, -0.5, 2000.5); + h_Event_nKaon_n = new TH1F("h_Event_nKaon_n", "Number of MC particles K^{-}; N_{MC} [1]; counts", 2001, -0.5, 2000.5); + h_Event_nProton_p = new TH1F("h_Event_nProton_p", "Number of MC particles p^{+}; N_{MC} [1]; counts", 2001, -0.5, 2000.5); + h_Event_nProton_n = new TH1F("h_Event_nProton_n", "Number of MC particles p^{-}; N_{MC} [1]; counts", 2001, -0.5, 2000.5); + h_Event_nElectron_p = new TH1F("h_Event_nElectron_p", "Number of MC particles e^{+}; N_{MC} [1]; counts", 2001, -0.5, 2000.5); + h_Event_nElectron_n = new TH1F("h_Event_nElectron_n", "Number of MC particles e^{-}; N_{MC} [1]; counts", 2001, -0.5, 2000.5); + + h_Event_nNeutron = new TH1F("h_Event_nNeutron", "Number of MC particles n; N_{MC} [1]; counts", 2001, -0.5, 2000.5); + h_Event_nGamma = new TH1F("h_Event_nGamma", "Number of MC particles #gamma; N_{MC} [1]; counts", 2001, -0.5, 2000.5); + + + // special + + h_Event_Q2 = new TH1F("h_Event_Q2", "Event Q^{2}; Q^{2} [GeV^{2}/c^{2}]; counts", 1000, 0.0, 10.0); + h_Event_x = new TH1F("h_Event_x", "Event x; x [1]; counts", nbins_x, logBinsArray_x); + h_Event_y = new TH1F("h_Event_y", "Event inelasticity y; y [1]; counts", 1000, 0.0, 1.0); + + h_Event_jets_Q2 = new TH1F("h_Event_jets_Q2", "Event with jets>=2 Q^{2}; Q^{2} [GeV^{2}/c^{2}]; counts", 1000, 0.0, 10.0); + h_Event_jets_x = new TH1F("h_Event_jets_x", "Event with jets>=2 x; x [1]; counts", nbins_x, logBinsArray_x); + h_Event_jets_y = new TH1F("h_Event_jets_y", "Event with jets>=2 inelasticity y; y [1]; counts", 1000, 0.0, 1.0); + + h_Event_nHCal_0_Q2 = new TH1F("h_Event_nHCal_0_Q2", "Event with 0 jets in nHCal Q^{2}; Q^{2} [GeV^{2}/c^{2}]; counts", 1000, 0.0, 10.0); + h_Event_nHCal_0_x = new TH1F("h_Event_nHCal_0_x", "Event with 0 jets in nHCal x; x [1]; counts", nbins_x, logBinsArray_x); + h_Event_nHCal_0_y = new TH1F("h_Event_nHCal_0_y", "Event with 0 jets in nHCal inelasticity y; y [1]; counts", 1000, 0.0, 1.0); + + h_Event_nHCal_1_Q2 = new TH1F("h_Event_nHCal_1_Q2", "Event with 1 jet in nHCal Q^{2}; Q^{2} [GeV^{2}/c^{2}]; counts", 1000, 0.0, 10.0); + h_Event_nHCal_1_x = new TH1F("h_Event_nHCal_1_x", "Event with 1 jet in nHCal x; x [1]; counts", nbins_x, logBinsArray_x); + h_Event_nHCal_1_y = new TH1F("h_Event_nHCal_1_y", "Event with 1 jet in nHCal inelasticity y; y [1]; counts", 1000, 0.0, 1.0); + + h_Event_nHCal_2_Q2 = new TH1F("h_Event_nHCal_2_Q2", "Event with 2 jets in nHCal Q^{2}; Q^{2} [GeV^{2}/c^{2}]; counts", 1000, 0.0, 10.0); + h_Event_nHCal_2_x = new TH1F("h_Event_nHCal_2_x", "Event with 2 jets in nHCal x; x [1]; counts", nbins_x, logBinsArray_x); + h_Event_nHCal_2_y = new TH1F("h_Event_nHCal_2_y", "Event with 2 jets in nHCal inelasticity y; y [1]; counts", 1000, 0.0, 1.0); + + h_Event_AllHCal_Q2 = new TH1F("h_Event_AllHCal_Q2", "Event with jets in any HCal Q^{2}; Q^{2} [GeV^{2}/c^{2}]; counts", 1000, 0.0, 10.0); + h_Event_AllHCal_x = new TH1F("h_Event_AllHCal_x", "Event with jets in any HCal x; x [1]; counts", nbins_x, logBinsArray_x); + h_Event_AllHCal_y = new TH1F("h_Event_AllHCal_y", "Event with jets in any HCal inelasticity y; y [1]; counts", 1000, 0.0, 1.0); + + + h_Event_JetMeas_nHCal_0_Q2 = new TH1F("h_Event_JetMeas_nHCal_0_Q2", "Event with 0 jets measured in nHCal Q^{2}; Q^{2} [GeV^{2}/c^{2}]; counts", 1000, 0.0, 10.0); + h_Event_JetMeas_nHCal_0_x = new TH1F("h_Event_JetMeas_nHCal_0_x", "Event with 0 jets measured in nHCal x; x [1]; counts", nbins_x, logBinsArray_x); + h_Event_JetMeas_nHCal_0_y = new TH1F("h_Event_JetMeas_nHCal_0_y", "Event with 0 jets measured in nHCal inelasticity y; y [1]; counts", 1000, 0.0, 1.0); + + h_Event_JetMeas_nHCal_1_Q2 = new TH1F("h_Event_JetMeas_nHCal_1_Q2", "Event with 1 jet measured in nHCal Q^{2}; Q^{2} [GeV^{2}/c^{2}]; counts", 1000, 0.0, 10.0); + h_Event_JetMeas_nHCal_1_x = new TH1F("h_Event_JetMeas_nHCal_1_x", "Event with 1 jet measured in nHCal x; x [1]; counts", nbins_x, logBinsArray_x); + h_Event_JetMeas_nHCal_1_y = new TH1F("h_Event_JetMeas_nHCal_1_y", "Event with 1 jet measured in nHCal inelasticity y; y [1]; counts", 1000, 0.0, 1.0); + + h_Event_JetMeas_nHCal_2_Q2 = new TH1F("h_Event_JetMeas_nHCal_2_Q2", "Event with 2 jets measured in nHCal Q^{2}; Q^{2} [GeV^{2}/c^{2}]; counts", 1000, 0.0, 10.0); + h_Event_JetMeas_nHCal_2_x = new TH1F("h_Event_JetMeas_nHCal_2_x", "Event with 2 jets measured in nHCal x; x [1]; counts", nbins_x, logBinsArray_x); + h_Event_JetMeas_nHCal_2_y = new TH1F("h_Event_JetMeas_nHCal_2_y", "Event with 2 jets measured in nHCal inelasticity y; y [1]; counts", 1000, 0.0, 1.0); + + h_Event_JetMeas_AllHCal_Q2 = new TH1F("h_Event_JetMeas_AllHCal_Q2", "Event with jets measured in any HCal Q^{2}; Q^{2} [GeV^{2}/c^{2}]; counts", 1000, 0.0, 10.0); + h_Event_JetMeas_AllHCal_x = new TH1F("h_Event_JetMeas_AllHCal_x", "Event with jets measured in any HCal x; x [1]; counts", nbins_x, logBinsArray_x); + h_Event_JetMeas_AllHCal_y = new TH1F("h_Event_JetMeas_AllHCal_y", "Event with jets measured in any HCal inelasticity y; y [1]; counts", 1000, 0.0, 1.0); + + + h_Event_JetMeas_no_nHCal_Q2 = new TH1F("h_Event_JetMeas_no_nHCal_Q2", "Event with jets measured without nHCal Q^{2}; Q^{2} [GeV^{2}/c^{2}]; counts", 1000, 0.0, 10.0); + h_Event_JetMeas_no_nHCal_x = new TH1F("h_Event_JetMeas_no_nHCal_x", "Event with jets measured without nHCal x; x [1]; counts", nbins_x, logBinsArray_x); + h_Event_JetMeas_no_nHCal_y = new TH1F("h_Event_JetMeas_no_nHCal_y", "Event with jets measured without nHCal inelasticity y; y [1]; counts", 1000, 0.0, 1.0); + + + + h_Event_HCal_jets = new TH2F("h_Event_HCal_jets", "Event with dijets in HCals; Jet #1; Jet #2; counts", 4, 0.0, 4.0, 4, 0.0, 4.0); + h_Event_HCal_jets_meas = new TH2F("h_Event_HCal_jets_meas", "Event with measured dijets in HCals; Jet #1; Jet #2; counts", 4, 0.0, 4.0, 4, 0.0, 4.0); + h_Event_HCal_jets_meas_no_nHCal = new TH2F("h_Event_HCal_jets_meas_no_nHCal", "Event with measured dijets in HCals without nHCal; Jet #1; Jet #2; counts", 4, 0.0, 4.0, 4, 0.0, 4.0); + + + h_Event_HCal_jets->GetXaxis()->SetBinLabel(1, "nHCal"); + h_Event_HCal_jets->GetXaxis()->SetBinLabel(2, "Barrel HCal"); + h_Event_HCal_jets->GetXaxis()->SetBinLabel(3, "LFHCAL"); + h_Event_HCal_jets->GetXaxis()->SetBinLabel(4, "Any"); + + h_Event_HCal_jets->GetYaxis()->SetBinLabel(1, "nHCal"); + h_Event_HCal_jets->GetYaxis()->SetBinLabel(2, "Barrel HCal"); + h_Event_HCal_jets->GetYaxis()->SetBinLabel(3, "LFHCAL"); + h_Event_HCal_jets->GetYaxis()->SetBinLabel(4, "Any"); + + + h_Event_HCal_jets_meas->GetXaxis()->SetBinLabel(1, "nHCal"); + h_Event_HCal_jets_meas->GetXaxis()->SetBinLabel(2, "Barrel HCal"); + h_Event_HCal_jets_meas->GetXaxis()->SetBinLabel(3, "LFHCAL"); + h_Event_HCal_jets_meas->GetXaxis()->SetBinLabel(4, "Any"); + + h_Event_HCal_jets_meas->GetYaxis()->SetBinLabel(1, "nHCal"); + h_Event_HCal_jets_meas->GetYaxis()->SetBinLabel(2, "Barrel HCal"); + h_Event_HCal_jets_meas->GetYaxis()->SetBinLabel(3, "LFHCAL"); + h_Event_HCal_jets_meas->GetYaxis()->SetBinLabel(4, "Any"); + + + h_Event_HCal_jets_meas_no_nHCal->GetXaxis()->SetBinLabel(1, "nHCal"); + h_Event_HCal_jets_meas_no_nHCal->GetXaxis()->SetBinLabel(2, "Barrel HCal"); + h_Event_HCal_jets_meas_no_nHCal->GetXaxis()->SetBinLabel(3, "LFHCAL"); + h_Event_HCal_jets_meas_no_nHCal->GetXaxis()->SetBinLabel(4, "Any"); + + h_Event_HCal_jets_meas_no_nHCal->GetYaxis()->SetBinLabel(1, "nHCal"); + h_Event_HCal_jets_meas_no_nHCal->GetYaxis()->SetBinLabel(2, "Barrel HCal"); + h_Event_HCal_jets_meas_no_nHCal->GetYaxis()->SetBinLabel(3, "LFHCAL"); + h_Event_HCal_jets_meas_no_nHCal->GetYaxis()->SetBinLabel(4, "Any"); + + + // full + h_Event_HCal_jets_meas_full = new TH2F("h_Event_HCal_jets_meas_full", "Event with measured dijets in HCals; Jet #1; Jet #2; counts", 6, 0.0, 6.0, 6, 0.0, 6.0); + + h_Event_HCal_jets_meas_full->GetXaxis()->SetBinLabel(1, "nHCal"); + h_Event_HCal_jets_meas_full->GetXaxis()->SetBinLabel(2, "nHCal+Barrel HCal"); + h_Event_HCal_jets_meas_full->GetXaxis()->SetBinLabel(3, "Barrel HCal"); + h_Event_HCal_jets_meas_full->GetXaxis()->SetBinLabel(4, "Barrel HCal+LFHCAL"); + h_Event_HCal_jets_meas_full->GetXaxis()->SetBinLabel(5, "LFHCAL"); + h_Event_HCal_jets_meas_full->GetXaxis()->SetBinLabel(6, "Any"); + + h_Event_HCal_jets_meas_full->GetYaxis()->SetBinLabel(1, "nHCal"); + h_Event_HCal_jets_meas_full->GetYaxis()->SetBinLabel(2, "nHCal+Barrel HCal"); + h_Event_HCal_jets_meas_full->GetYaxis()->SetBinLabel(3, "Barrel HCal"); + h_Event_HCal_jets_meas_full->GetYaxis()->SetBinLabel(4, "Barrel HCal+LFHCAL"); + h_Event_HCal_jets_meas_full->GetYaxis()->SetBinLabel(5, "LFHCAL"); + h_Event_HCal_jets_meas_full->GetYaxis()->SetBinLabel(6, "Any"); + + + + // Outgoing partons + h_Partons_status = new TH1F("h_Partons_status", "Outgoing parton status; status; counts", 401, -200.5, 200.5); + h_Partons_types = new TH1F("h_Partons_types", "Outgoing parton type; type; counts", 41, -0.5, 40.5); + h_Partons_types_anti = new TH1F("h_Partons_types_anti", "Outgoing parton type; type; counts", 41, -0.5, 40.5); + + h_Partons_types->GetXaxis()->SetBinLabel(h_Partons_types->GetXaxis()->FindBin(1), "d"); + h_Partons_types->GetXaxis()->SetBinLabel(h_Partons_types->GetXaxis()->FindBin(2), "u"); + h_Partons_types->GetXaxis()->SetBinLabel(h_Partons_types->GetXaxis()->FindBin(3), "s"); + h_Partons_types->GetXaxis()->SetBinLabel(h_Partons_types->GetXaxis()->FindBin(4), "c"); + h_Partons_types->GetXaxis()->SetBinLabel(h_Partons_types->GetXaxis()->FindBin(5), "b"); + h_Partons_types->GetXaxis()->SetBinLabel(h_Partons_types->GetXaxis()->FindBin(6), "t"); + //h_Partons_types->GetXaxis()->SetBinLabel(h_Partons_types->GetXaxis()->FindBin(7), "b^{#prime}"); + //h_Partons_types->GetXaxis()->SetBinLabel(h_Partons_types->GetXaxis()->FindBin(8), "t^{#prime}"); + + h_Partons_types->GetXaxis()->SetBinLabel(h_Partons_types->GetXaxis()->FindBin(11), "e^{-}"); + h_Partons_types->GetXaxis()->SetBinLabel(h_Partons_types->GetXaxis()->FindBin(12), "#nu_{e}"); + h_Partons_types->GetXaxis()->SetBinLabel(h_Partons_types->GetXaxis()->FindBin(13), "\mu^{-}"); + h_Partons_types->GetXaxis()->SetBinLabel(h_Partons_types->GetXaxis()->FindBin(14), "#nu_{#mu}"); + h_Partons_types->GetXaxis()->SetBinLabel(h_Partons_types->GetXaxis()->FindBin(15), "\tau^{-}"); + h_Partons_types->GetXaxis()->SetBinLabel(h_Partons_types->GetXaxis()->FindBin(16), "#nu_{#tau}"); + h_Partons_types->GetXaxis()->SetBinLabel(h_Partons_types->GetXaxis()->FindBin(17), "#tau^{#prime-}"); + h_Partons_types->GetXaxis()->SetBinLabel(h_Partons_types->GetXaxis()->FindBin(18), "#nu_{#tau^{#prime}}"); + + h_Partons_types->GetXaxis()->SetBinLabel(h_Partons_types->GetXaxis()->FindBin(21), "g"); + h_Partons_types->GetXaxis()->SetBinLabel(h_Partons_types->GetXaxis()->FindBin(22), "#gamma"); + h_Partons_types->GetXaxis()->SetBinLabel(h_Partons_types->GetXaxis()->FindBin(23), "Z^{0}"); + h_Partons_types->GetXaxis()->SetBinLabel(h_Partons_types->GetXaxis()->FindBin(24), "W^{+}"); + + + h_Partons_types_anti->GetXaxis()->SetBinLabel(h_Partons_types_anti->GetXaxis()->FindBin(1), "d"); + h_Partons_types_anti->GetXaxis()->SetBinLabel(h_Partons_types_anti->GetXaxis()->FindBin(2), "u"); + h_Partons_types_anti->GetXaxis()->SetBinLabel(h_Partons_types_anti->GetXaxis()->FindBin(3), "s"); + h_Partons_types_anti->GetXaxis()->SetBinLabel(h_Partons_types_anti->GetXaxis()->FindBin(4), "c"); + h_Partons_types_anti->GetXaxis()->SetBinLabel(h_Partons_types_anti->GetXaxis()->FindBin(5), "b"); + h_Partons_types_anti->GetXaxis()->SetBinLabel(h_Partons_types_anti->GetXaxis()->FindBin(6), "t"); + //h_Partons_types_anti->GetXaxis()->SetBinLabel(h_Partons_types_anti->GetXaxis()->FindBin(7), "b^{#prime}"); + //h_Partons_types_anti->GetXaxis()->SetBinLabel(h_Partons_types_anti->GetXaxis()->FindBin(8), "t^{#prime}"); + + h_Partons_types_anti->GetXaxis()->SetBinLabel(h_Partons_types_anti->GetXaxis()->FindBin(11), "e^{-}"); + h_Partons_types_anti->GetXaxis()->SetBinLabel(h_Partons_types_anti->GetXaxis()->FindBin(12), "#nu_{e}"); + h_Partons_types_anti->GetXaxis()->SetBinLabel(h_Partons_types_anti->GetXaxis()->FindBin(13), "#mu^{-}"); + h_Partons_types_anti->GetXaxis()->SetBinLabel(h_Partons_types_anti->GetXaxis()->FindBin(14), "#nu_{#mu}"); + h_Partons_types_anti->GetXaxis()->SetBinLabel(h_Partons_types_anti->GetXaxis()->FindBin(15), "#tau^{-}"); + h_Partons_types_anti->GetXaxis()->SetBinLabel(h_Partons_types_anti->GetXaxis()->FindBin(16), "#nu_{#tau}"); + h_Partons_types_anti->GetXaxis()->SetBinLabel(h_Partons_types_anti->GetXaxis()->FindBin(17), "#tau^{#prime-}"); + h_Partons_types_anti->GetXaxis()->SetBinLabel(h_Partons_types_anti->GetXaxis()->FindBin(18), "#nu_{#tau^{#prime}}"); + + h_Partons_types_anti->GetXaxis()->SetBinLabel(h_Partons_types_anti->GetXaxis()->FindBin(21), "g"); + h_Partons_types_anti->GetXaxis()->SetBinLabel(h_Partons_types_anti->GetXaxis()->FindBin(22), "#gamma"); + h_Partons_types_anti->GetXaxis()->SetBinLabel(h_Partons_types_anti->GetXaxis()->FindBin(23), "Z^{0}"); + h_Partons_types_anti->GetXaxis()->SetBinLabel(h_Partons_types_anti->GetXaxis()->FindBin(24), "W^{+}"); + + + h_Partons_eta = new TH2F("h_Partons_eta", "Outgoing partons #eta_{1} vs #eta_{2}; #eta_{1} [1]; #eta_{2} [1]; counts", 100, -5.0, 5.0, 100, -5.0, 5.0); + h_Partons_phi = new TH2F("h_Partons_phi", "Outgoing partons #phi_{1} vs #phi_{2}; #phi_{1} [1]; #phi_{2} [1]; counts", 200, -2.0*TMath::Pi(), 2.0*TMath::Pi(), 200, -2.0*TMath::Pi(), 2.0*TMath::Pi()); + + h_Partons_p = new TH2F("h_Partons_p", "Outgoing partons p_{1} vs p_{2}; p_{1} [GeV/c]; p_{2} [GeV/c]; counts", 200, 0.0, 20.0, 200, 0.0, 20.0); + h_Partons_pT = new TH2F("h_Partons_pT", "Outgoing partons p_{T}^{1} vs p_{T}^{2}; p_{T}^{1} [GeV/c]; p_{T}^{2} [GeV/c]; counts", 200, 0.0, 20.0, 200, 0.0, 20.0); + + h_Parton_y_eta1 = new TH2F("h_Parton_y_eta1", "Outgoing partons #eta_{1} vs. y; y [1]; #eta_{1} [1]; counts", 100, 0.0, 1.0, 100, -5.0, 5.0); + + h_Parton_eta_p = new TH2F("h_Parton_eta_p", "Outgoing parton p vs #eta; #eta [1]; p [GeV/c]; counts", 100, -5.0, 5.0, 200, 0.0, 20.0); + h_Parton_eta_pT = new TH2F("h_Parton_eta_pT", "Outgoing parton p_{T} vs #eta; #eta [1]; p_{T} [GeV/c]; counts", 100, -5.0, 5.0, 200, 0.0, 20.0); + h_Parton_eta_E = new TH2F("h_Parton_eta_E", "Outgoing parton E vs #eta; #eta [1]; E [GeV]; counts", 100, -5.0, 5.0, 200, 0.0, 20.0); + + h_Parton_x_eta = new TH2F("h_Parton_x_eta", "Outgoing parton #eta vs. x; x [1]; #eta [1]; counts", nbins_x, logBinsArray_x, 100, -5.0, 5.0); + h_Parton_y_eta = new TH2F("h_Parton_y_eta", "Outgoing partons #eta vs. y; y [1]; #eta [1]; counts", 100, 0.0, 1.0, 100, -5.0, 5.0); + + h_Parton_x_eta1 = new TH2F("h_Parton_x_eta1", "Outgoing partons #eta_{1} vs. x; x [1]; #eta_{1} [1]; counts", nbins_x, logBinsArray_x, 100, -5.0, 5.0); + h_Parton_x_eta2 = new TH2F("h_Parton_x_eta2", "Outgoing partons #eta_{2} vs. x; x [1]; #eta_{2} [1]; counts", nbins_x, logBinsArray_x, 100, -5.0, 5.0); + + h_Parton_y_eta1 = new TH2F("h_Parton_y_eta1", "Outgoing partons #eta_{1} vs. y; y [1]; #eta_{1} [1]; counts", 100, 0.0, 1.0, 100, -5.0, 5.0); + h_Parton_y_eta2 = new TH2F("h_Parton_y_eta2", "Outgoing partons #eta_{2} vs. y; y [1]; #eta_{2} [1]; counts", 100, 0.0, 1.0, 100, -5.0, 5.0); + + + // Particles + h_Particle_eta = new TH1F("h_Particle_eta", "MC particle #eta; #eta; counts", 200, -10.0, 10.0); + + h_Particle_eta_wE = new TH1F("h_Particle_eta_wE", "MC particle #eta, E weighed; #eta; counts", 200, -10.0, 10.0); + + h_Particle_eta_p = new TH2F("h_Particle_eta_p", "MC particles #eta vs. momentum; #eta [1]; p_{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_eta_pT = new TH2F("h_Particle_eta_pT", "MC particles #eta vs. momentum; #eta [1]; p_{T}^{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_eta_E = new TH2F("h_Particle_eta_E", "MC particles #eta vs. energy; #eta [1]; E_{MC} [GeV]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + + + // eta, momentum + h_Particle_pion_p_eta_p = new TH2F("h_Particle_Pion_p_eta_p", "MC particles #pi^{+} #eta vs. momentum; #eta [1]; p_{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_pion_n_eta_p = new TH2F("h_Particle_Pion_n_eta_p", "MC particles #pi^{-} #eta vs. momentum; #eta [1]; p_{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_Kaon_p_eta_p = new TH2F("h_Particle_Kaon_p_eta_p", "MC particles K^{+} #eta vs. momentum; #eta [1]; p_{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_Kaon_n_eta_p = new TH2F("h_Particle_Kaon_n_eta_p", "MC particles K^{-} #eta vs. momentum; #eta [1]; p_{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_proton_p_eta_p = new TH2F("h_Particle_Proton_p_eta_p", "MC particles p^{+} #eta vs. momentum; #eta [1]; p_{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_proton_n_eta_p = new TH2F("h_Particle_Proton_n_eta_p", "MC particles p^{-} #eta vs. momentum; #eta [1]; p_{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_Electron_p_eta_p = new TH2F("h_Particle_Electron_p_eta_p", "MC particles e^{+} #eta vs. momentum; #eta [1]; p_{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_Electron_n_eta_p = new TH2F("h_Particle_Electron_n_eta_p", "MC particles e^{-} #eta vs. momentum; #eta [1]; p_{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + + h_Particle_Neutron_eta_p = new TH2F("h_Particle_Neutron_eta_p", "MC particles n #eta vs. momentum; #eta [1]; p_{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_Gamma_eta_p = new TH2F("h_Particle_Gamma_eta_p", "MC particles #gamman #eta vs. momentum; #eta [1]; p_{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + + + // eta, transverse momentum pT + h_Particle_pion_p_eta_pT = new TH2F("h_Particle_Pion_p_eta_pT", "MC particles #pi^{+} #eta vs. momentum; #eta [1]; p_{T}^{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_pion_n_eta_pT = new TH2F("h_Particle_Pion_n_eta_pT", "MC particles #pi^{-} #eta vs. momentum; #eta [1]; p_{T}^{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_Kaon_p_eta_pT = new TH2F("h_Particle_Kaon_p_eta_pT", "MC particles K^{+} #eta vs. momentum; #eta [1]; p_{T}^{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_Kaon_n_eta_pT = new TH2F("h_Particle_Kaon_n_eta_pT", "MC particles K^{-} #eta vs. momentum; #eta [1]; p_{T}^{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_proton_p_eta_pT = new TH2F("h_Particle_Proton_p_eta_pT", "MC particles p^{+} #eta vs. momentum; #eta [1]; p_{T}^{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_proton_n_eta_pT = new TH2F("h_Particle_Proton_n_eta_pT", "MC particles p^{-} #eta vs. momentum; #eta [1]; p_{T}^{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_Electron_p_eta_pT = new TH2F("h_Particle_Electron_p_eta_pT", "MC particles e^{+} #eta vs. momentum; #eta [1]; p_{T}^{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_Electron_n_eta_pT = new TH2F("h_Particle_Electron_n_eta_pT", "MC particles e^{-} #eta vs. momentum; #eta [1]; p_{T}^{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + + h_Particle_Neutron_eta_pT = new TH2F("h_Particle_Neutron_eta_pT", "MC particles n #eta vs. momentum; #eta [1]; p_{T}^{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_Gamma_eta_pT = new TH2F("h_Particle_Gamma_eta_pT", "MC particles #gamman #eta vs. momentum; #eta [1]; p_{T}^{MC} [GeV/c]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + + + // eta, energy + h_Particle_Pion_p_eta_E = new TH2F("h_Particle_Pion_p_eta_E", "MC particles #pi^{+} #eta vs. energy; #eta [1]; E_{MC} [GeV]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_Pion_n_eta_E = new TH2F("h_Particle_Pion_n_eta_E", "MC particles #pi^{-} #eta vs. energy; #eta [1]; E_{MC} [GeV]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_Kaon_p_eta_E = new TH2F("h_Particle_Kaon_p_eta_E", "MC particles K^{+} #eta vs. energy; #eta [1]; E_{MC} [GeV]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_Kaon_n_eta_E = new TH2F("h_Particle_Kaon_n_eta_E", "MC particles K^{-} #eta vs. energy; #eta [1]; E_{MC} [GeV]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_Proton_p_eta_E = new TH2F("h_Particle_Proton_p_eta_E", "MC particles p^{+} #eta vs. energy; #eta [1]; E_{MC} [GeV]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_Proton_n_eta_E = new TH2F("h_Particle_Proton_n_eta_E", "MC particles p^{-} #eta vs. energy; #eta [1]; E_{MC} [GeV]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_Electron_p_eta_E = new TH2F("h_Particle_Electron_p_eta_E", "MC particles e^{+} #eta vs. energy; #eta [1]; E_{MC} [GeV]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_Electron_n_eta_E = new TH2F("h_Particle_Electron_n_eta_E", "MC particles e^{-} #eta vs. energy; #eta [1]; E_{MC} [GeV]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + + h_Particle_Neutron_eta_E = new TH2F("h_Particle_Neutron_eta_E", "MC particles n #eta vs. energy; #eta [1]; E_{MC} [GeV]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + h_Particle_Gamma_eta_E = new TH2F("h_Particle_Gamma_eta_E", "MC particles #gamma #eta vs. energy; #eta [1]; E_{MC} [GeV]; counts", 200, -10.0, 10.0, 500, 0.0, 50.0); + + // Jets + h_Jet_nPart = new TH1D("h_Jet_nPart", "Jet number of particles; N_{part} [1]; counts", 201, -0.5, 200.5); + h_Jet_mass = new TH1D("h_Jet_mass", "Jet mass; m [GeV/c^{2}]; counts", 2000, 0.0, 20.0); + h_Jet_charge = new TH1D("h_Jet_charge", "Jet charge; q [1]; counts", 101, -50.5, 50.5); + h_Jet_E = new TH1D("h_Jet_E", "Jet energy; E [GeV]; counts", 500, 0.0, 50.0); + h_Jet_p = new TH1D("h_Jet_p", "Jet momentum; p [GeV/c]; counts", 500, 0.0, 50.0); + h_Jet_pT = new TH1D("h_Jet_pT", "Jet transverse momentum; p_{T} [GeV/c]; counts", 500, 0.0, 50.0); + h_Jet_eta = new TH1D("h_Jet_eta", "Jet #eta; #eta [1]; counts", 200, -5.0, 5.0); + h_Jet_deta = new TH1D("h_Jet_deta", "Jet #Delta#eta; #Delta#eta [1]; counts", 400, -10.0, 10.0); + + h_Jets_eta = new TH2F("h_Jets_eta", "Jets #eta_{1} vs. #eta_{2}; #eta_{1} [1]; #eta_{2} [1]; counts", 200, -5.0, 5.0, 200, -5.0, 5.0); + h_Jets_phi = new TH2F("h_Jets_phi", "Jets #phi_{1} vs. #phi_{2}; #phi_{1} [1]; #phi_{2} [1]; counts", 200, -2.0*TMath::Pi(), 2.0*TMath::Pi(), 200, -2.0*TMath::Pi(), 2.0*TMath::Pi()); + h_Jets_p = new TH2F("h_Jets_p", "Jets p_{1} vs. p_{2}; p_{1} [GeV/c]; p_{2} [GeV/c]; counts", 200, 0.0, 20.0, 200, 0.0, 20.0); + h_Jets_pT = new TH2F("h_Jets_pT", "Jets p_{T}^{1} vs. p_{T}^{2}; p_{T}^{1} [GeV/c]; p_{T}^{2} [GeV/c]; counts", 200, 0.0, 20.0, 200, 0.0, 20.0); + h_Jets_E = new TH2F("h_Jets_E", "Jets E_{1} vs. E_{2}; E_{1} [GeV/c]; E_{2} [GeV]; counts", 200, 0.0, 20.0, 200, 0.0, 20.0); + + + // Jets measured + h_Jet_meas_nPart = new TH1D("h_Jet_meas_nPart", "Jet measured number of particles; N_{part} [1]; counts", 201, -0.5, 200.5); + h_Jet_meas_mass = new TH1D("h_Jet_meas_mass", "Jet measured mass; m [GeV/c^{2}]; counts", 2000, 0.0, 20.0); + h_Jet_meas_charge = new TH1D("h_Jet_meas_charge", "Jet measured charge; q [1]; counts", 101, -50.5, 50.5); + h_Jet_meas_E = new TH1D("h_Jet_meas_E", "Jet measured energy; E [GeV]; counts", 500, 0.0, 50.0); + h_Jet_meas_p = new TH1D("h_Jet_meas_p", "Jet measured momentum; p [GeV/c]; counts", 500, 0.0, 50.0); + h_Jet_meas_pT = new TH1D("h_Jet_meas_pT", "Jet measured transverse momentum; p_{T} [GeV/c]; counts", 500, 0.0, 50.0); + h_Jet_meas_eta = new TH1D("h_Jet_meas_eta", "Jet measured #eta; #eta [1]; counts", 200, -5.0, 5.0); + h_Jet_meas_deta = new TH1D("h_Jet_meas_deta", "Jet measured #Delta#eta; #Delta#eta [1]; counts", 400, -10.0, 10.0); + + h_Jets_meas_eta = new TH2F("h_Jets_meas_eta", "Jets measured #eta_{1} vs. #eta_{2}; #eta_{1} [1]; #eta_{2} [1]; counts", 200, -5.0, 5.0, 200, -5.0, 5.0); + h_Jets_meas_phi = new TH2F("h_Jets_meas_phi", "Jets measured #phi_{1} vs. #phi_{2}; #phi_{1} [1]; #phi_{2} [1]; counts", 200, -2.0*TMath::Pi(), 2.0*TMath::Pi(), 200, -2.0*TMath::Pi(), 2.0*TMath::Pi()); + h_Jets_meas_p = new TH2F("h_Jets_meas_p", "Jets measured p_{1} vs. p_{2}; p_{1} [GeV/c]; p_{2} [GeV/c]; counts", 200, 0.0, 20.0, 200, 0.0, 20.0); + h_Jets_meas_pT = new TH2F("h_Jets_meas_pT", "Jets measured p_{T}^{1} vs. p_{T}^{2}; p_{T}^{1} [GeV/c]; p_{T}^{2} [GeV/c]; counts", 200, 0.0, 20.0, 200, 0.0, 20.0); + h_Jets_meas_E = new TH2F("h_Jets_meas_E", "Jets measured E_{1} vs. E_{2}; E_{1} [GeV/c]; E_{2} [GeV]; counts", 200, 0.0, 20.0, 200, 0.0, 20.0); + + + h_Jet_bHCal_part_eta = new TH1D("h_Jet_bHCal_part_eta", "Jet in bHCal particle #eta; #eta [1]; counts", 200, -5.0, 5.0); + h_Jet_meas_bHCal_part_eta = new TH1D("h_Jet_meas_bHCal_part_eta", "Jet measured in bHCal particle #eta; #eta [1]; counts", 200, -5.0, 5.0); + + h_Jet_HCal_part_eta = new TH2D("h_Jet_HCal_part_eta", "Jet in HCals particle #eta; #eta [1]; HCal; counts", 200, -5.0, 5.0, 4, 0.0, 4.0); + h_Jet_meas_HCal_part_eta = new TH2D("h_Jet_meas_HCal_part_eta", "Jet measured in HCals particle #eta; #eta [1]; HCal; counts", 200, -5.0, 5.0, 4, 0.0, 4.0); + + h_Jet_HCal_part_eta->GetYaxis()->SetBinLabel(1, "no HCal"); + h_Jet_HCal_part_eta->GetYaxis()->SetBinLabel(2, "nHCal"); + h_Jet_HCal_part_eta->GetYaxis()->SetBinLabel(3, "bHCal"); + h_Jet_HCal_part_eta->GetYaxis()->SetBinLabel(4, "LFHCAL"); + + h_Jet_meas_HCal_part_eta->GetYaxis()->SetBinLabel(1, "no HCal"); + h_Jet_meas_HCal_part_eta->GetYaxis()->SetBinLabel(2, "nHCal"); + h_Jet_meas_HCal_part_eta->GetYaxis()->SetBinLabel(3, "bHCal"); + h_Jet_meas_HCal_part_eta->GetYaxis()->SetBinLabel(4, "LFHCAL"); + + + + // Jets measured without nHCal + h_Jet_meas_no_nHCal_nPart = new TH1D("h_Jet_meas_no_nHCal_nPart", "Jet measured without nHCal number of particles; N_{part} [1]; counts", 201, -0.5, 200.5); + h_Jet_meas_no_nHCal_mass = new TH1D("h_Jet_meas_no_nHCal_mass", "Jet measured without nHCal mass; m [GeV/c^{2}]; counts", 2000, 0.0, 20.0); + h_Jet_meas_no_nHCal_charge = new TH1D("h_Jet_meas_no_nHCal_charge", "Jet measured without nHCal charge; q [1]; counts", 101, -50.5, 50.5); + h_Jet_meas_no_nHCal_E = new TH1D("h_Jet_meas_no_nHCal_E", "Jet measured without nHCal energy; E [GeV]; counts", 500, 0.0, 50.0); + h_Jet_meas_no_nHCal_p = new TH1D("h_Jet_meas_no_nHCal_p", "Jet measured without nHCal momentum; p [GeV/c]; counts", 500, 0.0, 50.0); + h_Jet_meas_no_nHCal_pT = new TH1D("h_Jet_meas_no_nHCal_pT", "Jet measured without nHCal transverse momentum; p_{T} [GeV/c]; counts", 500, 0.0, 50.0); + h_Jet_meas_no_nHCal_eta = new TH1D("h_Jet_meas_no_nHCal_eta", "Jet measured without nHCal #eta; #eta [1]; counts", 200, -5.0, 5.0); + h_Jet_meas_no_nHCal_deta = new TH1D("h_Jet_meas_no_nHCal_deta", "Jet measured without nHCal #Delta#eta; #Delta#eta [1]; counts", 400, -10.0, 10.0); + + h_Jets_meas_no_nHCal_eta = new TH2F("h_Jets_meas_no_nHCal_eta", "Jets measured without nHCal #eta_{1} vs. #eta_{2}; #eta_{1} [1]; #eta_{2} [1]; counts", 200, -5.0, 5.0, 200, -5.0, 5.0); + h_Jets_meas_no_nHCal_phi = new TH2F("h_Jets_meas_no_nHCal_phi", "Jets measured without nHCal #phi_{1} vs. #phi_{2}; #phi_{1} [1]; #phi_{2} [1]; counts", 200, -2.0*TMath::Pi(), 2.0*TMath::Pi(), 200, -2.0*TMath::Pi(), 2.0*TMath::Pi()); + h_Jets_meas_no_nHCal_p = new TH2F("h_Jets_meas_no_nHCal_p", "Jets measured without nHCal p_{1} vs. p_{2}; p_{1} [GeV/c]; p_{2} [GeV/c]; counts", 200, 0.0, 20.0, 200, 0.0, 20.0); + h_Jets_meas_no_nHCal_pT = new TH2F("h_Jets_meas_no_nHCal_pT", "Jets measured without nHCal p_{T}^{1} vs. p_{T}^{2}; p_{T}^{1} [GeV/c]; p_{T}^{2} [GeV/c]; counts", 200, 0.0, 20.0, 200, 0.0, 20.0); + h_Jets_meas_no_nHCal_E = new TH2F("h_Jets_meas_no_nHCal_E", "Jets measured without nHCal E_{1} vs. E_{2}; E_{1} [GeV/c]; E_{2} [GeV]; counts", 200, 0.0, 20.0, 200, 0.0, 20.0); + + + // measured jets vs. partons + h_Jets_meas_Partons_eta = new TH2F("h_Jets_meas_Partons_eta", "Jet #eta vs parton #eta; #eta_{parton} [1]; #eta_{jet} [1]; counts", 100, -5.0, 5.0, 100, -5.0, 5.0); + h_Jets_meas_Partons_phi = new TH2F("h_Jets_meas_Partons_phi", "Jet #phi vs parton #phi; #phi_{parton} [1]; #phi_{jet} [1]; counts", 200, -2.0*TMath::Pi(), 2.0*TMath::Pi(), 200, -2.0*TMath::Pi(), 2.0*TMath::Pi()); + h_Jets_meas_Partons_E = new TH2F("h_Jets_meas_Partons_E", "Jet E vs parton E; E_{parton} [GeV]; E_{jet} [GeV]; counts", 100, 0.0, 50.0, 100, 0.0, 50.0); + + h_Jet_meas_Parton_eta1 = new TH2F("h_Jet_meas_Parton_eta1", "Jet #eta vs parton 1 #eta; #eta_{parton} [1]; #eta_{jet} [1]; counts", 100, -5.0, 5.0, 100, -5.0, 5.0); + h_Jet_meas_Parton_phi1 = new TH2F("h_Jet_meas_Parton_phi1", "Jet #phi vs parton 1 #phi; #phi_{parton} [1]; #phi_{jet} [1]; counts", 200, -2.0*TMath::Pi(), 2.0*TMath::Pi(), 200, -2.0*TMath::Pi(), 2.0*TMath::Pi()); + h_Jet_meas_Parton_E1 = new TH2F("h_Jet_meas_Parton_E1", "Jet E vs parton 1 E; E_{parton} [GeV]; E_{jet} [GeV]; counts", 100, 0.0, 50.0, 100, 0.0, 50.0); + + h_Jet_meas_Parton_eta2 = new TH2F("h_Jet_meas_Parton_eta2", "Jet #eta vs parton 2 #eta; #eta_{parton} [1]; #eta_{jet} [1]; counts", 100, -5.0, 5.0, 100, -5.0, 5.0); + h_Jet_meas_Parton_phi2 = new TH2F("h_Jet_meas_Parton_phi2", "Jet #phi vs parton 2 #phi; #phi_{parton} [1]; #phi_{jet} [1]; counts", 200, -2.0*TMath::Pi(), 2.0*TMath::Pi(), 200, -2.0*TMath::Pi(), 2.0*TMath::Pi()); + h_Jet_meas_Parton_E2 = new TH2F("h_Jet_meas_Parton_E2", "Jet E vs parton 2 E; E_{parton} [GeV]; E_{jet} [GeV]; counts", 100, 0.0, 50.0, 100, 0.0, 50.0); + + + // temp + const int nEtaBins = 7; + double EtaBins[nEtaBins+1] = {-10, -4.14, -1.18, -1.1, 1.1, 1.2, 4.2, 10}; + + hist_eta_energy_tmp = new TH1F("hist_eta_energy_tmp", "hist_eta_energy_tmp; #eta, E [GeV]", nEtaBins, EtaBins); + hist_eta_energy_denom_tmp = new TH1F("hist_eta_energy_denom_tmp", "hist_eta_energy_denom_tmp; #eta, E [GeV]", nEtaBins, EtaBins); + + + //hist_eta_energy_tmp->Sumw2(); + //hist_eta_energy_denom_tmp->Sumw2(); + + + return 1; +} + + +#endif /* HISTOGRAMSPYTHIA_H_ */ diff --git a/benchmarks/nhcal_dimuon_photoproduction/scripts/LinkDef.h b/benchmarks/nhcal_dimuon_photoproduction/scripts/LinkDef.h new file mode 100755 index 00000000..2502d63b --- /dev/null +++ b/benchmarks/nhcal_dimuon_photoproduction/scripts/LinkDef.h @@ -0,0 +1,19 @@ +// main144LinkDef.h is a part of the PYTHIA event generator. +// Copyright (C) 2024 Torbjorn Sjostrand. +// PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details. +// Please respect the MCnet Guidelines, see GUIDELINES for details. + +// Header used to generate a ROOT dictionary for the PYTHIA classes. +// Modified by Rene Brun and Axel Naumann to put the Pythia::event +// into a TTree. + +#ifdef __CINT__ +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; +#pragma link C++ namespace Pythia8; +#pragma link C++ class Pythia8::Event+; +#pragma link C++ class Pythia8::Particle+; +#pragma link C++ class Pythia8::Vec4+; +#pragma link C++ class PythiaEvent+; +#endif diff --git a/benchmarks/nhcal_dimuon_photoproduction/scripts/Makefile b/benchmarks/nhcal_dimuon_photoproduction/scripts/Makefile new file mode 100755 index 00000000..b0faa7b2 --- /dev/null +++ b/benchmarks/nhcal_dimuon_photoproduction/scripts/Makefile @@ -0,0 +1,40 @@ +#=============================================================================== +# File: Makefile +# +# Example Makefile for Pythia8 written for HF in STAR +# Author: Thomas Ullrich +# Last modified: September 9, 2008 +# +# This needs 3 packages installed: Pythia8, LHAPDF, and ROOT +# +# This is setup for Mac OS X 10.5.4 but should be easy to +# adapt for other (Unix) platforms. In principle changing the +# program name (PROGRAM), the location of PYTHIA (PYTHIAPATH), +# and the LHAPDF libraries (LHAPDFPATH) should do the job. +# Note that the environment variable ROOTSYS needs to be set. +# Otherwise define it here in the makefile. +#=============================================================================== + +PROGRAM = pythiaDiMuon +PYTHIAPATH = /opt/local +LHAPDFPATH = $(shell lhapdf-config --prefix) +ROOTSYS = $(shell root-config --prefix) + +FASTJETCXX = `fastjet-config --cxxflags` +FASTJETLIB = `fastjet-config --libs --plugins` + +DICTHEADERS = PythiaEvent.h LinkDef.h + +ROOTLIB = `root-config --libdir` + +CXX = g++ +CXXFLAGS = -O -W -Wall -m64 -std=c++20 $(FASTJETCXX) +CPPFLAGS = -I$(PYTHIAPATH)/include -I$(ROOTSYS)/include/root -I$(LHAPDFPATH)/include + +LDFLAGS = -L$(PYTHIAPATH)/lib -L$(ROOTLIB) -L$(LHAPDFPATH)/lib -lpythia8 -lHepMC3 -lstdc++ -lCore -lThread -lRIO -lNet -lHist -lMathCore -lGraf -lGraf3d -lGpad -lTree -lRint -lPostscript -lMatrix -lPhysics -lfreetype -lpthread -lm -ldl + +$(PROGRAM): $(PROGRAM).cpp DictOutput.cxx libPythiaEvent.so + $(CXX) $(CXXFLAGS) $(CPPFLAGS) $(PROGRAM).cpp DictOutput.cxx $(LDFLAGS) $(FASTJETLIB) libPythiaEvent.so -o $(PROGRAM) + +DictOutput.cxx : PythiaEvent.h + rootcint -f $@ -c $(DICTHEADERS) \ No newline at end of file diff --git a/benchmarks/nhcal_dimuon_photoproduction/scripts/PythiaEvent.h b/benchmarks/nhcal_dimuon_photoproduction/scripts/PythiaEvent.h new file mode 100755 index 00000000..c26368fe --- /dev/null +++ b/benchmarks/nhcal_dimuon_photoproduction/scripts/PythiaEvent.h @@ -0,0 +1,73 @@ +/* + * PythiaEvent.h + * + * Created on: 24 lis 2020 + * Author: Khaless + */ + +#ifndef PythiaEvent_H +#define PythiaEvent_H + +#include + +// ROOT +#include "TObject.h" + +#include "Pythia8/Pythia.h" +#include "Pythia8/Event.h" + +using namespace Pythia8; + +class PythiaEvent : public TObject { +public: + PythiaEvent() {} + //PythiaEvent(const PythiaEvent& ev); + ~PythiaEvent() {} + + // event + long eventId; + int nParticlesFinal; + float Q2; // momentum transfer + float x; // momentum fraction + float y; // inelasticity + + Particle scatteredEle; + + vector particles; + + +/* + // Upsilon candidate + float p; + float pt; + float y; + float eta; + float phi; + float m; + short charge; + + // electron 1 + float p1; + float pt1; + float eta1; + float phi1; + float nSigmaElectron1; + float nSigmaPion1; + + // electron 2 + float p2; + float pt2; + float eta2; + float phi2; + float nSigmaElectron2; + float nSigmaPion2;*/ + +private: + ClassDef(PythiaEvent,1); +}; + + + +#endif /* PythiaEvent_H */ + + diff --git a/benchmarks/nhcal_dimuon_photoproduction/scripts/dimuon_photoproduction_analysis.cxx b/benchmarks/nhcal_dimuon_photoproduction/scripts/dimuon_photoproduction_analysis.cxx new file mode 100644 index 00000000..c8cd76e5 --- /dev/null +++ b/benchmarks/nhcal_dimuon_photoproduction/scripts/dimuon_photoproduction_analysis.cxx @@ -0,0 +1,910 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include "TLorentzVector.h" + +#include "ROOT/RDataFrame.hxx" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "TROOT.h" +#include "TRandom.h" +#include "TH3.h" + +#include "DD4hep/Detector.h" +#include "DDRec/CellIDPositionConverter.h" +#include "DD4hep/Volumes.h" +#include "DD4hep/Objects.h" +#include "DD4hep/Shapes.h" +#include +#include +#include +#include +#include + +#include "podio/Frame.h" +#include "podio/CollectionBase.h" +#include "podio/ROOTReader.h" +#include "podio/CollectionIDTable.h" +#include "podio/ObjectID.h" + +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/MCParticle.h" + +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4hep/SimCalorimeterHit.h" + +#include "edm4hep/CalorimeterHit.h" +#include "edm4hep/CalorimeterHitCollection.h" + +#include "edm4eic/ClusterCollection.h" +#include "edm4eic/Cluster.h" +#include "edm4eic/ClusterData.h" + +#include "edm4eic/CalorimeterHit.h" +#include "edm4eic/CalorimeterHitCollection.h" + +#include "edm4eic/InclusiveKinematicsCollection.h" +#include "edm4eic/InclusiveKinematics.h" +#include "edm4hep/utils/kinematics.h" +#include "edm4hep/utils/vector_utils.h" +#include "edm4eic/vector_utils_legacy.h" + +#include "edm4eic/Track.h" +#include "edm4eic/TrackSegment.h" +#include "edm4eic/TrackPoint.h" +#include "edm4eic/TrackParameters.h" +#include "edm4eic/TrackSegmentCollection.h" + +#include "edm4eic/ReconstructedParticleCollection.h" +#include "edm4eic/ReconstructedParticle.h" + +#include "edm4eic/MCRecoCalorimeterHitAssociationCollection.h" +#include "edm4eic/MCRecoCalorimeterHitAssociation.h" +#include "edm4eic/MCRecoParticleAssociationCollection.h" +#include "edm4eic/MCRecoParticleAssociation.h" + +using namespace std; +using namespace ROOT; +using namespace TMath; +using namespace edm4hep; + +dd4hep::Detector* det = NULL; + +constexpr double ETA_MIN = -4.16, ETA_MAX = -1.16; + +inline bool inNHCal(double eta) {return (eta >= ETA_MIN && eta <= ETA_MAX);} + +inline string addPrefixAfterSlash(const string& path, + const string& prefix) { + const auto slash = path.find_last_of('/'); + if (slash == string::npos) + return prefix + path; + return path.substr(0, slash + 1) + prefix + path.substr(slash + 1); +} + +void MakeLogBins(double *array, int nbins, double binLo, double binHi) +{ + double logMin = log10(binLo); + double logMax = log10(binHi); + double binWidth = (logMax - logMin) / nbins; + + for (int i = 0; i <= nbins; ++i) { + array[i] = pow(10, logMin + i * binWidth); + } +} + +TLine* makeLine(double x1, double y1, double x2, double y2) { + TLine* l = new TLine(x1, y1, x2, y2); + l->SetLineColor(kRed); + l->SetLineStyle(2); + l->SetLineWidth(2); + l->Draw("same"); + return l; +} + +static TGraph* makeEffGraph_B(TH1D* h_sel, TH1D* h_all, + const char* name, + int minAll = 2, + bool logx = false) +{ + if (!h_sel || !h_all) return nullptr; + + std::unique_ptr h_eff((TH1D*)h_sel->Clone(Form("tmp_%s", name ? name : "eff"))); + h_eff->SetDirectory(nullptr); + h_eff->Sumw2(); + h_eff->Divide(h_sel, h_all, 1.0, 1.0, "B"); + + const int nb = h_all->GetNbinsX(); + + std::unique_ptr g_log; + std::unique_ptr g_lin; + + if (logx) g_log = std::make_unique(); + else g_lin = std::make_unique(); + + int ip = 0; + for (int b = 1; b <= nb; ++b) { + const double all = h_all->GetBinContent(b); + if (all < minAll) continue; + + const double xl = h_all->GetXaxis()->GetBinLowEdge(b); + const double xh = h_all->GetXaxis()->GetBinUpEdge(b); + + const double y = 100.0 * h_eff->GetBinContent(b); + const double ey = 100.0 * h_eff->GetBinError(b); + + if (logx) { + if (xl <= 0.0 || xh <= 0.0) continue; + const double x = sqrt(xl * xh); + const double exl = x - xl; + const double exh = xh - x; + + g_log->SetPoint(ip, x, y); + g_log->SetPointError(ip, exl, exh, ey, ey); + } else { + const double x = h_all->GetBinCenter(b); + const double ex = 0.5 * (xh - xl); + + g_lin->SetPoint(ip, x, y); + g_lin->SetPointError(ip, ex, ey); + } + ++ip; + } + + TGraph* g = logx ? (TGraph*)g_log.release() : (TGraph*)g_lin.release(); + g->SetName(Form("g_%s", name ? name : "eff")); + return g; +} + +inline void drawEffPanel(TH1D* h_all, + TH1D* const h_in[3], + const char* title, + const char* xlabel, + bool logx) +{ + gPad->SetGrid(); + if (logx) gPad->SetLogx(); + + auto* mg = new TMultiGraph(); + auto* leg = new TLegend(0.63, 0.7, 0.88, 0.88); + leg->SetBorderSize(0); + leg->SetFillStyle(0); + leg->SetTextSize(0.04); + + Color_t colors[3] = {kBlue, kRed, kBlack}; + Style_t markers[3] = {20, 21, 22}; + + int added = 0; + + for (int i = 0; i < 3; ++i) { + if (!h_all || !h_in[i]) continue; + auto* g = makeEffGraph_B(h_in[i], h_all, Form("eff_%d", i), 2, logx); + if (!g || g->GetN() == 0) continue; + + g->SetMarkerColor(colors[i]); + g->SetMarkerStyle(markers[i]); + g->SetMarkerSize(1.0); + g->SetLineColor(kBlack); + + mg->Add(g, "PE"); + leg->AddEntry(g, Form("%d mu-in nHCAL", i), "p"); + ++added; + } + + if (h_all && h_in[0] && h_in[1] && h_in[2]) { + std::unique_ptr h_sum((TH1D*)h_in[0]->Clone("h_sum")); + h_sum->SetDirectory(nullptr); + h_sum->Add(h_in[1]); + h_sum->Add(h_in[2]); + + auto* gsum = makeEffGraph_B(h_sum.get(), h_all, "eff_sum", 2, logx); + if (gsum && gsum->GetN() > 0) { + gsum->SetMarkerStyle(25); + gsum->SetMarkerColor(kMagenta + 1); + gsum->SetFillColor(kMagenta + 1); + gsum->SetLineColor(kBlack); + gsum->SetMarkerSize(1.0); + + mg->Add(gsum, "PE"); + leg->AddEntry(gsum, "All eligible", "p"); + ++added; + } + } + + mg->SetTitle(Form("%s;%s;geom. acc. [%%]", title ? title : "", xlabel ? xlabel : "")); + mg->Draw("A"); + + gPad->Update(); + if (mg->GetHistogram()) { + double ymax = mg->GetHistogram()->GetMaximum(); + mg->GetHistogram()->SetMaximum(ymax * 1.35); + } + gPad->Modified(); + gPad->Update(); + + if (logx) { + auto* ax = mg->GetXaxis(); + if (ax) { + double xmin = ax->GetXmin(), xmax = ax->GetXmax(); + if (xmin <= 0) xmin = 1e-12; + mg->GetXaxis()->SetLimits(xmin, xmax); + } + } + + leg->Draw(); +} + +inline TVector3 getPlasticDimensionsCM(dd4hep::Detector& det, + const dd4hep::DDSegmentation::BitFieldCoder* dec, + dd4hep::DDSegmentation::CellID cid, + int slice_idx, + int plastic_slice_value = 3) +{ + const double NaN = numeric_limits::quiet_NaN(); + TVector3 dims(NaN, NaN, NaN); + try { + if (!dec) throw runtime_error("decoder==nullptr"); + if (dec->get(cid, slice_idx) != plastic_slice_value) + throw runtime_error("cell is not plastic"); + + auto de = det.volumeManager().lookupDetElement(cid); + if (!de.isValid()) throw runtime_error("lookupDetElement failed"); + + auto pv = de.placement(); + if (!pv.isValid()) throw runtime_error("DetElement has no placement"); + + auto vol = pv.volume(); + if (!vol.isValid()) throw runtime_error("Invalid Volume"); + + auto* box = dynamic_cast(vol.solid().ptr()); + if (!box) throw runtime_error("Solid is not TGeoBBox"); + + dims.SetXYZ(2.0 * box->GetDX(), + 2.0 * box->GetDY(), + 2.0 * box->GetDZ()); + + //const double dz_cm = box->GetDZ(); + //const double thickness_cm = 2.0 * dz_cm; + //return thickness_cm; + return dims; + } catch (const exception& e) { + cerr << "[WARN] getPlasticThicknessMM: " << e.what() << " (cellID=" << cid << ")\n"; + return dims; + } +} + +inline double getPlasticCenterZ_cm( + const dd4hep::DDSegmentation::BitFieldCoder* decoder, + const dd4hep::rec::CellIDPositionConverter& cellid_converter, + dd4hep::DDSegmentation::CellID cid, + int slice_index, + int plastic_slice_value = 3) +{ + try { + if (!decoder) throw std::runtime_error("decoder==nullptr"); + if (decoder->get(cid, slice_index) != plastic_slice_value) + throw std::runtime_error("cell is not plastic (slice mismatch)"); + return cellid_converter.position(cid).z(); + } catch (const std::exception& e) { + cerr << "[WARN] getPlasticDimensionsCM: " << e.what() + << " (cellID=" << cid << ")\n"; + return std::numeric_limits::quiet_NaN();; + } +} + +struct GeomState { + double x, y, z; +}; + +inline double hypot2(double a, double b){ return sqrt(a*a + b*b); } +inline int sgnD(double v){ return (v>0) - (v<0); } + +inline pair +trackXYatZ_fromTwoStates(double x1,double y1,double z1, + double x2,double y2,double z2, + double zTarget) +{ + constexpr double EPSZ = 1e-12; + + const double dz = z2 - z1; + if (abs(dz) < EPSZ) { + const double d1 = abs(zTarget - z1); + const double d2 = abs(zTarget - z2); + if (d1 < d2) return {x1, y1}; + if (d2 < d1) return {x2, y2}; + return {0.5*(x1+x2), 0.5*(y1+y2)}; + } + + const double t = (zTarget - z1) / dz; + const double x = x1 + t * (x2 - x1); + const double y = y1 + t * (y2 - y1); + return {x, y}; +} + +inline pair +trackXYatZ(const GeomState& A, const GeomState& B, double zTarget){ + return trackXYatZ_fromTwoStates( + A.x, A.y, A.z, B.x, B.y, B.z, + zTarget + ); +} + + +int dimuon_photoproduction_analysis(const string& filename, string outname_pdf, string outname_png, TString compact_file) { + + gStyle->SetOptStat(0); + podio::ROOTReader reader; + reader.openFile(filename); + unsigned nEvents = reader.getEntries("events"); + cout << "Number of events: " << nEvents << endl; + + det = &(dd4hep::Detector::getInstance()); + det->fromCompact(compact_file.Data()); + det->volumeManager(); + det->apply("DD4hepVolumeManager", 0, 0); + + dd4hep::rec::CellIDPositionConverter cellid_converter(*det); + auto idSpec = det->readout("HcalEndcapNHits").idSpec(); + auto decoder = idSpec.decoder(); + const int slice_index = decoder->index("slice"); + if (slice_index < 0) { + cerr << "ERROR: 'slice' field not found in cell ID spec!" << endl; + return 1; + } + + //double thickness_cm = det->constantAsDouble("HcalEndcapNPolystyreneThickness")* 0.1; + + constexpr int NBINS = 60; + + constexpr double Q2_MIN = 1e-2; + constexpr double Q2_MAX = 1.0; + constexpr double X_MIN = 1e-6; + constexpr double X_MAX = 1e-3; + constexpr double Y_MIN = 0.0; + constexpr double Y_MAX = 1.0; + constexpr double W_MIN = 0.0; + constexpr double W_MAX = 160.0; + + constexpr int Z_NBINS = 100; + constexpr double Z_MIN_MM = -4500.0; + constexpr double Z_MAX_MM = -3600.0; + + constexpr int DXY_NBINS = 120; + constexpr double DXY_MIN_MM = -1000.0; + constexpr double DXY_MAX_MM = 1000.0; + + constexpr int DR_NBINS = 120; + constexpr double DR_MIN_MM = 0.0; + constexpr double DR_MAX_MM = 400.0; + + constexpr double P_MIN_GEV = 0.0; + constexpr double P_MAX_GEV = 25.0; + constexpr int P_NBINS = 50; + + constexpr double E_MIN_GEV = 2e-2; + constexpr double E_MAX_GEV = 10.0; + + constexpr int SIZE = 3; + constexpr double MIP_ENERGY_GEV = 0.002; + constexpr double E_CUTS[SIZE] = {1.5, 1.7, 2.0}; + constexpr double E_THRESH = E_CUTS[2]; + constexpr double LAYER_PROC[SIZE] = {0.6, 0.7, 0.8}; + + double t_cm; + double DR_CUTS_CM[SIZE] = {7.0, 10.0, 13.0}; + double DR_THRESH_MM = 10*DR_CUTS_CM[2]; + int LAYER_MAX = 10; + int LAYER_CUTS[SIZE] = {5, 6, 7}; + int LAYER_THRESH = LAYER_CUTS[0]; + + vector Xbins(NBINS+1), Q2bins(NBINS+1), Ebins(NBINS+1); + MakeLogBins(Q2bins.data(), NBINS, Q2_MIN, Q2_MAX); + MakeLogBins(Xbins.data(), NBINS, X_MIN, X_MAX); + MakeLogBins(Ebins.data(), NBINS, E_MIN_GEV, E_MAX_GEV); + + gStyle->SetTitleSize(0.06, "XYZ"); + gStyle->SetLabelSize(0.06, "XYZ"); + gStyle->SetPadLeftMargin(0.15); + gStyle->SetPadRightMargin(0.15); + gStyle->SetPadBottomMargin(0.15); + gStyle->SetPadTopMargin(0.10); + gROOT->ForceStyle(); + + TH2D* hEtaPt = new TH2D("hEtaPt", "Muon #eta vs p_{T};#eta;p_{T} [GeV]", 100, -6., 6., 100, 0., 7.); + TH2D* hx_Q2 = new TH2D("hx_Q2","Muon x vs Q^{2}; x; Q^{2} [GeV^{2}]", NBINS, Xbins.data(), NBINS, Q2bins.data()); + TH2D* hEta1_Eta2 = new TH2D("hEta1_Eta2", "Muon #eta_{+} vs #eta_{-}; #eta_{+} (PDG=-13); #eta_{-} (PDG=13)", 100, -6., 6., 100, -6., 6.); + + TH1D* hQ2_all = new TH1D("hQ2_all", "Q^{2}", NBINS, Q2bins.data()); + TH1D* hX_all = new TH1D("hX_all", "x", NBINS, Xbins.data()); + TH1D* hY_all = new TH1D("hY_all", "y", NBINS, Y_MIN, Y_MAX); + TH1D* hW_all = new TH1D("hW_all", "W", NBINS, W_MIN, W_MAX); + + TH1D* hQ2_in[3], *hX_in[3], *hY_in[3], *hW_in[3]; + for (int i = 0; i < 3; ++i) { + hQ2_in[i] = new TH1D(Form("hQ2_in_%d", i), "", NBINS, Q2bins.data()); + hX_in[i] = new TH1D(Form("hX_in_%d", i), "", NBINS, Xbins.data()); + hY_in[i] = new TH1D(Form("hY_in_%d", i), "", NBINS, Y_MIN, Y_MAX); + hW_in[i] = new TH1D(Form("hW_in_%d", i), "", NBINS, W_MIN, W_MAX); + } + + TH1D* hQ2_rec_all = new TH1D("hQ2_rec_all", "Q^{2} (rec)", NBINS, Q2bins.data()); + TH1D* hX_rec_all = new TH1D("hX_rec_all", "x (rec)", NBINS, Xbins.data()); + TH1D* hY_rec_all = new TH1D("hY_rec_all", "y (rec)", NBINS, Y_MIN, Y_MAX); + TH1D* hW_rec_all = new TH1D("hW_rec_all", "W (rec)", NBINS, W_MIN, W_MAX); + + TH1D* hQ2_rec_in[3], *hX_rec_in[3], *hY_rec_in[3], *hW_rec_in[3]; + for (int i = 0; i < 3; ++i) { + hQ2_rec_in[i] = new TH1D(Form("hQ2_rec_in_%d", i), "", NBINS, Q2bins.data()); + hX_rec_in[i] = new TH1D(Form("hX_rec_in_%d", i), "", NBINS, Xbins.data()); + hY_rec_in[i] = new TH1D(Form("hY_rec_in_%d", i), "", NBINS, Y_MIN, Y_MAX); + hW_rec_in[i] = new TH1D(Form("hW_rec_in_%d", i), "", NBINS, W_MIN, W_MAX); + } + + TH1D* hZ_proj = new TH1D("hZ_proj", "Muon track projections in nHCal; z [mm]; N_{proj}", NBINS, Z_MIN_MM, Z_MAX_MM); + + TH1D* hZ_hits = new TH1D("hZ_hits", "Reconstructed hit z in nHCal; z [mm]; N_{hits}", NBINS, Z_MIN_MM, Z_MAX_MM); + + TH3D* hDxDyZ_layer = new TH3D("hDxDyZ_layer","3D residuals (rec - proj): dx, dy vs z; dx [mm]; dy [mm]; z [mm]", NBINS, DXY_MIN_MM, DXY_MAX_MM, NBINS, DXY_MIN_MM, DXY_MAX_MM, NBINS, Z_MIN_MM, Z_MAX_MM); + + TH2D* hDxDy_all = new TH2D("hDxDy_all", + "Residuals (rec - proj): dx vs dy (all layers); dx [mm]; dy [mm]", + DXY_NBINS, DXY_MIN_MM, DXY_MAX_MM, DXY_NBINS, DXY_MIN_MM, DXY_MAX_MM); + + TH2D* hDrZ_layer = new TH2D("hDrZ_layer","Radial residual (rec - proj) vs z; dr [mm]; z [mm]", NBINS, DR_MIN_MM, DR_MAX_MM, NBINS, Z_MIN_MM, Z_MAX_MM); + + TH1D* hDr_all = new TH1D("hDr_all", + "Radial residual dr = #sqrt{dx^{2}+dy^{2}} (all layers); dr [mm]; N", + DR_NBINS, DR_MIN_MM, DR_MAX_MM); + + TH2D* hE_z = new TH2D("hE_z", "Reconstructed hit energy vs layer z; z_{layer} [mm]; E [GeV]", NBINS, Z_MIN_MM, Z_MAX_MM, NBINS, Ebins.data()); + + TH2D* hEsum_z = new TH2D("hEsum_z", "Layer energy sum vs z (reconstructed); z_{layer} [mm]; E_{sum} [GeV]", + NBINS, Z_MIN_MM, Z_MAX_MM, NBINS, Ebins.data()); + + TH1D* hE = new TH1D("hE", "Reconstructed hit energy; E [GeV]; N_{hits}", 10000, 0.0, 1.0); + + TH1D* hEsum = new TH1D("hEsum", "Reconstructed layer energy sum; E_{sum} [GeV]; N_{energy}", 10000, 0.0, 1.0); + + TH1D* hP_all_mu = new TH1D("hP_all_mu", "MC muon momentum (muons in nHCal acceptance); p_{MC} [GeV]; N", P_NBINS, P_MIN_GEV, P_MAX_GEV); + + TH1D* hP_pass_dr[3] = { + new TH1D(Form("hP_pass_dr%.1fcm",DR_CUTS_CM[0]), Form("Accepted (dr<%.1fcm); p_{MC} [GeV]; N",DR_CUTS_CM[0]), P_NBINS, P_MIN_GEV, P_MAX_GEV), + new TH1D(Form("hP_pass_dr%.1fcm",DR_CUTS_CM[1]), Form("Accepted (dr<%.1fcm); p_{MC} [GeV]; N",DR_CUTS_CM[1]), P_NBINS, P_MIN_GEV, P_MAX_GEV), + new TH1D(Form("hP_pass_dr%.1fcm",DR_CUTS_CM[2]), Form("Accepted (dr<%.1fcm); p_{MC} [GeV]; N",DR_CUTS_CM[2]), P_NBINS, P_MIN_GEV, P_MAX_GEV), + }; + + TH1D* hP_pass_ECut[3] = { + new TH1D(Form("hP_pass_E< %.1f MIP", E_CUTS[0]), Form("Accepted (E< %.1f MIP); p_{MC} [GeV]; N",E_CUTS[0]), P_NBINS, P_MIN_GEV, P_MAX_GEV), + new TH1D(Form("hP_pass_E< %.1f MIP", E_CUTS[1]), Form("Accepted (E< %.1f MIP); p_{MC} [GeV]; N",E_CUTS[1]), P_NBINS, P_MIN_GEV, P_MAX_GEV), + new TH1D(Form("hP_pass_E< %.1f MIP", E_CUTS[2]), Form("Accepted (E< %.1f MIP); p_{MC} [GeV]; N",E_CUTS[2]), P_NBINS, P_MIN_GEV, P_MAX_GEV), + }; + + TH1D* hP_pass_LayerCut[3] = { + new TH1D(Form("hP_pass_Layer< %d layers", LAYER_CUTS[0]), Form("Accepted (Layer < %d layers); p_{MC} [GeV]; N",LAYER_CUTS[0]), P_NBINS, P_MIN_GEV, P_MAX_GEV), + new TH1D(Form("hP_pass_Layer< %d layers", LAYER_CUTS[1]), Form("Accepted (Layer < %d layers); p_{MC} [GeV]; N",LAYER_CUTS[1]), P_NBINS, P_MIN_GEV, P_MAX_GEV), + new TH1D(Form("hP_pass_Layer< %d layers", LAYER_CUTS[2]), Form("Accepted (Layer < %d layers); p_{MC} [GeV]; N",LAYER_CUTS[2]), P_NBINS, P_MIN_GEV, P_MAX_GEV), + }; + + for (unsigned ev = 0; ev < nEvents; ev++) { + auto frameData = reader.readNextEntry(podio::Category::Event); + if (!frameData) continue; + podio::Frame frame(std::move(frameData)); + + auto& kinCol = frame.get("InclusiveKinematicsTruth"); + auto& mcCol = frame.get("MCParticles"); + auto& recParts = frame.get("ReconstructedParticles"); + auto& projSegs = frame.get("CalorimeterTrackProjections"); + auto& hcalRec = frame.get("HcalEndcapNRecHits"); + auto& assocCol = frame.get("ReconstructedParticleAssociations"); + + if (!kinCol.isValid() || !mcCol.isValid()) continue; + if (kinCol.empty()) continue; + + const auto& kin = kinCol.at(0); + double Q2 = kin.getQ2(), x = kin.getX(), y = kin.getY(), W = kin.getW(); + + bool m1_has_rec = false; + bool m2_has_rec = false; + + edm4hep::MCParticle m1, m2; + for (const auto& p : mcCol) { + if (abs(p.getPDG()) != 13 || p.getGeneratorStatus() == 0) continue; + if (p.getPDG() == 13) m1 = p; + else if (p.getPDG() == -13) m2 = p; + } + + TLorentzVector v1(m1.getMomentum().x, m1.getMomentum().y, m1.getMomentum().z, m1.getEnergy()); + TLorentzVector v2(m2.getMomentum().x, m2.getMomentum().y, m2.getMomentum().z, m2.getEnergy()); + + hEtaPt->Fill(v1.Eta(), v1.Pt()); + hEtaPt->Fill(v2.Eta(), v2.Pt()); + + hx_Q2->Fill(x,Q2); + hEta1_Eta2->Fill(v2.Eta(), v1.Eta()); + + hQ2_all->Fill(Q2); + hX_all->Fill(x); + hY_all->Fill(y); + hW_all->Fill(W); + + int nInNH = inNHCal(v1.Eta()) + inNHCal(v2.Eta()); + if (nInNH >= 0 && nInNH <= 2) { + hQ2_in[nInNH]->Fill(Q2); + hX_in[nInNH]->Fill(x); + hY_in[nInNH]->Fill(y); + hW_in[nInNH]->Fill(W); + } + + if(inNHCal(v1.Eta())) {hP_all_mu->Fill(v1.P());} + if(inNHCal(v2.Eta())){hP_all_mu->Fill(v2.P());} + + struct TaggedReco { edm4eic::ReconstructedParticle reco; int muTag; /*0=m1, 1=m2*/ }; + vector matchedRecos; + + auto find_associated_reco = [&](const edm4hep::MCParticle& mc)->void { + try { + if (!assocCol.isValid() || assocCol.empty()) return; + + const uint32_t mc_idx = mc.getObjectID().index; + + for (const auto& assoc : assocCol) { + if (assoc.getSimID() == mc_idx) { + uint32_t ridx = assoc.getRecID(); + if (!recParts.isValid() || ridx >= recParts.size()) continue; + auto reco = recParts.at(ridx); + if (reco.isAvailable()) { + matchedRecos.push_back(reco); + } + } + } + } catch (...) {} + }; + + if (!assocCol.isValid() || assocCol.empty()) continue; + if (m1.isAvailable() && abs(m1.getPDG())==13) find_associated_reco(m1, 0); + if (m2.isAvailable() && abs(m2.getPDG())==13) find_associated_reco(m2, 1); + + if (!recParts.isValid() || !projSegs.isValid() || !hcalRec.isValid() || recParts.empty() || projSegs.empty() || hcalRec.empty()) continue; + + set uniqueCentersZ10; + map layerData; + for (const auto& hit : hcalRec) { + double z = hit.getPosition().z; + double zBin = round(z); + uint64_t cid = hit.getCellID(); + double zCenter_mm = 10 * getPlasticCenterZ_cm(decoder, cellid_converter, cid, slice_index, /*plastic_slice_value=*/3); + // TVector3 dim = getPlasticDimensionsCM(*det, decoder, cid, slice_index, 3); + // printf("dim.x: %f, dim.y: %f, dim.z: %f", dim.x(), dim.y(), dim.z()); + // double dr = dim.x(); + // if(ev == 0){ + // DR_CUTS_CM[0] = dr; DR_CUTS_CM[1] = dr * 2; DR_CUTS_CM[2] = dr * 3; + // DR_THRESH_MM = DR_CUTS_CM[2]*10; + // } + if (!std::isnan(zCenter_mm)) { + int z10 = lround(zCenter_mm * 10.0); + uniqueCentersZ10.insert(z10); + } + + layerData[zBin] += hit.getEnergy(); + hZ_hits->Fill(z); + hE_z->Fill(z, hit.getEnergy()); + hE->Fill(hit.getEnergy()); + } + vector layerCentersZ; + layerCentersZ.reserve(uniqueCentersZ10.size()); + for (int z10 : uniqueCentersZ10) layerCentersZ.push_back(z10 / 10.0); + if(layerCentersZ.size() > LAYER_MAX) LAYER_MAX = layerCentersZ.size(); + + for(size_t n = 0; n < SIZE; n++) LAYER_CUTS[n] = static_cast(LAYER_MAX*LAYER_PROC[n]); + LAYER_THRESH = LAYER_CUTS[0]; + + for (const auto& [zValue, sumEnergy] : layerData) {hEsum_z->Fill(zValue, sumEnergy); hEsum->Fill(sumEnergy);} + + struct TaggedTrack { edm4eic::Track tr; int muTag; }; + vector allTracks; + for (const auto& R : matchedRecos) { + for (const auto& tr : R.reco.getTracks()) { + if (tr.isAvailable()) allTracks.push_back({tr, R.muTag}); + } + } + struct SegTag { edm4eic::TrackSegment seg; int muTag; }; + vector segsTagged; + for (const auto& seg : projSegs) { + auto linkedTr = seg.getTrack(); + if (!linkedTr.isAvailable()) continue; + for (const auto& TT : allTracks) { + if (linkedTr.getObjectID() == TT.tr.getObjectID()) { + segsTagged.push_back({seg, TT.muTag}); + break; + } + } + } + + for (const auto& ST : segsTagged) { + + vector segMinDistance(LAYER_MAX, std::numeric_limits::infinity()); + vector segHitEnergy(LAYER_MAX, std::numeric_limits::quiet_NaN()); + vector thickness_cm(LAYER_MAX, std::numeric_limits::quiet_NaN()); + vector count_DrCuts(SIZE, 0); + vector count_ECuts(SIZE, 0); + + GeomState A{}, B{}; + bool haveA = false, haveB = false; + + auto points = ST.seg.getPoints(); + + for (const auto& pt : points) { + if (pt.system != 113) continue; + hZ_proj->Fill(pt.position.z); + if (!haveA) { + A.x = pt.position.x; A.y = pt.position.y; A.z = pt.position.z; + haveA = true; + } else if (!haveB) { + B.x = pt.position.x; B.y = pt.position.y; B.z = pt.position.z; + haveB = true; + break; + } + } + + if (!haveA || !haveB) {continue;} + + for (int i = 0; i < LAYER_MAX; ++i) { + double best_dr_in_layer = 1e10; + double best_E_in_layer; + double partLayerEnergySum = 0; + double ratio_HitPartLayerEnergy = 0; + dd4hep::DDSegmentation::CellID best_cid_in_layer; + + double zc = layerCentersZ[i]; + auto [X, Y] = trackXYatZ(A, B, zc); + + for (const auto& hit : hcalRec) { + const auto& hp = hit.getPosition(); + + if (fabs(hp.z - zc) > 5.0 /*mm*/) continue; + + const double dx = X - hp.x; + const double dy = Y - hp.y; + const double dr = sqrt(dx*dx + dy*dy); + if (dr < 30*10 /*mm*/) partLayerEnergySum += hit.getEnergy(); + + hDxDyZ_layer->Fill(dx, dy, hp.z); + hDxDy_all->Fill(dx, dy); + hDrZ_layer->Fill(dr, hp.z); + hDr_all->Fill(dr); + + if (dr < best_dr_in_layer) { + best_dr_in_layer = dr; + best_E_in_layer = hit.getEnergy(); + best_cid_in_layer = hit.getCellID(); + } + } + + if (best_dr_in_layer < DR_THRESH_MM) { + segMinDistance[i] = best_dr_in_layer; + segHitEnergy[i] = best_E_in_layer; + thickness_cm[i] = getPlasticDimensionsCM(*det, decoder, best_cid_in_layer, slice_index, 3).z(); + t_cm = thickness_cm[0]; + } + ratio_HitPartLayerEnergy += segHitEnergy[i]/partLayerEnergySum; + if (std::isnan(ratio_HitPartLayerEnergy) || fabs(ratio_HitPartLayerEnergy - 1) >= 0.2) continue; + for(int j = 0; j < SIZE; ++j){ + if (segMinDistance[i] < DR_CUTS_CM[j] * 10 /*mm*/){++count_DrCuts[j];} + if (segHitEnergy[i] < (E_CUTS[j] * MIP_ENERGY_GEV * thickness_cm[i] * 100)){++count_ECuts[j];} + } + } + + if (ST.muTag == 0 && count_DrCuts[SIZE-1] > LAYER_THRESH) m1_has_rec = true; + if (ST.muTag == 1 && count_DrCuts[SIZE-1] > LAYER_THRESH) m2_has_rec = true; + + for (int j = 0; j < SIZE; ++j){ + const double p = ST.muTag ? v2.P() : v1.P(); + if (count_DrCuts[j] > LAYER_THRESH) hP_pass_dr[j]->Fill(p); + if (count_ECuts[j] > LAYER_THRESH) hP_pass_ECut[j]->Fill(p); + if (count_DrCuts[SIZE-1] > LAYER_CUTS[j]) hP_pass_LayerCut[j]->Fill(p); + } + } + + int nInNH_rec_local = 0; + if (m1_has_rec) ++nInNH_rec_local; + if (m2_has_rec) ++nInNH_rec_local; + + hQ2_rec_all->Fill(Q2); + hX_rec_all->Fill(x); + hY_rec_all->Fill(y); + hW_rec_all->Fill(W); + + if (nInNH_rec_local >= 0 && nInNH_rec_local <= 2) { + hQ2_rec_in[nInNH_rec_local]->Fill(Q2); + hX_rec_in[nInNH_rec_local]->Fill(x); + hY_rec_in[nInNH_rec_local]->Fill(y); + hW_rec_in[nInNH_rec_local]->Fill(W); + } + } + + TCanvas* c = new TCanvas("c", "Muon analysis", 1600, 1000); + c->Divide(2,2); + + c->cd(1); hEtaPt->Draw("COLZ"); + c->cd(2); gPad->SetLogx(); gPad->SetGrid(); gPad->SetLogy(); hx_Q2->Draw("COLZ"); + c->cd(3); hEta1_Eta2->Draw("COLZ"); + auto* lv1 = makeLine(ETA_MIN, -6, ETA_MIN, 6); + auto* lv2 = makeLine(ETA_MAX, -6, ETA_MAX, 6); + auto* lh1 = makeLine(hEta1_Eta2->GetXaxis()->GetXmin(), ETA_MIN, + hEta1_Eta2->GetXaxis()->GetXmax(), ETA_MIN); + auto* lh2 = makeLine(hEta1_Eta2->GetXaxis()->GetXmin(), ETA_MAX, + hEta1_Eta2->GetXaxis()->GetXmax(), ETA_MAX); + + c->SaveAs(outname_pdf.c_str()); + c->SaveAs(outname_png.c_str()); + + + TCanvas* c_sim = new TCanvas("c_sim", "Muon geometrical accuracy analysis", 1600, 1000); + c_sim->Divide(2, 2); + + c_sim->cd(1); drawEffPanel(hQ2_all, hQ2_in, "Geom. acc. vs Q^{2}", "Q^{2} [GeV^{2}]", true); + c_sim->cd(2); drawEffPanel(hX_all, hX_in, "Geom. acc. vs x", "x", true); + c_sim->cd(3); drawEffPanel(hY_all, hY_in, "Geom. acc. vs y", "y", false); + c_sim->cd(4); drawEffPanel(hW_all, hW_in, "Geom. acc. vs W", "W [GeV]", false); + + c_sim->SaveAs(addPrefixAfterSlash(outname_png, "sim_").c_str()); + c_sim->SaveAs(addPrefixAfterSlash(outname_pdf, "sim_").c_str()); + + TCanvas* c_rec = new TCanvas("c_rec", "Muon geometrical accuracy analysis rec", 1600, 1000); + c_rec->Divide(2, 2); + + c_rec->cd(1); drawEffPanel(hQ2_rec_all, hQ2_rec_in, "Geom. acc. vs Q^{2}", "Q^{2} [GeV^{2}]", true); + c_rec->cd(2); drawEffPanel(hX_rec_all, hX_rec_in, "Geom. acc. vs x", "x", true); + c_rec->cd(3); drawEffPanel(hY_rec_all, hY_rec_in, "Geom. acc. vs y", "y", false); + c_rec->cd(4); drawEffPanel(hW_rec_all, hW_rec_in, "Geom. acc. vs W", "W [GeV]", false); + + c_rec->SaveAs(addPrefixAfterSlash(outname_png, "rec_").c_str()); + c_rec->SaveAs(addPrefixAfterSlash(outname_pdf, "rec_").c_str()); + + TCanvas* c_hitTrack = new TCanvas("c_hitTrack", "Muon hit and track analysis", 1600, 1000); + c_hitTrack->Divide(2,2); + c_hitTrack->cd(1); gPad->SetGrid(); hZ_proj->Draw(); + c_hitTrack->cd(2); gPad->SetGrid(); hZ_hits->Draw(); + c_hitTrack->SaveAs(addPrefixAfterSlash(outname_png, "hitTrack_").c_str()); + c_hitTrack->SaveAs(addPrefixAfterSlash(outname_pdf, "hitTrack_").c_str()); + + TCanvas* c_hitTrackDistance = new TCanvas("c_hitTrackDistance", "Muon hit-track distance analysis", 1600, 1000); + c_hitTrackDistance->Divide(2,2); + c_hitTrackDistance->cd(1); gPad->SetGrid(); hDxDyZ_layer->Draw("COLZ"); + c_hitTrackDistance->cd(2); gPad->SetGrid(); hDxDy_all->Draw("COLZ"); + c_hitTrackDistance->cd(3); gPad->SetGrid(); hDrZ_layer->Draw("COLZ"); + c_hitTrackDistance->cd(4); gPad->SetGrid(); hDr_all->Draw("COLZ"); + c_hitTrackDistance->SaveAs(addPrefixAfterSlash(outname_png, "hitTrackDistance_").c_str()); + c_hitTrackDistance->SaveAs(addPrefixAfterSlash(outname_pdf, "hitTrackDistance_").c_str()); + + TCanvas* c_Eff = new TCanvas("c_Eff", "Muon efficiency analysis", 1600, 1000); + c_Eff->Divide(2,2); + c_Eff->cd(1); gPad->SetGrid(); hP_all_mu->SetLineColor(kBlack); hP_all_mu->Draw(); + c_Eff->cd(2); gPad->SetGrid(); + TH1D* hEff_dr[3]; + for (int idr=0; idr<3; ++idr) { + hEff_dr[idr] = (TH1D*)hP_pass_dr[idr]->Clone( + Form("hEff_dr%.1fcm", DR_CUTS_CM[idr]) + ); + hEff_dr[idr]->SetDirectory(nullptr); + hEff_dr[idr]->SetTitle(Form("Efficiency vs p_{MC} (Layers > %d); p_{MC} [GeV]; efficiency", LAYER_THRESH)); + hEff_dr[idr]->Divide(hP_pass_dr[idr], hP_all_mu, 1, 1, "B"); + } + + hEff_dr[0]->SetLineColor(kBlue); hEff_dr[0]->SetMinimum(0.0); hEff_dr[0]->SetMaximum(1.4); + hEff_dr[1]->SetLineColor(kRed); + hEff_dr[2]->SetLineColor(kGreen+2); + + hEff_dr[0]->Draw("HIST E"); + hEff_dr[1]->Draw("HIST E SAME"); + hEff_dr[2]->Draw("HIST E SAME"); + { auto leg_dr = new TLegend(0.60,0.70,0.88,0.88); + for(int idr=0; idr<3; ++idr) + {leg_dr->AddEntry(hEff_dr[idr],Form("dr < %.1f cm", DR_CUTS_CM[idr]),"l");} + leg_dr->Draw(); + } + + c_Eff->cd(3); gPad->SetGrid(); + TH1D* hEff_E[3]; + for (int ie=0; ie<3; ++ie) { + hEff_E[ie] = (TH1D*)hP_pass_ECut[ie]->Clone( + Form("hEff_E%.1fcm", E_CUTS[ie]) + ); + hEff_E[ie]->SetDirectory(nullptr); + hEff_E[ie]->SetTitle(Form("Efficiency vs p_{MC} (Layers > %d, dr < %.1f cm); p_{MC} [GeV]; efficiency",LAYER_THRESH, DR_THRESH_MM/10)); + hEff_E[ie]->Divide(hP_pass_ECut[ie], hP_all_mu, 1, 1, "B"); + } + + hEff_E[0]->SetLineColor(kBlue); hEff_E[0]->SetMinimum(0.0); hEff_E[0]->SetMaximum(1.4); + hEff_E[1]->SetLineColor(kRed); + hEff_E[2]->SetLineColor(kGreen+2); + + hEff_E[0]->Draw("HIST E "); + hEff_E[1]->Draw("HIST E SAME"); + hEff_E[2]->Draw("HIST E SAME"); + { auto leg_E = new TLegend(0.60,0.70,0.88,0.88); + for(int ie=0; ie<3; ++ie) + {leg_E->AddEntry(hEff_E[ie],Form("E < %.1f MIP", E_CUTS[ie]),"l");} + leg_E->Draw(); + } + + c_Eff->cd(4); gPad->SetGrid(); + TH1D* hEff_Layer[3]; + for (int il=0; il<3; ++il) { + hEff_Layer[il] = (TH1D*)hP_pass_LayerCut[il]->Clone( + Form("hEff_Layer%d", LAYER_CUTS[il]) + ); + hEff_Layer[il]->SetDirectory(nullptr); + hEff_Layer[il]->SetTitle(Form("Efficiency vs p_{MC} (dr < %.1f cm); p_{MC} [GeV]; efficiency", DR_THRESH_MM/10)); + hEff_Layer[il]->Divide(hP_pass_LayerCut[il], hP_all_mu, 1, 1, "B"); + } + + hEff_Layer[0]->SetLineColor(kBlue); hEff_Layer[0]->SetMinimum(0.0); hEff_Layer[0]->SetMaximum(1.4); + hEff_Layer[1]->SetLineColor(kRed); + hEff_Layer[2]->SetLineColor(kGreen+2); + + hEff_Layer[0]->Draw("HIST E "); + hEff_Layer[1]->Draw("HIST E SAME"); + hEff_Layer[2]->Draw("HIST E SAME"); + { auto leg_Layer = new TLegend(0.60,0.70,0.88,0.88); + for(int il=0; il<3; ++il) + {leg_Layer->AddEntry(hEff_Layer[il],Form("L > %d Layers", LAYER_CUTS[il]),"l");} + leg_Layer->Draw(); + } + c_Eff->SaveAs(addPrefixAfterSlash(outname_png, "matching_efficiency_").c_str()); + c_Eff->SaveAs(addPrefixAfterSlash(outname_pdf, "matching_efficiency_").c_str()); + + TCanvas* c_hEnergy = new TCanvas("c_hEnergy", "Muon hit energy analysis", 1600, 1000); + c_hEnergy->Divide(2,2); + + c_hEnergy->cd(1); gPad->SetGrid(); gPad->SetLogy(); hE_z->GetYaxis()->SetMoreLogLabels(); hE_z->Draw("COLZ"); + TProfile* pE_z = hE_z->ProfileX("pE_z"); pE_z->SetLineWidth(3); pE_z->SetLineColor(kRed); pE_z->SetMarkerColor(kRed); + pE_z->SetDirectory(nullptr); pE_z->Draw("SAME"); + + const double Ecut_GeV = MIP_ENERGY_GEV * t_cm * 100; + TLine* cut = new TLine(hE_z->GetXaxis()->GetXmin(), Ecut_GeV, hE_z->GetXaxis()->GetXmax(), Ecut_GeV); + cut->SetLineColor(kRed+1); cut->SetLineStyle(2); cut->SetLineWidth(2); + cut->Draw("SAME"); + auto* leg = new TLegend(0.58, 0.75, 0.84, 0.88); leg->SetTextSize(0.04);leg->AddEntry(cut, Form("E_{cut} = %.3f GeV", Ecut_GeV), "l"); + leg->Draw(); + + c_hEnergy->cd(2); gPad->SetGrid(); gPad->SetLogy(); hEsum_z->GetYaxis()->SetMoreLogLabels(); hEsum_z->Draw("COLZ"); + TProfile* pEsum_z = hEsum_z->ProfileX("pEsum_z"); pEsum_z->SetLineWidth(3); pEsum_z->SetLineColor(kRed); pEsum_z->SetMarkerColor(kRed); + pEsum_z->SetDirectory(nullptr); pEsum_z->Draw("SAME"); + c_hEnergy->cd(3); gPad->SetGrid(); hE->Draw(); + c_hEnergy->cd(4); gPad->SetGrid(); hEsum->Draw(); + + c_hEnergy->SaveAs(addPrefixAfterSlash(outname_png, "hit_energy_").c_str()); + c_hEnergy->SaveAs(addPrefixAfterSlash(outname_pdf, "hit_energy_").c_str()); + + delete c; + delete c_sim; + delete c_rec; + delete c_hitTrack; + delete c_hitTrackDistance; + delete c_Eff; + delete c_hEnergy; + + return 0; +} diff --git a/benchmarks/nhcal_dimuon_photoproduction/scripts/ep_DiMuon.cmnd b/benchmarks/nhcal_dimuon_photoproduction/scripts/ep_DiMuon.cmnd new file mode 100644 index 00000000..9fb98946 --- /dev/null +++ b/benchmarks/nhcal_dimuon_photoproduction/scripts/ep_DiMuon.cmnd @@ -0,0 +1,180 @@ +#============================================================================== +# STAR Heavy Flavor Tune 1.1 +# +# PYTHIA Version 8.1.08 +# Date: November 7, 2008 +# Last updated by: Thomas Ullrich +# +# This file contains commands to be read in for a Pythia8 run. +# Lines not beginning with a letter or digit are comments. +# Names are case-insensitive - but spellings-sensitive! +# +# Important notes: +# ================ +# For hard 2->2 processes, such as most heavy flavor producing one, +# there's divergence in PYTHIA that needs to be suppressed. Use +# these two lines in the code before you initialize PYTHIA to avoid +# any trouble: +# UserHooks *oniumUserHook = new SuppressSmallPT(); +# pythia.setUserHooksPtr(oniumUserHook); +# +#============================================================================== + +#------------------------------------------------------------------------------ +# Parameters that need to be set by whoever runs it. +# Note that they have no meaning unless restored and used +# in the user provided code (main program). +# This is not part of the star_hf tune (just convenient) +# Documentation: /htmldoc/MainProgramSettings.html +#------------------------------------------------------------------------------ +Main:numberOfEvents = 1000 ! number of events to generate +Next:numberShowEvent = 20 ! show event record +Next:numberCount = 0 ! number of events to print +Next:numberShowInfo = 100 ! show how far along run is this many times +Main:timesAllowErrors = 100 ! abort run after this many flawed events +Init:showChangedSettings = on ! print changed flags/modes/parameters +Init:showAllSettings = on ! print all flags/modes/parameters +#Init:showAllStatistics = on ! print statistics at the end + +#------------------------------------------------------------------------------ +# Colliding beams and center-of-mass energy +# Documentation: /htmldoc/BeamParameters.html +#------------------------------------------------------------------------------ +Beams:frameType = 2 +Beams:idA = 2212 ! proton +Beams:idB = 11 ! e- +Beams:eA = 275. ! RHIC nominal (GeV) +Beams:eB = 18. ! RHIC nominal (GeV) + +#------------------------------------------------------------------------------ +# Process Selection +#------------------------------------------------------------------------------ + + +# Enable equivalent photon approximation (EPA) for both beams +Photon:Q2Max = 1.0 ! Upper Q^2 limit for EPA photons (in GeV^2) +Photon:ProcessType = 4 ! 4 = direct-direct photons +#Photon:EPA = on + +PDF:beamA2gamma = on ! EPA photon flux from beam A +PDF:beamB2gamma = on ! EPA photon flux from beam B + +# Enable gamma-gamma -> mu+ mu− +PhotonCollision:gmgm2mumu = on + +# Optional: Turn off other QED or QCD backgrounds if you want exclusivity +PartonLevel:ISR = off +PartonLevel:FSR = off +HadronLevel:all = off + + +#------------------------------------------------------------------------------ +# K factor +# Multiply almost all cross sections by this common fix factor. +# This is usually no very useful. The data can be shifted up and down +# later anyhow as we please. +# Documentation: /htmldoc/CouplingsAndScales.html +#------------------------------------------------------------------------------ +# SigmaProcess:Kfactor = 3 + +#------------------------------------------------------------------------------ +# Scales (Ramona's suggestions) +# This sets the scale to settings typically for hard probes: +# mu_F = mu_R = 2*mT +# Documentation: /htmldoc/CouplingsAndScales.html +#------------------------------------------------------------------------------ +# SigmaProcess:renormScale2 = 3 +# SigmaProcess:factorScale2 = 3 +# SigmaProcess:renormMultFac = 2 ! 2mT +# SigmaProcess:factorMultFac = 2 ! 2mT + +#------------------------------------------------------------------------------ +# To limit particle production to a certain pthat range uncomment +# these lines. Use only when you 100% know what you are doing. +# It is extremely useful to split runs up in ptHat bins to generate +# statistics evenly in pt. Book keeping is important then (cross-sections, +# number of events) to compile the final complete spectra. +# Documentation: /htmldoc/PhaseSpaceCuts.html +#------------------------------------------------------------------------------ +# PhaseSpace:pTHatMin = 4 ! Minimal pT cut +# PhaseSpace:pTHatMax = 2 + +#------------------------------------------------------------------------------ +# Random Number +# Initialize random generator according to time. Otherwise multiple jobs +# will produce the same sequence (unless you pass a different seed every +# time which is not practical). +# Documentation: /htmldoc/RandomNumberSeed.html +#------------------------------------------------------------------------------ +Random:setSeed = on +Random:seed = 0 + +#------------------------------------------------------------------------------ +# PDF Selection: +# Note: you need LHAPDF to be installed. Pythia 8 only provides a +# minimal set to get started. +# The choice of PDF here is greatly motivated by: +# A.~Sherstnev and R.~S.~Thorne, arXiv:0807.2132 and arXiv:0711.2473v3 +# and W. Vogelsang (private communication) +# These are PDFs especially made for LO Monte-Carlo generators such +# as PYTHIA. +# The state-of-the-art NLO PDF is cteq66.LHgrid which can be used +# as an alternative (also from LHAPDF) but with the known problems +# that arise when using a NLO PDF in an LO simulator. +# Documentation: /htmldoc/PDFSelection.html +#------------------------------------------------------------------------------ +#PDF:useLHAPDF = off +#PDF:LHAPDFset = MRSTMCal.LHgrid +#PDF:extrapolateLHAPDF = on +#PDF:pSet = 8 +#PDF:beamA2gamma = on + +#------------------------------------------------------------------------------ +# Settings for the event generation process in the Pythia8 library +# Effect/Relevance of MI, ISR and FSR need to be checked. For sure +# the use more CPU and hence less events/s. +# If HadronLevel:Hadronize = off we end up with the pure c, b spectra +# (which might be useful at times) +# Documentation: /htmldoc/MasterSwitches.html +# Documentation: /htmldoc/MultipleInteractions.html +#------------------------------------------------------------------------------ +#PartonLevel:MPI = on ! multiple interactions +#PartonLevel:ISR = on ! initial-state radiation +#BeamRemnants:primordialKT = on ! primordial kt +#PartonLevel:FSR = on ! final-state radiation +#HadronLevel:Hadronize = off ! no hadronization use + +#------------------------------------------------------------------------------ +# Relative production ratio vector/pseudoscalar for charm and bottom mesons +# This was originally PARJ(13) where PARJ(13) = V/(PS+V) that is the +# vector meson fraction of primary charm+bottom mesons. +# Andre David (CERN/NA60) made an exhaustive study and found that the +# world data supports 0.6 while PYTHIA default was PARJ(13) = 3/4 = 0.75 +# from simple spin counting. +# In PYTHIA8 we now use V/PS not V/(PS+V) +# Documentation: /htmldoc/FlavourSelection.html +#------------------------------------------------------------------------------ +#StringFlav:mesonCvector = 1.5 ! same as PARJ(13)=0.6 +#StringFlav:mesonBvector = 3 ! leave at 0.75 + +#------------------------------------------------------------------------------ +# Heavy quark masses. +# Note that this should match with the ones used in the PDF. +# The masses are listed in the header of the refering PDF file. +# Documentation: /htmldoc/ParticleDataScheme.html +# Documentation: /htmldoc/ParticleData.html +#------------------------------------------------------------------------------ +#4:m0 = 1.43 +#5:m0 = 4.30 + +#------------------------------------------------------------------------------ +# Particle Decay limits +# When on, only particles with a decay within a volume limited by +# rho = sqrt(x^2 + y^2) < xyMax and |z| < zMax are decayed. +# The above xyMax, expressed in mm/c. +#------------------------------------------------------------------------------ +# ParticleDecays:limitCylinder = on +# ParticleDecays:xyMax = 600 +# ParticleDecays:zMax = 1000 + +# EOF diff --git a/benchmarks/nhcal_dimuon_photoproduction/scripts/libPythiaEvent.so b/benchmarks/nhcal_dimuon_photoproduction/scripts/libPythiaEvent.so new file mode 100755 index 00000000..5a83f2ae Binary files /dev/null and b/benchmarks/nhcal_dimuon_photoproduction/scripts/libPythiaEvent.so differ diff --git a/benchmarks/nhcal_dimuon_photoproduction/scripts/pythiaDiMuon.cpp b/benchmarks/nhcal_dimuon_photoproduction/scripts/pythiaDiMuon.cpp new file mode 100755 index 00000000..890dcf86 --- /dev/null +++ b/benchmarks/nhcal_dimuon_photoproduction/scripts/pythiaDiMuon.cpp @@ -0,0 +1,1490 @@ +//============================================================================== + +//============================================================================== + +#include "Pythia8/Pythia.h" +#include "TTree.h" +#include "TFile.h" +#include "TMath.h" + +#include +#include +#include + +#include "fastjet/ClusterSequence.hh" +#include + +#include "Pythia8Plugins/HepMC3.h" + +#include "PythiaEvent.h" +#include "HistogramsPythia.h" + +#define PR(x) std::cout << #x << " = " << (x) << std::endl; + +using namespace fastjet; +using namespace Pythia8; +using namespace std; + +double costhetastar(int, int, const Event&); +bool isInAcceptance(int, const Event&); // acceptance filter +int MakeEvent(Pythia *pythia, PythiaEvent *eventStore, int iev, bool writeTree = false); // event handler (analyze event and + // stores data in tuple) +int findFinalElectron(const Event& event); + +void FindPartonJet(vector jets, Particle parton1, Particle parton2, int &jetid1, int &jetid2); +double dR(fastjet::PseudoJet jet, Particle parton); + +int isInHcals(Particle part); +int isInHcals(double eta); +bool isInTracker(double eta); +//int isInTracker(Particle part); + +int FillHCals(TH2F *hist, TH1F *hEne, TH1F *hEneDenom, bool &anyHcal_jets); +int FillHCalsJets(TH2F *hist, vector jets); +int FillHCalsJetsShare(TH2F *hist, vector jets); + +void SortPairs(std::vector> &input, std::vector> &output); + +bool compFunc(pair a, pair b); +void PrintVec(std::vector> input); + +int main(int argc, char* argv[]) { + + if (argc != 3) { + cout << "Usage: " << argv[0] << " runcard outfile" << endl; + return 2; + } + char* runcard = argv[1]; + char* outfile = argv[2]; + //const char* xmlDB = "/users/PAS2524/lkosarz/Pythia/pythia8312/share/Pythia8/xmldoc"; + const char* xmlDB = "/opt/local/share/Pythia8/xmldoc"; + + bool WriteHepMC = true; + bool WriteTree = false; + + // https://pythia.org//latest-manual/examples/main345.html + + // + // Create instance of Pythia + // + //Pythia pythia(xmlDB); // the default parameters are read from xml files + // stored in the xmldoc directory. This includes + // particle data and decay definitions. + + Pythia *pythia = new Pythia("/opt/local/share/Pythia8/xmldoc"); + + + // Shorthand for (static) settings + Settings& settings = pythia->settings; + + // Read in runcard (should be star_hf_tune_v1.0.cmd) + pythia->readFile(runcard); + cout << "Runcard '" << runcard << "' loaded." << endl; + + + TFile *treeFile = new TFile(Form("%s_tree.root", outfile),"RECREATE"); + TFile *histFile = new TFile(Form("%s_hist.root", outfile),"RECREATE"); + + treeFile->cd(); + + TTree *eventTree = new TTree("eventTree", "Event tree"); + PythiaEvent *eventStore = new PythiaEvent(); + eventTree->Branch("eventTree", &eventStore, 32000, 1); + + + + + // + // Retrieve number of events and other parameters from the runcard. + // We need to deal with those settings ourself. Getting + // them through the runcard just avoids recompiling. + // + long maxNumberOfEvents = settings.mode("Main:numberOfEvents"); + //int nList = settings.mode("Main:numberToList"); + int nList = 10; + //int nShow = settings.mode("Main:timesToShow"); + int nShow = 10; + int maxErrors = settings.mode("Main:timesAllowErrors"); + //bool showCS = settings.flag("Main:showChangedSettings"); + //bool showAS = settings.flag("Main:showAllSettings"); + int pace = maxNumberOfEvents/nShow; + + + // Initialize Pythia, ready to go + pythia->init(); + HepMC3::WriterAscii hepmcWriter(Form("%s.hepmc3", outfile)); + HepMC3::Pythia8ToHepMC3 toHepMC; + + // List changed or all data + //if (showCS) settings.listChanged(); + //if (showAS) settings.listAll(); + settings.listChanged(); + settings.listAll(); + + + + //IO_GenEvent ascii_io(Form("%s.hepmc3", outfile)); + //HepMC::WriterRootTree WriterRootfile("file.root") + //GenEvent hepmcevt; + + + histFile->cd(); + CreateHistograms(); + + //-------------------------------------------------------------- + // Event loop + //-------------------------------------------------------------- + int ievent = 0; + int iErrors = 0; + + while (ievent < maxNumberOfEvents) { + + if (!pythia->next()) { + if (++iErrors < maxErrors) continue; + cout << "Error: too many errors in event generation - check your settings & code" << endl; + break; + } + + h_Events_cuts->Fill(0); + + h_XsecGen->Fill(pythia->info.sigmaGen()); + h_XsecGen_err->Fill(pythia->info.sigmaGen(), pythia->info.sigmaErr()); + + //-------------------- + + bool isDiffractive = pythia->info.isDiffractiveA() || pythia->info.isDiffractiveB(); + bool isHardDiffractive = pythia->info.isHardDiffractiveA() || pythia->info.isHardDiffractiveB(); + + //if (isDiffractive) cout<<"Diffractive"<Fill(1); + if (isHardDiffractive) h_Events_cuts->Fill(2); + + //if(!(isDiffractive || isHardDiffractive)) continue; + + int phMode = pythia->info.photonMode(); + + //------------------------- + + h_XsecSel->Fill(pythia->info.sigmaGen()); + h_XsecSel_err->Fill(pythia->info.sigmaGen(), pythia->info.sigmaErr()); + + MakeEvent(pythia, eventStore, ievent, WriteTree); // in MakeEvent we deal with the whole event and return + + if(WriteTree) + { + eventTree->Fill(); + eventStore->Clear(); + } + ievent++; + + if (ievent%pace == 0) { + cout << "# of events generated = " << ievent << endl; + } + + // List first few events. + if (ievent < nList) { + pythia->info.list(); + pythia->process.list(); + pythia->event.list(false, true, 3); + } + + + if (WriteHepMC) { + HepMC3::GenEvent evt; + toHepMC.fill_next_event(*pythia, &evt); + hepmcWriter.write_event(evt); + } + + //GenEvent* hepmcevt = new HepMC::GenEvent(); + //toHepMC.fill_next_event( *pythia.Pythia8(), hepmcevt); + //ascii_io << hepmcevt; + //delete hepmcevt; + + } + + //-------------------------------------------------------------- + // Finish up + //-------------------------------------------------------------- + //pythia.statistics(); + pythia->stat(); + + cout << "Writing files" << endl; + + histFile->cd(); + histFile->Write(); + + treeFile->cd(); + if(WriteTree) eventTree->Write(); + //treeFile->Write(); + + histFile->Close(); + treeFile->Close(); + hepmcWriter.close(); + + cout << "Finish!" << endl; + + return 0; +} + +// +// Event analysis +// +int MakeEvent(Pythia *pythia, PythiaEvent *eventStore, int iev, bool writeTree) +{ + Event &event = pythia->event; + + //cout<<"NEW EVENT ---------------------------"<info.isDiffractiveA()) h_Events_Diffractive->Fill(0); + if(pythia->info.isDiffractiveB()) h_Events_Diffractive->Fill(1); + if(pythia->info.isHardDiffractiveA()) h_Events_Diffractive->Fill(2); + if(pythia->info.isHardDiffractiveB()) h_Events_Diffractive->Fill(3); + + + hist_eta_energy_tmp->Reset(); + hist_eta_energy_denom_tmp->Reset(); + + //if (event[1].id() == 11) + //int eleid = findFinalElectron(event); + + // create a jet definition: + // a jet algorithm with a given radius parameter + //---------------------------------------------------------- + double R = 1.0; + double p = -1.0; + fastjet::JetDefinition jet_def(fastjet::antikt_algorithm, R); + //fastjet::JetDefinition jet_def(fastjet::ee_kt_algorithm); + //fastjet::JetDefinition jet_def(fastjet::ee_genkt_algorithm, R, p); + //fastjet::JetDefinition jet_def_meas(fastjet::ee_kt_algorithm); + + vector input_particles; + vector input_particles_meas; + vector input_particles_meas_no_nHCal; + + float EsumF = 0.0; + double EsumD = 0.0; + + bool has_nHcalActivity = false; + + int nPartFinal = 0; + int nPartonsOut = 0; + int iParton_1 = 0; + int iParton_2 = 0; + + bool Parton_1_inAcc = false; + bool Parton_2_inAcc = false; + + int nJetsInAcc = 0; + + int nPion_p = 0; + int nPion_n = 0; + int nKaon_p = 0; + int nKaon_n = 0; + int nProton_p = 0; + int nProton_n = 0; + int nElectron_p = 0; + int nElectron_n = 0; + + int nNeutron = 0; + int nGamma = 0; + + + // Four-momenta of proton, electron, virtual photon/Z^0/W^+-. + Vec4 pProton = event[1].p(); + Vec4 peIn = event[2].p(); + Vec4 pPhoton = event[4].p(); + + //Vec4 peIn = event[4].p(); + //Vec4 peOut = event[6].p(); + //Vec4 pPhoton = peIn - peOut; + + // Q2, W2, Bjorken x, y. + double Q2 = - pPhoton.m2Calc(); + double W2 = (pProton + pPhoton).m2Calc(); + double x = Q2 / (2. * pProton * pPhoton); + double y = (pProton * pPhoton) / (pProton * peIn); + + + for (int i = 0; i < event.size(); i++) { + + Particle part = event[i]; + + // count outgoing partons + if(part.status()==-23 || part.status()==-24) + { + nPartonsOut++; + + h_Partons_status->Fill(part.status()); + + h_Parton_eta_p->Fill(part.eta(), part.pAbs()); + h_Parton_eta_pT->Fill(part.eta(), part.pT()); + h_Parton_eta_E->Fill(part.eta(), part.e()); + + if(nPartonsOut == 1) iParton_1 = i; + if(nPartonsOut == 2) iParton_2 = i; + + if(part.id()>0) h_Partons_types->Fill(part.id()); + if(part.id()<0) h_Partons_types_anti->Fill(abs(part.id())); + + if((part.eta()>=-4.14 && part.eta()<=4.2) && nPartonsOut == 1) Parton_1_inAcc = true; + if((part.eta()>=-4.14 && part.eta()<=4.2) && nPartonsOut == 2) Parton_2_inAcc = true; + } + + if(nPartonsOut == 2) + { + h_Partons_eta->Fill(event[iParton_1].eta(), event[iParton_2].eta()); + h_Partons_phi->Fill(event[iParton_1].phi(), event[iParton_2].phi()); + h_Partons_p->Fill(event[iParton_1].pAbs(), event[iParton_2].pAbs()); + h_Partons_pT->Fill(event[iParton_1].pT(), event[iParton_2].pT()); + + h_Parton_x_eta->Fill(x, event[iParton_1].eta()); + h_Parton_x_eta->Fill(x, event[iParton_2].eta()); + + h_Parton_y_eta->Fill(y, event[iParton_1].eta()); + h_Parton_y_eta->Fill(y, event[iParton_2].eta()); + + h_Parton_x_eta1->Fill(x, event[iParton_1].eta()); + h_Parton_x_eta2->Fill(x, event[iParton_2].eta()); + + h_Parton_y_eta1->Fill(y, event[iParton_1].eta()); + h_Parton_y_eta2->Fill(y, event[iParton_2].eta()); + } + + // accept final particles only + if(!part.isFinal()) continue; + // ePIC acceptance only + if(part.eta()<-4.14 || part.eta()>4.2) continue; + if(!(11>part.status() || part.status()>19)) continue; // exclude beam particles + + nPartFinal++; + + + if(part.id() == 211) nPion_p++; + if(part.id() == -211) nPion_n++; + if(part.id() == 321) nKaon_p++; + if(part.id() == -321) nKaon_n++; + if(part.id() == 2212) nProton_p++; + if(part.id() == -2212) nProton_n++; + if(part.id() == -11) nElectron_p++; + if(part.id() == 11) nElectron_n++; + + if(part.id() == 2112) nNeutron++; + if(part.id() == 22) nGamma++; + + h_Particle_eta->Fill(part.eta()); + h_Particle_eta_wE->Fill(part.eta(), part.e()); + + h_Particle_eta_p->Fill(part.eta(), part.pAbs()); + h_Particle_eta_pT->Fill(part.eta(), part.pT()); + h_Particle_eta_E->Fill(part.eta(), part.e()); + + + // eta, momentum + if(part.id() == 211) h_Particle_pion_p_eta_p->Fill(part.eta(), part.pAbs()); + if(part.id() == -211) h_Particle_pion_n_eta_p->Fill(part.eta(), part.pAbs()); + if(part.id() == 321) h_Particle_Kaon_p_eta_p->Fill(part.eta(), part.pAbs()); + if(part.id() == -321) h_Particle_Kaon_n_eta_p->Fill(part.eta(), part.pAbs()); + if(part.id() == 2212) h_Particle_proton_p_eta_p->Fill(part.eta(), part.pAbs()); + if(part.id() == -2212) h_Particle_proton_n_eta_p->Fill(part.eta(), part.pAbs()); + if(part.id() == -11) h_Particle_Electron_p_eta_p->Fill(part.eta(), part.pAbs()); + if(part.id() == 11) h_Particle_Electron_n_eta_p->Fill(part.eta(), part.pAbs()); + + if(part.id() == 2112) h_Particle_Neutron_eta_p->Fill(part.eta(), part.pAbs()); + if(part.id() == 22) h_Particle_Gamma_eta_p->Fill(part.eta(), part.pAbs()); + + // eta, transverse momentum pT + if(part.id() == 211) h_Particle_pion_p_eta_pT->Fill(part.eta(), part.pT()); + if(part.id() == -211) h_Particle_pion_n_eta_pT->Fill(part.eta(), part.pT()); + if(part.id() == 321) h_Particle_Kaon_p_eta_pT->Fill(part.eta(), part.pT()); + if(part.id() == -321) h_Particle_Kaon_n_eta_pT->Fill(part.eta(), part.pT()); + if(part.id() == 2212) h_Particle_proton_p_eta_pT->Fill(part.eta(), part.pT()); + if(part.id() == -2212) h_Particle_proton_n_eta_pT->Fill(part.eta(), part.pT()); + if(part.id() == -11) h_Particle_Electron_p_eta_pT->Fill(part.eta(), part.pT()); + if(part.id() == 11) h_Particle_Electron_n_eta_pT->Fill(part.eta(), part.pT()); + + if(part.id() == 2112) h_Particle_Neutron_eta_pT->Fill(part.eta(), part.pT()); + if(part.id() == 22) h_Particle_Gamma_eta_pT->Fill(part.eta(), part.pT()); + + // eta, energy + if(part.id() == 211) h_Particle_Pion_p_eta_E->Fill(part.eta(), part.e()); + if(part.id() == -211) h_Particle_Pion_n_eta_E->Fill(part.eta(), part.e()); + if(part.id() == 321) h_Particle_Kaon_p_eta_E->Fill(part.eta(), part.e()); + if(part.id() == -321) h_Particle_Kaon_n_eta_E->Fill(part.eta(), part.e()); + if(part.id() == 2212) h_Particle_Proton_p_eta_E->Fill(part.eta(), part.e()); + if(part.id() == -2212) h_Particle_Proton_n_eta_E->Fill(part.eta(), part.e()); + if(part.id() == -11) h_Particle_Electron_p_eta_E->Fill(part.eta(), part.e()); + if(part.id() == 11) h_Particle_Electron_n_eta_E->Fill(part.eta(), part.e()); + + if(part.id() == 2112) h_Particle_Neutron_eta_E->Fill(part.eta(), part.e()); + if(part.id() == 22) h_Particle_Gamma_eta_E->Fill(part.eta(), part.e()); + + + // read in input particles + //---------------------------------------------------------- + input_particles.push_back(fastjet::PseudoJet(part.px(),part.py(),part.pz(),part.e())); + //cout<<"input_particles add = "<0) + { + if(part.isCharged() && isInTracker(part.eta())) + { + input_particles_meas.push_back(fastjet::PseudoJet(part.px(),part.py(),part.pz(),part.e())); + //cout<<"input_particles_meas add charged = "<Clone("hEne_Norm"); + + + float integral = hEne->Integral(); + + //cout<<"integral = "<Fill(a, b); + + if(b >= 0) + { + hist->Fill(a, 3); + } + if(a >= 0) + { + hist->Fill(3, b); + } + if(a >= 0 && b >= 0) + { + hist->Fill(3, 3); + } + + a = BinIdToCalo.at(sorted.at(0).first); + b = BinIdToCalo.at(sorted.at(1).first); + + //cout<<"a ="< jets) +{ + + + bool anyHcal_jets = false; + int nHCal_jets = 0; + int a = -1; + int b = -1; + + for (unsigned int i = 0; i < jets.size(); i++) { + + fastjet::PseudoJet jet = jets[i]; + + // nHCal + if(-4.14 < jet.eta() && jet.eta() < -1.18) + { + anyHcal_jets = true; + nHCal_jets++; + + if(i==0) a = 0; + if(i==1) b = 0; + } + + // bHCal + if(-1.1 < jet.eta() && jet.eta() < 1.1) + { + anyHcal_jets = true; + + if(i==0) a = 1; + if(i==1) b = 1; + } + + // LFHCal + if(1.2 < jet.eta() && jet.eta() < 4.2) + { + anyHcal_jets = true; + + if(i==0) a = 2; + if(i==1) b = 2; + } + + } + + hist->Fill(a, b); + + if(b >= 0) + { + hist->Fill(a, 3); + } + if(a >= 0) + { + hist->Fill(3, b); + } + if(a >= 0 && b >= 0) + { + hist->Fill(3, 3); + } + + //cout<<"a ="< jets) +{ + + + bool anyHcal_jets = false; + int nHCal_jets = 0; + int a = -1; + int b = -1; + + bool is_nHCal = false; + bool is_bHCal = false; + bool is_LFHCal = false; + + for (unsigned int i = 0; i < jets.size(); i++) { + + fastjet::PseudoJet jet = jets[i]; + + vector constituents = jet.constituents(); + + for (int j = 0; j < constituents.size(); ++j) { + + fastjet::PseudoJet con = constituents[j]; + + if(-4.14 < con.eta() && con.eta() < -1.18) is_nHCal = true; + if(-1.1 < con.eta() && con.eta() < 1.1) is_bHCal = true; + if(1.2 < con.eta() && con.eta() < 4.2) is_LFHCal = true; + } + + + //---------- + + // nHCal + if(is_nHCal && !is_bHCal && !is_LFHCal) + { + if(i==0) a = 0; + if(i==1) b = 0; + } + + // nHCal+bHCal + if(is_nHCal && is_bHCal && !is_LFHCal) + { + if(i==0) a = 1; + if(i==1) b = 1; + } + + // bHCal + if(!is_nHCal && is_bHCal && !is_LFHCal) + { + if(i==0) a = 2; + if(i==1) b = 2; + } + + // bHCal+LFHCal + if(!is_nHCal && is_bHCal && is_LFHCal) + { + if(i==0) a = 3; + if(i==1) b = 3; + } + + // LFHCal + if(!is_nHCal && !is_bHCal && is_LFHCal) + { + if(i==0) a = 4; + if(i==1) b = 4; + } + + } + + hist->Fill(a, b); + + if(b >= 0) + { + hist->Fill(a, 5); + } + if(a >= 0) + { + hist->Fill(5, b); + } + if(a >= 0 && b >= 0) + { + hist->Fill(5, 5); + } + + //cout<<"a ="<> &input, std::vector> &output) +{ + + //int maxId = 0; + //int maxVal = 0.0; + + std::vector> input_copy = input; + + //PrintVec(input_copy); + + sort(input_copy.rbegin(), input_copy.rend(), compFunc); + + //cout<<"Sorted:"< maxVal) + { + maxId = j+1; + maxVal = currVal; + } + } // j + + pair newPair; + newPair.first = maxId; + newPair.second = maxVal; + + + output.push_back(newPair); + input_copy.erase(maxId-1); + + } // i +*/ +} + + + +bool compFunc(pair a, pair b) +{ + return (a.second < b.second); +} + + +void PrintVec(std::vector> input) +{ + + for (int i = 0; i < input.size(); ++i) { + + int a = input.at(i).first; + float b = input.at(i).second; + + cout<\t"; + + } + cout< jets, Particle parton1, Particle parton2, int &jetid1, int &jetid2) +{ + + double r1_min = 1.0; + double r2_min = 1.0; + + int id1_min = -1; + int id2_min = -1; + + for (int i = 0; i < jets.size(); ++i) { + + fastjet::PseudoJet jet = jets[i]; + + double r1 = dR(jet, parton1); + double r2 = dR(jet, parton2); + + if(r1 [OUTDIR]} +OUTDIR=${2:-"$(pwd)/sim_input"} + +SCRIPT_DIR="$(dirname "$0")" + +mkdir -p "$OUTDIR" + +OUTBASE="${OUTDIR}/DiMuon_ep_18x275GeV.${INDEX}" +LOG="${OUTDIR}/DiMuon_ep_18x275GeV_${INDEX}.log" + +"${SCRIPT_DIR}/pythiaDiMuon" "${SCRIPT_DIR}/ep_DiMuon.cmnd" "$OUTBASE" | tee "$LOG" + diff --git a/benchmarks/nhcal_pion_rejection/Snakefile b/benchmarks/nhcal_pion_rejection/Snakefile new file mode 100644 index 00000000..1dec5388 --- /dev/null +++ b/benchmarks/nhcal_pion_rejection/Snakefile @@ -0,0 +1,76 @@ +rule nhcal_pion_rejection_simulate: + output: + "sim_output/nhcal_pion_rejection/sim_{DETECTOR_CONFIG}.{INDEX}.edm4hep.root", + params: + N_EVENTS=1000, + SEED=lambda wildcards: "1" + wildcards.INDEX, + DETECTOR_PATH=os.environ["DETECTOR_PATH"], + DETECTOR_CONFIG=lambda wildcards: wildcards.DETECTOR_CONFIG, + DD4HEP_HASH=get_spack_package_hash("dd4hep"), + NPSIM_HASH=get_spack_package_hash("npsim"), + wildcard_constraints: + INDEX = r"\d+", + cache: True + shell: + """ +exec npsim \ + --compactFile {params.DETECTOR_PATH}/{params.DETECTOR_CONFIG}.xml \ + --numberOfEvents {params.N_EVENTS} \ + --random.seed {params.SEED} \ + --enableGun \ + -v WARNING \ + --gun.particle pi- \ + --gun.thetaMin 120*degree \ + --gun.thetaMax 180*degree \ + --gun.distribution uniform \ + --gun.momentumMin "1*GeV" \ + --gun.momentumMax "18*GeV" \ + --outputFile {output} +""" + +rule nhcal_pion_rejection_recon: + input: + "sim_output/nhcal_pion_rejection/sim_{DETECTOR_CONFIG}.{INDEX}.edm4hep.root" + output: + "sim_output/nhcal_pion_rejection/sim_{DETECTOR_CONFIG}.{INDEX}.eicrecon.edm4hep.root" + params: + DETECTOR_CONFIG = lambda wildcards: wildcards.DETECTOR_CONFIG, + cache: True + shell: + """ +exec env DETECTOR_CONFIG={params.DETECTOR_CONFIG} \ + eicrecon {input} \ + -Ppodio:output_file={output} \ + +""" + +rule nhcal_pion_rejection_combine: + input: + lambda wildcards: expand( + "sim_output/nhcal_pion_rejection/sim_{DETECTOR_CONFIG}.{INDEX}.eicrecon.edm4hep.root", + DETECTOR_CONFIG=wildcards.DETECTOR_CONFIG, + INDEX=range(int(wildcards.N)), + ) + output: + temp("sim_output/nhcal_pion_rejection/sim_{DETECTOR_CONFIG}_{N}files.eicrecon.edm4hep.root"), + wildcard_constraints: + N = r"\d+", + shell: + """ +hadd -f {output} {input} +""" + +rule nhcal_pion_rejection_analysis: + input: + combined="sim_output/nhcal_pion_rejection/sim_{DETECTOR_CONFIG}_{N}files.eicrecon.edm4hep.root", + script="benchmarks/nhcal_pion_rejection/scripts/pion_rejection_analysis.cxx", + output: + png="results/nhcal_pion_rejection/analysis_{DETECTOR_CONFIG}_{N}files.png", + pdf="results/nhcal_pion_rejection/analysis_{DETECTOR_CONFIG}_{N}files.pdf", + params: + DETECTOR_PATH = os.environ["DETECTOR_PATH"], + DETECTOR_CONFIG = lambda wildcards: f"{wildcards.DETECTOR_CONFIG}.xml", + shell: + """ + root -l -b -q '{input.script}+O("{input.combined}","{output.pdf}","{output.png}","{params.DETECTOR_PATH}/{params.DETECTOR_CONFIG}")' +""" \ No newline at end of file diff --git a/benchmarks/nhcal_pion_rejection/config.yml b/benchmarks/nhcal_pion_rejection/config.yml new file mode 100644 index 00000000..f10648a8 --- /dev/null +++ b/benchmarks/nhcal_pion_rejection/config.yml @@ -0,0 +1,28 @@ +sim:nhcal_pion_rejection: + extends: .det_benchmark + stage: simulate + parallel: + matrix: + - INDEX_RANGE: ["0 1", "2 3", "4 5", "6 7", "8 9"] + script: + - snakemake $SNAKEMAKE_FLAGS --retries 1 --cores $MAX_CORES_PER_JOB $(for INDEX in $(seq $INDEX_RANGE); do echo sim_output/nhcal_pion_rejection/sim_epic_full.${INDEX}.edm4hep.root; done) + - snakemake $SNAKEMAKE_FLAGS --retries 1 --cores $MAX_CORES_PER_JOB $(for INDEX in $(seq $INDEX_RANGE); do echo sim_output/nhcal_pion_rejection/sim_epic_full.${INDEX}.eicrecon.edm4hep.root; done) + +bench:nhcal_pion_rejection_analysis: + extends: .det_benchmark + stage: benchmarks + needs: + - "sim:nhcal_pion_rejection" + script: + - snakemake $SNAKEMAKE_FLAGS --cores $MAX_CORES_PER_JOB results/nhcal_pion_rejection/analysis_epic_full_10files.pdf + +collect_results:nhcal_pion_rejection: + extends: .det_benchmark + stage: collect + needs: + - "bench:nhcal_pion_rejection_analysis" + script: + - ls -lrht + - mv results{,_save}/ + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/nhcal_pion_rejection/analysis_epic_full_10files.pdf + - mv results{_save,}/ diff --git a/benchmarks/nhcal_pion_rejection/scripts/pion_rejection_analysis.cxx b/benchmarks/nhcal_pion_rejection/scripts/pion_rejection_analysis.cxx new file mode 100644 index 00000000..5bd18b31 --- /dev/null +++ b/benchmarks/nhcal_pion_rejection/scripts/pion_rejection_analysis.cxx @@ -0,0 +1,792 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include "TLorentzVector.h" + +#include "ROOT/RDataFrame.hxx" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "TROOT.h" +#include "TRandom.h" +#include "TH3.h" + +#include "DD4hep/Detector.h" +#include "DDRec/CellIDPositionConverter.h" +#include "DD4hep/Volumes.h" +#include "DD4hep/Objects.h" +#include "DD4hep/Shapes.h" +#include +#include +#include +#include +#include + +#include "podio/Frame.h" +#include "podio/CollectionBase.h" +#include "podio/ROOTReader.h" +#include "podio/CollectionIDTable.h" +#include "podio/ObjectID.h" + +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/MCParticle.h" + +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4hep/SimCalorimeterHit.h" + +#include "edm4hep/CalorimeterHit.h" +#include "edm4hep/CalorimeterHitCollection.h" + +#include "edm4eic/ClusterCollection.h" +#include "edm4eic/Cluster.h" +#include "edm4eic/ClusterData.h" + +#include "edm4eic/CalorimeterHit.h" +#include "edm4eic/CalorimeterHitCollection.h" + +#include "edm4eic/InclusiveKinematicsCollection.h" +#include "edm4eic/InclusiveKinematics.h" +#include "edm4hep/utils/kinematics.h" +#include "edm4hep/utils/vector_utils.h" +#include "edm4eic/vector_utils_legacy.h" + +#include "edm4eic/Track.h" +#include "edm4eic/TrackSegment.h" +#include "edm4eic/TrackPoint.h" +#include "edm4eic/TrackParameters.h" +#include "edm4eic/TrackSegmentCollection.h" + +#include "edm4eic/ReconstructedParticleCollection.h" +#include "edm4eic/ReconstructedParticle.h" + +#include "edm4eic/MCRecoCalorimeterHitAssociationCollection.h" +#include "edm4eic/MCRecoCalorimeterHitAssociation.h" +#include "edm4eic/MCRecoParticleAssociationCollection.h" +#include "edm4eic/MCRecoParticleAssociation.h" + +using namespace std; +using namespace ROOT; +using namespace TMath; +using namespace edm4hep; + +dd4hep::Detector* det = NULL; + +constexpr double ETA_MIN = -4.14, ETA_MAX = -1.16; + +inline bool inNHCal(double eta) {return (eta >= ETA_MIN && eta <= ETA_MAX);} + +inline string addPrefixAfterSlash(const string& path, + const string& prefix) { + const auto slash = path.find_last_of('/'); + if (slash == string::npos) + return prefix + path; + return path.substr(0, slash + 1) + prefix + path.substr(slash + 1); +} + +void MakeLogBins(double *array, int nbins, double binLo, double binHi) +{ + double logMin = log10(binLo); + double logMax = log10(binHi); + double binWidth = (logMax - logMin) / nbins; + + for (int i = 0; i <= nbins; ++i) { + array[i] = pow(10, logMin + i * binWidth); + } +} + +TLine* makeLine(double x1, double y1, double x2, double y2) { + TLine* l = new TLine(x1, y1, x2, y2); + l->SetLineColor(kRed); + l->SetLineStyle(2); + l->SetLineWidth(2); + l->Draw("same"); + return l; +} + +static TGraph* makeEffGraph_B(TH1D* h_sel, TH1D* h_all, + const char* name, + int minAll = 2, + bool logx = false) +{ + if (!h_sel || !h_all) return nullptr; + + std::unique_ptr h_eff((TH1D*)h_sel->Clone(Form("tmp_%s", name ? name : "eff"))); + h_eff->SetDirectory(nullptr); + h_eff->Sumw2(); + h_eff->Divide(h_sel, h_all, 1.0, 1.0, "B"); + + const int nb = h_all->GetNbinsX(); + + std::unique_ptr g_log; + std::unique_ptr g_lin; + + if (logx) g_log = std::make_unique(); + else g_lin = std::make_unique(); + + int ip = 0; + for (int b = 1; b <= nb; ++b) { + const double all = h_all->GetBinContent(b); + if (all < minAll) continue; + + const double xl = h_all->GetXaxis()->GetBinLowEdge(b); + const double xh = h_all->GetXaxis()->GetBinUpEdge(b); + + const double y = 100.0 * h_eff->GetBinContent(b); + const double ey = 100.0 * h_eff->GetBinError(b); + + if (logx) { + if (xl <= 0.0 || xh <= 0.0) continue; + const double x = sqrt(xl * xh); + const double exl = x - xl; + const double exh = xh - x; + + g_log->SetPoint(ip, x, y); + g_log->SetPointError(ip, exl, exh, ey, ey); + } else { + const double x = h_all->GetBinCenter(b); + const double ex = 0.5 * (xh - xl); + + g_lin->SetPoint(ip, x, y); + g_lin->SetPointError(ip, ex, ey); + } + ++ip; + } + + TGraph* g = logx ? (TGraph*)g_log.release() : (TGraph*)g_lin.release(); + g->SetName(Form("g_%s", name ? name : "eff")); + return g; +} + +inline void drawEffPanel(TH1D* h_all, + TH1D* const h_in[3], + const char* title, + const char* xlabel, + bool logx) +{ + gPad->SetGrid(); + if (logx) gPad->SetLogx(); + + auto* mg = new TMultiGraph(); + auto* leg = new TLegend(0.63, 0.7, 0.88, 0.88); + leg->SetBorderSize(0); + leg->SetFillStyle(0); + leg->SetTextSize(0.04); + + Color_t colors[3] = {kBlue, kRed, kBlack}; + Style_t markers[3] = {20, 21, 22}; + + int added = 0; + + for (int i = 0; i < 3; ++i) { + if (!h_all || !h_in[i]) continue; + auto* g = makeEffGraph_B(h_in[i], h_all, Form("eff_%d", i), 2, logx); + if (!g || g->GetN() == 0) continue; + + g->SetMarkerColor(colors[i]); + g->SetMarkerStyle(markers[i]); + g->SetMarkerSize(1.0); + g->SetLineColor(kBlack); + + mg->Add(g, "PE"); + leg->AddEntry(g, Form("%d mu-in nHCAL", i), "p"); + ++added; + } + + if (h_all && h_in[0] && h_in[1] && h_in[2]) { + std::unique_ptr h_sum((TH1D*)h_in[0]->Clone("h_sum")); + h_sum->SetDirectory(nullptr); + h_sum->Add(h_in[1]); + h_sum->Add(h_in[2]); + + auto* gsum = makeEffGraph_B(h_sum.get(), h_all, "eff_sum", 2, logx); + if (gsum && gsum->GetN() > 0) { + gsum->SetMarkerStyle(25); + gsum->SetMarkerColor(kMagenta + 1); + gsum->SetFillColor(kMagenta + 1); + gsum->SetLineColor(kBlack); + gsum->SetMarkerSize(1.0); + + mg->Add(gsum, "PE"); + leg->AddEntry(gsum, "All eligible", "p"); + ++added; + } + } + + mg->SetTitle(Form("%s;%s;geom. acc. [%%]", title ? title : "", xlabel ? xlabel : "")); + mg->Draw("A"); + + gPad->Update(); + if (mg->GetHistogram()) { + double ymax = mg->GetHistogram()->GetMaximum(); + mg->GetHistogram()->SetMaximum(ymax * 1.35); + } + gPad->Modified(); + gPad->Update(); + + if (logx) { + auto* ax = mg->GetXaxis(); + if (ax) { + double xmin = ax->GetXmin(), xmax = ax->GetXmax(); + if (xmin <= 0) xmin = 1e-12; + mg->GetXaxis()->SetLimits(xmin, xmax); + } + } + + leg->Draw(); +} + +inline TVector3 getPlasticDimensionsCM(dd4hep::Detector& det, + const dd4hep::DDSegmentation::BitFieldCoder* dec, + dd4hep::DDSegmentation::CellID cid, + int slice_idx, + int plastic_slice_value = 3) +{ + const double NaN = numeric_limits::quiet_NaN(); + TVector3 dims(NaN, NaN, NaN); + try { + if (!dec) throw runtime_error("decoder==nullptr"); + int actual_value = dec->get(cid, slice_idx); + cerr << "slice_idx=" << slice_idx << " actual=" << actual_value + << " expected=" << plastic_slice_value << "\n"; + cerr << "types: actual=" << typeid(actual_value).name() + << " expected=" << typeid(plastic_slice_value).name() << "\n"; + cerr << "comparison: " << (actual_value != plastic_slice_value) << "\n"; + cerr << "hex: actual=0x" << hex << actual_value + << " expected=0x" << plastic_slice_value << dec << "\n"; + cerr << "type: " << typeid(dec->get(cid, slice_idx)).name() << "\n"; + if (actual_value != plastic_slice_value) + throw runtime_error("cell is not plastic"); + + auto de = det.volumeManager().lookupDetElement(cid); + if (!de.isValid()) throw runtime_error("lookupDetElement failed"); + + auto pv = de.placement(); + if (!pv.isValid()) throw runtime_error("DetElement has no placement"); + + auto vol = pv.volume(); + if (!vol.isValid()) throw runtime_error("Invalid Volume"); + + auto* box = dynamic_cast(vol.solid().ptr()); + if (!box) throw runtime_error("Solid is not TGeoBBox"); + + dims.SetXYZ(2.0 * box->GetDX(), + 2.0 * box->GetDY(), + 2.0 * box->GetDZ()); + return dims; + } catch (const exception& e) { + cerr << "[WARN] getPlasticThicknessCM: " << e.what() << " (cellID=" << cid << ")\n"; + return dims; + } +} + +inline double getPlasticCenterZ_cm( + const dd4hep::DDSegmentation::BitFieldCoder* decoder, + const dd4hep::rec::CellIDPositionConverter& cellid_converter, + dd4hep::DDSegmentation::CellID cid, + int slice_index, + int plastic_slice_value = 3) +{ + try { + if (!decoder) throw std::runtime_error("decoder==nullptr"); + if (decoder->get(cid, slice_index) != plastic_slice_value) + throw std::runtime_error("cell is not plastic (slice mismatch)"); + return cellid_converter.position(cid).z(); + } catch (const std::exception& e) { + cerr << "[WARN] getPlasticThicknessCM: " << e.what() + << " (cellID=" << cid << ")\n"; + return std::numeric_limits::quiet_NaN();; + } +} + +struct GeomState { + double x, y, z; +}; + +inline double hypot2(double a, double b){ return sqrt(a*a + b*b); } +inline int sgnD(double v){ return (v>0) - (v<0); } + +inline pair +trackXYatZ_fromTwoStates(double x1,double y1,double z1, + double x2,double y2,double z2, + double zTarget) +{ + constexpr double EPSZ = 1e-12; + + const double dz = z2 - z1; + if (abs(dz) < EPSZ) { + const double d1 = abs(zTarget - z1); + const double d2 = abs(zTarget - z2); + if (d1 < d2) return {x1, y1}; + if (d2 < d1) return {x2, y2}; + return {0.5*(x1+x2), 0.5*(y1+y2)}; + } + + const double t = (zTarget - z1) / dz; + const double x = x1 + t * (x2 - x1); + const double y = y1 + t * (y2 - y1); + return {x, y}; +} + +inline pair +trackXYatZ(const GeomState& A, const GeomState& B, double zTarget){ + return trackXYatZ_fromTwoStates( + A.x, A.y, A.z, B.x, B.y, B.z, + zTarget + ); +} + + +int pion_rejection_analysis(const string& filename, string outname_pdf, string outname_png, TString compact_file) { + + gStyle->SetOptStat(0); + podio::ROOTReader reader; + reader.openFile(filename); + unsigned nEvents = reader.getEntries("events"); + cout << "Number of events: " << nEvents << endl; + + det = &(dd4hep::Detector::getInstance()); + det->fromCompact(compact_file.Data()); + det->volumeManager(); + det->apply("DD4hepVolumeManager", 0, 0); + + dd4hep::rec::CellIDPositionConverter cellid_converter(*det); + auto idSpec = det->readout("HcalEndcapNHits").idSpec(); + auto decoder = idSpec.decoder(); + const int slice_index = decoder->index("slice"); + if (slice_index < 0) { + cerr << "ERROR: 'slice' field not found in cell ID spec!" << endl; + return 1; + } + + //double thickness_cm = det->constantAsDouble("HcalEndcapNPolystyreneThickness")* 0.1; + + constexpr int NBINS = 60; + + constexpr double Q2_MIN = 1e-2; + constexpr double Q2_MAX = 1.0; + constexpr double X_MIN = 1e-6; + constexpr double X_MAX = 1e-3; + constexpr double Y_MIN = 0.0; + constexpr double Y_MAX = 1.0; + constexpr double W_MIN = 0.0; + constexpr double W_MAX = 160.0; + + constexpr int Z_NBINS = 100; + constexpr double Z_MIN_MM = -4500.0; + constexpr double Z_MAX_MM = -3600.0; + + constexpr int DXY_NBINS = 120; + constexpr double DXY_MIN_MM = -1000.0; + constexpr double DXY_MAX_MM = 1000.0; + + constexpr int DR_NBINS = 120; + constexpr double DR_MIN_MM = 0.0; + constexpr double DR_MAX_MM = 400.0; + + constexpr double P_MIN_GEV = 0.0; + constexpr double P_MAX_GEV = 25.0; + constexpr int P_NBINS = 50; + + constexpr double E_MIN_GEV = 2e-2; + constexpr double E_MAX_GEV = 10.0; + + constexpr int SIZE = 3; + constexpr double MIP_ENERGY_GEV = 0.002; + constexpr double E_CUTS[SIZE] = {1.5, 1.7, 2.0}; + constexpr double E_THRESH = E_CUTS[2]; + constexpr double LAYER_PROC[SIZE] = {0.6, 0.7, 0.8}; + + double t_cm; + double DR_CUTS_CM[SIZE] = {7.0, 10.0, 13.0}; + double DR_THRESH_MM = 10*DR_CUTS_CM[2]; + int LAYER_MAX = 10; + int LAYER_CUTS[SIZE] = {5, 6, 7}; + int LAYER_THRESH = LAYER_CUTS[0]; + + vector Xbins(NBINS+1), Q2bins(NBINS+1), Ebins(NBINS+1); + MakeLogBins(Q2bins.data(), NBINS, Q2_MIN, Q2_MAX); + MakeLogBins(Xbins.data(), NBINS, X_MIN, X_MAX); + MakeLogBins(Ebins.data(), NBINS, E_MIN_GEV, E_MAX_GEV); + + gStyle->SetTitleSize(0.06, "XYZ"); + gStyle->SetLabelSize(0.06, "XYZ"); + gStyle->SetPadLeftMargin(0.15); + gStyle->SetPadRightMargin(0.15); + gStyle->SetPadBottomMargin(0.15); + gStyle->SetPadTopMargin(0.10); + gROOT->ForceStyle(); + + TH2D* hEtaPt = new TH2D("hEtaPt", "Pion #eta vs p_{T};#eta;p_{T} [GeV]", 100, -6., 6., 100, 0., 7.); + + TH1D* hZ_proj = new TH1D("hZ_proj", "Pion track projections in nHCal; z [mm]; N_{proj}", NBINS, Z_MIN_MM, Z_MAX_MM); + + TH1D* hZ_hits = new TH1D("hZ_hits", "Reconstructed hit z in nHCal; z [mm]; N_{hits}", NBINS, Z_MIN_MM, Z_MAX_MM); + + TH3D* hDxDyZ_layer = new TH3D("hDxDyZ_layer","3D residuals (rec - proj): dx, dy vs z; dx [mm]; dy [mm]; z [mm]", NBINS, DXY_MIN_MM, DXY_MAX_MM, NBINS, DXY_MIN_MM, DXY_MAX_MM, NBINS, Z_MIN_MM, Z_MAX_MM); + + TH2D* hDxDy_all = new TH2D("hDxDy_all", + "Residuals (rec - proj): dx vs dy (all layers); dx [mm]; dy [mm]", + DXY_NBINS, DXY_MIN_MM, DXY_MAX_MM, DXY_NBINS, DXY_MIN_MM, DXY_MAX_MM); + + TH2D* hDrZ_layer = new TH2D("hDrZ_layer","Radial residual (rec - proj) vs z; dr [mm]; z [mm]", NBINS, DR_MIN_MM, DR_MAX_MM, NBINS, Z_MIN_MM, Z_MAX_MM); + + TH1D* hDr_all = new TH1D("hDr_all", + "Radial residual dr = #sqrt{dx^{2}+dy^{2}} (all layers); dr [mm]; N", + DR_NBINS, DR_MIN_MM, DR_MAX_MM); + + TH2D* hE_z = new TH2D("hE_z", "Reconstructed hit energy vs layer z; z_{layer} [mm]; E [GeV]", NBINS, Z_MIN_MM, Z_MAX_MM, NBINS, Ebins.data()); + + TH2D* hEsum_z = new TH2D("hEsum_z", "Layer energy sum vs z (reconstructed); z_{layer} [mm]; E_{sum} [GeV]", + NBINS, Z_MIN_MM, Z_MAX_MM, NBINS, Ebins.data()); + + TH1D* hE = new TH1D("hE", "Reconstructed hit energy; E [GeV]; N_{hits}", 10000, 0.0, 1.0); + + TH1D* hEsum = new TH1D("hEsum", "Layer energy sum (reconstructed); E_{sum} [GeV]; N_{energy}", 10000, 0.0, 1.0); + + TH1D* hP_all_pion = new TH1D("hP_all_pion", "MC pion momentum (pions in nHCal acceptance); p_{MC} [GeV]; N_{pions}", P_NBINS, P_MIN_GEV, P_MAX_GEV); + + TH1D* hP_pass_dr[3] = { + new TH1D(Form("hP_pass_dr%.1fcm",DR_CUTS_CM[0]), Form("Accepted (dr<%.1fcm); p_{MC} [GeV]; N",DR_CUTS_CM[0]), P_NBINS, P_MIN_GEV, P_MAX_GEV), + new TH1D(Form("hP_pass_dr%.1fcm",DR_CUTS_CM[1]), Form("Accepted (dr<%.1fcm); p_{MC} [GeV]; N",DR_CUTS_CM[1]), P_NBINS, P_MIN_GEV, P_MAX_GEV), + new TH1D(Form("hP_pass_dr%.1fcm",DR_CUTS_CM[2]), Form("Accepted (dr<%.1fcm); p_{MC} [GeV]; N",DR_CUTS_CM[2]), P_NBINS, P_MIN_GEV, P_MAX_GEV), + }; + + TH1D* hP_pass_ECut[3] = { + new TH1D(Form("hP_pass_E< %.1f MIP", E_CUTS[0]), Form("Accepted (E< %.1f MIP); p_{MC} [GeV]; N",E_CUTS[0]), P_NBINS, P_MIN_GEV, P_MAX_GEV), + new TH1D(Form("hP_pass_E< %.1f MIP", E_CUTS[1]), Form("Accepted (E< %.1f MIP); p_{MC} [GeV]; N",E_CUTS[1]), P_NBINS, P_MIN_GEV, P_MAX_GEV), + new TH1D(Form("hP_pass_E< %.1f MIP", E_CUTS[2]), Form("Accepted (E< %.1f MIP); p_{MC} [GeV]; N",E_CUTS[2]), P_NBINS, P_MIN_GEV, P_MAX_GEV), + }; + + TH1D* hP_pass_LayerCut[3] = { + new TH1D(Form("hP_pass_Layer< %d layers", LAYER_CUTS[0]), Form("Accepted (Layer < %d layers); p_{MC} [GeV]; N",LAYER_CUTS[0]), P_NBINS, P_MIN_GEV, P_MAX_GEV), + new TH1D(Form("hP_pass_Layer< %d layers", LAYER_CUTS[1]), Form("Accepted (Layer < %d layers); p_{MC} [GeV]; N",LAYER_CUTS[1]), P_NBINS, P_MIN_GEV, P_MAX_GEV), + new TH1D(Form("hP_pass_Layer< %d layers", LAYER_CUTS[2]), Form("Accepted (Layer < %d layers); p_{MC} [GeV]; N",LAYER_CUTS[2]), P_NBINS, P_MIN_GEV, P_MAX_GEV), + }; + + for (unsigned ev = 0; ev < nEvents; ev++) { + auto frameData = reader.readNextEntry(podio::Category::Event); + if (!frameData) continue; + podio::Frame frame(std::move(frameData)); + + auto& mcCol = frame.get("MCParticles"); + auto& recParts = frame.get("ReconstructedParticles"); + auto& projSegs = frame.get("CalorimeterTrackProjections"); + auto& hcalRec = frame.get("HcalEndcapNRecHits"); + auto& assocCol = frame.get("ReconstructedParticleAssociations"); + + if (!mcCol.isValid()) continue; + + vector vPions; + vector vLorentzPions; + for (const auto& p : mcCol) { + if (p.getPDG() != -211 || p.getGeneratorStatus() == 0) continue; + vPions.push_back(p); + TLorentzVector v(p.getMomentum().x, p.getMomentum().y, p.getMomentum().z, p.getEnergy()); + vLorentzPions.push_back(v); + hEtaPt->Fill(v.Eta(), v.Pt()); + if(inNHCal(v.Eta())) {hP_all_pion->Fill(v.P());} + } + + vector matchedRecos; + auto find_associated_reco = [&](const edm4hep::MCParticle& mc)->void { + try { + if (!assocCol.isValid() || assocCol.empty()) return; + + const uint32_t mc_idx = mc.getObjectID().index; + + for (const auto& assoc : assocCol) { + if (assoc.getSimID() == mc_idx) { + uint32_t ridx = assoc.getRecID(); + if (!recParts.isValid() || ridx >= recParts.size()) continue; + auto reco = recParts.at(ridx); + if (reco.isAvailable()) { + matchedRecos.push_back(reco); + } + } + } + } catch (...) {} + }; + + if (!assocCol.isValid() || assocCol.empty()) continue; + for(const auto& p: vPions) if (p.getPDG() == -211) find_associated_reco(p); + + if (!recParts.isValid() || !projSegs.isValid() || !hcalRec.isValid() || recParts.empty() || projSegs.empty() || hcalRec.empty()) continue; + + set uniqueCentersZ10; + map layerData; + for (const auto& hit : hcalRec) { + double z = hit.getPosition().z; + double zBin = round(z); + uint64_t cid = hit.getCellID(); + double zCenter_mm = 10 * getPlasticCenterZ_cm(decoder, cellid_converter, cid, slice_index, /*plastic_slice_value=*/3); + if (!std::isnan(zCenter_mm)) { + int z10 = lround(zCenter_mm * 10.0); + uniqueCentersZ10.insert(z10); + } + layerData[zBin] += hit.getEnergy(); + hZ_hits->Fill(z); + hE_z->Fill(z, hit.getEnergy()); + hE->Fill(hit.getEnergy()); + } + + vector layerCentersZ; + layerCentersZ.reserve(uniqueCentersZ10.size()); + for (int z10 : uniqueCentersZ10) layerCentersZ.push_back(z10 / 10.0); + if(layerCentersZ.size() > LAYER_MAX) LAYER_MAX = layerCentersZ.size(); + + for(size_t n = 0; n < SIZE; n++) LAYER_CUTS[n] = static_cast(LAYER_MAX*LAYER_PROC[n]); + LAYER_THRESH = LAYER_CUTS[0]; + + for (const auto& [zValue, sumEnergy] : layerData) {hEsum_z->Fill(zValue, sumEnergy); hEsum->Fill(sumEnergy);} + + vector allTracks; + for (const auto& R : matchedRecos) { + for (const auto& tr : R.getTracks()) { + if (tr.isAvailable()) allTracks.push_back(tr); + } + } + vector segsTagged; + for (const auto& seg : projSegs) { + auto linkedTr = seg.getTrack(); + if (!linkedTr.isAvailable()) continue; + for (const auto& TT : allTracks) { + if (linkedTr.getObjectID() == TT.getObjectID()) { + segsTagged.push_back(seg); + break; + } + } + } + + for (size_t s = 0; s < segsTagged.size(); ++s) { + + vector segMinDistance(LAYER_MAX, std::numeric_limits::infinity()); + vector segHitEnergy(LAYER_MAX, std::numeric_limits::quiet_NaN()); + vector thickness_cm(LAYER_MAX, std::numeric_limits::quiet_NaN()); + vector count_DrCuts(SIZE, 0); + vector count_ECuts(SIZE, 0); + + GeomState A{}, B{}; + bool haveA = false, haveB = false; + + auto points = segsTagged[s].getPoints(); + + for (const auto& pt : points) { + if (pt.system != 113) continue; + hZ_proj->Fill(pt.position.z); + if (!haveA) { + A.x = pt.position.x; A.y = pt.position.y; A.z = pt.position.z; + haveA = true; + } else if (!haveB) { + B.x = pt.position.x; B.y = pt.position.y; B.z = pt.position.z; + haveB = true; + break; + } + } + + if (!haveA || !haveB) {continue;} + + for (size_t i = 0; i < LAYER_MAX; ++i) { + double best_dr_in_layer = 1e10; + double best_E_in_layer; + double partLayerEnergySum = 0; + double ratio_HitPartLayerEnergy = 0; + dd4hep::DDSegmentation::CellID best_cid_in_layer; + + double zc = layerCentersZ[i]; + auto [X, Y] = trackXYatZ(A, B, zc); + + for (const auto& hit : hcalRec) { + const auto& hp = hit.getPosition(); + + if (fabs(hp.z - zc) > 5.0 /*mm*/) continue; + + const double dx = X - hp.x; + const double dy = Y - hp.y; + const double dr = sqrt(dx*dx + dy*dy); + if(dr < 30*10 /*mm*/) partLayerEnergySum += hit.getEnergy(); + + hDxDyZ_layer->Fill(dx, dy, hp.z); + hDxDy_all->Fill(dx, dy); + hDrZ_layer->Fill(dr, hp.z); + hDr_all->Fill(dr); + + if (dr < best_dr_in_layer) { + best_dr_in_layer = dr; + best_E_in_layer = hit.getEnergy(); + best_cid_in_layer = hit.getCellID(); + } + } + + if (best_dr_in_layer < DR_THRESH_MM) { + segMinDistance[i] = best_dr_in_layer; + segHitEnergy[i] = best_E_in_layer; + thickness_cm[i] = getPlasticDimensionsCM(*det, decoder, best_cid_in_layer, slice_index, 3).z(); + t_cm = thickness_cm[0]; + } + ratio_HitPartLayerEnergy += segHitEnergy[i]/partLayerEnergySum; + if (std::isnan(ratio_HitPartLayerEnergy) || fabs(ratio_HitPartLayerEnergy - 1) >= 0.2) continue; + for(size_t j = 0; j < SIZE; ++j){ + if (segMinDistance[i] < DR_CUTS_CM[j] * 10 /*mm*/){++count_DrCuts[j];} + if (segHitEnergy[i] < (E_CUTS[j] * MIP_ENERGY_GEV * thickness_cm[i] * 100)){++count_ECuts[j];} + } + } + for (int j = 0; j < SIZE; ++j){ + if (count_DrCuts[j] > LAYER_THRESH) hP_pass_dr[j]->Fill(vLorentzPions[s].P()); + if (count_ECuts[j] > LAYER_THRESH) hP_pass_ECut[j]->Fill(vLorentzPions[s].P()); + if (count_DrCuts[SIZE-1] > LAYER_CUTS[j]) hP_pass_LayerCut[j]->Fill(vLorentzPions[s].P()); + } + } + } + + TCanvas* c = new TCanvas("c", "Pion analysis", 1600, 1000); + c->Divide(2,2); + + c->cd(1); hEtaPt->Draw("COLZ"); + + c->SaveAs(outname_pdf.c_str()); + c->SaveAs(outname_png.c_str()); + + TCanvas* c_hitTrack = new TCanvas("c_hitTrack", "Pion hit and track analysis", 1600, 1000); + c_hitTrack->Divide(2,2); + c_hitTrack->cd(1); gPad->SetGrid(); hZ_proj->Draw(); + c_hitTrack->cd(2); gPad->SetGrid(); hZ_hits->Draw(); + c_hitTrack->SaveAs(addPrefixAfterSlash(outname_png, "hitTrack_").c_str()); + c_hitTrack->SaveAs(addPrefixAfterSlash(outname_pdf, "hitTrack_").c_str()); + + TCanvas* c_hitTrackDistance = new TCanvas("c_hitTrackDistance", "Pion hit-track distance analysis", 1600, 1000); + c_hitTrackDistance->Divide(2,2); + c_hitTrackDistance->cd(1); gPad->SetGrid(); hDxDyZ_layer->Draw("COLZ"); + c_hitTrackDistance->cd(2); gPad->SetGrid(); hDxDy_all->Draw("COLZ"); + c_hitTrackDistance->cd(3); gPad->SetGrid(); hDrZ_layer->Draw("COLZ"); + c_hitTrackDistance->cd(4); gPad->SetGrid(); hDr_all->Draw("COLZ"); + c_hitTrackDistance->SaveAs(addPrefixAfterSlash(outname_png, "hitTrackDistance_").c_str()); + c_hitTrackDistance->SaveAs(addPrefixAfterSlash(outname_pdf, "hitTrackDistance_").c_str()); + + TCanvas* c_Eff = new TCanvas("c_Eff", "Pion efficiency analysis", 1600, 1000); + c_Eff->Divide(2,2); + c_Eff->cd(1); gPad->SetGrid(); hP_all_pion->SetLineColor(kBlack); hP_all_pion->Draw(); + c_Eff->cd(2); gPad->SetGrid(); + TH1D* hEff_dr[3]; + for (int idr=0; idr<3; ++idr) { + hEff_dr[idr] = (TH1D*)hP_pass_dr[idr]->Clone( + Form("hEff_dr%.1fcm", DR_CUTS_CM[idr]) + ); + hEff_dr[idr]->SetDirectory(nullptr); + hEff_dr[idr]->SetTitle(Form("Efficiency vs p_{MC} (Layers > %d); p_{MC} [GeV]; efficiency", LAYER_THRESH)); + hEff_dr[idr]->Divide(hP_pass_dr[idr], hP_all_pion, 1, 1, "B"); + } + + hEff_dr[0]->SetLineColor(kBlue); hEff_dr[0]->SetMinimum(0.0); hEff_dr[0]->SetMaximum(0.3); + hEff_dr[1]->SetLineColor(kRed); + hEff_dr[2]->SetLineColor(kGreen+2); + + hEff_dr[0]->Draw("HIST E"); + hEff_dr[1]->Draw("HIST E SAME"); + hEff_dr[2]->Draw("HIST E SAME"); + { auto leg_dr = new TLegend(0.60,0.70,0.88,0.88); + for(int idr=0; idr<3; ++idr) + {leg_dr->AddEntry(hEff_dr[idr],Form("dr < %.1f cm", DR_CUTS_CM[idr]),"l");} + leg_dr->Draw(); + } + + c_Eff->cd(3); gPad->SetGrid(); + TH1D* hEff_E[3]; + for (int ie=0; ie<3; ++ie) { + hEff_E[ie] = (TH1D*)hP_pass_ECut[ie]->Clone( + Form("hEff_E%.1fcm", E_CUTS[ie]) + ); + hEff_E[ie]->SetDirectory(nullptr); + hEff_E[ie]->SetTitle(Form("Efficiency vs p_{MC} (Layers > %d, dr < %.1f cm); p_{MC} [GeV]; efficiency",LAYER_THRESH, DR_THRESH_MM/10)); + hEff_E[ie]->Divide(hP_pass_ECut[ie], hP_all_pion, 1, 1, "B"); + } + + hEff_E[0]->SetLineColor(kBlue); hEff_E[0]->SetMinimum(0.0); hEff_E[0]->SetMaximum(0.3); + hEff_E[1]->SetLineColor(kRed); + hEff_E[2]->SetLineColor(kGreen+2); + + hEff_E[0]->Draw("HIST E"); + hEff_E[1]->Draw("HIST E SAME"); + hEff_E[2]->Draw("HIST E SAME"); + { auto leg_E = new TLegend(0.60,0.70,0.88,0.88); + for(int ie=0; ie<3; ++ie) + {leg_E->AddEntry(hEff_E[ie],Form("E < %.1f MIP", E_CUTS[ie]),"l");} + leg_E->Draw(); + } + + c_Eff->cd(4); gPad->SetGrid(); + TH1D* hEff_Layer[3]; + for (int il=0; il<3; ++il) { + hEff_Layer[il] = (TH1D*)hP_pass_LayerCut[il]->Clone( + Form("hEff_Layer%d", LAYER_CUTS[il]) + ); + hEff_Layer[il]->SetDirectory(nullptr); + hEff_Layer[il]->SetTitle(Form("Efficiency vs p_{MC} (dr < %.1f cm); p_{MC} [GeV]; efficiency", DR_THRESH_MM/10)); + hEff_Layer[il]->Divide(hP_pass_LayerCut[il], hP_all_pion, 1, 1, "B"); + } + + hEff_Layer[0]->SetLineColor(kBlue); hEff_Layer[0]->SetMinimum(0.0); hEff_Layer[0]->SetMaximum(0.3); + hEff_Layer[1]->SetLineColor(kRed); + hEff_Layer[2]->SetLineColor(kGreen+2); + + hEff_Layer[0]->Draw("HIST E"); + hEff_Layer[1]->Draw("HIST E SAME"); + hEff_Layer[2]->Draw("HIST E SAME"); + { auto leg_Layer = new TLegend(0.60,0.70,0.88,0.88); + for(int il=0; il<3; ++il) + {leg_Layer->AddEntry(hEff_Layer[il],Form("L > %d Layers", LAYER_CUTS[il]),"l");} + leg_Layer->Draw(); + } + c_Eff->SaveAs(addPrefixAfterSlash(outname_png, "matching_efficiency_").c_str()); + c_Eff->SaveAs(addPrefixAfterSlash(outname_pdf, "matching_efficiency_").c_str()); + + TCanvas* c_hEnergy = new TCanvas("c_hEnergy", "Pion hit energy analysis", 1600, 1000); + c_hEnergy->Divide(2,2); + + c_hEnergy->cd(1); gPad->SetGrid(); gPad->SetLogy(); hE_z->GetYaxis()->SetMoreLogLabels(); hE_z->Draw("COLZ"); + TProfile* pE_z = hE_z->ProfileX("pE_z"); pE_z->SetLineWidth(3); pE_z->SetLineColor(kRed); pE_z->SetMarkerColor(kRed); + pE_z->SetDirectory(nullptr); pE_z->Draw("SAME"); + + const double Ecut_GeV = MIP_ENERGY_GEV * t_cm * 100; + TLine* cut = new TLine(hE_z->GetXaxis()->GetXmin(), Ecut_GeV, hE_z->GetXaxis()->GetXmax(), Ecut_GeV); + cut->SetLineColor(kRed+1); cut->SetLineStyle(2); cut->SetLineWidth(2); + cut->Draw("SAME"); + auto* leg = new TLegend(0.58, 0.75, 0.84, 0.88); leg->SetTextSize(0.04);leg->AddEntry(cut, Form("E_{cut} = %.3f GeV", Ecut_GeV), "l"); + leg->Draw(); + + c_hEnergy->cd(2); gPad->SetGrid(); gPad->SetLogy(); hEsum_z->GetYaxis()->SetMoreLogLabels(); hEsum_z->Draw("COLZ"); + TProfile* pEsum_z = hEsum_z->ProfileX("pEsum_z"); pEsum_z->SetLineWidth(3); pEsum_z->SetLineColor(kRed); pEsum_z->SetMarkerColor(kRed); + pEsum_z->SetDirectory(nullptr); pEsum_z->Draw("SAME"); + c_hEnergy->cd(3); gPad->SetGrid(); hE->Draw(); + c_hEnergy->cd(4); gPad->SetGrid(); hEsum->Draw(); + + c_hEnergy->SaveAs(addPrefixAfterSlash(outname_png, "hit_energy_").c_str()); + c_hEnergy->SaveAs(addPrefixAfterSlash(outname_pdf, "hit_energy_").c_str()); + + delete c; + delete c_hitTrack; + delete c_hitTrackDistance; + delete c_Eff; + delete c_hEnergy; + + return 0; +} \ No newline at end of file diff --git a/benchmarks/nhcal_sampling_fraction/Snakefile b/benchmarks/nhcal_sampling_fraction/Snakefile new file mode 100644 index 00000000..27acc029 --- /dev/null +++ b/benchmarks/nhcal_sampling_fraction/Snakefile @@ -0,0 +1,70 @@ +def get_total_energy(particle, kinetic_energy): + """Convert kinetic energy to total energy by adding particle mass.""" + masses = { + "neutron": 0.94, + "pi-": 0.14, + "e-": 0.0005 + } + ke = float(kinetic_energy) + mass = masses[particle] + return f"{ke + mass:.2f}" + +rule nhcal_sampling_fraction_simulate: + output: + "sim_output/nhcal_sampling_fraction/{PARTICLE}/Ekin{ENERGY}GeV/sim.edm4hep.root", + params: + N_EVENTS=5000, + DETECTOR_PATH=os.environ["DETECTOR_PATH"], + DETECTOR_CONFIG="epic_backward_hcal_only_sampF.xml", + DD4HEP_HASH=get_spack_package_hash("dd4hep"), + NPSIM_HASH=get_spack_package_hash("npsim"), + TOTAL_ENERGY=lambda wildcards: get_total_energy(wildcards.PARTICLE, wildcards.ENERGY), + cache: True + shell: + """ +set -m +exec npsim \ + --compactFile {params.DETECTOR_PATH}/{params.DETECTOR_CONFIG} \ + --numberOfEvents {params.N_EVENTS} \ + --random.seed $RANDOM \ + --enableGun \ + -v WARNING \ + --gun.particle {wildcards.PARTICLE} \ + --gun.thetaMin 120*degree \ + --gun.thetaMax 180*degree \ + --gun.distribution uniform \ + --gun.energy "{params.TOTAL_ENERGY}*GeV" \ + --outputFile {output} +""" + + +rule nhcal_sampling_fraction_combine: + input: + lambda wildcards: expand( + "sim_output/nhcal_sampling_fraction/{PARTICLE}/Ekin{ENERGY}GeV/sim.edm4hep.root", + ENERGY=["0.5", "0.7", "1.0", "2.0", "5.0", "10.0"], + PARTICLE=["pi-", "neutron", "e-"], + ), + wildcard_constraints: + ENERGY=r"\d+" + output: + f"sim_output/nhcal_sampling_fraction/sim_combined.edm4hep.root", + shell: + """ +hadd -f {output} {input} +""" + +rule nhcal_sampling_fraction_analysis: + input: + combined="sim_output/nhcal_sampling_fraction/sim_combined.edm4hep.root", + script="benchmarks/nhcal_sampling_fraction/scripts/sampling_fraction_analysis.cxx", + output: + pdf="results/nhcal_sampling_fraction/hist_sampf_vs_Ehit.pdf", + png="results/nhcal_sampling_fraction/hist_sampf_vs_Ehit.png", + params: + DETECTOR_PATH=os.environ["DETECTOR_PATH"], + DETECTOR_CONFIG="epic_backward_hcal_only_sampF.xml" + shell: + """ + root -l -b -q '{input.script}+("{input.combined}","{output.pdf}","{output.png}","{params.DETECTOR_PATH}/{params.DETECTOR_CONFIG}")' +""" diff --git a/benchmarks/nhcal_sampling_fraction/config.yml b/benchmarks/nhcal_sampling_fraction/config.yml new file mode 100644 index 00000000..40697007 --- /dev/null +++ b/benchmarks/nhcal_sampling_fraction/config.yml @@ -0,0 +1,29 @@ +sim:nhcal_sampling_fraction: + extends: .det_benchmark + stage: simulate + parallel: + matrix: + - ENERGY: ["0.5", "0.7","1.0", "2.0", "5.0", "10.0"] + PARTICLE: ["pi-", "neutron", "e-"] + script: + - snakemake $SNAKEMAKE_FLAGS --cores $MAX_CORES_PER_JOB sim_output/nhcal_sampling_fraction/${PARTICLE}/Ekin${ENERGY}GeV/sim.edm4hep.root + +bench:nhcal_sampling_fraction: + extends: .det_benchmark + stage: benchmarks + needs: + - "sim:nhcal_sampling_fraction" + script: + - snakemake $SNAKEMAKE_FLAGS --cores 1 results/nhcal_sampling_fraction/hist_sampf_vs_Ehit.pdf + +collect_results:nhcal_sampling_fraction: + extends: .det_benchmark + stage: collect + needs: + - "bench:nhcal_sampling_fraction" + script: + - ls -lrht + - mv results{,_save}/ + - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/nhcal_sampling_fraction/hist_sampf_vs_Ehit.pdf + - mv results{_save,}/ + \ No newline at end of file diff --git a/benchmarks/nhcal_sampling_fraction/scripts/sampling_fraction_analysis.cxx b/benchmarks/nhcal_sampling_fraction/scripts/sampling_fraction_analysis.cxx new file mode 100644 index 00000000..a016ecb0 --- /dev/null +++ b/benchmarks/nhcal_sampling_fraction/scripts/sampling_fraction_analysis.cxx @@ -0,0 +1,289 @@ +#pragma cling load("edm4hep") +#pragma cling load("edm4hepDict") +#pragma cling load("edm4hepUtils") +#pragma cling load("edm4eic") +#pragma cling load("edm4eicDict") + +#include +#include +#include +#include +#include + +#include "ROOT/RDataFrame.hxx" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "TROOT.h" +#include "TRandom.h" +#include "TH3.h" +#include "TLorentzVector.h" +#include "TStyle.h" + +#include "DD4hep/Detector.h" +#include "DDRec/CellIDPositionConverter.h" + +#include +#include +#include "podio/ROOTReader.h" +#include "podio/CollectionIDTable.h" +#include "podio/ObjectID.h" + +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/MCParticle.h" + +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4hep/SimCalorimeterHit.h" + +#include "edm4hep/CalorimeterHitCollection.h" +#include "edm4hep/CalorimeterHit.h" + +#include "edm4eic/ClusterCollection.h" +#include "edm4eic/Cluster.h" + +#include "edm4eic/CalorimeterHitCollection.h" +#include "edm4eic/CalorimeterHit.h" + +#include + +using namespace std; +using namespace edm4hep; + +dd4hep::Detector* det = NULL; + +constexpr double ETA_MIN = -4.16, ETA_MAX = -1.16; + +inline bool inNHCal(double eta) {return (eta >= ETA_MIN && eta <= ETA_MAX);} + +inline string addPrefixAfterSlash(const string& path, + const string& prefix) { + const auto slash = path.find_last_of('/'); + const auto dot = path.find_last_of('.'); + + string extension = (dot != string::npos) ? path.substr(dot) : ""; + + if (slash == string::npos) + return prefix + extension; + + return path.substr(0, slash + 1) + prefix + extension; +} + +int sampling_fraction_analysis(const string &filename, string outname_pdf, string outname_png, TString compact_file) +{ + + podio::ROOTReader *reader = new podio::ROOTReader(); + reader->openFile(filename); + unsigned nEvents = reader->getEntries("events"); + cout << "Number of events: " << nEvents << endl; + + det = &(dd4hep::Detector::getInstance()); + det->fromCompact(compact_file.Data()); + det->volumeManager(); + det->apply("DD4hepVolumeManager", 0, 0); + + dd4hep::rec::CellIDPositionConverter cellid_converter(*det); + auto idSpec = det->readout("HcalEndcapNHits").idSpec(); + auto decoder = idSpec.decoder(); + const int slice_index = decoder->index("slice"); + if (slice_index < 0) { + cerr << "ERROR: 'slice' field not found in cell ID spec!" << endl; + return 1; + } + + constexpr int NBINS = 80; + + constexpr double E_MIN_GEV = 0.0; + constexpr double E_MAX_GEV = 12.0; + constexpr double SAMP_F_MIN = 0.0; + constexpr double SAMP_F_MAX = 1.0; + constexpr double SAMP_F_LOW = 0.2; + + gStyle->SetTitleSize(0.04, "XYZ"); + gStyle->SetLabelSize(0.04, "XYZ"); + gStyle->SetPadLeftMargin(0.20); + gStyle->SetPadRightMargin(0.15); + gStyle->SetPadBottomMargin(0.15); + gStyle->SetPadTopMargin(0.10); + + TH2D *h_sampF_e = new TH2D("h_sampF_e", "nHCal sampling fraction vs. energy (e-); E_{kin} [GeV]; sampling fraction from hits [E_{scint}/E_{hit}]; counts", + NBINS, E_MIN_GEV, E_MAX_GEV, NBINS, SAMP_F_MIN, SAMP_F_MAX); + + TH2D *h_sampF_n = new TH2D("h_sampF_n", "nHCal sampling fraction vs. energy (neutron); E_{kin} [GeV]; sampling fraction from hits [E_{scint}/E_{hit}]; counts", + NBINS, E_MIN_GEV, E_MAX_GEV, NBINS, SAMP_F_MIN, SAMP_F_MAX); + + TH2D *h_sampF_pi = new TH2D("h_sampF_pi", "nHCal sampling fraction vs. energy (#pi-); E_{kin} [GeV]; sampling fraction from hits [E_{scint}/E_{hit}]; counts", + NBINS, E_MIN_GEV, E_MAX_GEV, NBINS, SAMP_F_MIN, SAMP_F_MAX); + + TH2D *h_sampF_e_Ekin = new TH2D("h_sampF_e_Ekin", "nHCal sampling fraction vs. energy kin (e-); E_{kin} [GeV]; sampling fraction from hits [E_{scint}/E_{kin}]; counts", + NBINS, E_MIN_GEV, E_MAX_GEV, NBINS, SAMP_F_MIN, SAMP_F_LOW); + + TH2D *h_sampF_n_Ekin = new TH2D("h_sampF_n_Ekin", "nHCal sampling fraction vs. energy kin (neutron); E_{kin} [GeV]; sampling fraction from hits [E_{scint}/E_{kin}]; counts", + NBINS, E_MIN_GEV, E_MAX_GEV, NBINS, SAMP_F_MIN, SAMP_F_LOW); + + TH2D *h_sampF_pi_Ekin = new TH2D("h_sampF_pi_Ekin", "nHCal sampling fraction vs. energy kin (#pi-); E_{kin} [GeV]; sampling fraction from hits [E_{scint}/E_{kin}]; counts", + NBINS, E_MIN_GEV, E_MAX_GEV, NBINS, SAMP_F_MIN, SAMP_F_LOW); + + for (unsigned ev = 0; ev < nEvents; ev++) + { + double hit_Esum = 0; + double hit_scint_Esum = 0; + double singlePart_Ekin = 0; + + auto frameData = reader->readNextEntry(podio::Category::Event); + if (!frameData) + { + cerr << "Invalid FrameData at event " << ev << endl; + continue; + } + + podio::Frame frame(std::move(frameData)); + + auto& MCParticles_coll = frame.get("MCParticles"); + auto& SimCalorimeterHit_coll = frame.get("HcalEndcapNHits"); + + if (!SimCalorimeterHit_coll.isValid()) + { + cerr << "HcalEndcapNHits collection is not valid!" << endl; + } + if (!MCParticles_coll.isValid()) + { + cerr << "MCParticles collection is not valid!" << endl; + } + + edm4hep::MCParticle mcpart = MCParticles_coll.at(0); + TLorentzVector v(mcpart.getMomentum().x, mcpart.getMomentum().y, mcpart.getMomentum().z, mcpart.getEnergy()); + if(!inNHCal(v.Eta())) continue; + singlePart_Ekin = mcpart.getEnergy() - mcpart.getMass(); + int pdg = mcpart.getPDG(); + + for (const auto& hit : SimCalorimeterHit_coll) + { + const int hit_slice = decoder->get(hit.getCellID(), slice_index); + hit_Esum += hit.getEnergy(); + if(hit_slice == 3) hit_scint_Esum += hit.getEnergy(); + } + + if (pdg == 11) // e- + { + h_sampF_e->Fill(singlePart_Ekin, hit_scint_Esum/hit_Esum); + h_sampF_e_Ekin->Fill(singlePart_Ekin, hit_scint_Esum/singlePart_Ekin); + } + else if (pdg == -211) // pi- + { + h_sampF_pi->Fill(singlePart_Ekin, hit_scint_Esum/hit_Esum); + h_sampF_pi_Ekin->Fill(singlePart_Ekin, hit_scint_Esum/singlePart_Ekin); + } + else if (pdg == 2112) // neutron + { + h_sampF_n->Fill(singlePart_Ekin, hit_scint_Esum/hit_Esum); + h_sampF_n_Ekin->Fill(singlePart_Ekin, hit_scint_Esum/singlePart_Ekin); + } + } + + delete reader; + + h_sampF_e->Sumw2(); + h_sampF_e_Ekin->Sumw2(); + h_sampF_pi->Sumw2(); + h_sampF_pi_Ekin->Sumw2(); + h_sampF_n->Sumw2(); + h_sampF_n_Ekin->Sumw2(); + + TProfile* p_sampF_e = h_sampF_e->ProfileX(); + p_sampF_e ->SetTitle("nHCal sampling fraction vs. energy (e-); E_{kin} [GeV]; sampling fraction from hits [E_{scint}/E_{hit}]"); + TProfile* p_sampF_e_Ekin = h_sampF_e_Ekin->ProfileX(); + p_sampF_e_Ekin ->SetTitle("nHCal sampling fraction vs. energy kin (e-); E_{kin} [GeV]; f"); + TProfile* p_sampF_pi = h_sampF_pi->ProfileX(); + p_sampF_pi ->SetTitle("nHCal sampling fraction vs. energy (#pi-); E_{kin} [GeV]; sampling fraction from hits [E_{scint}/E_{hit}]"); + TProfile* p_sampF_pi_Ekin = h_sampF_pi_Ekin->ProfileX(); + p_sampF_pi_Ekin ->SetTitle("nHCal sampling fraction vs. energy kin (#pi-); E_{kin} [GeV]; f"); + TProfile* p_sampF_n = h_sampF_n->ProfileX(); + p_sampF_n ->SetTitle("nHCal sampling fraction vs. energy (neutron); E_{kin} [GeV]; sampling fraction from hits [E_{scint}/E_{hit}]"); + TProfile* p_sampF_n_Ekin = h_sampF_n_Ekin->ProfileX(); + p_sampF_n_Ekin ->SetTitle("nHCal sampling fraction vs. energy kin (neutron); E_{kin} [GeV]; f"); + + TH1D *h_e = p_sampF_e->ProjectionX("h_e"); + TH1D *h_pi = p_sampF_pi->ProjectionX("h_pi"); + + TH1D *h_e_over_pi = (TH1D*)h_e->Clone("h_e_over_pi"); + h_e_over_pi->SetTitle("e/h ratio;E_{kin} [GeV];e/h"); + h_e_over_pi->Divide(h_e, h_pi, 1, 1); + + + TCanvas *c_h = new TCanvas("canvas_h", "c_h", 1600, 800); + c_h->Divide(2,2); + c_h->cd(1); + h_sampF_e->Draw("COLZ"); + + c_h->cd(2); + h_sampF_pi->Draw("COLZ"); + + c_h->cd(3); + h_sampF_n->Draw("COLZ"); + + c_h->SaveAs(outname_pdf.c_str()); + c_h->SaveAs(outname_png.c_str()); + + TCanvas *c_p = new TCanvas("canvas_p", "c_p", 1600, 800); + c_p->Divide(2,2); + c_p->cd(1); + p_sampF_e->SetLineWidth(3); p_sampF_e->SetLineColor(kRed); p_sampF_e->SetMarkerColor(kRed); + p_sampF_e->Draw(); + c_p->cd(2); + p_sampF_pi->SetLineWidth(3); p_sampF_pi->SetLineColor(kRed); p_sampF_pi->SetMarkerColor(kRed); + p_sampF_pi->Draw(); + c_p->cd(3); + p_sampF_n->SetLineWidth(3); p_sampF_n->SetLineColor(kRed); p_sampF_n->SetMarkerColor(kRed); + p_sampF_n->Draw(); + c_p->cd(4); + h_e_over_pi->Draw(); + + c_p->SaveAs(addPrefixAfterSlash(outname_png, "prof_sampf_vs_Ehit").c_str()); + c_p->SaveAs(addPrefixAfterSlash(outname_pdf, "prof_sampf_vs_Ehit").c_str()); + + TH1D *h_e_Ekin = p_sampF_e_Ekin->ProjectionX("h_e_Ekin"); + TH1D *h_pi_Ekin = p_sampF_pi_Ekin->ProjectionX("h_pi_Ekin"); + + TH1D *h_e_over_pi_Ekin = (TH1D*)h_e_Ekin->Clone("h_e_over_pi_Ekin"); + h_e_over_pi_Ekin->SetTitle("e/#pi ratio;E_{kin} [GeV];e/#pi"); + h_e_over_pi_Ekin->Divide(h_e_Ekin, h_pi_Ekin, 1, 1); + + TCanvas *c_h_Ekin = new TCanvas("canvas_h_Ekin", "c_h_Ekin", 1600, 800); + c_h_Ekin->Divide(2,2); + c_h_Ekin->cd(1); + h_sampF_e_Ekin->Draw("COLZ"); + + c_h_Ekin->cd(2); + h_sampF_pi_Ekin->Draw("COLZ"); + + c_h_Ekin->cd(3); + h_sampF_n_Ekin->Draw("COLZ"); + + c_h_Ekin->SaveAs(addPrefixAfterSlash(outname_png, "hist_sampf_vs_Ekin").c_str()); + c_h_Ekin->SaveAs(addPrefixAfterSlash(outname_pdf, "hist_sampf_vs_Ekin").c_str()); + + TCanvas *c_p_Ekin = new TCanvas("canvas_p_Ekin", "c_p_Ekin", 1600, 800); + c_p_Ekin->Divide(2,2); + c_p_Ekin->cd(1); + p_sampF_e_Ekin->SetLineWidth(3); p_sampF_e_Ekin->SetLineColor(kRed); p_sampF_e_Ekin->SetMarkerColor(kRed); + p_sampF_e_Ekin->Draw(); + c_p_Ekin->cd(2); + p_sampF_pi_Ekin->SetLineWidth(3); p_sampF_pi_Ekin->SetLineColor(kRed); p_sampF_pi_Ekin->SetMarkerColor(kRed); + p_sampF_pi_Ekin->Draw(); + c_p_Ekin->cd(3); + p_sampF_n_Ekin->SetLineWidth(3); p_sampF_n_Ekin->SetLineColor(kRed); p_sampF_n_Ekin->SetMarkerColor(kRed); + p_sampF_n_Ekin->Draw(); + c_p_Ekin->cd(4); + h_e_over_pi_Ekin->Draw(); + + c_p_Ekin->SaveAs(addPrefixAfterSlash(outname_png, "prof_sampf_vs_Ekin").c_str()); + c_p_Ekin->SaveAs(addPrefixAfterSlash(outname_pdf, "prof_sampf_vs_Ekin").c_str()); + + return 0; +}