diff --git a/PWGLF/DataModel/LFAntinCexTables.h b/PWGLF/DataModel/LFAntinCexTables.h index de11af229c7..9e8f4ce1ca8 100644 --- a/PWGLF/DataModel/LFAntinCexTables.h +++ b/PWGLF/DataModel/LFAntinCexTables.h @@ -18,101 +18,101 @@ #ifndef PWGLF_DATAMODEL_LFANTINCEXTABLES_H_ #define PWGLF_DATAMODEL_LFANTINCEXTABLES_H_ -#include "Framework/AnalysisDataModel.h" #include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" namespace o2::aod { namespace antin_cex { // 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 +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); +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); +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); +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); +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(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); +DECLARE_SOA_COLUMN(AntipTrkNClsIts, antipTrkNClsIts, int16_t); -//Cuts Mask -DECLARE_SOA_COLUMN(SelMask, selMask, uint32_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(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(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(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(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); +DECLARE_SOA_COLUMN(PVtxZ, pVtxZ, float); // Proton ITS PID -DECLARE_SOA_COLUMN(PTrkItsNSigmaPr, pTrkItsNSigmaPr, float); -DECLARE_SOA_COLUMN(PTrkItsPidValid, pTrkItsPidValid, int8_t); -DECLARE_SOA_COLUMN(PTrkTgl, pTrkTgl, float); +DECLARE_SOA_COLUMN(PTrkItsNSigmaPr, pTrkItsNSigmaPr, float); +DECLARE_SOA_COLUMN(PTrkItsPidValid, pTrkItsPidValid, int8_t); +DECLARE_SOA_COLUMN(PTrkTgl, pTrkTgl, float); // Antiproton ITS PID -DECLARE_SOA_COLUMN(AntipTrkItsNSigmaPr, antipTrkItsNSigmaPr, float); -DECLARE_SOA_COLUMN(AntipTrkItsPidValid, antipTrkItsPidValid, int8_t); -DECLARE_SOA_COLUMN(AntipTrkTgl, antipTrkTgl, float); +DECLARE_SOA_COLUMN(AntipTrkItsNSigmaPr, antipTrkItsNSigmaPr, float); +DECLARE_SOA_COLUMN(AntipTrkItsPidValid, antipTrkItsPidValid, int8_t); +DECLARE_SOA_COLUMN(AntipTrkTgl, antipTrkTgl, float); } // namespace antin_cex // Table @@ -128,17 +128,15 @@ DECLARE_SOA_TABLE(AntinCexPairs, "AOD", "ANTINCEX", 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::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, antin_cex::PTrkItsNSigmaPr, antin_cex::PTrkItsPidValid, antin_cex::PTrkTgl, - antin_cex::AntipTrkItsNSigmaPr, antin_cex::AntipTrkItsPidValid, antin_cex::AntipTrkTgl - ); + antin_cex::AntipTrkItsNSigmaPr, antin_cex::AntipTrkItsPidValid, antin_cex::AntipTrkTgl); } // namespace o2::aod #endif // PWGLF_DATAMODEL_LFANTINCEXTABLES_H_ - diff --git a/PWGLF/TableProducer/Nuspex/nucleiAntineutronCex.cxx b/PWGLF/TableProducer/Nuspex/nucleiAntineutronCex.cxx index 5621b057b97..32a896db20c 100644 --- a/PWGLF/TableProducer/Nuspex/nucleiAntineutronCex.cxx +++ b/PWGLF/TableProducer/Nuspex/nucleiAntineutronCex.cxx @@ -15,24 +15,28 @@ /// \author Fabiola Lugo /// -#include +#include + +#include + +#include +#include +#include +#include #include -#include #include -#include -#include +#include #include -#include + #include -#include -#include #include -#include -#include #include -#include -#include -#include +#include + +#include +#include +#include +#include using namespace o2; using namespace o2::framework; @@ -41,90 +45,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.}}); @@ -135,7 +140,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.}}); @@ -146,57 +151,58 @@ 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.}}); + // ITS PID (protons / antiprotons, reconstructed tracks) - histos.add("pItsNsigmaPr", "ITS n#sigma (p hyp., proton);n#sigma_{ITS}(p);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("pItsNsigmaPr", "ITS n#sigma (p hyp., proton);n#sigma_{ITS}(p);Entries", kTH1F, {{100, -10., 10.}}); histos.add("apItsNsigmaPr", "ITS n#sigma (p hyp., antiproton);n#sigma_{ITS}(p);Entries", kTH1F, {{100, -10., 10.}}); - histos.add("pItsPidValid", "ITS PID valid flag (proton);PidValid;Entries", kTH1F, {{2, 0., 2.}}); + histos.add("pItsPidValid", "ITS PID valid flag (proton);PidValid;Entries", kTH1F, {{2, 0., 2.}}); histos.add("apItsPidValid", "ITS PID valid flag (antiproton);PidValid;Entries", kTH1F, {{2, 0., 2.}}); - histos.add("pTgl", "tgl (proton track);tgl;Entries", kTH1F, {{100, -2., 2.}}); + histos.add("pTgl", "tgl (proton track);tgl;Entries", kTH1F, {{100, -2., 2.}}); histos.add("apTgl", "tgl (antiproton track);tgl;Entries", kTH1F, {{100, -2., 2.}}); - histos.add("pItsNsigmaPr_bg", "ITS n#sigma (p hyp., proton);n#sigma_{ITS}(p);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("pItsNsigmaPr_bg", "ITS n#sigma (p hyp., proton);n#sigma_{ITS}(p);Entries", kTH1F, {{100, -10., 10.}}); histos.add("apItsNsigmaPr_bg", "ITS n#sigma (p hyp., antiproton);n#sigma_{ITS}(p);Entries", kTH1F, {{100, -10., 10.}}); - histos.add("pItsPidValid_bg", "ITS PID valid flag (proton);PidValid;Entries", kTH1F, {{2, 0., 2.}}); + histos.add("pItsPidValid_bg", "ITS PID valid flag (proton);PidValid;Entries", kTH1F, {{2, 0., 2.}}); histos.add("apItsPidValid_bg", "ITS PID valid flag (antiproton);PidValid;Entries", kTH1F, {{2, 0., 2.}}); - histos.add("pTgl_bg", "tgl (proton track);tgl;Entries", kTH1F, {{100, -2., 2.}}); + histos.add("pTgl_bg", "tgl (proton track);tgl;Entries", kTH1F, {{100, -2., 2.}}); histos.add("apTgl_bg", "tgl (antiproton track);tgl;Entries", kTH1F, {{100, -2., 2.}}); } - - 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(); @@ -205,9 +211,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()); @@ -215,7 +221,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()) { @@ -224,7 +231,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()) { @@ -233,7 +241,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()) { @@ -242,38 +251,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(); @@ -282,22 +294,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()) { @@ -305,23 +321,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()) { @@ -329,18 +348,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; @@ -354,21 +375,24 @@ 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; + if (!hasPrimaryMotherP) + continue; double protonVx = particle3.vx(); double protonVy = particle3.vy(); double protonVz = particle3.vz(); @@ -377,10 +401,10 @@ struct NucleiAntineutronCex { double pPz = particle3.pz(); double pE = particle3.e(); if (std::abs(protonVx - antipVx) < kVtxTol && std::abs(protonVy - antipVy) < kVtxTol && std::abs(protonVz - antipVz) < kVtxTol) { - //Same mother + // 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) { @@ -388,66 +412,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); @@ -455,15 +475,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); @@ -471,82 +495,86 @@ 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; - float pTrkItsNSigmaPr = -999.f; - int8_t pTrkItsPidValid = 0; - float pTrkTgl = 0.f; - + uint16_t apItsMap = 0; + float pTrkItsNSigmaPr = -999.f; + int8_t pTrkItsPidValid = 0; + float pTrkTgl = 0.f; + 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; - float antipTrkItsNSigmaPr = -999.f; + uint16_t pItsMap = 0; + float antipTrkItsNSigmaPr = -999.f; int8_t antipTrkItsPidValid = 0; - float antipTrkTgl = 0.f; - + float antipTrkTgl = 0.f; + o2::aod::ITSResponse itsResponse; - + 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(); antipTrkTgl = track.tgl(); const auto nsigmaITSantip = itsResponse.nSigmaITS(track); antipTrkItsNSigmaPr = static_cast(nsigmaITSantip); antipTrkItsPidValid = std::isfinite(nsigmaITSantip) ? 1 : 0; antipHasTrack = true; - apItsMap = static_cast(track.itsClusterMap()); + apItsMap = static_cast(track.itsClusterMap()); antipLayers = (apItsMap != 0); - if (layerCondition) antipLayers = true; + if (layerCondition) + antipLayers = true; if (motherPdg == -kNeutron) { histos.fill(HIST("apItsNsigmaPr"), antipTrkItsNSigmaPr); histos.fill(HIST("apItsPidValid"), antipTrkItsPidValid); @@ -557,25 +585,25 @@ struct NucleiAntineutronCex { histos.fill(HIST("apItsPidValid_bg"), antipTrkItsPidValid); histos.fill(HIST("apTgl_bg"), antipTrkTgl); } - } - else if (mc.globalIndex() == pId) { - pTrkP= track.p(); + } 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(); pTrkTgl = track.tgl(); const auto nsigmaITSp = - itsResponse.nSigmaITS(track); + itsResponse.nSigmaITS(track); pTrkItsNSigmaPr = static_cast(nsigmaITSp); pTrkItsPidValid = std::isfinite(nsigmaITSp) ? 1 : 0; 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 (motherPdg == -kNeutron) { histos.fill(HIST("pItsNsigmaPr"), pTrkItsNSigmaPr); histos.fill(HIST("pItsPidValid"), pTrkItsPidValid); @@ -588,50 +616,57 @@ struct NucleiAntineutronCex { } } } - 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) { - //TrackParCov - auto trP = makeTPCovFromAOD(*pRow); + // 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); @@ -639,12 +674,15 @@ 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; @@ -652,32 +690,32 @@ 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; for (int i = 0; i < static_cast(rLayers.size()); ++i) { @@ -687,8 +725,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); @@ -696,8 +734,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); @@ -706,108 +743,107 @@ 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, - - pTrkItsNSigmaPr, - pTrkItsPidValid, - pTrkTgl, - - antipTrkItsNSigmaPr, - antipTrkItsPidValid, - antipTrkTgl - ); + 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, + + pTrkItsNSigmaPr, + pTrkItsPidValid, + pTrkTgl, + + antipTrkItsNSigmaPr, + antipTrkItsPidValid, + antipTrkTgl); } } // ==== end DCAFitter2 ==== @@ -821,4 +857,3 @@ WorkflowSpec defineDataProcessing(ConfigContext const& ctx) { return WorkflowSpec{adaptAnalysisTask(ctx)}; } -