diff --git a/PWGLF/Tasks/Nuspex/hyperhelium4sigmaAnalysis.cxx b/PWGLF/Tasks/Nuspex/hyperhelium4sigmaAnalysis.cxx index bdbb7c5911b..22ee0a41d0f 100644 --- a/PWGLF/Tasks/Nuspex/hyperhelium4sigmaAnalysis.cxx +++ b/PWGLF/Tasks/Nuspex/hyperhelium4sigmaAnalysis.cxx @@ -10,10 +10,11 @@ // or submit itself to any jurisdiction. // /// \file hyperhelium4sigmaAnalysis.cxx -/// \brief Simple check for injected hyper-helium4sigma (H4S) in MC productions +/// \brief QA and analysis task for hyper-helium4sigma (He4S) /// \author Yuanzhe Wang #include +#include #include "Framework/runDataProcessing.h" #include "Framework/AnalysisTask.h" @@ -30,10 +31,22 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; -using CollisionsFull = soa::Join; +using CollisionsFull = soa::Join; using FullTracksExtIU = soa::Join; using MCLabeledTracksIU = soa::Join; +namespace +{ +constexpr int kITSLayers = 7; +constexpr int kITSInnerBarrelLayers = 3; +// constexpr int kITSOuterBarrelLayers = 4; +std::shared_ptr hMotherCounter; +std::shared_ptr hDauAlphaCounter; +std::shared_ptr hDauTritonCounter; +std::shared_ptr hDauProtonCounter; +std::shared_ptr hDauPionCounter; +} // namespace + //-------------------------------Check the decay channel of H4S------------------------------- enum Channel { k2body = 0, // helium4, pion0 @@ -43,11 +56,16 @@ enum Channel { }; template -Channel getDecayChannelH4S(TMCParticle const& particle) +Channel getDecayChannelH4S(TMCParticle const& particle, std::vector& list) { if (std::abs(particle.pdgCode()) != o2::constants::physics::Pdg::kHyperHelium4Sigma) { return kNDecayChannel; } + + // list: charged, 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; @@ -55,36 +73,47 @@ Channel getDecayChannelH4S(TMCParticle const& particle) 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(); } if (mcDaughter.pdgCode() == PDG_t::kNeutron) { haveNeuteron = true; + list[2] = mcDaughter.globalIndex(); } if (mcDaughter.pdgCode() == -PDG_t::kNeutron) { haveAntiNeuteron = true; + list[2] = mcDaughter.globalIndex(); } if (mcDaughter.pdgCode() == PDG_t::kPiPlus) { havePionPlus = true; + list[1] = mcDaughter.globalIndex(); } if (mcDaughter.pdgCode() == -PDG_t::kPiPlus) { havePionMinus = true; + list[1] = mcDaughter.globalIndex(); } if (mcDaughter.pdgCode() == PDG_t::kPi0) { havePion0 = true; + list[2] = mcDaughter.globalIndex(); } } @@ -98,12 +127,31 @@ Channel getDecayChannelH4S(TMCParticle const& particle) return kNDecayChannel; } + //-------------------------------------------------------------- +// check if the mcparticle is daughter of hyperhelium4sigma +template +bool isDaughterTrack(TMCParticle const& mcparticle, int pdgcode = o2::constants::physics::Pdg::kHyperHelium4Sigma) +{ + if (!mcparticle.has_mothers()) { + return false; + } + for (const auto& particleMother : mcparticle.template mothers_as()) { + if (std::abs(particleMother.pdgCode()) == pdgcode) { + return true; + } + } + return false; +} + +//-------------------------------------------------------------- +// analysis task for hyperhelium4sigma 2-body decay struct Hyperhelium4sigmaAnalysis { // Histograms are defined with HistogramRegistry HistogramRegistry registry{"registry", {}}; // Configurable for event selection + Configurable doEventCut{"doEventCut", true, "Apply event selection"}; Configurable cutzvertex{"cutzvertex", 10.0f, "Accepted z-vertex range (cm)"}; Configurable cutNSigmaAl{"cutNSigmaAl", 5, "NSigmaTPCAlpha"}; @@ -115,42 +163,131 @@ struct Hyperhelium4sigmaAnalysis { const AxisSpec nSigmaAxis{120, -6.f, 6.f, "n#sigma_{#alpha}"}; const AxisSpec massAxis{100, 3.85, 4.25, "m (GeV/#it{c}^{2})"}; - registry.add("hVertexZRec", "hVertexZRec", {HistType::kTH1F, {vertexZAxis}}); + registry.add("hEventCounter", "hEventCounter", HistType::kTH1F, {{2, 0, 2}}); + registry.add("hVertexZCollision", "hVertexZCollision", HistType::kTH1F, {vertexZAxis}); - registry.add("h2MassHyperhelium4sigmaPt", "h2MassHyperhelium4sigmaPt", {HistType::kTH2F, {ptAxis, massAxis}}); - registry.add("h2NSigmaAlPt", "h2NSigmaAlPt", {HistType::kTH2F, {ptAxis, nSigmaAxis}}); + if (doprocessData == true) { + registry.add("hCandidateCounter", "hCandidateCounter", HistType::kTH1F, {{3, 0, 3}}); + } + + if (doprocessMC == true) { + registry.add("hCandidateCounter", "hCandidateCounter", HistType::kTH1F, {{6, 0, 6}}); + registry.add("hDiffSVx", ";;#Delta x (cm)", HistType::kTH1F, {{200, -10, 10}}); + registry.add("hDiffSVy", ";;#Delta y (cm)", HistType::kTH1F, {{200, -10, 10}}); + registry.add("hDiffSVz", ";;#Delta z (cm)", HistType::kTH1F, {{200, -10, 10}}); + 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("h2MassHyperhelium4sigmaPt", "h2MassHyperhelium4sigmaPt", HistType::kTH2F, {{ptAxis, massAxis}}); + registry.add("h2NSigmaAlPt", "h2NSigmaAlPt", HistType::kTH2F, {{ptAxis, nSigmaAxis}}); } - void process(soa::Join::iterator const& collision, - aod::KinkCands const& KinkCands, FullTracksExtIU const&) + void processData(CollisionsFull const& collisions, aod::KinkCands const& KinkCands, FullTracksExtIU const&) { - if (std::abs(collision.posZ()) > cutzvertex || !collision.sel8()) { - return; + for (const auto& collision : collisions) { + registry.fill(HIST("hEventCounter"), 0); + if (doEventCut && (std::abs(collision.posZ()) > cutzvertex || !collision.sel8())) { + continue; + } + registry.fill(HIST("hEventCounter"), 1); + registry.fill(HIST("hVertexZCollision"), collision.posZ()); } - registry.fill(HIST("hVertexZRec"), collision.posZ()); + for (const auto& kinkCand : KinkCands) { + registry.fill(HIST("hCandidateCounter"), 0); + auto collision = kinkCand.collision_as(); + if (doEventCut && (std::abs(collision.posZ()) > cutzvertex || !collision.sel8())) { + continue; + } + registry.fill(HIST("hCandidateCounter"), 1); + auto dauTrack = kinkCand.trackDaug_as(); if (std::abs(dauTrack.tpcNSigmaAl()) > cutNSigmaAl) { 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}); + registry.fill(HIST("hCandidateCounter"), 2); registry.fill(HIST("h2MassHyperhelium4sigmaPt"), kinkCand.mothSign() * kinkCand.ptMoth(), invMass); registry.fill(HIST("h2NSigmaAlPt"), kinkCand.mothSign() * kinkCand.ptDaug(), dauTrack.tpcNSigmaAl()); } } + PROCESS_SWITCH(Hyperhelium4sigmaAnalysis, processData, "process data", true); + + void processMC(CollisionsFull const& collisions, aod::KinkCands const& KinkCands, MCLabeledTracksIU const&, aod::McParticles const&) + { + std::vector dauIDList; + + for (const auto& collision : collisions) { + registry.fill(HIST("hEventCounter"), 0); + if (doEventCut && (!collision.selection_bit(aod::evsel::kIsTriggerTVX) || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder) || std::abs(collision.posZ()) > cutzvertex)) { + continue; + } + registry.fill(HIST("hEventCounter"), 1); + registry.fill(HIST("hVertexZCollision"), collision.posZ()); + } + + for (const auto& kinkCand : KinkCands) { + registry.fill(HIST("hCandidateCounter"), 0); + auto collision = kinkCand.collision_as(); + if (doEventCut && (!collision.selection_bit(aod::evsel::kIsTriggerTVX) || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder) || std::abs(collision.posZ()) > cutzvertex)) { + continue; + } + registry.fill(HIST("hCandidateCounter"), 1); + + auto motherTrack = kinkCand.trackMoth_as(); + auto dauTrack = kinkCand.trackDaug_as(); + if (!motherTrack.has_mcParticle() || !dauTrack.has_mcParticle()) { + continue; + } + registry.fill(HIST("hCandidateCounter"), 2); + + auto mcMotherTrack = motherTrack.mcParticle_as(); + auto mcDauTrack = dauTrack.mcParticle_as(); + auto dChannel = getDecayChannelH4S(mcMotherTrack, dauIDList); + if (dChannel != k2body) { + continue; + } + registry.fill(HIST("hCandidateCounter"), 3); + + if (dauIDList[0] != mcDauTrack.globalIndex()) { + continue; + } + registry.fill(HIST("hCandidateCounter"), 4); + + if (std::abs(dauTrack.tpcNSigmaAl()) > cutNSigmaAl) { + continue; + } + registry.fill(HIST("hCandidateCounter"), 5); + + 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}); + registry.fill(HIST("hDiffSVx"), kinkCand.xDecVtx() - mcDauTrack.vx()); + registry.fill(HIST("hDiffSVy"), kinkCand.yDecVtx() - mcDauTrack.vy()); + registry.fill(HIST("hDiffSVz"), kinkCand.zDecVtx() - mcDauTrack.vz()); + 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("h2MassHyperhelium4sigmaPt"), kinkCand.mothSign() * kinkCand.ptMoth(), invMass); + registry.fill(HIST("h2NSigmaAlPt"), kinkCand.mothSign() * kinkCand.ptDaug(), dauTrack.tpcNSigmaAl()); + } + } + PROCESS_SWITCH(Hyperhelium4sigmaAnalysis, processMC, "process MC", false); }; //-------------------------------------------------------------- // check the performance of mcparticle struct Hyperhelium4sigmaQa { - // Basic checks - HistogramRegistry registry{"registry", {}}; + + HistogramRegistry genQAHist{"genQAHist", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + HistogramRegistry recoQAHist{"recoQAHist", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; ConfigurableAxis ptBins{"ptBins", {200, 0.f, 10.f}, "Binning for #it{p}_{T} (GeV/#it{c})"}; ConfigurableAxis ctBins{"ctBins", {100, 0.f, 25.f}, "Binning for c#it{t} (cm)"}; ConfigurableAxis rigidityBins{"rigidityBins", {200, -10.f, 10.f}, "Binning for #it{p}_{T} (GeV/#it{c})"}; ConfigurableAxis nsigmaBins{"nsigmaBins", {120, -6.f, 6.f}, "Binning for n sigma"}; - ConfigurableAxis invMassBins{"invMassBins", {100, 3.85, 4.15f}, "Binning for invariant mass (GeV/#it{c}^{2})"}; + ConfigurableAxis invMassBins{"invMassBins", {100, 3.85f, 4.15f}, "Binning for invariant mass (GeV/#it{c}^{2})"}; + ConfigurableAxis radiusBins{"radiusBins", {40, 0.f, 40.f}, "Binning for radius in xy plane (cm)"}; void init(InitContext&) { @@ -160,49 +297,82 @@ struct Hyperhelium4sigmaQa { const AxisSpec rigidityAxis{rigidityBins, "p/z (GeV/#it{c})"}; const AxisSpec nsigmaAxis{nsigmaBins, "TPC n#sigma"}; const AxisSpec invMassAxis{invMassBins, "Inv Mass (GeV/#it{c}^{2})"}; - - auto hCollCounter = registry.add("hCollCounter", "hCollCounter", HistType::kTH1F, {{2, 0.0f, 2.0f}}); - registry.get(HIST("hCollCounter"))->GetXaxis()->SetBinLabel(1, "Reconstructed Collisions"); - registry.get(HIST("hCollCounter"))->GetXaxis()->SetBinLabel(2, "Selected"); - auto hMcCollCounter = registry.add("hMcCollCounter", "hMcCollCounter", HistType::kTH1F, {{2, 0.0f, 2.0f}}); - registry.get(HIST("hMcCollCounter"))->GetXaxis()->SetBinLabel(1, "MC Collisions"); - registry.get(HIST("hMcCollCounter"))->GetXaxis()->SetBinLabel(2, "Reconstructed"); - - auto hGenHyperHelium4SigmaCounter = registry.add("hGenHyperHelium4SigmaCounter", "", HistType::kTH1F, {{10, 0.f, 10.f}}); - registry.get(HIST("hGenHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(1, "H4S All"); - registry.get(HIST("hGenHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(2, "Matter"); - registry.get(HIST("hGenHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(3, "AntiMatter"); - registry.get(HIST("hGenHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(4, "#alpha + #pi^{0}"); - registry.get(HIST("hGenHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(5, "#bar{#alpha} + #pi^{0}"); - registry.get(HIST("hGenHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(6, "t + p + #pi^{0}"); - registry.get(HIST("hGenHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(7, "#bar{t} + #bar{p} + #pi^{0}"); - registry.get(HIST("hGenHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(8, "t + n + #pi^{+}"); - registry.get(HIST("hGenHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(9, "#bar{t} + #bar{n} + #pi^{+}"); - registry.get(HIST("hGenHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(10, "Unexpected"); - - auto hEvtSelectedHyperHelium4SigmaCounter = registry.add("hEvtSelectedHyperHelium4SigmaCounter", "", HistType::kTH1F, {{2, 0.f, 2.f}}); - registry.get(HIST("hEvtSelectedHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(1, "Generated"); - registry.get(HIST("hEvtSelectedHyperHelium4SigmaCounter"))->GetXaxis()->SetBinLabel(2, "Survived"); - - registry.add("hGenHyperHelium4SigmaP", "", HistType::kTH1F, {ptAxis}); - registry.add("hGenHyperHelium4SigmaPt", "", HistType::kTH1F, {ptAxis}); - registry.add("hGenHyperHelium4SigmaCt", "", HistType::kTH1F, {ctAxis}); - registry.add("hMcRecoInvMass", "", HistType::kTH1F, {invMassAxis}); - - registry.add("hDauHelium4TPCNSigma", "", HistType::kTH2F, {rigidityAxis, nsigmaAxis}); - registry.add("hDauTritonTPCNSigma", "", HistType::kTH2F, {rigidityAxis, nsigmaAxis}); - registry.add("hDauProtonTPCNSigma", "", HistType::kTH2F, {rigidityAxis, nsigmaAxis}); + const AxisSpec diffPtAxis{200, -10.f, 10.f, "#Delta p_{T} (GeV/#it{c})"}; + const AxisSpec radiusAxis{radiusBins, "R (cm)"}; + + auto hCollCounter = genQAHist.add("hCollCounter", "hCollCounter", HistType::kTH1F, {{2, 0.0f, 2.0f}}); + hCollCounter->GetXaxis()->SetBinLabel(1, "Reconstructed Collisions"); + hCollCounter->GetXaxis()->SetBinLabel(2, "Selected"); + auto hMcCollCounter = genQAHist.add("hMcCollCounter", "hMcCollCounter", HistType::kTH1F, {{2, 0.0f, 2.0f}}); + hMcCollCounter->GetXaxis()->SetBinLabel(1, "MC Collisions"); + hMcCollCounter->GetXaxis()->SetBinLabel(2, "Reconstructed"); + + auto hGenHyperHelium4SigmaCounter = genQAHist.add("hGenHyperHelium4SigmaCounter", "", HistType::kTH1F, {{10, 0.f, 10.f}}); + hGenHyperHelium4SigmaCounter->GetXaxis()->SetBinLabel(1, "H4S 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, "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, {ptAxis}); + genQAHist.add("hGenHyperHelium4SigmaPt", "", HistType::kTH1F, {ptAxis}); + genQAHist.add("hGenHyperHelium4SigmaCt", "", HistType::kTH1F, {ctAxis}); + genQAHist.add("hMcRecoInvMass", "", HistType::kTH1F, {invMassAxis}); + + // efficiency/criteria studies for tracks which are true candidates + hMotherCounter = recoQAHist.add("hMotherCounter", "", HistType::kTH1F, {{9, 0.f, 9.f}}); + hMotherCounter->GetXaxis()->SetBinLabel(1, "Generated"); + hMotherCounter->GetXaxis()->SetBinLabel(2, "Reconstructed"); + hMotherCounter->GetXaxis()->SetBinLabel(3, "eta"); + hMotherCounter->GetXaxis()->SetBinLabel(4, "has collision"); + hMotherCounter->GetXaxis()->SetBinLabel(5, "ITSonly"); + hMotherCounter->GetXaxis()->SetBinLabel(6, "ITS hits"); + hMotherCounter->GetXaxis()->SetBinLabel(7, "ITS IR"); + hMotherCounter->GetXaxis()->SetBinLabel(8, "ITS chi2"); + hMotherCounter->GetXaxis()->SetBinLabel(9, "pt"); + recoQAHist.add("hTrueMotherRVsDiffPt", ";#Delta p_{T} (GeV/#it{c});R (cm);", HistType::kTH2F, {diffPtAxis, radiusAxis}); + recoQAHist.add("hTrueMotherRVsDiffPz", ";#Delta p_{z} (GeV/#it{c});R (cm);", HistType::kTH2F, {diffPtAxis, radiusAxis}); + recoQAHist.add("hGoodMotherRVsDiffPt", ";#Delta p_{T} (GeV/#it{c});R (cm);", HistType::kTH2F, {diffPtAxis, radiusAxis}); + recoQAHist.add("hGoodMotherRVsDiffPz", ";#Delta p_{z} (GeV/#it{c});R (cm);", HistType::kTH2F, {diffPtAxis, radiusAxis}); + + hDauAlphaCounter = recoQAHist.add("hDauAlphaCounter", "", HistType::kTH1F, {{7, 0.f, 7.f}}); + hDauTritonCounter = recoQAHist.add("hDauTritonCounter", "", HistType::kTH1F, {{7, 0.f, 7.f}}); + hDauProtonCounter = recoQAHist.add("hDauProtonCounter", "", HistType::kTH1F, {{7, 0.f, 7.f}}); + hDauPionCounter = recoQAHist.add("hDauPionCounter", "", HistType::kTH1F, {{7, 0.f, 7.f}}); + + recoQAHist.add("hDauAlphaTPCNSigma", "", HistType::kTH2F, {rigidityAxis, nsigmaAxis}); + recoQAHist.add("hDauTritonTPCNSigma", "", HistType::kTH2F, {rigidityAxis, nsigmaAxis}); + recoQAHist.add("hDauProtonTPCNSigma", "", HistType::kTH2F, {rigidityAxis, nsigmaAxis}); + recoQAHist.add("hDauPionTPCNSigma", "", HistType::kTH2F, {rigidityAxis, nsigmaAxis}); } } - // Configurable eventSel8Cut{"eventSel8Cut", true, "flag to enable event sel8 selection"}; + Configurable skipRejectedEvents{"skipRejectedEvents", false, "Flag to skip events that fail event selection cuts"}; Configurable mcEventCut{"mcEventCut", true, "flag to enable mc event selection: kIsTriggerTVX and kNoTimeFrameBorder"}; Configurable eventPosZCut{"eventPosZCut", true, "flag to enable event posZ selection"}; Configurable maxPosZ{"maxPosZ", 10.f, "max pv posZ for event selection"}; + Configurable etaMax{"etaMax", 1., "eta cut for tracks"}; + Configurable minPtMoth{"minPtMoth", 0.5, "Minimum pT of the hypercandidate"}; + Configurable tpcPidNsigmaCut{"tpcPidNsigmaCut", 5, "tpcPidNsigmaCut"}; + Configurable nTPCClusMinDaug{"nTPCClusMinDaug", 80, "daug NTPC clusters cut"}; + Configurable itsMaxChi2{"itsMaxChi2", 36, "max chi2 for ITS"}; + Configurable minRatioTPCNCls{"minRatioTPCNCls", 0.8, "min ratio of TPC crossed rows to findable clusters"}; + Preslice permcCollision = o2::aod::mcparticle::mcCollisionId; + // construct index array from mcParticle to track std::vector mcPartIndices; + template void setTrackIDForMC(aod::McParticles const& particlesMC, TTrackTable const& tracks) { @@ -225,6 +395,82 @@ struct Hyperhelium4sigmaQa { } } + // qa for mother track selection + template + bool motherTrackCheck(const TTrack& track, const std::shared_ptr hist) + { + hist->Fill(1); + + if (std::abs(track.eta()) > etaMax) { + return false; + } + hist->Fill(2); + + if (!track.has_collision()) { + return false; + } + hist->Fill(3); + + if (!track.hasITS() || track.hasTPC() || track.hasTOF()) { + return false; + } + hist->Fill(4); + + if (track.itsNCls() >= kITSLayers - 1) { + return false; + } + hist->Fill(5); + + if (track.itsNClsInnerBarrel() != kITSInnerBarrelLayers) { + return false; + } + hist->Fill(6); + + if (track.itsChi2NCl() >= itsMaxChi2) { + return false; + } + hist->Fill(7); + + if (track.pt() <= minPtMoth) { + return false; + } + hist->Fill(8); + + return true; + } + + // qa for daughter track selection + template + void daughterTrackCheck(const TTrack& track, const std::shared_ptr hist, float tpcNSigma) + { + hist->Fill(1); + + if (std::abs(track.eta()) > etaMax) { + return; + } + hist->Fill(2); + + if (!track.hasITS() || !track.hasTPC()) { + return; + } + hist->Fill(3); + + if (track.itsNClsInnerBarrel() != 0 && track.itsNCls() > kITSInnerBarrelLayers && track.tpcNClsCrossedRows() <= minRatioTPCNCls * track.tpcNClsFindable() && track.tpcNClsFound() <= nTPCClusMinDaug) { + return; + } + hist->Fill(4); + + if (std::abs(tpcNSigma) > tpcPidNsigmaCut) { + return; + } + hist->Fill(5); + + if (track.hasTOF()) { + return; + } + hist->Fill(6); + } + void processData(o2::aod::Collisions const&) { // dummy process function; @@ -233,29 +479,34 @@ struct Hyperhelium4sigmaQa { void processMC(aod::McCollisions const& mcCollisions, aod::McParticles const& particlesMC, o2::soa::Join const& collisions, MCLabeledTracksIU const& tracks) { + + // check mcparticles setTrackIDForMC(particlesMC, tracks); std::vector selectedEvents(collisions.size()); + std::vector dauIDList; int nevts = 0; for (const auto& collision : collisions) { - registry.fill(HIST("hCollCounter"), 0.5); + genQAHist.fill(HIST("hCollCounter"), 0.5); if (mcEventCut && (!collision.selection_bit(aod::evsel::kIsTriggerTVX) || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder))) { continue; } if (eventPosZCut && std::abs(collision.posZ()) > maxPosZ) { // 10cm continue; } - registry.fill(HIST("hCollCounter"), 1.5); + genQAHist.fill(HIST("hCollCounter"), 1.5); selectedEvents[nevts++] = collision.mcCollision_as().globalIndex(); } selectedEvents.resize(nevts); for (const auto& mcCollision : mcCollisions) { - registry.fill(HIST("hMcCollCounter"), 0.5); + genQAHist.fill(HIST("hMcCollCounter"), 0.5); const auto evtReconstructedAndSelected = std::find(selectedEvents.begin(), selectedEvents.end(), mcCollision.globalIndex()) != selectedEvents.end(); if (evtReconstructedAndSelected) { // Check that the event is reconstructed and that the reconstructed events pass the selection - registry.fill(HIST("hMcCollCounter"), 1.5); + genQAHist.fill(HIST("hMcCollCounter"), 1.5); } else { - // continue; + if (skipRejectedEvents) { + continue; + } } const auto& dparticlesMC = particlesMC.sliceBy(permcCollision, mcCollision.globalIndex()); @@ -264,47 +515,55 @@ struct Hyperhelium4sigmaQa { bool isMatter; if (mcparticle.pdgCode() == o2::constants::physics::Pdg::kHyperHelium4Sigma) { - registry.fill(HIST("hGenHyperHelium4SigmaCounter"), 1.5); + genQAHist.fill(HIST("hGenHyperHelium4SigmaCounter"), 1.5); isMatter = true; } else if (mcparticle.pdgCode() == -o2::constants::physics::Pdg::kHyperHelium4Sigma) { - registry.fill(HIST("hGenHyperHelium4SigmaCounter"), 2.5); + genQAHist.fill(HIST("hGenHyperHelium4SigmaCounter"), 2.5); isMatter = false; } else { continue; } - registry.fill(HIST("hGenHyperHelium4SigmaCounter"), 0.5); - registry.fill(HIST("hEvtSelectedHyperHelium4SigmaCounter"), 0.5); + genQAHist.fill(HIST("hGenHyperHelium4SigmaCounter"), 0.5); + genQAHist.fill(HIST("hEvtSelectedHyperHelium4SigmaCounter"), 0.5); if (evtReconstructedAndSelected) { - registry.fill(HIST("hEvtSelectedHyperHelium4SigmaCounter"), 1.5); + genQAHist.fill(HIST("hEvtSelectedHyperHelium4SigmaCounter"), 1.5); + } + + auto dChannel = getDecayChannelH4S(mcparticle, dauIDList); + if (dChannel == kNDecayChannel) { + genQAHist.fill(HIST("hGenHyperHelium4SigmaCounter"), 9.5); + continue; } + // bool isAllDauReconstructed = mcPartIndices[dauIDList[0]] != -1 && (dChannel == k2body ? true : mcPartIndices[dauIDList[1]] != -1); + + // check for mother tracks + recoQAHist.fill(HIST("hMotherCounter"), 0); + double svPos[3] = {-999, -999, -999}; - double dauHelium4Mom[3] = {-999, -999, -999}; + double dauAlphaMom[3] = {-999, -999, -999}; double dauTritonMom[3] = {-999, -999, -999}; double dauProtonMom[3] = {-999, -999, -999}; double dauNeuteronMom[3] = {-999, -999, -999}; double dauChargedPionMom[3] = {-999, -999, -999}; double dauPion0Mom[3] = {-999, -999, -999}; - auto dChannel = getDecayChannelH4S(mcparticle); - if (dChannel == kNDecayChannel) { - registry.fill(HIST("hGenHyperHelium4SigmaCounter"), 9.5); - continue; - } for (const auto& mcparticleDaughter : mcparticle.daughters_as()) { if (std::abs(mcparticleDaughter.pdgCode()) == o2::constants::physics::Pdg::kAlpha) { - dauHelium4Mom[0] = mcparticleDaughter.px(); - dauHelium4Mom[1] = mcparticleDaughter.py(); - dauHelium4Mom[2] = mcparticleDaughter.pz(); + dauAlphaMom[0] = mcparticleDaughter.px(); + dauAlphaMom[1] = mcparticleDaughter.py(); + dauAlphaMom[2] = mcparticleDaughter.pz(); // get SV position for 2body decay svPos[0] = mcparticleDaughter.vx(); svPos[1] = mcparticleDaughter.vy(); svPos[2] = mcparticleDaughter.vz(); + recoQAHist.fill(HIST("hDauAlphaCounter"), 0); if (mcPartIndices[mcparticleDaughter.globalIndex()] != -1) { auto track = tracks.rawIteratorAt(mcPartIndices[mcparticleDaughter.globalIndex()]); - registry.fill(HIST("hDauHelium4TPCNSigma"), track.p() * track.sign(), track.tpcNSigmaAl()); + recoQAHist.fill(HIST("hDauAlphaTPCNSigma"), track.p() * track.sign(), track.tpcNSigmaAl()); + daughterTrackCheck(track, hDauAlphaCounter, track.tpcNSigmaAl()); } } else if (std::abs(mcparticleDaughter.pdgCode()) == o2::constants::physics::Pdg::kTriton) { dauTritonMom[0] = mcparticleDaughter.px(); @@ -316,18 +575,22 @@ struct Hyperhelium4sigmaQa { svPos[1] = mcparticleDaughter.vy(); svPos[2] = mcparticleDaughter.vz(); + recoQAHist.fill(HIST("hDauTritonCounter"), 0); if (mcPartIndices[mcparticleDaughter.globalIndex()] != -1) { auto track = tracks.rawIteratorAt(mcPartIndices[mcparticleDaughter.globalIndex()]); - registry.fill(HIST("hDauTritonTPCNSigma"), track.p() * track.sign(), track.tpcNSigmaTr()); + recoQAHist.fill(HIST("hDauTritonTPCNSigma"), track.p() * track.sign(), track.tpcNSigmaTr()); + daughterTrackCheck(track, hDauTritonCounter, track.tpcNSigmaTr()); } } else if (std::abs(mcparticleDaughter.pdgCode()) == PDG_t::kProton) { dauProtonMom[0] = mcparticleDaughter.px(); dauProtonMom[1] = mcparticleDaughter.py(); dauProtonMom[2] = mcparticleDaughter.pz(); + recoQAHist.fill(HIST("hDauProtonCounter"), 0); if (mcPartIndices[mcparticleDaughter.globalIndex()] != -1) { auto track = tracks.rawIteratorAt(mcPartIndices[mcparticleDaughter.globalIndex()]); - registry.fill(HIST("hDauProtonTPCNSigma"), track.p() * track.sign(), track.tpcNSigmaPr()); + recoQAHist.fill(HIST("hDauProtonTPCNSigma"), track.p() * track.sign(), track.tpcNSigmaPr()); + daughterTrackCheck(track, hDauProtonCounter, track.tpcNSigmaPr()); } } else if (std::abs(mcparticleDaughter.pdgCode()) == PDG_t::kNeutron) { dauNeuteronMom[0] = mcparticleDaughter.px(); @@ -337,6 +600,13 @@ struct Hyperhelium4sigmaQa { dauChargedPionMom[0] = mcparticleDaughter.px(); dauChargedPionMom[1] = mcparticleDaughter.py(); dauChargedPionMom[2] = mcparticleDaughter.pz(); + + recoQAHist.fill(HIST("hDauPionCounter"), 0); + if (mcPartIndices[mcparticleDaughter.globalIndex()] != -1) { + auto track = tracks.rawIteratorAt(mcPartIndices[mcparticleDaughter.globalIndex()]); + recoQAHist.fill(HIST("hDauPionTPCNSigma"), track.p() * track.sign(), track.tpcNSigmaPi()); + daughterTrackCheck(track, hDauPionCounter, track.tpcNSigmaPi()); + } } else if (mcparticleDaughter.pdgCode() == PDG_t::kPi0) { dauPion0Mom[0] = mcparticleDaughter.px(); dauPion0Mom[1] = mcparticleDaughter.py(); @@ -344,35 +614,47 @@ struct Hyperhelium4sigmaQa { } } - registry.fill(HIST("hGenHyperHelium4SigmaP"), mcparticle.p()); - registry.fill(HIST("hGenHyperHelium4SigmaPt"), mcparticle.pt()); + genQAHist.fill(HIST("hGenHyperHelium4SigmaP"), mcparticle.p()); + genQAHist.fill(HIST("hGenHyperHelium4SigmaPt"), mcparticle.pt()); double ct = RecoDecay::sqrtSumOfSquares(svPos[0] - mcparticle.vx(), svPos[1] - mcparticle.vy(), svPos[2] - mcparticle.vz()) * o2::constants::physics::MassHyperHelium4Sigma / mcparticle.p(); - registry.fill(HIST("hGenHyperHelium4SigmaCt"), ct); + genQAHist.fill(HIST("hGenHyperHelium4SigmaCt"), ct); + + if (mcPartIndices[mcparticle.globalIndex()] != -1) { + auto motherTrack = tracks.rawIteratorAt(mcPartIndices[mcparticle.globalIndex()]); + bool isGoodMother = motherTrackCheck(motherTrack, hMotherCounter); + double svR = RecoDecay::sqrtSumOfSquares(svPos[0], svPos[1]); + recoQAHist.fill(HIST("hTrueMotherRVsDiffPt"), mcparticle.pt() - 2 * motherTrack.pt(), svR); + recoQAHist.fill(HIST("hTrueMotherRVsDiffPz"), mcparticle.pz() - 2 * motherTrack.pz(), svR); + if (isGoodMother) { + recoQAHist.fill(HIST("hGoodMotherRVsDiffPt"), mcparticle.pt() - 2 * motherTrack.pt(), svR); + recoQAHist.fill(HIST("hGoodMotherRVsDiffPz"), mcparticle.pz() - 2 * motherTrack.pz(), svR); + } + } if (dChannel == k2body) { if (isMatter) { - registry.fill(HIST("hGenHyperHelium4SigmaCounter"), 3.5); + genQAHist.fill(HIST("hGenHyperHelium4SigmaCounter"), 3.5); } else { - registry.fill(HIST("hGenHyperHelium4SigmaCounter"), 4.5); + genQAHist.fill(HIST("hGenHyperHelium4SigmaCounter"), 4.5); } - double hyperHelium4SigmaMCMass = RecoDecay::m(std::array{std::array{dauHelium4Mom[0], dauHelium4Mom[1], dauHelium4Mom[2]}, std::array{dauPion0Mom[0], dauPion0Mom[1], dauPion0Mom[2]}}, std::array{o2::constants::physics::MassAlpha, o2::constants::physics::MassPi0}); - registry.fill(HIST("hMcRecoInvMass"), hyperHelium4SigmaMCMass); + double hyperHelium4SigmaMCMass = RecoDecay::m(std::array{std::array{dauAlphaMom[0], dauAlphaMom[1], dauAlphaMom[2]}, std::array{dauPion0Mom[0], dauPion0Mom[1], dauPion0Mom[2]}}, std::array{o2::constants::physics::MassAlpha, o2::constants::physics::MassPi0}); + genQAHist.fill(HIST("hMcRecoInvMass"), hyperHelium4SigmaMCMass); } else if (dChannel == k3body_p) { if (isMatter) { - registry.fill(HIST("hGenHyperHelium4SigmaCounter"), 5.5); + genQAHist.fill(HIST("hGenHyperHelium4SigmaCounter"), 5.5); } else { - registry.fill(HIST("hGenHyperHelium4SigmaCounter"), 6.5); + genQAHist.fill(HIST("hGenHyperHelium4SigmaCounter"), 6.5); } double hyperHelium4SigmaMCMass = RecoDecay::m(std::array{std::array{dauTritonMom[0], dauTritonMom[1], dauTritonMom[2]}, std::array{dauProtonMom[0], dauProtonMom[1], dauProtonMom[2]}, std::array{dauPion0Mom[0], dauPion0Mom[1], dauPion0Mom[2]}}, std::array{o2::constants::physics::MassTriton, o2::constants::physics::MassProton, o2::constants::physics::MassPi0}); - registry.fill(HIST("hMcRecoInvMass"), hyperHelium4SigmaMCMass); + genQAHist.fill(HIST("hMcRecoInvMass"), hyperHelium4SigmaMCMass); } else if (dChannel == k3body_n) { if (isMatter) { - registry.fill(HIST("hGenHyperHelium4SigmaCounter"), 7.5); + genQAHist.fill(HIST("hGenHyperHelium4SigmaCounter"), 7.5); } else { - registry.fill(HIST("hGenHyperHelium4SigmaCounter"), 8.5); + genQAHist.fill(HIST("hGenHyperHelium4SigmaCounter"), 8.5); } double hyperHelium4SigmaMCMass = RecoDecay::m(std::array{std::array{dauTritonMom[0], dauTritonMom[1], dauTritonMom[2]}, std::array{dauNeuteronMom[0], dauNeuteronMom[1], dauNeuteronMom[2]}, std::array{dauChargedPionMom[0], dauChargedPionMom[1], dauChargedPionMom[2]}}, std::array{o2::constants::physics::MassTriton, o2::constants::physics::MassNeutron, o2::constants::physics::MassPionCharged}); - registry.fill(HIST("hMcRecoInvMass"), hyperHelium4SigmaMCMass); + genQAHist.fill(HIST("hMcRecoInvMass"), hyperHelium4SigmaMCMass); } } }