diff --git a/external/InstallAnalysisTree.cmake b/external/InstallAnalysisTree.cmake index 3ca9c626ebaf83300668bb15f35198c8a6f5d732..d3d20d6454568781ecb42aafcdeccf73e774e2fa 100644 --- a/external/InstallAnalysisTree.cmake +++ b/external/InstallAnalysisTree.cmake @@ -1,4 +1,4 @@ -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") diff --git a/macro/analysis/common/cuts/cbm_cuts.h b/macro/analysis/common/cuts/cbm_cuts.h new file mode 100644 index 0000000000000000000000000000000000000000..7b930dbee960ff6507221c8f2362dcb5aa2f4382 --- /dev/null +++ b/macro/analysis/common/cuts/cbm_cuts.h @@ -0,0 +1,33 @@ + +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 diff --git a/macro/analysis/common/qa/run_analysistree_qa.C b/macro/analysis/common/qa/run_analysistree_qa.C index ec6d9f8b1aaacecd8a1e0f3a192ac83575521568..d3633886c69bf2e04128ebb9e6e864d0df79dd76 100644 --- a/macro/analysis/common/qa/run_analysistree_qa.C +++ b/macro/analysis/common/qa/run_analysistree_qa.C @@ -1,10 +1,12 @@ //#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