diff --git a/PWGMM/Mult/Tasks/dndetaMFTPbPb.cxx b/PWGMM/Mult/Tasks/dndetaMFTPbPb.cxx index e7c6b9c5b5d..4380a7471c6 100644 --- a/PWGMM/Mult/Tasks/dndetaMFTPbPb.cxx +++ b/PWGMM/Mult/Tasks/dndetaMFTPbPb.cxx @@ -20,6 +20,7 @@ #include "bestCollisionTable.h" #include "Common/CCDB/ctpRateFetcher.h" +#include "Common/Core/fwdtrackUtilities.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/CollisionAssociationTables.h" #include "Common/DataModel/EventSelection.h" @@ -72,10 +73,13 @@ using namespace o2::aod::rctsel; auto static constexpr CminCharge = 3.f; auto static constexpr CintZero = 0; -auto static constexpr CintTwo = 2; auto static constexpr Czero = 0.f; auto static constexpr Cninety = 90.f; auto static constexpr ConeHeighty = 180.f; +auto static constexpr CminAccFT0A = 3.5f; +auto static constexpr CmaxAccFT0A = 4.9f; +auto static constexpr CminAccFT0C = -3.3f; +auto static constexpr CmaxAccFT0C = -2.1f; enum TrkSel { trkSelAll, @@ -106,6 +110,11 @@ enum AmbTrkType { nAmbTrkType }; +std::unordered_map mapVtxXrec; +std::unordered_map mapVtxYrec; +std::unordered_map mapVtxZrec; +std::unordered_map mapMcCollIdPerRecColl; + struct DndetaMFTPbPb { SliceCache cache; @@ -141,6 +150,7 @@ struct DndetaMFTPbPb { Configurable cfgRemoveReassigned{"cfgRemoveReassigned", false, "Remove reassgined tracks"}; Configurable cfgUseTrackParExtra{"cfgUseTrackParExtra", false, "Use table with refitted track parameters"}; Configurable cfgUseInelgt0{"cfgUseInelgt0", false, "Use INEL > 0 condition"}; + Configurable cfgUseInelgt0wMFT{"cfgUseInelgt0wMFT", false, "Use INEL > 0 condition with MFT acceptance"}; Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; Configurable cfgDCAtype{"cfgDCAtype", 2, "DCA coordinate type [0: DCA-X, 1: DCA-Y, 2: DCA-XY]"}; } gConf; @@ -217,6 +227,9 @@ struct DndetaMFTPbPb { Service ccdb; Configurable ccdbNoLaterThan{"ccdbNoLaterThan", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object"}; Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable cfgApplyZShiftFromCCDB{"cfgApplyZShiftFromCCDB", false, "flag to apply z shift from CCDB"}; + Configurable cfgZShiftPath{"cfgZShiftPath", "Users/m/mcoquet/ZShift", "CCDB path for z shift to apply to forward tracks"}; + Configurable cfgManualZShift{"cfgManualZShift", 0.0f, "manual z-shift for propagation of global muon to PV"}; int mRunNumber{-1}; uint64_t mSOR{0}; @@ -228,6 +241,7 @@ struct DndetaMFTPbPb { float bZ = 0; // Magnetic field for MFT static constexpr double CcenterMFT[3] = {0, 0, -61.4}; // Field at center of MFT + float mZShift = 0; // z-vertex shift o2::parameters::GRPMagField* grpmag = nullptr; @@ -908,27 +922,24 @@ struct DndetaMFTPbPb { if (doprocessCollAssocMC) { // tracks not associated to any collision - hCollAssoc[0] = qaregistry.add("TrackToColl/hNonAssocTracks", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis}); + hCollAssoc[0] = qaregistry.add("TrackToColl/hNonAssocTracks", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{X}_{vtx}^{reco}#minus#it{X}_{vtx}^{gen} (cm);#it{Y}_{vtx}^{reco}#minus#it{Y}_{vtx}^{gen} (cm);#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis, deltaZAxis, deltaZAxis}); // tracks associasted to a collision - hCollAssoc[1] = qaregistry.add("TrackToColl/hAssocTracks", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis}); + hCollAssoc[1] = qaregistry.add("TrackToColl/hAssocTracks", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{X}_{vtx}^{reco}#minus#it{X}_{vtx}^{gen} (cm);#it{Y}_{vtx}^{reco}#minus#it{Y}_{vtx}^{gen} (cm);#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis, deltaZAxis, deltaZAxis}); // tracks associated to the correct collision considering only first reco collision (based on the MC collision index) - hCollAssoc[2] = qaregistry.add("TrackToColl/hGoodAssocTracks", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis}); + hCollAssoc[2] = qaregistry.add("TrackToColl/hGoodAssocTracks", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{X}_{vtx}^{reco}#minus#it{X}_{vtx}^{gen} (cm);#it{Y}_{vtx}^{reco}#minus#it{Y}_{vtx}^{gen} (cm);#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis, deltaZAxis, deltaZAxis}); // tracks associated to the correct collision considering all ambiguous reco collisions (based on the MC collision index) - hCollAssoc[3] = qaregistry.add("TrackToColl/hGoodAssocTracksAmb", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis}); + hCollAssoc[3] = qaregistry.add("TrackToColl/hGoodAssocTracksAmb", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{X}_{vtx}^{reco}#minus#it{X}_{vtx}^{gen} (cm);#it{Y}_{vtx}^{reco}#minus#it{Y}_{vtx}^{gen} (cm);#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis, deltaZAxis, deltaZAxis}); qaregistry.add("TrackToColl/histFracTracksFakeMcColl", "Fraction of tracks originating from fake collision; fraction; entries", {HistType::kTH1F, {{101, 0., 1.01}}}); qaregistry.add("TrackToColl/histFracGoodTracks", "Fraction of tracks originating from the correct collision; fraction; entries", {HistType::kTH1F, {{101, 0., 1.01}}}); qaregistry.add("TrackToColl/histAmbTrackNumColls", "Number of collisions associated to an ambiguous track; no. collisions; entries", {HistType::kTH1F, {{30, -0.5, 29.5}}}); + qaregistry.add("TrackToColl/histAmbTrackZvtxRMS", "RMS of #it{Z}^{reco} of collisions associated to a track; RMS(#it{Z}^{reco}) (cm); entries", {HistType::kTH1F, {{100, 0., 0.5}}}); } if (doprocessReAssocMC) { - // tracks not associated to any collision - hReAssoc[0] = qaregistry.add("ReAssoc/hNonAssocTracks", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{X}_{vtx}^{reco}#minus#it{X}_{vtx}^{gen} (cm);#it{Y}_{vtx}^{reco}#minus#it{Y}_{vtx}^{gen} (cm);#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis, deltaZAxis, deltaZAxis}); - // tracks associasted to a collision - hReAssoc[1] = qaregistry.add("ReAssoc/hAssocTracks", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{X}_{vtx}^{reco}#minus#it{X}_{vtx}^{gen} (cm);#it{Y}_{vtx}^{reco}#minus#it{Y}_{vtx}^{gen} (cm);#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis, deltaZAxis, deltaZAxis}); - // tracks associated to the correct collision considering only first reco collision (based on the MC collision index) - hReAssoc[2] = qaregistry.add("ReAssoc/hGoodAssocTracks", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{X}_{vtx}^{reco}#minus#it{X}_{vtx}^{gen} (cm);#it{Y}_{vtx}^{reco}#minus#it{Y}_{vtx}^{gen} (cm);#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis, deltaZAxis, deltaZAxis}); - // tracks associated to the correct collision considering all ambiguous reco collisions (based on the MC collision index) - hReAssoc[3] = qaregistry.add("ReAssoc/hGoodAssocTracksAmb", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{X}_{vtx}^{reco}#minus#it{X}_{vtx}^{gen} (cm);#it{Y}_{vtx}^{reco}#minus#it{Y}_{vtx}^{gen} (cm);#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis, deltaZAxis, deltaZAxis}); + hReAssoc[0] = qaregistry.add("ReAssoc/hAssocBestTrue", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{X}_{vtx}^{reco}#minus#it{X}_{vtx}^{gen} (cm);#it{Y}_{vtx}^{reco}#minus#it{Y}_{vtx}^{gen} (cm);#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis, deltaZAxis, deltaZAxis}); + hReAssoc[1] = qaregistry.add("ReAssoc/hAssocBestWrong", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{X}_{vtx}^{reco}#minus#it{X}_{vtx}^{gen} (cm);#it{Y}_{vtx}^{reco}#minus#it{Y}_{vtx}^{gen} (cm);#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis, deltaZAxis, deltaZAxis}); + hReAssoc[2] = qaregistry.add("ReAssoc/hReAssocBestTrue", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{X}_{vtx}^{reco}#minus#it{X}_{vtx}^{gen} (cm);#it{Y}_{vtx}^{reco}#minus#it{Y}_{vtx}^{gen} (cm);#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis, deltaZAxis, deltaZAxis}); + hReAssoc[3] = qaregistry.add("ReAssoc/hReAssocBestWrong", ";#it{p}_{T}^{reco} (GeV/#it{c});#it{#eta}^{reco};#it{X}_{vtx}^{reco}#minus#it{X}_{vtx}^{gen} (cm);#it{Y}_{vtx}^{reco}#minus#it{Y}_{vtx}^{gen} (cm);#it{Z}_{vtx}^{reco}#minus#it{Z}_{vtx}^{gen} (cm)", HistType::kTHnSparseF, {ptAxis, etaAxis, deltaZAxis, deltaZAxis, deltaZAxis}); } if (doprocessEfficiencyInclusive) { @@ -1073,12 +1084,11 @@ struct DndetaMFTPbPb { } } - if (doprocessReassocEfficiencyCentFT0C) { - auto hNevt = registry.add("Events/hNReEffColls", "Number of generated and reconstructed MC collisions", HistType::kTH1F, {{3, 0.5, 3.5}}); - hNevt->GetXaxis()->SetBinLabel(1, "Reconstructed collisions"); - hNevt->GetXaxis()->SetBinLabel(2, "Generated collisions"); + if (doprocessReassocEfficiency) { + registry.add("Events/hNReEffColls", "Number of times generated collisions are reconstructed; N; Counts", HistType::kTH1F, {{10, -0.5, 9.5}}); + registry.add("Events/hNchTVX", "; status;", HistType::kTH2F, {multAxis, {2, 0, 2}}); - registry.add("Events/Centrality/ReEffStatus", ";status;centrality", HistType::kTH2F, {{8, 0.5, 8.5}, centralityAxis}); + registry.add("Events/Centrality/ReEffStatus", ";status;centrality", HistType::kTH2F, {{9, 0.5, 9.5}, centralityAxis}); auto hstat = registry.get(HIST("Events/Centrality/ReEffStatus")); hstat->GetXaxis()->SetBinLabel(1, "All tracks"); hstat->GetXaxis()->SetBinLabel(2, "Ambiguous tracks"); @@ -1086,8 +1096,11 @@ struct DndetaMFTPbPb { hstat->GetXaxis()->SetBinLabel(4, "Extra tracks"); hstat->GetXaxis()->SetBinLabel(5, "Originally correctly reassgined"); hstat->GetXaxis()->SetBinLabel(6, "Correctly reassigned"); - hstat->GetXaxis()->SetBinLabel(7, "Not reassigned"); - hstat->GetXaxis()->SetBinLabel(8, "Correctly assigned true"); + hstat->GetXaxis()->SetBinLabel(7, "Not reassigned (reassigned)"); + hstat->GetXaxis()->SetBinLabel(8, "Not reassigned"); + hstat->GetXaxis()->SetBinLabel(9, "Correctly assigned true"); + + registry.add({"AmbTracks/hVtxzMCrec", " ; Z_{vtx} (cm)", {HistType::kTH1F, {zAxis}}}); registry.add({"AmbTracks/DCAXY", " ; DCA_{XY} (cm)", {HistType::kTH1F, {dcaxyAxis}}}); registry.add({"AmbTracks/DCAZ", " ; DCA_{Z} (cm)", {HistType::kTH1F, {dcazAxis}}}); @@ -1110,6 +1123,13 @@ struct DndetaMFTPbPb { registry.add({"AmbTracks/Centrality/THnDCAxyBestFalse", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm)", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis, centralityAxis}}}); registry.add({"AmbTracks/Centrality/THnDCAxyBestTrueOrigAssoc", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm)", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis, centralityAxis}}}); registry.add({"AmbTracks/Centrality/THnDCAxyBestTrueOrigReAssoc", "; p_{T} (GeV/c); #eta; Z_{vtx} (cm); DCA_{XY} (cm); DCA_{Z} (cm)", {HistType::kTHnSparseF, {ptAxis, etaAxis, zAxis, dcaxyAxis, dcazAxis, centralityAxis}}}); + + registry.add({"AmbTracks/BestGenDxyz", "; #Delta X (cm); #Delta Y (cm); #Delta Z (cm); DCA_{XY} (cm); DCA_{Z} (cm)", {HistType::kTHnSparseF, {deltaZAxis, deltaZAxis, deltaZAxis, dcaxyAxis, dcazAxis}}}); + registry.add({"AmbTracks/BestPrimDxyz", "; #Delta X (cm); #Delta Y (cm); #Delta Z (cm); DCA_{XY} (cm); DCA_{Z} (cm)", {HistType::kTHnSparseF, {deltaZAxis, deltaZAxis, deltaZAxis, dcaxyAxis, dcazAxis}}}); + registry.add({"AmbTracks/BestTrueDxyz", "; #Delta X (cm); #Delta Y (cm); #Delta Z (cm); DCA_{XY} (cm); DCA_{Z} (cm)", {HistType::kTHnSparseF, {deltaZAxis, deltaZAxis, deltaZAxis, dcaxyAxis, dcazAxis}}}); + registry.add({"AmbTracks/BestFalseDxyz", "; #Delta X (cm); #Delta Y (cm); #Delta Z (cm); DCA_{XY} (cm); DCA_{Z} (cm)", {HistType::kTHnSparseF, {deltaZAxis, deltaZAxis, deltaZAxis, dcaxyAxis, dcazAxis}}}); + registry.add({"AmbTracks/BestTrueOrigReAssocDxyz", "; #Delta X (cm); #Delta Y (cm); #Delta Z (cm); DCA_{XY} (cm); DCA_{Z} (cm)", {HistType::kTHnSparseF, {deltaZAxis, deltaZAxis, deltaZAxis, dcaxyAxis, dcazAxis}}}); + registry.add({"AmbTracks/BestTrueOrigAssocDxyz", "; #Delta X (cm); #Delta Y (cm); #Delta Z (cm); DCA_{XY} (cm); DCA_{Z} (cm)", {HistType::kTHnSparseF, {deltaZAxis, deltaZAxis, deltaZAxis, dcaxyAxis, dcazAxis}}}); } if (doprocessEventAndSignalLossCentFT0C) { @@ -1184,6 +1204,22 @@ struct DndetaMFTPbPb { using FiltParticles = soa::Filtered; + /// \brief RMS calculation + /// \param vec vector of values to compute RMS + template + T computeRMS(std::vector& vec) + { + T sum = std::accumulate(vec.begin(), vec.end(), 0.0); + T mean = sum / vec.size(); + + std::vector diff(vec.size()); + std::transform(vec.begin(), vec.end(), diff.begin(), [mean](T x) { return x - mean; }); + T sqSum = std::inner_product(diff.begin(), diff.end(), diff.begin(), 0.0); + T stdev = std::sqrt(sqSum / vec.size()); + + return stdev; + } + void initCCDB(ExtBCs::iterator const& bc) { if (mRunNumber == bc.runNumber()) { @@ -1200,6 +1236,20 @@ struct DndetaMFTPbPb { o2::field::MagneticField* field = static_cast(TGeoGlobalMagField::Instance()->GetField()); bZ = field->getBz(CcenterMFT); LOG(info) << "The field at the center of the MFT is bZ = " << bZ; + + if (cfgApplyZShiftFromCCDB) { + auto* zShift = ccdb->getForTimeStamp>(cfgZShiftPath, bc.timestamp()); + if (zShift != nullptr && !zShift->empty()) { + LOGF(info, "reading z shift %f from %s", (*zShift)[0], cfgZShiftPath.value); + mZShift = (*zShift)[0]; + } else { + LOGF(info, "z shift is not found in ccdb path %s. set to 0 cm", cfgZShiftPath.value); + mZShift = 0; + } + } else { + LOGF(info, "z shift is manually set to %f cm", cfgManualZShift.value); + mZShift = cfgManualZShift; + } } template @@ -1501,6 +1551,46 @@ struct DndetaMFTPbPb { return nCharged; } + template + bool isInelGt0wMft(P const& particles) + { + int nChrgMc = 0; + int nChrgFT0A = 0; + int nChrgFT0C = 0; + for (auto const& particle : particles) { + if (!isChrgParticle(particle.pdgCode())) { + continue; + } + if (!particle.isPhysicalPrimary()) { + continue; + } + // trigger TVX + if (particle.eta() > CminAccFT0A && particle.eta() < CmaxAccFT0A) { + nChrgFT0A++; + } + if (particle.eta() > CminAccFT0C && particle.eta() < CmaxAccFT0C) { + nChrgFT0C++; + } + // acceptance MFT + if (particle.eta() < trackCuts.minEta || particle.eta() > trackCuts.maxEta) { + continue; + } + nChrgMc++; + } + + if (nChrgFT0A == CintZero || nChrgFT0C == CintZero) { + registry.fill(HIST("Events/hNchTVX"), nChrgMc, 0.5); + return false; + } + registry.fill(HIST("Events/hNchTVX"), nChrgMc, 1.5); + + if (nChrgMc == CintZero) { + return false; + } + + return true; + } + template bool isParticleSelected(P const& particle) { @@ -2934,37 +3024,51 @@ struct DndetaMFTPbPb { qaregistry.fill(HIST("TrackToColl/histFracTracksFakeMcColl"), fracFake); } - for (auto const& track : tracks) { + for (const auto& track : tracks) { uint index = uint(track.collisionId() >= 0); if (track.has_mcParticle()) { // auto particle = track.mcParticle_as(); - auto particle = track.mcParticle(); - bool isAmbiguous = (track.compatibleCollIds().size() != 1); + const auto& particle = track.mcParticle(); + bool isAmbiguous = (track.compatibleCollIds().size() > 1); if (isAmbiguous) { qaregistry.fill(HIST("TrackToColl/histAmbTrackNumColls"), track.compatibleCollIds().size()); + + std::vector ambVtxZ{}; + for (const auto& collIdx : track.compatibleCollIds()) { + const auto& ambColl = collisions.rawIteratorAt(collIdx); + ambVtxZ.push_back(ambColl.posZ()); + } + if (!ambVtxZ.empty()) { + qaregistry.fill(HIST("TrackToColl/histAmbTrackZvtxRMS"), computeRMS(ambVtxZ)); + } } + float deltaX = -999.f; + float deltaY = -999.f; float deltaZ = -999.f; if (index) { - auto collision = track.collision_as(); - auto mcCollision = particle.mcCollision_as(); + const auto& collision = track.collision_as(); + const auto& mcCollision = particle.mcCollision_as(); + deltaX = collision.posX() - mcCollision.posX(); + deltaY = collision.posY() - mcCollision.posY(); deltaZ = collision.posZ() - mcCollision.posZ(); if (collision.has_mcCollision() && collision.mcCollisionId() == particle.mcCollisionId()) { - hCollAssoc[index + 1]->Fill(track.pt(), track.eta(), deltaZ); + hCollAssoc[index + 1]->Fill(track.pt(), track.eta(), deltaX, deltaY, deltaZ); } else { if (isAmbiguous) { for (const auto& collIdx : track.compatibleCollIds()) { - auto ambCollision = collisions.rawIteratorAt(collIdx); - if (ambCollision.has_mcCollision() && ambCollision.mcCollisionId() == particle.mcCollisionId()) { - hCollAssoc[index + 2]->Fill(track.pt(), track.eta(), deltaZ); + auto ambColl = collisions.rawIteratorAt(collIdx); + if (ambColl.has_mcCollision() && ambColl.mcCollisionId() == particle.mcCollisionId()) { + hCollAssoc[index + 2]->Fill(track.pt(), track.eta(), deltaX, deltaY, deltaZ); + // hCollAssoc[index + 2]->Fill(track.pt(), track.eta(), ambColl.posX() - mcCollision.posX(), ambColl.posY() - mcCollision.posY(), ambColl.posZ() - mcCollision.posZ()); break; } } } } + hCollAssoc[index]->Fill(track.pt(), track.eta(), deltaX, deltaY, deltaZ); } - hCollAssoc[index]->Fill(track.pt(), track.eta(), deltaZ); } else { - hCollAssoc[index]->Fill(track.pt(), track.eta(), -999.f); + hCollAssoc[index]->Fill(track.pt(), track.eta(), -999.f, -999.f, -999.f); } } } @@ -2979,50 +3083,136 @@ struct DndetaMFTPbPb { PROCESS_SWITCH(DndetaMFTPbPb, processCollAssocMC, "Process collision-association information, requires extra table from TrackToCollisionAssociation task (fillTableOfCollIdsPerTrack=true)", false); template - void processCheckReAssocMC(C const& /*collisions*/, + void processCheckReAssocMC(C const& collisions, soa::SmallGroups const& besttracks, FiltMcMftTracks const& /*tracks*/, FiltParticles const& /*particles*/, aod::McCollisions const& /*mccollisions*/ ) { - for (auto const& track : besttracks) { - uint index = uint(track.bestCollisionId() >= 0); // assigned - if (!isBestTrackSelected(track)) { - continue; + const auto& nRecoColls = collisions.size(); + // Generated evets with >= 1 reco collisions + if (nRecoColls > CintZero) { + + mapVtxXrec.clear(); + mapVtxYrec.clear(); + mapVtxZrec.clear(); + mapMcCollIdPerRecColl.clear(); + mapVtxXrec.reserve(collisions.size()); + mapVtxYrec.reserve(collisions.size()); + mapVtxZrec.reserve(collisions.size()); + mapMcCollIdPerRecColl.reserve(collisions.size()); + + auto maxNcontributors = -1; + auto bestCollIndex = -1; + for (auto const& collision : collisions) { + if (!isGoodEvent(collision)) { + continue; + } + if (!collision.has_mcCollision()) { + continue; + } + if (maxNcontributors < collision.numContrib()) { + maxNcontributors = collision.numContrib(); + bestCollIndex = collision.globalIndex(); + mapVtxXrec.emplace(collision.globalIndex(), collision.posX()); + mapVtxYrec.emplace(collision.globalIndex(), collision.posY()); + mapVtxZrec.emplace(collision.globalIndex(), collision.posZ()); + mapMcCollIdPerRecColl.emplace(collision.globalIndex(), collision.mcCollisionId()); + } else { + continue; + } } - auto itrack = track.mfttrack_as(); - if (gConf.cfgRemoveReassigned) { - if (itrack.collisionId() != track.bestCollisionId()) { + for (const auto& collision : collisions) { + if (!isGoodEvent(collision)) { + continue; + } + if (!collision.has_mcCollision()) { + continue; + } + // Select collisions with the largest number of contributors + if (bestCollIndex != collision.globalIndex()) { continue; } - } - if (!isTrackSelected(itrack)) { - continue; - } - if (itrack.has_mcParticle()) { - auto particle = itrack.mcParticle_as(); - float deltaX = -999.f; - float deltaY = -999.f; - float deltaZ = -999.f; - if (index) { - auto collision = itrack.collision_as(); - auto mcCollision = particle.mcCollision_as(); - deltaX = collision.posX() - mcCollision.posX(); - deltaY = collision.posY() - mcCollision.posY(); - deltaZ = collision.posZ() - mcCollision.posZ(); - if (collision.has_mcCollision() && collision.mcCollisionId() == particle.mcCollisionId()) { - hReAssoc[index + 1]->Fill(itrack.pt(), itrack.eta(), deltaX, deltaY, deltaZ); - } else { - hReAssoc[index + 2]->Fill(itrack.pt(), itrack.eta(), deltaX, deltaY, deltaZ); + auto perCollisionASample = besttracks.sliceBy(perColU, collision.globalIndex()); + for (auto const& atrack : perCollisionASample) { + if (!isBestTrackSelected(atrack)) { + continue; + } + auto itrack = atrack.template mfttrack_as(); + if (!isTrackSelected(itrack)) { + continue; + } + float phi = itrack.phi(); + o2::math_utils::bringTo02Pi(phi); + if (phi < Czero || TwoPI < phi) { + continue; + } + if (!itrack.has_collision()) { + continue; + } + if (gConf.cfgRemoveReassigned) { + if (itrack.collisionId() != atrack.bestCollisionId()) { + continue; + } + } + if (itrack.has_mcParticle()) { + auto particle = itrack.template mcParticle_as(); + auto collision = itrack.template collision_as(); + auto mcCollision = particle.template mcCollision_as(); + + if (eventCuts.useZVtxCutMC && (std::abs(mcCollision.posZ()) >= eventCuts.maxZvtx)) { + continue; + } + + float deltaX = -999.f; + float deltaY = -999.f; + float deltaZ = -999.f; + + if (mapVtxXrec.find(atrack.bestCollisionId()) == mapVtxXrec.end()) { + continue; + } + if (mapVtxYrec.find(atrack.bestCollisionId()) == mapVtxYrec.end()) { + continue; + } + if (mapVtxZrec.find(atrack.bestCollisionId()) == mapVtxZrec.end()) { + continue; + } + if (mapMcCollIdPerRecColl.find(atrack.bestCollisionId()) == mapMcCollIdPerRecColl.end()) { + continue; + } + const float vtxXbest = mapVtxXrec.find(atrack.bestCollisionId())->second; + const float vtxYbest = mapVtxYrec.find(atrack.bestCollisionId())->second; + const float vtxZbest = mapVtxZrec.find(atrack.bestCollisionId())->second; + // LOGP(info, "\t ---> \t .... \t vtxZrec: {} - collision.posZ(): {}", vtxZrec, collision.posZ()); + const float mcCollIdRec = mapMcCollIdPerRecColl.find(atrack.bestCollisionId())->second; + // LOGP(info, "\t ---> \t .... \t mcCollIdRec: {} - bestMCCol: {}", mcCollIdRec, bestMCCol); + deltaX = vtxXbest - mcCollision.posX(); + deltaY = vtxYbest - mcCollision.posY(); + deltaZ = vtxZbest - mcCollision.posZ(); + + if (itrack.collisionId() == atrack.bestCollisionId()) { // associated + if (collision.has_mcCollision() && mcCollIdRec == particle.mcCollisionId()) { + hReAssoc[0]->Fill(itrack.pt(), itrack.eta(), deltaX, deltaY, deltaZ); + } else { + hReAssoc[1]->Fill(itrack.pt(), itrack.eta(), deltaX, deltaY, deltaZ); + } + } else { + if (collision.has_mcCollision() && mcCollIdRec == particle.mcCollisionId()) { + hReAssoc[2]->Fill(itrack.pt(), itrack.eta(), deltaX, deltaY, deltaZ); + } else { + hReAssoc[3]->Fill(itrack.pt(), itrack.eta(), deltaX, deltaY, deltaZ); + } + } } } - hReAssoc[index]->Fill(itrack.pt(), itrack.eta(), deltaX, deltaY, deltaZ); - } else { - hReAssoc[index]->Fill(itrack.pt(), itrack.eta(), -999.f, -999.f, -999.f); } + mapVtxXrec.clear(); + mapVtxYrec.clear(); + mapVtxZrec.clear(); + mapMcCollIdPerRecColl.clear(); } } @@ -3266,8 +3456,8 @@ struct DndetaMFTPbPb { } else { if (isAmbiguous) { for (const auto& collIdx : track.compatibleCollIds()) { - auto ambCollision = collisions.rawIteratorAt(collIdx); - if (ambCollision.has_mcCollision() && ambCollision.mcCollisionId() == particle.mcCollisionId()) { + auto ambColl = collisions.rawIteratorAt(collIdx); + if (ambColl.has_mcCollision() && ambColl.mcCollisionId() == particle.mcCollisionId()) { if (!particle.isPhysicalPrimary()) { // Secondaries (weak decays and material) if constexpr (has_reco_cent) { registry.fill(HIST("Tracks/Centrality/THnGenSecAmb"), particle.pt(), particle.eta(), particle.mcCollision().posZ(), crec); @@ -3649,11 +3839,10 @@ struct DndetaMFTPbPb { Preslice compCollPerCol = o2::aod::fwdtrack::collisionId; /// @brief process function to check reassociation efficiency - template - void processReassocEfficiency(typename soa::Join const& collisions, + void processReassocEfficiency(CollsMCExtra::iterator const& mcCollision, + soa::SmallGroups const& collisions, MftTracksWCollsMC const& tracks, - MC const& mcCollisions, - aod::McParticles const& /*particles*/, + aod::McParticles const& particles, ExtBCs const& bcs) { if (bcs.size() == 0) { @@ -3662,158 +3851,203 @@ struct DndetaMFTPbPb { auto bc = bcs.begin(); initCCDB(bc); - registry.fill(HIST("Events/hNReEffColls"), 1.f, collisions.size()); - registry.fill(HIST("Events/hNReEffColls"), 2.f, mcCollisions.size()); - - for (const auto& collision : collisions) { - if (!isGoodEvent(collision)) { - continue; - } - float crec = getRecoCent(collision); - registry.fill(HIST("Events/Centrality/ReEffStatus"), 1, crec); - - if (!collision.has_mcCollision()) { - continue; - } + // At least one generated primary in MFT acceptance + TVX triggered collisions + if (gConf.cfgUseInelgt0wMFT && !isInelGt0wMft(particles)) { + return; + } + if (eventCuts.useZVtxCutMC && (std::abs(mcCollision.posZ()) >= eventCuts.maxZvtx)) { + return; + } - std::array dcaInfOrig; - std::array dcaInfo; - double bestDCA[2]; - auto trkPerColl = tracks.sliceBy(compCollPerCol, collision.globalIndex()); + const auto& nRecoColls = collisions.size(); + registry.fill(HIST("Events/hNReEffColls"), 1.f, nRecoColls); - for (auto const& track : trkPerColl) { - dcaInfOrig[0] = 999.f; // original DCAx from propagation - dcaInfOrig[1] = 999.f; // original DCAy from propagation - dcaInfOrig[2] = 999.f; // original DCAz from propagation - dcaInfo[0] = 999.f; // calcualted DCAxy - dcaInfo[1] = 999.f; // calculated DCAz - same as original - bestDCA[0] = 999.f; // minimal DCAxy - bestDCA[1] = 999.f; // minimal DCAz + // Generated evets with >= 1 reco collisions + if (nRecoColls > CintZero) { - if (!isTrackSelected(track)) { + auto maxNcontributors = -1; + auto bestCollIndex = -1; + auto crec = -1; + for (auto const& collision : collisions) { + if (!isGoodEvent(collision)) { continue; } - - auto bestCol = track.collisionId(); - uint index = uint(track.collisionId() >= 0); - auto ids = track.compatibleCollIds(); - bool isAmbiguous = (ids.size() > 1); - - if (!track.has_mcParticle()) { + if (maxNcontributors < collision.numContrib()) { + maxNcontributors = collision.numContrib(); + bestCollIndex = collision.globalIndex(); + crec = getRecoCent(collision); + } else { continue; } - auto particle = track.template mcParticle_as(); - if (!isChrgParticle(particle.pdgCode())) { + } + + registry.fill(HIST("AmbTracks/hVtxzMCrec"), mcCollision.posZ()); + + for (const auto& collision : collisions) { + if (!isGoodEvent(collision)) { continue; } - if (particle.eta() <= trackCuts.minEta || particle.eta() >= trackCuts.maxEta) { + // Select collisions with the largest number of contributors + if (bestCollIndex != collision.globalIndex()) { continue; } - if (gConf.cfgUseParticleSel && !isParticleSelected(particle)) { + + registry.fill(HIST("Events/Centrality/ReEffStatus"), 1, crec); + + if (!collision.has_mcCollision()) { continue; } - int bestMCCol = -1; - std::vector v1; // Temporary null vector for the computation of the covariance matrix - SMatrix55 tcovs(v1.begin(), v1.end()); - SMatrix5 tpars(track.x(), track.y(), track.phi(), track.tgl(), track.signed1Pt()); - o2::track::TrackParCovFwd trackPar{track.z(), tpars, tcovs, track.chi2()}; + std::array dcaInfOrig; + std::array dcaInfo; + double bestDCA[2]; + auto trkPerColl = tracks.sliceBy(compCollPerCol, collision.globalIndex()); + + for (auto const& track : trkPerColl) { + dcaInfOrig[0] = 999.f; // original DCAx from propagation + dcaInfOrig[1] = 999.f; // original DCAy from propagation + dcaInfOrig[2] = 999.f; // original DCAz from propagation + dcaInfo[0] = 999.f; // calcualted DCAxy + dcaInfo[1] = 999.f; // calculated DCAz - same as original + bestDCA[0] = 999.f; // minimal DCAxy + bestDCA[1] = 999.f; // minimal DCAz + + if (!isTrackSelected(track)) { + continue; + } - if (index) { - auto mcCollision = particle.template mcCollision_as(); - if (eventCuts.useZDiffCut) { - if (std::abs(collision.posZ() - mcCollision.posZ()) > eventCuts.maxZvtxDiff) { + auto bestCol = track.has_collision() ? track.collisionId() : -1; + uint index = uint(track.collisionId() >= 0); + bool isAmbiguous = (track.compatibleCollIds().size() > 1); + + if (!track.has_mcParticle()) { + continue; + } + auto particle = track.mcParticle_as(); + if (!isChrgParticle(particle.pdgCode())) { + continue; + } + if (particle.eta() <= trackCuts.minEta || particle.eta() >= trackCuts.maxEta) { + continue; + } + if (gConf.cfgUseParticleSel && !isParticleSelected(particle)) { + continue; + } + + if (gConf.cfgRemoveReassigned) { + if (track.compatibleCollIds().empty() || (track.compatibleCollIds().size() == 1 && bestCol == track.compatibleCollIds()[0])) { continue; } } - if (isAmbiguous) { - for (const auto& collIdx : track.compatibleCollIds()) { - auto ambCollision = collisions.rawIteratorAt(collIdx); - if (!ambCollision.has_mcCollision()) { + int bestMCCol = -1; + o2::track::TrackParCovFwd trackPar = o2::aod::fwdtrackutils::getTrackParCovFwdShift(track, mZShift); + + if (index) { + auto mcCollision = particle.mcCollision_as(); + if (eventCuts.useZDiffCut) { + if (std::abs(collision.posZ() - mcCollision.posZ()) > eventCuts.maxZvtxDiff) { continue; } + } - trackPar.propagateToDCAhelix(bZ, {ambCollision.posX(), ambCollision.posY(), ambCollision.posZ()}, dcaInfOrig); + if (isAmbiguous) { + for (const auto& collIdx : track.compatibleCollIds()) { + auto ambColl = collisions.rawIteratorAt(collIdx); + if (!ambColl.has_mcCollision()) { + continue; + } - if (gConf.cfgDCAtype == 0) { - dcaInfo[0] = dcaInfOrig[0]; - } else if (gConf.cfgDCAtype == 1) { - dcaInfo[0] = dcaInfOrig[1]; - } else if (gConf.cfgDCAtype == CintTwo) { - dcaInfo[0] = std::sqrt(dcaInfOrig[0] * dcaInfOrig[0] + dcaInfOrig[1] * dcaInfOrig[1]); - } - dcaInfo[1] = dcaInfOrig[2]; + trackPar.propagateToDCAhelix(bZ, {ambColl.posX(), ambColl.posY(), ambColl.posZ()}, dcaInfOrig); - // LOGP(info, "-> track {}: {}", track.globalIndex(), dcaInfo[0]); - registry.fill(HIST("AmbTracks/DCAXY"), dcaInfo[0]); - registry.fill(HIST("AmbTracks/DCAZ"), dcaInfo[1]); + if (gConf.cfgDCAtype == 0) { + dcaInfo[0] = dcaInfOrig[0]; + } else if (gConf.cfgDCAtype == 1) { + dcaInfo[0] = dcaInfOrig[1]; + } else if (gConf.cfgDCAtype == 2) { + dcaInfo[0] = std::sqrt(dcaInfOrig[0] * dcaInfOrig[0] + dcaInfOrig[1] * dcaInfOrig[1]); + } + dcaInfo[1] = dcaInfOrig[2]; - if ((std::abs(dcaInfo[0]) < std::abs(bestDCA[0])) && (std::abs(dcaInfo[1]) < std::abs(bestDCA[1]))) { - bestCol = ambCollision.globalIndex(); - bestMCCol = ambCollision.mcCollisionId(); - bestDCA[0] = dcaInfo[0]; - bestDCA[1] = dcaInfo[1]; - } - } + if ((std::abs(dcaInfo[0]) < std::abs(bestDCA[0])) && (std::abs(dcaInfo[1]) < std::abs(bestDCA[1]))) { + bestCol = ambColl.globalIndex(); + bestMCCol = ambColl.mcCollisionId(); + bestDCA[0] = dcaInfo[0]; + bestDCA[1] = dcaInfo[1]; + } - registry.fill(HIST("AmbTracks/DCAXYBest"), bestDCA[0]); - registry.fill(HIST("AmbTracks/DCAZBest"), bestDCA[1]); - registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestGen"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); + // LOGP(info, "-> track {}: {}", track.globalIndex(), dcaInfo[0]); + registry.fill(HIST("AmbTracks/DCAXY"), dcaInfo[0]); + registry.fill(HIST("AmbTracks/DCAZ"), dcaInfo[1]); + } - if (particle.isPhysicalPrimary()) { - registry.fill(HIST("AmbTracks/DCAXYBestPrim"), bestDCA[0]); - registry.fill(HIST("AmbTracks/DCAZBestPrim"), bestDCA[1]); - registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestGenPrim"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); - } + registry.fill(HIST("AmbTracks/DCAXYBest"), bestDCA[0]); + registry.fill(HIST("AmbTracks/DCAZBest"), bestDCA[1]); + registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestGen"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); + registry.fill(HIST("AmbTracks/BestGenDxyz"), collision.posX() - mcCollision.posX(), collision.posY() - mcCollision.posY(), collision.posZ() - mcCollision.posZ(), bestDCA[0], bestDCA[1]); - auto mcCollID = particle.mcCollisionId(); - registry.fill(HIST("Events/Centrality/ReEffStatus"), 2, crec); - - if (!track.has_collision()) { - registry.fill(HIST("Events/Centrality/ReEffStatus"), 4, crec); - if (bestMCCol == mcCollID) // correctly assigned to bestCol - { - registry.fill(HIST("Events/Centrality/ReEffStatus"), 6, crec); - registry.fill(HIST("AmbTracks/DCAXYBestTrue"), bestDCA[0]); - registry.fill(HIST("AmbTracks/DCAZBestTrue"), bestDCA[1]); - registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestTrue"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); - } else { - registry.fill(HIST("AmbTracks/DCAXYBestFalse"), bestDCA[0]); - registry.fill(HIST("AmbTracks/DCAZBestFalse"), bestDCA[1]); - registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestFalse"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); - } - } else if (track.collisionId() != bestCol) { // reassgined - auto collOrig = track.template collision_as>(); - registry.fill(HIST("Events/Centrality/ReEffStatus"), 3, crec); - if (bestMCCol == mcCollID) // correctly reassigned - { - registry.fill(HIST("Events/Centrality/ReEffStatus"), 6, crec); - registry.fill(HIST("AmbTracks/DCAXYBestTrue"), bestDCA[0]); - registry.fill(HIST("AmbTracks/DCAZBestTrue"), bestDCA[1]); - registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestTrue"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); - } else { // uncorrectly reassigned - registry.fill(HIST("AmbTracks/DCAXYBestFalse"), bestDCA[0]); - registry.fill(HIST("AmbTracks/DCAZBestFalse"), bestDCA[1]); - registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestFalse"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); + if (particle.isPhysicalPrimary()) { + registry.fill(HIST("AmbTracks/DCAXYBestPrim"), bestDCA[0]); + registry.fill(HIST("AmbTracks/DCAZBestPrim"), bestDCA[1]); + registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestGenPrim"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); + registry.fill(HIST("AmbTracks/BestPrimDxyz"), collision.posX() - mcCollision.posX(), collision.posY() - mcCollision.posY(), collision.posZ() - mcCollision.posZ(), bestDCA[0], bestDCA[1]); } - if (collOrig.mcCollisionId() == mcCollID) { // initially correctly assigned - registry.fill(HIST("Events/Centrality/ReEffStatus"), 5, crec); - registry.fill(HIST("AmbTracks/DCAXYBestTrueOrigReAssoc"), bestDCA[0]); - registry.fill(HIST("AmbTracks/DCAZBestTrueOrigReAssoc"), bestDCA[1]); - registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestTrueOrigReAssoc"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); - } - } else { // the track has a collision and track.collisionId() == bestCol - if (track.collisionId() != bestCol) { - // LOGP(info, "-> track id {}: bestCollid {}", track.collisionId(), bestCol); - } - registry.fill(HIST("Events/Centrality/ReEffStatus"), 7, crec); - if (bestMCCol == mcCollID) // correctly assigned - { + + auto mcCollID = particle.mcCollisionId(); + registry.fill(HIST("Events/Centrality/ReEffStatus"), 2, crec); + + if (!track.has_collision()) { + registry.fill(HIST("Events/Centrality/ReEffStatus"), 4, crec); + if (bestMCCol == mcCollID) // correctly assigned to bestCol + { + registry.fill(HIST("Events/Centrality/ReEffStatus"), 6, crec); + registry.fill(HIST("AmbTracks/DCAXYBestTrue"), bestDCA[0]); + registry.fill(HIST("AmbTracks/DCAZBestTrue"), bestDCA[1]); + registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestTrue"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); + registry.fill(HIST("AmbTracks/BestTrueDxyz"), collision.posX() - mcCollision.posX(), collision.posY() - mcCollision.posY(), collision.posZ() - mcCollision.posZ(), bestDCA[0], bestDCA[1]); + } else { + registry.fill(HIST("AmbTracks/DCAXYBestFalse"), bestDCA[0]); + registry.fill(HIST("AmbTracks/DCAZBestFalse"), bestDCA[1]); + registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestFalse"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); + registry.fill(HIST("AmbTracks/BestFalseDxyz"), collision.posX() - mcCollision.posX(), collision.posY() - mcCollision.posY(), collision.posZ() - mcCollision.posZ(), bestDCA[0], bestDCA[1]); + } + } else if (track.collisionId() != bestCol) { // reassgined + auto collOrig = collisions.rawIteratorAt(track.collisionId()); + registry.fill(HIST("Events/Centrality/ReEffStatus"), 3, crec); + if (bestMCCol == mcCollID) // correctly reassigned + { + registry.fill(HIST("Events/Centrality/ReEffStatus"), 6, crec); + registry.fill(HIST("AmbTracks/DCAXYBestTrue"), bestDCA[0]); + registry.fill(HIST("AmbTracks/DCAZBestTrue"), bestDCA[1]); + registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestTrue"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); + registry.fill(HIST("AmbTracks/BestTrueDxyz"), collision.posX() - mcCollision.posX(), collision.posY() - mcCollision.posY(), collision.posZ() - mcCollision.posZ(), bestDCA[0], bestDCA[1]); + } else { // uncorrectly reassigned + registry.fill(HIST("AmbTracks/DCAXYBestFalse"), bestDCA[0]); + registry.fill(HIST("AmbTracks/DCAZBestFalse"), bestDCA[1]); + registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestFalse"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); + registry.fill(HIST("AmbTracks/BestFalseDxyz"), collision.posX() - mcCollision.posX(), collision.posY() - mcCollision.posY(), collision.posZ() - mcCollision.posZ(), bestDCA[0], bestDCA[1]); + } + if (collOrig.mcCollisionId() == mcCollID) { // initially correctly assigned + registry.fill(HIST("Events/Centrality/ReEffStatus"), 5, crec); + registry.fill(HIST("AmbTracks/DCAXYBestTrueOrigReAssoc"), bestDCA[0]); + registry.fill(HIST("AmbTracks/DCAZBestTrueOrigReAssoc"), bestDCA[1]); + registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestTrueOrigReAssoc"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); + registry.fill(HIST("AmbTracks/BestTrueOrigReAssocDxyz"), collision.posX() - mcCollision.posX(), collision.posY() - mcCollision.posY(), collision.posZ() - mcCollision.posZ(), bestDCA[0], bestDCA[1]); + } + } else { // not reassigned - the track has a collision and track.collisionId() == bestCol + if (track.collisionId() != bestCol) { + registry.fill(HIST("Events/Centrality/ReEffStatus"), 7, crec); + // LOGP(info, "-> track id {}: bestCollid {}", track.collisionId(), bestCol); + } registry.fill(HIST("Events/Centrality/ReEffStatus"), 8, crec); - registry.fill(HIST("AmbTracks/DCAXYBestTrueOrigAssoc"), bestDCA[0]); - registry.fill(HIST("AmbTracks/DCAZBestTrueOrigAssoc"), bestDCA[1]); - registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestTrueOrigAssoc"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); + if (bestMCCol == mcCollID) // correctly assigned + { + registry.fill(HIST("Events/Centrality/ReEffStatus"), 9, crec); + registry.fill(HIST("AmbTracks/DCAXYBestTrueOrigAssoc"), bestDCA[0]); + registry.fill(HIST("AmbTracks/DCAZBestTrueOrigAssoc"), bestDCA[1]); + registry.fill(HIST("AmbTracks/Centrality/THnDCAxyBestTrueOrigAssoc"), particle.pt(), particle.eta(), mcCollision.posZ(), bestDCA[0], bestDCA[1], crec); + registry.fill(HIST("AmbTracks/BestTrueOrigAssocDxyz"), collision.posX() - mcCollision.posX(), collision.posY() - mcCollision.posY(), collision.posZ() - mcCollision.posZ(), bestDCA[0], bestDCA[1]); + } } } } @@ -3822,16 +4056,7 @@ struct DndetaMFTPbPb { } } - void processReassocEfficiencyCentFT0C(soa::Join const& collisions, - MftTracksWCollsMC const& tracks, - aod::McCollisions const& mccollisions, - aod::McParticles const& particles, - ExtBCs const& bcs) - { - processReassocEfficiency(collisions, tracks, mccollisions, particles, bcs); - } - - PROCESS_SWITCH(DndetaMFTPbPb, processReassocEfficiencyCentFT0C, "Process collision-reassociation efficiency based on track propagation DCA 3D (in FT0C centrality bins)", false); + PROCESS_SWITCH(DndetaMFTPbPb, processReassocEfficiency, "Process collision-reassociation efficiency based on track propagation DCA 3D (in FT0C centrality bins)", false); /// @brief process function to calculate signal loss based on MC void processEventAndSignalLossCentFT0C(CollsMCExtra::iterator const& mcCollision,