From 87bc031562bb185078b1e19a65d8a9dcd14acc95 Mon Sep 17 00:00:00 2001 From: Mattia Faggin Date: Fri, 29 Jan 2021 14:46:34 +0100 Subject: [PATCH] PWGHF: first implementation of Xic->pKpi filtering --- .../include/AnalysisCore/HFConfigurables.h | 6 + .../HFCandidateSelectionTables.h | 11 + .../AnalysisDataModel/HFSecondaryVertex.h | 38 +- .../Tasks/PWGHF/HFCandidateCreator3Prong.cxx | 16 + .../Tasks/PWGHF/HFTrackIndexSkimsCreator.cxx | 23 +- .../Tasks/PWGHF/HFXicCandidateSelector.cxx | 447 ++++++++++++++++++ Analysis/Tasks/PWGHF/taskXic.cxx | 166 +++++++ 7 files changed, 705 insertions(+), 2 deletions(-) create mode 100644 Analysis/Tasks/PWGHF/HFXicCandidateSelector.cxx create mode 100644 Analysis/Tasks/PWGHF/taskXic.cxx diff --git a/Analysis/Core/include/AnalysisCore/HFConfigurables.h b/Analysis/Core/include/AnalysisCore/HFConfigurables.h index 34d0d1c8c14c3..50faa19f0f50b 100644 --- a/Analysis/Core/include/AnalysisCore/HFConfigurables.h +++ b/Analysis/Core/include/AnalysisCore/HFConfigurables.h @@ -53,6 +53,12 @@ class HFTrackIndexSkimsCreatorConfigs double mInvMassDsToPiKKMax = 2.15; //original value 2.2 double mCPADsToPiKKMin = 0.5; double mDecLenDsToPiKKMin = 0.; + // 3-prong cuts - XicToPKPi + double mPtXicToPKPiMin = 1.; // + double mInvMassXicToPKPiMin = 2.25; // + double mInvMassXicToPKPiMax = 2.70; // + double mCPAXicToPKPiMin = 0.5; + double mDecLenXicToPKPiMin = 0.; private: ClassDef(HFTrackIndexSkimsCreatorConfigs, 1); diff --git a/Analysis/DataModel/include/AnalysisDataModel/HFCandidateSelectionTables.h b/Analysis/DataModel/include/AnalysisDataModel/HFCandidateSelectionTables.h index 12eb12a04586c..b2ef5557bfa3f 100644 --- a/Analysis/DataModel/include/AnalysisDataModel/HFCandidateSelectionTables.h +++ b/Analysis/DataModel/include/AnalysisDataModel/HFCandidateSelectionTables.h @@ -39,4 +39,15 @@ DECLARE_SOA_COLUMN(IsSelJpsiToEE, isSelJpsiToEE, int); } // namespace hf_selcandidate_jpsi DECLARE_SOA_TABLE(HFSelJpsiToEECandidate, "AOD", "HFSELJPSICAND", hf_selcandidate_jpsi::IsSelJpsiToEE); } // namespace o2::aod + +namespace o2::aod +{ +namespace hf_selcandidate_xic +{ +DECLARE_SOA_COLUMN(IsSelXicpKpi, isSelXicpKpi, int); +DECLARE_SOA_COLUMN(IsSelXicpiKp, isSelXicpiKp, int); +} // namespace hf_selcandidate_xic +DECLARE_SOA_TABLE(HFSelXicpKpiCandidate, "AOD", "HFSELXICCAND", hf_selcandidate_xic::IsSelXicpKpi, hf_selcandidate_xic::IsSelXicpiKp); +} // namespace o2::aod + #endif // O2_ANALYSIS_HFCANDIDATESELECTIONTABLES_H_ diff --git a/Analysis/DataModel/include/AnalysisDataModel/HFSecondaryVertex.h b/Analysis/DataModel/include/AnalysisDataModel/HFSecondaryVertex.h index b4df0f9597b23..0d9e3185f19b3 100644 --- a/Analysis/DataModel/include/AnalysisDataModel/HFSecondaryVertex.h +++ b/Analysis/DataModel/include/AnalysisDataModel/HFSecondaryVertex.h @@ -55,6 +55,7 @@ DECLARE_SOA_COLUMN(JpsiToEEFlag, jpsiToEEFlag, uint8_t); DECLARE_SOA_COLUMN(DPlusPiKPiFlag, dPlusPiKPiFlag, uint8_t); DECLARE_SOA_COLUMN(LcPKPiFlag, lcPKPiFlag, uint8_t); DECLARE_SOA_COLUMN(DsKKPiFlag, dsKKPiFlag, uint8_t); +DECLARE_SOA_COLUMN(XicPKPiFlag, xicPKPiFlag, uint8_t); } // namespace hf_track_index DECLARE_SOA_TABLE(HfTrackIndexProng2, "AOD", "HFTRACKIDXP2", @@ -75,7 +76,8 @@ DECLARE_SOA_TABLE(HfTrackIndexProng3, "AOD", "HFTRACKIDXP3", DECLARE_SOA_TABLE(HfCutStatusProng3, "AOD", "HFCUTSTATUSP3", hf_track_index::DPlusPiKPiFlag, hf_track_index::LcPKPiFlag, - hf_track_index::DsKKPiFlag); + hf_track_index::DsKKPiFlag, + hf_track_index::XicPKPiFlag); // general decay properties namespace hf_cand @@ -312,6 +314,7 @@ DECLARE_SOA_DYNAMIC_COLUMN(MaxNormalisedDeltaIP, maxNormalisedDeltaIP, [](float // MC matching result: // - ±DPlusToPiKPi: D± → π± K∓ π± // - ±LcToPKPi: Λc± → p± K∓ π± +// - ±XicToPKPi: Ξc± → p± K∓ π± DECLARE_SOA_COLUMN(FlagMCMatchRec, flagMCMatchRec, int8_t); // reconstruction level DECLARE_SOA_COLUMN(FlagMCMatchGen, flagMCMatchGen, int8_t); // generator level @@ -319,6 +322,7 @@ DECLARE_SOA_COLUMN(FlagMCMatchGen, flagMCMatchGen, int8_t); // generator level enum DecayType { DPlusToPiKPi = 0, LcToPKPi, DsToPiKK, + XicToPKPi, N3ProngDecays }; //always keep N3ProngDecays at the end // functions for specific particles @@ -380,6 +384,38 @@ auto InvMassLcpiKp(const T& candidate) { return candidate.m(array{RecoDecay::getMassPDG(kPiPlus), RecoDecay::getMassPDG(kKPlus), RecoDecay::getMassPDG(kProton)}); } + +// Ξc± → p± K∓ π± +template +auto CtXic(const T& candidate) +{ + return candidate.ct(RecoDecay::getMassPDG(4232)); +} + +template +auto YXic(const T& candidate) +{ + return candidate.y(RecoDecay::getMassPDG(4232)); +} + +template +auto EXic(const T& candidate) +{ + return candidate.e(RecoDecay::getMassPDG(4232)); +} + +template +auto InvMassXicpKpi(const T& candidate) +{ + return candidate.m(array{RecoDecay::getMassPDG(kProton), RecoDecay::getMassPDG(kKPlus), RecoDecay::getMassPDG(kPiPlus)}); +} + +template +auto InvMassXicpiKp(const T& candidate) +{ + return candidate.m(array{RecoDecay::getMassPDG(kPiPlus), RecoDecay::getMassPDG(kKPlus), RecoDecay::getMassPDG(kProton)}); +} + } // namespace hf_cand_prong3 // 3-prong decay candidate table diff --git a/Analysis/Tasks/PWGHF/HFCandidateCreator3Prong.cxx b/Analysis/Tasks/PWGHF/HFCandidateCreator3Prong.cxx index b56d8451fa0d8..547a21032b0fc 100644 --- a/Analysis/Tasks/PWGHF/HFCandidateCreator3Prong.cxx +++ b/Analysis/Tasks/PWGHF/HFCandidateCreator3Prong.cxx @@ -173,6 +173,14 @@ struct HFCandidateCreator3ProngMC { } } + // Ξc± → p± K∓ π± + if (result == N3ProngDecays) { + //Printf("Checking Ξc± → p± K∓ π±"); + if (RecoDecay::getMatchedMCRec(particlesMC, std::move(arrayDaughters), 4232, array{+kProton, -kKPlus, +kPiPlus}, true, &sign) > -1) { + result = sign * XicToPKPi; + } + } + rowMCMatchRec(result); } @@ -195,6 +203,14 @@ struct HFCandidateCreator3ProngMC { } } + // Ξc± → p± K∓ π± + if (result == N3ProngDecays) { + //Printf("Checking Λc± → p± K∓ π±"); + if (RecoDecay::isMatchedMCGen(particlesMC, particle, 4122, array{+kProton, -kKPlus, +kPiPlus}, true, &sign)) { + result = sign * LcToPKPi; + } + } + rowMCMatchGen(result); } } diff --git a/Analysis/Tasks/PWGHF/HFTrackIndexSkimsCreator.cxx b/Analysis/Tasks/PWGHF/HFTrackIndexSkimsCreator.cxx index 935f6ad29eff4..f82656e846dd8 100644 --- a/Analysis/Tasks/PWGHF/HFTrackIndexSkimsCreator.cxx +++ b/Analysis/Tasks/PWGHF/HFTrackIndexSkimsCreator.cxx @@ -183,7 +183,8 @@ struct HFTrackIndexSkimsCreator { {"hNCand3ProngVsNTracks", "3-prong candidates preselected;# of selected tracks;# of candidates;entries", {HistType::kTH2F, {{2500, 0., 25000.}, {5000, 0., 500000.}}}}, {"hmassDPlusToPiKPi", "D+ candidates;inv. mass (#pi K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}}, {"hmassLcToPKPi", "Lc candidates;inv. mass (p K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}}, - {"hmassDsToPiKK", "Ds candidates;inv. mass (K K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}}}}; + {"hmassDsToPiKK", "Ds candidates;inv. mass (K K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}}, + {"hmassXicToPKPi", "Xic candidates;inv. mass (p K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 2., 3.}}}}}}; Filter filterSelectTracks = (aod::hf_seltrack::isSelProng > 0); using SelectedTracks = soa::Filtered>; @@ -278,6 +279,12 @@ struct HFTrackIndexSkimsCreator { cut3ProngCPACandMin[DsToPiKK] = configs->mCPADsToPiKKMin; cut3ProngDecLenCandMin[DsToPiKK] = configs->mDecLenDsToPiKKMin; + cut3ProngPtCandMin[XicToPKPi] = configs->mPtXicToPKPiMin; + cut3ProngInvMassCandMin[XicToPKPi] = configs->mInvMassXicToPKPiMin; + cut3ProngInvMassCandMax[XicToPKPi] = configs->mInvMassXicToPKPiMax; + cut3ProngCPACandMin[XicToPKPi] = configs->mCPAXicToPKPiMin; + cut3ProngDecLenCandMin[XicToPKPi] = configs->mDecLenXicToPKPiMin; + bool cutStatus2Prong[n2ProngDecays][nCuts2Prong]; bool cutStatus3Prong[n3ProngDecays][nCuts3Prong]; int nCutStatus2ProngBit = TMath::Power(2, nCuts2Prong) - 1; //bit value for selection status for each 2 prongs candidate where each selection is one bit and they are all set it 1 @@ -295,11 +302,13 @@ struct HFTrackIndexSkimsCreator { arr3Mass1[DPlusToPiKPi] = array{massPi, massK, massPi}; arr3Mass1[LcToPKPi] = array{massProton, massK, massPi}; arr3Mass1[DsToPiKK] = array{massK, massK, massPi}; + arr3Mass1[XicToPKPi] = array{massProton, massK, massPi}; array, n3ProngDecays> arr3Mass2; arr3Mass2[DPlusToPiKPi] = array{massPi, massK, massPi}; arr3Mass2[LcToPKPi] = array{massPi, massK, massProton}; arr3Mass2[DsToPiKK] = array{massPi, massK, massK}; + arr3Mass2[XicToPKPi] = array{massPi, massK, massProton}; double mass2ProngHypo1[n2ProngDecays]; double mass2ProngHypo2[n2ProngDecays]; @@ -654,6 +663,9 @@ struct HFTrackIndexSkimsCreator { if (n3 == DsToPiKK) { registry.get(HIST("hmassDsToPiKK"))->Fill(mass3ProngHypo1[n3]); } + if (n3 == XicToPKPi) { + registry.get(HIST("hmassXicToPKPi"))->Fill(mass3ProngHypo1[n3]); + } } if ((cut3ProngInvMassCandMin[n3] < 0. && cut3ProngInvMassCandMax[n3] <= 0.) || (mass3ProngHypo2[n3] >= cut3ProngInvMassCandMin[n3] && mass3ProngHypo2[n3] < cut3ProngInvMassCandMax[n3])) { mass3ProngHypo2[n3] = RecoDecay::M(arr3Mom, arr3Mass2[n3]); @@ -663,6 +675,9 @@ struct HFTrackIndexSkimsCreator { if (n3 == DsToPiKK) { registry.get(HIST("hmassDsToPiKK"))->Fill(mass3ProngHypo2[n3]); } + if (n3 == XicToPKPi) { + registry.get(HIST("hmassXicToPKPi"))->Fill(mass3ProngHypo2[n3]); + } } } } @@ -821,6 +836,9 @@ struct HFTrackIndexSkimsCreator { if (n3 == DsToPiKK) { registry.get(HIST("hmassDsToPiKK"))->Fill(mass3ProngHypo1[n3]); } + if (n3 == XicToPKPi) { + registry.get(HIST("hmassXicToPKPi"))->Fill(mass3ProngHypo1[n3]); + } } if ((cut3ProngInvMassCandMin[n3] < 0. && cut3ProngInvMassCandMax[n3] <= 0.) || (mass3ProngHypo2[n3] >= cut3ProngInvMassCandMin[n3] && mass3ProngHypo2[n3] < cut3ProngInvMassCandMax[n3])) { mass3ProngHypo2[n3] = RecoDecay::M(arr3Mom, arr3Mass2[n3]); @@ -830,6 +848,9 @@ struct HFTrackIndexSkimsCreator { if (n3 == DsToPiKK) { registry.get(HIST("hmassDsToPiKK"))->Fill(mass3ProngHypo2[n3]); } + if (n3 == XicToPKPi) { + registry.get(HIST("hmassXicToPKPi"))->Fill(mass3ProngHypo2[n3]); + } } } } diff --git a/Analysis/Tasks/PWGHF/HFXicCandidateSelector.cxx b/Analysis/Tasks/PWGHF/HFXicCandidateSelector.cxx new file mode 100644 index 0000000000000..1b2f3585d1bc9 --- /dev/null +++ b/Analysis/Tasks/PWGHF/HFXicCandidateSelector.cxx @@ -0,0 +1,447 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file HFXicCandidateSelector.cxx +/// \brief Xic->pKpi selection task. +/// \note Inspired from HFXicCandidateSelector.cxx +/// +/// \author Mattia Faggin , University and INFN PADOVA + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "AnalysisDataModel/HFSecondaryVertex.h" +#include "AnalysisDataModel/HFCandidateSelectionTables.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::aod::hf_cand_prong3; + +static const int npTBins = 10; +static const int nCutVars = 8; +//temporary until 2D array in configurable is solved - then move to json +// m ptp ptk ptpi DCA sigmavtx dlenght cosp +constexpr double cuts[npTBins][nCutVars] = {{0.400, 0.4, 0.4, 0.4, 0.05, 0.09, 0.005, 0.}, /* pt<1 */ + {0.400, 0.4, 0.4, 0.4, 0.05, 0.09, 0.005, 0.}, /* 1 hfSelXicCandidate; + + Configurable d_FilterPID{"d_FilterPID", true, "Bool to use or not the PID at filtering level"}; + + Configurable d_pTCandMin{"d_pTCandMin", 0., "Lower bound of candidate pT"}; + Configurable d_pTCandMax{"d_pTCandMax", 36., "Upper bound of candidate pT"}; + + Configurable d_pidTPCMinpT{"d_pidTPCMinpT", 0.15, "Lower bound of track pT for TPC PID"}; + Configurable d_pidTPCMaxpT{"d_pidTPCMaxpT", 1., "Upper bound of track pT for TPC PID"}; + Configurable d_pidTOFMinpT{"d_pidTOFMinpT", 0.5, "Lower bound of track pT for TOF PID"}; + Configurable d_pidTOFMaxpT{"d_pidTOFMaxpT", 4., "Upper bound of track pT for TOF PID"}; + + Configurable d_TPCNClsFindablePIDCut{"d_TPCNClsFindablePIDCut", 70., "Lower bound of TPC findable clusters for good PID"}; + Configurable d_nSigmaTPC{"d_nSigmaTPC", 3., "Nsigma cut on TPC only"}; + Configurable d_nSigmaTPCCombined{"d_nSigmaTPCCombined", 5., "Nsigma cut on TPC combined with TOF"}; + Configurable d_nSigmaTOF{"d_nSigmaTOF", 3., "Nsigma cut on TOF only"}; + Configurable d_nSigmaTOFCombined{"d_nSigmaTOFCombined", 5., "Nsigma cut on TOF combined with TPC"}; + + /// Gets corresponding pT bin from cut file array + /// \param candpT is the pT of the candidate + /// \return corresponding bin number of array + template + int getpTBin(T candpT) + { + double pTBins[npTBins + 1] = {0, 1., 2., 3., 4., 5., 6., 8., 12., 24., 36.}; + if (candpT < pTBins[0] || candpT >= pTBins[npTBins]) { + return -1; + } + for (int i = 0; i < npTBins; i++) { + if (candpT < pTBins[i + 1]) { + return i; + } + } + return -1; + } + + /// Selection on goodness of daughter tracks + /// \note should be applied at candidate selection + /// \param track is daughter track + /// \return true if track is good + template + bool daughterSelection(const T& track) + { + if (track.charge() == 0) { + return false; + } + if (track.tpcNClsFound() == 0) { + return false; //is it clusters findable or found - need to check + } + return true; + } + + /// Conjugate independent toplogical cuts + /// \param hfCandProng3 is candidate + /// \return true if candidate passes all cuts + template + bool selectionTopol(const T& hfCandProng3) + { + auto candpT = hfCandProng3.pt(); + int pTBin = getpTBin(candpT); + if (pTBin == -1) { + return false; + } + + if (candpT < d_pTCandMin || candpT >= d_pTCandMax) { + return false; //check that the candidate pT is within the analysis range + } + + if (hfCandProng3.cpa() <= cuts[pTBin][7]) { + return false; //cosine of pointing angle + } + + /* if (hfCandProng3.chi2PCA() > cuts[pTBin][5]) { //candidate DCA + return false; + }*/ + + if (hfCandProng3.decayLength() <= cuts[pTBin][6]) { + return false; + } + return true; + } + + /// Conjugate dependent toplogical cuts + /// \param hfCandProng3 is candidate + /// \param trackProton is the track with the proton hypothesis + /// \param trackPion is the track with the pion hypothesis + /// \param trackKaon is the track with the kaon hypothesis + /// \return true if candidate passes all cuts for the given Conjugate + template + bool selectionTopolConjugate(const T1& hfCandProng3, const T2& trackProton, const T2& trackKaon, const T2& trackPion) + { + + auto candpT = hfCandProng3.pt(); + int pTBin = getpTBin(candpT); + if (pTBin == -1) { + return false; + } + + if (trackProton.pt() < cuts[pTBin][1] || trackKaon.pt() < cuts[pTBin][2] || trackPion.pt() < cuts[pTBin][3]) { + return false; //cut on daughter pT + } + + if (trackProton.globalIndex() == hfCandProng3.index0Id()) { + if (TMath::Abs(InvMassXicpKpi(hfCandProng3) - RecoDecay::getMassPDG(4232)) > cuts[pTBin][0]) { + return false; + } + } else { + if (TMath::Abs(InvMassXicpiKp(hfCandProng3) - RecoDecay::getMassPDG(4232)) > cuts[pTBin][0]) { + return false; + } + } + + return true; + } + + /// Check if track is ok for TPC PID + /// \param track is the track + /// \note function to be expanded + /// \return true if track is ok for TPC PID + template + bool validTPCPID(const T& track) + { + if (TMath::Abs(track.pt()) < d_pidTPCMinpT || TMath::Abs(track.pt()) >= d_pidTPCMaxpT) { + return false; + } + //if (track.TPCNClsFindable() < d_TPCNClsFindablePIDCut) return false; + return true; + } + + /// Check if track is ok for TOF PID + /// \param track is the track + /// \note function to be expanded + /// \return true if track is ok for TOF PID + template + bool validTOFPID(const T& track) + { + if (TMath::Abs(track.pt()) < d_pidTOFMinpT || TMath::Abs(track.pt()) >= d_pidTOFMaxpT) { + return false; + } + return true; + } + + /// Check if track is compatible with given TPC Nsigma cut for a given flavour hypothesis + /// \param track is the track + /// \param nPDG is the flavour hypothesis PDG number + /// \param nSigmaCut is the nsigma threshold to test against + /// \note nPDG=2212 proton nPDG=211 pion nPDG=321 kaon + /// \return true if track satisfies TPC PID hypothesis for given Nsigma cut + template + bool selectionPIDTPC(const T& track, int nPDG, int nSigmaCut) + { + double nSigma = 99.; //arbitarily large value + nPDG = TMath::Abs(nPDG); + if (nPDG == 2212) { + nSigma = track.tpcNSigmaPr(); + } else if (nPDG == 321) { + nSigma = track.tpcNSigmaKa(); + } else if (nPDG == 211) { + nSigma = track.tpcNSigmaPi(); + } else { + return false; + } + + return nSigma < nSigmaCut; + } + + /// Check if track is compatible with given TOF NSigma cut for a given flavour hypothesis + /// \param track is the track + /// \param nPDG is the flavour hypothesis PDG number + /// \param nSigmaCut is the nSigma threshold to test against + /// \note nPDG=2212 proton nPDG=211 pion nPDG=321 kaon + /// \return true if track satisfies TOF PID hypothesis for given NSigma cut + template + bool selectionPIDTOF(const T& track, int nPDG, int nSigmaCut) + { + double nSigma = 99.; //arbitarily large value + nPDG = TMath::Abs(nPDG); + if (nPDG == 2212) { + nSigma = track.tofNSigmaPr(); + } else if (nPDG == 321) { + nSigma = track.tofNSigmaKa(); + } else if (nPDG == 211) { + nSigma = track.tofNSigmaPi(); + } else { + return false; + } + + return nSigma < nSigmaCut; + } + + /// PID selection on daughter track + /// \param track is the daughter track + /// \param nPDG is the PDG code of the flavour hypothesis + /// \note nPDG=2212 nPDG=211 pion nPDG=321 kaon + /// \return 1 if successful PID match, 0 if successful PID rejection, -1 if no PID info + template + int selectionPID(const T& track, int nPDG) + { + int statusTPC = -1; + int statusTOF = -1; + + if (validTPCPID(track)) { + if (!selectionPIDTPC(track, nPDG, d_nSigmaTPC)) { + if (!selectionPIDTPC(track, nPDG, d_nSigmaTPCCombined)) { + statusTPC = 0; //rejected by PID + } else { + statusTPC = 1; //potential to be acceepted if combined with TOF + } + } else { + statusTPC = 2; //positive PID + } + } else { + statusTPC = -1; //no PID info + } + + if (validTOFPID(track)) { + if (!selectionPIDTOF(track, nPDG, d_nSigmaTOF)) { + if (!selectionPIDTOF(track, nPDG, d_nSigmaTOFCombined)) { + statusTOF = 0; //rejected by PID + } else { + statusTOF = 1; //potential to be acceepted if combined with TOF + } + } else { + statusTOF = 2; //positive PID + } + } else { + statusTOF = -1; //no PID info + } + + if (statusTPC == 2 || statusTOF == 2) { + return 1; //what if we have 2 && 0 ? + } else if (statusTPC == 1 && statusTOF == 1) { + return 1; + } else if (statusTPC == 0 || statusTOF == 0) { + return 0; + } else { + return -1; + } + + /* + /// alternative: conservative PID + + /// TPC + if (validTPCPID(track)) { + if (!selectionPIDTPC(track, nPDG, d_nSigmaTPC)) { + statusTPC = 0; //rejected by PID + } else { + statusTPC = 2; //positive PID + } + } + else{ + statusTPC = -1; //no PID info + } + + /// TOF + if (validTOFPID(track)) { + if (!selectionPIDTOF(track, nPDG, d_nSigmaTOF)) { + statusTOF = 0; //rejected by PID + } else { + statusTOF = 2; //positive PID + } + } + else statusTOF = -1; //no PID info + + /// Conservative PID (no need to combine TPC and TOF) + /// case 1. if both PID dectectors info is missing: keep for analysis + /// case 2. if one of the PID detector info is missing: consider the other detector only + /// * the detector accepts the track: keep + /// * the detector does not accept the track: reject + /// case 3. if either at least TPC or TOF accept the track: keep + /// (even if the other would reject it) + if( statusTPC==-1 && statusTOF==-1 ) return -1; // case 1. + else{ + if( statusTPC==0 && statusTOF==0 ) return 0; + else if( statusTPC==-1 && statusTOF==0 ) return 0; + else if( statusTPC==0 && statusTOF==-1 ) return 0; + else return 1; + } + */ + } + + void process(aod::HfCandProng3 const& hfCandProng3s, aod::BigTracksPID const& tracks) + { + int statusXicpKpi, statusXicpiKp; // final selection flag : 0-rejected 1-accepted + bool topolXicpKpi, topolXicpiKp; + int pidXicPKPi, pidXicPiKP, kaonNeg, prot1, prot2, pion1, pion2; // proton, kaonMinus, pionPlus; + + for (auto& hfCandProng3 : hfCandProng3s) { //looping over 3 prong candidates + + statusXicpKpi = 0; + statusXicpiKp = 0; + + if (!(hfCandProng3.hfflag() & 1 << XicToPKPi)) { + hfSelXicCandidate(statusXicpKpi, statusXicpiKp); + continue; + } + + auto trackPos1 = hfCandProng3.index0_as(); //positive daughter (negative for the antiparticles) + auto trackNeg1 = hfCandProng3.index1_as(); //negative daughter (positive for the antiparticles) + auto trackPos2 = hfCandProng3.index2_as(); //positive daughter (negative for the antiparticles) + + topolXicpKpi = true; + topolXicpiKp = true; + pidXicPKPi = -1; + pidXicPiKP = -1; + kaonNeg = -1; + prot1 = -1; + prot2 = -1; + pion1 = -1; + pion2 = -1; + //proton = -1; + //kaonMinus = -1; + //pionPlus = -1; + + // daughter track validity selection + if (!daughterSelection(trackPos1) || !daughterSelection(trackNeg1) || !daughterSelection(trackPos2)) { + hfSelXicCandidate(statusXicpKpi, statusXicpiKp); + continue; + } + + //implement filter bit 4 cut - should be done before this task at the track selection level + + //conjugate independent topological selection + if (!selectionTopol(hfCandProng3)) { + hfSelXicCandidate(statusXicpKpi, statusXicpiKp); + continue; + } + + //conjugate dependent topplogical selection for Xic + + topolXicpKpi = selectionTopolConjugate(hfCandProng3, trackPos1, trackNeg1, trackPos2); + topolXicpiKp = selectionTopolConjugate(hfCandProng3, trackPos2, trackNeg1, trackPos1); + + if (!topolXicpKpi && !topolXicpiKp) { + hfSelXicCandidate(statusXicpKpi, statusXicpiKp); + continue; + } + + /* + proton = selectionPID(trackPos1, 2212); + kaonMinus = selectionPID(trackNeg1, 321); + pionPlus = selectionPID(trackPos2, 211); + + if (proton == 0 || kaonMinus == 0 || pionPlus == 0) { + pidLc = 0; //exclude Lc + } + if (proton == 1 && kaonMinus == 1 && pionPlus == 1) { + pidLc = 1; //accept Lc + } + + if (pidLc == 0) { + hfSelLcCandidate(statusLcpKpi, statusLcpiKp); + continue; + } + + if ((pidLc == -1 || pidLc == 1) && topolLcpKpi) { + statusLcpKpi = 1; //identified as Lc + } + if ((pidLc == -1 || pidLc == 1) && topolLcpiKp) { + statusLcpiKp = 1; //identified as Lc + } + */ + + kaonNeg = selectionPID(trackNeg1, 321); + prot1 = selectionPID(trackPos1, 2212); + prot2 = selectionPID(trackPos2, 2212); + pion1 = selectionPID(trackPos1, 211); + pion2 = selectionPID(trackPos2, 211); + + if (!d_FilterPID) { + /// PID not applied at filtering level + /// set the PID flags to 1 + pidXicPKPi = 1; + pidXicPiKP = 1; + } else { + /// PID applied at filtering level + /// Accept only candidates recognised either at least as pKpi or as piKp + if (kaonNeg == 1) { + if (std::abs(prot1) == 1 && std::abs(pion2) == 1) { + pidXicPKPi = 1; + } + if (std::abs(pion1) == 1 && std::abs(prot2) == 1) { + pidXicPiKP = 1; + } + } + } + + if ((pidXicPKPi == 1) && topolXicpKpi) { + statusXicpKpi = 1; //identified as Xic->pKpi + } + if ((pidXicPiKP == 1) && topolXicpiKp) { + statusXicpiKp = 1; //identified as Xic->piKp + } + + hfSelXicCandidate(statusXicpKpi, statusXicpiKp); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const&) +{ + return WorkflowSpec{ + adaptAnalysisTask("hf-xic-candidate-selector")}; +} diff --git a/Analysis/Tasks/PWGHF/taskXic.cxx b/Analysis/Tasks/PWGHF/taskXic.cxx new file mode 100644 index 0000000000000..21ac204ac30b1 --- /dev/null +++ b/Analysis/Tasks/PWGHF/taskXic.cxx @@ -0,0 +1,166 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file taskXic.cxx +/// \brief Xic± analysis task +/// \note Inspired from taskLc.cxx +/// +/// \author Mattia Faggin , University and INFN PADOVA + +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "AnalysisDataModel/HFSecondaryVertex.h" +#include "AnalysisDataModel/HFCandidateSelectionTables.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::aod::hf_cand_prong3; +using namespace o2::framework::expressions; + +void customize(std::vector& workflowOptions) +{ + ConfigParamSpec optionDoMC{"doMC", VariantType::Bool, false, {"Fill MC histograms."}}; + workflowOptions.push_back(optionDoMC); +} + +#include "Framework/runDataProcessing.h" + +/// Xic± analysis task +struct TaskXic { + HistogramRegistry registry{ + "registry", + {{"hmass", "3-prong candidates;inv. mass (p K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 1.6, 3.1}}}}, + {"hptcand", "3-prong candidates;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, + {"hptprong0", "3-prong candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, + {"hptprong1", "3-prong candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, + {"hptprong2", "3-prong candidates;prong 2 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, + {"hdeclength", "3-prong candidates;decay length (cm);entries", {HistType::kTH1F, {{200, 0., 2.}}}}, + {"hd0Prong0", "3-prong candidates;prong 0 DCAxy to prim. vertex (cm);entries", {HistType::kTH1F, {{100, -1., 1.}}}}, + {"hd0Prong1", "3-prong candidates;prong 1 DCAxy to prim. vertex (cm);entries", {HistType::kTH1F, {{100, -1., 1.}}}}, + {"hd0Prong2", "3-prong candidates;prong 1 DCAxy to prim. vertex (cm);entries", {HistType::kTH1F, {{100, -1., 1.}}}}, + {"hCt", "3-prong candidates;proper lifetime (#Lambda_{c}) * #it{c} (cm);entries", {HistType::kTH1F, {{120, -20., 100.}}}}, + {"hCPA", "3-prong candidates;cosine of pointing angle;entries", {HistType::kTH1F, {{110, -1.1, 1.1}}}}, + {"hEta", "3-prong candidates;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -2., 2.}}}}, + {"hselectionstatus", "3-prong candidates;selection status;entries", {HistType::kTH1F, {{5, -0.5, 4.5}}}}, + {"hImpParErr", "3-prong candidates;impact parameter error (cm);entries", {HistType::kTH1F, {{100, -1., 1.}}}}, + {"hDecLenErr", "3-prong candidates;decay length error (cm);entries", {HistType::kTH1F, {{100, 0., 1.}}}}, + {"hdca2", "3-prong candidates;prong DCA to sec. vertex (cm);entries", {HistType::kTH1F, {{100, 0., 1.}}}}}}; + + Configurable d_selectionFlagXic{"d_selectionFlagXic", 1, "Selection Flag for Xic"}; + Configurable cutEtaCandMax{"cutEtaCandMax", -1., "max. cand. pseudorapidity"}; + + Filter filterSelectCandidates = (aod::hf_selcandidate_xic::isSelXicpKpi >= d_selectionFlagXic || aod::hf_selcandidate_xic::isSelXicpiKp >= d_selectionFlagXic); + + //void process(aod::HfCandProng3 const& candidates) + void process(soa::Filtered> const& candidates) + { + for (auto& candidate : candidates) { + if (!(candidate.hfflag() & 1 << XicToPKPi)) { + continue; + } + if (cutEtaCandMax >= 0. && std::abs(candidate.eta()) > cutEtaCandMax) { + //Printf("Candidate: eta rejection: %g", candidate.eta()); + continue; + } + if (candidate.isSelXicpKpi() >= d_selectionFlagXic) { + registry.fill(HIST("hmass"), InvMassXicpKpi(candidate)); + } + if (candidate.isSelXicpiKp() >= d_selectionFlagXic) { + registry.fill(HIST("hmass"), InvMassXicpiKp(candidate)); + } + registry.fill(HIST("hptcand"), candidate.pt()); + registry.fill(HIST("hptprong0"), candidate.ptProng0()); + registry.fill(HIST("hptprong1"), candidate.ptProng1()); + registry.fill(HIST("hptprong2"), candidate.ptProng2()); + registry.fill(HIST("hdeclength"), candidate.decayLength()); + registry.fill(HIST("hd0Prong0"), candidate.impactParameter0()); + registry.fill(HIST("hd0Prong1"), candidate.impactParameter1()); + registry.fill(HIST("hd0Prong2"), candidate.impactParameter2()); + registry.fill(HIST("hCt"), CtXic(candidate)); + registry.fill(HIST("hCPA"), candidate.cpa()); + registry.fill(HIST("hEta"), candidate.eta()); + registry.fill(HIST("hselectionstatus"), candidate.isSelXicpKpi()); + registry.fill(HIST("hselectionstatus"), candidate.isSelXicpiKp()); + registry.fill(HIST("hImpParErr"), candidate.errorImpactParameter0()); + registry.fill(HIST("hImpParErr"), candidate.errorImpactParameter1()); + registry.fill(HIST("hImpParErr"), candidate.errorImpactParameter2()); + registry.fill(HIST("hDecLenErr"), candidate.errorDecayLength()); + registry.fill(HIST("hDecLenErr"), candidate.chi2PCA()); + } + } +}; + +/// Fills MC histograms. +struct TaskXicMC { + HistogramRegistry registry{ + "registry", + {{"hPtRecSig", "3-prong candidates (rec. matched);#it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, + {"hPtRecBg", "3-prong candidates (rec. unmatched);#it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, + {"hPtGen", "3-prong candidates (gen. matched);#it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, + {"hCPARecSig", "3-prong candidates (rec. matched);cosine of pointing angle;entries", {HistType::kTH1F, {{110, -1.1, 1.1}}}}, + {"hCPARecBg", "3-prong candidates (rec. unmatched);cosine of pointing angle;entries", {HistType::kTH1F, {{110, -1.1, 1.1}}}}, + {"hEtaRecSig", "3-prong candidates (rec. matched);#it{#eta};entries", {HistType::kTH1F, {{100, -2., 2.}}}}, + {"hEtaRecBg", "3-prong candidates (rec. unmatched);#it{#eta};entries", {HistType::kTH1F, {{100, -2., 2.}}}}, + {"hEtaGen", "3-prong candidates (gen. matched);#it{#eta};entries", {HistType::kTH1F, {{100, -2., 2.}}}}}}; + + Configurable d_selectionFlagXic{"d_selectionFlagXic", 1, "Selection Flag for Xic"}; + Configurable d_selectionFlagXicbar{"d_selectionFlagXicbar", 1, "Selection Flag for Xicbar"}; + Configurable cutEtaCandMax{"cutEtaCandMax", -1., "max. cand. pseudorapidity"}; + + Filter filterSelectCandidates = (aod::hf_selcandidate_xic::isSelXicpKpi >= d_selectionFlagXic || aod::hf_selcandidate_xic::isSelXicpiKp >= d_selectionFlagXic); + + void process(soa::Filtered> const& candidates, + soa::Join const& particlesMC) + { + // MC rec. + //Printf("MC Candidates: %d", candidates.size()); + for (auto& candidate : candidates) { + if (!(candidate.hfflag() & 1 << XicToPKPi)) { + continue; + } + if (cutEtaCandMax >= 0. && std::abs(candidate.eta()) > cutEtaCandMax) { + //Printf("MC Rec.: eta rejection: %g", candidate.eta()); + continue; + } + if (std::abs(candidate.flagMCMatchRec()) == XicToPKPi) { + registry.fill(HIST("hPtRecSig"), candidate.pt()); + registry.fill(HIST("hCPARecSig"), candidate.cpa()); + registry.fill(HIST("hEtaRecSig"), candidate.eta()); + } else { + registry.fill(HIST("hPtRecBg"), candidate.pt()); + registry.fill(HIST("hCPARecBg"), candidate.cpa()); + registry.fill(HIST("hEtaRecBg"), candidate.eta()); + } + } + // MC gen. + //Printf("MC Particles: %d", particlesMC.size()); + for (auto& particle : particlesMC) { + if (cutEtaCandMax >= 0. && std::abs(particle.eta()) > cutEtaCandMax) { + //Printf("MC Gen.: eta rejection: %g", particle.eta()); + continue; + } + if (std::abs(particle.flagMCMatchGen()) == XicToPKPi) { + registry.fill(HIST("hPtGen"), particle.pt()); + registry.fill(HIST("hEtaGen"), particle.eta()); + } + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow{ + adaptAnalysisTask("hf-task-xic")}; + const bool doMC = cfgc.options().get("doMC"); + if (doMC) { + workflow.push_back(adaptAnalysisTask("hf-task-xic-mc")); + } + return workflow; +}