From e5ce99eb1b5b907155392ade24eec73c8b556b31 Mon Sep 17 00:00:00 2001 From: FabiolaLP Date: Wed, 3 Sep 2025 16:38:50 -0600 Subject: [PATCH 1/2] PWGLF/Nuspex: add nuclei-antineutroncex workflow --- PWGLF/Tasks/Nuspex/AntineutronCEX.cxx | 773 ++++++++++++++++++++++++++ PWGLF/Tasks/Nuspex/CMakeLists.txt | 5 + 2 files changed, 778 insertions(+) create mode 100644 PWGLF/Tasks/Nuspex/AntineutronCEX.cxx diff --git a/PWGLF/Tasks/Nuspex/AntineutronCEX.cxx b/PWGLF/Tasks/Nuspex/AntineutronCEX.cxx new file mode 100644 index 00000000000..79da961662b --- /dev/null +++ b/PWGLF/Tasks/Nuspex/AntineutronCEX.cxx @@ -0,0 +1,773 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* +Study for the antineutron CEX interaction background +author: Fabiola Lugo +*/ + +#include "Framework/AnalysisTask.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 +//ROOT +#include "TTree.h" +#include "TVector3.h" +#include "TMath.h" +#include + +using namespace o2; +using namespace o2::framework; + +struct AntineutronTask { + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + OutputObj cexTree{"cex_tree"}; // CEX + OutputObj bgTree{"bg_tree"}; // background + // === Branch buffers === + struct PairBuffers { + // MC (pair) + double mc_pair_p, mc_pair_pt, mc_pair_pz; + double mc_dplane; + double mc_angle_deg; + double mc_vtx_x, mc_vtx_y, mc_vtx_z; + // TRK (pair, fitter) + double trk_pair_p, trk_pair_pt, trk_pair_pz; + double trk_angle_deg; + double trk_vtxfit_dcaPair; + double trk_vtxfit_R; + double trk_vtxfit_distToPV; + double trk_vtxfit_secVtx_x, trk_vtxfit_secVtx_y, trk_vtxfit_secVtx_z; + // fit quality + double vtxfit_chi2; + int vtxfit_status; + int nCand; + double vtxfit_dX, vtxfit_dY, vtxfit_dZ; // (fit − MC) + double vtxfit_d3D; + // proton trk + double p_trk_p, p_trk_px, p_trk_py, p_trk_pz, p_trk_eta, p_trk_tpcSignal; + int p_trk_nClsITS; + // antiproton trk + double antip_trk_p, antip_trk_px, antip_trk_py, antip_trk_pz, antip_trk_eta, antip_trk_tpcSignal; + int antip_trk_nClsITS; + // Meta + int mother_pdg, colId, p_id, antip_id; + } buf; + + 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("antip_p", "Total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("antip_px", "p_{x};p_{x} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("antip_py", "p_{y};p_{y} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("antip_pz", "p_{z};p_{z} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("antip_eta","Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); + histos.add("antip_p_ITScuts", "Momentum with ITS cuts;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + + // Primary protons + histos.add("p_p", "Total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("p_px", "p_{x};p_{x} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("p_py", "p_{y};p_{y} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("p_pz", "p_{z};p_{z} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("p_eta","Pseudorapidity;#eta;Entries", kTH1F, {{100, -10., 10.}}); + histos.add("p_p_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.}}); + histos.add("primmom_test", "Secondary particles;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + + // CEX pair from antineutron (MC) + histos.add("cex_pairmc_p", "CEX pair total momentum;|p| (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("cex_pairmc_pt", "CEX pair p_{T};p_{T} (GeV/c);Entries", kTH1F, {{100, 0., 10.}}); + histos.add("cex_pairmc_pz", "CEX pair p_{z};p_{z} (GeV/c);Entries", kTH1F, {{100, -10., 10.}}); + histos.add("cex_pairmc_dplane","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("cex_pairmc_pITScuts","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_pairmc_dplane","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("cex_pairtrk_p", "Pair momentum (tracks);|p| (GeV/c);Entries", kTH1F, {{120, 0., 12.}}); + histos.add("cex_pairtrk_pt", "Pair p_{T} (tracks);p_{T} (GeV/c);Entries", kTH1F, {{120, 0., 12.}}); + histos.add("cex_pairtrk_pz", "Pair p_{z} (tracks);p_{z} (GeV/c);Entries", kTH1F, {{120, -12., 12.}}); + histos.add("cex_pairtrk_vtxfit_dcaPair", "DCA between tracks at PCA;DCA (cm);Entries", kTH1F, {{200, 0., 10.}}); + histos.add("cex_pairtrk_vtxfit_R", "Secondary-vertex radius (PCA);R (cm);Entries", kTH1F, {{200, 0., 60.}}); + histos.add("cex_pairtrk_vtxfit_distToPV", "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_pairtrk_vtxfit_dcaPair", "DCA between tracks at PCA;DCA (cm);Entries", kTH1F, {{200, 0., 10.}}); + histos.add("cexbg_pairtrk_vtxfit_R", "Secondary-vertex radius (PCA);R (cm);Entries", kTH1F, {{200, 0., 60.}}); + histos.add("cexbg_pairtrk_vtxfit_distToPV", "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("vtxfit_chi2", "DCAFitter2 #chi^{2};#chi^{2};Entries", kTH1F, {{200, 0., 100.}}); + histos.add("vtxfit_status", "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.}}); + + //Trees + cexTree.setObject(new TTree("cex_pairs", "CEX (antineutron daughters)")); + bgTree.setObject(new TTree("bg_pairs", "Background pairs")); + // Helper + auto defineBranches = [&](o2::framework::OutputObj& t) { + // MC + t->Branch("mc_pair_p", &buf.mc_pair_p); + t->Branch("mc_pair_pt", &buf.mc_pair_pt); + t->Branch("mc_pair_pz", &buf.mc_pair_pz); + t->Branch("mc_dplane", &buf.mc_dplane); + t->Branch("mc_angle_deg", &buf.mc_angle_deg); + t->Branch("mc_vtx_x", &buf.mc_vtx_x); + t->Branch("mc_vtx_y", &buf.mc_vtx_y); + t->Branch("mc_vtx_z", &buf.mc_vtx_z); + // TRK (par) + t->Branch("trk_pair_p", &buf.trk_pair_p); + t->Branch("trk_pair_pt", &buf.trk_pair_pt); + t->Branch("trk_pair_pz", &buf.trk_pair_pz); + t->Branch("trk_angle_deg",&buf.trk_angle_deg); + t->Branch("trk_vtxfit_dcaPair", &buf.trk_vtxfit_dcaPair); + t->Branch("trk_vtxfit_R", &buf.trk_vtxfit_R); + t->Branch("trk_vtxfit_distToPV", &buf.trk_vtxfit_distToPV); + t->Branch("trk_vtxfit_secVtx_x", &buf.trk_vtxfit_secVtx_x); + t->Branch("trk_vtxfit_secVtx_y", &buf.trk_vtxfit_secVtx_y); + t->Branch("trk_vtxfit_secVtx_z", &buf.trk_vtxfit_secVtx_z); + // Fit/residuales + t->Branch("vtxfit_chi2", &buf.vtxfit_chi2); + t->Branch("vtxfit_status", &buf.vtxfit_status); + t->Branch("nCand", &buf.nCand); + t->Branch("vtxfit_dX", &buf.vtxfit_dX); + t->Branch("vtxfit_dY", &buf.vtxfit_dY); + t->Branch("vtxfit_dZ", &buf.vtxfit_dZ); + t->Branch("vtxfit_d3D", &buf.vtxfit_d3D); + // Tracks individuales + t->Branch("p_trk_p", &buf.p_trk_p); + t->Branch("p_trk_px", &buf.p_trk_px); + t->Branch("p_trk_py", &buf.p_trk_py); + t->Branch("p_trk_pz", &buf.p_trk_pz); + t->Branch("p_trk_eta", &buf.p_trk_eta); + t->Branch("p_trk_tpcSignal", &buf.p_trk_tpcSignal); + t->Branch("p_trk_nClsITS", &buf.p_trk_nClsITS); + t->Branch("antip_trk_p", &buf.antip_trk_p); + t->Branch("antip_trk_px", &buf.antip_trk_px); + t->Branch("antip_trk_py", &buf.antip_trk_py); + t->Branch("antip_trk_pz", &buf.antip_trk_pz); + t->Branch("antip_trk_eta", &buf.antip_trk_eta); + t->Branch("antip_trk_tpcSignal", &buf.antip_trk_tpcSignal); + t->Branch("antip_trk_nClsITS", &buf.antip_trk_nClsITS); + // Meta + t->Branch("mother_pdg", &buf.mother_pdg); + t->Branch("colId", &buf.colId); + t->Branch("p_id", &buf.p_id); + t->Branch("antip_id", &buf.antip_id); + }; + + defineBranches(cexTree); + defineBranches(bgTree); + } + //Check available tables in the AOD, specifically TracksIU, TracksCovIU + using TracksWCovMc = soa::Join; + + 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 pvtx_x = 0; + double pvtx_y = 0; + double pvtx_z = 0; + for (auto const& col : cols) { + const auto colId= col.globalIndex(); + if (std::isfinite(col.posX()) && std::isfinite(col.posY()) && std::isfinite(col.posZ())) { + pvtx_x = col.posX(); + pvtx_y = col.posY(); + pvtx_z = col.posZ(); + histos.fill(HIST("hVx"), pvtx_x); + histos.fill(HIST("hVy"), pvtx_y); + histos.fill(HIST("hVz"), pvtx_z); + } + + for (auto& particle : particles) { + const auto procEnum = particle.getProcess(); + const bool isSecondaryFromMaterial = (!particle.producedByGenerator()) && (procEnum == kPHadronic || procEnum == kPHInhelastic); + if (particle.mcCollisionId() != colId) continue; + + // Primary antineutrons + if (particle.pdgCode() == -2112 && 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 (TMath::Abs(particle.eta()) < 1.5 && TMath::Abs(particle.vz())< 5.3) histos.fill(HIST("antin_p_ITScuts"),particle.p()); + } + // Primary neutrons + if (particle.pdgCode() == 2112 && 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 (TMath::Abs(particle.eta()) < 1.5 && TMath::Abs(particle.vz())< 5.3) histos.fill(HIST("n_p_ITScuts"),particle.p()); + } + // Primary antiprotons + if (particle.pdgCode() == -2212 && particle.isPhysicalPrimary()) { + histos.fill(HIST("antip_p"), particle.p()); + histos.fill(HIST("antip_px"), particle.px()); + histos.fill(HIST("antip_py"), particle.py()); + histos.fill(HIST("antip_pz"), particle.pz()); + histos.fill(HIST("antip_eta"), particle.eta()); + if (TMath::Abs(particle.eta()) < 1.5 && TMath::Abs(particle.vz())< 5.3) histos.fill(HIST("antip_p_ITScuts"),particle.p()); + } + // Primary protons + if (particle.pdgCode() == 2212 && particle.isPhysicalPrimary()) { + histos.fill(HIST("p_p"), particle.p()); + histos.fill(HIST("p_px"), particle.px()); + histos.fill(HIST("p_py"), particle.py()); + histos.fill(HIST("p_pz"), particle.pz()); + histos.fill(HIST("p_eta"), particle.eta()); + if (TMath::Abs(particle.eta()) < 1.5 && TMath::Abs(particle.vz())< 5.3) histos.fill(HIST("p_p_ITScuts"),particle.p()); + } + + if (particle.pdgCode() == -2212 && isSecondaryFromMaterial && !particle.mothersIds().empty()) { + for (auto mother : particle.mothers_as()) { + if ((mother.isPhysicalPrimary())) histos.fill(HIST("primmom_test"), particle.p()); + } + } + + //Seconday antiprotons from material + if (particle.pdgCode() == -2212 && isSecondaryFromMaterial && !particle.mothersIds().empty()) { + histos.fill(HIST("antip_test"), particle.p()); + //Primary mother + for (auto mother : particle.mothers_as()) { + if (!(mother.isPhysicalPrimary())) continue; + double mother_pt = mother.pt(); + double mother_pz = mother.pz(); + double mother_vz = mother.vz(); + double mother_p = mother.p(); + double mother_eta = mother.eta(); + int mother_pdg = mother.pdgCode(); + + double antip_vx = particle.vx(); + double antip_vy = particle.vy(); + double antip_vz = particle.vz(); + double antip_p = particle.p(); + double antip_px = particle.px(); + double antip_py = particle.py(); + double antip_pz = particle.pz(); + double antip_e = particle.e(); + double antip_eta = particle.eta(); + int antip_id = particle.globalIndex(); + + //Selection conditions + const double R = TMath::Sqrt(antip_vx*antip_vx + antip_vy*antip_vy); + //Config for ITS + //if(3.9<=R && R<=43.0 && TMath::Abs(antip_vz)<=48.9){ + //Config for ITS2 + if(2.2<=R && R<=39.0 && TMath::Abs(antip_vz)<=39.0){ + if (TMath::Abs(mother_eta) < 1.5 && TMath::Abs(mother_vz) < 5.3) + { + //Pion minus veto + bool pion_m = false; + for (auto& particle1 : particles) { + const auto proc1Enum = particle1.getProcess(); + const bool isSecondaryFromMaterial1 = (!particle1.producedByGenerator()) && (proc1Enum == kPHadronic || proc1Enum == kPHInhelastic); + if (particle1.mcCollisionId() != colId) continue; + if (particle1.pdgCode() != -211 || !isSecondaryFromMaterial1 || particle1.mothersIds().empty()) continue; + bool hasPrimaryMother_pim = false; + for (auto mother : particle1.mothers_as()) { + if (mother.isPhysicalPrimary()) { + hasPrimaryMother_pim = true; + break; + } + } + if (!hasPrimaryMother_pim) continue; + double pim_vx = particle1.vx(); + double pim_vy = particle1.vy(); + double pim_vz = particle1.vz(); + if(pim_vx == antip_vx && pim_vy == antip_vy && pim_vz == antip_vz){ + pion_m = true; + break; + } + } + //Pion plus veto + bool pion_p = false; + for (auto& particle2 : particles) { + if (particle2.mcCollisionId() != colId) continue; + const auto proc2Enum = particle2.getProcess(); + const bool isSecondaryFromMaterial2 = (!particle2.producedByGenerator()) && (proc2Enum == kPHadronic || proc2Enum == kPHInhelastic); + if (particle2.pdgCode() != 211 || !isSecondaryFromMaterial2 || particle2.mothersIds().empty()) continue; + bool hasPrimaryMother_pip = false; + for (auto mother : particle2.mothers_as()) { + if (mother.isPhysicalPrimary()) { + hasPrimaryMother_pip = true; + break; + } + } + if (!hasPrimaryMother_pip) continue; + double pip_vx = particle2.vx(); + double pip_vy = particle2.vy(); + double pip_vz = particle2.vz(); + if(pip_vx == antip_vx && pip_vy == antip_vy && pip_vz == antip_vz){ + pion_p = true; + break; + } + } + if (pion_p == false && pion_m == false){ + // CEX selection + double dplane = 10; + double dplane_tmp = 0; + double p = 0; + double p_tmp = 0; + double pcex_px = 0; + double pcex_py = 0; + double pcex_pz = 0; + double pcex_p = 0; + double e = 0; + double e_tmp = 0; + double pcex_e = 0; + int k_plane = -1; + int k_e = -1; + int k_p = -1; + //Secondary proton from material + for (auto& particle3 : particles) { + if (particle3.mcCollisionId() != colId) continue; + const auto proc3Enum = particle3.getProcess(); + const bool isSecondaryFromMaterial3 = (!particle3.producedByGenerator()) && (proc3Enum == kPHadronic || proc3Enum == kPHInhelastic); + if (particle3.pdgCode() != 2212 || !isSecondaryFromMaterial3 || particle3.mothersIds().empty()) continue; + bool hasPrimaryMother_p = false; + for (auto mother : particle3.mothers_as()) { + if (mother.isPhysicalPrimary()) { + hasPrimaryMother_p = true; + break; + } + } + if (!hasPrimaryMother_p) continue; + double p_vx = particle3.vx(); + double p_vy = particle3.vy(); + double p_vz = particle3.vz(); + double p_p = particle3.p(); + double p_px = particle3.px(); + double p_py = particle3.py(); + double p_pz = particle3.pz(); + double p_e = particle3.e(); + double p_eta = particle3.eta(); + + if(p_vx == antip_vx && p_vy == antip_vy && p_vz == antip_vz){ + bool shareMother = false; + { + const auto& momsAp = particle.mothersIds(); //antiproton + const auto& momsP = particle3.mothersIds(); //proton + for (auto ida : momsAp) { + for (auto idp : momsP) { + if (ida == idp) { + shareMother = true; + break; + } + } + if (shareMother) break; + } + } + if (!shareMother) continue; + //dplane_tmp = (p_py*antip_pz - p_pz*antip_py)*(pvtx_x-antip_vx) + (p_pz*antip_px - p_px*antip_pz)*(pvtx_y-antip_vy) + (p_px*antip_py - p_py*antip_px)*(pvtx_z-antip_vz); + double nx = (p_py*antip_pz - p_pz*antip_py); + double ny = (p_pz*antip_px - p_px*antip_pz); + double nz = (p_px*antip_py - p_py*antip_px); + double rx = (pvtx_x - antip_vx); + double ry = (pvtx_y - antip_vy); + double rz = (pvtx_z - antip_vz); + double denom = nx*nx + ny*ny + nz*nz; + if (denom > 0.) { + dplane_tmp = TMath::Abs(nx*rx + ny*ry + nz*rz) / TMath::Sqrt(denom); + } else { + dplane_tmp = 1e9; + } + if(TMath::Abs(dplane_tmp) < TMath::Abs(dplane)) + { + k_plane = particle3.globalIndex(); + dplane = dplane_tmp; + } + + e_tmp = antip_e + p_e; + if(TMath::Abs(e_tmp) > TMath::Abs(e)) + { + k_e = particle3.globalIndex(); + e = e_tmp; + pcex_e = p_e; + } + + p_tmp = TMath::Sqrt(TMath::Power((p_px+antip_px),2)+TMath::Power((p_py+antip_py),2)+TMath::Power((p_pz+antip_pz),2)); + if(TMath::Abs(p_tmp) > TMath::Abs(p)) + { + k_p = particle3.globalIndex(); + p = p_tmp; + pcex_p = p_p; + pcex_px = p_px; + pcex_py = p_py; + pcex_pz = p_pz; + } + if(k_plane == k_e && k_plane == k_p && k_plane >= 0 ){ + int p_id = k_plane; + TVector3 pVecProton = TVector3(pcex_px, pcex_py, pcex_pz); + TVector3 pVecAntiproton = TVector3(antip_px, antip_py, antip_pz); + TVector3 total_mc_pVec = pVecProton + pVecAntiproton; + double cex_pairmc_p = total_mc_pVec.Mag(); + double cex_pairmc_pt = total_mc_pVec.Pt(); + double cex_pairmc_pz = pcex_pz+antip_pz; + double mcangleRad = pVecProton.Angle(pVecAntiproton); + double mcangleDeg = mcangleRad * TMath::RadToDeg(); + + //Antineutron mother + if (mother_pdg == -2112) { + //CEX pair + histos.fill(HIST("cex_pairmc_p"), cex_pairmc_p); + histos.fill(HIST("cex_pairmc_pt"), cex_pairmc_pt); + histos.fill(HIST("cex_pairmc_pz"), cex_pairmc_pz); + histos.fill(HIST("cex_pairmc_dplane"), dplane); + histos.fill(HIST("cex_pairmc_angle"), mcangleDeg); + histos.fill(HIST("cex_pairmc_vtx"), antip_vx, antip_vy); + histos.fill(HIST("cex_pairmc_vtxz"), antip_vz); + if (TMath::Abs(mother_eta) < 0.9 && TMath::Abs(mother_vz)< 5.3) histos.fill(HIST("cex_pairmc_pITScuts"),cex_pairmc_p); + //CEX pair normalized + if (mother_p != 0) histos.fill(HIST("cexn_pairmc_p"), cex_pairmc_p /mother_p); + if (mother_pt != 0) histos.fill(HIST("cexn_pairmc_pt"), cex_pairmc_pt/mother_pt); + if (mother_pz != 0) histos.fill(HIST("cexn_pairmc_pz"), cex_pairmc_pz/mother_pz); + } + //BG mother + if (mother_pdg != -2112) { + //CEX pair + histos.fill(HIST("cexbg_pairmc_p"), cex_pairmc_p); + histos.fill(HIST("cexbg_pairmc_pt"), cex_pairmc_pt); + histos.fill(HIST("cexbg_pairmc_pz"), cex_pairmc_pz); + histos.fill(HIST("cexbg_pairmc_dplane"), dplane); + histos.fill(HIST("cexbg_pairmc_angle"), mcangleDeg); + histos.fill(HIST("cexbg_pairmc_vtx"), antip_vx, antip_vy); + histos.fill(HIST("cexbg_pairmc_vtxz"), antip_vz); + if (TMath::Abs(mother_eta) < 0.9 && TMath::Abs(mother_vz)< 5.3) histos.fill(HIST("cexbg_pairmc_pITScuts"),cex_pairmc_p); + } + + //Detector signal + bool antip_layers = false; + bool antip_hasTrack = false; + double antip_trk_px = 0.; + double antip_trk_py = 0.; + double antip_trk_pz = 0.; + double antip_trk_p = 0.; + double antip_trk_eta = 0.; + double antip_trk_tpcSignal = 0; + //int antip_trk_nClsTPC = 0; + int antip_trk_nClsITS = 0; + + bool p_layers = false; + bool p_hasTrack = false; + double p_trk_px = 0.; + double p_trk_py = 0.; + double p_trk_pz = 0.; + double p_trk_p = 0.; + double p_trk_eta = 0.; + double p_trk_tpcSignal = 0; + //int p_trk_nClsTPC = 0; + int p_trk_nClsITS = 0; + + for (auto& track : tracks) { + if (!track.has_mcParticle()) continue; + const auto& mc = track.mcParticle(); + 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 layer_condition = (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 layer_condition = (!hitIB) && hitOuter && (nITS >= 2); + + if (mc.globalIndex() == antip_id) { + antip_trk_p = track.p(); + antip_trk_px = track.px(); + antip_trk_py = track.py(); + antip_trk_pz = track.pz(); + antip_trk_eta = track.eta(); + antip_trk_tpcSignal = track.tpcSignal(); + //antip_trk_nClsTPC = track.tpcNCls(); + antip_trk_nClsITS = track.itsNCls(); + antip_hasTrack = true; + if (layer_condition) antip_layers = true; + } + else if (mc.globalIndex() == p_id) { + p_trk_p = track.p(); + p_trk_px = track.px(); + p_trk_py = track.py(); + p_trk_pz = track.pz(); + p_trk_eta = track.eta(); + p_trk_tpcSignal = track.tpcSignal(); + //p_trk_nClsTPC = track.tpcNCls(); + p_trk_nClsITS = track.itsNCls(); + p_hasTrack = true; + if (layer_condition) p_layers = true; + } + } + + if (p_hasTrack == true && antip_hasTrack == true){ + TVector3 pVecProton_trk(p_trk_px, p_trk_py, p_trk_pz); + TVector3 AntipVecProton_trk(antip_trk_px, antip_trk_py, antip_trk_pz); + TVector3 total_trk_pVec = pVecProton_trk + AntipVecProton_trk; + double trkangleRad = AntipVecProton_trk.Angle(pVecProton_trk); + double trkangleDeg = trkangleRad * TMath::RadToDeg(); + if(mother_pdg == -2112) histos.fill(HIST("cex_pairtrk_angle"), trkangleDeg); + if(mother_pdg != -2112) histos.fill(HIST("cexbg_pairtrk_angle"), trkangleDeg); + + // ==== Secondary vertex via central DCA vertexer (DCAFitter2) ==== + using o2::vertexing::DCAFitter2; + constexpr float BZ_TESLA = 0.5f; + + DCAFitter2 fitter(/*bz=*/BZ_TESLA, /*useAbsDCA=*/true, /*propagateToPCA=*/true); + fitter.setBz(BZ_TESLA); + + //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 (auto& tr : tracks) { + if (!tr.has_mcParticle()) continue; + const auto& mc = tr.mcParticle(); + if (mc.globalIndex() == antip_id) apRow = tr; + if (mc.globalIndex() == p_id) 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("vtxfit_status"), (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 (mother_pdg == -2112) histos.fill(HIST("cex_pairtrk_vtxfit_dcaPair"), dcaPair); + if (mother_pdg != -2112) histos.fill(HIST("cexbg_pairtrk_vtxfit_dcaPair"), dcaPair); + + //Cuts + //if (dcaPair > 10.0) continue; + //if (trkangleDeg > 60) continue; + //if (TMath::Abs(dplane) > 2) continue; + if (antip_layers == true && p_layers == true){ + double cex_pairtrk_p = total_trk_pVec.Mag(); + double cex_pairtrk_pt = total_trk_pVec.Pt(); + double cex_pairtrk_pz = p_trk_pz + antip_trk_pz; + const double radius = std::hypot(secX, secY); + const double dxPV = secX - pvtx_x; + const double dyPV = secY - pvtx_y; + const double dzPV = secZ - pvtx_z; + const double distToPrimary = std::sqrt(dxPV*dxPV + dyPV*dyPV + dzPV*dzPV); + + if(mother_pdg == -2112){ + histos.fill(HIST("cex_pairtrk_p"), cex_pairtrk_p); + histos.fill(HIST("cex_pairtrk_pt"), cex_pairtrk_pt); + histos.fill(HIST("cex_pairtrk_pz"), cex_pairtrk_pz); + histos.fill(HIST("cex_pairtrk_vtxfit_R"), radius); + histos.fill(HIST("cex_pairtrk_vtxfit_distToPV"), 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"), cex_pairtrk_p); + histos.fill(HIST("cexbg_pairtrk_pt"), cex_pairtrk_pt); + histos.fill(HIST("cexbg_pairtrk_pz"), cex_pairtrk_pz); + histos.fill(HIST("cexbg_pairtrk_vtxfit_R"), radius); + histos.fill(HIST("cexbg_pairtrk_vtxfit_distToPV"), 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("vtxfit_chi2"), chi2); + const double dX = secX - antip_vx; + const double dY = secY - antip_vy; + const double dZ = secZ - antip_vz; + 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); + + // Fill Trees + buf.mc_pair_p = cex_pairmc_p; + buf.mc_pair_pt = cex_pairmc_pt; + buf.mc_pair_pz = cex_pairmc_pz; + buf.mc_dplane = dplane; + buf.mc_angle_deg = mcangleDeg; + buf.mc_vtx_x = antip_vx; + buf.mc_vtx_y = antip_vy; + buf.mc_vtx_z = antip_vz; + buf.trk_pair_p = cex_pairtrk_p; + buf.trk_pair_pt = cex_pairtrk_pt; + buf.trk_pair_pz = cex_pairtrk_pz; + buf.trk_angle_deg = trkangleDeg; + buf.trk_vtxfit_dcaPair = dcaPair; + buf.trk_vtxfit_R = radius; + buf.trk_vtxfit_distToPV = distToPrimary; + buf.trk_vtxfit_secVtx_x = secX; + buf.trk_vtxfit_secVtx_y = secY; + buf.trk_vtxfit_secVtx_z = secZ; + buf.vtxfit_chi2 = chi2; + buf.vtxfit_status = (int)status; + buf.nCand = nCand; + buf.vtxfit_dX = secX - antip_vx; + buf.vtxfit_dY = secY - antip_vy; + buf.vtxfit_dZ = secZ - antip_vz; + buf.vtxfit_d3D = std::sqrt(buf.vtxfit_dX*buf.vtxfit_dX + buf.vtxfit_dY*buf.vtxfit_dY + buf.vtxfit_dZ*buf.vtxfit_dZ); + buf.p_trk_p = p_trk_p; + buf.p_trk_px = p_trk_px; + buf.p_trk_py = p_trk_py; + buf.p_trk_pz = p_trk_pz; + buf.p_trk_eta = p_trk_eta; + buf.p_trk_tpcSignal = p_trk_tpcSignal; + buf.p_trk_nClsITS = p_trk_nClsITS; + buf.antip_trk_p = antip_trk_p; + buf.antip_trk_px = antip_trk_px; + buf.antip_trk_py = antip_trk_py; + buf.antip_trk_pz = antip_trk_pz; + buf.antip_trk_eta = antip_trk_eta; + buf.antip_trk_tpcSignal = antip_trk_tpcSignal; + buf.antip_trk_nClsITS = antip_trk_nClsITS; + buf.mother_pdg = mother_pdg; + buf.colId = colId; + buf.p_id = p_id; + buf.antip_id = antip_id; + + if (mother_pdg == -2112) cexTree->Fill(); + else bgTree->Fill(); + } + } + } + // ==== end DCAFitter2 ==== + } + } + } + } + } + } + } + break; + } + } + } + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& ctx) +{ + return WorkflowSpec{adaptAnalysisTask(ctx)}; +} + diff --git a/PWGLF/Tasks/Nuspex/CMakeLists.txt b/PWGLF/Tasks/Nuspex/CMakeLists.txt index 38d5704f683..9e230d5c66b 100644 --- a/PWGLF/Tasks/Nuspex/CMakeLists.txt +++ b/PWGLF/Tasks/Nuspex/CMakeLists.txt @@ -134,6 +134,11 @@ o2physics_add_dpl_workflow(nuclei-from-hypertriton-map PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(nuclei-antineutroncex + SOURCES AntineutronCEX.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + if(FastJet_FOUND) o2physics_add_dpl_workflow(angular-correlations-in-jets SOURCES angularCorrelationsInJets.cxx From 6d9fbaa446b98f56653953031c2b692702372716 Mon Sep 17 00:00:00 2001 From: FabiolaLP Date: Sun, 21 Sep 2025 22:44:53 -0600 Subject: [PATCH 2/2] PWGLF: add nuclei-antineutron-cex workflow and ANTINCEX AOD table --- PWGLF/DataModel/LFAntinCexTables.h | 90 +++ PWGLF/TableProducer/Nuspex/CMakeLists.txt | 5 + .../Nuspex/nucleiAntineutronCex.cxx | 708 ++++++++++++++++++ 3 files changed, 803 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..80822e3f155 --- /dev/null +++ b/PWGLF/DataModel/LFAntinCexTables.h @@ -0,0 +1,90 @@ +#ifndef PWGLF_DATAMODEL_LFANTINCEXTABLES_H_ +#define PWGLF_DATAMODEL_LFANTINCEXTABLES_H_ + +#include "Framework/AnalysisDataModel.h" +#include "Framework/ASoAHelpers.h" + +namespace o2::aod +{ +namespace antinCexNS +{ +// Etiqueta/metadata básica +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 (par) +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 (par, del 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); + +// Calidad del fit y residuales (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); + +// Track del protón +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); + +// Track del antiprotón +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); + +// (Opcional pero útil) Máscara de selección por cortes +DECLARE_SOA_COLUMN(SelMask, selMask, uint32_t); // bits de cortes superados +} // namespace antinCexNS + +// Tabla única para CEX y BG (flag isCex) +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 +) + +} // 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..a282078865d --- /dev/null +++ b/PWGLF/TableProducer/Nuspex/nucleiAntineutronCex.cxx @@ -0,0 +1,708 @@ +/// \file nucleiAntineutronCex.cxx +/// \brief Analysis task for antineutron charge-exchange (n̄ + n → p + p̄) +/// \author Fabiola Lugo + +#include "Framework/AnalysisTask.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 +//ROOT +#include "TVector3.h" +#include "TMath.h" +#include +#include +#include "CommonConstants/MathConstants.h" +#include "TPDGCode.h" +// +#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; + + 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; + + 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; + 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; + 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); + + 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 + ); + } + } + // ==== end DCAFitter2 ==== + } + } + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& ctx) +{ + return WorkflowSpec{adaptAnalysisTask(ctx)}; +} +