Skip to content
Snippets Groups Projects

Simple e/pi separation

Closed Whitney Armstrong requested to merge e_pi_sep into master
Compare and
3 files
+ 215
6
Compare changes
  • Side-by-side
  • Inline
Files
3
 
////////////////////////////////////////
 
// Read reconstruction ROOT output file
 
// Plot variables
 
////////////////////////////////////////
 
 
#include "ROOT/RDataFrame.hxx"
 
#include <iostream>
 
 
#include "dd4pod/Geant4ParticleCollection.h"
 
#include "dd4pod/CalorimeterHitCollection.h"
 
 
#include "TCanvas.h"
 
#include "TStyle.h"
 
#include "TMath.h"
 
#include "TH1.h"
 
#include "TF1.h"
 
#include "TH1D.h"
 
 
R__LOAD_LIBRARY(libfmt.so)
 
#include "fmt/core.h"
 
#include "DD4hep/Detector.h"
 
#include "DDG4/Geant4Data.h"
 
#include "DDRec/CellIDPositionConverter.h"
 
//#include "DDRec/SurfaceManager.h"
 
//#include "DDRec/Surface.h"
 
 
using ROOT::RDataFrame;
 
using namespace ROOT::VecOps;
 
 
void e_pi_separation(const char* input_fname =
 
"sim_output/sim_emcal_barrel_uniform_electrons.root",
 
const char* input_fname2 =
 
"sim_output/sim_emcal_barrel_uniform_pions.root") {
 
//// 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_fname,input_fname2});
 
 
std::string detector_path = "";
 
std::string detector_name = "topside";
 
if(std::getenv("DETECTOR_PATH")) {
 
detector_path = std::getenv("DETECTOR_PATH");
 
}
 
if(std::getenv("JUGGLER_DETECTOR")) {
 
detector_name = std::getenv("JUGGLER_DETECTOR");
 
}
 
 
dd4hep::Detector& detector = dd4hep::Detector::getInstance();
 
detector.fromCompact(fmt::format("{}/{}.xml", detector_path,detector_name));
 
//dd4hep::rec::CellIDPositionConverter cellid_converter(detector);
 
 
auto decoder = detector.readout("EcalBarrelHits").idSpec().decoder();
 
fmt::print("{}\n", decoder->fieldDescription());
 
auto layer_index = decoder->index("layer");
 
fmt::print(" layer index is {}.\n", layer_index);
 
 
 
// Thrown Energy [GeV]
 
auto is_electron = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
 
auto p = input[2];
 
if(p.pdgID == 11) return true;
 
return false;
 
};
 
auto Ethr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
 
auto p = input[2];
 
auto energy = TMath::Sqrt(p.psx * p.psx + p.psy * p.psy + p.psz * p.psz + p.mass * p.mass);
 
return energy;
 
};
 
 
// Number of hits
 
auto nhits = [] (const std::vector<dd4pod::CalorimeterHitData>& evt) {return (int) evt.size(); };
 
 
// Energy deposition [GeV]
 
auto Esim = [](const std::vector<dd4pod::CalorimeterHitData>& evt) {
 
auto total_edep = 0.0;
 
for (const auto& i: evt)
 
total_edep += i.energyDeposit;
 
return total_edep;
 
};
 
auto Esim_front = [=](const std::vector<dd4pod::CalorimeterHitData>& evt) {
 
auto total_edep = 0.0;
 
for (const auto& i: evt) {
 
fmt::print("cell id {}, layer {}\n",i.cellID, decoder->get(i.cellID, layer_index));
 
if( decoder->get(i.cellID, layer_index) < 4 ){
 
total_edep += i.energyDeposit;
 
}
 
}
 
return total_edep;
 
};
 
 
// Sampling fraction = Esampling / Ethrown
 
auto fsam = [](const double sampled, const double thrown) {
 
return sampled / thrown;
 
};
 
 
// Define variables
 
auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"})
 
.Define("is_electron", is_electron, {"mcparticles"})
 
.Define("nhits", nhits, {"EcalBarrelHits"})
 
.Define("Esim", Esim, {"EcalBarrelHits"})
 
.Define("Esim_front", Esim_front, {"EcalBarrelHits"})
 
.Define("Eratio", "Esim_front/Esim")
 
.Define("fsam", fsam, {"Esim", "Ethr"});
 
auto d2 = d1.Filter("is_electron");
 
 
// Define Histograms
 
auto hEthr = d1.Histo1D(
 
{"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 7.5},
 
"Ethr");
 
auto hNhits =
 
d1.Histo1D({"hNhits", "Number of hits per events; Number of hits; Events",
 
100, 0.0, 2000.0},
 
"nhits");
 
auto hEsim = d1.Histo1D(
 
{"hEsim", "Energy Deposit; Energy Deposit [GeV]; Events", 100, 0.0, 0.2},
 
"Esim");
 
 
auto hEsim_front = d1.Histo1D(
 
{"hEsim_front", "; Front Energy Deposit [GeV]; Events", 100, 0.0, 0.2},
 
"Esim_front");
 
 
auto hfsam = d1.Histo1D(
 
{"hfsam", "Sampling Fraction; Sampling Fraction; Events", 100, 0.0, 0.1},
 
"fsam");
 
 
auto hEratio = d1.Histo1D(
 
{"Eratio", ";E_front/E_tot; Events", 100, 0.0, 0.4},
 
"Eratio");
 
auto hEratio2 = d2.Histo1D(
 
{"Eratio2", ";E_front/E_tot; Events", 100, 0.0, 0.4},
 
"Eratio");
 
 
// Event Counts
 
auto nevents_thrown = d1.Count();
 
std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
 
 
// Draw Histograms
 
{
 
TCanvas* c1 = new TCanvas("c1", "c1", 700, 500);
 
c1->SetLogy(1);
 
auto h = hEthr->DrawCopy();
 
//h->GetYaxis()->SetTitleOffset(1.4);
 
h->SetLineWidth(2);
 
h->SetLineColor(kBlue);
 
c1->SaveAs("results/emcal_barrel_electrons_Ethr.png");
 
c1->SaveAs("results/emcal_barrel_electrons_Ethr.pdf");
 
}
 
std::cout << "derp1\n";
 
 
{
 
TCanvas* c2 = new TCanvas("c2", "c2", 700, 500);
 
c2->SetLogy(1);
 
auto h = hNhits->DrawCopy();
 
//h->GetYaxis()->SetTitleOffset(1.4);
 
h->SetLineWidth(2);
 
h->SetLineColor(kBlue);
 
c2->SaveAs("results/emcal_barrel_electrons_nhits.png");
 
c2->SaveAs("results/emcal_barrel_electrons_nhits.pdf");
 
}
 
 
{
 
TCanvas* c3 = new TCanvas("c3", "c3", 700, 500);
 
c3->SetLogy(1);
 
auto h = hEsim->DrawCopy();
 
//h->GetYaxis()->SetTitleOffset(1.4);
 
h->SetLineWidth(2);
 
h->SetLineColor(kBlue);
 
h = hEsim_front->DrawCopy("same");
 
h->SetLineWidth(2);
 
h->SetLineColor(2);
 
c3->SaveAs("results/emcal_barrel_electrons_Esim.png");
 
c3->SaveAs("results/emcal_barrel_electrons_Esim.pdf");
 
}
 
 
std::cout << "derp4\n";
 
{
 
TCanvas* c4 = new TCanvas("c4", "c4", 700, 500);
 
c4->SetLogy(1);
 
auto h = hEratio->DrawCopy();
 
//h->GetYaxis()->SetTitleOffset(1.4);
 
h->SetLineWidth(2);
 
h->SetLineColor(kBlue);
 
h = hEratio2->DrawCopy("same");
 
h->SetLineWidth(2);
 
h->SetLineColor(2);
 
//h->Fit("gaus", "", "", 0.01, 0.1);
 
//h->GetFunction("gaus")->SetLineWidth(2);
 
//h->GetFunction("gaus")->SetLineColor(kRed);
 
c4->SaveAs("results/emcal_barrel_Eratio.png");
 
c4->SaveAs("results/emcal_barrel_Eratio.pdf");
 
}
 
}
Loading