From b7d95913a5ac46a60df61000c9b32435b8608933 Mon Sep 17 00:00:00 2001 From: Alexandre Bigot Date: Fri, 24 Mar 2023 17:39:40 +0100 Subject: [PATCH 1/5] Change structure: loop over aod::Collisions, over D candidates attached to this collision and over tracks attached to this collision --- PWGHF/TableProducer/candidateCreatorB0.cxx | 245 ++++++++++++--------- 1 file changed, 136 insertions(+), 109 deletions(-) diff --git a/PWGHF/TableProducer/candidateCreatorB0.cxx b/PWGHF/TableProducer/candidateCreatorB0.cxx index cd05108995a..4e27bf2a0de 100644 --- a/PWGHF/TableProducer/candidateCreatorB0.cxx +++ b/PWGHF/TableProducer/candidateCreatorB0.cxx @@ -19,6 +19,7 @@ #include "Framework/AnalysisTask.h" #include "DCAFitter/DCAFitterN.h" #include "Common/Core/trackUtilities.h" +#include "Common/DataModel/CollisionAssociation.h" #include "ReconstructionDataFormats/DCA.h" #include "ReconstructionDataFormats/V0.h" #include "PWGHF/DataModel/CandidateReconstructionTables.h" @@ -76,9 +77,19 @@ struct HfCandidateCreatorB0 { using TracksWithSel = soa::Join; + // Filter filterSelectCollisions = (aod::hf_sel_collision::whyRejectColl == 0); // select only triggered events + // using SelectedCollisions = soa::Filtered>; Filter filterSelectTracks = (!usePionIsGlobalTrackWoDCA) || ((usePionIsGlobalTrackWoDCA) && requireGlobalTrackWoDCAInFilter()); Filter filterSelectCandidates = (aod::hf_sel_candidate_dplus::isSelDplusToPiKPi >= selectionFlagD); + using TracksWithSelFiltered = soa::Filtered; + + using CandsDFiltered = soa::Filtered>; + Preslice candsDPerCollision = aod::track_association::collisionId; + + //using TrackAssocSelFiltered = soa::Filtered; // soa::Filtered>; + Preslice trackIndicesPerCollision = aod::track_association::collisionId; + OutputObj hMassDToPiKPi{TH1F("hMassDToPiKPi", "D^{#minus} candidates;inv. mass (p^{#minus} K^{#plus} #pi^{#minus}) (GeV/#it{c}^{2});entries", 500, 0., 5.)}; OutputObj hPtD{TH1F("hPtD", "D^{#minus} candidates;D^{#minus} candidate #it{p}_{T} (GeV/#it{c});entries", 100, 0., 10.)}; OutputObj hPtPion{TH1F("hPtPion", "#pi^{#plus} candidates;#pi^{#plus} candidate #it{p}_{T} (GeV/#it{c});entries", 100, 0., 10.)}; @@ -107,10 +118,12 @@ struct HfCandidateCreatorB0 { return true; } - void process(aod::Collision const& collision, - soa::Filtered> const& candsD, + // TODO: filter Collisions table to only keep triggered events ? + void process(aod::Collisions const& collisions, + CandsDFiltered const& candsD, + aod::TrackAssoc const& trackIndices, TracksWithSel const&, - soa::Filtered const& tracksPion) + TracksWithSelFiltered const&) { // FIXME: store B0 creator configurable (until https://alice.its.cern.ch/jira/browse/O2-3582 solved) if (!isHfCandB0ConfigFilled) { @@ -140,117 +153,131 @@ struct HfCandidateCreatorB0 { df3.setUseAbsDCA(useAbsDCA); df3.setWeightedFinalPCA(useWeightedFinalPCA); - // loop over D candidates - for (const auto& candD : candsD) { - hMassDToPiKPi->Fill(invMassDplusToPiKPi(candD), candD.pt()); - hPtD->Fill(candD.pt()); - hCPAD->Fill(candD.cpa()); - - // track0 <-> pi, track1 <-> K, track2 <-> pi - auto track0 = candD.prong0_as(); - auto track1 = candD.prong1_as(); - auto track2 = candD.prong2_as(); - auto trackParCov0 = getTrackParCov(track0); - auto trackParCov1 = getTrackParCov(track1); - auto trackParCov2 = getTrackParCov(track2); - - // reconstruct 3-prong secondary vertex (D±) - if (df3.process(trackParCov0, trackParCov1, trackParCov2) == 0) { - continue; - } + static int ncol = 0; - const auto& secondaryVertex = df3.getPCACandidate(); - trackParCov0.propagateTo(secondaryVertex[0], bz); - trackParCov1.propagateTo(secondaryVertex[0], bz); - trackParCov2.propagateTo(secondaryVertex[0], bz); - - // D∓ → π∓ K± π∓ - array pVecpiK = {track0.px() + track1.px(), track0.py() + track1.py(), track0.pz() + track1.pz()}; - array pVecD = {pVecpiK[0] + track2.px(), pVecpiK[1] + track2.py(), pVecpiK[2] + track2.pz()}; - auto trackParCovPiK = o2::dataformats::V0(df3.getPCACandidatePos(), pVecpiK, df3.calcPCACovMatrixFlat(), - trackParCov0, trackParCov1, {0, 0}, {0, 0}); - auto trackParCovD = o2::dataformats::V0(df3.getPCACandidatePos(), pVecD, df3.calcPCACovMatrixFlat(), - trackParCovPiK, trackParCov2, {0, 0}, {0, 0}); - - int index0D = track0.globalIndex(); - int index1D = track1.globalIndex(); - int index2D = track2.globalIndex(); - - // loop over pions - for (const auto& trackPion : tracksPion) { - // minimum pT selection - if (trackPion.pt() < ptPionMin || !isSelectedTrackDCA(trackPion)) { - continue; - } - // reject pions that are D daughters - if (trackPion.globalIndex() == index0D || trackPion.globalIndex() == index1D || trackPion.globalIndex() == index2D) { - continue; - } - // reject pi and D with same sign - if (trackPion.sign() * track0.sign() > 0) { - continue; - } + for (const auto& collision : collisions) { + auto primaryVertex = getPrimaryVertex(collision); + auto covMatrixPV = primaryVertex.getCov(); - hPtPion->Fill(trackPion.pt()); - array pVecPion = {trackPion.px(), trackPion.py(), trackPion.pz()}; - auto trackParCovPi = getTrackParCov(trackPion); - // --------------------------------- - // reconstruct the 2-prong B0 vertex - if (df2.process(trackParCovD, trackParCovPi) == 0) { + if (ncol % 10000 == 0) { + LOG(info) << ncol << " collisions parsed"; + } + ncol++; + + auto thisCollId = collision.globalIndex(); + auto candsDThisColl = candsD.sliceBy(candsDPerCollision, thisCollId); + + for (const auto& candD : candsDThisColl) { // start loop over filtered D candidates indices as associated to this collision in candidateCreator3Prong.cxx + hMassDToPiKPi->Fill(invMassDplusToPiKPi(candD), candD.pt()); + hPtD->Fill(candD.pt()); + hCPAD->Fill(candD.cpa()); + + // track0 <-> pi, track1 <-> K, track2 <-> pi + auto track0 = candD.prong0_as(); + auto track1 = candD.prong1_as(); + auto track2 = candD.prong2_as(); + auto trackParCov0 = getTrackParCov(track0); + auto trackParCov1 = getTrackParCov(track1); + auto trackParCov2 = getTrackParCov(track2); + + // reconstruct 3-prong secondary vertex (D±) + if (df3.process(trackParCov0, trackParCov1, trackParCov2) == 0) { continue; } - // calculate invariant mass - massDPi = RecoDecay::m(array{pVecD, pVecPion}, array{massD, massPi}); - if (std::abs(massDPi - massB0) > invMassWindowB0) { - continue; - } - hMassB0ToDPi->Fill(massDPi); - - // calculate relevant properties - const auto& secondaryVertexB0 = df2.getPCACandidate(); - auto chi2PCA = df2.getChi2AtPCACandidate(); - auto covMatrixPCA = df2.calcPCACovMatrixFlat(); - - // must be called before getTrack query - df2.propagateTracksToVertex(); - // track.getPxPyPzGlo(pVec) modifies pVec of track - df2.getTrack(0).getPxPyPzGlo(pVecD); - df2.getTrack(1).getPxPyPzGlo(pVecPion); - - auto primaryVertex = getPrimaryVertex(collision); - auto covMatrixPV = primaryVertex.getCov(); - o2::dataformats::DCA impactParameter0; - o2::dataformats::DCA impactParameter1; - trackParCovD.propagateToDCA(primaryVertex, bz, &impactParameter0); - trackParCovPi.propagateToDCA(primaryVertex, bz, &impactParameter1); - - hCovSVXX->Fill(covMatrixPCA[0]); - hCovPVXX->Fill(covMatrixPV[0]); - - // get uncertainty of the decay length - double phi, theta; - // getPointDirection modifies phi and theta - getPointDirection(array{collision.posX(), collision.posY(), collision.posZ()}, secondaryVertexB0, phi, theta); - auto errorDecayLength = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta)); - auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.)); - - int hfFlag = BIT(hf_cand_b0::DecayType::B0ToDPi); - - // fill the candidate table for the B0 here: - rowCandidateBase(collision.globalIndex(), - collision.posX(), collision.posY(), collision.posZ(), - secondaryVertexB0[0], secondaryVertexB0[1], secondaryVertexB0[2], - errorDecayLength, errorDecayLengthXY, - chi2PCA, - pVecD[0], pVecD[1], pVecD[2], - pVecPion[0], pVecPion[1], pVecPion[2], - impactParameter0.getY(), impactParameter1.getY(), - std::sqrt(impactParameter0.getSigmaY2()), std::sqrt(impactParameter1.getSigmaY2()), - candD.globalIndex(), trackPion.globalIndex(), - hfFlag); - } // pi loop - } // D loop + const auto& secondaryVertex = df3.getPCACandidate(); + trackParCov0.propagateTo(secondaryVertex[0], bz); + trackParCov1.propagateTo(secondaryVertex[0], bz); + trackParCov2.propagateTo(secondaryVertex[0], bz); + + // D∓ → π∓ K± π∓ + array pVecpiK = {track0.px() + track1.px(), track0.py() + track1.py(), track0.pz() + track1.pz()}; + array pVecD = {pVecpiK[0] + track2.px(), pVecpiK[1] + track2.py(), pVecpiK[2] + track2.pz()}; + auto trackParCovPiK = o2::dataformats::V0(df3.getPCACandidatePos(), pVecpiK, df3.calcPCACovMatrixFlat(), + trackParCov0, trackParCov1, {0, 0}, {0, 0}); + auto trackParCovD = o2::dataformats::V0(df3.getPCACandidatePos(), pVecD, df3.calcPCACovMatrixFlat(), + trackParCovPiK, trackParCov2, {0, 0}, {0, 0}); + + int index0D = track0.globalIndex(); + int index1D = track1.globalIndex(); + int index2D = track2.globalIndex(); + + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + for (const auto& trackId : trackIdsThisCollision) { // start loop over track indices associated to this collision + auto trackPion = trackId.track_as(); // TODO: test with TracksWithSelFiltereds + + // minimum pT selection + if (trackPion.pt() < ptPionMin || !isSelectedTrackDCA(trackPion)) { + continue; + } + // reject pions that are D daughters + if (trackPion.globalIndex() == index0D || trackPion.globalIndex() == index1D || trackPion.globalIndex() == index2D) { + continue; + } + // reject pi and D with same sign + if (trackPion.sign() * track0.sign() > 0) { + continue; + } + + hPtPion->Fill(trackPion.pt()); + array pVecPion = {trackPion.px(), trackPion.py(), trackPion.pz()}; + auto trackParCovPi = getTrackParCov(trackPion); + // --------------------------------- + // reconstruct the 2-prong B0 vertex + if (df2.process(trackParCovD, trackParCovPi) == 0) { + continue; + } + + // calculate invariant mass + massDPi = RecoDecay::m(array{pVecD, pVecPion}, array{massD, massPi}); + if (std::abs(massDPi - massB0) > invMassWindowB0) { + continue; + } + hMassB0ToDPi->Fill(massDPi); + + // calculate relevant properties + const auto& secondaryVertexB0 = df2.getPCACandidate(); + auto chi2PCA = df2.getChi2AtPCACandidate(); + auto covMatrixPCA = df2.calcPCACovMatrixFlat(); + + // must be called before getTrack query + df2.propagateTracksToVertex(); + // track.getPxPyPzGlo(pVec) modifies pVec of track + df2.getTrack(0).getPxPyPzGlo(pVecD); + df2.getTrack(1).getPxPyPzGlo(pVecPion); + + o2::dataformats::DCA impactParameter0; + o2::dataformats::DCA impactParameter1; + trackParCovD.propagateToDCA(primaryVertex, bz, &impactParameter0); + trackParCovPi.propagateToDCA(primaryVertex, bz, &impactParameter1); + + hCovSVXX->Fill(covMatrixPCA[0]); + hCovPVXX->Fill(covMatrixPV[0]); + + // get uncertainty of the decay length + double phi, theta; + // getPointDirection modifies phi and theta + getPointDirection(array{collision.posX(), collision.posY(), collision.posZ()}, secondaryVertexB0, phi, theta); + auto errorDecayLength = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta)); + auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.)); + + int hfFlag = BIT(hf_cand_b0::DecayType::B0ToDPi); + + // fill the candidate table for the B0 here: + rowCandidateBase(collision.globalIndex(), + collision.posX(), collision.posY(), collision.posZ(), + secondaryVertexB0[0], secondaryVertexB0[1], secondaryVertexB0[2], + errorDecayLength, errorDecayLengthXY, + chi2PCA, + pVecD[0], pVecD[1], pVecD[2], + pVecPion[0], pVecPion[1], pVecPion[2], + impactParameter0.getY(), impactParameter1.getY(), + std::sqrt(impactParameter0.getSigmaY2()), std::sqrt(impactParameter1.getSigmaY2()), + candD.globalIndex(), trackPion.globalIndex(), + hfFlag); + } // pi loop + } // D loop + } // collision loop } // process }; // struct From 8fa08c3f201b64fa077ea077da102159d0d30713 Mon Sep 17 00:00:00 2001 From: Alexandre Bigot Date: Mon, 27 Mar 2023 15:07:38 +0200 Subject: [PATCH 2/5] Repropagate 3prong daughters to this collision if needed + Apply fixes --- PWGHF/TableProducer/candidateCreatorB0.cxx | 121 +++++++++++++-------- 1 file changed, 75 insertions(+), 46 deletions(-) diff --git a/PWGHF/TableProducer/candidateCreatorB0.cxx b/PWGHF/TableProducer/candidateCreatorB0.cxx index 4e27bf2a0de..21b92cf940f 100644 --- a/PWGHF/TableProducer/candidateCreatorB0.cxx +++ b/PWGHF/TableProducer/candidateCreatorB0.cxx @@ -77,8 +77,6 @@ struct HfCandidateCreatorB0 { using TracksWithSel = soa::Join; - // Filter filterSelectCollisions = (aod::hf_sel_collision::whyRejectColl == 0); // select only triggered events - // using SelectedCollisions = soa::Filtered>; Filter filterSelectTracks = (!usePionIsGlobalTrackWoDCA) || ((usePionIsGlobalTrackWoDCA) && requireGlobalTrackWoDCAInFilter()); Filter filterSelectCandidates = (aod::hf_sel_candidate_dplus::isSelDplusToPiKPi >= selectionFlagD); @@ -87,7 +85,6 @@ struct HfCandidateCreatorB0 { using CandsDFiltered = soa::Filtered>; Preslice candsDPerCollision = aod::track_association::collisionId; - //using TrackAssocSelFiltered = soa::Filtered; // soa::Filtered>; Preslice trackIndicesPerCollision = aod::track_association::collisionId; OutputObj hMassDToPiKPi{TH1F("hMassDToPiKPi", "D^{#minus} candidates;inv. mass (p^{#minus} K^{#plus} #pi^{#minus}) (GeV/#it{c}^{2});entries", 500, 0., 5.)}; @@ -118,7 +115,6 @@ struct HfCandidateCreatorB0 { return true; } - // TODO: filter Collisions table to only keep triggered events ? void process(aod::Collisions const& collisions, CandsDFiltered const& candsD, aod::TrackAssoc const& trackIndices, @@ -180,38 +176,70 @@ struct HfCandidateCreatorB0 { auto trackParCov1 = getTrackParCov(track1); auto trackParCov2 = getTrackParCov(track2); + std::array pVec0; // = {track0.px(), track0.py(), track0.pz()}; + std::array pVec1; // = {track1.px(), track1.py(), track1.pz()}; + std::array pVec2; // = {track2.px(), track2.py(), track2.pz()}; + + o2::dataformats::DCA dca0; + // dca0.set(track0.dcaXY(), track0.dcaZ()); + o2::dataformats::DCA dca1; + // dca1.set(track1.dcaXY(), track1.dcaZ()); + o2::dataformats::DCA dca2; + // dca2.set(track2.dcaXY(), track2.dcaZ()); + + // repropagate tracks to this collision if needed + if (track0.collisionId() != thisCollId) { + trackParCov0.propagateToDCA(primaryVertex, bz, &dca0); + } + + if (track1.collisionId() != thisCollId) { + trackParCov1.propagateToDCA(primaryVertex, bz, &dca1); + } + + if (track2.collisionId() != thisCollId) { + trackParCov2.propagateToDCA(primaryVertex, bz, &dca2); + } + + // --------------------------------- // reconstruct 3-prong secondary vertex (D±) if (df3.process(trackParCov0, trackParCov1, trackParCov2) == 0) { continue; } - const auto& secondaryVertex = df3.getPCACandidate(); - trackParCov0.propagateTo(secondaryVertex[0], bz); - trackParCov1.propagateTo(secondaryVertex[0], bz); - trackParCov2.propagateTo(secondaryVertex[0], bz); + const auto& secondaryVertexD = df3.getPCACandidate(); + // propagate the 3 prongs to the secondary vertex + trackParCov0.propagateTo(secondaryVertexD[0], bz); + trackParCov1.propagateTo(secondaryVertexD[0], bz); + trackParCov2.propagateTo(secondaryVertexD[0], bz); + + // update pVec of tracks + getPxPyPz(trackParCov0, pVec0); // FIXME: or df3.getTrack(0).getPxPyPzGlo(pVec0) ? + getPxPyPz(trackParCov1, pVec1); + getPxPyPz(trackParCov2, pVec2); // D∓ → π∓ K± π∓ - array pVecpiK = {track0.px() + track1.px(), track0.py() + track1.py(), track0.pz() + track1.pz()}; - array pVecD = {pVecpiK[0] + track2.px(), pVecpiK[1] + track2.py(), pVecpiK[2] + track2.pz()}; - auto trackParCovPiK = o2::dataformats::V0(df3.getPCACandidatePos(), pVecpiK, df3.calcPCACovMatrixFlat(), + array pVecPiK = RecoDecay::pVec(pVec0, pVec1); + array pVecD = RecoDecay::pVec(pVec0, pVec1, pVec2); + auto trackParCovPiK = o2::dataformats::V0(df3.getPCACandidatePos(), pVecPiK, df3.calcPCACovMatrixFlat(), trackParCov0, trackParCov1, {0, 0}, {0, 0}); auto trackParCovD = o2::dataformats::V0(df3.getPCACandidatePos(), pVecD, df3.calcPCACovMatrixFlat(), trackParCovPiK, trackParCov2, {0, 0}, {0, 0}); - int index0D = track0.globalIndex(); - int index1D = track1.globalIndex(); - int index2D = track2.globalIndex(); + int indexTrack0 = track0.globalIndex(); + int indexTrack1 = track1.globalIndex(); + int indexTrack2 = track2.globalIndex(); auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + for (const auto& trackId : trackIdsThisCollision) { // start loop over track indices associated to this collision - auto trackPion = trackId.track_as(); // TODO: test with TracksWithSelFiltereds - + auto trackPion = trackId.track_as(); + // minimum pT selection if (trackPion.pt() < ptPionMin || !isSelectedTrackDCA(trackPion)) { continue; } // reject pions that are D daughters - if (trackPion.globalIndex() == index0D || trackPion.globalIndex() == index1D || trackPion.globalIndex() == index2D) { + if (trackPion.globalIndex() == indexTrack0 || trackPion.globalIndex() == indexTrack1 || trackPion.globalIndex() == indexTrack2) { continue; } // reject pi and D with same sign @@ -222,37 +250,38 @@ struct HfCandidateCreatorB0 { hPtPion->Fill(trackPion.pt()); array pVecPion = {trackPion.px(), trackPion.py(), trackPion.pz()}; auto trackParCovPi = getTrackParCov(trackPion); + // --------------------------------- // reconstruct the 2-prong B0 vertex if (df2.process(trackParCovD, trackParCovPi) == 0) { continue; } - // calculate invariant mass - massDPi = RecoDecay::m(array{pVecD, pVecPion}, array{massD, massPi}); - if (std::abs(massDPi - massB0) > invMassWindowB0) { - continue; - } - hMassB0ToDPi->Fill(massDPi); - // calculate relevant properties const auto& secondaryVertexB0 = df2.getPCACandidate(); auto chi2PCA = df2.getChi2AtPCACandidate(); auto covMatrixPCA = df2.calcPCACovMatrixFlat(); + hCovSVXX->Fill(covMatrixPCA[0]); + hCovPVXX->Fill(covMatrixPV[0]); - // must be called before getTrack query + // propagate D and Pi to the B0 vertex df2.propagateTracksToVertex(); // track.getPxPyPzGlo(pVec) modifies pVec of track - df2.getTrack(0).getPxPyPzGlo(pVecD); - df2.getTrack(1).getPxPyPzGlo(pVecPion); + df2.getTrack(0).getPxPyPzGlo(pVecD); // momentum of D at the B0 vertex + df2.getTrack(1).getPxPyPzGlo(pVecPion); // momentum of Pi at the B0 vertex - o2::dataformats::DCA impactParameter0; - o2::dataformats::DCA impactParameter1; - trackParCovD.propagateToDCA(primaryVertex, bz, &impactParameter0); - trackParCovPi.propagateToDCA(primaryVertex, bz, &impactParameter1); + // calculate invariant mass and apply selection + massDPi = RecoDecay::m(array{pVecD, pVecPion}, array{massD, massPi}); + if (std::abs(massDPi - massB0) > invMassWindowB0) { + continue; + } + hMassB0ToDPi->Fill(massDPi); - hCovSVXX->Fill(covMatrixPCA[0]); - hCovPVXX->Fill(covMatrixPV[0]); + // compute impact parameters of D and Pi + o2::dataformats::DCA dcaD; + o2::dataformats::DCA dcaPion; + trackParCovD.propagateToDCA(primaryVertex, bz, &dcaD); + trackParCovPi.propagateToDCA(primaryVertex, bz, &dcaPion); // get uncertainty of the decay length double phi, theta; @@ -265,21 +294,21 @@ struct HfCandidateCreatorB0 { // fill the candidate table for the B0 here: rowCandidateBase(collision.globalIndex(), - collision.posX(), collision.posY(), collision.posZ(), - secondaryVertexB0[0], secondaryVertexB0[1], secondaryVertexB0[2], - errorDecayLength, errorDecayLengthXY, - chi2PCA, - pVecD[0], pVecD[1], pVecD[2], - pVecPion[0], pVecPion[1], pVecPion[2], - impactParameter0.getY(), impactParameter1.getY(), - std::sqrt(impactParameter0.getSigmaY2()), std::sqrt(impactParameter1.getSigmaY2()), - candD.globalIndex(), trackPion.globalIndex(), - hfFlag); + collision.posX(), collision.posY(), collision.posZ(), + secondaryVertexB0[0], secondaryVertexB0[1], secondaryVertexB0[2], + errorDecayLength, errorDecayLengthXY, + chi2PCA, + pVecD[0], pVecD[1], pVecD[2], + pVecPion[0], pVecPion[1], pVecPion[2], + dcaD.getY(), dcaPion.getY(), + std::sqrt(dcaD.getSigmaY2()), std::sqrt(dcaPion.getSigmaY2()), + candD.globalIndex(), trackPion.globalIndex(), + hfFlag); } // pi loop } // D loop - } // collision loop - } // process -}; // struct + } // collision loop + } // process +}; // struct /// Extends the base table with expression columns and performs MC matching. struct HfCandidateCreatorB0Expressions { From 7d9af704328d59307073a2e1a8a2692740944f2f Mon Sep 17 00:00:00 2001 From: Alexandre Bigot Date: Wed, 29 Mar 2023 13:43:20 +0200 Subject: [PATCH 3/5] PR comments --- PWGHF/TableProducer/candidateCreatorB0.cxx | 39 ++++++++++------------ 1 file changed, 18 insertions(+), 21 deletions(-) diff --git a/PWGHF/TableProducer/candidateCreatorB0.cxx b/PWGHF/TableProducer/candidateCreatorB0.cxx index 21b92cf940f..5e7b2065c36 100644 --- a/PWGHF/TableProducer/candidateCreatorB0.cxx +++ b/PWGHF/TableProducer/candidateCreatorB0.cxx @@ -77,11 +77,7 @@ struct HfCandidateCreatorB0 { using TracksWithSel = soa::Join; - Filter filterSelectTracks = (!usePionIsGlobalTrackWoDCA) || ((usePionIsGlobalTrackWoDCA) && requireGlobalTrackWoDCAInFilter()); Filter filterSelectCandidates = (aod::hf_sel_candidate_dplus::isSelDplusToPiKPi >= selectionFlagD); - - using TracksWithSelFiltered = soa::Filtered; - using CandsDFiltered = soa::Filtered>; Preslice candsDPerCollision = aod::track_association::collisionId; @@ -118,8 +114,7 @@ struct HfCandidateCreatorB0 { void process(aod::Collisions const& collisions, CandsDFiltered const& candsD, aod::TrackAssoc const& trackIndices, - TracksWithSel const&, - TracksWithSelFiltered const&) + TracksWithSel const&) { // FIXME: store B0 creator configurable (until https://alice.its.cern.ch/jira/browse/O2-3582 solved) if (!isHfCandB0ConfigFilled) { @@ -156,7 +151,7 @@ struct HfCandidateCreatorB0 { auto covMatrixPV = primaryVertex.getCov(); if (ncol % 10000 == 0) { - LOG(info) << ncol << " collisions parsed"; + LOG(debug) << ncol << " collisions parsed"; } ncol++; @@ -176,16 +171,13 @@ struct HfCandidateCreatorB0 { auto trackParCov1 = getTrackParCov(track1); auto trackParCov2 = getTrackParCov(track2); - std::array pVec0; // = {track0.px(), track0.py(), track0.pz()}; - std::array pVec1; // = {track1.px(), track1.py(), track1.pz()}; - std::array pVec2; // = {track2.px(), track2.py(), track2.pz()}; + std::array pVec0 = {track0.px(), track0.py(), track0.pz()}; + std::array pVec1 = {track1.px(), track1.py(), track1.pz()}; + std::array pVec2 = {track2.px(), track2.py(), track2.pz()}; - o2::dataformats::DCA dca0; - // dca0.set(track0.dcaXY(), track0.dcaZ()); - o2::dataformats::DCA dca1; - // dca1.set(track1.dcaXY(), track1.dcaZ()); - o2::dataformats::DCA dca2; - // dca2.set(track2.dcaXY(), track2.dcaZ()); + auto dca0 = o2::dataformats::DCA(track0.dcaXY(), track0.dcaZ()); + auto dca1 = o2::dataformats::DCA(track1.dcaXY(), track1.dcaZ()); + auto dca2 = o2::dataformats::DCA(track2.dcaXY(), track2.dcaZ()); // repropagate tracks to this collision if needed if (track0.collisionId() != thisCollId) { @@ -213,9 +205,9 @@ struct HfCandidateCreatorB0 { trackParCov2.propagateTo(secondaryVertexD[0], bz); // update pVec of tracks - getPxPyPz(trackParCov0, pVec0); // FIXME: or df3.getTrack(0).getPxPyPzGlo(pVec0) ? - getPxPyPz(trackParCov1, pVec1); - getPxPyPz(trackParCov2, pVec2); + df3.getTrack(0).getPxPyPzGlo(pVec0); + df3.getTrack(1).getPxPyPzGlo(pVec1); + df3.getTrack(2).getPxPyPzGlo(pVec2); // D∓ → π∓ K± π∓ array pVecPiK = RecoDecay::pVec(pVec0, pVec1); @@ -232,7 +224,12 @@ struct HfCandidateCreatorB0 { auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); for (const auto& trackId : trackIdsThisCollision) { // start loop over track indices associated to this collision - auto trackPion = trackId.track_as(); + auto trackPion = trackId.track_as(); + + // check isGlobalTrackWoDCA status for pions if wanted + if (usePionIsGlobalTrackWoDCA && !trackPion.isGlobalTrackWoDCA()) { + continue; + } // minimum pT selection if (trackPion.pt() < ptPionMin || !isSelectedTrackDCA(trackPion)) { @@ -293,7 +290,7 @@ struct HfCandidateCreatorB0 { int hfFlag = BIT(hf_cand_b0::DecayType::B0ToDPi); // fill the candidate table for the B0 here: - rowCandidateBase(collision.globalIndex(), + rowCandidateBase(thisCollId, collision.posX(), collision.posY(), collision.posZ(), secondaryVertexB0[0], secondaryVertexB0[1], secondaryVertexB0[2], errorDecayLength, errorDecayLengthXY, From ebff8ed7a5c0b5a1f18dcffadc6ca57ab7b9e725 Mon Sep 17 00:00:00 2001 From: Alexandre Bigot Date: Wed, 29 Mar 2023 14:02:24 +0200 Subject: [PATCH 4/5] Modify initialization of dca for D daughters --- PWGHF/TableProducer/candidateCreatorB0.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/PWGHF/TableProducer/candidateCreatorB0.cxx b/PWGHF/TableProducer/candidateCreatorB0.cxx index 5e7b2065c36..ded4b2e779b 100644 --- a/PWGHF/TableProducer/candidateCreatorB0.cxx +++ b/PWGHF/TableProducer/candidateCreatorB0.cxx @@ -175,9 +175,9 @@ struct HfCandidateCreatorB0 { std::array pVec1 = {track1.px(), track1.py(), track1.pz()}; std::array pVec2 = {track2.px(), track2.py(), track2.pz()}; - auto dca0 = o2::dataformats::DCA(track0.dcaXY(), track0.dcaZ()); - auto dca1 = o2::dataformats::DCA(track1.dcaXY(), track1.dcaZ()); - auto dca2 = o2::dataformats::DCA(track2.dcaXY(), track2.dcaZ()); + auto dca0 = o2::dataformats::DCA(track0.dcaXY(), track0.dcaZ(), track0.cYY(), track0.cZY(), track0.cZZ()); + auto dca1 = o2::dataformats::DCA(track1.dcaXY(), track1.dcaZ(), track1.cYY(), track1.cZY(), track1.cZZ()); + auto dca2 = o2::dataformats::DCA(track2.dcaXY(), track2.dcaZ(), track2.cYY(), track2.cZY(), track2.cZZ()); // repropagate tracks to this collision if needed if (track0.collisionId() != thisCollId) { From 2af990b0d18780233c5871f99caa5f4f8c340a2a Mon Sep 17 00:00:00 2001 From: Alexandre Bigot Date: Wed, 29 Mar 2023 14:21:04 +0200 Subject: [PATCH 5/5] Get bz from CCDB --- PWGHF/TableProducer/candidateCreatorB0.cxx | 48 ++++++++++++++++++++-- 1 file changed, 44 insertions(+), 4 deletions(-) diff --git a/PWGHF/TableProducer/candidateCreatorB0.cxx b/PWGHF/TableProducer/candidateCreatorB0.cxx index ded4b2e779b..153164ecd84 100644 --- a/PWGHF/TableProducer/candidateCreatorB0.cxx +++ b/PWGHF/TableProducer/candidateCreatorB0.cxx @@ -24,6 +24,7 @@ #include "ReconstructionDataFormats/V0.h" #include "PWGHF/DataModel/CandidateReconstructionTables.h" #include "PWGHF/DataModel/CandidateSelectionTables.h" +#include "PWGHF/Utils/utilsBfieldCCDB.h" using namespace o2; using namespace o2::aod; @@ -52,7 +53,7 @@ struct HfCandidateCreatorB0 { Produces rowCandidateConfig; // vertexing - Configurable bz{"bz", 5., "magnetic field"}; + // Configurable bz{"bz", 5., "magnetic field"}; Configurable propagateToPCA{"propagateToPCA", true, "create tracks version propagated to PCA"}; Configurable useAbsDCA{"useAbsDCA", false, "Minimise abs. distance rather than chi2"}; Configurable useWeightedFinalPCA{"useWeightedFinalPCA", false, "Recalculate vertex position using track covariances, effective only if useAbsDCA is true"}; @@ -69,11 +70,24 @@ struct HfCandidateCreatorB0 { Configurable selectionFlagD{"selectionFlagD", 1, "Selection Flag for D"}; // FIXME: store B0 creator configurable (until https://alice.its.cern.ch/jira/browse/O2-3582 solved) bool isHfCandB0ConfigFilled = false; + // magnetic field setting from CCDB + Configurable isRun2{"isRun2", false, "enable Run 2 or Run 3 GRP objects for magnetic field"}; + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable ccdbPathLut{"ccdbPathLut", "GLO/Param/MatLUT", "Path for LUT parametrization"}; + Configurable ccdbPathGeo{"ccdbPathGeo", "GLO/Config/GeometryAligned", "Path of the geometry file"}; + Configurable ccdbPathGrp{"ccdbPathGrp", "GLO/GRP/GRP", "Path of the grp file (Run 2)"}; + Configurable ccdbPathGrpMag{"ccdbPathGrpMag", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object (Run 3)"}; + + Service ccdb; + o2::base::MatLayerCylSet* lut; + o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrLUT; + int runNumber; double massPi = RecoDecay::getMassPDG(kPiPlus); double massD = RecoDecay::getMassPDG(pdg::Code::kDMinus); double massB0 = RecoDecay::getMassPDG(pdg::Code::kB0); double massDPi{0.}; + double bz{0.}; using TracksWithSel = soa::Join; @@ -91,6 +105,18 @@ struct HfCandidateCreatorB0 { OutputObj hCovPVXX{TH1F("hCovPVXX", "2-prong candidates;XX element of cov. matrix of prim. vtx. position (cm^{2});entries", 100, 0., 1.e-4)}; OutputObj hCovSVXX{TH1F("hCovSVXX", "2-prong candidates;XX element of cov. matrix of sec. vtx. position (cm^{2});entries", 100, 0., 0.2)}; + void init(InitContext const&) + { + ccdb->setURL(ccdbUrl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->get(ccdbPathLut)); + if (!o2::base::GeometryManager::isGeometryLoaded()) { + ccdb->get(ccdbPathGeo); + } + runNumber = 0; + } + /// Single-track cuts for pions on dcaXY /// \param track is a track /// \return true if track passes all cuts @@ -114,7 +140,8 @@ struct HfCandidateCreatorB0 { void process(aod::Collisions const& collisions, CandsDFiltered const& candsD, aod::TrackAssoc const& trackIndices, - TracksWithSel const&) + TracksWithSel const&, + aod::BCsWithTimestamps const&) { // FIXME: store B0 creator configurable (until https://alice.its.cern.ch/jira/browse/O2-3582 solved) if (!isHfCandB0ConfigFilled) { @@ -124,7 +151,7 @@ struct HfCandidateCreatorB0 { } // Initialise fitter for B vertex (2-prong vertex filter) o2::vertexing::DCAFitterN<2> df2; - df2.setBz(bz); + // df2.setBz(bz); df2.setPropagateToPCA(propagateToPCA); df2.setMaxR(maxR); df2.setMaxDZIni(maxDZIni); @@ -135,7 +162,7 @@ struct HfCandidateCreatorB0 { // Initial fitter to redo D-vertex to get extrapolated daughter tracks (3-prong vertex filter) o2::vertexing::DCAFitterN<3> df3; - df3.setBz(bz); + // df3.setBz(bz); df3.setPropagateToPCA(propagateToPCA); df3.setMaxR(maxR); df3.setMaxDZIni(maxDZIni); @@ -155,6 +182,19 @@ struct HfCandidateCreatorB0 { } ncol++; + /// Set the magnetic field from ccdb. + /// The static instance of the propagator was already modified in the HFTrackIndexSkimCreator, + /// but this is not true when running on Run2 data/MC already converted into AO2Ds. + auto bc = collision.bc_as(); + if (runNumber != bc.runNumber()) { + LOG(info) << ">>>>>>>>>>>> Current run number: " << runNumber; + initCCDB(bc, runNumber, ccdb, isRun2 ? ccdbPathGrp : ccdbPathGrpMag, lut, isRun2); + bz = o2::base::Propagator::Instance()->getNominalBz(); + LOG(info) << ">>>>>>>>>>>> Magnetic field: " << bz; + } + df2.setBz(bz); + df3.setBz(bz); + auto thisCollId = collision.globalIndex(); auto candsDThisColl = candsD.sliceBy(candsDPerCollision, thisCollId);