diff --git a/PWGLF/DataModel/Reduced3BodyTables.h b/PWGLF/DataModel/Reduced3BodyTables.h index a11081cd848..9f963951cf0 100644 --- a/PWGLF/DataModel/Reduced3BodyTables.h +++ b/PWGLF/DataModel/Reduced3BodyTables.h @@ -21,18 +21,18 @@ namespace o2::aod { -DECLARE_SOA_TABLE(ReducedCollisions, "AOD", "REDUCEDCOLLISION", //! reduced collision table (same structure as the original collision table) - o2::soa::Index<>, collision::BCId, +DECLARE_SOA_TABLE(ReducedCollisions, "AOD", "REDCOLLISION", //! reduced collision table (same structure as the original collision table) + o2::soa::Index<>, collision::PosX, collision::PosY, collision::PosZ, collision::CovXX, collision::CovXY, collision::CovYY, collision::CovXZ, collision::CovYZ, collision::CovZZ, collision::Flags, collision::Chi2, collision::NumContrib, collision::CollisionTime, collision::CollisionTimeRes, bc::RunNumber); -DECLARE_SOA_TABLE(ReducedPVMults, "AOD", "REDUCEDPVMULT", //! Multiplicity from the PV contributors, joinable with reducedCollisions +DECLARE_SOA_TABLE(ReducedPVMults, "AOD", "REDPVMULT", //! Multiplicity from the PV contributors, joinable with reducedCollisions mult::MultNTracksPV); -DECLARE_SOA_TABLE(ReducedCentFT0Cs, "AOD", "REDUCEDCENTFT0C", //! Reduced Run 3 FT0C centrality table, joinable with reducedCollisions +DECLARE_SOA_TABLE(ReducedCentFT0Cs, "AOD", "REDCENTFT0C", //! Reduced Run 3 FT0C centrality table, joinable with reducedCollisions cent::CentFT0C); namespace reducedtracks3body @@ -189,7 +189,7 @@ DECLARE_SOA_COLUMN(TOFNSigmaDe, tofNSigmaDe, float); //! Nsigma separation with } // namespace reducedtracks3body -DECLARE_SOA_TABLE_FULL(StoredReducedTracksIU, "ReducedTracks_IU", "AOD", "REDUCEDTRACK_IU", //! On disk version of the track parameters at inner most update (e.g. ITS) as it comes from the tracking +DECLARE_SOA_TABLE_FULL(StoredReducedTracksIU, "ReducedTracks_IU", "AOD", "REDTRACK_IU", //! On disk version of the track parameters at inner most update (e.g. ITS) as it comes from the tracking o2::soa::Index<>, reducedtracks3body::CollisionId, reducedtracks3body::X, reducedtracks3body::Alpha, reducedtracks3body::Y, reducedtracks3body::Z, reducedtracks3body::Snp, reducedtracks3body::Tgl, @@ -224,7 +224,7 @@ DECLARE_SOA_TABLE_FULL(StoredReducedTracksIU, "ReducedTracks_IU", "AOD", "REDUCE reducedtracks3body::ITSClsSizeInLayer, reducedtracks3body::TPCCrossedRowsOverFindableCls); -DECLARE_SOA_EXTENDED_TABLE(ReducedTracksIU, StoredReducedTracksIU, "EXREDUCEDTRACK_IU", 0, //! Track parameters at inner most update (e.g. ITS) as it comes from the tracking +DECLARE_SOA_EXTENDED_TABLE(ReducedTracksIU, StoredReducedTracksIU, "EXREDTRACK_IU", 0, //! Track parameters at inner most update (e.g. ITS) as it comes from the tracking reducedtracks3body::Pt, reducedtracks3body::P, reducedtracks3body::Eta, @@ -258,7 +258,7 @@ DECLARE_SOA_INDEX_COLUMN_FULL(Track2, track2, int, ReducedTracksIU, "_2"); DECLARE_SOA_INDEX_COLUMN_FULL(Collision, collision, int, ReducedCollisions, ""); //! Collision index } // namespace reduceddecay3body -DECLARE_SOA_TABLE(ReducedDecay3Bodys, "AOD", "REDUCEDDECAY3BODY", //! reduced 3-body decay table +DECLARE_SOA_TABLE(ReducedDecay3Bodys, "AOD", "REDDECAY3BODY", //! reduced 3-body decay table o2::soa::Index<>, reduceddecay3body::CollisionId, reduceddecay3body::Track0Id, reduceddecay3body::Track1Id, reduceddecay3body::Track2Id); using ReducedDecay3BodysLinked = soa::Join; diff --git a/PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx b/PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx index 53a44016ddb..09f2b7ba530 100644 --- a/PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx +++ b/PWGLF/TableProducer/Nuspex/decay3bodybuilder.cxx @@ -35,8 +35,6 @@ #include "PWGLF/DataModel/Reduced3BodyTables.h" #include "PWGLF/DataModel/Vtx3BodyTables.h" #include "PWGLF/DataModel/pidTOFGeneric.h" -#include "Common/Core/TrackSelection.h" -#include "Common/DataModel/TrackSelectionTables.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/PIDResponse.h" #include "Common/Core/PID/PIDTOF.h" @@ -80,6 +78,9 @@ using TrackExtPIDIUwithEvTimes = soa::Join; +using ReducedCollisionsMults = soa::Join; +using ReducedCollisionsMultsCents = soa::Join; + struct vtxCandidate { int track0Id; int track1Id; @@ -170,8 +171,6 @@ struct decay3bodyBuilder { std::vector kfPionDaughters; std::vector kfDeuteronDaughters; std::vector kf3bodyVtxCandidates; - std::vector kfCandidates; - // kfCandidate candidate; // Configurables Configurable d_UseAbsDCA{"d_UseAbsDCA", true, "Use Abs DCAs"}; @@ -325,10 +324,8 @@ struct decay3bodyBuilder { // KF event mixing using BinningTypeKF = ColumnBinningPolicy; - BinningTypeKF binningOnPosAndMult{{kfparticleConfigurations.binsVtxZ, kfparticleConfigurations.binsMultiplicity}, true}; // Filters and slices - // Filter collisionFilter = (aod::evsel::sel8 == true && nabs(aod::collision::posZ) < 10.f); Preslice perCollision = o2::aod::decay3body::collisionId; Preslice perReducedCollision = o2::aod::reduceddecay3body::collisionId; @@ -418,7 +415,7 @@ struct decay3bodyBuilder { registry.add("hBachelorTOFNSigmaDe", "", HistType::kTH2F, {{40, -10.0f, 10.0f, "p/z (GeV/c)"}, {40, -10.0f, 10.0f, "TOF n#sigma"}}); } - if (doprocessRun3withKFParticle == true || doprocessRun3withKFParticleEM == true || doprocessRun3withKFParticleStrangenessTracking == true || doprocessRun3withKFParticleTrackRotation == true) { + if (doprocessRun3withKFParticle == true || doprocessRun3withKFParticleEM == true || doprocessRun3withKFParticleStrangenessTracking == true || doprocessRun3withKFParticleTrackRotation == true || doprocessRun3withKFParticleReduced == true || doprocessRun3withKFParticleReducedEM == true) { auto hEventCounterZorro = registry.add("Counters/hEventCounterZorro", "hEventCounterZorro", HistType::kTH1D, {{2, -0.5, 1.5}}); hEventCounterZorro->GetXaxis()->SetBinLabel(1, "Zorro before evsel"); hEventCounterZorro->GetXaxis()->SetBinLabel(2, "Zorro after evsel"); @@ -486,7 +483,7 @@ struct decay3bodyBuilder { registry.add("QA/Event/hVtxCovYZ", "hVtxCovYZ", HistType::kTH1F, {{200, -0.0001f, 0.0001f, "PV cov(YZ) (cm^{2})"}}); } - if (doprocessRun3withKFParticleEM == true) { + if (doprocessRun3withKFParticleEM == true || doprocessRun3withKFParticleReducedEM == true) { auto hPairCounterMixing = registry.add("QA/EM/hPairCounterMixing", "hPairCounterMixing", HistType::kTH1F, {{3, 0.0f, 3.0f}}); hPairCounterMixing->GetXaxis()->SetBinLabel(1, "total"); hPairCounterMixing->GetXaxis()->SetBinLabel(2, "sel8"); @@ -495,8 +492,15 @@ struct decay3bodyBuilder { auto hCombinationCounterMixing = registry.add("QA/EM/hCombinationCounterMixing", "hCombinationCounterMixing", HistType::kTH1F, {{3, 0.0f, 3.0f}}); hCombinationCounterMixing->GetXaxis()->SetBinLabel(1, "total"); hCombinationCounterMixing->GetXaxis()->SetBinLabel(2, "bach sign/ID"); - hCombinationCounterMixing->GetXaxis()->SetBinLabel(3, "bach pT"); + hCombinationCounterMixing->GetXaxis()->SetBinLabel(3, "bach pT or radius, phi"); hCombinationCounterMixing->LabelsOption("v"); + + registry.add("QA/EM/hRadius1", "hRadius1", HistType::kTH1F, {{200, 0.0f, 20.0f, "Radius (cm)"}}); + registry.add("QA/EM/hRadius2", "hRadius2", HistType::kTH1F, {{200, 0.0f, 20.0f, "Radius (cm)"}}); + registry.add("QA/EM/hPhi1", "hPhi1", HistType::kTH1F, {{360, 0.0f, 360.0f, "#phi (degree)"}}); + registry.add("QA/EM/hPhi2", "hPhi2", HistType::kTH1F, {{360, 0.0f, 360.0f, "#phi (degree)"}}); + registry.add("QA/EM/hDeltaRadius", "hDeltaRadius", HistType::kTH1F, {{200, 0.0f, 10.0f, "#Delta Radius (cm)"}}); + registry.add("QA/EM/hDeltaPhi", "hDeltaPhi", HistType::kTH1F, {{360, 0.0f, 360.0f, "#Delta #phi (degree)"}}); } } @@ -620,6 +624,30 @@ struct decay3bodyBuilder { bachelorTOFPID.SetParams(mRespParamsV2); } + void initCCDBReduced(const int& runNumber) + { + o2::parameters::GRPMagField* grpmag = ccdb->getForRun(grpmagPath, runNumber); + if (!grpmag) { + LOG(fatal) << "Got nullptr from CCDB for path " << grpmagPath << " of object GRPMagField and " << grpPath << " of object GRPObject for run number " << runNumber; + } + o2::base::Propagator::initFieldFromGRP(grpmag); + // Fetch magnetic field from ccdb for current collision + d_bz = o2::base::Propagator::Instance()->getNominalBz(); + LOG(info) << "Retrieved GRP for run number " << runNumber << " with magnetic field of " << d_bz << " kZG"; + + // Set magnetic field for KF vertexing +#ifdef HomogeneousField + KFParticle::SetField(d_bz); +#endif + // Set field for DCAfitter + fitter3body.setBz(d_bz); + + if (useMatCorrType == 2) { + // setMatLUT only after magfield has been initalized + o2::base::Propagator::Instance()->setMatLUT(lut); + } + } + //------------------------------------------------------------------ //-------------------- DCA fitter reconstruction ------------------- //------------------------------------------------------------------ @@ -792,6 +820,22 @@ struct decay3bodyBuilder { return true; } + template + double getTOFnSigma(TCollision const& collision, TTrack const& track, bool isEventMixing) + { + // TOF PID of deuteron (set motherhyp correctly) + double tofNSigmaDeuteron = -999; + if (track.has_collision() && track.hasTOF()) { + if (isEventMixing) { + tofNSigmaDeuteron = bachelorTOFPID.GetTOFNSigma(track, collision, collision); + } else { + auto originalcol = track.template collision_as(); + tofNSigmaDeuteron = bachelorTOFPID.GetTOFNSigma(track, originalcol, collision); + } + } + return tofNSigmaDeuteron; + } + //------------------------------------------------------------------ // function to fill candidate table template @@ -916,7 +960,7 @@ struct decay3bodyBuilder { //------------------------------------------------------------------ // 3body candidate builder with KFParticle template - void buildVtx3BodyDataTableKFParticle(TCollision const& collision, TTrack const& trackPos, TTrack const& trackNeg, TTrack const& trackBach, int64_t decay3bodyID, int bachelorcharge, int nRotations) + void buildVtx3BodyDataTableKFParticle(TCollision const& collision, TTrack const& trackPos, TTrack const& trackNeg, TTrack const& trackBach, int64_t decay3bodyID, int bachelorcharge, int nRotations, double tofNSigmaDeuteron) { gROOT->SetBatch(true); gRandom->SetSeed(42); @@ -1013,25 +1057,13 @@ struct decay3bodyBuilder { registry.fill(HIST("Counters/hVtx3BodyCounterKFParticle"), kKfVtxEta); // number of TPC clusters - int tpcNClProton; - int tpcNClPion; - int tpcNClDeuteron = trackBach.tpcNClsFound(); - float tpcChi2NCl = trackBach.tpcChi2NCl(); if (trackBach.tpcNClsFound() <= kfparticleConfigurations.mintpcNClsBach) { return; } - if (isMatter) { - tpcNClPion = trackNeg.tpcNClsFound(); - tpcNClProton = trackPos.tpcNClsFound(); - if ((kfparticleConfigurations.useTPCforPion && trackNeg.tpcNClsFound() <= kfparticleConfigurations.mintpcNClsPion) || trackPos.tpcNClsFound() <= kfparticleConfigurations.mintpcNClsProton) { - return; - } - } else if (!isMatter) { - tpcNClPion = trackPos.tpcNClsFound(); - tpcNClProton = trackNeg.tpcNClsFound(); - if ((kfparticleConfigurations.useTPCforPion && trackPos.tpcNClsFound() <= kfparticleConfigurations.mintpcNClsPion) || trackNeg.tpcNClsFound() <= kfparticleConfigurations.mintpcNClsProton) { - return; - } + if (isMatter && ((kfparticleConfigurations.useTPCforPion && trackNeg.tpcNClsFound() <= kfparticleConfigurations.mintpcNClsPion) || trackPos.tpcNClsFound() <= kfparticleConfigurations.mintpcNClsProton)) { + return; + } else if (!isMatter && ((kfparticleConfigurations.useTPCforPion && trackPos.tpcNClsFound() <= kfparticleConfigurations.mintpcNClsPion) || trackNeg.tpcNClsFound() <= kfparticleConfigurations.mintpcNClsProton)) { + return; } registry.fill(HIST("Counters/hVtx3BodyCounterKFParticle"), kKfVtxTPCNcls); @@ -1074,17 +1106,6 @@ struct decay3bodyBuilder { registry.fill(HIST("Counters/hVtx3BodyCounterKFParticle"), kKfVtxTPCPID); LOG(debug) << "Basic track selections done."; - // TOF PID of deuteron (set motherhyp correctly) - double tofNSigmaDeuteron = -999; - if (trackBach.has_collision() && trackBach.hasTOF()) { - if (isEventMixing) { - tofNSigmaDeuteron = bachelorTOFPID.GetTOFNSigma(trackBach, collision, collision); - } else { - auto originalcol = trackBach.template collision_as(); - tofNSigmaDeuteron = bachelorTOFPID.GetTOFNSigma(trackBach, originalcol, collision); - } - } - // Average ITS cluster size of deuteron track double averageClusterSizeDeuteron(0); int nCls(0); @@ -1209,13 +1230,6 @@ struct decay3bodyBuilder { kfpDeuteron = createKFParticleFromTrackParCov(trackParCovBach, trackBach.sign() * bachelorcharge, constants::physics::MassDeuteron); LOG(debug) << "KFParticle created from deuteron track."; float kfpvPos[3] = {kfpv.GetX(), kfpv.GetY(), kfpv.GetZ()}; - // Check phi angle of "default" candidate daughter tracks - float phiDeuteron_noRotation; - float phiDeuteronSigma_noRotation; - float phiProton_noRotation; - float phiProtonSigma_noRotation; - kfpDeuteron.GetPhi(phiDeuteron_noRotation, phiDeuteronSigma_noRotation); - kfpProton.GetPhi(phiProton_noRotation, phiProtonSigma_noRotation); // Construct vertex /// BRIEF: Case 1: fully uncorrelated bkg via rotation of proton track @@ -1361,14 +1375,6 @@ struct decay3bodyBuilder { trackedClSize = !fTrackedClSizeVector.empty() ? fTrackedClSizeVector[decay3bodyID] : 0; } - // get deuteron and proton phi after rotation - float phiDeuteron_afterRotation; - float phiDeuteronSigma_afterRotation; - float phiProton_afterRotation; - float phiProtonSigma_afterRotation; - kfDeuteronDaughters[i].GetPhi(phiDeuteron_afterRotation, phiDeuteronSigma_afterRotation); - kfProtonDaughters[i].GetPhi(phiProton_afterRotation, phiProtonSigma_afterRotation); - // candidate filling kfCandidate candidate; candidate.collisionID = collision.globalIndex(); @@ -1406,8 +1412,6 @@ struct decay3bodyBuilder { candidate.chi2topoNDF = chi2topoNDF; candidate.ctau = KFHtPV.GetLifeTime(); candidate.trackedClSize = trackedClSize; - candidate.DeltaPhiRotDeuteron = phiDeuteron_noRotation - phiDeuteron_afterRotation; - candidate.DeltaPhiRotProton = phiProton_noRotation - phiProton_afterRotation; // V0 candidate.massV0 = massV0; candidate.chi2massV0 = chi2massV0; @@ -1425,11 +1429,6 @@ struct decay3bodyBuilder { candidate.tpcInnerParam[0] = tpcInnerParamProton; candidate.tpcInnerParam[1] = tpcInnerParamPion; candidate.tpcInnerParam[2] = tpcInnerParamDeuteron; - // daughter track quality - candidate.tpcNClDaughters[0] = tpcNClProton; - candidate.tpcNClDaughters[1] = tpcNClPion; - candidate.tpcNClDaughters[2] = tpcNClDeuteron; - candidate.tpcChi2NClDeuteron = tpcChi2NCl; // daughter DCAs with KF candidate.DCAdaughterToPV[0] = kfProtonDaughters[i].GetDistanceFromVertex(kfpv); candidate.DCAdaughterToPV[1] = kfPionDaughters[i].GetDistanceFromVertex(kfpv); @@ -1518,13 +1517,21 @@ struct decay3bodyBuilder { registry.fill(HIST("hEventCounter"), 0.5, collisions.size()); + int lastRunNumber = -1; + for (const auto& d3body : decay3bodys) { auto t0 = d3body.template track0_as(); auto t1 = d3body.template track1_as(); auto t2 = d3body.template track2_as(); auto collision = d3body.template collision_as(); - auto bc = collision.bc_as(); - initCCDB(bc); + // auto bc = collision.bc_as(); + // initCCDB(bc); + // set magnetic field only when run number changes + if (collision.runNumber() != lastRunNumber) { + initCCDBReduced(collision.runNumber()); + lastRunNumber = collision.runNumber(); // Update the last run number + LOG(debug) << "CCDB initialized for run " << lastRunNumber; + } fillVtxCand(collision, t0, t1, t2, d3body.globalIndex(), bachelorcharge, t2.tofNSigmaDe()); } @@ -1691,9 +1698,9 @@ struct decay3bodyBuilder { auto trackPos = vtx3body.template track0_as(); auto trackNeg = vtx3body.template track1_as(); auto trackBach = vtx3body.template track2_as(); - buildVtx3BodyDataTableKFParticle(collision, trackPos, trackNeg, trackBach, vtx3body.globalIndex(), bachelorcharge, 1); + buildVtx3BodyDataTableKFParticle(collision, trackPos, trackNeg, trackBach, vtx3body.globalIndex(), bachelorcharge, 1 /*nRotations*/, getTOFnSigma(collision, trackBach, false /*isEventMixing*/)); + LOG(debug) << "End of processKFParticle."; } - LOG(debug) << "End of processKFParticle."; } } PROCESS_SWITCH(decay3bodyBuilder, processRun3withKFParticle, "Produce KFParticle decay3body tables", false); @@ -1711,9 +1718,8 @@ struct decay3bodyBuilder { void processRun3withKFParticleEM(ColwithEvTimesMults const& collisions, TrackExtPIDIUwithEvTimes const& tracksIU, aod::Decay3Bodys const& decay3bodys, aod::BCsWithTimestamps const&) { - kf3bodyVtxCandidates.clear(); - auto tuple = std::make_tuple(decay3bodys, tracksIU); + BinningTypeKF binningOnPosAndMult{{kfparticleConfigurations.binsVtxZ, kfparticleConfigurations.binsMultiplicity}, true}; Pair pair{binningOnPosAndMult, kfparticleConfigurations.nEvtMixing, -1, collisions, tuple, &cache}; // indicates that under/overflow (-1) to be ignored for (auto& [c1, decays3body, c2, tracks] : pair) { @@ -1748,7 +1754,7 @@ struct decay3bodyBuilder { continue; } registry.fill(HIST("QA/EM/hCombinationCounterMixing"), 2.5); - buildVtx3BodyDataTableKFParticle(c1, trackPos, trackNeg, track, -1, bachelorcharge, 1); + buildVtx3BodyDataTableKFParticle(c1, trackPos, trackNeg, track, -1 /*vtx3bodyID*/, bachelorcharge, 1 /*nRotations*/, getTOFnSigma(c1, trackBach, true /*isEventMixing*/)); LOG(debug) << "buildVtx3BodyDataTableKFParticle called."; } } @@ -1757,8 +1763,6 @@ struct decay3bodyBuilder { void processRun3withKFParticleTrackRotation(ColwithEvTimes const& collisions, TrackExtPIDIUwithEvTimes const&, aod::Decay3Bodys const& decay3bodys, aod::BCsWithTimestamps const&) { - kf3bodyVtxCandidates.clear(); - // same process as processRun3withKFParticle but with nRotations set according to configurable for (const auto& collision : collisions) { // event selection @@ -1784,12 +1788,131 @@ struct decay3bodyBuilder { auto trackPos = vtx3body.template track0_as(); auto trackNeg = vtx3body.template track1_as(); auto trackBach = vtx3body.template track2_as(); - buildVtx3BodyDataTableKFParticle(collision, trackPos, trackNeg, trackBach, vtx3body.globalIndex(), bachelorcharge, kfparticleConfigurations.nRotations); + buildVtx3BodyDataTableKFParticle(collision, trackPos, trackNeg, trackBach, vtx3body.globalIndex(), bachelorcharge, kfparticleConfigurations.nRotations, getTOFnSigma(collision, trackBach, false /*isEventMixing*/)); } LOG(debug) << "End of processKFParticle."; } } PROCESS_SWITCH(decay3bodyBuilder, processRun3withKFParticleTrackRotation, "Produce KFParticle decay3body tables with rotated tracks", false); + + void processRun3withKFParticleReduced(aod::ReducedCollisions const& collisions, aod::ReducedTracksIU const&, aod::ReducedDecay3Bodys const& decay3bodys) + { + int lastRunNumber = -1; + + for (const auto& collision : collisions) { + // set magnetic field only when run number changes + if (collision.runNumber() != lastRunNumber) { + initCCDBReduced(collision.runNumber()); + lastRunNumber = collision.runNumber(); // Update the last run number + LOG(debug) << "CCDB initialized for run " << lastRunNumber; + } + + // event selection + registry.fill(HIST("Counters/hEventCounterKFParticle"), 2.5); + registry.fill(HIST("QA/Event/hAllSelEventsVtxZ"), collision.posZ()); + + // slice Decay3Body table by collision + const uint64_t collIdx = collision.globalIndex(); + auto Decay3BodyTable_thisCollision = decay3bodys.sliceBy(perReducedCollision, collIdx); + for (auto& vtx3body : Decay3BodyTable_thisCollision) { + auto trackPos = vtx3body.template track0_as(); + auto trackNeg = vtx3body.template track1_as(); + auto trackBach = vtx3body.template track2_as(); + buildVtx3BodyDataTableKFParticle(collision, trackPos, trackNeg, trackBach, vtx3body.globalIndex(), bachelorcharge, 1 /*nRotations*/, trackBach.tofNSigmaDe()); + } + LOG(debug) << "End of processKFParticleDerived."; + } + } + PROCESS_SWITCH(decay3bodyBuilder, processRun3withKFParticleReduced, "Produce KFParticle decay3body tables from derived decay3body data", false); + + void processRun3withKFParticleReducedEM(ReducedCollisionsMults const& collisions, aod::ReducedTracksIU const&, aod::ReducedDecay3Bodys const& decay3bodys) + { + auto tuple = std::make_tuple(decay3bodys); + BinningTypeKF binningOnPosAndMult{{kfparticleConfigurations.binsVtxZ, kfparticleConfigurations.binsMultiplicity}, true}; + SameKindPair pair{binningOnPosAndMult, kfparticleConfigurations.nEvtMixing, -1, collisions, tuple, &cache}; // indicates that under/overflow (-1) to be ignored + + int lastRunNumber = -1; + + for (auto& [c1, decays3body1, c2, decays3body2] : pair) { + registry.fill(HIST("QA/EM/hPairCounterMixing"), 0.5); + + // event selection already applied in reducer task + + // set magnetic field only when run number changes + if (c1.runNumber() != lastRunNumber) { + initCCDBReduced(c1.runNumber()); + lastRunNumber = c1.runNumber(); // Update the last run number + LOG(debug) << "CCDB initialized for run " << lastRunNumber; + } + + for (auto& [decay3body1, decay3body2] : soa::combinations(soa::CombinationsFullIndexPolicy(decays3body1, decays3body2))) { + auto trackPos1 = decay3body1.template track0_as(); + auto trackNeg1 = decay3body1.template track1_as(); + auto trackBach1 = decay3body1.template track2_as(); + auto trackPos2 = decay3body2.template track0_as(); + auto trackNeg2 = decay3body2.template track1_as(); + auto trackBach2 = decay3body2.template track2_as(); + + registry.fill(HIST("QA/EM/hCombinationCounterMixing"), 0.5); + + // ---------- selections bachelor track ---------- + if ((trackBach1.sign() > 0 && !(trackBach2.sign() > 0)) || (trackBach1.sign() < 0 && !(trackBach2.sign() < 0)) || trackBach1.globalIndex() == trackBach2.globalIndex()) { // only combine if trackBach2 has correct sign and is not same as trackBach1 + continue; + } + registry.fill(HIST("QA/EM/hCombinationCounterMixing"), 1.5); + + // ---------- check radius and phi of decay3bodys ---------- + auto trackParCovPos1 = getTrackParCov(trackPos1); + auto trackParCovNeg1 = getTrackParCov(trackNeg1); + auto trackParCovBach1 = getTrackParCov(trackBach1); + auto trackParCovPos2 = getTrackParCov(trackPos2); + auto trackParCovNeg2 = getTrackParCov(trackNeg2); + auto trackParCovBach2 = getTrackParCov(trackBach2); + // create KFParticle objects from tracks + KFParticle kfpProton1, kfpPion1, kfpDeuteron1; + if (trackBach1.sign() > 0) { + kfpProton1 = createKFParticleFromTrackParCov(trackParCovPos1, trackPos1.sign(), constants::physics::MassProton); + kfpPion1 = createKFParticleFromTrackParCov(trackParCovNeg1, trackNeg1.sign(), constants::physics::MassPionCharged); + } else if (!(trackBach1.sign() > 0)) { + kfpProton1 = createKFParticleFromTrackParCov(trackParCovNeg1, trackNeg1.sign(), constants::physics::MassProton); + kfpPion1 = createKFParticleFromTrackParCov(trackParCovPos1, trackPos1.sign(), constants::physics::MassPionCharged); + } + kfpDeuteron1 = createKFParticleFromTrackParCov(trackParCovBach1, trackBach1.sign() * bachelorcharge, constants::physics::MassDeuteron); + KFParticle kfpProton2, kfpPion2, kfpDeuteron2; + if (trackBach2.sign() > 0) { + kfpProton2 = createKFParticleFromTrackParCov(trackParCovPos2, trackPos2.sign(), constants::physics::MassProton); + kfpPion2 = createKFParticleFromTrackParCov(trackParCovNeg2, trackNeg2.sign(), constants::physics::MassPionCharged); + } else if (!(trackBach2.sign() > 0)) { + kfpProton2 = createKFParticleFromTrackParCov(trackParCovNeg2, trackNeg2.sign(), constants::physics::MassProton); + kfpPion2 = createKFParticleFromTrackParCov(trackParCovPos2, trackPos2.sign(), constants::physics::MassPionCharged); + } + kfpDeuteron1 = createKFParticleFromTrackParCov(trackParCovBach2, trackBach2.sign() * bachelorcharge, constants::physics::MassDeuteron); + // fit vertices + KFParticle KFHt1, KFHt2; + fit3bodyVertex(kfpProton1, kfpPion1, kfpDeuteron1, KFHt1); + fit3bodyVertex(kfpProton2, kfpPion2, kfpDeuteron2, KFHt2); + + // ---------- select common radius and phi region ---------- + auto radius1 = std::sqrt(KFHt1.GetX() * KFHt1.GetX() + KFHt1.GetY() * KFHt1.GetY()); + auto radius2 = std::sqrt(KFHt2.GetX() * KFHt2.GetX() + KFHt2.GetY() * KFHt2.GetY()); + registry.fill(HIST("QA/EM/hRadius1"), radius1); + registry.fill(HIST("QA/EM/hRadius2"), radius2); + registry.fill(HIST("QA/EM/hPhi1"), radius1); + registry.fill(HIST("QA/EM/hPhi2"), radius2); + registry.fill(HIST("QA/EM/hDeltaRadius"), std::abs(radius1 - radius2)); + registry.fill(HIST("QA/EM/hDeltaPhi"), std::abs(KFHt1.GetPhi() - KFHt2.GetPhi()) * (180.0 / TMath::Pi())); + if (std::abs(KFHt1.GetPhi() - KFHt2.GetPhi()) * (180.0 / TMath::Pi()) > 10 || std::abs(radius1 - radius2) > 2) { + continue; + } + registry.fill(HIST("QA/EM/hCombinationCounterMixing"), 2.5); + + // ---------- do candidate analysis ---------- + buildVtx3BodyDataTableKFParticle(c1, trackPos1, trackNeg1, trackBach2, -1 /*vtx3bodyID*/, bachelorcharge, 1 /*nRotations*/, trackBach2.tofNSigmaDe()); + LOG(debug) << "buildVtx3BodyDataTableKFParticle called."; + } + } // end pairing loop + } + PROCESS_SWITCH(decay3bodyBuilder, processRun3withKFParticleReducedEM, "Produce KFParticle event mixing decay3body tables from derived decay3body data", false); }; // build link from decay3body -> vtx3body @@ -1833,27 +1956,34 @@ struct kfdecay3bodyDataLinkBuilder { void init(InitContext const&) {} - void processDoNotBuildLink(aod::Collisions::iterator const&) - { - // dummy process function - } - PROCESS_SWITCH(kfdecay3bodyDataLinkBuilder, processDoNotBuildLink, "Do not build data link table.", false); - - // build Decay3Body -> KFDecay3BodyData link table - void processBuildLink(aod::Decay3Bodys const& decay3bodytable, aod::KFVtx3BodyDatas const& vtxdatatable) + template + void buildDataLink(TDecay3Bodys const& decay3bodytable, TVtx3BodyDatas const& vtxdatatable) { std::vector lIndices; lIndices.reserve(decay3bodytable.size()); for (int ii = 0; ii < decay3bodytable.size(); ii++) lIndices[ii] = -1; for (auto& vtxdata : vtxdatatable) { - lIndices[vtxdata.decay3bodyId()] = vtxdata.globalIndex(); + if (vtxdata.decay3bodyId() != -1) { + lIndices[vtxdata.decay3bodyId()] = vtxdata.globalIndex(); + } } for (int ii = 0; ii < decay3bodytable.size(); ii++) { kfvtxdataLink(lIndices[ii]); } } - PROCESS_SWITCH(kfdecay3bodyDataLinkBuilder, processBuildLink, "Build data link table.", true); + + void processStandard(aod::Decay3Bodys const& decay3bodytable, aod::KFVtx3BodyDatas const& vtxdatatable) + { + buildDataLink(decay3bodytable, vtxdatatable); // build Decay3Body -> KFDecay3BodyData link table + } + PROCESS_SWITCH(kfdecay3bodyDataLinkBuilder, processStandard, "Build data link table.", true); + + void processReduced(aod::ReducedDecay3Bodys const& decay3bodytable, aod::KFVtx3BodyDatas const& vtxdatatable) + { + buildDataLink(decay3bodytable, vtxdatatable); // build ReducedDecay3Body -> KFDecay3BodyData link table + } + PROCESS_SWITCH(kfdecay3bodyDataLinkBuilder, processReduced, "Build data link table for reduced data.", true); }; struct decay3bodyLabelBuilder { diff --git a/PWGLF/TableProducer/Nuspex/reduced3bodyCreator.cxx b/PWGLF/TableProducer/Nuspex/reduced3bodyCreator.cxx index e6d38c83bc9..99ab175eca0 100644 --- a/PWGLF/TableProducer/Nuspex/reduced3bodyCreator.cxx +++ b/PWGLF/TableProducer/Nuspex/reduced3bodyCreator.cxx @@ -17,7 +17,6 @@ #include #include #include -#include #include #include "Framework/runDataProcessing.h" @@ -192,6 +191,29 @@ struct reduced3bodyCreator { bachelorTOFPID.SetParams(mRespParamsV2); } + template + void fillTrackTable(TTrack const& daughter, double tofNSigmaTrack, auto collisionIndex) + { + reducedFullTracksPIDIU( + // TrackIU + collisionIndex, + daughter.x(), daughter.alpha(), + daughter.y(), daughter.z(), daughter.snp(), daughter.tgl(), + daughter.signed1Pt(), + // TracksCovIU + daughter.sigmaY(), daughter.sigmaZ(), daughter.sigmaSnp(), daughter.sigmaTgl(), daughter.sigma1Pt(), + daughter.rhoZY(), daughter.rhoSnpY(), daughter.rhoSnpZ(), daughter.rhoTglY(), daughter.rhoTglZ(), + daughter.rhoTglSnp(), daughter.rho1PtY(), daughter.rho1PtZ(), daughter.rho1PtSnp(), daughter.rho1PtTgl(), + // TracksExtra + daughter.tpcInnerParam(), daughter.flags(), daughter.itsClusterSizes(), + daughter.tpcNClsFindable(), daughter.tpcNClsFindableMinusFound(), daughter.tpcNClsFindableMinusCrossedRows(), + daughter.trdPattern(), daughter.tpcChi2NCl(), daughter.tofChi2(), + daughter.tpcSignal(), daughter.tofExpMom(), + // PID + daughter.tpcNSigmaPr(), daughter.tpcNSigmaPi(), daughter.tpcNSigmaDe(), + tofNSigmaTrack); + } + void process(ColwithEvTimesMultsCents const& collisions, TrackExtPIDIUwithEvTimes const&, aod::Decay3Bodys const& decay3bodys, aod::BCsWithTimestamps const&) { @@ -246,7 +268,6 @@ struct reduced3bodyCreator { if (collision.globalIndex() != lastCollisionID) { int runNumber = bc.runNumber(); reducedCollisions( - collision.bcId(), collision.posX(), collision.posY(), collision.posZ(), collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ(), collision.flags(), collision.chi2(), collision.numContrib(), @@ -258,6 +279,9 @@ struct reduced3bodyCreator { lastCollisionID = collision.globalIndex(); } + // Precompute collision index + const auto collisionIndex = reducedCollisions.lastIndex(); + // Save daughter tracks const auto daughter0 = d3body.template track0_as(); const auto daughter1 = d3body.template track1_as(); @@ -270,34 +294,13 @@ struct reduced3bodyCreator { // ---------------------------------------------- // save reduced track table with decay3body daughters - daughterTracks.push_back(daughter0); - daughterTracks.push_back(daughter1); - daughterTracks.push_back(daughter2); - for (int i = 0; i < 3; i++) { - double tofNSigmaTrack = (i == 2) ? tofNSigmaBach : -999.; - reducedFullTracksPIDIU( - // TrackIU - // reducedTrackID + i, - reducedCollisions.lastIndex(), - daughterTracks[i].x(), daughterTracks[i].alpha(), - daughterTracks[i].y(), daughterTracks[i].z(), daughterTracks[i].snp(), daughterTracks[i].tgl(), - daughterTracks[i].signed1Pt(), - // TracksCovIU - daughterTracks[i].sigmaY(), daughterTracks[i].sigmaZ(), daughterTracks[i].sigmaSnp(), daughterTracks[i].sigmaTgl(), daughterTracks[i].sigma1Pt(), - daughterTracks[i].rhoZY(), daughterTracks[i].rhoSnpY(), daughterTracks[i].rhoSnpZ(), daughterTracks[i].rhoTglY(), daughterTracks[i].rhoTglZ(), - daughterTracks[i].rhoTglSnp(), daughterTracks[i].rho1PtY(), daughterTracks[i].rho1PtZ(), daughterTracks[i].rho1PtSnp(), daughterTracks[i].rho1PtTgl(), - // TracksExtra - daughterTracks[i].tpcInnerParam(), daughterTracks[i].flags(), daughterTracks[i].itsClusterSizes(), - daughterTracks[i].tpcNClsFindable(), daughterTracks[i].tpcNClsFindableMinusFound(), daughterTracks[i].tpcNClsFindableMinusCrossedRows(), - daughterTracks[i].trdPattern(), daughterTracks[i].tpcChi2NCl(), daughterTracks[i].tofChi2(), - daughterTracks[i].tpcSignal(), daughterTracks[i].tofExpMom(), - // PID - daughterTracks[i].tpcNSigmaPr(), daughterTracks[i].tpcNSigmaPi(), daughterTracks[i].tpcNSigmaDe(), - tofNSigmaTrack); - } + fillTrackTable(daughter0, -999, collisionIndex); + fillTrackTable(daughter1, -999, collisionIndex); + fillTrackTable(daughter2, tofNSigmaBach, collisionIndex); // save reduced decay3body table - reducedDecay3Bodys(reducedCollisions.lastIndex(), reducedFullTracksPIDIU.lastIndex() - 2, reducedFullTracksPIDIU.lastIndex() - 1, reducedFullTracksPIDIU.lastIndex()); + const auto trackStartIndex = reducedFullTracksPIDIU.lastIndex(); + reducedDecay3Bodys(collisionIndex, trackStartIndex - 2, trackStartIndex - 1, trackStartIndex); } registry.fill(HIST("hEventCounter"), 3.5, reducedCollisions.lastIndex() + 1);