diff --git a/PWGLF/Tasks/Strangeness/kinkAnalysis.cxx b/PWGLF/Tasks/Strangeness/kinkAnalysis.cxx index f01a9b3d11e..1f49710ba99 100644 --- a/PWGLF/Tasks/Strangeness/kinkAnalysis.cxx +++ b/PWGLF/Tasks/Strangeness/kinkAnalysis.cxx @@ -13,6 +13,7 @@ /// \author everyone #include +#include #include "Math/Vector4D.h" #include "Common/Core/CollisionAssociation.h" @@ -46,7 +47,6 @@ #include "ReconstructionDataFormats/Track.h" -#include using namespace std; using namespace o2; @@ -93,7 +93,8 @@ struct kinkAnalysis { Pion, Xi, OmegaToL, - OmegaToXi }; + OmegaToXi, + Hypertriton }; Configurable cfgCutEta{"cfgCutEta", 0.8f, "Eta range for tracks"}; Configurable cfgBz{"cfgBz", -999, "bz field, -999 is automatic"}; @@ -160,7 +161,7 @@ struct kinkAnalysis { float mMother, mChargedDaughter, mNeutralDaughter; - using CompleteTracks = soa::Join; + using CompleteTracks = soa::Join; using CompleteCollisions = soa::Join; struct TrackCand { @@ -172,7 +173,7 @@ struct kinkAnalysis { float angleCutFunction(int particleName, float x) { - float mSigmaMinus = 1.197449, mSigmaPlus = 1.189449, mKaon = 0.493677, mPion = 0.13957018, mXi = 1.321, mOmega = 1.67245; + float mSigmaMinus = 1.197449, mSigmaPlus = 1.189449, mKaon = 0.493677, mPion = 0.13957018, mXi = 1.321, mOmega = 1.67245, mHypertriton = 2.99131; float par1 = mSigmaMinus; /// Default value is for SigmaMinus float par2 = 0.8; float par3 = M_PI; @@ -211,9 +212,14 @@ struct kinkAnalysis { par1 = mPion; par2 = 0.2731374; break; + + case Hypertriton: + par1 = mHypertriton; + par2 = 0.07; + break; } - if ((particleName == SigmaMinus) || (particleName == SigmaPlusToPi) || (particleName == SigmaPlusToProton) || (particleName == Xi) || (particleName == OmegaToXi) || (particleName == OmegaToL)) + if ((particleName == SigmaMinus) || (particleName == SigmaPlusToPi) || (particleName == SigmaPlusToProton) || (particleName == Xi) || (particleName == OmegaToXi) || (particleName == OmegaToL) || (particleName == Hypertriton)) return par1 * (par2 / (sqrt((x * x) * (1 - (par2 * par2)) - ((par1 * par1) * (par2 * par2))))) * (180. / par3) + 1; return ((atan(par1 * par2 * (1.0 / (sqrt((x * x) * (1.0 - (par1 * par1)) - (par1 * par1) * (par2 * par2)))))) * 180.) / par3; @@ -271,6 +277,13 @@ struct kinkAnalysis { particlePdgCode = -3112; break; + case Hypertriton: + if (cfgMotherCharge == -1) + particlePdgCode = 1010010030; + else + particlePdgCode = -1010010030; + break; + default: particlePdgCode = 3112; } @@ -284,7 +297,11 @@ struct kinkAnalysis { const AxisSpec axisqT{1000, 0.0, 1.0, "q_{T}"}; const AxisSpec axisRmother{500, 0.0, 50.0, "#it{R}_{mother} (cm)"}; const AxisSpec axisDCAdaugh{1000, 0.0, 10.0, "#it{DCA}_{daughter} (cm)"}; - const AxisSpec axisSigmaMass{300, cfgLowerHistLimit, cfgUpperHistLimit, "#it{M}_{inv} (Gev/#it{c^{2}})"}; + + const AxisSpec axisSigmaMass{1000, 1.1, 1.4, "#it{M}_{inv} (Gev/#it{c^{2}})"}; + const AxisSpec hypAxisMass{100, 2.9, 3.6, "#it{M}_{inv} (Gev/#it{c^{2}})"}; + const AxisSpec axisPtHyp{100, 0., 15, "#it{p}_{T} (GeV/#it{c})"}; + const AxisSpec axisdtime{10000, 0, 50, "#delta t"}; const AxisSpec axisPt{500, 0., 50., "#it{p}_{T} (GeV/#it{c})"}; const AxisSpec axisdphi{360, 0, 360, "#delta phi"}; @@ -318,6 +335,13 @@ struct kinkAnalysis { histos.add("hPtMinusRecMcTrthM", "hPtMinusRecMcTrthM", kTH2F, {axisSigmaMass, axisPt}); histos.add("hPtMinusRecMcTrthelse", "hPtMinusRecMcTrthelse", kTH2F, {axisSigmaMass, axisPt}); + histos.add("hHypMass", "hHypMass", kTH1F, {hypAxisMass}); + histos.add("hHypMassMC", "hHypMassMC", kTH1F, {hypAxisMass}); + histos.add("hHypMassPt", "hHypMassPt", kTH2F, {hypAxisMass, axisPtHyp}); + + histos.add("hNSigmaTrVsPt", "hNSigmaTrVsPt", kTH2F, {axisPtHyp, {100, -5, 5, "nSigmaTPC"}}); + histos.add("hpRes", "hpRes", kTH2F, {axisPtHyp, {100, -0.5, 0.5, "p_{T} Res"}}); + lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->get(cfgLUTpath)); ft2.setMaxChi2(5); ft2.setUseAbsDCA(true); @@ -379,7 +403,6 @@ struct kinkAnalysis { TrackCand trForpool; std::array, 4> pools; // pools of positive and negative seeds sorted in min VtxID - std::vector selected(tracks.size(), 0u); std::vector globalBCvector; @@ -494,6 +517,7 @@ struct kinkAnalysis { const int poolCh1 = chargeM < 0 ? 1 : 0; const int poolCh2 = chargeD < 0 ? 1 : 0; int motherPdg = 999; + float mcMotherPt = 0; int daughterPdg = 777; switch (particleName) { @@ -538,6 +562,12 @@ struct kinkAnalysis { mChargedDaughter = 0.13957039; mNeutralDaughter = 1.31486; break; + + case Hypertriton: + mMother = 2.99131; + mChargedDaughter = 2.808921; + mNeutralDaughter = 0.1349766; + break; } o2::dataformats::VertexBase primaryVertex; @@ -574,6 +604,7 @@ struct kinkAnalysis { if ((seedM.mcParticleIdx != -1) && partTable) { auto mcParticle = partTable->rawIteratorAt(seedM.mcParticleIdx); motherPdg = mcParticle.pdgCode(); + mcMotherPt = mcParticle.pt(); } if (trackM.has_collision()) { @@ -612,6 +643,15 @@ struct kinkAnalysis { auto mcParticle = partTable->rawIteratorAt(seedD.mcParticleIdx); daughterPdg = mcParticle.pdgCode(); } + + bool isDaughter = false; + + if (cfgIsMC && (particleName == Hypertriton)) { + int tritPDG = 1000010030; + if (abs(daughterPdg) == tritPDG) + isDaughter = true; + } + PionTr = getTrackParCov(trackDgh); SigmaTr2.getPxPyPzGlo(sigmaPin); @@ -628,6 +668,10 @@ struct kinkAnalysis { if (trackDgh.tpcNSigmaKa() > cfgNsigmaTPCdaughter) continue; } + if (particleName == Hypertriton) { + if (trackDgh.tpcNSigmaTr() > cfgNsigmaTPCdaughter) + continue; + } gpu::gpustd::array dcaInfo; o2::base::Propagator::Instance()->propagateToDCABxByBz({primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, PionTr, 2.f, static_cast(cfgMaterialCorrection.value), &dcaInfo); @@ -686,6 +730,9 @@ struct kinkAnalysis { neutronPabs = sqrt(pow((sigmaPDC[2] - pionPDC[2]), 2) + pow((sigmaPDC[1] - pionPDC[1]), 2) + pow((sigmaPDC[0] - pionPDC[0]), 2)); neutronM = sqrt((sigmaE - pionE) * (sigmaE - pionE) - neutronPabs * neutronPabs); + if (abs(neutronM - mNeutralDaughter) / mNeutralDaughter > 0.1) + continue; + if ((particleName == SigmaMinus) || (particleName == SigmaPlusToPi)) { if ((theta * radToDeg > (angleCutFunction(particleName, sigmaPabsDC))) && (sigmaPabsDC > 1.6)) continue; @@ -720,12 +767,25 @@ struct kinkAnalysis { continue; } + if (particleName == Hypertriton) { + // if (sigmaPabsDC < 2.5) + // if (theta * radToDeg < 2.5) + // continue; + if (sigmaPabsDC > 0.6) + if (theta * radToDeg > angleCutFunction(particleName, sigmaPabsDC)) + continue; + } + if (!cfgIsMC) if (sigmaPt < 1.6) continue; - if (theta * radToDeg < 0.5) + if (particleName != Hypertriton) { + if (theta * radToDeg < 0.5) + continue; + } else if (theta * radToDeg < 0.2) { continue; + } if ((qT < qTlower) || (qT > qTupper)) continue; @@ -740,6 +800,21 @@ struct kinkAnalysis { mass = sqrt((neutronE + pionE) * (neutronE + pionE) - sigmaPabsDC * sigmaPabsDC); + if (particleName == Hypertriton) { + histos.fill(HIST("hHypMass"), mass); + histos.fill(HIST("hHypMassPt"), mass, sigmaPt); + + histos.fill(HIST("hNSigmaTrVsPt"), sigmaPt, trackDgh.tpcNSigmaTr()); + if (isDaughter) + histos.fill(HIST("hHypMassMC"), mass); + } + + if (((chargeM == -1) && (chargeD == -1)) || ((chargeM == 1) && (chargeD == 1))) { + if (cfgIsMC) { + histos.fill(HIST("hcodes"), motherPdg, daughterPdg); + } + } + if ((chargeM == -1) && (chargeD == -1)) { if (cfgIsMC) { histos.fill(HIST("hcodes"), motherPdg, daughterPdg); @@ -749,7 +824,9 @@ struct kinkAnalysis { } else if ((motherPdg == particlePdgCode || motherPdg == -3222) && (daughterPdg != -211)) { histos.fill(HIST("hptMtrue"), sigmaPt, PionTr.getPt()); histos.fill(HIST("hPtMinusRecMcTrthM"), mass, sigmaPt); - } else { // if ((motherPdg != particlePdgCode)&&(daughterPdg!=-211)) { + } else if ((motherPdg == particlePdgCode || motherPdg == -1010010030) && (daughterPdg == -1000010030 || daughterPdg == 1000010030)) { + histos.fill(HIST("hpRes"), sigmaPt, (mcMotherPt - sigmaPt) / mcMotherPt); + } else { histos.fill(HIST("hptMDelse"), sigmaPt, PionTr.getPt()); histos.fill(HIST("hPtMinusRecMcTrthelse"), mass, sigmaPt); }