diff --git a/benchmarks/imaging_ecal/analysis/emcal_barrel_pion_rejection_analysis.cxx b/benchmarks/imaging_ecal/analysis/emcal_barrel_pion_rejection_analysis.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b07f5728b38685f19831b33fde5c01f905d3e2f1
--- /dev/null
+++ b/benchmarks/imaging_ecal/analysis/emcal_barrel_pion_rejection_analysis.cxx
@@ -0,0 +1,946 @@
+////////////////////////////////////////
+// Read reconstruction ROOT output file
+// Plot variables
+////////////////////////////////////////
+
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+#include <algorithm>
+#include <string>
+
+#include "dd4pod/CalorimeterHitCollection.h"
+#include "dd4pod/Geant4ParticleCollection.h"
+
+#include "eicd/CalorimeterHitCollection.h"
+#include "eicd/CalorimeterHitData.h"
+#include "eicd/ClusterCollection.h"
+#include "eicd/ClusterData.h"
+#include "eicd/ClusterLayerData.h"
+
+#include "eicd/RawCalorimeterHitCollection.h"
+#include "eicd/RawCalorimeterHitData.h"
+
+#include "common_bench/benchmark.h"
+#include "common_bench/mt.h"
+#include "common_bench/util.h"
+
+#include <boost/range/combine.hpp>
+
+#include "DD4hep/IDDescriptor.h"
+#include "DD4hep/Readout.h"
+#include "DD4hep/Segmentations.h"
+
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TMath.h"
+#include "TH1.h"
+#include "TF1.h"
+#include "TH1D.h"
+#include "TH2D.h"
+#include "TFitResult.h"
+#include "TLegend.h"
+#include "TString.h"
+#include "TGraph.h"
+#include "TGraph2D.h"
+#include "TGraphErrors.h"
+#include "TObject.h"
+#include "TLine.h"
+#include "TError.h"
+#include "TKey.h"
+#include "TList.h"
+#include "TCollection.h"
+#include "TSystem.h"
+#include "TROOT.h"
+
+R__LOAD_LIBRARY(libfmt.so)
+#include "fmt/core.h"
+#include "DD4hep/Detector.h"
+#include "DDG4/Geant4Data.h"
+#include "DDRec/CellIDPositionConverter.h"
+#include <nlohmann/json.hpp>
+
+using ROOT::RDataFrame;
+using namespace ROOT::VecOps;
+
+void emcal_barrel_pion_rejection_analysis(
+ const char* input_fname1 = "results/rec_emcal_barrel_pion_rej_electron.root",
+ const char* input_fname2 = "results/rec_emcal_barrel_pion_rej_piminus.root"
+ //const char* input_fname1 = "../rec_electron_4k.root",
+ //const char* input_fnamE1 = "../rec_piminus_4k.root"//16
+ //const char* input_fname1 = "../rec_emcal_barrel_electron_new4k.root",
+ //const char* input_fnamE1 = "../rec_emcal_barrel_piminus_new16k.root"//16
+ //const char* input_fname1 = "../rec_emcal_barrel_pion_rej_electron_new_E1.root",
+ //const char* input_fname2 = "../rec_emcal_barrel_pion_rej_piminus_new_E1.root"//100k
+ )
+{
+
+ // Error Ignore Level Set
+ //gErrorIgnoreLevel = kFatal;
+
+ // Setting for graphs
+ gROOT->SetStyle("Plain");
+ gStyle->SetOptFit(1);
+ gStyle->SetLineWidth(2);
+ gStyle->SetPadTickX(1);
+ gStyle->SetPadTickY(1);
+ gStyle->SetPadGridX(1);
+ gStyle->SetPadGridY(1);
+ gStyle->SetPadLeftMargin(0.14);
+ gStyle->SetPadRightMargin(0.14);
+
+ ROOT::EnableImplicitMT();
+ ROOT::RDataFrame d0("events", {input_fname1, input_fname2});
+
+ // Environment Variables
+ std::string detector_path = "";
+ std::string detector_name = "athena";//athena
+ if(std::getenv("DETECTOR_PATH")) {
+ detector_path = std::getenv("DETECTOR_PATH");
+ }
+ if(std::getenv("JUGGLER_DETECTOR")) {
+ detector_name = std::getenv("JUGGLER_DETECTOR");
+ }
+
+ // MCParticles Functions/////////////////////////////////////////////////////////////////////////////////////
+
+ // Filter function to get electrons
+ auto is_electron = [](std::vector<dd4pod::Geant4ParticleData> const& input){
+ return (input[2].pdgID == 11 ? true : false);
+ };
+
+ // Filter function to get just negative pions
+ auto is_piMinus = [](std::vector<dd4pod::Geant4ParticleData> const& input){
+ return (input[2].pdgID == -211 ? true : false);
+ };
+
+ // Thrown Energy [GeV]
+ auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
+ return TMath::Sqrt(input[2].ps.x*input[2].ps.x + input[2].ps.y*input[2].ps.y + input[2].ps.z*input[2].ps.z + input[2].mass*input[2].mass);
+ };
+
+ // Thrown Momentum [GeV]
+ auto Pthr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
+ return TMath::Sqrt(input[2].ps.x*input[2].ps.x + input[2].ps.y*input[2].ps.y + input[2].ps.z*input[2].ps.z);
+ };
+
+ // Thrown Eta
+ auto Eta = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
+ double E = TMath::Sqrt(input[2].ps.x*input[2].ps.x + input[2].ps.y*input[2].ps.y + input[2].ps.z*input[2].ps.z + input[2].mass*input[2].mass);
+ return 0.5*TMath::Log((E + input[2].ps.z) / (E - input[2].ps.z));
+ };
+
+ // Thrown pT [GeV]
+ auto pT = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
+ return TMath::Sqrt(input[2].ps.x*input[2].ps.x + input[2].ps.y*input[2].ps.y);
+ };
+
+ // RecoEcalBarrelImagingHits Functions/////////////////////////////////////////////////////////////////////////////////////
+ int layerNumImg;
+ // Energy deposition [GeV]
+ auto ERecoImg = [](const std::vector<eic::CalorimeterHitData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ total_edep += i.energy;
+ }
+ return total_edep;
+ };
+
+ // Nhits in layer
+ auto nhitsImgLayer = [&layerNumImg](const std::vector<eic::CalorimeterHitData>& evt) {
+ int nhitsTot;
+ for (const auto& i: evt){
+ if (i.layer == layerNumImg){nhitsTot++;}
+ }
+ return nhitsTot;
+ };
+
+ // Energy deposition at a layer
+ auto ERecoImgLayer = [&layerNumImg](const std::vector<eic::CalorimeterHitData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ if (i.layer == layerNumImg){total_edep += i.energy;}
+ }
+ return total_edep;
+ };
+ //###########################################################################################
+ // Energy deposition at a layer
+ auto ERecoImgLayer1 = [](const std::vector<eic::CalorimeterHitData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ if (i.layer == 1){total_edep += i.energy;}
+ }
+ return total_edep;
+ };
+ auto ERecoImgLayer2 = [](const std::vector<eic::CalorimeterHitData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ if (i.layer == 2){total_edep += i.energy;}
+ }
+ return total_edep;
+ };
+ auto ERecoImgLayer3 = [](const std::vector<eic::CalorimeterHitData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ if (i.layer == 3){total_edep += i.energy;}
+ }
+ return total_edep;
+ };
+ auto ERecoImgLayer4 = [](const std::vector<eic::CalorimeterHitData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ if (i.layer == 4){total_edep += i.energy;}
+ }
+ return total_edep;
+ };
+ auto ERecoImgLayer5 = [](const std::vector<eic::CalorimeterHitData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ if (i.layer == 5){total_edep += i.energy;}
+ }
+ return total_edep;
+ };
+ auto ERecoImgLayer6 = [](const std::vector<eic::CalorimeterHitData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ if (i.layer == 6){total_edep += i.energy;}
+ }
+ return total_edep;
+ };
+ //###########################################################################################
+
+ // Energy deposition at a layer / Energy of all layers
+ auto ERecoImgLayerOverERecoImg = [&layerNumImg](const std::vector<eic::CalorimeterHitData>& evt, const double energy) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ if (i.layer == layerNumImg){total_edep += i.energy;}
+ }
+ return total_edep / energy;
+ };
+
+ // Nhits / Energy Deposition in that layer
+ auto nhitsImgLayerOverERecoImg = [&layerNumImg](const std::vector<eic::CalorimeterHitData>& evt) {
+ double total_edep = 0.0;
+ int nhits;
+ for (const auto& i: evt){
+ if (i.layer == layerNumImg){total_edep += i.energy; nhits++;}
+ }
+ return (double)nhits / total_edep;
+ };
+
+ // EcalBarrelImagingClustersLayers Functions////////////////////////////////////////////////////////////////////////////////
+ // Number of hits
+ auto nhitsClusLayerImg = [] (const std::vector<eic::ClusterLayerData>& evt) {
+ int nhitsTot = 0;
+ for (const auto &i : evt){
+ nhitsTot += i.nhits;
+ }
+ return nhitsTot;
+ };
+
+ // Number of clusters
+ auto nClusLayerImg = [] (const std::vector<eic::ClusterLayerData>& evt) {return (int) evt.size(); };
+
+ // Energy deposition in cluster [GeV]
+ auto EClusLayerImg = [](const std::vector<eic::ClusterLayerData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ total_edep += i.energy;
+ }
+ return total_edep;
+ };
+
+ // Max Energy deposition in cluster [GeV]
+ auto EClusLayerMaxImg = [](const std::vector<eic::ClusterLayerData>& evt) {
+ double max = 0.0;
+ for (const auto& i: evt){
+ if (i.energy > max){max = i.energy;}
+ }
+ return max;
+ };
+
+ // Min Energy deposition in cluster [GeV]
+ auto EClusLayerMax2Img = [](const std::vector<eic::ClusterLayerData>& evt) {
+ double max1 = 0.0;
+ double max2 = 0.0;
+ for (const auto& i: evt){
+ if (i.energy > max1){max2 = max1; max1 = i.energy;}
+ }
+ return max2;
+ };
+
+ // EcalBarrelImagingClusters Functions////////////////////////////////////////////////////////////////////
+ // May want to look at position data as well
+
+ // Energy deposition [GeV]
+ auto EClusImg = [](const std::vector<eic::ClusterData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ total_edep += i.energy;
+ }
+ return total_edep;
+ };
+
+ // Number of energy deposits
+ auto nhitsClusImg = [](const std::vector<eic::ClusterData>& evt) {
+ int nhitsTot = 0;
+ for (const auto& i: evt){
+ nhitsTot += i.nhits;
+ }
+ return (double)nhitsTot;
+ };
+
+ // Mean number of energy deposits
+ auto nhitsClusMeanImg = [](const std::vector<eic::ClusterData>& evt) {
+ int nhitsTot = 0;
+ int count = 0;
+ for (const auto& i: evt){
+ nhitsTot += i.nhits;
+ count++;
+ }
+ return (double)nhitsTot / (double)count;
+ };
+
+ // EcalBarrelScFiHitsReco Functions/////////////////////////////////////////////////////////////////////////////////////
+ int layerNumScFi;
+ // Energy deposition [GeV]
+ auto ERecoScFi = [](const std::vector<eic::CalorimeterHitData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ total_edep += i.energy;
+ }
+ return total_edep;
+ };
+
+ // Nhits in layer
+ auto nhitsScFiLayer = [&layerNumScFi](const std::vector<eic::CalorimeterHitData>& evt) {
+ int nhits;
+ for (const auto& i: evt){
+ if (i.layer == layerNumScFi){nhits++;}
+ }
+ return nhits;
+ };
+
+ // Energy deposition at a layer
+ auto ERecoScFiLayer = [&layerNumScFi](const std::vector<eic::CalorimeterHitData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ if (i.layer == layerNumScFi){total_edep += i.energy;}
+ }
+ return total_edep;
+ };
+
+ // STD Deviation
+ auto ERecoSTDScFi = [](const std::vector<eic::CalorimeterHitData>& evt, const double total_energy) {
+ double std = 0.0;
+ double count = 0.0;
+ for (const auto& i: evt){
+ std += std::pow(i.energy - total_energy, 2.0);
+ count++;
+ }
+ //std::printf("%e %f\n", std, count);
+ if (count == 0){count = 1;}
+ return std::sqrt(std / count);
+ };
+ //########################################################################################
+ auto ERecoScFiLayer1 = [](const std::vector<eic::CalorimeterHitData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ if (i.layer == 1){total_edep += i.energy;}
+ }
+ return total_edep;
+ };
+ auto ERecoScFiLayer2 = [](const std::vector<eic::CalorimeterHitData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ if (i.layer == 2){total_edep += i.energy;}
+ }
+ return total_edep;
+ };
+ auto ERecoScFiLayer3 = [](const std::vector<eic::CalorimeterHitData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ if (i.layer == 3){total_edep += i.energy;}
+ }
+ return total_edep;
+ };
+ auto ERecoScFiLayer4 = [](const std::vector<eic::CalorimeterHitData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ if (i.layer == 4){total_edep += i.energy;}
+ }
+ return total_edep;
+ };
+ auto ERecoScFiLayer5 = [](const std::vector<eic::CalorimeterHitData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ if (i.layer == 5){total_edep += i.energy;}
+ }
+ return total_edep;
+ };
+ auto ERecoScFiLayer6 = [](const std::vector<eic::CalorimeterHitData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ if (i.layer == 6){total_edep += i.energy;}
+ }
+ return total_edep;
+ };
+ auto ERecoScFiLayer7 = [](const std::vector<eic::CalorimeterHitData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ if (i.layer == 7){total_edep += i.energy;}
+ }
+ return total_edep;
+ };
+ auto ERecoScFiLayer8 = [](const std::vector<eic::CalorimeterHitData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ if (i.layer == 8){total_edep += i.energy;}
+ }
+ return total_edep;
+ };
+ auto ERecoScFiLayer9 = [](const std::vector<eic::CalorimeterHitData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ if (i.layer == 9){total_edep += i.energy;}
+ }
+ return total_edep;
+ };
+ //#########################################################################################
+
+ // Energy deposition at a layer / Energy of all layers
+ auto ERecoScFiLayerOverERecoScFi = [&layerNumScFi](const std::vector<eic::CalorimeterHitData>& evt, const double energy) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ if (i.layer == layerNumScFi){total_edep += i.energy;}
+ }
+ return total_edep / energy;
+ };
+
+ // Energy deposition at layer 4 / Energy of all layers
+ auto ERecoScFiLayerOverERecoScFi4 = [](const std::vector<eic::CalorimeterHitData>& evt, const double energy) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ if (i.layer == 4){total_edep += i.energy;}
+ }
+ return total_edep / energy;
+ };
+
+ // Nhits / Energy Deposition in that layer
+ auto nhitsScFiLayerOverERecoScFi = [&layerNumScFi](const std::vector<eic::CalorimeterHitData>& evt) {
+ double total_edep = 0.0;
+ int nhits;
+ for (const auto& i: evt){
+ if (i.layer == layerNumScFi){total_edep += i.energy; nhits++;}
+ }
+ return (double)nhits / total_edep;
+ };
+
+ // EcalBarrelScFiClusters Functions////////////////////////////////////////////////////////////////////
+ // May want to look at position data as well
+
+ // Energy deposition [GeV]
+ auto EClusScFi = [](const std::vector<eic::ClusterData>& evt) {
+ double total_edep = 0.0;
+ for (const auto& i: evt){
+ total_edep += i.energy;
+ }
+ return total_edep;
+ };
+
+ // Number of energy deposits
+ auto nhitsClusScFi = [](const std::vector<eic::ClusterData>& evt) {
+ int nhitsTot = 0;
+ for (const auto& i: evt){
+ nhitsTot += i.nhits;
+ }
+ return (double)nhitsTot;
+ };
+
+ // Mean number of energy deposits
+ auto nhitsClusMeanScFi = [](const std::vector<eic::ClusterData>& evt) {
+ int nhitsTot = 0;
+ int count = 0;
+ for (const auto& i: evt){
+ nhitsTot += i.nhits;
+ count++;
+ }
+ return (double)nhitsTot / (double)count;
+ };
+
+ std::string std_str = "TMath::Sqrt(TMath::Power(ERecoScFiLayer1-ERecoScFi/9.0,2)+";
+ std_str+= "TMath::Power(ERecoScFiLayer2-ERecoScFi/9.0,2)+TMath::Power(ERecoScFiLayer3-ERecoScFi/9.0,2)+";
+ std_str+= "TMath::Power(ERecoScFiLayer4-ERecoScFi/9.0,2)+TMath::Power(ERecoScFiLayer5-ERecoScFi/9.0,2)+";
+ std_str+= "TMath::Power(ERecoScFiLayer6-ERecoScFi/9.0,2)+TMath::Power(ERecoScFiLayer7-ERecoScFi/9.0,2)+";
+ std_str+= "TMath::Power(ERecoScFiLayer8-ERecoScFi/9.0,2)+TMath::Power(ERecoScFiLayer9-ERecoScFi/9.0,2))/TMath::Sqrt(9)";
+
+ // Define variables
+ auto d1 = d0
+ .Define("Ethr", Ethr, {"mcparticles"})
+ .Define("Pthr", Pthr, {"mcparticles"})
+ .Define("Eta", Eta, {"mcparticles"})
+ .Define("pT", pT, {"mcparticles"})
+
+ .Define("ERecoImg", ERecoImg, {"RecoEcalBarrelImagingHits"})
+ .Define("nhitsImgLayer", nhitsImgLayer, {"RecoEcalBarrelImagingHits"})
+ .Define("ERecoImgLayer", ERecoImgLayer, {"RecoEcalBarrelImagingHits"})
+ .Define("ERecoImgLayerOverERecoImg", ERecoImgLayerOverERecoImg, {"RecoEcalBarrelImagingHits", "ERecoImg"})
+ .Define("nhitsImgLayerOverERecoImg", nhitsImgLayerOverERecoImg, {"RecoEcalBarrelImagingHits"})
+
+ .Define("nhitsClusLayerImg", nhitsClusLayerImg, {"EcalBarrelImagingClustersLayers"})
+ .Define("nClusLayerImg", nClusLayerImg, {"EcalBarrelImagingClustersLayers"})
+ .Define("EClusLayerImg", EClusLayerImg, {"EcalBarrelImagingClustersLayers"})
+ .Define("EClusLayerMaxImg", EClusLayerMaxImg, {"EcalBarrelImagingClustersLayers"})
+ .Define("EClusLayerMax2Img", EClusLayerMax2Img, {"EcalBarrelImagingClustersLayers"})
+ .Define("EClusLayerMaxOverECluster", "EClusLayerMaxImg/EClusLayerImg")
+ .Define("EClusLayerMax2OverECluster", "EClusLayerMax2Img/EClusLayerImg")
+
+ .Define("EClusImg", EClusImg, {"EcalBarrelImagingClusters"})
+ .Define("nhitsClusImg", nhitsClusImg, {"EcalBarrelImagingClusters"})
+ .Define("nhitsClusMeanImg", nhitsClusMeanImg, {"EcalBarrelImagingClusters"})
+
+ .Define("ERecoScFi", ERecoScFi, {"EcalBarrelScFiHitsReco"})
+ .Define("nhitsScFiLayer", nhitsScFiLayer, {"EcalBarrelScFiHitsReco"})
+ .Define("ERecoScFiLayer", ERecoScFiLayer, {"EcalBarrelScFiHitsReco"})
+ .Define("ERecoSTDScFi", ERecoSTDScFi, {"EcalBarrelScFiHitsReco", "ERecoScFi"})
+ .Define("ERecoScFiLayerOverERecoScFi", ERecoScFiLayerOverERecoScFi, {"EcalBarrelScFiHitsReco", "ERecoScFi"})
+ .Define("ERecoScFiLayerOverERecoScFi4", ERecoScFiLayerOverERecoScFi4, {"EcalBarrelScFiHitsReco", "ERecoScFi"})
+ .Define("nhitsScFiLayerOverERecoScFi", nhitsScFiLayerOverERecoScFi, {"EcalBarrelScFiHitsReco"})
+
+ .Define("EClusScFi", EClusScFi, {"EcalBarrelScFiClusters"})
+ .Define("nhitsClusScFi", nhitsClusScFi, {"EcalBarrelScFiClusters"})
+ .Define("nhitsClusMeanScFi", nhitsClusMeanScFi, {"EcalBarrelScFiClusters"})
+ .Define("EClusDiff", "EClusScFi-EClusImg")
+ .Define("EClusRes", "EClusDiff/EClusScFi")
+ .Define("nhitsDiff", "nhitsClusScFi-nhitsClusImg")
+
+ .Define("EClusImgOverP", "EClusImg/Pthr")
+ .Define("EClusScFiOverP", "EClusScFi/Pthr")
+ .Define("ERecoImgLayerOverP", "ERecoImgLayer/Pthr")
+ .Define("ERecoScFiLayerOverP", "ERecoScFiLayer/Pthr")
+
+
+ .Define("ERecoDiff", "ERecoScFi-ERecoImg")
+ .Define("ERecoScFiLayer1", ERecoScFiLayer1, {"EcalBarrelScFiHitsReco"})
+ .Define("ERecoScFiLayer2", ERecoScFiLayer2, {"EcalBarrelScFiHitsReco"})
+ .Define("ERecoScFiLayer3", ERecoScFiLayer3, {"EcalBarrelScFiHitsReco"})
+ .Define("ERecoScFiLayer4", ERecoScFiLayer4, {"EcalBarrelScFiHitsReco"})
+ .Define("ERecoScFiLayer5", ERecoScFiLayer5, {"EcalBarrelScFiHitsReco"})
+ .Define("ERecoScFiLayer6", ERecoScFiLayer6, {"EcalBarrelScFiHitsReco"})
+ .Define("ERecoScFiLayer7", ERecoScFiLayer7, {"EcalBarrelScFiHitsReco"})
+ .Define("ERecoScFiLayer8", ERecoScFiLayer8, {"EcalBarrelScFiHitsReco"})
+ .Define("ERecoScFiLayer9", ERecoScFiLayer9, {"EcalBarrelScFiHitsReco"})
+ .Define("ERecoScFiLayerSTD", std_str)
+ .Define("ERecoScFiLayer10verP", "ERecoScFiLayer1/Pthr")
+ .Define("ERecoScFiLayer2OverP", "ERecoScFiLayer2/Pthr")
+ .Define("ERecoScFiLayer3OverP", "ERecoScFiLayer3/Pthr")
+ .Define("ERecoScFiLayer4OverP", "ERecoScFiLayer4/Pthr")
+ .Define("ERecoScFiLayer5OverP", "ERecoScFiLayer5/Pthr")
+ .Define("ERecoScFiLayer6OverP", "ERecoScFiLayer6/Pthr")
+ .Define("ERecoScFiLayer7OverP", "ERecoScFiLayer7/Pthr")
+ .Define("ERecoScFiLayer8OverP", "ERecoScFiLayer8/Pthr")
+ .Define("ERecoScFiLayer9OverP", "ERecoScFiLayer9/Pthr")
+ ;
+
+ // Preliminary Cut
+ int binNum = 50;
+ std::string cut = "nhitsClusImg>0&&nhitsClusScFi>0";
+
+ // Preliminary Plots
+ // Distributions eventually used for the cuts
+ std::vector<std::string> colList = {"EClusImg", "EClusScFi", "nhitsClusImg", "nhitsDiff", "nhitsClusScFi",
+ "ERecoScFiLayerOverERecoScFi4", "EClusRes", "EClusImgOverP", "EClusScFiOverP",
+ "ERecoScFiLayer3OverP", "ERecoScFiLayer4OverP", "ERecoScFiLayer5OverP","ERecoScFiLayer6OverP"};
+ std::vector<std::vector<double>> colRange = {{0,23},{0,19},{0,570},{-300,300},{0,280},{0.065,0.45},{-0.5,1},{0,1.5},{0,1.5},{0,2.5e-2},{0,2.5e-2},{0,2.5e-2},{0,2.5e-2},};
+
+ auto d_temp_ele = d1.Filter(is_electron, {"mcparticles"});
+ auto d_temp_pim = d1.Filter(is_piMinus, {"mcparticles"});
+
+ int rangeCount = 0;
+ for (auto && name : colList){
+ double min = colRange[rangeCount][0];
+ double max = colRange[rangeCount][1];
+
+ auto he = d_temp_ele.Filter(cut.c_str()).Histo1D({"he", name.c_str(), binNum, min, max},name.c_str());
+ auto hp = d_temp_pim.Filter(cut.c_str()).Histo1D({"hp", name.c_str(), binNum, min, max},name.c_str());
+
+ auto c = new TCanvas("c", "c", 700, 500);
+ hp->GetYaxis()->SetTitleOffset(1.4);
+ he->SetLineWidth(2);
+ he->SetLineColor(kBlue);
+
+ hp->SetLineWidth(2);
+ hp->SetLineColor(kRed);
+ hp->DrawClone();
+ he->DrawClone("same");
+
+ auto leng = new TLegend(0.7, 0.7, 0.9, 0.9);
+ leng->AddEntry(he.GetPtr(),("e^{-}(" + std::to_string((int)he->Integral()) + ")").c_str(),"l");
+ leng->AddEntry(hp.GetPtr(),("#pi^{-}(" + std::to_string((int)hp->Integral()) + ")").c_str(),"l");
+ leng->Draw();
+ c->SaveAs((fmt::format("results/emcal_barrel_preliminary_testing_{}.png", name)).c_str());
+ rangeCount++;
+ }
+
+ // Gather Cuts//////////////////////////////////////////////////////////////////////////////////
+ // The cuts were generated and tested using the code below. Currently the resultant code/cuts are just saved.
+ std::vector<std::string> cutVec;
+ string finalCut = "(";
+ double deltag = 50;
+/*
+ for (double g = -1; g < 18; g+=deltag){
+ // Gather Primary Cut
+ double fitMin = 100;
+ double fitMax = 200;
+ fitMin = 0;
+ fitMax = 100;
+ string energyCut = fmt::format("EClusScFi>={}&&EClusScFi<{}", g, g + deltag);
+ auto he_cut = d_temp_ele.Filter(energyCut.c_str()).Histo1D({"", "", binNum, fitMin, fitMax}, "nhitsClusScFi");
+ auto he_cut_copy = he_cut->DrawCopy();
+ he_cut_copy->Fit("gaus", "", "", fitMin, fitMax);
+ double* res_cut = he_cut_copy->GetFunction("gaus")->GetParameters();
+ string cut2 = fmt::format("nhitsClusScFi<{}", res_cut[1] + 2.0*res_cut[2]);
+ cut2 ="nhitsClusScFi<80";// < 2Gev
+
+ string cut3 = cut3 = cut;//Changed for test
+
+ // Gather Teritiary Cut
+ // Note a > 0.03 cut only removes ~3 events.
+ // A > 0.1 cut yields a 0.955 electrion eff, which on top of the previous eff would yield 0.955*0.95 ~0.91
+ // At the moment it is set to the > 0.03 cut for further testing
+ string cut4 = fmt::format("ERecoScFiLayerOverERecoScFi4>{}", 0.03);
+ finalCut += fmt::format("(EClusScFi>={}&&EClusScFi<{}&&{}&&{}&&{})||", g, g + deltag, cut2, cut3, cut4);
+ cutVec.push_back("nhitsClusImg>0");
+ cutVec.push_back("nhitsClusScFi>0");
+ cutVec.push_back(cut2);
+ cutVec.push_back(cut3);
+ cutVec.push_back("EClusScFi>-1");//cut4
+ }
+ finalCut.pop_back();
+ finalCut.pop_back();
+ finalCut+="||(EClusScFi>=2&&EClusScFi<4&&ERecoScFiLayerOverERecoScFi4>0.09&&nhitsClusImg>20)||(EClusScFi>=4&&EClusScFi<6&&ERecoScFiLayerOverERecoScFi4>0.08)";
+ finalCut+=")&&" + cut;
+ //finalCut+="&&EClusRes<0.8";
+ finalCut+="&&nhitsDiff<100";
+ cutVec.push_back("(EClusScFi>=2&&EClusScFi<4&&ERecoScFiLayerOverERecoScFi4>0.09&&nhitsClusImg>20)");
+ cutVec.push_back("(EClusScFi>=4&&EClusScFi<6&&ERecoScFiLayerOverERecoScFi4>0.08)");
+ cutVec.push_back("EClusRes<1.5");//0.8 for E >2
+ cutVec.push_back("nhitsDiff<100");
+
+ // E/p Cuts
+ std::vector<double> val = {2.5, 3, 4};//334
+ std::string newcut;
+ int count = 3;
+ for (auto &&num : val){
+ std::string colName = (fmt::format("ERecoScFiLayer{}OverP", count)).c_str();
+ auto he_cut = d_temp_ele.Filter(finalCut.c_str()).Histo1D({"", "", binNum, 0, 4e-2}, colName);//2.5e-2 for >2
+ auto he_cut_copy = he_cut->DrawCopy();
+
+ he_cut_copy->Fit("gaus", "", "", 0, 2.5e-2);
+ double* res_cut = he_cut_copy->GetFunction("gaus")->GetParameters();
+ newcut += (fmt::format("{}>={}", colName, res_cut[1] - num*res_cut[2])).c_str();
+ newcut += "&&";
+ count++;
+ //if (colName == "ERecoScFiLayer3OverP"){num=3;}
+ cutVec.push_back(fmt::format("{}>={}", colName, res_cut[1] - num*res_cut[2]));
+ }
+ newcut.pop_back();
+ newcut.pop_back();
+ finalCut+= "&&" + newcut;
+ finalCut+= "&&ERecoScFiLayer2OverP>0.005";// only for 1 GeV
+ finalCut+= "&&nhitsDiff<60";// only for 1 GeV
+ finalCut+= "&&nhitsClusMeanImg>1";// only for 1 GeV
+*/
+
+ // All Cuts/////////////////////////////////////////////////////////////////////////////////////////
+ finalCut="(";
+ finalCut+="(EClusScFi<0.5||";
+ finalCut+="(EClusScFi>=0.5&&EClusScFi<1.5&&nhitsClusScFi<80&&ERecoScFiLayerOverERecoScFi4>0.03)||";
+ finalCut+="(EClusScFi>=1.5&&EClusScFi<4&&ERecoScFiLayerOverERecoScFi4>0.05&&nhitsClusImg>20)||";
+ finalCut+="(EClusScFi>=4&&EClusScFi<6&&ERecoScFiLayerOverERecoScFi4>0.08)||";
+ finalCut+="(EClusScFi>6))&&";
+ finalCut+="nhitsClusImg>0&&nhitsClusScFi>0&&";
+ finalCut+="ERecoScFiLayer3OverP>=0.0024095049104228966&&";
+ finalCut+="ERecoScFiLayer2OverP>0.005&&";
+ finalCut+="nhitsDiff<60&&";
+ finalCut+="nhitsClusMeanImg>1";
+ finalCut+=")";
+ finalCut+="||";
+ finalCut+="(EClusScFi>3&&(";
+ finalCut+="(nhitsClusScFi<292.16412543952&&";
+ finalCut+="nhitsClusImg>0&&";
+ finalCut+="nhitsClusScFi>0&&";
+ finalCut+="ERecoScFiLayerOverERecoScFi4>0.03)||";
+ finalCut+="(EClusScFi>=2&&EClusScFi<4&&";
+ finalCut+="ERecoScFiLayerOverERecoScFi4>0.09&&";
+ finalCut+="nhitsClusImg>20)||";
+ finalCut+="(EClusScFi>=4&&EClusScFi<6&&ERecoScFiLayerOverERecoScFi4>0.08))&&EClusRes<0.8&&nhitsDiff<100&&ERecoScFiLayer3OverP>=0.0033134401319620185&&";
+ finalCut+="ERecoScFiLayer4OverP>=0.005112020043084812&&ERecoScFiLayer5OverP>=0.0010761376551889999";
+ finalCut+="&&((EClusScFi<5 && nhitsDiff < 80)||(EClusScFi>=5 && nhitsDiff<60))";
+ finalCut+="&&EClusScFiOverP>0.7";
+ finalCut+="&&((EClusScFi>4 && nhitsDiff < 40) || EClusScFi <=4)";
+ finalCut+="&&((EClusScFi>8 && nhitsDiff < -40) || EClusScFi <=8)";
+ finalCut+=")";
+
+ // Filtered distributions//////////////////////////////////////////////////////////////////////////////////////
+ auto d_ele = d_temp_ele.Filter(finalCut.c_str());
+ auto d_pim = d_temp_pim.Filter(finalCut.c_str());
+ int ele_count = *d_ele.Count();
+ int pim_count = *d_pim.Count();
+ int ele_count_pass_hits = *d_temp_ele.Filter("nhitsClusImg>0&&nhitsClusScFi>0").Count();
+ int pim_count_pass_hits = *d_temp_pim.Filter("nhitsClusImg>0&&nhitsClusScFi>0").Count();
+
+ ///////////////////////////////////////////////
+ //Generate Energy specific Plots
+/*
+ {
+ auto he1 = d_ele.Histo1D({"he", "EClusScFi/P : E = 5 GeV", binNum, 0.6, 1.2}, "EClusScFiOverP");
+ auto he_cut_copy = he1->DrawCopy();
+
+ he_cut_copy->Fit("gaus", "", "", 0.7, 1.2);
+ double* res_cut = he_cut_copy->GetFunction("gaus")->GetParameters();
+ double d = 6.5;
+ std::string newcut = (fmt::format("{}>={}", "EClusScFiOverP", res_cut[1] - d*res_cut[2])).c_str();
+ //newcut = "EClusScFi>-1";
+
+ auto he = d_ele.Filter(newcut.c_str()).Histo1D({"he", "EClusScFi/P : E = 5 GeV", binNum, 0.6, 1.2}, "EClusScFiOverP");
+ auto hp = d_pim.Filter(newcut.c_str()).Histo1D({"hp", "EClusScFi/P : E = 5 GeV", binNum, 0.6, 1.2}, "EClusScFiOverP");
+
+ auto c = new TCanvas("c", "c", 700, 500);
+ hp->GetYaxis()->SetTitleOffset(1.4);
+ he->SetLineWidth(2);
+ he->SetLineColor(kBlue);
+
+ hp->SetLineWidth(2);
+ hp->SetLineColor(kRed);
+ he->DrawClone();
+ hp->DrawClone("same");
+
+ gStyle->SetLegendTextSize(0.03);
+ auto leng = new TLegend(0.2, 0.7, 0.5, 0.9);
+ ele_count = *d_ele.Filter(newcut.c_str()).Count();
+ pim_count = *d_pim.Filter(newcut.c_str()).Count();
+ std::string text = std::to_string(100000.0 / pim_count);
+ text.resize(4);
+ leng->AddEntry((TObject*)0, ("Pion Rejection = " + text).c_str(), "");
+ text = std::to_string(ele_count / 10000.0);
+ text.resize(4);
+ leng->AddEntry((TObject*)0, ("#varepsilon_{e} = " + text).c_str(), "");
+ leng->AddEntry(he.GetPtr(),("e^{-}(" + std::to_string(ele_count) + ")").c_str(),"l");
+ leng->AddEntry(hp.GetPtr(),("#pi^{-}(" + std::to_string(pim_count) + ")").c_str(),"l");
+ leng->Draw();
+ c->SaveAs("results/emcal_barrel_CutComp_E1_EClusScFiOverP.png");
+ }
+*/
+ //////////////////////////////////////////////
+ //Testing sigmas to get 95% electron efficiency
+/*
+ {
+ auto he_cut = d_ele.Histo1D({"", "", binNum, 0.6, 1.2}, "EClusScFiOverP");
+ auto he_cut_copy = he_cut->DrawCopy();
+
+ he_cut_copy->Fit("gaus", "", "", 0.7, 1.2);
+ double* res_cut = he_cut_copy->GetFunction("gaus")->GetParameters();
+ for (double d = 5.5; d < 7.5; d+=0.5){
+ std::string newcut = (fmt::format("{}>={}", "EClusScFiOverP", res_cut[1] - d*res_cut[2])).c_str();
+ ele_count = *d_ele.Filter(newcut.c_str()).Count();
+ pim_count = *d_pim.Filter(newcut.c_str()).Count();
+ std::printf("E1 %fsigma : E %d P %d\n", d, ele_count, pim_count);
+ }
+ exit(0);
+ }
+ ////////////////////////////////////////////
+*/
+
+
+ //Final Plots ////////////////////////////////////////////////////////////////////////////////////////////////
+ rangeCount = 0;
+ for (auto && name : colList){
+ double min = colRange[rangeCount][0];
+ double max = colRange[rangeCount][1];
+ double mean = *d_ele.Mean(name);
+ double std = *d_ele.StdDev(name);
+ min = mean - 2.5*std;
+ max = mean + 2.5*std;
+
+ auto he = d_ele.Histo1D({"#e^{-}", name.c_str(), binNum, min, max},name.c_str());
+ auto hp = d_pim.Histo1D({"#pi^{-}", name.c_str(), binNum, min, max},name.c_str());
+
+ auto c = new TCanvas("c", "c", 700, 500);
+ hp->GetYaxis()->SetTitleOffset(1.4);
+ he->SetLineWidth(2);
+ he->SetLineColor(kBlue);
+
+ hp->SetLineWidth(2);
+ hp->SetLineColor(kRed);
+ he->DrawClone();
+ hp->DrawClone("same");
+
+ auto leng = new TLegend(0.7, 0.7, 0.9, 0.9);
+ leng->AddEntry(he.GetPtr(),"e^{-}","l");
+ leng->AddEntry(hp.GetPtr(),"#pi^{-}","l");
+ leng->Draw();
+ c->SaveAs((fmt::format("results/emcal_barrel_final_{}.png", name)).c_str());
+ rangeCount++;
+ }
+ ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+ // Pion Rejection Plot
+ auto tg = new TGraphErrors();
+ std::vector<std::vector<double>> EBins = {{0,2}, {2, 4}, {4, 6}, {6, 9}, {9, 12}, {12, 15}, {15, 18}};
+ std::vector<double> binEdges = {0, 2, 4, 6, 9, 12, 15, 18};
+ double effEle[7];
+ double effPim[7];
+ double rejRatios[7];
+
+ auto he_uncut = d1.Filter(is_electron, {"mcparticles"}).Histo1D({"he_uncut", "P", 7, &binEdges[0]}, "Pthr");
+ auto he_cut = d_ele.Histo1D({"he_cut", "P", 7, &binEdges[0]}, "Pthr");
+ auto hp_uncut = d1.Filter(is_piMinus, {"mcparticles"}).Histo1D({"hp_uncut", "P", 7, &binEdges[0]}, "Pthr");
+ auto hp_cut = d_pim.Histo1D({"hp_cut", "P", 7, &binEdges[0]}, "Pthr");
+
+ double ele_eff_total = (double)*d_ele.Count() / (double)*d1.Filter(is_electron, {"mcparticles"}).Count();
+
+ he_cut->Divide(he_uncut.GetPtr());//Effienciency
+ hp_uncut->Divide(hp_cut.GetPtr());//Rejection power
+
+ for (int i = 1; i < 8; i++){
+ effEle[i - 1] = he_cut->GetBinContent(i);
+ effPim[i - 1] = hp_uncut->GetBinContent(i);
+ rejRatios[i - 1] = effPim[i-1];
+
+ tg->SetPoint(i-1, 0.5*(binEdges[i-1] + binEdges[i]), effPim[i-1]);
+ tg->SetPointError(i-1, 0.5*(binEdges[i] - binEdges[i-1]), hp_uncut ->GetBinError(i));
+ }
+
+ tg->SetTitle(("#pi Rejection Power (#varepsilon_{e} = " + std::to_string(ele_eff_total) + ")").c_str());
+ tg->GetXaxis()->SetTitle("p [GeV]");
+ tg->GetYaxis()->SetTitle("R_{#pi}");
+ tg->SetMarkerColor(kBlue);
+ tg->SetMarkerStyle(20);
+
+ auto cp = new TCanvas("cp", "cp");
+ tg->DrawClone("ap");
+ cp->SaveAs("results/emcal_barrel_pion_rej_RatioRej.png");
+ cp->Clear();
+ auto cp_log = new TCanvas("cp_log", "cp_log");
+ cp_log->SetLogy();
+ cp_log->SetLogx();
+ tg->DrawClone("ap");
+ cp_log->SaveAs("results/emcal_barrel_pion_rej_RatioRej_log.png");
+
+ // Writing out benchmarks
+ //Tests
+ std::string test_tag = "Barrel_emcal_pion_rejection";
+ //TODO: Change test_tag to something else
+ std:string detectorEle = "Barrel_emcal";
+ double target = 1;
+
+ int energyInt = 2;
+ // E 2
+ common_bench::Test pion_rejection_E1{
+ {{"name", fmt::format("{}_E{}", test_tag, energyInt)},
+ {"title", "Pion Rejection"},
+ {"description", fmt::format("Pion rejection with E = {}", energyInt)},
+ {"quantity", "e-/pi-"},
+ {"cut", finalCut},
+ {"e- efficiency", std::to_string(effEle[0])},
+ {"pi- efficiency", std::to_string(effPim[0])},
+ {"target", std::to_string(target)}
+ }
+ };
+ target <= rejRatios[0] ? pion_rejection_E1.pass(rejRatios[0]) : pion_rejection_E1.fail(rejRatios[0]);
+ energyInt = 4;
+ // E 4
+ common_bench::Test pion_rejection_E4{
+ {{"name", fmt::format("{}_E{}", test_tag, energyInt)},
+ {"title", "Pion Rejection"},
+ {"description", fmt::format("Pion rejection with E = {}", energyInt)},
+ {"quantity", "e-/pi-"},
+ {"cut", finalCut},
+ {"e- efficiency", std::to_string(effEle[1])},
+ {"pi- efficiency", std::to_string(effPim[1])},
+ {"target", std::to_string(target)}
+ }
+ };
+ target <= rejRatios[1] ? pion_rejection_E4.pass(rejRatios[1]) : pion_rejection_E4.fail(rejRatios[1]);
+ energyInt = 6;
+ // E 6
+ common_bench::Test pion_rejection_E6{
+ {{"name", fmt::format("{}_E{}", test_tag, energyInt)},
+ {"title", "Pion Rejection"},
+ {"description", fmt::format("Pion rejection with E = {}", energyInt)},
+ {"quantity", "e-/pi-"},
+ {"cut", finalCut},
+ {"e- efficiency", std::to_string(effEle[2])},
+ {"pi- efficiency", std::to_string(effPim[2])},
+ {"target", std::to_string(target)}
+ }
+ };
+ target <= rejRatios[2] ? pion_rejection_E6.pass(rejRatios[2]) : pion_rejection_E6.fail(rejRatios[2]);
+ energyInt = 9;
+ // E 9
+ common_bench::Test pion_rejection_E9{
+ {{"name", fmt::format("{}_E{}", test_tag, energyInt)},
+ {"title", "Pion Rejection"},
+ {"description", fmt::format("Pion rejection with E = {}", energyInt)},
+ {"quantity", "e-/pi-"},
+ {"cut", finalCut},
+ {"e- efficiency", std::to_string(effEle[3])},
+ {"pi- efficiency", std::to_string(effPim[3])},
+ {"target", std::to_string(target)}
+ }
+ };
+ target <= rejRatios[3] ? pion_rejection_E9.pass(rejRatios[3]) : pion_rejection_E9.fail(rejRatios[3]);
+ energyInt = 12;
+ // E 12
+ common_bench::Test pion_rejection_E12{
+ {{"name", fmt::format("{}_E{}", test_tag, energyInt)},
+ {"title", "Pion Rejection"},
+ {"description", fmt::format("Pion rejection with E = {}", energyInt)},
+ {"quantity", "e-/pi-"},
+ {"cut", finalCut},
+ {"e- efficiency", std::to_string(effEle[4])},
+ {"pi- efficiency", std::to_string(effPim[4])},
+ {"target", std::to_string(target)}
+ }
+ };
+ target <= rejRatios[4] ? pion_rejection_E12.pass(rejRatios[4]) : pion_rejection_E12.fail(rejRatios[4]);
+ energyInt = 15;
+ // E 15
+ common_bench::Test pion_rejection_E15{
+ {{"name", fmt::format("{}_E{}", test_tag, energyInt)},
+ {"title", "Pion Rejection"},
+ {"description", fmt::format("Pion rejection with E = {}", energyInt)},
+ {"quantity", "e-/pi-"},
+ {"cut", finalCut},
+ {"e- efficiency", std::to_string(effEle[5])},
+ {"pi- efficiency", std::to_string(effPim[5])},
+ {"target", std::to_string(target)}
+ }
+ };
+ target <= rejRatios[5] ? pion_rejection_E15.pass(rejRatios[5]) : pion_rejection_E15.fail(rejRatios[5]);
+ energyInt = 18;
+ // E 18
+ common_bench::Test pion_rejection_E18{
+ {{"name", fmt::format("{}_E{}", test_tag, energyInt)},
+ {"title", "Pion Rejection"},
+ {"description", fmt::format("Pion rejection with E = {}", energyInt)},
+ {"quantity", "e-/pi-"},
+ {"cut", finalCut},
+ {"e- efficiency", std::to_string(effEle[6])},
+ {"pi- efficiency", std::to_string(effPim[6])},
+ {"target", std::to_string(target)}
+ }
+ };
+ target <= rejRatios[6] ? pion_rejection_E18.pass(rejRatios[6]) : pion_rejection_E18.fail(rejRatios[6]);
+
+ // Writing out all tests
+ common_bench::write_test({pion_rejection_E1,
+ pion_rejection_E4,
+ pion_rejection_E6,
+ pion_rejection_E9,
+ pion_rejection_E12,
+ pion_rejection_E15,
+ pion_rejection_E18},
+ fmt::format("results/{}_pion_rej.json", detectorEle));
+
+}
diff --git a/benchmarks/imaging_ecal/config.yml b/benchmarks/imaging_ecal/config.yml
index 0eb3d818efca2d832cb8f42eea745a82ea9e384a..44e321a0f21bd5c890423d3958fd83125bf2c41c 100644
--- a/benchmarks/imaging_ecal/config.yml
+++ b/benchmarks/imaging_ecal/config.yml
@@ -26,3 +26,23 @@ imaging_ecal_pion0:
script:
- bash benchmarks/imaging_ecal/run_imcal_pion0.sh -t imcal_barrel_pion0 -p "pion0" -n 100
+imaging_ecal_pion_rejection:
+ extends: .rec_benchmark
+ timeout: 48 hours
+ stage: run
+ script:
+ - compile_analyses.py imaging_ecal
+ - bash benchmarks/imaging_ecal/run_emcal_barrel_pion_rej.sh -t emcal_barrel_pion_rej_electron -p "electron" -n 100
+ - bash benchmarks/imaging_ecal/run_emcal_barrel_pion_rej.sh -t emcal_barrel_pion_rej_piminus -p "pion-" -n 100
+
+imaging_ecal_pion_rejection:bench:
+ extends: .rec_benchmark
+ timeout: 1 hours
+ stage: benchmarks2
+ needs:
+ - ["imaging_ecal_pion_rejection"]
+ script:
+ - ls -lhtR
+ - compile_analyses.py imaging_ecal
+ - root -b -q benchmarks/imaging_ecal/analysis/emcal_barrel_pion_rejection_analysis.cxx+
+ #- bash run_pion_rej_analysis.sh
diff --git a/benchmarks/imaging_ecal/options/hybrid_cluster.py b/benchmarks/imaging_ecal/options/hybrid_cluster.py
index 223b61696e2e457084fc9939b6fe51fbcec8f765..f4d3388fbf72d76f9903dd68922d800abab89420 100644
--- a/benchmarks/imaging_ecal/options/hybrid_cluster.py
+++ b/benchmarks/imaging_ecal/options/hybrid_cluster.py
@@ -19,11 +19,13 @@ from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco
with open('config/emcal_barrel_calibration.json') as f:
calib_data = json.load(f)['electron']
-print(calib_data)
+#print(calib_data)
kwargs = dict()
kwargs['img_sf'] = float(calib_data['sampling_fraction_img'])
kwargs['scfi_sf'] = float(calib_data['sampling_fraction_scfi'])
+#kwargs['img_sf'] = 0.0060795
+#kwargs['scfi_sf'] = 0.10343
# input arguments through environment variables
kwargs['input'] = os.environ['CB_EMCAL_SIM_FILE']
@@ -56,14 +58,15 @@ imcaldaq = dict(
pedestalMean=400,
pedestalSigma=50) # 50/32767*3 MeV ~ 5 keV
+
imcaldigi = CalHitDigi('imcal_digi',
- # OutputLevel=DEBUG,
+ OutputLevel=DEBUG,
inputHitCollection='EcalBarrelHits',
outputHitCollection='DigiEcalBarrelImagingHits',
energyResolutions=[0., 0.02, 0.],
**imcaldaq)
imcalreco = ImagingPixelReco('imcal_reco',
- # OutputLevel=DEBUG,
+ #OutputLevel=DEBUG,
inputHitCollection=imcaldigi.outputHitCollection,
outputHitCollection='RecoEcalBarrelImagingHits',
readoutClass='EcalBarrelHits',
@@ -72,7 +75,7 @@ imcalreco = ImagingPixelReco('imcal_reco',
**imcaldaq)
imcalcluster = ImagingTopoCluster('imcal_cluster',
- # OutputLevel=DEBUG,
+ #OutputLevel=DEBUG,
inputHitCollection=imcalreco.outputHitCollection,
outputProtoClusterCollection='EcalBarrelImagingProtoClusters',
localDistXY=[2.*mm, 2*mm],
@@ -80,7 +83,7 @@ imcalcluster = ImagingTopoCluster('imcal_cluster',
neighbourLayersRange=2,
sectorDist=3.*cm)
clusterreco = ImagingClusterReco('imcal_clreco',
- # OutputLevel=DEBUG,
+ #OutputLevel=DEBUG,
inputHitCollection=imcalcluster.inputHitCollection,
inputProtoClusterCollection=imcalcluster.outputProtoClusterCollection,
outputLayerCollection='EcalBarrelImagingClustersLayers',
@@ -112,7 +115,7 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
# merge hits in different layer (projection to local x-y plane)
scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger",
- # OutputLevel=DEBUG,
+ #OutputLevel=DEBUG,
inputHitCollection=scfi_barrel_reco.outputHitCollection,
outputHitCollection="EcalBarrelScFiGridReco",
fields=["fiber"],
@@ -120,7 +123,7 @@ scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger",
readoutClass="EcalBarrelScFiHits")
scfi_barrel_cl = IslandCluster("scfi_barrel_cl",
- # OutputLevel=DEBUG,
+ #OutputLevel=DEBUG,
inputHitCollection=scfi_barrel_merger.outputHitCollection,
outputProtoClusterCollection="EcalBarrelScFiProtoClusters",
splitCluster=False,
diff --git a/benchmarks/imaging_ecal/run_emcal_barrel_pion_rej.sh b/benchmarks/imaging_ecal/run_emcal_barrel_pion_rej.sh
new file mode 100755
index 0000000000000000000000000000000000000000..c4adc4013a3296bada2d724cb2c1f74249e222ef
--- /dev/null
+++ b/benchmarks/imaging_ecal/run_emcal_barrel_pion_rej.sh
@@ -0,0 +1,173 @@
+#!/bin/bash
+
+print_env.sh
+
+function print_the_help {
+ echo "USAGE: ${0} -n <nevents> -t <nametag> -p <particle> "
+ echo " OPTIONS: "
+ echo " -n,--nevents Number of events"
+ echo " -t,--nametag name tag"
+ echo " -p,--particle particle type"
+ echo " allowed types: pion0, pion+, pion-, kaon0, kaon+, kaon-, proton, neutron, electron, positron, photon"
+ exit
+}
+
+POSITIONAL=()
+while [[ $# -gt 0 ]]
+do
+ key="$1"
+
+ case $key in
+ -h|--help)
+ shift # past argument
+ print_the_help
+ ;;
+ -t|--nametag)
+ nametag="$2"
+ shift # past argument
+ shift # past value
+ ;;
+ -p|--particle)
+ particle="$2"
+ shift # past argument
+ shift # past value
+ ;;
+ -n|--nevents)
+ export CB_EMCAL_NUMEV="$2"
+ shift # past argument
+ shift # past value
+ ;;
+ *) # unknown option
+ #POSITIONAL+=("$1") # save it in an array for later
+ echo "unknown option $1"
+ print_the_help
+ shift # past argument
+ ;;
+ esac
+done
+set -- "${POSITIONAL[@]}" # restore positional parameters
+
+export CB_EMCAL_COMPACT_PATH=${DETECTOR_PATH}/${JUGGLER_DETECTOR}.xml
+
+if [[ ! -n "${CB_EMCAL_NUMEV}" ]] ; then
+ export CB_EMCAL_NUMEV=1000
+fi
+
+if [[ ! -n "${CB_EMCAL_IEV}" ]] ; then
+ export CB_EMCAL_IEV="0, 1"
+fi
+
+if [[ ! -n "${CB_EMCAL_ENERGY}" ]] ; then
+ export CB_EMCAL_ENERGY=5.0
+fi
+
+if [[ ! -n "${CB_EMCAL_SAMP_FRAC}" ]] ; then
+ export CB_EMCAL_SAMP_FRAC=0.014
+fi
+
+export CB_EMCAL_NAME_TAG="${nametag}"
+export CB_EMCAL_GEN_FILE="${CB_EMCAL_NAME_TAG}.hepmc"
+
+export CB_EMCAL_SIM_FILE="sim_${CB_EMCAL_NAME_TAG}.root"
+export CB_EMCAL_REC_FILE="rec_${CB_EMCAL_NAME_TAG}.root"
+
+echo "CB_EMCAL_NUMEV = ${CB_EMCAL_NUMEV}"
+echo "CB_EMCAL_COMPACT_PATH = ${CB_EMCAL_COMPACT_PATH}"
+
+# Generate the input events
+python benchmarks/imaging_ecal/scripts/gen_particles.py ${CB_EMCAL_GEN_FILE} -n ${CB_EMCAL_NUMEV}\
+ --angmin 60 --angmax 120 --pmin 0.05 --pmax 18.0 --particles="${particle}"
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running script: generating input events"
+ exit 1
+fi
+
+ls -lh ${CB_EMCAL_GEN_FILE}
+
+# Run geant4 simulations
+npsim --runType batch \
+ -v WARNING \
+ --part.minimalKineticEnergy "0.5*MeV" \
+ --numberOfEvents ${CB_EMCAL_NUMEV} \
+ --compactFile ${CB_EMCAL_COMPACT_PATH} \
+ --inputFiles ${CB_EMCAL_GEN_FILE} \
+ --outputFile ${CB_EMCAL_SIM_FILE}
+
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running npdet"
+ exit 1
+fi
+
+echo "Simulation file"
+rootls -t "${CB_EMCAL_SIM_FILE}"
+
+# Directory for plots
+mkdir -p results
+
+CB_EMCAL_OPTION_DIR=benchmarks/imaging_ecal/options
+# Run Juggler
+xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+ gaudirun.py ${CB_EMCAL_OPTION_DIR}/hybrid_cluster.py
+# gaudirun.py ${CB_EMCAL_OPTION_DIR}/imaging_topocluster.py
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running juggler"
+ exit 1
+fi
+
+# Plot clusters first
+FULL_CAL_SCRIPT_DIR=benchmarks/clustering/scripts
+python ${FULL_CAL_SCRIPT_DIR}/cluster_plots.py ${CB_EMCAL_SIM_FILE} ${CB_EMCAL_REC_FILE} -o results \
+ --collections "EcalBarrelImagingClusters, EcalBarrelScFiClusters"
+
+# check required python modules
+python -m pip install -r benchmarks/imaging_ecal/requirements.txt
+CB_EMCAL_SCRIPT_DIR=benchmarks/imaging_ecal/scripts
+IFS=',' read -ra ADDR <<< "$CB_EMCAL_IEV"
+for iev in "${ADDR[@]}"; do
+ if [[ $iev -ge "${CB_EMCAL_NUMEV}" ]] ; then
+ continue
+ fi
+
+ python ${CB_EMCAL_SCRIPT_DIR}/draw_cluster_layers.py \
+ ${CB_EMCAL_REC_FILE} -e ${iev} --topo-size=1.0 --compact=${CB_EMCAL_COMPACT_PATH} -o results/${particle}
+ if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running analysis script: draw_cluster_layers"
+ exit 1
+ fi
+
+ python ${CB_EMCAL_SCRIPT_DIR}/draw_cluster.py \
+ ${CB_EMCAL_REC_FILE} -e ${iev} --topo-size=2.0 --compact=${CB_EMCAL_COMPACT_PATH} -o results/${particle}
+ if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running analysis script: draw_cluster"
+ exit 1
+ fi
+done
+
+python ${CB_EMCAL_SCRIPT_DIR}/energy_profile.py \
+ ${CB_EMCAL_REC_FILE} --type=EM --energy=${CB_EMCAL_ENERGY} -o results/${particle} \
+ --save=results/profile.csv --color=royalblue
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running analysis script: energy_profile"
+ exit 1
+fi
+
+echo "Reconstruction File"
+rootls -t "${CB_EMCAL_REC_FILE}"
+
+root_filesize=$(stat --format=%s "${CB_EMCAL_REC_FILE}")
+if [[ "${CB_EMCAL_NUMEV}" -lt "500" ]] ; then
+ # file must be less than 10 MB to upload
+ if [[ "${root_filesize}" -lt "10000000" ]] ; then
+ cp ${CB_EMCAL_REC_FILE} results/.
+ fi
+fi
+
+pwd
+echo "Current directory"
+ls
+
+#echo "Results directory"
+#ls benchmarks/imaging_ecal/results/
+
+
+