first commit of LFHCAL benchmarks

.gitlab-ci.yml

@@ -132,6 +132,7 @@ include:
   - 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_photon/config.yml'

Snakefile

@@ -7,6 +7,7 @@ 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/insert_muon/Snakefile"
 include: "benchmarks/zdc_lambda/Snakefile"
 include: "benchmarks/zdc_lyso/Snakefile"

benchmarks/lfhcal/LFHCAL_Performance.C

@@ -0,0 +1,156 @@
+// Code to extract the Tracking Performances
+// Shyam Kumar; INFN Bari, Italy
+// [email protected]; [email protected]
+
+#include "TGraphErrors.h"
+#include "TF1.h"
+#include "TRandom.h"
+#include "TCanvas.h"
+#include "TLegend.h"
+#include "TMath.h"
+#include "TVector3.h"
+
+#define mpi 0.139  // 1.864 GeV/c^2
+
+void LFHCAL_Performance(TString filename="tracking_output",TString particle="pi-", double mom=0.1, Double_t pTcut = 0.0, TString name = "")
+{
+
+  // style of the plot
+   gStyle->SetPalette(1);
+   gStyle->SetOptTitle(1);
+   gStyle->SetTitleOffset(.85,"X");gStyle->SetTitleOffset(.85,"Y");
+   gStyle->SetTitleSize(.05,"X");gStyle->SetTitleSize(.05,"Y");
+   gStyle->SetLabelSize(.04,"X");gStyle->SetLabelSize(.04,"Y");
+   gStyle->SetHistLineWidth(2);
+   gStyle->SetOptFit(1);
+   gStyle->SetOptStat(1);
+
+   TString dir = "";
+   TString dist_dir_mom = "mom_resol";
+
+   bool debug=true;	  
+  // Tree with reconstructed tracks
+   const int nbins_eta = 5;
+   int nfiles = 100; 
+   double eta[nbins_eta+1]={1.2,1.5,2,2.5,3,3.5};
+   double pt[nbins_eta+1]={0.5,1.0,2.0,5.0,10.0,20.1};
+   TH1D *histp[nbins_eta]; 
+
+
+   for (int i=0; i<nbins_eta; i++){
+   histp[i] = new TH1D(Form("hist_etabin%d",i),Form("hist_etabin%d",i),600,-1,1);
+   histp[i]->SetTitle(Form("%1.1f < #eta < %1.1f && p = %1.1f ",eta[i],eta[i+1],mom));
+   histp[i]->SetName(Form("hist_mom_%1.1f_%1.1f_pmax_%1.1f",mom,eta[i],eta[i+1]));
+   }
+
+   TFile* file = TFile::Open(filename.Data());
+   if (!file) {printf("file not found !!!"); return;}
+   TTreeReader myReader("events", file); // name of tree and file
+   if (debug) cout<<"Filename: "<<file->GetName()<<"\t NEvents: "<<myReader.GetEntries()<<endl;
+
+   // MC and Reco information 
+   TTreeReaderArray<Float_t> charge(myReader, "MCParticles.charge"); 
+   TTreeReaderArray<Double_t> vx_mc(myReader, "MCParticles.vertex.x"); 
+   TTreeReaderArray<Double_t> vy_mc(myReader, "MCParticles.vertex.y"); 
+   TTreeReaderArray<Double_t> vz_mc(myReader, "MCParticles.vertex.z");
+   TTreeReaderArray<Float_t> px_mc(myReader, "MCParticles.momentum.x"); 
+   TTreeReaderArray<Float_t> py_mc(myReader, "MCParticles.momentum.y"); 
+   TTreeReaderArray<Float_t> pz_mc(myReader, "MCParticles.momentum.z");
+   TTreeReaderArray<Int_t> status(myReader, "MCParticles.generatorStatus"); 
+   TTreeReaderArray<Int_t> pdg(myReader, "MCParticles.PDG");
+
+   TTreeReaderArray<Float_t> pe_lc(myReader, "LFHCALClusters.energy"); 
+   TTreeReaderArray<Float_t> px_lc(myReader, "LFHCALClusters.position.x"); 
+   TTreeReaderArray<Float_t> py_lc(myReader, "LFHCALClusters.position.y"); 
+   TTreeReaderArray<Float_t> pz_lc(myReader, "LFHCALClusters.position.z"); 
+   TTreeReaderArray<Float_t> pe_ec(myReader, "EcalEndcapPClusters.energy"); 
+   TTreeReaderArray<Float_t> px_ec(myReader, "EcalEndcapPClusters.position.x"); 
+   TTreeReaderArray<Float_t> py_ec(myReader, "EcalEndcapPClusters.position.y"); 
+   TTreeReaderArray<Float_t> pz_ec(myReader, "EcalEndcapPClusters.position.z"); 
+
+  int count =0;
+  int matchId = 1; // Always matched track assigned the index 0 
+  while (myReader.Next()) 
+    {
+      std::cout << "events = " << count++ << std::endl;
+      for (int j = 0; j < pdg.GetSize(); ++j)
+	{
+	  if (status[j] !=1 && pdg.GetSize()!=1) continue;
+	  Double_t pzmc = pz_mc[j];  
+	  Double_t ptmc = sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]); 
+	  Double_t pmc = sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]+pz_mc[j]*pz_mc[j]); // 1./(q/p); similar to prec
+	  Double_t etamc = -1.0*TMath::Log(TMath::Tan((TMath::ACos(pzmc/fabs(pmc)))/2));
+	  Double_t phimc = TMath::ATan2(py_mc[j],px_mc[j]);
+	  std::cout << "neutron p=" << pmc << " pt=" << ptmc << std::endl;
+
+	  if (fabs(ptmc) < pTcut) continue;
+
+	  float l_px_tot=0;
+	  float l_py_tot=0;
+	  float l_pz_tot=0;
+	  float l_e_tot=0;
+
+	  std::cout << "LFHCAL nclus=" << px_lc.GetSize() << " ECAL nclus=" << pe_ec.GetSize() << std::endl;
+	  for (int jl = 0;jl<px_lc.GetSize();jl++)
+	    {
+	      float e = pe_lc[jl];
+	      TVector3 v(px_lc[jl],py_lc[jl],pz_lc[jl]);
+	      v.Print();
+	      float eta = v.PseudoRapidity();
+	      float phi = v.Phi();
+	      float pt = e/cosh(eta);
+	      std::cout << "LFHCAL clus: e=" << e << " eta=" << eta << " pt=" << pt << std::endl;
+	      l_e_tot += e;
+	      l_px_tot += pt*cos(phi);
+	      l_py_tot += pt*sin(phi);
+	      l_pz_tot += pt*sinh(eta);
+	    }
+
+	  float e_px_tot=0;
+	  float e_py_tot=0;
+	  float e_pz_tot=0;
+	  float e_e_tot=0;
+
+	  for (int je = 0;je<px_ec.GetSize();je++)
+	    {
+	      float e = pe_ec[je];
+	      TVector3 v(px_ec[je],py_ec[je],pz_ec[je]);
+	      float eta = v.PseudoRapidity();
+	      float phi = v.Phi();
+	      float pt = e/cosh(eta);
+	      std::cout << "ECAL clus: e=" << e << " eta=" << eta << " pt=" << pt << std::endl;
+	      e_e_tot += e;
+	      e_px_tot += pt*cos(phi);
+	      e_py_tot += pt*sin(phi);
+	      e_pz_tot += pt*sinh(eta);
+	    }
+
+	  std::cout << "LFHCAL e=" <<l_e_tot << " ECAL e=" << e_e_tot << std::endl;
+	  float px_tot = l_px_tot+e_px_tot;
+	  float py_tot = l_py_tot+e_py_tot;
+	  float pz_tot = l_pz_tot+e_pz_tot;
+	  float e_tot = l_e_tot+e_e_tot;
+
+	  float prec = sqrt(px_tot*px_tot+py_tot*py_tot+pz_tot*pz_tot);
+	  float ptrec = sqrt(px_tot*px_tot+py_tot*py_tot);
+	  float pzrec = pz_tot;
+
+	  float p_resol = (e_tot-pmc)/pmc;
+	  std::cout << "p_resol = " << p_resol << std::endl;
+	  for (int ibin=0; ibin<nbins_eta; ++ibin){ 
+	    if(etamc>eta[ibin] && etamc<eta[ibin+1]) histp[ibin]->Fill(p_resol); 
+	  }
+	} // Generated Tracks  
+
+    }// event loop ends    
+
+   TFile *fout_mom = new TFile(Form("%s/mom/lfhcal_mom_%1.1f_%s_%s.root",particle.Data(),mom,dist_dir_mom.Data(),particle.Data()),"recreate");
+   fout_mom->cd();
+   for (int ibin=0; ibin<nbins_eta; ++ibin) histp[ibin]->Write();
+   fout_mom->Close();
+
+}
+
+
+
+

benchmarks/lfhcal/README.md

@@ -0,0 +1,11 @@
+Detector Benchmark for LFHCAL
+===============================
+
+## Overview
+This benchmark generates events with a range of particle species at fixed energies, all aimed at, or forward of, the LFHCAL, and produces basic plots of clustering performance, both in isolation and combined with the FEMCAL.
+
+
+## Contacts
+[Peter Steinberg]([email protected])
+
+

benchmarks/lfhcal/Snakefile

@@ -0,0 +1,145 @@
+rule lfhcal_sim:
+    input:
+        steering_file="EPIC/EVGEN/SINGLE/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.steer",
+        warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root",
+    output:
+        "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.edm4hep.root",
+    log:
+        "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.edm4hep.root.log",
+    wildcard_constraints:
+        ENERGY="[0-9]+[kMG]eV",
+	PARTICLE="(neutron|pi-|gamma)",
+        PHASE_SPACE="3to50deg",
+        INDEX="\d{4}",
+    params:
+        N_EVENTS=1000
+    shell:
+        """
+ddsim \
+  --runType batch \
+  --enableGun \
+  --steeringFile "{input.steering_file}" \
+  --random.seed 1{wildcards.INDEX} \
+  --filter.tracker edep0 \
+  -v INFO \
+  --numberOfEvents {params.N_EVENTS} \
+  --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
+  --outputFile {output}
+"""
+
+
+rule lfhcal_recon:
+    input:
+        "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.edm4hep.root",
+    output:
+        "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.eicrecon.tree.edm4eic.root",
+    log:
+        "sim_output/lfhcal/{DETECTOR_CONFIG}/{PARTICLE}/{ENERGY}/{PHASE_SPACE}/{PARTICLE}_{ENERGY}_{PHASE_SPACE}.{INDEX}.eicrecon.tree.edm4eic.root.log",
+    wildcard_constraints:
+        INDEX="\d{4}",
+    shell: """
+env DETECTOR_CONFIG={wildcards.DETECTOR_CONFIG} \
+  eicrecon {input} -Ppodio:output_file={output} \
+  -Ppodio:output_collections=MCParticles,ReconstructedParticles,LFHCALTruthClusters,LFHCALClusters,LFHCALHits,EcalEndcapPTruthClusters,EcalEndcapPClusters,EcalEndcapPHits
+"""
+
+rule lfhcal_at_momentum:
+    input:
+        script="benchmarks/lfhcal/LFHCAL_Performance.C",
+        # TODO pass as a file list?
+        sim=lambda wildcards:
+          expand(
+              "sim_output/lfhcal/{DETECTOR_CONFIG}/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{{PARTICLE}}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root",
+              DETECTOR_CONFIG="epic_craterlake",
+              ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV",
+              PHASE_SPACE=["3to50deg"],
+              INDEX=range(1),
+          )
+          if wildcards.CAMPAIGN == "local" else
+          ancient(expand(
+              "EPIC/RECO/{{CAMPAIGN}}/epic_craterlake/SINGLE/{{PARTICLE}}/{ENERGY}/{PHASE_SPACE}/{{PARTICLE}}_{ENERGY}_{PHASE_SPACE}.{INDEX:04d}.eicrecon.tree.edm4eic.root",
+              DETECTOR_CONFIG="epic_craterlake",
+              ENERGY=f"{float(wildcards.MOMENTUM):.0f}GeV" if float(wildcards.MOMENTUM) >= 1 else f"{float(wildcards.MOMENTUM) * 1000:.0f}MeV",
+              PHASE_SPACE=["3to50deg"],
+              INDEX=range(1),
+          )),
+    output:
+        "{CAMPAIGN}/{PARTICLE}/mom/lfhcal_mom_{MOMENTUM}_mom_resol_{PARTICLE}.root",
+        combined_root=temp("{CAMPAIGN}/lfhcal_sim_{MOMENTUM}_{PARTICLE}.root"),
+    shell:
+        """
+hadd {output.combined_root} {input.sim}
+cd {wildcards.CAMPAIGN}
+root -l -b -q ../{input.script}'("../{output.combined_root}", "{wildcards.PARTICLE}", {wildcards.MOMENTUM}, 0.15)'
+"""
+
+rule lfhcal_summary_at_eta:
+    input:
+        expand(
+            [
+                "{{CAMPAIGN}}/{{PARTICLE}}/mom/lfhcal_mom_{MOMENTUM:.1f}_mom_resol_{{PARTICLE}}.root",
+            ],
+            MOMENTUM=[0.5, 1.0, 2.0, 5.0, 10.0, 20.0],
+        ),
+        script="benchmarks/lfhcal/doCompare_widebins_mom.C",
+    output:
+        "{CAMPAIGN}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{PARTICLE}_{ETA_MIN}_eta_{ETA_MAX}.png",
+        "{CAMPAIGN}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{PARTICLE}_{ETA_MIN}_eta_{ETA_MAX}.root",
+    shell:
+        """
+if [[ "{wildcards.CAMPAIGN}" == "local" ]]; then
+	set +e
+	EPIC_VERSION="${{DETECTOR_VERSION:-}}"
+	EICRECON_VERSION="$(eicrecon -v | sed -n -e 's/.*\(v[0-9\.]\+\).*/\\1/p')"
+	# Legacy detection
+	: ${{EPIC_VERSION:="$(echo $DETECTOR_PATH | sed -n -e 's/.*epic-\([^-/]\+\).*/\\1/p')"}}
+	set -e
+
+	echo "ePIC version: $EPIC_VERSION"
+	echo "EICrecon version: $EICRECON_VERSION"
+	EXTRA_LEGEND="ePIC $EPIC_VERSION / EICrecon $EICRECON_VERSION"
+else
+	EXTRA_LEGEND="ePIC Simulation {wildcards.CAMPAIGN}"
+fi
+cd {wildcards.CAMPAIGN}
+root -l -b -q ../{input.script}'("{wildcards.PARTICLE}", {wildcards.ETA_MIN}, {wildcards.ETA_MAX}, 1., true, "'"$EXTRA_LEGEND"'")'
+"""
+
+LFHCAL_ETA_BINS = [1.2,1.5,2,2.5,3,3.5]
+
+rule lfhcal:
+    input:
+        lambda wildcards: expand(
+            [
+	      "{{CAMPAIGN}}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{PARTICLE}_{ETA_BIN}.png",
+	      "{{CAMPAIGN}}/Final_Results/{PARTICLE}/mom/lfhcal_mom_resol_{PARTICLE}_{ETA_BIN}.root",
+            ],
+	    ETA_BIN=[f"{eta_min:.1f}_eta_{eta_max:.1f}" for eta_min, eta_max in zip(LFHCAL_ETA_BINS[:-1], LFHCAL_ETA_BINS[1:])],
+	    PARTICLE=["neutron", "pi-", "gamma"] if wildcards.CAMPAIGN == "local" else ["pi-"],
+        )
+    output:
+        directory("results/lfhcal/{CAMPAIGN}/")
+    shell:
+        """
+mkdir {output}
+cp {input} {output}
+"""
+
+
+rule lfhcal_local:
+    input:
+        "results/lfhcal/local",
+
+
+rule lfhcal_campaigns:
+    input:
+        expand(
+            "results/lfhcal/{CAMPAIGN}",
+            CAMPAIGN=[
+                "23.10.0",
+                "23.11.0",
+                "23.12.0",
+                "24.03.1",
+                "24.04.0",
+            ],
+        )

benchmarks/lfhcal/config.yml

@@ -0,0 +1,48 @@
+sim:lfhcal:
+  extends: .det_benchmark
+  stage: simulate
+  parallel:
+    matrix:
+      - PARTICLE: ["neutron", "pi-"]
+        MOMENTUM: ["500MeV", "1GeV", "2GeV", "5GeV", "10GeV", "20GeV"]
+  script:
+    - |
+      snakemake --cache --cores 5 \
+        sim_output/lfhcal/epic_craterlake/${PARTICLE}/${MOMENTUM}/3to50deg/${PARTICLE}_${MOMENTUM}_3to50deg.0000.eicrecon.tree.edm4eic.root
+
+bench:lfhcal:
+  extends: .det_benchmark
+  stage: benchmarks
+  needs:
+    - ["sim:lfhcal"]
+  script:
+    - snakemake --cores 3 lfhcal_local
+
+collect_results:lfhcal:
+  extends: .det_benchmark
+  stage: collect
+  needs:
+    - "bench:lfhcal"
+  script:
+    - ls -lrht
+    - mv results{,_save}/ # move results directory out of the way to preserve it
+    - snakemake --cores 1 --delete-all-output lfhcal_local
+    - mv results{_save,}/
+
+bench:lfhcal_campaigns:
+  extends: .det_benchmark
+  stage: benchmarks
+  #when: manual
+  script:
+    - snakemake --cores 1 lfhcal_campaigns
+
+collect_results:lfhcal_campaigns:
+  extends: .det_benchmark
+  stage: collect
+  needs:
+    - "bench:lfhcal_campaigns"
+  script:
+    - ls -lrht
+    - mv results{,_save}/ # move results directory out of the way to preserve it
+    - snakemake --cores 1 --delete-all-output lfhcal_campaigns
+    - mv results{_save,}/