diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 62211ceb..fbedbf8d 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -111,25 +111,27 @@ common:detector: - runner_system_failure include: - - local: 'benchmarks/diffractive_vm/config.yml' + #- local: 'benchmarks/diffractive_vm/config.yml' - local: 'benchmarks/dis/config.yml' #- local: 'benchmarks/dvmp/config.yml' - - local: 'benchmarks/dvcs/config.yml' - - local: 'benchmarks/tcs/config.yml' - - local: 'benchmarks/u_omega/config.yml' + #- local: 'benchmarks/dvcs/config.yml' + #- local: 'benchmarks/tcs/config.yml' + #- local: 'benchmarks/u_omega/config.yml' + - local: 'benchmarks/u_rho/config.yml' - local: 'benchmarks/single/config.yml' - - local: 'benchmarks/backgrounds/config.yml' + #- local: 'benchmarks/backgrounds/config.yml' summary: stage: finish needs: - - "diffractive_vm:results" + #- "diffractive_vm:results" - "dis:results" - - "dvcs:results" - - "tcs:results" - - "u_omega:results" + #- "dvcs:results" + #- "tcs:results" + #- "u_omega:results" + - "u_channel_rho:results" - "single:results" - - "backgrounds:results" + #- "backgrounds:results" script: - collect_benchmarks.py - find results -print | sort | tee summary.txt diff --git a/Snakefile b/Snakefile index 7389e85f..93cbf7ee 100644 --- a/Snakefile +++ b/Snakefile @@ -18,3 +18,4 @@ root -l -b -q -e '.L {input}+' """ include: "benchmarks/diffractive_vm/Snakefile" +include: "benchmarks/u_rho/Snakefile" diff --git a/benchmarks/u_rho/BenchmarkPlotsExplained.pdf b/benchmarks/u_rho/BenchmarkPlotsExplained.pdf new file mode 100644 index 00000000..bdae335b Binary files /dev/null and b/benchmarks/u_rho/BenchmarkPlotsExplained.pdf differ diff --git a/benchmarks/u_rho/Snakefile b/benchmarks/u_rho/Snakefile new file mode 100644 index 00000000..e6323337 --- /dev/null +++ b/benchmarks/u_rho/Snakefile @@ -0,0 +1,133 @@ +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"], +) + + +rule uchannelrho_compile: + input: + ROOT_BUILD_DIR_PREFIX + "benchmarks/u_rho/analysis/uchannelrho_cxx.so", + + +rule uchannelrho_compile_manual: + input: + tar=HTTPRemoteProvider().remote("https://github.com/tectonic-typesetting/tectonic/releases/download/tectonic%400.15.0/tectonic-0.15.0-x86_64-unknown-linux-musl.tar.gz"), + cls=workflow.source_path("bench.cls"), + tex=workflow.source_path("bench.tex"), + output: + temp("tectonic"), + cls_tmp=temp("bench.cls"), + pdf="results/bench.pdf", + shell: """ +tar zxf {input.tar} +cp {input.cls} {output.cls_tmp} # copy to local directory +./tectonic {input.tex} --outdir="$(dirname {output.pdf})" +""" + + +rule uchannelrho_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: + "../../sim_output/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 uchannelrho_campaign_raw_get: + input: + S3.remote(f"eictest/EPIC/EVGEN/EXCLUSIVE/UCHANNEL_RHO/10x100/rho_10x100_uChannel_Q2of0to10_hiDiv.hepmc3.tree.root"), + output: + "../../sim_output/rho_10x100_uChannel_Q2of0to10_hiDiv.hepmc3.tree.root", + shell: + """ +echo "Getting afterburned event-generator file" +ln {input} {output} +ls ../../sim_output/* +""" + +rule uchannelrho_simulate: + input: + "../../sim_output/rho_10x100_uChannel_Q2of0to10_hiDiv.hepmc3.tree.root", + output: + "../../sim_output/rho_10x100_uChannel_Q2of0to10_hiDiv.hepmc3.tree.detectorsim.root", + params: + N_EVENTS=100 + shell: + """ +echo "Simulating detector response!" +bash benchmarks/u_rho/simulate.sh {input} {output} {params.N_EVENTS} +""" + + +rule uchannelrho_analysis: + input: + script="benchmarks/u_rho/analysis/uchannelrho.cxx", + #script_compiled=ROOT_BUILD_DIR_PREFIX + "benchmarks/u_rho/analysis/uchannelrho_cxx.so", + data="../../sim_output/campaign_24.07.0_rho_10x100_uChannel_Q2of0to10_hiDiv_{INDEX}_eicrecon.edm4eic.root", + output: + plots="../../sim_output/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 uchannelrho_combine: + input: + #lambda wildcards: [f"../../sim_output/campaign_24.07.0_{ix:04d}_eicrecon.edm4eic/plots.root" for ix in range(int(wildcards.NUM_FILES))], + lambda wildcards: expand( + "../../sim_output/campaign_24.07.0_{INDEX:04d}_eicrecon.edm4eic/plots.root", + INDEX=range(int(wildcards.N)), + ), + wildcard_constraints: + N="\d+", + output: + "../../sim_output/campaign_24.07.0_combined_{N}files_eicrecon.edm4eic.plots.root", + shell: + """ +hadd {output} {input} +""" + + +ruleorder: uchannelrho_combine > uchannelrho_analysis + + +rule uchannelrho_plots: + input: + script="benchmarks/u_rho/macros/plot_rho_physics_benchmark.C", + plots="../../sim_output/campaign_24.07.0_combined_{N}files_eicrecon.edm4eic.plots.root", + output: + "../../sim_output/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}")' +cat benchmark_output/*.json +#cp {input.plots}_figures/*.pdf figures/ +""" + + +# Couple examples of invocation: + +rule uchannelrho_run_over_a_campaign: + input: + "sim_output/figures/plots_benchmark_rho_dsigmadt.pdf", + message: + "See output in {input[0]}" + diff --git a/benchmarks/u_rho/analysis/pleaseIncludeMe.h b/benchmarks/u_rho/analysis/pleaseIncludeMe.h new file mode 100644 index 00000000..eb05ba11 --- /dev/null +++ b/benchmarks/u_rho/analysis/pleaseIncludeMe.h @@ -0,0 +1,116 @@ +#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 "nlohmann/json.hpp" +#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 + +auto getNtrk(const std::vector& parts) +{ + std::vector mult; + int n=0; + for(auto& i1 : parts){ + if(i1.charge!=0) n++; + } + mult.push_back( n ); + return mult; +} +auto getNtrkMC(const std::vector& parts) +{ + std::vector mult; + int n=0; + for(auto& i1 : parts){ + if(i1.charge!=0 && i1.generatorStatus==1) n++; + } + mult.push_back( n ); + return mult; +} +auto momenta_from_chargedparticles(const std::vector& parts) { + std::vector momenta; + for(auto& i1 : parts){ + TVector3 trk(i1.momentum.x,i1.momentum.y,i1.momentum.z); + if(i1.charge!=0) momenta.push_back(trk); + } + return momenta; +} +auto momenta_from_mcparticles(const std::vector& parts) { + std::vector momenta; + for(auto& i1 : parts){ + TVector3 trk(i1.momentum.x,i1.momentum.y,i1.momentum.z); + if(i1.charge!=0 && i1.generatorStatus==1) momenta.push_back(trk); + } + return momenta; +} +auto pt_resolution(const std::vector& mcs, + const std::vector& recos){ + + std::vector resolution; + for(auto& i1: recos){ + TVector3 trk(i1.momentum.x,i1.momentum.y,i1.momentum.z); + if(i1.charge==0) continue; + double minR=99; + TVector3 matchMCtrk(-99,-99,-99); + for(auto& i2 : mcs){ + TVector3 trkMC(i2.momentum.x,i2.momentum.y,i2.momentum.z); + if(i2.charge!=0 && i2.generatorStatus==1){ + if(trk.DeltaR(trkMC) < minR ){ + minR=trk.DeltaR(trkMC); + matchMCtrk=trkMC; + } + } + } + double res= (matchMCtrk.Perp()-trk.Perp()) / matchMCtrk.Perp(); + resolution.push_back( res ); + + } + + return resolution; +} +auto getPt(const std::vector& tracks) +{ + std::vector Pt; + for(auto& i1 : tracks){Pt.push_back(i1.Perp());} + return Pt; +} +auto getEta(const std::vector& tracks) +{ + std::vector eta; + for(auto& i1 : tracks){eta.push_back(i1.Eta());} + return eta; +} +auto getPhi(const std::vector& tracks) +{ + std::vector Phi; + for(auto& i1 : tracks){Phi.push_back(i1.Phi());} + return Phi; +} diff --git a/benchmarks/u_rho/analysis/uchannelrho.cxx b/benchmarks/u_rho/analysis/uchannelrho.cxx new file mode 100644 index 00000000..b787421a --- /dev/null +++ b/benchmarks/u_rho/analysis/uchannelrho.cxx @@ -0,0 +1,567 @@ +#include "pleaseIncludeMe.h" + +auto giveme_t_method_L(TLorentzVector eIn, + TLorentzVector eOut, + TLorentzVector pIn, + TLorentzVector vmOut) +{ + TLorentzVector aInVec(pIn.Px(),pIn.Py(),pIn.Pz(),sqrt(pIn.Px()*pIn.Px() + pIn.Py()*pIn.Py() + pIn.Pz()*pIn.Pz() + MASS_PROTON*MASS_PROTON) ); + double method_L = 0; + TLorentzVector a_beam_scattered = aInVec-(vmOut+eOut-eIn); + double p_Aplus = a_beam_scattered.E()+a_beam_scattered.Pz(); + double p_TAsquared = TMath::Power(a_beam_scattered.Pt(),2); + double p_Aminus = (MASS_PROTON*MASS_PROTON + p_TAsquared) / p_Aplus; + TLorentzVector a_beam_scattered_corr; + a_beam_scattered_corr.SetPxPyPzE(a_beam_scattered.Px(),a_beam_scattered.Py(),(p_Aplus-p_Aminus)/2., (p_Aplus+p_Aminus)/2. ); + method_L = -(a_beam_scattered_corr-aInVec).Mag2(); + + return method_L; +} + +auto giveme_u(TLorentzVector pIn, TLorentzVector vmOut){ + double uvalue = -(vmOut-pIn).Mag2(); + return uvalue; +} + +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_endx_array = {tree_reader, "MCParticles.endpoint.x"}; +TTreeReaderArray mc_endy_array = {tree_reader, "MCParticles.endpoint.y"}; +TTreeReaderArray mc_endz_array = {tree_reader, "MCParticles.endpoint.z"}; +TTreeReaderArray mc_mass_array = {tree_reader, "MCParticles.mass"}; +TTreeReaderArray mc_pdg_array = {tree_reader, "MCParticles.PDG"}; + +//Reconstructed EcalEndcapNClusters +TTreeReaderArray em_energy_array = {tree_reader, "EcalEndcapNClusters.energy"}; +TTreeReaderArray em_x_array = {tree_reader, "EcalEndcapNClusters.position.x"}; +TTreeReaderArray em_y_array = {tree_reader, "EcalEndcapNClusters.position.y"}; +TTreeReaderArray emhits_x_array = {tree_reader, "EcalEndcapNRecHits.position.x"}; +TTreeReaderArray emhits_y_array = {tree_reader, "EcalEndcapNRecHits.position.y"}; +TTreeReaderArray emhits_energy_array = {tree_reader, "EcalEndcapNRecHits.energy"}; + +TTreeReaderArray em_rec_id_array = {tree_reader, "EcalEndcapNClusterAssociations.recID"}; +TTreeReaderArray em_sim_id_array = {tree_reader, "EcalEndcapNClusterAssociations.simID"}; + +//B0 Tracker hits +TTreeReaderArray reco_x_array = {tree_reader, "MCParticles.endpoint.x"};//= {tree_reader, "B0TrackerHits.position.x"}; +TTreeReaderArray reco_y_array = {tree_reader, "MCParticles.endpoint.y"};//= {tree_reader, "B0TrackerHits.position.y"}; +TTreeReaderArray reco_z_array = {tree_reader, "MCParticles.endpoint.z"};//= {tree_reader, "B0TrackerHits.position.z"}; + +// Reconstructed particles pz array for each reconstructed particle +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 rec_id = {tree_reader, "ReconstructedChargedParticleAssociations.recID"}; +TTreeReaderArray sim_id = {tree_reader, "ReconstructedChargedParticleAssociations.simID"}; + +TTreeReaderArray reco_type = {tree_reader,"ReconstructedChargedParticles.type"}; + +TTreeReaderArray reco_PDG = {tree_reader,"ReconstructedChargedParticles.PDG"}; + +TString output_name_dir = outputfile; +cout << "Output file = " << output_name_dir << endl; +TFile* output = new TFile(output_name_dir,"RECREATE"); + +//events +TH1D* h_Q2_e = new TH1D("h_Q2_e",";Q^{2}_{e,MC}",100,0,20); +TH1D* h_y_e = new TH1D("h_y_e",";y_{e,MC}",100,0,1); +TH1D* h_energy_MC = new TH1D("h_energy_MC",";E_{MC} (GeV)",100,0,20); +TH1D* h_t_MC = new TH1D("h_t_MC",";t_{MC}; counts",1000,0,5); + +TH1D* h_Q2REC_e = new TH1D("h_Q2REC_e",";Q^{2}_{e,REC}",100,0,20); +TH1D* h_yREC_e = new TH1D("h_yREC_e",";y_{e,REC}",100,0,1); +TH1D* h_energy_REC = new TH1D("h_energy_REC",";E_{REC} (GeV)",100,0,20); +TH1D* h_trk_energy_REC = new TH1D("h_trk_energy_REC",";E_{REC} (GeV)",100,0,20); +TH1D* h_trk_Epz_REC = new TH1D("h_trk_Epz_REC",";E - p_{z} (GeV)",200,0,50); + +//track +TH1D* h_eta = new TH1D("h_eta",";#eta",100,-5,5); +TH2D* h_trk_energy_res = new TH2D("h_trk_energy_res",";E_{MC} (GeV); E_{MC}-E_{REC}/E_{MC} track-base ",100,0,20,1000,-1,1); +TH2D* h_trk_Pt_res = new TH2D("h_trk_Pt_res",";p_{T,MC} (GeV); P_{T,MC}-P_{T,REC}/P_{T,MC} track-base ",100,0,15,1000,-1,1); +TH1D* h_Epz_REC = new TH1D("h_Epz_REC",";E - p_{z} (GeV)",200,0,50); +TH2D* h_Acceptance_REC = new TH2D("h_Acceptance_REC",";#eta; p_{T}(GeV/c REC)",500,0,10,500,0,10); +TH2D* h_Acceptance_angular_REC = new TH2D("h_Acceptance_angular_REC",";#phi (radians); #theta (mrad)",100,0,6.4,200,0,100); +TH2D* h_Acceptance_angular_RECPI = new TH2D("h_Acceptance_angular_RECPI",";#phi (radians); #theta (mrad)",100,0,6.4,200,0,100); +TH2D* h_Acceptance_angular_RECPIP = new TH2D("h_Acceptance_angular_RECPIP",";#phi (radians); #theta (mrad)",100,0,6.4,200,0,100); +TH2D* h_Acceptance_angular_RECPIM = new TH2D("h_Acceptance_angular_RECPIM",";#phi (radians); #theta (mrad)",100,0,6.4,200,0,100); +TH2D* h_Acceptance_xy_RECPI = new TH2D("h_Acceptance_xy_RECPI",";x(mm); y(mm)",1000,-1000,1000,1000,-1000,1000); +TH2D* h_Acceptance_angular_MC = new TH2D("h_Acceptance_angular_MC" ,";#phi (radians); #theta (mrad)",500,0,6.4,500,0,100); +TH2D* h_etavseta_MC = new TH2D("h_etavseta_MC",";#eta(#pi^{+} MC); #eta(#pi^{-} MC)",250,0,10,250,0,10); +TH2D* h_etavseta_REC = new TH2D("h_etavseta_REC",";#eta(#pi^{+} REC); #eta(#pi^{-} REC)",250,0,10,250,0,10); + +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); + +TH2D* h_RecoMomPi = new TH2D("h_RecoMomPi", ";p (GeV/c) MC;p (GeV/c) reco.",100,0,100,100,0,100); +TH2D* h_RecoMomPip = new TH2D("h_RecoMomPip",";p (GeV/c) MC;p (GeV/c) reco.",100,0,100,100,0,100); +TH2D* h_RecoMomPim = new TH2D("h_RecoMomPim",";p (GeV/c) MC;p (GeV/c) reco.",100,0,100,100,0,100); + +TH2D* h_RecoTransMomPi = new TH2D("h_RecoTransMomPi", ";p_{T} (GeV/c) MC;p_{T} (GeV/c) reco.",100,0,1.6,100,0,1.6); +TH2D* h_RecoTransMomPip = new TH2D("h_RecoTransMomPip",";p_{T} (GeV/c) MC;p_{T} (GeV/c) reco.",100,0,1.6,100,0,1.6); +TH2D* h_RecoTransMomPim = new TH2D("h_RecoTransMomPim",";p_{T} (GeV/c) MC;p_{T} (GeV/c) reco.",100,0,1.6,100,0,1.6); + +TH1D* h_PiMomRecoQuality = new TH1D("h_PiMomRecoQuality",";(p_{reco}-p_{MC})/p_{MC};counts",500,-1,1); +TH1D* h_PipMomRecoQuality = new TH1D("h_PipMomRecoQuality",";(p_{reco}-p_{MC})/p_{MC};counts",500,-1,1); +TH1D* h_PimMomRecoQuality = new TH1D("h_PimMomRecoQuality",";(p_{reco}-p_{MC})/p_{MC};counts",500,-1,1); +TH1D* h_PiTransMomRecoQuality = new TH1D("h_PiTransMomRecoQuality",";(p_{T,reco}-p_{T,MC})/p_{T,MC};counts",500,-1,1); +TH1D* h_PipTransMomRecoQuality = new TH1D("h_PipTransMomRecoQuality",";(p_{T,reco}-p_{T,MC})/p_{T,MC};counts",500,-1,1); +TH1D* h_PimTransMomRecoQuality = new TH1D("h_PimTransMomRecoQuality",";(p_{T,reco}-p_{T,MC})/p_{T,MC};counts",500,-1,1); + +//VM & t +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); +TH1D* h_VM_mass_REC_justpionsB0 = new TH1D("h_VM_mass_REC_justpionsB0",";mass (GeV)",200,0,4); +TH1D* h_VM_mass_REC_notjustpionsB0 = new TH1D("h_VM_mass_REC_notjustpionsB0",";mass (GeV)",200,0,4); +TH1D* h_VM_pt_REC = new TH1D("h_VM_pt_REC",";p_{T} (GeV/c)",200,0,2); +TH2D* h_VM_res = new TH2D("h_VM_res",";p_{T,MC} (GeV); p_{T,MC}-E_{T,REC}/p_{T,MC}",100,0,2,1000,-1,1); +TH1D* h_t_REC = new TH1D("h_t_REC",";t_{REC} (GeV^{2}); counts",1000,0,5); +TH1D* h_t_trk_REC = new TH1D("h_t_trk_REC",";t_{REC}(GeV^{2}) track-base; counts",1000,0,5); +TH1D* h_t_combo_REC = new TH1D("h_t_combo_REC",";t_{combo,REC}(GeV^{2}); counts",1000,0,5); +TH2D* h_t_res = new TH2D("h_t_res",";t_{MC} (GeV^{2}); t_{MC}-t_{REC}/t_{MC}",1000,0,5,1000,-10,10); +TH2D* h_trk_t_res = new TH2D("h_trk_t_res",";t_{MC} (GeV^{2}); t_{MC}-t_{REC}/t_{MC} track-base",1000,0,5,1000,-10,10); +TH2D* h_t_2D = new TH2D("h_t_2D",";t_{MC} (GeV^{2}); t_{REC} (GeV^{2}) track-base",1000,0,5,1000,0,5); +TH2D* h_t_REC_2D = new TH2D("h_t_REC_2D",";t_{trk,REC} (GeV^{2}); t_{EEMC,REC} (GeV^{2})",1000,0,5,1000,0,5); +TH2D* h_t_RECMC_2D = new TH2D("h_t_RECMC_2D",";t_{MC} (GeV^{2}); t_{trk,REC} / t_{EEMC,REC} ",1000,0,5,200,-10,10); +TH2D* h_VM_endpointXY_MC = new TH2D("h_VM_endpointXY_MC",";x (cm); y (cm)",1000,-300,300,1000,-300,300); +TH1D* h_VM_endpointZ_MC = new TH1D("h_VM_endpointZ_MC",";z (cm); counts",1000,0,1000); +TH2D* h_u_REC_2D = new TH2D("h_u_REC_2D",";-#it{u}_{MC} (GeV^{2}); -#it{u}_{REC} (GeV^{2})",1000,0,5,1000,0,5); +double umin = -1.0; +double umax = 5.0; +int n_ubins = 100; +double u_binwidth = (umax-umin)/((double)n_ubins); +TH1D* h_u_REC = new TH1D("h_u_REC", ";-#it{u}_{REC} (GeV^{2}); dN/d|#it{u}| (GeV^{-2} scaled)",n_ubins,umin,umax); +TH1D* h_u_REC_justpions = new TH1D("h_u_REC_justpions", ";-#it{u}_{REC} (GeV^{2}); dN/d|#it{u}| (GeV^{-2} scaled)",n_ubins,umin,umax); +TH1D* h_u_REC_justpionsB0 = new TH1D("h_u_REC_justpionsB0", ";-#it{u}_{REC} (GeV^{2}); dN/d|#it{u}| (GeV^{-2} scaled)",n_ubins,umin,umax); +TH1D* h_u_REC_notjustpionsB0 = new TH1D("h_u_REC_notjustpionsB0", ";-#it{u}_{REC} (GeV^{2}); dN/d|#it{u}| (GeV^{-2} scaled)",n_ubins,umin,umax); +TH1D* h_u_MC = new TH1D("h_u_MC", ";-#it{u}_{MC} (GeV^{2}); dN/d|#it{u}| (GeV^{-2} scaled)",n_ubins,umin,umax); + +//energy clus +TH2D* h_emClus_position_REC = new TH2D("h_emClus_position_REC",";x (mm);y (mm)",80,-800,800,80,-800,800); +TH2D* h_emHits_position_REC = new TH2D("h_emHits_position_REC",";x (mm);y (mm)",80,-800,800,80,-800,800); +TH2D* h_energy_res = new TH2D("h_energy_res",";E_{MC} (GeV); E_{MC}-E_{REC}/E_{MC} emcal",100,0,20,1000,-1,1); +TH1D* h_energy_calibration_REC = new TH1D("h_energy_calibration_REC",";E (GeV)",200,0,2); +TH1D* h_EoverP_REC = new TH1D("h_EoverP_REC",";E/p",200,0,2); +TH1D* h_ClusOverHit_REC = new TH1D("h_ClusOverHit_REC",";cluster energy / new cluster energy",200,0,2); + +// +TH1D* h_numPositiveTracks = new TH1D("h_numPositiveTracks",";number of positive tracks;counts",10,-0.5,9.5); + +cout<<"There are "<GetEntries()<<" events!!!!!!!"<GetEntries()); +while (tree_reader.Next()) { + + /* + Beam particles + */ + TLorentzVector ebeam(0,0,0,0); + TLorentzVector pbeam(0,0,0,0); + + TLorentzVector vmMC(0,0,0,0); + TLorentzVector piplusMC(0,0,0,0); + TLorentzVector piminusMC(0,0,0,0); + + //MC level + TLorentzVector scatMC(0,0,0,0); + unsigned int mc_elect_index=-1; + double maxPt=-99.; + for(unsigned int imc=0;imcmaxPt){ + maxPt=mctrk.Perp(); + mc_elect_index=imc; + scatMC.SetVectM(mctrk,mc_mass_array[imc]); + } + if(mc_pdg_array[imc]==211 && mc_genStatus_array[imc]==1){ + piplusMC.SetVectM(mctrk,MASS_PION); + h_VM_endpointXY_MC->Fill(mc_endx_array[imc],mc_endy_array[imc]); + h_VM_endpointZ_MC->Fill(mc_endz_array[imc]); + double phi = mctrk.Phi()>0 ? mctrk.Phi() : mctrk.Phi()+6.2831853; + h_Acceptance_angular_MC->Fill(phi,1000.0*mctrk.Theta()); + } + if(mc_pdg_array[imc]==-211 && mc_genStatus_array[imc]==1){ + piminusMC.SetVectM(mctrk,MASS_PION); + h_VM_endpointXY_MC->Fill(mc_endx_array[imc],mc_endy_array[imc]); + h_VM_endpointZ_MC->Fill(mc_endz_array[imc]); + double phi = mctrk.Phi()>0 ? mctrk.Phi() : mctrk.Phi()+6.2831853; + h_Acceptance_angular_MC->Fill(phi,1000.0*mctrk.Theta()); + } + } + vmMC=piplusMC+piminusMC; + h_etavseta_MC->Fill(piplusMC.Eta(),piminusMC.Eta()); + //protection. + if(ebeam.E()==pbeam.E() && ebeam.E()==0) { + std::cout << "problem with MC incoming beams" << std::endl; + continue; + } + TLorentzVector qbeam=ebeam-scatMC; + double Q2=-(qbeam).Mag2(); + double pq=pbeam.Dot(qbeam); + double y= pq/pbeam.Dot(ebeam); + + //MC level phase space cut + //if(Q2<1.||Q2>10.) continue; + //if(y<0.01||y>0.85) continue; + + h_Q2_e->Fill(Q2); + h_y_e->Fill(y); + h_energy_MC->Fill(scatMC.E()); + + double t_MC=0.; + if(vmMC.E()!=0) + { + double method_E = -(qbeam-vmMC).Mag2(); + t_MC=method_E; + h_t_MC->Fill( method_E ); + 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()); + } + + //rec level + //leading cluster + double maxEnergy=-99.; + double xpos=-999.; + double ypos=-999.; + for(unsigned int iclus=0;iclusmaxEnergy){ + maxEnergy=em_energy_array[iclus]; + xpos=em_x_array[iclus]; + ypos=em_y_array[iclus]; + } + } + //leading hit energy + double maxHitEnergy=0.01;//threshold 10 MeV + double xhitpos=-999.; + double yhitpos=-999.; + int hit_index=-1; + for(unsigned int ihit=0;ihitmaxHitEnergy){ + maxHitEnergy=emhits_energy_array[ihit]; + xhitpos=emhits_x_array[ihit]; + yhitpos=emhits_y_array[ihit]; + hit_index=ihit; + } + } + //sum over all 3x3 towers around the leading tower + double xClus=xhitpos*maxHitEnergy; + double yClus=yhitpos*maxHitEnergy; + for(unsigned int ihit=0;ihit0.01) { + maxHitEnergy+=hitenergy;//clustering around leading tower 3 crystal = 60mm. + xClus+=x*hitenergy; + yClus+=y*hitenergy; + } + } + + h_ClusOverHit_REC->Fill( maxEnergy / maxHitEnergy ); + //weighted average cluster position. + xClus = xClus/maxHitEnergy; + yClus = yClus/maxHitEnergy; + double radius=sqrt(xClus*xClus+yClus*yClus); + //if( radius>550. ) continue; //geometric acceptance cut + //4.4% energy calibration. + double clusEnergy=1.044*maxHitEnergy; + + h_energy_REC->Fill(clusEnergy); + //ratio of reco / truth Energy + h_energy_calibration_REC->Fill( clusEnergy / scatMC.E() ); + //energy resolution + double res= (scatMC.E()-clusEnergy)/scatMC.E(); + h_energy_res->Fill(scatMC.E(), res); + h_emClus_position_REC->Fill(xpos,ypos);//default clustering position + h_emHits_position_REC->Fill(xClus,yClus); //self clustering position + + //association of rec level scat' e + int rec_elect_index=-1; + for(unsigned int i=0;imaxP){ + //track-base 4 vector + maxP=trk.Mag(); + scatREC.SetVectM(trk,MASS_ELECTRON); + + //use emcal energy to define 4 vector + double p = sqrt(clusEnergy*clusEnergy- MASS_ELECTRON*MASS_ELECTRON ); + double eta=scatREC.Eta(); + double phi=scatREC.Phi(); + double pt = TMath::Sin(scatREC.Theta())*p; + scatClusEREC.SetPtEtaPhiM(pt,eta,phi,MASS_ELECTRON); + } + } + //loop over track again; + int numpositivetracks = 0; + int failed = 0; + + for(unsigned int itrk=0;itrkpi+pi- daughters; + h_eta->Fill(trk.Eta()); + //if(fabs(trk.Eta())<3.0){ + if(reco_charge_array[itrk]>0){ + numpositivetracks++; + if ((sim_id[itrk - failed]==4 || sim_id[itrk - failed]==5) && reco_charge_array[itrk - failed]==1){ + // if(reco_PDG[itrk]==211){ + piplusREC.SetVectM(trk,MASS_PION); + isPiPlusFound=true; + } + if(sim_id[itrk - failed]==6){ + protonREC.SetVectM(trk,MASS_PROTON); + 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;} + //} + double pt = sqrt(reco_px_array[itrk]*reco_px_array[itrk] + reco_py_array[itrk]*reco_py_array[itrk]); + h_Acceptance_REC->Fill(fabs(trk.Eta()),pt); + double phi = trk.Phi()>0 ? trk.Phi() : trk.Phi()+6.2831853; + h_Acceptance_angular_REC->Fill(phi,1000.0*trk.Theta()); + if(sim_id[itrk - failed]==4 || sim_id[itrk - failed]==5) h_Acceptance_angular_RECPI->Fill(phi,1000.0*trk.Theta()); + if((sim_id[itrk - failed]==4 || sim_id[itrk - failed]==5) && reco_charge_array[itrk - failed]==1) h_Acceptance_angular_RECPIP->Fill(phi,1000.0*trk.Theta()); + if((sim_id[itrk - failed]==4 || sim_id[itrk - failed]==5) && reco_charge_array[itrk - failed]==-1) h_Acceptance_angular_RECPIM->Fill(phi,1000.0*trk.Theta()); + } + } + h_numPositiveTracks->Fill(numpositivetracks); + //loop over B0 hits + for(unsigned int ihit=0; ihitFill(hit.X(),hit.Y()); + } + //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 + 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); + // + if(isPiPlusFound) h_RecoMomPi ->Fill(piplusMC.P() ,piplusREC.P() ); + if(isPiMinusFound) h_RecoMomPi ->Fill(piminusMC.P(),piminusREC.P()); + if(isPiPlusFound) h_RecoMomPip->Fill(piplusMC.P() ,piplusREC.P() ); + if(isPiMinusFound) h_RecoMomPim->Fill(piminusMC.P(),piminusREC.P()); + // + if(isPiPlusFound) h_RecoTransMomPi ->Fill(piplusMC.Pt() ,piplusREC.Pt() ); + if(isPiMinusFound) h_RecoTransMomPi ->Fill(piminusMC.Pt(),piminusREC.Pt()); + if(isPiPlusFound) h_RecoTransMomPip->Fill(piplusMC.Pt() ,piplusREC.Pt() ); + if(isPiMinusFound) h_RecoTransMomPim->Fill(piminusMC.Pt(),piminusREC.Pt()); + // + if(isPiPlusFound){ + h_PiMomRecoQuality->Fill((piplusREC.P()-piplusMC.P())/piplusMC.P()); + h_PipMomRecoQuality->Fill((piplusREC.P()-piplusMC.P())/piplusMC.P()); + h_PiTransMomRecoQuality->Fill((piplusREC.Pt()-piplusMC.Pt())/piplusMC.Pt()); + h_PipTransMomRecoQuality->Fill((piplusREC.Pt()-piplusMC.Pt())/piplusMC.Pt()); + } + if(isPiMinusFound){ + h_PiMomRecoQuality->Fill((piminusREC.P()-piminusMC.P())/piminusMC.P()); + h_PimMomRecoQuality->Fill((piminusREC.P()-piminusMC.P())/piminusMC.P()); + h_PiTransMomRecoQuality->Fill((piminusREC.Pt()-piminusMC.Pt())/piminusMC.Pt()); + h_PimTransMomRecoQuality->Fill((piminusREC.Pt()-piminusMC.Pt())/piminusMC.Pt()); + } + + //track-base e' energy REC; + h_trk_energy_REC->Fill(scatMCmatchREC.E()); + + //track-base e' energy resolution; + res= (scatMC.E()-scatMCmatchREC.E())/scatMC.E(); + h_trk_energy_res->Fill(scatMC.E(), res); + + //track-base e' pt resolution; + res= (scatMC.Pt()-scatMCmatchREC.Pt())/scatMC.Pt(); + h_trk_Pt_res->Fill(scatMC.Pt(), res); + + //track-base Epz scat' e + double EpzREC= (scatMCmatchREC+hfs).E() - (scatMCmatchREC+hfs).Pz(); + h_trk_Epz_REC->Fill( EpzREC ); + + //EEMC cluster Epz scat' e + EpzREC= (scatClusEREC+hfs).E() - (scatClusEREC+hfs).Pz(); + h_Epz_REC->Fill( EpzREC ); + + //E over p + double EoverP=scatClusEREC.E() / scatMCmatchREC.P(); + h_EoverP_REC->Fill( EoverP ); + + //cluster-base DIS kine; + TLorentzVector qbeamREC=ebeam-scatClusEREC; + double Q2REC=-(qbeamREC).Mag2(); + double pqREC=pbeam.Dot(qbeamREC); + double yREC= pqREC/pbeam.Dot(ebeam); + h_Q2REC_e->Fill(Q2REC); + h_yREC_e->Fill(yREC); + + //Event selection: + //if( EpzREC<27||EpzREC>40 ) continue; + //if( EoverP<0.8||EoverP>1.18 ) continue; + + //MC level phase space cut + //if(Q2REC<1.||Q2REC>10.) continue; + //if(yREC<0.01||yREC>0.85) continue; + + //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(piplusREC.Eta()>3.9 && piminusREC.Eta()>3.9) h_VM_mass_REC_notjustpionsB0->Fill(rho_mass); + if(isPiMinusFound && isPiPlusFound){ + h_VM_mass_REC_justpions->Fill(rho_mass); + h_etavseta_REC->Fill(piplusREC.Eta(),piminusREC.Eta()); + if(piplusREC.Eta()>3.9 && piminusREC.Eta()>3.9) h_VM_mass_REC_justpionsB0->Fill(rho_mass); + } + h_VM_pt_REC->Fill(vmREC.Pt()); + + //2 versions: track and energy cluster: + double t_trk_REC = giveme_t_method_L(ebeam,scatMCmatchREC,pbeam,vmREC); + double t_REC = giveme_t_method_L(ebeam,scatClusEREC,pbeam,vmREC); + h_t_trk_REC->Fill( t_trk_REC ); + h_t_REC->Fill( t_REC ); + h_t_REC_2D->Fill(t_trk_REC,t_REC); + double u_REC = giveme_u(pbeam,vmREC); + double u_REC_fromproton = giveme_u(pbeam,vmRECfromproton); + double u_MC = giveme_u(pbeam,vmMC); + //4vector of VM; + if(piplusREC.E()!=0. && piminusREC.E()!=0.){ + h_u_REC->Fill(u_REC); + h_u_REC_justpions->Fill(u_REC); + } + if(protonRECasifpion.E()!=0. && piminusREC.E()!=0.){ + h_u_REC->Fill(u_REC_fromproton); + } + h_u_REC_2D->Fill(u_MC,u_REC); + /*h_u_REC->Fill(u_REC); + if(piplusREC.Eta()>3.9 && piminusREC.Eta()>3.9) h_u_REC_notjustpionsB0->Fill(u_REC); + if(isPiMinusFound && isPiPlusFound){ + h_u_REC_justpions->Fill(u_REC); + if(piplusREC.Eta()>3.9 && piminusREC.Eta()>3.9) h_u_REC_justpionsB0->Fill(u_REC); + } + */ + h_u_MC->Fill(u_MC); + if( (t_trk_REC/t_REC) > 0.5 && (t_trk_REC/t_REC) < 1.5 ){ + h_t_combo_REC->Fill( (t_trk_REC+t_REC)/2. );//w=1./(fabs(1.0-(t_trk_REC/t_REC))) + } + h_t_RECMC_2D->Fill(t_MC,t_trk_REC/t_REC); + + //t track resolution + res= (t_MC-t_trk_REC)/t_MC; + h_trk_t_res->Fill(t_MC, res); + + //t EEMC resolution; + res= (t_MC-t_REC)/t_MC; + h_t_res->Fill(t_MC, res); + + //2D t + h_t_2D->Fill(t_MC,t_trk_REC); + + //VM pt resolution; + res= (vmMC.Pt()-vmREC.Pt())/vmMC.Pt(); + h_VM_res->Fill(vmMC.Pt(), res); + +} +h_u_REC->Scale(1.0/u_binwidth); +h_u_REC_justpions->Scale(1.0/u_binwidth); +h_u_REC_justpionsB0->Scale(1.0/u_binwidth); +h_u_REC_notjustpionsB0->Scale(1.0/u_binwidth); +h_u_MC->Scale(1.0/u_binwidth); + +output->Write(); +output->Close(); + +return 0; +} + diff --git a/benchmarks/u_rho/analyze.sh b/benchmarks/u_rho/analyze.sh new file mode 100644 index 00000000..388d0035 --- /dev/null +++ b/benchmarks/u_rho/analyze.sh @@ -0,0 +1,24 @@ +#!/bin/bash +source strict-mode.sh + + +source benchmarks/u_rho/setup.config $* + +OUTPUT_PLOTS_DIR=sim_output/nocampaign +mkdir -p ${OUTPUT_PLOTS_DIR} +# Analyze +command time -v \ +root -l -b -q "benchmarks/u_rho/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/u_rho/macros/plot_rho_physics_benchmark.C(\"${OUTPUT_PLOTS_DIR}/plots.root\")" +cat benchmark_output/*.json diff --git a/benchmarks/u_rho/bench.cls b/benchmarks/u_rho/bench.cls new file mode 100644 index 00000000..13e49564 --- /dev/null +++ b/benchmarks/u_rho/bench.cls @@ -0,0 +1,166 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% BENCH.CLS % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\NeedsTeXFormat{LaTeX2e} +\LoadClass[11pt]{article} +\ProvidesClass{bench} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Making the Title, with definitions for date, abstract etc % +% This part written by David Cassel % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\def\@aabuffer{} +\def\author #1{\expandafter\def\expandafter\@aabuffer\expandafter +{\@aabuffer \small\rm #1\relax \par}} +\def\address#1{\expandafter\def\expandafter\@aabuffer\expandafter +{\@aabuffer \small\it #1\relax +\\ +\Photo +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +\par\vspace{1em}}} + +\def\maketitle{ +\begin{center} + {\bf \@title \par} + \vskip 2em % Vertical space after title. + \@aabuffer\relax +\end{center} \par +\gdef\@aabuffer{} +} + +\def\abstracts#1{ +\begin{center} +{\begin{minipage}{5.2truein} + \footnotesize + \parindent=0pt #1\par + \end{minipage}}\end{center} + \vskip 2em \par} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Margins, textwidths, indentations etc % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\fussy +\flushbottom +\parindent 0.25in +%\oddsidemargin 4mm +%\evensidemargin 4mm +%\topmargin=0mm +%\headheight=0mm +%\headsep=0mm +%\footskip=5mm +%\textheight = 240mm +%\textwidth = 160mm +%\RequirePackage[a4paper,hmargin=1.5cm,vmargin={1.5cm,1.5cm}]{geometry} +\RequirePackage[a4paper,margin=2.5cm]{geometry} +\RequirePackage[english]{babel} +\RequirePackage{graphicx,url} +\RequirePackage[colorlinks=true,urlcolor=blue,linkcolor=black,citecolor=black]{hyperref} +% for BibTeX - sorted numerical labels by order of first citation. +\bibliographystyle{unsrt} + +\def\section{\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus + -.2ex}{2.3ex plus .2ex}{\bf }} +\def\subsection{\@startsection{subsection}{2}{\z@}{-3.25ex plus -1ex minus + -.2ex}{1.5ex plus .2ex}{\it }} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% alpha footnotes, no running heads and silly citations. % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\def\thefootnote{\alph{footnote}} +\def\@makefnmark{{$\!^{\@thefnmark}$}} + +\pagestyle{empty} + + +\renewenvironment{thebibliography}[1] + {\begin{list}{\arabic{enumi}.} + {\usecounter{enumi}\setlength{\parsep}{0pt} + \setlength{\itemsep}{0pt} + \settowidth + {\labelwidth}{#1.}\sloppy}}{\end{list}} + +%--------------------------------------------------------------------------- +%FOLLOWING THREE COMMANDS ARE FOR `LIST' COMMAND. +\topsep=0in\parsep=0in\itemsep=0in + +\newcounter{arabiclistc} +\newenvironment{arabiclist} + {\setcounter{arabiclistc}{0} + \begin{list}{\arabic{arabiclistc}} + {\usecounter{arabiclistc} + \setlength{\parsep}{0pt} + \setlength{\itemsep}{0pt}}}{\end{list}} + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% ACKNOWLEDGEMENT: this portion is from John Hershberger % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\def\@citex[#1]#2{\if@filesw\immediate\write\@auxout + {\string\citation{#2}}\fi +\def\@citea{}\@cite{\@for\@citeb:=#2\do + {\@citea\def\@citea{,}\@ifundefined + {b@\@citeb}{{\bf ?}\@warning + {Citation `\@citeb' on page \thepage \space undefined}} + {\csname b@\@citeb\endcsname}}}{#1}} + +\newif\if@cghi +\def\cite{\@cghitrue\@ifnextchar [{\@tempswatrue + \@citex}{\@tempswafalse\@citex[]}} +\def\citelow{\@cghifalse\@ifnextchar [{\@tempswatrue + \@citex}{\@tempswafalse\@citex[]}} +\def\@cite#1#2{{$\!^{#1}$\if@tempswa\typeout + {IJCGA warning: optional citation argument + ignored: `#2'} \fi}} +\newcommand{\citeup}{\cite} + +\setcounter{secnumdepth}{2} + +\def\baselinestretch{1.0} +\ifx\selectfont\undefined +\@normalsize\else\let\glb@currsize=\relax\selectfont +\fi + +\ifx\selectfont\undefined +\def\@singlespacing{% +\def\baselinestretch{1}\ifx\@currsize\normalsize\@normalsize\else\@currsize\fi% +} +\else +\def\@singlespacing{\def\baselinestretch{1}\let\glb@currsize=\relax\selectfont} +\fi + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Footnote size table and figure captions % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +\long\def\@caption#1[#2]#3{\par\addcontentsline{\csname + ext@#1\endcsname}{#1}{\protect\numberline{\csname + the#1\endcsname}{\ignorespaces #2}}\begingroup + \@parboxrestore + \footnotesize + \expandafter\let\expandafter\@tempa\csname @make#1caption\endcsname + \ifx\@tempa\relax\let\@tempa\@makecaption\fi + \@tempa{\csname fnum@#1\endcsname}{\ignorespaces #3}\par + \endgroup} +% +% Here is the content of a .sty file containing definitions using the above +% hook. +% +% The following is the same as the \@makecaption in book.sty: +\long\def\@makefigurecaption#1#2{% + \vskip 10pt + \setbox\@tempboxa\hbox{#1 -- {\footnotesize #2}}% + \ifdim \wd\@tempboxa >\hsize #1 -- {\footnotesize #2}\par \else + \hbox to\hsize{\hfil\box\@tempboxa\hfil} % + \fi + } +% + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% END OF FILE bench.cls % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/benchmarks/u_rho/bench.tex b/benchmarks/u_rho/bench.tex new file mode 100644 index 00000000..a7712191 --- /dev/null +++ b/benchmarks/u_rho/bench.tex @@ -0,0 +1,55 @@ +%====================================================================% +% BENCH.TEX % +% Written by Zachary Sweger % +%====================================================================% + %\documentclass[final,3p]{elsarticle} +\documentclass{bench} + + + +% A useful Journal macro +\def\Journal#1#2#3#4{{#1} {\bf #2}, #3 (#4)} + +\NewDocumentCommand{\codeword}{v}{% +\texttt{\textcolor{black}{#1}}% +} + +% Some other macros used in the sample text +\def\st{\scriptstyle} +\def\sst{\scriptscriptstyle} +\def\epp{\epsilon^{\prime}} +\def\vep{\varepsilon} +\def\vp{{\bf p}} +\def\be{\begin{equation}} +\def\ee{\end{equation}} +\def\bea{\begin{eqnarray}} +\def\eea{\end{eqnarray}} +\def\CPbar{\hbox{{\rm CP}\hskip-1.80em{/}}} + + +\begin{document} +\title{$u$-channel $\rho^0$ Benchmark Figures} + + +\maketitle + +\codeword{benchmark_rho_mass.pdf}: +This figure shows the reconstruction of the $\rho^0$ mass. The \textbf{black} histogram is the invariant mass of each MC $\pi^+\pi^-$ pair after being processed by the afterburner. The \textbf{\textcolor{blue}{blue}} histogram is the invariant mass of reconstructed $\pi^+\pi^-$ pairs with no cuts on acceptance. PDG codes were used to select pions, although this PID is unrealistic. In the absence of PID, the $\rho^0$ will be reconstructed from each oppositely-charged track. The dominant combinatorial background from this approach comes from pairing protons with the $\pi^-$. This $m_{p\pi^-}$ background is shown by the red histogram. The sum of the signal and background is shown in +\textbf{\textcolor{magenta}{magenta}}. + + +\codeword{benchmark_rho_mass_cuts.pdf}: +This figure shows $\rho^0$ mass reconstruction for events in which both MC-level pions should be within the B0 acceptance (9$<\theta<$13 mrad with respect to the hadron beam pipe). The \textbf{black} histogram is the invariant mass of each MC $\pi^+\pi^-$ pair which passes this $\theta$ cut after being processed by the afterburner. The \textbf{\textcolor{magenta}{magenta}} histogram is the invariant mass of reconstructed $\pi^+\pi^-$ pairs for these same events. PDG codes were used to select pions. The \textbf{\textcolor{magenta}{magenta}} and \textbf{black} distributions were integrated over $0.6 +#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" + +#include "common_bench/benchmark.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/u_rho/macros/plot_rho_physics_benchmark.C b/benchmarks/u_rho/macros/plot_rho_physics_benchmark.C new file mode 100644 index 00000000..f868ccbb --- /dev/null +++ b/benchmarks/u_rho/macros/plot_rho_physics_benchmark.C @@ -0,0 +1,745 @@ +#include "RiceStyle.h" +using namespace std; + +int setbenchstatus(double eff){ + ///////////// Set benchmark status! + // create our test definition + // test_tag + common_bench::Test rho_reco_eff_test{ + { + {"name", "rho_reconstruction_efficiency"}, + {"title", "rho Reconstruction Efficiency for rho -> pi+pi- in the B0"}, + {"description", "u-channel rho->pi+pi- reconstruction efficiency " + "when both pions should be within B0 acceptance"}, + {"quantity", "efficiency"}, + {"target", "0.9"} + } + }; //these 2 need to be consistent + double eff_target = 0.9; //going to find a way to use the same variable + + if(eff<0 || eff>1){ + rho_reco_eff_test.error(-1); + }else if(eff > eff_target){ + rho_reco_eff_test.pass(eff); + }else{ + rho_reco_eff_test.fail(eff); + } + + // write out our test data + common_bench::write_test(rho_reco_eff_test, "./benchmark_output/u_rho_eff.json"); + return 0; +} + +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}"; + //t distribution + TH1D* h_t_MC = (TH1D*) file->Get("h_t_MC"); + TH1D* h_t_REC = (TH1D*) file->Get("h_t_REC"); + TH1D* h_t_trk_REC = (TH1D*) file->Get("h_t_trk_REC"); + TH1D* h_t_combo_REC = (TH1D*) file->Get("h_t_combo_REC"); + //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"); + //dN/du distribution + TH1D* h_dNdu_MC = (TH1D*) file->Get("h_u_MC"); + TH1D* h_dNdu_REC = (TH1D*) file->Get("h_u_REC"); + TH1D* h_dNdu_REC_justpions = (TH1D*) file->Get("h_u_REC_justpions"); + //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"); + //reco quality + TH2D* h_RecoMomPi = (TH2D*) file->Get("h_RecoMomPi"); + TH2D* h_RecoMomPip= (TH2D*) file->Get("h_RecoMomPip"); + TH2D* h_RecoMomPim= (TH2D*) file->Get("h_RecoMomPim"); + TH2D* h_RecoTransMomPi = (TH2D*) file->Get("h_RecoTransMomPi"); + TH2D* h_RecoTransMomPip= (TH2D*) file->Get("h_RecoTransMomPip"); + TH2D* h_RecoTransMomPim= (TH2D*) file->Get("h_RecoTransMomPim"); + + + TCanvas* c1 = new TCanvas("c1","c1",1,1,600,600); + gPad->SetLogy(1); + gPad->SetTicks(); + gPad->SetLeftMargin(0.15); + gPad->SetBottomMargin(0.15); + gPad->SetRightMargin(0.01); + TH1D* base1 = makeHist("base1", "", "|#it{t} | (GeV^{2})", "dN/d|#it{t} | (GeV^{-2}) ", 100,0,5.0,kBlack); + base1->GetYaxis()->SetRangeUser(8e-2, 8e5); + base1->GetXaxis()->SetTitleColor(kBlack); + fixedFontHist1D(base1,1.,1.2); + base1->GetYaxis()->SetTitleSize(base1->GetYaxis()->GetTitleSize()*1.5); + base1->GetXaxis()->SetTitleSize(base1->GetXaxis()->GetTitleSize()*1.5); + base1->GetYaxis()->SetLabelSize(base1->GetYaxis()->GetLabelSize()*1.5); + base1->GetXaxis()->SetLabelSize(base1->GetXaxis()->GetLabelSize()*1.5); + base1->GetXaxis()->SetNdivisions(4,4,0); + base1->GetYaxis()->SetNdivisions(5,5,0); + base1->Draw(); + + h_t_MC->Draw("same"); + + h_t_REC->SetMarkerStyle(20); + h_t_REC->Draw("PEsame"); + + h_t_trk_REC->SetFillColorAlpha(kBlue,0.4); + h_t_trk_REC->SetFillStyle(1001); + h_t_trk_REC->SetMarkerStyle(24); + h_t_trk_REC->SetMarkerColor(kBlue); + // h_t_trk_REC->Draw("PE3same"); + + h_t_combo_REC->SetFillColorAlpha(kRed,0.4); + h_t_combo_REC->SetFillStyle(1001); + h_t_combo_REC->SetMarkerStyle(24); + h_t_combo_REC->SetMarkerColor(kRed); + // h_t_combo_REC->Draw("PE3same"); + + TLatex* r42 = new TLatex(0.18, 0.91, "ep 10#times100 GeV"); + r42->SetNDC(); + r42->SetTextSize(22); + r42->SetTextFont(43); + r42->SetTextColor(kBlack); + r42->Draw("same"); + + TLatex* r43 = new TLatex(0.9,0.91, "#bf{EPIC}"); + r43->SetNDC(); + //r43->SetTextSize(0.04); + r43->SetTextSize(22); + r43->Draw("same"); + + TLatex* r44 = new TLatex(0.53, 0.78, "10^{-3}2 GeV"); + r44->SetNDC(); + r44->SetTextSize(20); + r44->SetTextFont(43); + r44->SetTextColor(kBlack); + r44->Draw("same"); + + 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); + r44_2->Draw("same"); + + TLegend *w7 = new TLegend(0.58,0.68,0.93,0.76); + w7->SetLineColor(kWhite); + w7->SetFillColor(0); + w7->SetTextSize(17); + w7->SetTextFont(45); + w7->AddEntry(h_t_MC, "eSTARlight "+vm_label+" MC ", "L"); + w7->AddEntry(h_t_REC, "eSTARlight "+vm_label+" RECO ", "P"); + w7->Draw("same"); + + //c1->Print("./sim_output/figures/benchmark_rho_dsigmadt.pdf"); + //TString figure1name = figure_directory+"/benchmark_rho_dsigmadt.pdf"; + //c1->Print(figure1name); + + 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_REC_justpions->SetFillColorAlpha(kViolet+10,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 figure2name = figure_directory+"/benchmark_rho_mass.pdf"; + c2->Print(figure2name); + + ///////////////////// Figure 3 + TCanvas* c3 = new TCanvas("c3","c3",1,1,600,600); + gPad->SetTicks(); + gPad->SetLogy(1); + gPad->SetLeftMargin(0.18); + gPad->SetBottomMargin(0.18); + gPad->SetTopMargin(0.10); + gPad->SetRightMargin(0.01); + TH1D* base3 = makeHist("base3", "", "-#it{u} (GeV^{2})", "dN/d#it{u} (GeV^{-2} scaled)", 100,-0.25,3.05,kBlack); + base3->GetYaxis()->SetRangeUser(0.5, 100*(h_dNdu_MC->GetMaximum())); + base3->GetXaxis()->SetTitleColor(kBlack); + fixedFontHist1D(base3,1.,1.2); + base3->GetYaxis()->SetTitleSize(base3->GetYaxis()->GetTitleSize()*1.5); + base3->GetXaxis()->SetTitleSize(base3->GetXaxis()->GetTitleSize()*1.5); + base3->GetYaxis()->SetLabelSize(base3->GetYaxis()->GetLabelSize()*1.5); + base3->GetXaxis()->SetLabelSize(base3->GetXaxis()->GetLabelSize()*1.5); + base3->GetYaxis()->SetTitleOffset(1.2); + base3->GetXaxis()->SetNdivisions(4,4,0); + base3->GetYaxis()->SetNdivisions(5,5,0); + base3->Draw(); + + h_dNdu_MC->SetLineColor(kBlack); + h_dNdu_REC->SetLineColor(kMagenta); + h_dNdu_REC_justpions->SetLineColor(kBlue); + h_dNdu_MC->SetLineWidth(2); + h_dNdu_REC->SetLineWidth(2); + h_dNdu_REC_justpions->SetLineWidth(2); + + h_dNdu_MC->Draw("HIST E same"); + h_dNdu_REC->Draw("HIST E same"); + h_dNdu_REC_justpions->Draw("HIST E same"); + + r42->Draw("same"); + r43->Draw("same"); + r44->Draw("same"); + r44_2->Draw("same"); + + TF1* fit_mc = new TF1("fit_mc", "[0]*exp([1]*x)",0.2, 1.2); + TF1* fit_rec = new TF1("fit_rec", "[0]*exp([1]*x)",0.2, 1.2); + TF1* fit_recpi = new TF1("fit_recpi", "[0]*exp([1]*x)",0.2, 1.2); + fit_mc->SetLineColor(kBlack); + fit_rec->SetLineColor(kMagenta); + fit_recpi->SetLineColor(kBlue); + fit_mc->SetLineWidth(3); + fit_rec->SetLineWidth(3); + fit_recpi->SetLineWidth(3); + TFitResultPtr r1 = h_dNdu_MC->Fit("fit_mc","RLMSIN"); + TFitResultPtr r2 = h_dNdu_REC->Fit("fit_rec","RLMSIN"); + TFitResultPtr r3 = h_dNdu_REC_justpions->Fit("fit_recpi","RLMSIN"); + double alpha_mc = r1->Parameter(1); + double alpha_rec = r2->Parameter(1); + double alpha_recpi = r3->Parameter(1); + fit_mc->Draw("SAME"); + fit_rec->Draw("SAME"); + fit_recpi->Draw("SAME"); + + TH1I* hwhite = new TH1I("hwhite","hwhite",1,0,1); + hwhite->SetLineColor(kWhite); + TLegend *w9 = new TLegend(0.53,0.61,0.86,0.76); + w9->SetLineColor(kWhite); + w9->SetFillColor(0); + w9->SetTextSize(17); + w9->SetTextFont(45); + w9->AddEntry(hwhite, "~exp[#alpha(-#it{u})]", "L"); + w9->AddEntry(h_dNdu_MC, vm_label+Form(" MC, #alpha=%.1fGeV^{-2}",alpha_mc), "L"); + w9->AddEntry(h_dNdu_REC, vm_label+Form(" reco. #alpha=%.1fGeV^{-2}",alpha_rec), "L"); + w9->AddEntry(h_dNdu_REC_justpions, vm_label+Form(" reco. (#pi^{#minus}#pi^{#plus}), #alpha=%.1fGeV^{-2}",alpha_recpi), "L"); + w9->Draw("same"); + + TString figure3name = figure_directory+"/benchmark_rho_dNdu.pdf"; + c3->Print(figure3name); + + ///////////////////// 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 figure4name = figure_directory+"/benchmark_rho_mass_cuts.pdf"; + c4->Print(figure4name); + + ///////////////////// 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(0.07); + 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 figure5name = figure_directory+"/benchmark_rho_efficiencies.pdf"; + c5->Print(figure5name); + + ///////////////////// Figure 6 + TCanvas* c6 = new TCanvas("c6","c6",1,1,700,560); + TPad* p6 = new TPad("p6","Pad5",0.,0.,1.,1.); + p6->Divide(3,2,0,0); + p6->Draw(); + gStyle->SetPalette(kBlueRedYellow); + gStyle->SetOptStat(0); + + h_RecoMomPi ->GetXaxis()->SetLabelSize(h_RecoMomPi ->GetXaxis()->GetLabelSize()*1.8); + h_RecoMomPip ->GetXaxis()->SetLabelSize(h_RecoMomPip ->GetXaxis()->GetLabelSize()*1.8); + h_RecoMomPim ->GetXaxis()->SetLabelSize(h_RecoMomPim ->GetXaxis()->GetLabelSize()*1.8); + h_RecoMomPi ->GetYaxis()->SetLabelSize(h_RecoMomPi ->GetYaxis()->GetLabelSize()*1.8); + h_RecoMomPim ->GetZaxis()->SetLabelSize(h_RecoMomPim ->GetZaxis()->GetLabelSize()*0.5); + h_RecoMomPim ->GetZaxis()->SetTitleSize(h_RecoMomPim ->GetZaxis()->GetTitleSize()*0.5); + h_RecoTransMomPi ->GetXaxis()->SetLabelSize(h_RecoTransMomPi ->GetXaxis()->GetLabelSize()*1.8); + h_RecoTransMomPip->GetXaxis()->SetLabelSize(h_RecoTransMomPip->GetXaxis()->GetLabelSize()*1.8); + h_RecoTransMomPim->GetXaxis()->SetLabelSize(h_RecoTransMomPim->GetXaxis()->GetLabelSize()*1.8); + h_RecoTransMomPi ->GetYaxis()->SetLabelSize(h_RecoTransMomPi ->GetYaxis()->GetLabelSize()*1.8); + h_RecoTransMomPim->GetZaxis()->SetLabelSize(h_RecoTransMomPim->GetZaxis()->GetLabelSize()*0.5); + h_RecoTransMomPim->GetZaxis()->SetTitleSize(h_RecoTransMomPim->GetZaxis()->GetTitleSize()*0.5); + + fixedFontHist1D(h_RecoMomPi,1.,1.2); + fixedFontHist1D(h_RecoMomPip,1.,1.2); + fixedFontHist1D(h_RecoMomPim,1.,1.2); + fixedFontHist1D(h_RecoTransMomPi,1.,1.2); + fixedFontHist1D(h_RecoTransMomPip,1.,1.2); + fixedFontHist1D(h_RecoTransMomPim,1.,1.2); + + double maxz = h_RecoMomPi->GetMaximum(); + double maxzt = h_RecoTransMomPi->GetMaximum(); + + p6->cd(1); + TVirtualPad* p61 = p6->cd(1); + p61->SetLogz(); + p61->SetTopMargin(0.08); + p61->SetRightMargin(0); + p61->SetLeftMargin(0.21); + p61->SetBottomMargin(0.2); + h_RecoMomPi->GetXaxis()->SetRangeUser(0,99); + h_RecoMomPi->GetYaxis()->SetRangeUser(0,99); + h_RecoMomPi->GetZaxis()->SetRangeUser(1,maxz); + h_RecoMomPi->GetXaxis()->SetNdivisions(5); + h_RecoMomPi->GetYaxis()->SetNdivisions(5); + h_RecoMomPi->SetContour(99); + h_RecoMomPi->Draw("COLZ"); + TBox* box1 = new TBox(7,70,32,95); + box1->SetLineColor(kBlack); + box1->SetFillColor(kWhite); + box1->Draw("l same"); + TLatex* pilabelz = new TLatex(0.3, 0.75, "#pi^{#pm}"); + pilabelz->SetNDC(); + pilabelz->SetTextSize(40); + pilabelz->SetTextFont(43); + pilabelz->SetTextColor(kBlack); + pilabelz->Draw("same"); + r44fig5c->Draw("same"); + + p6->cd(2); + TVirtualPad* p62 = p6->cd(2); + p62->SetLogz(); + p62->SetTopMargin(0.08); + p62->SetRightMargin(0); + p62->SetLeftMargin(0); + p62->SetBottomMargin(0.2); + h_RecoMomPip->GetXaxis()->SetRangeUser(1,99); + h_RecoMomPip->GetYaxis()->SetRangeUser(0,99); + h_RecoMomPip->GetZaxis()->SetRangeUser(1,maxz); + h_RecoMomPip->GetXaxis()->SetNdivisions(5); + h_RecoMomPip->GetYaxis()->SetNdivisions(5); + h_RecoMomPip->SetContour(99); + h_RecoMomPip->Draw("COLZ"); + TBox* box2 = new TBox(9,70,30,95); + box2->SetLineColor(kBlack); + box2->SetFillColor(kWhite); + box2->Draw("l same"); + TLatex* piplabelz = new TLatex(0.11, 0.75, "#pi^{#plus}"); + piplabelz->SetNDC(); + piplabelz->SetTextSize(40); + piplabelz->SetTextFont(43); + piplabelz->SetTextColor(kBlack); + piplabelz->Draw("same"); + r44fig5a->Draw("same"); + + p6->cd(3); + TVirtualPad* p63 = p6->cd(3); + p63->SetLogz(); + p63->SetTopMargin(0.08); + p63->SetRightMargin(0.2); + p63->SetLeftMargin(0); + p63->SetBottomMargin(0.2); + h_RecoMomPim->SetTitle(";p (GeV/c) MC;;counts"); + h_RecoMomPim->GetXaxis()->SetRangeUser(1,99); + h_RecoMomPim->GetYaxis()->SetRangeUser(0,99); + h_RecoMomPim->GetZaxis()->SetRangeUser(1,maxz); + h_RecoMomPim->GetXaxis()->SetNdivisions(5); + h_RecoMomPim->GetYaxis()->SetNdivisions(5); + h_RecoMomPim->SetContour(99); + h_RecoMomPim->Draw("COLZ"); + TBox* box3 = new TBox(12,70,40,95); + box3->SetLineColor(kBlack); + box3->SetFillColor(kWhite); + box3->Draw("l same"); + TLatex* pimlabelz = new TLatex(0.12, 0.75, "#pi^{#minus}"); + pimlabelz->SetNDC(); + pimlabelz->SetTextSize(40); + pimlabelz->SetTextFont(43); + pimlabelz->SetTextColor(kBlack); + pimlabelz->Draw("same"); + r43fig5->Draw("same"); + r44fig5b->Draw("same"); + + p6->cd(4); + TVirtualPad* p64 = p6->cd(4); + p64->SetLogz(); + p64->SetTopMargin(0.05); + p64->SetRightMargin(0); + p64->SetLeftMargin(0.21); + p64->SetBottomMargin(0.2); + h_RecoTransMomPi->GetXaxis()->SetRangeUser(0.01,1.49); + h_RecoTransMomPi->GetYaxis()->SetRangeUser(0,1.49); + h_RecoTransMomPi->GetZaxis()->SetRangeUser(1,maxzt); + h_RecoTransMomPi->GetXaxis()->SetNdivisions(5); + h_RecoTransMomPi->GetYaxis()->SetNdivisions(5); + h_RecoTransMomPi->SetContour(99); + h_RecoTransMomPi->Draw("COLZ"); + TBox* box4 = new TBox(0.15,1.15,0.5,1.45); + box4->SetLineColor(kBlack); + box4->SetFillColor(kWhite); + box4->Draw("l same"); + TLatex* pilabelaa = new TLatex(0.3, 0.8, "#pi^{#pm}"); + pilabelaa->SetNDC(); + pilabelaa->SetTextSize(40); + pilabelaa->SetTextFont(43); + pilabelaa->SetTextColor(kBlack); + pilabelaa->Draw("same"); + + p6->cd(5); + TVirtualPad* p65 = p6->cd(5); + p65->SetLogz(); + p65->SetTopMargin(0.05); + p65->SetRightMargin(0); + p65->SetLeftMargin(0); + p65->SetBottomMargin(0.2); + h_RecoTransMomPip->GetXaxis()->SetRangeUser(0.02,1.49); + h_RecoTransMomPip->GetYaxis()->SetRangeUser(0,1.49); + h_RecoTransMomPip->GetZaxis()->SetRangeUser(1,maxzt); + h_RecoTransMomPip->GetXaxis()->SetNdivisions(5); + h_RecoTransMomPip->GetYaxis()->SetNdivisions(5); + h_RecoTransMomPip->SetContour(99); + h_RecoTransMomPip->Draw("COLZ"); + TBox* box5 = new TBox(0.1,1.15,0.4,1.45); + box5->SetLineColor(kBlack); + box5->SetFillColor(kWhite); + box5->Draw("l same"); + TLatex* piplabelaz = new TLatex(0.1, 0.8, "#pi^{#plus}"); + piplabelaz->SetNDC(); + piplabelaz->SetTextSize(40); + piplabelaz->SetTextFont(43); + piplabelaz->SetTextColor(kBlack); + piplabelaz->Draw("same"); + + p6->cd(6); + TVirtualPad* p66 = p6->cd(6); + p66->SetLogz(); + p66->SetTopMargin(0.05); + p66->SetRightMargin(0.2); + p66->SetLeftMargin(0); + p66->SetBottomMargin(0.2); + h_RecoTransMomPim->SetTitle(";p_{T} (GeV/c) MC;;counts"); + h_RecoTransMomPim->GetXaxis()->SetRangeUser(0.02,1.49); + h_RecoTransMomPim->GetYaxis()->SetRangeUser(0,1.49); + h_RecoTransMomPim->GetZaxis()->SetRangeUser(1,maxz); + h_RecoTransMomPim->GetXaxis()->SetNdivisions(5); + h_RecoTransMomPim->GetYaxis()->SetNdivisions(5); + h_RecoTransMomPim->SetContour(99); + h_RecoTransMomPim->Draw("COLZ"); + TBox* box6 = new TBox(0.13,1.15,0.5,1.45); + box6->SetLineColor(kBlack); + box6->SetFillColor(kWhite); + box6->Draw("l same"); + TLatex* pimlabelaz = new TLatex(0.1, 0.8, "#pi^{#minus}"); + pimlabelaz->SetNDC(); + pimlabelaz->SetTextSize(40); + pimlabelaz->SetTextFont(43); + pimlabelaz->SetTextColor(kBlack); + pimlabelaz->Draw("same"); + + TString figure6name = figure_directory+"/benchmark_rho_recoquality.pdf"; + c6->Print(figure6name); + + double rhorecoeff = thiseff/100.0; + setbenchstatus(rhorecoeff); + +} diff --git a/benchmarks/u_rho/reconstruct.sh b/benchmarks/u_rho/reconstruct.sh new file mode 100644 index 00000000..2b7a4c10 --- /dev/null +++ b/benchmarks/u_rho/reconstruct.sh @@ -0,0 +1,26 @@ +#!/bin/bash +source strict-mode.sh + +source benchmarks/u_rho/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/u_rho/setup.config b/benchmarks/u_rho/setup.config new file mode 100644 index 00000000..5b091d4b --- /dev/null +++ b/benchmarks/u_rho/setup.config @@ -0,0 +1,14 @@ +#!/bin/bash +source strict-mode.sh + +USE_SIMULATION_CAMPAIGN=true + +N_EVENTS=600 + +FILE_BASE=sim_output/rho_10x100_uChannel_Q2of0to10_hiDiv.hepmc3.tree +#INPUT_FILE=${FILE_BASE}.root +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/u_rho/simulate.sh b/benchmarks/u_rho/simulate.sh new file mode 100644 index 00000000..ce22d3c1 --- /dev/null +++ b/benchmarks/u_rho/simulate.sh @@ -0,0 +1,24 @@ +#!/bin/bash +source strict-mode.sh + +source benchmarks/u_rho/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