diff --git a/benchmarks/barrel_ecal/scripts/e_pi_separation.cxx b/benchmarks/barrel_ecal/scripts/e_pi_separation.cxx new file mode 100644 index 0000000000000000000000000000000000000000..8b002c646725a264257f6ffe68652c942720a5d0 --- /dev/null +++ b/benchmarks/barrel_ecal/scripts/e_pi_separation.cxx @@ -0,0 +1,156 @@ +//////////////////////////////////////// +// 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" + +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(); + std::cout << decoder->fieldDescription() << "\n"; + + + + + // Thrown Energy [GeV] + 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; + }; + + // 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("nhits", nhits, {"EcalBarrelHits"}) + .Define("Esim", Esim, {"EcalBarrelHits"}) + .Define("fsam", fsam, {"Esim", "Ethr"}); + + // 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, 1.0}, + "Esim"); + auto hfsam = d1.Histo1D( + {"hfsam", "Sampling Fraction; Sampling Fraction; Events", 100, 0.0, 0.1}, + "fsam"); + + // 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); + 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 = hfsam->DrawCopy(); + //h->GetYaxis()->SetTitleOffset(1.4); + h->SetLineWidth(2); + h->SetLineColor(kBlue); + //h->Fit("gaus", "", "", 0.01, 0.1); + //h->GetFunction("gaus")->SetLineWidth(2); + //h->GetFunction("gaus")->SetLineColor(kRed); + c4->SaveAs("results/emcal_barrel_electrons_fsam.png"); + c4->SaveAs("results/emcal_barrel_electrons_fsam.pdf"); + } +}