Benchmarks for nhcal

.gitlab-ci.yml

@@ -124,64 +124,71 @@ get_data:
       - 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/ecal_gaps/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/barrel_hcal/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/timing/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/backgrounds/config.yml'
+  # - local: 'benchmarks/backwards_ecal/config.yml'
+  # - local: 'benchmarks/beamline/config.yml'
+  # - local: 'benchmarks/calo_pid/config.yml'
+  # - local: 'benchmarks/ecal_gaps/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/barrel_hcal/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/timing/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'
+
 deploy_results:
   allow_failure: true
   stage: deploy
   needs:
-    - "collect_results:backgrounds"
-    - "collect_results:backwards_ecal"
-    - "collect_results:barrel_ecal"
-    - "collect_results:barrel_hcal"
-    - "collect_results:beamline"
-    - "collect_results:calo_pid"
-    - "collect_results:ecal_gaps"
-    - "collect_results:lfhcal"
-    - "collect_results:material_scan"
-    - "collect_results:pid"
-    - "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:backgrounds"
+    # - "collect_results:backwards_ecal"
+    # - "collect_results:barrel_ecal"
+    # - "collect_results:barrel_hcal"
+    # - "collect_results:beamline"
+    # - "collect_results:calo_pid"
+    # - "collect_results:ecal_gaps"
+    # - "collect_results:lfhcal"
+    # - "collect_results:material_scan"
+    # - "collect_results:pid"
+    # - "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"
   script:
     - snakemake $SNAKEMAKE_FLAGS --cores 1 results/metadata.json
     - find results -print | sort | tee summary.txt

Snakefile

@@ -30,28 +30,31 @@ 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/ecal_gaps/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/backgrounds/Snakefile"
+# include: "benchmarks/backwards_ecal/Snakefile"
+# include: "benchmarks/barrel_ecal/Snakefile"
+# include: "benchmarks/beamline/Snakefile"
+# include: "benchmarks/calo_pid/Snakefile"
+# include: "benchmarks/ecal_gaps/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"
 
 use_s3 = config["remote_provider"].lower() == "s3"
 use_xrootd = config["remote_provider"].lower() == "xrootd"

benchmarks/nhcal_acceptance/Snakefile

@@ -0,0 +1,57 @@
+import os
+
+DETECTOR_PATH = os.environ["DETECTOR_PATH"]
+
+rule nhcal_acceptance_simulate:
+    output:
+        "sim_output/nhcal_acceptance/E{ENERGY}GeV/sim_{DETECTOR_CONFIG}.{INDEX}.edm4hep.root",
+    params:
+        N_EVENTS=10000,
+    shell:
+        """
+set -m
+mkdir -p sim_output/nhcal_acceptance/E{wildcards.ENERGY}GeV
+exec ddsim \
+    --compactFile {DETECTOR_PATH}/{wildcards.DETECTOR_CONFIG}.xml \
+    --numberOfEvents {params.N_EVENTS} \
+    --random.seed $RANDOM \
+    --enableGun \
+    -v WARNING \
+    --gun.particle pi- \
+    --gun.thetaMin 120*degree \
+    --gun.thetaMax 180*degree \
+    --gun.distribution uniform \
+    --gun.energy "{wildcards.ENERGY}*GeV" \
+    --outputFile {output}
+"""
+
+
+rule nhcal_acceptance_combine:
+    input:
+        lambda wildcards: expand(
+            "sim_output/nhcal_acceptance/E{ENERGY}GeV/sim_{DETECTOR_CONFIG}.{INDEX:02d}.edm4hep.root", 
+            DETECTOR_CONFIG=wildcards.DETECTOR_CONFIG,
+            ENERGY=wildcards.ENERGY,
+            INDEX=range(int(wildcards.N)),
+        ),
+    wildcard_constraints:
+        N=r"\d+",
+        ENERGY=r"\d+"
+    output:
+        f"sim_output/nhcal_acceptance/sim_{{DETECTOR_CONFIG}}_E{{ENERGY}}GeV_combined_{{N}}files.edm4hep.root",
+    shell:
+        """
+hadd -f {output} {input} 
+"""
+
+rule nhcal_acceptance_analysis:
+    input:
+        combined="sim_output/nhcal_acceptance/sim_{DETECTOR_CONFIG}_E{ENERGY}GeV_combined_{N}files.edm4hep.root",
+        script="benchmarks/nhcal_acceptance/scripts/acceptance_analysis.cxx",
+    output:
+        png=f"results/nhcal_acceptance/analysis_{{DETECTOR_CONFIG}}_E{{ENERGY}}GeV_combined_{{N}}files.png",
+        pdf=f"results/nhcal_acceptance/analysis_{{DETECTOR_CONFIG}}_E{{ENERGY}}GeV_combined_{{N}}files.pdf",
+    shell:
+        """
+    root -l -b -q '{input.script}("{input.combined}","{output.pdf}","{output.png}")'
+"""

benchmarks/nhcal_acceptance/config.yml

@@ -0,0 +1,34 @@
+sim:nhcal_acceptance:
+  extends: .det_benchmark
+  stage: simulate
+  parallel:
+    matrix:
+      - ENERGY: ["1GeV", "5GeV", "10GeV"]
+        INDEX: ["00","01","02","03","04","05","06","07","08","09"]
+  script:
+    - snakemake --cores 5 sim_output/nhcal_acceptance/E${ENERGY}/sim_epic_backward_hcal_only.${INDEX}.edm4hep.root
+
+bench:nhcal_acceptance_analysis:
+  extends: .det_benchmark
+  stage: benchmarks
+  needs:
+    - "sim:nhcal_acceptance"
+  parallel:
+    matrix:
+      - ENERGY: ["1GeV", "5GeV", "10GeV"]
+  script:
+    - snakemake --cores 1 results/nhcal_acceptance/analysis_epic_backward_hcal_only_E${ENERGY}_combined_10files.pdf
+
+collect_results:nhcal_acceptance:
+  extends: .det_benchmark
+  stage: collect
+  needs: 
+    - "bench:nhcal_acceptance_analysis"
+  parallel:
+    matrix:
+      - ENERGY: ["1GeV", "5GeV", "10GeV"]
+  script:
+    - ls -lrht
+    - mv results{,_save}/ 
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output results/nhcal_acceptance/analysis_epic_backward_hcal_only_E${ENERGY}_combined_10files.pdf
+    - mv results{_save,}/

benchmarks/nhcal_acceptance/scripts/acceptance_analysis.cxx

@@ -0,0 +1,96 @@
+#include <TChain.h>
+#include <TTreeReader.h>
+#include <TTreeReaderArray.h>
+#include <TTreeReaderValue.h>
+#include <TH2D.h>
+#include <TCanvas.h>
+#include <TFile.h>
+#include <TMath.h>
+#include <set>
+#include <iostream>
+#include <TString.h>
+
+using namespace std;
+
+int acceptance_analysis(TString filename, string outname_pdf, string outname_png) 
+{
+    TChain *chain = new TChain("events");
+    chain->Add(filename);
+
+    TTreeReader reader(chain);
+
+    TTreeReaderArray<int> mc_pdg(reader, "MCParticles.PDG");
+    TTreeReaderArray<int> mc_genStatus(reader, "MCParticles.generatorStatus");
+    TTreeReaderArray<double> mc_px(reader, "MCParticles.momentum.x");
+    TTreeReaderArray<double> mc_py(reader, "MCParticles.momentum.y");
+    TTreeReaderArray<double> mc_pz(reader, "MCParticles.momentum.z");
+
+    TTreeReaderArray<int> contrib_particle_idx(reader, "_HcalEndcapNHitsContributions_particle.index");
+    TTreeReaderArray<unsigned int> contrib_particle_cid(reader, "_HcalEndcapNHitsContributions_particle.collectionID");
+
+    int nEtaBins = 100;
+    int nPhiBins = 100;
+    double etaMin = -5, etaMax = 0;
+
+    TH2D* hEtaPhiAll = new TH2D("hEtaPhiAll", "All #pi- (status==1); #eta[1]; #phi[rad]",
+                                nEtaBins, etaMin, etaMax, nPhiBins, -TMath::Pi(), TMath::Pi());
+
+    TH2D* hEtaPhiDetected = new TH2D("hEtaPhiDetected", "#pi- detected in nHCal; #eta[1]; #phi[rad]",
+                                     nEtaBins, etaMin, etaMax, nPhiBins, -TMath::Pi(), TMath::Pi());
+
+    while (reader.Next()) 
+    {
+        map<int, pair<double, double>> pi_minus_eta_phi; 
+        set<int> detected;
+        for (size_t i = 0; i < mc_pdg.GetSize(); ++i) 
+        {
+            if (mc_pdg[i] == -211 && mc_genStatus[i] == 1) 
+            {
+                float px = mc_px[i];
+                float py = mc_py[i];
+                float pz = mc_pz[i];
+
+                float p = sqrt(px * px + py * py + pz * pz);
+                float eta = 0.5 * log((p + pz) / (p - pz + 1e-8));
+                float phi = atan2(py, px);
+
+                hEtaPhiAll->Fill(eta, phi);
+                pi_minus_eta_phi[i] = make_pair(eta, phi);
+            }
+        }
+
+        for (size_t i = 0; i < contrib_particle_idx.GetSize(); i++) {
+            int idx = contrib_particle_idx[i];
+            if (pi_minus_eta_phi.count(idx)) {
+                detected.insert(idx);
+            }
+        }
+
+        for (auto idx : detected) {
+            auto [eta, phi] = pi_minus_eta_phi[idx];
+            hEtaPhiDetected->Fill(eta, phi);
+        }
+    }
+
+    TH2D* hAcceptance = (TH2D*)hEtaPhiAll->Clone("hAcceptance");
+    hAcceptance->Divide(hEtaPhiDetected);
+    hAcceptance->SetTitle("#pi- detected/All");
+    hAcceptance->SetMinimum(0);
+    hAcceptance->SetMaximum(1);
+
+    TCanvas *canvas = new TCanvas("canvas", "pi- All", 1600, 600);
+    canvas->Divide(3,1);
+    canvas->cd(1);
+    hEtaPhiAll->Draw("COLZ");
+    canvas->cd(2);
+    hEtaPhiDetected->Draw("COLZ");
+    canvas->cd(3);
+    hAcceptance->Draw("COLZ");
+
+    canvas->SaveAs(outname_pdf.c_str());
+    canvas->SaveAs(outname_png.c_str());
+
+    return 0;
+}
+
+