README.md

@@ -11,21 +11,20 @@ Reconstruction Benchmarks for the EIC
 Here we setup to use our local build of the `juggler` library.
 First set some environment variables.
 ```
-export JUGGLER_INSTALL_PREFIX=$HOME/stow/juggler
+export JUGGLER_INSTALL_PREFIX=/usr/local
 export JUGGLER_DETECTOR=athena   # athena is the default
 export BEAMLINE_CONFIG=ip6       # ip6 is the default
 ```
 
 
+
 ```
-git@eicweb.phy.anl.gov:EIC/benchmarks/reconstruction_benchmarks.git && cd reconstruction_benchmarks
+git clone https://eicweb.phy.anl.gov/EIC/benchmarks/reconstruction_benchmarks.git && cd reconstruction_benchmarks
 git clone https://eicweb.phy.anl.gov/EIC/benchmarks/common_bench.git setup
 source setup/bin/env.sh && ./setup/bin/install_common.sh
 source .local/bin/env.sh && build_detector.sh
 mkdir_local_data_link sim_output
-mkdir -p results
-mkdir -p config
-
+mkdir -p results config
 ```
 
 

benchmarks/imaging_ecal/analysis/emcal_barrel_pion_rejection_analysis.cxx

+////////////////////////////////////////
+// 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));
+
+}

benchmarks/imaging_ecal/config.yml

@@ -26,3 +26,24 @@ imaging_ecal_pion0:
   script:
     - bash benchmarks/imaging_ecal/run_imcal_pion0.sh -t imcal_barrel_pion0 -p "pion0" -n 100
 
+## SJJ removed as these are causing the CI to fail
+#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

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,

benchmarks/imaging_ecal/run_emcal_barrel_pion_rej.sh

+#!/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/
+
+
+

benchmarks/imaging_ecal/scripts/draw_cluster_layers.py

@@ -110,7 +110,7 @@ if __name__ == '__main__':
     else:
         dfc = df[(df['eta'] <= df['eta'].quantile(0.95)) & (df['eta'] >= df['eta'].quantile(0.05)) &
                  (df['phi'] <= df['phi'].quantile(0.95)) & (df['phi'] >= df['phi'].quantile(0.05))]
-        vec = np.average(dfc[['x', 'y', 'z']].values, axis=0, weights=dfc['energy'].values)
+        vec = np.average(dfc[['x', 'y', 'z']].values, axis=0, weights= dfc['energy'].values if np.sum(dfc['energy'].values) > 0 else None)
     vec = vec/np.linalg.norm(vec)
 
     # particle line from (0, 0, 0) to the inner Ecal surface

benchmarks/track_finding/options/track_reconstruction.py

@@ -5,11 +5,17 @@ from GaudiKernel import SystemOfUnits as units
 
 detector_name = "topside"
 if "JUGGLER_DETECTOR" in os.environ :
-  detector_name = str(os.environ["JUGGLER_DETECTOR"])
+    detector_name = str(os.environ["JUGGLER_DETECTOR"])
 
 detector_path = ""
 if "DETECTOR_PATH" in os.environ :
-  detector_path = str(os.environ["DETECTOR_PATH"])
+    detector_path = str(os.environ["DETECTOR_PATH"])
+
+detector_version = 'default'
+if "JUGGLER_DETECTOR_VERSION" in os.environ:
+    env_version = str(os.environ["JUGGLER_DETECTOR_VERSION"])
+    if 'acadia' in env_version:
+        detector_version = 'acadia'
 
 # todo add checks
 input_sim_file  = str(os.environ["JUGGLER_SIM_FILE"])
@@ -40,8 +46,13 @@ from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
 
 algorithms = [ ]
 
+input_collections = ['mcparticles','TrackerEndcapHits','TrackerBarrelHits','VertexBarrelHits','GEMTrackerEndcapHits']
+if 'acadia' in detector_version:
+    input_collections.append('VertexEndcapHits')
+else:
+    input_collections.append('MPGDTrackerBarrelHits')
 podioinput = PodioInput("PodioReader", 
-                        collections=["mcparticles","TrackerEndcapHits","TrackerBarrelHits","VertexBarrelHits","VertexEndcapHits","GEMTrackerEndcapHits"])
+                        collections=input_collections)#, OutputLevel=DEBUG)
 algorithms.append( podioinput )
 
 trk_b_digi = TrackerDigi("trk_b_digi", 
@@ -61,11 +72,18 @@ vtx_b_digi = TrackerDigi("vtx_b_digi",
         timeResolution=8)
 algorithms.append( vtx_b_digi )
 
-vtx_ec_digi = TrackerDigi("vtx_ec_digi", 
-        inputHitCollection="VertexEndcapHits",
-        outputHitCollection="VertexEndcapRawHits",
-        timeResolution=8)
-algorithms.append( vtx_ec_digi )
+if 'acadia' in detector_version:
+    vtx_ec_digi = TrackerDigi("vtx_ec_digi", 
+            inputHitCollection="VertexEndcapHits",
+            outputHitCollection="VertexEndcapRawHits",
+            timeResolution=8)
+    algorithms.append( vtx_ec_digi )
+else:
+    mm_b_digi = TrackerDigi("mm_b_digi", 
+            inputHitCollection="MPGDTrackerBarrelHits",
+            outputHitCollection="MPGDTrackerBarrelRawHits",
+            timeResolution=8)
+    algorithms.append( mm_b_digi )
 
 gem_ec_digi = TrackerDigi("gem_ec_digi",
         inputHitCollection="GEMTrackerEndcapHits",
@@ -89,23 +107,35 @@ vtx_b_reco = TrackerHitReconstruction("vtx_b_reco",
         outputHitCollection="VertexBarrelRecHits")
 algorithms.append( vtx_b_reco )
 
-vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
-        inputHitCollection = vtx_ec_digi.outputHitCollection,
-        outputHitCollection="VertexEndcapRecHits")
-algorithms.append( vtx_ec_reco )
+if 'acadia' in detector_version:
+    vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
+            inputHitCollection = vtx_ec_digi.outputHitCollection,
+            outputHitCollection="VertexEndcapRecHits")
+    algorithms.append( vtx_ec_reco )
+else:
+    mm_b_reco = TrackerHitReconstruction("mm_b_reco",
+            inputHitCollection = mm_b_digi.outputHitCollection,
+            outputHitCollection="MPGDTrackerBarrelRecHits")
+    algorithms.append( mm_b_reco )
+
 
 gem_ec_reco = TrackerHitReconstruction("gem_ec_reco",
         inputHitCollection=gem_ec_digi.outputHitCollection,
         outputHitCollection="GEMTrackerEndcapRecHits")
 algorithms.append(gem_ec_reco)
 
+input_tracking_hits = [
+    str(trk_b_reco.outputHitCollection),
+    str(trk_ec_reco.outputHitCollection),
+    str(vtx_b_reco.outputHitCollection),
+    str(gem_ec_reco.outputHitCollection) ]
+if 'acadia' in detector_version:
+    input_tracking_hits.append(str(vtx_ec_reco.outputHitCollection))
+else:
+    input_tracking_hits.append(str(mm_b_reco.outputHitCollection))
+
 trk_hit_col = TrackingHitsCollector("trk_hit_col",
-        inputTrackingHits=[
-            str(trk_b_reco.outputHitCollection),
-            str(trk_ec_reco.outputHitCollection),
-            str(vtx_b_reco.outputHitCollection),
-            str(vtx_ec_reco.outputHitCollection),
-            str(gem_ec_reco.outputHitCollection) ],
+        inputTrackingHits=input_tracking_hits,
         trackingHits="trackingHits",
         OutputLevel=DEBUG)
 algorithms.append( trk_hit_col )

benchmarks/track_fitting/options/track_reconstruction.py

@@ -5,11 +5,17 @@ from GaudiKernel import SystemOfUnits as units
 
 detector_name = "athena"
 if "JUGGLER_DETECTOR" in os.environ :
-  detector_name = str(os.environ["JUGGLER_DETECTOR"])
+    detector_name = str(os.environ["JUGGLER_DETECTOR"])
 
 detector_path = ""
 if "DETECTOR_PATH" in os.environ :
-  detector_path = str(os.environ["DETECTOR_PATH"])
+    detector_path = str(os.environ["DETECTOR_PATH"])
+
+detector_version = 'default'
+if "JUGGLER_DETECTOR_VERSION" in os.environ:
+    env_version = str(os.environ["JUGGLER_DETECTOR_VERSION"])
+    if 'acadia' in env_version:
+        detector_version = 'acadia'
 
 # todo add checks
 input_sim_file  = str(os.environ["JUGGLER_SIM_FILE"])
@@ -40,8 +46,13 @@ from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
 
 algorithms = [ ]
 
+input_collections = ['mcparticles','TrackerEndcapHits','TrackerBarrelHits','VertexBarrelHits','GEMTrackerEndcapHits']
+if 'acadia' in detector_version:
+    input_collections.append('VertexEndcapHits')
+else:
+    input_collections.append('MPGDTrackerBarrelHits')
 podioinput = PodioInput("PodioReader", 
-                        collections=["mcparticles","TrackerEndcapHits","TrackerBarrelHits","VertexBarrelHits","VertexEndcapHits","GEMTrackerEndcapHits"])#, OutputLevel=DEBUG)
+                        collections=input_collections)#, OutputLevel=DEBUG)
 algorithms.append( podioinput )
 
 trk_b_digi = TrackerDigi("trk_b_digi", 
@@ -61,11 +72,18 @@ vtx_b_digi = TrackerDigi("vtx_b_digi",
         timeResolution=8)
 algorithms.append( vtx_b_digi )
 
-vtx_ec_digi = TrackerDigi("vtx_ec_digi", 
-        inputHitCollection="VertexEndcapHits",
-        outputHitCollection="VertexEndcapRawHits",
-        timeResolution=8)
-algorithms.append( vtx_ec_digi )
+if 'acadia' in detector_version:
+    vtx_ec_digi = TrackerDigi("vtx_ec_digi", 
+            inputHitCollection="VertexEndcapHits",
+            outputHitCollection="VertexEndcapRawHits",
+            timeResolution=8)
+    algorithms.append( vtx_ec_digi )
+else:
+    mm_b_digi = TrackerDigi("mm_b_digi", 
+            inputHitCollection="MPGDTrackerBarrelHits",
+            outputHitCollection="MPGDTrackerBarrelRawHits",
+            timeResolution=8)
+    algorithms.append( mm_b_digi )
 
 gem_ec_digi = TrackerDigi("gem_ec_digi",
         inputHitCollection="GEMTrackerEndcapHits",
@@ -89,24 +107,36 @@ vtx_b_reco = TrackerHitReconstruction("vtx_b_reco",
         outputHitCollection="VertexBarrelRecHits")
 algorithms.append( vtx_b_reco )
 
-vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
-        inputHitCollection = vtx_ec_digi.outputHitCollection,
-        outputHitCollection="VertexEndcapRecHits")
-algorithms.append( vtx_ec_reco )
+if 'acadia' in detector_version:
+    vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
+            inputHitCollection = vtx_ec_digi.outputHitCollection,
+            outputHitCollection="VertexEndcapRecHits")
+    algorithms.append( vtx_ec_reco )
+else:
+    mm_b_reco = TrackerHitReconstruction("mm_b_reco",
+            inputHitCollection = mm_b_digi.outputHitCollection,
+            outputHitCollection="MPGDTrackerBarrelRecHits")
+    algorithms.append( mm_b_reco )
 
 gem_ec_reco = TrackerHitReconstruction("gem_ec_reco",
         inputHitCollection=gem_ec_digi.outputHitCollection,
         outputHitCollection="GEMTrackerEndcapRecHits")
 algorithms.append(gem_ec_reco)
 
+input_tracking_hits = [
+    str(trk_b_reco.outputHitCollection),
+    str(trk_ec_reco.outputHitCollection),
+    str(vtx_b_reco.outputHitCollection),
+    str(gem_ec_reco.outputHitCollection) ]
+if 'acadia' in detector_version:
+    input_tracking_hits.append(str(vtx_ec_reco.outputHitCollection))
+else:
+    input_tracking_hits.append(str(mm_b_reco.outputHitCollection))
+
 trk_hit_col = TrackingHitsCollector("trk_hit_col",
-        inputTrackingHits=[
-            str(vtx_b_reco.outputHitCollection),
-            str(vtx_ec_reco.outputHitCollection),
-            str(trk_b_reco.outputHitCollection),
-            str(trk_ec_reco.outputHitCollection),
-            str(gem_ec_reco.outputHitCollection) ],
-        trackingHits="trackingHits")
+        inputTrackingHits=input_tracking_hits,
+        trackingHits="trackingHits",
+        OutputLevel=DEBUG)
 algorithms.append( trk_hit_col )
 
 # Hit Source linker 

benchmarks/tracking/options/track_reconstruction.py

@@ -5,18 +5,36 @@ from GaudiKernel import SystemOfUnits as units
 
 detector_name = "topside"
 if "JUGGLER_DETECTOR" in os.environ :
-  detector_name = str(os.environ["JUGGLER_DETECTOR"])
+    detector_name = str(os.environ["JUGGLER_DETECTOR"])
 
 detector_path = ""
 if "DETECTOR_PATH" in os.environ :
-  detector_path = str(os.environ["DETECTOR_PATH"])
+    detector_path = str(os.environ["DETECTOR_PATH"])
+
+detector_version = 'default'
+if "JUGGLER_DETECTOR_VERSION" in os.environ:
+    env_version = str(os.environ["JUGGLER_DETECTOR_VERSION"])
+    if 'acadia' in env_version:
+        detector_version = 'acadia'
 
 # todo add checks
 input_sim_file  = str(os.environ["JUGGLER_SIM_FILE"])
 output_rec_file = str(os.environ["JUGGLER_REC_FILE"])
 n_events = str(os.environ["JUGGLER_N_EVENTS"])
 
-geo_service  = GeoSvc("GeoSvc", detectors=["{}/{}.xml".format(detector_path,detector_name)], OutputLevel=WARNING)
+## note: old version of material map is called material-maps.XXX, new version is materials-map.XXX
+##       these names are somewhat inconsistent, and should probably all be renamed to 'material-map.XXX' 
+##       FIXME
+if detector_version == 'acadia':
+    geo_service = GeoSvc("GeoSvc",
+                         detectors=["{}/{}.xml".format(detector_path,detector_name)],
+                         materials="config/material-maps.json",
+                         OutputLevel=WARNING)
+else:
+    geo_service = GeoSvc("GeoSvc",
+                         detectors=["{}/{}.xml".format(detector_path,detector_name)],
+                         materials="calibrations/materials-map.cbor",
+                         OutputLevel=WARNING)
 podioevent   = EICDataSvc("EventDataSvc", inputs=[input_sim_file], OutputLevel=WARNING)
 
 from Configurables import PodioInput
@@ -33,15 +51,20 @@ from Configurables import Jug__Reco__TrackParamVertexClusterInit as TrackParamVe
 
 from Configurables import Jug__Reco__TrackFindingAlgorithm as TrackFindingAlgorithm
 from Configurables import Jug__Reco__ParticlesFromTrackFit as ParticlesFromTrackFit
-from Configurables import Jug__Reco__TrajectoryFromTrackFit as TrajectoryFromTrackFit
+#from Configurables import Jug__Reco__TrajectoryFromTrackFit as TrajectoryFromTrackFit
 
 from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
 
 
 algorithms = [ ]
 
+input_collections = ['mcparticles','TrackerEndcapHits','TrackerBarrelHits','VertexBarrelHits','GEMTrackerEndcapHits']
+if 'acadia' in detector_version:
+    input_collections.append('VertexEndcapHits')
+else:
+    input_collections.append('MPGDTrackerBarrelHits')
 podioinput = PodioInput("PodioReader", 
-                        collections=["mcparticles","TrackerEndcapHits","TrackerBarrelHits","VertexBarrelHits","VertexEndcapHits","GEMTrackerEndcapHits"])#, OutputLevel=DEBUG)
+                        collections=input_collections)#, OutputLevel=DEBUG)
 algorithms.append( podioinput )
 
 trk_b_digi = TrackerDigi("trk_b_digi", 
@@ -61,11 +84,18 @@ vtx_b_digi = TrackerDigi("vtx_b_digi",
         timeResolution=8)
 algorithms.append( vtx_b_digi )
 
-vtx_ec_digi = TrackerDigi("vtx_ec_digi", 
-        inputHitCollection="VertexEndcapHits",
-        outputHitCollection="VertexEndcapRawHits",
-        timeResolution=8)
-algorithms.append( vtx_ec_digi )
+if 'acadia' in detector_version:
+    vtx_ec_digi = TrackerDigi("vtx_ec_digi", 
+            inputHitCollection="VertexEndcapHits",
+            outputHitCollection="VertexEndcapRawHits",
+            timeResolution=8)
+    algorithms.append( vtx_ec_digi )
+else:
+    mm_b_digi = TrackerDigi("mm_b_digi", 
+            inputHitCollection="MPGDTrackerBarrelHits",
+            outputHitCollection="MPGDTrackerBarrelRawHits",
+            timeResolution=8)
+    algorithms.append( mm_b_digi )
 
 gem_ec_digi = TrackerDigi("gem_ec_digi",
         inputHitCollection="GEMTrackerEndcapHits",
@@ -89,25 +119,35 @@ vtx_b_reco = TrackerHitReconstruction("vtx_b_reco",
         outputHitCollection="VertexBarrelRecHits")
 algorithms.append( vtx_b_reco )
 
-vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
-        inputHitCollection = vtx_ec_digi.outputHitCollection,
-        outputHitCollection="VertexEndcapRecHits")
-algorithms.append( vtx_ec_reco )
+if 'acadia' in detector_version:
+    vtx_ec_reco = TrackerHitReconstruction("vtx_ec_reco",
+            inputHitCollection = vtx_ec_digi.outputHitCollection,
+            outputHitCollection="VertexEndcapRecHits")
+    algorithms.append( vtx_ec_reco )
+else:
+    mm_b_reco = TrackerHitReconstruction("mm_b_reco",
+            inputHitCollection = mm_b_digi.outputHitCollection,
+            outputHitCollection="MPGDTrackerBarrelRecHits")
+    algorithms.append( mm_b_reco )
 
 gem_ec_reco = TrackerHitReconstruction("gem_ec_reco",
         inputHitCollection=gem_ec_digi.outputHitCollection,
         outputHitCollection="GEMTrackerEndcapRecHits")
 algorithms.append(gem_ec_reco)
 
+input_tracking_hits = [
+    str(trk_b_reco.outputHitCollection),
+    str(trk_ec_reco.outputHitCollection),
+    str(vtx_b_reco.outputHitCollection),
+    str(gem_ec_reco.outputHitCollection) ]
+if 'acadia' in detector_version:
+    input_tracking_hits.append(str(vtx_ec_reco.outputHitCollection))
+else:
+    input_tracking_hits.append(str(mm_b_reco.outputHitCollection))
+
 trk_hit_col = TrackingHitsCollector("trk_hit_col",
-        inputTrackingHits=[
-            str(trk_b_reco.outputHitCollection),
-            str(trk_ec_reco.outputHitCollection),
-            str(vtx_b_reco.outputHitCollection),
-            str(vtx_ec_reco.outputHitCollection),
-            str(gem_ec_reco.outputHitCollection) ],
-        trackingHits="trackingHits",
-        OutputLevel=DEBUG)
+        inputTrackingHits=input_tracking_hits,
+        trackingHits="trackingHits")
 algorithms.append( trk_hit_col )
 
 # Hit Source linker 
@@ -115,23 +155,21 @@ sourcelinker = TrackerSourceLinker("sourcelinker",
         inputHitCollection=trk_hit_col.trackingHits,
         outputSourceLinks="TrackSourceLinks",
         outputMeasurements="TrackMeasurements",
-        OutputLevel=DEBUG)
+        OutputLevel=VERBOSE)
 algorithms.append( sourcelinker )
 
 ## Track param init
 truth_trk_init = TrackParamTruthInit("truth_trk_init",
         inputMCParticles="mcparticles",
         outputInitialTrackParameters="InitTrackParams")
-        #OutputLevel=DEBUG)
 algorithms.append( truth_trk_init )
 
 # Tracking algorithms
 trk_find_alg = TrackFindingAlgorithm("trk_find_alg",
         inputSourceLinks = sourcelinker.outputSourceLinks,
         inputMeasurements = sourcelinker.outputMeasurements,
-        inputInitialTrackParameters= "InitTrackParams",#"InitTrackParamsFromClusters", 
+        inputInitialTrackParameters= "InitTrackParams",
         outputTrajectories="trajectories")
-        #OutputLevel=DEBUG)
 algorithms.append( trk_find_alg )
 
 parts_from_fit = ParticlesFromTrackFit("parts_from_fit",
@@ -141,29 +179,11 @@ parts_from_fit = ParticlesFromTrackFit("parts_from_fit",
         #OutputLevel=DEBUG)
 algorithms.append( parts_from_fit )
 
-trajs_from_fit = TrajectoryFromTrackFit("trajs_from_fit",
-inputTrajectories = trk_find_alg.outputTrajectories,
-outputTrajectoryParameters = "outputTrajectoryParameters")
+#trajs_from_fit = TrajectoryFromTrackFit("trajs_from_fit",
+#inputTrajectories = trk_find_alg.outputTrajectories,
+#outputTrajectoryParameters = "outputTrajectoryParameters")
         #OutputLevel=DEBUG)
-algorithms.append(trajs_from_fit)
-
-#types = []
-## this printout is useful to check that the type information is passed to python correctly
-#print("---------------------------------------\n")
-#print("---\n# List of input and output types by class")
-#for configurable in sorted([ PodioInput, EICDataSvc, PodioOutput,
-#                             TrackerHitReconstruction,ExampleCaloDigi, 
-#                             UFSDTrackerDigi, TrackerSourceLinker,
-#                             PodioOutput],
-#                           key=lambda c: c.getType()):
-#    print("\"{}\":".format(configurable.getType()))
-#    props = configurable.getDefaultProperties()
-#    for propname, prop in sorted(props.items()):
-#        print(" prop name: {}".format(propname))
-#        if isinstance(prop, DataHandleBase):
-#            types.append(prop.type())
-#            print("  {}: \"{}\"".format(propname, prop.type()))
-#print("---")
+#algorithms.append(trajs_from_fit)
 
 out = PodioOutput("out", filename=output_rec_file)
 out.outputCommands = ["keep *", 
@@ -183,4 +203,3 @@ ApplicationMgr(
     ExtSvc = [podioevent,geo_service],
     OutputLevel=WARNING
  )
-

benchmarks/tracking/scripts/tracking_performance.py

@@ -186,7 +186,7 @@ if __name__ == '__main__':
                              weights=np.repeat(1./float(rec_p.shape[0]), rec_p.shape[0]), ec='k')
     nbins = hbins.shape[0] - 1
     ax.set_ylabel('Normalized Counts / {:d} Bins'.format(nbins), fontsize=20)
-    ax.set_xlabel('$\delta p$ (%)', fontsize=20)
+    ax.set_xlabel(r'$\delta$ p (%)', fontsize=20)
 
     # theta resolution
     ax = axs.flat[2]

bin/get_calibrations

@@ -4,7 +4,7 @@ branch=${1:-master}
 detector_benchmarks=https://eicweb.phy.anl.gov/EIC/benchmarks/detector_benchmarks/-/jobs/artifacts/${branch}/raw/
 
 mkdir -p config
-for i in results/emcal_barrel_calibration.json ; do
+for i in results/emcal_barrel_calibration.json results/material-maps.json ; do
   curl --fail -sL ${detector_benchmarks}/${i}?job=deploy_results --output config/$(basename ${i})
   if [[ "$?" -ne "0" ]] ; then
     if find ${DETECTOR_PATH} -name $(basename $i) ; then