I'm beefing up my benchmark!

Snakefile

@@ -44,3 +44,4 @@ ddsim \
 include: "benchmarks/diffractive_vm/Snakefile"
 include: "benchmarks/dis/Snakefile"
 include: "benchmarks/demp/Snakefile"
+include: "benchmarks/your_benchmark/Snakefile"

benchmarks/your_benchmark/Snakefile

@@ -0,0 +1,71 @@
+import os
+from snakemake.remote.S3 import RemoteProvider as S3RemoteProvider
+from snakemake.remote.HTTP import RemoteProvider as HTTPRemoteProvider
+
+S3 = S3RemoteProvider(
+    endpoint_url="https://dtn01.sdcc.bnl.gov:9000",
+    access_key_id=os.environ["S3_ACCESS_KEY"],
+    secret_access_key=os.environ["S3_SECRET_KEY"],
+)
+
+# Set environment mode (local or eicweb)
+ENV_MODE = os.getenv("ENV_MODE", "local")  # Defaults to "local" if not set
+# Output directory based on environment
+OUTPUT_DIR = "../../sim_output/" if ENV_MODE == "eicweb" else "sim_output/"
+# Benchmark directory based on environment
+BENCH_DIR = "benchmarks/your_benchmark/" if ENV_MODE == "eicweb" else "./"
+
+rule your_benchmark_campaign_reco_get:
+    input:
+        lambda wildcards: S3.remote(f"eictest/EPIC/RECO/24.07.0/epic_craterlake/EXCLUSIVE/UCHANNEL_RHO/10x100/rho_10x100_uChannel_Q2of0to10_hiDiv.{wildcards.INDEX}.eicrecon.tree.edm4eic.root"),
+    output:
+        f"{OUTPUT_DIR}campaign_24.07.0_rho_10x100_uChannel_Q2of0to10_hiDiv_{{INDEX}}_eicrecon.edm4eic.root",
+    shell:
+        """
+echo "Getting file for INDEX {wildcards.INDEX}"
+ln {input} {output}
+"""
+
+rule your_benchmark_analysis:
+    input:
+        script=f"{BENCH_DIR}analysis/uchannelrho.cxx",
+        data=f"{OUTPUT_DIR}campaign_24.07.0_rho_10x100_uChannel_Q2of0to10_hiDiv_{{INDEX}}_eicrecon.edm4eic.root",
+    output:
+        plots=f"{OUTPUT_DIR}campaign_24.07.0_{{INDEX}}_eicrecon.edm4eic/plots.root",
+    shell:
+        """
+mkdir -p $(dirname "{output.plots}")
+root -l -b -q '{input.script}+("{input.data}","{output.plots}")'
+"""
+
+rule your_benchmark_combine:
+    input:
+        lambda wildcards: expand(
+           f"{OUTPUT_DIR}campaign_24.07.0_{{INDEX:04d}}_eicrecon.edm4eic/plots.root",
+           INDEX=range(int(wildcards.N)),
+        ),	
+    wildcard_constraints:
+        N="\d+",
+    output:
+        f"{OUTPUT_DIR}campaign_24.07.0_combined_{{N}}files_eicrecon.edm4eic.plots.root",
+    shell:
+        """
+hadd {output} {input}
+"""
+
+rule your_benchmark_plots:
+    input:
+        script=f"{BENCH_DIR}macros/plot_rho_physics_benchmark.C",
+        plots=f"{OUTPUT_DIR}campaign_24.07.0_combined_{{N}}files_eicrecon.edm4eic.plots.root",
+    output:
+        f"{OUTPUT_DIR}campaign_24.07.0_combined_{{N}}files_eicrecon.edm4eic.plots_figures/benchmark_rho_mass.pdf",
+    shell:
+        """
+if [ ! -d "{input.plots}_figures" ]; then
+    mkdir "{input.plots}_figures"
+    echo "{input.plots}_figures directory created successfully."
+else
+    echo "{input.plots}_figures directory already exists."
+fi
+root -l -b -q '{input.script}("{input.plots}")'
+"""

benchmarks/your_benchmark/analysis/uchannelrho.cxx

@@ -0,0 +1,202 @@
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <string>
+#include <vector>
+#include <algorithm>
+#include <utility>
+
+#include "ROOT/RDataFrame.hxx"
+#include <TH1D.h>
+#include <TFitResult.h>
+#include <TRandom3.h>
+#include <TCanvas.h>
+#include <TSystem.h>
+#include "TFile.h"
+#include "TLorentzVector.h"
+#include "TLorentzRotation.h"
+#include "TVector2.h"
+#include "TVector3.h"
+
+#include "fmt/color.h"
+#include "fmt/core.h"
+
+#include "edm4eic/InclusiveKinematicsData.h"
+#include "edm4eic/ReconstructedParticleData.h"
+#include "edm4hep/MCParticleData.h"
+
+#define PI            3.1415926
+#define MASS_ELECTRON 0.00051
+#define MASS_PROTON   0.93827
+#define MASS_PION     0.13957
+#define MASS_KAON     0.493667
+#define MASS_AU197    183.45406466643374
+
+int uchannelrho(TString rec_file="input.root", TString outputfile="output.root"){	
+	if (gSystem->AccessPathName(rec_file.Data()) != 0) {
+	   // File does not exist
+	   cout<<Form("File %s does not exist.", rec_file.Data())<<endl;
+	   return 0;
+	}
+
+	// read our configuration	
+	auto tree = new TChain("events");
+	TString name_of_input = (TString) rec_file;
+	std::cout << "Input file = " << name_of_input << endl;
+	tree->Add(name_of_input);
+	TTreeReader tree_reader(tree);       // !the tree reader
+
+	//MC-level track attributes
+	TTreeReaderArray<int>   mc_genStatus_array = {tree_reader, "MCParticles.generatorStatus"};
+	TTreeReaderArray<float> mc_px_array = {tree_reader, "MCParticles.momentum.x"};
+	TTreeReaderArray<float> mc_py_array = {tree_reader, "MCParticles.momentum.y"};
+	TTreeReaderArray<float> mc_pz_array = {tree_reader, "MCParticles.momentum.z"};
+	TTreeReaderArray<int>   mc_pdg_array= {tree_reader, "MCParticles.PDG"};
+
+	//Reco-level track attributes
+	TTreeReaderArray<float> reco_px_array = {tree_reader, "ReconstructedChargedParticles.momentum.x"};
+	TTreeReaderArray<float> reco_py_array = {tree_reader, "ReconstructedChargedParticles.momentum.y"};
+	TTreeReaderArray<float> reco_pz_array = {tree_reader, "ReconstructedChargedParticles.momentum.z"};
+	TTreeReaderArray<float> reco_charge_array = {tree_reader, "ReconstructedChargedParticles.charge"};
+	TTreeReaderArray<int>   reco_type     = {tree_reader,"ReconstructedChargedParticles.type"};
+
+	TTreeReaderArray<unsigned int> rec_id = {tree_reader, "ReconstructedChargedParticleAssociations.recID"};
+	TTreeReaderArray<unsigned int> sim_id = {tree_reader, "ReconstructedChargedParticleAssociations.simID"};
+
+	TString output_name_dir = outputfile;
+	cout << "Output file = " << output_name_dir << endl;
+	TFile* output = new TFile(output_name_dir,"RECREATE");
+
+	//Pion reconstruction efficiency histograms
+	TProfile2D* h_effEtaPtPi  = new TProfile2D("h_effEtaPtPi",";#eta; p_{T}(GeV/c)",50,4,6,30,0,1.6);
+	TProfile2D* h_effEtaPtPip = new TProfile2D("h_effEtaPtPip",";#eta; p_{T}(GeV/c)",50,4,6,30,0,1.6);
+	TProfile2D* h_effEtaPtPim = new TProfile2D("h_effEtaPtPim",";#eta; p_{T}(GeV/c)",50,4,6,30,0,1.6);
+	TProfile2D* h_effPhiEtaPi  = new TProfile2D("h_effPhiEtaPi",";#phi (rad);#eta",50,0,6.4,50,4,6);
+	TProfile2D* h_effPhiEtaPip = new TProfile2D("h_effPhiEtaPip",";#phi (rad);#eta",50,0,6.4,50,4,6);
+	TProfile2D* h_effPhiEtaPim = new TProfile2D("h_effPhiEtaPim",";#phi (rad);#eta",50,0,6.4,50,4,6);
+
+	//rho vector meson mass
+	TH1D* h_VM_mass_MC = new TH1D("h_VM_mass_MC",";mass (GeV)",200,0,4);
+	TH1D* h_VM_mass_MC_etacut = new TH1D("h_VM_mass_MC_etacut",";mass (GeV)",200,0,4);
+	TH1D* h_VM_mass_REC = new TH1D("h_VM_mass_REC",";mass (GeV)",200,0,4);
+	TH1D* h_VM_mass_REC_etacut = new TH1D("h_VM_mass_REC_etacut",";mass (GeV)",200,0,4);
+	TH1D* h_VM_mass_REC_justpions = new TH1D("h_VM_mass_REC_justpions",";mass (GeV)",200,0,4);
+
+	cout<<"There are "<<tree->GetEntries()<<" events!!!!!!!"<<endl;
+	tree_reader.SetEntriesRange(0, tree->GetEntries());
+	while (tree_reader.Next()) {
+
+		TLorentzVector vmMC(0,0,0,0);
+		TLorentzVector piplusMC(0,0,0,0);
+		TLorentzVector piminusMC(0,0,0,0);
+
+		//MC level
+		for(unsigned int imc=0;imc<mc_px_array.GetSize();imc++){
+			TVector3 mctrk(mc_px_array[imc],mc_py_array[imc],mc_pz_array[imc]);	
+			if(mc_pdg_array[imc]!=11) mctrk.RotateY(0.025);//Rotate all non-electrons to hadron beam coordinate system
+			if(mc_genStatus_array[imc]!=1) continue;
+			if(mc_pdg_array[imc]==211 && mc_genStatus_array[imc]==1){ 
+			  piplusMC.SetVectM(mctrk,MASS_PION);  
+			}
+		  if(mc_pdg_array[imc]==-211 && mc_genStatus_array[imc]==1){ 
+				piminusMC.SetVectM(mctrk,MASS_PION); 
+			}
+		}
+
+		vmMC=piplusMC+piminusMC;
+		if(vmMC.E()!=0){
+			h_VM_mass_MC->Fill(vmMC.M());
+			if(piplusMC.Theta()>0.009  && piplusMC.Theta()<0.013 && 
+			    piminusMC.Theta()>0.009 && piminusMC.Theta()<0.013 ) h_VM_mass_MC_etacut->Fill(vmMC.M());
+		}
+
+		//4-vector for proton reconstructed as if it were a pi+
+		TLorentzVector protonRECasifpion(0,0,0,0);
+		//pion 4-vectors
+		TLorentzVector piplusREC(0,0,0,0);
+		TLorentzVector piminusREC(0,0,0,0);
+		//rho reconstruction 4-vector
+		TLorentzVector vmREC(0,0,0,0);
+		//rho reconstruction from mis-identified proton 4-vector
+		TLorentzVector vmRECfromproton(0,0,0,0);
+
+		bool isPiMinusFound = false;
+		bool isPiPlusFound = false;
+		bool isProtonFound = false;
+
+		//track loop
+		int numpositivetracks = 0;
+		int failed = 0;
+		for(unsigned int itrk=0;itrk<reco_pz_array.GetSize();itrk++){
+			TVector3 trk(reco_px_array[itrk],reco_py_array[itrk],reco_pz_array[itrk]);
+
+			//  Rotate in order to account for crossing angle 
+			//  and express coordinates in hadron beam pipe frame
+			//  This is just a patch, not a final solution.
+			trk.RotateY(0.025);
+
+			if(reco_type[itrk] == -1){ 
+				failed++;
+				continue;
+			}
+
+
+		  if(reco_charge_array[itrk]>0){ 
+				numpositivetracks++; 
+			  if ((sim_id[itrk - failed]==4 || sim_id[itrk - failed]==5) && reco_charge_array[itrk - failed]==1){
+			    piplusREC.SetVectM(trk,MASS_PION); 
+			    isPiPlusFound=true;
+			  }
+		     if(sim_id[itrk - failed]==6){
+		     	protonRECasifpion.SetVectM(trk,MASS_PION);
+		     	isProtonFound=true; 
+		     }
+			}
+		  if(reco_charge_array[itrk]<0){ 
+		  	piminusREC.SetVectM(trk,MASS_PION); 
+		  	if((sim_id[itrk - failed]==4 || sim_id[itrk - failed]==5) && reco_charge_array[itrk - failed]==-1)	isPiMinusFound=true;
+		  }
+
+		}
+
+		//4vector of VM;
+		if(piplusREC.E()!=0. && piminusREC.E()!=0.){
+			vmREC=piplusREC+piminusREC;
+		}
+		if(protonRECasifpion.E()!=0. && piminusREC.E()!=0.){
+		  vmRECfromproton=protonRECasifpion+piminusREC;
+		}
+
+		//pion reconstruction efficiency
+		double thispipeff = (isPiPlusFound) ? 1.0 : 0.0;
+		double thispimeff = (isPiMinusFound) ? 1.0 : 0.0;
+		h_effEtaPtPi ->Fill(piplusMC.Eta(), piplusMC.Pt(), thispipeff);
+		h_effEtaPtPi ->Fill(piminusMC.Eta(),piminusMC.Pt(),thispimeff);
+		h_effEtaPtPip->Fill(piplusMC.Eta(), piplusMC.Pt(), thispipeff);
+		h_effEtaPtPim->Fill(piminusMC.Eta(),piminusMC.Pt(),thispimeff);
+		//
+		double thispipphi = piplusMC.Phi()>0  ? piplusMC.Phi()  : piplusMC.Phi()+6.2831853;
+		double thispimphi = piminusMC.Phi()>0 ? piminusMC.Phi() : piminusMC.Phi()+6.2831853;
+		if(piplusMC.Pt() >0.2) h_effPhiEtaPi ->Fill(thispipphi, piplusMC.Eta(), thispipeff);
+		if(piminusMC.Pt()>0.2) h_effPhiEtaPi ->Fill(thispimphi,piminusMC.Eta(),thispimeff);
+		if(piplusMC.Pt() >0.2) h_effPhiEtaPip->Fill(thispipphi, piplusMC.Eta(), thispipeff);
+		if(piminusMC.Pt()>0.2) h_effPhiEtaPim->Fill(thispimphi,piminusMC.Eta(),thispimeff);
+		//
+
+		//VM rec
+		if(vmREC.E()==0 && vmRECfromproton.E()==0) continue;
+		double rho_mass = vmREC.M();
+		double rho_mass_fromproton = vmRECfromproton.M();
+		h_VM_mass_REC->Fill(rho_mass);
+		h_VM_mass_REC->Fill(rho_mass_fromproton);
+		if(piplusMC.Theta()>0.009  && piplusMC.Theta()<0.013 &&
+		                  piminusMC.Theta()>0.009 && piminusMC.Theta()<0.013 ) h_VM_mass_REC_etacut->Fill(vmREC.M());
+		if(isPiMinusFound && isPiPlusFound){ 
+		  h_VM_mass_REC_justpions->Fill(rho_mass);
+		}
+	}
+	output->Write();
+	output->Close();
+
+	return 0;
+}

benchmarks/your_benchmark/analyze.sh

@@ -0,0 +1,22 @@
+#!/bin/bash
+source strict-mode.sh
+source benchmarks/your_benchmark/setup.config $*
+
+OUTPUT_PLOTS_DIR=sim_output/nocampaign
+mkdir -p ${OUTPUT_PLOTS_DIR}
+# Analyze
+command time -v \
+root -l -b -q "benchmarks/your_benchmark/analysis/uchannelrho.cxx+(\"${REC_FILE}\",\"${OUTPUT_PLOTS_DIR}/plots.root\")"
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR analysis failed"
+  exit 1
+fi
+
+if [ ! -d "${OUTPUT_PLOTS_DIR}/plots_figures" ]; then
+    mkdir "${OUTPUT_PLOTS_DIR}/plots_figures"
+    echo "${OUTPUT_PLOTS_DIR}/plots_figures directory created successfully."
+else
+    echo "${OUTPUT_PLOTS_DIR}/plots_figures directory already exists."
+fi
+root -l -b -q "benchmarks/your_benchmark/macros/plot_rho_physics_benchmark.C(\"${OUTPUT_PLOTS_DIR}/plots.root\")"
+cat benchmark_output/*.json

benchmarks/your_benchmark/config.yml

@@ -7,11 +7,44 @@ your_benchmark:compile:
 your_benchmark:simulate:
   extends: .phy_benchmark
   stage: simulate
+  needs: ["common:setup"]
+  timeout: 10 hour
   script:
     - echo "I will simulate detector response here!"
+    - config_file=benchmarks/your_benchmark/setup.config
+    - source $config_file
+    - if [ "$USE_SIMULATION_CAMPAIGN" = true ] ; then
+    -     echo "Using simulation campaign so skipping this step!"
+    - else
+    -     echo "Grabbing raw events from S3 and running Geant4"
+    -     bash benchmarks/your_benchmark/simulate.sh
+    -     echo "Geant4 simulations done! Starting eicrecon now!"
+    -     bash benchmarks/your_benchmark/reconstruct.sh
+    - fi
+    - echo "Finished simulating detector response"
+  retry:
+    max: 2
+    when:
+      - runner_system_failure
 
 your_benchmark:results:
   extends: .phy_benchmark
   stage: collect
+  needs:
+    - ["your_benchmark:simulate"]
   script:
     - echo "I will collect results here!"
+    - mkdir -p results/your_benchmark
+    - mkdir -p benchmark_output
+    - config_file=benchmarks/your_benchmark/setup.config
+    - source $config_file
+    - if [ "$USE_SIMULATION_CAMPAIGN" = true ] ; then
+    -     echo "Using simulation campaign!"
+    -     snakemake --cores 2 ../../sim_output/campaign_24.07.0_combined_45files_eicrecon.edm4eic.plots_figures/benchmark_rho_mass.pdf
+    -     cp ../../sim_output/campaign_24.07.0_combined_45files_eicrecon.edm4eic.plots_figures/*.pdf results/your_benchmark/
+    - else
+    -     echo "Not using simulation campaign!"
+    -     bash benchmarks/your_benchmark/analyze.sh
+    -     cp sim_output/nocampaign/plots_figures/*.pdf results/your_benchmark/
+    - fi
+    - echo "Finished copying!"