diff --git a/Detectors/Upgrades/ITS3/CMakeLists.txt b/Detectors/Upgrades/ITS3/CMakeLists.txt index 73ad4b9d53e37..5e40e59ad0068 100644 --- a/Detectors/Upgrades/ITS3/CMakeLists.txt +++ b/Detectors/Upgrades/ITS3/CMakeLists.txt @@ -19,3 +19,4 @@ add_subdirectory(base) add_subdirectory(workflow) add_subdirectory(reconstruction) add_subdirectory(macros) +add_subdirectory(study) diff --git a/Detectors/Upgrades/ITS3/macros/test/CMakeLists.txt b/Detectors/Upgrades/ITS3/macros/test/CMakeLists.txt index cb6812445283c..6b274e764f276 100644 --- a/Detectors/Upgrades/ITS3/macros/test/CMakeLists.txt +++ b/Detectors/Upgrades/ITS3/macros/test/CMakeLists.txt @@ -12,7 +12,6 @@ its3_add_macro(CheckDigitsITS3.C) its3_add_macro(CheckClustersITS3.C) its3_add_macro(CheckTracksITS3.C) -its3_add_macro(CheckDCA.C) its3_add_macro(CreateDictionariesITS3.C) its3_add_macro(buildMatBudLUT.C) its3_add_macro(CheckHits.C) diff --git a/Detectors/Upgrades/ITS3/macros/test/CheckClustersITS3.C b/Detectors/Upgrades/ITS3/macros/test/CheckClustersITS3.C index f245a047377ae..5e56321c7676d 100644 --- a/Detectors/Upgrades/ITS3/macros/test/CheckClustersITS3.C +++ b/Detectors/Upgrades/ITS3/macros/test/CheckClustersITS3.C @@ -283,7 +283,7 @@ void CheckClustersITS3(const std::string& clusfile = "o2clus_its.root", nt.Draw("cgy:cgx>>h_cgy_vs_cgx_OB(1000, -50, 50, 1000, -50, 50)", "id >= 3456", "colz"); canvCgXCgY->cd(4); nt.Draw("cgy:cgz>>h_cgy_vs_cgz_OB(1000, -100, 100, 1000, -50, 50)", "id >= 3456", "colz"); - canvCgXCgY->SaveAs("it3clusters_y_vs_x_vs_z.pdf"); + canvCgXCgY->SaveAs("it3clusters_y_vs_x_vs_z.png"); auto canvdXdZ = new TCanvas("canvdXdZ", "", 1600, 800); canvdXdZ->Divide(2, 2); @@ -295,7 +295,7 @@ void CheckClustersITS3(const std::string& clusfile = "o2clus_its.root", nt.Draw("dx:dz>>h_dx_vs_dz_IB_z(1000, -0.01, 0.01, 1000, -0.01, 0.01)", "id < 3456 && abs(cgz) < 2", "colz"); canvdXdZ->cd(4)->SetLogz(); nt.Draw("dx:dz>>h_dx_vs_dz_OB_z(1000, -0.01, 0.01, 1000, -0.01, 0.01)", "id >= 3456 && abs(cgz) < 2", "colz"); - canvdXdZ->SaveAs("it3clusters_dx_vs_dz.pdf"); + canvdXdZ->SaveAs("it3clusters_dx_vs_dz.png"); auto canvCHXZ = new TCanvas("canvCHXZ", "", 1600, 1600); canvCHXZ->Divide(2, 2); @@ -307,7 +307,7 @@ void CheckClustersITS3(const std::string& clusfile = "o2clus_its.root", nt.Draw("(cgz-hgz)*10000:eta>>h_chz_IB(101,-1.4,1.4,101,-50,50)", "id<3456", "prof"); canvCHXZ->cd(4); nt.Draw("(cgz-hgz)*10000:eta>>h_chz_OB(101,-1.4,1.4,101,-50,50)", "id>=3456", "prof"); - canvCgXCgY->SaveAs("it3clusters_xz_eta.pdf"); + canvCgXCgY->SaveAs("it3clusters_xz_eta.png"); auto c1 = new TCanvas("p1", "pullX"); c1->cd(); diff --git a/Detectors/Upgrades/ITS3/macros/test/CheckDCA.C b/Detectors/Upgrades/ITS3/macros/test/CheckDCA.C deleted file mode 100644 index b2872431384f1..0000000000000 --- a/Detectors/Upgrades/ITS3/macros/test/CheckDCA.C +++ /dev/null @@ -1,965 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. - -/// \file CheckDCA.C -/// \brief Simple macro to check ITS3 impact parameter resolution - -#if !defined(__CLING__) || defined(__ROOTCLING__) -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include "DataFormatsITS/TrackITS.h" -#include "DataFormatsTPC/TrackTPC.h" -#include "DetectorsBase/Propagator.h" -#include "Field/MagneticField.h" -#include "ITSBase/GeometryTGeo.h" -#include "DetectorsBase/Propagator.h" -#include "ReconstructionDataFormats/TrackTPCITS.h" -#include "ReconstructionDataFormats/Vertex.h" -#include "ReconstructionDataFormats/DCA.h" -#include "DataFormatsITSMFT/CompCluster.h" -#include "SimulationDataFormat/MCCompLabel.h" -#include "SimulationDataFormat/MCEventHeader.h" -#include "SimulationDataFormat/MCTrack.h" -#include "SimulationDataFormat/MCTruthContainer.h" -#include "SimulationDataFormat/MCUtils.h" -#include "SimulationDataFormat/TrackReference.h" -#include "Steer/MCKinematicsReader.h" - -#include -#include -#include -#include -#include -#include -#include -#endif - -namespace fs = std::filesystem; - -constexpr auto mMatCorr{o2::base::Propagator::MatCorrType::USEMatCorrNONE}; -constexpr float mMaxStep{2}; - -constexpr float rapMax{0.9}; - -std::vector find_dirs(fs::path const& dir, std::function filter, std::optional> sort = std::nullopt) -{ - std::vector result; - if (fs::exists(dir)) { // Find Dirs matching filter - for (auto const& entry : fs::recursive_directory_iterator(dir, fs::directory_options::follow_directory_symlink)) { - if (fs::is_directory(entry) && filter(entry)) { - result.emplace_back(entry); - } - } - } - if (sort) { // Optionally sort paths - std::sort(result.begin(), result.end(), *sort); - } - return result; -} - -void CheckDCA(const std::string& collisioncontextFileName = "collisioncontext.root", - const std::string& tpcTracksFileName = "tpctracks.root", - const std::string& itsTracksFileName = "o2trac_its.root", - const std::string& itstpcTracksFileName = "o2match_itstpc.root", - const std::string& magFileName = "o2sim_grp.root") -{ - gROOT->SetBatch(); - gStyle->SetOptStat(0); - gStyle->SetPalette(kRainBow); - gStyle->SetPadLeftMargin(0.16); - gStyle->SetPadTickX(1); - gStyle->SetPadTickY(1); - gErrorIgnoreLevel = 2001; // suppress warnings - ProcInfo_t procInfo; - - const int nPtBins = 35; - const int nPtBinsEff = 39; - double ptLimits[nPtBins] = {0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2., 2.2, 2.5, 3., 4., 5., 6., 8., 10., 15., 20.}; - double ptLimitsEff[nPtBinsEff] = {0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2., 2.2, 2.5, 3., 4., 5., 6., 8., 10., 15., 20.}; - - const std::regex tf_pattern(R"(tf\d+)"); - auto tf_matcher = [&tf_pattern](fs::path const& p) -> bool { - return std::regex_search(p.string(), tf_pattern); - }; - auto tf_sorter = [&tf_pattern](fs::path const& a, fs::path const& b) -> bool { - const auto &as = a.string(), &bs = b.string(); - std::smatch am, bm; - if (std::regex_search(as, am, tf_pattern) && std::regex_search(bs, bm, tf_pattern)) { - return std::stoi(am.str().substr(2)) < std::stoi(bm.str().substr(2)); - } else { - LOGP(fatal, "TF Regex matching failed"); - return false; - } - }; - - const int nSpecies = 4; - std::array pdgCodes{11, 211, 321, 2212}; - auto fGaus = new TF1("fGaus", "gaus", -200., 200.); - std::map partNames = { - {11, "Electrons"}, - {211, "Pions"}, - {321, "Kaons"}, - {2212, "Protons"}}; - std::map colors{{11, kOrange + 7}, {211, kRed + 1}, {321, kAzure + 4}, {2212, kGreen + 2}}; - /// ITS - std::map hDcaxyResAllLayersITS = { - {11, new TH1F("hDcaxyResElectronsAllLayersITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}, - {211, new TH1F("hDcaxyResPionsAllLayersITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}, - {321, new TH1F("hDcaxyResKaonsAllLayersITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}, - {2212, new TH1F("hDcaxyResProtonsAllLayersITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}}; - std::map hDcazResAllLayersITS = { - {11, new TH1F("hDcazResElectronsAllLayersITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}, - {211, new TH1F("hDcazResPionsAllLayersITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}, - {321, new TH1F("hDcazResKaonsAllLayersITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}, - {2212, new TH1F("hDcazResProtonsAllLayersITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}}; - std::map hPtResAllLayersITS = { - {11, new TH1F("hPtResElectronsAllLayersITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(#Delta#it{p}_{T}/#it{p}_{T})", nPtBins - 1, ptLimits)}, - {211, new TH1F("hPtResPionsAllLayersITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(#Delta#it{p}_{T}/#it{p}_{T})", nPtBins - 1, ptLimits)}, - {321, new TH1F("hPtResKaonsAllLayersITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(#Delta#it{p}_{T}/#it{p}_{T})", nPtBins - 1, ptLimits)}, - {2212, new TH1F("hPtResProtonsAllLayersITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(#Delta#it{p}_{T}/#it{p}_{T})", nPtBins - 1, ptLimits)}}; - std::map hDcaxyResNoFirstLayerITS = { - {11, new TH1F("hDcaxyResElectronsNoFirstLayerITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}, - {211, new TH1F("hDcaxyResPionsNoFirstLayerITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}, - {321, new TH1F("hDcaxyResKaonsNoFirstLayerITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}, - {2212, new TH1F("hDcaxyResProtonsNoFirstLayerITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}}; - std::map hDcazResNoFirstLayerITS = { - {11, new TH1F("hDcazResElectronsNoFirstLayerITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}, - {211, new TH1F("hDcazResPionsNoFirstLayerITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}, - {321, new TH1F("hDcazResKaonsNoFirstLayerITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}, - {2212, new TH1F("hDcazResProtonsNoFirstLayerITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}}; - std::map hDcaxyReskAnyITS = { - {11, new TH1F("hDcaxyResElectronskAnyITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}, - {211, new TH1F("hDcaxyResPionskAnyITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}, - {321, new TH1F("hDcaxyResKaonskAnyITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}, - {2212, new TH1F("hDcaxyResProtonskAnyITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}}; - std::map hDcazReskAnyITS = { - {11, new TH1F("hDcazResElectronskAnyITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}, - {211, new TH1F("hDcazResPionskAnyITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}, - {321, new TH1F("hDcazResKaonskAnyITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}, - {2212, new TH1F("hDcazResProtonskAnyITS", "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}}; - - std::map hDcaxyVsPtAllLayersITS = { - {11, new TH2F("hDcaxyVsPtElectronsAllLayersITS", "ITS Electrons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {211, new TH2F("hDcaxyVsPtPionsAllLayersITS", "ITS Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {321, new TH2F("hDcaxyVsPtKaonsAllLayersITS", "ITS Kaons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {2212, new TH2F("hDcaxyVsPtProtonsAllLayersITS", "ITS Protons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}}; - std::map hDcazVsPtAllLayersITS = { - {11, new TH2F("hDcazVsPtElectronsAllLayersITS", "ITS Electrons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {211, new TH2F("hDcazVsPtPionsAllLayersITS", "ITS Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {321, new TH2F("hDcazVsPtKaonsAllLayersITS", "ITS Kaons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {2212, new TH2F("hDcazVsPtProtonsAllLayersITS", "ITS Protons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}}; - std::map hDcaxyVsPhiAllLayersITS = { - {11, new TH2F("hDcaxyVsPhiElectronsAllLayersITS", "ITS Electrons (>2 Gev);#varphi (rad);#sigma(DCA_{#it{xy}}) (#mum)", 100, 0.f, 2 * TMath::Pi(), 1000, -500, 500)}, - {211, new TH2F("hDcaxyVsPhiPionsAllLayersITS", "ITS Pions (>2 Gev);#varphi (rad);#sigma(DCA_{#it{xy}}) (#mum)", 100, 0.f, 2 * TMath::Pi(), 1000, -500, 500)}, - {321, new TH2F("hDcaxyVsPhiKaonsAllLayersITS", "ITS Kaons (>2 Gev);#varphi (rad);#sigma(DCA_{#it{xy}}) (#mum)", 100, 0.f, 2 * TMath::Pi(), 1000, -500, 500)}, - {2212, new TH2F("hDcaxyVsPhiProtonsAllLayersITS", "ITS Protons (>2 Gev);#varphi (rad);#sigma(DCA_{#it{xy}}) (#mum)", 100, 0.f, 2 * TMath::Pi(), 1000, -500, 500)}}; - std::map hDcazVsPhiAllLayersITS = { - {11, new TH2F("hDcazVsPhiElectronsAllLayersITS", "ITS Electrons (>2 Gev);#varphi (rad);#sigma(DCA_{#it{z}}) (#mum)", 100, 0.f, 2 * TMath::Pi(), 1000, -500, 500)}, - {211, new TH2F("hDcazVsPhiPionsAllLayersITS", "ITS Pions (>2 Gev);#varphi (rad);#sigma(DCA_{#it{z}}) (#mum)", 100, 0.f, 2 * TMath::Pi(), 1000, -500, 500)}, - {321, new TH2F("hDcazVsPhiKaonsAllLayersITS", "ITS Kaons (>2 Gev);#varphi (rad);#sigma(DCA_{#it{z}}) (#mum)", 100, 0.f, 2 * TMath::Pi(), 1000, -500, 500)}, - {2212, new TH2F("hDcazVsPhiProtonsAllLayersITS", "ITS Protons (>2 Gev);#varphi (rad);#sigma(DCA_{#it{z}}) (#mum)", 100, 0.f, 2 * TMath::Pi(), 1000, -500, 500)}}; - std::map hDcaxyVsPtNoFirstLayerITS = { - {11, new TH2F("hDcaxyVsPtElectronsNoFirstLayerITS", "ITS Electrons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {211, new TH2F("hDcaxyVsPtPionsNoFirstLayerITS", "ITS Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {321, new TH2F("hDcaxyVsPtKaonsNoFirstLayerITS", "ITS Kaons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {2212, new TH2F("hDcaxyVsPtProtonsNoFirstLayerITS", "ITS Protons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}}; - std::map hDcazVsPtNoFirstLayerITS = { - {11, new TH2F("hDcazVsPtElectronsNoFirstLayerITS", "ITS Electrons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {211, new TH2F("hDcazVsPtPionsNoFirstLayerITS", "ITS Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {321, new TH2F("hDcazVsPtKaonsNoFirstLayerITS", "ITS Kaons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {2212, new TH2F("hDcazVsPtProtonsNoFirstLayerITS", "ITS Protons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}}; - std::map hDcazVsPtkAnyITS = { - {11, new TH2F("hDcazVsPtElectronskAnyITS ", "ITS Electrons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {211, new TH2F("hDcazVsPtPionskAnyITS", "ITS Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {321, new TH2F("hDcazVsPtKaonskAnyITS", "ITS Kaons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {2212, new TH2F("hDcazVsPtProtonskAnyITS", "ITS Protons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}}; - std::map hDcaxyVsPtkAnyITS = { - {11, new TH2F("hDcaxyVsPtElectronskAnyITS", "ITS Electrons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {211, new TH2F("hDcaxyVsPtPionskAnyITS", "ITS Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {321, new TH2F("hDcaxyVsPtKaonskAnyITS", "ITS Kaons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {2212, new TH2F("hDcaxyVsPtProtonskAnyITS", "ITS Protons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}}; - std::map hDeltaPtVsPtAllLayersITS = { - {11, new TH2F("hDeltaPtVsPtElectronsAllLayersITS", "ITS Electrons;#it{p}_{T} (GeV/#it{c});#sigma(#Delta#it{p}_{T}/#it{p}_{T})", nPtBins - 1, ptLimits, 200, -0.2, 0.2)}, - {211, new TH2F("hDeltaPtVsPtPionsAllLayersITS", "ITS Pions;#it{p}_{T} (GeV/#it{c});#sigma(#Delta#it{p}_{T}/#it{p}_{T})", nPtBins - 1, ptLimits, 200, -0.2, 0.2)}, - {321, new TH2F("hDeltaPtVsPtKaonsAllLayersITS", "ITS Kaons;#it{p}_{T} (GeV/#it{c});#sigma(#Delta#it{p}_{T}/#it{p}_{T})", nPtBins - 1, ptLimits, 200, -0.2, 0.2)}, - {2212, new TH2F("hDeltaPtVsPtProtonsAllLayersITS", "ITS Protons;#it{p}_{T} (GeV/#it{c});#sigma(#Delta#it{p}_{T}/#it{p}_{T})", nPtBins - 1, ptLimits, 200, -0.2, 0.2)}}; - // ITS-TPC - std::map hDcaxyResAllLayersITSTPC = { - {11, new TH1F("hDcaxyResElectronsAllLayersITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}, - {211, new TH1F("hDcaxyResPionsAllLayersITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}, - {321, new TH1F("hDcaxyResKaonsAllLayersITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}, - {2212, new TH1F("hDcaxyResProtonsAllLayersITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}}; - std::map hDcazResAllLayersITSTPC = { - {11, new TH1F("hDcazResElectronsAllLayersITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}, - {211, new TH1F("hDcazResPionsAllLayersITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}, - {321, new TH1F("hDcazResKaonsAllLayersITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}, - {2212, new TH1F("hDcazResProtonsAllLayersITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}}; - std::map hPtResAllLayersITSTPC = { - {11, new TH1F("hPtResElectronsAllLayersITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(#Delta#it{p}_{T}/#it{p}_{T})", nPtBins - 1, ptLimits)}, - {211, new TH1F("hPtResPionsAllLayersITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(#Delta#it{p}_{T}/#it{p}_{T})", nPtBins - 1, ptLimits)}, - {321, new TH1F("hPtResKaonsAllLayersITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(#Delta#it{p}_{T}/#it{p}_{T})", nPtBins - 1, ptLimits)}, - {2212, new TH1F("hPtResProtonsAllLayersITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(#Delta#it{p}_{T}/#it{p}_{T})", nPtBins - 1, ptLimits)}}; - std::map hDcaxyResNoFirstLayerITSTPC = { - {11, new TH1F("hDcaxyResElectronsNoFirstLayerITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}, - {211, new TH1F("hDcaxyResPionsNoFirstLayerITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}, - {321, new TH1F("hDcaxyResKaonsNoFirstLayerITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}, - {2212, new TH1F("hDcaxyResProtonsNoFirstLayerITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}}; - std::map hDcazResNoFirstLayerITSTPC = { - {11, new TH1F("hDcazResElectronsNoFirstLayerITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}, - {211, new TH1F("hDcazResPionsNoFirstLayerITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}, - {321, new TH1F("hDcazResKaonsNoFirstLayerITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}, - {2212, new TH1F("hDcazResProtonsNoFirstLayerITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}}; - std::map hDcaxyReskAnyITSTPC = { - {11, new TH1F("hDcaxyResElectronskAnyITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}, - {211, new TH1F("hDcaxyResPionskAnyITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}, - {321, new TH1F("hDcaxyResKaonskAnyITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}, - {2212, new TH1F("hDcaxyResProtonskAnyITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits)}}; - std::map hDcazReskAnyITSTPC = { - {11, new TH1F("hDcazResElectronskAnyITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}, - {211, new TH1F("hDcazResPionskAnyITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}, - {321, new TH1F("hDcazResKaonskAnyITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}, - {2212, new TH1F("hDcazResProtonskAnyITSTPC", "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits)}}; - - std::map hDcaxyVsPtAllLayersITSTPC = { - {11, new TH2F("hDcaxyVsPtElectronsAllLayersITSTPC", "ITS-TPC Electrons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {211, new TH2F("hDcaxyVsPtPionsAllLayersITSTPC", "ITS-TPC Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {321, new TH2F("hDcaxyVsPtKaonsAllLayersITSTPC", "ITS-TPC Kaons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {2212, new TH2F("hDcaxyVsPtProtonsAllLayersITSTPC", "ITS-TPC Protons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}}; - std::map hDcazVsPtAllLayersITSTPC = { - {11, new TH2F("hDcazVsPtElectronsAllLayersITSTPC", "ITS-TPC Electrons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {211, new TH2F("hDcazVsPtPionsAllLayersITSTPC", "ITS-TPC Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {321, new TH2F("hDcazVsPtKaonsAllLayersITSTPC", "ITS-TPC Kaons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {2212, new TH2F("hDcazVsPtProtonsAllLayersITSTPC", "ITS-TPC Protons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}}; - std::map hDcaxyVsPhiAllLayersITSTPC = { - {11, new TH2F("hDcaxyVsPhiElectronsAllLayersITSTPC", "ITS-TPC Electrons (>2 Gev);#varphi (rad);#sigma(DCA_{#it{xy}}) (#mum)", 100, 0.f, 2 * TMath::Pi(), 1000, -500, 500)}, - {211, new TH2F("hDcaxyVsPhiPionsAllLayersITSTPC", "ITS-TPC Pions (>2 Gev);#varphi (rad);#sigma(DCA_{#it{xy}}) (#mum)", 100, 0.f, 2 * TMath::Pi(), 1000, -500, 500)}, - {321, new TH2F("hDcaxyVsPhiKaonsAllLayersITSTPC", "ITS-TPC Kaons (>2 Gev);#varphi (rad);#sigma(DCA_{#it{xy}}) (#mum)", 100, 0.f, 2 * TMath::Pi(), 1000, -500, 500)}, - {2212, new TH2F("hDcaxyVsPhiProtonsAllLayersITSTPC", "ITS-TPC Protons (>2 Gev);#varphi (rad);#sigma(DCA_{#it{xy}}) (#mum)", 100, 0.f, 2 * TMath::Pi(), 1000, -500, 500)}}; - std::map hDcazVsPhiAllLayersITSTPC = { - {11, new TH2F("hDcazVsPhiElectronsAllLayersITSTPC", "ITS-TPC Electrons (>2 Gev);#varphi (rad);#sigma(DCA_{#it{z}}) (#mum)", 100, 0.f, 2 * TMath::Pi(), 1000, -500, 500)}, - {211, new TH2F("hDcazVsPhiPionsAllLayersITSTPC", "ITS-TPC Pions (>2 Gev);#varphi (rad);#sigma(DCA_{#it{z}}) (#mum)", 100, 0.f, 2 * TMath::Pi(), 1000, -500, 500)}, - {321, new TH2F("hDcazVsPhiKaonsAllLayersITSTPC", "ITS-TPC Kaons (>2 Gev);#varphi (rad);#sigma(DCA_{#it{z}}) (#mum)", 100, 0.f, 2 * TMath::Pi(), 1000, -500, 500)}, - {2212, new TH2F("hDcazVsPhiProtonsAllLayersITSTPC", "ITS-TPC Protons (>2 Gev);#varphi (rad);#sigma(DCA_{#it{z}}) (#mum)", 100, 0.f, 2 * TMath::Pi(), 1000, -500, 500)}}; - std::map hDcaxyVsPtNoFirstLayerITSTPC = { - {11, new TH2F("hDcaxyVsPtElectronsNoFirstLayerITSTPC", "ITS-TPC Electrons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {211, new TH2F("hDcaxyVsPtPionsNoFirstLayerITSTPC", "ITS-TPC Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {321, new TH2F("hDcaxyVsPtKaonsNoFirstLayerITSTPC", "ITS-TPC Kaons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {2212, new TH2F("hDcaxyVsPtProtonsNoFirstLayerITSTPC", "ITS-TPC Protons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}}; - std::map hDcazVsPtNoFirstLayerITSTPC = { - {11, new TH2F("hDcazVsPtElectronsNoFirstLayerITSTPC", "ITS-TPC Electrons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {211, new TH2F("hDcazVsPtPionsNoFirstLayerITSTPC", "ITS-TPC Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {321, new TH2F("hDcazVsPtKaonsNoFirstLayerITSTPC", "ITS-TPC Kaons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {2212, new TH2F("hDcazVsPtProtonsNoFirstLayerITSTPC", "ITS-TPC Protons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}}; - std::map hDcazVsPtkAnyITSTPC = { - {11, new TH2F("hDcazVsPtElectronskAnyITSTPC", "ITS-TPC Electrons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {211, new TH2F("hDcazVsPtPionskAnyITSTPC", "ITS-TPC Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {321, new TH2F("hDcazVsPtKaonskAnyITSTPC", "ITS-TPC Kaons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {2212, new TH2F("hDcazVsPtProtonskAnyITSTPC", "ITS-TPC Protons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}}; - std::map hDcaxyVsPtkAnyITSTPC = { - {11, new TH2F("hDcaxyVsPtElectronskAnyITSTPC", "ITS-TPC Electrons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {211, new TH2F("hDcaxyVsPtPionskAnyITSTPC", "ITS-TPC Pions;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {321, new TH2F("hDcaxyVsPtKaonskAnyITSTPC", "ITS-TPC Kaons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}, - {2212, new TH2F("hDcaxyVsPtProtonskAnyITSTPC", "ITS-TPC Protons;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)", nPtBins - 1, ptLimits, 1000, -500, 500)}}; - std::map hDeltaPtVsPtAllLayersITSTPC = { - {11, new TH2F("hDeltaPtVsPtElectronsAllLayersITSTPC", "ITS-TPC Electrons;#it{p}_{T} (GeV/#it{c});#sigma(#Delta#it{p}_{T}/#it{p}_{T})", nPtBins - 1, ptLimits, 200, -0.2, 0.2)}, - {211, new TH2F("hDeltaPtVsPtPionsAllLayersITSTPC", "ITS-TPC Pions;#it{p}_{T} (GeV/#it{c});#sigma(#Delta#it{p}_{T}/#it{p}_{T})", nPtBins - 1, ptLimits, 200, -0.2, 0.2)}, - {321, new TH2F("hDeltaPtVsPtKaonsAllLayersITSTPC", "ITS-TPC Kaons;#it{p}_{T} (GeV/#it{c});#sigma(#Delta#it{p}_{T}/#it{p}_{T})", nPtBins - 1, ptLimits, 200, -0.2, 0.2)}, - {2212, new TH2F("hDeltaPtVsPtProtonsAllLayersITSTPC", "ITS-TPC Protons;#it{p}_{T} (GeV/#it{c});#sigma(#Delta#it{p}_{T}/#it{p}_{T})", nPtBins - 1, ptLimits, 200, -0.2, 0.2)}}; - - o2::dataformats::VertexBase collision; - o2::dataformats::DCA impactParameter; - - const auto origWD{fs::current_path()}; - const auto tfDirs = find_dirs(fs::current_path(), tf_matcher, tf_sorter); - for (const auto& tfDir : tfDirs) { - LOGP(info, "Analysing {:?}", tfDir.c_str()); - fs::current_path(tfDir); - - // MC Information - o2::steer::MCKinematicsReader mcReader; - if (!mcReader.initFromDigitContext(collisioncontextFileName)) { - LOGP(error, "Cannot init MC reader in {:?}", tfDir.c_str()); - continue; - } - - // Magnetic field and Propagator - float bz{-999}; - static bool initOnce{false}; - if (!initOnce) { - initOnce = true; - o2::base::Propagator::initFieldFromGRP(magFileName); - bz = o2::base::Propagator::Instance()->getNominalBz(); - } - - LOGP(info, "Loading ITS Tracks"); - auto fITSTracks = TFile::Open(itsTracksFileName.c_str(), "READ"); - auto tITSTracks = fITSTracks->Get("o2sim"); - std::vector* itsTracks{nullptr}; - tITSTracks->SetBranchAddress("ITSTrack", &itsTracks); - std::vector* itsTrkLab{nullptr}; - tITSTracks->SetBranchAddress("ITSTrackMCTruth", &itsTrkLab); - - for (Long64_t iEntry{0}; tITSTracks->LoadTree(iEntry) >= 0; ++iEntry) { - tITSTracks->GetEntry(iEntry); - for (size_t iTrk{0}; iTrk < itsTracks->size(); ++iTrk) { - auto trk = itsTracks->at(iTrk); - const auto& lbl = itsTrkLab->at(iTrk); - if (!lbl.isValid()) { - continue; - } - - const auto& mcEvent = mcReader.getMCEventHeader(lbl.getSourceID(), lbl.getEventID()); - const auto& mcTrack = mcReader.getTrack(lbl); - if (!mcTrack->isPrimary() || !(std::abs(mcTrack->GetEta()) < rapMax)) { - continue; - } - auto pdg = std::abs(mcTrack->GetPdgCode()); - if (pdg != 11 && pdg != 211 && pdg != 321 && pdg != 2212) { - continue; - } - - collision.setXYZ(mcEvent.GetX(), mcEvent.GetY(), mcEvent.GetZ()); - if (!o2::base::Propagator::Instance()->propagateToDCA(collision, trk, bz, mMaxStep, mMatCorr, &impactParameter)) { - continue; - } - - auto ptReco = trk.getPt(); - auto ptGen = mcTrack->GetPt(); - auto deltaPt = (1. / ptReco - 1. / ptGen) / (1. / ptGen); - auto dcaXY = impactParameter.getY() * 10000; - auto dcaZ = impactParameter.getZ() * 10000; - auto phiReco = trk.getPhi(); - - if (trk.getNumberOfClusters() == 7) { - hDcaxyVsPtAllLayersITS[pdg]->Fill(ptGen, dcaXY); - hDcazVsPtAllLayersITS[pdg]->Fill(ptGen, dcaZ); - hDeltaPtVsPtAllLayersITS[pdg]->Fill(ptGen, deltaPt); - if (ptGen > 2.) { - hDcaxyVsPhiAllLayersITS[pdg]->Fill(phiReco, dcaXY); - hDcazVsPhiAllLayersITS[pdg]->Fill(phiReco, dcaZ); - } - } else if (!trk.hasHitOnLayer(0)) { - hDcaxyVsPtNoFirstLayerITS[pdg]->Fill(ptGen, dcaXY); - hDcazVsPtNoFirstLayerITS[pdg]->Fill(ptGen, dcaZ); - } else { - hDcaxyVsPtkAnyITS[pdg]->Fill(ptGen, dcaXY); - hDcazVsPtkAnyITS[pdg]->Fill(ptGen, dcaZ); - } - } - } - - LOGP(info, "Loading ITS-TPC Tracks"); - auto fITSTPCTracks = TFile::Open(itstpcTracksFileName.c_str(), "READ"); - auto tITSTPCTracks = fITSTPCTracks->Get("matchTPCITS"); - std::vector* itstpcTracks{nullptr}; - tITSTPCTracks->SetBranchAddress("TPCITS", &itstpcTracks); - std::vector* itstpcTrkLab{nullptr}; - tITSTPCTracks->SetBranchAddress("MatchMCTruth", &itstpcTrkLab); - // TPC Tracks - auto fTPCTracks = TFile::Open(tpcTracksFileName.c_str(), "READ"); - auto tTPCTracks = fTPCTracks->Get("tpcrec"); - std::vector* tpcTracks{nullptr}; - tTPCTracks->SetBranchAddress("TPCTracks", &tpcTracks); - std::vector* tpcTrkLab{nullptr}; - tTPCTracks->SetBranchAddress("TPCTracksMCTruth", &tpcTrkLab); - for (Long64_t iEntry{0}; tITSTPCTracks->LoadTree(iEntry) >= 0; ++iEntry) { - tITSTPCTracks->GetEntry(iEntry); - tITSTracks->GetEntry(iEntry); - tTPCTracks->GetEntry(iEntry); - for (size_t iTrk{0}; iTrk < itstpcTracks->size(); ++iTrk) { - auto trk = itstpcTracks->at(iTrk); - const auto& lbl = itstpcTrkLab->at(iTrk); - - const auto& trkITS = itsTracks->at(trk.getRefITS().getIndex()); - const auto& trkITSLbl = itsTrkLab->at(trk.getRefITS().getIndex()); - const auto& trkTPC = tpcTracks->at(trk.getRefTPC().getIndex()); - const auto& trkTPCLbl = tpcTrkLab->at(trk.getRefTPC().getIndex()); - if (!lbl.isValid() || trkITSLbl != trkTPCLbl) { - continue; - } - - const auto& mcEvent = mcReader.getMCEventHeader(lbl.getSourceID(), lbl.getEventID()); - const auto& mcTrack = mcReader.getTrack(lbl); - if (!mcTrack->isPrimary() || !(std::abs(mcTrack->GetEta()) < rapMax)) { - continue; - } - - auto pdg = std::abs(mcTrack->GetPdgCode()); - if (pdg != 11 && pdg != 211 && pdg != 321 && pdg != 2212) { - continue; - } - - collision.setXYZ(mcEvent.GetX(), mcEvent.GetY(), mcEvent.GetZ()); - if (!o2::base::Propagator::Instance()->propagateToDCA(collision, trk, bz, mMaxStep, mMatCorr, &impactParameter)) { - continue; - } - - auto ptReco = trk.getPt(); - auto ptGen = mcTrack->GetPt(); - auto deltaPt = (1. / ptReco - 1. / ptGen) / (1. / ptGen); - auto dcaXY = impactParameter.getY() * 10000; - auto dcaZ = impactParameter.getZ() * 10000; - auto phiReco = trk.getPhi(); - - if (trkITS.getNumberOfClusters() == 7) { - hDcaxyVsPtAllLayersITSTPC[pdg]->Fill(ptGen, dcaXY); - hDcazVsPtAllLayersITSTPC[pdg]->Fill(ptGen, dcaZ); - hDeltaPtVsPtAllLayersITSTPC[pdg]->Fill(ptGen, deltaPt); - if (ptGen > 2.) { - hDcaxyVsPhiAllLayersITSTPC[pdg]->Fill(phiReco, dcaXY); - hDcazVsPhiAllLayersITSTPC[pdg]->Fill(phiReco, dcaZ); - } - } else if (!trkITS.hasHitOnLayer(0) && !trkITS.hasHitOnLayer(1) && !trkITS.hasHitOnLayer(2)) { - hDcaxyVsPtNoFirstLayerITSTPC[pdg]->Fill(ptGen, dcaXY); - hDcazVsPtNoFirstLayerITSTPC[pdg]->Fill(ptGen, dcaZ); - } else { - hDcaxyVsPtkAnyITSTPC[pdg]->Fill(ptGen, dcaXY); - hDcazVsPtkAnyITSTPC[pdg]->Fill(ptGen, dcaZ); - } - } - } - - delete itsTracks; - delete itsTrkLab; - delete tpcTracks; - delete tpcTrkLab; - delete itstpcTracks; - delete itstpcTrkLab; - delete tITSTracks; - delete tTPCTracks; - delete tITSTPCTracks; - delete fITSTracks; - delete fTPCTracks; - delete fITSTPCTracks; - - gSystem->GetProcInfo(&procInfo); - LOGF(info, "MemVirtual (%ld), MemResident (%ld)", procInfo.fMemVirtual, procInfo.fMemResident); - LOGP(info, "Done with {:?}", tfDir.c_str()); - if (procInfo.fMemResident > 200'000'000) { - LOGP(error, "Exceeding 200GBs stopping!"); - break; - } - } - LOGP(info, "Restoring original CWD to {:?}", origWD.c_str()); - fs::current_path(origWD); // restore original wd - - LOGP(info, "Projecting Plots"); - TH1* hProj; - const char* fitOpt{"QWMER"}; - /* const char* fitOpt{"Q"}; */ - std::map lProjITS = { - {11, new TList()}, - {211, new TList()}, - {321, new TList()}, - {2212, new TList()}, - }; - std::map lProjITSTPC = { - {11, new TList()}, - {211, new TList()}, - {321, new TList()}, - {2212, new TList()}, - }; - for (const auto& pdgCode : pdgCodes) { - for (auto iPt{0}; iPt < nPtBins; ++iPt) { - // ITS - auto ptMin = hDcaxyVsPtAllLayersITS[pdgCode]->GetXaxis()->GetBinLowEdge(iPt + 1); - float minFit = (ptMin < 1.) ? -200. : -50.; - float maxFit = (ptMin < 1.) ? 200. : 50.; - - hProj = hDeltaPtVsPtAllLayersITS[pdgCode]->ProjectionY(Form("hProjDeltaPtAll%d__%dITS", pdgCode, iPt), iPt + 1, iPt + 1); - hProj->Fit("fGaus", fitOpt); - lProjITS[pdgCode]->Add(hProj); - hPtResAllLayersITS[pdgCode]->SetBinContent(iPt + 1, fGaus->GetParameter(2)); - hPtResAllLayersITS[pdgCode]->SetBinError(iPt + 1, fGaus->GetParError(2)); - - hProj = hDcaxyVsPtAllLayersITS[pdgCode]->ProjectionY(Form("hProjDcaxyAll%d__%dITS", pdgCode, iPt), iPt + 1, iPt + 1); - hProj->Fit("fGaus", fitOpt, "", minFit, maxFit); - lProjITS[pdgCode]->Add(hProj); - hDcaxyResAllLayersITS[pdgCode]->SetBinContent(iPt + 1, fGaus->GetParameter(2)); - hDcaxyResAllLayersITS[pdgCode]->SetBinError(iPt + 1, fGaus->GetParError(2)); - - hProj = hDcazVsPtAllLayersITS[pdgCode]->ProjectionY(Form("hProjDcazAll%d__%dITS", pdgCode, iPt), iPt + 1, iPt + 1); - hProj->Fit("fGaus", fitOpt, "", minFit, maxFit); - lProjITS[pdgCode]->Add(hProj); - hDcazResAllLayersITS[pdgCode]->SetBinContent(iPt + 1, fGaus->GetParameter(2)); - hDcazResAllLayersITS[pdgCode]->SetBinError(iPt + 1, fGaus->GetParError(2)); - - hProj = hDcaxyVsPtNoFirstLayerITS[pdgCode]->ProjectionY(Form("hProjDcaxyNoFirst%d__%dITS", pdgCode, iPt), iPt + 1, iPt + 1); - hProj->Fit("fGaus", fitOpt, "", minFit, maxFit); - lProjITS[pdgCode]->Add(hProj); - hDcaxyResNoFirstLayerITS[pdgCode]->SetBinContent(iPt + 1, fGaus->GetParameter(2)); - hDcaxyResNoFirstLayerITS[pdgCode]->SetBinError(iPt + 1, fGaus->GetParError(2)); - - hProj = hDcazVsPtNoFirstLayerITS[pdgCode]->ProjectionY(Form("hProjDcazNoFirst%d__%dITS", pdgCode, iPt), iPt + 1, iPt + 1); - hProj->Fit("fGaus", fitOpt, "", minFit, maxFit); - lProjITS[pdgCode]->Add(hProj); - hDcazResNoFirstLayerITS[pdgCode]->SetBinContent(iPt + 1, fGaus->GetParameter(2)); - hDcazResNoFirstLayerITS[pdgCode]->SetBinError(iPt + 1, fGaus->GetParError(2)); - - hProj = hDcaxyVsPtkAnyITS[pdgCode]->ProjectionY(Form("hProjDcaxyAny%d__%dITS", pdgCode, iPt), iPt + 1, iPt + 1); - hProj->Fit("fGaus", fitOpt, "", minFit, maxFit); - lProjITS[pdgCode]->Add(hProj); - hDcaxyReskAnyITS[pdgCode]->SetBinContent(iPt + 1, fGaus->GetParameter(2)); - hDcaxyReskAnyITS[pdgCode]->SetBinError(iPt + 1, fGaus->GetParError(2)); - - hProj = hDcazVsPtkAnyITS[pdgCode]->ProjectionY(Form("hProjDcazAny%d__%dITS", pdgCode, iPt), iPt + 1, iPt + 1); - hProj->Fit("fGaus", fitOpt, "", minFit, maxFit); - lProjITS[pdgCode]->Add(hProj); - hDcazReskAnyITS[pdgCode]->SetBinContent(iPt + 1, fGaus->GetParameter(2)); - hDcazReskAnyITS[pdgCode]->SetBinError(iPt + 1, fGaus->GetParError(2)); - - // ITS-TPC - hProj = hDeltaPtVsPtAllLayersITSTPC[pdgCode]->ProjectionY(Form("hProjDeltaPtAll%d__%dITSTPC", pdgCode, iPt), iPt + 1, iPt + 1); - hProj->Fit("fGaus", fitOpt); - lProjITSTPC[pdgCode]->Add(hProj); - hPtResAllLayersITSTPC[pdgCode]->SetBinContent(iPt + 1, fGaus->GetParameter(2)); - hPtResAllLayersITSTPC[pdgCode]->SetBinError(iPt + 1, fGaus->GetParError(2)); - - ptMin = hDcaxyVsPtAllLayersITSTPC[pdgCode]->GetXaxis()->GetBinLowEdge(iPt + 1); - minFit = (ptMin < 1.) ? -200. : -50.; - maxFit = (ptMin < 1.) ? 200. : 50.; - - hProj = hDcaxyVsPtAllLayersITSTPC[pdgCode]->ProjectionY(Form("hProjDcaxyAll%d__%dITSTPC", pdgCode, iPt), iPt + 1, iPt + 1); - hProj->Fit("fGaus", fitOpt, "", minFit, maxFit); - lProjITSTPC[pdgCode]->Add(hProj); - hDcaxyResAllLayersITSTPC[pdgCode]->SetBinContent(iPt + 1, fGaus->GetParameter(2)); - hDcaxyResAllLayersITSTPC[pdgCode]->SetBinError(iPt + 1, fGaus->GetParError(2)); - - hProj = hDcazVsPtAllLayersITSTPC[pdgCode]->ProjectionY(Form("hProjDcazAll%d__%dITSTPC", pdgCode, iPt), iPt + 1, iPt + 1); - hProj->Fit("fGaus", fitOpt, "", minFit, maxFit); - lProjITSTPC[pdgCode]->Add(hProj); - hDcazResAllLayersITSTPC[pdgCode]->SetBinContent(iPt + 1, fGaus->GetParameter(2)); - hDcazResAllLayersITSTPC[pdgCode]->SetBinError(iPt + 1, fGaus->GetParError(2)); - - hProj = hDcaxyVsPtNoFirstLayerITSTPC[pdgCode]->ProjectionY(Form("hProjDcaxyNoFirst%d__%dITSTPC", pdgCode, iPt), iPt + 1, iPt + 1); - hProj->Fit("fGaus", fitOpt, "", minFit, maxFit); - lProjITSTPC[pdgCode]->Add(hProj); - hDcaxyResNoFirstLayerITSTPC[pdgCode]->SetBinContent(iPt + 1, fGaus->GetParameter(2)); - hDcaxyResNoFirstLayerITSTPC[pdgCode]->SetBinError(iPt + 1, fGaus->GetParError(2)); - - hProj = hDcazVsPtNoFirstLayerITSTPC[pdgCode]->ProjectionY(Form("hProjDcazNoFirst%d__%dITSTPC", pdgCode, iPt), iPt + 1, iPt + 1); - hProj->Fit("fGaus", fitOpt, "", minFit, maxFit); - lProjITSTPC[pdgCode]->Add(hProj); - hDcazResNoFirstLayerITSTPC[pdgCode]->SetBinContent(iPt + 1, fGaus->GetParameter(2)); - hDcazResNoFirstLayerITSTPC[pdgCode]->SetBinError(iPt + 1, fGaus->GetParError(2)); - - hProj = hDcaxyVsPtkAnyITSTPC[pdgCode]->ProjectionY(Form("hProjDcaxAnty%d__%dITSTPC", pdgCode, iPt), iPt + 1, iPt + 1); - hProj->Fit("fGaus", fitOpt, "", minFit, maxFit); - lProjITSTPC[pdgCode]->Add(hProj); - hDcaxyReskAnyITSTPC[pdgCode]->SetBinContent(iPt + 1, fGaus->GetParameter(2)); - hDcaxyReskAnyITSTPC[pdgCode]->SetBinError(iPt + 1, fGaus->GetParError(2)); - - hProj = hDcazVsPtkAnyITSTPC[pdgCode]->ProjectionY(Form("hProjDcazAny%d__%dITSTPC", pdgCode, iPt), iPt + 1, iPt + 1); - hProj->Fit("fGaus", fitOpt, "", minFit, maxFit); - lProjITSTPC[pdgCode]->Add(hProj); - hDcazReskAnyITSTPC[pdgCode]->SetBinContent(iPt + 1, fGaus->GetParameter(2)); - hDcazReskAnyITSTPC[pdgCode]->SetBinError(iPt + 1, fGaus->GetParError(2)); - } - } - - // Style - LOGP(info, "Styling Plots"); - for (const auto& pdgCode : pdgCodes) { - // ITS - hPtResAllLayersITS[pdgCode]->SetLineWidth(2); - hPtResAllLayersITS[pdgCode]->SetMarkerColor(colors[pdgCode]); - hPtResAllLayersITS[pdgCode]->SetLineColor(colors[pdgCode]); - hPtResAllLayersITS[pdgCode]->SetMarkerStyle(kFullCircle); - - hDcaxyResAllLayersITS[pdgCode]->SetLineWidth(2); - hDcaxyResAllLayersITS[pdgCode]->SetMarkerColor(colors[pdgCode]); - hDcaxyResAllLayersITS[pdgCode]->SetLineColor(colors[pdgCode]); - hDcaxyResAllLayersITS[pdgCode]->SetMarkerStyle(kFullCircle); - - hDcazResAllLayersITS[pdgCode]->SetLineWidth(2); - hDcazResAllLayersITS[pdgCode]->SetMarkerColor(colors[pdgCode]); - hDcazResAllLayersITS[pdgCode]->SetLineColor(colors[pdgCode]); - hDcazResAllLayersITS[pdgCode]->SetMarkerStyle(kFullCircle); - - hDcaxyResNoFirstLayerITS[pdgCode]->SetLineWidth(2); - hDcaxyResNoFirstLayerITS[pdgCode]->SetMarkerColor(colors[pdgCode]); - hDcaxyResNoFirstLayerITS[pdgCode]->SetLineColor(colors[pdgCode]); - hDcaxyResNoFirstLayerITS[pdgCode]->SetMarkerStyle(kFullCircle); - - hDcazResNoFirstLayerITS[pdgCode]->SetLineWidth(2); - hDcazResNoFirstLayerITS[pdgCode]->SetMarkerColor(colors[pdgCode]); - hDcazResNoFirstLayerITS[pdgCode]->SetLineColor(colors[pdgCode]); - hDcazResNoFirstLayerITS[pdgCode]->SetMarkerStyle(kFullCircle); - - hDcaxyReskAnyITS[pdgCode]->SetLineWidth(2); - hDcaxyReskAnyITS[pdgCode]->SetMarkerColor(colors[pdgCode]); - hDcaxyReskAnyITS[pdgCode]->SetLineColor(colors[pdgCode]); - hDcaxyReskAnyITS[pdgCode]->SetMarkerStyle(kFullCircle); - - hDcazReskAnyITS[pdgCode]->SetLineWidth(2); - hDcazReskAnyITS[pdgCode]->SetMarkerColor(colors[pdgCode]); - hDcazReskAnyITS[pdgCode]->SetLineColor(colors[pdgCode]); - hDcazReskAnyITS[pdgCode]->SetMarkerStyle(kFullCircle); - - // ITS-TPC - hPtResAllLayersITSTPC[pdgCode]->SetLineWidth(2); - hPtResAllLayersITSTPC[pdgCode]->SetMarkerColor(colors[pdgCode]); - hPtResAllLayersITSTPC[pdgCode]->SetLineColor(colors[pdgCode]); - hPtResAllLayersITSTPC[pdgCode]->SetMarkerStyle(kOpenCircle); - - hDcaxyResAllLayersITSTPC[pdgCode]->SetLineWidth(2); - hDcaxyResAllLayersITSTPC[pdgCode]->SetMarkerColor(colors[pdgCode]); - hDcaxyResAllLayersITSTPC[pdgCode]->SetLineColor(colors[pdgCode]); - hDcaxyResAllLayersITSTPC[pdgCode]->SetMarkerStyle(kOpenCircle); - - hDcazResAllLayersITSTPC[pdgCode]->SetLineWidth(2); - hDcazResAllLayersITSTPC[pdgCode]->SetMarkerColor(colors[pdgCode]); - hDcazResAllLayersITSTPC[pdgCode]->SetLineColor(colors[pdgCode]); - hDcazResAllLayersITSTPC[pdgCode]->SetMarkerStyle(kOpenCircle); - - hDcaxyResNoFirstLayerITSTPC[pdgCode]->SetLineWidth(2); - hDcaxyResNoFirstLayerITSTPC[pdgCode]->SetMarkerColor(colors[pdgCode]); - hDcaxyResNoFirstLayerITSTPC[pdgCode]->SetLineColor(colors[pdgCode]); - hDcaxyResNoFirstLayerITSTPC[pdgCode]->SetMarkerStyle(kOpenCircle); - - hDcazResNoFirstLayerITSTPC[pdgCode]->SetLineWidth(2); - hDcazResNoFirstLayerITSTPC[pdgCode]->SetMarkerColor(colors[pdgCode]); - hDcazResNoFirstLayerITSTPC[pdgCode]->SetLineColor(colors[pdgCode]); - hDcazResNoFirstLayerITSTPC[pdgCode]->SetMarkerStyle(kOpenCircle); - - hDcaxyReskAnyITSTPC[pdgCode]->SetLineWidth(2); - hDcaxyReskAnyITSTPC[pdgCode]->SetMarkerColor(colors[pdgCode]); - hDcaxyReskAnyITSTPC[pdgCode]->SetLineColor(colors[pdgCode]); - hDcaxyReskAnyITSTPC[pdgCode]->SetMarkerStyle(kOpenCircle); - - hDcazReskAnyITSTPC[pdgCode]->SetLineWidth(2); - hDcazReskAnyITSTPC[pdgCode]->SetMarkerColor(colors[pdgCode]); - hDcazReskAnyITSTPC[pdgCode]->SetLineColor(colors[pdgCode]); - hDcazReskAnyITSTPC[pdgCode]->SetMarkerStyle(kOpenCircle); - } - - /// Output - LOGP(info, "Writing final output"); - // ITS - auto canvPtDeltaITS = new TCanvas("canvPtDeltaITS", "", 1500, 500); - canvPtDeltaITS->Divide(nSpecies, 1); - canvPtDeltaITS->cd(1)->SetLogz(); - hDeltaPtVsPtAllLayersITS[11]->Draw("colz"); - canvPtDeltaITS->cd(2)->SetLogz(); - hDeltaPtVsPtAllLayersITS[211]->Draw("colz"); - canvPtDeltaITS->cd(3)->SetLogz(); - hDeltaPtVsPtAllLayersITS[321]->Draw("colz"); - canvPtDeltaITS->cd(4)->SetLogz(); - hDeltaPtVsPtAllLayersITS[2212]->Draw("colz"); - - auto canvDcaVsPtITS = new TCanvas("canvDcaVsPtITS", "", 1500, 1000); - canvDcaVsPtITS->Divide(nSpecies, 2); - canvDcaVsPtITS->cd(1)->SetLogz(); - hDcaxyVsPtAllLayersITS[11]->Draw("colz"); - canvDcaVsPtITS->cd(2)->SetLogz(); - hDcaxyVsPtAllLayersITS[211]->Draw("colz"); - canvDcaVsPtITS->cd(3)->SetLogz(); - hDcaxyVsPtAllLayersITS[321]->Draw("colz"); - canvDcaVsPtITS->cd(4)->SetLogz(); - hDcaxyVsPtAllLayersITS[2212]->Draw("colz"); - canvDcaVsPtITS->cd(5)->SetLogz(); - hDcazVsPtAllLayersITS[11]->Draw("colz"); - canvDcaVsPtITS->cd(6)->SetLogz(); - hDcazVsPtAllLayersITS[211]->Draw("colz"); - canvDcaVsPtITS->cd(7)->SetLogz(); - hDcazVsPtAllLayersITS[321]->Draw("colz"); - canvDcaVsPtITS->cd(8)->SetLogz(); - hDcazVsPtAllLayersITS[2212]->Draw("colz"); - - auto canvPtResITS = new TCanvas("canvPtResITS", "", 500, 500); - canvPtResITS->DrawFrame(ptLimits[0], 0., ptLimits[nPtBins - 1], 0.2, "ITS;#it{p}_{T} (GeV/#it{c});#sigma(#Delta#it{p}_{T}/#it{p}_{T})"); - for (const auto& pdgCode : pdgCodes) { - hPtResAllLayersITS[pdgCode]->Draw("same"); - } - - auto canvDcaxyResITS = new TCanvas("canvDcaxyResITS", "", 500, 500); - canvDcaxyResITS->DrawFrame(ptLimits[0], 1., ptLimits[nPtBins - 1], 1.e3, "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)"); - canvDcaxyResITS->SetLogx(); - canvDcaxyResITS->SetLogy(); - canvDcaxyResITS->SetGrid(); - for (const auto& pdgCode : pdgCodes) { - hDcaxyResAllLayersITS[pdgCode]->Draw("same"); - } - - auto canvDcazResITS = new TCanvas("canvDcazResITS", "", 500, 500); - canvDcazResITS->DrawFrame(ptLimits[0], 1., ptLimits[nPtBins - 1], 1.e3, "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)"); - canvDcazResITS->SetLogx(); - canvDcazResITS->SetLogy(); - canvDcazResITS->SetGrid(); - for (const auto& pdgCode : pdgCodes) { - hDcazResAllLayersITS[pdgCode]->Draw("same"); - } - - auto canvDcaxyResNoFirstLayerITS = new TCanvas("canvDcaxyResNoFirstLayerITS", "", 500, 500); - canvDcaxyResNoFirstLayerITS->DrawFrame(ptLimits[0], 1., ptLimits[nPtBins - 1], 1.e3, "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)"); - canvDcaxyResNoFirstLayerITS->SetLogx(); - canvDcaxyResNoFirstLayerITS->SetLogy(); - canvDcaxyResNoFirstLayerITS->SetGrid(); - for (const auto& pdgCode : pdgCodes) { - hDcaxyResNoFirstLayerITS[pdgCode]->Draw("same"); - } - - auto canvDcazResNoFirstLayerITS = new TCanvas("canvDcazResNoFirstLayerITS", "", 500, 500); - canvDcazResNoFirstLayerITS->DrawFrame(ptLimits[0], 1., ptLimits[nPtBins - 1], 1.e3, "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)"); - canvDcazResNoFirstLayerITS->SetLogx(); - canvDcazResNoFirstLayerITS->SetLogy(); - canvDcazResNoFirstLayerITS->SetGrid(); - for (const auto& pdgCode : pdgCodes) { - hDcazResNoFirstLayerITS[pdgCode]->Draw("same"); - } - - auto canvDcaxyReskAnyITS = new TCanvas("canvDcaxyReskAnyITS", "", 500, 500); - canvDcaxyReskAnyITS->DrawFrame(ptLimits[0], 1., ptLimits[nPtBins - 1], 1.e3, "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)"); - canvDcaxyReskAnyITS->SetLogx(); - canvDcaxyReskAnyITS->SetLogy(); - canvDcaxyReskAnyITS->SetGrid(); - for (const auto& pdgCode : pdgCodes) { - hDcaxyReskAnyITS[pdgCode]->Draw("same"); - } - - auto canvDcazReskAnyITS = new TCanvas("canvDcazReskAnyITS", "", 500, 500); - canvDcazReskAnyITS->DrawFrame(ptLimits[0], 1., ptLimits[nPtBins - 1], 1.e3, "ITS;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)"); - canvDcazReskAnyITS->SetLogx(); - canvDcazReskAnyITS->SetLogy(); - canvDcazReskAnyITS->SetGrid(); - for (const auto& pdgCode : pdgCodes) { - hDcazReskAnyITS[pdgCode]->Draw("same"); - } - - // ITS-TPC - auto canvPtDeltaITSTPC = new TCanvas("canvPtDeltaITSTPC", "", 1500, 500); - canvPtDeltaITSTPC->Divide(nSpecies, 1); - canvPtDeltaITSTPC->cd(1)->SetLogz(); - hDeltaPtVsPtAllLayersITSTPC[11]->Draw("colz"); - canvPtDeltaITSTPC->cd(2)->SetLogz(); - hDeltaPtVsPtAllLayersITSTPC[211]->Draw("colz"); - canvPtDeltaITSTPC->cd(3)->SetLogz(); - hDeltaPtVsPtAllLayersITSTPC[321]->Draw("colz"); - canvPtDeltaITSTPC->cd(4)->SetLogz(); - hDeltaPtVsPtAllLayersITSTPC[2212]->Draw("colz"); - - auto canvDcaVsPtITSTPC = new TCanvas("canvDcaVsPtITSTPC", "", 1500, 1000); - canvDcaVsPtITSTPC->Divide(nSpecies, 2); - canvDcaVsPtITSTPC->cd(1)->SetLogz(); - hDcaxyVsPtAllLayersITSTPC[11]->Draw("colz"); - canvDcaVsPtITSTPC->cd(2)->SetLogz(); - hDcaxyVsPtAllLayersITSTPC[211]->Draw("colz"); - canvDcaVsPtITSTPC->cd(3)->SetLogz(); - hDcaxyVsPtAllLayersITSTPC[321]->Draw("colz"); - canvDcaVsPtITSTPC->cd(4)->SetLogz(); - hDcaxyVsPtAllLayersITSTPC[2212]->Draw("colz"); - canvDcaVsPtITSTPC->cd(5)->SetLogz(); - hDcazVsPtAllLayersITSTPC[11]->Draw("colz"); - canvDcaVsPtITSTPC->cd(6)->SetLogz(); - hDcazVsPtAllLayersITSTPC[211]->Draw("colz"); - canvDcaVsPtITSTPC->cd(7)->SetLogz(); - hDcazVsPtAllLayersITSTPC[321]->Draw("colz"); - canvDcaVsPtITSTPC->cd(8)->SetLogz(); - hDcazVsPtAllLayersITSTPC[2212]->Draw("colz"); - - auto canvPtResITSTPC = new TCanvas("canvPtResITSTPC", "", 500, 500); - canvPtResITSTPC->DrawFrame(ptLimits[0], 0., ptLimits[nPtBins - 1], 0.2, "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(#Delta#it{p}_{T}/#it{p}_{T})"); - for (const auto& pdgCode : pdgCodes) { - hPtResAllLayersITSTPC[pdgCode]->Draw("same"); - } - - auto canvDcaxyResITSTPC = new TCanvas("canvDcaxyResITSTPC", "", 500, 500); - canvDcaxyResITSTPC->DrawFrame(ptLimits[0], 1., ptLimits[nPtBins - 1], 1.e3, "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)"); - canvDcaxyResITSTPC->SetLogx(); - canvDcaxyResITSTPC->SetLogy(); - canvDcaxyResITSTPC->SetGrid(); - for (const auto& pdgCode : pdgCodes) { - hDcaxyResAllLayersITSTPC[pdgCode]->Draw("same"); - } - - auto canvDcazResITSTPC = new TCanvas("canvDcazResITSTPC", "", 500, 500); - canvDcazResITSTPC->DrawFrame(ptLimits[0], 1., ptLimits[nPtBins - 1], 1.e3, "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)"); - canvDcazResITSTPC->SetLogx(); - canvDcazResITSTPC->SetLogy(); - canvDcazResITSTPC->SetGrid(); - for (const auto& pdgCode : pdgCodes) { - hDcazResAllLayersITSTPC[pdgCode]->Draw("same"); - } - - auto canvDcaxyResNoFirstLayerITSTPC = new TCanvas("canvDcaxyResNoFirstLayerITSTPC", "", 500, 500); - canvDcaxyResNoFirstLayerITSTPC->DrawFrame(ptLimits[0], 1., ptLimits[nPtBins - 1], 1.e3, "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)"); - canvDcaxyResNoFirstLayerITSTPC->SetLogx(); - canvDcaxyResNoFirstLayerITSTPC->SetLogy(); - canvDcaxyResNoFirstLayerITSTPC->SetGrid(); - for (const auto& pdgCode : pdgCodes) { - hDcaxyResNoFirstLayerITSTPC[pdgCode]->Draw("same"); - } - - auto canvDcazResNoFirstLayerITSTPC = new TCanvas("canvDcazResNoFirstLayerITSTPC", "", 500, 500); - canvDcazResNoFirstLayerITSTPC->DrawFrame(ptLimits[0], 1., ptLimits[nPtBins - 1], 1.e3, "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)"); - canvDcazResNoFirstLayerITSTPC->SetLogx(); - canvDcazResNoFirstLayerITSTPC->SetLogy(); - canvDcazResNoFirstLayerITSTPC->SetGrid(); - for (const auto& pdgCode : pdgCodes) { - hDcazResNoFirstLayerITSTPC[pdgCode]->Draw("same"); - } - - auto canvDcaxyReskAnyITSTPC = new TCanvas("canvDcaxyReskAnyITSTPC", "", 500, 500); - canvDcaxyReskAnyITSTPC->DrawFrame(ptLimits[0], 1., ptLimits[nPtBins - 1], 1.e3, "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)"); - canvDcaxyReskAnyITSTPC->SetLogx(); - canvDcaxyReskAnyITSTPC->SetLogy(); - canvDcaxyReskAnyITSTPC->SetGrid(); - for (const auto& pdgCode : pdgCodes) { - hDcaxyReskAnyITSTPC[pdgCode]->Draw("same"); - } - - auto canvDcazReskAnyITSTPC = new TCanvas("canvDcazReskAnyITSTPC", "", 500, 500); - canvDcazReskAnyITSTPC->DrawFrame(ptLimits[0], 1., ptLimits[nPtBins - 1], 1.e3, "ITS-TPC;#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)"); - canvDcazReskAnyITSTPC->SetLogx(); - canvDcazReskAnyITSTPC->SetLogy(); - canvDcazReskAnyITSTPC->SetGrid(); - for (const auto& pdgCode : pdgCodes) { - hDcazReskAnyITSTPC[pdgCode]->Draw("same"); - } - - // Compare ITS-TPC resolution; - auto canvDcaxyResComp = new TCanvas("canvDcaxyResAllLayersComp", "", 500, 500); - canvDcaxyResComp->DrawFrame(ptLimits[0], 1., ptLimits[nPtBins - 1], 1.e3, "ITS vs. ITS-TPC (all layers);#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum)"); - canvDcaxyResComp->SetLogx(); - canvDcaxyResComp->SetLogy(); - canvDcaxyResComp->SetGrid(); - hDcaxyResAllLayersITS[211]->Draw("same"); - hDcaxyResAllLayersITSTPC[211]->Draw("same"); - gPad->BuildLegend(0.8, 0.8, 0.94, 0.94); - - auto canvDcazResComp = new TCanvas("canvDcazResAllLayersComp", "", 500, 500); - canvDcazResComp->DrawFrame(ptLimits[0], 1., ptLimits[nPtBins - 1], 1.e3, "ITS vs. ITS-TPC (all layers);#it{p}_{T} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum)"); - canvDcazResComp->SetLogx(); - canvDcazResComp->SetLogy(); - canvDcazResComp->SetGrid(); - hDcazResAllLayersITS[211]->Draw("same"); - hDcazResAllLayersITSTPC[211]->Draw("same"); - gPad->BuildLegend(0.8, 0.8, 0.94, 0.94); - - auto canvPtResComp = new TCanvas("canvPtResAllLayersComp", "", 500, 500); - canvPtResComp->DrawFrame(ptLimits[0], 0., ptLimits[nPtBins - 1], 0.2, "ITS vs. ITS-TPC (all layers);#it{p}_{T} (GeV/#it{c});#sigma(#Delta#it{p}_{T}/#it{p}_{T})"); - canvPtResComp->SetLogx(); - canvPtResComp->SetGrid(); - hPtResAllLayersITS[211]->Draw("same"); - hPtResAllLayersITSTPC[211]->Draw("same"); - gPad->BuildLegend(0.8, 0.8, 0.94, 0.94); - - auto canvDcaPtComp = new TCanvas("canvDcaPtAllLayersComp", "", 500, 500); - canvDcaPtComp->Divide(2, 2); - canvDcaPtComp->cd(1); - hDcaxyVsPtAllLayersITS[211]->Draw(); - canvDcaPtComp->cd(2); - hDcazVsPtAllLayersITS[211]->Draw(); - canvDcaPtComp->cd(3); - hDcaxyVsPtAllLayersITSTPC[211]->Draw(); - canvDcaPtComp->cd(4); - hDcazVsPtAllLayersITSTPC[211]->Draw(); - - auto canvDcaPhiComp = new TCanvas("canvDcaPhiAllLayersComp", "", 500, 500); - canvDcaPhiComp->Divide(2, 2); - canvDcaPhiComp->cd(1); - hDcaxyVsPhiAllLayersITS[211]->Draw(); - canvDcaPhiComp->cd(2); - hDcazVsPhiAllLayersITS[211]->Draw(); - canvDcaPhiComp->cd(3); - hDcaxyVsPhiAllLayersITSTPC[211]->Draw(); - canvDcaPhiComp->cd(4); - hDcazVsPhiAllLayersITSTPC[211]->Draw(); - - // Write - TFile outFile("checkDCA.root", "RECREATE"); - outFile.mkdir("ITS"); - outFile.cd("ITS"); - gDirectory->WriteTObject(canvPtResITS); - gDirectory->WriteTObject(canvDcaxyResITS); - gDirectory->WriteTObject(canvDcazResITS); - gDirectory->WriteTObject(canvDcazResNoFirstLayerITS); - gDirectory->WriteTObject(canvDcaxyResNoFirstLayerITS); - gDirectory->WriteTObject(canvDcaxyReskAnyITS); - gDirectory->WriteTObject(canvDcazReskAnyITS); - gDirectory->WriteTObject(canvPtDeltaITS); - gDirectory->WriteTObject(canvDcaVsPtITS); - - outFile.mkdir("ITS-TPC"); - outFile.cd("ITS-TPC"); - gDirectory->WriteTObject(canvPtResITSTPC); - gDirectory->WriteTObject(canvDcaxyResITSTPC); - gDirectory->WriteTObject(canvDcazResITSTPC); - gDirectory->WriteTObject(canvDcazResNoFirstLayerITSTPC); - gDirectory->WriteTObject(canvDcaxyResNoFirstLayerITSTPC); - gDirectory->WriteTObject(canvDcaxyReskAnyITSTPC); - gDirectory->WriteTObject(canvDcazReskAnyITSTPC); - gDirectory->WriteTObject(canvPtDeltaITSTPC); - gDirectory->WriteTObject(canvDcaVsPtITSTPC); - - outFile.mkdir("Compare"); - outFile.cd("Compare"); - gDirectory->WriteTObject(canvDcaxyResComp); - gDirectory->WriteTObject(canvDcazResComp); - gDirectory->WriteTObject(canvPtResComp); - gDirectory->WriteTObject(canvDcaPtComp); - gDirectory->WriteTObject(canvDcaPhiComp); - - for (const auto& pdgCode : pdgCodes) { - const char* dirName = partNames[pdgCode].c_str(); - auto dir = outFile.mkdir(dirName); - outFile.cd(dirName); - - gDirectory->mkdir("ITS"); - gDirectory->cd("ITS"); - gDirectory->WriteTObject(hDeltaPtVsPtAllLayersITS[pdgCode]); - gDirectory->WriteTObject(hDcaxyVsPtAllLayersITS[pdgCode]); - gDirectory->WriteTObject(hDcazVsPtAllLayersITS[pdgCode]); - gDirectory->WriteTObject(hDcazResAllLayersITS[pdgCode]); - gDirectory->WriteTObject(hDcaxyResAllLayersITS[pdgCode]); - gDirectory->WriteTObject(hDcaxyVsPtNoFirstLayerITS[pdgCode]); - gDirectory->WriteTObject(hDcazVsPtNoFirstLayerITS[pdgCode]); - gDirectory->WriteTObject(hDcazResNoFirstLayerITS[pdgCode]); - gDirectory->WriteTObject(hDcaxyResNoFirstLayerITS[pdgCode]); - gDirectory->WriteTObject(hDcaxyVsPhiAllLayersITS[pdgCode]); - gDirectory->WriteTObject(hDcazVsPhiAllLayersITS[pdgCode]); - gDirectory->mkdir("projections"); - gDirectory->cd("projections"); - for (TObject* obj : *lProjITS[pdgCode]) { - obj->Write(); - } - - dir->cd(); - gDirectory->mkdir("ITS-TPC"); - gDirectory->cd("ITS-TPC"); - gDirectory->WriteTObject(hDeltaPtVsPtAllLayersITSTPC[pdgCode]); - gDirectory->WriteTObject(hDcaxyVsPtAllLayersITSTPC[pdgCode]); - gDirectory->WriteTObject(hDcazVsPtAllLayersITSTPC[pdgCode]); - gDirectory->WriteTObject(hDcazResAllLayersITSTPC[pdgCode]); - gDirectory->WriteTObject(hDcaxyResAllLayersITSTPC[pdgCode]); - gDirectory->WriteTObject(hDcaxyVsPtNoFirstLayerITSTPC[pdgCode]); - gDirectory->WriteTObject(hDcazVsPtNoFirstLayerITSTPC[pdgCode]); - gDirectory->WriteTObject(hDcazResNoFirstLayerITSTPC[pdgCode]); - gDirectory->WriteTObject(hDcaxyResNoFirstLayerITSTPC[pdgCode]); - gDirectory->WriteTObject(hDcaxyVsPhiAllLayersITSTPC[pdgCode]); - gDirectory->WriteTObject(hDcazVsPhiAllLayersITSTPC[pdgCode]); - gDirectory->mkdir("projections"); - gDirectory->cd("projections"); - for (TObject* obj : *lProjITSTPC[pdgCode]) { - obj->Write(); - } - } -} diff --git a/Detectors/Upgrades/ITS3/macros/test/CheckDigitsITS3.C b/Detectors/Upgrades/ITS3/macros/test/CheckDigitsITS3.C index 240b1bd344af5..82578cc406f0c 100644 --- a/Detectors/Upgrades/ITS3/macros/test/CheckDigitsITS3.C +++ b/Detectors/Upgrades/ITS3/macros/test/CheckDigitsITS3.C @@ -40,7 +40,7 @@ #include "fairlogger/Logger.h" #endif -void CheckDigitsITS3(std::string digifile = "it3digits.root", std::string hitfile = "o2sim_HitsIT3.root", std::string inputGeom = "", bool batch = false) +void CheckDigitsITS3(bool readFromFile = false, std::string digifile = "it3digits.root", std::string hitfile = "o2sim_HitsIT3.root", std::string inputGeom = "", bool batch = false) { gROOT->SetBatch(batch); gStyle->SetPalette(kRainBow); @@ -53,176 +53,211 @@ void CheckDigitsITS3(std::string digifile = "it3digits.root", std::string hitfil using o2::itsmft::SegmentationAlpide; std::array mMosaixSegmentations{0, 1, 2}; - TFile* f = TFile::Open("CheckDigits.root", "recreate"); - TNtuple* nt = new TNtuple("ntd", "digit ntuple", "id:x:y:z:rowD:colD:rowH:colH:xlH:zlH:xlcH:zlcH:dx:dz"); + TFile* f{nullptr}; + TNtuple* nt{nullptr}; + if (!readFromFile) { + f = TFile::Open("CheckDigits.root", "recreate"); + nt = new TNtuple("ntd", "digit ntuple", "id:x:y:z:rowD:colD:rowH:colH:xlH:zlH:xlcH:zlcH:dx:dz:etaH"); - // Geometry - o2::base::GeometryManager::loadGeometry(inputGeom); - auto* gman = o2::its::GeometryTGeo::Instance(); - gman->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); + // Geometry + o2::base::GeometryManager::loadGeometry(inputGeom); + auto* gman = o2::its::GeometryTGeo::Instance(); + gman->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::L2G)); - // Hits - TFile* hitFile = TFile::Open(hitfile.data()); - TTree* hitTree = (TTree*)hitFile->Get("o2sim"); - int nevH = hitTree->GetEntries(); // hits are stored as one event per entry - std::vector*> hitArray(nevH, nullptr); - std::vector> mc2hitVec(nevH); + // Hits + TFile* hitFile = TFile::Open(hitfile.data()); + TTree* hitTree = (TTree*)hitFile->Get("o2sim"); + int nevH = hitTree->GetEntries(); // hits are stored as one event per entry + std::vector*> hitArray(nevH, nullptr); + std::vector> mc2hitVec(nevH); - // Digits - TFile* digFile = TFile::Open(digifile.data()); - TTree* digTree = (TTree*)digFile->Get("o2sim"); + // Digits + TFile* digFile = TFile::Open(digifile.data()); + TTree* digTree = (TTree*)digFile->Get("o2sim"); - std::vector* digArr = nullptr; - digTree->SetBranchAddress("IT3Digit", &digArr); + std::vector* digArr = nullptr; + digTree->SetBranchAddress("IT3Digit", &digArr); - o2::dataformats::IOMCTruthContainerView* plabels = nullptr; - digTree->SetBranchAddress("IT3DigitMCTruth", &plabels); + o2::dataformats::IOMCTruthContainerView* plabels = nullptr; + digTree->SetBranchAddress("IT3DigitMCTruth", &plabels); - int nevD = digTree->GetEntries(); // digits in cont. readout may be grouped as few events per entry + int nevD = digTree->GetEntries(); // digits in cont. readout may be grouped as few events per entry - int nDigitReadIB{0}, nDigitReadOB{0}; - int nDigitFilledIB{0}, nDigitFilledOB{0}; + int nDigitReadIB{0}, nDigitReadOB{0}; + int nDigitFilledIB{0}, nDigitFilledOB{0}; - // Get Read Out Frame arrays - std::vector* ROFRecordArrray = nullptr; - digTree->SetBranchAddress("IT3DigitROF", &ROFRecordArrray); - std::vector& ROFRecordArrrayRef = *ROFRecordArrray; + // Get Read Out Frame arrays + std::vector* ROFRecordArrray = nullptr; + digTree->SetBranchAddress("IT3DigitROF", &ROFRecordArrray); + std::vector& ROFRecordArrrayRef = *ROFRecordArrray; - std::vector* MC2ROFRecordArrray = nullptr; - digTree->SetBranchAddress("IT3DigitMC2ROF", &MC2ROFRecordArrray); - std::vector& MC2ROFRecordArrrayRef = *MC2ROFRecordArrray; + std::vector* MC2ROFRecordArrray = nullptr; + digTree->SetBranchAddress("IT3DigitMC2ROF", &MC2ROFRecordArrray); + std::vector& MC2ROFRecordArrrayRef = *MC2ROFRecordArrray; - digTree->GetEntry(0); + digTree->GetEntry(0); - int nROFRec = (int)ROFRecordArrrayRef.size(); - std::vector mcEvMin(nROFRec, hitTree->GetEntries()); - std::vector mcEvMax(nROFRec, -1); - o2::dataformats::ConstMCTruthContainer labels; - plabels->copyandflatten(labels); - delete plabels; + int nROFRec = (int)ROFRecordArrrayRef.size(); + std::vector mcEvMin(nROFRec, hitTree->GetEntries()); + std::vector mcEvMax(nROFRec, -1); + o2::dataformats::ConstMCTruthContainer labels; + plabels->copyandflatten(labels); + delete plabels; - LOGP(debug, "Build min and max MC events used by each ROF"); - for (int imc = MC2ROFRecordArrrayRef.size(); imc--;) { - const auto& mc2rof = MC2ROFRecordArrrayRef[imc]; - /* LOGP(debug, "MCRecord: {}", mc2rof.asString()); */ + LOGP(debug, "Build min and max MC events used by each ROF"); + for (int imc = MC2ROFRecordArrrayRef.size(); imc--;) { + const auto& mc2rof = MC2ROFRecordArrrayRef[imc]; + /* LOGP(debug, "MCRecord: {}", mc2rof.asString()); */ - if (mc2rof.rofRecordID < 0) { - continue; // this MC event did not contribute to any ROF - } + if (mc2rof.rofRecordID < 0) { + continue; // this MC event did not contribute to any ROF + } - for (int irfd = mc2rof.maxROF - mc2rof.minROF + 1; irfd--;) { + for (int irfd = mc2rof.maxROF - mc2rof.minROF + 1; irfd--;) { - int irof = mc2rof.rofRecordID + irfd; + int irof = mc2rof.rofRecordID + irfd; - if (irof >= nROFRec) { - LOGP(error, "ROF={} from MC2ROF record is >= N ROFs={}", irof, nROFRec); - } - if (mcEvMin[irof] > imc) { - mcEvMin[irof] = imc; - } - if (mcEvMax[irof] < imc) { - mcEvMax[irof] = imc; + if (irof >= nROFRec) { + LOGP(error, "ROF={} from MC2ROF record is >= N ROFs={}", irof, nROFRec); + } + if (mcEvMin[irof] > imc) { + mcEvMin[irof] = imc; + } + if (mcEvMax[irof] < imc) { + mcEvMax[irof] = imc; + } } } - } - LOGP(debug, "LOOP on: ROFRecord array"); - unsigned int rofIndex = 0; - unsigned int rofNEntries = 0; - for (unsigned int iROF = 0; iROF < ROFRecordArrrayRef.size(); iROF++) { + LOGP(debug, "LOOP on: ROFRecord array"); + unsigned int rofIndex = 0; + unsigned int rofNEntries = 0; + for (unsigned int iROF = 0; iROF < ROFRecordArrrayRef.size(); iROF++) { - rofIndex = ROFRecordArrrayRef[iROF].getFirstEntry(); - rofNEntries = ROFRecordArrrayRef[iROF].getNEntries(); + rofIndex = ROFRecordArrrayRef[iROF].getFirstEntry(); + rofNEntries = ROFRecordArrrayRef[iROF].getNEntries(); - // >> read and map MC events contributing to this ROF - for (int im = mcEvMin[iROF]; im <= mcEvMax[iROF]; im++) { + // >> read and map MC events contributing to this ROF + for (int im = mcEvMin[iROF]; im <= mcEvMax[iROF]; im++) { - if (hitArray[im] == nullptr) { + if (hitArray[im] == nullptr) { - hitTree->SetBranchAddress("IT3Hit", &hitArray[im]); - hitTree->GetEntry(im); + hitTree->SetBranchAddress("IT3Hit", &hitArray[im]); + hitTree->GetEntry(im); - auto& mc2hit = mc2hitVec[im]; + auto& mc2hit = mc2hitVec[im]; - for (size_t ih = hitArray[im]->size(); ih--;) { - const auto& hit = (*hitArray[im])[ih]; - uint64_t key = (uint64_t(hit.GetTrackID()) << 32) + hit.GetDetectorID(); - mc2hit.emplace(key, ih); + for (size_t ih = hitArray[im]->size(); ih--;) { + const auto& hit = (*hitArray[im])[ih]; + uint64_t key = (uint64_t(hit.GetTrackID()) << 32) + hit.GetDetectorID(); + mc2hit.emplace(key, ih); + } } } - } - - LOGP(debug, " `-> LOOP on: Digits array(size={}) starting at ROFIndex={} to {}", digArr->size(), rofIndex, rofIndex + rofNEntries); - for (unsigned int iDigit = rofIndex; iDigit < rofIndex + rofNEntries; iDigit++) { - int ix = (*digArr)[iDigit].getRow(), iz = (*digArr)[iDigit].getColumn(); - auto chipID = (*digArr)[iDigit].getChipIndex(); - auto layer = its3::constants::detID::getDetID2Layer(chipID); - bool isIB{its3::constants::detID::isDetITS3(chipID)}; - float x{0.f}, y{0.f}, z{0.f}; - (isIB) ? ++nDigitReadIB : ++nDigitReadOB; - - if (isIB) { - // ITS3 IB - float xFlat{0.f}, yFlat{0.f}; - mMosaixSegmentations[layer].detectorToLocal(ix, iz, xFlat, z); - mMosaixSegmentations[layer].flatToCurved(xFlat, 0., x, y); - } else { - // ITS2 OB - SegmentationAlpide::detectorToLocal(ix, iz, x, z); - } - o2::math_utils::Point3D locD(x, y, z); - auto lab = (labels.getLabels(iDigit))[0]; - int trID = lab.getTrackID(); - if (!lab.isValid()) { // not noise - continue; - } - - // get MC info - uint64_t key = (uint64_t(trID) << 32) + chipID; - const auto* mc2hit = &mc2hitVec[lab.getEventID()]; - const auto& hitEntry = mc2hit->find(key); - if (hitEntry == mc2hit->end()) { - LOGP(debug, "Failed to find MC hit entry for Tr {} chipID {}", trID, chipID); - continue; - } + LOGP(debug, " `-> LOOP on: Digits array(size={}) starting at ROFIndex={} to {}", digArr->size(), rofIndex, rofIndex + rofNEntries); + for (unsigned int iDigit = rofIndex; iDigit < rofIndex + rofNEntries; iDigit++) { + int ix = (*digArr)[iDigit].getRow(), iz = (*digArr)[iDigit].getColumn(); + auto chipID = (*digArr)[iDigit].getChipIndex(); + auto layer = its3::constants::detID::getDetID2Layer(chipID); + bool isIB{its3::constants::detID::isDetITS3(chipID)}; + float x{0.f}, y{0.f}, z{0.f}; + (isIB) ? ++nDigitReadIB : ++nDigitReadOB; + + if (isIB) { + // ITS3 IB + float xFlat{0.f}, yFlat{0.f}; + mMosaixSegmentations[layer].detectorToLocal(ix, iz, xFlat, z); + mMosaixSegmentations[layer].flatToCurved(xFlat, 0., x, y); + } else { + // ITS2 OB + SegmentationAlpide::detectorToLocal(ix, iz, x, z); + } - auto gloD = gman->getMatrixL2G(chipID)(locD); // convert to global - - ////// HITS - Hit& hit = (*hitArray[lab.getEventID()])[hitEntry->second]; - - auto xyzLocE = gman->getMatrixL2G(chipID) ^ (hit.GetPos()); // inverse conversion from global to local - auto xyzLocS = gman->getMatrixL2G(chipID) ^ (hit.GetPosStart()); - o2::math_utils::Vector3D xyzLocM; - xyzLocM.SetCoordinates(0.5f * (xyzLocE.X() + xyzLocS.X()), 0.5f * (xyzLocE.Y() + xyzLocS.Y()), 0.5f * (xyzLocE.Z() + xyzLocS.Z())); - float xlc = 0., zlc = 0.; - int row = 0, col = 0; - - if (isIB) { - float xFlat{0.}, yFlat{0.}; - mMosaixSegmentations[layer].curvedToFlat(xyzLocM.X(), xyzLocM.Y(), xFlat, yFlat); - xyzLocM.SetCoordinates(xFlat, yFlat, xyzLocM.Z()); - mMosaixSegmentations[layer].curvedToFlat(locD.X(), locD.Y(), xFlat, yFlat); - locD.SetCoordinates(xFlat, yFlat, locD.Z()); - if (auto v1 = !mMosaixSegmentations[layer].localToDetector(xyzLocM.X(), xyzLocM.Z(), row, col), - v2 = !mMosaixSegmentations[layer].detectorToLocal(row, col, xlc, zlc); - v1 || v2) { + o2::math_utils::Point3D locD(x, y, z); + auto lab = (labels.getLabels(iDigit))[0]; + int trID = lab.getTrackID(); + if (!lab.isValid()) { // not noise continue; } - } else { - if (auto v1 = !SegmentationAlpide::localToDetector(xyzLocM.X(), xyzLocM.Z(), row, col), - v2 = !SegmentationAlpide::detectorToLocal(row, col, xlc, zlc); - v1 || v2) { + + // get MC info + uint64_t key = (uint64_t(trID) << 32) + chipID; + const auto* mc2hit = &mc2hitVec[lab.getEventID()]; + const auto& hitEntry = mc2hit->find(key); + if (hitEntry == mc2hit->end()) { + LOGP(debug, "Failed to find MC hit entry for Tr {} chipID {}", trID, chipID); continue; } - } - nt->Fill(chipID, gloD.X(), gloD.Y(), gloD.Z(), ix, iz, row, col, xyzLocM.X(), xyzLocM.Z(), xlc, zlc, xyzLocM.X() - locD.X(), xyzLocM.Z() - locD.Z()); + auto gloD = gman->getMatrixL2G(chipID)(locD); // convert to global + + ////// HITS + Hit& hit = (*hitArray[lab.getEventID()])[hitEntry->second]; + // mean local position of the hit + auto locH = gman->getMatrixL2G(chipID) ^ (hit.GetPos()); // inverse conversion from global to local + auto locHsta = gman->getMatrixL2G(chipID) ^ (hit.GetPosStart()); + o2::math_utils::Point3D locHmid; + float x0, y0, z0, dltx, dlty, dltz, r; + if (isIB) { + float xFlat{0.}, yFlat{0.}; + mMosaixSegmentations[layer].curvedToFlat(locH.X(), locH.Y(), xFlat, yFlat); + locH.SetCoordinates(xFlat, yFlat, locH.Z()); + mMosaixSegmentations[layer].curvedToFlat(locHsta.X(), locHsta.Y(), xFlat, yFlat); + locHsta.SetCoordinates(xFlat, yFlat, locHsta.Z()); + x0 = locHsta.X(); + dltx = locH.X() - x0; + y0 = locHsta.Y(); + dlty = locH.Y() - y0; + z0 = locHsta.Z(); + dltz = locH.Z() - z0; + r = (o2::its3::constants::pixelarray::pixels::apts::responseYShift - y0) / dlty; + } else { + x0 = locHsta.X(); + dltx = locH.X() - x0; + y0 = locHsta.Y(); + dlty = locH.Y() - y0; + z0 = locHsta.Z(); + dltz = locH.Z() - z0; + r = (0.5 * (o2::itsmft::SegmentationAlpide::SensorLayerThickness - o2::itsmft::SegmentationAlpide::SensorLayerThicknessEff) - y0) / dlty; + } + locHmid.SetXYZ(x0 + r * dltx, y0 + r * dlty, z0 + r * dltz); + auto gloHmid = gman->getMatrixL2G(chipID) * locHmid; + float theta = std::acos(gloHmid.Z() / gloHmid.Rho()); + float eta = -std::log(std::tan(theta / 2.f)); + + float xlc = 0., zlc = 0.; + int row = 0, col = 0; + + if (isIB) { + float xFlat{0.}, yFlat{0.}; + mMosaixSegmentations[layer].curvedToFlat(locD.X(), locD.Y(), xFlat, yFlat); + locD.SetCoordinates(xFlat, yFlat, locD.Z()); + if (auto v1 = !mMosaixSegmentations[layer].localToDetector(locHmid.X(), locHmid.Z(), row, col), + v2 = !mMosaixSegmentations[layer].detectorToLocal(row, col, xlc, zlc); + v1 || v2) { + continue; + } + } else { + if (auto v1 = !SegmentationAlpide::localToDetector(locHmid.X(), locHmid.Z(), row, col), + v2 = !SegmentationAlpide::detectorToLocal(row, col, xlc, zlc); + v1 || v2) { + continue; + } + } + + nt->Fill(chipID, gloD.X(), gloD.Y(), gloD.Z(), ix, iz, row, col, locHmid.X(), locHmid.Z(), xlc, zlc, locHmid.X() - locD.X(), locHmid.Z() - locD.Z(), eta); - (isIB) ? ++nDigitFilledIB : ++nDigitFilledOB; - } // end loop on digits array - } // end loop on ROFRecords array + (isIB) ? ++nDigitFilledIB : ++nDigitFilledOB; + } // end loop on digits array + } // end loop on ROFRecords array + + Info("EXIT", "read %d filled %d in IB\n", nDigitReadIB, nDigitFilledIB); + Info("EXIT", "read %d filled %d in OB\n", nDigitReadOB, nDigitFilledOB); + } else { + f = TFile::Open("CheckDigits.root", "Open"); + nt = f->Get("ntd"); + } auto canvXY = new TCanvas("canvXY", "", 1600, 1600); canvXY->Divide(2, 2); @@ -234,7 +269,7 @@ void CheckDigitsITS3(std::string digifile = "it3digits.root", std::string hitfil nt->Draw("y:x>>h_y_vs_x_OB(1000, -50, 50, 1000, -50, 50)", "id >= 3456", "colz"); canvXY->cd(4); nt->Draw("y:z>>h_y_vs_z_OB(1000, -100, 100, 1000, -50, 50)", "id >= 3456", "colz"); - canvXY->SaveAs("it3digits_y_vs_x_vs_z.pdf"); + canvXY->SaveAs("it3digits_y_vs_x_vs_z.png"); auto canvdXdZ = new TCanvas("canvdXdZ", "", 1600, 800); canvdXdZ->Divide(2, 2); @@ -258,10 +293,8 @@ void CheckDigitsITS3(std::string digifile = "it3digits.root", std::string hitfil h = (TH2F*)gPad->GetPrimitive("h_dx_vs_dz_OB_z"); Info("OB |z|<2", "RMS(dx)=%.1f mu", h->GetRMS(2) * 1e4); Info("OB |z|<2", "RMS(dz)=%.1f mu", h->GetRMS(1) * 1e4); - canvdXdZ->SaveAs("it3digits_dx_vs_dz.pdf"); + canvdXdZ->SaveAs("it3digits_dx_vs_dz.png"); f->Write(); f->Close(); - Info("EXIT", "read %d filled %d in IB\n", nDigitReadIB, nDigitFilledIB); - Info("EXIT", "read %d filled %d in OB\n", nDigitReadOB, nDigitFilledOB); } diff --git a/Detectors/Upgrades/ITS3/study/CMakeLists.txt b/Detectors/Upgrades/ITS3/study/CMakeLists.txt new file mode 100644 index 0000000000000..4bb1cbca7dcb0 --- /dev/null +++ b/Detectors/Upgrades/ITS3/study/CMakeLists.txt @@ -0,0 +1,37 @@ +# Copyright 2019-2020 CERN and copyright holders of ALICE O2. +# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +# All rights not expressly granted are reserved. +# +# This software is distributed under the terms of the GNU General Public +# License v3 (GPL Version 3), copied verbatim in the file "COPYING". +# +# In applying this license CERN does not waive the privileges and immunities +# granted to it by virtue of its status as an Intergovernmental Organization +# or submit itself to any jurisdiction. + +# add_compile_options(-O0 -g -fPIC -fno-omit-frame-pointer) + +o2_add_library(ITS3TrackingStudy + TARGETVARNAME targetName + SOURCES src/ITS3TrackingStudyParam.cxx + src/TrackingStudy.cxx + src/ParticleInfoExt.cxx + PUBLIC_LINK_LIBRARIES O2::ITS3Workflow + O2::GlobalTracking + O2::GlobalTrackingWorkflowReaders + O2::GlobalTrackingWorkflowHelpers + O2::DataFormatsGlobalTracking + O2::DetectorsVertexing + O2::SimulationDataFormat) + +o2_target_root_dictionary(ITS3TrackingStudy + HEADERS include/ITS3TrackingStudy/ITS3TrackingStudyParam.h + include/ITS3TrackingStudy/ParticleInfoExt.h + LINKDEF src/ITS3TrackingStudyLinkDef.h) + +o2_add_executable(study-workflow + COMPONENT_NAME its3-tracking + SOURCES src/its3-tracking-study-workflow.cxx + PUBLIC_LINK_LIBRARIES O2::ITS3TrackingStudy) + +add_subdirectory(macros) diff --git a/Detectors/Upgrades/ITS3/study/include/ITS3TrackingStudy/ITS3TrackingStudyParam.h b/Detectors/Upgrades/ITS3/study/include/ITS3TrackingStudy/ITS3TrackingStudyParam.h new file mode 100644 index 0000000000000..2e718622daa90 --- /dev/null +++ b/Detectors/Upgrades/ITS3/study/include/ITS3TrackingStudy/ITS3TrackingStudyParam.h @@ -0,0 +1,49 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef O2_TRACKING_STUDY_CONFIG_H +#define O2_TRACKING_STUDY_CONFIG_H +#include "CommonUtils/ConfigurableParam.h" +#include "CommonUtils/ConfigurableParamHelper.h" + +#include "DetectorsBase/Propagator.h" + +namespace o2::its3::study +{ + +struct ITS3TrackingStudyParam : o2::conf::ConfigurableParamHelper { + /// general track selection + float maxChi2{36}; + float maxEta{1.0}; + float minPt{0.1}; + float maxPt{1e2}; + /// PV selection + int minPVCont{5}; + /// ITS track selection + int minITSCls{7}; + bool refitITS{true}; // refit ITS track including the PV + /// TPC track selection + int minTPCCls{110}; + + // propagator + o2::base::PropagatorImpl::MatCorrType CorrType = o2::base::PropagatorImpl::MatCorrType::USEMatCorrLUT; + + /// studies + bool doDCA = true; + bool doDCARefit = true; + bool doPull = true; + bool doMC = false; + O2ParamDef(ITS3TrackingStudyParam, "ITS3TrackingStudyParam"); +}; + +} // namespace o2::its3::study + +#endif diff --git a/Detectors/Upgrades/ITS3/study/include/ITS3TrackingStudy/ParticleInfoExt.h b/Detectors/Upgrades/ITS3/study/include/ITS3TrackingStudy/ParticleInfoExt.h new file mode 100644 index 0000000000000..c66068418377d --- /dev/null +++ b/Detectors/Upgrades/ITS3/study/include/ITS3TrackingStudy/ParticleInfoExt.h @@ -0,0 +1,42 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef ALICEO2_PARTICLEINFO_EXT_H +#define ALICEO2_PARTICLEINFO_EXT_H + +#include "ReconstructionDataFormats/Track.h" +#include "ReconstructionDataFormats/GlobalTrackID.h" +#include "SimulationDataFormat/MCTrack.h" + +namespace o2::its3::study +{ + +struct ParticleInfoExt { + // cluster info + uint8_t clusters{0}; + uint8_t fakeClusters{0}; + // reco info + uint8_t isReco{0}; + uint8_t isFake{0}; + // matching info + uint8_t recoTracks; + uint8_t fakeTracks; + // reco track + track::TrackParCov recoTrack; + // mc info + MCTrack mcTrack; + + ClassDefNV(ParticleInfoExt, 1); +}; + +} // namespace o2::its3::study + +#endif diff --git a/Detectors/Upgrades/ITS3/study/include/ITS3TrackingStudy/TrackingStudy.h b/Detectors/Upgrades/ITS3/study/include/ITS3TrackingStudy/TrackingStudy.h new file mode 100644 index 0000000000000..065629058fd32 --- /dev/null +++ b/Detectors/Upgrades/ITS3/study/include/ITS3TrackingStudy/TrackingStudy.h @@ -0,0 +1,25 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef O2_ITS3_TRACKING_STUDY_H +#define O2_ITS3_TRACKING_STUDY_H + +#include "ReconstructionDataFormats/GlobalTrackID.h" +#include "Framework/DataProcessorSpec.h" + +namespace o2::its3::study +{ + +o2::framework::DataProcessorSpec getTrackingStudySpec(o2::dataformats::GlobalTrackID::mask_t srcTracks, o2::dataformats::GlobalTrackID::mask_t srcClus, bool useMC); + +} // namespace o2::its3::study + +#endif diff --git a/Detectors/Upgrades/ITS3/study/macros/CMakeLists.txt b/Detectors/Upgrades/ITS3/study/macros/CMakeLists.txt new file mode 100644 index 0000000000000..aaf763888c5e0 --- /dev/null +++ b/Detectors/Upgrades/ITS3/study/macros/CMakeLists.txt @@ -0,0 +1,18 @@ +# Copyright 2019-2020 CERN and copyright holders of ALICE O2. +# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +# All rights not expressly granted are reserved. +# +# This software is distributed under the terms of the GNU General Public +# License v3 (GPL Version 3), copied verbatim in the file "COPYING". +# +# In applying this license CERN does not waive the privileges and immunities +# granted to it by virtue of its status as an Intergovernmental Organization +# or submit itself to any jurisdiction. + +o2_add_test_root_macro(PlotDCA.C + PUBLIC_LINK_LIBRARIES O2::ITS3TrackingStudy + LABELS its COMPILE_ONLY) + +o2_add_test_root_macro(PlotPulls.C + PUBLIC_LINK_LIBRARIES O2::ITS3TrackingStudy + LABELS its COMPILE_ONLY) diff --git a/Detectors/Upgrades/ITS3/study/macros/PlotDCA.C b/Detectors/Upgrades/ITS3/study/macros/PlotDCA.C new file mode 100644 index 0000000000000..ac92fa491c1ac --- /dev/null +++ b/Detectors/Upgrades/ITS3/study/macros/PlotDCA.C @@ -0,0 +1,190 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file PlotDCA.C +/// \brief Simple macro to plot ITS3 impact parameter resolution + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include + +#include +#include +#include +#include +#include +#include + +#include "ReconstructionDataFormats/GlobalTrackID.h" +#include "ReconstructionDataFormats/Track.h" +#include "ReconstructionDataFormats/DCA.h" +#include "SimulationDataFormat/MCTrack.h" +#endif + +using GTrackID = o2::dataformats::GlobalTrackID; + +static std::string SanitizeSourceName(std::string_view raw) +{ + std::string s(raw); + s.erase(std::remove(s.begin(), s.end(), '-'), s.end()); + return s; +} + +void PlotDCA(const char* fName = "its3TrackStudy.root") +{ + TH1::SetDefaultSumw2(); + std::unique_ptr inFile(TFile::Open(fName)); + auto tree = inFile->Get("dca"); + + int src; // track type + tree->SetBranchAddress("src", &src); + o2::dataformats::DCA* dca{nullptr}; + tree->SetBranchAddress("dca2MC", &dca); + o2::track::TrackParCov* trk{nullptr}; + tree->SetBranchAddress("trk", &trk); + o2::MCTrack* mcTrk{nullptr}; + tree->SetBranchAddress("mcTrk", &mcTrk); + + const int nPtBins = 35; + const double ptLimits[nPtBins] = {0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2., 2.2, 2.5, 3., 4., 5., 6., 8., 10., 15., 20.}; + const int yDCABins{1000}; + const float yDCARange{500}; + const int nSpecies = 5; + std::array pdgCodes{-1, 11, 211, 321, 2212}; + auto fGaus = new TF1("fGaus", "gaus", -200., 200.); + std::map partNames = { + {-1, "All"}, + {11, "Electrons"}, + {211, "Pions"}, + {321, "Kaons"}, + {2212, "Protons"}}; + + std::map> hMapDCAxyVsPtAllLayers; // species -> [src, {dca}] + std::map> hMapResDCAxyVsPtAllLayers; // species -> [src, {dca}] + std::map> hMapDCAzVsPtAllLayers; // species -> [src, {dca}] + std::map> hMapResDCAzVsPtAllLayers; // species -> [src, {dca}] + std::map> hMapDeltaPtVsPtAllLayers; // species -> [src, {dca}] + std::map> hMapResPtVsPtAllLayers; // species -> [src, {dca}] + for (const auto& [sPDG, sName] : partNames) { + std::map histsDCAxy, histsDCAz, histsDeltaPt; + std::map histsResDCAxy, histsResDCAz, histsResDeltaPt; + + for (int cis = 0; cis < GTrackID::NSources; ++cis) { + const auto cdm = GTrackID::getSourceDetectorsMask(cis); + if (!cdm[GTrackID::ITS]) { + continue; // keep same logic as original + } + + const std::string srcRaw = GTrackID::getSourceName(cis); + const std::string src = SanitizeSourceName(srcRaw); + + histsDCAxy[cis] = new TH2F(Form("hDCAxyVsPtAllLayers_%s_%s", sName.c_str(), src.c_str()), Form("%s;#it{p}_{T,MC} (GeV/#it{c});DCA_{#it{xy}} (#mum);entries", srcRaw.c_str()), nPtBins - 1, ptLimits, yDCABins, -yDCARange, yDCARange); + + histsResDCAxy[cis] = new TH1F(Form("hResDCAxyVsPtAllLayers_%s_%s", sName.c_str(), src.c_str()), Form("%s;#it{p}_{T,MC} (GeV/#it{c});#sigma(DCA_{#it{xy}}) (#mum);entries", srcRaw.c_str()), nPtBins - 1, ptLimits); + + histsDCAz[cis] = new TH2F(Form("hDCAzVsPtAllLayers_%s_%s", sName.c_str(), src.c_str()), Form("%s;#it{p}_{T,MC} (GeV/#it{c});DCA_{#it{z}} (#mum);entries", srcRaw.c_str()), nPtBins - 1, ptLimits, yDCABins, -yDCARange, yDCARange); + + histsResDCAz[cis] = new TH1F(Form("hResDCAzVsPtAllLayers_%s_%s", sName.c_str(), src.c_str()), Form("%s;#it{p}_{T,MC} (GeV/#it{c});#sigma(DCA_{#it{z}}) (#mum);entries", srcRaw.c_str()), nPtBins - 1, ptLimits); + + histsDeltaPt[cis] = new TH2F(Form("hDeltaPtVsPtAllLayers_%s_%s", sName.c_str(), src.c_str()), Form("%s;#it{p}_{T,MC} (GeV/#it{c});#Delta_{#it{p}_{T}}/#it{p}_{T};entries", srcRaw.c_str()), nPtBins - 1, ptLimits, 200, -0.2, 0.2); + + histsResDeltaPt[cis] = new TH1F(Form("hResDeltaPtVsPtAllLayers_%s_%s", sName.c_str(), src.c_str()), Form("%s;#it{p}_{T,MC} (GeV/#it{c});#sigma(#Delta#it{p}_{T}/#it{p}_{T});entries", srcRaw.c_str()), nPtBins - 1, ptLimits); + } + + hMapDCAxyVsPtAllLayers[sPDG] = std::move(histsDCAxy); + hMapResDCAxyVsPtAllLayers[sPDG] = std::move(histsResDCAxy); + hMapDCAzVsPtAllLayers[sPDG] = std::move(histsDCAz); + hMapResDCAzVsPtAllLayers[sPDG] = std::move(histsResDCAz); + hMapDeltaPtVsPtAllLayers[sPDG] = std::move(histsDeltaPt); + hMapResPtVsPtAllLayers[sPDG] = std::move(histsResDeltaPt); + } + + for (int iEntry = 0; tree->LoadTree(iEntry) >= 0; ++iEntry) { + tree->GetEntry(iEntry); + if (!mcTrk->isPrimary()) { + continue; + } + auto pdg = std::abs(mcTrk->GetPdgCode()); + if (pdg != 11 && pdg != 211 && pdg != 321 && pdg != 2212) { + continue; + } + auto ptReco = trk->getPt(); + auto ptGen = mcTrk->GetPt(); + auto deltaPt = (1. / ptReco - 1. / ptGen) / (1. / ptGen); + auto dcaXY = dca->getY() * 10000.; + auto dcaZ = dca->getZ() * 10000.; + auto phiReco = trk->getPhi(); + + for (int spe : {-1, pdg}) { + hMapDeltaPtVsPtAllLayers[spe][src]->Fill(ptGen, deltaPt); + hMapDCAxyVsPtAllLayers[spe][src]->Fill(ptGen, dcaXY); + hMapDCAzVsPtAllLayers[spe][src]->Fill(ptGen, dcaZ); + } + } + + const char* fitOpt{"QWMER"}; + for (const auto& [sPDG, sName] : partNames) { + for (int cis = 0; cis < GTrackID::NSources; cis++) { + const auto cdm = GTrackID::getSourceDetectorsMask(cis); + if (!cdm[GTrackID::ITS]) { + continue; + } + for (auto iPt{0}; iPt < nPtBins; ++iPt) { + auto ptMin = hMapDCAxyVsPtAllLayers[sPDG][cis]->GetXaxis()->GetBinLowEdge(iPt + 1); + float minFit = (ptMin < 1.) ? -200. : -50.; + float maxFit = (ptMin < 1.) ? 200. : 50.; + auto doProjection = [&](auto& hIn, auto& hOut, bool useRange = true) { + auto hProj = hIn[sPDG][cis]->ProjectionY(Form("%s_%d", hOut[sPDG][cis]->GetName(), iPt), iPt + 1, iPt + 1); + if (hProj->GetEntries() < 100) { + return; + } + if (useRange) { + hProj->Fit("fGaus", fitOpt, "", minFit, maxFit); + } else { + hProj->Fit("fGaus", fitOpt); + } + hOut[sPDG][cis]->SetBinContent(iPt + 1, fGaus->GetParameter(2)); + hOut[sPDG][cis]->SetBinError(iPt + 1, fGaus->GetParError(2)); + }; + + doProjection(hMapDeltaPtVsPtAllLayers, hMapResPtVsPtAllLayers, false); + doProjection(hMapDCAxyVsPtAllLayers, hMapResDCAxyVsPtAllLayers); + doProjection(hMapDCAzVsPtAllLayers, hMapResDCAzVsPtAllLayers); + } + } + } + + TFile outFile("plotDCA.root", "RECREATE"); + for (const auto& [sPDG, sName] : partNames) { + outFile.mkdir(sName.c_str()); + outFile.cd(sName.c_str()); + for (int cis = 0; cis < GTrackID::NSources; cis++) { + const auto cdm = GTrackID::getSourceDetectorsMask(cis); + if (!cdm[GTrackID::ITS]) { + continue; + } + const std::string srcRaw = GTrackID::getSourceName(cis); + const std::string src = SanitizeSourceName(srcRaw); + gDirectory->mkdir(src.c_str()); + gDirectory->cd(src.c_str()); + + hMapDCAxyVsPtAllLayers[sPDG][cis]->Write(); + hMapResDCAxyVsPtAllLayers[sPDG][cis]->Write(); + hMapDCAzVsPtAllLayers[sPDG][cis]->Write(); + hMapResDCAzVsPtAllLayers[sPDG][cis]->Write(); + hMapDeltaPtVsPtAllLayers[sPDG][cis]->Write(); + hMapResPtVsPtAllLayers[sPDG][cis]->Write(); + + outFile.cd(sName.c_str()); + } + + outFile.cd(); + } +} diff --git a/Detectors/Upgrades/ITS3/study/macros/PlotPulls.C b/Detectors/Upgrades/ITS3/study/macros/PlotPulls.C new file mode 100644 index 0000000000000..371a94cda0e70 --- /dev/null +++ b/Detectors/Upgrades/ITS3/study/macros/PlotPulls.C @@ -0,0 +1,176 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file PlotDCA.C +/// \brief Simple macro to plot ITS3 impact parameter resolution + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "ReconstructionDataFormats/GlobalTrackID.h" +#include "ReconstructionDataFormats/Track.h" +#include "ReconstructionDataFormats/DCA.h" +#include "SimulationDataFormat/MCTrack.h" +#endif + +// chi2 PDF with amplitude A, degrees of freedom k, scale s +Double_t chi2_pdf(Double_t* x, Double_t* par) +{ + const Double_t xx = x[0]; + const Double_t A = par[0]; + const Double_t k = par[1]; + const Double_t s = par[2]; + if (xx <= 0.0 || k <= 0.0 || s <= 0.0) { + return 0.0; + } + const Double_t coef = 1.0 / (TMath::Power(2.0 * s, k * 0.5) * TMath::Gamma(k * 0.5)); + return A * coef * TMath::Power(xx, (k * 0.5) - 1.0) * TMath::Exp(-xx / (2.0 * s)); +} + +void PlotPulls(const char* fName = "its3TrackStudy.root") +{ + TH1::SetDefaultSumw2(); + std::unique_ptr inFile(TFile::Open(fName)); + auto tree = inFile->Get("pull"); + + uint8_t src; // track type + tree->SetBranchAddress("src", &src); + o2::track::TrackParCov* trk{nullptr}; + tree->SetBranchAddress("trk", &trk); + o2::track::TrackPar* mcTrk{nullptr}; + tree->SetBranchAddress("mcTrk", &mcTrk); + o2::MCTrack* part{nullptr}; + tree->SetBranchAddress("mcPart", &part); + const int nPtBins = 35; + const double ptLimits[nPtBins] = {0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2., 2.2, 2.5, 3., 4., 5., 6., 8., 10., 15., 20.}; + const int yBins{100}, yRange{5}; + const char* pNames[5] = {"Y", "Z", "Snp", "Tgl", "Q2Pt"}; + auto fGaus = new TF1("fGaus", "[0]*exp(-0.5*((x-[1])/[2])**2)", -3., 3.); + + std::array pulls{ + new TH2F("hPullY", "", nPtBins - 1, ptLimits, yBins, -yRange, yRange), + new TH2F("hPullZ", "", nPtBins - 1, ptLimits, yBins, -yRange, yRange), + new TH2F("hPullSnp", "", nPtBins - 1, ptLimits, yBins, -yRange, yRange), + new TH2F("hPullTgl", "", nPtBins - 1, ptLimits, yBins, -yRange, yRange), + new TH2F("hPullQ2Pt", "", nPtBins - 1, ptLimits, yBins, -yRange, yRange)}; + + std::array means{ + new TH1F("hPullYMean", "", nPtBins - 1, ptLimits), + new TH1F("hPullZMean", "", nPtBins - 1, ptLimits), + new TH1F("hPullSnpMean", "", nPtBins - 1, ptLimits), + new TH1F("hPullTglMean", "", nPtBins - 1, ptLimits), + new TH1F("hPullQ2PtMean", "", nPtBins - 1, ptLimits)}; + + std::array sigmas{ + new TH1F("hPullYSigma", "", nPtBins - 1, ptLimits), + new TH1F("hPullZSigma", "", nPtBins - 1, ptLimits), + new TH1F("hPullSnpSigma", "", nPtBins - 1, ptLimits), + new TH1F("hPullTglSigma", "", nPtBins - 1, ptLimits), + new TH1F("hPullQ2PtSigma", "", nPtBins - 1, ptLimits)}; + + auto calcMahalanobisDist2 = [&](const auto* trk, const auto* mc) -> float { + o2::math_utils::SMatrix> cov; + cov(o2::track::kY, o2::track::kY) = trk->getSigmaY2(); + cov(o2::track::kZ, o2::track::kY) = trk->getSigmaZY(); + cov(o2::track::kZ, o2::track::kZ) = trk->getSigmaZ2(); + cov(o2::track::kSnp, o2::track::kY) = trk->getSigmaSnpY(); + cov(o2::track::kSnp, o2::track::kZ) = trk->getSigmaSnpZ(); + cov(o2::track::kSnp, o2::track::kSnp) = trk->getSigmaSnp2(); + cov(o2::track::kTgl, o2::track::kY) = trk->getSigmaTglY(); + cov(o2::track::kTgl, o2::track::kZ) = trk->getSigmaTglZ(); + cov(o2::track::kTgl, o2::track::kSnp) = trk->getSigmaTglSnp(); + cov(o2::track::kTgl, o2::track::kTgl) = trk->getSigmaTgl2(); + cov(o2::track::kQ2Pt, o2::track::kY) = trk->getSigma1PtY(); + cov(o2::track::kQ2Pt, o2::track::kZ) = trk->getSigma1PtZ(); + cov(o2::track::kQ2Pt, o2::track::kSnp) = trk->getSigma1PtSnp(); + cov(o2::track::kQ2Pt, o2::track::kTgl) = trk->getSigma1PtTgl(); + cov(o2::track::kQ2Pt, o2::track::kQ2Pt) = trk->getSigma1Pt2(); + if (!cov.Invert()) { + return -1.f; + } + o2::math_utils::SVector trkPar(trk->getParams(), o2::track::kNParams), mcPar(mc->getParams(), o2::track::kNParams); + auto res = trkPar - mcPar; + return ROOT::Math::Similarity(cov, res); + }; + + auto hMahDist2 = new TH1F("hMahDist2", ";Mahalanobis distance 2;n. entries", 100, 0, 10); + + auto getIndex = [](int i) -> int { return i * (i + 3) / 2; }; + + for (int iEntry = 0; tree->LoadTree(iEntry) >= 0; ++iEntry) { + tree->GetEntry(iEntry); + if (src != o2::dataformats::GlobalTrackID::ITS || std::abs(part->GetPdgCode()) != 211) { + continue; + } + for (int i{0}; i < o2::track::kNParams; ++i) { + pulls[i]->Fill(part->GetPt(), (trk->getParam(i) - mcTrk->getParam(i)) / std::sqrt(trk->getCov()[getIndex(i)])); + } + if (part->GetPt() >= 1.0 && part->GetPt() < 2) { + if (auto dist = calcMahalanobisDist2(trk, mcTrk); dist >= 0.) { + hMahDist2->Fill(dist); + } + } + } + + std::vector projs; + const char* fitOpt{"QWMERSB"}; + for (int i{0}; i < o2::track::kNParams; ++i) { + for (auto iPt{0}; iPt < nPtBins; ++iPt) { + auto hProj = pulls[i]->ProjectionY(Form("%s_%d", pulls[i]->GetName(), iPt), iPt + 1, iPt + 1); + hProj->SetName(Form("p%s_pt%d", pNames[i], iPt)); + hProj->SetTitle(Form("Pull %s #it{p}_{T}#in[%.2f, %.2f)", pNames[i], ptLimits[iPt], ptLimits[iPt + 1])); + projs.push_back(hProj); + if (hProj->GetEntries() < 100) { + return; + } + fGaus->SetParameter(1, 0); + fGaus->SetParameter(2, 1); + auto fRes = hProj->Fit(fGaus, fitOpt); + if (fRes->IsValid() && fGaus->GetParameter(2) > 0) { + means[i]->SetBinContent(iPt + 1, fGaus->GetParameter(1)); + means[i]->SetBinError(iPt + 1, fGaus->GetParError(1)); + sigmas[i]->SetBinContent(iPt + 1, fGaus->GetParameter(2)); + sigmas[i]->SetBinError(iPt + 1, fGaus->GetParError(2)); + } + } + } + + hMahDist2->Scale(1. / hMahDist2->Integral("width")); + TF1* fchi2Fit = new TF1("fchi2_fit", chi2_pdf, 0.1, 6, 3); + fchi2Fit->SetParNames("A", "k", "s"); + fchi2Fit->SetParameter(0, 1); + fchi2Fit->SetParameter(1, 5); + fchi2Fit->SetParameter(2, 1); + auto fitres = hMahDist2->Fit(fchi2Fit, "RMQS"); + fitres->Print(); + + TFile outFile("plotPulls.root", "RECREATE"); + for (int i{0}; i < o2::track::kNParams; ++i) { + pulls[i]->Write(); + means[i]->Write(); + sigmas[i]->Write(); + } + for (const auto& p : projs) { + p->Write(); + } + hMahDist2->Write(); + fchi2Fit->Write(); +} diff --git a/Detectors/Upgrades/ITS3/study/src/ITS3TrackingStudyLinkDef.h b/Detectors/Upgrades/ITS3/study/src/ITS3TrackingStudyLinkDef.h new file mode 100644 index 0000000000000..182ffd858629c --- /dev/null +++ b/Detectors/Upgrades/ITS3/study/src/ITS3TrackingStudyLinkDef.h @@ -0,0 +1,23 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifdef __CLING__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class o2::its3::study::ITS3TrackingStudyParam + ; +#pragma link C++ class o2::conf::ConfigurableParamHelper < o2::its3::study::ITS3TrackingStudyParam> + ; + +#pragma link C++ class o2::its3::study::ParticleInfoExt + ; + +#endif diff --git a/Detectors/Upgrades/ITS3/study/src/ITS3TrackingStudyParam.cxx b/Detectors/Upgrades/ITS3/study/src/ITS3TrackingStudyParam.cxx new file mode 100644 index 0000000000000..00bb800e65f8c --- /dev/null +++ b/Detectors/Upgrades/ITS3/study/src/ITS3TrackingStudyParam.cxx @@ -0,0 +1,13 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "ITS3TrackingStudy/ITS3TrackingStudyParam.h" +O2ParamImpl(o2::its3::study::ITS3TrackingStudyParam); diff --git a/Detectors/Upgrades/ITS3/study/src/ParticleInfoExt.cxx b/Detectors/Upgrades/ITS3/study/src/ParticleInfoExt.cxx new file mode 100644 index 0000000000000..aa5edbf408270 --- /dev/null +++ b/Detectors/Upgrades/ITS3/study/src/ParticleInfoExt.cxx @@ -0,0 +1,13 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "ITS3TrackingStudy/ParticleInfoExt.h" +ClassImp(o2::its3::study::ParticleInfoExt); diff --git a/Detectors/Upgrades/ITS3/study/src/TrackingStudy.cxx b/Detectors/Upgrades/ITS3/study/src/TrackingStudy.cxx new file mode 100644 index 0000000000000..cb1d7f381983d --- /dev/null +++ b/Detectors/Upgrades/ITS3/study/src/TrackingStudy.cxx @@ -0,0 +1,841 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include +#include + +#include +#include + +#include "CommonUtils/TreeStreamRedirector.h" +#include "DataFormatsGlobalTracking/RecoContainer.h" +#include "DataFormatsGlobalTracking/RecoContainerCreateTracksVariadic.h" +#include "DataFormatsITSMFT/Digit.h" +#include "ITSMFTSimulation/Hit.h" +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsBase/GRPGeomHelper.h" +#include "DetectorsBase/Propagator.h" +#include "DetectorsCommonDataFormats/DetID.h" +#include "DetectorsCommonDataFormats/SimTraits.h" +#include "DetectorsVertexing/PVertexer.h" +#include "Framework/CCDBParamSpec.h" +#include "Framework/DeviceSpec.h" +#include "Framework/Task.h" +#include "ITSBase/GeometryTGeo.h" +#include "ITS3Base/SpecsV2.h" +#include "ITS3Reconstruction/TopologyDictionary.h" +#include "ITS3Reconstruction/IOUtils.h" +#include "ITS3TrackingStudy/ITS3TrackingStudyParam.h" +#include "ITS3TrackingStudy/ParticleInfoExt.h" +#include "ReconstructionDataFormats/DCA.h" +#include "ReconstructionDataFormats/GlobalTrackID.h" +#include "ReconstructionDataFormats/PrimaryVertex.h" +#include "ReconstructionDataFormats/PrimaryVertexExt.h" +#include "ReconstructionDataFormats/VtxTrackRef.h" +#include "SimulationDataFormat/MCEventLabel.h" +#include "SimulationDataFormat/MCUtils.h" +#include "Steer/MCKinematicsReader.h" + +namespace o2::its3::study +{ + +using namespace o2::framework; +using DetID = o2::detectors::DetID; +using DataRequest = o2::globaltracking::DataRequest; +using PVertex = o2::dataformats::PrimaryVertex; +using GTrackID = o2::dataformats::GlobalTrackID; +using VtxTrackID = o2::dataformats::VtxTrackIndex; +using T2VMap = std::unordered_map; + +class TrackingStudySpec : public Task +{ + public: + TrackingStudySpec(std::shared_ptr dr, std::shared_ptr gr, GTrackID::mask_t src, bool useMC) + : mDataRequest(dr), mGGCCDBRequest(gr), mTracksSrc(src), mUseMC(useMC) {} + ~TrackingStudySpec() final = default; + void init(InitContext& ic) final; + void run(ProcessingContext& pc) final; + void endOfStream(EndOfStreamContext& ec) final; + void finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) final; + + private: + void process(o2::globaltracking::RecoContainer& recoData); + void updateTimeDependentParams(ProcessingContext& pc); + std::vector> prepareITSClusters(const o2::globaltracking::RecoContainer& data) const; + bool selectTrack(GTrackID trkID, o2::globaltracking::RecoContainer& recoData, bool checkMCTruth = true) const; + T2VMap buildT2V(o2::globaltracking::RecoContainer& recoData, bool includeCont = false, bool requireMCMatch = true) const; + bool refitITSPVTrack(o2::globaltracking::RecoContainer& recoData, o2::track::TrackParCov& trFit, GTrackID gidx); + void doDCAStudy(o2::globaltracking::RecoContainer& recoData); + void doDCARefitStudy(o2::globaltracking::RecoContainer& recoData); + void doPullStudy(o2::globaltracking::RecoContainer& recoData); + void doMCStudy(o2::globaltracking::RecoContainer& recoData); + + struct TrackCounter { + TrackCounter() = default; + + void operator+=(int src) + { + if (src >= 0 && src < static_cast(mSuccess.size())) { + ++mSuccess[src]; + } + } + + void operator-=(int src) + { + if (src >= 0 && src < static_cast(mirrors.size())) { + ++mirrors[src]; + } + } + + void operator&=(int src) + { + if (src >= 0 && src < static_cast(mRejected.size())) { + ++mRejected[src]; + } + } + + void print() const + { + LOGP(info, "\t\t\tSuccess / Error / Rejected"); + for (int cis = 0; cis < GTrackID::NSources; ++cis) { + const auto cdm = GTrackID::getSourceDetectorsMask(cis); + if (cdm[DetID::ITS]) { + LOGP(info, "\t{:{}}\t{} / {} / {}", GTrackID::getSourceName(cis), 15, mSuccess[cis], mirrors[cis], mRejected[cis]); + } + } + } + + void reset() + { + mSuccess.fill(0); + mirrors.fill(0); + mRejected.fill(0); + } + + std::array mSuccess{}; + std::array mirrors{}; + std::array mRejected{}; + }; + TrackCounter mTrackCounter; + + std::unique_ptr mDBGOut; + std::shared_ptr mDataRequest; + std::shared_ptr mGGCCDBRequest; + bool mUseMC{false}; + GTrackID::mask_t mTracksSrc; + o2::vertexing::PVertexer mVertexer; + o2::steer::MCKinematicsReader mcReader; // reader of MC information + const o2::its3::TopologyDictionary* mITSDict = nullptr; // cluster patterns dictionary +}; + +void TrackingStudySpec::init(InitContext& ic) +{ + o2::base::GRPGeomHelper::instance().setRequest(mGGCCDBRequest); + int lane = ic.services().get().inputTimesliceId; + int maxLanes = ic.services().get().maxInputTimeslices; + std::string dbgnm = maxLanes == 1 ? "its3TrackStudy.root" : fmt::format("its3TrackStudy_{}.root", lane); + mDBGOut = std::make_unique(dbgnm.c_str(), "recreate"); + + if (mUseMC && !mcReader.initFromDigitContext(o2::base::NameConf::getCollisionContextFileName())) { + LOGP(fatal, "initialization of MCKinematicsReader failed"); + } +} + +void TrackingStudySpec::run(ProcessingContext& pc) +{ + o2::globaltracking::RecoContainer recoData; + recoData.collectData(pc, *mDataRequest); + updateTimeDependentParams(pc); + process(recoData); +} + +void TrackingStudySpec::updateTimeDependentParams(ProcessingContext& pc) +{ + o2::base::GRPGeomHelper::instance().checkUpdates(pc); + if (static bool initOnceDone{false}; !initOnceDone) { // this params need to be queried only once + initOnceDone = true; + auto grp = o2::base::GRPGeomHelper::instance().getGRPECS(); + mVertexer.init(); + o2::its::GeometryTGeo::Instance()->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::L2G, o2::math_utils::TransformType::T2G)); + } +} + +void TrackingStudySpec::endOfStream(EndOfStreamContext& ec) +{ + mDBGOut.reset(); +} + +void TrackingStudySpec::finaliseCCDB(ConcreteDataMatcher& matcher, void* obj) +{ + if (o2::base::GRPGeomHelper::instance().finaliseCCDB(matcher, obj)) { + return; + } + if (matcher == ConcreteDataMatcher("IT3", "CLUSDICT", 0)) { + LOG(info) << "cluster dictionary updated"; + mITSDict = (const o2::its3::TopologyDictionary*)obj; + return; + } +} + +void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData) +{ + const auto& conf = ITS3TrackingStudyParam::Instance(); + if (conf.doDCA) { + doDCAStudy(recoData); + } + if (conf.doDCARefit) { + doDCARefitStudy(recoData); + } + if (mUseMC && conf.doPull) { + doPullStudy(recoData); + } + if (mUseMC && conf.doMC) { + doMCStudy(recoData); + } +} + +std::vector> TrackingStudySpec::prepareITSClusters(const o2::globaltracking::RecoContainer& data) const +{ + std::vector> itscl; + const auto& clusITS = data.getITSClusters(); + if (clusITS.size()) { + const auto& patterns = data.getITSClustersPatterns(); + itscl.reserve(clusITS.size()); + auto pattIt = patterns.begin(); + o2::its3::ioutils::convertCompactClusters(clusITS, pattIt, itscl, mITSDict); + } + return std::move(itscl); +} + +bool TrackingStudySpec::selectTrack(GTrackID trkID, o2::globaltracking::RecoContainer& recoData, bool checkMCTruth) const +{ + const auto& conf = ITS3TrackingStudyParam::Instance(); + if (!trkID.includesDet(GTrackID::ITS)) { + return false; + } + if (!recoData.isTrackSourceLoaded(trkID.getSource())) { + return false; + } + auto contributorsGID = recoData.getSingleDetectorRefs(trkID); + if (!contributorsGID[GTrackID::ITS].isIndexSet()) { // we need of course ITS + return false; + } + // ITS specific + const auto& itsTrk = recoData.getITSTrack(contributorsGID[GTrackID::ITS]); + if (itsTrk.getChi2() > conf.maxChi2 || itsTrk.getNClusters() < conf.minITSCls) { + return false; + } + // TPC specific + if (contributorsGID[GTrackID::TPC].isIndexSet()) { + const auto& tpcTrk = recoData.getTPCTrack(contributorsGID[GTrackID::TPC]); + if (tpcTrk.getNClusters() < conf.minTPCCls) { + return false; + } + } + // general + const auto& gTrk = recoData.getTrackParam(trkID); + if (gTrk.getPt() < conf.minPt || gTrk.getPt() > conf.maxPt) { + return false; + } + if (std::abs(gTrk.getEta()) > conf.maxEta) { + return false; + } + if (mUseMC && checkMCTruth) { + const auto& itsLbl = recoData.getTrackMCLabel(contributorsGID[GTrackID::ITS]); + if (!itsLbl.isValid()) { + return false; + } + if (contributorsGID[GTrackID::TPC].isIndexSet()) { + const auto& tpcLbl = recoData.getTrackMCLabel(contributorsGID[GTrackID::TPC]); + if (itsLbl != tpcLbl) { + return false; + } + } + if (contributorsGID[GTrackID::TRD].isIndexSet()) { + // TODO + } + if (contributorsGID[GTrackID::TOF].isIndexSet()) { + const auto& tofLbls = recoData.getTOFClustersMCLabels()->getLabels(contributorsGID[GTrackID::TOF]); + for (const auto& lbl : tofLbls) { + if (lbl.isValid()) { + return true; + } + } + } + } + return true; +} + +T2VMap TrackingStudySpec::buildT2V(o2::globaltracking::RecoContainer& recoData, bool includeCont, bool requireMCMatch) const +{ + // build track->vertex assoc., maybe including contributor tracks + const auto& conf = ITS3TrackingStudyParam::Instance(); + auto pvvec = recoData.getPrimaryVertices(); + auto trackIndex = recoData.getPrimaryVertexMatchedTracks(); // Global ID's for associated tracks + auto vtxRefs = recoData.getPrimaryVertexMatchedTrackRefs(); // references from vertex to these track IDs + auto nv = vtxRefs.size() - 1; // last entry is for unassigned tracks, ignore them + T2VMap t2v; + for (size_t iv = 0; iv < nv; ++iv) { + const auto& pv = pvvec[iv]; + if (pv.getNContributors() - 1 < conf.minPVCont) { + continue; + } + if (requireMCMatch) { + auto pvl = recoData.getPrimaryVertexMCLabel(iv); + } + const auto& vtxRef = vtxRefs[iv]; + int it = vtxRef.getFirstEntry(), itLim = it + vtxRef.getEntries(); + for (; it < itLim; it++) { + const auto& tvid = trackIndex[it]; + if (tvid.isAmbiguous()) { + continue; + } + if (!recoData.isTrackSourceLoaded(tvid.getSource())) { + continue; + } + if (mUseMC && requireMCMatch) { + const auto& pvlbl = recoData.getPrimaryVertexMCLabel(iv); + if (pvlbl.getEventID() != recoData.getTrackMCLabel(tvid).getEventID()) { + continue; + } + } + t2v[tvid] = iv; + if (includeCont) { + auto contributorsGID = recoData.getSingleDetectorRefs(tvid); + for (int cis = 0; cis < GTrackID::NSources; cis++) { + const auto cdm = GTrackID::getSourceDetectorsMask(cis); + if (!recoData.isTrackSourceLoaded(cis) || !cdm[DetID::ITS] || !contributorsGID[cis].isIndexSet()) { + continue; + } + if (mUseMC && requireMCMatch) { + const auto& pvlbl = recoData.getPrimaryVertexMCLabel(iv); + if (pvlbl.getEventID() != recoData.getTrackMCLabel(contributorsGID[cis]).getEventID()) { + continue; + } + } + t2v[contributorsGID[cis]] = iv; + } + } + } + } + return std::move(t2v); +} + +bool TrackingStudySpec::refitITSPVTrack(o2::globaltracking::RecoContainer& recoData, o2::track::TrackParCov& trFit, GTrackID gidx) +{ + if (gidx.getSource() != GTrackID::ITS) { + return false; + } + static auto pvvec = recoData.getPrimaryVertices(); + static auto t2v = buildT2V(recoData, true, true); + static const auto itsClusters = prepareITSClusters(recoData); + static std::vector itsTracksROF; + if (static bool done{false}; !done) { + done = true; + const auto& itsTracksROFRec = recoData.getITSTracksROFRecords(); + itsTracksROF.resize(recoData.getITSTracks().size()); + for (unsigned irf = 0, cnt = 0; irf < itsTracksROFRec.size(); irf++) { + int ntr = itsTracksROFRec[irf].getNEntries(); + for (int itr = 0; itr < ntr; itr++) { + itsTracksROF[cnt++] = irf; + } + } + } + auto prop = o2::base::Propagator::Instance(); + const auto& conf = ITS3TrackingStudyParam::Instance(); + std::array, 8> clArr{}; + std::array clAlpha{}; + const auto trkIn = recoData.getTrackParam(gidx); + const auto trkOut = recoData.getTrackParamOut(gidx); + const auto& itsTrOrig = recoData.getITSTrack(gidx); + int ncl = itsTrOrig.getNumberOfClusters(), rof = itsTracksROF[gidx.getIndex()]; + const auto& itsTrackClusRefs = recoData.getITSTracksClusterRefs(); + int clEntry = itsTrOrig.getFirstClusterEntry(); + const auto propagator = o2::base::Propagator::Instance(); + // convert PV to a fake cluster in the track DCA frame + const auto& pv = pvvec[t2v[gidx]]; + auto trkPV = trkIn; + if (!prop->propagateToDCA(pv, trkPV, prop->getNominalBz(), 2.0, conf.CorrType)) { + mTrackCounter -= gidx.getSource(); + return false; + } + // create base cluster from the PV, with the alpha corresponding to the track at DCA + float cosAlp = NAN, sinAlp = NAN; + o2::math_utils::sincos(trkPV.getAlpha(), sinAlp, cosAlp); + // vertex position rotated to track frame + clArr[0].setXYZ(pv.getX() * cosAlp + pv.getY() * sinAlp, -pv.getX() * sinAlp + pv.getY() * cosAlp, pv.getZ()); + clArr[0].setSigmaY2(0.5 * (pv.getSigmaX2() + pv.getSigmaY2())); + clArr[0].setSigmaZ2(pv.getSigmaZ2()); + clAlpha[0] = trkPV.getAlpha(); + for (int icl = 0; icl < ncl; ++icl) { // ITS clusters are referred in layer decreasing order + clArr[ncl - icl] = itsClusters[itsTrackClusRefs[clEntry + icl]]; + clAlpha[ncl - icl] = o2::its::GeometryTGeo::Instance()->getSensorRefAlpha(clArr[ncl - icl].getSensorID()); + } + // start refit + trFit = trkOut; + trFit.resetCovariance(1'000); + float chi2{0}; + for (int icl = ncl; icl >= 0; --icl) { // go backwards + if (!trFit.rotate(clAlpha[icl]) || !prop->propagateToX(trFit, clArr[icl].getX(), prop->getNominalBz(), 0.85, 2.0, conf.CorrType)) { + mTrackCounter -= gidx.getSource(); + return false; + } + chi2 += trFit.getPredictedChi2(clArr[icl]); + if (!trFit.update(clArr[icl])) { + mTrackCounter -= gidx.getSource(); + return false; + } + } + // chi2 < conf.maxChi2; should I cut here? + return true; +}; + +void TrackingStudySpec::doDCAStudy(o2::globaltracking::RecoContainer& recoData) +{ + /// analyse DCA of impact parameter for different track types + LOGP(info, "Doing DCA study"); + mTrackCounter.reset(); + const auto& conf = ITS3TrackingStudyParam::Instance(); + auto prop = o2::base::Propagator::Instance(); + TStopwatch sw; + sw.Start(); + int nDCAFits{0}, nDCAFitsFail{0}; + auto pvvec = recoData.getPrimaryVertices(); + auto trackIndex = recoData.getPrimaryVertexMatchedTracks(); // Global ID's for associated tracks + auto vtxRefs = recoData.getPrimaryVertexMatchedTrackRefs(); // references from vertex to these track IDs + auto nv = vtxRefs.size() - 1; // last entry is for unassigned tracks, ignore them + auto& stream = (*mDBGOut) << "dca"; + for (int iv = 0; iv < nv; iv++) { + const auto& pv = pvvec[iv]; + const auto& vtref = vtxRefs[iv]; + for (int is = 0; is < GTrackID::NSources; is++) { + const auto dm = GTrackID::getSourceDetectorsMask(is); + if (!recoData.isTrackSourceLoaded(is) || !dm[DetID::ITS]) { + mTrackCounter &= is; + continue; + } + int idMin = vtref.getFirstEntryOfSource(is), idMax = idMin + vtref.getEntriesOfSource(is); + for (int i = idMin; i < idMax; i++) { + const auto vid = trackIndex[i]; + if (!vid.isPVContributor()) { + mTrackCounter &= vid.getSource(); + continue; + } + + // we fit each different sub-track type, that include ITS, e.g. + // ITS,ITS-TPC,ITS-TPC-TRD,ITS-TPC-TOF,ITS-TPC-TRD-TOF + auto contributorsGID = recoData.getSingleDetectorRefs(vid); + for (int cis = 0; cis < GTrackID::NSources && cis <= is; cis++) { + const auto cdm = GTrackID::getSourceDetectorsMask(cis); + if (!recoData.isTrackSourceLoaded(cis) || !cdm[DetID::ITS] || !contributorsGID[cis].isIndexSet()) { + mTrackCounter &= cis; + continue; + } + if (!selectTrack(contributorsGID[cis], recoData)) { + mTrackCounter &= vid.getSource(); + continue; + } + + o2::dataformats::DCA dcaInfo; + const auto& trk = recoData.getTrackParam(contributorsGID[cis]); + auto trkRefit = trk; + // for ITS standalone tracks instead of having the trk at the pv we refit with the pv + if (conf.refitITS && cis == GTrackID::ITS && !refitITSPVTrack(recoData, trkRefit, contributorsGID[cis])) { + mTrackCounter -= cis; + continue; + } else { + trkRefit.invalidate(); + }; + + auto trkDCA = trk; + if (!prop->propagateToDCABxByBz(pv, trkDCA, 2.f, conf.CorrType, &dcaInfo)) { + mTrackCounter -= cis; + ++nDCAFitsFail; + continue; + } + + stream << "src=" << cis + << "pv=" << pv + << "trk=" << trk + << "trkRefit=" << trkRefit + << "trkAtPV=" << trkDCA + << "dca=" << dcaInfo; + + if (mUseMC) { + const auto& lbl = recoData.getTrackMCLabel(contributorsGID[cis]); + lbl.print(); + o2::dataformats::DCA dcaInfoMC; + const auto& eve = mcReader.getMCEventHeader(lbl.getSourceID(), lbl.getEventID()); + o2::dataformats::VertexBase mcEve; + mcEve.setPos({(float)eve.GetX(), (float)eve.GetY(), (float)eve.GetZ()}); + auto trkC = trk; + if (!prop->propagateToDCABxByBz(mcEve, trkC, 2.f, conf.CorrType, &dcaInfoMC)) { + mTrackCounter -= cis; + ++nDCAFitsFail; + continue; + } + const auto& mcTrk = mcReader.getTrack(lbl); + if (mcTrk == nullptr) { + LOGP(fatal, "mcTrk is null did selection fail?"); + } + stream << "mcTrk=" << *mcTrk + << "dca2MC=" << dcaInfoMC + << "lbl=" << lbl; + } + stream << "\n"; + + ++nDCAFits; + mTrackCounter += cis; + } + } + } + } + sw.Stop(); + LOGP(info, "doDCAStudy: accepted {} fits, failed {} (in {:.2f} seconds)", nDCAFits, nDCAFitsFail, sw.RealTime()); + mTrackCounter.print(); +} + +void TrackingStudySpec::doDCARefitStudy(o2::globaltracking::RecoContainer& recoData) +{ + /// analyse DCA of impact parameter for different track types while refitting the PV without the cand track + LOGP(info, "Doing DCARefit study"); + mTrackCounter.reset(); + const auto& conf = ITS3TrackingStudyParam::Instance(); + auto prop = o2::base::Propagator::Instance(); + TStopwatch sw; + sw.Start(); + + // build track->vertex assoc. + auto pvvec = recoData.getPrimaryVertices(); + auto vtxRefs = recoData.getPrimaryVertexMatchedTrackRefs(); // references from vertex to these track IDs + auto nv = vtxRefs.size() - 1; // last entry is for unassigned tracks, ignore them + auto t2v = buildT2V(recoData); + std::vector> v2t; + v2t.resize(nv); + auto creator = [&](const auto& trk, GTrackID trkID, float _t0, float terr) -> bool { + if constexpr (!isBarrelTrack()) { + mTrackCounter &= trkID.getSource(); + return false; + } + if (!trkID.includesDet(GTrackID::ITS)) { + mTrackCounter &= trkID.getSource(); + return false; + } + // general + if constexpr (isBarrelTrack()) { + if (trk.getPt() < conf.minPt || trk.getPt() > conf.maxPt) { + mTrackCounter &= trkID.getSource(); + return false; + } + if (std::abs(trk.getEta()) > conf.maxEta) { + mTrackCounter &= trkID.getSource(); + return false; + } + if (!t2v.contains(trkID)) { + mTrackCounter &= trkID.getSource(); + return false; + } + if (!selectTrack(trkID, recoData, mUseMC)) { + mTrackCounter &= trkID.getSource(); + return false; + } + } + v2t[t2v[trkID]].push_back(trkID); + return true; + }; + recoData.createTracksVariadic(creator); + + int nDCAFits{0}, nDCAFitsFail{0}; + auto& stream = (*mDBGOut) << "dcaRefit"; + for (size_t iv = 0; iv < nv; ++iv) { + const auto& pv = pvvec[iv]; + const auto& trkIDs = v2t[iv]; + if (trkIDs.size() - 1 < conf.minPVCont) { + continue; + } + std::vector trks; + trks.reserve(trkIDs.size()); + for (const auto& trkID : trkIDs) { + trks.push_back(recoData.getTrackParam(trkID)); + } + + if (!mVertexer.prepareVertexRefit(trks, pv)) { + continue; + } + std::vector trkMask(trkIDs.size(), true); + for (size_t it{0}; it < trkMask.size(); ++it) { + trkMask[it] = false; // mask current track from pv refit + if (it != 0) { + trkMask[it - 1] = true; // unmask previoustrack from pv refit + } + auto pvRefit = mVertexer.refitVertex(trkMask, pv); + if (pvRefit.getChi2() < 0) { + trkMask[it] = true; + continue; + } + + // check DCA both for refitted and original PV + o2::dataformats::DCA dcaInfo; + auto trkC = trks[it]; + if (!prop->propagateToDCABxByBz(pv, trkC, 2.f, conf.CorrType, &dcaInfo)) { + mTrackCounter -= trkIDs[it].getSource(); + ++nDCAFitsFail; + continue; + } + o2::dataformats::DCA dcaInfoRefit; + auto trkCRefit = trks[it]; + if (!prop->propagateToDCABxByBz(pv, trkCRefit, 2.f, conf.CorrType, &dcaInfoRefit)) { + mTrackCounter -= trkIDs[it].getSource(); + ++nDCAFitsFail; + continue; + } + + stream << "src=" << trkIDs[it].getSource() + << "pv=" << pv + << "trkAtPV=" << trkC + << "dca=" << dcaInfo + << "pvRefit=" << pvRefit + << "trkAtPVRefit=" << trkC + << "dcaRefit=" << dcaInfoRefit; + if (mUseMC) { + const auto& mcTrk = mcReader.getTrack(recoData.getTrackMCLabel(trkIDs[it])); + if (mcTrk == nullptr) { + LOGP(fatal, "mcTrk is null did selection fail?"); + } + stream << "mcTrk=" << *mcTrk; + } + stream << "\n"; + ++nDCAFits; + mTrackCounter += trkIDs[it].getSource(); + } + } + sw.Stop(); + LOGP(info, "doDCARefitStudy: accepted {} fits, failed {} (in {:.2f} seconds)", nDCAFits, nDCAFitsFail, sw.RealTime()); + mTrackCounter.print(); +} + +void TrackingStudySpec::doPullStudy(o2::globaltracking::RecoContainer& recoData) +{ + // check track pulls compared to mc generation + LOGP(info, "Doing Pull study"); + mTrackCounter.reset(); + TStopwatch sw; + sw.Start(); + int nPulls{0}, nPullsFail{0}; + auto prop = o2::base::Propagator::Instance(); + const auto& conf = ITS3TrackingStudyParam::Instance(); + + auto checkInTrack = [&](GTrackID trkID) { + if (!selectTrack(trkID, recoData)) { + mTrackCounter &= trkID.getSource(); + return; + } + const auto mcTrk = mcReader.getTrack(recoData.getTrackMCLabel(trkID)); + if (!mcTrk) { + return; + } + auto trk = recoData.getTrackParam(trkID); + + // for ITS standalone tracks we add the PV as an additional measurement point + if (conf.refitITS && trkID.getSource() == GTrackID::ITS && !refitITSPVTrack(recoData, trk, trkID)) { + mTrackCounter -= trkID.getSource(); + ++nPullsFail; + return; + } + + std::array xyz{(float)mcTrk->GetStartVertexCoordinatesX(), (float)mcTrk->GetStartVertexCoordinatesY(), (float)mcTrk->GetStartVertexCoordinatesZ()}, + pxyz{(float)mcTrk->GetStartVertexMomentumX(), (float)mcTrk->GetStartVertexMomentumY(), (float)mcTrk->GetStartVertexMomentumZ()}; + TParticlePDG* pPDG = TDatabasePDG::Instance()->GetParticle(mcTrk->GetPdgCode()); + if (!pPDG) { + mTrackCounter -= trkID.getSource(); + ++nPullsFail; + return; + } + o2::track::TrackPar mcTrkO2(xyz, pxyz, TMath::Nint(pPDG->Charge() / 3), false); + // propagate it to the alpha/X of the reconstructed track + if (!mcTrkO2.rotate(trk.getAlpha()) || !prop->PropagateToXBxByBz(mcTrkO2, trk.getX())) { + mTrackCounter -= trkID.getSource(); + ++nPullsFail; + return; + } + const auto contTrk = recoData.getSingleDetectorRefs(trkID); + const auto& itsTrk = recoData.getITSTrack(contTrk[GTrackID::ITS]); + + (*mDBGOut) + << "pull" + << "src=" << trkID.getSource() + << "itsTrk=" << itsTrk + << "mcTrk=" << mcTrkO2 + << "mcPart=" << mcTrk + << "trk=" << trk + << "\n"; + ++nPulls; + mTrackCounter += trkID.getSource(); + }; + + for (size_t iTrk{0}; iTrk < recoData.getITSTracks().size(); ++iTrk) { + checkInTrack(GTrackID(iTrk, GTrackID::ITS)); + } + for (size_t iTrk{0}; iTrk < recoData.getTPCITSTracks().size(); ++iTrk) { + checkInTrack(GTrackID(iTrk, GTrackID::ITSTPC)); + } + for (size_t iTrk{0}; iTrk < recoData.getITSTPCTRDTracksMCLabels().size(); ++iTrk) { + checkInTrack(GTrackID(iTrk, GTrackID::ITSTPCTRD)); + } + for (size_t iTrk{0}; iTrk < recoData.getITSTPCTOFMatches().size(); ++iTrk) { + checkInTrack(GTrackID(iTrk, GTrackID::ITSTPCTOF)); + } + for (size_t iTrk{0}; iTrk < recoData.getITSTPCTRDTOFMatches().size(); ++iTrk) { + checkInTrack(GTrackID(iTrk, GTrackID::ITSTPCTRDTOF)); + } + sw.Stop(); + LOGP(info, "doPullStudy: accepted {} pulls; rejected {} (in {:.2f} seconds)", nPulls, nPullsFail, sw.RealTime()); + mTrackCounter.print(); +} + +void TrackingStudySpec::doMCStudy(o2::globaltracking::RecoContainer& recoData) +{ + LOGP(info, "Doing MC study"); + mTrackCounter.reset(); + TStopwatch sw; + sw.Start(); + int nTracks{0}; + + const int iSrc{0}; + const int nev = mcReader.getNEvents(iSrc); + std::unordered_map info; + + LOGP(info, "** Filling particle table ... "); + for (int iEve{0}; iEve < nev; ++iEve) { + const auto& mcTrks = mcReader.getTracks(iSrc, iEve); + for (int iTrk{0}; iTrk < mcTrks.size(); ++iTrk) { + const auto& mcTrk = mcTrks[iTrk]; + const auto pdg = mcTrk.GetPdgCode(); + if (o2::O2DatabasePDG::Instance()->GetParticle(pdg) == nullptr) { + continue; + } + const auto apdg = std::abs(pdg); + if (apdg != 11 && apdg != 211 && apdg != 321 && apdg != 2212) { + continue; + } + o2::MCCompLabel lbl(iTrk, iEve, iSrc); + auto& part = info[lbl]; + part.mcTrack = mcTrk; + } + } + LOGP(info, "** Creating particle/clusters correspondence ... "); + const auto& clusters = recoData.getITSClusters(); + const auto& clustersMCLCont = recoData.getITSClustersMCLabels(); + for (auto iCluster{0}; iCluster < clusters.size(); ++iCluster) { + auto labs = clustersMCLCont->getLabels(iCluster); + for (auto& lab : labs) { + if (!lab.isValid() || lab.getSourceID() != 0 || !lab.isCorrect()) { + continue; + } + int trackID = 0, evID = 0, srcID = 0; + bool fake = false; + lab.get(trackID, evID, srcID, fake); + auto& cluster = clusters[iCluster]; + auto layer = o2::its::GeometryTGeo::Instance()->getLayer(cluster.getSensorID()); + auto& part = info[{trackID, evID, srcID}]; + part.clusters |= (1 << layer); + if (fake) { + part.fakeClusters |= (1 << layer); + } + } + } + LOGP(info, "** Analysing tracks ... "); + auto accountLbl = [&](const globaltracking::RecoContainer::GlobalIDSet& contributorsGID, DetID::ID det) { + if (contributorsGID[det].isIndexSet()) { + const auto& lbl = recoData.getTrackMCLabel(contributorsGID[det]); + if (lbl.isValid()) { + o2::MCCompLabel iLbl(lbl.getTrackID(), lbl.getEventID(), lbl.getSourceID()); + if (info.contains(iLbl)) { + auto& part = info[iLbl]; + SETBIT(part.recoTracks, det); + if (lbl.isFake()) { + SETBIT(part.fakeTracks, det); + } + } + } + } + }; + auto creator = [&](const auto& trk, GTrackID trkID, float _t0, float terr) -> bool { + if constexpr (!isBarrelTrack()) { + return false; + } + if (!trkID.includesDet(GTrackID::ITS)) { + return false; + } + // general + auto contributorsGID = recoData.getSingleDetectorRefs(trkID); + if (!contributorsGID[GTrackID::ITS].isIndexSet()) { // we need of course ITS + return false; + } + const auto& gLbl = recoData.getTrackMCLabel(trkID); + if (!gLbl.isValid()) { + return false; + } + o2::MCCompLabel iLbl(gLbl.getTrackID(), gLbl.getEventID(), gLbl.getSourceID()); + if (!info.contains(iLbl)) { + return false; + } + auto& part = info[iLbl]; + part.recoTrack = recoData.getTrackParam(trkID); + + accountLbl(contributorsGID, DetID::ITS); + accountLbl(contributorsGID, DetID::TPC); + accountLbl(contributorsGID, DetID::TRD); + accountLbl(contributorsGID, DetID::TOF); + + ++nTracks; + return true; + }; + recoData.createTracksVariadic(creator); + + LOGP(info, "Streaming output to tree"); + for (const auto& [_, part] : info) { + (*mDBGOut) << "mc" + << "part=" << part + << "\n"; + } + + sw.Stop(); + LOGP(info, "doMCStudy: accounted {} MCParticles and {} tracks (in {:.2f} seconds)", info.size(), nTracks, sw.RealTime()); +} + +DataProcessorSpec getTrackingStudySpec(GTrackID::mask_t srcTracks, GTrackID::mask_t srcClusters, bool useMC) +{ + std::vector outputs; + auto dataRequest = std::make_shared(); + + dataRequest->requestTracks(srcTracks, useMC); + dataRequest->requestIT3Clusters(useMC); + dataRequest->requestClusters(srcClusters, useMC); + dataRequest->requestPrimaryVertices(useMC); + auto ggRequest = std::make_shared(false, // orbitResetTime + true, // GRPECS=true + true, // GRPLHCIF + true, // GRPMagField + true, // askMatLUT + o2::base::GRPGeomRequest::Aligned, // geometry + dataRequest->inputs, + true); + + return DataProcessorSpec{ + .name = "its3-track-study", + .inputs = dataRequest->inputs, + .outputs = outputs, + .algorithm = AlgorithmSpec{adaptFromTask(dataRequest, ggRequest, srcTracks, useMC)}, + .options = {}}; +} + +} // namespace o2::its3::study diff --git a/Detectors/Upgrades/ITS3/study/src/its3-tracking-study-workflow.cxx b/Detectors/Upgrades/ITS3/study/src/its3-tracking-study-workflow.cxx new file mode 100644 index 0000000000000..e0a0aea1c368a --- /dev/null +++ b/Detectors/Upgrades/ITS3/study/src/its3-tracking-study-workflow.cxx @@ -0,0 +1,73 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "ITS3TrackingStudy/TrackingStudy.h" +#include "ReconstructionDataFormats/GlobalTrackID.h" +#include "DetectorsCommonDataFormats/DetID.h" +#include "CommonUtils/ConfigurableParam.h" +#include "Framework/CompletionPolicy.h" +#include "Framework/ConfigParamSpec.h" +#include "Framework/CompletionPolicyHelpers.h" +#include "Framework/CallbacksPolicy.h" +#include "DetectorsBase/DPLWorkflowUtils.h" +#include "GlobalTrackingWorkflowHelpers/InputHelper.h" +#include "DetectorsRaw/HBFUtilsInitializer.h" + +using namespace o2::framework; +using GID = o2::dataformats::GlobalTrackID; +using DetID = o2::detectors::DetID; + +// ------------------------------------------------------------------ +void customize(std::vector& policies) +{ + o2::raw::HBFUtilsInitializer::addNewTimeSliceCallback(policies); +} + +// we need to add workflow options before including Framework/runDataProcessing +void customize(std::vector& workflowOptions) +{ + // option allowing to set parameters + std::vector options{ + {"disable-mc", o2::framework::VariantType::Bool, false, {"disable MC propagation"}}, + {"track-sources", VariantType::String, std::string{GID::ALL}, {"comma-separated list of track sources to use"}}, + {"cluster-sources", VariantType::String, "ITS,TRD,TOF", {"comma-separated list of cluster sources to use"}}, + {"disable-root-input", VariantType::Bool, false, {"disable root-files input reader"}}, + {"configKeyValues", VariantType::String, "", {"Semicolon separated key=value strings ..."}}}; + o2::raw::HBFUtilsInitializer::addConfigOption(options); + std::swap(workflowOptions, options); +} + +// ------------------------------------------------------------------ + +#include "Framework/runDataProcessing.h" + +WorkflowSpec defineDataProcessing(ConfigContext const& configcontext) +{ + WorkflowSpec specs; + + GID::mask_t allowedSourcesTrc = GID::getSourcesMask("ITS,TPC,TRD,TOF,ITS-TPC,ITS-TPC-TRD,ITS-TPC-TOF,ITS-TPC-TRD-TOF"); + GID::mask_t allowedSourcesClus = GID::getSourcesMask("ITS,TPC,TRD,TOF"); + + o2::conf::ConfigurableParam::updateFromString(configcontext.options().get("configKeyValues")); + auto useMC = !configcontext.options().get("disable-mc"); + + GID::mask_t srcTrc = allowedSourcesTrc & GID::getSourcesMask(configcontext.options().get("track-sources")); + GID::mask_t srcCls = allowedSourcesClus & GID::getSourcesMask(configcontext.options().get("cluster-sources")); + + o2::globaltracking::InputHelper::addInputSpecs(configcontext, specs, srcCls, srcTrc, srcTrc, useMC); + o2::globaltracking::InputHelper::addInputSpecsPVertex(configcontext, specs, useMC); + + specs.emplace_back(o2::its3::study::getTrackingStudySpec(srcTrc, srcCls, useMC)); + + o2::raw::HBFUtilsInitializer hbfIni(configcontext, specs); + + return std::move(specs); +}