From 3914fd1c136b362f5154f26a1c535893b4f76bc2 Mon Sep 17 00:00:00 2001 From: FabiolaLP Date: Fri, 26 Sep 2025 20:10:41 -0600 Subject: [PATCH 1/4] Nuspex: add antineutron CEX table producer & datamodel; update Nuspex CMake; clang-format --- PWGLF/DataModel/LFAntinCexTables.h | 130 +++ PWGLF/TableProducer/Nuspex/CMakeLists.txt | 5 + .../Nuspex/nucleiAntineutronCex.cxx | 815 ++++++++++++++++++ 3 files changed, 950 insertions(+) create mode 100644 PWGLF/DataModel/LFAntinCexTables.h create mode 100644 PWGLF/TableProducer/Nuspex/nucleiAntineutronCex.cxx diff --git a/PWGLF/DataModel/LFAntinCexTables.h b/PWGLF/DataModel/LFAntinCexTables.h new file mode 100644 index 00000000000..63a68dbc9a1 --- /dev/null +++ b/PWGLF/DataModel/LFAntinCexTables.h @@ -0,0 +1,130 @@ +// 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 LFAntinCexTables.h +/// \brief Slim tables for nucleiAntineutronCex +/// \author Fabiola Lugo +/// + +#ifndef PWGLF_DATAMODEL_LFANTINCEXTABLES_H_ +#define PWGLF_DATAMODEL_LFANTINCEXTABLES_H_ + +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" + +namespace o2::aod +{ +namespace antinCexNS +{ +// Metadata +DECLARE_SOA_COLUMN(IsCex, isCex, bool); // 1=CEX (from antin), 0=BG +DECLARE_SOA_COLUMN(MotherPdg, motherPdg, int32_t); // mother PDG +DECLARE_SOA_COLUMN(ColId, colId, int32_t); // mcCollisionId +DECLARE_SOA_COLUMN(PId, pId, int32_t); // proton MC id +DECLARE_SOA_COLUMN(AntipId, antipId, int32_t); // antiproton MC id + +// MC (pair) +DECLARE_SOA_COLUMN(McPairP, mcPairP, float); +DECLARE_SOA_COLUMN(McPairPt, mcPairPt, float); +DECLARE_SOA_COLUMN(McPairPz, mcPairPz, float); +DECLARE_SOA_COLUMN(McDplane, mcDplane, float); +DECLARE_SOA_COLUMN(McAngleDeg, mcAngleDeg, float); +DECLARE_SOA_COLUMN(McVtxX, mcVtxX, float); +DECLARE_SOA_COLUMN(McVtxY, mcVtxY, float); +DECLARE_SOA_COLUMN(McVtxZ, mcVtxZ, float); + +// Tracks (pair, fitter) +DECLARE_SOA_COLUMN(TrkPairP, trkPairP, float); +DECLARE_SOA_COLUMN(TrkPairPt, trkPairPt, float); +DECLARE_SOA_COLUMN(TrkPairPz, trkPairPz, float); +DECLARE_SOA_COLUMN(TrkAngleDeg, trkAngleDeg, float); +DECLARE_SOA_COLUMN(TrkVtxfitDcaPair, trkVtxfitDcaPair, float); +DECLARE_SOA_COLUMN(TrkVtxfitR, trkVtxfitR, float); +DECLARE_SOA_COLUMN(TrkVtxfitDistToPv, trkVtxfitDistToPv, float); +DECLARE_SOA_COLUMN(TrkVtxfitSecVtxX, trkVtxfitSecVtxX, float); +DECLARE_SOA_COLUMN(TrkVtxfitSecVtxY, trkVtxfitSecVtxY, float); +DECLARE_SOA_COLUMN(TrkVtxfitSecVtxZ, trkVtxfitSecVtxZ, float); + +// Fit quality (fit − MC) +DECLARE_SOA_COLUMN(VtxfitChi2, vtxfitChi2, float); +DECLARE_SOA_COLUMN(VtxfitStatus, vtxfitStatus, int32_t); +DECLARE_SOA_COLUMN(NCand, nCand, int32_t); +DECLARE_SOA_COLUMN(VtxfitDX, vtxfitDX, float); +DECLARE_SOA_COLUMN(VtxfitDY, vtxfitDY, float); +DECLARE_SOA_COLUMN(VtxfitDZ, vtxfitDZ, float); +DECLARE_SOA_COLUMN(VtxfitD3D, vtxfitD3D, float); + +// Proton track +DECLARE_SOA_COLUMN(PTrkP, pTrkP, float); +DECLARE_SOA_COLUMN(PTrkPx, pTrkPx, float); +DECLARE_SOA_COLUMN(PTrkPy, pTrkPy, float); +DECLARE_SOA_COLUMN(PTrkPz, pTrkPz, float); +DECLARE_SOA_COLUMN(PTrkEta, pTrkEta, float); +DECLARE_SOA_COLUMN(PTrkTpcSignal, pTrkTpcSignal, float); +DECLARE_SOA_COLUMN(PTrkNClsIts, pTrkNClsIts, int16_t); + +// Antiproton track +DECLARE_SOA_COLUMN(AntipTrkP, antipTrkP, float); +DECLARE_SOA_COLUMN(AntipTrkPx, antipTrkPx, float); +DECLARE_SOA_COLUMN(AntipTrkPy, antipTrkPy, float); +DECLARE_SOA_COLUMN(AntipTrkPz, antipTrkPz, float); +DECLARE_SOA_COLUMN(AntipTrkEta, antipTrkEta, float); +DECLARE_SOA_COLUMN(AntipTrkTpcSignal, antipTrkTpcSignal, float); +DECLARE_SOA_COLUMN(AntipTrkNClsIts, antipTrkNClsIts, int16_t); + +// Cuts Mask +DECLARE_SOA_COLUMN(SelMask, selMask, uint32_t); + +DECLARE_SOA_COLUMN(PairPointingAngleDeg, pairPointingAngleDeg, float); +DECLARE_SOA_COLUMN(PairPBalance, pairPBalance, float); +DECLARE_SOA_COLUMN(PairPtBalance, pairPtBalance, float); +DECLARE_SOA_COLUMN(PairQ, pairQ, float); + +DECLARE_SOA_COLUMN(DPairP, dPairP, float); +DECLARE_SOA_COLUMN(DPairPt, dPairPt, float); +DECLARE_SOA_COLUMN(DPairPz, dPairPz, float); +DECLARE_SOA_COLUMN(DOpenAngle, dOpenAngle, float); + +DECLARE_SOA_COLUMN(SVNearestLayerId, svNearestLayerId, int16_t); +DECLARE_SOA_COLUMN(SVDeltaRtoLayer, svDeltaRtoLayer, float); + +DECLARE_SOA_COLUMN(PTrkItsHitMap, pTrkItsHitMap, uint16_t); +DECLARE_SOA_COLUMN(APTrkItsHitMap, apTrkItsHitMap, uint16_t); +DECLARE_SOA_COLUMN(PLayersOK, pLayersOK, int8_t); +DECLARE_SOA_COLUMN(APLayersOK, apLayersOK, int8_t); + +DECLARE_SOA_COLUMN(PvtxZ, pvtxZ, float); +} // namespace antinCexNS + +// Table +DECLARE_SOA_TABLE(antinCexPairs, "AOD", "ANTINCEX", + antinCexNS::IsCex, + antinCexNS::MotherPdg, antinCexNS::ColId, antinCexNS::PId, antinCexNS::AntipId, + antinCexNS::McPairP, antinCexNS::McPairPt, antinCexNS::McPairPz, + antinCexNS::McDplane, antinCexNS::McAngleDeg, antinCexNS::McVtxX, antinCexNS::McVtxY, antinCexNS::McVtxZ, + antinCexNS::TrkPairP, antinCexNS::TrkPairPt, antinCexNS::TrkPairPz, antinCexNS::TrkAngleDeg, + antinCexNS::TrkVtxfitDcaPair, antinCexNS::TrkVtxfitR, antinCexNS::TrkVtxfitDistToPv, + antinCexNS::TrkVtxfitSecVtxX, antinCexNS::TrkVtxfitSecVtxY, antinCexNS::TrkVtxfitSecVtxZ, + antinCexNS::VtxfitChi2, antinCexNS::VtxfitStatus, antinCexNS::NCand, + antinCexNS::VtxfitDX, antinCexNS::VtxfitDY, antinCexNS::VtxfitDZ, antinCexNS::VtxfitD3D, + antinCexNS::PTrkP, antinCexNS::PTrkPx, antinCexNS::PTrkPy, antinCexNS::PTrkPz, antinCexNS::PTrkEta, antinCexNS::PTrkTpcSignal, antinCexNS::PTrkNClsIts, + antinCexNS::AntipTrkP, antinCexNS::AntipTrkPx, antinCexNS::AntipTrkPy, antinCexNS::AntipTrkPz, antinCexNS::AntipTrkEta, antinCexNS::AntipTrkTpcSignal, antinCexNS::AntipTrkNClsIts, + antinCexNS::SelMask, + antinCexNS::PairPointingAngleDeg, antinCexNS::PairPBalance, antinCexNS::PairPtBalance, antinCexNS::PairQ, + antinCexNS::DPairP, antinCexNS::DPairPt, antinCexNS::DPairPz, antinCexNS::DOpenAngle, + antinCexNS::SVNearestLayerId, antinCexNS::SVDeltaRtoLayer, + antinCexNS::PTrkItsHitMap, antinCexNS::APTrkItsHitMap, antinCexNS::PLayersOK, antinCexNS::APLayersOK, + antinCexNS::PvtxZ); + +} // namespace o2::aod + +#endif // PWGLF_DATAMODEL_LFANTINCEXTABLES_H_ diff --git a/PWGLF/TableProducer/Nuspex/CMakeLists.txt b/PWGLF/TableProducer/Nuspex/CMakeLists.txt index 98dac784da5..127b4d465de 100644 --- a/PWGLF/TableProducer/Nuspex/CMakeLists.txt +++ b/PWGLF/TableProducer/Nuspex/CMakeLists.txt @@ -103,3 +103,8 @@ o2physics_add_dpl_workflow(he3-lambda-analysis SOURCES he3LambdaAnalysis.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::EventFilteringUtils COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(nuclei-antineutron-cex + SOURCES nucleiAntineutronCex.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) diff --git a/PWGLF/TableProducer/Nuspex/nucleiAntineutronCex.cxx b/PWGLF/TableProducer/Nuspex/nucleiAntineutronCex.cxx new file mode 100644 index 00000000000..0c96de014ce --- /dev/null +++ b/PWGLF/TableProducer/Nuspex/nucleiAntineutronCex.cxx @@ -0,0 +1,815 @@ +// 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 nucleiAntineutronCex.cxx +/// \brief Analysis task for antineutron detection through cex interactions +/// \author Fabiola Lugo +/// + +#include "DCAFitter/DCAFitterN.h" +#include "DetectorsBase/Propagator.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/Logger.h" +#include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/TrackParametrization.h" + +#include "TMCProcess.h" + +#include +// ROOT +#include "CommonConstants/MathConstants.h" + +#include "TMath.h" +#include "TPDGCode.h" +#include "TVector3.h" + +#include +#include +// +#include "PWGLF/DataModel/LFAntinCexTables.h" + +using namespace o2; +using namespace o2::framework; +using o2::constants::math::Rad2Deg; + +struct NucleiAntineutronCex { + // Slicing per colision + Preslice perMcByColl = aod::mcparticle::mcCollisionId; + // Check available tables in the AOD, specifically TracksIU, TracksCovIU + using TracksWCovMc = soa::Join; + + // === Cut values === + static constexpr double kIts2MinR = 2.2; // ITS2 min radius [cm] + static constexpr double kIts2MaxR = 39.0; // ITS2 max radius [cm] + static constexpr double kIts2MaxVz = 39.0; // ITS2 max |vz| [cm] + static constexpr double kAccMaxEta = 1.2; // acceptance |eta| + static constexpr double kAccMaxVz = 5.3; // acceptance |vz| [cm] + static constexpr double kStrictEta = 0.9; // tighter eta cut + static constexpr double kInitDplane = 10.0; // init dplane + static constexpr double kHuge = 1e9; // fallback for bad denom + static constexpr int kMinItsHits = 2; + static constexpr double kVtxTol = 1e-4; + + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + Produces outPairs; + + void init(InitContext const&) + { + // Primary vertex + histos.add("hVx", "Primary vertex X;X (cm);Entries", kTH1F, {{100, -5., 5.}}); + histos.add("hVy", "Primary vertex Y;Y (cm);Entries", kTH1F, {{100, -5., 5.}}); + histos.add("hVz", "Primary vertex Z;Z (cm);Entries", kTH1F, {{200, -20., 20.}}); + + // Primary antineutrons + histos.add("antin_p", "Total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("antin_px", "p_{x};p_{x} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("antin_py", "p_{y};p_{y} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("antin_pz", "p_{z};p_{z} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("antin_eta", "Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); + histos.add("antin_p_ITScuts", "Momentum with ITS cuts;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + + // Primary neutrons + histos.add("n_p", "Total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("n_px", "p_{x};p_{x} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("n_py", "p_{y};p_{y} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("n_pz", "p_{z};p_{z} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("n_eta", "Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); + histos.add("n_p_ITScuts", "Momentum with ITS cuts;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + + // Primary antiprotons + histos.add("antipP", "Total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("antipPx", "p_{x};p_{x} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("antipPy", "p_{y};p_{y} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("antipPz", "p_{z};p_{z} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("antipEta", "Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); + histos.add("antipP_ITScuts", "Momentum with ITS cuts;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + + // Primary protons + histos.add("pP", "Total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("pPx", "p_{x};p_{x} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("pPy", "p_{y};p_{y} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("pPz", "p_{z};p_{z} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("pEta", "Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); + histos.add("pP_ITScuts", "Momentum with ITS cuts;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + + // test (MC) + histos.add("antip_test", "Secondary antiprotons;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + + // CEX pair from antineutron (MC) + histos.add("cexPairMcP", "CEX pair total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("cexPairMcPt", "CEX pair p_{T};p_{T} (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("cexPairMcPz", "CEX pair p_{z};p_{z} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("cex_pairmcDplane", "CEX pair d_{plane};d_{plane} (cm);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("cex_pairmc_angle", "Pair opening angle;Angle (°);Entries", kTH1F, {{180, 0., 180.}}); + histos.add("cex_pairmc_vtx", "MC CEX pair vertex;X (cm);Y (cm)", kTH2F, {{100, -50., 50.}, {100, -50., 50.}}); + histos.add("cex_pairmc_vtxz", "MC secondary vertex Z;Z (cm);Entries", kTH1F, {{200, -60., 60.}}); + histos.add("cexPairMcPITScuts", "CEX pair momentum (ITS cuts);|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + + // CEX pair normalized to antineutron (MC) + histos.add("cexn_pairmc_p", "Pair p / antineutron p;p/p_{#bar{n}};Entries", kTH1F, {{100, 0., 2.}}); + histos.add("cexn_pairmc_pt", "Pair p_{T} / antineutron p_{T};p_{T}/p_{T,#bar{n}};Entries", kTH1F, {{100, 0., 2.}}); + histos.add("cexn_pairmc_pz", "Pair p_{z} / antineutron p_{z};p_{z}/p_{z,#bar{n}};Entries", kTH1F, {{100, -2., 2.}}); + + // BG pair (not from antineutron) (MC) + histos.add("cexbg_pairmc_p", "Background pair momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("cexbg_pairmc_pt", "Background pair p_{T};p_{T} (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("cexbg_pairmc_pz", "Background pair p_{z};p_{z} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("cexbg_pairmcDplane", "Background d_{plane};d_{plane} (cm);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("cexbg_pairmc_angle", "Background opening angle;Angle (°);Entries", kTH1F, {{180, 0., 180.}}); + histos.add("cexbg_pairmc_vtx", "Background pair vertex;X (cm);Y (cm)", kTH2F, {{100, -50., 50.}, {100, -50., 50.}}); + histos.add("cexbg_pairmc_vtxz", "Background secondary vertex Z;Z (cm);Entries", kTH1F, {{200, -60., 60.}}); + histos.add("cexbg_pairmc_pITScuts", "Background momentum (ITS cuts);|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + + // CEX pair from antineutron (TRK) + histos.add("cex_pairtrk_angle", "Pair opening angle (tracks);Angle (°);Entries", kTH1F, {{180, 0., 180.}}); + histos.add("cexPairTrkP", "Pair momentum (tracks);|p| (GeV/c);Entries", kTH1F, {{120, 0., 12.}}); + histos.add("cexPairTrkPt", "Pair p_{T} (tracks);p_{T} (GeV/c);Entries", kTH1F, {{120, 0., 12.}}); + histos.add("cexPairTrkPz", "Pair p_{z} (tracks);p_{z} (GeV/c);Entries", kTH1F, {{120, -12., 12.}}); + histos.add("cex_pairtrkVtxfitDcaPair", "DCA between tracks at PCA;DCA (cm);Entries", kTH1F, {{200, 0., 10.}}); + histos.add("cex_pairtrkVtxfitR", "Secondary-vertex radius (PCA);R (cm);Entries", kTH1F, {{200, 0., 60.}}); + histos.add("cex_pairtrkVtxfitDistToPv", "Distance from secondary vertex to PV;dist (cm);Entries", kTH1F, {{240, 0., 120.}}); + histos.add("cex_pairtrk_vtxfit_secVtxXY", "Secondary vertex (PCA);X (cm);Y (cm)", kTH2F, {{200, -60., 60.}, {200, -60., 60.}}); + histos.add("cex_pairtrk_vtxfit_secVtxZ", "Secondary vertex Z (PCA);Z (cm);Entries", kTH1F, {{240, -60., 60.}}); + + // BG pair (not from antineutron) (TRK) + histos.add("cexbg_pairtrk_angle", "Background opening angle (tracks);Angle (°);Entries", kTH1F, {{180, 0., 180.}}); + histos.add("cexbg_pairtrk_p", "Pair momentum (tracks);|p| (GeV/c);Entries", kTH1F, {{120, 0., 12.}}); + histos.add("cexbg_pairtrk_pt", "Pair p_{T} (tracks);p_{T} (GeV/c);Entries", kTH1F, {{120, 0., 12.}}); + histos.add("cexbg_pairtrk_pz", "Pair p_{z} (tracks);p_{z} (GeV/c);Entries", kTH1F, {{120, -12., 12.}}); + histos.add("cexbg_pairtrkVtxfitDcaPair", "DCA between tracks at PCA;DCA (cm);Entries", kTH1F, {{200, 0., 10.}}); + histos.add("cexbg_pairtrkVtxfitR", "Secondary-vertex radius (PCA);R (cm);Entries", kTH1F, {{200, 0., 60.}}); + histos.add("cexbg_pairtrkVtxfitDistToPv", "Distance from secondary vertex to PV;dist (cm);Entries", kTH1F, {{240, 0., 120.}}); + histos.add("cexbg_pairtrk_vtxfit_secVtxXY", "Secondary vertex (PCA);X (cm);Y (cm)", kTH2F, {{200, -60., 60.}, {200, -60., 60.}}); + histos.add("cexbg_pairtrk_vtxfit_secVtxZ", "Secondary vertex Z (PCA);Z (cm);Entries", kTH1F, {{240, -60., 60.}}); + + // Vertex fit (DCAFitter2 / PCA) + histos.add("vtxfitChi2", "DCAFitter2 #chi^{2};#chi^{2};Entries", kTH1F, {{200, 0., 100.}}); + histos.add("vtxfitStatus", "Fit status (0=OK);code;Entries", kTH1I, {{10, 0., 10.}}); + histos.add("vtxfit_mc_dX", "SV residual X (fit - MC);#Delta X (cm);Entries", kTH1F, {{400, -20., 20.}}); + histos.add("vtxfit_mc_dY", "SV residual Y (fit - MC);#Delta Y (cm);Entries", kTH1F, {{400, -20., 20.}}); + histos.add("vtxfit_mc_dZ", "SV residual Z (fit - MC);#Delta Z (cm);Entries", kTH1F, {{400, -20., 20.}}); + histos.add("vtxfit_mc_d3D", "SV distance |fit - MC|;#Delta r (cm);Entries", kTH1F, {{300, 0., 30.}}); + } + + static o2::track::TrackParCov makeTPCovFromAOD(const TracksWCovMc::iterator& tr) + { + using o2::track::TrackParCov; + TrackParCov tpcov; + // Local state: x, alpha, y, z, snp, tgl, q/pt + float par[5] = {tr.y(), tr.z(), tr.snp(), tr.tgl(), tr.signed1Pt()}; + tpcov.set(tr.x(), tr.alpha(), par); + + // Covariance matrix (15 terms) in O2 order + std::array cov = { + tr.cYY(), tr.cZY(), tr.cZZ(), + tr.cSnpY(), tr.cSnpZ(), tr.cSnpSnp(), + tr.cTglY(), tr.cTglZ(), tr.cTglSnp(), + tr.cTglTgl(), tr.c1PtY(), tr.c1PtZ(), + tr.c1PtSnp(), tr.c1PtTgl(), tr.c1Pt21Pt2()}; + tpcov.setCov(cov); + return tpcov; + } + + void process(aod::McCollisions const& cols, aod::McParticles const& particles, TracksWCovMc const& tracks) + { + double pvtxX = 0; + double pvtxY = 0; + double pvtxZ = 0; + for (auto const& col : cols) { + const auto colId = col.globalIndex(); + auto mcPartsThis = particles.sliceBy(perMcByColl, colId); + + if (std::isfinite(col.posX()) && std::isfinite(col.posY()) && std::isfinite(col.posZ())) { + pvtxX = col.posX(); + pvtxY = col.posY(); + pvtxZ = col.posZ(); + histos.fill(HIST("hVx"), pvtxX); + histos.fill(HIST("hVy"), pvtxY); + histos.fill(HIST("hVz"), pvtxZ); + } + + for (const auto& particle : mcPartsThis) { + + // Primary antineutrons + if (particle.pdgCode() == -kNeutron && particle.isPhysicalPrimary()) { + histos.fill(HIST("antin_p"), particle.p()); + histos.fill(HIST("antin_px"), particle.px()); + histos.fill(HIST("antin_py"), particle.py()); + histos.fill(HIST("antin_pz"), particle.pz()); + histos.fill(HIST("antin_eta"), particle.eta()); + if (std::abs(particle.eta()) < kAccMaxEta && std::abs(particle.vz()) < kAccMaxVz) + histos.fill(HIST("antin_p_ITScuts"), particle.p()); + } + // Primary neutrons + if (particle.pdgCode() == kNeutron && particle.isPhysicalPrimary()) { + histos.fill(HIST("n_p"), particle.p()); + histos.fill(HIST("n_px"), particle.px()); + histos.fill(HIST("n_py"), particle.py()); + histos.fill(HIST("n_pz"), particle.pz()); + histos.fill(HIST("n_eta"), particle.eta()); + if (std::abs(particle.eta()) < kAccMaxEta && std::abs(particle.vz()) < kAccMaxVz) + histos.fill(HIST("n_p_ITScuts"), particle.p()); + } + // Primary antiprotons + if (particle.pdgCode() == -kProton && particle.isPhysicalPrimary()) { + histos.fill(HIST("antipP"), particle.p()); + histos.fill(HIST("antipPx"), particle.px()); + histos.fill(HIST("antipPy"), particle.py()); + histos.fill(HIST("antipPz"), particle.pz()); + histos.fill(HIST("antipEta"), particle.eta()); + if (std::abs(particle.eta()) < kAccMaxEta && std::abs(particle.vz()) < kAccMaxVz) + histos.fill(HIST("antipP_ITScuts"), particle.p()); + } + // Primary protons + if (particle.pdgCode() == kProton && particle.isPhysicalPrimary()) { + histos.fill(HIST("pP"), particle.p()); + histos.fill(HIST("pPx"), particle.px()); + histos.fill(HIST("pPy"), particle.py()); + histos.fill(HIST("pPz"), particle.pz()); + histos.fill(HIST("pEta"), particle.eta()); + if (std::abs(particle.eta()) < kAccMaxEta && std::abs(particle.vz()) < kAccMaxVz) + histos.fill(HIST("pP_ITScuts"), particle.p()); + } + + // Seconday antiprotons from material + const auto procEnum = particle.getProcess(); + const bool isSecondaryFromMaterial = (!particle.producedByGenerator()) && (procEnum == kPHadronic || procEnum == kPHInhelastic); + if (particle.pdgCode() != -kProton || !isSecondaryFromMaterial || particle.mothersIds().empty()) + continue; + histos.fill(HIST("antip_test"), particle.p()); + + // Primary mother + bool hasPrimaryMotherAntip = false; + double motherPt = 0.0; + double motherPz = 0.0; + double motherVz = 0.0; + double motherP = 0.0; + double motherEta = 0.0; + int motherPdg = 0; + + for (const auto& mother : particle.mothers_as()) { + if (mother.isPhysicalPrimary()) { + hasPrimaryMotherAntip = true; + motherPt = mother.pt(); + motherPz = mother.pz(); + motherVz = mother.vz(); + motherP = mother.p(); + motherEta = mother.eta(); + motherPdg = mother.pdgCode(); + break; + } + } + if (!hasPrimaryMotherAntip) + continue; + + double antipVx = particle.vx(); + double antipVy = particle.vy(); + double antipVz = particle.vz(); + double antipP = particle.p(); + double antipPx = particle.px(); + double antipPy = particle.py(); + double antipPz = particle.pz(); + double antipE = particle.e(); + double antipEta = particle.eta(); + int antipId = particle.globalIndex(); + + // Selection conditions: Produced in the ITS + const double R = std::sqrt(antipVx * antipVx + antipVy * antipVy); + // Config for ITS + // if(3.9<=R && R<=43.0 && std::abs(antipVz)<=48.9){ + // Config for ITS2 + if (R < kIts2MinR || R > kIts2MaxR || std::abs(antipVz) > kIts2MaxVz) + continue; + if (std::abs(motherEta) >= kAccMaxEta || std::abs(motherVz) >= kAccMaxVz) + continue; + + // Pion minus veto + bool pionMinus = false; + for (const auto& particle1 : mcPartsThis) { + const auto proc1Enum = particle1.getProcess(); + const bool isSecondaryFromMaterial1 = (!particle1.producedByGenerator()) && (proc1Enum == kPHadronic || proc1Enum == kPHInhelastic); + if (particle1.mcCollisionId() != colId) + continue; + if (particle1.pdgCode() != kPiMinus || !isSecondaryFromMaterial1 || particle1.mothersIds().empty()) + continue; + bool hasPrimaryMotherPim = false; + for (const auto& mother : particle1.mothers_as()) { + if (mother.isPhysicalPrimary()) { + hasPrimaryMotherPim = true; + break; + } + } + if (!hasPrimaryMotherPim) + continue; + double pimVx = particle1.vx(); + double pimVy = particle1.vy(); + double pimVz = particle1.vz(); + if (std::abs(pimVx - antipVx) < kVtxTol && std::abs(pimVy - antipVy) < kVtxTol && std::abs(pimVz - antipVz) < kVtxTol) { + pionMinus = true; + break; + } + } + + // Pion plus veto + bool pionPlus = false; + for (const auto& particle2 : mcPartsThis) { + if (particle2.mcCollisionId() != colId) + continue; + const auto proc2Enum = particle2.getProcess(); + const bool isSecondaryFromMaterial2 = (!particle2.producedByGenerator()) && (proc2Enum == kPHadronic || proc2Enum == kPHInhelastic); + if (particle2.pdgCode() != kPiPlus || !isSecondaryFromMaterial2 || particle2.mothersIds().empty()) + continue; + bool hasPrimaryMotherPip = false; + for (const auto& mother : particle2.mothers_as()) { + if (mother.isPhysicalPrimary()) { + hasPrimaryMotherPip = true; + break; + } + } + if (!hasPrimaryMotherPip) + continue; + double pipVx = particle2.vx(); + double pipVy = particle2.vy(); + double pipVz = particle2.vz(); + if (std::abs(pipVx - antipVx) < kVtxTol && std::abs(pipVy - antipVy) < kVtxTol && std::abs(pipVz - antipVz) < kVtxTol) { + pionPlus = true; + break; + } + } + + if (pionPlus || pionMinus) + continue; + + // CEX selection + double dplane = kInitDplane; + double dplaneTmp = 0; + double p = 0; + double pTmp = 0; + double pcexPx = 0; + double pcexPy = 0; + double pcexPz = 0; + double pcexP = 0; + double e = 0; + double eTmp = 0; + double pcexE = 0; + int k_plane = -1; + int k_e = -1; + int k_p = -1; + + // Secondary proton from material + for (const auto& particle3 : mcPartsThis) { + if (particle3.mcCollisionId() != colId) + continue; + const auto proc3Enum = particle3.getProcess(); + const bool isSecondaryFromMaterial3 = (!particle3.producedByGenerator()) && (proc3Enum == kPHadronic || proc3Enum == kPHInhelastic); + if (particle3.pdgCode() != kProton || !isSecondaryFromMaterial3 || particle3.mothersIds().empty()) + continue; + bool hasPrimaryMotherP = false; + for (const auto& mother : particle3.mothers_as()) { + if (mother.isPhysicalPrimary()) { + hasPrimaryMotherP = true; + break; + } + } + if (!hasPrimaryMotherP) + continue; + double pVx = particle3.vx(); + double pVy = particle3.vy(); + double pVz = particle3.vz(); + double pP = particle3.p(); + double pPx = particle3.px(); + double pPy = particle3.py(); + double pPz = particle3.pz(); + double pE = particle3.e(); + double pEta = particle3.eta(); + if (std::abs(pVx - antipVx) < kVtxTol && std::abs(pVy - antipVy) < kVtxTol && std::abs(pVz - antipVz) < kVtxTol) { + // Same mother + bool shareMother = false; + const auto& momsAp = particle.mothersIds(); // antiproton + const auto& momsP = particle3.mothersIds(); // proton + for (const auto& ida : momsAp) { + for (const auto& idp : momsP) { + if (ida == idp) { + shareMother = true; + break; + } + } + if (shareMother) + break; + } + if (!shareMother) + continue; + + // CEX proton selection + // dplaneTmp = (pPy*antipPz - pPz*antipPy)*(pvtxX-antipVx) + (pPz*antipPx - pPx*antipPz)*(pvtxY-antipVy) + (pPx*antipPy - pPy*antipPx)*(pvtxZ-antipVz); + double nx = (pPy * antipPz - pPz * antipPy); + double ny = (pPz * antipPx - pPx * antipPz); + double nz = (pPx * antipPy - pPy * antipPx); + double rx = (pvtxX - antipVx); + double ry = (pvtxY - antipVy); + double rz = (pvtxZ - antipVz); + double denom = nx * nx + ny * ny + nz * nz; + if (denom > 0.) { + dplaneTmp = std::abs(nx * rx + ny * ry + nz * rz) / std::sqrt(denom); + } else { + dplaneTmp = kHuge; + } + if (std::abs(dplaneTmp) < std::abs(dplane)) { + k_plane = particle3.globalIndex(); + dplane = dplaneTmp; + } + + eTmp = antipE + pE; + if (std::abs(eTmp) > std::abs(e)) { + k_e = particle3.globalIndex(); + e = eTmp; + pcexE = pE; + } + + pTmp = std::sqrt(std::pow((pPx + antipPx), 2) + std::pow((pPy + antipPy), 2) + std::pow((pPz + antipPz), 2)); + if (std::abs(pTmp) > std::abs(p)) { + k_p = particle3.globalIndex(); + p = pTmp; + pcexP = pP; + pcexPx = pPx; + pcexPy = pPy; + pcexPz = pPz; + } + } + } + + if (k_plane == k_e && k_plane == k_p && k_plane >= 0) { + int pId = k_plane; + TVector3 pVecProton = TVector3(pcexPx, pcexPy, pcexPz); + TVector3 pVecAntiproton = TVector3(antipPx, antipPy, antipPz); + TVector3 total_mc_pVec = pVecProton + pVecAntiproton; + double cexPairMcP = total_mc_pVec.Mag(); + double cexPairMcPt = total_mc_pVec.Pt(); + double cexPairMcPz = pcexPz + antipPz; + double mcangleRad = pVecProton.Angle(pVecAntiproton); + double mcangleDeg = mcangleRad * Rad2Deg; + + // Antineutron mother + if (motherPdg == -kNeutron) { + // CEX pair + histos.fill(HIST("cexPairMcP"), cexPairMcP); + histos.fill(HIST("cexPairMcPt"), cexPairMcPt); + histos.fill(HIST("cexPairMcPz"), cexPairMcPz); + histos.fill(HIST("cex_pairmcDplane"), dplane); + histos.fill(HIST("cex_pairmc_angle"), mcangleDeg); + histos.fill(HIST("cex_pairmc_vtx"), antipVx, antipVy); + histos.fill(HIST("cex_pairmc_vtxz"), antipVz); + if (std::abs(motherEta) < kStrictEta && std::abs(motherVz) < kAccMaxVz) + histos.fill(HIST("cexPairMcPITScuts"), cexPairMcP); + // CEX pair normalized + if (motherP != 0) + histos.fill(HIST("cexn_pairmc_p"), cexPairMcP / motherP); + if (motherPt != 0) + histos.fill(HIST("cexn_pairmc_pt"), cexPairMcPt / motherPt); + if (motherPz != 0) + histos.fill(HIST("cexn_pairmc_pz"), cexPairMcPz / motherPz); + } + // BG mother + if (motherPdg != -kNeutron) { + // CEX pair + histos.fill(HIST("cexbg_pairmc_p"), cexPairMcP); + histos.fill(HIST("cexbg_pairmc_pt"), cexPairMcPt); + histos.fill(HIST("cexbg_pairmc_pz"), cexPairMcPz); + histos.fill(HIST("cexbg_pairmcDplane"), dplane); + histos.fill(HIST("cexbg_pairmc_angle"), mcangleDeg); + histos.fill(HIST("cexbg_pairmc_vtx"), antipVx, antipVy); + histos.fill(HIST("cexbg_pairmc_vtxz"), antipVz); + if (std::abs(motherEta) < kStrictEta && std::abs(motherVz) < kAccMaxVz) + histos.fill(HIST("cexbg_pairmc_pITScuts"), cexPairMcP); + } + + // Detector signal + bool antipLayers = false; + bool antipHasTrack = false; + double antipTrkPx = 0.; + double antipTrkPy = 0.; + double antipTrkPz = 0.; + double antipTrkP = 0.; + double antipTrkEta = 0.; + double antipTrkTpcSignal = 0; + // int antip_trk_nClsTPC = 0; + int antipTrkNClsIts = 0; + uint16_t apItsMap = 0; + UShort_t apNClsIts = 0; + + bool pLayers = false; + bool pHasTrack = false; + double pTrkPx = 0.; + double pTrkPy = 0.; + double pTrkPz = 0.; + double pTrkP = 0.; + double pTrkEta = 0.; + double pTrkTpcSignal = 0; + // int p_trk_nClsTPC = 0; + int pTrkNClsIts = 0; + uint16_t pItsMap = 0; + UShort_t pNClsIts = 0; + + for (const auto& track : tracks) { + if (!track.has_mcParticle()) + continue; + const auto& mc = track.mcParticle(); + if (mc.mcCollisionId() != colId) + continue; + uint8_t itsMap = track.itsClusterMap(); + // Config for ITS1 + /*bool hitSPD = (itsMap & 0x3) != 0; // bits 0 (SPD L1) & 1 (SPD L2) + bool hitSDD = (itsMap & 0xC) != 0; // bits 2–3 + bool hitSSD = (itsMap & 0x30) != 0; // bits 4–5 + bool layerCondition = (hitSDD || hitSSD) && !hitSPD;*/ + // Config for ITS2 + bool hitL0 = (itsMap & (1u << 0)) != 0; + bool hitL1 = (itsMap & (1u << 1)) != 0; + bool hitL2 = (itsMap & (1u << 2)) != 0; + bool hitL3 = (itsMap & (1u << 3)) != 0; + bool hitL4 = (itsMap & (1u << 4)) != 0; + bool hitL5 = (itsMap & (1u << 5)) != 0; + bool hitL6 = (itsMap & (1u << 6)) != 0; + bool hitIB = (hitL0 || hitL1 || hitL2); + bool hitOuter = (hitL3 || hitL4 || hitL5 || hitL6); + int nITS = track.itsNCls(); + bool layerCondition = (!hitIB) && hitOuter && (nITS >= kMinItsHits); + + if (mc.globalIndex() == antipId) { + antipTrkP = track.p(); + antipTrkPx = track.px(); + antipTrkPy = track.py(); + antipTrkPz = track.pz(); + antipTrkEta = track.eta(); + antipTrkTpcSignal = track.tpcSignal(); + // antip_trk_nClsTPC = track.tpcNCls(); + antipTrkNClsIts = track.itsNCls(); + antipHasTrack = true; + apItsMap = static_cast(track.itsClusterMap()); + antipLayers = (apItsMap != 0); + if (layerCondition) + antipLayers = true; + } else if (mc.globalIndex() == pId) { + pTrkP = track.p(); + pTrkPx = track.px(); + pTrkPy = track.py(); + pTrkPz = track.pz(); + pTrkEta = track.eta(); + pTrkTpcSignal = track.tpcSignal(); + // p_trk_nClsTPC = track.tpcNCls(); + pTrkNClsIts = track.itsNCls(); + pHasTrack = true; + pItsMap = static_cast(track.itsClusterMap()); + pLayers = (pItsMap != 0); + if (layerCondition) + pLayers = true; + } + } + if (!(pHasTrack && antipHasTrack)) + continue; + + TVector3 pVecProton_trk(pTrkPx, pTrkPy, pTrkPz); + TVector3 AntipVecProton_trk(antipTrkPx, antipTrkPy, antipTrkPz); + TVector3 total_trk_pVec = pVecProton_trk + AntipVecProton_trk; + double trkangleRad = AntipVecProton_trk.Angle(pVecProton_trk); + double trkangleDeg = trkangleRad * Rad2Deg; + if (motherPdg == -kNeutron) + histos.fill(HIST("cex_pairtrk_angle"), trkangleDeg); + if (motherPdg != -kNeutron) + histos.fill(HIST("cexbg_pairtrk_angle"), trkangleDeg); + + // ==== Secondary vertex via central DCA vertexer (DCAFitter2) ==== + using o2::vertexing::DCAFitter2; + constexpr float kBzTesla = 0.5f; + DCAFitter2 fitter(/*bz=*/kBzTesla, /*useAbsDCA=*/true, /*propagateToPCA=*/true); + fitter.setBz(kBzTesla); + // float bz = o2::base::Propagator::Instance()->getNominalBz(); // en kGauss + // DCAFitter2 fitter(bz, /*useAbsDCA=*/true, /*propagateToPCA=*/true); + fitter.setMaxR(45.f); // cm + fitter.setMaxDZIni(4.f); // cm + fitter.setMaxDXYIni(4.f); // cm + fitter.setMaxChi2(50.f); + fitter.setPropagateToPCA(true); + std::optional pRow, apRow; + for (const auto& tr : tracks) { + if (!tr.has_mcParticle()) + continue; + const auto& mc = tr.mcParticle(); + if (mc.globalIndex() == antipId) + apRow = tr; + if (mc.globalIndex() == pId) + pRow = tr; + if (pRow && apRow) + break; + } + if (!(pRow && apRow)) { + } else { + // TrackParCov + auto trP = makeTPCovFromAOD(*pRow); + auto trAP = makeTPCovFromAOD(*apRow); + int nCand = fitter.process(trP, trAP); + auto status = fitter.getFitStatus(); + histos.fill(HIST("vtxfitStatus"), (int)status); + if (nCand > 0 && (status == DCAFitter2::FitStatus::Converged || status == DCAFitter2::FitStatus::MaxIter)) { + // Secondary vertex (commom PCA) [x,y,z] cm + auto vtx = fitter.getPCACandidatePos(); + const double secX = vtx[0]; + const double secY = vtx[1]; + const double secZ = vtx[2]; + // DCA of the pair in the PCA (equivalent to minDCA) + fitter.propagateTracksToVertex(); + auto tp0 = fitter.getTrackParamAtPCA(0); + auto tp1 = fitter.getTrackParamAtPCA(1); + const auto p0 = tp0.getXYZGlo(); + const auto p1 = tp1.getXYZGlo(); + const double x0 = p0.X(), y0 = p0.Y(), z0 = p0.Z(); + const double x1 = p1.X(), y1 = p1.Y(), z1 = p1.Z(); + const double dcaPair = std::sqrt((x0 - x1) * (x0 - x1) + (y0 - y1) * (y0 - y1) + (z0 - z1) * (z0 - z1)); + + if (motherPdg == -kNeutron) + histos.fill(HIST("cex_pairtrkVtxfitDcaPair"), dcaPair); + if (motherPdg != -kNeutron) + histos.fill(HIST("cexbg_pairtrkVtxfitDcaPair"), dcaPair); + + // Cuts + // if (dcaPair > 10.0) continue; + // if (trkangleDeg > 60) continue; + // if (std::abs(dplane) > 2) continue; + if (!(antipLayers && pLayers)) + continue; + double cexPairTrkP = total_trk_pVec.Mag(); + double cexPairTrkPt = total_trk_pVec.Pt(); + double cexPairTrkPz = pTrkPz + antipTrkPz; + const double radius = std::hypot(secX, secY); + const double dxPV = secX - pvtxX; + const double dyPV = secY - pvtxY; + const double dzPV = secZ - pvtxZ; + const double distToPrimary = std::sqrt(dxPV * dxPV + dyPV * dyPV + dzPV * dzPV); + + const TVector3 pv2sv(secX - pvtxX, secY - pvtxY, secZ - pvtxZ); + const double pairPointingAngleDeg = pv2sv.Angle(total_trk_pVec) * Rad2Deg; + + const double pP = pVecProton_trk.Mag(); + const double pAP = AntipVecProton_trk.Mag(); + const double ptP = pVecProton_trk.Pt(); + const double ptAP = AntipVecProton_trk.Pt(); + + const double denomP = std::max(1e-9, pP + pAP); + const double denomPt = std::max(1e-9, ptP + ptAP); + + const float pairPBalance = std::abs(pP - pAP) / denomP; + const float pairPtBalance = std::abs(ptP - ptAP) / denomPt; + + const float pairQ = (pVecProton_trk - AntipVecProton_trk).Mag(); + + // Trk - MC + const float dPairP = cexPairTrkP - cexPairMcP; + const float dPairPt = cexPairTrkPt - cexPairMcPt; + const float dPairPz = cexPairTrkPz - cexPairMcPz; + const float dOpenAngle = trkangleDeg - mcangleDeg; + + // closest ITS layer: Radius need to be checked + static const std::array Rlayers = {2.2, 2.8, 3.6, 19.6, 24.0, 29.0, 35.0}; + int16_t svNearestLayerId = -1; + float svDeltaRtoLayer = 1e9f; + for (int i = 0; i < (int)Rlayers.size(); ++i) { + const float dR = std::abs((float)radius - (float)Rlayers[i]); + if (dR < svDeltaRtoLayer) { + svDeltaRtoLayer = dR; + svNearestLayerId = i; + } + } + + if (motherPdg == -kNeutron) { + histos.fill(HIST("cexPairTrkP"), cexPairTrkP); + histos.fill(HIST("cexPairTrkPt"), cexPairTrkPt); + histos.fill(HIST("cexPairTrkPz"), cexPairTrkPz); + histos.fill(HIST("cex_pairtrkVtxfitR"), radius); + histos.fill(HIST("cex_pairtrkVtxfitDistToPv"), distToPrimary); + histos.fill(HIST("cex_pairtrk_vtxfit_secVtxXY"), secX, secY); + histos.fill(HIST("cex_pairtrk_vtxfit_secVtxZ"), secZ); + } else { + histos.fill(HIST("cexbg_pairtrk_p"), cexPairTrkP); + histos.fill(HIST("cexbg_pairtrk_pt"), cexPairTrkPt); + histos.fill(HIST("cexbg_pairtrk_pz"), cexPairTrkPz); + histos.fill(HIST("cexbg_pairtrkVtxfitR"), radius); + histos.fill(HIST("cexbg_pairtrkVtxfitDistToPv"), distToPrimary); + histos.fill(HIST("cexbg_pairtrk_vtxfit_secVtxXY"), secX, secY); + histos.fill(HIST("cexbg_pairtrk_vtxfit_secVtxZ"), secZ); + } + + const float chi2 = fitter.getChi2AtPCACandidate(); + histos.fill(HIST("vtxfitChi2"), chi2); + const double dX = secX - antipVx; + const double dY = secY - antipVy; + const double dZ = secZ - antipVz; + const double d3D = std::sqrt(dX * dX + dY * dY + dZ * dZ); + histos.fill(HIST("vtxfit_mc_dX"), dX); + histos.fill(HIST("vtxfit_mc_dY"), dY); + histos.fill(HIST("vtxfit_mc_dZ"), dZ); + histos.fill(HIST("vtxfit_mc_d3D"), d3D); + + const bool isCex = (motherPdg == -kNeutron); + + const float vtxfitDX = secX - antipVx; + const float vtxfitDY = secY - antipVy; + const float vtxfitDZ = secZ - antipVz; + const float vtxfitD3D = std::sqrt(vtxfitDX * vtxfitDX + vtxfitDY * vtxfitDY + vtxfitDZ * vtxfitDZ); + + const uint32_t selMask = 0u; + + outPairs( + isCex, + motherPdg, + colId, + pId, + antipId, + + cexPairMcP, + cexPairMcPt, + cexPairMcPz, + dplane, + mcangleDeg, + antipVx, + antipVy, + antipVz, + + cexPairTrkP, + cexPairTrkPt, + cexPairTrkPz, + trkangleDeg, + dcaPair, + radius, + distToPrimary, + secX, + secY, + secZ, + + chi2, + static_cast(status), + nCand, + vtxfitDX, + vtxfitDY, + vtxfitDZ, + vtxfitD3D, + + pTrkP, + pTrkPx, + pTrkPy, + pTrkPz, + pTrkEta, + pTrkTpcSignal, + pTrkNClsIts, + + antipTrkP, + antipTrkPx, + antipTrkPy, + antipTrkPz, + antipTrkEta, + antipTrkTpcSignal, + antipTrkNClsIts, + + selMask, + + pairPointingAngleDeg, + pairPBalance, + pairPtBalance, + pairQ, + + dPairP, + dPairPt, + dPairPz, + dOpenAngle, + + svNearestLayerId, + svDeltaRtoLayer, + + pItsMap, + apItsMap, + (int8_t)(pLayers ? 1 : 0), + (int8_t)(antipLayers ? 1 : 0), + + pvtxZ); + } + } + // ==== end DCAFitter2 ==== + } + } + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& ctx) +{ + return WorkflowSpec{adaptAnalysisTask(ctx)}; +} From 49249fbf77a18d447dc6725f6733d429953960d9 Mon Sep 17 00:00:00 2001 From: FabiolaLP Date: Tue, 30 Sep 2025 21:20:24 -0600 Subject: [PATCH 2/4] Fixed linter comments --- PWGLF/DataModel/LFAntinCexTables.h | 48 +++++++++---------- .../Nuspex/nucleiAntineutronCex.cxx | 27 +++++------ 2 files changed, 36 insertions(+), 39 deletions(-) diff --git a/PWGLF/DataModel/LFAntinCexTables.h b/PWGLF/DataModel/LFAntinCexTables.h index 63a68dbc9a1..f49eab6ad80 100644 --- a/PWGLF/DataModel/LFAntinCexTables.h +++ b/PWGLF/DataModel/LFAntinCexTables.h @@ -23,7 +23,7 @@ namespace o2::aod { -namespace antinCexNS +namespace antin_cex { // Metadata DECLARE_SOA_COLUMN(IsCex, isCex, bool); // 1=CEX (from antin), 0=BG @@ -95,35 +95,35 @@ DECLARE_SOA_COLUMN(DPairPz, dPairPz, float); DECLARE_SOA_COLUMN(DOpenAngle, dOpenAngle, float); DECLARE_SOA_COLUMN(SVNearestLayerId, svNearestLayerId, int16_t); -DECLARE_SOA_COLUMN(SVDeltaRtoLayer, svDeltaRtoLayer, float); +DECLARE_SOA_COLUMN(SVDeltaRToLayer, svDeltaRToLayer, float); DECLARE_SOA_COLUMN(PTrkItsHitMap, pTrkItsHitMap, uint16_t); DECLARE_SOA_COLUMN(APTrkItsHitMap, apTrkItsHitMap, uint16_t); -DECLARE_SOA_COLUMN(PLayersOK, pLayersOK, int8_t); -DECLARE_SOA_COLUMN(APLayersOK, apLayersOK, int8_t); +DECLARE_SOA_COLUMN(PLayersOk, pLayersOk, int8_t); +DECLARE_SOA_COLUMN(APLayersOk, apLayersOk, int8_t); -DECLARE_SOA_COLUMN(PvtxZ, pvtxZ, float); -} // namespace antinCexNS +DECLARE_SOA_COLUMN(PVtxZ, pVtxZ, float); +} // namespace antin_cex // Table -DECLARE_SOA_TABLE(antinCexPairs, "AOD", "ANTINCEX", - antinCexNS::IsCex, - antinCexNS::MotherPdg, antinCexNS::ColId, antinCexNS::PId, antinCexNS::AntipId, - antinCexNS::McPairP, antinCexNS::McPairPt, antinCexNS::McPairPz, - antinCexNS::McDplane, antinCexNS::McAngleDeg, antinCexNS::McVtxX, antinCexNS::McVtxY, antinCexNS::McVtxZ, - antinCexNS::TrkPairP, antinCexNS::TrkPairPt, antinCexNS::TrkPairPz, antinCexNS::TrkAngleDeg, - antinCexNS::TrkVtxfitDcaPair, antinCexNS::TrkVtxfitR, antinCexNS::TrkVtxfitDistToPv, - antinCexNS::TrkVtxfitSecVtxX, antinCexNS::TrkVtxfitSecVtxY, antinCexNS::TrkVtxfitSecVtxZ, - antinCexNS::VtxfitChi2, antinCexNS::VtxfitStatus, antinCexNS::NCand, - antinCexNS::VtxfitDX, antinCexNS::VtxfitDY, antinCexNS::VtxfitDZ, antinCexNS::VtxfitD3D, - antinCexNS::PTrkP, antinCexNS::PTrkPx, antinCexNS::PTrkPy, antinCexNS::PTrkPz, antinCexNS::PTrkEta, antinCexNS::PTrkTpcSignal, antinCexNS::PTrkNClsIts, - antinCexNS::AntipTrkP, antinCexNS::AntipTrkPx, antinCexNS::AntipTrkPy, antinCexNS::AntipTrkPz, antinCexNS::AntipTrkEta, antinCexNS::AntipTrkTpcSignal, antinCexNS::AntipTrkNClsIts, - antinCexNS::SelMask, - antinCexNS::PairPointingAngleDeg, antinCexNS::PairPBalance, antinCexNS::PairPtBalance, antinCexNS::PairQ, - antinCexNS::DPairP, antinCexNS::DPairPt, antinCexNS::DPairPz, antinCexNS::DOpenAngle, - antinCexNS::SVNearestLayerId, antinCexNS::SVDeltaRtoLayer, - antinCexNS::PTrkItsHitMap, antinCexNS::APTrkItsHitMap, antinCexNS::PLayersOK, antinCexNS::APLayersOK, - antinCexNS::PvtxZ); +DECLARE_SOA_TABLE(AntinCexPairs, "AOD", "ANTINCEX", + antin_cex::IsCex, + antin_cex::MotherPdg, antin_cex::ColId, antin_cex::PId, antin_cex::AntipId, + antin_cex::McPairP, antin_cex::McPairPt, antin_cex::McPairPz, + antin_cex::McDplane, antin_cex::McAngleDeg, antin_cex::McVtxX, antin_cex::McVtxY, antin_cex::McVtxZ, + antin_cex::TrkPairP, antin_cex::TrkPairPt, antin_cex::TrkPairPz, antin_cex::TrkAngleDeg, + antin_cex::TrkVtxfitDcaPair, antin_cex::TrkVtxfitR, antin_cex::TrkVtxfitDistToPv, + antin_cex::TrkVtxfitSecVtxX, antin_cex::TrkVtxfitSecVtxY, antin_cex::TrkVtxfitSecVtxZ, + antin_cex::VtxfitChi2, antin_cex::VtxfitStatus, antin_cex::NCand, + antin_cex::VtxfitDX, antin_cex::VtxfitDY, antin_cex::VtxfitDZ, antin_cex::VtxfitD3D, + antin_cex::PTrkP, antin_cex::PTrkPx, antin_cex::PTrkPy, antin_cex::PTrkPz, antin_cex::PTrkEta, antin_cex::PTrkTpcSignal, antin_cex::PTrkNClsIts, + antin_cex::AntipTrkP, antin_cex::AntipTrkPx, antin_cex::AntipTrkPy, antin_cex::AntipTrkPz, antin_cex::AntipTrkEta, antin_cex::AntipTrkTpcSignal, antin_cex::AntipTrkNClsIts, + antin_cex::SelMask, + antin_cex::PairPointingAngleDeg, antin_cex::PairPBalance, antin_cex::PairPtBalance, antin_cex::PairQ, + antin_cex::DPairP, antin_cex::DPairPt, antin_cex::DPairPz, antin_cex::DOpenAngle, + antin_cex::SVNearestLayerId, antin_cex::SVDeltaRToLayer, + antin_cex::PTrkItsHitMap, antin_cex::APTrkItsHitMap, antin_cex::PLayersOk, antin_cex::APLayersOk, + antin_cex::PVtxZ); } // namespace o2::aod diff --git a/PWGLF/TableProducer/Nuspex/nucleiAntineutronCex.cxx b/PWGLF/TableProducer/Nuspex/nucleiAntineutronCex.cxx index 0c96de014ce..232208917a9 100644 --- a/PWGLF/TableProducer/Nuspex/nucleiAntineutronCex.cxx +++ b/PWGLF/TableProducer/Nuspex/nucleiAntineutronCex.cxx @@ -25,6 +25,7 @@ #include "TMCProcess.h" +#include #include // ROOT #include "CommonConstants/MathConstants.h" @@ -61,7 +62,7 @@ struct NucleiAntineutronCex { static constexpr double kVtxTol = 1e-4; HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - Produces outPairs; + Produces outPairs; void init(InitContext const&) { @@ -619,7 +620,7 @@ struct NucleiAntineutronCex { auto trAP = makeTPCovFromAOD(*apRow); int nCand = fitter.process(trP, trAP); auto status = fitter.getFitStatus(); - histos.fill(HIST("vtxfitStatus"), (int)status); + histos.fill(HIST("vtxfitStatus"), static_cast(status)); if (nCand > 0 && (status == DCAFitter2::FitStatus::Converged || status == DCAFitter2::FitStatus::MaxIter)) { // Secondary vertex (commom PCA) [x,y,z] cm auto vtx = fitter.getPCACandidatePos(); @@ -641,10 +642,6 @@ struct NucleiAntineutronCex { if (motherPdg != -kNeutron) histos.fill(HIST("cexbg_pairtrkVtxfitDcaPair"), dcaPair); - // Cuts - // if (dcaPair > 10.0) continue; - // if (trkangleDeg > 60) continue; - // if (std::abs(dplane) > 2) continue; if (!(antipLayers && pLayers)) continue; double cexPairTrkP = total_trk_pVec.Mag(); @@ -681,12 +678,12 @@ struct NucleiAntineutronCex { // closest ITS layer: Radius need to be checked static const std::array Rlayers = {2.2, 2.8, 3.6, 19.6, 24.0, 29.0, 35.0}; int16_t svNearestLayerId = -1; - float svDeltaRtoLayer = 1e9f; - for (int i = 0; i < (int)Rlayers.size(); ++i) { - const float dR = std::abs((float)radius - (float)Rlayers[i]); - if (dR < svDeltaRtoLayer) { - svDeltaRtoLayer = dR; - svNearestLayerId = i; + float SVDeltaRToLayer = 1e9f; + for (int i = 0; i < static_cast(Rlayers.size()); ++i) { + const float dR = static_cast(std::abs(radius - Rlayers[i])); + if (dR < SVDeltaRToLayer) { + SVDeltaRToLayer = dR; + svNearestLayerId = static_cast(i); } } @@ -792,12 +789,12 @@ struct NucleiAntineutronCex { dOpenAngle, svNearestLayerId, - svDeltaRtoLayer, + SVDeltaRToLayer, pItsMap, apItsMap, - (int8_t)(pLayers ? 1 : 0), - (int8_t)(antipLayers ? 1 : 0), + static_cast(pLayers ? 1 : 0), + static_cast(antipLayers ? 1 : 0), pvtxZ); } From 72a299c5f4c1c45a78311911570b8324d5175ff6 Mon Sep 17 00:00:00 2001 From: FabiolaLP Date: Sat, 4 Oct 2025 20:04:25 -0600 Subject: [PATCH 3/4] chore: remove unused vars and tidy selectors to satisfy -Werror --- .../Nuspex/nucleiAntineutronCex.cxx | 676 ++++++++---------- 1 file changed, 316 insertions(+), 360 deletions(-) diff --git a/PWGLF/TableProducer/Nuspex/nucleiAntineutronCex.cxx b/PWGLF/TableProducer/Nuspex/nucleiAntineutronCex.cxx index 232208917a9..023f1283b57 100644 --- a/PWGLF/TableProducer/Nuspex/nucleiAntineutronCex.cxx +++ b/PWGLF/TableProducer/Nuspex/nucleiAntineutronCex.cxx @@ -15,27 +15,23 @@ /// \author Fabiola Lugo /// -#include "DCAFitter/DCAFitterN.h" -#include "DetectorsBase/Propagator.h" -#include "Framework/AnalysisDataModel.h" +#include #include "Framework/AnalysisTask.h" -#include "Framework/Logger.h" #include "Framework/runDataProcessing.h" +#include "Framework/Logger.h" +#include "Framework/AnalysisDataModel.h" +#include "DCAFitter/DCAFitterN.h" #include "ReconstructionDataFormats/TrackParametrization.h" - +#include "DetectorsBase/Propagator.h" #include "TMCProcess.h" - -#include #include -// ROOT -#include "CommonConstants/MathConstants.h" - -#include "TMath.h" -#include "TPDGCode.h" +//ROOT #include "TVector3.h" - -#include +#include "TMath.h" #include +#include +#include "CommonConstants/MathConstants.h" +#include "TPDGCode.h" // #include "PWGLF/DataModel/LFAntinCexTables.h" @@ -46,91 +42,90 @@ using o2::constants::math::Rad2Deg; struct NucleiAntineutronCex { // Slicing per colision Preslice perMcByColl = aod::mcparticle::mcCollisionId; - // Check available tables in the AOD, specifically TracksIU, TracksCovIU + //Check available tables in the AOD, specifically TracksIU, TracksCovIU using TracksWCovMc = soa::Join; - + // === Cut values === - static constexpr double kIts2MinR = 2.2; // ITS2 min radius [cm] - static constexpr double kIts2MaxR = 39.0; // ITS2 max radius [cm] - static constexpr double kIts2MaxVz = 39.0; // ITS2 max |vz| [cm] - static constexpr double kAccMaxEta = 1.2; // acceptance |eta| - static constexpr double kAccMaxVz = 5.3; // acceptance |vz| [cm] - static constexpr double kStrictEta = 0.9; // tighter eta cut - static constexpr double kInitDplane = 10.0; // init dplane - static constexpr double kHuge = 1e9; // fallback for bad denom - static constexpr int kMinItsHits = 2; - static constexpr double kVtxTol = 1e-4; - + static constexpr double kIts2MinR = 2.2; // ITS2 min radius [cm] + static constexpr double kIts2MaxR = 39.0; // ITS2 max radius [cm] + static constexpr double kIts2MaxVz = 39.0; // ITS2 max |vz| [cm] + static constexpr double kAccMaxEta = 1.2; // acceptance |eta| + static constexpr double kAccMaxVz = 5.3; // acceptance |vz| [cm] + static constexpr double kStrictEta = 0.9; // tighter eta cut + static constexpr double kInitDplane = 10.0; // init dplane + static constexpr double kHuge = 1e9; // fallback for bad denom + static constexpr int kMinItsHits = 2; + static constexpr double kVtxTol = 1e-4; + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; Produces outPairs; - - void init(InitContext const&) - { + + void init(InitContext const&) { // Primary vertex histos.add("hVx", "Primary vertex X;X (cm);Entries", kTH1F, {{100, -5., 5.}}); histos.add("hVy", "Primary vertex Y;Y (cm);Entries", kTH1F, {{100, -5., 5.}}); histos.add("hVz", "Primary vertex Z;Z (cm);Entries", kTH1F, {{200, -20., 20.}}); - + // Primary antineutrons histos.add("antin_p", "Total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); histos.add("antin_px", "p_{x};p_{x} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); histos.add("antin_py", "p_{y};p_{y} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); histos.add("antin_pz", "p_{z};p_{z} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); - histos.add("antin_eta", "Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); + histos.add("antin_eta","Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); histos.add("antin_p_ITScuts", "Momentum with ITS cuts;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); - + // Primary neutrons - histos.add("n_p", "Total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("n_p", "Total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); histos.add("n_px", "p_{x};p_{x} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); histos.add("n_py", "p_{y};p_{y} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); histos.add("n_pz", "p_{z};p_{z} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); - histos.add("n_eta", "Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); + histos.add("n_eta","Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); histos.add("n_p_ITScuts", "Momentum with ITS cuts;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); - + // Primary antiprotons - histos.add("antipP", "Total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("antipP", "Total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); histos.add("antipPx", "p_{x};p_{x} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); histos.add("antipPy", "p_{y};p_{y} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); histos.add("antipPz", "p_{z};p_{z} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); - histos.add("antipEta", "Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); + histos.add("antipEta","Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); histos.add("antipP_ITScuts", "Momentum with ITS cuts;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); - + // Primary protons histos.add("pP", "Total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); histos.add("pPx", "p_{x};p_{x} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); histos.add("pPy", "p_{y};p_{y} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); histos.add("pPz", "p_{z};p_{z} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); - histos.add("pEta", "Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); + histos.add("pEta","Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); histos.add("pP_ITScuts", "Momentum with ITS cuts;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); - + // test (MC) histos.add("antip_test", "Secondary antiprotons;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); - + // CEX pair from antineutron (MC) histos.add("cexPairMcP", "CEX pair total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); histos.add("cexPairMcPt", "CEX pair p_{T};p_{T} (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); histos.add("cexPairMcPz", "CEX pair p_{z};p_{z} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); - histos.add("cex_pairmcDplane", "CEX pair d_{plane};d_{plane} (cm);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("cex_pairmcDplane","CEX pair d_{plane};d_{plane} (cm);Entries", kTH1F, {{100, 0., 10.}}); histos.add("cex_pairmc_angle", "Pair opening angle;Angle (°);Entries", kTH1F, {{180, 0., 180.}}); histos.add("cex_pairmc_vtx", "MC CEX pair vertex;X (cm);Y (cm)", kTH2F, {{100, -50., 50.}, {100, -50., 50.}}); histos.add("cex_pairmc_vtxz", "MC secondary vertex Z;Z (cm);Entries", kTH1F, {{200, -60., 60.}}); - histos.add("cexPairMcPITScuts", "CEX pair momentum (ITS cuts);|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); - + histos.add("cexPairMcPITScuts","CEX pair momentum (ITS cuts);|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + // CEX pair normalized to antineutron (MC) histos.add("cexn_pairmc_p", "Pair p / antineutron p;p/p_{#bar{n}};Entries", kTH1F, {{100, 0., 2.}}); histos.add("cexn_pairmc_pt", "Pair p_{T} / antineutron p_{T};p_{T}/p_{T,#bar{n}};Entries", kTH1F, {{100, 0., 2.}}); histos.add("cexn_pairmc_pz", "Pair p_{z} / antineutron p_{z};p_{z}/p_{z,#bar{n}};Entries", kTH1F, {{100, -2., 2.}}); - + // BG pair (not from antineutron) (MC) histos.add("cexbg_pairmc_p", "Background pair momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); histos.add("cexbg_pairmc_pt", "Background pair p_{T};p_{T} (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); histos.add("cexbg_pairmc_pz", "Background pair p_{z};p_{z} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); - histos.add("cexbg_pairmcDplane", "Background d_{plane};d_{plane} (cm);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("cexbg_pairmcDplane","Background d_{plane};d_{plane} (cm);Entries", kTH1F, {{100, 0., 10.}}); histos.add("cexbg_pairmc_angle", "Background opening angle;Angle (°);Entries", kTH1F, {{180, 0., 180.}}); histos.add("cexbg_pairmc_vtx", "Background pair vertex;X (cm);Y (cm)", kTH2F, {{100, -50., 50.}, {100, -50., 50.}}); histos.add("cexbg_pairmc_vtxz", "Background secondary vertex Z;Z (cm);Entries", kTH1F, {{200, -60., 60.}}); - histos.add("cexbg_pairmc_pITScuts", "Background momentum (ITS cuts);|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); - + histos.add("cexbg_pairmc_pITScuts","Background momentum (ITS cuts);|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + // CEX pair from antineutron (TRK) histos.add("cex_pairtrk_angle", "Pair opening angle (tracks);Angle (°);Entries", kTH1F, {{180, 0., 180.}}); histos.add("cexPairTrkP", "Pair momentum (tracks);|p| (GeV/c);Entries", kTH1F, {{120, 0., 12.}}); @@ -141,7 +136,7 @@ struct NucleiAntineutronCex { histos.add("cex_pairtrkVtxfitDistToPv", "Distance from secondary vertex to PV;dist (cm);Entries", kTH1F, {{240, 0., 120.}}); histos.add("cex_pairtrk_vtxfit_secVtxXY", "Secondary vertex (PCA);X (cm);Y (cm)", kTH2F, {{200, -60., 60.}, {200, -60., 60.}}); histos.add("cex_pairtrk_vtxfit_secVtxZ", "Secondary vertex Z (PCA);Z (cm);Entries", kTH1F, {{240, -60., 60.}}); - + // BG pair (not from antineutron) (TRK) histos.add("cexbg_pairtrk_angle", "Background opening angle (tracks);Angle (°);Entries", kTH1F, {{180, 0., 180.}}); histos.add("cexbg_pairtrk_p", "Pair momentum (tracks);|p| (GeV/c);Entries", kTH1F, {{120, 0., 12.}}); @@ -152,44 +147,43 @@ struct NucleiAntineutronCex { histos.add("cexbg_pairtrkVtxfitDistToPv", "Distance from secondary vertex to PV;dist (cm);Entries", kTH1F, {{240, 0., 120.}}); histos.add("cexbg_pairtrk_vtxfit_secVtxXY", "Secondary vertex (PCA);X (cm);Y (cm)", kTH2F, {{200, -60., 60.}, {200, -60., 60.}}); histos.add("cexbg_pairtrk_vtxfit_secVtxZ", "Secondary vertex Z (PCA);Z (cm);Entries", kTH1F, {{240, -60., 60.}}); - + // Vertex fit (DCAFitter2 / PCA) histos.add("vtxfitChi2", "DCAFitter2 #chi^{2};#chi^{2};Entries", kTH1F, {{200, 0., 100.}}); histos.add("vtxfitStatus", "Fit status (0=OK);code;Entries", kTH1I, {{10, 0., 10.}}); histos.add("vtxfit_mc_dX", "SV residual X (fit - MC);#Delta X (cm);Entries", kTH1F, {{400, -20., 20.}}); histos.add("vtxfit_mc_dY", "SV residual Y (fit - MC);#Delta Y (cm);Entries", kTH1F, {{400, -20., 20.}}); histos.add("vtxfit_mc_dZ", "SV residual Z (fit - MC);#Delta Z (cm);Entries", kTH1F, {{400, -20., 20.}}); - histos.add("vtxfit_mc_d3D", "SV distance |fit - MC|;#Delta r (cm);Entries", kTH1F, {{300, 0., 30.}}); + histos.add("vtxfit_mc_d3D", "SV distance |fit - MC|;#Delta r (cm);Entries", kTH1F, {{300, 0., 30.}}); } - - static o2::track::TrackParCov makeTPCovFromAOD(const TracksWCovMc::iterator& tr) - { + + static o2::track::TrackParCov makeTPCovFromAOD(const TracksWCovMc::iterator& tr) { using o2::track::TrackParCov; TrackParCov tpcov; // Local state: x, alpha, y, z, snp, tgl, q/pt float par[5] = {tr.y(), tr.z(), tr.snp(), tr.tgl(), tr.signed1Pt()}; tpcov.set(tr.x(), tr.alpha(), par); - + // Covariance matrix (15 terms) in O2 order std::array cov = { - tr.cYY(), tr.cZY(), tr.cZZ(), - tr.cSnpY(), tr.cSnpZ(), tr.cSnpSnp(), - tr.cTglY(), tr.cTglZ(), tr.cTglSnp(), - tr.cTglTgl(), tr.c1PtY(), tr.c1PtZ(), - tr.c1PtSnp(), tr.c1PtTgl(), tr.c1Pt21Pt2()}; + tr.cYY(), tr.cZY(), tr.cZZ(), + tr.cSnpY(), tr.cSnpZ(), tr.cSnpSnp(), + tr.cTglY(), tr.cTglZ(), tr.cTglSnp(), + tr.cTglTgl(), tr.c1PtY(), tr.c1PtZ(), + tr.c1PtSnp(), tr.c1PtTgl(), tr.c1Pt21Pt2() + }; tpcov.setCov(cov); return tpcov; } - - void process(aod::McCollisions const& cols, aod::McParticles const& particles, TracksWCovMc const& tracks) - { + + void process(aod::McCollisions const& cols, aod::McParticles const& particles, TracksWCovMc const& tracks) { double pvtxX = 0; double pvtxY = 0; double pvtxZ = 0; for (auto const& col : cols) { - const auto colId = col.globalIndex(); + const auto colId= col.globalIndex(); auto mcPartsThis = particles.sliceBy(perMcByColl, colId); - + if (std::isfinite(col.posX()) && std::isfinite(col.posY()) && std::isfinite(col.posZ())) { pvtxX = col.posX(); pvtxY = col.posY(); @@ -198,9 +192,9 @@ struct NucleiAntineutronCex { histos.fill(HIST("hVy"), pvtxY); histos.fill(HIST("hVz"), pvtxZ); } - + for (const auto& particle : mcPartsThis) { - + // Primary antineutrons if (particle.pdgCode() == -kNeutron && particle.isPhysicalPrimary()) { histos.fill(HIST("antin_p"), particle.p()); @@ -208,8 +202,7 @@ struct NucleiAntineutronCex { histos.fill(HIST("antin_py"), particle.py()); histos.fill(HIST("antin_pz"), particle.pz()); histos.fill(HIST("antin_eta"), particle.eta()); - if (std::abs(particle.eta()) < kAccMaxEta && std::abs(particle.vz()) < kAccMaxVz) - histos.fill(HIST("antin_p_ITScuts"), particle.p()); + if (std::abs(particle.eta()) < kAccMaxEta && std::abs(particle.vz())< kAccMaxVz) histos.fill(HIST("antin_p_ITScuts"),particle.p()); } // Primary neutrons if (particle.pdgCode() == kNeutron && particle.isPhysicalPrimary()) { @@ -218,8 +211,7 @@ struct NucleiAntineutronCex { histos.fill(HIST("n_py"), particle.py()); histos.fill(HIST("n_pz"), particle.pz()); histos.fill(HIST("n_eta"), particle.eta()); - if (std::abs(particle.eta()) < kAccMaxEta && std::abs(particle.vz()) < kAccMaxVz) - histos.fill(HIST("n_p_ITScuts"), particle.p()); + if (std::abs(particle.eta()) < kAccMaxEta && std::abs(particle.vz())< kAccMaxVz) histos.fill(HIST("n_p_ITScuts"),particle.p()); } // Primary antiprotons if (particle.pdgCode() == -kProton && particle.isPhysicalPrimary()) { @@ -228,8 +220,7 @@ struct NucleiAntineutronCex { histos.fill(HIST("antipPy"), particle.py()); histos.fill(HIST("antipPz"), particle.pz()); histos.fill(HIST("antipEta"), particle.eta()); - if (std::abs(particle.eta()) < kAccMaxEta && std::abs(particle.vz()) < kAccMaxVz) - histos.fill(HIST("antipP_ITScuts"), particle.p()); + if (std::abs(particle.eta()) < kAccMaxEta && std::abs(particle.vz())< kAccMaxVz) histos.fill(HIST("antipP_ITScuts"),particle.p()); } // Primary protons if (particle.pdgCode() == kProton && particle.isPhysicalPrimary()) { @@ -238,71 +229,62 @@ struct NucleiAntineutronCex { histos.fill(HIST("pPy"), particle.py()); histos.fill(HIST("pPz"), particle.pz()); histos.fill(HIST("pEta"), particle.eta()); - if (std::abs(particle.eta()) < kAccMaxEta && std::abs(particle.vz()) < kAccMaxVz) - histos.fill(HIST("pP_ITScuts"), particle.p()); + if (std::abs(particle.eta()) < kAccMaxEta && std::abs(particle.vz())< kAccMaxVz) histos.fill(HIST("pP_ITScuts"),particle.p()); } - - // Seconday antiprotons from material + + //Seconday antiprotons from material const auto procEnum = particle.getProcess(); const bool isSecondaryFromMaterial = (!particle.producedByGenerator()) && (procEnum == kPHadronic || procEnum == kPHInhelastic); - if (particle.pdgCode() != -kProton || !isSecondaryFromMaterial || particle.mothersIds().empty()) - continue; + if (particle.pdgCode() != -kProton || !isSecondaryFromMaterial || particle.mothersIds().empty()) continue; histos.fill(HIST("antip_test"), particle.p()); - - // Primary mother + + //Primary mother bool hasPrimaryMotherAntip = false; double motherPt = 0.0; double motherPz = 0.0; double motherVz = 0.0; - double motherP = 0.0; + double motherP = 0.0; double motherEta = 0.0; - int motherPdg = 0; - + int motherPdg = 0; + for (const auto& mother : particle.mothers_as()) { - if (mother.isPhysicalPrimary()) { + if (mother.isPhysicalPrimary()){ hasPrimaryMotherAntip = true; motherPt = mother.pt(); motherPz = mother.pz(); motherVz = mother.vz(); - motherP = mother.p(); + motherP= mother.p(); motherEta = mother.eta(); motherPdg = mother.pdgCode(); break; } } - if (!hasPrimaryMotherAntip) - continue; - + if (!hasPrimaryMotherAntip) continue; + double antipVx = particle.vx(); double antipVy = particle.vy(); double antipVz = particle.vz(); - double antipP = particle.p(); double antipPx = particle.px(); double antipPy = particle.py(); double antipPz = particle.pz(); double antipE = particle.e(); - double antipEta = particle.eta(); int antipId = particle.globalIndex(); - - // Selection conditions: Produced in the ITS - const double R = std::sqrt(antipVx * antipVx + antipVy * antipVy); - // Config for ITS - // if(3.9<=R && R<=43.0 && std::abs(antipVz)<=48.9){ - // Config for ITS2 - if (R < kIts2MinR || R > kIts2MaxR || std::abs(antipVz) > kIts2MaxVz) - continue; - if (std::abs(motherEta) >= kAccMaxEta || std::abs(motherVz) >= kAccMaxVz) - continue; - - // Pion minus veto + + //Selection conditions: Produced in the ITS + const double R = std::sqrt(antipVx*antipVx + antipVy*antipVy); + //Config for ITS + //if(3.9<=R && R<=43.0 && std::abs(antipVz)<=48.9){ + //Config for ITS2 + if (R < kIts2MinR || R > kIts2MaxR || std::abs(antipVz) > kIts2MaxVz) continue; + if (std::abs(motherEta) >= kAccMaxEta || std::abs(motherVz) >= kAccMaxVz) continue; + + //Pion minus veto bool pionMinus = false; for (const auto& particle1 : mcPartsThis) { const auto proc1Enum = particle1.getProcess(); const bool isSecondaryFromMaterial1 = (!particle1.producedByGenerator()) && (proc1Enum == kPHadronic || proc1Enum == kPHInhelastic); - if (particle1.mcCollisionId() != colId) - continue; - if (particle1.pdgCode() != kPiMinus || !isSecondaryFromMaterial1 || particle1.mothersIds().empty()) - continue; + if (particle1.mcCollisionId() != colId) continue; + if (particle1.pdgCode() != kPiMinus || !isSecondaryFromMaterial1 || particle1.mothersIds().empty()) continue; bool hasPrimaryMotherPim = false; for (const auto& mother : particle1.mothers_as()) { if (mother.isPhysicalPrimary()) { @@ -310,26 +292,23 @@ struct NucleiAntineutronCex { break; } } - if (!hasPrimaryMotherPim) - continue; + if (!hasPrimaryMotherPim) continue; double pimVx = particle1.vx(); double pimVy = particle1.vy(); double pimVz = particle1.vz(); - if (std::abs(pimVx - antipVx) < kVtxTol && std::abs(pimVy - antipVy) < kVtxTol && std::abs(pimVz - antipVz) < kVtxTol) { + if(std::abs(pimVx - antipVx) < kVtxTol && std::abs(pimVy - antipVy) < kVtxTol && std::abs(pimVz - antipVz) < kVtxTol){ pionMinus = true; break; } } - - // Pion plus veto + + //Pion plus veto bool pionPlus = false; for (const auto& particle2 : mcPartsThis) { - if (particle2.mcCollisionId() != colId) - continue; + if (particle2.mcCollisionId() != colId) continue; const auto proc2Enum = particle2.getProcess(); const bool isSecondaryFromMaterial2 = (!particle2.producedByGenerator()) && (proc2Enum == kPHadronic || proc2Enum == kPHInhelastic); - if (particle2.pdgCode() != kPiPlus || !isSecondaryFromMaterial2 || particle2.mothersIds().empty()) - continue; + if (particle2.pdgCode() != kPiPlus || !isSecondaryFromMaterial2 || particle2.mothersIds().empty()) continue; bool hasPrimaryMotherPip = false; for (const auto& mother : particle2.mothers_as()) { if (mother.isPhysicalPrimary()) { @@ -337,20 +316,18 @@ struct NucleiAntineutronCex { break; } } - if (!hasPrimaryMotherPip) - continue; + if (!hasPrimaryMotherPip) continue; double pipVx = particle2.vx(); double pipVy = particle2.vy(); double pipVz = particle2.vz(); - if (std::abs(pipVx - antipVx) < kVtxTol && std::abs(pipVy - antipVy) < kVtxTol && std::abs(pipVz - antipVz) < kVtxTol) { + if(std::abs(pipVx - antipVx) < kVtxTol && std::abs(pipVy - antipVy) < kVtxTol && std::abs(pipVz - antipVz) < kVtxTol){ pionPlus = true; break; } } - - if (pionPlus || pionMinus) - continue; - + + if (pionPlus || pionMinus) continue; + // CEX selection double dplane = kInitDplane; double dplaneTmp = 0; @@ -359,45 +336,38 @@ struct NucleiAntineutronCex { double pcexPx = 0; double pcexPy = 0; double pcexPz = 0; - double pcexP = 0; double e = 0; double eTmp = 0; - double pcexE = 0; int k_plane = -1; int k_e = -1; int k_p = -1; - - // Secondary proton from material + + //Secondary proton from material for (const auto& particle3 : mcPartsThis) { - if (particle3.mcCollisionId() != colId) - continue; + if (particle3.mcCollisionId() != colId) continue; const auto proc3Enum = particle3.getProcess(); const bool isSecondaryFromMaterial3 = (!particle3.producedByGenerator()) && (proc3Enum == kPHadronic || proc3Enum == kPHInhelastic); - if (particle3.pdgCode() != kProton || !isSecondaryFromMaterial3 || particle3.mothersIds().empty()) - continue; - bool hasPrimaryMotherP = false; + if (particle3.pdgCode() != kProton || !isSecondaryFromMaterial3 || particle3.mothersIds().empty()) continue; + bool hasPrimaryMotherP= false; for (const auto& mother : particle3.mothers_as()) { if (mother.isPhysicalPrimary()) { - hasPrimaryMotherP = true; + hasPrimaryMotherP= true; break; } } - if (!hasPrimaryMotherP) - continue; + if (!hasPrimaryMotherP) continue; double pVx = particle3.vx(); double pVy = particle3.vy(); double pVz = particle3.vz(); - double pP = particle3.p(); double pPx = particle3.px(); double pPy = particle3.py(); double pPz = particle3.pz(); double pE = particle3.e(); - double pEta = particle3.eta(); - if (std::abs(pVx - antipVx) < kVtxTol && std::abs(pVy - antipVy) < kVtxTol && std::abs(pVz - antipVz) < kVtxTol) { - // Same mother + if(std::abs(pVx - antipVx) < kVtxTol && std::abs(pVy - antipVy) < kVtxTol && std::abs(pVz - antipVz) < kVtxTol){ + //Same mother bool shareMother = false; - const auto& momsAp = particle.mothersIds(); // antiproton - const auto& momsP = particle3.mothersIds(); // proton + const auto& momsAp = particle.mothersIds(); //antiproton + const auto& momsP = particle3.mothersIds(); //proton for (const auto& ida : momsAp) { for (const auto& idp : momsP) { if (ida == idp) { @@ -405,64 +375,66 @@ struct NucleiAntineutronCex { break; } } - if (shareMother) - break; + if (shareMother) break; } - if (!shareMother) - continue; - - // CEX proton selection - // dplaneTmp = (pPy*antipPz - pPz*antipPy)*(pvtxX-antipVx) + (pPz*antipPx - pPx*antipPz)*(pvtxY-antipVy) + (pPx*antipPy - pPy*antipPx)*(pvtxZ-antipVz); - double nx = (pPy * antipPz - pPz * antipPy); - double ny = (pPz * antipPx - pPx * antipPz); - double nz = (pPx * antipPy - pPy * antipPx); + if (!shareMother) continue; + + //CEX proton selection + //dplaneTmp = (pPy*antipPz - pPz*antipPy)*(pvtxX-antipVx) + (pPz*antipPx - pPx*antipPz)*(pvtxY-antipVy) + (pPx*antipPy - pPy*antipPx)*(pvtxZ-antipVz); + double nx = (pPy*antipPz - pPz*antipPy); + double ny = (pPz*antipPx - pPx*antipPz); + double nz = (pPx*antipPy - pPy*antipPx); double rx = (pvtxX - antipVx); double ry = (pvtxY - antipVy); double rz = (pvtxZ - antipVz); - double denom = nx * nx + ny * ny + nz * nz; - if (denom > 0.) { - dplaneTmp = std::abs(nx * rx + ny * ry + nz * rz) / std::sqrt(denom); - } else { + double denom = nx*nx + ny*ny + nz*nz; + if (denom > 0.) + { + dplaneTmp = std::abs(nx*rx + ny*ry + nz*rz) / std::sqrt(denom); + } + else + { dplaneTmp = kHuge; } - if (std::abs(dplaneTmp) < std::abs(dplane)) { + if(std::abs(dplaneTmp) < std::abs(dplane)) + { k_plane = particle3.globalIndex(); dplane = dplaneTmp; } - + eTmp = antipE + pE; - if (std::abs(eTmp) > std::abs(e)) { + if(std::abs(eTmp) > std::abs(e)) + { k_e = particle3.globalIndex(); e = eTmp; - pcexE = pE; } - - pTmp = std::sqrt(std::pow((pPx + antipPx), 2) + std::pow((pPy + antipPy), 2) + std::pow((pPz + antipPz), 2)); - if (std::abs(pTmp) > std::abs(p)) { + + pTmp = std::sqrt(std::pow((pPx+antipPx),2)+std::pow((pPy+antipPy),2)+std::pow((pPz+antipPz),2)); + if(std::abs(pTmp) > std::abs(p)) + { k_p = particle3.globalIndex(); p = pTmp; - pcexP = pP; - pcexPx = pPx; - pcexPy = pPy; - pcexPz = pPz; + pcexPx = pPx; + pcexPy = pPy; + pcexPz = pPz; } } } - - if (k_plane == k_e && k_plane == k_p && k_plane >= 0) { + + if(k_plane == k_e && k_plane == k_p && k_plane >= 0 ){ int pId = k_plane; TVector3 pVecProton = TVector3(pcexPx, pcexPy, pcexPz); TVector3 pVecAntiproton = TVector3(antipPx, antipPy, antipPz); TVector3 total_mc_pVec = pVecProton + pVecAntiproton; double cexPairMcP = total_mc_pVec.Mag(); double cexPairMcPt = total_mc_pVec.Pt(); - double cexPairMcPz = pcexPz + antipPz; + double cexPairMcPz = pcexPz+antipPz; double mcangleRad = pVecProton.Angle(pVecAntiproton); double mcangleDeg = mcangleRad * Rad2Deg; - - // Antineutron mother + + //Antineutron mother if (motherPdg == -kNeutron) { - // CEX pair + //CEX pair histos.fill(HIST("cexPairMcP"), cexPairMcP); histos.fill(HIST("cexPairMcPt"), cexPairMcPt); histos.fill(HIST("cexPairMcPz"), cexPairMcPz); @@ -470,19 +442,15 @@ struct NucleiAntineutronCex { histos.fill(HIST("cex_pairmc_angle"), mcangleDeg); histos.fill(HIST("cex_pairmc_vtx"), antipVx, antipVy); histos.fill(HIST("cex_pairmc_vtxz"), antipVz); - if (std::abs(motherEta) < kStrictEta && std::abs(motherVz) < kAccMaxVz) - histos.fill(HIST("cexPairMcPITScuts"), cexPairMcP); - // CEX pair normalized - if (motherP != 0) - histos.fill(HIST("cexn_pairmc_p"), cexPairMcP / motherP); - if (motherPt != 0) - histos.fill(HIST("cexn_pairmc_pt"), cexPairMcPt / motherPt); - if (motherPz != 0) - histos.fill(HIST("cexn_pairmc_pz"), cexPairMcPz / motherPz); + if (std::abs(motherEta) < kStrictEta && std::abs(motherVz)< kAccMaxVz ) histos.fill(HIST("cexPairMcPITScuts"),cexPairMcP); + //CEX pair normalized + if (motherP != 0) histos.fill(HIST("cexn_pairmc_p"), cexPairMcP /motherP); + if (motherPt != 0) histos.fill(HIST("cexn_pairmc_pt"), cexPairMcPt/motherPt); + if (motherPz != 0) histos.fill(HIST("cexn_pairmc_pz"), cexPairMcPz/motherPz); } - // BG mother + //BG mother if (motherPdg != -kNeutron) { - // CEX pair + //CEX pair histos.fill(HIST("cexbg_pairmc_p"), cexPairMcP); histos.fill(HIST("cexbg_pairmc_pt"), cexPairMcPt); histos.fill(HIST("cexbg_pairmc_pz"), cexPairMcPz); @@ -490,144 +458,132 @@ struct NucleiAntineutronCex { histos.fill(HIST("cexbg_pairmc_angle"), mcangleDeg); histos.fill(HIST("cexbg_pairmc_vtx"), antipVx, antipVy); histos.fill(HIST("cexbg_pairmc_vtxz"), antipVz); - if (std::abs(motherEta) < kStrictEta && std::abs(motherVz) < kAccMaxVz) - histos.fill(HIST("cexbg_pairmc_pITScuts"), cexPairMcP); + if (std::abs(motherEta) < kStrictEta && std::abs(motherVz)< kAccMaxVz) histos.fill(HIST("cexbg_pairmc_pITScuts"),cexPairMcP); } - - // Detector signal + + //Detector signal bool antipLayers = false; bool antipHasTrack = false; double antipTrkPx = 0.; double antipTrkPy = 0.; double antipTrkPz = 0.; - double antipTrkP = 0.; + double antipTrkP= 0.; double antipTrkEta = 0.; double antipTrkTpcSignal = 0; - // int antip_trk_nClsTPC = 0; + //int antip_trk_nClsTPC = 0; int antipTrkNClsIts = 0; - uint16_t apItsMap = 0; - UShort_t apNClsIts = 0; - + uint16_t apItsMap = 0; + bool pLayers = false; bool pHasTrack = false; double pTrkPx = 0.; double pTrkPy = 0.; double pTrkPz = 0.; - double pTrkP = 0.; + double pTrkP= 0.; double pTrkEta = 0.; double pTrkTpcSignal = 0; - // int p_trk_nClsTPC = 0; + //int p_trk_nClsTPC = 0; int pTrkNClsIts = 0; - uint16_t pItsMap = 0; - UShort_t pNClsIts = 0; - + uint16_t pItsMap = 0; + for (const auto& track : tracks) { - if (!track.has_mcParticle()) - continue; + if (!track.has_mcParticle()) continue; const auto& mc = track.mcParticle(); - if (mc.mcCollisionId() != colId) - continue; + if (mc.mcCollisionId() != colId) continue; uint8_t itsMap = track.itsClusterMap(); - // Config for ITS1 + //Config for ITS1 /*bool hitSPD = (itsMap & 0x3) != 0; // bits 0 (SPD L1) & 1 (SPD L2) bool hitSDD = (itsMap & 0xC) != 0; // bits 2–3 bool hitSSD = (itsMap & 0x30) != 0; // bits 4–5 bool layerCondition = (hitSDD || hitSSD) && !hitSPD;*/ - // Config for ITS2 - bool hitL0 = (itsMap & (1u << 0)) != 0; - bool hitL1 = (itsMap & (1u << 1)) != 0; - bool hitL2 = (itsMap & (1u << 2)) != 0; - bool hitL3 = (itsMap & (1u << 3)) != 0; - bool hitL4 = (itsMap & (1u << 4)) != 0; - bool hitL5 = (itsMap & (1u << 5)) != 0; - bool hitL6 = (itsMap & (1u << 6)) != 0; + //Config for ITS2 + bool hitL0 = (itsMap & (1u<<0)) != 0; + bool hitL1 = (itsMap & (1u<<1)) != 0; + bool hitL2 = (itsMap & (1u<<2)) != 0; + bool hitL3 = (itsMap & (1u<<3)) != 0; + bool hitL4 = (itsMap & (1u<<4)) != 0; + bool hitL5 = (itsMap & (1u<<5)) != 0; + bool hitL6 = (itsMap & (1u<<6)) != 0; bool hitIB = (hitL0 || hitL1 || hitL2); bool hitOuter = (hitL3 || hitL4 || hitL5 || hitL6); int nITS = track.itsNCls(); bool layerCondition = (!hitIB) && hitOuter && (nITS >= kMinItsHits); - + if (mc.globalIndex() == antipId) { - antipTrkP = track.p(); + antipTrkP= track.p(); antipTrkPx = track.px(); antipTrkPy = track.py(); antipTrkPz = track.pz(); antipTrkEta = track.eta(); antipTrkTpcSignal = track.tpcSignal(); - // antip_trk_nClsTPC = track.tpcNCls(); + //antip_trk_nClsTPC = track.tpcNCls(); antipTrkNClsIts = track.itsNCls(); antipHasTrack = true; - apItsMap = static_cast(track.itsClusterMap()); + apItsMap = static_cast(track.itsClusterMap()); antipLayers = (apItsMap != 0); - if (layerCondition) - antipLayers = true; - } else if (mc.globalIndex() == pId) { - pTrkP = track.p(); + if (layerCondition) antipLayers = true; + } + else if (mc.globalIndex() == pId) { + pTrkP= track.p(); pTrkPx = track.px(); pTrkPy = track.py(); pTrkPz = track.pz(); pTrkEta = track.eta(); pTrkTpcSignal = track.tpcSignal(); - // p_trk_nClsTPC = track.tpcNCls(); + //p_trk_nClsTPC = track.tpcNCls(); pTrkNClsIts = track.itsNCls(); pHasTrack = true; - pItsMap = static_cast(track.itsClusterMap()); - pLayers = (pItsMap != 0); - if (layerCondition) - pLayers = true; + pItsMap = static_cast(track.itsClusterMap()); + pLayers = (pItsMap != 0); + if (layerCondition) pLayers = true; } } - if (!(pHasTrack && antipHasTrack)) - continue; - + if (!(pHasTrack && antipHasTrack)) continue; + TVector3 pVecProton_trk(pTrkPx, pTrkPy, pTrkPz); TVector3 AntipVecProton_trk(antipTrkPx, antipTrkPy, antipTrkPz); TVector3 total_trk_pVec = pVecProton_trk + AntipVecProton_trk; double trkangleRad = AntipVecProton_trk.Angle(pVecProton_trk); double trkangleDeg = trkangleRad * Rad2Deg; - if (motherPdg == -kNeutron) - histos.fill(HIST("cex_pairtrk_angle"), trkangleDeg); - if (motherPdg != -kNeutron) - histos.fill(HIST("cexbg_pairtrk_angle"), trkangleDeg); - + if(motherPdg == -kNeutron) histos.fill(HIST("cex_pairtrk_angle"), trkangleDeg); + if(motherPdg != -kNeutron) histos.fill(HIST("cexbg_pairtrk_angle"), trkangleDeg); + // ==== Secondary vertex via central DCA vertexer (DCAFitter2) ==== using o2::vertexing::DCAFitter2; constexpr float kBzTesla = 0.5f; DCAFitter2 fitter(/*bz=*/kBzTesla, /*useAbsDCA=*/true, /*propagateToPCA=*/true); fitter.setBz(kBzTesla); - // float bz = o2::base::Propagator::Instance()->getNominalBz(); // en kGauss - // DCAFitter2 fitter(bz, /*useAbsDCA=*/true, /*propagateToPCA=*/true); - fitter.setMaxR(45.f); // cm - fitter.setMaxDZIni(4.f); // cm - fitter.setMaxDXYIni(4.f); // cm + //float bz = o2::base::Propagator::Instance()->getNominalBz(); // en kGauss + //DCAFitter2 fitter(bz, /*useAbsDCA=*/true, /*propagateToPCA=*/true); + fitter.setMaxR(45.f); // cm + fitter.setMaxDZIni(4.f); // cm + fitter.setMaxDXYIni(4.f); // cm fitter.setMaxChi2(50.f); fitter.setPropagateToPCA(true); std::optional pRow, apRow; for (const auto& tr : tracks) { - if (!tr.has_mcParticle()) - continue; + if (!tr.has_mcParticle()) continue; const auto& mc = tr.mcParticle(); - if (mc.globalIndex() == antipId) - apRow = tr; - if (mc.globalIndex() == pId) - pRow = tr; - if (pRow && apRow) - break; + if (mc.globalIndex() == antipId) apRow = tr; + if (mc.globalIndex() == pId) pRow = tr; + if (pRow && apRow) break; } if (!(pRow && apRow)) { - } else { - // TrackParCov - auto trP = makeTPCovFromAOD(*pRow); + } + else { + //TrackParCov + auto trP = makeTPCovFromAOD(*pRow); auto trAP = makeTPCovFromAOD(*apRow); int nCand = fitter.process(trP, trAP); auto status = fitter.getFitStatus(); histos.fill(HIST("vtxfitStatus"), static_cast(status)); if (nCand > 0 && (status == DCAFitter2::FitStatus::Converged || status == DCAFitter2::FitStatus::MaxIter)) { - // Secondary vertex (commom PCA) [x,y,z] cm + //Secondary vertex (commom PCA) [x,y,z] cm auto vtx = fitter.getPCACandidatePos(); const double secX = vtx[0]; const double secY = vtx[1]; const double secZ = vtx[2]; - // DCA of the pair in the PCA (equivalent to minDCA) + //DCA of the pair in the PCA (equivalent to minDCA) fitter.propagateTracksToVertex(); auto tp0 = fitter.getTrackParamAtPCA(0); auto tp1 = fitter.getTrackParamAtPCA(1); @@ -635,15 +591,12 @@ struct NucleiAntineutronCex { const auto p1 = tp1.getXYZGlo(); const double x0 = p0.X(), y0 = p0.Y(), z0 = p0.Z(); const double x1 = p1.X(), y1 = p1.Y(), z1 = p1.Z(); - const double dcaPair = std::sqrt((x0 - x1) * (x0 - x1) + (y0 - y1) * (y0 - y1) + (z0 - z1) * (z0 - z1)); - - if (motherPdg == -kNeutron) - histos.fill(HIST("cex_pairtrkVtxfitDcaPair"), dcaPair); - if (motherPdg != -kNeutron) - histos.fill(HIST("cexbg_pairtrkVtxfitDcaPair"), dcaPair); - - if (!(antipLayers && pLayers)) - continue; + const double dcaPair = std::sqrt((x0-x1)*(x0-x1) + (y0-y1)*(y0-y1) + (z0-z1)*(z0-z1)); + + if (motherPdg == -kNeutron) histos.fill(HIST("cex_pairtrkVtxfitDcaPair"), dcaPair); + if (motherPdg != -kNeutron) histos.fill(HIST("cexbg_pairtrkVtxfitDcaPair"), dcaPair); + + if (!(antipLayers && pLayers)) continue; double cexPairTrkP = total_trk_pVec.Mag(); double cexPairTrkPt = total_trk_pVec.Pt(); double cexPairTrkPz = pTrkPz + antipTrkPz; @@ -651,34 +604,34 @@ struct NucleiAntineutronCex { const double dxPV = secX - pvtxX; const double dyPV = secY - pvtxY; const double dzPV = secZ - pvtxZ; - const double distToPrimary = std::sqrt(dxPV * dxPV + dyPV * dyPV + dzPV * dzPV); - + const double distToPrimary = std::sqrt(dxPV*dxPV + dyPV*dyPV + dzPV*dzPV); + const TVector3 pv2sv(secX - pvtxX, secY - pvtxY, secZ - pvtxZ); const double pairPointingAngleDeg = pv2sv.Angle(total_trk_pVec) * Rad2Deg; - - const double pP = pVecProton_trk.Mag(); - const double pAP = AntipVecProton_trk.Mag(); - const double ptP = pVecProton_trk.Pt(); + + const double pP = pVecProton_trk.Mag(); + const double pAP = AntipVecProton_trk.Mag(); + const double ptP = pVecProton_trk.Pt(); const double ptAP = AntipVecProton_trk.Pt(); - - const double denomP = std::max(1e-9, pP + pAP); + + const double denomP = std::max(1e-9, pP + pAP); const double denomPt = std::max(1e-9, ptP + ptAP); - - const float pairPBalance = std::abs(pP - pAP) / denomP; + + const float pairPBalance = std::abs(pP - pAP) / denomP; const float pairPtBalance = std::abs(ptP - ptAP) / denomPt; - + const float pairQ = (pVecProton_trk - AntipVecProton_trk).Mag(); - - // Trk - MC - const float dPairP = cexPairTrkP - cexPairMcP; - const float dPairPt = cexPairTrkPt - cexPairMcPt; - const float dPairPz = cexPairTrkPz - cexPairMcPz; - const float dOpenAngle = trkangleDeg - mcangleDeg; - + + //Trk - MC + const float dPairP = cexPairTrkP - cexPairMcP; + const float dPairPt = cexPairTrkPt - cexPairMcPt; + const float dPairPz = cexPairTrkPz - cexPairMcPz; + const float dOpenAngle = trkangleDeg - mcangleDeg; + // closest ITS layer: Radius need to be checked - static const std::array Rlayers = {2.2, 2.8, 3.6, 19.6, 24.0, 29.0, 35.0}; + static const std::array Rlayers = {2.2, 2.8, 3.6, 19.6, 24.0, 29.0, 35.0}; int16_t svNearestLayerId = -1; - float SVDeltaRToLayer = 1e9f; + float SVDeltaRToLayer = 1e9f; for (int i = 0; i < static_cast(Rlayers.size()); ++i) { const float dR = static_cast(std::abs(radius - Rlayers[i])); if (dR < SVDeltaRToLayer) { @@ -686,8 +639,8 @@ struct NucleiAntineutronCex { svNearestLayerId = static_cast(i); } } - - if (motherPdg == -kNeutron) { + + if(motherPdg == -kNeutron){ histos.fill(HIST("cexPairTrkP"), cexPairTrkP); histos.fill(HIST("cexPairTrkPt"), cexPairTrkPt); histos.fill(HIST("cexPairTrkPz"), cexPairTrkPz); @@ -695,7 +648,8 @@ struct NucleiAntineutronCex { histos.fill(HIST("cex_pairtrkVtxfitDistToPv"), distToPrimary); histos.fill(HIST("cex_pairtrk_vtxfit_secVtxXY"), secX, secY); histos.fill(HIST("cex_pairtrk_vtxfit_secVtxZ"), secZ); - } else { + } + else{ histos.fill(HIST("cexbg_pairtrk_p"), cexPairTrkP); histos.fill(HIST("cexbg_pairtrk_pt"), cexPairTrkPt); histos.fill(HIST("cexbg_pairtrk_pz"), cexPairTrkPz); @@ -704,99 +658,100 @@ struct NucleiAntineutronCex { histos.fill(HIST("cexbg_pairtrk_vtxfit_secVtxXY"), secX, secY); histos.fill(HIST("cexbg_pairtrk_vtxfit_secVtxZ"), secZ); } - + const float chi2 = fitter.getChi2AtPCACandidate(); histos.fill(HIST("vtxfitChi2"), chi2); const double dX = secX - antipVx; const double dY = secY - antipVy; const double dZ = secZ - antipVz; - const double d3D = std::sqrt(dX * dX + dY * dY + dZ * dZ); + const double d3D = std::sqrt(dX*dX + dY*dY + dZ*dZ); histos.fill(HIST("vtxfit_mc_dX"), dX); histos.fill(HIST("vtxfit_mc_dY"), dY); histos.fill(HIST("vtxfit_mc_dZ"), dZ); histos.fill(HIST("vtxfit_mc_d3D"), d3D); - + const bool isCex = (motherPdg == -kNeutron); - + const float vtxfitDX = secX - antipVx; const float vtxfitDY = secY - antipVy; const float vtxfitDZ = secZ - antipVz; - const float vtxfitD3D = std::sqrt(vtxfitDX * vtxfitDX + vtxfitDY * vtxfitDY + vtxfitDZ * vtxfitDZ); - + const float vtxfitD3D = std::sqrt(vtxfitDX*vtxfitDX + vtxfitDY*vtxfitDY + vtxfitDZ*vtxfitDZ); + const uint32_t selMask = 0u; - + outPairs( - isCex, - motherPdg, - colId, - pId, - antipId, - - cexPairMcP, - cexPairMcPt, - cexPairMcPz, - dplane, - mcangleDeg, - antipVx, - antipVy, - antipVz, - - cexPairTrkP, - cexPairTrkPt, - cexPairTrkPz, - trkangleDeg, - dcaPair, - radius, - distToPrimary, - secX, - secY, - secZ, - - chi2, - static_cast(status), - nCand, - vtxfitDX, - vtxfitDY, - vtxfitDZ, - vtxfitD3D, - - pTrkP, - pTrkPx, - pTrkPy, - pTrkPz, - pTrkEta, - pTrkTpcSignal, - pTrkNClsIts, - - antipTrkP, - antipTrkPx, - antipTrkPy, - antipTrkPz, - antipTrkEta, - antipTrkTpcSignal, - antipTrkNClsIts, - - selMask, - - pairPointingAngleDeg, - pairPBalance, - pairPtBalance, - pairQ, - - dPairP, - dPairPt, - dPairPz, - dOpenAngle, - - svNearestLayerId, - SVDeltaRToLayer, - - pItsMap, - apItsMap, - static_cast(pLayers ? 1 : 0), - static_cast(antipLayers ? 1 : 0), - - pvtxZ); + isCex, + motherPdg, + colId, + pId, + antipId, + + cexPairMcP, + cexPairMcPt, + cexPairMcPz, + dplane, + mcangleDeg, + antipVx, + antipVy, + antipVz, + + cexPairTrkP, + cexPairTrkPt, + cexPairTrkPz, + trkangleDeg, + dcaPair, + radius, + distToPrimary, + secX, + secY, + secZ, + + chi2, + static_cast(status), + nCand, + vtxfitDX, + vtxfitDY, + vtxfitDZ, + vtxfitD3D, + + pTrkP, + pTrkPx, + pTrkPy, + pTrkPz, + pTrkEta, + pTrkTpcSignal, + pTrkNClsIts, + + antipTrkP, + antipTrkPx, + antipTrkPy, + antipTrkPz, + antipTrkEta, + antipTrkTpcSignal, + antipTrkNClsIts, + + selMask, + + pairPointingAngleDeg, + pairPBalance, + pairPtBalance, + pairQ, + + dPairP, + dPairPt, + dPairPz, + dOpenAngle, + + svNearestLayerId, + SVDeltaRToLayer, + + pItsMap, + apItsMap, + static_cast(pLayers ? 1 : 0), + static_cast(antipLayers ? 1 : 0), + + pvtxZ + ); } } // ==== end DCAFitter2 ==== @@ -810,3 +765,4 @@ WorkflowSpec defineDataProcessing(ConfigContext const& ctx) { return WorkflowSpec{adaptAnalysisTask(ctx)}; } + From eaeb9e30d7b9db06ca1a7c1afec5d94b7cb8a1be Mon Sep 17 00:00:00 2001 From: FabiolaLP Date: Sun, 5 Oct 2025 21:11:31 -0600 Subject: [PATCH 4/4] style: apply git-clang-format --- .../Nuspex/nucleiAntineutronCex.cxx | 706 +++++++++--------- 1 file changed, 368 insertions(+), 338 deletions(-) diff --git a/PWGLF/TableProducer/Nuspex/nucleiAntineutronCex.cxx b/PWGLF/TableProducer/Nuspex/nucleiAntineutronCex.cxx index 023f1283b57..9fea0a7e487 100644 --- a/PWGLF/TableProducer/Nuspex/nucleiAntineutronCex.cxx +++ b/PWGLF/TableProducer/Nuspex/nucleiAntineutronCex.cxx @@ -15,25 +15,26 @@ /// \author Fabiola Lugo /// -#include +#include "PWGLF/DataModel/LFAntinCexTables.h" + +#include "CommonConstants/MathConstants.h" +#include "DCAFitter/DCAFitterN.h" +#include "DetectorsBase/Propagator.h" +#include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" -#include "Framework/runDataProcessing.h" #include "Framework/Logger.h" -#include "Framework/AnalysisDataModel.h" -#include "DCAFitter/DCAFitterN.h" +#include "Framework/runDataProcessing.h" #include "ReconstructionDataFormats/TrackParametrization.h" -#include "DetectorsBase/Propagator.h" + #include "TMCProcess.h" -#include -//ROOT -#include "TVector3.h" #include "TMath.h" -#include -#include -#include "CommonConstants/MathConstants.h" #include "TPDGCode.h" -// -#include "PWGLF/DataModel/LFAntinCexTables.h" +#include "TVector3.h" + +#include +#include +#include +#include using namespace o2; using namespace o2::framework; @@ -42,90 +43,91 @@ using o2::constants::math::Rad2Deg; struct NucleiAntineutronCex { // Slicing per colision Preslice perMcByColl = aod::mcparticle::mcCollisionId; - //Check available tables in the AOD, specifically TracksIU, TracksCovIU + // Check available tables in the AOD, specifically TracksIU, TracksCovIU using TracksWCovMc = soa::Join; - + // === Cut values === - static constexpr double kIts2MinR = 2.2; // ITS2 min radius [cm] - static constexpr double kIts2MaxR = 39.0; // ITS2 max radius [cm] - static constexpr double kIts2MaxVz = 39.0; // ITS2 max |vz| [cm] - static constexpr double kAccMaxEta = 1.2; // acceptance |eta| - static constexpr double kAccMaxVz = 5.3; // acceptance |vz| [cm] - static constexpr double kStrictEta = 0.9; // tighter eta cut - static constexpr double kInitDplane = 10.0; // init dplane - static constexpr double kHuge = 1e9; // fallback for bad denom - static constexpr int kMinItsHits = 2; - static constexpr double kVtxTol = 1e-4; - + static constexpr double kIts2MinR = 2.2; // ITS2 min radius [cm] + static constexpr double kIts2MaxR = 39.0; // ITS2 max radius [cm] + static constexpr double kIts2MaxVz = 39.0; // ITS2 max |vz| [cm] + static constexpr double kAccMaxEta = 1.2; // acceptance |eta| + static constexpr double kAccMaxVz = 5.3; // acceptance |vz| [cm] + static constexpr double kStrictEta = 0.9; // tighter eta cut + static constexpr double kInitDplane = 10.0; // init dplane + static constexpr double kHuge = 1e9; // fallback for bad denom + static constexpr int kMinItsHits = 2; + static constexpr double kVtxTol = 1e-4; + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; Produces outPairs; - - void init(InitContext const&) { + + void init(InitContext const&) + { // Primary vertex histos.add("hVx", "Primary vertex X;X (cm);Entries", kTH1F, {{100, -5., 5.}}); histos.add("hVy", "Primary vertex Y;Y (cm);Entries", kTH1F, {{100, -5., 5.}}); histos.add("hVz", "Primary vertex Z;Z (cm);Entries", kTH1F, {{200, -20., 20.}}); - + // Primary antineutrons histos.add("antin_p", "Total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); histos.add("antin_px", "p_{x};p_{x} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); histos.add("antin_py", "p_{y};p_{y} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); histos.add("antin_pz", "p_{z};p_{z} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); - histos.add("antin_eta","Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); + histos.add("antin_eta", "Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); histos.add("antin_p_ITScuts", "Momentum with ITS cuts;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); - + // Primary neutrons - histos.add("n_p", "Total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("n_p", "Total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); histos.add("n_px", "p_{x};p_{x} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); histos.add("n_py", "p_{y};p_{y} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); histos.add("n_pz", "p_{z};p_{z} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); - histos.add("n_eta","Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); + histos.add("n_eta", "Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); histos.add("n_p_ITScuts", "Momentum with ITS cuts;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); - + // Primary antiprotons - histos.add("antipP", "Total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("antipP", "Total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); histos.add("antipPx", "p_{x};p_{x} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); histos.add("antipPy", "p_{y};p_{y} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); histos.add("antipPz", "p_{z};p_{z} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); - histos.add("antipEta","Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); + histos.add("antipEta", "Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); histos.add("antipP_ITScuts", "Momentum with ITS cuts;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); - + // Primary protons histos.add("pP", "Total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); histos.add("pPx", "p_{x};p_{x} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); histos.add("pPy", "p_{y};p_{y} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); histos.add("pPz", "p_{z};p_{z} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); - histos.add("pEta","Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); + histos.add("pEta", "Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); histos.add("pP_ITScuts", "Momentum with ITS cuts;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); - + // test (MC) histos.add("antip_test", "Secondary antiprotons;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); - + // CEX pair from antineutron (MC) histos.add("cexPairMcP", "CEX pair total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); histos.add("cexPairMcPt", "CEX pair p_{T};p_{T} (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); histos.add("cexPairMcPz", "CEX pair p_{z};p_{z} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); - histos.add("cex_pairmcDplane","CEX pair d_{plane};d_{plane} (cm);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("cex_pairmcDplane", "CEX pair d_{plane};d_{plane} (cm);Entries", kTH1F, {{100, 0., 10.}}); histos.add("cex_pairmc_angle", "Pair opening angle;Angle (°);Entries", kTH1F, {{180, 0., 180.}}); histos.add("cex_pairmc_vtx", "MC CEX pair vertex;X (cm);Y (cm)", kTH2F, {{100, -50., 50.}, {100, -50., 50.}}); histos.add("cex_pairmc_vtxz", "MC secondary vertex Z;Z (cm);Entries", kTH1F, {{200, -60., 60.}}); - histos.add("cexPairMcPITScuts","CEX pair momentum (ITS cuts);|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); - + histos.add("cexPairMcPITScuts", "CEX pair momentum (ITS cuts);|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + // CEX pair normalized to antineutron (MC) histos.add("cexn_pairmc_p", "Pair p / antineutron p;p/p_{#bar{n}};Entries", kTH1F, {{100, 0., 2.}}); histos.add("cexn_pairmc_pt", "Pair p_{T} / antineutron p_{T};p_{T}/p_{T,#bar{n}};Entries", kTH1F, {{100, 0., 2.}}); histos.add("cexn_pairmc_pz", "Pair p_{z} / antineutron p_{z};p_{z}/p_{z,#bar{n}};Entries", kTH1F, {{100, -2., 2.}}); - + // BG pair (not from antineutron) (MC) histos.add("cexbg_pairmc_p", "Background pair momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); histos.add("cexbg_pairmc_pt", "Background pair p_{T};p_{T} (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); histos.add("cexbg_pairmc_pz", "Background pair p_{z};p_{z} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); - histos.add("cexbg_pairmcDplane","Background d_{plane};d_{plane} (cm);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("cexbg_pairmcDplane", "Background d_{plane};d_{plane} (cm);Entries", kTH1F, {{100, 0., 10.}}); histos.add("cexbg_pairmc_angle", "Background opening angle;Angle (°);Entries", kTH1F, {{180, 0., 180.}}); histos.add("cexbg_pairmc_vtx", "Background pair vertex;X (cm);Y (cm)", kTH2F, {{100, -50., 50.}, {100, -50., 50.}}); histos.add("cexbg_pairmc_vtxz", "Background secondary vertex Z;Z (cm);Entries", kTH1F, {{200, -60., 60.}}); - histos.add("cexbg_pairmc_pITScuts","Background momentum (ITS cuts);|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); - + histos.add("cexbg_pairmc_pITScuts", "Background momentum (ITS cuts);|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + // CEX pair from antineutron (TRK) histos.add("cex_pairtrk_angle", "Pair opening angle (tracks);Angle (°);Entries", kTH1F, {{180, 0., 180.}}); histos.add("cexPairTrkP", "Pair momentum (tracks);|p| (GeV/c);Entries", kTH1F, {{120, 0., 12.}}); @@ -136,7 +138,7 @@ struct NucleiAntineutronCex { histos.add("cex_pairtrkVtxfitDistToPv", "Distance from secondary vertex to PV;dist (cm);Entries", kTH1F, {{240, 0., 120.}}); histos.add("cex_pairtrk_vtxfit_secVtxXY", "Secondary vertex (PCA);X (cm);Y (cm)", kTH2F, {{200, -60., 60.}, {200, -60., 60.}}); histos.add("cex_pairtrk_vtxfit_secVtxZ", "Secondary vertex Z (PCA);Z (cm);Entries", kTH1F, {{240, -60., 60.}}); - + // BG pair (not from antineutron) (TRK) histos.add("cexbg_pairtrk_angle", "Background opening angle (tracks);Angle (°);Entries", kTH1F, {{180, 0., 180.}}); histos.add("cexbg_pairtrk_p", "Pair momentum (tracks);|p| (GeV/c);Entries", kTH1F, {{120, 0., 12.}}); @@ -147,43 +149,44 @@ struct NucleiAntineutronCex { histos.add("cexbg_pairtrkVtxfitDistToPv", "Distance from secondary vertex to PV;dist (cm);Entries", kTH1F, {{240, 0., 120.}}); histos.add("cexbg_pairtrk_vtxfit_secVtxXY", "Secondary vertex (PCA);X (cm);Y (cm)", kTH2F, {{200, -60., 60.}, {200, -60., 60.}}); histos.add("cexbg_pairtrk_vtxfit_secVtxZ", "Secondary vertex Z (PCA);Z (cm);Entries", kTH1F, {{240, -60., 60.}}); - + // Vertex fit (DCAFitter2 / PCA) histos.add("vtxfitChi2", "DCAFitter2 #chi^{2};#chi^{2};Entries", kTH1F, {{200, 0., 100.}}); histos.add("vtxfitStatus", "Fit status (0=OK);code;Entries", kTH1I, {{10, 0., 10.}}); histos.add("vtxfit_mc_dX", "SV residual X (fit - MC);#Delta X (cm);Entries", kTH1F, {{400, -20., 20.}}); histos.add("vtxfit_mc_dY", "SV residual Y (fit - MC);#Delta Y (cm);Entries", kTH1F, {{400, -20., 20.}}); histos.add("vtxfit_mc_dZ", "SV residual Z (fit - MC);#Delta Z (cm);Entries", kTH1F, {{400, -20., 20.}}); - histos.add("vtxfit_mc_d3D", "SV distance |fit - MC|;#Delta r (cm);Entries", kTH1F, {{300, 0., 30.}}); + histos.add("vtxfit_mc_d3D", "SV distance |fit - MC|;#Delta r (cm);Entries", kTH1F, {{300, 0., 30.}}); } - - static o2::track::TrackParCov makeTPCovFromAOD(const TracksWCovMc::iterator& tr) { + + static o2::track::TrackParCov makeTPCovFromAOD(const TracksWCovMc::iterator& tr) + { using o2::track::TrackParCov; TrackParCov tpcov; // Local state: x, alpha, y, z, snp, tgl, q/pt float par[5] = {tr.y(), tr.z(), tr.snp(), tr.tgl(), tr.signed1Pt()}; tpcov.set(tr.x(), tr.alpha(), par); - + // Covariance matrix (15 terms) in O2 order std::array cov = { - tr.cYY(), tr.cZY(), tr.cZZ(), - tr.cSnpY(), tr.cSnpZ(), tr.cSnpSnp(), - tr.cTglY(), tr.cTglZ(), tr.cTglSnp(), - tr.cTglTgl(), tr.c1PtY(), tr.c1PtZ(), - tr.c1PtSnp(), tr.c1PtTgl(), tr.c1Pt21Pt2() - }; + tr.cYY(), tr.cZY(), tr.cZZ(), + tr.cSnpY(), tr.cSnpZ(), tr.cSnpSnp(), + tr.cTglY(), tr.cTglZ(), tr.cTglSnp(), + tr.cTglTgl(), tr.c1PtY(), tr.c1PtZ(), + tr.c1PtSnp(), tr.c1PtTgl(), tr.c1Pt21Pt2()}; tpcov.setCov(cov); return tpcov; } - - void process(aod::McCollisions const& cols, aod::McParticles const& particles, TracksWCovMc const& tracks) { + + void process(aod::McCollisions const& cols, aod::McParticles const& particles, TracksWCovMc const& tracks) + { double pvtxX = 0; double pvtxY = 0; double pvtxZ = 0; for (auto const& col : cols) { - const auto colId= col.globalIndex(); + const auto colId = col.globalIndex(); auto mcPartsThis = particles.sliceBy(perMcByColl, colId); - + if (std::isfinite(col.posX()) && std::isfinite(col.posY()) && std::isfinite(col.posZ())) { pvtxX = col.posX(); pvtxY = col.posY(); @@ -192,9 +195,9 @@ struct NucleiAntineutronCex { histos.fill(HIST("hVy"), pvtxY); histos.fill(HIST("hVz"), pvtxZ); } - + for (const auto& particle : mcPartsThis) { - + // Primary antineutrons if (particle.pdgCode() == -kNeutron && particle.isPhysicalPrimary()) { histos.fill(HIST("antin_p"), particle.p()); @@ -202,7 +205,8 @@ struct NucleiAntineutronCex { histos.fill(HIST("antin_py"), particle.py()); histos.fill(HIST("antin_pz"), particle.pz()); histos.fill(HIST("antin_eta"), particle.eta()); - if (std::abs(particle.eta()) < kAccMaxEta && std::abs(particle.vz())< kAccMaxVz) histos.fill(HIST("antin_p_ITScuts"),particle.p()); + if (std::abs(particle.eta()) < kAccMaxEta && std::abs(particle.vz()) < kAccMaxVz) + histos.fill(HIST("antin_p_ITScuts"), particle.p()); } // Primary neutrons if (particle.pdgCode() == kNeutron && particle.isPhysicalPrimary()) { @@ -211,7 +215,8 @@ struct NucleiAntineutronCex { histos.fill(HIST("n_py"), particle.py()); histos.fill(HIST("n_pz"), particle.pz()); histos.fill(HIST("n_eta"), particle.eta()); - if (std::abs(particle.eta()) < kAccMaxEta && std::abs(particle.vz())< kAccMaxVz) histos.fill(HIST("n_p_ITScuts"),particle.p()); + if (std::abs(particle.eta()) < kAccMaxEta && std::abs(particle.vz()) < kAccMaxVz) + histos.fill(HIST("n_p_ITScuts"), particle.p()); } // Primary antiprotons if (particle.pdgCode() == -kProton && particle.isPhysicalPrimary()) { @@ -220,7 +225,8 @@ struct NucleiAntineutronCex { histos.fill(HIST("antipPy"), particle.py()); histos.fill(HIST("antipPz"), particle.pz()); histos.fill(HIST("antipEta"), particle.eta()); - if (std::abs(particle.eta()) < kAccMaxEta && std::abs(particle.vz())< kAccMaxVz) histos.fill(HIST("antipP_ITScuts"),particle.p()); + if (std::abs(particle.eta()) < kAccMaxEta && std::abs(particle.vz()) < kAccMaxVz) + histos.fill(HIST("antipP_ITScuts"), particle.p()); } // Primary protons if (particle.pdgCode() == kProton && particle.isPhysicalPrimary()) { @@ -229,38 +235,41 @@ struct NucleiAntineutronCex { histos.fill(HIST("pPy"), particle.py()); histos.fill(HIST("pPz"), particle.pz()); histos.fill(HIST("pEta"), particle.eta()); - if (std::abs(particle.eta()) < kAccMaxEta && std::abs(particle.vz())< kAccMaxVz) histos.fill(HIST("pP_ITScuts"),particle.p()); + if (std::abs(particle.eta()) < kAccMaxEta && std::abs(particle.vz()) < kAccMaxVz) + histos.fill(HIST("pP_ITScuts"), particle.p()); } - - //Seconday antiprotons from material + + // Seconday antiprotons from material const auto procEnum = particle.getProcess(); const bool isSecondaryFromMaterial = (!particle.producedByGenerator()) && (procEnum == kPHadronic || procEnum == kPHInhelastic); - if (particle.pdgCode() != -kProton || !isSecondaryFromMaterial || particle.mothersIds().empty()) continue; + if (particle.pdgCode() != -kProton || !isSecondaryFromMaterial || particle.mothersIds().empty()) + continue; histos.fill(HIST("antip_test"), particle.p()); - - //Primary mother + + // Primary mother bool hasPrimaryMotherAntip = false; double motherPt = 0.0; double motherPz = 0.0; double motherVz = 0.0; - double motherP = 0.0; + double motherP = 0.0; double motherEta = 0.0; - int motherPdg = 0; - + int motherPdg = 0; + for (const auto& mother : particle.mothers_as()) { - if (mother.isPhysicalPrimary()){ + if (mother.isPhysicalPrimary()) { hasPrimaryMotherAntip = true; motherPt = mother.pt(); motherPz = mother.pz(); motherVz = mother.vz(); - motherP= mother.p(); + motherP = mother.p(); motherEta = mother.eta(); motherPdg = mother.pdgCode(); break; } } - if (!hasPrimaryMotherAntip) continue; - + if (!hasPrimaryMotherAntip) + continue; + double antipVx = particle.vx(); double antipVy = particle.vy(); double antipVz = particle.vz(); @@ -269,22 +278,26 @@ struct NucleiAntineutronCex { double antipPz = particle.pz(); double antipE = particle.e(); int antipId = particle.globalIndex(); - - //Selection conditions: Produced in the ITS - const double R = std::sqrt(antipVx*antipVx + antipVy*antipVy); - //Config for ITS - //if(3.9<=R && R<=43.0 && std::abs(antipVz)<=48.9){ - //Config for ITS2 - if (R < kIts2MinR || R > kIts2MaxR || std::abs(antipVz) > kIts2MaxVz) continue; - if (std::abs(motherEta) >= kAccMaxEta || std::abs(motherVz) >= kAccMaxVz) continue; - - //Pion minus veto + + // Selection conditions: Produced in the ITS + const double r = std::sqrt(antipVx * antipVx + antipVy * antipVy); + // Config for ITS + // if(3.9<=r && r<=43.0 && std::abs(antipVz)<=48.9){ + // Config for ITS2 + if (r < kIts2MinR || r > kIts2MaxR || std::abs(antipVz) > kIts2MaxVz) + continue; + if (std::abs(motherEta) >= kAccMaxEta || std::abs(motherVz) >= kAccMaxVz) + continue; + + // Pion minus veto bool pionMinus = false; for (const auto& particle1 : mcPartsThis) { const auto proc1Enum = particle1.getProcess(); const bool isSecondaryFromMaterial1 = (!particle1.producedByGenerator()) && (proc1Enum == kPHadronic || proc1Enum == kPHInhelastic); - if (particle1.mcCollisionId() != colId) continue; - if (particle1.pdgCode() != kPiMinus || !isSecondaryFromMaterial1 || particle1.mothersIds().empty()) continue; + if (particle1.mcCollisionId() != colId) + continue; + if (particle1.pdgCode() != kPiMinus || !isSecondaryFromMaterial1 || particle1.mothersIds().empty()) + continue; bool hasPrimaryMotherPim = false; for (const auto& mother : particle1.mothers_as()) { if (mother.isPhysicalPrimary()) { @@ -292,23 +305,26 @@ struct NucleiAntineutronCex { break; } } - if (!hasPrimaryMotherPim) continue; + if (!hasPrimaryMotherPim) + continue; double pimVx = particle1.vx(); double pimVy = particle1.vy(); double pimVz = particle1.vz(); - if(std::abs(pimVx - antipVx) < kVtxTol && std::abs(pimVy - antipVy) < kVtxTol && std::abs(pimVz - antipVz) < kVtxTol){ + if (std::abs(pimVx - antipVx) < kVtxTol && std::abs(pimVy - antipVy) < kVtxTol && std::abs(pimVz - antipVz) < kVtxTol) { pionMinus = true; break; } } - - //Pion plus veto + + // Pion plus veto bool pionPlus = false; for (const auto& particle2 : mcPartsThis) { - if (particle2.mcCollisionId() != colId) continue; + if (particle2.mcCollisionId() != colId) + continue; const auto proc2Enum = particle2.getProcess(); const bool isSecondaryFromMaterial2 = (!particle2.producedByGenerator()) && (proc2Enum == kPHadronic || proc2Enum == kPHInhelastic); - if (particle2.pdgCode() != kPiPlus || !isSecondaryFromMaterial2 || particle2.mothersIds().empty()) continue; + if (particle2.pdgCode() != kPiPlus || !isSecondaryFromMaterial2 || particle2.mothersIds().empty()) + continue; bool hasPrimaryMotherPip = false; for (const auto& mother : particle2.mothers_as()) { if (mother.isPhysicalPrimary()) { @@ -316,18 +332,20 @@ struct NucleiAntineutronCex { break; } } - if (!hasPrimaryMotherPip) continue; + if (!hasPrimaryMotherPip) + continue; double pipVx = particle2.vx(); double pipVy = particle2.vy(); double pipVz = particle2.vz(); - if(std::abs(pipVx - antipVx) < kVtxTol && std::abs(pipVy - antipVy) < kVtxTol && std::abs(pipVz - antipVz) < kVtxTol){ + if (std::abs(pipVx - antipVx) < kVtxTol && std::abs(pipVy - antipVy) < kVtxTol && std::abs(pipVz - antipVz) < kVtxTol) { pionPlus = true; break; } } - - if (pionPlus || pionMinus) continue; - + + if (pionPlus || pionMinus) + continue; + // CEX selection double dplane = kInitDplane; double dplaneTmp = 0; @@ -341,33 +359,36 @@ struct NucleiAntineutronCex { int k_plane = -1; int k_e = -1; int k_p = -1; - - //Secondary proton from material + + // Secondary proton from material for (const auto& particle3 : mcPartsThis) { - if (particle3.mcCollisionId() != colId) continue; + if (particle3.mcCollisionId() != colId) + continue; const auto proc3Enum = particle3.getProcess(); const bool isSecondaryFromMaterial3 = (!particle3.producedByGenerator()) && (proc3Enum == kPHadronic || proc3Enum == kPHInhelastic); - if (particle3.pdgCode() != kProton || !isSecondaryFromMaterial3 || particle3.mothersIds().empty()) continue; - bool hasPrimaryMotherP= false; + if (particle3.pdgCode() != kProton || !isSecondaryFromMaterial3 || particle3.mothersIds().empty()) + continue; + bool hasPrimaryMotherP = false; for (const auto& mother : particle3.mothers_as()) { if (mother.isPhysicalPrimary()) { - hasPrimaryMotherP= true; + hasPrimaryMotherP = true; break; } } - if (!hasPrimaryMotherP) continue; - double pVx = particle3.vx(); - double pVy = particle3.vy(); - double pVz = particle3.vz(); + if (!hasPrimaryMotherP) + continue; + double protonVx = particle3.vx(); + double protonVy = particle3.vy(); + double protonVz = particle3.vz(); double pPx = particle3.px(); double pPy = particle3.py(); double pPz = particle3.pz(); double pE = particle3.e(); - if(std::abs(pVx - antipVx) < kVtxTol && std::abs(pVy - antipVy) < kVtxTol && std::abs(pVz - antipVz) < kVtxTol){ - //Same mother + if (std::abs(protonVx - antipVx) < kVtxTol && std::abs(protonVy - antipVy) < kVtxTol && std::abs(protonVz - antipVz) < kVtxTol) { + // Same mother bool shareMother = false; - const auto& momsAp = particle.mothersIds(); //antiproton - const auto& momsP = particle3.mothersIds(); //proton + const auto& momsAp = particle.mothersIds(); // antiproton + const auto& momsP = particle3.mothersIds(); // proton for (const auto& ida : momsAp) { for (const auto& idp : momsP) { if (ida == idp) { @@ -375,66 +396,62 @@ struct NucleiAntineutronCex { break; } } - if (shareMother) break; + if (shareMother) + break; } - if (!shareMother) continue; - - //CEX proton selection - //dplaneTmp = (pPy*antipPz - pPz*antipPy)*(pvtxX-antipVx) + (pPz*antipPx - pPx*antipPz)*(pvtxY-antipVy) + (pPx*antipPy - pPy*antipPx)*(pvtxZ-antipVz); - double nx = (pPy*antipPz - pPz*antipPy); - double ny = (pPz*antipPx - pPx*antipPz); - double nz = (pPx*antipPy - pPy*antipPx); + if (!shareMother) + continue; + + // CEX proton selection + // dplaneTmp = (pPy*antipPz - pPz*antipPy)*(pvtxX-antipVx) + (pPz*antipPx - pPx*antipPz)*(pvtxY-antipVy) + (pPx*antipPy - pPy*antipPx)*(pvtxZ-antipVz); + double nx = (pPy * antipPz - pPz * antipPy); + double ny = (pPz * antipPx - pPx * antipPz); + double nz = (pPx * antipPy - pPy * antipPx); double rx = (pvtxX - antipVx); double ry = (pvtxY - antipVy); double rz = (pvtxZ - antipVz); - double denom = nx*nx + ny*ny + nz*nz; - if (denom > 0.) - { - dplaneTmp = std::abs(nx*rx + ny*ry + nz*rz) / std::sqrt(denom); - } - else - { + double denom = nx * nx + ny * ny + nz * nz; + if (denom > 0.) { + dplaneTmp = std::abs(nx * rx + ny * ry + nz * rz) / std::sqrt(denom); + } else { dplaneTmp = kHuge; } - if(std::abs(dplaneTmp) < std::abs(dplane)) - { + if (std::abs(dplaneTmp) < std::abs(dplane)) { k_plane = particle3.globalIndex(); dplane = dplaneTmp; } - + eTmp = antipE + pE; - if(std::abs(eTmp) > std::abs(e)) - { + if (std::abs(eTmp) > std::abs(e)) { k_e = particle3.globalIndex(); e = eTmp; } - - pTmp = std::sqrt(std::pow((pPx+antipPx),2)+std::pow((pPy+antipPy),2)+std::pow((pPz+antipPz),2)); - if(std::abs(pTmp) > std::abs(p)) - { + + pTmp = std::sqrt(std::pow((pPx + antipPx), 2) + std::pow((pPy + antipPy), 2) + std::pow((pPz + antipPz), 2)); + if (std::abs(pTmp) > std::abs(p)) { k_p = particle3.globalIndex(); p = pTmp; - pcexPx = pPx; - pcexPy = pPy; - pcexPz = pPz; + pcexPx = pPx; + pcexPy = pPy; + pcexPz = pPz; } } } - - if(k_plane == k_e && k_plane == k_p && k_plane >= 0 ){ + + if (k_plane == k_e && k_plane == k_p && k_plane >= 0) { int pId = k_plane; TVector3 pVecProton = TVector3(pcexPx, pcexPy, pcexPz); TVector3 pVecAntiproton = TVector3(antipPx, antipPy, antipPz); TVector3 total_mc_pVec = pVecProton + pVecAntiproton; double cexPairMcP = total_mc_pVec.Mag(); double cexPairMcPt = total_mc_pVec.Pt(); - double cexPairMcPz = pcexPz+antipPz; + double cexPairMcPz = pcexPz + antipPz; double mcangleRad = pVecProton.Angle(pVecAntiproton); double mcangleDeg = mcangleRad * Rad2Deg; - - //Antineutron mother + + // Antineutron mother if (motherPdg == -kNeutron) { - //CEX pair + // CEX pair histos.fill(HIST("cexPairMcP"), cexPairMcP); histos.fill(HIST("cexPairMcPt"), cexPairMcPt); histos.fill(HIST("cexPairMcPz"), cexPairMcPz); @@ -442,15 +459,19 @@ struct NucleiAntineutronCex { histos.fill(HIST("cex_pairmc_angle"), mcangleDeg); histos.fill(HIST("cex_pairmc_vtx"), antipVx, antipVy); histos.fill(HIST("cex_pairmc_vtxz"), antipVz); - if (std::abs(motherEta) < kStrictEta && std::abs(motherVz)< kAccMaxVz ) histos.fill(HIST("cexPairMcPITScuts"),cexPairMcP); - //CEX pair normalized - if (motherP != 0) histos.fill(HIST("cexn_pairmc_p"), cexPairMcP /motherP); - if (motherPt != 0) histos.fill(HIST("cexn_pairmc_pt"), cexPairMcPt/motherPt); - if (motherPz != 0) histos.fill(HIST("cexn_pairmc_pz"), cexPairMcPz/motherPz); + if (std::abs(motherEta) < kStrictEta && std::abs(motherVz) < kAccMaxVz) + histos.fill(HIST("cexPairMcPITScuts"), cexPairMcP); + // CEX pair normalized + if (motherP != 0) + histos.fill(HIST("cexn_pairmc_p"), cexPairMcP / motherP); + if (motherPt != 0) + histos.fill(HIST("cexn_pairmc_pt"), cexPairMcPt / motherPt); + if (motherPz != 0) + histos.fill(HIST("cexn_pairmc_pz"), cexPairMcPz / motherPz); } - //BG mother + // BG mother if (motherPdg != -kNeutron) { - //CEX pair + // CEX pair histos.fill(HIST("cexbg_pairmc_p"), cexPairMcP); histos.fill(HIST("cexbg_pairmc_pt"), cexPairMcPt); histos.fill(HIST("cexbg_pairmc_pz"), cexPairMcPz); @@ -458,132 +479,141 @@ struct NucleiAntineutronCex { histos.fill(HIST("cexbg_pairmc_angle"), mcangleDeg); histos.fill(HIST("cexbg_pairmc_vtx"), antipVx, antipVy); histos.fill(HIST("cexbg_pairmc_vtxz"), antipVz); - if (std::abs(motherEta) < kStrictEta && std::abs(motherVz)< kAccMaxVz) histos.fill(HIST("cexbg_pairmc_pITScuts"),cexPairMcP); + if (std::abs(motherEta) < kStrictEta && std::abs(motherVz) < kAccMaxVz) + histos.fill(HIST("cexbg_pairmc_pITScuts"), cexPairMcP); } - - //Detector signal + + // Detector signal bool antipLayers = false; bool antipHasTrack = false; double antipTrkPx = 0.; double antipTrkPy = 0.; double antipTrkPz = 0.; - double antipTrkP= 0.; + double antipTrkP = 0.; double antipTrkEta = 0.; double antipTrkTpcSignal = 0; - //int antip_trk_nClsTPC = 0; + // int antip_trk_nClsTPC = 0; int antipTrkNClsIts = 0; - uint16_t apItsMap = 0; - + uint16_t apItsMap = 0; + bool pLayers = false; bool pHasTrack = false; double pTrkPx = 0.; double pTrkPy = 0.; double pTrkPz = 0.; - double pTrkP= 0.; + double pTrkP = 0.; double pTrkEta = 0.; double pTrkTpcSignal = 0; - //int p_trk_nClsTPC = 0; + // int p_trk_nClsTPC = 0; int pTrkNClsIts = 0; - uint16_t pItsMap = 0; - + uint16_t pItsMap = 0; + for (const auto& track : tracks) { - if (!track.has_mcParticle()) continue; + if (!track.has_mcParticle()) + continue; const auto& mc = track.mcParticle(); - if (mc.mcCollisionId() != colId) continue; + if (mc.mcCollisionId() != colId) + continue; uint8_t itsMap = track.itsClusterMap(); - //Config for ITS1 + // Config for ITS1 /*bool hitSPD = (itsMap & 0x3) != 0; // bits 0 (SPD L1) & 1 (SPD L2) bool hitSDD = (itsMap & 0xC) != 0; // bits 2–3 bool hitSSD = (itsMap & 0x30) != 0; // bits 4–5 bool layerCondition = (hitSDD || hitSSD) && !hitSPD;*/ - //Config for ITS2 - bool hitL0 = (itsMap & (1u<<0)) != 0; - bool hitL1 = (itsMap & (1u<<1)) != 0; - bool hitL2 = (itsMap & (1u<<2)) != 0; - bool hitL3 = (itsMap & (1u<<3)) != 0; - bool hitL4 = (itsMap & (1u<<4)) != 0; - bool hitL5 = (itsMap & (1u<<5)) != 0; - bool hitL6 = (itsMap & (1u<<6)) != 0; + // Config for ITS2 + bool hitL0 = (itsMap & (1u << 0)) != 0; + bool hitL1 = (itsMap & (1u << 1)) != 0; + bool hitL2 = (itsMap & (1u << 2)) != 0; + bool hitL3 = (itsMap & (1u << 3)) != 0; + bool hitL4 = (itsMap & (1u << 4)) != 0; + bool hitL5 = (itsMap & (1u << 5)) != 0; + bool hitL6 = (itsMap & (1u << 6)) != 0; bool hitIB = (hitL0 || hitL1 || hitL2); bool hitOuter = (hitL3 || hitL4 || hitL5 || hitL6); int nITS = track.itsNCls(); bool layerCondition = (!hitIB) && hitOuter && (nITS >= kMinItsHits); - + if (mc.globalIndex() == antipId) { - antipTrkP= track.p(); + antipTrkP = track.p(); antipTrkPx = track.px(); antipTrkPy = track.py(); antipTrkPz = track.pz(); antipTrkEta = track.eta(); antipTrkTpcSignal = track.tpcSignal(); - //antip_trk_nClsTPC = track.tpcNCls(); + // antip_trk_nClsTPC = track.tpcNCls(); antipTrkNClsIts = track.itsNCls(); antipHasTrack = true; - apItsMap = static_cast(track.itsClusterMap()); + apItsMap = static_cast(track.itsClusterMap()); antipLayers = (apItsMap != 0); - if (layerCondition) antipLayers = true; - } - else if (mc.globalIndex() == pId) { - pTrkP= track.p(); + if (layerCondition) + antipLayers = true; + } else if (mc.globalIndex() == pId) { + pTrkP = track.p(); pTrkPx = track.px(); pTrkPy = track.py(); pTrkPz = track.pz(); pTrkEta = track.eta(); pTrkTpcSignal = track.tpcSignal(); - //p_trk_nClsTPC = track.tpcNCls(); + // p_trk_nClsTPC = track.tpcNCls(); pTrkNClsIts = track.itsNCls(); pHasTrack = true; - pItsMap = static_cast(track.itsClusterMap()); - pLayers = (pItsMap != 0); - if (layerCondition) pLayers = true; + pItsMap = static_cast(track.itsClusterMap()); + pLayers = (pItsMap != 0); + if (layerCondition) + pLayers = true; } } - if (!(pHasTrack && antipHasTrack)) continue; - + if (!(pHasTrack && antipHasTrack)) + continue; + TVector3 pVecProton_trk(pTrkPx, pTrkPy, pTrkPz); TVector3 AntipVecProton_trk(antipTrkPx, antipTrkPy, antipTrkPz); TVector3 total_trk_pVec = pVecProton_trk + AntipVecProton_trk; double trkangleRad = AntipVecProton_trk.Angle(pVecProton_trk); double trkangleDeg = trkangleRad * Rad2Deg; - if(motherPdg == -kNeutron) histos.fill(HIST("cex_pairtrk_angle"), trkangleDeg); - if(motherPdg != -kNeutron) histos.fill(HIST("cexbg_pairtrk_angle"), trkangleDeg); - + if (motherPdg == -kNeutron) + histos.fill(HIST("cex_pairtrk_angle"), trkangleDeg); + if (motherPdg != -kNeutron) + histos.fill(HIST("cexbg_pairtrk_angle"), trkangleDeg); + // ==== Secondary vertex via central DCA vertexer (DCAFitter2) ==== using o2::vertexing::DCAFitter2; constexpr float kBzTesla = 0.5f; DCAFitter2 fitter(/*bz=*/kBzTesla, /*useAbsDCA=*/true, /*propagateToPCA=*/true); fitter.setBz(kBzTesla); - //float bz = o2::base::Propagator::Instance()->getNominalBz(); // en kGauss - //DCAFitter2 fitter(bz, /*useAbsDCA=*/true, /*propagateToPCA=*/true); - fitter.setMaxR(45.f); // cm - fitter.setMaxDZIni(4.f); // cm - fitter.setMaxDXYIni(4.f); // cm + // float bz = o2::base::Propagator::Instance()->getNominalBz(); // en kGauss + // DCAFitter2 fitter(bz, /*useAbsDCA=*/true, /*propagateToPCA=*/true); + fitter.setMaxR(45.f); // cm + fitter.setMaxDZIni(4.f); // cm + fitter.setMaxDXYIni(4.f); // cm fitter.setMaxChi2(50.f); fitter.setPropagateToPCA(true); std::optional pRow, apRow; for (const auto& tr : tracks) { - if (!tr.has_mcParticle()) continue; + if (!tr.has_mcParticle()) + continue; const auto& mc = tr.mcParticle(); - if (mc.globalIndex() == antipId) apRow = tr; - if (mc.globalIndex() == pId) pRow = tr; - if (pRow && apRow) break; - } - if (!(pRow && apRow)) { + if (mc.globalIndex() == antipId) + apRow = tr; + if (mc.globalIndex() == pId) + pRow = tr; + if (pRow && apRow) + break; } - else { - //TrackParCov - auto trP = makeTPCovFromAOD(*pRow); + if (pRow && apRow) { + // TrackParCov + auto trP = makeTPCovFromAOD(*pRow); auto trAP = makeTPCovFromAOD(*apRow); int nCand = fitter.process(trP, trAP); auto status = fitter.getFitStatus(); histos.fill(HIST("vtxfitStatus"), static_cast(status)); if (nCand > 0 && (status == DCAFitter2::FitStatus::Converged || status == DCAFitter2::FitStatus::MaxIter)) { - //Secondary vertex (commom PCA) [x,y,z] cm + // Secondary vertex (commom PCA) [x,y,z] cm auto vtx = fitter.getPCACandidatePos(); const double secX = vtx[0]; const double secY = vtx[1]; const double secZ = vtx[2]; - //DCA of the pair in the PCA (equivalent to minDCA) + // DCA of the pair in the PCA (equivalent to minDCA) fitter.propagateTracksToVertex(); auto tp0 = fitter.getTrackParamAtPCA(0); auto tp1 = fitter.getTrackParamAtPCA(1); @@ -591,56 +621,59 @@ struct NucleiAntineutronCex { const auto p1 = tp1.getXYZGlo(); const double x0 = p0.X(), y0 = p0.Y(), z0 = p0.Z(); const double x1 = p1.X(), y1 = p1.Y(), z1 = p1.Z(); - const double dcaPair = std::sqrt((x0-x1)*(x0-x1) + (y0-y1)*(y0-y1) + (z0-z1)*(z0-z1)); - - if (motherPdg == -kNeutron) histos.fill(HIST("cex_pairtrkVtxfitDcaPair"), dcaPair); - if (motherPdg != -kNeutron) histos.fill(HIST("cexbg_pairtrkVtxfitDcaPair"), dcaPair); - - if (!(antipLayers && pLayers)) continue; + const double dcaPair = std::sqrt((x0 - x1) * (x0 - x1) + (y0 - y1) * (y0 - y1) + (z0 - z1) * (z0 - z1)); + + if (motherPdg == -kNeutron) + histos.fill(HIST("cex_pairtrkVtxfitDcaPair"), dcaPair); + if (motherPdg != -kNeutron) + histos.fill(HIST("cexbg_pairtrkVtxfitDcaPair"), dcaPair); + + if (!(antipLayers && pLayers)) + continue; double cexPairTrkP = total_trk_pVec.Mag(); double cexPairTrkPt = total_trk_pVec.Pt(); double cexPairTrkPz = pTrkPz + antipTrkPz; const double radius = std::hypot(secX, secY); - const double dxPV = secX - pvtxX; - const double dyPV = secY - pvtxY; - const double dzPV = secZ - pvtxZ; - const double distToPrimary = std::sqrt(dxPV*dxPV + dyPV*dyPV + dzPV*dzPV); - + const double dxPv = secX - pvtxX; + const double dyPv = secY - pvtxY; + const double dzPv = secZ - pvtxZ; + const double distToPrimary = std::sqrt(dxPv * dxPv + dyPv * dyPv + dzPv * dzPv); + const TVector3 pv2sv(secX - pvtxX, secY - pvtxY, secZ - pvtxZ); const double pairPointingAngleDeg = pv2sv.Angle(total_trk_pVec) * Rad2Deg; - - const double pP = pVecProton_trk.Mag(); - const double pAP = AntipVecProton_trk.Mag(); - const double ptP = pVecProton_trk.Pt(); + + const double pP = pVecProton_trk.Mag(); + const double pAP = AntipVecProton_trk.Mag(); + const double ptP = pVecProton_trk.Pt(); const double ptAP = AntipVecProton_trk.Pt(); - - const double denomP = std::max(1e-9, pP + pAP); + + const double denomP = std::max(1e-9, pP + pAP); const double denomPt = std::max(1e-9, ptP + ptAP); - - const float pairPBalance = std::abs(pP - pAP) / denomP; + + const float pairPBalance = std::abs(pP - pAP) / denomP; const float pairPtBalance = std::abs(ptP - ptAP) / denomPt; - + const float pairQ = (pVecProton_trk - AntipVecProton_trk).Mag(); - - //Trk - MC - const float dPairP = cexPairTrkP - cexPairMcP; - const float dPairPt = cexPairTrkPt - cexPairMcPt; - const float dPairPz = cexPairTrkPz - cexPairMcPz; - const float dOpenAngle = trkangleDeg - mcangleDeg; - - // closest ITS layer: Radius need to be checked - static const std::array Rlayers = {2.2, 2.8, 3.6, 19.6, 24.0, 29.0, 35.0}; + + // Trk - MC + const float dPairP = cexPairTrkP - cexPairMcP; + const float dPairPt = cexPairTrkPt - cexPairMcPt; + const float dPairPz = cexPairTrkPz - cexPairMcPz; + const float dOpenAngle = trkangleDeg - mcangleDeg; + + // Closest ITS layer: Radius need to be checked + static const std::array rLayers = {2.2, 2.8, 3.6, 19.6, 24.0, 29.0, 35.0}; int16_t svNearestLayerId = -1; - float SVDeltaRToLayer = 1e9f; - for (int i = 0; i < static_cast(Rlayers.size()); ++i) { - const float dR = static_cast(std::abs(radius - Rlayers[i])); - if (dR < SVDeltaRToLayer) { - SVDeltaRToLayer = dR; + float svDeltaRToLayer = 1e9f; + for (int i = 0; i < static_cast(rLayers.size()); ++i) { + const float dR = static_cast(std::abs(radius - rLayers[i])); + if (dR < svDeltaRToLayer) { + svDeltaRToLayer = dR; svNearestLayerId = static_cast(i); } } - - if(motherPdg == -kNeutron){ + + if (motherPdg == -kNeutron) { histos.fill(HIST("cexPairTrkP"), cexPairTrkP); histos.fill(HIST("cexPairTrkPt"), cexPairTrkPt); histos.fill(HIST("cexPairTrkPz"), cexPairTrkPz); @@ -648,8 +681,7 @@ struct NucleiAntineutronCex { histos.fill(HIST("cex_pairtrkVtxfitDistToPv"), distToPrimary); histos.fill(HIST("cex_pairtrk_vtxfit_secVtxXY"), secX, secY); histos.fill(HIST("cex_pairtrk_vtxfit_secVtxZ"), secZ); - } - else{ + } else { histos.fill(HIST("cexbg_pairtrk_p"), cexPairTrkP); histos.fill(HIST("cexbg_pairtrk_pt"), cexPairTrkPt); histos.fill(HIST("cexbg_pairtrk_pz"), cexPairTrkPz); @@ -658,100 +690,99 @@ struct NucleiAntineutronCex { histos.fill(HIST("cexbg_pairtrk_vtxfit_secVtxXY"), secX, secY); histos.fill(HIST("cexbg_pairtrk_vtxfit_secVtxZ"), secZ); } - + const float chi2 = fitter.getChi2AtPCACandidate(); histos.fill(HIST("vtxfitChi2"), chi2); - const double dX = secX - antipVx; - const double dY = secY - antipVy; - const double dZ = secZ - antipVz; - const double d3D = std::sqrt(dX*dX + dY*dY + dZ*dZ); - histos.fill(HIST("vtxfit_mc_dX"), dX); - histos.fill(HIST("vtxfit_mc_dY"), dY); - histos.fill(HIST("vtxfit_mc_dZ"), dZ); - histos.fill(HIST("vtxfit_mc_d3D"), d3D); - + const double dx = secX - antipVx; + const double dy = secY - antipVy; + const double dz = secZ - antipVz; + const double d3d = std::sqrt(dx * dx + dy * dy + dz * dz); + histos.fill(HIST("vtxfit_mc_dX"), dx); + histos.fill(HIST("vtxfit_mc_dY"), dy); + histos.fill(HIST("vtxfit_mc_dZ"), dz); + histos.fill(HIST("vtxfit_mc_d3D"), d3d); + const bool isCex = (motherPdg == -kNeutron); - + const float vtxfitDX = secX - antipVx; const float vtxfitDY = secY - antipVy; const float vtxfitDZ = secZ - antipVz; - const float vtxfitD3D = std::sqrt(vtxfitDX*vtxfitDX + vtxfitDY*vtxfitDY + vtxfitDZ*vtxfitDZ); - + const float vtxfitD3D = std::sqrt(vtxfitDX * vtxfitDX + vtxfitDY * vtxfitDY + vtxfitDZ * vtxfitDZ); + const uint32_t selMask = 0u; - + outPairs( - isCex, - motherPdg, - colId, - pId, - antipId, - - cexPairMcP, - cexPairMcPt, - cexPairMcPz, - dplane, - mcangleDeg, - antipVx, - antipVy, - antipVz, - - cexPairTrkP, - cexPairTrkPt, - cexPairTrkPz, - trkangleDeg, - dcaPair, - radius, - distToPrimary, - secX, - secY, - secZ, - - chi2, - static_cast(status), - nCand, - vtxfitDX, - vtxfitDY, - vtxfitDZ, - vtxfitD3D, - - pTrkP, - pTrkPx, - pTrkPy, - pTrkPz, - pTrkEta, - pTrkTpcSignal, - pTrkNClsIts, - - antipTrkP, - antipTrkPx, - antipTrkPy, - antipTrkPz, - antipTrkEta, - antipTrkTpcSignal, - antipTrkNClsIts, - - selMask, - - pairPointingAngleDeg, - pairPBalance, - pairPtBalance, - pairQ, - - dPairP, - dPairPt, - dPairPz, - dOpenAngle, - - svNearestLayerId, - SVDeltaRToLayer, - - pItsMap, - apItsMap, - static_cast(pLayers ? 1 : 0), - static_cast(antipLayers ? 1 : 0), - - pvtxZ - ); + isCex, + motherPdg, + colId, + pId, + antipId, + + cexPairMcP, + cexPairMcPt, + cexPairMcPz, + dplane, + mcangleDeg, + antipVx, + antipVy, + antipVz, + + cexPairTrkP, + cexPairTrkPt, + cexPairTrkPz, + trkangleDeg, + dcaPair, + radius, + distToPrimary, + secX, + secY, + secZ, + + chi2, + static_cast(status), + nCand, + vtxfitDX, + vtxfitDY, + vtxfitDZ, + vtxfitD3D, + + pTrkP, + pTrkPx, + pTrkPy, + pTrkPz, + pTrkEta, + pTrkTpcSignal, + pTrkNClsIts, + + antipTrkP, + antipTrkPx, + antipTrkPy, + antipTrkPz, + antipTrkEta, + antipTrkTpcSignal, + antipTrkNClsIts, + + selMask, + + pairPointingAngleDeg, + pairPBalance, + pairPtBalance, + pairQ, + + dPairP, + dPairPt, + dPairPz, + dOpenAngle, + + svNearestLayerId, + svDeltaRToLayer, + + pItsMap, + apItsMap, + static_cast(pLayers ? 1 : 0), + static_cast(antipLayers ? 1 : 0), + + pvtxZ); } } // ==== end DCAFitter2 ==== @@ -765,4 +796,3 @@ WorkflowSpec defineDataProcessing(ConfigContext const& ctx) { return WorkflowSpec{adaptAnalysisTask(ctx)}; } -