Skip to content
Snippets Groups Projects
Commit 6a53b574 authored by Marshall Scott's avatar Marshall Scott Committed by Jihee Kim
Browse files

Resolve "barrel ecal pi minus benchmark"

parent 5e26fad9
No related tags found
1 merge request!129Resolve "barrel ecal pi minus benchmark"
Hide whitespace changes
Inline Side-by-side
benchmarks/imaging_ecal/analysis/emcal_barrel_pion_rejection_analysis.cxx 0 → 100644
+ 946
− 0
View file @ 6a53b574
////////////////////////////////////////
// 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
+ 20
− 0
View file @ 6a53b574
...@@ -26,3 +26,23 @@ imaging_ecal_pion0: ...@@ -26,3 +26,23 @@ imaging_ecal_pion0:
script: script:
- bash benchmarks/imaging_ecal/run_imcal_pion0.sh -t imcal_barrel_pion0 -p "pion0" -n 100 - 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
benchmarks/imaging_ecal/options/hybrid_cluster.py
+ 10
− 7
View file @ 6a53b574
...@@ -19,11 +19,13 @@ from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco ...@@ -19,11 +19,13 @@ from Configurables import Jug__Reco__ImagingClusterReco as ImagingClusterReco
with open('config/emcal_barrel_calibration.json') as f: with open('config/emcal_barrel_calibration.json') as f:
calib_data = json.load(f)['electron'] calib_data = json.load(f)['electron']
print(calib_data) #print(calib_data)
kwargs = dict() kwargs = dict()
kwargs['img_sf'] = float(calib_data['sampling_fraction_img']) kwargs['img_sf'] = float(calib_data['sampling_fraction_img'])
kwargs['scfi_sf'] = float(calib_data['sampling_fraction_scfi']) kwargs['scfi_sf'] = float(calib_data['sampling_fraction_scfi'])
#kwargs['img_sf'] = 0.0060795
#kwargs['scfi_sf'] = 0.10343
# input arguments through environment variables # input arguments through environment variables
kwargs['input'] = os.environ['CB_EMCAL_SIM_FILE'] kwargs['input'] = os.environ['CB_EMCAL_SIM_FILE']
...@@ -56,14 +58,15 @@ imcaldaq = dict( ...@@ -56,14 +58,15 @@ imcaldaq = dict(
pedestalMean=400, pedestalMean=400,
pedestalSigma=50) # 50/32767*3 MeV ~ 5 keV pedestalSigma=50) # 50/32767*3 MeV ~ 5 keV
imcaldigi = CalHitDigi('imcal_digi', imcaldigi = CalHitDigi('imcal_digi',
# OutputLevel=DEBUG, OutputLevel=DEBUG,
inputHitCollection='EcalBarrelHits', inputHitCollection='EcalBarrelHits',
outputHitCollection='DigiEcalBarrelImagingHits', outputHitCollection='DigiEcalBarrelImagingHits',
energyResolutions=[0., 0.02, 0.], energyResolutions=[0., 0.02, 0.],
**imcaldaq) **imcaldaq)
imcalreco = ImagingPixelReco('imcal_reco', imcalreco = ImagingPixelReco('imcal_reco',
# OutputLevel=DEBUG, #OutputLevel=DEBUG,
inputHitCollection=imcaldigi.outputHitCollection, inputHitCollection=imcaldigi.outputHitCollection,
outputHitCollection='RecoEcalBarrelImagingHits', outputHitCollection='RecoEcalBarrelImagingHits',
readoutClass='EcalBarrelHits', readoutClass='EcalBarrelHits',
...@@ -72,7 +75,7 @@ imcalreco = ImagingPixelReco('imcal_reco', ...@@ -72,7 +75,7 @@ imcalreco = ImagingPixelReco('imcal_reco',
**imcaldaq) **imcaldaq)
imcalcluster = ImagingTopoCluster('imcal_cluster', imcalcluster = ImagingTopoCluster('imcal_cluster',
# OutputLevel=DEBUG, #OutputLevel=DEBUG,
inputHitCollection=imcalreco.outputHitCollection, inputHitCollection=imcalreco.outputHitCollection,
outputProtoClusterCollection='EcalBarrelImagingProtoClusters', outputProtoClusterCollection='EcalBarrelImagingProtoClusters',
localDistXY=[2.*mm, 2*mm], localDistXY=[2.*mm, 2*mm],
...@@ -80,7 +83,7 @@ imcalcluster = ImagingTopoCluster('imcal_cluster', ...@@ -80,7 +83,7 @@ imcalcluster = ImagingTopoCluster('imcal_cluster',
neighbourLayersRange=2, neighbourLayersRange=2,
sectorDist=3.*cm) sectorDist=3.*cm)
clusterreco = ImagingClusterReco('imcal_clreco', clusterreco = ImagingClusterReco('imcal_clreco',
# OutputLevel=DEBUG, #OutputLevel=DEBUG,
inputHitCollection=imcalcluster.inputHitCollection, inputHitCollection=imcalcluster.inputHitCollection,
inputProtoClusterCollection=imcalcluster.outputProtoClusterCollection, inputProtoClusterCollection=imcalcluster.outputProtoClusterCollection,
outputLayerCollection='EcalBarrelImagingClustersLayers', outputLayerCollection='EcalBarrelImagingClustersLayers',
...@@ -112,7 +115,7 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco", ...@@ -112,7 +115,7 @@ scfi_barrel_reco = CalHitReco("scfi_barrel_reco",
# merge hits in different layer (projection to local x-y plane) # merge hits in different layer (projection to local x-y plane)
scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger", scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger",
# OutputLevel=DEBUG, #OutputLevel=DEBUG,
inputHitCollection=scfi_barrel_reco.outputHitCollection, inputHitCollection=scfi_barrel_reco.outputHitCollection,
outputHitCollection="EcalBarrelScFiGridReco", outputHitCollection="EcalBarrelScFiGridReco",
fields=["fiber"], fields=["fiber"],
...@@ -120,7 +123,7 @@ scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger", ...@@ -120,7 +123,7 @@ scfi_barrel_merger = CalHitsMerger("scfi_barrel_merger",
readoutClass="EcalBarrelScFiHits") readoutClass="EcalBarrelScFiHits")
scfi_barrel_cl = IslandCluster("scfi_barrel_cl", scfi_barrel_cl = IslandCluster("scfi_barrel_cl",
# OutputLevel=DEBUG, #OutputLevel=DEBUG,
inputHitCollection=scfi_barrel_merger.outputHitCollection, inputHitCollection=scfi_barrel_merger.outputHitCollection,
outputProtoClusterCollection="EcalBarrelScFiProtoClusters", outputProtoClusterCollection="EcalBarrelScFiProtoClusters",
splitCluster=False, splitCluster=False,
......
benchmarks/imaging_ecal/run_emcal_barrel_pion_rej.sh 0 → 100755
+ 173
− 0
View file @ 6a53b574
#!/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/
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment