Skip to content
Snippets Groups Projects
Commit 87c352f1 authored by Sylvester Joosten's avatar Sylvester Joosten
Browse files

Prepare canyonlands

parent 4573b716
No related branches found
No related tags found
1 merge request!96Prepare canyonlands
Hide whitespace changes
Inline Side-by-side
.clang-format
+ 5
− 0
View file @ 87c352f1
......@@ -7,4 +7,9 @@ Standard: Cpp11
#SpaceBeforeParens: ControlStatements
SpaceAfterControlStatementKeyword: true
PointerBindsToType: true
IncludeBlocks: Preserve
UseTab: Never
ColumnLimit: 120
NamespaceIndentation: Inner
AlignConsecutiveAssignments: true
...
.rootlogon.C
+ 2
− 0
View file @ 87c352f1
......@@ -2,6 +2,8 @@
// top-level include-dir
gROOT->ProcessLine(".include include");
R__LOAD_LIBRARY(fmt)
// Setting for Graphs
gROOT->SetStyle("Plain");
gStyle->SetOptFit(1);
......
README.md
+ 1
− 1
View file @ 87c352f1
......@@ -35,7 +35,7 @@ export BEAMLINE_CONFIG=ip6 # ip6 is the default
```
git@eicweb.phy.anl.gov:EIC/benchmarks/detector_benchmarks.git && cd detector_benchmarks
git clone https://eicweb.phy.anl.gov/EIC/benchmarks/detector_benchmarks.git && cd detector_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
......
benchmarks/pid/config.yml
+ 25
− 4
View file @ 87c352f1
......@@ -13,8 +13,30 @@ bench:mrich:
stage: benchmarks
needs: ["sim:backward"]
script:
- mkdir -p results/pid/backward/mrich/
- root -t -x -q -b "benchmarks/pid/scripts/mrich_analysis.cxx+(\"sim_output/sim_pid_backward_${PARTICLE}_5GeV.root\", \"${PARTICLE}\")"
- |
if [[ ${JUGGLER_DETECTOR_VERSION} =~ acadia ]]; then
echo "Performing MRICH benchmarks for acadia"
mkdir -p results/pid/backward/mrich/
root -t -x -q -b "benchmarks/pid/scripts/mrich_analysis.cxx+(\"sim_output/sim_pid_backward_${PARTICLE}_5GeV.root\", \"${PARTICLE}\")"
else
echo "Skipping MRICH Benchmarks for non-acadia detector version ${JUGGLER_DETECTOR_VERSION}"
fi
parallel:
matrix:
- PARTICLE: ["e-", "pi+", "proton"]
bench:erich:
extends: .det_benchmark
stage: benchmarks
needs: ["sim:backward"]
script:
- |
if [[ ${JUGGLER_DETECTOR_VERSION} =~ acadia ]]; then
echo "Skipping ERICH benchmarks for acadia"
else
echo "Doing MRICH Benchmarks for ${JUGGLER_DETECTOR_VERSION}"
echo "FIXME: needs implementing"
fi
parallel:
matrix:
- PARTICLE: ["e-", "pi+", "proton"]
......@@ -22,8 +44,7 @@ bench:mrich:
collect_results:pid:
extends: .det_benchmark
stage: collect
needs:
- ["bench:mrich"]
#needs:
script:
- ls -lrht
benchmarks/tracking_detectors/scripts/sim_track_hits.cxx
+ 113
− 162
View file @ 87c352f1
#include "ROOT/RDataFrame.hxx"
#include "TCanvas.h"
#include "TLegend.h"
#include "TH1D.h"
#include "TLegend.h"
#include "TProfile.h"
#include <cstdlib>
#include <iostream>
R__LOAD_LIBRARY(libeicd.so)
R__LOAD_LIBRARY(libDD4pod.so)
#include <fmt/format.h>
#include "dd4pod/Geant4ParticleCollection.h"
#include "dd4pod/TrackerHitCollection.h"
#include "dd4pod/TrackerHitData.h"
#include "dd4pod/Geant4ParticleCollection.h"
#include "eicd/TrackParametersCollection.h"
#include "eicd/ClusterCollection.h"
#include "eicd/ClusterData.h"
#include "eicd/TrackParametersCollection.h"
#include "eicd/TrackerHitCollection.h"
using ROOT::RDataFrame;
......@@ -23,13 +26,12 @@ using namespace ROOT::VecOps;
auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
std::vector<double> result;
for (size_t i = 0; i < in.size(); ++i) {
result.push_back(std::abs(1.0/(in[i].qOverP)));
result.push_back(std::abs(1.0 / (in[i].qOverP)));
}
return result;
};
std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){
std::vector<float> pt(std::vector<dd4pod::Geant4ParticleData> const& in) {
std::vector<float> result;
for (size_t i = 0; i < in.size(); ++i) {
result.push_back(std::sqrt(in[i].ps.x * in[i].ps.x + in[i].ps.y * in[i].ps.y));
......@@ -40,14 +42,14 @@ std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){
auto momentum = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
std::vector<double> result;
for (size_t i = 0; i < in.size(); ++i) {
result.push_back(in[i].E());
result.push_back(in[i].E());
}
return result;
};
auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) {
std::vector<double> result;
for (size_t i = 0; i < in.size(); ++i) {
result.push_back(in[i].Theta()*180/M_PI);
result.push_back(in[i].Theta() * 180 / M_PI);
}
return result;
};
......@@ -75,173 +77,122 @@ auto delta_p_over_p = [](const std::vector<double>& tracks, const std::vector<do
std::vector<double> res;
for (const auto& p1 : thrown) {
for (const auto& p2 : tracks) {
res.push_back((p1 - p2)/p1);
res.push_back((p1 - p2) / p1);
}
}
return res;
};
int sim_track_hits(const char* fname = "sim_track_hits.root")
{
// Create dataframe with extra definitions for hit multiplicities
ROOT::RDF::RNode add_subsystems(ROOT::RDF::RNode df, std::vector<std::pair<std::string, std::string>> hitcols) {
if (hitcols.empty()) {
return df;
}
const auto [name, collection] = hitcols.back();
std::cout << " --> registering subsystem " << collection << "\n";
auto df2 = df.Define("N_" + name, [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size(); }, {collection});
hitcols.pop_back();
return add_subsystems(df2, hitcols);
};
int sim_track_hits(const char* fname = "sim_track_hits.root") {
ROOT::EnableImplicitMT();
ROOT::RDataFrame df("events", fname);
// detect detector setup
std::string detector = "default";
if (const char* env_detector = std::getenv("JUGGLER_DETECTOR_VERSION")) {
if (detector.find("acadia") != std::string::npos) {
detector = "acadia";
} else if (detector.find("canyonlands") != std::string::npos) {
detector = "canyonlands";
}
}
std::cout << "sim_track_hits: detector set to " << detector << std::endl;
// minimal hit collection setup
std::vector<std::pair<std::string, std::string>> hitCollections{{"vtx_barrel", "VertexBarrelHits"},
{"trk_barrel", "TrackerBarrelHits"},
{"trk_endcap", "TrackerEndcapHits"},
{"gem_endcap", "GEMTrackerEndcapHits"}};
// append extra hit collections based on detector setup
if (detector == "acadia") {
hitCollections.push_back({"vtx_endcap", "VertexEndcapHits"});
} else if (detector == "canyonlands" || detector == "default") {
hitCollections.push_back({"mm_barrel", "MPGDTrackerBarrelHits"});
}
auto df0 = df.Define("isThrown", "mcparticles.genStatus == 1")
.Define("thrownParticles", "mcparticles[isThrown]")
.Define("thrownP", fourvec, {"thrownParticles"})
.Define("p_thrown", momentum, {"thrownP"})
.Define("theta_thrown", theta, {"thrownP"})
.Define("theta0", "theta_thrown[0]")
.Define("N_VertexBarrelHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"VertexBarrelHits"})
.Define("N_VertexEndcapHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"VertexEndcapHits"})
.Define("N_BarrelHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"TrackerBarrelHits"})
.Define("N_EndcapHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"TrackerEndcapHits"})
;
auto hBarrel_x_vs_y = df0.Histo2D(
{"hBarrel_x_vs_y", "; x ; y ", 100, -600, 600, 100, -600, 600},
"TrackerBarrelHits.position.x", "TrackerBarrelHits.position.y");
auto hEndcap_x_vs_y = df0.Histo2D(
{"hEndcap_x_vs_y", "; x ; y ", 100, -600, 600, 100, -600, 600},
"TrackerEndcapHits.position.x", "TrackerEndcapHits.position.y");
auto hVertexBarrel_x_vs_y = df0.Histo2D(
{"hVertexBarrel_x_vs_y", "; x ; y ", 100, -600, 600, 100, -600, 600},
"VertexBarrelHits.position.x", "VertexBarrelHits.position.y");
auto hVertexEndcap_x_vs_y = df0.Histo2D(
{"hVertexEndcap_x_vs_y", "; x ; y ", 100, -600, 600, 100, -600, 600},
"VertexEndcapHits.position.x", "VertexEndcapHits.position.y");
auto hBarrel_x_vs_z = df0.Histo2D(
{"hBarrel_x_vs_z", "; z ; x ", 200, -1200, 1200, 100, -900, 900},
"TrackerBarrelHits.position.z", "TrackerBarrelHits.position.x");
auto hEndcap_x_vs_z = df0.Histo2D(
{"hEndcap_x_vs_z", "; z ; x ", 200, -1200, 1200, 100, -900, 900},
"TrackerEndcapHits.position.z", "TrackerEndcapHits.position.x");
auto hVertexBarrel_x_vs_z = df0.Histo2D(
{"hVertexBarrel_x_vs_z", "; z ; x ", 200, -1200, 1200, 100, -900, 900},
"VertexBarrelHits.position.z", "VertexBarrelHits.position.x");
auto hVertexEndcap_x_vs_z = df0.Histo2D(
{"hVertexEndcap_x_vs_z", "; z ; x ", 200, -1200, 1200, 100, -900, 900},
"VertexEndcapHits.position.z", "VertexEndcapHits.position.x");
auto hBarrel_y_vs_z = df0.Histo2D(
{"hBarrel_y_vs_z", "; z ; x ", 200, -1200, 1200, 100, -900, 900},
"TrackerBarrelHits.position.z", "TrackerBarrelHits.position.y");
auto hEndcap_y_vs_z = df0.Histo2D(
{"hEndcap_y_vs_z", "; z ; x ", 200, -1200, 1200, 100, -900, 900},
"TrackerEndcapHits.position.z", "TrackerEndcapHits.position.y");
auto hVertexBarrel_y_vs_z = df0.Histo2D(
{"hVertexBarrel_y_vs_z", "; z ; x ", 200, -1200, 1200, 100, -900, 900},
"VertexBarrelHits.position.z", "VertexBarrelHits.position.y");
auto hVertexEndcap_y_vs_z = df0.Histo2D(
{"hVertexEndcap_y_vs_z", "; z ; x ", 200, -1200, 1200, 100, -900, 900},
"VertexEndcapHits.position.z", "VertexEndcapHits.position.y");
auto hBarrel_N_vs_theta = df0.Histo1D({"hBarrel_N_vs_theta", "; #theta [deg.]", 20, 0, 180 }, "theta0", "N_BarrelHits");
auto hEndcap_N_vs_theta = df0.Histo1D({"hEndcap_N_vs_theta", "; #theta [deg.]", 20, 0, 180 }, "theta0", "N_EndcapHits");
auto hVertexBarrel_N_vs_theta = df0.Histo1D({"hVertexBarrel_N_vs_theta", "; #theta [deg.]", 20, 0, 180 }, "theta0", "N_VertexBarrelHits");
auto hVertexEndcap_N_vs_theta = df0.Histo1D({"hVertexEndcap_N_vs_theta", "; #theta [deg.]", 20, 0, 180 }, "theta0", "N_VertexEndcapHits");
auto hBarrel_Nhits = df0.Histo1D({"hBarrel_Nhits", "; #theta [deg.]", 20, 0, 20 }, "N_BarrelHits");
auto hEndcap_Nhits = df0.Histo1D({"hEndcap_Nhits", "; #theta [deg.]", 20, 0, 20 }, "N_EndcapHits");
auto hVertexBarrel_Nhits = df0.Histo1D({"hVertexBarrel_Nhits", "; #theta [deg.]", 20, 0, 20 }, "N_VertexBarrelHits");
auto hVertexEndcap_Nhits = df0.Histo1D({"hVertexEndcap_Nhits", "; #theta [deg.]", 20, 0, 20 }, "N_VertexEndcapHits");
auto hBarrel_Ntheta = df0.Histo1D({"hBarrel_Ntheta", "; #theta [deg.]", 20, 0, 180 }, "theta0");
auto hEndcap_Ntheta = df0.Histo1D({"hEndcap_Ntheta", "; #theta [deg.]", 20, 0, 180 }, "theta0");
auto hVertexBarrel_Ntheta = df0.Histo1D({"hVertexBarrel_Ntheta", "; #theta [deg.]", 20, 0, 180 }, "theta0");
auto hVertexEndcap_Ntheta = df0.Histo1D({"hVertexEndcap_Ntheta", "; #theta [deg.]", 20, 0, 180 }, "theta0");
auto c = new TCanvas();
auto hs = new THStack("n_hits","; #theta ");
auto h1 = (TH1D*) hBarrel_N_vs_theta->Clone();
auto h2 = (TH1D*) hBarrel_Ntheta->Clone();
h1->Divide(h2);
hs->Add(h1);
h1 = (TH1D*) hEndcap_N_vs_theta->Clone();
h2 = (TH1D*) hEndcap_Ntheta->Clone();
h1->Divide(h2);
h1->SetLineColor(2);
hs->Add(h1);
h1 = (TH1D*) hVertexEndcap_N_vs_theta->Clone();
h2 = (TH1D*) hVertexEndcap_Ntheta->Clone();
h1->Divide(h2);
h1->SetLineColor(4);
hs->Add(h1);
h1 = (TH1D*) hVertexBarrel_N_vs_theta->Clone();
h2 = (TH1D*) hVertexBarrel_Ntheta->Clone();
h1->Divide(h2);
h1->SetLineColor(8);
hs->Add(h1);
hs->Draw("nostack, hist");
c->BuildLegend();
c->SaveAs("results/tracking_detectors/sim_track_hits_n_hits_vs_theta.png");
c->SaveAs("results/tracking_detectors/sim_track_hits_n_hits_vs_theta.pdf");
c = new TCanvas();
hs = new THStack("theta","; #theta ");
h1 = (TH1D*) hBarrel_N_vs_theta->Clone();
h2 = (TH1D*) hBarrel_Ntheta->Clone();
//h1->Divide(h2);
hs->Add(h2);
h1 = (TH1D*) hEndcap_N_vs_theta->Clone();
h2 = (TH1D*) hEndcap_Ntheta->Clone();
//h1->Divide(h2);
h1->SetLineColor(2);
h2->SetLineColor(2);
hs->Add(h2);
//h1 = (TH1D*) hVtxBarrel_vs_theta->Clone();
//h1->SetLineColor(4);
//h1->SetFillStyle(3001);
//h1->SetFillColor(4);
//hs->Add(h1);
hs->Draw("nostack hist");
c->BuildLegend();
c->SaveAs("results/tracking_detectors/sim_track_hits_theta.png");
c->SaveAs("results/tracking_detectors/sim_track_hits_theta.pdf");
c = new TCanvas();
hs = new THStack("hits","; hits ");
h1 = (TH1D*) hBarrel_Nhits->Clone();
hs->Add(h1);
h1 = (TH1D*) hEndcap_Nhits->Clone();
h1->SetLineColor(2);
h2->SetLineColor(2);
hs->Add(h2);
//h1 = (TH1D*) hVtxBarrel_Nhits->Clone();
//h1->SetLineColor(4);
//h1->SetFillStyle(3001);
//h1->SetFillColor(4);
//hs->Add(h1);
//hs->Draw("nostack hist");
c->BuildLegend();
c->SaveAs("results/tracking_detectors/sim_track_hits_nhits.png");
c->SaveAs("results/tracking_detectors/sim_track_hits_nhits.pdf");
c = new TCanvas();
hBarrel_x_vs_y->DrawCopy("colz");
c->SaveAs("results/tracking_detectors/sim_track_hits_trkBarrel_xy.png");
c->SaveAs("results/tracking_detectors/sim_track_hits_trkBarrel_xy.pdf");
c = new TCanvas();
hEndcap_x_vs_y->DrawCopy("colz");
hVertexEndcap_x_vs_y->DrawCopy("colz same");
//hBarrel_x_vs_y->DrawCopy("colz same");
//hvtxEndcap_x_vs_y->DrawCopy("colz same");
c->SaveAs("results/tracking_detectors/sim_track_hits_Hits_xy.png");
c->SaveAs("results/tracking_detectors/sim_track_hits_Hits_xy.pdf");
//hAllHits_x_vs_z->DrawCopy("colz");
hBarrel_x_vs_z->DrawCopy("colz");
hEndcap_x_vs_z->DrawCopy("colz same");
c->SaveAs("results/tracking_detectors/sim_track_hits_Hits_xz.png");
c->SaveAs("results/tracking_detectors/sim_track_hits_Hits_xz.pdf");
hBarrel_y_vs_z->DrawCopy("colz");
hEndcap_y_vs_z->DrawCopy("colz same");
c->SaveAs("results/tracking_detectors/sim_track_hits_Hits_yz.png");
c->SaveAs("results/tracking_detectors/sim_track_hits_Hits_yz.pdf");
.Define("theta0", "theta_thrown[0]");
auto df1 = add_subsystems(df0, hitCollections);
// get 1D and 2D histogram definitions
std::vector<decltype(df1.Histo2D({}, "", ""))> h2D;
std::vector<decltype(df1.Histo1D(""))> hnth;
std::vector<decltype(df1.Histo1D(""))> hn;
for (const auto& [name, col] : hitCollections) {
hnth.push_back(
df1.Histo1D({(name + "_nhits_vs_theta").c_str(), "; #theta [deg.]; ", 20, 0, 180}, "theta0", "N_" + name));
hn.push_back(df1.Histo1D({(name + "_nhits").c_str(), "; Nhits; #", 20, 0, 20}, "N_" + name));
h2D.push_back(df1.Histo2D({(name + "_x_vs_y").c_str(), "; x ; y ", 100, -600, 600, 100, -600, 600},
col + ".position.x", col + ".position.y"));
h2D.push_back(df1.Histo2D({(name + "_x_vs_z").c_str(), "; z ; x ", 100, -1200, 1200, 100, -900, 900},
col + ".position.z", col + ".position.x"));
}
auto hth = df1.Histo1D({"theta0", "; #theta [deg.]", 20, 0, 180}, "theta0");
// Draw multiplicity versus theta
{
TCanvas c;
THStack hs{"n_hits", "; #theta [deg.] "};
int idx = 1;
for (auto& h : hnth) {
h->Divide(&*hth);
h->SetLineColor(idx++);
hs.Add((TH1D*)h->Clone());
}
hs.Draw("nostack, hist");
c.BuildLegend();
c.SaveAs("results/tracking_detectors/sim_track_hits_n_hits_vs_theta.png");
c.SaveAs("results/tracking_detectors/sim_track_hits_n_hits_vs_theta.pdf");
}
// Draw nhits
{
TCanvas c;
THStack hs{"n_hits", "; NHits "};
int idx = 1;
for (auto& h : hn) {
h->SetLineColor(idx++);
hs.Add((TH1D*)h->Clone());
}
hs.Draw("nostack, hist");
c.BuildLegend();
c.SaveAs("results/tracking_detectors/sim_track_hits_n_hits.png");
c.SaveAs("results/tracking_detectors/sim_track_hits_n_hits.pdf");
}
// Draw total multiplicity versus theta
{
TCanvas c;
hth->DrawClone();
c.BuildLegend();
c.SaveAs("results/tracking_detectors/sim_track_hits_theta.png");
c.SaveAs("results/tracking_detectors/sim_track_hits_theta.pdf");
}
for (auto& h : h2D) {
TCanvas c;
h->DrawClone("colz");
c.SaveAs(fmt::format("results/tracking_detectors/sim_track_hits_{}.png", h->GetName()).c_str());
c.SaveAs(fmt::format("results/tracking_detectors/sim_track_hits_{}.pdf", h->GetName()).c_str());
}
return 0;
}
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