diff --git a/PWGCF/Femto/Core/femtoUtils.h b/PWGCF/Femto/Core/femtoUtils.h index fa4b9a23466..c623e161f2b 100644 --- a/PWGCF/Femto/Core/femtoUtils.h +++ b/PWGCF/Femto/Core/femtoUtils.h @@ -30,6 +30,7 @@ #include #include #include +#include #include namespace o2::analysis::femto @@ -83,7 +84,7 @@ inline double getPdgMass(int pdgCode) { // use this function instead of TDatabasePDG to return masses defined in the PhysicsConstants.h header // this approach saves a lot of memory and important partilces like deuteron are missing in TDatabasePDG anyway - double mass = 0.f; + double mass = 0.; // add new particles if necessary here switch (std::abs(pdgCode)) { case kPiPlus: @@ -102,10 +103,10 @@ inline double getPdgMass(int pdgCode) mass = o2::constants::physics::MassPhi; break; case kRho770_0: - mass = 775.26; // not defined in O2? + mass = 0.77526; // not defined in O2? break; case kRho770Plus: - mass = 775.11; // not defined in O2? + mass = 0.77511; // not defined in O2? break; case o2::constants::physics::Pdg::kK0Star892: mass = o2::constants::physics::MassK0Star892; @@ -138,7 +139,7 @@ inline double getPdgMass(int pdgCode) mass = o2::constants::physics::MassOmegaMinus; break; default: - LOG(warn) << "PDG code is not suppored. Return 0..."; + LOG(warn) << "PDG code " << pdgCode << " is not suppored. Return 0... "; } return mass; } @@ -154,11 +155,15 @@ concept HasQvectors = requires(T col) { }; template -concept HasEventShape = requires(T col) { - col.qvec(); - col.eventPlaneAngle(); +concept HasEventShapeRow = requires(T row) { + row.qvec(); + row.eventPlaneAngle(); }; +/// accepts either a row/iterator or a table (Filtered> etc.) +template +concept HasEventShape = HasEventShapeRow> || (requires { typename std::decay_t::iterator; } && HasEventShapeRow::iterator>); + /// Recalculate pT for Kinks (Sigmas) using kinematic constraints inline float calcPtnew(float pxMother, float pyMother, float pzMother, float pxDaughter, float pyDaughter, float pzDaughter) { diff --git a/PWGCF/Femto/Core/pairBuilder.h b/PWGCF/Femto/Core/pairBuilder.h index 3998732a004..ea4963b964a 100644 --- a/PWGCF/Femto/Core/pairBuilder.h +++ b/PWGCF/Femto/Core/pairBuilder.h @@ -202,8 +202,8 @@ class PairTrackTrackBuilder } } - template - void processMixedEvent(T1 const& cols, T2& trackTable, T3& partition1, T4& partition2, T5& cache, T6& binsVtxMult, T7& binsVtxCent, T8& binsVtxMultCent) + template + void processMixedEvent(T1 const& cols, T2& trackTable, T3& partition1, T4& partition2, T5& cache, T6& binsVtxMult, T7& binsVtxCent, T8& binsVtxMultCent, T9& binsVtxCentEventPlaneAngle) { if (mSameSpecies) { @@ -217,6 +217,13 @@ class PairTrackTrackBuilder case static_cast(pairhistmanager::kVtxMultCent): pairprocesshelpers::processMixedEvent(cols, partition1, partition1, trackTable, cache, binsVtxMultCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPcMe); break; + case static_cast(pairhistmanager::kVtxCentEventPlaneAngle): + if constexpr (utils::HasEventShape) { + pairprocesshelpers::processMixedEvent(cols, partition1, partition1, trackTable, cache, binsVtxCentEventPlaneAngle, mMixingDepth, mPairHistManagerMe, mCprMe, mPcMe); + } else { + LOG(fatal) << "Mixing policy kVtxCentEventPlaneAngle requires a collision table with event-shape columns. Breaking..."; + } + break; default: LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; } @@ -231,14 +238,21 @@ class PairTrackTrackBuilder case static_cast(pairhistmanager::kVtxMultCent): pairprocesshelpers::processMixedEvent(cols, partition1, partition2, trackTable, cache, binsVtxMultCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPcMe); break; + case static_cast(pairhistmanager::kVtxCentEventPlaneAngle): + if constexpr (utils::HasEventShape) { + pairprocesshelpers::processMixedEvent(cols, partition1, partition2, trackTable, cache, binsVtxCentEventPlaneAngle, mMixingDepth, mPairHistManagerMe, mCprMe, mPcMe); + } else { + LOG(fatal) << "Mixing policy kVtxCentEventPlaneAngle requires a collision table with event-shape columns. Breaking..."; + } + break; default: LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; } } } - template - void processMixedEvent(T1 const& cols, T2 const& mcCols, T3& trackTable, T4& partition1, T5& partition2, T6 const& mcParticles, T7 const& mcMothers, T8 const& mcPartonicMothers, T9& cache, T10& binsVtxMult, T11& binsVtxCent, T12& binsVtxMultCent) + template + void processMixedEvent(T1 const& cols, T2 const& mcCols, T3& trackTable, T4& partition1, T5& partition2, T6 const& mcParticles, T7 const& mcMothers, T8 const& mcPartonicMothers, T9& cache, T10& binsVtxMult, T11& binsVtxCent, T12& binsVtxMultCent, T13& binsVtxCentEventPlaneAngle) { if (mSameSpecies) { switch (mMixingPolicy) { @@ -251,6 +265,13 @@ class PairTrackTrackBuilder case static_cast(pairhistmanager::kVtxMultCent): pairprocesshelpers::processMixedEvent(cols, mcCols, partition1, partition1, trackTable, mcParticles, mcMothers, mcPartonicMothers, cache, binsVtxMultCent, mMixingDepth, mPairHistManagerMe, mTrackCleaner1, mTrackCleaner1, mCprMe, mPcMe); break; + case static_cast(pairhistmanager::kVtxCentEventPlaneAngle): + if constexpr (utils::HasEventShape) { + pairprocesshelpers::processMixedEvent(cols, mcCols, partition1, partition1, trackTable, mcParticles, mcMothers, mcPartonicMothers, cache, binsVtxCentEventPlaneAngle, mMixingDepth, mPairHistManagerMe, mTrackCleaner1, mTrackCleaner1, mCprMe, mPcMe); + } else { + LOG(fatal) << "Mixing policy kVtxCentEventPlaneAngle requires a collision table with event-shape columns. Breaking..."; + } + break; default: LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; } @@ -265,6 +286,13 @@ class PairTrackTrackBuilder case static_cast(pairhistmanager::kVtxMultCent): pairprocesshelpers::processMixedEvent(cols, mcCols, partition1, partition2, trackTable, mcParticles, mcMothers, mcPartonicMothers, cache, binsVtxMultCent, mMixingDepth, mPairHistManagerMe, mTrackCleaner1, mTrackCleaner2, mCprMe, mPcMe); break; + case static_cast(pairhistmanager::kVtxCentEventPlaneAngle): + if constexpr (utils::HasEventShape) { + pairprocesshelpers::processMixedEvent(cols, mcCols, partition1, partition2, trackTable, mcParticles, mcMothers, mcPartonicMothers, cache, binsVtxCentEventPlaneAngle, mMixingDepth, mPairHistManagerMe, mTrackCleaner1, mTrackCleaner2, mCprMe, mPcMe); + } else { + LOG(fatal) << "Mixing policy kVtxCentEventPlaneAngle requires a collision table with event-shape columns. Breaking..."; + } + break; default: LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; } diff --git a/PWGCF/Femto/Core/pairHistManager.h b/PWGCF/Femto/Core/pairHistManager.h index acedd9e59e8..8554c676b15 100644 --- a/PWGCF/Femto/Core/pairHistManager.h +++ b/PWGCF/Femto/Core/pairHistManager.h @@ -156,6 +156,10 @@ enum PairHist { kQside, kQlong, kQoutQsideQlong, + + // event shape enginerring + kQoutQsideQlongEventPlaneAngleQvector, + kPairHistogramLast }; @@ -163,6 +167,7 @@ enum MixingPolicy { kVtxMult, kVtxCent, kVtxMultCent, + kVtxCentEventPlaneAngle, kMixingPolicyLast }; @@ -172,6 +177,7 @@ struct ConfMixing : o2::framework::ConfigurableGroup { o2::framework::ConfigurableAxis multBins{"multBins", {o2::framework::VARIABLE_WIDTH, 0.0f, 4.0f, 8.0f, 12.0f, 16.0f, 20.0f, 24.0f, 28.0f, 32.0f, 36.0f, 40.0f, 44.0f, 48.0f, 52.0f, 56.0f, 60.0f, 64.0f, 68.0f, 72.0f, 76.0f, 80.0f, 84.0f, 88.0f, 92.0f, 96.0f, 100.0f, 200.0f}, "Mixing bins - multiplicity"}; o2::framework::ConfigurableAxis centBins{"centBins", {o2::framework::VARIABLE_WIDTH, 0.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f, 100.0f}, "Mixing bins - centrality"}; o2::framework::ConfigurableAxis vtxBins{"vtxBins", {o2::framework::VARIABLE_WIDTH, -10.0f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - z-vertex"}; + o2::framework::ConfigurableAxis eventPlaneAngle{"eventPlaneAngle", {10, 0.f, 1.f * o2::constants::math::TwoPI}, "Mixing bins - event plane angle"}; o2::framework::Configurable depth{"depth", 5, "Number of events for mixing"}; o2::framework::Configurable policy{"policy", 0, "Binning policy for mixing (alywas in combination with z-vertex) -> 0: multiplicity, -> 1: centrality, -> 2: both"}; o2::framework::Configurable sameSpecies{"sameSpecies", false, "Enable if particle 1 and particle 2 are the same"}; @@ -228,6 +234,9 @@ struct ConfPairBinning : o2::framework::ConfigurableGroup { o2::framework::ConfigurableAxis qout{"qout", {{300, -1.5f, 1.5f}}, "q_{out} (GeV/c) in LCMS"}; o2::framework::ConfigurableAxis qside{"qside", {{300, -1.5f, 1.5f}}, "q_{side} (GeV/c) in LCMS"}; o2::framework::ConfigurableAxis qlong{"qlong", {{300, -1.5f, 1.5f}}, "q_{long} (GeV/c) in LCMS"}; + o2::framework::Configurable plotEventShape{"plotEventShape", false, "(Reco/Mc) Enable 5D (q_out, q_side, q_long, event plane anglke, qvector) histogram"}; + o2::framework::ConfigurableAxis eventPlaneAngle{"eventPlaneAngle", {{10, 0.f, 1.f * o2::constants::math::TwoPI}}, "event plane angle"}; + o2::framework::ConfigurableAxis qvector{"qvector", {{o2::framework::VARIABLE_WIDTH, 0.50f, 68.50f, 100.50f, 126.50f, 151.50f, 176.50f, 203.50f, 232.50f, 269.50f, 322.50f, 833.50f}}, "qvector"}; o2::framework::Configurable plotSH{"plotSH", false, "(Reco) Enable spherical-harmonics decomposition of the pair momentum-difference vector"}; o2::framework::Configurable shLMax{"shLMax", 2, "Maximum l for SH decomposition (0..5). FemtoUniverse hard-codes 1."}; o2::framework::Configurable shFrame{"shFrame", 1, "SH reference frame/variable: 0=LCMS non-identical (k*), 1=LCMS identical (qinv, FemtoUniverse default), 2=PRF (q_PRF, matches FemtoUniverse isIdenPRF=true)"}; @@ -354,6 +363,7 @@ constexpr std::array, kPairHistogramLast> {kTrueQoutVsQout, o2::framework::HistType::kTH2F, "hTrueQoutVsQout", "q_{out,True} vs q_{out}; q_{out,True} (GeV/#it{c}); q_{out} (GeV/#it{c})"}, {kTrueQsideVsQside, o2::framework::HistType::kTH2F, "hTrueQsideVsQside", "q_{side,True} vs q_{side}; q_{side,True} (GeV/#it{c}); q_{side} (GeV/#it{c})"}, {kTrueQlongVsQlong, o2::framework::HistType::kTH2F, "hTrueQlongVsQlong", "q_{long,True} vs q_{long}; q_{long,True} (GeV/#it{c}); q_{long} (GeV/#it{c})"}, + {kQoutQsideQlongEventPlaneAngleQvector, o2::framework::HistType::kTHnSparseF, "hQoutQsideQlongEventPlaneAngleQvector", "Event shape enginering; q_{out} (GeV/#it{c}); q_{side} (GeV/#it{c}); q_{long} (GeV/#it{c}); #varphi_{EP}; q-vector;"}, }}; // NOLINTNEXTLINE(cppcoreguidelines-macro-usage) @@ -403,7 +413,8 @@ constexpr std::array, kPairHistogramLast> {kQout, {(confAnalysis).qout}}, \ {kQside, {(confAnalysis).qside}}, \ {kQlong, {(confAnalysis).qlong}}, \ - {kQoutQsideQlong, {(confAnalysis).qout, (confAnalysis).qside, (confAnalysis).qlong}}, + {kQoutQsideQlong, {(confAnalysis).qout, (confAnalysis).qside, (confAnalysis).qlong}}, \ + {kQoutQsideQlongEventPlaneAngleQvector, {(confAnalysis).qout, (confAnalysis).qside, (confAnalysis).qlong, (confAnalysis).eventPlaneAngle, (confAnalysis).qvector}}, // mixing-qa entries are independent of reco vs mc-truth status — both the reco // analysis path and the pure mc-truth path need them whenever kSe/kMe is set @@ -578,6 +589,7 @@ class PairHistManager mPlotDalitz = ConfPairBinning.plotDalitz.value; mPlotDeltaEtaDeltaPhi = ConfPairBinning.plotDeltaEtaDeltaPhi.value; mPlotBertschPratt = ConfPairBinning.plotBertschPratt.value; + mPlotEventShape = ConfPairBinning.plotEventShape.value; mPlotSH = ConfPairBinning.plotSH.value; mShUseCent = ConfPairBinning.shUseCent.value; @@ -658,81 +670,10 @@ class PairHistManager mAbsCharge2 = std::abs(chargeAbsParticle2); } - template - void setPair(T1 const& particle1, T2 const& particle2, T3 const& trackTable) - { - // if one of the particles has a mass getter (like lambda), we cache the value for the filling later - // otherwise set it to the pdg mass - if constexpr (utils::HasMass) { - mRecoMass1 = particle1.mass(); - } else { - mRecoMass1 = mPdgMass1; - } - if constexpr (utils::HasMass) { - mRecoMass2 = particle2.mass(); - } else { - mRecoMass2 = mPdgMass2; - } - - // get mass for 4-vectors - double mass1 = 0.f; - double mass2 = 0.f; - if (mUsePdgMass) { - mass1 = mPdgMass1; - mass2 = mPdgMass2; - } else { - mass1 = mRecoMass1; - mass2 = mRecoMass2; - } - - // pt in track table is calculated from 1/signedPt from the original track table - // in case of He with Z=2, we have to rescale the pt with the absolute charge - mParticle1 = ROOT::Math::PtEtaPhiMVector(mAbsCharge1 * particle1.pt(), particle1.eta(), particle1.phi(), mass1); - mParticle2 = ROOT::Math::PtEtaPhiMVector(mAbsCharge2 * particle2.pt(), particle2.eta(), particle2.phi(), mass2); - - // set kT - mKt = getKt(mParticle1, mParticle2); - - // set mT - mMt = getMt(mParticle1, mParticle2); - - // set Minv - mMassInv = getMinv(mParticle1, mParticle2); - - // set kstar - mKstar = getKstar(mParticle1, mParticle2); - - if (mPlotBertschPratt) { - std::tie(mQout, mQside, mQlong) = computeBertschPrattLCMS(mParticle1, mParticle2); - } - - if (mPlotSH) { - std::tie(mShKv, mShOut, mShSide, mShLong) = computeShKinematics(mParticle1, mParticle2); - } - - if (mPlotDeltaEtaDeltaPhi) { - mDeltaEta = particle1.eta() - particle2.eta(); - mDeltaPhi = RecoDecay::constrainAngle(particle1.phi() - particle2.phi(), -o2::constants::math::PIHalf); - } - - if (mPlotDalitz) { - if constexpr (modes::isEqual(particleType1, modes::Particle::kTrack) && (modes::isEqual(particleType2, modes::Particle::kV0) || modes::isEqual(particleType2, modes::Particle::kTwoTrackResonance) || modes::isEqual(particleType2, modes::Particle::kCharmHadron)) && - requires(T2 p) { p.posDauId(); p.negDauId(); }) { - auto posDaughter = trackTable.rawIteratorAt(particle2.posDauId() - trackTable.offset()); - auto negDaughter = trackTable.rawIteratorAt(particle2.negDauId() - trackTable.offset()); - ROOT::Math::PtEtaPhiMVector posDau4v = ROOT::Math::PtEtaPhiMVector(posDaughter.pt(), posDaughter.eta(), posDaughter.phi(), mPdgMassPosDau2); - ROOT::Math::PtEtaPhiMVector negDau4v = ROOT::Math::PtEtaPhiMVector(negDaughter.pt(), negDaughter.eta(), negDaughter.phi(), mPdgMassNegDau2); - mMassTot2 = (mParticle1 + posDau4v + negDau4v).M2(); - mMass12 = (mParticle1 + posDau4v).M2(); - mMass13 = (mParticle1 + negDau4v).M2(); - } - } - } - template void setPair(T1 const& particle1, T2 const& particle2, T3 const& trackTable, T4 const& col) { - setPair(particle1, particle2, trackTable); + setPairKinematics(particle1, particle2, trackTable, col); mMult = col.mult(); mCent = col.cent(); } @@ -740,7 +681,7 @@ class PairHistManager template void setPair(T1 const& particle1, T2 const& particle2, T3 const& trackTable, T4 const& col1, T5 const& col2) { - setPair(particle1, particle2, trackTable); + setPairKinematics(particle1, particle2, trackTable, col1); mMult = 0.5f * (col1.mult() + col2.mult()); // if mixing with multiplicity, should be in the same mixing bin mCent = 0.5f * (col1.cent() + col2.cent()); // if mixing with centrality, should be in the same mixing bin } @@ -1047,6 +988,11 @@ class PairHistManager mHistogramRegistry->add(analysisDir + getHistNameV2(kQlong, HistTable), getHistDesc(kQlong, HistTable), getHistType(kQlong, HistTable), {Specs.at(kQlong)}); mHistogramRegistry->add(analysisDir + getHistNameV2(kQoutQsideQlong, HistTable), getHistDesc(kQoutQsideQlong, HistTable), getHistType(kQoutQsideQlong, HistTable), {Specs.at(kQoutQsideQlong)}); } + + if (mPlotEventShape) { + mHistogramRegistry->add(analysisDir + getHistNameV2(kQoutQsideQlongEventPlaneAngleQvector, HistTable), getHistDesc(kQoutQsideQlongEventPlaneAngleQvector, HistTable), getHistType(kQoutQsideQlongEventPlaneAngleQvector, HistTable), {Specs.at(kQoutQsideQlongEventPlaneAngleQvector)}); + } + if (mPlotSH) { const int nJM = (mShLMax + 1) * (mShLMax + 1); const int nCent = static_cast(mShCentEdges.size()) - 1; // n edges -> n - 1 bins @@ -1225,6 +1171,89 @@ class PairHistManager } } + template + void setPairKinematics(T1 const& particle1, T2 const& particle2, T3 const& trackTable, T4 const& col) + { + // if one of the particles has a mass getter (like lambda), we cache the value for the filling later + // otherwise set it to the pdg mass + if constexpr (utils::HasMass) { + mRecoMass1 = particle1.mass(); + } else { + mRecoMass1 = mPdgMass1; + } + if constexpr (utils::HasMass) { + mRecoMass2 = particle2.mass(); + } else { + mRecoMass2 = mPdgMass2; + } + + // get mass for 4-vectors + double mass1 = 0.f; + double mass2 = 0.f; + if (mUsePdgMass) { + mass1 = mPdgMass1; + mass2 = mPdgMass2; + } else { + mass1 = mRecoMass1; + mass2 = mRecoMass2; + } + + // pt in track table is calculated from 1/signedPt from the original track table + // in case of He with Z=2, we have to rescale the pt with the absolute charge + mParticle1 = ROOT::Math::PtEtaPhiMVector(mAbsCharge1 * particle1.pt(), particle1.eta(), particle1.phi(), mass1); + mParticle2 = ROOT::Math::PtEtaPhiMVector(mAbsCharge2 * particle2.pt(), particle2.eta(), particle2.phi(), mass2); + + // set kT + mKt = getKt(mParticle1, mParticle2); + + // set mT + mMt = getMt(mParticle1, mParticle2); + + // set Minv + mMassInv = getMinv(mParticle1, mParticle2); + + // set kstar + mKstar = getKstar(mParticle1, mParticle2); + + if (mPlotBertschPratt) { + std::tie(mQout, mQside, mQlong) = computeBertschPrattLCMS(mParticle1, mParticle2); + } + + if (mPlotEventShape) { + if (!mPlotBertschPratt) { + std::tie(mQout, mQside, mQlong) = computeBertschPrattLCMS(mParticle1, mParticle2); + } + if constexpr (utils::HasEventShape) { + mPairPhiFromEventPlaneAngle = getPairPhiFromEventPlaneAngle(mParticle1, mParticle2, col); + mQvector = col.qvec(); + } else { + LOG(fatal) << "plotEventShape is enabled but the collision table has no event-shape columns. Breaking..."; + } + } + + if (mPlotSH) { + std::tie(mShKv, mShOut, mShSide, mShLong) = computeShKinematics(mParticle1, mParticle2); + } + + if (mPlotDeltaEtaDeltaPhi) { + mDeltaEta = particle1.eta() - particle2.eta(); + mDeltaPhi = RecoDecay::constrainAngle(particle1.phi() - particle2.phi(), -o2::constants::math::PIHalf); + } + + if (mPlotDalitz) { + if constexpr (modes::isEqual(particleType1, modes::Particle::kTrack) && (modes::isEqual(particleType2, modes::Particle::kV0) || modes::isEqual(particleType2, modes::Particle::kTwoTrackResonance) || modes::isEqual(particleType2, modes::Particle::kCharmHadron)) && + requires(T2 p) { p.posDauId(); p.negDauId(); }) { + auto posDaughter = trackTable.rawIteratorAt(particle2.posDauId() - trackTable.offset()); + auto negDaughter = trackTable.rawIteratorAt(particle2.negDauId() - trackTable.offset()); + ROOT::Math::PtEtaPhiMVector posDau4v = ROOT::Math::PtEtaPhiMVector(posDaughter.pt(), posDaughter.eta(), posDaughter.phi(), mPdgMassPosDau2); + ROOT::Math::PtEtaPhiMVector negDau4v = ROOT::Math::PtEtaPhiMVector(negDaughter.pt(), negDaughter.eta(), negDaughter.phi(), mPdgMassNegDau2); + mMassTot2 = (mParticle1 + posDau4v + negDau4v).M2(); + mMass12 = (mParticle1 + posDau4v).M2(); + mMass13 = (mParticle1 + negDau4v).M2(); + } + } + } + void fillAnalysis() { if (mPlot1d) { @@ -1327,6 +1356,9 @@ class PairHistManager mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQlong, HistTable)), mQlong); mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQoutQsideQlong, HistTable)), mQout, mQside, mQlong); } + if (mPlotEventShape) { + mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQoutQsideQlongEventPlaneAngleQvector, HistTable)), mQout, mQside, mQlong, mPairPhiFromEventPlaneAngle, mQvector); + } if (mPlotSH) { const float shCentValue = mShUseCent ? mCent : mMult; const int iCent = findShBin(shCentValue, mShCentEdges); @@ -1611,6 +1643,14 @@ class PairHistManager static_cast(fDKOut), static_cast(fDKSide), static_cast(fDKLong)}; } + template + float getPairPhiFromEventPlaneAngle(const ROOT::Math::PtEtaPhiMVector& part1, const ROOT::Math::PtEtaPhiMVector& part2, const T1& col) + { + const ROOT::Math::PtEtaPhiMVector trackSum = part1 + part2; + float phi = RecoDecay::constrainAngle(trackSum.Phi() - col.eventPlaneAngle()); + return phi; + } + o2::framework::HistogramRegistry* mHistogramRegistry = nullptr; bool mUsePdgMass = true; double mPdgMass1 = 0.; @@ -1702,6 +1742,8 @@ class PairHistManager bool mPlotBertschPratt = false; + bool mPlotEventShape = false; + float mQout = 0.f; float mQside = 0.f; float mQlong = 0.f; @@ -1709,6 +1751,9 @@ class PairHistManager float mDeltaEta = 0.f; float mDeltaPhi = 0.f; + float mQvector = 0.f; + float mPairPhiFromEventPlaneAngle = 0.f; + // Spherical harmonics bool mPlotSH = false; bool mShUseCent = false; diff --git a/PWGCF/Femto/Tasks/femtoPairTrackTrack.cxx b/PWGCF/Femto/Tasks/femtoPairTrackTrack.cxx index 495ed345efd..818db97e7eb 100644 --- a/PWGCF/Femto/Tasks/femtoPairTrackTrack.cxx +++ b/PWGCF/Femto/Tasks/femtoPairTrackTrack.cxx @@ -26,6 +26,7 @@ #include "PWGCF/Femto/Core/trackHistManager.h" #include "PWGCF/Femto/DataModel/FemtoTables.h" +#include #include #include #include @@ -54,6 +55,11 @@ struct FemtoPairTrackTrack { using FilteredFemtoCollisionsWithLabel = o2::soa::Filtered; using FilteredFemtoCollisionWithLabel = FilteredFemtoCollisionsWithLabel::iterator; + // collision table with event shape, same filter as standard collisions + using FemtoCollisionsWithEventShape = o2::soa::Join; + using FilteredFemtoCollisionsWithEventShape = o2::soa::Filtered; + using FilteredFemtoCollisionWithEventShape = FilteredFemtoCollisionsWithEventShape::iterator; + using FemtoTracks = o2::soa::Join; // for analysis which require particles at high pt, add tof mass so sidebands can be used @@ -111,19 +117,22 @@ struct FemtoPairTrackTrack { std::vector defaultVtxBins{10, -10, 10}; std::vector defaultMultBins{50, 0, 200}; std::vector defaultCentBins{10, 0, 100}; + std::vector defaultEventPlaneAngleBins{10, 0, o2::constants::math::TwoPI}; o2::framework::ColumnBinningPolicy mixBinsVtxMult{{defaultVtxBins, defaultMultBins}, true}; o2::framework::ColumnBinningPolicy mixBinsVtxCent{{defaultVtxBins, defaultCentBins}, true}; o2::framework::ColumnBinningPolicy mixBinsVtxMultCent{{defaultVtxBins, defaultMultBins, defaultCentBins}, true}; + o2::framework::ColumnBinningPolicy mixBinsVtxCentEventPlaneAngle{{defaultVtxBins, defaultCentBins, defaultEventPlaneAngleBins}, true}; pairhistmanager::ConfMixing confMixing; o2::framework::HistogramRegistry hRegistry{"FemtoTrackTrack", {}, o2::framework::OutputObjHandlingPolicy::AnalysisObject}; void init(o2::framework::InitContext&) { - if ((static_cast(doprocessSameEvent) + static_cast(doprocessSameEventWithMass) + static_cast(doprocessSameEventMc)) > 1 || (static_cast(doprocessMixedEvent) + static_cast(doprocessMixedEventWithMass) + static_cast(doprocessMixedEventMc)) > 1) { + if ((static_cast(doprocessSameEvent) + static_cast(doprocessSameEventWithMass) + static_cast(doprocessSameEventWithEventShape) + static_cast(doprocessSameEventMc)) > 1 || + (static_cast(doprocessMixedEvent) + static_cast(doprocessMixedEventWithMass) + static_cast(doprocessMixedEventWithEventShape) + static_cast(doprocessMixedEventMc)) > 1) { LOG(fatal) << "More than 1 same or mixed event process function is activated. Breaking..."; } - bool processData = doprocessSameEvent || doprocessMixedEvent || doprocessSameEventWithMass || doprocessMixedEventWithMass; + bool processData = doprocessSameEvent || doprocessMixedEvent || doprocessSameEventWithMass || doprocessMixedEventWithMass || doprocessSameEventWithEventShape || doprocessMixedEventWithEventShape; bool processMc = doprocessSameEventMc || doprocessMixedEventMc; if (processData && processMc) { LOG(fatal) << "Both data and mc processing is activated. Breaking..."; @@ -134,6 +143,7 @@ struct FemtoPairTrackTrack { mixBinsVtxMult = {{confMixing.vtxBins.value, confMixing.multBins.value}, true}; mixBinsVtxCent = {{confMixing.vtxBins.value, confMixing.centBins.value}, true}; mixBinsVtxMultCent = {{confMixing.vtxBins.value, confMixing.multBins.value, confMixing.centBins.value}, true}; + mixBinsVtxCentEventPlaneAngle = {{confMixing.vtxBins.value, confMixing.centBins.value, confMixing.eventPlaneAngle.value}, true}; // setup histogram specs std::map> colHistSpec; @@ -171,6 +181,12 @@ struct FemtoPairTrackTrack { } PROCESS_SWITCH(FemtoPairTrackTrack, processSameEventWithMass, "Enable processing same event processing (with track masses)", false); + void processSameEventWithEventShape(FilteredFemtoCollisionWithEventShape const& col, FemtoTracks const& tracks) + { + pairTrackTrackBuilder.processSameEvent(col, tracks, trackPartition1, trackPartition2, cache); + } + PROCESS_SWITCH(FemtoPairTrackTrack, processSameEventWithEventShape, "Enable processing same event processing with event shape information", false); + void processSameEventMc(FilteredFemtoCollisionWithLabel const& col, o2::aod::FMcCols const& mcCols, FemtoTracksWithLabel const& tracks, FemtoMcParticlesWithLabel const& mcParticles, o2::aod::FMcMothers const& mcMothers, o2::aod::FMcPartMoths const& mcPartonicMothers) { pairTrackTrackBuilder.processSameEvent(col, mcCols, tracks, trackWithLabelPartition1, trackWithLabelPartition2, mcParticles, mcMothers, mcPartonicMothers, cache); @@ -179,19 +195,25 @@ struct FemtoPairTrackTrack { void processMixedEvent(FilteredFemtoCollisions const& cols, FemtoTracks const& tracks) { - pairTrackTrackBuilder.processMixedEvent(cols, tracks, trackPartition1, trackPartition2, cache, mixBinsVtxMult, mixBinsVtxCent, mixBinsVtxMultCent); + pairTrackTrackBuilder.processMixedEvent(cols, tracks, trackPartition1, trackPartition2, cache, mixBinsVtxMult, mixBinsVtxCent, mixBinsVtxMultCent, mixBinsVtxCentEventPlaneAngle); } PROCESS_SWITCH(FemtoPairTrackTrack, processMixedEvent, "Enable processing mixed event processing", true); void processMixedEventWithMass(FilteredFemtoCollisions const& cols, FemtoTracksWithMass const& tracks) { - pairTrackTrackBuilder.processMixedEvent(cols, tracks, trackWithMassPartition1, trackWithMassPartition2, cache, mixBinsVtxMult, mixBinsVtxCent, mixBinsVtxMultCent); + pairTrackTrackBuilder.processMixedEvent(cols, tracks, trackWithMassPartition1, trackWithMassPartition2, cache, mixBinsVtxMult, mixBinsVtxCent, mixBinsVtxMultCent, mixBinsVtxCentEventPlaneAngle); } PROCESS_SWITCH(FemtoPairTrackTrack, processMixedEventWithMass, "Enable processing mixed event processing (with track masses)", false); + void processMixedEventWithEventShape(FilteredFemtoCollisionsWithEventShape const& cols, FemtoTracks const& tracks) + { + pairTrackTrackBuilder.processMixedEvent(cols, tracks, trackPartition1, trackPartition2, cache, mixBinsVtxMult, mixBinsVtxCent, mixBinsVtxMultCent, mixBinsVtxCentEventPlaneAngle); + } + PROCESS_SWITCH(FemtoPairTrackTrack, processMixedEventWithEventShape, "Enable processing mixed event processing with event shape information", false); + void processMixedEventMc(FilteredFemtoCollisionsWithLabel const& cols, o2::aod::FMcCols const& mcCols, FemtoTracksWithLabel const& tracks, FemtoMcParticlesWithLabel const& mcParticles, o2::aod::FMcMothers const& mcMothers, o2::aod::FMcPartMoths const& mcPartonicMothers) { - pairTrackTrackBuilder.processMixedEvent(cols, mcCols, tracks, trackWithLabelPartition1, trackWithLabelPartition2, mcParticles, mcMothers, mcPartonicMothers, cache, mixBinsVtxMult, mixBinsVtxCent, mixBinsVtxMultCent); + pairTrackTrackBuilder.processMixedEvent(cols, mcCols, tracks, trackWithLabelPartition1, trackWithLabelPartition2, mcParticles, mcMothers, mcPartonicMothers, cache, mixBinsVtxMult, mixBinsVtxCent, mixBinsVtxMultCent, mixBinsVtxCentEventPlaneAngle); } PROCESS_SWITCH(FemtoPairTrackTrack, processMixedEventMc, "Enable processing mixed event processing", false); };