Skip to content
Snippets Groups Projects
Commit 78e69c26 authored by Marshall Scott's avatar Marshall Scott
Browse files

rebasing

parent 747958ea
No related branches found
No related tags found
1 merge request!43Resolve "Add detector name to results folder for benchmarks"
...@@ -3,6 +3,10 @@ ...@@ -3,6 +3,10 @@
// M. Scott 05/2021 // M. Scott 05/2021
////////////////////////// //////////////////////////
#include <cctype>
#include "TH1.h"
// Returns particle pdgID and mass in [GeV]
std::tuple <int, double> extract_particle_parameters(std::string particle_name) { std::tuple <int, double> extract_particle_parameters(std::string particle_name) {
if (particle_name == "electron") return std::make_tuple(11, 0.51099895e-3); if (particle_name == "electron") return std::make_tuple(11, 0.51099895e-3);
if (particle_name == "photon") return std::make_tuple(22, 0.0); if (particle_name == "photon") return std::make_tuple(22, 0.0);
...@@ -15,4 +19,18 @@ std::tuple <int, double> extract_particle_parameters(std::string particle_name) ...@@ -15,4 +19,18 @@ std::tuple <int, double> extract_particle_parameters(std::string particle_name)
std::cout << "wrong particle name" << std::endl; std::cout << "wrong particle name" << std::endl;
abort(); abort();
}
// Returns Environment Variables
std::string getEnvVar( std::string const & key ){
char * val = getenv( key.c_str() );
return val == NULL ? std::string("") : std::string(val);
}
// Added detetcor name to title
void addDetectorName(std::string name, TH1 * inhist){
std::string newName = inhist -> GetTitle();
for (auto& x : name){
x = toupper(x);
}
inhist -> SetTitle((name + " " + newName).c_str());
} }
\ No newline at end of file
...@@ -16,11 +16,15 @@ ...@@ -16,11 +16,15 @@
#include "TF1.h" #include "TF1.h"
#include "TH1D.h" #include "TH1D.h"
#include "emcal_barrel_common_functions.h"
#include <cctype>
using ROOT::RDataFrame; using ROOT::RDataFrame;
using namespace ROOT::VecOps; using namespace ROOT::VecOps;
void emcal_barrel_electrons_analysis(const char* input_fname = "sim_output/sim_emcal_barrel_uniform_electrons.root") void emcal_barrel_electrons_analysis(const char* input_fname = "sim_output/sim_emcal_barrel_uniform_electrons.root")
{ {
input_fname = "temp_pions_electrons.root";
// Setting for graphs // Setting for graphs
gROOT->SetStyle("Plain"); gROOT->SetStyle("Plain");
gStyle->SetOptFit(1); gStyle->SetOptFit(1);
...@@ -35,6 +39,16 @@ void emcal_barrel_electrons_analysis(const char* input_fname = "sim_output/sim_e ...@@ -35,6 +39,16 @@ void emcal_barrel_electrons_analysis(const char* input_fname = "sim_output/sim_e
ROOT::EnableImplicitMT(); ROOT::EnableImplicitMT();
ROOT::RDataFrame d0("events", input_fname); ROOT::RDataFrame d0("events", input_fname);
// Environment Variables
std::string detector_path = "";
std::string detector_name = "topside";
if(std::getenv("DETECTOR_PATH")) {
detector_path = std::getenv("DETECTOR_PATH");
}
if(std::getenv("JUGGLER_DETECTOR")) {
detector_name = std::getenv("JUGGLER_DETECTOR");
}
// Thrown Energy [GeV] // Thrown Energy [GeV]
auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) { auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
auto p = input[2]; auto p = input[2];
...@@ -79,6 +93,10 @@ void emcal_barrel_electrons_analysis(const char* input_fname = "sim_output/sim_e ...@@ -79,6 +93,10 @@ void emcal_barrel_electrons_analysis(const char* input_fname = "sim_output/sim_e
{"hfsam", "Sampling Fraction; Sampling Fraction; Events", 100, 0.0, 0.1}, {"hfsam", "Sampling Fraction; Sampling Fraction; Events", 100, 0.0, 0.1},
"fsam"); "fsam");
addDetectorName(detector_name, hEthr.GetPtr());
addDetectorName(detector_name, hEsim.GetPtr());
addDetectorName(detector_name, hfsam.GetPtr());
// Event Counts // Event Counts
auto nevents_thrown = d1.Count(); auto nevents_thrown = d1.Count();
std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n"; std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
......
...@@ -51,6 +51,16 @@ void emcal_barrel_particles_analysis(std::string particle_name = "electron") ...@@ -51,6 +51,16 @@ void emcal_barrel_particles_analysis(std::string particle_name = "electron")
std::string input_fname = fmt::format("sim_output/sim_emcal_barrel_{}.root", particle_name); std::string input_fname = fmt::format("sim_output/sim_emcal_barrel_{}.root", particle_name);
ROOT::RDataFrame d0("events", input_fname); ROOT::RDataFrame d0("events", input_fname);
// Environment Variables
std::string detector_path = "";
std::string detector_name = "topside";
if(std::getenv("DETECTOR_PATH")) {
detector_path = std::getenv("DETECTOR_PATH");
}
if(std::getenv("JUGGLER_DETECTOR")) {
detector_name = std::getenv("JUGGLER_DETECTOR");
}
// Thrown Energy [GeV] // Thrown Energy [GeV]
auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) { auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
auto p = input[2]; auto p = input[2];
...@@ -94,6 +104,10 @@ void emcal_barrel_particles_analysis(std::string particle_name = "electron") ...@@ -94,6 +104,10 @@ void emcal_barrel_particles_analysis(std::string particle_name = "electron")
{"hfsam", "Sampling Fraction; Sampling Fraction; Events", 200, 0.0, 0.05}, {"hfsam", "Sampling Fraction; Sampling Fraction; Events", 200, 0.0, 0.05},
"fsam"); "fsam");
addDetectorName(detector_name, hEthr.GetPtr());
addDetectorName(detector_name, hEsim.GetPtr());
addDetectorName(detector_name, hfsam.GetPtr());
// Event Counts // Event Counts
auto nevents_thrown = d1.Count(); auto nevents_thrown = d1.Count();
std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n"; std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
......
...@@ -24,7 +24,7 @@ R__LOAD_LIBRARY(libfmt.so) ...@@ -24,7 +24,7 @@ R__LOAD_LIBRARY(libfmt.so)
#include "DD4hep/Detector.h" #include "DD4hep/Detector.h"
#include "DDG4/Geant4Data.h" #include "DDG4/Geant4Data.h"
#include "DDRec/CellIDPositionConverter.h" #include "DDRec/CellIDPositionConverter.h"
#include "emcal_barrel_common_functions.h"
using ROOT::RDataFrame; using ROOT::RDataFrame;
using namespace ROOT::VecOps; using namespace ROOT::VecOps;
...@@ -45,7 +45,7 @@ void emcal_barrel_pions_electrons_analysis(const char* input_fname1 = "sim_outpu ...@@ -45,7 +45,7 @@ void emcal_barrel_pions_electrons_analysis(const char* input_fname1 = "sim_outpu
ROOT::EnableImplicitMT(); ROOT::EnableImplicitMT();
ROOT::RDataFrame d0("events", {input_fname1, input_fname2}); ROOT::RDataFrame d0("events", {input_fname1, input_fname2});
//Using the detector layers // Environment Variables
std::string detector_path = ""; std::string detector_path = "";
std::string detector_name = "topside"; std::string detector_name = "topside";
if(std::getenv("DETECTOR_PATH")) { if(std::getenv("DETECTOR_PATH")) {
...@@ -55,6 +55,7 @@ void emcal_barrel_pions_electrons_analysis(const char* input_fname1 = "sim_outpu ...@@ -55,6 +55,7 @@ void emcal_barrel_pions_electrons_analysis(const char* input_fname1 = "sim_outpu
detector_name = std::getenv("JUGGLER_DETECTOR"); detector_name = std::getenv("JUGGLER_DETECTOR");
} }
// Using the detector layers
dd4hep::Detector& detector = dd4hep::Detector::getInstance(); dd4hep::Detector& detector = dd4hep::Detector::getInstance();
detector.fromCompact(fmt::format("{}/{}.xml", detector_path,detector_name)); detector.fromCompact(fmt::format("{}/{}.xml", detector_path,detector_name));
//dd4hep::rec::CellIDPositionConverter cellid_converter(detector); //dd4hep::rec::CellIDPositionConverter cellid_converter(detector);
...@@ -156,48 +157,66 @@ void emcal_barrel_pions_electrons_analysis(const char* input_fname1 = "sim_outpu ...@@ -156,48 +157,66 @@ void emcal_barrel_pions_electrons_analysis(const char* input_fname1 = "sim_outpu
auto d_pim_rej = d_pim.Filter("Esim_front > " + to_string(rejectCut)); auto d_pim_rej = d_pim.Filter("Esim_front > " + to_string(rejectCut));
// Define Histograms // Define Histograms
auto hEthr = d1.Histo1D({"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 7.5}, "Ethr"); auto hEthr = d1.Histo1D({"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 7.5}, "Ethr");
auto hNhits = d1.Histo1D({"hNhits", "Number of hits per events; Number of hits; Events", 100, 0.0, 2000.0}, "nhits"); auto hNhits = d1.Histo1D({"hNhits", "Number of hits per events; Number of hits; Events", 100, 0.0, 2000.0}, "nhits");
auto hEsim = d1.Histo1D({"hEsim", "Energy Deposit; Energy Deposit [GeV]; Events", 10, 0.0, 0.05}, "Esim"); auto hEsim = d1.Histo1D({"hEsim", "Energy Deposit; Energy Deposit [GeV]; Events", 10, 0.0, 0.05}, "Esim");
auto hfsam = d1.Histo1D({"hfsam", "Sampling Fraction; Sampling Fraction; Events", 100, 0.0, 0.1}, "fsam"); auto hfsam = d1.Histo1D({"hfsam", "Sampling Fraction; Sampling Fraction; Events", 100, 0.0, 0.1}, "fsam");
auto hEsim_front = d1.Histo1D({"hEsim_front", "Energy Deposit Front; Energy Deposit [GeV]; Events", 10, 0.0, 0.05}, "Esim_front"); auto hEsim_front = d1.Histo1D({"hEsim_front", "Energy Deposit Front; Energy Deposit [GeV]; Events", 10, 0.0, 0.05}, "Esim_front");
addDetectorName(detector_name, hEthr.GetPtr());
addDetectorName(detector_name, hNhits.GetPtr());
addDetectorName(detector_name, hfsam.GetPtr());
addDetectorName(detector_name, hEsim_front.GetPtr());
auto hEsim_ele = d_ele.Histo1D({"hEsim_ele", "Energy Deposit Electron; Energy Deposit [GeV]; Events", 10, 0.0, 0.05}, "Esim"); auto hEsim_ele = d_ele.Histo1D({"hEsim_ele", "Energy Deposit Electron; Energy Deposit [GeV]; Events", 10, 0.0, 0.05}, "Esim");
auto hEsim_ele_front = d_ele.Histo1D({"hEsim_ele_front", "Energy Deposit Front Electron; Energy Deposit [GeV]; Events", 10, 0.0, 0.05}, "Esim_front"); auto hEsim_ele_front = d_ele.Histo1D({"hEsim_ele_front", "Energy Deposit Front Electron; Energy Deposit [GeV]; Events", 10, 0.0, 0.05}, "Esim_front");
auto hEsim_pim = d_pim.Histo1D({"hEsim_pim", "Energy Deposit Electron; Energy Deposit [GeV]; Events", 10, 0.0, 0.05}, "Esim"); auto hEsim_pim = d_pim.Histo1D({"hEsim_pim", "Energy Deposit Electron; Energy Deposit [GeV]; Events", 10, 0.0, 0.05}, "Esim");
auto hEsim_pim_front = d_pim.Histo1D({"hEsim_pim_front", "Energy Deposit Front Pion-; Energy Deposit [GeV]; Events", 10, 0.0, 0.05}, "Esim_front"); auto hEsim_pim_front = d_pim.Histo1D({"hEsim_pim_front", "Energy Deposit Front Pion-; Energy Deposit [GeV]; Events", 10, 0.0, 0.05}, "Esim_front");
addDetectorName(detector_name, hEsim_ele.GetPtr());
addDetectorName(detector_name, hEsim_pim_front.GetPtr());
addDetectorName(detector_name, hEsim_pim.GetPtr());
addDetectorName(detector_name, hEsim_pim_front.GetPtr());
auto hEsim_ele_front_rej = d_ele_rej.Histo1D({"hEsim_ele_front_rej", "Energy Deposit Front Electron; Energy Deposit [GeV]; Events", 10, 0.0, 0.05}, "Esim_front"); auto hEsim_ele_front_rej = d_ele_rej.Histo1D({"hEsim_ele_front_rej", "Energy Deposit Front Electron; Energy Deposit [GeV]; Events", 10, 0.0, 0.05}, "Esim_front");
auto hEsim_pim_front_rej = d_pim_rej.Histo1D({"hEsim_pim_front_rej", "Energy Deposit Front Pion-; Energy Deposit [GeV]; Events", 10, 0.0, 0.05}, "Esim_front"); auto hEsim_pim_front_rej = d_pim_rej.Histo1D({"hEsim_pim_front_rej", "Energy Deposit Front Pion-; Energy Deposit [GeV]; Events", 10, 0.0, 0.05}, "Esim_front");
addDetectorName(detector_name, hEsim_ele_front_rej.GetPtr());
addDetectorName(detector_name, hEsim_pim_front_rej.GetPtr());
auto hEpvp_ele = d_ele.Histo2D({"hEpvp_ele", "Energy Deposit 1st 4 Layers/P vs P; E/P; P [GeV]", 20, 0.0, 0.01, 20, 1.0, 7.0}, "EOverP", "Pthr"); auto hEpvp_ele = d_ele.Histo2D({"hEpvp_ele", "Energy Deposit 1st 4 Layers/P vs P; E/P; P [GeV]", 20, 0.0, 0.01, 20, 1.0, 7.0}, "EOverP", "Pthr");
auto hEpvp_pim = d_pim.Histo2D({"hEpvp_pim", "Energy Deposit 1st 4 Layers/P vs P; E/P; P [GeV]", 20, 0.0, 0.01, 20, 1.0, 7.0}, "EOverP", "Pthr"); auto hEpvp_pim = d_pim.Histo2D({"hEpvp_pim", "Energy Deposit 1st 4 Layers/P vs P; E/P; P [GeV]", 20, 0.0, 0.01, 20, 1.0, 7.0}, "EOverP", "Pthr");
addDetectorName(detector_name, hEpvp_ele.GetPtr());
addDetectorName(detector_name, hEpvp_pim.GetPtr());
TH1D* hElePurity_initial = (TH1D *)hEsim_ele -> Clone(); TH1D* hElePurity_initial = (TH1D *)hEsim_ele -> Clone();
hElePurity_initial -> Divide(hEsim.GetPtr()); hElePurity_initial -> Divide(hEsim.GetPtr());
hElePurity_initial -> SetTitle("Electron/Pion Rejection"); hElePurity_initial -> SetTitle("Electron/Pion Rejection");
addDetectorName(detector_name, hElePurity_initial);
TH1D* hElePurity_ele = (TH1D *)hEsim_ele_front -> Clone(); TH1D* hElePurity_ele = (TH1D *)hEsim_ele_front -> Clone();
hElePurity_ele -> Divide(hEsim_front.GetPtr()); hElePurity_ele -> Divide(hEsim_front.GetPtr());
hElePurity_ele -> SetTitle("Electron/Pion Rejection : Electron"); hElePurity_ele -> SetTitle("Electron/Pion Rejection : Electron");
addDetectorName(detector_name, hElePurity_ele);
TH1D* hElePurity_pim = (TH1D *)hEsim_pim_front -> Clone(); TH1D* hElePurity_pim = (TH1D *)hEsim_pim_front -> Clone();
hElePurity_pim -> Divide(hEsim_front.GetPtr()); hElePurity_pim -> Divide(hEsim_front.GetPtr());
hElePurity_pim -> SetTitle("Electron/Pion Rejection : Pion Minus"); hElePurity_pim -> SetTitle("Electron/Pion Rejection : Pion Minus");
addDetectorName(detector_name, hElePurity_pim);
// Rejection // Rejection
TH1D* hElePurity_rej = (TH1D *)hEsim_ele_front_rej -> Clone(); TH1D* hElePurity_rej = (TH1D *)hEsim_ele_front_rej -> Clone();
hElePurity_rej -> Add(hEsim_pim_front_rej.GetPtr()); hElePurity_rej -> Add(hEsim_pim_front_rej.GetPtr());
hElePurity_rej -> SetTitle("Electron/Pion Rejection"); hElePurity_rej -> SetTitle("Electron/Pion Rejection");
addDetectorName(detector_name, hElePurity_rej);
TH1D* hElePurity_ele_rej = (TH1D *)hEsim_ele_front_rej -> Clone(); TH1D* hElePurity_ele_rej = (TH1D *)hEsim_ele_front_rej -> Clone();
hElePurity_ele_rej -> Divide(hElePurity_rej); hElePurity_ele_rej -> Divide(hElePurity_rej);
hElePurity_ele_rej -> SetTitle("Electron/Pion Rejection : Electron"); hElePurity_ele_rej -> SetTitle("Electron/Pion Rejection : Electron");
addDetectorName(detector_name, hElePurity_ele_rej);
TH1D* hElePurity_pim_rej = (TH1D *)hEsim_pim_front_rej -> Clone(); TH1D* hElePurity_pim_rej = (TH1D *)hEsim_pim_front_rej -> Clone();
hElePurity_pim_rej -> Divide(hElePurity_rej); hElePurity_pim_rej -> Divide(hElePurity_rej);
hElePurity_pim_rej -> SetTitle("Electron/Pion Rejection : Pi Minus"); hElePurity_pim_rej -> SetTitle("Electron/Pion Rejection : Pi Minus");
addDetectorName(detector_name, hElePurity_pim_rej);
// Event Counts // Event Counts
auto nevents_thrown = d1.Count(); auto nevents_thrown = d1.Count();
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment