diff --git a/PWGHF/DataModel/CandidateReconstructionTables.h b/PWGHF/DataModel/CandidateReconstructionTables.h index 50dc192d66d..54f17c4c478 100644 --- a/PWGHF/DataModel/CandidateReconstructionTables.h +++ b/PWGHF/DataModel/CandidateReconstructionTables.h @@ -1394,6 +1394,9 @@ namespace hf_cand_toxipi { // Data processing results: DECLARE_SOA_INDEX_COLUMN(Collision, collision); +DECLARE_SOA_COLUMN(XPv, xPv, float); +DECLARE_SOA_COLUMN(YPv, yPv, float); +DECLARE_SOA_COLUMN(ZPv, zPv, float); DECLARE_SOA_COLUMN(XDecayVtxOmegac, xDecayVtxOmegac, float); DECLARE_SOA_COLUMN(YDecayVtxOmegac, yDecayVtxOmegac, float); DECLARE_SOA_COLUMN(ZDecayVtxOmegac, zDecayVtxOmegac, float); @@ -1459,7 +1462,6 @@ DECLARE_SOA_INDEX_COLUMN_FULL(PrimaryPi, primaryPi, int, Tracks, "_primarypi"); DECLARE_SOA_COLUMN(ImpactParOmegacXY, impactParOmegacXY, float); DECLARE_SOA_COLUMN(ImpactParOmegacZ, impactParOmegacZ, float); DECLARE_SOA_COLUMN(InvMassLambda, invMassLambda, double); -DECLARE_SOA_COLUMN(InvMassAntiLambda, invMassAntiLambda, double); DECLARE_SOA_COLUMN(InvMassCascade, invMassCascade, double); DECLARE_SOA_COLUMN(InvMassOmegac, invMassOmegac, double); DECLARE_SOA_COLUMN(CosPAV0, cosPAV0, double); @@ -1479,13 +1481,15 @@ DECLARE_SOA_COLUMN(EtaPiFromOme, etaPiFromOme, double); DECLARE_SOA_COLUMN(EtaOmegac, etaOmegac, double); DECLARE_SOA_COLUMN(EtaCascade, etaCascade, double); DECLARE_SOA_COLUMN(EtaV0, etaV0, double); -DECLARE_SOA_COLUMN(DcaXYToPVPrimaryPi, dcaXYToPVPrimaryPi, double); -DECLARE_SOA_COLUMN(DcaXYToPVV0Dau0, dcaXYToPVV0dau0, double); -DECLARE_SOA_COLUMN(DcaXYToPVV0Dau1, dcaXYToPVV0dau1, double); -DECLARE_SOA_COLUMN(DcaXYToPVCascDau, dcaXYToPVCascdau, double); -DECLARE_SOA_COLUMN(DcaCascDau, dcaCascDau, double); -DECLARE_SOA_COLUMN(DcaV0Dau, dcaV0Dau, double); -DECLARE_SOA_COLUMN(DcaOmegacDau, dcaOmegacDau, double); +DECLARE_SOA_COLUMN(DcaXYToPvV0Dau0, dcaXYToPvV0Dau0, float); +DECLARE_SOA_COLUMN(DcaXYToPvV0Dau1, dcaXYToPvV0Dau1, float); +DECLARE_SOA_COLUMN(DcaXYToPvCascDau, dcaXYToPvCascDau, float); +DECLARE_SOA_COLUMN(DcaZToPvV0Dau0, dcaZToPvV0Dau0, float); +DECLARE_SOA_COLUMN(DcaZToPvV0Dau1, dcaZToPvV0Dau1, float); +DECLARE_SOA_COLUMN(DcaZToPvCascDau, dcaZToPvCascDau, float); +DECLARE_SOA_COLUMN(DcaCascDau, dcaCascDau, float); +DECLARE_SOA_COLUMN(DcaV0Dau, dcaV0Dau, float); +DECLARE_SOA_COLUMN(DcaOmegacDau, dcaOmegacDau, float); // MC matching result: DECLARE_SOA_COLUMN(FlagMcMatchRec, flagMcMatchRec, int8_t); // reconstruction level @@ -1502,7 +1506,7 @@ enum DecayType { DecayToXiPi = 0, // declare dedicated Omegac and Xic to Xi Pi candidate table DECLARE_SOA_TABLE(HfCandToXiPi, "AOD", "HFCANDTOXIPI", o2::soa::Index<>, - hf_cand_toxipi::CollisionId, collision::PosX, collision::PosY, collision::PosZ, + hf_cand_toxipi::CollisionId, hf_cand_toxipi::XPv, hf_cand_toxipi::YPv, hf_cand_toxipi::ZPv, hf_cand_toxipi::XDecayVtxOmegac, hf_cand_toxipi::YDecayVtxOmegac, hf_cand_toxipi::ZDecayVtxOmegac, hf_cand_toxipi::XDecayVtxCascade, hf_cand_toxipi::YDecayVtxCascade, hf_cand_toxipi::ZDecayVtxCascade, hf_cand_toxipi::XDecayVtxV0, hf_cand_toxipi::YDecayVtxV0, hf_cand_toxipi::ZDecayVtxV0, @@ -1523,12 +1527,13 @@ DECLARE_SOA_TABLE(HfCandToXiPi, "AOD", "HFCANDTOXIPI", hf_cand_toxipi::ErrImpactParCascXY, hf_cand_toxipi::ErrImpactParPrimaryPiXY, hf_cand_toxipi::ErrImpactParV0XY, hf_cand_toxipi::V0Id, v0data::PosTrackId, v0data::NegTrackId, hf_cand_toxipi::CascadeId, hf_cand_toxipi::PrimaryPiId, cascdata::BachelorId, hf_cand_toxipi::ImpactParOmegacXY, hf_cand_toxipi::ImpactParOmegacZ, - hf_cand_toxipi::InvMassLambda, hf_cand_toxipi::InvMassAntiLambda, hf_cand_toxipi::InvMassCascade, hf_cand_toxipi::InvMassOmegac, + hf_cand_toxipi::InvMassLambda, hf_cand_toxipi::InvMassCascade, hf_cand_toxipi::InvMassOmegac, hf_cand_toxipi::CosPAV0, hf_cand_toxipi::CosPAOmegac, hf_cand_toxipi::CosPACasc, hf_cand_toxipi::CosPAXYV0, hf_cand_toxipi::CosPAXYOmegac, hf_cand_toxipi::CosPAXYCasc, hf_cand_toxipi::CTauOmegac, hf_cand_toxipi::CTauCascade, hf_cand_toxipi::CTauV0, hf_cand_toxipi::CTauXic, hf_cand_toxipi::EtaV0PosDau, hf_cand_toxipi::EtaV0NegDau, hf_cand_toxipi::EtaPiFromCasc, hf_cand_toxipi::EtaPiFromOme, hf_cand_toxipi::EtaOmegac, hf_cand_toxipi::EtaCascade, hf_cand_toxipi::EtaV0, - hf_cand_toxipi::DcaXYToPVPrimaryPi, hf_cand_toxipi::DcaXYToPVV0Dau0, hf_cand_toxipi::DcaXYToPVV0Dau1, hf_cand_toxipi::DcaXYToPVCascDau, + hf_cand_toxipi::DcaXYToPvV0Dau0, hf_cand_toxipi::DcaXYToPvV0Dau1, hf_cand_toxipi::DcaXYToPvCascDau, + hf_cand_toxipi::DcaZToPvV0Dau0, hf_cand_toxipi::DcaZToPvV0Dau1, hf_cand_toxipi::DcaZToPvCascDau, hf_cand_toxipi::DcaCascDau, hf_cand_toxipi::DcaV0Dau, hf_cand_toxipi::DcaOmegacDau, hf_track_index::HFflag); // table with results of reconstruction level MC matching diff --git a/PWGHF/TableProducer/candidateCreatorToXiPi.cxx b/PWGHF/TableProducer/candidateCreatorToXiPi.cxx index db749503cb7..8a4bc59c428 100644 --- a/PWGHF/TableProducer/candidateCreatorToXiPi.cxx +++ b/PWGHF/TableProducer/candidateCreatorToXiPi.cxx @@ -17,6 +17,7 @@ #include "Common/Core/trackUtilities.h" #include "Common/Core/RecoDecay.h" #include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/CollisionAssociation.h" #include "DataFormatsParameters/GRPMagField.h" #include "DataFormatsParameters/GRPObject.h" #include "DetectorsBase/Propagator.h" @@ -50,6 +51,8 @@ using namespace o2::aod::hf_cand_toxipi; struct HfCandidateCreatorToXiPi { Produces rowCandidate; + Configurable doPvRefit{"doPvRefit", false, "set to true if you do PV refit in trackIndexSkimCreator.cxx"}; + Configurable propagateToPCA{"propagateToPCA", false, "create tracks version propagated to PCA"}; Configurable useAbsDCA{"useAbsDCA", true, "Minimise abs. distance rather than chi2"}; Configurable useWeightedFinalPCA{"useWeightedFinalPCA", true, "Recalculate vertex position using track covariances, effective only if useAbsDCA is true"}; @@ -75,14 +78,19 @@ struct HfCandidateCreatorToXiPi { o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrLUT; int runNumber; - // filter to use only HF selected collisions - Filter filterSelectCollisions = (aod::hf_sel_collision::whyRejectColl == 0); - using SelectedCollisions = soa::Filtered>; - using MyTracks = soa::Join; - using MyCascTable = soa::Join; + using MyTracks = soa::Join; + using FilteredHfTrackAssocSel = soa::Filtered>; + using MyCascTable = soa::Join; using MyV0Table = soa::Join; + Filter filterSelectCollisions = (aod::hf_sel_collision::whyRejectColl == 0); // filter to use only HF selected collisions + Filter filterSelectTrackIds = (aod::hf_sel_track::isSelProng >= 4); + + Preslice tracksPerCollision = aod::track::collisionId; // needed for PV refit + Preslice trackIndicesPerCollision = aod::track_association::collisionId; // aod::hf_track_association::collisionId + Preslice cascadesPerCollision = aod::cascdata::collisionId; + OutputObj hInvMassOmegac{TH1F("hInvMassOmegac", "Omegac invariant mass;inv mass;entries", 500, 2.2, 3.1)}; void init(InitContext const&) @@ -94,251 +102,307 @@ struct HfCandidateCreatorToXiPi { runNumber = 0; } - void process(SelectedCollisions::iterator const& collision, + void process(SelectedCollisions const& collisions, aod::BCsWithTimestamps const& bcWithTimeStamps, - MyCascTable const& cascades, MyTracks const& tracks, + FilteredHfTrackAssocSel const& trackIndices, + MyCascTable const& cascades, MyV0Table const&, aod::V0sLinked const&) { - // set the magnetic field from CCDB - auto bc = collision.bc_as(); - initCCDB(bc, runNumber, ccdb, isRun2 ? ccdbPathGrp : ccdbPathGrpMag, lut, isRun2); - auto magneticField = o2::base::Propagator::Instance()->getNominalBz(); // z component - - // 2-prong vertex fitter to build the omegac vertex - o2::vertexing::DCAFitterN<2> df; - df.setBz(magneticField); - df.setPropagateToPCA(propagateToPCA); - df.setMaxR(maxR); - df.setMaxDZIni(maxDZIni); - df.setMaxDXYIni(maxDXYIni); - df.setMinParamChange(minParamChange); - df.setMinRelChi2Change(minRelChi2Change); - df.setMaxChi2(maxChi2); - df.setUseAbsDCA(useAbsDCA); - df.setWeightedFinalPCA(useWeightedFinalPCA); - df.setRefitWithMatCorr(refitWithMatCorr); - - double massPionFromPDG = RecoDecay::getMassPDG(kPiPlus); // pdg code 211 - double massLambdaFromPDG = RecoDecay::getMassPDG(kLambda0); // pdg code 3122 - double massXiFromPDG = RecoDecay::getMassPDG(kXiMinus); // pdg code 3312 - double massOmegacFromPDG = RecoDecay::getMassPDG(kOmegaC0); // pdg code 4332 - double massXicFromPDG = RecoDecay::getMassPDG(kXiCZero); // pdg code 4132 - - // loop over cascades reconstructed by cascadebuilder.cxx - for (auto const& casc : cascades) { - - if (collision.globalIndex() != casc.collisionId()) { // check to be further processed when the problem of ambiguous tracks will be solved - continue; - } - - //----------------accessing particles in the decay chain------------- - // cascade daughter - charged particle - // int indexTrackXiDauCharged = casc.bachelorId(); // pion <- xi index from cascade table (not used) - auto trackXiDauCharged = casc.bachelor_as(); // pion <- xi track from MyTracks table - // cascade daughter - V0 - if (!casc.v0_as().has_v0Data()) { // check that V0 data are stored - continue; - } - auto v0 = casc.v0_as(); - auto v0Element = v0.v0Data_as(); // V0 element from LF table containing V0 info - // V0 positive daughter - auto trackV0Dau0 = v0Element.posTrack_as(); // p <- V0 track (positive track) from MyTracks table - // V0 negative daughter - auto trackV0Dau1 = v0Element.negTrack_as(); // pion <- V0 track (negative track) from MyTracks table - - // check that particles come from the same collision - if (rejDiffCollTrack) { // check to be further processed when the problem of ambiguous tracks will be solved - if (trackV0Dau0.collisionId() != trackV0Dau1.collisionId()) { + for (const auto& collision : collisions) { + + // set the magnetic field from CCDB + auto bc = collision.bc_as(); + initCCDB(bc, runNumber, ccdb, isRun2 ? ccdbPathGrp : ccdbPathGrpMag, lut, isRun2); + auto magneticField = o2::base::Propagator::Instance()->getNominalBz(); // z component + + // 2-prong vertex fitter to build the omegac vertex + o2::vertexing::DCAFitterN<2> df; + df.setBz(magneticField); + df.setPropagateToPCA(propagateToPCA); + df.setMaxR(maxR); + df.setMaxDZIni(maxDZIni); + df.setMaxDXYIni(maxDXYIni); + df.setMinParamChange(minParamChange); + df.setMinRelChi2Change(minRelChi2Change); + df.setMaxChi2(maxChi2); + df.setUseAbsDCA(useAbsDCA); + df.setWeightedFinalPCA(useWeightedFinalPCA); + df.setRefitWithMatCorr(refitWithMatCorr); + + double massPionFromPDG = RecoDecay::getMassPDG(kPiPlus); // pdg code 211 + double massLambdaFromPDG = RecoDecay::getMassPDG(kLambda0); // pdg code 3122 + double massXiFromPDG = RecoDecay::getMassPDG(kXiMinus); // pdg code 3312 + double massOmegacFromPDG = RecoDecay::getMassPDG(kOmegaC0); // pdg code 4332 + double massXicFromPDG = RecoDecay::getMassPDG(kXiCZero); // pdg code 4132 + + // loop over cascades reconstructed by cascadebuilder.cxx + auto thisCollId = collision.globalIndex(); + auto groupedCascades = cascades.sliceBy(cascadesPerCollision, thisCollId); + + for (auto const& casc : groupedCascades) { + + //----------------accessing particles in the decay chain------------- + // cascade daughter - charged particle + // int indexTrackXiDauCharged = casc.bachelorId(); // pion <- xi index from cascade table (not used) + auto trackXiDauCharged = casc.bachelor_as(); // pion <- xi track from MyTracks table + // cascade daughter - V0 + if (!casc.v0_as().has_v0Data()) { // check that V0 data are stored continue; } - if (trackXiDauCharged.collisionId() != trackV0Dau0.collisionId()) { - continue; + auto v0 = casc.v0_as(); + auto v0Element = v0.v0Data_as(); // V0 element from LF table containing V0 info + // V0 positive daughter + auto trackV0Dau0 = v0Element.posTrack_as(); // p <- V0 track (positive track) from MyTracks table + // V0 negative daughter + auto trackV0Dau1 = v0Element.negTrack_as(); // pion <- V0 track (negative track) from MyTracks table + + // check that particles come from the same collision + if (rejDiffCollTrack) { + if (trackV0Dau0.collisionId() != trackV0Dau1.collisionId()) { + continue; + } + if (trackXiDauCharged.collisionId() != trackV0Dau0.collisionId()) { + continue; + } } - } - - //--------------------------reconstruct V0 track--------------------------- - // pseudorapidity - double pseudorapV0PosDau = trackV0Dau0.eta(); - double pseudorapV0NegDau = trackV0Dau1.eta(); - - // pion & p <- V0 tracks - auto trackParCovV0Dau0 = getTrackParCov(trackV0Dau0); - auto trackParCovV0Dau1 = getTrackParCov(trackV0Dau1); - // info from LF table - std::array pVecV0 = {v0Element.px(), v0Element.py(), v0Element.pz()}; // pVec stands for vector containing the 3-momentum components - std::array vertexV0 = {v0Element.x(), v0Element.y(), v0Element.z()}; - const std::array covVtxV0 = {v0Element.positionCovMat()[0], v0Element.positionCovMat()[1], v0Element.positionCovMat()[2], v0Element.positionCovMat()[3], v0Element.positionCovMat()[4], v0Element.positionCovMat()[5]}; - - std::array pVecV0Dau0 = {casc.pxpos(), casc.pypos(), casc.pzpos()}; - std::array pVecV0Dau1 = {casc.pxneg(), casc.pyneg(), casc.pzneg()}; - - // create V0 track - auto trackV0 = o2::dataformats::V0(vertexV0, pVecV0, covVtxV0, trackParCovV0Dau0, trackParCovV0Dau1, {0, 0}, {0, 0}); - auto trackV0Copy = trackV0; - - //-----------------------------reconstruct cascade track----------------------------- - // pseudorapidity - double pseudorapPiFromCas = trackXiDauCharged.eta(); + //--------------------------reconstruct V0 track--------------------------- + // pseudorapidity + double pseudorapV0PosDau = trackV0Dau0.eta(); + double pseudorapV0NegDau = trackV0Dau1.eta(); + + // pion & p <- V0 tracks + auto trackParCovV0Dau0 = getTrackParCov(trackV0Dau0); + auto trackParCovV0Dau1 = getTrackParCov(trackV0Dau1); + + // info from LF table + std::array pVecV0 = {casc.pxlambda(), casc.pylambda(), casc.pzlambda()}; // pVec stands for vector containing the 3-momentum components + std::array vertexV0 = {casc.xlambda(), casc.ylambda(), casc.zlambda()}; + std::array covV0 = {0.}; + constexpr int MomInd[6] = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component + for (int i = 0; i < 6; i++) { + covV0[MomInd[i]] = v0Element.momentumCovMat()[i]; + covV0[i] = v0Element.positionCovMat()[i]; + } + // create V0 track + auto trackV0 = o2::track::TrackParCov(vertexV0, pVecV0, covV0, 0, true); + trackV0.setAbsCharge(0); + trackV0.setPID(o2::track::PID::Lambda); + + std::array pVecV0Dau0 = {casc.pxpos(), casc.pypos(), casc.pzpos()}; + std::array pVecV0Dau1 = {casc.pxneg(), casc.pyneg(), casc.pzneg()}; + + auto trackV0Copy = trackV0; + + //-----------------------------reconstruct cascade track----------------------------- + // pseudorapidity + double pseudorapPiFromCas = trackXiDauCharged.eta(); + + // pion <- casc track to be processed with DCAfitter + auto trackParCovXiDauCharged = getTrackParCov(trackXiDauCharged); + + // info from LF table + std::array vertexCasc = {casc.x(), casc.y(), casc.z()}; + std::array pVecCasc = {casc.px(), casc.py(), casc.pz()}; + std::array covCasc = {0.}; + for (int i = 0; i < 6; i++) { + covCasc[MomInd[i]] = casc.momentumCovMat()[i]; + covCasc[i] = casc.positionCovMat()[i]; + } + // create cascade track + o2::track::TrackParCov trackCasc; + if (trackXiDauCharged.sign() > 0) { + trackCasc = o2::track::TrackParCov(vertexCasc, pVecCasc, covCasc, 1, true); + } else if (trackXiDauCharged.sign() < 0) { + trackCasc = o2::track::TrackParCov(vertexCasc, pVecCasc, covCasc, -1, true); + } else { + continue; + } + trackCasc.setAbsCharge(1); + trackCasc.setPID(o2::track::PID::XiMinus); - // info from LF table - std::array vertexCasc = {casc.x(), casc.y(), casc.z()}; - std::array pVecCasc = {casc.px(), casc.py(), casc.pz()}; - const std::array covVtxCasc = {casc.positionCovMat()[0], casc.positionCovMat()[1], casc.positionCovMat()[2], casc.positionCovMat()[3], casc.positionCovMat()[4], casc.positionCovMat()[5]}; + std::array pVecPionFromCasc = {casc.pxbach(), casc.pybach(), casc.pzbach()}; - std::array pVecPionFromCasc = {casc.pxbach(), casc.pybach(), casc.pzbach()}; + auto trackCascCopy = trackCasc; - // pion <- casc track to be processed with DCAfitter - auto trackParVarXiDauCharged = getTrackParCov(trackXiDauCharged); + //-------------------combining cascade and pion tracks-------------------------- + auto groupedTrackIndices = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + for (auto const& trackIndexPion : groupedTrackIndices) { - // create cascade track - auto trackCasc = o2::dataformats::V0(vertexCasc, pVecCasc, covVtxCasc, trackV0, trackParVarXiDauCharged, {0, 0}, {0, 0}); - auto trackCascCopy = trackCasc; + auto trackPion = trackIndexPion.track_as(); - //-------------------combining cascade and pion tracks-------------------------- - for (auto const& trackPion : tracks) { + if ((rejDiffCollTrack) && (trackXiDauCharged.collisionId() != trackPion.collisionId())) { + continue; + } - if ((rejDiffCollTrack) && (trackXiDauCharged.collisionId() != trackPion.collisionId())) { // check to be further processed when the problem of ambiguous tracks will be solved - continue; - } + // ask for opposite sign daughters (omegac daughters) + if (trackPion.sign() * trackXiDauCharged.sign() >= 0) { + continue; + } - // ask for opposite sign daughters (omegac daughters) - if (trackPion.sign() * trackXiDauCharged.sign() >= 0) { - continue; - } + // check not to take the same particle twice in the decay chain + if (trackPion.globalIndex() == trackXiDauCharged.globalIndex() || trackPion.globalIndex() == trackV0Dau0.globalIndex() || trackPion.globalIndex() == trackV0Dau1.globalIndex() || trackPion.globalIndex() == casc.globalIndex()) { + continue; + } - // check not to take the same particle twice in the decay chain - if (trackPion.globalIndex() == trackXiDauCharged.globalIndex() || trackPion.globalIndex() == trackV0Dau0.globalIndex() || trackPion.globalIndex() == trackV0Dau1.globalIndex() || trackPion.globalIndex() == casc.globalIndex()) { - continue; - } + // pseudorapidity + double pseudorapPiFromOme = trackPion.eta(); - // pseudirapidity - double pseudorapPiFromOme = trackPion.eta(); + // primary pion track to be processed with DCAFitter + auto trackParVarPi = getTrackParCov(trackPion); + auto trackParVarPiCopy = trackParVarPi; - // primary pion track to be processed with DCAFitter - auto trackParVarPi = getTrackParCov(trackPion); - auto trackParVarPiCopy = trackParVarPi; + // reconstruct omegac with DCAFitter + int nVtxFromFitterOmegac = df.process(trackCasc, trackParVarPi); + if (nVtxFromFitterOmegac == 0) { + continue; + } + auto vertexOmegacFromFitter = df.getPCACandidate(); + auto chi2PCAOmegac = df.getChi2AtPCACandidate(); + std::array pVecCascAsD; + std::array pVecPionFromOmegac; + df.propagateTracksToVertex(); + if (!df.isPropagateTracksToVertexDone()) { + continue; + } + df.getTrack(0).getPxPyPzGlo(pVecCascAsD); + df.getTrack(1).getPxPyPzGlo(pVecPionFromOmegac); + std::array pVecOmegac = {pVecCascAsD[0] + pVecPionFromOmegac[0], pVecCascAsD[1] + pVecPionFromOmegac[1], pVecCascAsD[2] + pVecPionFromOmegac[2]}; + + std::array coordVtxOmegac = df.getPCACandidatePos(); + std::array covVtxOmegac = df.calcPCACovMatrixFlat(); + + // create omegac track + o2::track::TrackParCov trackOmegac = df.createParentTrackParCov(); + trackOmegac.setAbsCharge(0); + + // DCAxy (computed with propagateToDCABxByBz method) + float dcaxyV0Dau0 = trackV0Dau0.dcaXY(); + float dcaxyV0Dau1 = trackV0Dau1.dcaXY(); + float dcaxyPiFromCasc = trackXiDauCharged.dcaXY(); + + // DCAz (computed with propagateToDCABxByBz method) + float dcazV0Dau0 = trackV0Dau0.dcaZ(); + float dcazV0Dau1 = trackV0Dau1.dcaZ(); + float dcazPiFromCasc = trackXiDauCharged.dcaZ(); + + // primary vertex of the collision + auto primaryVertex = getPrimaryVertex(collision); // get the associated covariance matrix with auto covMatrixPV = primaryVertex.getCov(); + std::array pvCoord = {collision.posX(), collision.posY(), collision.posZ()}; + + if (doPvRefit && ((trackPion.pvRefitSigmaX2() != 1e10f) || (trackPion.pvRefitSigmaY2() != 1e10f) || (trackPion.pvRefitSigmaZ2() != 1e10f))) { // if I asked for PV refit in trackIndexSkimCreator.cxx + pvCoord[0] = trackPion.pvRefitX(); + pvCoord[1] = trackPion.pvRefitY(); + pvCoord[2] = trackPion.pvRefitZ(); + + // o2::dataformats::VertexBase Pvtx; + primaryVertex.setX(trackPion.pvRefitX()); + primaryVertex.setY(trackPion.pvRefitY()); + primaryVertex.setZ(trackPion.pvRefitZ()); + primaryVertex.setCov(trackPion.pvRefitSigmaX2(), trackPion.pvRefitSigmaXY(), trackPion.pvRefitSigmaY2(), trackPion.pvRefitSigmaXZ(), trackPion.pvRefitSigmaYZ(), trackPion.pvRefitSigmaZ2()); + + o2::dataformats::DCA impactParameterV0Dau0; + o2::dataformats::DCA impactParameterV0Dau1; + o2::dataformats::DCA impactParameterPiFromCasc; + o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackParCovV0Dau0, 2.f, matCorr, &impactParameterV0Dau0); + o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackParCovV0Dau1, 2.f, matCorr, &impactParameterV0Dau1); + o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackParCovXiDauCharged, 2.f, matCorr, &impactParameterPiFromCasc); + dcaxyV0Dau0 = impactParameterV0Dau0.getY(); + dcaxyV0Dau1 = impactParameterV0Dau1.getY(); + dcaxyPiFromCasc = impactParameterPiFromCasc.getY(); + dcazV0Dau0 = impactParameterV0Dau0.getZ(); + dcazV0Dau1 = impactParameterV0Dau1.getZ(); + dcazPiFromCasc = impactParameterPiFromCasc.getZ(); + } - // reconstruct omegac with DCAFitter - int nVtxFromFitterOmegac = df.process(trackCasc, trackParVarPi); - if (nVtxFromFitterOmegac == 0) { - continue; - } - auto vertexOmegacFromFitter = df.getPCACandidate(); - auto chi2PCAOmegac = df.getChi2AtPCACandidate(); - std::array pVecCascAsD; - std::array pVecPionFromOmegac; - df.propagateTracksToVertex(); - if (!df.isPropagateTracksToVertexDone()) { - continue; - } - df.getTrack(0).getPxPyPzGlo(pVecCascAsD); - df.getTrack(1).getPxPyPzGlo(pVecPionFromOmegac); - std::array pVecOmegac = {pVecCascAsD[0] + pVecPionFromOmegac[0], pVecCascAsD[1] + pVecPionFromOmegac[1], pVecCascAsD[2] + pVecPionFromOmegac[2]}; - - std::array coordVtxOmegac = df.getPCACandidatePos(); - std::array covVtxOmegac = df.calcPCACovMatrixFlat(); - - // create omegac track - auto trackOmegac = o2::dataformats::V0(coordVtxOmegac, pVecOmegac, covVtxOmegac, trackCasc, trackParVarPi, {0, 0}, {0, 0}); - - // impact parameter omegac - o2::dataformats::DCA impactParameterOmegac; - auto primaryVertex = getPrimaryVertex(collision); // get the associated covariance matrix with auto covMatrixPV = primaryVertex.getCov(); - trackOmegac.propagateToDCA(primaryVertex, magneticField, &impactParameterOmegac); - - // impact parameter - o2::dataformats::DCA impactParameterCasc; - o2::dataformats::DCA impactParameterPrimaryPi; - o2::dataformats::DCA impactParameterV0; - trackCascCopy.propagateToDCA(primaryVertex, magneticField, &impactParameterCasc); - trackParVarPiCopy.propagateToDCA(primaryVertex, magneticField, &impactParameterPrimaryPi); - trackV0Copy.propagateToDCA(primaryVertex, magneticField, &impactParameterV0); - - // DCAxy - double dcaxyPrimaryPi = trackPion.dcaXY(); - double dcaxyV0Dau0 = trackV0Dau0.dcaXY(); - double dcaxyV0Dau1 = trackV0Dau1.dcaXY(); - double dcaxyCascDau = trackXiDauCharged.dcaXY(); - - // invariant mass under the hypothesis of particles ID corresponding to the decay chain - double mLambda = v0Element.mLambda(); // from LF table, V0 mass under lambda hypothesis - double mAntiLambda = v0Element.mAntiLambda(); // from LF table, V0 mass under anti-lambda hypothesis - double mCasc = casc.mXi(); - const std::array arrMassOmegac = {massXiFromPDG, massPionFromPDG}; - double mOmegac = RecoDecay::m(std::array{pVecCascAsD, pVecPionFromOmegac}, arrMassOmegac); - - // computing cosPA - double cpaV0 = RecoDecay::cpa(vertexCasc, vertexV0, pVecV0); - double cpaOmegac = RecoDecay::cpa(std::array{collision.posX(), collision.posY(), collision.posZ()}, coordVtxOmegac, pVecOmegac); - double cpaCasc = RecoDecay::cpa(coordVtxOmegac, vertexCasc, pVecCasc); - double cpaxyV0 = RecoDecay::cpaXY(vertexCasc, vertexV0, pVecV0); - double cpaxyOmegac = RecoDecay::cpaXY(std::array{collision.posX(), collision.posY(), collision.posZ()}, coordVtxOmegac, pVecOmegac); - double cpaxyCasc = RecoDecay::cpaXY(coordVtxOmegac, vertexCasc, pVecCasc); - - // computing decay length and ctau - double decLenOmegac = RecoDecay::distance(std::array{collision.posX(), collision.posY(), collision.posZ()}, coordVtxOmegac); - double decLenCascade = RecoDecay::distance(coordVtxOmegac, vertexCasc); - double decLenV0 = RecoDecay::distance(vertexCasc, vertexV0); - double ctOmegac = RecoDecay::ct(pVecOmegac, decLenOmegac, massOmegacFromPDG); - double ctXic = RecoDecay::ct(pVecOmegac, decLenOmegac, massXicFromPDG); - double ctCascade = RecoDecay::ct(pVecCasc, decLenCascade, massXiFromPDG); - double ctV0 = RecoDecay::ct(pVecV0, decLenV0, massLambdaFromPDG); - - // computing eta - double pseudorapOmegac = RecoDecay::eta(pVecOmegac); - double pseudorapCascade = RecoDecay::eta(pVecCasc); - double pseudorapV0 = RecoDecay::eta(pVecV0); - - // DCA between daughters - double dcaCascDau = casc.dcacascdaughters(); - double dcaV0Dau = casc.dcaV0daughters(); - double dcaOmegacDau = std::sqrt(df.getChi2AtPCACandidate()); - - // set hfFlag - int hfFlag = 1 << DecayType::DecayToXiPi; - - // fill test histograms - - hInvMassOmegac->Fill(mOmegac); - - // fill the table - rowCandidate(collision.globalIndex(), - collision.posX(), collision.posY(), collision.posZ(), - vertexOmegacFromFitter[0], vertexOmegacFromFitter[1], vertexOmegacFromFitter[2], - vertexCasc[0], vertexCasc[1], vertexCasc[2], - vertexV0[0], vertexV0[1], vertexV0[2], - trackXiDauCharged.sign(), - chi2PCAOmegac, covVtxOmegac[0], covVtxOmegac[1], covVtxOmegac[2], covVtxOmegac[3], covVtxOmegac[4], covVtxOmegac[5], - covVtxV0[0], covVtxV0[1], covVtxV0[2], covVtxV0[3], covVtxV0[4], covVtxV0[5], - covVtxCasc[0], covVtxCasc[1], covVtxCasc[2], covVtxCasc[3], covVtxCasc[4], covVtxCasc[5], - pVecOmegac[0], pVecOmegac[1], pVecOmegac[2], - pVecCasc[0], pVecCasc[1], pVecCasc[2], - pVecPionFromOmegac[0], pVecPionFromOmegac[1], pVecPionFromOmegac[2], - pVecV0[0], pVecV0[1], pVecV0[2], - pVecPionFromCasc[0], pVecPionFromCasc[1], pVecPionFromCasc[2], - pVecV0Dau0[0], pVecV0Dau0[1], pVecV0Dau0[2], - pVecV0Dau1[0], pVecV0Dau1[1], pVecV0Dau1[2], - impactParameterCasc.getY(), impactParameterPrimaryPi.getY(), - impactParameterCasc.getZ(), impactParameterPrimaryPi.getZ(), - impactParameterV0.getY(), impactParameterV0.getZ(), - std::sqrt(impactParameterCasc.getSigmaY2()), std::sqrt(impactParameterPrimaryPi.getSigmaY2()), std::sqrt(impactParameterV0.getSigmaY2()), - v0Element.globalIndex(), v0Element.posTrackId(), v0Element.negTrackId(), - casc.globalIndex(), trackPion.globalIndex(), trackXiDauCharged.globalIndex(), - impactParameterOmegac.getY(), impactParameterOmegac.getZ(), - mLambda, mAntiLambda, mCasc, mOmegac, - cpaV0, cpaOmegac, cpaCasc, cpaxyV0, cpaxyOmegac, cpaxyCasc, - ctOmegac, ctCascade, ctV0, ctXic, - pseudorapV0PosDau, pseudorapV0NegDau, pseudorapPiFromCas, pseudorapPiFromOme, - pseudorapOmegac, pseudorapCascade, pseudorapV0, - dcaxyPrimaryPi, dcaxyV0Dau0, dcaxyV0Dau1, dcaxyCascDau, - dcaCascDau, dcaV0Dau, dcaOmegacDau, hfFlag); - - } // loop over pions - } // loop over candidates + // impact parameters + o2::dataformats::DCA impactParameterCasc; + o2::dataformats::DCA impactParameterPrimaryPi; + o2::dataformats::DCA impactParameterV0; + o2::dataformats::DCA impactParameterOmegac; + o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackCascCopy, 2.f, matCorr, &impactParameterCasc); + o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackParVarPiCopy, 2.f, matCorr, &impactParameterPrimaryPi); + o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackV0Copy, 2.f, matCorr, &impactParameterV0); + o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVertex, trackOmegac, 2.f, matCorr, &impactParameterOmegac); + + // invariant mass under the hypothesis of particles ID corresponding to the decay chain + double mLambda = casc.mLambda(); // from LF table, V0 mass under lambda hypothesis + double mCasc = casc.mXi(); + const std::array arrMassOmegac = {massXiFromPDG, massPionFromPDG}; + double mOmegac = RecoDecay::m(std::array{pVecCascAsD, pVecPionFromOmegac}, arrMassOmegac); + + // computing cosPA + double cpaV0 = RecoDecay::cpa(vertexCasc, vertexV0, pVecV0); + double cpaOmegac = RecoDecay::cpa(pvCoord, coordVtxOmegac, pVecOmegac); + double cpaCasc = RecoDecay::cpa(coordVtxOmegac, vertexCasc, pVecCasc); + double cpaxyV0 = RecoDecay::cpaXY(vertexCasc, vertexV0, pVecV0); + double cpaxyOmegac = RecoDecay::cpaXY(pvCoord, coordVtxOmegac, pVecOmegac); + double cpaxyCasc = RecoDecay::cpaXY(coordVtxOmegac, vertexCasc, pVecCasc); + + // computing decay length and ctau + double decLenOmegac = RecoDecay::distance(pvCoord, coordVtxOmegac); + double decLenCascade = RecoDecay::distance(coordVtxOmegac, vertexCasc); + double decLenV0 = RecoDecay::distance(vertexCasc, vertexV0); + double ctOmegac = RecoDecay::ct(pVecOmegac, decLenOmegac, massOmegacFromPDG); + double ctXic = RecoDecay::ct(pVecOmegac, decLenOmegac, massXicFromPDG); + double ctCascade = RecoDecay::ct(pVecCasc, decLenCascade, massXiFromPDG); + double ctV0 = RecoDecay::ct(pVecV0, decLenV0, massLambdaFromPDG); + + // computing eta + double pseudorapOmegac = RecoDecay::eta(pVecOmegac); + double pseudorapCascade = RecoDecay::eta(pVecCasc); + double pseudorapV0 = RecoDecay::eta(pVecV0); + + // DCA between daughters + float dcaCascDau = casc.dcacascdaughters(); + float dcaV0Dau = casc.dcaV0daughters(); + float dcaOmegacDau = std::sqrt(df.getChi2AtPCACandidate()); + + // set hfFlag + int hfFlag = 1 << DecayType::DecayToXiPi; + + // fill test histograms + hInvMassOmegac->Fill(mOmegac); + + // fill the table + rowCandidate(collision.globalIndex(), + pvCoord[0], pvCoord[1], pvCoord[2], + vertexOmegacFromFitter[0], vertexOmegacFromFitter[1], vertexOmegacFromFitter[2], + vertexCasc[0], vertexCasc[1], vertexCasc[2], + vertexV0[0], vertexV0[1], vertexV0[2], + trackXiDauCharged.sign(), + chi2PCAOmegac, covVtxOmegac[0], covVtxOmegac[1], covVtxOmegac[2], covVtxOmegac[3], covVtxOmegac[4], covVtxOmegac[5], + covV0[0], covV0[1], covV0[2], covV0[3], covV0[4], covV0[5], + covCasc[0], covCasc[1], covCasc[2], covCasc[3], covCasc[4], covCasc[5], + pVecOmegac[0], pVecOmegac[1], pVecOmegac[2], + pVecCasc[0], pVecCasc[1], pVecCasc[2], + pVecPionFromOmegac[0], pVecPionFromOmegac[1], pVecPionFromOmegac[2], + pVecV0[0], pVecV0[1], pVecV0[2], + pVecPionFromCasc[0], pVecPionFromCasc[1], pVecPionFromCasc[2], + pVecV0Dau0[0], pVecV0Dau0[1], pVecV0Dau0[2], + pVecV0Dau1[0], pVecV0Dau1[1], pVecV0Dau1[2], + impactParameterCasc.getY(), impactParameterPrimaryPi.getY(), + impactParameterCasc.getZ(), impactParameterPrimaryPi.getZ(), + impactParameterV0.getY(), impactParameterV0.getZ(), + std::sqrt(impactParameterCasc.getSigmaY2()), std::sqrt(impactParameterPrimaryPi.getSigmaY2()), std::sqrt(impactParameterV0.getSigmaY2()), + v0Element.globalIndex(), v0Element.posTrackId(), v0Element.negTrackId(), + casc.globalIndex(), trackPion.globalIndex(), trackXiDauCharged.globalIndex(), + impactParameterOmegac.getY(), impactParameterOmegac.getZ(), + mLambda, mCasc, mOmegac, + cpaV0, cpaOmegac, cpaCasc, cpaxyV0, cpaxyOmegac, cpaxyCasc, + ctOmegac, ctCascade, ctV0, ctXic, + pseudorapV0PosDau, pseudorapV0NegDau, pseudorapPiFromCas, pseudorapPiFromOme, + pseudorapOmegac, pseudorapCascade, pseudorapV0, + dcaxyV0Dau0, dcaxyV0Dau1, dcaxyPiFromCasc, + dcazV0Dau0, dcazV0Dau1, dcazPiFromCasc, + dcaCascDau, dcaV0Dau, dcaOmegacDau, hfFlag); + + } // loop over pions + } // loop over cascades + } // close loop collisions } // end of process }; // end of struct diff --git a/PWGHF/TableProducer/candidateSelectorToXiPi.cxx b/PWGHF/TableProducer/candidateSelectorToXiPi.cxx index 22ea4d13384..85cea7ed51d 100644 --- a/PWGHF/TableProducer/candidateSelectorToXiPi.cxx +++ b/PWGHF/TableProducer/candidateSelectorToXiPi.cxx @@ -52,9 +52,16 @@ struct HfCandidateSelectorToXiPi { Configurable etaTrackMax{"etaTrackMax", 0.8, "Max absolute value of eta"}; Configurable ptPiFromCascMin{"ptPiFromCascMin", 0.15, "Min pT pi <- casc"}; Configurable ptPiFromOmeMin{"ptPiFromOmeMin", 0.2, "Min pT pi <- omegac"}; - Configurable dcaXYPriPiMin{"dcaXYPriPiMin", 0., "Min dcaxy primary pi track to PV"}; - Configurable dcaXYPriPiMax{"dcaXYPriPiMax", 10., "Max dcaxy primary pi track to PV"}; - Configurable impactParameterZPriPiMax{"impactParameterZPriPiMax", 10., "Max impact parameter primary pi track to PV in Z direction"}; + + Configurable impactParameterXYPriPiMin{"impactParameterXYPriPiMin", 0., "Min dcaxy primary pi track to PV"}; + Configurable impactParameterXYPriPiMax{"impactParameterXYPriPiMax", 10., "Max dcaxy primary pi track to PV"}; + Configurable impactParameterZPriPiMin{"impactParameterZPriPiMin", 0., "Min dcaz primary pi track to PV"}; + Configurable impactParameterZPriPiMax{"impactParameterZPriPiMax", 10., "Max dcaz primary pi track to PV"}; + + Configurable impactParameterXYCascMin{"impactParameterXYCascMin", 0., "Min dcaxy cascade track to PV"}; + Configurable impactParameterXYCascMax{"impactParameterXYCascMax", 10., "Max dcaxy cascade track to PV"}; + Configurable impactParameterZCascMin{"impactParameterZCascMin", 0., "Min dcaz cascade track to PV"}; + Configurable impactParameterZCascMax{"impactParameterZCascMax", 10., "Max dcaz cascade track to PV"}; Configurable ptCandMin{"ptCandMin", 0., "Lower bound of candidate pT"}; Configurable ptCandMax{"ptCandMax", 50., "Upper bound of candidate pT"}; @@ -97,7 +104,7 @@ struct HfCandidateSelectorToXiPi { Configurable nClustersItsMin{"nClustersItsMin", 3, "Minimum number of ITS clusters requirement for pi <- Omegac"}; Configurable nClustersItsInnBarrMin{"nClustersItsInnBarrMin", 1, "Minimum number of ITS clusters in inner barrel requirement for pi <- Omegac"}; - using MyTrackInfo = soa::Join; + using MyTrackInfo = aod::BigTracksPIDExtended; OutputObj hxVertexOmegac{TH1F("hxVertexOmegac", "x Omegac vertex;xVtx;entries", 500, -10, 10)}; OutputObj hInvMassOmegac{TH1F("hInvMassOmegac", "Omegac invariant mass;inv mass;entries", 500, 2.2, 3.1)}; @@ -131,8 +138,7 @@ struct HfCandidateSelectorToXiPi { double massLambdaFromPDG = RecoDecay::getMassPDG(kLambda0); double massXiFromPDG = RecoDecay::getMassPDG(kXiMinus); - int collIdCasc = -999; - int collIdPrimaryPi = -999; + int collId = -999; // looping over omegac candidates for (auto const& candidate : candidates) { @@ -201,11 +207,19 @@ struct HfCandidateSelectorToXiPi { resultSelections = false; } - // cut on primary pion dcaXY and impact parameter - if ((candidate.dcaXYToPVPrimaryPi() < dcaXYPriPiMin) || (candidate.dcaXYToPVPrimaryPi() > dcaXYPriPiMax)) { + // cut on primary pion dcaXY and dcaZ + if ((candidate.impactParPrimaryPiXY() < impactParameterXYPriPiMin) || (candidate.impactParPrimaryPiXY() > impactParameterXYPriPiMax)) { + resultSelections = false; + } + if ((candidate.impactParPrimaryPiZ() < impactParameterZPriPiMin) || (candidate.impactParPrimaryPiZ() > impactParameterZPriPiMax)) { + resultSelections = false; + } + + // cut on cascade dcaXY and dcaZ + if ((candidate.impactParCascXY() < impactParameterXYCascMin) || (candidate.impactParCascXY() > impactParameterXYCascMax)) { resultSelections = false; } - if (candidate.impactParPrimaryPiZ() > impactParameterZPriPiMax) { + if ((candidate.impactParCascZ() < impactParameterZCascMin) || (candidate.impactParCascZ() > impactParameterZCascMax)) { resultSelections = false; } @@ -311,12 +325,7 @@ struct HfCandidateSelectorToXiPi { bool statusInvMassCascade = false; bool statusInvMassOmegac = false; - double invMassLambda = 0; - if (signDecay < 0) { - invMassLambda = candidate.invMassLambda(); - } else if (signDecay > 0) { - invMassLambda = candidate.invMassAntiLambda(); - } + double invMassLambda = candidate.invMassLambda(); double invMassCascade = candidate.invMassCascade(); double invMassOmegac = candidate.invMassOmegac(); @@ -388,16 +397,9 @@ struct HfCandidateSelectorToXiPi { hCTauOmegac->Fill(candidate.ctauOmegac()); hCTauXic->Fill(candidate.ctauXic()); - if (trackPiFromCasc.collisionId() != collIdCasc) { + if (candidate.collisionId() != collId) { hNEventsSaved->Fill(0.5); - collIdCasc = trackPiFromCasc.collisionId(); - } - if (trackPiFromOmeg.collisionId() != collIdPrimaryPi) { - hNEventsSaved->Fill(1.5); - collIdPrimaryPi = trackPiFromOmeg.collisionId(); - } - if (trackPiFromOmeg.collisionId() != trackPiFromCasc.collisionId()) { - hNEventsSaved->Fill(2.5); + collId = trackPiFromCasc.collisionId(); } } }