diff --git a/benchmarks/demp/Snakefile b/benchmarks/demp/Snakefile
index 920a0e91b948138bdef979da2c4997b4a3c27488..a397c3200bb591aa36b1135cfb86023f5ec07aba 100644
--- a/benchmarks/demp/Snakefile
+++ b/benchmarks/demp/Snakefile
@@ -1,12 +1,23 @@
import os
-
+#Compile the analysis and plotting scripts
rule demp_compile:
input:
ROOT_BUILD_DIR_PREFIX + "benchmarks/demp/analysis/demp_analysis_cxx.so",
ROOT_BUILD_DIR_PREFIX + "benchmarks/demp/analysis/demp_plots_cxx.so",
+#Process the simulated files based on the user-defined campaign
+rule demp_campaign_reco_get:
+ input:
+ lambda wildcards: f"EPIC/RECO/{wildcards.RELEASE_TAG}/{wildcards.DETECTOR_CONFIG}/EXCLUSIVE/DEMP/DEMPgen-1.1.0/{wildcards.EBEAM}x{wildcards.PBEAM}/pi+/DEMPgen-1.1.0_{wildcards.EBEAM}x{wildcards.PBEAM}_pi+_1.000{wildcards.INDEX}.eicrecon.tree.edm4eic.root",
+ output:
+ temp("reco/{DETECTOR_CONFIG}/campaign_{RELEASE_TAG}_demp_{EBEAM}on{PBEAM}_{INDEX}.edm4eic.root"),
+ shell:
+ """
+mv {input} {output}
+"""
+#Process the afterburned files through the simulations
rule demp_sim:
input:
warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root",
@@ -27,11 +38,11 @@ ddsim \
-v WARNING \
--numberOfEvents {params.N_EVENTS} \
--compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
- --inputFiles root://dtn-eic.jlab.org//work/eic2/EPIC/EVGEN/EXCLUSIVE/DEMP/{wildcards.EBEAM}on{wildcards.PBEAM}/eic_DEMPGen_{wildcards.EBEAM}on{wildcards.PBEAM}_ip6_pi+_1B_{wildcards.INDEX}.hepmc3.tree.root \
+ --inputFiles root://dtn-eic.jlab.org//work/eic2/EPIC/EVGEN/EXCLUSIVE/DEMP/DEMPgen-1.1.0/{wildcards.EBEAM}x{wildcards.PBEAM}/pi+/DEMPgen-1.1.0_{wildcards.EBEAM}x{wildcards.PBEAM}_pi+_{wildcards.INDEX}.hepmc3.tree.root \
--outputFile {output}
"""
-
+#Process the files produced in the previous step through eicrecon
rule demp_reco:
input:
"sim/{DETECTOR_CONFIG}/demp_{EBEAM}on{PBEAM}_{INDEX}.edm4hep.root",
@@ -46,16 +57,17 @@ rule demp_reco:
DETECTOR_CONFIG={wildcards.DETECTOR_CONFIG} eicrecon {input} -Ppodio:output_file={output}
"""
-
+#Process the files (either from the campaign or eicrecon) through the analysis script
rule demp_analysis:
input:
script="benchmarks/demp/analysis/demp_analysis.cxx",
script_compiled=ROOT_BUILD_DIR_PREFIX + "benchmarks/demp/analysis/demp_analysis_cxx.so",
- data="reco/{DETECTOR_CONFIG}/demp_{EBEAM}on{PBEAM}_{INDEX}.edm4eic.root",
+ data="reco/{DETECTOR_CONFIG}/{PREFIX}demp_{EBEAM}on{PBEAM}_{INDEX}.edm4eic.root",
output:
- config="results/{DETECTOR_CONFIG}/demp/demp_{EBEAM}on{PBEAM}_{INDEX}/config.json",
- hists="results/{DETECTOR_CONFIG}/demp/demp_{EBEAM}on{PBEAM}_{INDEX}/hists.root",
+ config="results/{DETECTOR_CONFIG}/demp/{PREFIX}demp_{EBEAM}on{PBEAM}_{INDEX}/config.json",
+ hists="results/{DETECTOR_CONFIG}/demp/{PREFIX}demp_{EBEAM}on{PBEAM}_{INDEX}/hists.root",
wildcard_constraints:
+ PREFIX= ".*",
EBEAM="\d+",
PBEAM="\d+",
INDEX="\d+",
@@ -73,14 +85,15 @@ EOF
root -l -b -q '{input.script}+("{output.config}")'
"""
-
+#Merge all the files produced in the previous step
rule demp_combine:
input:
- lambda wildcards: [f"results/{wildcards.DETECTOR_CONFIG}/demp/demp_{wildcards.EBEAM}on{wildcards.PBEAM}_{ix}/hists.root" for ix in range(1,int(wildcards.NUM_FILES)+1)],
+ lambda wildcards: [f"results/{wildcards.DETECTOR_CONFIG}/demp/{wildcards.PREFIX}demp_{wildcards.EBEAM}on{wildcards.PBEAM}_{ix}/hists.root" for ix in range(1,int(wildcards.NUM_FILES)+1)],
output:
- config="results/{DETECTOR_CONFIG}/demp/demp_{EBEAM}on{PBEAM}_combined_{NUM_FILES}/config.json",
- hists="results/{DETECTOR_CONFIG}/demp/demp_{EBEAM}on{PBEAM}_combined_{NUM_FILES}/hists.root",
+ config="results/{DETECTOR_CONFIG}/demp/{PREFIX}demp_{EBEAM}on{PBEAM}_combined_{NUM_FILES}/config.json",
+ hists="results/{DETECTOR_CONFIG}/demp/{PREFIX}demp_{EBEAM}on{PBEAM}_combined_{NUM_FILES}/hists.root",
wildcard_constraints:
+ PREFIX= ".*",
EBEAM="\d+",
PBEAM="\d+",
NUM_FILES="\d+",
@@ -92,21 +105,23 @@ cat > {output.config} <<EOF
"detector": "{wildcards.DETECTOR_CONFIG}",
"ebeam": {wildcards.EBEAM},
"pbeam": {wildcards.PBEAM},
+ "nfiles": {wildcards.NUM_FILES},
"output_prefix": "$(dirname "{output.hists}")/plots"
}}
EOF
hadd {output.hists} {input}
"""
-
+#Process the merged file through the plotting script
rule demp_plots:
input:
script="benchmarks/demp/analysis/demp_plots.cxx",
script_compiled=ROOT_BUILD_DIR_PREFIX + "benchmarks/demp/analysis/demp_plots_cxx.so",
- config="results/{DETECTOR_CONFIG}/demp/demp_{EBEAM}on{PBEAM}_combined_{NUM_FILES}/config.json",
+ config="results/{DETECTOR_CONFIG}/demp/{PREFIX}demp_{EBEAM}on{PBEAM}_combined_{NUM_FILES}/config.json",
output:
- "results/{DETECTOR_CONFIG}/demp/demp_{EBEAM}on{PBEAM}_combined_{NUM_FILES}/plots.pdf"
+ "results/{DETECTOR_CONFIG}/demp/{PREFIX}demp_{EBEAM}on{PBEAM}_combined_{NUM_FILES}/plots.pdf"
wildcard_constraints:
+ PREFIX= ".*",
EBEAM="\d+",
PBEAM="\d+",
NUM_FILES="\d+",
@@ -115,10 +130,16 @@ rule demp_plots:
root -l -b -q '{input.script}+("{input.config}")'
"""
-# Example of invocation
+#Example of invocation
rule demp_run_locally:
input:
"results/" + os.environ["DETECTOR_CONFIG"] + "/demp/demp_5on41_combined_5/plots.pdf",
message:
"See output in {input[0]}"
+rule demp_run_campaign:
+ input:
+ "results/epic_craterlake/demp/campaign_24.06.0_demp_5on41_combined_5/plots.pdf",
+ message:
+ "See output in {input[0]}"
+
diff --git a/benchmarks/demp/analysis/demp_analysis.cxx b/benchmarks/demp/analysis/demp_analysis.cxx
index b57beea181d8546044fb673ad6986da8460a712e..3095c2b0355d84f3b893121ae18999904eb12813 100644
--- a/benchmarks/demp/analysis/demp_analysis.cxx
+++ b/benchmarks/demp/analysis/demp_analysis.cxx
@@ -20,9 +20,10 @@
void demp_analysis(const std::string& config_name)
{
-
- //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------
- // read our configuration
+
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
+ // Read our configuration
std::ifstream config_file{config_name};
nlohmann::json config;
config_file >> config;
@@ -40,19 +41,23 @@ void demp_analysis(const std::string& config_name)
fmt::print(" - output prefix for histograms: {}\n", output_prefix);
fmt::print(" - ebeam: {}\n", ebeam);
fmt::print(" - pbeam: {}\n", pbeam);
-
- //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------
+
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
// Set output file for the histograms
std::string output_name_hists = fmt::format("{}.root", output_prefix);
cout << "Output file for histograms = " << output_name_hists << endl;
TFile* ofile = new TFile(output_name_hists.c_str(), "RECREATE");
- //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
// Set up input file chain
TChain *mychain = new TChain("events");
mychain->Add(rec_file.c_str());
- //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------
+
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
// Initialize reader
TTreeReader tree_reader(mychain);
@@ -88,8 +93,9 @@ void demp_analysis(const std::string& config_name)
TTreeReaderArray<float> neutPosY_hcal(tree_reader, "HcalFarForwardZDCClusters.position.y");
TTreeReaderArray<float> neutPosZ_hcal(tree_reader, "HcalFarForwardZDCClusters.position.z");
TTreeReaderArray<float> neutEng_hcal(tree_reader, "HcalFarForwardZDCClusters.energy");
-
- //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------
+
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
// Define Histograms
TH2* eTruthw_Thetap = new TH2D("eTruthw_Thetap","e' truth #theta vs P; #theta (deg); P (GeV/c); Rate/bin (Hz)",100,120,165,100,4,6);
@@ -108,6 +114,7 @@ void demp_analysis(const std::string& config_name)
TH2* nRecw_Thetap = new TH2D("nRecw_Thetap","n rec #theta vs P; #theta (Deg); P (GeV/c); Rate/bin (Hz)",100,0.8,2.0,100,0,60);
TH2* nRecw_Thetaphi = new TH2D("nRecw_Thetaphi","n rec #theta vs #phi; #theta (mRad); #phi (deg); Rate/bin (Hz)",100,15.0,35.0,100,-200,200);
TH2* nRecw_rot_Thetaphi = new TH2D("nRecw_rot_Thetaphi","n rec #theta* vs #phi* around p axis; #theta* (mRad); #phi* (deg); Rate/bin (Hz)",100,0.0,10.0,100,-200,200);
+ TH2* nRecw_rot_PosXY = new TH2D("nRecw_rot_PosXY","n X vs Y around proton axis at Z = 35 m ( #theta* < 4.0 mRad, E > 10 GeV ); x (mm); y (mm); Rate/bin (Hz)",100,-200,200,100,-200,200);
TH2* nRecw_rot_Thetap = new TH2D("nRecw_rot_Thetap","n rec #theta* vs P around p axis ( #theta* < 4.0 mRad, E > 10 GeV ); #theta* (mRad); P (GeV/c); Rate/bin (Hz)",100,0.0,4.0,100,5,60);
TH1* nRec_en = new TH1D("nRec_en", "n rec E for all clusters ( #theta* < 4.0 mRad ); E (GeV); Rate (Hz)", 100, 0.0, 60);
@@ -120,10 +127,10 @@ void demp_analysis(const std::string& config_name)
n_ThetaDiff->SetLineWidth(2);
TH1* n_PhiDiff = new TH1D("n_PhiDiff", " #phi*_{pMiss_rec} - #phi*_{ZDC}; #phi*_{pMiss_rec} - #phi*_{ZDC}(Deg); Rate (Hz)", 100, -50, 50);
n_PhiDiff->SetLineWidth(2);
- TH2* n_ThetaPhiDiff = new TH2D("n_ThetaPhiDiff", "#theta*_{pMiss_rec} - #theta*_{ZDC} vs #phi*_{pMiss_rec} - #phi*_{ZDC}; #theta*_{pMiss_rec} - #theta*_{ZDC} (Deg); #phi*_{pMiss_rec} - #phi*_{ZDC} (Deg); Rate/bin (Hz)",100, -0.3, 1.5, 100, -50, 50);
+ TH2* n_ThetaPhiDiff = new TH2D("n_ThetaPhiDiff", "#theta*_{pMiss_rec} - #theta*_{ZDC} vs #phi*_{pMiss_rec} - #phi*_{ZDC}; #theta*_{pMiss_rec} - #theta*_{ZDC} (Deg); #phi*_{pMiss_rec} - #phi*_{ZDC} (Deg); Rate/bin (Hz)",100, -1.0, 1.0, 100, -75, 75);
TH2* pMissRecw_Thetaphi = new TH2D("pMissRecw_Thetaphi", "pMiss rec #theta vs #phi; #theta (mRad); #phi (deg); Rate/bin (Hz)",100,15.0,35.0,100,-200,200);
TH2* pMissRecw_rot_Thetaphi = new TH2D("pMissRecw_rot_Thetaphi", "pMiss rec #theta* vs #phi* around p axis; #theta* (mRad); #phi* (deg); Rate/bin (Hz)",100,0.0,10.0,100,-200,200);
- TH2* n_TruthRecw_ThetaPhiDiff = new TH2D("n_TruthRecw_ThetaPhiDiff", " #theta*_{n_MC} - #theta*_{n_rec} vs #phi*_{n_MC} - #phi*_{n_rec}; #theta*_{n_MC} - #theta*_{n_rec} (Deg); #phi*_{n_MC} - #phi*_{n_rec} (Deg); Rate/bin (Hz)",100, -0.2, 0.2, 100, -20, 20);
+ TH2* n_TruthRecw_ThetaPhiDiff = new TH2D("n_TruthRecw_ThetaPhiDiff", " #theta*_{n_MC} - #theta*_{n_rec} vs #phi*_{n_MC} - #phi*_{n_rec}; #theta*_{n_MC} - #theta*_{n_rec} (Deg); #phi*_{n_MC} - #phi*_{n_rec} (Deg); Rate/bin (Hz)",100, -0.2, 0.2, 100, -25, 25);
// Absolute difference -t plots
TH1* htw_t1 = new TH1D("htw_t1", "-t_{rec, alt_rec, rec_pT, rec_corr} - -t_{truth} Distribution; #Delta -t (GeV^{2}); Rate (Hz) ", 100, -0.1,0.1);
@@ -216,8 +223,9 @@ void demp_analysis(const std::string& config_name)
// Neutrons HCal
TH2* nRecw_Thetap_hcal = new TH2D("nRecw_Thetap_hcal","n rec #theta vs P for 1 cluster events; #theta (Deg); P (GeV/c); Rate/bin (Hz)",100,0.8,2.0,100,0,40);
TH2* nRecw_rot_Thetap_hcal = new TH2D("nRecw_rot_Thetap_hcal","n rec #theta* vs P around p axis for 1 cluster events; #theta* (mRad); P (GeV/c); Rate/bin (Hz)",100,0,4.0,100,5,40);
-
- //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------
+
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
//Defining the four vectors
ROOT::Math::PxPyPzEVector elec_beam; // initialized the 4 vector for electron beam
ROOT::Math::PxPyPzEVector prot_beam; // initialized the 4 vector for proton beam
@@ -261,8 +269,10 @@ void demp_analysis(const std::string& config_name)
double weight, partEng; // weight and energy of the particles
double Q2_truth, W_truth, y_truth, t_truth, t_alttruth; // Truth kinematic variables
double Q2_rec, W_rec, y_rec, t_rec, t_altrec, t_recpT, t_reccorr; // Reconstructed kinematic variables
+ double neutPosX, neutPosY; // neutron position
- //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
// Defining initial colliding beams
double eMass = 0.000510998950; //electron beam
double eEng = ebeam;
@@ -281,7 +291,8 @@ void demp_analysis(const std::string& config_name)
double p_p3 = p_pmag*cos(c_a);
prot_beam.SetPxPyPzE(p_p1, p_p2, p_p3, pEng);
- //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
bool x,y,z; // x,y, and z are for reconstructed electron, pion, and neutron
while(tree_reader.Next()) { // Loop over events
@@ -292,8 +303,9 @@ void demp_analysis(const std::string& config_name)
vector<string> weight_value = weight_map[0].second;
weight = std::stod( *(weight_value.begin()) );
- //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------
-
+
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
for(unsigned int i=0; i<partGenStat.GetSize(); i++) { // Loop over thrown particles
partEng = sqrt(pow(partMomX[i],2) + pow(partMomY[i],2) + pow(partMomZ[i],2) + pow(partMass[i],2)); // Energy of all Monte Carlo particles
@@ -323,8 +335,9 @@ void demp_analysis(const std::string& config_name)
} // for over thrown particles
- //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------
-
+
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
for(unsigned int i=0; i<trackPdg.GetSize(); i++) { // Loop over reconstructed particles
// if(trackPdg[i] == 11) { // Look at electron
if(trackCharge[i] == -1 && trackMomZ[i] < 0) {
@@ -348,8 +361,8 @@ void demp_analysis(const std::string& config_name)
}// for over reconstructed particles
- //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------
-
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
for(unsigned int i=0; i<neutEng.GetSize(); i++) { // Loop over zdc neutrons
neut_rec.SetPxPyPzE(neutMomX[i],neutMomY[i],neutMomZ[i], neutEng[i]);
@@ -366,14 +379,18 @@ void demp_analysis(const std::string& config_name)
if(neut_rot_rec.E()>10.0){ // neutron energy cut
z = true;
nRecw_rot_Thetap -> Fill(neut_rot_rec.Theta()*1000., neut_rot_rec.P(), weight);
+
+ neutPosX = 35000 * sin(neut_rot_rec.Theta()) * cos(neut_rot_rec.Phi()); // neutron position at r = z = 35.0 m
+ neutPosY = 35000 * sin(neut_rot_rec.Theta()) * sin(neut_rot_rec.Phi());
+ nRecw_rot_PosXY -> Fill(neutPosX, neutPosY, weight);
}
}
}// for over zdc neutrons
- //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------
-
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
for(unsigned int i=0; i<neutEng.GetSize(); i++) { // Loop over zdc neutrons in HCal
hcal_clus_size = neutEng_hcal.GetSize(); //ZDC HCal cluster size -> No. of clusters in ZDC
@@ -399,9 +416,9 @@ void demp_analysis(const std::string& config_name)
}// for over zdc neutrons in HCal
- //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
// Truth kinematic variables
-
virtphoton_truth = (elec_beam - elec_mc);
Q2_truth = -1*(virtphoton_truth.mag2());
W_truth = ((virtphoton_truth + prot_beam).mag());
@@ -421,9 +438,9 @@ void demp_analysis(const std::string& config_name)
}
}
- //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
// Reconstructed kinematic variables
-
if (x == true && y == true && z == true ){ // if e', pi, and neutron are in coincidence
virtphoton_rec = (elec_beam - elec_rec);
@@ -509,7 +526,8 @@ void demp_analysis(const std::string& config_name)
cout<<"truth_neutron = "<<count2<<endl;
- //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
ofile->Write(); // Write histograms to file
ofile->Close(); // Close output file
diff --git a/benchmarks/demp/analysis/demp_plots.cxx b/benchmarks/demp/analysis/demp_plots.cxx
index ba7317ecef5778da72f26c3186d65e7b459e4a73..e9a7d3f565584f596413dce76da4859b3fc5fdfe 100644
--- a/benchmarks/demp/analysis/demp_plots.cxx
+++ b/benchmarks/demp/analysis/demp_plots.cxx
@@ -16,11 +16,70 @@
#include "nlohmann/json.hpp"
+//--------------------------------------------------------------------------------------------------------------------------------------------
+
+//Plotting style for histograms
+void plots_2D(TH2D* hist, const int& nfiles = 1, const double font = 42, const double tsize = 0.04){
+
+ hist->SetStats(1);
+
+ hist->GetXaxis()->SetTitleFont(font); //x-axis
+ hist->GetXaxis()->SetTitleSize(tsize);
+ hist->GetXaxis()->SetTitleOffset(1.3);
+ hist->GetXaxis()->SetLabelFont(font);
+ hist->GetXaxis()->SetLabelSize(tsize);
+ hist->GetXaxis()->SetLabelOffset(0.01);
+
+ hist->GetYaxis()->SetTitleFont(font); //y-axis
+ hist->GetYaxis()->SetTitleSize(tsize);
+ hist->GetYaxis()->SetTitleOffset(1.3);
+ hist->GetYaxis()->SetLabelFont(font);
+ hist->GetYaxis()->SetLabelSize(tsize);
+ hist->GetYaxis()->SetLabelOffset(0.01);
+
+ hist->GetZaxis()->SetTitleFont(font); //z-axis
+ hist->GetZaxis()->SetTitleSize(0.036);
+ hist->GetZaxis()->SetTitleOffset(1.05);
+ hist->GetZaxis()->SetLabelFont(font);
+ hist->GetZaxis()->SetLabelSize(0.03);
+ hist->GetZaxis()->SetLabelOffset(0.01);
+
+ hist->Scale(1.0/nfiles);
+ hist->Draw("colz");
+
+}
+
+void plots_1D(TH1D* hist, const int& nfiles = 1, const double font = 42, const double tsize = 0.04){
+
+ hist->SetStats(1);
+
+ hist->GetXaxis()->SetTitleFont(font); //x-axis
+ hist->GetXaxis()->SetTitleSize(tsize);
+ hist->GetXaxis()->SetTitleOffset(1.3);
+ hist->GetXaxis()->SetLabelFont(font);
+ hist->GetXaxis()->SetLabelSize(tsize);
+ hist->GetXaxis()->SetLabelOffset(0.01);
+
+ hist->GetYaxis()->SetTitleFont(font); //y-axis
+ hist->GetYaxis()->SetTitleSize(tsize);
+ hist->GetYaxis()->SetTitleOffset(1.3);
+ hist->GetYaxis()->SetLabelFont(font);
+ hist->GetYaxis()->SetLabelSize(tsize);
+ hist->GetYaxis()->SetLabelOffset(0.01);
+
+ hist->Scale(1.0/nfiles);
+ hist->Draw("HIST");
+
+}
+
+//--------------------------------------------------------------------------------------------------------------------------------------------
+
void demp_plots(const std::string& config_name)
{
- //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------
- // read our configuration
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
+ // Read our configuration
std::ifstream config_file{config_name};
nlohmann::json config;
config_file >> config;
@@ -30,6 +89,7 @@ void demp_plots(const std::string& config_name)
const std::string output_prefix = config["output_prefix"];
const int ebeam = config["ebeam"];
const int pbeam = config["pbeam"];
+ const int nfiles = config["nfiles"];
fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green),
"Running DEMP analysis...\n");
@@ -38,8 +98,10 @@ void demp_plots(const std::string& config_name)
fmt::print(" - output prefix for plots: {}\n", output_prefix);
fmt::print(" - ebeam: {}\n", ebeam);
fmt::print(" - pbeam: {}\n", pbeam);
-
- //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------
+ fmt::print(" - nfiles: {}\n", nfiles);
+
+ //-------------------------------------------------------------------------------------------------------------------------------------------
+
// Read file with histograms
TFile* file = new TFile(hists_file.c_str());
@@ -57,6 +119,7 @@ void demp_plots(const std::string& config_name)
TH2D* piRecw_Thetaphi = (TH2D*) file->Get("piRecw_Thetaphi");
TH2D* nRecw_Thetaphi = (TH2D*) file->Get("nRecw_Thetaphi");
TH2D* nRecw_rot_Thetaphi = (TH2D*) file->Get("nRecw_rot_Thetaphi");
+ TH2D* nRecw_rot_PosXY = (TH2D*) file->Get("nRecw_rot_PosXY");
TH2D* nRecw_Thetap_hcal = (TH2D*) file->Get("nRecw_Thetap_hcal");
TH2D* nRecw_rot_Thetap_hcal = (TH2D*) file->Get("nRecw_rot_Thetap_hcal");
TH1D* nRec_en = (TH1D*) file->Get("nRec_en");
@@ -109,179 +172,194 @@ void demp_plots(const std::string& config_name)
TH1D* piTruthw_P_Uncut = (TH1D*) file->Get("piTruthw_P_Uncut");
TH1D* piRecw_P_Cut = (TH1D*) file->Get("piRecw_P_Cut");
- //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
// Make plots and save to PDF file
std::cout<<"Making plots..."<<std::endl;
-
- gStyle->SetPadRightMargin(0.125); // left space on right side
- gStyle->SetStatX(0.895); // move the stat bax on left or right side
- gStyle->SetStatY(0.90); // move the stat bax on up or down side
+
+ gStyle->SetPadRightMargin(0.125);
+ gStyle->SetPadTopMargin(0.1);
+ gStyle->SetOptTitle(1); //setting title
+ gStyle->SetTitleAlign(1);
+ gStyle->SetTitleX(0.35);
+ gStyle->SetTitleY(0.92);
+ gStyle->SetTitleW(0.32);
+ gStyle->SetTitleH(0.065);
+ gStyle->SetOptStat(1); //setting stat box
+ gStyle->SetStatX(0.875);
+ gStyle->SetStatY(0.9);
+ gStyle->SetStatW(0.14);
+ gStyle->SetStatH(0.11);
gStyle->SetPalette(55);
- gStyle->SetLineStyleString(2,"[12 12]");
gStyle->SetHistLineWidth(2);
TCanvas *c1 = new TCanvas("c1"); //truth information
- c1->SetLogz();
- eTruthw_Thetap->Draw("colz");
+ c1->SetLogz();
+ plots_2D(eTruthw_Thetap, nfiles);
TCanvas *c2 = new TCanvas("c2");
c2->SetLogz();
- piTruthw_Thetap->Draw("colz");
-
+ plots_2D(piTruthw_Thetap, nfiles);
+
TCanvas *c3 = new TCanvas("c3");
c3->SetLogz();
- nTruthw_Thetaphi->Draw("colz");
-
+ plots_2D(nTruthw_Thetaphi, nfiles);
+
TCanvas *c3a = new TCanvas("c3a");
c3a->SetLogz();
- nTruthw_rot_Thetaphi->Draw("colz");
+ plots_2D(nTruthw_rot_Thetaphi, nfiles);
TCanvas *c3b = new TCanvas("c3b");
c3b->SetLogz();
- nTruthw_Thetap->Draw("colz");
+ plots_2D(nTruthw_Thetap, nfiles);
TCanvas *c3c = new TCanvas("c3c");
c3c->SetLogz();
- nTruthw_rot_Thetap->Draw("colz");
+ plots_2D(nTruthw_rot_Thetap, nfiles);
TCanvas *c4 = new TCanvas("c4"); //reconstructed information
c4->SetLogz();
- eRecw_Thetap->Draw("colz");
+ plots_2D(eRecw_Thetap, nfiles);
TCanvas *c4a = new TCanvas("c4a");
c4a->SetLogz();
- eRecw_Thetaphi->Draw("colz");
+ plots_2D(eRecw_Thetaphi, nfiles);
TCanvas *c5 = new TCanvas("c5");
c5->SetLogz();
- piRecw_Thetap->Draw("colz");
+ plots_2D(piRecw_Thetap, nfiles);
TCanvas *c5_a = new TCanvas("c5_a");
c5_a->SetLogz();
- piRecw_Thetaphi->Draw("colz");
+ plots_2D(piRecw_Thetap, nfiles);
TCanvas *c5a = new TCanvas("c5a");
c5a->SetLogz();
- nRecw_Thetaphi->Draw("colz");
+ plots_2D(nRecw_Thetaphi, nfiles);
TCanvas *c5b = new TCanvas("c5b");
c5b->SetLogz();
- nRecw_rot_Thetaphi->Draw("colz");
+ plots_2D(nRecw_rot_Thetaphi, nfiles);
TCanvas *c5c = new TCanvas("c5c");
c5c->SetLogz();
- nRecw_Thetap_hcal->Draw("colz");
+ plots_2D(nRecw_Thetap_hcal, nfiles);
TCanvas *c5d = new TCanvas("c5d");
c5d->SetLogz();
- nRecw_rot_Thetap_hcal->Draw("colz");
+ plots_2D(nRecw_rot_Thetap_hcal, nfiles);
TCanvas *c6a = new TCanvas("c6a");
- nRec_en->Draw("HIST");
+ plots_1D(nRec_en, nfiles);
TCanvas *c6b = new TCanvas("c6b");
- nRec_clus->Draw();
+ plots_1D(nRec_clus);
TCanvas *c7 = new TCanvas("c7");
c7->SetLogz();
- nRecw_Thetap->Draw("colz");
+ plots_2D(nRecw_Thetap, nfiles);
TCanvas *c7a = new TCanvas("c7a");
c7a->SetLogz();
- nRecw_rot_Thetap->Draw("colz");
-
+ plots_2D(nRecw_rot_Thetap, nfiles);
+
+ TCanvas *c7b = new TCanvas("c7b");
+ c7b->SetLogz();
+ plots_2D(nRecw_rot_PosXY, nfiles);
+
TCanvas *c8 = new TCanvas("c8"); // -t reconstruction plots
c8->SetLogz();
- htw_rec1->Draw("colz");
+ plots_2D(htw_rec1, nfiles);
TCanvas *c8a = new TCanvas("c8a");
c8a->SetLogz();
- htwz_rec1->Draw("colz");
+ plots_2D(htwz_rec1, nfiles);
TCanvas *c9 = new TCanvas("c9");
c9->SetLogz();
- htw_rec2->Draw("colz");
+ plots_2D(htw_rec2, nfiles);
TCanvas *c9a = new TCanvas("c9a");
c9a->SetLogz();
- htwz_rec2->Draw("colz");
+ plots_2D(htwz_rec2, nfiles);
TCanvas *c10 = new TCanvas("c10");
c10->SetLogz();
- htw_rec3->Draw("colz");
+ plots_2D(htw_rec3, nfiles);
TCanvas *c10a = new TCanvas("c10a");
c10a->SetLogz();
- htwz_rec3->Draw("colz");
+ plots_2D(htwz_rec3, nfiles);
TCanvas *c11 = new TCanvas("c11");
c11->SetLogz();
- htw_rec4->Draw("colz");
+ plots_2D(htw_rec4, nfiles);
TCanvas *c11a = new TCanvas("c11a");
c11a->SetLogz();
- htwz_rec4->Draw("colz");
+ plots_2D(htwz_rec4, nfiles);
TCanvas *c12 = new TCanvas("c12"); // Resolution plots
- htw_res_e->Draw("HIST");
+ plots_1D(htw_res_e, nfiles);
TCanvas *c13 = new TCanvas("c13");
- htw_res_pi->Draw("HIST");
+ plots_1D(htw_res_pi, nfiles);
TCanvas *c14a = new TCanvas("c14a");
- htw_res_n1->Draw("HIST");
+ plots_1D(htw_res_n1, nfiles);
TCanvas *c14b = new TCanvas("c14b");
- htw_res_n2->Draw("HIST");
+ plots_1D(htw_res_n2, nfiles);
TCanvas *c14c = new TCanvas("c14c");
- htw_res_n3->Draw("HIST");
+ plots_1D(htw_res_n3, nfiles);
TCanvas *c14d = new TCanvas("c14d");
- htw_res_n4->Draw("HIST");
+ plots_1D(htw_res_n4, nfiles);
TCanvas *c15a = new TCanvas("c15a"); //t-resoutions
- htw_res1->Draw("HIST");
+ plots_1D(htw_res1, nfiles);
TCanvas *c15b = new TCanvas("c15b");
- htw_res2->Draw("HIST");
+ plots_1D(htw_res2, nfiles);
TCanvas *c15c = new TCanvas("c15c");
- htw_res3->Draw("HIST");
+ plots_1D(htw_res3, nfiles);
TCanvas *c15d = new TCanvas("c15d");
- htw_res4->Draw("HIST");
+ plots_1D(htw_res4, nfiles);
TCanvas *c15e = new TCanvas("c15e");
- htw_res5->Draw("HIST");
+ plots_1D(htw_res5, nfiles);
TCanvas *c15f = new TCanvas("c15f");
- htw_res6->Draw("HIST");
+ plots_1D(htw_res6, nfiles);
TCanvas *c16a = new TCanvas("c16a"); // Neutron theta-phi plots
- n_ThetaDiff->Draw("HIST");
+ plots_1D(n_ThetaDiff, nfiles);
TCanvas *c16b = new TCanvas("c16b");
- n_PhiDiff->Draw("HIST");
+ plots_1D(n_PhiDiff, nfiles);
TCanvas *c16c = new TCanvas("c16c");
c16c->SetLogz();
- n_ThetaPhiDiff->Draw("colz");
+ plots_2D(n_ThetaPhiDiff, nfiles);
TCanvas *c16d = new TCanvas("c16d");
c16d->SetLogz();
- pMissRecw_Thetaphi->Draw("colz");
+ plots_2D(pMissRecw_Thetaphi, nfiles);
TCanvas *c16e = new TCanvas("c16e");
c16e->SetLogz();
- pMissRecw_rot_Thetaphi->Draw("colz");
+ plots_2D(pMissRecw_rot_Thetaphi, nfiles);
TCanvas *c16f = new TCanvas("c16f");
c16f->SetLogz();
- n_TruthRecw_ThetaPhiDiff->Draw("colz");
+ plots_2D(n_TruthRecw_ThetaPhiDiff, nfiles);
TCanvas *c17 = new TCanvas("c17"); // Absolute difference -t plots
- htw_t4->Draw("HIST");
+ htw_t1->Scale(1.0/nfiles), htw_t2->Scale(1.0/nfiles), htw_t3->Scale(1.0/nfiles), htw_t4->Scale(1.0/nfiles); // special case
+ plots_1D(htw_t4);
htw_t3->Draw("HIST SAME");
htw_t2->Draw("HIST SAME");
htw_t1->Draw("HIST SAME");
@@ -296,51 +374,52 @@ void demp_plots(const std::string& config_name)
TCanvas *c18 = new TCanvas("c18"); // Efficiency plots
Q2_t_DetEff-> GetZaxis()->SetRangeUser(0.0,1.0);
- Q2_t_DetEff->Draw("colz");
+ plots_2D(Q2_t_DetEff);
TCanvas *c18a = new TCanvas("c18a");
- Q2_t_DetEff_Cut->Draw("colz");
+ plots_2D(Q2_t_DetEff_Cut, nfiles);
TCanvas *c18b = new TCanvas("c18b");
- Q2_t_DetEff_Uncut->Draw("colz");
+ plots_2D(Q2_t_DetEff_Uncut, nfiles);
TCanvas *c19 = new TCanvas("c19"); // elec eta eff.
- eEff_Eta->Draw("HIST");
+ plots_1D(eEff_Eta);
TCanvas *c19a = new TCanvas("c19a");
- eTruthw_Eta_Uncut->Draw("HIST");
+ plots_1D(eTruthw_Eta_Uncut, nfiles);
TCanvas *c19b = new TCanvas("c19b");
- eRecw_Eta_Cut->Draw("HIST");
+ plots_1D(eRecw_Eta_Cut, nfiles);
TCanvas *c20 = new TCanvas("c20"); // elec mom eff.
- eEff_P->Draw("HIST");
+ plots_1D(eEff_P);
TCanvas *c20a = new TCanvas("c20a");
- eTruthw_P_Uncut->Draw("HIST");
+ plots_1D(eTruthw_P_Uncut, nfiles);
TCanvas *c20b = new TCanvas("c20b");
- eRecw_P_Cut->Draw("HIST");
+ plots_1D(eRecw_P_Cut, nfiles);
TCanvas *c21 = new TCanvas("c21"); // pi eta eff.
- piEff_Eta->Draw("HIST");
+ plots_1D(piEff_Eta);
TCanvas *c21a = new TCanvas("c21a");
- piTruthw_Eta_Uncut->Draw("HIST");
+ plots_1D(piTruthw_Eta_Uncut, nfiles);
TCanvas *c21b = new TCanvas("c21b");
- piRecw_Eta_Cut->Draw("HIST");
+ plots_1D(piRecw_Eta_Cut, nfiles);
TCanvas *c22 = new TCanvas("c22"); //pi mom eff.
- piEff_P->Draw("HIST");
+ plots_1D(piEff_P);
TCanvas *c22a = new TCanvas("c22a");
- piTruthw_P_Uncut->Draw("HIST");
+ plots_1D(piTruthw_P_Uncut, nfiles);
TCanvas *c22b = new TCanvas("c22b");
- piRecw_P_Cut->Draw("HIST");
+ plots_1D(piRecw_P_Cut, nfiles);
+
+ //--------------------------------------------------------------------------------------------------------------------------------------------
- //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------
// Print plots to pdf file
c1->Print(fmt::format("{}.pdf[", output_prefix).c_str());
c1->Print(fmt::format("{}.pdf", output_prefix).c_str());
@@ -361,6 +440,7 @@ void demp_plots(const std::string& config_name)
c6a->Print(fmt::format("{}.pdf", output_prefix).c_str());
c7->Print(fmt::format("{}.pdf", output_prefix).c_str());
c7a->Print(fmt::format("{}.pdf", output_prefix).c_str());
+ c7b->Print(fmt::format("{}.pdf", output_prefix).c_str());
c8->Print(fmt::format("{}.pdf", output_prefix).c_str());
c8a->Print(fmt::format("{}.pdf", output_prefix).c_str());
c9->Print(fmt::format("{}.pdf", output_prefix).c_str());