Commit e62205fe authored by Viktor's avatar Viktor Committed by Florian Uhlig
Browse files

add some default cuts for QA

parent 43a7c5cf
set(ANALYSISTREE_VERSION 16de36dc21744aa553a9470084af9c34da0d379b) #v2.1.0
set(ANALYSISTREE_VERSION e5a1f89993ecce314371d64bc921485ac4888511) #v2.1.0
set(ANALYSISTREE_SRC_URL "https://github.com/HeavyIonAnalysis/AnalysisTree.git")
set(ANALYSISTREE_DESTDIR "${CMAKE_BINARY_DIR}/external/ANALYSISTREE-prefix")
......
AnalysisTree::Cuts* GetCbmEventCuts(const std::string& branch, std::string name="CbmGoodEvent") {
AnalysisTree::SimpleCut vtx_x_cut = AnalysisTree::RangeCut(branch+".vtx_x", -0.5, 0.5);
AnalysisTree::SimpleCut vtx_y_cut = AnalysisTree::RangeCut(branch+".vtx_y", -0.5, 0.5);
AnalysisTree::SimpleCut vtx_z_cut = AnalysisTree::RangeCut(branch+".vtx_z", -0.1, 0.1);
AnalysisTree::SimpleCut vtx_chi2_cut = AnalysisTree::RangeCut(branch+".vtx_chi2", 0.8, 1.7);
auto* event_cuts = new AnalysisTree::Cuts(std::move(name), {vtx_x_cut, vtx_y_cut, vtx_z_cut, vtx_chi2_cut});
return event_cuts;
}
AnalysisTree::Cuts* GetCbmTrackCuts(const std::string& branch, std::string name="CbmGoodVertexTrack") {
AnalysisTree::SimpleCut vtx_chi2_track_cut = AnalysisTree::RangeCut(branch+".vtx_chi2", 0, 3);
AnalysisTree::SimpleCut nhits_cut = AnalysisTree::RangeCut(branch+".nhits", 4, 100);
AnalysisTree::SimpleCut chi2_cut({branch + ".chi2", branch + ".ndf"}, [](std::vector<double> par){ return par[0]/par[1] < 3; });
AnalysisTree::SimpleCut eta_cut = AnalysisTree::RangeCut(branch+".eta", 0.2, 6);
auto* vertex_tracks_cuts = new AnalysisTree::Cuts(std::move(name), {vtx_chi2_track_cut, nhits_cut, chi2_cut, eta_cut});
return vertex_tracks_cuts;
}
AnalysisTree::Cuts* GetCbmTofHitsCuts(const std::string& branch, std::string name="CbmGoodTofHit") {
AnalysisTree::SimpleCut tof_cuts_dx = AnalysisTree::RangeCut(branch+".dx", -2, 2);
AnalysisTree::SimpleCut tof_cuts_dy = AnalysisTree::RangeCut(branch+ ".dy", -1, 1);
auto* tof_cuts = new AnalysisTree::Cuts(std::move(name), {tof_cuts_dx, tof_cuts_dy});
return tof_cuts;
}
AnalysisTree::Cuts* GetCbmMcTracksCuts(const std::string& branch, std::string name="CbmMcPrimaryTrack") {
auto* sim_tracks_cut = new AnalysisTree::Cuts(std::move(name), {AnalysisTree::EqualsCut(branch+".mother_id", -1)});
return sim_tracks_cut;
}
\ No newline at end of file
//#include <AnalysisTree/TaskManager.hpp>
//#include "AnalysisTree/Cuts.hpp"
//
//#include "src/EntryConfig.hpp"
//#include "src/Task.hpp"
//#include "src/Utils.hpp"
#include "../cuts/cbm_cuts.h"
typedef AnalysisTree::QA::EntryConfig::PlotType PlotType;
using AnalysisTree::QA::gNbins;
......@@ -30,20 +32,8 @@ const std::string tof_hits = "TofHits";
const std::string trd_tracks = "TrdTracks";
const std::string rich_rings = "RichRings";
int main(int argc, char** argv) {
if (argc <= 1) {
std::cout << "Not enough arguments! Please use:" << std::endl;
std::cout << " ./cbm_qa filelist" << std::endl;
return -1;
}
const std::string filelist = argv[1];
run_analysistree_qa(filelist);
return 0;
}
void run_analysistree_qa(std::string filelist, bool is_single_file) {
void run_analysistree_qa(std::string filelist, bool is_single_file)
{
if (is_single_file) {
std::ofstream fl("fl_temp.txt");
fl << filelist << "\n";
......@@ -56,6 +46,10 @@ void run_analysistree_qa(std::string filelist, bool is_single_file) {
auto* task = new QA::Task;
task->SetOutputFileName("cbm_qa.root");
task->SetEventCuts(GetCbmEventCuts("RecEventHeader"));
task->AddBranchCut(GetCbmMcTracksCuts("SimParticles"));
task->AddBranchCut(GetCbmTrackCuts("RecTracks"));
RichRingsQA(*task);
TrdTracksQA(*task);
// KFPFTracksQA(*task);
......@@ -80,61 +74,39 @@ void run_analysistree_qa(std::string filelist, bool is_single_file) {
}
}
void TrdTracksQA(QA::Task& task) {
void TrdTracksQA(QA::Task& task)
{
AddTrackQA(&task, trd_tracks);
AddTracksMatchQA(&task, trd_tracks, rec_tracks);
task.AddH1({"TRD total energy loss, keV",
{trd_tracks, "energy_loss_total"},
{gNbins, 0, 1}});
task.AddH1({"TRD energy loss in 1st station, keV",
{trd_tracks, "energy_loss_0"},
{gNbins, 0, 1}});
task.AddH1({"TRD energy loss in 2nd station",
{trd_tracks, "energy_loss_1"},
{gNbins, 0, 1}});
task.AddH1({"TRD energy loss in 3rd station",
{trd_tracks, "energy_loss_2"},
{gNbins, 0, 1}});
task.AddH1({"TRD energy loss in 4th station",
{trd_tracks, "energy_loss_3"},
{gNbins, 0, 1}});
task.AddH1({"TRD energy loss in 1st station, keV", {trd_tracks, "energy_loss_0"}, {gNbins, 0, 1}});
task.AddH1({"TRD energy loss in 2nd station", {trd_tracks, "energy_loss_1"}, {gNbins, 0, 1}});
task.AddH1({"TRD energy loss in 3rd station", {trd_tracks, "energy_loss_2"}, {gNbins, 0, 1}});
task.AddH1({"TRD energy loss in 4th station", {trd_tracks, "energy_loss_3"}, {gNbins, 0, 1}});
task.AddH1({"Number of hits in TRD", {trd_tracks, "n_hits"}, {6, 0, 6}});
task.AddH1({"PID ANN", {trd_tracks, "pid_ann"}, {gNbins, -1, 1}});
// task.AddH1({"PID WKN", {trd_tracks, "pid_wkn"}, {gNbins, -1, 1}});
task.AddH1(
{"PID Likelihood, e^{-}", {trd_tracks, "pid_like_e"}, {gNbins, -1, 1}});
task.AddH1(
{"PID Likelihood, #pi", {trd_tracks, "pid_like_pi"}, {gNbins, -1, 1}});
task.AddH1(
{"PID Likelihood, K", {trd_tracks, "pid_like_k"}, {gNbins, -1, 1}});
task.AddH1(
{"PID Likelihood, p", {trd_tracks, "pid_like_p"}, {gNbins, -1, 1}});
task.AddH1(
{"PID Likelihood, #mu", {trd_tracks, "pid_like_mu"}, {gNbins, -1, 1}});
task.AddH1({"PID Likelihood, e^{-}", {trd_tracks, "pid_like_e"}, {gNbins, -1, 1}});
task.AddH1({"PID Likelihood, #pi", {trd_tracks, "pid_like_pi"}, {gNbins, -1, 1}});
task.AddH1({"PID Likelihood, K", {trd_tracks, "pid_like_k"}, {gNbins, -1, 1}});
task.AddH1({"PID Likelihood, p", {trd_tracks, "pid_like_p"}, {gNbins, -1, 1}});
task.AddH1({"#chi^{2}/NDF", {trd_tracks, "chi2_ov_ndf"}, {gNbins, 0, 30}});
task.AddH1({"p_{T}^{last} (GeV/c)", {trd_tracks, "pT_out"}, {gNbins, 0, 10}});
task.AddH1({"p^{last} (GeV/c)", {trd_tracks, "p_out"}, {gNbins, 0, 10}});
}
void RichRingsQA(QA::Task& task) {
task.AddH1(
{"RICh ring x-position (cm)", {rich_rings, "x"}, {gNbins, -100, 100}});
task.AddH1(
{"RICh ring y-position (cm)", {rich_rings, "y"}, {gNbins, -250, 250}});
task.AddH2(
{"RICh ring x-position (cm)", {rich_rings, "x"}, {gNbins, -100, 100}},
{"RICh ring y-position (cm)", {rich_rings, "y"}, {gNbins, -250, 250}});
void RichRingsQA(QA::Task& task)
{
task.AddH1({"RICh ring x-position (cm)", {rich_rings, "x"}, {gNbins, -100, 100}});
task.AddH1({"RICh ring y-position (cm)", {rich_rings, "y"}, {gNbins, -250, 250}});
task.AddH2({"RICh ring x-position (cm)", {rich_rings, "x"}, {gNbins, -100, 100}},
{"RICh ring y-position (cm)", {rich_rings, "y"}, {gNbins, -250, 250}});
task.AddH1({"Ring radius", {rich_rings, "radius"}, {gNbins, 0, 10}});
task.AddH1({"n_hits", {rich_rings, "n_hits"}, {gNbins, 0, 60}});
task.AddH1(
{"n_hits_on_ring", {rich_rings, "n_hits_on_ring"}, {gNbins, 0, 100}});
task.AddH1({"n_hits_on_ring", {rich_rings, "n_hits_on_ring"}, {gNbins, 0, 100}});
task.AddH1({"axis_a", {rich_rings, "axis_a"}, {gNbins, 0, 10}});
task.AddH1({"axis_b", {rich_rings, "axis_b"}, {gNbins, 0, 10}});
......@@ -142,16 +114,13 @@ void RichRingsQA(QA::Task& task) {
task.AddH1({"phi_ellipse", {rich_rings, "phi_ellipse"}, {gNbins, -3.2, 3.2}});
}
void VertexTracksQA(QA::Task& task) {
void VertexTracksQA(QA::Task& task)
{
AddTrackQA(&task, rec_tracks);
if (!sim_particles.empty()) {
AddTracksMatchQA(&task, rec_tracks, sim_particles);
}
if (!sim_particles.empty()) { AddTracksMatchQA(&task, rec_tracks, sim_particles); }
Variable chi2_over_ndf(
"chi2_ndf",
{{rec_tracks, "chi2"}, {rec_tracks, "ndf"}},
[](std::vector<double>& var) { return var.at(0) / var.at(1); });
Variable chi2_over_ndf("chi2_ndf", {{rec_tracks, "chi2"}, {rec_tracks, "ndf"}},
[](std::vector<double>& var) { return var.at(0) / var.at(1); });
task.AddH1({"DCA_{x}, cm", {rec_tracks, "dcax"}, {gNbins, -1, 1}});
task.AddH1({"DCA_{y}, cm", {rec_tracks, "dcay"}, {gNbins, -1, 1}});
......@@ -165,7 +134,8 @@ void VertexTracksQA(QA::Task& task) {
{"DCA_{y}, cm", {rec_tracks, "dcay"}, {gNbins, -1, 1}});
}
void TofHitsQA(QA::Task& task) {
void TofHitsQA(QA::Task& task)
{
task.AddH1({"TOF hit x-position (cm)", {tof_hits, "x"}, {gNbins, -600, 600}});
task.AddH1({"TOF hit y-position (cm)", {tof_hits, "y"}, {gNbins, -400, 400}});
task.AddH1({"TOF hit z-position (cm)", {tof_hits, "z"}, {gNbins, 600, 800}});
......@@ -173,42 +143,52 @@ void TofHitsQA(QA::Task& task) {
task.AddH2({"TOF hit x-position (cm)", {tof_hits, "x"}, {gNbins, -600, 600}},
{"TOF hit y-position (cm)", {tof_hits, "y"}, {gNbins, -600, 600}});
Variable qp_sts("qp_reco",
{{rec_tracks, "q"}, {rec_tracks, "p"}},
Variable qp_sts("qp_reco", {{rec_tracks, "q"}, {rec_tracks, "p"}},
[](std::vector<double>& qp) { return qp.at(0) * qp.at(1); });
task.AddH2({"sign(q)*p, GeV/c", qp_sts, {gNbins, -10, 10}},
{"m^{2}, GeV^{2}/c^{2}", {tof_hits, "mass2"}, {gNbins, -1, 2}});
task.AddH2({"sign(q)*p STS, GeV/c", qp_sts, {gNbins, -10, 10}},
{"sign(q)*p TOF, GeV/c", {tof_hits, "qp_tof"}, {gNbins, -10, 10}});
SimpleCut sc_protons = EqualsCut("SimParticles.pid", 2212);
SimpleCut sc_prim = EqualsCut("SimParticles.mother_id", -1);
auto prim = new Cuts("mc_primary", {sc_prim});
auto protons = new Cuts("mc_protons", {sc_protons});
auto prim_protons = new Cuts("mc_primary_protons", {sc_protons, sc_prim});
auto prim_pions = new Cuts("mc_primary_pions_pos", {EqualsCut("SimParticles.pid", 211), sc_prim});
auto prim_kaons = new Cuts("mc_primary_kaons_pos", {EqualsCut("SimParticles.pid", 321), sc_prim});
std::vector<Cuts*> cuts = {nullptr, prim, protons, prim_protons, prim_pions, prim_kaons};
for(auto cut : cuts){
task.AddH2({"p_{MC}, GeV/c", {"SimParticles", "p"}, {250, 0, 10}},{"m^{2}, GeV^{2}/c^{2}", {"TofHits", "mass2"}, {500, -1, 2}}, cut);
}
}
void SimParticlesQA(QA::Task& task) {
void SimParticlesQA(QA::Task& task)
{
AddParticleQA(&task, sim_particles);
// task.AddH1({"PDG code", {sim_particles, "pid"}, {8001, -3999.5, 4000.5}});
}
void SimEventHeaderQA(QA::Task& task) {
task.AddH1(
{"x_{vertex}^{MC} (cm)", {sim_event_header, "vtx_x"}, {gNbins, -1, 1}});
task.AddH1(
{"y_{vertex}^{MC} (cm)", {sim_event_header, "vtx_y"}, {gNbins, -1, 1}});
task.AddH1(
{"z_{vertex}^{MC} (cm)", {sim_event_header, "vtx_z"}, {gNbins, -1, 1}});
void SimEventHeaderQA(QA::Task& task)
{
task.AddH1({"x_{vertex}^{MC} (cm)", {sim_event_header, "vtx_x"}, {gNbins, -1, 1}});
task.AddH1({"y_{vertex}^{MC} (cm)", {sim_event_header, "vtx_y"}, {gNbins, -1, 1}});
task.AddH1({"z_{vertex}^{MC} (cm)", {sim_event_header, "vtx_z"}, {gNbins, -1, 1}});
task.AddH1({"b (fm)", {sim_event_header, "b"}, {gNbins, 0, 20}});
task.AddH1({"#Psi_{RP}", {sim_event_header, "psi_RP"}, {gNbins, 0, 6.5}});
task.AddH2(
{"x_{vertex}^{MC} (cm)", {sim_event_header, "vtx_x"}, {gNbins, -1, 1}},
{"y_{vertex}^{MC} (cm)", {sim_event_header, "vtx_y"}, {gNbins, -1, 1}});
task.AddH2({"x_{vertex}^{MC} (cm)", {sim_event_header, "vtx_x"}, {gNbins, -1, 1}},
{"y_{vertex}^{MC} (cm)", {sim_event_header, "vtx_y"}, {gNbins, -1, 1}});
}
void RecEventHeaderQA(QA::Task& task) {
void RecEventHeaderQA(QA::Task& task)
{
task.AddH1({"x_{vertex} (cm)", {rec_event_header, "vtx_x"}, {gNbins, -1, 1}});
task.AddH1({"y_{vertex} (cm)", {rec_event_header, "vtx_y"}, {gNbins, -1, 1}});
task.AddH1({"z_{vertex} (cm)", {rec_event_header, "vtx_z"}, {gNbins, -1, 1}});
task.AddH1(
{"#chi^{2}_{vertex fit}", {rec_event_header, "vtx_chi2"}, {gNbins, 0, 5}});
task.AddH1({"#chi^{2}_{vertex fit}", {rec_event_header, "vtx_chi2"}, {gNbins, 0, 5}});
task.AddH1({"E_{PSD} (GeV)", {rec_event_header, "Epsd"}, {gNbins, 0, 60}});
task.AddH1({"M_{tracks}", {rec_event_header, "M"}, {800, 0, 800}});
......@@ -225,73 +205,47 @@ void RecEventHeaderQA(QA::Task& task) {
{"z_{vertex} (cm)", {sim_event_header, "vtx_z"}, {gNbins, -1, 1}});
}
void EfficiencyMaps(QA::Task& task) {
void EfficiencyMaps(QA::Task& task)
{
const float y_beam = 1.62179f; // TODO from DataHeader
const float p_mass = 0.938;
const float pi_mass = 0.14;
Variable proton_y("y-y_{beam}",
{{rec_tracks, "p"}, {rec_tracks, "pz"}},
[y_beam, p_mass](std::vector<double>& var) {
const float e = sqrt(p_mass * p_mass + var[0] * var[0]);
return 0.5 * log((e + var[1]) / (e - var[1]));
});
Variable pion_y("y-y_{beam}",
{{rec_tracks, "p"}, {rec_tracks, "pz"}},
[y_beam, pi_mass](std::vector<double>& var) {
const float e = sqrt(pi_mass * pi_mass + var[0] * var[0]);
return 0.5 * log((e + var[1]) / (e - var[1]));
});
Cuts* mc_protons =
new Cuts("McProtons", {EqualsCut({sim_particles + ".pid"}, 2212)});
Cuts* mc_pions_neg =
new Cuts("McPionsNeg", {EqualsCut({sim_particles + ".pid"}, -211)});
Cuts* mc_pions_pos =
new Cuts("McPionsPos", {EqualsCut({sim_particles + ".pid"}, -211)});
Variable proton_y("y-y_{beam}", {{rec_tracks, "p"}, {rec_tracks, "pz"}}, [y_beam, p_mass](std::vector<double>& var) {
const float e = sqrt(p_mass * p_mass + var[0] * var[0]);
return 0.5 * log((e + var[1]) / (e - var[1]));
});
Variable pion_y("y-y_{beam}", {{rec_tracks, "p"}, {rec_tracks, "pz"}}, [y_beam, pi_mass](std::vector<double>& var) {
const float e = sqrt(pi_mass * pi_mass + var[0] * var[0]);
return 0.5 * log((e + var[1]) / (e - var[1]));
});
Cuts* mc_protons = new Cuts("McProtons", {EqualsCut({sim_particles + ".pid"}, 2212)});
Cuts* mc_pions_neg = new Cuts("McPionsNeg", {EqualsCut({sim_particles + ".pid"}, -211)});
Cuts* mc_pions_pos = new Cuts("McPionsPos", {EqualsCut({sim_particles + ".pid"}, -211)});
task.AddH2({"#it{y}_{Lab}", {sim_particles, "rapidity"}, {gNbins, -1, 5}},
{"p_{T}, GeV/c", {sim_particles, "pT"}, {gNbins, 0, 2}},
mc_protons);
task.AddH2({"#it{y}_{Lab}", proton_y, {gNbins, -1, 5}},
{"p_{T}, GeV/c", {rec_tracks, "pT"}, {gNbins, 0, 2}},
{"p_{T}, GeV/c", {sim_particles, "pT"}, {gNbins, 0, 2}}, mc_protons);
task.AddH2({"#it{y}_{Lab}", proton_y, {gNbins, -1, 5}}, {"p_{T}, GeV/c", {rec_tracks, "pT"}, {gNbins, 0, 2}},
mc_protons);
task.AddH2({"#it{y}_{Lab}", {sim_particles, "rapidity"}, {gNbins, -1, 5}},
{"p_{T}, GeV/c", {sim_particles, "pT"}, {gNbins, 0, 2}},
mc_pions_neg);
task.AddH2({"#it{y}_{Lab}", pion_y, {gNbins, -1, 5}},
{"p_{T}, GeV/c", {rec_tracks, "pT"}, {gNbins, 0, 2}},
{"p_{T}, GeV/c", {sim_particles, "pT"}, {gNbins, 0, 2}}, mc_pions_neg);
task.AddH2({"#it{y}_{Lab}", pion_y, {gNbins, -1, 5}}, {"p_{T}, GeV/c", {rec_tracks, "pT"}, {gNbins, 0, 2}},
mc_pions_neg);
task.AddH2({"#it{y}_{Lab}", {sim_particles, "rapidity"}, {gNbins, -1, 5}},
{"p_{T}, GeV/c", {sim_particles, "pT"}, {gNbins, 0, 2}},
mc_pions_pos);
task.AddH2({"#it{y}_{Lab}", pion_y, {gNbins, -1, 5}},
{"p_{T}, GeV/c", {rec_tracks, "pT"}, {gNbins, 0, 2}},
{"p_{T}, GeV/c", {sim_particles, "pT"}, {gNbins, 0, 2}}, mc_pions_pos);
task.AddH2({"#it{y}_{Lab}", pion_y, {gNbins, -1, 5}}, {"p_{T}, GeV/c", {rec_tracks, "pT"}, {gNbins, 0, 2}},
mc_pions_pos);
}
void KFPFTracksQA(QA::Task& task) {
const std::vector<std::string> field_par_names_ {
"cx0", "cx1", "cx2", "cy1", "cy2", "cz0", "cz1", "cz2"};
const std::vector<std::string> cov_names_ {"cov1",
"cov2",
"cov3",
"cov4",
"cov5",
"cov6",
"cov7",
"cov8",
"cov9",
"cov10",
"cov11",
"cov12",
"cov13",
"cov14",
"cov15"};
void KFPFTracksQA(QA::Task& task)
{
const std::vector<std::string> field_par_names_ {"cx0", "cx1", "cx2", "cy1", "cy2", "cz0", "cz1", "cz2"};
const std::vector<std::string> cov_names_ {"cov1", "cov2", "cov3", "cov4", "cov5", "cov6", "cov7", "cov8",
"cov9", "cov10", "cov11", "cov12", "cov13", "cov14", "cov15"};
// KFPF tracks
task.AddH1({"x, cm", {rec_tracks, "x"}, {gNbins, -10, 10}});
task.AddH1({"y, cm", {rec_tracks, "y"}, {gNbins, -10, 10}});
......@@ -309,3 +263,16 @@ void KFPFTracksQA(QA::Task& task) {
task.AddH1({"cy0", {rec_tracks, "cy0"}, {gNbins, -12, -8}});
task.AddH1({"z0", {rec_tracks, "z0"}, {gNbins, 0, 40}});
}
int main(int argc, char** argv)
{
if (argc <= 1) {
std::cout << "Not enough arguments! Please use:" << std::endl;
std::cout << " ./cbm_qa filelist" << std::endl;
return -1;
}
const std::string filelist = argv[1];
run_analysistree_qa(filelist);
return 0;
}
\ No newline at end of file
Supports Markdown
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