From 0a7a98c7530c6069632329a5ade54e0084aacb1d Mon Sep 17 00:00:00 2001 From: Daiki Sekihata Date: Fri, 26 Jun 2026 22:08:30 +0200 Subject: [PATCH] [PWGEM/Dilepton] update studyDCAFitter.cxx --- PWGEM/Dilepton/DataModel/lmeeMLTables.h | 48 +++-- PWGEM/Dilepton/Tasks/CMakeLists.txt | 2 +- PWGEM/Dilepton/Tasks/studyDCAFitter.cxx | 247 ++++++++++++++++++++++-- 3 files changed, 267 insertions(+), 30 deletions(-) diff --git a/PWGEM/Dilepton/DataModel/lmeeMLTables.h b/PWGEM/Dilepton/DataModel/lmeeMLTables.h index 0f85851075b..6ab6a8d70fe 100644 --- a/PWGEM/Dilepton/DataModel/lmeeMLTables.h +++ b/PWGEM/Dilepton/DataModel/lmeeMLTables.h @@ -346,25 +346,39 @@ using EMMLLCascPair = EMMLLCascPairs::iterator; namespace emmldilepton { -DECLARE_SOA_INDEX_COLUMN(EMMLEvent, emmlevent); //! index to event table -DECLARE_SOA_COLUMN(Signed1Pt1, signed1Pt1, float); //! q/pt of lepton1 at PV -DECLARE_SOA_COLUMN(Eta1, eta1, float); //! eta of lepton1 at PV -DECLARE_SOA_COLUMN(ImpParXY1, impParXY1, float); //! impact parameter for lepton1 in XY plane -DECLARE_SOA_COLUMN(ImpParZ1, impParZ1, float); //! impact parameter for lepton1 in Z plane -DECLARE_SOA_COLUMN(ImpParCYY1, impParCYY1, float); //! sigma of impact parameter for lepton1 in XY -DECLARE_SOA_COLUMN(ImpParCZY1, impParCZY1, float); //! sigma of impact parameter for lepton1, correlaion term -DECLARE_SOA_COLUMN(ImpParCZZ1, impParCZZ1, float); //! sigma of impact parameter for lepton1 in Z +DECLARE_SOA_INDEX_COLUMN(EMMLEvent, emmlevent); //! index to event table +DECLARE_SOA_COLUMN(Signed1Pt1, signed1Pt1, float); //! q/pt of lepton1 at PV +DECLARE_SOA_COLUMN(Eta1, eta1, float); //! eta of lepton1 at PV +DECLARE_SOA_COLUMN(ImpParXY1, impParXY1, float); //! impact parameter for lepton1 in XY plane +DECLARE_SOA_COLUMN(ImpParZ1, impParZ1, float); //! impact parameter for lepton1 in Z plane +DECLARE_SOA_COLUMN(ImpParCYY1, impParCYY1, float); //! sigma of impact parameter for lepton1 in XY +DECLARE_SOA_COLUMN(ImpParCZY1, impParCZY1, float); //! sigma of impact parameter for lepton1, correlaion term +DECLARE_SOA_COLUMN(ImpParCZZ1, impParCZZ1, float); //! sigma of impact parameter for lepton1 in Z + +DECLARE_SOA_COLUMN(UnbiasedImpParXY1, unbiasedImpParXY1, float); //! impact parameter for lepton1 in XY plane +DECLARE_SOA_COLUMN(UnbiasedImpParZ1, unbiasedImpParZ1, float); //! impact parameter for lepton1 in Z plane +DECLARE_SOA_COLUMN(UnbiasedImpParCYY1, unbiasedImpParCYY1, float); //! sigma of impact parameter for lepton1 in XY +DECLARE_SOA_COLUMN(UnbiasedImpParCZY1, unbiasedImpParCZY1, float); //! sigma of impact parameter for lepton1, correlaion term +DECLARE_SOA_COLUMN(UnbiasedImpParCZZ1, unbiasedImpParCZZ1, float); //! sigma of impact parameter for lepton1 in Z + DECLARE_SOA_COLUMN(IsCorrectCollision1, isCorrectCollision1, bool); //! lepton1 is associated to correct collision. DECLARE_SOA_COLUMN(IsReassociated1, isReassociated1, bool); //! lepton1 is reassociated. DECLARE_SOA_COLUMN(PdgCodeMother1, pdgCodeMother1, int); //! pdg code of mother of lepton1 -DECLARE_SOA_COLUMN(Signed1Pt2, signed1Pt2, float); //! q/pt of lepton2 at PV -DECLARE_SOA_COLUMN(Eta2, eta2, float); //! eta of lepton1 at PV -DECLARE_SOA_COLUMN(ImpParXY2, impParXY2, float); //! impact parameter for lepton2 in XY plane -DECLARE_SOA_COLUMN(ImpParZ2, impParZ2, float); //! impact parameter for lepton2 in Z plane -DECLARE_SOA_COLUMN(ImpParCYY2, impParCYY2, float); //! sigma of impact parameter for lepton2 in XY -DECLARE_SOA_COLUMN(ImpParCZY2, impParCZY2, float); //! sigma of impact parameter for lepton2, correlaion term -DECLARE_SOA_COLUMN(ImpParCZZ2, impParCZZ2, float); //! sigma of impact parameter for lepton2 in Z +DECLARE_SOA_COLUMN(Signed1Pt2, signed1Pt2, float); //! q/pt of lepton2 at PV +DECLARE_SOA_COLUMN(Eta2, eta2, float); //! eta of lepton1 at PV +DECLARE_SOA_COLUMN(ImpParXY2, impParXY2, float); //! impact parameter for lepton2 in XY plane +DECLARE_SOA_COLUMN(ImpParZ2, impParZ2, float); //! impact parameter for lepton2 in Z plane +DECLARE_SOA_COLUMN(ImpParCYY2, impParCYY2, float); //! sigma of impact parameter for lepton2 in XY +DECLARE_SOA_COLUMN(ImpParCZY2, impParCZY2, float); //! sigma of impact parameter for lepton2, correlaion term +DECLARE_SOA_COLUMN(ImpParCZZ2, impParCZZ2, float); //! sigma of impact parameter for lepton2 in Z + +DECLARE_SOA_COLUMN(UnbiasedImpParXY2, unbiasedImpParXY2, float); //! impact parameter for lepton2 in XY plane +DECLARE_SOA_COLUMN(UnbiasedImpParZ2, unbiasedImpParZ2, float); //! impact parameter for lepton2 in Z plane +DECLARE_SOA_COLUMN(UnbiasedImpParCYY2, unbiasedImpParCYY2, float); //! sigma of impact parameter for lepton2 in XY +DECLARE_SOA_COLUMN(UnbiasedImpParCZY2, unbiasedImpParCZY2, float); //! sigma of impact parameter for lepton2, correlaion term +DECLARE_SOA_COLUMN(UnbiasedImpParCZZ2, unbiasedImpParCZZ2, float); //! sigma of impact parameter for lepton2 in Z + DECLARE_SOA_COLUMN(IsCorrectCollision2, isCorrectCollision2, bool); //! lepton is associated to correct collision. DECLARE_SOA_COLUMN(IsReassociated2, isReassociated2, bool); //! lepton2 is reassociated. DECLARE_SOA_COLUMN(PdgCodeMother2, pdgCodeMother2, int); //! pdg code of mother of lepton1 @@ -392,6 +406,10 @@ DECLARE_SOA_TABLE(EMMLDielectronsAtSV, "AOD", "EMMLEESV", //! emmldilepton::EMMLEventId, emmldilepton::Signed1Pt1, emmldilepton::Eta1, emmldilepton::ImpParXY1, emmldilepton::ImpParZ1, emmldilepton::ImpParCYY1, emmldilepton::ImpParCZY1, emmldilepton::ImpParCZZ1, emmldilepton::IsCorrectCollision1, emmldilepton::IsReassociated1, emmldilepton::PdgCodeMother1, emmldilepton::Signed1Pt2, emmldilepton::Eta2, emmldilepton::ImpParXY2, emmldilepton::ImpParZ2, emmldilepton::ImpParCYY2, emmldilepton::ImpParCZY2, emmldilepton::ImpParCZZ2, emmldilepton::IsCorrectCollision2, emmldilepton::IsReassociated2, emmldilepton::PdgCodeMother2, + + emmldilepton::UnbiasedImpParXY1, emmldilepton::UnbiasedImpParZ1, emmldilepton::UnbiasedImpParCYY1, emmldilepton::UnbiasedImpParCZY1, emmldilepton::UnbiasedImpParCZZ1, + emmldilepton::UnbiasedImpParXY2, emmldilepton::UnbiasedImpParZ2, emmldilepton::UnbiasedImpParCYY2, emmldilepton::UnbiasedImpParCZY2, emmldilepton::UnbiasedImpParCZZ2, + emmldilepton::Mass, emmldilepton::Pt, emmldilepton::Rapidity, emmldilepton::Chi2PCA, emmldilepton::CPA, emmldilepton::CPAXY, emmldilepton::CPARZ, diff --git a/PWGEM/Dilepton/Tasks/CMakeLists.txt b/PWGEM/Dilepton/Tasks/CMakeLists.txt index 8ac8b8fba4a..bb0adb4d478 100644 --- a/PWGEM/Dilepton/Tasks/CMakeLists.txt +++ b/PWGEM/Dilepton/Tasks/CMakeLists.txt @@ -152,7 +152,7 @@ o2physics_add_dpl_workflow(mc-particle-predictions-otf o2physics_add_dpl_workflow(study-dcafitter SOURCES studyDCAFitter.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::DCAFitter O2Physics::AnalysisCore + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsVertexing O2::DCAFitter O2Physics::AnalysisCore COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(evaluate-acceptance diff --git a/PWGEM/Dilepton/Tasks/studyDCAFitter.cxx b/PWGEM/Dilepton/Tasks/studyDCAFitter.cxx index 832f3685bf2..40a5946c206 100644 --- a/PWGEM/Dilepton/Tasks/studyDCAFitter.cxx +++ b/PWGEM/Dilepton/Tasks/studyDCAFitter.cxx @@ -35,6 +35,7 @@ #include #include #include +#include // for PV refit #include #include #include @@ -45,6 +46,7 @@ #include #include #include +#include // for PV refit #include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) #include @@ -56,6 +58,7 @@ #include #include #include +#include #include #include @@ -70,8 +73,19 @@ using namespace o2::aod::pwgem::dilepton::utils::pairutil; struct studyDCAFitter { using MyCollisions = soa::Join; using MyTracks = soa::Join; + using MyTrack = MyTracks::iterator; using MyBCs = soa::Join; + // struct unbiasedDCA { + // int globalIndex{0}; + // float dcaXY{999.f}; + // float dcaZ{999.f}; + // float cYY{999.f}; + // float cZZ{999.f}; + // float cZY{999.f}; + // bool isOK{false}; + // }; + Produces eventTable; Produces dileptonTable; @@ -233,6 +247,11 @@ struct studyDCAFitter { fRegistry.add("Event/hCentFT0M", "hCentFT0M;centrality FT0M (%)", kTH1F, {{110, 0, 110}}, false); fRegistry.add("Event/hCentFT0CvsMultNTracksPV", "hCentFT0CvsMultNTracksPV;centrality FT0C (%);N_{track} to PV", kTH2F, {{110, 0, 110}, {600, 0, 6000}}, false); fRegistry.add("Event/hMultFT0CvsMultNTracksPV", "hMultFT0CvsMultNTracksPV;mult. FT0C;N_{track} to PV", kTH2F, {{60, 0, 60000}, {600, 0, 6000}}, false); + fRegistry.add("Event/refitPV/hNContrib", "hNContrib;before;after", kTH2F, {{1001, -0.5, 1000.5}, {1001, -0.5, 1000.5}}, false); + fRegistry.add("Event/refitPV/hChi2", "hChi2;before;after", kTH2F, {{100, 0, 1000}, {100, 0, 1000}}, false); + fRegistry.add("Event/refitPV/hDeltaXvsNContrib", "hDeltaXvsNContrib;numContrib after;X_{after} - X_{before} (cm)", kTH2F, {{1001, -0.5, 1000.5}, {200, -0.01, 0.01}}, false); + fRegistry.add("Event/refitPV/hDeltaYvsNContrib", "hDeltaYvsNContrib;numContrib after;Y_{after} - Y_{before} (cm)", kTH2F, {{1001, -0.5, 1000.5}, {200, -0.01, 0.01}}, false); + fRegistry.add("Event/refitPV/hDeltaZvsNContrib", "hDeltaZvsNContrib;numContrib after;Z_{after} - Z_{before} (cm)", kTH2F, {{1001, -0.5, 1000.5}, {200, -0.01, 0.01}}, false); const o2::framework::AxisSpec axis_mass{ConfMllBins, "m_{ll} (GeV/c^{2})"}; const o2::framework::AxisSpec axis_pt{ConfPtllBins, "p_{T,ee} (GeV/c)"}; @@ -510,13 +529,156 @@ struct studyDCAFitter { } // end of HFee } - template - void runSVFinder(TCollision const& collision, TTrack const& t1, TTrack const& t2, TMCParticles const& mcParticles) + template + o2::dataformats::VertexBase refitPV(TCollision const& collision, std::vector const& vecPvContributorTrackParCov, std::vector const& vecPvContributorGlobalId, std::vector const& tracksToRemove) + { + std::vector vecPvRefitContributorUsed(vecPvContributorGlobalId.size(), true); + + // build the VertexBase to initialize the vertexer + o2::dataformats::VertexBase primVtx; + primVtx.setX(collision.posX()); + primVtx.setY(collision.posY()); + primVtx.setZ(collision.posZ()); + primVtx.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()); + + o2::vertexing::PVertexer vertexer; + o2::conf::ConfigurableParam::updateFromString("pvertexer.useMeanVertexConstraint=false"); // remove diamond constraint + vertexer.init(); + const bool pvRefitDoable = vertexer.prepareVertexRefit(vecPvContributorTrackParCov, primVtx); + if (!pvRefitDoable) { + LOG(info) << "Not enough tracks accepted for the refit"; // this should not happen by definition, because nPV>=2 is required in the reconstruction. I analyze the reconstructed collisions in AO2D. + return primVtx; + } + + for (const auto& trackToRemove : tracksToRemove) { + const auto trackIterator = std::find(vecPvContributorGlobalId.begin(), vecPvContributorGlobalId.end(), trackToRemove.globalIndex()); // track global index + if (trackIterator != vecPvContributorGlobalId.end()) { + const int entry = std::distance(vecPvContributorGlobalId.begin(), trackIterator); // this track contributed to this collision: let's do the refit without it + vecPvRefitContributorUsed[entry] = false; // remove the track from the PV refitting + } + } + + const auto primVtxRefitted = vertexer.refitVertexFull(vecPvRefitContributorUsed, primVtx); // vertex refit + // LOG(info) << "refit for track with global index " << static_cast(trackToRemove.globalIndex()) << " " << primVtxRefitted.asString(); + if (primVtxRefitted.getChi2() < 0) { + LOG(info) << "---> Refitted vertex has bad chi2 = " << primVtxRefitted.getChi2(); + return primVtx; + } + + fRegistry.fill(HIST("Event/refitPV/hNContrib"), collision.numContrib(), primVtxRefitted.getNContributors()); + fRegistry.fill(HIST("Event/refitPV/hChi2"), collision.chi2(), primVtxRefitted.getChi2()); + fRegistry.fill(HIST("Event/refitPV/hDeltaXvsNContrib"), primVtxRefitted.getNContributors(), primVtxRefitted.getX() - primVtx.getX()); + fRegistry.fill(HIST("Event/refitPV/hDeltaYvsNContrib"), primVtxRefitted.getNContributors(), primVtxRefitted.getY() - primVtx.getY()); + fRegistry.fill(HIST("Event/refitPV/hDeltaZvsNContrib"), primVtxRefitted.getNContributors(), primVtxRefitted.getZ() - primVtx.getZ()); + + o2::dataformats::VertexBase primVtxBaseRecalc; + primVtxBaseRecalc.setX(primVtxRefitted.getX()); + primVtxBaseRecalc.setY(primVtxRefitted.getY()); + primVtxBaseRecalc.setZ(primVtxRefitted.getZ()); + primVtxBaseRecalc.setCov(primVtxRefitted.getSigmaX2(), primVtxRefitted.getSigmaXY(), primVtxRefitted.getSigmaY2(), primVtxRefitted.getSigmaXZ(), primVtxRefitted.getSigmaYZ(), primVtxRefitted.getSigmaZ2()); + + vecPvRefitContributorUsed.clear(); + vecPvRefitContributorUsed.shrink_to_fit(); + // return primVtxRefitted; + return primVtxBaseRecalc; + } + + // template + // bool refitPV(TCollision const& collision, std::vector const& vecPvContributorTrackParCov, std::vector const& vecPvContributorGlobalId, std::vector const& tracksToRemove, unbiasedDCA& candidate) + // { + // std::vector vecPvRefitContributorUsed(vecPvContributorGlobalId.size(), true); + + // // build the VertexBase to initialize the vertexer + // o2::dataformats::VertexBase primVtx; + // primVtx.setX(collision.posX()); + // primVtx.setY(collision.posY()); + // primVtx.setZ(collision.posZ()); + // primVtx.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()); + + // o2::vertexing::PVertexer vertexer; + // o2::conf::ConfigurableParam::updateFromString("pvertexer.useMeanVertexConstraint=false"); // remove diamond constraint + // vertexer.init(); + // const bool pvRefitDoable = vertexer.prepareVertexRefit(vecPvContributorTrackParCov, primVtx); + // if (!pvRefitDoable) { + // LOG(info) << "Not enough tracks accepted for the refit"; // this should not happen by definition, because nPV>=2 is required in the reconstruction. I analyze the reconstructed collisions in AO2D. + // candidate.isOK = false; + // return false; + // } + + // // PV refitting, if the tracks contributed to this at the beginning + // o2::dataformats::VertexBase primVtxBaseRecalc; + // bool recalcImpPar = false; + // if (pvRefitDoable) { + // recalcImpPar = true; + + // for (const auto& trackToRemove : tracksToRemove) { + // const auto trackIterator = std::find(vecPvContributorGlobalId.begin(), vecPvContributorGlobalId.end(), trackToRemove.globalIndex()); /// track global index + // if (trackIterator != vecPvContributorGlobalId.end()) { + // const int entry = std::distance(vecPvContributorGlobalId.begin(), trackIterator); // this track contributed to this collision: let's do the refit without it + // vecPvRefitContributorUsed[entry] = false; // remove the track from the PV refitting + // } + // } + + // const auto primVtxRefitted = vertexer.refitVertexFull(vecPvRefitContributorUsed, primVtx); // vertex refit + + // // LOG(info) << "refit for track with global index " << static_cast(trackToRemove.globalIndex()) << " " << primVtxRefitted.asString(); + // if (primVtxRefitted.getChi2() < 0) { + // LOG(info) << "---> Refitted vertex has bad chi2 = " << primVtxRefitted.getChi2(); + // candidate.isOK = false; + // recalcImpPar = false; + // return false; + // } + + // // LOGF(info, "collision.numContrib() = %d, primVtxRefitted.getNContributors() = %d", collision.numContrib(), primVtxRefitted.getNContributors()); + // fRegistry.fill(HIST("Event/refitPV/hNContrib"), collision.numContrib(), primVtxRefitted.getNContributors()); + // fRegistry.fill(HIST("Event/refitPV/hChi2"), collision.chi2(), primVtxRefitted.getChi2()); + // fRegistry.fill(HIST("Event/refitPV/hDeltaXvsNContrib"), primVtxRefitted.getNContributors(), primVtxRefitted.getX() - primVtx.getX()); + // fRegistry.fill(HIST("Event/refitPV/hDeltaYvsNContrib"), primVtxRefitted.getNContributors(), primVtxRefitted.getY() - primVtx.getY()); + // fRegistry.fill(HIST("Event/refitPV/hDeltaZvsNContrib"), primVtxRefitted.getNContributors(), primVtxRefitted.getZ() - primVtx.getZ()); + + // // vecPvRefitContributorUsed[entry] = true; // restore the track for the next PV refitting (probably not necessary here) + + // if (recalcImpPar) { + // // fill the newly calculated PV + // primVtxBaseRecalc.setX(primVtxRefitted.getX()); + // primVtxBaseRecalc.setY(primVtxRefitted.getY()); + // primVtxBaseRecalc.setZ(primVtxRefitted.getZ()); + // primVtxBaseRecalc.setCov(primVtxRefitted.getSigmaX2(), primVtxRefitted.getSigmaXY(), primVtxRefitted.getSigmaY2(), primVtxRefitted.getSigmaXZ(), primVtxRefitted.getSigmaYZ(), primVtxRefitted.getSigmaZ2()); + // } + + // // updated value after PV recalculation + // if (recalcImpPar) { + // for (const auto& trackToRemove : tracksToRemove) { + // auto trackParCov = getTrackParCov(trackToRemove); + // o2::dataformats::DCA impactParameter; + // impactParameter.set(999, 999, 999, 999, 999); + // trackParCov.setPID(o2::track::PID::Electron); + + // // if (o2::base::Propagator::Instance()->propagateToDCABxByBz(primVtxBaseRecalc, trackParCov, 2.f, matCorr, &impactParameter)) { + // if (o2::base::Propagator::Instance()->propagateToDCABxByBz(primVtxRefitted, trackParCov, 2.f, matCorr, &impactParameter)) { + // candidate.dcaXY = impactParameter.getY(); + // candidate.dcaZ = impactParameter.getZ(); + // candidate.cYY = trackParCov.getSigmaY2(); + // candidate.cZZ = trackParCov.getSigmaZ2(); + // candidate.cZY = trackParCov.getSigmaZY(); + // // LOGF(info, "trackToRemove.globalIndex() = %d, candidate.dcaXY = %f, candidate.dcaZ = %f, candidate.cYY = %.16f, candidate.cZZ = %.16f, candidate.cZY = %.16f", trackToRemove.globalIndex(), candidate.dcaXY, candidate.dcaZ, candidate.cYY, candidate.cZZ, candidate.cZY); + // candidate.isOK = true; + // } + // } + // } + // } /// end 'if (doPvRefit && pvRefitDoable)' + + // vecPvRefitContributorUsed.clear(); + // vecPvRefitContributorUsed.shrink_to_fit(); + // return true; + // } + + template + void runSVFinder(TCollision const& collision, TTrack const& t1, TTrack const& t2, TMCParticles const& mcParticles, TRefittedPV const& refittedPV) { auto trackParCov1 = getTrackParCov(t1); trackParCov1.setPID(o2::track::PID::Electron); mDcaInfoCov.set(999, 999, 999, 999, 999); - trackParCov1.setPID(o2::track::PID::Electron); o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, trackParCov1, 2.f, matCorr, &mDcaInfoCov); float dcaXY1 = mDcaInfoCov.getY(); float dcaZ1 = mDcaInfoCov.getZ(); @@ -526,10 +688,17 @@ struct studyDCAFitter { float signed1Pt1 = trackParCov1.getQ2Pt(); // at PV float eta1 = trackParCov1.getEta(); // at PV + mDcaInfoCov.set(999, 999, 999, 999, 999); + o2::base::Propagator::Instance()->propagateToDCABxByBz(refittedPV, trackParCov1, 2.f, matCorr, &mDcaInfoCov); + float unbiased_dcaXY1 = mDcaInfoCov.getY(); + float unbiased_dcaZ1 = mDcaInfoCov.getZ(); + float unbiased_CYY1 = trackParCov1.getSigmaY2(); + float unbiased_CZY1 = trackParCov1.getSigmaZY(); + float unbiased_CZZ1 = trackParCov1.getSigmaZ2(); + auto trackParCov2 = getTrackParCov(t2); trackParCov2.setPID(o2::track::PID::Electron); mDcaInfoCov.set(999, 999, 999, 999, 999); - trackParCov2.setPID(o2::track::PID::Electron); o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, trackParCov2, 2.f, matCorr, &mDcaInfoCov); float dcaXY2 = mDcaInfoCov.getY(); float dcaZ2 = mDcaInfoCov.getZ(); @@ -539,7 +708,16 @@ struct studyDCAFitter { float signed1Pt2 = trackParCov2.getQ2Pt(); // at PV float eta2 = trackParCov2.getEta(); // at PV - std::array pVtx = {collision.posX(), collision.posY(), collision.posZ()}; + mDcaInfoCov.set(999, 999, 999, 999, 999); + o2::base::Propagator::Instance()->propagateToDCABxByBz(refittedPV, trackParCov2, 2.f, matCorr, &mDcaInfoCov); + float unbiased_dcaXY2 = mDcaInfoCov.getY(); + float unbiased_dcaZ2 = mDcaInfoCov.getZ(); + float unbiased_CYY2 = trackParCov2.getSigmaY2(); + float unbiased_CZY2 = trackParCov2.getSigmaZY(); + float unbiased_CZZ2 = trackParCov2.getSigmaZ2(); + + // std::array pVtx = {collision.posX(), collision.posY(), collision.posZ()}; + std::array pVtx = {refittedPV.getX(), refittedPV.getY(), refittedPV.getZ()}; std::array svpos = {0.}; // secondary vertex position std::array pvec0 = {0.}; std::array pvec1 = {0.}; @@ -575,11 +753,12 @@ struct studyDCAFitter { return; } - float lxy = std::sqrt(std::pow(svpos[0] - collision.posX(), 2) + std::pow(svpos[1] - collision.posY(), 2)); // in cm - float lz = svpos[2] - collision.posZ(); // in cm + float lxy = std::sqrt(std::pow(svpos[0] - refittedPV.getX(), 2) + std::pow(svpos[1] - refittedPV.getY(), 2)); // in cm + float lz = svpos[2] - refittedPV.getZ(); // in cm float lxyz = std::sqrt(std::pow(lxy, 2) + std::pow(lz, 2)); - auto primaryVertex = getPrimaryVertex(collision); + // auto primaryVertex = getPrimaryVertex(collision); + auto primaryVertex = refittedPV; std::array covVtx = fitter.calcPCACovMatrixFlat(); double phi{}, theta{}; getPointDirection(std::array{primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, svpos, phi, theta); @@ -772,6 +951,8 @@ struct studyDCAFitter { dileptonTable(eventTable.lastIndex() + 1, signed1Pt1, eta1, dcaXY1, dcaZ1, CYY1, CZY1, CZZ1, isCorrectCollision1, isReassociated1, pdgCodeMother1, signed1Pt2, eta2, dcaXY2, dcaZ2, CYY2, CZY2, CZZ2, isCorrectCollision2, isReassociated2, pdgCodeMother2, + unbiased_dcaXY1, unbiased_dcaZ1, unbiased_CYY1, unbiased_CZY1, unbiased_CZZ1, + unbiased_dcaXY2, unbiased_dcaZ2, unbiased_CYY2, unbiased_CZY2, unbiased_CZZ2, meeAtSV, pteeAtSV, yeeAtSV, chi2PCA, cpa, cpaXY, cpaRZ, @@ -801,6 +982,8 @@ struct studyDCAFitter { fillEventHistograms(collision); auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, collision.globalIndex()); + std::vector electronIds; + std::vector positronIds; electronIds.reserve(trackIdsThisCollision.size()); positronIds.reserve(trackIdsThisCollision.size()); mVtx.setPos({collision.posX(), collision.posY(), collision.posZ()}); @@ -851,11 +1034,38 @@ struct studyDCAFitter { } } // end of track loop for electron selection + auto pvContributors_per_collision = pvContributors->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache); + std::vector vecPvContributorTrackParCov; + std::vector vecPvContributorGlobalId; + vecPvContributorTrackParCov.reserve(pvContributors_per_collision.size()); + vecPvContributorGlobalId.reserve(pvContributors_per_collision.size()); + + for (const auto& track : pvContributors_per_collision) { + vecPvContributorGlobalId.push_back(track.globalIndex()); + vecPvContributorTrackParCov.push_back(getTrackParCov(track)); + } + // LOGF(info, "collision.numContrib() = %d, pvContributors_per_collision.size() = %d, vecPvContributorGlobalId.size() = %d", collision.numContrib(), pvContributors_per_collision.size(), vecPvContributorGlobalId.size()); // these 3 numbers must be consistent. + + // std::unordered_map map_unbiasedDCA; // track.globalIndex() -> unbiasedDCA + // for (const auto& posId : positronIds) { + // unbiasedDCA candidate; + // auto pos = tracks.rawIteratorAt(posId); + // refitPV(collision, vecPvContributorTrackParCov, vecPvContributorGlobalId, std::vector{pos}, candidate); + // map_unbiasedDCA[pos.globalIndex()] = candidate; + // } + // for (const auto& eleId : electronIds) { + // unbiasedDCA candidate; + // auto ele = tracks.rawIteratorAt(eleId); + // refitPV(collision, vecPvContributorTrackParCov, vecPvContributorGlobalId, std::vector{ele}, candidate); + // map_unbiasedDCA[ele.globalIndex()] = candidate; + // } + for (const auto& posId : positronIds) { auto pos = tracks.rawIteratorAt(posId); for (const auto& eleId : electronIds) { auto ele = tracks.rawIteratorAt(eleId); - runSVFinder<0>(collision, pos, ele, mcParticles); + const auto& refittedPV = refitPV(collision, vecPvContributorTrackParCov, vecPvContributorGlobalId, std::vector{pos, ele}); + runSVFinder<0>(collision, pos, ele, mcParticles, refittedPV); runPairingAtPV<0>(pos, ele, mcParticles); } // end of electron loop } // end of positron loop @@ -864,7 +1074,8 @@ struct studyDCAFitter { auto pos1 = tracks.rawIteratorAt(positronIds[i]); for (size_t j = i + 1; j < positronIds.size(); j++) { auto pos2 = tracks.rawIteratorAt(positronIds[j]); - runSVFinder<1>(collision, pos1, pos2, mcParticles); + const auto& refittedPV = refitPV(collision, vecPvContributorTrackParCov, vecPvContributorGlobalId, std::vector{pos1, pos2}); + runSVFinder<1>(collision, pos1, pos2, mcParticles, refittedPV); runPairingAtPV<1>(pos1, pos2, mcParticles); } // end of positron loop } // end of positron loop @@ -873,7 +1084,8 @@ struct studyDCAFitter { auto ele1 = tracks.rawIteratorAt(electronIds[i]); for (size_t j = i + 1; j < electronIds.size(); j++) { auto ele2 = tracks.rawIteratorAt(electronIds[j]); - runSVFinder<2>(collision, ele1, ele2, mcParticles); + const auto& refittedPV = refitPV(collision, vecPvContributorTrackParCov, vecPvContributorGlobalId, std::vector{ele1, ele2}); + runSVFinder<2>(collision, ele1, ele2, mcParticles, refittedPV); runPairingAtPV<2>(ele1, ele2, mcParticles); } // end of electron loop } // end of electron loop @@ -887,6 +1099,14 @@ struct studyDCAFitter { positronIds.clear(); positronIds.shrink_to_fit(); + vecPvContributorGlobalId.clear(); + vecPvContributorTrackParCov.clear(); + + vecPvContributorGlobalId.shrink_to_fit(); + vecPvContributorTrackParCov.shrink_to_fit(); + + // map_unbiasedDCA.clear(); + } // end of collision loop } @@ -897,13 +1117,12 @@ struct studyDCAFitter { Filter collisionFilter_centrality = (eventCut.cfgCentMin < o2::aod::cent::centFT0M && o2::aod::cent::centFT0M < eventCut.cfgCentMax) || (eventCut.cfgCentMin < o2::aod::cent::centFT0A && o2::aod::cent::centFT0A < eventCut.cfgCentMax) || (eventCut.cfgCentMin < o2::aod::cent::centFT0C && o2::aod::cent::centFT0C < eventCut.cfgCentMax); using FilteredMyCollisions = soa::Filtered; + Partition pvContributors = ((o2::aod::track::flags & static_cast(o2::aod::track::PVContributor)) == static_cast(o2::aod::track::PVContributor)); + Preslice trackIndicesPerCollision = aod::track_association::collisionId; // Partition posTracks = o2::aod::track::signed1Pt > 0.f; // Partition negTracks = o2::aod::track::signed1Pt < 0.f; - std::vector electronIds; - std::vector positronIds; - void processMC(FilteredMyCollisions const& collisions, MyBCs const& bcs, MyTracks const& tracks, aod::TrackAssoc const& trackIndices, aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles) { run(bcs, collisions, tracks, trackIndices, mcCollisions, mcParticles);