Skip to content
Snippets Groups Projects
Commit d9a49def authored by Wouter Deconinck's avatar Wouter Deconinck
Browse files

ZDC analysis of reconstructed quantities

parent 3495a7dd
No related branches found
No related tags found
1 merge request!255ZDC analysis of reconstructed quantities
Hide whitespace changes
Inline Side-by-side
Showing
with 362 additions and 551 deletions
.gitlab-ci.yml
+ 1
− 1
View file @ d9a49def
......@@ -88,7 +88,7 @@ include:
final_report:
stage: finish
needs: ["ecal_collect", "tracking_central_electrons", "clustering:results", "track_finding:collect", "track_fitting:collect"]
needs: ["ecal_collect", "tracking_central_electrons", "clustering:results", "track_finding:collect", "track_fitting:collect", "far_forward:collect"]
script:
# disabled while we address ACTS issues
#- mkdir -p results/views && cd results/views && bash ../../bin/download_views
......
benchmarks/far_forward/analysis/analysis_zdc_neutrons.cxx deleted 100644 → 0
+ 0
− 156
View file @ 3495a7dd
////////////////////////////////////////
// Read ROOT output file
// Plot variables
// Jihee Kim 09/2021
////////////////////////////////////////
#include "ROOT/RDataFrame.hxx"
#include <iostream>
#include <vector>
#include "edm4hep/MCParticleCollection.h"
#include "edm4hep/SimCalorimeterHitCollection.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TGaxis.h"
#include "TPaveStats.h"
#include "TMath.h"
#include "TH1.h"
#include "TF1.h"
#include "TH1D.h"
#include "TVector3.h"
R__LOAD_LIBRARY(libeicd.so)
#include "fmt/core.h"
#include "DD4hep/Detector.h"
#include "DDG4/Geant4Data.h"
#include "DDRec/CellIDPositionConverter.h"
using ROOT::RDataFrame;
void analysis_zdc_neutrons(const char* input_fname = "sim_zdc_uniform_neutron.edm4hep.root")
{
// Setting for graphs
TStyle* kStyle = new TStyle("kStyle","Kim's Style");
kStyle->SetOptStat("emr");
kStyle->SetOptTitle(0);
kStyle->SetOptFit(1101);
kStyle->SetStatColor(0);
kStyle->SetStatW(0.15);
kStyle->SetStatH(0.10);
kStyle->SetStatX(0.9);
kStyle->SetStatY(0.9);
kStyle->SetStatBorderSize(1);
kStyle->SetLabelSize(0.045,"xyz");
kStyle->SetTitleSize(0.045,"xyz");
kStyle->SetTitleOffset(1.2,"y");
kStyle->SetLineWidth(2);
kStyle->SetTitleFont(42,"xyz");
kStyle->SetLabelFont(42,"xyz");
kStyle->SetCanvasDefW(700);
kStyle->SetCanvasDefH(500);
kStyle->SetCanvasColor(0);
kStyle->SetPadTickX(1);
kStyle->SetPadTickY(1);
kStyle->SetPadGridX(1);
kStyle->SetPadGridY(1);
kStyle->SetPadLeftMargin(0.1);
kStyle->SetPadRightMargin(0.1);
kStyle->SetPadTopMargin(0.1);
kStyle->SetPadBottomMargin(0.1);
TGaxis::SetMaxDigits(3);
gROOT->SetStyle("kStyle");
ROOT::EnableImplicitMT();
ROOT::RDataFrame d0("events", input_fname);
// Thrown Energy [GeV]
auto Ethr = [](const std::vector<edm4hep::MCParticleData> & input) {
auto p = input[2];
auto energy = TMath::Sqrt(p.momentum.x * p.momentum.x + p.momentum.y * p.momentum.y + p.momentum.z * p.momentum.z + p.mass * p.mass);
return energy;
};
// Theta [mrad]
auto Thetathr = [](const std::vector<edm4hep::MCParticleData> & input) {
auto p = input[2];
auto theta = std::atan2(std::hypot(p.momentum.x, p.momentum.y), p.momentum.z);
return theta*1000.0;
};
// Phi [rad]
auto Phithr = [](const std::vector<edm4hep::MCParticleData> & input) {
auto p = input[2];
auto phi = std::atan2(p.momentum.y, p.momentum.x);
return phi;
};
// Eta
auto Etathr = [](const std::vector<edm4hep::MCParticleData> & input) {
auto p = input[2];
auto theta = std::atan2(std::hypot(p.momentum.x, p.momentum.y), p.momentum.z);
auto eta = -std::log(std::tan(theta/2));
return eta;
};
// Define variables
auto d1 = d0.Define("Ethr", Ethr, {"MCParticles"})
.Define("Thetathr", Thetathr, {"MCParticles"})
.Define("Phithr", Phithr, {"MCParticles"})
.Define("Etathr", Etathr, {"MCParticles"})
;
// Define Histograms
auto hEthr = d1.Histo1D({"hEthr", "Thrown Energy; E_{thr} [GeV]; Events", 100, 0.0, 0.3}, "Ethr");
auto hThetathr = d1.Histo1D({"hThetathr", "Thrown Theta; #theta_{thr} [mrad]; Events", 100, 0.0, 5.0}, "Thetathr");
auto hPhithr = d1.Histo1D({"hPhithr", "Thrown Phi; #phi_{thr} [rad]; Events", 100, 0.0, 5.0}, "Phithr");
auto hEtathr = d1.Histo1D({"hEtathr", "Thrown Eta; #eta_{thr}; Events", 100, 6.0, 15.0}, "Etathr");
auto hPhiThetathr = d1.Histo2D({"hPhiThetathr", "Thrown Phi vs Theta; #theta_{thr} [mrad]; #phi_{thr}", 20, 0.0, 5.0, 20, 0.0, 5.0},
"Thetathr", "Phithr");
// Event Counts
auto nevents_thrown = d1.Count();
std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
// Draw Histograms
{
TCanvas* c1 = new TCanvas("c1", "c1");
auto h = hEthr->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c1->SaveAs("results/far_forward/zdc/zdc_neutron_Ethr.png");
c1->SaveAs("results/far_forward/zdc/zdc_neutron_Ethr.pdf");
}
{
TCanvas* c2 = new TCanvas("c2", "c2");
auto h = hThetathr->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c2->SaveAs("results/far_forward/zdc/zdc_neutron_Thetathr.png");
c2->SaveAs("results/far_forward/zdc/zdc_neutron_Thetathr.pdf");
}
{
TCanvas* c3 = new TCanvas("c3", "c3");
auto h = hPhithr->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c3->SaveAs("results/far_forward/zdc/zdc_neutron_Phithr.png");
c3->SaveAs("results/far_forward/zdc/zdc_neutron_Phithr.pdf");
}
{
TCanvas* c4 = new TCanvas("c4", "c4");
auto h = hEtathr->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c4->SaveAs("results/far_forward/zdc/zdc_neutron_Etathr.png");
c4->SaveAs("results/far_forward/zdc/zdc_neutron_Etathr.pdf");
}
{
TCanvas* c5 = new TCanvas("c5", "c5");
auto h = hPhiThetathr->DrawCopy("COLZ");
c5->SaveAs("results/far_forward/zdc/zdc_neutron_PhiThetathr.png");
c5->SaveAs("results/far_forward/zdc/zdc_neutron_PhiThetathr.pdf");
}
}
benchmarks/far_forward/analysis/analysis_zdc_particles.cxx 0 → 100644
+ 306
− 0
View file @ d9a49def
////////////////////////////////////////
// Read ROOT output file
// Plot variables
// Jihee Kim 09/2021
////////////////////////////////////////
#include "ROOT/RDataFrame.hxx"
#include <iostream>
#include <vector>
#include "edm4hep/MCParticleCollection.h"
#include "edm4hep/SimCalorimeterHitCollection.h"
#include "eicd/ClusterCollection.h"
#include "eicd/ClusterData.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TGaxis.h"
#include "TPaveStats.h"
#include "TMath.h"
#include "TH1.h"
#include "TF1.h"
#include "TH1D.h"
#include "TVector3.h"
R__LOAD_LIBRARY(libeicd.so)
#include "fmt/core.h"
#include "DD4hep/Detector.h"
#include "DDG4/Geant4Data.h"
#include "DDRec/CellIDPositionConverter.h"
using ROOT::RDataFrame;
void analysis_zdc_particles(
const std::string& input_fname = "sim_zdc_uniform.edm4hep.root",
const std::string& results_path = "results/far_forward/zdc/"
) {
// Setting for graphs
TStyle* kStyle = new TStyle("kStyle","Kim's Style");
kStyle->SetOptStat("emr");
kStyle->SetOptTitle(0);
kStyle->SetOptFit(1101);
kStyle->SetStatColor(0);
kStyle->SetStatW(0.15);
kStyle->SetStatH(0.10);
kStyle->SetStatX(0.9);
kStyle->SetStatY(0.9);
kStyle->SetStatBorderSize(1);
kStyle->SetLabelSize(0.045,"xyz");
kStyle->SetTitleSize(0.045,"xyz");
kStyle->SetTitleOffset(1.2,"y");
kStyle->SetLineWidth(2);
kStyle->SetTitleFont(42,"xyz");
kStyle->SetLabelFont(42,"xyz");
kStyle->SetCanvasDefW(700);
kStyle->SetCanvasDefH(500);
kStyle->SetCanvasColor(0);
kStyle->SetPadTickX(1);
kStyle->SetPadTickY(1);
kStyle->SetPadGridX(1);
kStyle->SetPadGridY(1);
kStyle->SetPadLeftMargin(0.1);
kStyle->SetPadRightMargin(0.1);
kStyle->SetPadTopMargin(0.1);
kStyle->SetPadBottomMargin(0.1);
TGaxis::SetMaxDigits(3);
gROOT->SetStyle("kStyle");
ROOT::EnableImplicitMT();
ROOT::RDataFrame d0("events", input_fname);
// Thrown Energy [GeV]
auto E_gen = [](const std::vector<edm4hep::MCParticleData> & input) {
auto p = input[2];
auto energy = TMath::Sqrt(p.momentum.x * p.momentum.x + p.momentum.y * p.momentum.y + p.momentum.z * p.momentum.z + p.mass * p.mass);
return energy;
};
// Theta [rad]
auto Theta_gen = [](const std::vector<edm4hep::MCParticleData> & input) {
auto p = input[2];
auto theta = std::atan2(std::hypot(p.momentum.x, p.momentum.y), p.momentum.z);
return theta;
};
auto Thetarec = [](const std::vector<eicd::ClusterData> & input) {
std::vector<float> output;
for (const auto &c: input) {
auto r = c.position;
auto theta = std::atan2(std::hypot(r.x, r.y), r.z);
output.push_back(theta);
}
return output;
};
// Phi [rad]
auto Phi_gen = [](const std::vector<edm4hep::MCParticleData> & input) {
auto p = input[2];
auto phi = std::atan2(p.momentum.y, p.momentum.x);
if (phi < 0) phi += 2.*M_PI;
return phi;
};
auto Phirec = [](const std::vector<eicd::ClusterData> & input) {
std::vector<float> output;
for (const auto &c: input) {
auto r = c.position;
auto phi = std::atan2(r.y, r.x);
if (phi < 0) phi += 2.*M_PI;
output.push_back(phi);
}
return output;
};
// Eta
auto Eta_gen = [](const std::vector<edm4hep::MCParticleData> & input) {
auto p = input[2];
auto theta = std::atan2(std::hypot(p.momentum.x, p.momentum.y), p.momentum.z);
auto eta = -std::log(std::tan(theta/2));
return eta;
};
auto Etarec = [](const std::vector<eicd::ClusterData> & input) {
std::vector<float> output;
for (const auto &c: input) {
auto r = c.position;
auto theta = std::atan2(std::hypot(r.x, r.y), r.z);
auto eta = -std::log(std::tan(theta/2));
output.push_back(eta);
}
return output;
};
// Define variables
auto d1 = d0
// Truth
.Define("E_gen", E_gen, {"MCParticles"})
.Define("Theta_gen", Theta_gen, {"MCParticles"})
.Define("Phi_gen", Phi_gen, {"MCParticles"})
.Define("Eta_gen", Eta_gen, {"MCParticles"})
// Ecal
.Define("E_Ecal", {"ZDCEcalClusters.energy"})
.Define("Theta_Ecal", Thetarec, {"ZDCEcalClusters"})
.Define("Phi_Ecal", Phirec, {"ZDCEcalClusters"})
.Define("Eta_Ecal", Etarec, {"ZDCEcalClusters"})
// HCal
.Define("E_Hcal", {"ZDCHcalClusters.energy"})
.Define("Theta_Hcal", Thetarec, {"ZDCHcalClusters"})
.Define("Phi_Hcal", Phirec, {"ZDCHcalClusters"})
.Define("Eta_Hcal", Etarec, {"ZDCHcalClusters"})
;
// Histogram boundaries
const auto E_nbins = 100;
const auto E_max = 50.0;
const auto theta_nbins = 50;
const auto theta_max = 0.050;
const auto phi_nbins = 90;
const auto phi_min = 0.0;
const auto phi_max = 2.*M_PI;
const auto eta_nbins = 40;
const auto eta_min = 6.0;
const auto eta_max = 10.0;
// Define Histograms
auto hE_gen = d1.Histo1D({"hE_gen", "Thrown Energy; E_{gen} [GeV]; Events", E_nbins, 0.0, E_max}, "E_gen");
auto hTheta_gen = d1.Histo1D({"hTheta_gen", "Thrown Theta; #theta_{gen} [rad]; Events", theta_nbins, 0.0, theta_max}, "Theta_gen");
auto hPhi_gen = d1.Histo1D({"hPhi_gen", "Thrown Phi; #phi_{gen} [rad]; Events", phi_nbins, phi_min, phi_max}, "Phi_gen");
auto hEta_gen = d1.Histo1D({"hEta_gen", "Thrown Eta; #eta_{gen}; Events", eta_nbins, eta_min, eta_max}, "Eta_gen");
auto hPhiTheta_gen = d1.Histo2D({"hPhiTheta_gen", "Thrown Phi vs Theta; #theta_{gen} [rad]; #phi_{gen}", theta_nbins, 0.0, theta_max, phi_nbins, phi_min, phi_max},
"Theta_gen", "Phi_gen");
auto hE_Ecal = d1.Histo1D({"hE_Ecal", "Reconstructed Ecal Energy; E_{rec} [GeV]; Events", E_nbins, 0.0, E_max}, "E_Ecal");
auto hTheta_Ecal = d1.Histo1D({"hTheta_Ecal", "Reconstructed Ecal Theta; #theta_{rec} [rad]; Events", theta_nbins, 0.0, theta_max}, "Theta_Ecal");
auto hPhi_Ecal = d1.Histo1D({"hPhi_Ecal", "Reconstructed Ecal Phi; #phi_{rec} [rad]; Events", phi_nbins, phi_min, phi_max}, "Phi_Ecal");
auto hEta_Ecal = d1.Histo1D({"hEta_Ecal", "Reconstructed Ecal Eta; #eta_{rec}; Events", eta_nbins, eta_min, eta_max}, "Eta_Ecal");
auto hPhiTheta_Ecal = d1.Histo2D({"hPhiTheta_Ecal", "Reconstructed Ecal Phi vs Theta; #theta_{rec} [rad]; #phi_{rec}", theta_nbins, 0.0, theta_max, phi_nbins, phi_min, phi_max},
"Theta_Ecal", "Phi_Ecal");
auto hE_Hcal = d1.Histo1D({"hE_Hcal", "Reconstructed Hcal Energy; E_{rec} [GeV]; Events", E_nbins, 0.0, E_max}, "E_Hcal");
auto hTheta_Hcal = d1.Histo1D({"hTheta_Hcal", "Reconstructed Hcal Theta; #theta_{rec} [rad]; Events", theta_nbins, 0.0, theta_max}, "Theta_Hcal");
auto hPhi_Hcal = d1.Histo1D({"hPhi_Hcal", "Reconstructed Hcal Phi; #phi_{rec} [rad]; Events", phi_nbins, phi_min, phi_max}, "Phi_Hcal");
auto hEta_Hcal = d1.Histo1D({"hEta_Hcal", "Reconstructed Hcal Eta; #eta_{rec}; Events", eta_nbins, eta_min, eta_max}, "Eta_Hcal");
auto hPhiTheta_Hcal = d1.Histo2D({"hPhiTheta_Hcal", "Reconstructed Hcal Phi vs Theta; #theta_{rec} [rad]; #phi_{rec}", theta_nbins, 0.0, theta_max, phi_nbins, phi_min, phi_max},
"Theta_Hcal", "Phi_Hcal");
// Event Counts
auto nevents_thrown = d1.Count();
std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
// Draw Histograms
{
TCanvas* c1 = new TCanvas("c1", "c1");
auto h = hE_gen->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c1->SaveAs(TString(results_path + "/zdc_E_gen.png"));
c1->SaveAs(TString(results_path + "/zdc_E_gen.pdf"));
}
{
TCanvas* c2 = new TCanvas("c2", "c2");
auto h = hTheta_gen->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c2->SaveAs(TString(results_path + "/zdc_Theta_gen.png"));
c2->SaveAs(TString(results_path + "/zdc_Theta_gen.pdf"));
}
{
TCanvas* c3 = new TCanvas("c3", "c3");
auto h = hPhi_gen->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c3->SaveAs(TString(results_path + "/zdc_Phi_gen.png"));
c3->SaveAs(TString(results_path + "/zdc_Phi_gen.pdf"));
}
{
TCanvas* c4 = new TCanvas("c4", "c4");
auto h = hEta_gen->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c4->SaveAs(TString(results_path + "/zdc_Eta_gen.png"));
c4->SaveAs(TString(results_path + "/zdc_Eta_gen.pdf"));
}
{
TCanvas* c5 = new TCanvas("c5", "c5");
auto h = hPhiTheta_gen->DrawCopy("COLZ");
c5->SaveAs(TString(results_path + "/zdc_PhiTheta_gen.png"));
c5->SaveAs(TString(results_path + "/zdc_PhiTheta_gen.pdf"));
}
// Draw Histograms Ecal
{
TCanvas* c1 = new TCanvas("c1", "c1");
auto h = hE_Ecal->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c1->SaveAs(TString(results_path + "/zdc_E_Ecal.png"));
c1->SaveAs(TString(results_path + "/zdc_E_Ecal.pdf"));
}
{
TCanvas* c2 = new TCanvas("c2", "c2");
auto h = hTheta_Ecal->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c2->SaveAs(TString(results_path + "/zdc_Theta_Ecal.png"));
c2->SaveAs(TString(results_path + "/zdc_Theta_Ecal.pdf"));
}
{
TCanvas* c3 = new TCanvas("c3", "c3");
auto h = hPhi_Ecal->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c3->SaveAs(TString(results_path + "/zdc_Phi_Ecal.png"));
c3->SaveAs(TString(results_path + "/zdc_Phi_Ecal.pdf"));
}
{
TCanvas* c4 = new TCanvas("c4", "c4");
auto h = hEta_Ecal->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c4->SaveAs(TString(results_path + "/zdc_Eta_Ecal.png"));
c4->SaveAs(TString(results_path + "/zdc_Eta_Ecal.pdf"));
}
{
TCanvas* c5 = new TCanvas("c5", "c5");
auto h = hPhiTheta_Ecal->DrawCopy("COLZ");
c5->SaveAs(TString(results_path + "/zdc_PhiTheta_Ecal.png"));
c5->SaveAs(TString(results_path + "/zdc_PhiTheta_Ecal.pdf"));
}
// Draw Histograms Hcal
{
TCanvas* c1 = new TCanvas("c1", "c1");
auto h = hE_Hcal->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c1->SaveAs(TString(results_path + "/zdc_E_Hcal.png"));
c1->SaveAs(TString(results_path + "/zdc_E_Hcal.pdf"));
}
{
TCanvas* c2 = new TCanvas("c2", "c2");
auto h = hTheta_Hcal->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c2->SaveAs(TString(results_path + "/zdc_Theta_Hcal.png"));
c2->SaveAs(TString(results_path + "/zdc_Theta_Hcal.pdf"));
}
{
TCanvas* c3 = new TCanvas("c3", "c3");
auto h = hPhi_Hcal->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c3->SaveAs(TString(results_path + "/zdc_Phi_Hcal.png"));
c3->SaveAs(TString(results_path + "/zdc_Phi_Hcal.pdf"));
}
{
TCanvas* c4 = new TCanvas("c4", "c4");
auto h = hEta_Hcal->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c4->SaveAs(TString(results_path + "/zdc_Eta_Hcal.png"));
c4->SaveAs(TString(results_path + "/zdc_Eta_Hcal.pdf"));
}
{
TCanvas* c5 = new TCanvas("c5", "c5");
auto h = hPhiTheta_Hcal->DrawCopy("COLZ");
c5->SaveAs(TString(results_path + "/zdc_PhiTheta_Hcal.png"));
c5->SaveAs(TString(results_path + "/zdc_PhiTheta_Hcal.pdf"));
}
}
benchmarks/far_forward/analysis/analysis_zdc_photons.cxx deleted 100644 → 0
+ 0
− 156
View file @ 3495a7dd
////////////////////////////////////////
// Read ROOT output file
// Plot variables
// Jihee Kim 09/2021
////////////////////////////////////////
#include "ROOT/RDataFrame.hxx"
#include <iostream>
#include <vector>
#include "edm4hep/MCParticleCollection.h"
#include "edm4hep/CalorimeterHitCollection.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TGaxis.h"
#include "TPaveStats.h"
#include "TMath.h"
#include "TH1.h"
#include "TF1.h"
#include "TH1D.h"
#include "TVector3.h"
R__LOAD_LIBRARY(libeicd.so)
#include "fmt/core.h"
#include "DD4hep/Detector.h"
#include "DDG4/Geant4Data.h"
#include "DDRec/CellIDPositionConverter.h"
using ROOT::RDataFrame;
void analysis_zdc_photons(const char* input_fname = "sim_zdc_uniform_photon.edm4hep.root")
{
// Setting for graphs
TStyle* kStyle = new TStyle("kStyle","Kim's Style");
kStyle->SetOptStat("emr");
kStyle->SetOptTitle(0);
kStyle->SetOptFit(1101);
kStyle->SetStatColor(0);
kStyle->SetStatW(0.15);
kStyle->SetStatH(0.10);
kStyle->SetStatX(0.9);
kStyle->SetStatY(0.9);
kStyle->SetStatBorderSize(1);
kStyle->SetLabelSize(0.045,"xyz");
kStyle->SetTitleSize(0.045,"xyz");
kStyle->SetTitleOffset(1.2,"y");
kStyle->SetLineWidth(2);
kStyle->SetTitleFont(42,"xyz");
kStyle->SetLabelFont(42,"xyz");
kStyle->SetCanvasDefW(700);
kStyle->SetCanvasDefH(500);
kStyle->SetCanvasColor(0);
kStyle->SetPadTickX(1);
kStyle->SetPadTickY(1);
kStyle->SetPadGridX(1);
kStyle->SetPadGridY(1);
kStyle->SetPadLeftMargin(0.1);
kStyle->SetPadRightMargin(0.1);
kStyle->SetPadTopMargin(0.1);
kStyle->SetPadBottomMargin(0.1);
TGaxis::SetMaxDigits(3);
gROOT->SetStyle("kStyle");
ROOT::EnableImplicitMT();
ROOT::RDataFrame d0("events", input_fname);
// Thrown Energy [GeV]
auto Ethr = [](const std::vector<edm4hep::MCParticleData> & input) {
auto p = input[2];
auto energy = TMath::Sqrt(p.momentum.x * p.momentum.x + p.momentum.y * p.momentum.y + p.momentum.z * p.momentum.z + p.mass * p.mass);
return energy;
};
// Theta [mrad]
auto Thetathr = [](const std::vector<edm4hep::MCParticleData> & input) {
auto p = input[2];
auto theta = std::atan2(std::hypot(p.momentum.x, p.momentum.y), p.momentum.z);
return theta*1000.0;
};
// Phi [rad]
auto Phithr = [](const std::vector<edm4hep::MCParticleData> & input) {
auto p = input[2];
auto phi = std::atan2(p.momentum.y, p.momentum.x);
return phi;
};
// Eta
auto Etathr = [](const std::vector<edm4hep::MCParticleData> & input) {
auto p = input[2];
auto theta = std::atan2(std::hypot(p.momentum.x, p.momentum.y), p.momentum.z);
auto eta = -std::log(std::tan(theta/2));
return eta;
};
// Define variables
auto d1 = d0.Define("Ethr", Ethr, {"MCParticles"})
.Define("Thetathr", Thetathr, {"MCParticles"})
.Define("Phithr", Phithr, {"MCParticles"})
.Define("Etathr", Etathr, {"MCParticles"})
;
// Define Histograms
auto hEthr = d1.Histo1D({"hEthr", "Thrown Energy; E_{thr} [GeV]; Events", 100, 0.0, 0.3}, "Ethr");
auto hThetathr = d1.Histo1D({"hThetathr", "Thrown Theta; #theta_{thr} [mrad]; Events", 100, 0.0, 5.0}, "Thetathr");
auto hPhithr = d1.Histo1D({"hPhithr", "Thrown Phi; #phi_{thr} [rad]; Events", 100, 0.0, 5.0}, "Phithr");
auto hEtathr = d1.Histo1D({"hEtathr", "Thrown Eta; #eta_{thr}; Events", 100, 6.0, 15.0}, "Etathr");
auto hPhiThetathr = d1.Histo2D({"hPhiThetathr", "Thrown Phi vs Theta; #theta_{thr} [mrad]; #phi_{thr}", 20, 0.0, 5.0, 20, 0.0, 5.0},
"Thetathr", "Phithr");
// Event Counts
auto nevents_thrown = d1.Count();
std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
// Draw Histograms
{
TCanvas* c1 = new TCanvas("c1", "c1");
auto h = hEthr->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c1->SaveAs("results/far_forward/zdc/zdc_photon_Ethr.png");
c1->SaveAs("results/far_forward/zdc/zdc_photon_Ethr.pdf");
}
{
TCanvas* c2 = new TCanvas("c2", "c2");
auto h = hThetathr->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c2->SaveAs("results/far_forward/zdc/zdc_photon_Thetathr.png");
c2->SaveAs("results/far_forward/zdc/zdc_photon_Thetathr.pdf");
}
{
TCanvas* c3 = new TCanvas("c3", "c3");
auto h = hPhithr->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c3->SaveAs("results/far_forward/zdc/zdc_photon_Phithr.png");
c3->SaveAs("results/far_forward/zdc/zdc_photon_Phithr.pdf");
}
{
TCanvas* c4 = new TCanvas("c4", "c4");
auto h = hEtathr->DrawCopy();
h->SetLineWidth(3);
h->SetLineColor(kBlue);
c4->SaveAs("results/far_forward/zdc/zdc_photon_Etathr.png");
c4->SaveAs("results/far_forward/zdc/zdc_photon_Etathr.pdf");
}
{
TCanvas* c5 = new TCanvas("c5", "c5");
auto h = hPhiThetathr->DrawCopy("COLZ");
c5->SaveAs("results/far_forward/zdc/zdc_photon_PhiThetathr.png");
c5->SaveAs("results/far_forward/zdc/zdc_photon_PhiThetathr.pdf");
}
}
benchmarks/far_forward/analysis/gen_zdc_neutrons.cxx → benchmarks/far_forward/analysis/gen_zdc_particles.cxx
+ 29
− 12
View file @ d9a49def
//////////////////////////////////////////////////////////////
// ZDC detector
// Single Neutron dataset
// Single Particle dataset
// J. Kim 09/2021
//////////////////////////////////////////////////////////////
#include "HepMC3/GenEvent.h"
......@@ -14,10 +14,12 @@
#include <math.h>
#include <random>
#include <TRandom.h>
#include <Math/Vector3D.h>
#include <Math/RotationY.h>
using namespace HepMC3;
void gen_zdc_neutrons(int n_events = 1e6, double e_start = 20.0, double e_end = 140.0, const char* out_fname = "zdc_neutron.hepmc") {
void gen_zdc_particles(int n_events = 1e6, std::string particle = "neutron", double p_start = 20.0, double p_end = 140.0, const std::string& out_fname = "zdc_neutron.hepmc") {
WriterAscii hepmc_output(out_fname);
int events_parsed = 0;
GenEvent evt(Units::GEV, Units::MM);
......@@ -25,9 +27,12 @@ void gen_zdc_neutrons(int n_events = 1e6, double e_start = 20.0, double e_end =
// Random number generator
TRandom* r1 = new TRandom();
// Set crossing angle [rad]
double crossing_angle = -0.025;
// Set polar angle range [rad]
double theta_min = 0.0;
double theta_max = 0.004;
double theta_max = 0.015;
for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
// FourVector(px,py,pz,e,pdgid,status)
......@@ -38,22 +43,34 @@ void gen_zdc_neutrons(int n_events = 1e6, double e_start = 20.0, double e_end =
GenParticlePtr p2 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
// Define momentum
Double_t p = r1->Uniform(e_start, e_end);
Double_t p_mag = r1->Uniform(p_start, p_end);
Double_t phi = r1->Uniform(0.0, 2.0 * M_PI);
Double_t theta = acos(r1->Uniform(cos(theta_max), cos(theta_min)));
Double_t px = p * std::cos(phi) * std::sin(theta);
Double_t py = p * std::sin(phi) * std::sin(theta);
Double_t pz = p * std::cos(theta);
Double_t p0_x = p_mag * std::cos(phi) * std::sin(theta);
Double_t p0_y = p_mag * std::sin(phi) * std::sin(theta);
Double_t p0_z = p_mag * std::cos(theta);
// Rotate the vector in Y by crossing angle -25 mrad when particles are being generated
//ROOT::Math::XYZVector p0 = {px,py,pz};
//ROOT::Math::RotationY r(-0.025);
//auto p_rot = r*p0;
// Rotate the vector in Y by crossing angle when particles are being generated
ROOT::Math::XYZVector p0{p0_x, p0_y, p0_z};
ROOT::Math::RotationY r(crossing_angle);
auto p = r * p0;
auto px = p.X();
auto py = p.Y();
auto pz = p.Z();
// FourVector(px,py,pz,e,pdgid,status)
// status - type 1 is final state
// pdgid 22 - photon massless
// pdgid 2112 - neutron 939.565 MeV/c^2
GenParticlePtr p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p + (0.939565 * 0.939565))), 2112, 1);
GenParticlePtr p3;
if (particle == "neutron") {
p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, std::hypot(p_mag, 0.939565)), 2112, 1);
} else if (particle == "photon") {
p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, p_mag), 22, 1);
} else {
std::cout << "Particle type " << particle << " not recognized!" << std::endl;
exit(-1);
}
// Create vertex
GenVertexPtr v1 = std::make_shared<GenVertex>();
......
benchmarks/far_forward/analysis/gen_zdc_photons.cxx deleted 100644 → 0
+ 0
− 78
View file @ 3495a7dd
//////////////////////////////////////////////////////////////
// ZDC detector
// Single Neutron dataset
// J. Kim 09/2021
//////////////////////////////////////////////////////////////
#include "HepMC3/GenEvent.h"
#include "HepMC3/Print.h"
#include "HepMC3/ReaderAscii.h"
#include "HepMC3/WriterAscii.h"
#include <TMath.h>
#include <cmath>
#include <iostream>
#include <math.h>
#include <random>
#include <TRandom.h>
using namespace HepMC3;
void gen_zdc_photons(int n_events = 1e6, double e_start = 20.0, double e_end = 140.0, const char* out_fname = "zdc_photon.hepmc") {
WriterAscii hepmc_output(out_fname);
int events_parsed = 0;
GenEvent evt(Units::GEV, Units::MM);
// Random number generator
TRandom* r1 = new TRandom();
// Set polar angle range [rad]
double theta_min = 0.0;
double theta_max = 0.004;
for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
// FourVector(px,py,pz,e,pdgid,status)
// status - type 4 is beam
// pdgid 11 - electron
// pdgid 2212 - proton
GenParticlePtr p1 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
GenParticlePtr p2 = std::make_shared<GenParticle>(FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4);
// Define momentum
Double_t p = r1->Uniform(e_start, e_end);
Double_t phi = r1->Uniform(0.0, 2.0 * M_PI);
Double_t theta = acos(r1->Uniform(cos(theta_max), cos(theta_min)));
Double_t px = p * std::cos(phi) * std::sin(theta);
Double_t py = p * std::sin(phi) * std::sin(theta);
Double_t pz = p * std::cos(theta);
// Rotate the vector in Y by crossing angle -25 mrad when particles are being generated
//ROOT::Math::XYZVector p0 = {px,py,pz};
//ROOT::Math::RotationY r(-0.025);
//auto p_rot = r*p0;
// FourVector(px,py,pz,e,pdgid,status)
// status - type 1 is final state
// pdgid 22 - photon
GenParticlePtr p3 = std::make_shared<GenParticle>(FourVector(px, py, pz, sqrt(p * p)), 22, 1);
// Create vertex
GenVertexPtr v1 = std::make_shared<GenVertex>();
v1->add_particle_in(p1);
v1->add_particle_in(p2);
v1->add_particle_out(p3);
evt.add_vertex(v1);
if (events_parsed == 0) {
std::cout << "First event: " << std::endl;
Print::listing(evt);
}
hepmc_output.write_event(evt);
if (events_parsed % 100 == 0) {
std::cout << "Event: " << events_parsed << std::endl;
}
evt.clear();
}
hepmc_output.close();
std::cout << "Events parsed and written: " << events_parsed << std::endl;
}
benchmarks/far_forward/config.yml
+ 11
− 7
View file @ d9a49def
......@@ -11,16 +11,20 @@ B0_far_forward_protons:
script:
- bash benchmarks/far_forward/far_forward_protons.sh
ZDC_far_forward_neutrons:
ZDC_far_forward:
extends: .rec_benchmark
stage: run
needs: ["far_forward:compile"]
script:
- bash benchmarks/far_forward/run_zdc_neutrons.sh
- bash benchmarks/far_forward/run_zdc.sh --particle $PARTICLE
parallel:
matrix:
- PARTICLE: ["neutron", "photon"]
ZDC_far_forward_photons:
extends: .rec_benchmark
stage: run
needs: ["far_forward:compile"]
far_forward:collect:
stage: collect
needs:
- ["B0_far_forward_protons", "ZDC_far_forward"]
script:
- bash benchmarks/far_forward/run_zdc_photons.sh
- echo "Done collecting artifacts."
benchmarks/far_forward/run_zdc_neutrons.sh → benchmarks/far_forward/run_zdc.sh
+ 15
− 7
View file @ d9a49def
#!/bin/bash
function print_the_help {
echo "USAGE: ${0} [--rec-only] "
echo "USAGE: ${0} [--rec-only] [--ana-only] [--particle=neutron] "
echo " OPTIONS: "
echo " --rec-only only run gaudi reconstruction part"
echo " --ana-only only run analysis macros part"
exit
}
PARTICLE="neutron"
REC_ONLY=
ANALYSIS_ONLY=
......@@ -28,6 +30,11 @@ do
ANALYSIS_ONLY=1
shift # past value
;;
--particle)
PARTICLE="$2"
shift # past argument
shift # past value
;;
*) # unknown option
#POSITIONAL+=("$1") # save it in an array for later
echo "unknown option $1"
......@@ -69,19 +76,22 @@ fi
export DETECTOR_COMPACT_PATH=${DETECTOR_PATH}/${JUGGLER_DETECTOR_CONFIG}.xml
export JUGGLER_FILE_NAME_TAG="zdc_neutrons"
export JUGGLER_FILE_NAME_TAG="zdc_${PARTICLE}"
export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.edm4hep.root"
export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
export RESULTS_PATH="results/far_forward/zdc/${PARTICLE}"
mkdir -p ${RESULTS_PATH}
echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
if [[ -z "${REC_ONLY}" && -z "${ANALYSIS_ONLY}" ]] ;
then
echo "Generating input events"
root -b -q "benchmarks/far_forward/analysis/gen_zdc_neutrons.cxx(${JUGGLER_N_EVENTS}, ${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
root -b -q "benchmarks/far_forward/analysis/gen_zdc_particles.cxx(${JUGGLER_N_EVENTS}, \"${PARTICLE}\", ${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running script: generating input events"
exit 1
......@@ -113,11 +123,10 @@ then
fi
fi
mkdir -p results/far_forward
mkdir -p results/far_forward/zdc
rootls -t ${JUGGLER_REC_FILE}
echo "Running analysis root scripts"
root -b -q "benchmarks/far_forward/analysis/analysis_zdc_neutrons.cxx+(\"${JUGGLER_REC_FILE}\")"
root -b -q "benchmarks/far_forward/analysis/analysis_zdc_particles.cxx+(\"${JUGGLER_REC_FILE}\", \"${RESULTS_PATH}\")"
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running root analysis script"
exit 1
......@@ -131,4 +140,3 @@ if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then
cp ${JUGGLER_REC_FILE} results/far_forward/zdc/.
fi
fi
benchmarks/far_forward/run_zdc_photons.sh deleted 100644 → 0
+ 0
− 134
View file @ 3495a7dd
#!/bin/bash
function print_the_help {
echo "USAGE: ${0} [--rec-only] "
echo " OPTIONS: "
echo " --rec-only only run gaudi reconstruction part"
exit
}
REC_ONLY=
ANALYSIS_ONLY=
POSITIONAL=()
while [[ $# -gt 0 ]]
do
key="$1"
case $key in
-h|--help)
shift # past argument
print_the_help
;;
--rec-only)
REC_ONLY=1
shift # past value
;;
--ana-only)
ANALYSIS_ONLY=1
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
print_env.sh
## To run the reconstruction, we need the following global variables:
## - JUGGLER_INSTALL_PREFIX: Install prefix for Juggler (simu/recon)
## - JUGGLER_DETECTOR: the detector package we want to use for this benchmark
## - JUGGLER_DETECTOR_VERSION: the detector package we want to use for this benchmark
## - DETECTOR_PATH: full path to the detector definitions
##
## You can ready options/env.sh for more in-depth explanations of the variables
## and how they can be controlled.
export DETECTOR_PATH=${DETECTOR_PATH}
if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then
export JUGGLER_N_EVENTS=100
fi
if [[ ! -n "${E_start}" ]] ; then
export E_start=0.05
fi
if [[ ! -n "${E_end}" ]] ; then
export E_end=0.25
fi
if [[ ! -n "${ZDC_SCI_SAMP_FRAC}" ]] ; then
export ZDC_PbSCI_SAMP_FRAC=1.0
fi
export DETECTOR_COMPACT_PATH=${DETECTOR_PATH}/${JUGGLER_DETECTOR_CONFIG}.xml
export JUGGLER_FILE_NAME_TAG="zdc_photons"
export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc"
export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.edm4hep.root"
export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root"
echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}"
echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}"
if [[ -z "${REC_ONLY}" && -z "${ANALYSIS_ONLY}" ]] ;
then
echo "Generating input events"
root -b -q "benchmarks/far_forward/analysis/gen_zdc_photons.cxx(${JUGGLER_N_EVENTS}, ${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running script: generating input events"
exit 1
fi
echo "Running Geant4 simulation"
ddsim --runType batch \
--part.minimalKineticEnergy 0.5*MeV \
--filter.tracker edep0 \
-v WARNING \
--numberOfEvents ${JUGGLER_N_EVENTS} \
--compactFile ${DETECTOR_PATH}/${JUGGLER_DETECTOR_CONFIG}.xml \
--inputFiles ${JUGGLER_FILE_NAME_TAG}.hepmc \
--outputFile ${JUGGLER_SIM_FILE}
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running npdet: Geant4 simulation"
exit 1
fi
fi
rootls -t ${JUGGLER_SIM_FILE}
if [[ -z "${ANALYSIS_ONLY}" ]] ;
then
gaudirun.py benchmarks/far_forward/options/zdc_reconstruction.py
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running juggler"
exit 1
fi
fi
mkdir -p results/far_forward
mkdir -p results/far_forward/zdc
echo "Running analysis root scripts"
root -b -q "benchmarks/far_forward/analysis/analysis_zdc_photons.cxx+(\"${JUGGLER_REC_FILE}\")"
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running root analysis script"
exit 1
fi
root_filesize=$(stat --format=%s "${JUGGLER_REC_FILE}")
if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then
# file must be less than 10 MB to upload
if [[ "${root_filesize}" -lt "10000000" ]] ; then
cp ${JUGGLER_SIM_FILE} results/far_forward/zdc/.
cp ${JUGGLER_REC_FILE} results/far_forward/zdc/.
fi
fi
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment