From 8f8b8afb5c8d4cce53883f16a2d18875e0a0370c Mon Sep 17 00:00:00 2001 From: abhas-19 Date: Mon, 14 Oct 2024 19:58:01 +0530 Subject: [PATCH] I'm beefing up my benchmark! --- Snakefile | 1 + benchmarks/your_benchmark/Snakefile | 58 ++ .../your_benchmark/analysis/uchannelrho.cxx | 202 ++++++ benchmarks/your_benchmark/analyze.sh | 22 + benchmarks/your_benchmark/config.yml | 34 + benchmarks/your_benchmark/macros/RiceStyle.h | 684 ++++++++++++++++++ .../macros/plot_rho_physics_benchmark.C | 381 ++++++++++ benchmarks/your_benchmark/reconstruct.sh | 26 + benchmarks/your_benchmark/setup.config | 15 + benchmarks/your_benchmark/simulate.sh | 24 + 10 files changed, 1447 insertions(+) create mode 100644 benchmarks/your_benchmark/Snakefile create mode 100644 benchmarks/your_benchmark/analysis/uchannelrho.cxx create mode 100644 benchmarks/your_benchmark/analyze.sh create mode 100644 benchmarks/your_benchmark/macros/RiceStyle.h create mode 100644 benchmarks/your_benchmark/macros/plot_rho_physics_benchmark.C create mode 100644 benchmarks/your_benchmark/reconstruct.sh create mode 100644 benchmarks/your_benchmark/setup.config create mode 100644 benchmarks/your_benchmark/simulate.sh diff --git a/Snakefile b/Snakefile index 4a4b743b..4a6c17ab 100644 --- a/Snakefile +++ b/Snakefile @@ -44,3 +44,4 @@ ddsim \ include: "benchmarks/diffractive_vm/Snakefile" include: "benchmarks/dis/Snakefile" include: "benchmarks/demp/Snakefile" +include: "benchmarks/your_benchmark/Snakefile" diff --git a/benchmarks/your_benchmark/Snakefile b/benchmarks/your_benchmark/Snakefile new file mode 100644 index 00000000..2e24b4df --- /dev/null +++ b/benchmarks/your_benchmark/Snakefile @@ -0,0 +1,58 @@ +import os +from snakemake.remote.HTTP import RemoteProvider as HTTPRemoteProvider + +# 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: + output: + f"{OUTPUT_DIR}rho_10x100_uChannel_Q2of0to10_hiDiv.{{INDEX}}.eicrecon.tree.edm4eic.root", + shell: """ +xrdcp root://dtn-eic.jlab.org//work/eic2/EPIC/RECO/24.07.0/epic_craterlake/EXCLUSIVE/UCHANNEL_RHO/10x100/rho_10x100_uChannel_Q2of0to10_hiDiv.{wildcards.INDEX}.eicrecon.tree.edm4eic.root {output} +""" +rule your_benchmark_analysis: + input: + script=f"{BENCH_DIR}analysis/uchannelrho.cxx", + data=f"{OUTPUT_DIR}rho_10x100_uChannel_Q2of0to10_hiDiv.{{INDEX}}.eicrecon.tree.edm4eic.root", + output: + plots=f"{OUTPUT_DIR}campaign_24.07.0_{{INDEX}}.eicrecon.tree.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.tree.edm4eic/plots.root", + INDEX=range(int(wildcards.N)), + ), + wildcard_constraints: + N="\d+", + output: + f"{OUTPUT_DIR}campaign_24.07.0_combined_{{N}}files.eicrecon.tree.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.tree.edm4eic.plots.root", + output: + f"{OUTPUT_DIR}campaign_24.07.0_combined_{{N}}files.eicrecon.tree.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}")' +""" + diff --git a/benchmarks/your_benchmark/analysis/uchannelrho.cxx b/benchmarks/your_benchmark/analysis/uchannelrho.cxx new file mode 100644 index 00000000..91637837 --- /dev/null +++ b/benchmarks/your_benchmark/analysis/uchannelrho.cxx @@ -0,0 +1,202 @@ +#include +#include +#include +#include +#include +#include +#include + +#include "ROOT/RDataFrame.hxx" +#include +#include +#include +#include +#include +#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< mc_px_array = {tree_reader, "MCParticles.momentum.x"}; + TTreeReaderArray mc_py_array = {tree_reader, "MCParticles.momentum.y"}; + TTreeReaderArray mc_pz_array = {tree_reader, "MCParticles.momentum.z"}; + TTreeReaderArray mc_pdg_array= {tree_reader, "MCParticles.PDG"}; + + //Reco-level track attributes + TTreeReaderArray reco_px_array = {tree_reader, "ReconstructedChargedParticles.momentum.x"}; + TTreeReaderArray reco_py_array = {tree_reader, "ReconstructedChargedParticles.momentum.y"}; + TTreeReaderArray reco_pz_array = {tree_reader, "ReconstructedChargedParticles.momentum.z"}; + TTreeReaderArray reco_charge_array = {tree_reader, "ReconstructedChargedParticles.charge"}; + TTreeReaderArray reco_type = {tree_reader,"ReconstructedChargedParticles.type"}; + + TTreeReaderArray rec_id = {tree_reader, "ReconstructedChargedParticleAssociations.recID"}; + TTreeReaderArray 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 "<GetEntries()<<" events!!!!!!!"<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;imcFill(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;itrk0){ + 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; +} diff --git a/benchmarks/your_benchmark/analyze.sh b/benchmarks/your_benchmark/analyze.sh new file mode 100644 index 00000000..94f5d9d7 --- /dev/null +++ b/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 diff --git a/benchmarks/your_benchmark/config.yml b/benchmarks/your_benchmark/config.yml index 52aa9c84..bc861b27 100644 --- a/benchmarks/your_benchmark/config.yml +++ b/benchmarks/your_benchmark/config.yml @@ -7,13 +7,47 @@ 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 script: - echo "I will collect results here!" + needs: + - ["your_benchmark:simulate"] + script: + - 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.tree.edm4eic.plots_figures/benchmark_rho_mass.pdf + - cp ../../sim_output/campaign_24.07.0_combined_45files.eicrecon.tree.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!" diff --git a/benchmarks/your_benchmark/macros/RiceStyle.h b/benchmarks/your_benchmark/macros/RiceStyle.h new file mode 100644 index 00000000..6f9862e9 --- /dev/null +++ b/benchmarks/your_benchmark/macros/RiceStyle.h @@ -0,0 +1,684 @@ +#include +#include +#include +#include + +#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 makeMultiPad(const int num){//we only have 4,6,8 options for now + + cout << "num: "<< num << endl; + vector 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;iGetNbinsX();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;ibinGetNbinsX();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"); + + } + +} + diff --git a/benchmarks/your_benchmark/macros/plot_rho_physics_benchmark.C b/benchmarks/your_benchmark/macros/plot_rho_physics_benchmark.C new file mode 100644 index 00000000..7348acb7 --- /dev/null +++ b/benchmarks/your_benchmark/macros/plot_rho_physics_benchmark.C @@ -0,0 +1,381 @@ +#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}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; ibinGetNbinsX(); 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.6SetNDC(); + 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}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); +} diff --git a/benchmarks/your_benchmark/reconstruct.sh b/benchmarks/your_benchmark/reconstruct.sh new file mode 100644 index 00000000..286f5a8e --- /dev/null +++ b/benchmarks/your_benchmark/reconstruct.sh @@ -0,0 +1,26 @@ +#!/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} + diff --git a/benchmarks/your_benchmark/setup.config b/benchmarks/your_benchmark/setup.config new file mode 100644 index 00000000..3d165e35 --- /dev/null +++ b/benchmarks/your_benchmark/setup.config @@ -0,0 +1,15 @@ +#!/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 diff --git a/benchmarks/your_benchmark/simulate.sh b/benchmarks/your_benchmark/simulate.sh new file mode 100644 index 00000000..2cc8e023 --- /dev/null +++ b/benchmarks/your_benchmark/simulate.sh @@ -0,0 +1,24 @@ +#!/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 +