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

+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

+#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

+#!/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!"

benchmarks/your_benchmark/macros/RiceStyle.h

+#include <iostream>
+#include <vector>
+#include <sstream>
+#include <string>
+
+#include "TGaxis.h"
+#include "TString.h"
+#include "TF1.h"
+#include "TH1.h"
+#include "TH2.h"
+#include "TH3.h"
+#include "TMath.h"
+#include "TTree.h"
+#include "TChain.h"
+#include "TFile.h"
+#include "TCanvas.h"
+#include "TSystem.h"
+#include "TROOT.h"
+#include "TGraph.h"
+#include "TGraphErrors.h"
+#include "TGraphAsymmErrors.h"
+#include "TMultiGraph.h"
+#include "TCanvas.h"
+#include "TPad.h"
+#include "TLegend.h"
+#include "TLatex.h"
+#include "TLine.h"
+#include "TAxis.h"
+#include "TGraph.h"
+#include "TGraphErrors.h"
+#include "TLorentzVector.h"
+#include "TFitResult.h"
+#include "TFitResultPtr.h"
+#include "TMatrixDSym.h"
+#include "TMatrixD.h"
+#include "TArrow.h"
+
+void RiceStyle(){
+
+std::cout << "Welcome to Rice Heavy Ion group!! " << std::endl;
+
+}
+
+//make Canvas
+
+TCanvas* makeCanvas(const char* name, const char *title, bool doLogx = false, bool doLogy = false ){
+	
+	// Start with a canvas
+	TCanvas *canvas = new TCanvas(name,title, 1, 1 ,600 ,600 );
+	// General overall stuff
+	canvas->SetFillColor      (0);
+	canvas->SetBorderMode     (0);
+	canvas->SetBorderSize     (10);
+	// Set margins to reasonable defaults
+	canvas->SetLeftMargin     (0.13);
+	canvas->SetRightMargin    (0.10);
+	canvas->SetTopMargin      (0.10);
+	canvas->SetBottomMargin   (0.13);
+	// Setup a frame which makes sense
+	canvas->SetFrameFillStyle (0);
+	canvas->SetFrameLineStyle (0);
+	canvas->SetFrameBorderMode(0);
+	canvas->SetFrameBorderSize(10);
+	canvas->SetFrameFillStyle (0);
+	canvas->SetFrameLineStyle (0);
+	canvas->SetFrameBorderMode(0);
+	canvas->SetFrameBorderSize(10);
+
+	if( doLogy == true ) gPad->SetLogy(1);
+	if( doLogx == true ) gPad->SetLogx(1);
+	
+	gPad->SetTicks();
+
+	return canvas;
+}
+
+TCanvas* makeMultiCanvas(const char* name, 
+						 const char* title,
+						 int nRows,
+						 int nColumns
+){
+
+	double ratio = nRows/nColumns;
+
+	TCanvas* canvas = new TCanvas( name, title, 1, 1, 400*nRows, 400*nColumns );
+	canvas->SetFillColor      (0);
+	canvas->SetBorderMode     (0);
+	canvas->SetBorderSize     (10);
+	// Set margins to reasonable defaults
+	canvas->SetLeftMargin     (0.30);
+	canvas->SetRightMargin    (0.10);
+	canvas->SetTopMargin      (0.10);
+	canvas->SetBottomMargin   (0.30);
+	// Setup a frame which makes sense
+	canvas->SetFrameFillStyle (0);
+	canvas->SetFrameLineStyle (0);
+	canvas->SetFrameBorderMode(0);
+	canvas->SetFrameBorderSize(10);
+	canvas->SetFrameFillStyle (0);
+	canvas->SetFrameLineStyle (0);
+	canvas->SetFrameBorderMode(0);
+	canvas->SetFrameBorderSize(10);
+ 	
+ 	canvas->Divide(nRows,nColumns,0.01,0.01);
+
+ 	gPad->SetLeftMargin(0.3);
+ 	gPad->SetBottomMargin(0.3);
+	return canvas;
+
+}
+
+void saveCanvas(TCanvas* c, TString dir, TString filename)
+{
+   TDatime* date = new TDatime();
+   //c->Print(Form("../%s/%s_%d.eps",dir.Data(),filename.Data(),date->GetDate()));
+   //c->Print(Form("../%s/%s_%d.gif",dir.Data(),filename.Data(),date->GetDate()));
+   c->Print(Form("../%s/%s_%d.pdf",dir.Data(),filename.Data(),date->GetDate()));
+   //c->Print(Form("../%s/%s_%d.C",dir.Data(),filename.Data(),date->GetDate()));
+   c->Print(Form("../%s/%s_%d.png",dir.Data(),filename.Data(),date->GetDate()));
+}
+
+void initSubPad(TPad* pad, int i)
+{
+  //printf("Pad: %p, index: %d\n",pad,i);
+
+  pad->cd(i);
+  TPad *tmpPad = (TPad*) pad->GetPad(i);
+  tmpPad->SetLeftMargin  (0.20);
+  tmpPad->SetTopMargin   (0.06);
+  tmpPad->SetRightMargin (0.08);
+  tmpPad->SetBottomMargin(0.15);
+  return;
+}
+
+vector<TPad*> makeMultiPad(const int num){//we only have 4,6,8 options for now
+
+	cout << "num: "<< num << endl;
+	vector<TPad*> temp;
+
+	TPad* pad[ num ];
+
+	double setting1[] = {0.0, 0.35, 0.56, 1.0};
+	double setting2[] = {0.0, 0.35, 0.40, 0.7, 1.0 };
+	double setting3[] = {0.0, 0.35, 0.3, 0.533, 0.766, 1.0};
+
+ 	if ( num == 4 ){
+
+		pad[0] = new TPad("pad20", "pad20",setting1[0], setting1[1], setting1[2], setting1[3]);
+		pad[1] = new TPad("pad21", "pad21",setting1[2], setting1[1], setting1[3], setting1[3]);
+		pad[2] = new TPad("pad28", "pad28",setting1[0], setting1[0], setting1[2], setting1[1]);
+		pad[3] = new TPad("pad29", "pad29",setting1[2], setting1[0], setting1[3], setting1[1]);
+
+		for(int i = 0; i < num; i++){
+
+			pad[i]->SetLeftMargin(0.0);
+			pad[i]->SetRightMargin(0);
+			pad[i]->SetTopMargin(0.0);
+			pad[i]->SetBottomMargin(0);
+			pad[i]->Draw();
+			gPad->SetTicks();
+
+		}
+
+		pad[0]->SetLeftMargin(0.265);
+		pad[2]->SetLeftMargin(0.265);
+
+		pad[1]->SetRightMargin(0.05);
+		pad[3]->SetRightMargin(0.05);
+
+		pad[0]->SetTopMargin(0.02);
+		pad[1]->SetTopMargin(0.02);
+
+		pad[2]->SetBottomMargin(0.3);
+		pad[3]->SetBottomMargin(0.3);
+
+	}
+	else if( num == 6 ){
+
+	  	pad[0] = new TPad("pad10", "pad10",setting2[0], setting2[1], setting2[2], setting2[4]);
+	  	pad[1] = new TPad("pad11", "pad11",setting2[2], setting2[1], setting2[3], setting2[4]);
+	  	pad[2] = new TPad("pad12", "pad12",setting2[3], setting2[1], setting2[4], setting2[4]);
+
+	  	pad[3] = new TPad("pad18", "pad18",  setting2[0], setting2[0], setting2[2],  setting2[1]);
+	  	pad[4] = new TPad("pad19", "pad19",  setting2[2], setting2[0], setting2[3],  setting2[1]);
+	  	pad[5] = new TPad("pad110", "pad110",setting2[3], setting2[0], setting2[4],  setting2[1]);
+
+		for(int i = 0; i < num; i++){
+
+			pad[i]->SetLeftMargin(0.0);
+			pad[i]->SetRightMargin(0);
+			pad[i]->SetTopMargin(0.0);
+			pad[i]->SetBottomMargin(0);
+			pad[i]->SetTicks();
+			pad[i]->Draw();
+
+		}
+
+		pad[0]->SetLeftMargin(0.265);
+		pad[3]->SetLeftMargin(0.265);
+
+		pad[2]->SetRightMargin(0.05);
+		pad[5]->SetRightMargin(0.05);
+
+		pad[0]->SetTopMargin(0.02);
+		pad[1]->SetTopMargin(0.02);
+		pad[2]->SetTopMargin(0.02);
+
+		pad[3]->SetBottomMargin(0.30);
+		pad[4]->SetBottomMargin(0.30);
+		pad[5]->SetBottomMargin(0.30);
+
+	}
+	else if( num == 8 ){
+
+		pad[0] = new TPad("pad10", "pad10",setting3[0], setting3[1], setting3[2], setting3[5]);
+		pad[1] = new TPad("pad11", "pad11",setting3[2], setting3[1], setting3[3], setting3[5]);		
+  		pad[2] = new TPad("pad12", "pad12",setting3[3], setting3[1], setting3[4], setting3[5]);
+		pad[3] = new TPad("pad13", "pad13",setting3[4], setting3[1], setting3[5], setting3[5]);
+
+		pad[4] = new TPad("pad18", "pad18",  setting3[0],  setting3[0], setting3[2], setting3[1]);
+		pad[5] = new TPad("pad19", "pad19",  setting3[2],  setting3[0], setting3[3], setting3[1]);
+		pad[6] = new TPad("pad110", "pad110",setting3[3],  setting3[0], setting3[4], setting3[1]);
+		pad[7] = new TPad("pad111", "pad111",setting3[4],  setting3[0], setting3[5], setting3[1]);
+
+		for( int i = 0; i < num; i++ ){
+
+			pad[i]->SetLeftMargin(0.0);
+			pad[i]->SetRightMargin(0);
+			pad[i]->SetTopMargin(0.0);
+			pad[i]->SetBottomMargin(0);
+			pad[i]->SetTicks();
+			pad[i]->Draw();		
+		}
+
+		pad[0]->SetLeftMargin(0.265);
+		pad[4]->SetLeftMargin(0.265);
+
+		pad[3]->SetRightMargin(0.05);
+		pad[7]->SetRightMargin(0.05);
+
+		pad[0]->SetTopMargin(0.05);
+		pad[1]->SetTopMargin(0.05);
+		pad[2]->SetTopMargin(0.05);
+		pad[3]->SetTopMargin(0.05);
+
+		pad[4]->SetBottomMargin(0.30);
+		pad[5]->SetBottomMargin(0.30);
+		pad[6]->SetBottomMargin(0.30);
+		pad[7]->SetBottomMargin(0.30);
+	}
+
+	for( int i = 0; i < num; i++){
+
+		temp.push_back( pad[i] );
+	}
+
+	return temp;
+}
+
+TH1D* makeHist(const char*name, const char*title, const char*xtit, const char*ytit, const int nBins, const double lower, const double higher, EColor color = kBlack ){
+
+	TH1D* temp = new TH1D(name, title, nBins, lower, higher);
+	
+	temp->SetMarkerSize(1.0);
+	temp->SetMarkerStyle(20);
+	temp->SetMarkerColor(color);
+	temp->SetLineColor(color);
+	temp->SetStats(kFALSE);
+
+	temp->GetXaxis()->SetTitle( xtit );
+	temp->GetXaxis()->SetTitleSize(0.05);
+	temp->GetXaxis()->SetTitleFont(42);
+	temp->GetXaxis()->SetTitleOffset(1.25);
+	temp->GetXaxis()->SetLabelSize(0.05);
+	temp->GetXaxis()->SetLabelOffset(0.01);
+	temp->GetXaxis()->SetLabelFont(42);
+	temp->GetXaxis()->SetLabelColor(kBlack);
+	temp->GetXaxis()->CenterTitle();
+
+	temp->GetYaxis()->SetTitle( ytit );
+	temp->GetYaxis()->SetTitleSize(0.05);
+	temp->GetYaxis()->SetTitleFont(42);
+	temp->GetYaxis()->SetTitleOffset(1.4);
+	temp->GetYaxis()->SetLabelSize(0.05);
+	temp->GetYaxis()->SetLabelOffset(0.01);
+	temp->GetYaxis()->SetLabelFont(42);
+	temp->GetYaxis()->SetLabelColor(kBlack);
+	temp->GetYaxis()->CenterTitle();
+
+	return temp;
+}
+
+TH1D* makeHistDifferentBins(const char*name, const char*title, const char*xtit, const char*ytit, const int nBins, double bins[], EColor color = kBlack ){
+
+	TH1D* temp = new TH1D(name, title, nBins, bins);
+	
+	temp->SetMarkerSize(1.0);
+	temp->SetMarkerStyle(20);
+	temp->SetMarkerColor(color);
+	temp->SetLineColor(color);
+	temp->SetStats(kFALSE);
+
+	temp->GetXaxis()->SetTitle( xtit );
+	temp->GetXaxis()->SetTitleSize(0.05);
+	temp->GetXaxis()->SetTitleFont(42);
+	temp->GetXaxis()->SetTitleOffset(1.25);
+	temp->GetXaxis()->SetLabelSize(0.05);
+	temp->GetXaxis()->SetLabelOffset(0.01);
+	temp->GetXaxis()->SetLabelFont(42);
+	temp->GetXaxis()->SetLabelColor(kBlack);
+	temp->GetXaxis()->CenterTitle();
+
+	temp->GetYaxis()->SetTitle( ytit );
+	temp->GetYaxis()->SetTitleSize(0.05);
+	temp->GetYaxis()->SetTitleFont(42);
+	temp->GetYaxis()->SetTitleOffset(1.4);
+	temp->GetYaxis()->SetLabelSize(0.05);
+	temp->GetYaxis()->SetLabelOffset(0.01);
+	temp->GetYaxis()->SetLabelFont(42);
+	temp->GetYaxis()->SetLabelColor(kBlack);
+	temp->GetYaxis()->CenterTitle();
+
+	return temp;
+}
+
+void fixedFontHist1D(TH1 * h, Float_t xoffset=1.5, Float_t yoffset=2.3)
+{
+  h->SetLabelFont(43,"X");
+  h->SetLabelFont(43,"Y");
+  //h->SetLabelOffset(0.01);
+  h->SetLabelSize(16);
+  h->SetTitleFont(43);
+  h->SetTitleSize(20);
+  h->SetLabelSize(15,"Y");
+  h->SetTitleFont(43,"Y");
+  h->SetTitleSize(20,"Y");
+  h->SetTitleOffset(xoffset,"X");
+  h->SetTitleOffset(yoffset,"Y");
+  h->GetXaxis()->CenterTitle();
+  h->GetYaxis()->CenterTitle();
+  
+}
+
+TH2D* make2DHist( const char*name, 
+			      const char*title, 
+			      const char*xtit, 
+			      const char*ytit,
+			      const int nxbins,
+			      const double xlow, 
+			      const double xhigh,
+			      const int nybins,
+			      const double ylow, 
+			      const double yhigh
+){
+
+	TH2D* temp2D = new TH2D(name, title, nxbins, xlow, xhigh, nybins, ylow, yhigh);
+
+	temp2D->SetMarkerSize(1.0);
+	temp2D->SetMarkerStyle(20);
+	temp2D->SetMarkerColor(kBlack);
+	temp2D->SetLineColor(kBlack);
+	temp2D->SetStats(kFALSE);
+
+	temp2D->GetXaxis()->SetTitle( xtit );
+	temp2D->GetXaxis()->SetTitleSize(0.04);
+	temp2D->GetXaxis()->SetTitleFont(42);
+	temp2D->GetXaxis()->SetTitleOffset(1.4);
+	temp2D->GetXaxis()->SetLabelSize(0.04);
+	temp2D->GetXaxis()->SetLabelOffset(0.01);
+	temp2D->GetXaxis()->SetLabelFont(42);
+	temp2D->GetXaxis()->SetLabelColor(kBlack);
+
+	temp2D->GetYaxis()->SetTitle( ytit );
+	temp2D->GetYaxis()->SetTitleSize(0.04);
+	temp2D->GetYaxis()->SetTitleFont(42);
+	temp2D->GetYaxis()->SetTitleOffset(1.7);
+	temp2D->GetYaxis()->SetLabelSize(0.04);
+	temp2D->GetYaxis()->SetLabelOffset(0.01);
+	temp2D->GetYaxis()->SetLabelFont(42);
+	temp2D->GetYaxis()->SetLabelColor(kBlack);
+
+	return temp2D;
+
+}
+
+void fixedFontHist(TH2D * h, Float_t xoffset=0.9, Float_t yoffset=2.7)
+{
+  h->SetLabelFont(43,"X");
+  h->SetLabelFont(43,"Y");
+  //h->SetLabelOffset(0.01);
+  h->SetLabelSize(16);
+  h->SetTitleFont(43);
+  h->SetTitleSize(20);
+  h->SetLabelSize(15,"Y");
+  h->SetTitleFont(43,"Y");
+  h->SetTitleSize(17,"Y");
+  h->SetTitleOffset(xoffset,"X");
+  h->SetTitleOffset(yoffset,"Y");
+  h->GetXaxis()->CenterTitle();
+  h->GetYaxis()->CenterTitle();
+}
+void make_dNdX( TH1D* hist ){
+
+	for(int i=0;i<hist->GetNbinsX();i++){
+		double value = hist->GetBinContent(i+1);
+		double error = hist->GetBinError(i+1);
+		double binwidth = hist->GetBinWidth(i+1);
+
+		hist->SetBinContent(i+1, value / binwidth );
+		hist->SetBinError(i+1, error / binwidth );
+	}
+}
+
+double calColError(double Ea, double Eb, double Sa, double Sb){
+
+	double temp = Ea/Eb;
+	double temp2 = (Sa*Sa)/(Ea*Ea) + (Sb*Sb)/(Eb*Eb);
+	double temp3 = (2.*Sa*Sa)/(Ea*Eb);
+
+	return temp*(sqrt(TMath::Abs(temp2-temp3)) );
+}
+
+TH1D* make_systematicRatio(TH1D* hist1, TH1D* hist2){
+
+	TH1D* hist_ratio = (TH1D*) hist1->Clone("hist_ratio");
+
+	if( hist1->GetNbinsX() != hist2->GetNbinsX() ){
+		std::cout << "Not compatible binning, abort!" << std::endl;
+		return 0;
+	}
+
+	for(int ibin=0;ibin<hist1->GetNbinsX();ibin++){
+		double value_a = hist1->GetBinContent(ibin+1);
+		double error_a = hist1->GetBinError(ibin+1);
+
+		double value_b = hist2->GetBinContent(ibin+1);
+		double error_b = hist2->GetBinError(ibin+1);
+
+		hist_ratio->SetBinContent(ibin+1, value_a / value_b );
+		hist_ratio->SetBinError(ibin+1, calColError(value_a, value_b, error_a, error_b) );
+	}
+
+	return hist_ratio;
+}
+
+TLegend* makeLegend(){
+
+	TLegend *w2 = new TLegend(0.65,0.15,0.90,0.45);
+	w2->SetLineColor(kWhite);
+	w2->SetFillColor(0);
+	return w2;
+
+}
+
+TGraphAsymmErrors* makeEfficiency(TH1D* hist1, TH1D* hist2, const char*Draw = "cp", EColor color = kBlack  ){
+
+	TGraphAsymmErrors* temp = new TGraphAsymmErrors();
+	temp->Divide( hist1, hist2, Draw );
+	temp->SetMarkerStyle(20);
+	temp->SetMarkerColor(color);
+	temp->SetLineColor(color);
+
+	return temp;
+
+}
+
+TLatex* makeLatex(const char* txt,  double x, double y){
+
+	TLatex* r = new TLatex(x, y, txt);
+	r->SetTextSize(0.05);
+	r->SetNDC();
+	return r;
+
+}
+
+void drawBox(TH1D* hist1, double sys, bool doPercentage, double xe= 0.05){
+
+	TBox* temp_box[100];
+	int bins = hist1->GetNbinsX();
+
+	for(int deta = 0; deta < bins; deta++){
+
+		double value = hist1->GetBinContent(deta+1);
+		double bincenter = hist1->GetBinCenter(deta+1);
+		double sys_temp = 0.;
+
+		if( doPercentage ) sys_temp = value*sys;
+		else sys_temp = sys;
+
+		temp_box[deta] = new TBox(bincenter-xe,value-sys_temp,bincenter+xe,value+sys_temp);
+	    temp_box[deta]->SetFillColorAlpha(kGray+2,0.4);
+	    temp_box[deta]->SetFillStyle(1001);
+		temp_box[deta]->SetLineWidth(0);
+    	temp_box[deta]->Draw("same");
+    }
+
+}
+
+void drawBoxRatio(TH1D* hist1, TH1D* hist2, double sys, bool doPercentage){
+
+	TBox* temp_box[100];
+	double xe = 0.08;
+
+	int bins = hist1->GetNbinsX();
+
+	for(int deta = 0; deta < bins; deta++){
+
+		if(deta > 6) continue;
+
+		double factor = hist2->GetBinContent(deta+1); 
+		double value = hist1->GetBinContent(deta+1);
+		double bincenter = hist1->GetBinCenter(deta+1);
+
+		if( doPercentage ) sys = sqrt((value*0.045)*(value*0.045));
+		else sys = sys;
+
+		double sys_temp;
+		sys_temp = sys;
+	
+		temp_box[deta] = new TBox(bincenter-xe,value-sys_temp,bincenter+xe,value+sys_temp);
+	    temp_box[deta]->SetFillColorAlpha(kGray+2,0.4);
+	    temp_box[deta]->SetFillStyle(1001);
+		temp_box[deta]->SetLineWidth(0);
+    	temp_box[deta]->Draw("same");
+
+    }
+
+}
+
+void drawBoxTGraphRatio(TGraphErrors* gr1, int bins, double sys, bool doPercentage){
+
+	double xe[11];
+	TBox* box1[11];
+
+    for(int mult = 0; mult < bins; mult++){
+
+    	if( mult < 6 ){
+
+    		xe[mult] = 15*log(1.1*(mult+1));
+    		if(mult == 0) xe[mult] = 10;
+    	}
+  		if( mult ==  6) xe[mult] = 37;
+  		if( mult ==  7) xe[mult] = 50;
+  		if( mult ==  8) xe[mult] = 50;
+  		if( mult ==  9) xe[mult] = 63;
+  		if( mult ==  10) xe[mult] = 73;
+    	
+    	
+    	double x1;
+    	double value1;
+    	gr1->GetPoint(mult, x1, value1);
+
+    	double ye = sys;
+  		if( doPercentage ) ye = sqrt((value1*0.045)*(value1*0.045));
+  		else ye = sys;
+
+    	box1[mult] = new TBox(x1-xe[mult],value1-ye,x1+xe[mult],value1+ye);
+	 	box1[mult]->SetFillColorAlpha(kGray+2,0.4);
+        box1[mult]->SetFillStyle(1001);
+    	box1[mult]->SetLineWidth(0);
+    	box1[mult]->SetLineColor(kRed);
+        box1[mult]->Draw("SAME");
+
+    }
+
+}
+
+void drawBoxTGraph(TGraphErrors* gr1, int bins, double sys, bool doPercentage, bool doConstantWidth){
+
+	double xe[100];
+	TBox* box1[100];
+
+    for(int mult = 0; mult < bins; mult++){
+
+    	if(!doConstantWidth){
+	    	if( mult < 6 ){
+
+	    		xe[mult] = 15*log(1.1*(mult+1));
+	    		if(mult == 0) xe[mult] = 10;
+	    	}
+	  		if( mult ==  6) xe[mult] = 37;
+	  		if( mult ==  7) xe[mult] = 50;
+	  		if( mult ==  8) xe[mult] = 50;
+	  		if( mult ==  9) xe[mult] = 63;
+	  		if( mult ==  10) xe[mult] = 73;
+  		}
+  		else{ xe[mult] = 0.02;}
+    	
+    	
+    	double x1;
+    	double value1;
+    	gr1->GetPoint(mult, x1, value1);
+
+    	double ye = sys;
+  		if( doPercentage ) ye = value1 * sys;
+  		else ye = sys;
+
+    	box1[mult] = new TBox(x1-xe[mult],value1-ye,x1+xe[mult],value1+ye);
+	 	box1[mult]->SetFillColorAlpha(kGray+2,0.4);
+        box1[mult]->SetFillStyle(1001);
+    	box1[mult]->SetLineWidth(0);
+    	box1[mult]->SetLineColor(kBlack);
+        box1[mult]->Draw("LSAME");
+
+    }
+}
+
+void drawBoxTGraph_alt(TGraphErrors* gr1, int bins, double sys, bool doPercentage, bool doConstantWidth){
+
+	double xe[100];
+	TBox* box1[100];
+
+    for(int mult = 0; mult < bins; mult++){
+
+    	if(!doConstantWidth){
+	    	if( mult < 6 ){
+
+	    		xe[mult] = 15*log(1.1*(mult+1));
+	    		if(mult == 0) xe[mult] = 10;
+	    	}
+	  		if( mult ==  6) xe[mult] = 37;
+	  		if( mult ==  7) xe[mult] = 50;
+	  		if( mult ==  8) xe[mult] = 50;
+	  		if( mult ==  9) xe[mult] = 63;
+	  		if( mult ==  10) xe[mult] = 73;
+  		}
+  		else{ xe[mult] = 0.005;}
+    	
+    	
+    	double x1;
+    	double value1;
+    	gr1->GetPoint(mult, x1, value1);
+
+    	double ye = sys;
+  		if( doPercentage ) ye = value1 * sys;
+  		else ye = sys;
+
+    	box1[mult] = new TBox(x1-xe[mult],value1-ye,x1+xe[mult],value1+ye);
+	 	box1[mult]->SetFillColorAlpha(kGray+2,0.4);
+        box1[mult]->SetFillStyle(1001);
+    	box1[mult]->SetLineWidth(0);
+    	box1[mult]->SetLineColor(kRed);
+        box1[mult]->Draw("SAME");
+
+    }
+}
+
+
+void drawBoxTGraphDiff(TGraphErrors* gr1, TGraphErrors* gr2, int bins, double sys, bool doPercentage){
+
+	double xe[11];
+	TBox* box1[11];
+	TBox* box2[11];
+
+    for(int mult = 0; mult < bins; mult++){
+
+    	xe[mult] = 10*log(1.1*(mult+1));
+    	if(mult == 0) xe[mult] = 7;
+
+    	double x1;
+    	double value1;
+    	gr1->GetPoint(mult, x1, value1);
+
+    	double x2;
+    	double value2;
+    	gr2->GetPoint(mult, x2, value2);
+
+    	double value = value2 - value1;
+
+    	double ye = sys;
+  		if( doPercentage ) ye = sqrt((value*0.045)*(value*0.045)+sys*sys);
+  		else ye = sys;
+
+    	box1[mult] = new TBox(x1-xe[mult],value-ye,x1+xe[mult],value+ye);
+	 	box1[mult]->SetFillColorAlpha(kGray+2,0.4);
+        box1[mult]->SetFillStyle(1001);
+    	box1[mult]->SetLineWidth(0);
+    	box1[mult]->SetLineColor(kRed);
+        box1[mult]->Draw("SAME");
+
+    }
+
+}
+	

benchmarks/your_benchmark/macros/plot_rho_physics_benchmark.C

+#include "RiceStyle.h"
+
+using namespace std;
+
+void plot_rho_physics_benchmark(TString filename="./sim_output/plot_combined.root"){
+  Ssiz_t dotPosition = filename.Last('.');
+  TString figure_directory = filename(0, dotPosition);
+  figure_directory += "_figures";	
+  
+  TFile* file = new TFile(filename);
+  TString vm_label="#rho^{0}";
+  TString daug_label="#pi^{#plus}#pi^{#minus}";
+  //mass distribution
+  TH1D* h_VM_mass_MC = (TH1D*) file->Get("h_VM_mass_MC");
+  TH1D* h_VM_mass_REC = (TH1D*) file->Get("h_VM_mass_REC");
+  TH1D* h_VM_mass_REC_justpions = (TH1D*) file->Get("h_VM_mass_REC_justpions");
+  //mass distribution within B0
+  TH1D* h_VM_mass_MC_etacut = (TH1D*) file->Get("h_VM_mass_MC_etacut");
+  TH1D* h_VM_mass_REC_etacut = (TH1D*) file->Get("h_VM_mass_REC_etacut");
+  //efficiencies
+  TProfile2D* h_effEtaPtPi = (TProfile2D*) file->Get("h_effEtaPtPi");
+  TProfile2D* h_effEtaPtPip = (TProfile2D*) file->Get("h_effEtaPtPip"); 
+  TProfile2D* h_effEtaPtPim = (TProfile2D*) file->Get("h_effEtaPtPim"); 
+  TProfile2D* h_effPhiEtaPi = (TProfile2D*) file->Get("h_effPhiEtaPi"); 
+  TProfile2D* h_effPhiEtaPip = (TProfile2D*) file->Get("h_effPhiEtaPip"); 
+  TProfile2D* h_effPhiEtaPim = (TProfile2D*) file->Get("h_effPhiEtaPim"); 	
+
+
+  TLatex* r42 = new TLatex(0.18, 0.91, "ep 10#times100 GeV");
+  r42->SetNDC();
+  r42->SetTextSize(22);
+  r42->SetTextFont(43);
+  r42->SetTextColor(kBlack);
+
+  TLatex* r43 = new TLatex(0.9,0.91, "#bf{EPIC}");
+  r43->SetNDC();
+  //r43->SetTextSize(0.04);
+  r43->SetTextSize(22);
+
+  TLatex* r44 = new TLatex(0.53, 0.78, "10^{-3}<Q^{2}<10 GeV^{2}, W>2 GeV");
+  r44->SetNDC();
+  r44->SetTextSize(20);
+  r44->SetTextFont(43);
+  r44->SetTextColor(kBlack);
+
+  TLatex* r44_2 = new TLatex(0.5, 0.83, ""+vm_label+" #rightarrow "+daug_label+" eSTARlight");
+  r44_2->SetNDC();
+  r44_2->SetTextSize(30);
+  r44_2->SetTextFont(43);
+  r44_2->SetTextColor(kBlack);
+
+  TCanvas* c2 = new TCanvas("c2","c2",1,1,600,600);
+  gPad->SetTicks();
+  gPad->SetLeftMargin(0.18);
+  gPad->SetBottomMargin(0.18);
+  gPad->SetTopMargin(0.10);
+  gPad->SetRightMargin(0.01);
+  TH1D* base2 = makeHist("base2", "", "#pi^{#plus}#pi^{#minus} inv. mass (GeV)", "counts", 100,0.05,2.05,kBlack);
+  base2->GetYaxis()->SetRangeUser(0.5, 1.2*(h_VM_mass_MC->GetMaximum()));
+  base2->GetXaxis()->SetTitleColor(kBlack);
+  fixedFontHist1D(base2,1.,1.2);
+  base2->GetYaxis()->SetTitleSize(base2->GetYaxis()->GetTitleSize()*1.5);
+  base2->GetXaxis()->SetTitleSize(base2->GetXaxis()->GetTitleSize()*1.5);
+  base2->GetYaxis()->SetLabelSize(base2->GetYaxis()->GetLabelSize()*1.5);
+  base2->GetXaxis()->SetLabelSize(base2->GetXaxis()->GetLabelSize()*1.5);
+  base2->GetXaxis()->SetNdivisions(4,4,0);
+  base2->GetYaxis()->SetNdivisions(5,5,0);
+  base2->GetYaxis()->SetTitleOffset(1.3);
+  base2->Draw();
+
+  TH1D* h_VM_mass_REC_justprotons = (TH1D*)h_VM_mass_REC->Clone("h_VM_mass_REC_justprotons");
+  for(int ibin=1; ibin<h_VM_mass_REC_justprotons->GetNbinsX(); ibin++){
+    h_VM_mass_REC_justprotons->SetBinContent(ibin,h_VM_mass_REC_justprotons->GetBinContent(ibin) - h_VM_mass_REC_justpions->GetBinContent(ibin));
+  }
+
+  h_VM_mass_MC->SetFillColorAlpha(kBlack,0.2);
+  h_VM_mass_REC->SetFillColorAlpha(kMagenta,0.2);
+  h_VM_mass_MC->SetLineColor(kBlack);
+  h_VM_mass_REC->SetLineColor(kMagenta);
+  h_VM_mass_REC_justpions->SetLineColor(kViolet+10);
+  h_VM_mass_REC_justprotons->SetLineColor(kRed);
+  h_VM_mass_MC->SetLineWidth(2);
+  h_VM_mass_REC->SetLineWidth(2);
+  h_VM_mass_REC_justpions->SetLineWidth(2);
+  h_VM_mass_REC_justprotons->SetLineWidth(2);
+
+  h_VM_mass_REC->Scale(3.0);
+  h_VM_mass_REC_justpions->Scale(3.0);
+  h_VM_mass_REC_justprotons->Scale(3.0);
+
+  h_VM_mass_MC->Draw("HIST E same");
+  h_VM_mass_REC->Draw("HIST E same");
+  h_VM_mass_REC_justpions->Draw("HIST same");
+  h_VM_mass_REC_justprotons->Draw("HIST same");
+
+  r42->Draw("same");
+  r43->Draw("same");
+  r44->Draw("same");
+  r44_2->Draw("same");
+
+  TLegend *w8 = new TLegend(0.58,0.63,0.93,0.76);
+  w8->SetLineColor(kWhite);
+  w8->SetFillColor(0);
+  w8->SetTextSize(17);
+  w8->SetTextFont(45);
+  w8->AddEntry(h_VM_mass_MC, ""+vm_label+" MC ", "L");
+  w8->AddEntry(h_VM_mass_REC, vm_label+" reco.#times3", "L");
+  w8->AddEntry(h_VM_mass_REC_justpions, vm_label+" reco.#times3 (#pi^{#minus}#pi^{#plus})", "L");
+  w8->AddEntry(h_VM_mass_REC_justprotons, vm_label+" reco.#times3 (#pi^{#minus}p)", "L");
+  w8->Draw("same");
+
+  TString figure1name = figure_directory+"/benchmark_rho_mass.pdf";
+  c2->Print(figure1name);
+
+
+  ///////////////////// Figure 4
+  TCanvas* c4 = new TCanvas("c4","c4",1,1,600,600);
+  gPad->SetTicks();
+  gPad->SetLeftMargin(0.18);
+  gPad->SetBottomMargin(0.18);
+  gPad->SetTopMargin(0.10);
+  gPad->SetRightMargin(0.01);
+  TH1D* base4 = makeHist("base4", "", "#pi^{#plus}#pi^{#minus} inv. mass (GeV)", "counts", 100,0.05,2.05,kBlack);
+  base4->GetYaxis()->SetRangeUser(0.5, 1.2*(h_VM_mass_MC_etacut->GetMaximum()));
+  base4->GetXaxis()->SetTitleColor(kBlack);
+  fixedFontHist1D(base4,1.,1.2);
+  base4->GetYaxis()->SetTitleSize(base4->GetYaxis()->GetTitleSize()*1.5);
+  base4->GetXaxis()->SetTitleSize(base4->GetXaxis()->GetTitleSize()*1.5);
+  base4->GetYaxis()->SetLabelSize(base4->GetYaxis()->GetLabelSize()*1.5);
+  base4->GetXaxis()->SetLabelSize(base4->GetXaxis()->GetLabelSize()*1.5);
+  base4->GetXaxis()->SetNdivisions(4,4,0);
+  base4->GetYaxis()->SetNdivisions(5,5,0);
+  base4->GetYaxis()->SetTitleOffset(1.3);
+  base4->Draw();
+
+  h_VM_mass_MC_etacut->SetFillColorAlpha(kBlack,0.2);
+  h_VM_mass_REC_etacut->SetFillColorAlpha(kMagenta,0.2);
+  h_VM_mass_MC_etacut->SetLineColor(kBlack);
+  h_VM_mass_REC_etacut->SetLineColor(kMagenta);
+  h_VM_mass_MC_etacut->SetLineWidth(2);
+  h_VM_mass_REC_etacut->SetLineWidth(2);
+
+  h_VM_mass_MC_etacut->Draw("HIST E same");
+  h_VM_mass_REC_etacut->Draw("HIST E same");
+
+  double minbineff = h_VM_mass_MC_etacut->FindBin(0.6);
+  double maxbineff = h_VM_mass_MC_etacut->FindBin(1.0);
+  double thiseff = 100.0*(1.0*h_VM_mass_REC_etacut->Integral(minbineff,maxbineff))/(1.0*h_VM_mass_MC_etacut->Integral(minbineff,maxbineff));
+
+  r42->Draw("same");
+  r43->Draw("same");
+  r44->Draw("same");
+  r44_2->Draw("same");
+
+  TLegend *w10 = new TLegend(0.58,0.62,0.93,0.7);
+  w10->SetLineColor(kWhite);
+  w10->SetFillColor(0);
+  w10->SetTextSize(17);
+  w10->SetTextFont(45);
+  w10->AddEntry(h_VM_mass_MC_etacut, vm_label+" MC ", "L");
+  w10->AddEntry(h_VM_mass_REC_etacut, vm_label+" reco. (#pi^{#minus}#pi^{#plus})", "L");
+  w10->Draw("same");
+
+  TLatex* anglelabel = new TLatex(0.59, 0.73, "9<#theta_{#pi^{#pm},MC}<13 mrad");
+  anglelabel->SetNDC();
+  anglelabel->SetTextSize(20);
+  anglelabel->SetTextFont(43);
+  anglelabel->SetTextColor(kBlack);
+  anglelabel->Draw("same");
+
+  TLatex* efflabel = new TLatex(0.59, 0.55, "reco. eff (0.6<m^{2}<1 GeV^{2})");
+  efflabel->SetNDC();
+  efflabel->SetTextSize(20);
+  efflabel->SetTextFont(43);
+  efflabel->SetTextColor(kBlack);
+  efflabel->Draw("same");
+  TLatex* effnlabel = new TLatex(0.59, 0.51, Form("          = %.0f%%",thiseff));
+  effnlabel->SetNDC();
+  effnlabel->SetTextSize(20);
+  effnlabel->SetTextFont(43);
+  effnlabel->SetTextColor(kBlack);
+  effnlabel->Draw("same");
+
+  TString figure2name = figure_directory+"/benchmark_rho_mass_cuts.pdf";
+  c4->Print(figure2name);
+
+  ///////////////////// Figure 5
+  TCanvas* c5 = new TCanvas("c5","c5",1,1,700,560);
+  TPad* p5 = new TPad("p5","Pad5",0.,0.,1.,1.);
+  p5->Divide(3,2,0,0);
+  p5->Draw();
+  gStyle->SetPalette(kBlueRedYellow);
+  gStyle->SetOptStat(0);	
+
+  h_effEtaPtPi  ->GetXaxis()->SetLabelSize(h_effEtaPtPi  ->GetXaxis()->GetLabelSize()*1.8);
+  h_effEtaPtPip ->GetXaxis()->SetLabelSize(h_effEtaPtPip ->GetXaxis()->GetLabelSize()*1.8);
+  h_effEtaPtPim ->GetXaxis()->SetLabelSize(h_effEtaPtPim ->GetXaxis()->GetLabelSize()*1.8);
+  h_effEtaPtPi  ->GetYaxis()->SetLabelSize(h_effEtaPtPi  ->GetYaxis()->GetLabelSize()*1.8);
+  h_effEtaPtPim ->GetZaxis()->SetLabelSize(h_effEtaPtPim ->GetZaxis()->GetLabelSize()*0.5);
+  h_effEtaPtPim ->GetZaxis()->SetTitleSize(h_effEtaPtPim ->GetZaxis()->GetTitleSize()*0.5);
+  h_effPhiEtaPi ->GetXaxis()->SetLabelSize(h_effPhiEtaPi ->GetXaxis()->GetLabelSize()*1.8);
+  h_effPhiEtaPip->GetXaxis()->SetLabelSize(h_effPhiEtaPip->GetXaxis()->GetLabelSize()*1.8);
+  h_effPhiEtaPim->GetXaxis()->SetLabelSize(h_effPhiEtaPim->GetXaxis()->GetLabelSize()*1.8);
+  h_effPhiEtaPi ->GetYaxis()->SetLabelSize(h_effPhiEtaPi ->GetYaxis()->GetLabelSize()*1.8);
+  h_effPhiEtaPim->GetZaxis()->SetLabelSize(h_effPhiEtaPim->GetZaxis()->GetLabelSize()*0.5);
+  h_effPhiEtaPim->GetZaxis()->SetTitleSize(h_effPhiEtaPim->GetZaxis()->GetTitleSize()*0.5);
+
+  fixedFontHist1D(h_effEtaPtPi,1.,1.2);
+  fixedFontHist1D(h_effEtaPtPip,1.,1.2);
+  fixedFontHist1D(h_effEtaPtPim,1.,1.2);
+  fixedFontHist1D(h_effPhiEtaPi,1.,1.2);
+  fixedFontHist1D(h_effPhiEtaPip,1.,1.2);
+  fixedFontHist1D(h_effPhiEtaPim,1.,1.2);
+
+  p5->cd(1);
+  TVirtualPad* p51 = p5->cd(1);
+  p51->SetTopMargin(0.08);
+  p51->SetRightMargin(0);
+  p51->SetLeftMargin(0.21);
+  p51->SetBottomMargin(0.2);
+  h_effEtaPtPi->GetXaxis()->SetRangeUser(3.9,6.05);
+  h_effEtaPtPi->GetYaxis()->SetRangeUser(0,1.7);
+  h_effEtaPtPi->GetZaxis()->SetRangeUser(0,1);
+  h_effEtaPtPi->GetXaxis()->SetNdivisions(5);
+  h_effEtaPtPi->GetYaxis()->SetNdivisions(5);
+  h_effEtaPtPi->SetContour(99);
+  h_effEtaPtPi->Draw("COLZ");
+  TLatex* pilabel = new TLatex(0.81, 0.75, "#pi^{#pm}");
+  pilabel->SetNDC();
+  pilabel->SetTextSize(40);
+  pilabel->SetTextFont(43);
+  pilabel->SetTextColor(kBlack);
+  pilabel->Draw("same");
+  TLatex* r44fig5c = new TLatex(0.21, 0.93, "ep 10#times100 GeV         #rho^{0}#rightarrow#pi^{#plus}#pi^{#minus}");
+  r44fig5c->SetNDC();
+  r44fig5c->SetTextSize(15);
+  r44fig5c->SetTextFont(43);
+  r44fig5c->SetTextColor(kBlack);
+  r44fig5c->Draw("same");
+
+  p5->cd(2);
+  TVirtualPad* p52 = p5->cd(2);
+  p52->SetTopMargin(0.08);
+  p52->SetRightMargin(0);
+  p52->SetLeftMargin(0);
+  p52->SetBottomMargin(0.2);
+  h_effEtaPtPip->GetXaxis()->SetRangeUser(4.05,6.05);
+  h_effEtaPtPip->GetYaxis()->SetRangeUser(0,1.7);
+  h_effEtaPtPip->GetZaxis()->SetRangeUser(0,1);
+  h_effEtaPtPip->GetXaxis()->SetNdivisions(5);
+  h_effEtaPtPip->GetYaxis()->SetNdivisions(5);
+  h_effEtaPtPip->SetContour(99);
+  h_effEtaPtPip->Draw("COLZ");
+  TLatex* piplabel = new TLatex(0.81, 0.75, "#pi^{#plus}");
+  piplabel->SetNDC();
+  piplabel->SetTextSize(40);
+  piplabel->SetTextFont(43);
+  piplabel->SetTextColor(kBlack);
+  piplabel->Draw("same");
+  TLatex* r44fig5a = new TLatex(0.01, 0.93, "eSTARlight        10^{-3}<Q^{2}<10 GeV^{2}");
+  r44fig5a->SetNDC();
+  r44fig5a->SetTextSize(15);
+  r44fig5a->SetTextFont(43);
+  r44fig5a->SetTextColor(kBlack);
+  r44fig5a->Draw("same");
+
+  p5->cd(3);
+  TVirtualPad* p53 = p5->cd(3);
+  p53->SetTopMargin(0.08);
+  p53->SetRightMargin(0.2);
+  p53->SetLeftMargin(0);
+  p53->SetBottomMargin(0.2);
+  h_effEtaPtPim->SetTitle(";#eta;;efficiency");
+  h_effEtaPtPim->GetXaxis()->SetRangeUser(4.05,6.05);
+  h_effEtaPtPim->GetYaxis()->SetRangeUser(0,1.7);
+  h_effEtaPtPim->GetZaxis()->SetRangeUser(0,1);
+  h_effEtaPtPim->GetXaxis()->SetNdivisions(5);
+  h_effEtaPtPim->GetYaxis()->SetNdivisions(5);
+  h_effEtaPtPim->SetContour(99);
+  h_effEtaPtPim->Draw("COLZ");
+  TLatex* pimlabel = new TLatex(0.62, 0.75, "#pi^{#minus}");
+  pimlabel->SetNDC();
+  pimlabel->SetTextSize(40);
+  pimlabel->SetTextFont(43);
+  pimlabel->SetTextColor(kBlack);
+  pimlabel->Draw("same");
+  TLatex* r43fig5 = new TLatex(0.66,0.93, "#bf{EPIC}");
+  r43fig5->SetNDC();
+  r43fig5->SetTextSize(15);
+  r43fig5->SetTextFont(43);
+  r43fig5->SetTextColor(kBlack);
+  r43fig5->Draw("same");
+  TLatex* r44fig5b = new TLatex(0.01, 0.93, "W>2 GeV");
+  r44fig5b->SetNDC();
+  r44fig5b->SetTextSize(15);
+  r44fig5b->SetTextFont(43);
+  r44fig5b->SetTextColor(kBlack);
+  r44fig5b->Draw("same");
+
+  p5->cd(4);
+  TVirtualPad* p54 = p5->cd(4);
+  p54->SetTopMargin(0.05);
+  p54->SetRightMargin(0);
+  p54->SetLeftMargin(0.2);
+  p54->SetBottomMargin(0.21);
+  h_effPhiEtaPi->GetXaxis()->SetRangeUser(0,6.2);
+  h_effPhiEtaPi->GetYaxis()->SetRangeUser(4,6);
+  h_effPhiEtaPi->GetZaxis()->SetRangeUser(0,1);
+  h_effPhiEtaPi->GetXaxis()->SetNdivisions(8);
+  h_effPhiEtaPi->GetYaxis()->SetNdivisions(5);
+  h_effPhiEtaPi->SetContour(99);
+  h_effPhiEtaPi->Draw("COLZ");
+  TLatex* pilabela = new TLatex(0.3, 0.82, "#pi^{#pm}");
+  TLatex* pilabelb = new TLatex(0.5, 0.84, "(p_{T}>0.2 GeV/c)");
+  pilabela->SetNDC();
+  pilabelb->SetNDC();
+  pilabela->SetTextSize(40);
+  pilabelb->SetTextSize(15);
+  pilabela->SetTextFont(43);
+  pilabelb->SetTextFont(43);
+  pilabela->SetTextColor(kWhite);
+  pilabelb->SetTextColor(kWhite);
+  pilabela->Draw("same");
+  pilabelb->Draw("same");
+
+  p5->cd(5);
+  TVirtualPad* p55 = p5->cd(5);
+  p55->SetTopMargin(0.05);
+  p55->SetRightMargin(0);
+  p55->SetLeftMargin(0);
+  p55->SetBottomMargin(0.2);
+  h_effPhiEtaPip->GetXaxis()->SetRangeUser(0.15,6.2);
+  h_effPhiEtaPip->GetYaxis()->SetRangeUser(4,6);
+  h_effPhiEtaPip->GetZaxis()->SetRangeUser(0,1);
+  h_effPhiEtaPip->GetXaxis()->SetNdivisions(8);
+  h_effPhiEtaPip->GetYaxis()->SetNdivisions(5);
+  h_effPhiEtaPip->SetContour(99);
+  h_effPhiEtaPip->Draw("COLZ");
+  TLatex* piplabela = new TLatex(0.2, 0.82, "#pi^{#plus}");
+  TLatex* piplabelb = new TLatex(0.4, 0.84, "(p_{T}>0.2 GeV/c)");
+  piplabela->SetNDC();
+  piplabelb->SetNDC();
+  piplabela->SetTextSize(40);
+  piplabelb->SetTextSize(15);
+  piplabela->SetTextFont(43);
+  piplabelb->SetTextFont(43);
+  piplabela->SetTextColor(kWhite);
+  piplabelb->SetTextColor(kWhite);
+  piplabela->Draw("same");
+  piplabelb->Draw("same");
+
+  p5->cd(6);
+  TVirtualPad* p56 = p5->cd(6);
+  p56->SetTopMargin(0.05);
+  p56->SetRightMargin(0.2);
+  p56->SetLeftMargin(0);
+  p56->SetBottomMargin(0.2);
+  h_effPhiEtaPim->SetTitle(";#phi (rad);;efficiency");
+  h_effPhiEtaPim->GetXaxis()->SetRangeUser(0.15,6.2);
+  h_effPhiEtaPim->GetYaxis()->SetRangeUser(4,6);
+  h_effPhiEtaPim->GetZaxis()->SetRangeUser(0,1);
+  h_effPhiEtaPim->GetXaxis()->SetNdivisions(8);
+  h_effPhiEtaPim->GetYaxis()->SetNdivisions(5);
+  h_effPhiEtaPim->SetContour(99);
+  h_effPhiEtaPim->Draw("COLZ");
+  TLatex* pimlabela = new TLatex(0.1, 0.82, "#pi^{#minus}");
+  TLatex* pimlabelb = new TLatex(0.25, 0.84, "(p_{T}>0.2 GeV/c)");
+  pimlabela->SetNDC();
+  pimlabelb->SetNDC();
+  pimlabela->SetTextSize(40);
+  pimlabelb->SetTextSize(15);
+  pimlabela->SetTextFont(43);
+  pimlabelb->SetTextFont(43);
+  pimlabela->SetTextColor(kWhite);
+  pimlabelb->SetTextColor(kWhite);
+  pimlabela->Draw("same");
+  pimlabelb->Draw("same");
+
+  TString figure3name = figure_directory+"/benchmark_rho_efficiencies.pdf";
+  c5->Print(figure3name);
+}

benchmarks/your_benchmark/reconstruct.sh

+#!/bin/bash
+source strict-mode.sh
+source benchmarks/your_benchmark/setup.config $*
+
+# Reconstruct
+if [ ${RECO} == "eicrecon" ] ; then
+  eicrecon ${OUTPUT_FILE} -Ppodio:output_file=${REC_FILE}
+  if [[ "$?" -ne "0" ]] ; then
+    echo "ERROR running eicrecon"
+    exit 1
+  fi
+fi
+
+if [[ ${RECO} == "juggler" ]] ; then
+  gaudirun.py options/reconstruction.py || [ $? -eq 4 ]
+  if [ "$?" -ne "0" ] ; then
+    echo "ERROR running juggler"
+    exit 1
+  fi
+fi
+
+if [ -f jana.dot ] ; then cp jana.dot ${REC_FILE_BASE}.dot ; fi
+
+#rootls -t ${REC_FILE_BASE}.tree.edm4eic.root
+rootls -t ${REC_FILE}

benchmarks/your_benchmark/setup.config

+#!/bin/bash
+source strict-mode.sh
+
+export ENV_MODE=eicweb
+
+USE_SIMULATION_CAMPAIGN=true
+
+N_EVENTS=100
+
+FILE_BASE=sim_output/rho_10x100_uChannel_Q2of0to10_hiDiv.hepmc3.tree
+INPUT_FILE=root://dtn-eic.jlab.org//work/eic2/EPIC/EVGEN/EXCLUSIVE/UCHANNEL_RHO/10x100/rho_10x100_uChannel_Q2of0to10_hiDiv.hepmc3.tree.root
+OUTPUT_FILE=${FILE_BASE}.detectorsim.root
+
+REC_FILE_BASE=${FILE_BASE}.detectorsim.edm4eic
+REC_FILE=${REC_FILE_BASE}.root
+

benchmarks/your_benchmark/simulate.sh

+#!/bin/bash
+source strict-mode.sh
+source benchmarks/your_benchmark/setup.config $*
+
+if [ -f ${INPUT_FILE} ]; then
+  echo "ERROR: Input simulation file does ${INPUT_FILE} not exist."
+else
+  echo "GOOD: Input simulation file ${INPUT_FILE} exists!"
+fi
+
+# Simulate
+ddsim --runType batch \
+      -v WARNING \
+      --numberOfEvents ${N_EVENTS} \
+      --part.minimalKineticEnergy 100*GeV  \
+      --filter.tracker edep0 \
+      --compactFile ${DETECTOR_PATH}/${DETECTOR_CONFIG}.xml \
+      --inputFiles ${INPUT_FILE} \
+      --outputFile  ${OUTPUT_FILE}
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running ddsim"
+  exit 1
+fi