diff --git a/PWGLF/DataModel/LFHyperNucleiKinkTables.h b/PWGLF/DataModel/LFHyperNucleiKinkTables.h index 797a1d3c109..4ffd6f4ec46 100644 --- a/PWGLF/DataModel/LFHyperNucleiKinkTables.h +++ b/PWGLF/DataModel/LFHyperNucleiKinkTables.h @@ -24,6 +24,7 @@ namespace o2::aod namespace hyperkink { +DECLARE_SOA_COLUMN(MagPolarity, magPolarity, int8_t); //! Magnetic field polarity DECLARE_SOA_COLUMN(XPV, xPV, float); //! Primary vertex of the candidate (x direction) DECLARE_SOA_COLUMN(YPV, yPV, float); //! Primary vertex of the candidate (y direction) DECLARE_SOA_COLUMN(ZPV, zPV, float); //! Primary vertex of the candidate (z direction) @@ -82,6 +83,7 @@ DECLARE_SOA_COLUMN(UpdatePzMothPV, updatePzMothPV, float); //! updated pz o DECLARE_SOA_TABLE(HypKinkCand, "AOD", "HYPKINKCANDS", o2::soa::Index<>, + hyperkink::MagPolarity, hyperkink::XPV, hyperkink::YPV, hyperkink::ZPV, hyperkink::XSV, hyperkink::YSV, hyperkink::ZSV, hyperkink::IsMatter, @@ -96,6 +98,7 @@ DECLARE_SOA_TABLE(HypKinkCand, "AOD", "HYPKINKCANDS", DECLARE_SOA_TABLE(MCHypKinkCand, "AOD", "MCHYPKINKCANDS", o2::soa::Index<>, + hyperkink::MagPolarity, hyperkink::XPV, hyperkink::YPV, hyperkink::ZPV, hyperkink::XSV, hyperkink::YSV, hyperkink::ZSV, hyperkink::IsMatter, diff --git a/PWGLF/TableProducer/Common/kinkBuilder.cxx b/PWGLF/TableProducer/Common/kinkBuilder.cxx index fe1cdc61a87..eb68fe85b39 100644 --- a/PWGLF/TableProducer/Common/kinkBuilder.cxx +++ b/PWGLF/TableProducer/Common/kinkBuilder.cxx @@ -91,6 +91,7 @@ struct kinkCandidate { struct kinkBuilder { enum PartType { kSigmaMinus = 0, + kHypertriton, kHyperhelium4sigma }; Produces outputDataTable; @@ -161,6 +162,11 @@ struct kinkBuilder { mothMass = o2::constants::physics::MassSigmaMinus; chargedDauMass = o2::constants::physics::MassPiMinus; neutDauMass = o2::constants::physics::MassNeutron; + } else if (hypoMoth == kHypertriton) { + charge = 1; + mothMass = o2::constants::physics::MassHyperTriton; + chargedDauMass = o2::constants::physics::MassTriton; + neutDauMass = o2::constants::physics::MassPi0; } else if (hypoMoth == kHyperhelium4sigma) { charge = 2; mothMass = o2::constants::physics::MassHyperHelium4; @@ -199,6 +205,8 @@ struct kinkBuilder { AxisSpec massAxis(100, 1.1, 1.4, "m (GeV/#it{c}^{2})"); if (hypoMoth == kSigmaMinus) { massAxis = AxisSpec{100, 1.1, 1.4, "m (GeV/#it{c}^{2})"}; + } else if (hypoMoth == kHypertriton) { + massAxis = AxisSpec{100, 2.94, 3.2, "m (GeV/#it{c}^{2})"}; } else if (hypoMoth == kHyperhelium4sigma) { massAxis = AxisSpec{100, 3.85, 4.25, "m (GeV/#it{c}^{2})"}; } diff --git a/PWGLF/TableProducer/Nuspex/CMakeLists.txt b/PWGLF/TableProducer/Nuspex/CMakeLists.txt index 66c91c16392..2f1295a8f4f 100644 --- a/PWGLF/TableProducer/Nuspex/CMakeLists.txt +++ b/PWGLF/TableProducer/Nuspex/CMakeLists.txt @@ -94,8 +94,8 @@ o2physics_add_dpl_workflow(nuclei-flow-trees PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DetectorsBase O2Physics::EventFilteringUtils COMPONENT_NAME Analysis) -o2physics_add_dpl_workflow(hyperhelium4sigma-reco-task - SOURCES hyperhelium4sigmaRecoTask.cxx +o2physics_add_dpl_workflow(hyperkink-reco-task + SOURCES hyperkinkRecoTask.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) diff --git a/PWGLF/TableProducer/Nuspex/hyperhelium4sigmaRecoTask.cxx b/PWGLF/TableProducer/Nuspex/hyperkinkRecoTask.cxx similarity index 53% rename from PWGLF/TableProducer/Nuspex/hyperhelium4sigmaRecoTask.cxx rename to PWGLF/TableProducer/Nuspex/hyperkinkRecoTask.cxx index 300fe513e8e..e465cee83fc 100644 --- a/PWGLF/TableProducer/Nuspex/hyperhelium4sigmaRecoTask.cxx +++ b/PWGLF/TableProducer/Nuspex/hyperkinkRecoTask.cxx @@ -9,8 +9,8 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. // -/// \file hyperhelium4sigmaRecoTask.cxx -/// \brief QA and analysis task for hyper-helium4sigma (HyperHe4S) +/// \file hyperkinkRecoTask.cxx +/// \brief QA and analysis task for kink decay of hypernuclei /// \author Yuanzhe Wang #include "PWGLF/DataModel/LFHyperNucleiKinkTables.h" @@ -48,22 +48,15 @@ using MCLabeledCollisionsFull = soa::Join; using MCLabeledTracksIU = soa::Join; -enum Channel { - k2body = 0, // helium4, pion0 - k3body_p, // triton, proton, pion0 - k3body_n, // triton, neutron, pion+ - kNDecayChannel +enum PartType { + kHypertriton = 0, + kHyperhelium4sigma }; enum DaughterType { - kDaugAlpha = 0, - kDaugTriton, - kDaugProton, - kDaugChargedPion, - kDaugNeutron, - kDaugPion0, - kNDaughterType, - kNChargedDaughterType = kDaugNeutron + kDaugCharged = 0, + kDaugNeutral, + kNDaughterType }; namespace @@ -72,104 +65,169 @@ constexpr std::array LayerRadii{2.33959f, 3.14076f, 3.91924f, 19.6213f constexpr int kITSLayers = 7; constexpr int kITSInnerBarrelLayers = 3; // constexpr int kITSOuterBarrelLayers = 4; -constexpr std::array kDaugghterPDG = { - o2::constants::physics::Pdg::kAlpha, - o2::constants::physics::Pdg::kTriton, - PDG_t::kProton, - PDG_t::kPiPlus, - PDG_t::kNeutron, - PDG_t::kPi0}; const std::array covPosSV{6.4462712107237135f, 0.1309793068144521f, 6.626654155592929f, -0.4510297694023185f, 0.16996629627762413f, 4.109195981415627f}; std::shared_ptr hMothCounter; -std::shared_ptr hMoth2BCounter; -std::shared_ptr hDaugCounter[kNChargedDaughterType]; -std::shared_ptr hDaugTPCNSigma[kNChargedDaughterType]; +std::shared_ptr hDaugCounter; +std::shared_ptr hDaugTPCNSigma; std::shared_ptr hRecoMothCounter; -std::shared_ptr hRecoDaugAlphaCounter; +std::shared_ptr hRecoDaugCounter; } // namespace //-------------------------------------------------------------- -// Check the decay channel of hyperhelium4sigma -template -Channel getDecayChannelHe4S(TMCParticle const& particle, std::vector& list) -{ - if (std::abs(particle.pdgCode()) != o2::constants::physics::Pdg::kHyperHelium4Sigma) { - return kNDecayChannel; - } +struct H3LDecay { + enum Channel { + k2bodyNeutral = 0, // triton, pion0 + k2bodyCharged, + k3bodyCharged, + kNChannel + }; - // list: charged (alpha or triton), charged or empty, neutral - list.clear(); - list.resize(3, -1); - - bool haveAlpha = false, haveTriton = false, haveProton = false, haveNeuteron = false; - bool haveAntiAlpha = false, haveAntiTriton = false, haveAntiProton = false, haveAntiNeuteron = false; - bool havePionPlus = false, havePionMinus = false, havePion0 = false; - for (const auto& mcDaughter : particle.template daughters_as()) { - if (mcDaughter.pdgCode() == o2::constants::physics::Pdg::kAlpha) { - haveAlpha = true; - list[0] = mcDaughter.globalIndex(); - } - if (mcDaughter.pdgCode() == -o2::constants::physics::Pdg::kAlpha) { - haveAntiAlpha = true; - list[0] = mcDaughter.globalIndex(); - } - if (mcDaughter.pdgCode() == o2::constants::physics::Pdg::kTriton) { - haveTriton = true; - list[0] = mcDaughter.globalIndex(); - } - if (mcDaughter.pdgCode() == -o2::constants::physics::Pdg::kTriton) { - haveAntiTriton = true; - list[0] = mcDaughter.globalIndex(); - } - if (mcDaughter.pdgCode() == PDG_t::kProton) { - haveProton = true; - list[1] = mcDaughter.globalIndex(); - } - if (mcDaughter.pdgCode() == -PDG_t::kProton) { - haveAntiProton = true; - list[1] = mcDaughter.globalIndex(); + template + static Channel getDecayChannel(TMCParticle const& particle, std::vector& list) + { + if (std::abs(particle.pdgCode()) != o2::constants::physics::Pdg::kHyperTriton) { + return kNChannel; } - if (mcDaughter.pdgCode() == PDG_t::kNeutron) { - haveNeuteron = true; - list[2] = mcDaughter.globalIndex(); + + list.clear(); + list.resize(2, -1); + + bool haveHelium3 = false, haveAntiHelium3 = false, haveDeuteron = false, haveAntiDeuteron = false; + bool haveProton = false, haveAntiProton = false, havePionPlus = false, havePionMinus = false; + bool haveTriton = false, haveAntiTriton = false, havePion0 = false; + for (const auto& mcDaughter : particle.template daughters_as()) { + if (mcDaughter.pdgCode() == o2::constants::physics::Pdg::kTriton) { + haveTriton = true; + list[0] = mcDaughter.globalIndex(); + } + if (mcDaughter.pdgCode() == -o2::constants::physics::Pdg::kTriton) { + haveAntiTriton = true; + list[0] = mcDaughter.globalIndex(); + } + if (mcDaughter.pdgCode() == PDG_t::kPi0) { + havePion0 = true; + list[1] = mcDaughter.globalIndex(); + } + if (mcDaughter.pdgCode() == o2::constants::physics::Pdg::kHelium3) { + haveHelium3 = true; + } + if (mcDaughter.pdgCode() == -o2::constants::physics::Pdg::kHelium3) { + haveAntiHelium3 = true; + } + if (mcDaughter.pdgCode() == o2::constants::physics::Pdg::kDeuteron) { + haveDeuteron = true; + } + if (mcDaughter.pdgCode() == -o2::constants::physics::Pdg::kDeuteron) { + haveAntiDeuteron = true; + } + if (mcDaughter.pdgCode() == PDG_t::kProton) { + haveProton = true; + } + if (mcDaughter.pdgCode() == -PDG_t::kProton) { + haveAntiProton = true; + } + if (mcDaughter.pdgCode() == PDG_t::kPiPlus) { + havePionPlus = true; + } + if (mcDaughter.pdgCode() == -PDG_t::kPiPlus) { + havePionMinus = true; + } } - if (mcDaughter.pdgCode() == -PDG_t::kNeutron) { - haveAntiNeuteron = true; - list[2] = mcDaughter.globalIndex(); + + if ((haveTriton && havePion0) || (haveAntiTriton && havePion0)) { + return H3LDecay::k2bodyNeutral; + } else if ((haveHelium3 && havePion0) || (haveAntiHelium3 && havePion0)) { + return H3LDecay::k2bodyCharged; + } else if ((haveDeuteron && haveProton && havePionPlus) || (haveAntiDeuteron && haveAntiProton && havePionMinus)) { + return H3LDecay::k3bodyCharged; + } else { + return kNChannel; } - if (mcDaughter.pdgCode() == PDG_t::kPiPlus) { - havePionPlus = true; - list[1] = mcDaughter.globalIndex(); + } +}; + +//-------------------------------------------------------------- +// Check the decay channel of hyperhelium4sigma +struct He4SDecay { + enum Channel { + k2body = 0, // helium4, pion0 + k3body_p, // triton, proton, pion0 + k3body_n, // triton, neutron, pion+ + kNChannel + }; + + template + static Channel getDecayChannel(TMCParticle const& particle, std::vector& list) + { + if (std::abs(particle.pdgCode()) != o2::constants::physics::Pdg::kHyperHelium4Sigma) { + return kNChannel; } - if (mcDaughter.pdgCode() == -PDG_t::kPiPlus) { - havePionMinus = true; - list[1] = mcDaughter.globalIndex(); + + list.clear(); + list.resize(2, -1); + + bool haveAlpha = false, haveTriton = false, haveProton = false, haveNeuteron = false; + bool haveAntiAlpha = false, haveAntiTriton = false, haveAntiProton = false, haveAntiNeuteron = false; + bool havePionPlus = false, havePionMinus = false, havePion0 = false; + for (const auto& mcDaughter : particle.template daughters_as()) { + if (mcDaughter.pdgCode() == o2::constants::physics::Pdg::kAlpha) { + haveAlpha = true; + list[0] = mcDaughter.globalIndex(); + } + if (mcDaughter.pdgCode() == -o2::constants::physics::Pdg::kAlpha) { + haveAntiAlpha = true; + list[0] = mcDaughter.globalIndex(); + } + if (mcDaughter.pdgCode() == o2::constants::physics::Pdg::kTriton) { + haveTriton = true; + } + if (mcDaughter.pdgCode() == -o2::constants::physics::Pdg::kTriton) { + haveAntiTriton = true; + } + if (mcDaughter.pdgCode() == PDG_t::kProton) { + haveProton = true; + } + if (mcDaughter.pdgCode() == -PDG_t::kProton) { + haveAntiProton = true; + } + if (mcDaughter.pdgCode() == PDG_t::kNeutron) { + haveNeuteron = true; + } + if (mcDaughter.pdgCode() == -PDG_t::kNeutron) { + haveAntiNeuteron = true; + } + if (mcDaughter.pdgCode() == PDG_t::kPiPlus) { + havePionPlus = true; + } + if (mcDaughter.pdgCode() == -PDG_t::kPiPlus) { + havePionMinus = true; + } + if (mcDaughter.pdgCode() == PDG_t::kPi0) { + havePion0 = true; + list[1] = mcDaughter.globalIndex(); + } } - if (mcDaughter.pdgCode() == PDG_t::kPi0) { - havePion0 = true; - list[2] = mcDaughter.globalIndex(); + + if ((haveAlpha && havePion0) || (haveAntiAlpha && havePion0)) { + return He4SDecay::k2body; + } else if ((haveTriton && haveProton && havePion0) || (haveAntiTriton && haveAntiProton && havePion0)) { + return He4SDecay::k3body_p; + } else if ((haveTriton && haveNeuteron && havePionPlus) || (haveAntiTriton && haveAntiNeuteron && havePionMinus)) { + return He4SDecay::k3body_n; } - } - if ((haveAlpha && havePion0) || (haveAntiAlpha && havePion0)) { - return k2body; - } else if ((haveTriton && haveProton && havePion0) || (haveAntiTriton && haveAntiProton && havePion0)) { - return k3body_p; - } else if ((haveTriton && haveNeuteron && havePionPlus) || (haveAntiTriton && haveAntiNeuteron && havePionMinus)) { - return k3body_n; + return kNChannel; } - - return kNDecayChannel; -} +}; //-------------------------------------------------------------- // Extract track parameters from a mcparticle, use global coordinates as the local one template o2::track::TrackParametrization getTrackParFromMC(const T& mcparticle, int charge = 1) { - int sign = mcparticle.pdgCode() > 0 ? 1 : -1; // ok for hyperhelium4sigma + int sign = mcparticle.pdgCode() > 0 ? 1 : -1; // ok for hypernuclei TrackPrecision snp = mcparticle.py() / (mcparticle.pt() + 1.e-10f); TrackPrecision tgl = mcparticle.pz() / (mcparticle.pt() + 1.e-10f); std::array arraypar = {mcparticle.vy(), mcparticle.vz(), snp, @@ -201,25 +259,38 @@ void setTrackIDForMC(std::vector& mcPartIndices, aod::McParticles const } } +//-------------------------------------------------------------- +// get ITSNSigma for daughter track +template +float getITSNSigma(const TTrack& track, o2::aod::ITSResponse& itsResponse, o2::track::PID partType) +{ + float nSigma = -999.f; + switch (partType) { + case o2::track::PID::Alpha: + nSigma = itsResponse.nSigmaITS(track); + break; + case o2::track::PID::Triton: + nSigma = itsResponse.nSigmaITS(track); + break; + default: + break; + } + return nSigma; +} + //-------------------------------------------------------------- // get TPCNSigma for daughter track template -float getTPCNSigma(const TTrack& track, const int daughterType) +float getTPCNSigma(const TTrack& track, o2::track::PID partType) { float nSigma = -999.f; - switch (daughterType) { - case kDaugAlpha: + switch (partType) { + case o2::track::PID::Alpha: nSigma = track.tpcNSigmaAl(); break; - case kDaugTriton: + case o2::track::PID::Triton: nSigma = track.tpcNSigmaTr(); break; - case kDaugProton: - nSigma = track.tpcNSigmaPr(); - break; - case kDaugChargedPion: - nSigma = track.tpcNSigmaPi(); - break; default: break; } @@ -256,7 +327,7 @@ bool refitMotherTrack(const TTrack& track, o2::track::TrackParametrizationWithEr } //-------------------------------------------------------------- -struct Hyphe4sCandidate { +struct HypKinkCandidate { bool isMatter = false; @@ -293,8 +364,8 @@ struct Hyphe4sCandidate { }; //-------------------------------------------------------------- -// analysis task for hyperhelium4sigma 2-body decay -struct Hyperhelium4sigmaRecoTask { +// analysis task for 2-body kink decay +struct HyperkinkRecoTask { Produces outputDataTable; Produces outputMCTable; @@ -302,15 +373,16 @@ struct Hyperhelium4sigmaRecoTask { Service ccdb; o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrLUT; - std::vector mcHe4sIndices; + std::vector mcPartIndices; // Histograms are defined with HistogramRegistry HistogramRegistry registry{"registry", {}}; + Configurable hypoMoth{"hypoMoth", kHypertriton, "Mother particle hypothesis"}; // Configurable for event selection Configurable doEventCut{"doEventCut", true, "Apply event selection"}; Configurable maxZVertex{"maxZVertex", 10.0f, "Accepted z-vertex range (cm)"}; - Configurable cutNSigmaAl{"cutNSigmaAl", 5, "NSigmaTPCAlpha"}; + Configurable cutNSigmaDaug{"cutNSigmaDaug", 5, "TPC NSigmaTPC cut for daughter tracks"}; // CCDB options Configurable inputBz{"inputBz", -999, "bz field, -999 is automatic"}; @@ -324,15 +396,37 @@ struct Hyperhelium4sigmaRecoTask { o2::aod::ITSResponse itsResponse; - void init(InitContext const&) + float massChargedDaug = 999.f; + float massNeutralDaug = 999.f; + int pdgMoth = 0; + std::array pdgDaug = {o2::constants::physics::Pdg::kTriton, PDG_t::kPi0}; // pdgcode of charged (0) and neutral (1) daughter particles + o2::track::PID pidTypeDaug = o2::track::PID::Triton; + + void init(InitContext&) { + if (hypoMoth == kHypertriton) { + massChargedDaug = o2::constants::physics::MassTriton; + massNeutralDaug = o2::constants::physics::MassPi0; + pdgMoth = o2::constants::physics::Pdg::kHyperTriton; + pdgDaug[kDaugCharged] = o2::constants::physics::Pdg::kTriton; + pdgDaug[kDaugNeutral] = PDG_t::kPi0; + pidTypeDaug = o2::track::PID::Triton; + } else if (hypoMoth == kHyperhelium4sigma) { + massChargedDaug = o2::constants::physics::MassAlpha; + massNeutralDaug = o2::constants::physics::MassPi0; + pdgMoth = o2::constants::physics::Pdg::kHyperHelium4Sigma; + pdgDaug[kDaugCharged] = o2::constants::physics::Pdg::kAlpha; + pdgDaug[kDaugNeutral] = PDG_t::kPi0; + pidTypeDaug = o2::track::PID::Alpha; + } else { + LOG(fatal) << "Unknown mother particle hypothesis"; + } + // Axes const AxisSpec vertexZAxis{100, -15., 15., "vtx_{Z} [cm]"}; const AxisSpec ptAxis{50, -10, 10, "#it{p}_{T} (GeV/#it{c})"}; - const AxisSpec nSigmaAxis{120, -6.f, 6.f, "n#sigma_{#alpha}"}; + const AxisSpec nSigmaAxis{120, -6.f, 6.f, "n#sigma"}; const AxisSpec massAxis{100, 3.85, 4.25, "m (GeV/#it{c}^{2})"}; - const AxisSpec diffPxAxis{200, -10.f, 10.f, "#Delta #it{p}_{x} (GeV/#it{c})"}; - const AxisSpec diffPyAxis{200, -10.f, 10.f, "#Delta #it{p}_{y} (GeV/#it{c})"}; const AxisSpec diffPtAxis{200, -10.f, 10.f, "#Delta #it{p}_{T} (GeV/#it{c})"}; const AxisSpec diffPzAxis{200, -10.f, 10.f, "#Delta #it{p}_{z} (GeV/#it{c})"}; const AxisSpec radiusAxis{40, 0.f, 40.f, "R (cm)"}; @@ -349,24 +443,20 @@ struct Hyperhelium4sigmaRecoTask { registry.add("hDiffSVy", ";#Delta y (cm);", HistType::kTH1F, {{200, -10, 10}}); registry.add("hDiffSVz", ";#Delta z (cm);", HistType::kTH1F, {{200, -10, 10}}); registry.add("h2RecSVRVsTrueSVR", ";Reconstruced SV R (cm);True SV R (cm);", HistType::kTH2F, {radiusAxis, radiusAxis}); - registry.add("h2TrueMotherDiffPxVsRecSVR", ";Reconstruced SV R (cm);#Delta #it{p}_{T} (GeV/#it{c});", HistType::kTH2F, {radiusAxis, diffPxAxis}); - registry.add("h2TrueMotherDiffPyVsRecSVR", ";Reconstruced SV R (cm);#Delta #it{p}_{T} (GeV/#it{c});", HistType::kTH2F, {radiusAxis, diffPyAxis}); registry.add("h2TrueMotherDiffPtVsRecSVR", ";Reconstruced SV R (cm);#Delta #it{p}_{T} (GeV/#it{c});", HistType::kTH2F, {radiusAxis, diffPtAxis}); - registry.add("h2TrueMotherDiffPzVsRecSVR", ";Reconstruced SV R (cm);#Delta #it{p}_{z} (GeV/#it{c});", HistType::kTH2F, {radiusAxis, diffPzAxis}); - registry.add("h2TrueMotherDiffTglVsRecSVR", ";Reconstruced SV R (cm);#Delta tan#lambda;", HistType::kTH2F, {radiusAxis, {200, -0.1, 0.1}}); registry.add("h2TrueMotherDiffEtaVsRecSVR", ";Reconstruced SV R (cm);#Delta #eta;", HistType::kTH2F, {radiusAxis, {200, -0.1, 0.1}}); registry.add("hDiffDauPx", ";#Delta p_{x} (GeV/#it{c}); ", HistType::kTH1D, {{200, -10, 10}}); registry.add("hDiffDauPy", ";#Delta p_{y} (GeV/#it{c}); ", HistType::kTH1D, {{200, -10, 10}}); registry.add("hDiffDauPz", ";#Delta p_{z} (GeV/#it{c}); ", HistType::kTH1D, {{200, -10, 10}}); registry.add("h2TrueSignalMassPt", "h2TrueSignalMassPt", HistType::kTH2F, {{ptAxis, massAxis}}); - registry.add("h2TrueSignalNSigmaAlPt", "h2TrueSignalNSigmaAlPt", HistType::kTH2F, {{ptAxis, nSigmaAxis}}); + registry.add("h2TrueDaugTPCNSigmaPt", "h2TrueDaugTPCNSigmaPt", HistType::kTH2F, {{ptAxis, nSigmaAxis}}); registry.add("hDCAXYMothToRecSV", "hDCAXYMothToRecSV", HistType::kTH1F, {{200, -10, 10}}); registry.add("hDCAZMothToRecSV", "hDCAZMothToRecSV", HistType::kTH1F, {{200, -10, 10}}); } - registry.add("h2MassHyperhelium4sigmaPt", "h2MassHyperhelium4sigmaPt", HistType::kTH2F, {{ptAxis, massAxis}}); - registry.add("h2NSigmaAlPt", "h2NSigmaAlPt", HistType::kTH2F, {{ptAxis, nSigmaAxis}}); + registry.add("h2MothMassPt", "h2MothMassPt", HistType::kTH2F, {{ptAxis, massAxis}}); + registry.add("h2DaugTPCNSigmaPt", "h2DaugTPCNSigmaPt", HistType::kTH2F, {{ptAxis, nSigmaAxis}}); ccdb->setURL(ccdbPath); ccdb->setCaching(true); @@ -394,32 +484,32 @@ struct Hyperhelium4sigmaRecoTask { } template - void fillCandidate(Hyphe4sCandidate& hyphe4sCand, TCollision const& collision, TKindCandidate const& kinkCand, TTrack const& trackMoth, TTrack const& trackDaug) + void fillCandidate(HypKinkCandidate& hypkinkCand, TCollision const& collision, TKindCandidate const& kinkCand, TTrack const& trackMoth, TTrack const& trackDaug) { - hyphe4sCand.isMatter = kinkCand.mothSign() > 0; - hyphe4sCand.posPV[0] = collision.posX(); - hyphe4sCand.posPV[1] = collision.posY(); - hyphe4sCand.posPV[2] = collision.posZ(); - hyphe4sCand.posSV[0] = kinkCand.xDecVtx() + collision.posX(); - hyphe4sCand.posSV[1] = kinkCand.yDecVtx() + collision.posY(); - hyphe4sCand.posSV[2] = kinkCand.zDecVtx() + collision.posZ(); - - hyphe4sCand.momMothSV[0] = kinkCand.pxMoth(); - hyphe4sCand.momMothSV[1] = kinkCand.pyMoth(); - hyphe4sCand.momMothSV[2] = kinkCand.pzMoth(); - hyphe4sCand.momDaugSV[0] = kinkCand.pxDaug(); - hyphe4sCand.momDaugSV[1] = kinkCand.pyDaug(); - hyphe4sCand.momDaugSV[2] = kinkCand.pzDaug(); - - hyphe4sCand.dcaXYMothPv = kinkCand.dcaMothPv(); - hyphe4sCand.dcaXYDaugPv = kinkCand.dcaDaugPv(); - hyphe4sCand.dcaKinkTopo = kinkCand.dcaKinkTopo(); - - fillCandidateRecoMoth(hyphe4sCand, collision, trackMoth); - - hyphe4sCand.itsClusterSizeDaug = trackDaug.itsClusterSizes(); - hyphe4sCand.nSigmaTPCDaug = trackDaug.tpcNSigmaAl(); - hyphe4sCand.nSigmaITSDaug = itsResponse.nSigmaITS(trackDaug); + hypkinkCand.isMatter = kinkCand.mothSign() > 0; + hypkinkCand.posPV[0] = collision.posX(); + hypkinkCand.posPV[1] = collision.posY(); + hypkinkCand.posPV[2] = collision.posZ(); + hypkinkCand.posSV[0] = kinkCand.xDecVtx() + collision.posX(); + hypkinkCand.posSV[1] = kinkCand.yDecVtx() + collision.posY(); + hypkinkCand.posSV[2] = kinkCand.zDecVtx() + collision.posZ(); + + hypkinkCand.momMothSV[0] = kinkCand.pxMoth(); + hypkinkCand.momMothSV[1] = kinkCand.pyMoth(); + hypkinkCand.momMothSV[2] = kinkCand.pzMoth(); + hypkinkCand.momDaugSV[0] = kinkCand.pxDaug(); + hypkinkCand.momDaugSV[1] = kinkCand.pyDaug(); + hypkinkCand.momDaugSV[2] = kinkCand.pzDaug(); + + hypkinkCand.dcaXYMothPv = kinkCand.dcaMothPv(); + hypkinkCand.dcaXYDaugPv = kinkCand.dcaDaugPv(); + hypkinkCand.dcaKinkTopo = kinkCand.dcaKinkTopo(); + + fillCandidateRecoMoth(hypkinkCand, collision, trackMoth); + + hypkinkCand.itsClusterSizeDaug = trackDaug.itsClusterSizes(); + hypkinkCand.nSigmaTPCDaug = getTPCNSigma(trackDaug, pidTypeDaug); + hypkinkCand.nSigmaITSDaug = getITSNSigma(trackDaug, itsResponse, pidTypeDaug); int lastLayerMoth = 0; for (int i = 6; i >= 0; i--) { @@ -432,17 +522,17 @@ struct Hyperhelium4sigmaRecoTask { o2::base::Propagator::Instance()->PropagateToXBxByBz(trackparMother, LayerRadii[lastLayerMoth]); std::array posLastHit{-999.f}; trackparMother.getXYZGlo(posLastHit); - hyphe4sCand.lastPosMoth[0] = posLastHit[0]; - hyphe4sCand.lastPosMoth[1] = posLastHit[1]; - hyphe4sCand.lastPosMoth[2] = posLastHit[2]; + hypkinkCand.lastPosMoth[0] = posLastHit[0]; + hypkinkCand.lastPosMoth[1] = posLastHit[1]; + hypkinkCand.lastPosMoth[2] = posLastHit[2]; } template - void fillCandidateRecoMoth(Hyphe4sCandidate& hyphe4sCand, TCollision const& collision, TTrack const& trackMoth) + void fillCandidateRecoMoth(HypKinkCandidate& hypkinkCand, TCollision const& collision, TTrack const& trackMoth) { - hyphe4sCand.isMothReco = true; - hyphe4sCand.chi2ITSMoth = trackMoth.itsChi2NCl(); - hyphe4sCand.itsClusterSizeMoth = trackMoth.itsClusterSizes(); + hypkinkCand.isMothReco = true; + hypkinkCand.chi2ITSMoth = trackMoth.itsChi2NCl(); + hypkinkCand.itsClusterSizeMoth = trackMoth.itsClusterSizes(); auto motherTrackPar = getTrackParCov(trackMoth); o2::dataformats::VertexBase primaryVtx = {{collision.posX(), collision.posY(), collision.posZ()}, {collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()}}; @@ -450,34 +540,34 @@ struct Hyperhelium4sigmaRecoTask { if (o2::base::Propagator::Instance()->propagateToDCABxByBz(primaryVtx, motherTrackPar, 2.f, o2::base::Propagator::MatCorrType::USEMatCorrLUT)) { motherTrackPar.getPxPyPzGlo(pMotherPv); } - hyphe4sCand.momMothPV[0] = pMotherPv[0]; - hyphe4sCand.momMothPV[1] = pMotherPv[1]; - hyphe4sCand.momMothPV[2] = pMotherPv[2]; + hypkinkCand.momMothPV[0] = pMotherPv[0]; + hypkinkCand.momMothPV[1] = pMotherPv[1]; + hypkinkCand.momMothPV[2] = pMotherPv[2]; std::array updatePMotherPv = {-999.f}; if (motherTrackPar.update(primaryVtx)) { motherTrackPar.getPxPyPzGlo(updatePMotherPv); } - hyphe4sCand.updateMomMothPV[0] = updatePMotherPv[0]; - hyphe4sCand.updateMomMothPV[1] = updatePMotherPv[1]; - hyphe4sCand.updateMomMothPV[2] = updatePMotherPv[2]; + hypkinkCand.updateMomMothPV[0] = updatePMotherPv[0]; + hypkinkCand.updateMomMothPV[1] = updatePMotherPv[1]; + hypkinkCand.updateMomMothPV[2] = updatePMotherPv[2]; } template - void fillCandidateMCInfo(Hyphe4sCandidate& hyphe4sCand, TMCParticle const& mcMothTrack, TMCParticle const& mcDauTrack, TMCParticle const& mcNeutDauTrack) + void fillCandidateMCInfo(HypKinkCandidate& hypkinkCand, TMCParticle const& mcMothTrack, TMCParticle const& mcDaugTrack, TMCParticle const& mcNeutDauTrack) { - hyphe4sCand.truePosSV[0] = mcDauTrack.vx(); - hyphe4sCand.truePosSV[1] = mcDauTrack.vy(); - hyphe4sCand.truePosSV[2] = mcDauTrack.vz(); - hyphe4sCand.trueMomMothPV[0] = mcMothTrack.px(); - hyphe4sCand.trueMomMothPV[1] = mcMothTrack.py(); - hyphe4sCand.trueMomMothPV[2] = mcMothTrack.pz(); - hyphe4sCand.trueMomMothSV[0] = mcDauTrack.px() + mcNeutDauTrack.px(); - hyphe4sCand.trueMomMothSV[1] = mcDauTrack.py() + mcNeutDauTrack.py(); - hyphe4sCand.trueMomMothSV[2] = mcDauTrack.pz() + mcNeutDauTrack.pz(); - hyphe4sCand.trueMomDaugSV[0] = mcDauTrack.px(); - hyphe4sCand.trueMomDaugSV[1] = mcDauTrack.py(); - hyphe4sCand.trueMomDaugSV[2] = mcDauTrack.pz(); + hypkinkCand.truePosSV[0] = mcDaugTrack.vx(); + hypkinkCand.truePosSV[1] = mcDaugTrack.vy(); + hypkinkCand.truePosSV[2] = mcDaugTrack.vz(); + hypkinkCand.trueMomMothPV[0] = mcMothTrack.px(); + hypkinkCand.trueMomMothPV[1] = mcMothTrack.py(); + hypkinkCand.trueMomMothPV[2] = mcMothTrack.pz(); + hypkinkCand.trueMomMothSV[0] = mcDaugTrack.px() + mcNeutDauTrack.px(); + hypkinkCand.trueMomMothSV[1] = mcDaugTrack.py() + mcNeutDauTrack.py(); + hypkinkCand.trueMomMothSV[2] = mcDaugTrack.pz() + mcNeutDauTrack.pz(); + hypkinkCand.trueMomDaugSV[0] = mcDaugTrack.px(); + hypkinkCand.trueMomDaugSV[1] = mcDaugTrack.py(); + hypkinkCand.trueMomDaugSV[2] = mcDaugTrack.pz(); } void processData(CollisionsFull const& collisions, aod::KinkCands const& KinkCands, FullTracksExtIU const&, aod::BCs const&) @@ -500,19 +590,20 @@ struct Hyperhelium4sigmaRecoTask { registry.fill(HIST("hCandidateCounter"), 1); auto dauTrack = kinkCand.trackDaug_as(); - if (std::abs(dauTrack.tpcNSigmaAl()) > cutNSigmaAl) { + float tpcNSigmaDaug = getTPCNSigma(dauTrack, pidTypeDaug); + if (std::abs(tpcNSigmaDaug) > cutNSigmaDaug) { continue; } - float invMass = RecoDecay::m(std::array{std::array{kinkCand.pxDaug(), kinkCand.pyDaug(), kinkCand.pzDaug()}, std::array{kinkCand.pxDaugNeut(), kinkCand.pyDaugNeut(), kinkCand.pzDaugNeut()}}, std::array{o2::constants::physics::MassAlpha, o2::constants::physics::MassPi0}); + float invMass = RecoDecay::m(std::array{std::array{kinkCand.pxDaug(), kinkCand.pyDaug(), kinkCand.pzDaug()}, std::array{kinkCand.pxDaugNeut(), kinkCand.pyDaugNeut(), kinkCand.pzDaugNeut()}}, std::array{massChargedDaug, massNeutralDaug}); registry.fill(HIST("hCandidateCounter"), 2); - registry.fill(HIST("h2MassHyperhelium4sigmaPt"), kinkCand.mothSign() * kinkCand.ptMoth(), invMass); - registry.fill(HIST("h2NSigmaAlPt"), kinkCand.mothSign() * kinkCand.ptDaug(), dauTrack.tpcNSigmaAl()); + registry.fill(HIST("h2MothMassPt"), kinkCand.mothSign() * kinkCand.ptMoth(), invMass); + registry.fill(HIST("h2DaugTPCNSigmaPt"), kinkCand.mothSign() * kinkCand.ptDaug(), tpcNSigmaDaug); auto bc = collision.bc_as(); initCCDB(bc); auto motherTrack = kinkCand.trackMoth_as(); - Hyphe4sCandidate hyphe4sCand; - fillCandidate(hyphe4sCand, collision, kinkCand, motherTrack, dauTrack); + HypKinkCandidate hypkinkCand; + fillCandidate(hypkinkCand, collision, kinkCand, motherTrack, dauTrack); o2::dataformats::VertexBase primaryVtx = {{collision.posX(), collision.posY(), collision.posZ()}, {collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()}}; o2::dataformats::VertexBase secondaryVtx = {{kinkCand.xDecVtx() + collision.posX(), kinkCand.yDecVtx() + collision.posY(), kinkCand.zDecVtx() + collision.posZ()}, {covPosSV[0], covPosSV[1], covPosSV[2], covPosSV[3], covPosSV[4], covPosSV[5]}}; @@ -526,30 +617,31 @@ struct Hyperhelium4sigmaRecoTask { } outputDataTable( - hyphe4sCand.posPV[0], hyphe4sCand.posPV[1], hyphe4sCand.posPV[2], - hyphe4sCand.posSV[0], hyphe4sCand.posSV[1], hyphe4sCand.posSV[2], - hyphe4sCand.isMatter, - hyphe4sCand.lastPosMoth[0], hyphe4sCand.lastPosMoth[1], hyphe4sCand.lastPosMoth[2], - hyphe4sCand.momMothSV[0], hyphe4sCand.momMothSV[1], hyphe4sCand.momMothSV[2], + mBz > 0 ? 1 : -1, + hypkinkCand.posPV[0], hypkinkCand.posPV[1], hypkinkCand.posPV[2], + hypkinkCand.posSV[0], hypkinkCand.posSV[1], hypkinkCand.posSV[2], + hypkinkCand.isMatter, + hypkinkCand.lastPosMoth[0], hypkinkCand.lastPosMoth[1], hypkinkCand.lastPosMoth[2], + hypkinkCand.momMothSV[0], hypkinkCand.momMothSV[1], hypkinkCand.momMothSV[2], refitPPV[0], refitPPV[1], refitPPV[2], refitPSV[0], refitPSV[1], refitPSV[2], - hyphe4sCand.momDaugSV[0], hyphe4sCand.momDaugSV[1], hyphe4sCand.momDaugSV[2], - hyphe4sCand.dcaXYMothPv, hyphe4sCand.dcaXYDaugPv, hyphe4sCand.dcaKinkTopo, - hyphe4sCand.chi2ITSMoth, hyphe4sCand.itsClusterSizeMoth, hyphe4sCand.itsClusterSizeDaug, - hyphe4sCand.nSigmaTPCDaug, hyphe4sCand.nSigmaITSDaug); + hypkinkCand.momDaugSV[0], hypkinkCand.momDaugSV[1], hypkinkCand.momDaugSV[2], + hypkinkCand.dcaXYMothPv, hypkinkCand.dcaXYDaugPv, hypkinkCand.dcaKinkTopo, + hypkinkCand.chi2ITSMoth, hypkinkCand.itsClusterSizeMoth, hypkinkCand.itsClusterSizeDaug, + hypkinkCand.nSigmaTPCDaug, hypkinkCand.nSigmaITSDaug); } } - PROCESS_SWITCH(Hyperhelium4sigmaRecoTask, processData, "process data", true); + PROCESS_SWITCH(HyperkinkRecoTask, processData, "process data", true); void processMC(MCLabeledCollisionsFull const& collisions, aod::KinkCands const& KinkCands, MCLabeledTracksIU const& tracks, aod::McParticles const& particlesMC, aod::McCollisions const& mcCollisions, aod::BCs const&) { - mcHe4sIndices.clear(); + mcPartIndices.clear(); std::vector mcPartIndices; setTrackIDForMC(mcPartIndices, particlesMC, tracks); std::vector isReconstructedMCCollisions(mcCollisions.size(), false); std::vector isSelectedMCCollisions(mcCollisions.size(), false); std::vector isGoodCollisions(collisions.size(), false); - std::vector dauIDList(3, -1); + std::vector dauIDList(2, -1); for (const auto& collision : collisions) { isReconstructedMCCollisions[collision.mcCollisionId()] = true; @@ -567,18 +659,25 @@ struct Hyperhelium4sigmaRecoTask { auto motherTrack = kinkCand.trackMoth_as(); auto dauTrack = kinkCand.trackDaug_as(); - bool isTrueSignal = false; + bool isKinkSignal = false; if (motherTrack.has_mcParticle() && dauTrack.has_mcParticle()) { - auto mcMotherTrack = motherTrack.mcParticle_as(); - auto mcDauTrack = dauTrack.mcParticle_as(); - auto dChannel = getDecayChannelHe4S(mcMotherTrack, dauIDList); - if (dChannel == k2body && dauIDList[0] == mcDauTrack.globalIndex()) { - isTrueSignal = true; + auto mcMothTrack = motherTrack.mcParticle_as(); + auto mcDaugTrack = dauTrack.mcParticle_as(); + if (hypoMoth == kHypertriton) { + auto dChannel = H3LDecay::getDecayChannel(mcMothTrack, dauIDList); + if (dChannel == H3LDecay::k2bodyNeutral && dauIDList[0] == mcDaugTrack.globalIndex()) { + isKinkSignal = true; + } + } else if (hypoMoth == kHyperhelium4sigma) { + auto dChannel = He4SDecay::getDecayChannel(mcMothTrack, dauIDList); + if (dChannel == He4SDecay::k2body && dauIDList[0] == mcDaugTrack.globalIndex()) { + isKinkSignal = true; + } } } registry.fill(HIST("hCandidateCounter"), 0); - if (isTrueSignal) { + if (isKinkSignal) { registry.fill(HIST("hTrueCandidateCounter"), 0); } auto collision = kinkCand.collision_as(); @@ -586,25 +685,26 @@ struct Hyperhelium4sigmaRecoTask { continue; } registry.fill(HIST("hCandidateCounter"), 1); - if (isTrueSignal) { + if (isKinkSignal) { registry.fill(HIST("hTrueCandidateCounter"), 1); } - if (std::abs(dauTrack.tpcNSigmaAl()) > cutNSigmaAl) { + float tpcNSigmaDaug = getTPCNSigma(dauTrack, pidTypeDaug); + if (std::abs(tpcNSigmaDaug) > cutNSigmaDaug) { continue; } - float invMass = RecoDecay::m(std::array{std::array{kinkCand.pxDaug(), kinkCand.pyDaug(), kinkCand.pzDaug()}, std::array{kinkCand.pxDaugNeut(), kinkCand.pyDaugNeut(), kinkCand.pzDaugNeut()}}, std::array{o2::constants::physics::MassAlpha, o2::constants::physics::MassPi0}); + float invMass = RecoDecay::m(std::array{std::array{kinkCand.pxDaug(), kinkCand.pyDaug(), kinkCand.pzDaug()}, std::array{kinkCand.pxDaugNeut(), kinkCand.pyDaugNeut(), kinkCand.pzDaugNeut()}}, std::array{massChargedDaug, massNeutralDaug}); registry.fill(HIST("hCandidateCounter"), 2); - if (isTrueSignal) { + if (isKinkSignal) { registry.fill(HIST("hTrueCandidateCounter"), 2); } - registry.fill(HIST("h2MassHyperhelium4sigmaPt"), kinkCand.mothSign() * kinkCand.ptMoth(), invMass); - registry.fill(HIST("h2NSigmaAlPt"), kinkCand.mothSign() * kinkCand.ptDaug(), dauTrack.tpcNSigmaAl()); + registry.fill(HIST("h2MothMassPt"), kinkCand.mothSign() * kinkCand.ptMoth(), invMass); + registry.fill(HIST("h2DaugTPCNSigmaPt"), kinkCand.mothSign() * kinkCand.ptDaug(), tpcNSigmaDaug); auto bc = collision.bc_as(); initCCDB(bc); - Hyphe4sCandidate hyphe4sCand; - fillCandidate(hyphe4sCand, collision, kinkCand, motherTrack, dauTrack); + HypKinkCandidate hypkinkCand; + fillCandidate(hypkinkCand, collision, kinkCand, motherTrack, dauTrack); std::array posDecVtx = {kinkCand.xDecVtx() + collision.posX(), kinkCand.yDecVtx() + collision.posY(), kinkCand.zDecVtx() + collision.posZ()}; @@ -619,79 +719,74 @@ struct Hyperhelium4sigmaRecoTask { motherTrackPar.getPxPyPzGlo(refitPPV); } - // qa for true signal - if (isTrueSignal) { - auto mcMotherTrack = motherTrack.mcParticle_as(); - auto mcDauTrack = dauTrack.mcParticle_as(); - auto mcNeutTrack = particlesMC.rawIteratorAt(dauIDList[2]); + // QA, store mcInfo for true signals + if (isKinkSignal) { + auto mcMothTrack = motherTrack.mcParticle_as(); + auto mcDaugTrack = dauTrack.mcParticle_as(); + auto mcNeutTrack = particlesMC.rawIteratorAt(dauIDList[1]); float recSVR = std::sqrt(posDecVtx[0] * posDecVtx[0] + posDecVtx[1] * posDecVtx[1]); - registry.fill(HIST("hDiffSVx"), posDecVtx[0] - mcDauTrack.vx()); - registry.fill(HIST("hDiffSVy"), posDecVtx[1] - mcDauTrack.vy()); - registry.fill(HIST("hDiffSVz"), posDecVtx[2] - mcDauTrack.vz()); - registry.fill(HIST("h2RecSVRVsTrueSVR"), recSVR, std::hypot(mcDauTrack.vx(), mcDauTrack.vy())); - registry.fill(HIST("h2TrueMotherDiffPtVsRecSVR"), recSVR, mcMotherTrack.pt() - kinkCand.ptMoth()); - registry.fill(HIST("h2TrueMotherDiffPzVsRecSVR"), recSVR, mcMotherTrack.pz() - kinkCand.pzMoth()); - registry.fill(HIST("h2TrueMotherDiffTglVsRecSVR"), recSVR, mcMotherTrack.pz() / mcMotherTrack.pt() - motherTrack.tgl()); - registry.fill(HIST("h2TrueMotherDiffEtaVsRecSVR"), recSVR, mcMotherTrack.eta() - motherTrack.eta()); - registry.fill(HIST("hDiffDauPx"), kinkCand.pxDaug() - mcDauTrack.px()); - registry.fill(HIST("hDiffDauPy"), kinkCand.pyDaug() - mcDauTrack.py()); - registry.fill(HIST("hDiffDauPz"), kinkCand.pzDaug() - mcDauTrack.pz()); + registry.fill(HIST("hDiffSVx"), posDecVtx[0] - mcDaugTrack.vx()); + registry.fill(HIST("hDiffSVy"), posDecVtx[1] - mcDaugTrack.vy()); + registry.fill(HIST("hDiffSVz"), posDecVtx[2] - mcDaugTrack.vz()); + registry.fill(HIST("h2RecSVRVsTrueSVR"), recSVR, std::hypot(mcDaugTrack.vx(), mcDaugTrack.vy())); + registry.fill(HIST("h2TrueMotherDiffPtVsRecSVR"), recSVR, mcMothTrack.pt() - kinkCand.ptMoth()); + registry.fill(HIST("h2TrueMotherDiffEtaVsRecSVR"), recSVR, mcMothTrack.eta() - motherTrack.eta()); + registry.fill(HIST("hDiffDauPx"), kinkCand.pxDaug() - mcDaugTrack.px()); + registry.fill(HIST("hDiffDauPy"), kinkCand.pyDaug() - mcDaugTrack.py()); + registry.fill(HIST("hDiffDauPz"), kinkCand.pzDaug() - mcDaugTrack.pz()); registry.fill(HIST("h2TrueSignalMassPt"), kinkCand.mothSign() * kinkCand.ptMoth(), invMass); - registry.fill(HIST("h2TrueSignalNSigmaAlPt"), kinkCand.mothSign() * kinkCand.ptDaug(), dauTrack.tpcNSigmaAl()); + registry.fill(HIST("h2TrueDaugTPCNSigmaPt"), kinkCand.mothSign() * kinkCand.ptDaug(), tpcNSigmaDaug); - hyphe4sCand.isSignal = true; - hyphe4sCand.isSignalReco = true; - hyphe4sCand.isCollReco = true; - hyphe4sCand.isSurvEvSelection = true; - fillCandidateMCInfo(hyphe4sCand, mcMotherTrack, mcDauTrack, mcNeutTrack); - mcHe4sIndices.push_back(mcMotherTrack.globalIndex()); + hypkinkCand.isSignal = true; + hypkinkCand.isSignalReco = true; + hypkinkCand.isCollReco = true; + hypkinkCand.isSurvEvSelection = true; + fillCandidateMCInfo(hypkinkCand, mcMothTrack, mcDaugTrack, mcNeutTrack); + mcPartIndices.push_back(mcMothTrack.globalIndex()); std::array dcaInfo; - auto mcMotherTrackPar = getTrackParFromMC(mcMotherTrack, 2); - o2::base::Propagator::Instance()->propagateToDCABxByBz({posDecVtx[0], posDecVtx[1], posDecVtx[2]}, mcMotherTrackPar, 2.f, matCorr, &dcaInfo); + auto mcMothTrackPar = getTrackParFromMC(mcMothTrack, 2); + o2::base::Propagator::Instance()->propagateToDCABxByBz({posDecVtx[0], posDecVtx[1], posDecVtx[2]}, mcMothTrackPar, 2.f, matCorr, &dcaInfo); registry.fill(HIST("hDCAXYMothToRecSV"), dcaInfo[0]); registry.fill(HIST("hDCAZMothToRecSV"), dcaInfo[1]); - std::array mcPMotherSV = {-999.f, -999.f, -999.f}; - mcMotherTrackPar.getPxPyPzGlo(mcPMotherSV); - registry.fill(HIST("h2TrueMotherDiffPxVsRecSVR"), recSVR, mcPMotherSV[0] - kinkCand.pxMoth()); - registry.fill(HIST("h2TrueMotherDiffPyVsRecSVR"), recSVR, mcPMotherSV[1] - kinkCand.pyMoth()); } outputMCTable( - hyphe4sCand.posPV[0], hyphe4sCand.posPV[1], hyphe4sCand.posPV[2], - hyphe4sCand.posSV[0], hyphe4sCand.posSV[1], hyphe4sCand.posSV[2], - hyphe4sCand.isMatter, - hyphe4sCand.lastPosMoth[0], hyphe4sCand.lastPosMoth[1], hyphe4sCand.lastPosMoth[2], - hyphe4sCand.momMothSV[0], hyphe4sCand.momMothSV[1], hyphe4sCand.momMothSV[2], + mBz > 0 ? 1 : -1, + hypkinkCand.posPV[0], hypkinkCand.posPV[1], hypkinkCand.posPV[2], + hypkinkCand.posSV[0], hypkinkCand.posSV[1], hypkinkCand.posSV[2], + hypkinkCand.isMatter, + hypkinkCand.lastPosMoth[0], hypkinkCand.lastPosMoth[1], hypkinkCand.lastPosMoth[2], + hypkinkCand.momMothSV[0], hypkinkCand.momMothSV[1], hypkinkCand.momMothSV[2], refitPPV[0], refitPPV[1], refitPPV[2], refitPSV[0], refitPSV[1], refitPSV[2], - hyphe4sCand.momDaugSV[0], hyphe4sCand.momDaugSV[1], hyphe4sCand.momDaugSV[2], - hyphe4sCand.dcaXYMothPv, hyphe4sCand.dcaXYDaugPv, hyphe4sCand.dcaKinkTopo, - hyphe4sCand.chi2ITSMoth, hyphe4sCand.itsClusterSizeMoth, hyphe4sCand.itsClusterSizeDaug, - hyphe4sCand.nSigmaTPCDaug, hyphe4sCand.nSigmaITSDaug, - hyphe4sCand.isSignal, hyphe4sCand.isSignalReco, hyphe4sCand.isCollReco, hyphe4sCand.isSurvEvSelection, - hyphe4sCand.truePosSV[0], hyphe4sCand.truePosSV[1], hyphe4sCand.truePosSV[2], - hyphe4sCand.trueMomMothPV[0], hyphe4sCand.trueMomMothPV[1], hyphe4sCand.trueMomMothPV[2], - hyphe4sCand.trueMomMothSV[0], hyphe4sCand.trueMomMothSV[1], hyphe4sCand.trueMomMothSV[2], - hyphe4sCand.trueMomDaugSV[0], hyphe4sCand.trueMomDaugSV[1], hyphe4sCand.trueMomDaugSV[2], - hyphe4sCand.isMothReco, hyphe4sCand.momMothPV[0], hyphe4sCand.momMothPV[1], hyphe4sCand.momMothPV[2], - hyphe4sCand.updateMomMothPV[0], hyphe4sCand.updateMomMothPV[1], hyphe4sCand.updateMomMothPV[2]); + hypkinkCand.momDaugSV[0], hypkinkCand.momDaugSV[1], hypkinkCand.momDaugSV[2], + hypkinkCand.dcaXYMothPv, hypkinkCand.dcaXYDaugPv, hypkinkCand.dcaKinkTopo, + hypkinkCand.chi2ITSMoth, hypkinkCand.itsClusterSizeMoth, hypkinkCand.itsClusterSizeDaug, + hypkinkCand.nSigmaTPCDaug, hypkinkCand.nSigmaITSDaug, + hypkinkCand.isSignal, hypkinkCand.isSignalReco, hypkinkCand.isCollReco, hypkinkCand.isSurvEvSelection, + hypkinkCand.truePosSV[0], hypkinkCand.truePosSV[1], hypkinkCand.truePosSV[2], + hypkinkCand.trueMomMothPV[0], hypkinkCand.trueMomMothPV[1], hypkinkCand.trueMomMothPV[2], + hypkinkCand.trueMomMothSV[0], hypkinkCand.trueMomMothSV[1], hypkinkCand.trueMomMothSV[2], + hypkinkCand.trueMomDaugSV[0], hypkinkCand.trueMomDaugSV[1], hypkinkCand.trueMomDaugSV[2], + hypkinkCand.isMothReco, hypkinkCand.momMothPV[0], hypkinkCand.momMothPV[1], hypkinkCand.momMothPV[2], + hypkinkCand.updateMomMothPV[0], hypkinkCand.updateMomMothPV[1], hypkinkCand.updateMomMothPV[2]); } - // fill hyperhelium4sigma signals which are not reconstructed + // fill kink signals which are not reconstructed for (auto const& mcparticle : particlesMC) { - auto dChannel = getDecayChannelHe4S(mcparticle, dauIDList); - if (dChannel != k2body) { + auto dChannel = He4SDecay::getDecayChannel(mcparticle, dauIDList); + if (dChannel != He4SDecay::k2body) { continue; } - if (std::find(mcHe4sIndices.begin(), mcHe4sIndices.end(), mcparticle.globalIndex()) != mcHe4sIndices.end()) { + if (std::find(mcPartIndices.begin(), mcPartIndices.end(), mcparticle.globalIndex()) != mcPartIndices.end()) { continue; } - Hyphe4sCandidate hyphe4sCand; - auto mcDauTrack = particlesMC.rawIteratorAt(dauIDList[0]); - auto mcNeutTrack = particlesMC.rawIteratorAt(dauIDList[2]); - fillCandidateMCInfo(hyphe4sCand, mcparticle, mcDauTrack, mcNeutTrack); + HypKinkCandidate hypkinkCand; + auto mcDaugTrack = particlesMC.rawIteratorAt(dauIDList[0]); + auto mcNeutTrack = particlesMC.rawIteratorAt(dauIDList[1]); + fillCandidateMCInfo(hypkinkCand, mcparticle, mcDaugTrack, mcNeutTrack); if (mcPartIndices[mcparticle.globalIndex()] != -1) { auto mothTrack = tracks.rawIteratorAt(mcPartIndices[mcparticle.globalIndex()]); @@ -699,11 +794,12 @@ struct Hyperhelium4sigmaRecoTask { auto collision = mothTrack.collision_as(); auto bc = collision.template bc_as(); initCCDB(bc); - fillCandidateRecoMoth(hyphe4sCand, collision, mothTrack); + fillCandidateRecoMoth(hypkinkCand, collision, mothTrack); } } outputMCTable( + mBz > 0 ? 1 : -1, -1, -1, -1, -1, -1, -1, -1, @@ -716,20 +812,22 @@ struct Hyperhelium4sigmaRecoTask { -1, -1, -1, -1, -1, true, false, isReconstructedMCCollisions[mcparticle.mcCollisionId()], isSelectedMCCollisions[mcparticle.mcCollisionId()], - hyphe4sCand.truePosSV[0], hyphe4sCand.truePosSV[1], hyphe4sCand.truePosSV[2], - hyphe4sCand.trueMomMothPV[0], hyphe4sCand.trueMomMothPV[1], hyphe4sCand.trueMomMothPV[2], - hyphe4sCand.trueMomMothSV[0], hyphe4sCand.trueMomMothSV[1], hyphe4sCand.trueMomMothSV[2], - hyphe4sCand.trueMomDaugSV[0], hyphe4sCand.trueMomDaugSV[1], hyphe4sCand.trueMomDaugSV[2], - hyphe4sCand.isMothReco, hyphe4sCand.momMothPV[0], hyphe4sCand.momMothPV[1], hyphe4sCand.momMothPV[2], - hyphe4sCand.updateMomMothPV[0], hyphe4sCand.updateMomMothPV[1], hyphe4sCand.updateMomMothPV[2]); + hypkinkCand.truePosSV[0], hypkinkCand.truePosSV[1], hypkinkCand.truePosSV[2], + hypkinkCand.trueMomMothPV[0], hypkinkCand.trueMomMothPV[1], hypkinkCand.trueMomMothPV[2], + hypkinkCand.trueMomMothSV[0], hypkinkCand.trueMomMothSV[1], hypkinkCand.trueMomMothSV[2], + hypkinkCand.trueMomDaugSV[0], hypkinkCand.trueMomDaugSV[1], hypkinkCand.trueMomDaugSV[2], + hypkinkCand.isMothReco, hypkinkCand.momMothPV[0], hypkinkCand.momMothPV[1], hypkinkCand.momMothPV[2], + hypkinkCand.updateMomMothPV[0], hypkinkCand.updateMomMothPV[1], hypkinkCand.updateMomMothPV[2]); } } - PROCESS_SWITCH(Hyperhelium4sigmaRecoTask, processMC, "process MC", false); + PROCESS_SWITCH(HyperkinkRecoTask, processMC, "process MC", false); }; //-------------------------------------------------------------- // check the performance of mcparticle -struct Hyperhelium4sigmaQa { +struct HyperkinkQa { + + Configurable hypoMoth{"hypoMoth", kHypertriton, "Mother particle hypothesis"}; HistogramRegistry genQAHist{"genQAHist", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; HistogramRegistry recoQAHist{"recoQAHist", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; @@ -743,11 +841,38 @@ struct Hyperhelium4sigmaQa { o2::aod::ITSResponse itsResponse; + float massMoth = 999.f; + float massChargedDaug = 999.f; + float massNeutralDaug = 999.f; + int pdgMoth = 0; + std::array pdgDaug = {o2::constants::physics::Pdg::kTriton, PDG_t::kPi0}; // pdgcode of charged (0) and neutral (1) daughter particles + o2::track::PID pidTypeDaug = o2::track::PID::Triton; + void init(InitContext&) { if (doprocessMC == true) { itsResponse.setMCDefaultParameters(); + if (hypoMoth == kHypertriton) { + massMoth = o2::constants::physics::MassHyperTriton; + massChargedDaug = o2::constants::physics::MassTriton; + massNeutralDaug = o2::constants::physics::MassPi0; + pdgMoth = o2::constants::physics::Pdg::kHyperTriton; + pdgDaug[kDaugCharged] = o2::constants::physics::Pdg::kTriton; + pdgDaug[kDaugNeutral] = PDG_t::kPi0; + pidTypeDaug = o2::track::PID::Triton; + } else if (hypoMoth == kHyperhelium4sigma) { + massMoth = o2::constants::physics::MassHyperHelium4Sigma; + massChargedDaug = o2::constants::physics::MassAlpha; + massNeutralDaug = o2::constants::physics::MassPi0; + pdgMoth = o2::constants::physics::Pdg::kHyperHelium4Sigma; + pdgDaug[kDaugCharged] = o2::constants::physics::Pdg::kAlpha; + pdgDaug[kDaugNeutral] = PDG_t::kPi0; + pidTypeDaug = o2::track::PID::Alpha; + } else { + LOG(fatal) << "Unknown mother particle hypothesis"; + } + const AxisSpec pAxis{ptBins, "#it{p} (GeV/#it{c})"}; const AxisSpec ptAxis{ptBins, "#it{p}_{T} (GeV/#it{c})"}; const AxisSpec ctAxis{ctBins, "c#it{t} (cm)"}; @@ -767,62 +892,62 @@ struct Hyperhelium4sigmaQa { hMcCollCounter->GetXaxis()->SetBinLabel(1, "MC Collisions"); hMcCollCounter->GetXaxis()->SetBinLabel(2, "Reconstructed"); - auto hGenHyperHelium4SigmaCounter = genQAHist.add("hGenHyperHelium4SigmaCounter", "", HistType::kTH1F, {{11, 0.f, 11.f}}); - hGenHyperHelium4SigmaCounter->GetXaxis()->SetBinLabel(1, "He4S All"); - hGenHyperHelium4SigmaCounter->GetXaxis()->SetBinLabel(2, "Matter"); - hGenHyperHelium4SigmaCounter->GetXaxis()->SetBinLabel(3, "AntiMatter"); - hGenHyperHelium4SigmaCounter->GetXaxis()->SetBinLabel(4, "#alpha + #pi^{0}"); - hGenHyperHelium4SigmaCounter->GetXaxis()->SetBinLabel(5, "#bar{#alpha} + #pi^{0}"); - hGenHyperHelium4SigmaCounter->GetXaxis()->SetBinLabel(6, "t + p + #pi^{0}"); - hGenHyperHelium4SigmaCounter->GetXaxis()->SetBinLabel(7, "#bar{t} + #bar{p} + #pi^{0}"); - hGenHyperHelium4SigmaCounter->GetXaxis()->SetBinLabel(8, "t + n + #pi^{+}"); - hGenHyperHelium4SigmaCounter->GetXaxis()->SetBinLabel(9, "#bar{t} + #bar{n} + #pi^{+}"); - hGenHyperHelium4SigmaCounter->GetXaxis()->SetBinLabel(10, "Tracks found"); - hGenHyperHelium4SigmaCounter->GetXaxis()->SetBinLabel(11, "Unexpected"); - - auto hEvtSelectedHyperHelium4SigmaCounter = genQAHist.add("hEvtSelectedHyperHelium4SigmaCounter", "", HistType::kTH1F, {{2, 0.f, 2.f}}); - hEvtSelectedHyperHelium4SigmaCounter->GetXaxis()->SetBinLabel(1, "Generated"); - hEvtSelectedHyperHelium4SigmaCounter->GetXaxis()->SetBinLabel(2, "Survived"); - - genQAHist.add("hGenHyperHelium4SigmaP", "", HistType::kTH1F, {pAxis}); - genQAHist.add("hGenHyperHelium4SigmaPt", "", HistType::kTH1F, {ptAxis}); - genQAHist.add("hGenHyperHelium4SigmaCt", "", HistType::kTH1F, {ctAxis}); + if (hypoMoth == kHypertriton) { + auto hGenHyperMothCounter = genQAHist.add("hGenHyperMothCounter", "", HistType::kTH1F, {{10, 0.f, 10.f}}); + hGenHyperMothCounter->GetXaxis()->SetBinLabel(1, "H3L All"); + hGenHyperMothCounter->GetXaxis()->SetBinLabel(2, "Matter"); + hGenHyperMothCounter->GetXaxis()->SetBinLabel(3, "AntiMatter"); + hGenHyperMothCounter->GetXaxis()->SetBinLabel(4, "t + #pi^{0}"); + hGenHyperMothCounter->GetXaxis()->SetBinLabel(5, "#bar{t} + #pi^{0}"); + hGenHyperMothCounter->GetXaxis()->SetBinLabel(6, "he3 + #pi^{-}"); + hGenHyperMothCounter->GetXaxis()->SetBinLabel(7, "#bar{he3} + #pi^{+}"); + hGenHyperMothCounter->GetXaxis()->SetBinLabel(8, "d + p + #pi^{-}"); + hGenHyperMothCounter->GetXaxis()->SetBinLabel(9, "#bar{d} + #bar{p} + #pi^{+}"); + hGenHyperMothCounter->GetXaxis()->SetBinLabel(10, "Others"); + } else if (hypoMoth == kHyperhelium4sigma) { + auto hGenHyperMothCounter = genQAHist.add("hGenHyperMothCounter", "", HistType::kTH1F, {{10, 0.f, 10.f}}); + hGenHyperMothCounter->GetXaxis()->SetBinLabel(1, "He4S All"); + hGenHyperMothCounter->GetXaxis()->SetBinLabel(2, "Matter"); + hGenHyperMothCounter->GetXaxis()->SetBinLabel(3, "AntiMatter"); + hGenHyperMothCounter->GetXaxis()->SetBinLabel(4, "#alpha + #pi^{0}"); + hGenHyperMothCounter->GetXaxis()->SetBinLabel(5, "#bar{#alpha} + #pi^{0}"); + hGenHyperMothCounter->GetXaxis()->SetBinLabel(6, "t + p + #pi^{0}"); + hGenHyperMothCounter->GetXaxis()->SetBinLabel(7, "#bar{t} + #bar{p} + #pi^{0}"); + hGenHyperMothCounter->GetXaxis()->SetBinLabel(8, "t + n + #pi^{+}"); + hGenHyperMothCounter->GetXaxis()->SetBinLabel(9, "#bar{t} + #bar{n} + #pi^{+}"); + hGenHyperMothCounter->GetXaxis()->SetBinLabel(10, "Others"); + } + + auto hEvtSelectedHyperMothCounter = genQAHist.add("hEvtSelectedHyperMothCounter", "", HistType::kTH1F, {{2, 0.f, 2.f}}); + hEvtSelectedHyperMothCounter->GetXaxis()->SetBinLabel(1, "Generated"); + hEvtSelectedHyperMothCounter->GetXaxis()->SetBinLabel(2, "Survived"); + + genQAHist.add("hGenHyperMothP", "", HistType::kTH1F, {pAxis}); + genQAHist.add("hGenHyperMothPt", "", HistType::kTH1F, {ptAxis}); + genQAHist.add("hGenHyperMothCt", "", HistType::kTH1F, {ctAxis}); genQAHist.add("hMcRecoInvMass", "", HistType::kTH1F, {invMassAxis}); // efficiency/criteria studies for tracks which are true candidates hMothCounter = recoQAHist.add("hMothCounter", "", HistType::kTH1F, {{9, 0.f, 9.f}}); - hMoth2BCounter = recoQAHist.add("hMoth2BCounter", "", HistType::kTH1F, {{9, 0.f, 9.f}}); - for (const auto& hist : {hMothCounter, hMoth2BCounter}) { - hist->GetXaxis()->SetBinLabel(1, "Generated"); - hist->GetXaxis()->SetBinLabel(2, "Reconstructed"); - hist->GetXaxis()->SetBinLabel(3, "eta"); - hist->GetXaxis()->SetBinLabel(4, "has collision"); - hist->GetXaxis()->SetBinLabel(5, "ITSonly"); - hist->GetXaxis()->SetBinLabel(6, "ITS hits"); - hist->GetXaxis()->SetBinLabel(7, "ITS IR"); - hist->GetXaxis()->SetBinLabel(8, "ITS chi2"); - hist->GetXaxis()->SetBinLabel(9, "pt"); - } + hMothCounter->GetXaxis()->SetBinLabel(1, "Generated"); + hMothCounter->GetXaxis()->SetBinLabel(2, "Reconstructed"); + hMothCounter->GetXaxis()->SetBinLabel(3, "eta"); + hMothCounter->GetXaxis()->SetBinLabel(4, "has collision"); + hMothCounter->GetXaxis()->SetBinLabel(5, "ITSonly"); + hMothCounter->GetXaxis()->SetBinLabel(6, "ITS hits"); + hMothCounter->GetXaxis()->SetBinLabel(7, "ITS IR"); + hMothCounter->GetXaxis()->SetBinLabel(8, "ITS chi2"); + hMothCounter->GetXaxis()->SetBinLabel(9, "pt"); recoQAHist.add("h2TrueMotherDiffPtVsTrueSVR", ";Decay Vertex R (cm);#Delta p_{T} (GeV/#it{c});", HistType::kTH2F, {svRadiuAxis, diffPtAxis}); - recoQAHist.add("h2TrueMotherDiffPzVsTrueSVR", ";Decay Vertex R (cm);#Delta p_{z} (GeV/#it{c});", HistType::kTH2F, {svRadiuAxis, diffPzAxis}); - recoQAHist.add("h2TrueMotherDiffTglVsTrueSVR", ";Decay Vertex R (cm);#Delta tgl;", HistType::kTH2F, {svRadiuAxis, {200, -1.f, 1.f}}); recoQAHist.add("h2TrueMotherDiffEtaVsTrueSVR", ";Decay Vertex R (cm);#Delta #eta;", HistType::kTH2F, {svRadiuAxis, {200, -1.f, 1.f}}); recoQAHist.add("h2GoodMotherDiffPtVsTrueSVR", ";Decay Vertex R (cm);#Delta p_{T} (GeV/#it{c});", HistType::kTH2F, {svRadiuAxis, diffPtAxis}); - recoQAHist.add("h2GoodMotherDiffPzVsTrueSVR", ";Decay Vertex R (cm);#Delta p_{z} (GeV/#it{c});", HistType::kTH2F, {svRadiuAxis, diffPzAxis}); - - hDaugCounter[kDaugAlpha] = recoQAHist.add("hDaugAlphaCounter", "", HistType::kTH1F, {{9, 0.f, 9.f}}); - hDaugCounter[kDaugTriton] = recoQAHist.add("hDaugTritonCounter", "", HistType::kTH1F, {{9, 0.f, 9.f}}); - hDaugCounter[kDaugProton] = recoQAHist.add("hDaugProtonCounter", "", HistType::kTH1F, {{9, 0.f, 9.f}}); - hDaugCounter[kDaugChargedPion] = recoQAHist.add("hDaugPionCounter", "", HistType::kTH1F, {{9, 0.f, 9.f}}); - hDaugTPCNSigma[kDaugAlpha] = recoQAHist.add("hDaugAlphaTPCNSigma", "", HistType::kTH2F, {rigidityAxis, nsigmaAxis}); - hDaugTPCNSigma[kDaugTriton] = recoQAHist.add("hDaugTritonTPCNSigma", "", HistType::kTH2F, {rigidityAxis, nsigmaAxis}); - hDaugTPCNSigma[kDaugProton] = recoQAHist.add("hDaugProtonTPCNSigma", "", HistType::kTH2F, {rigidityAxis, nsigmaAxis}); - hDaugTPCNSigma[kDaugChargedPion] = recoQAHist.add("hDaugPionTPCNSigma", "", HistType::kTH2F, {rigidityAxis, nsigmaAxis}); - recoQAHist.add("hDaugAlphaITSNSigmaCheck", "", HistType::kTH2F, {rigidityAxis, itsnsigmaAxis}); + + hDaugCounter = recoQAHist.add("hDaugCounter", "", HistType::kTH1F, {{9, 0.f, 9.f}}); + hDaugTPCNSigma = recoQAHist.add("hDaugTPCNSigma", "", HistType::kTH2F, {rigidityAxis, nsigmaAxis}); hRecoMothCounter = recoQAHist.add("hRecoMothCounter", "", HistType::kTH1F, {{9, 0.f, 9.f}}); - hRecoDaugAlphaCounter = recoQAHist.add("hRecoDaugAlphaCounter", "", HistType::kTH1F, {{9, 0.f, 9.f}}); - for (const auto& hist : {hDaugCounter[kDaugAlpha], hDaugCounter[kDaugTriton], hDaugCounter[kDaugProton], hDaugCounter[kDaugChargedPion], hRecoMothCounter, hRecoDaugAlphaCounter}) { + hRecoDaugCounter = recoQAHist.add("hRecoDaugCounter", "", HistType::kTH1F, {{9, 0.f, 9.f}}); + for (const auto& hist : {hDaugCounter, hRecoMothCounter, hRecoDaugCounter}) { hist->GetXaxis()->SetBinLabel(1, "Generated"); hist->GetXaxis()->SetBinLabel(2, "Reconstructed"); hist->GetXaxis()->SetBinLabel(3, "eta"); @@ -836,21 +961,20 @@ struct Hyperhelium4sigmaQa { recoQAHist.add("hMothIsPVContributer", "", HistType::kTH1F, {{2, 0.f, 2.f}}); recoQAHist.add("hMothITSCls", "", HistType::kTH1F, {{8, 0.f, 8.f}}); - recoQAHist.add("hDaugAlphaIsPVContributer", "", HistType::kTH1F, {{2, 0.f, 2.f}}); - recoQAHist.add("hDaugAlphaITSCls", "", HistType::kTH1F, {{8, 0.f, 8.f}}); - recoQAHist.add("hDaugAlphaITSNSigma", "", HistType::kTH2F, {rigidityAxis, itsnsigmaAxis}); - recoQAHist.add("hReco2BDauAlphaPVsITSNSigma", "", HistType::kTH2F, {rigidityAxis, itsnsigmaAxis}); - recoQAHist.add("hReco2BCandidateCount", "", HistType::kTH1F, {{4, 0.f, 4.f}}); + recoQAHist.add("hDaugIsPVContributer", "", HistType::kTH1F, {{2, 0.f, 2.f}}); + recoQAHist.add("hDaugITSCls", "", HistType::kTH1F, {{8, 0.f, 8.f}}); + recoQAHist.add("hDaugITSNSigma", "", HistType::kTH2F, {rigidityAxis, itsnsigmaAxis}); + recoQAHist.add("hRecoDaugPVsITSNSigma", "", HistType::kTH2F, {rigidityAxis, itsnsigmaAxis}); + recoQAHist.add("hRecoCandidateCount", "", HistType::kTH1F, {{4, 0.f, 4.f}}); } } Configurable skipRejectedEvents{"skipRejectedEvents", false, "Flag to skip events that fail event selection cuts"}; Configurable doEventCut{"doEventCut", true, "Apply event selection"}; Configurable maxZVertex{"maxZVertex", 10.0f, "Accepted z-vertex range (cm)"}; - Configurable only2BodyDecay{"only2BodyDecay", true, "Only consider 2-body decays for hyperhelium4sigma"}; Configurable etaMax{"etaMax", 1., "eta cut for tracks"}; - Configurable minPtMoth{"minPtMoth", 0.5, "Minimum pT/z of the hyperhelium4sigma candidate"}; + Configurable minPtMoth{"minPtMoth", 0.5, "Minimum pT/z of the mother track"}; Configurable tpcPidNsigmaCut{"tpcPidNsigmaCut", 5, "tpcPidNsigmaCut"}; Configurable nTPCClusMinDaug{"nTPCClusMinDaug", 80, "daug NTPC clusters cut"}; Configurable itsMaxChi2{"itsMaxChi2", 36, "max chi2 for ITS"}; @@ -858,7 +982,7 @@ struct Hyperhelium4sigmaQa { Preslice permcCollision = o2::aod::mcparticle::mcCollisionId; - // qa for mother track selection + // QA for mother track selection template bool motherTrackCheck(const TTrack& track, const std::shared_ptr hist) { @@ -950,14 +1074,14 @@ struct Hyperhelium4sigmaQa { { // dummy process function; } - PROCESS_SWITCH(Hyperhelium4sigmaQa, processData, "process data", true); + PROCESS_SWITCH(HyperkinkQa, processData, "process data", true); void processMC(aod::McCollisions const& mcCollisions, aod::McParticles const& particlesMC, MCLabeledCollisionsFull const& collisions, MCLabeledTracksIU const& tracks, aod::BCs const&) { std::vector mcPartIndices; setTrackIDForMC(mcPartIndices, particlesMC, tracks); std::vector isSelectedMCCollisions(mcCollisions.size(), false); - std::vector dauIDList(3, -1); + std::vector dauIDList(2, -1); for (const auto& collision : collisions) { genQAHist.fill(HIST("hCollCounter"), 0.5); if (doEventCut && (!collision.selection_bit(aod::evsel::kIsTriggerTVX) || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder) || std::abs(collision.posZ()) > maxZVertex)) { @@ -980,66 +1104,75 @@ struct Hyperhelium4sigmaQa { const auto& dparticlesMC = particlesMC.sliceBy(permcCollision, mcCollision.globalIndex()); for (const auto& mcparticle : dparticlesMC) { - - bool isMatter; - if (mcparticle.pdgCode() == o2::constants::physics::Pdg::kHyperHelium4Sigma) { - isMatter = true; - } else if (mcparticle.pdgCode() == -o2::constants::physics::Pdg::kHyperHelium4Sigma) { - isMatter = false; - } else { + if (std::abs(mcparticle.pdgCode()) != pdgMoth) { continue; } - - auto dChannel = getDecayChannelHe4S(mcparticle, dauIDList); - if (dChannel == kNDecayChannel) { - genQAHist.fill(HIST("hGenHyperHelium4SigmaCounter"), 10.5); - continue; + bool isMatter = mcparticle.pdgCode() > 0; + genQAHist.fill(HIST("hGenHyperMothCounter"), 0.5); + genQAHist.fill(HIST("hGenHyperMothCounter"), isMatter ? 1.5 : 2.5); + + // QA for decay channels + bool isKinkSignal = false; + if (hypoMoth == kHypertriton) { + auto dChannel = H3LDecay::getDecayChannel(mcparticle, dauIDList); + if (dChannel == H3LDecay::k2bodyNeutral) { + genQAHist.fill(HIST("hGenHyperMothCounter"), isMatter ? 3.5 : 4.5); + isKinkSignal = true; + } else if (dChannel == H3LDecay::k2bodyCharged) { + genQAHist.fill(HIST("hGenHyperMothCounter"), isMatter ? 5.5 : 6.5); + } else if (dChannel == H3LDecay::k3bodyCharged) { + genQAHist.fill(HIST("hGenHyperMothCounter"), isMatter ? 7.5 : 8.5); + } else if (dChannel == H3LDecay::kNChannel) { + genQAHist.fill(HIST("hGenHyperMothCounter"), 9.5); + continue; + } + } else if (hypoMoth == kHyperhelium4sigma) { + auto dChannel = He4SDecay::getDecayChannel(mcparticle, dauIDList); + if (dChannel == He4SDecay::k2body) { + genQAHist.fill(HIST("hGenHyperMothCounter"), isMatter ? 3.5 : 4.5); + isKinkSignal = true; + } else if (dChannel == He4SDecay::k3body_p) { + genQAHist.fill(HIST("hGenHyperMothCounter"), isMatter ? 5.5 : 6.5); + } else if (dChannel == He4SDecay::k3body_n) { + genQAHist.fill(HIST("hGenHyperMothCounter"), isMatter ? 7.5 : 8.5); + } else if (dChannel == He4SDecay::kNChannel) { + genQAHist.fill(HIST("hGenHyperMothCounter"), 9.5); + continue; + } } - // qa for mother tracks - if (dChannel == k2body) { - recoQAHist.fill(HIST("hMoth2BCounter"), 0); - } else { - if (only2BodyDecay) { - continue; // skip 3-body decays - } + if (!isKinkSignal) { + continue; } recoQAHist.fill(HIST("hMothCounter"), 0); - genQAHist.fill(HIST("hGenHyperHelium4SigmaCounter"), 0.5); - genQAHist.fill(HIST("hGenHyperHelium4SigmaCounter"), isMatter ? 1.5 : 2.5); - genQAHist.fill(HIST("hEvtSelectedHyperHelium4SigmaCounter"), 0.5); + genQAHist.fill(HIST("hEvtSelectedHyperMothCounter"), 0.5); if (isSelectedMCCollisions[mcCollision.globalIndex()]) { - genQAHist.fill(HIST("hEvtSelectedHyperHelium4SigmaCounter"), 1.5); + genQAHist.fill(HIST("hEvtSelectedHyperMothCounter"), 1.5); } float svPos[3] = {-999, -999, -999}; std::vector> dauMom(kNDaughterType, std::vector(3, -999.0f)); for (const auto& mcparticleDaughter : mcparticle.daughters_as()) { for (int type = 0; type < kNDaughterType; type++) { - if (std::abs(mcparticleDaughter.pdgCode()) == kDaugghterPDG[type]) { + if (std::abs(mcparticleDaughter.pdgCode()) == pdgDaug[type]) { dauMom[type][0] = mcparticleDaughter.px(); dauMom[type][1] = mcparticleDaughter.py(); dauMom[type][2] = mcparticleDaughter.pz(); - if (type <= kDaugTriton) { + if (type == kDaugCharged) { svPos[0] = mcparticleDaughter.vx(); svPos[1] = mcparticleDaughter.vy(); svPos[2] = mcparticleDaughter.vz(); - } - if (type < kNChargedDaughterType) { - hDaugCounter[type]->Fill(0.f); // if daughter track is reconstructed - if (type <= kNChargedDaughterType && mcPartIndices[mcparticleDaughter.globalIndex()] != -1) { + if (mcPartIndices[mcparticleDaughter.globalIndex()] != -1) { + hDaugCounter->Fill(0.f); auto track = tracks.rawIteratorAt(mcPartIndices[mcparticleDaughter.globalIndex()]); - float tpcNSigma = getTPCNSigma(track, type); - daughterTrackCheck(track, hDaugCounter[type], tpcNSigma); + float tpcNSigma = getTPCNSigma(track, pidTypeDaug); + daughterTrackCheck(track, hDaugCounter, tpcNSigma); if (track.hasTPC()) { - hDaugTPCNSigma[type]->Fill(track.p() * track.sign(), tpcNSigma); - } - if (type == kDaugAlpha && track.itsNCls() > kITSLayers - 2) { - recoQAHist.fill(HIST("hDaugAlphaITSNSigmaCheck"), track.p() * track.sign(), itsResponse.nSigmaITS(track)); + hDaugTPCNSigma->Fill(track.p() * track.sign(), tpcNSigma); } } } @@ -1047,94 +1180,69 @@ struct Hyperhelium4sigmaQa { } } - genQAHist.fill(HIST("hGenHyperHelium4SigmaP"), mcparticle.p()); - genQAHist.fill(HIST("hGenHyperHelium4SigmaPt"), mcparticle.pt()); - float ct = RecoDecay::sqrtSumOfSquares(svPos[0] - mcparticle.vx(), svPos[1] - mcparticle.vy(), svPos[2] - mcparticle.vz()) * o2::constants::physics::MassHyperHelium4Sigma / mcparticle.p(); - genQAHist.fill(HIST("hGenHyperHelium4SigmaCt"), ct); + genQAHist.fill(HIST("hGenHyperMothP"), mcparticle.p()); + genQAHist.fill(HIST("hGenHyperMothPt"), mcparticle.pt()); + float ct = RecoDecay::sqrtSumOfSquares(svPos[0] - mcparticle.vx(), svPos[1] - mcparticle.vy(), svPos[2] - mcparticle.vz()) * massMoth / mcparticle.p(); + genQAHist.fill(HIST("hGenHyperMothCt"), ct); + float hypermothMCMass = RecoDecay::m(std::array{std::array{dauMom[kDaugCharged][0], dauMom[kDaugCharged][1], dauMom[kDaugCharged][2]}, std::array{dauMom[kDaugNeutral][0], dauMom[kDaugNeutral][1], dauMom[kDaugNeutral][2]}}, std::array{massChargedDaug, massNeutralDaug}); + genQAHist.fill(HIST("hMcRecoInvMass"), hypermothMCMass); // if mother track is reconstructed if (mcPartIndices[mcparticle.globalIndex()] != -1) { auto motherTrack = tracks.rawIteratorAt(mcPartIndices[mcparticle.globalIndex()]); bool isGoodMother = motherTrackCheck(motherTrack, hMothCounter); - if (dChannel == k2body) { - motherTrackCheck(motherTrack, hMoth2BCounter); - } float svR = RecoDecay::sqrtSumOfSquares(svPos[0], svPos[1]); float diffpt = mcparticle.pt() - 2 * motherTrack.pt(); - float diffpz = mcparticle.pz() - 2 * motherTrack.pz(); recoQAHist.fill(HIST("h2TrueMotherDiffPtVsTrueSVR"), svR, diffpt); - recoQAHist.fill(HIST("h2TrueMotherDiffPzVsTrueSVR"), svR, diffpz); - recoQAHist.fill(HIST("h2TrueMotherDiffTglVsTrueSVR"), svR, mcparticle.pz() / mcparticle.pt() - motherTrack.tgl()); recoQAHist.fill(HIST("h2TrueMotherDiffEtaVsTrueSVR"), svR, mcparticle.eta() - motherTrack.eta()); if (isGoodMother) { recoQAHist.fill(HIST("h2GoodMotherDiffPtVsTrueSVR"), svR, diffpt); - recoQAHist.fill(HIST("h2GoodMotherDiffPzVsTrueSVR"), svR, diffpz); } // if mother track and charged daughters are all reconstructed - bool isDauReconstructed = mcPartIndices[dauIDList[0]] != -1 && (dChannel == k2body ? true : mcPartIndices[dauIDList[1]] != -1); + bool isDauReconstructed = mcPartIndices[dauIDList[0]] != -1; if (isDauReconstructed) { - genQAHist.fill(HIST("hGenHyperHelium4SigmaCounter"), 9.5); - - // qa for bc matching for reconstructed tracks - if (dChannel == k2body) { - auto daughterTrack = tracks.rawIteratorAt(mcPartIndices[dauIDList[0]]); - bool isMoth = motherTrackCheck(motherTrack, hRecoMothCounter); - bool isDaug = daughterTrackCheck(daughterTrack, hRecoDaugAlphaCounter, daughterTrack.tpcNSigmaAl()); - - recoQAHist.fill(HIST("hReco2BCandidateCount"), 0.5); - recoQAHist.fill(HIST("hRecoMothCounter"), 0.5); - recoQAHist.fill(HIST("hMothITSCls"), motherTrack.itsNCls()); - recoQAHist.fill(HIST("hRecoDaugAlphaCounter"), 0.5); - recoQAHist.fill(HIST("hMothIsPVContributer"), motherTrack.isPVContributor() ? 1.5 : 0.5); - recoQAHist.fill(HIST("hDaugAlphaIsPVContributer"), daughterTrack.isPVContributor() ? 1.5 : 0.5); - - float itsNSigma = itsResponse.nSigmaITS(daughterTrack); - if (daughterTrack.hasITS()) { - recoQAHist.fill(HIST("hDaugAlphaITSNSigma"), daughterTrack.sign() * daughterTrack.p(), itsNSigma); - recoQAHist.fill(HIST("hDaugAlphaITSCls"), daughterTrack.itsNCls()); - } + auto daughterTrack = tracks.rawIteratorAt(mcPartIndices[dauIDList[0]]); + bool isMoth = motherTrackCheck(motherTrack, hRecoMothCounter); + bool isDaug = daughterTrackCheck(daughterTrack, hRecoDaugCounter, getTPCNSigma(daughterTrack, pidTypeDaug)); + + recoQAHist.fill(HIST("hRecoCandidateCount"), 0.5); + recoQAHist.fill(HIST("hRecoMothCounter"), 0.5); + recoQAHist.fill(HIST("hMothITSCls"), motherTrack.itsNCls()); + recoQAHist.fill(HIST("hRecoDaugCounter"), 0.5); + recoQAHist.fill(HIST("hMothIsPVContributer"), motherTrack.isPVContributor() ? 1.5 : 0.5); + recoQAHist.fill(HIST("hDaugIsPVContributer"), daughterTrack.isPVContributor() ? 1.5 : 0.5); + + float itsNSigma = getITSNSigma(daughterTrack, itsResponse, pidTypeDaug); + if (daughterTrack.hasITS()) { + recoQAHist.fill(HIST("hDaugITSNSigma"), daughterTrack.sign() * daughterTrack.p(), itsNSigma); + recoQAHist.fill(HIST("hDaugITSCls"), daughterTrack.itsNCls()); + } - if (motherTrack.has_collision() && daughterTrack.has_collision()) { - recoQAHist.fill(HIST("hReco2BCandidateCount"), 1.5); - if (motherTrack.collisionId() == daughterTrack.collisionId()) { - recoQAHist.fill(HIST("hReco2BCandidateCount"), 2.5); - } + if (motherTrack.has_collision() && daughterTrack.has_collision()) { + recoQAHist.fill(HIST("hRecoCandidateCount"), 1.5); + if (motherTrack.collisionId() == daughterTrack.collisionId()) { + recoQAHist.fill(HIST("hRecoCandidateCount"), 2.5); } + } - if (isMoth && isDaug) { - recoQAHist.fill(HIST("hReco2BCandidateCount"), 3.5); - recoQAHist.fill(HIST("hReco2BDauAlphaPVsITSNSigma"), daughterTrack.sign() * daughterTrack.p(), itsNSigma); - } + if (isMoth && isDaug) { + recoQAHist.fill(HIST("hRecoCandidateCount"), 3.5); + recoQAHist.fill(HIST("hRecoDaugPVsITSNSigma"), daughterTrack.sign() * daughterTrack.p(), itsNSigma); } } } - - // qa for branching ratios and invariant mass - if (dChannel == k2body) { - genQAHist.fill(HIST("hGenHyperHelium4SigmaCounter"), isMatter ? 3.5 : 4.5); - float hyperHelium4SigmaMCMass = RecoDecay::m(std::array{std::array{dauMom[kDaugAlpha][0], dauMom[kDaugAlpha][1], dauMom[kDaugAlpha][2]}, std::array{dauMom[kDaugPion0][0], dauMom[kDaugPion0][1], dauMom[kDaugPion0][2]}}, std::array{o2::constants::physics::MassAlpha, o2::constants::physics::MassPi0}); - genQAHist.fill(HIST("hMcRecoInvMass"), hyperHelium4SigmaMCMass); - } else if (dChannel == k3body_p) { - genQAHist.fill(HIST("hGenHyperHelium4SigmaCounter"), isMatter ? 5.5 : 6.5); - float hyperHelium4SigmaMCMass = RecoDecay::m(std::array{std::array{dauMom[kDaugTriton][0], dauMom[kDaugTriton][1], dauMom[kDaugTriton][2]}, std::array{dauMom[kDaugProton][0], dauMom[kDaugProton][1], dauMom[kDaugProton][2]}, std::array{dauMom[kDaugPion0][0], dauMom[kDaugPion0][1], dauMom[kDaugPion0][2]}}, std::array{o2::constants::physics::MassTriton, o2::constants::physics::MassProton, o2::constants::physics::MassPi0}); - genQAHist.fill(HIST("hMcRecoInvMass"), hyperHelium4SigmaMCMass); - } else if (dChannel == k3body_n) { - genQAHist.fill(HIST("hGenHyperHelium4SigmaCounter"), isMatter ? 7.5 : 8.5); - float hyperHelium4SigmaMCMass = RecoDecay::m(std::array{std::array{dauMom[kDaugTriton][0], dauMom[kDaugTriton][1], dauMom[kDaugTriton][2]}, std::array{dauMom[kDaugNeutron][0], dauMom[kDaugNeutron][1], dauMom[kDaugNeutron][2]}, std::array{dauMom[kDaugChargedPion][0], dauMom[kDaugChargedPion][1], dauMom[kDaugChargedPion][2]}}, std::array{o2::constants::physics::MassTriton, o2::constants::physics::MassNeutron, o2::constants::physics::MassPiPlus}); - genQAHist.fill(HIST("hMcRecoInvMass"), hyperHelium4SigmaMCMass); - } } } } - PROCESS_SWITCH(Hyperhelium4sigmaQa, processMC, "do QA for MC prodcutions", false); + PROCESS_SWITCH(HyperkinkQa, processMC, "do QA for MC prodcutions", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask(cfgc), - adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), }; }