diff --git a/Analysis/DataModel/include/Analysis/HFCandidateSelectionTables.h b/Analysis/DataModel/include/Analysis/HFCandidateSelectionTables.h index 9507e32113e4b..fa5ffde065dc4 100644 --- a/Analysis/DataModel/include/Analysis/HFCandidateSelectionTables.h +++ b/Analysis/DataModel/include/Analysis/HFCandidateSelectionTables.h @@ -21,4 +21,13 @@ DECLARE_SOA_COLUMN(IsSelD0bar, isSelD0bar, int); DECLARE_SOA_TABLE(HFSelD0Candidate, "AOD", "HFSELD0CAND", hf_selcandidate::IsSelD0, hf_selcandidate::IsSelD0bar); } // namespace o2::aod +namespace o2::aod +{ +namespace hf_selcandidate_lc +{ +DECLARE_SOA_COLUMN(IsSelLc, isSelLc, int); +} // namespace hf_selcandidate_lc +DECLARE_SOA_TABLE(HFSelLcCandidate, "AOD", "HFSELLcCAND", hf_selcandidate_lc::IsSelLc); +} // namespace o2::aod + #endif // O2_ANALYSIS_HFCANDIDATESELECTIONTABLES_H_ diff --git a/Analysis/DataModel/include/Analysis/HFSecondaryVertex.h b/Analysis/DataModel/include/Analysis/HFSecondaryVertex.h index 4e4437a8daa85..7637fd5571754 100644 --- a/Analysis/DataModel/include/Analysis/HFSecondaryVertex.h +++ b/Analysis/DataModel/include/Analysis/HFSecondaryVertex.h @@ -280,6 +280,32 @@ auto InvMassDPlus(const T& candidate) { return candidate.m(array{RecoDecay::getMassPDG(kPiPlus), RecoDecay::getMassPDG(kKPlus), RecoDecay::getMassPDG(kPiPlus)}); } + +// Lc+ → p K- π+ + +template +auto CtLc(const T& candidate) +{ + return candidate.ct(RecoDecay::getMassPDG(4122)); +} + +template +auto YLc(const T& candidate) +{ + return candidate.y(RecoDecay::getMassPDG(4122)); +} + +template +auto ELc(const T& candidate) +{ + return candidate.e(RecoDecay::getMassPDG(4122)); +} + +template +auto InvMassLc(const T& candidate) +{ + return candidate.m(array{RecoDecay::getMassPDG(kProton), RecoDecay::getMassPDG(kKPlus), RecoDecay::getMassPDG(kPiPlus)}); +} } // namespace hf_cand_prong3 // 3-prong decay candidate table diff --git a/Analysis/Tasks/PWGHF/CMakeLists.txt b/Analysis/Tasks/PWGHF/CMakeLists.txt index f37b15d1ad06a..934b0800e2bb2 100644 --- a/Analysis/Tasks/PWGHF/CMakeLists.txt +++ b/Analysis/Tasks/PWGHF/CMakeLists.txt @@ -28,6 +28,11 @@ o2_add_dpl_workflow(hf-d0-candidate-selector PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing COMPONENT_NAME Analysis) +o2_add_dpl_workflow(hf-lc-candidate-selector + SOURCES HFLcCandidateSelector.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing + COMPONENT_NAME Analysis) + o2_add_dpl_workflow(hf-task-d0 SOURCES taskD0.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing @@ -37,3 +42,8 @@ o2_add_dpl_workflow(hf-task-dplus SOURCES taskDPlus.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing COMPONENT_NAME Analysis) + +o2_add_dpl_workflow(hf-task-lc + SOURCES taskLc.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing + COMPONENT_NAME Analysis) diff --git a/Analysis/Tasks/PWGHF/HFLcCandidateSelector.cxx b/Analysis/Tasks/PWGHF/HFLcCandidateSelector.cxx new file mode 100644 index 0000000000000..c4865dd141d19 --- /dev/null +++ b/Analysis/Tasks/PWGHF/HFLcCandidateSelector.cxx @@ -0,0 +1,356 @@ +// 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 HFD0CandidateSelector.cxx +/// \brief Lc->pKpi selection task. +/// +/// \author Luigi Dello Stritto , CERN +/// \author Nima Zardoshti , CERN + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Analysis/HFSecondaryVertex.h" +#include "Analysis/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.5, 0.1, 0.1, 0.1, 0.05, 0.09, 0.005, 0.}, /* pt<1 */ + {0.5, 0.1, 0.1, 0.1, 0.05, 0.09, 0.005, 0.}, /* 1 hfSelLcCandidate; + + 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.}; + for (int i = 0; i < npTBins; i++) { + if (candpT >= pTBins[i] && 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.dca() > cuts[pTBin][4]) return false; //candidate DCA + if (hfCandProng3.chi2PCA() > cuts[pTBin][5]) { //candidate DCA + return false; + } + + if (hfCandProng3.decayLength() * hfCandProng3.decayLength() < cuts[pTBin][6] * 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 + /// \note trackPion = positive and trackKaon = negative for D0 selection and inverse for D0bar + /// \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; + } + + //invariant mass cut + if (TMath::Abs(InvMassLc(hfCandProng3) - RecoDecay::getMassPDG(4122)) > cuts[pTBin][0]) { + return false; + } + + if (TMath::Abs(trackProton.pt()) < TMath::Abs(cuts[pTBin][1]) || TMath::Abs(trackKaon.pt()) < TMath::Abs(cuts[pTBin][2]) || TMath::Abs(trackPion.pt()) < TMath::Abs(cuts[pTBin][3])) { + return false; //cut on daughter pT + } + + /* if (TMath::Sqrt( trackProton.dcaPrim0()*trackProton.dcaPrim0() + trackKaon.dcaPrim0()*trackKaon.dcaPrim0() + trackPion.dcaPrim0()*trackPion.dcaPrim0() ) > cuts[pTBin][5]) { + return false; //cut on daughter dca - need to add secondary vertex constraint here + }*/ + + 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 = 0.0; + nPDG = TMath::Abs(nPDG); + if (nPDG == 2212) { + nSigma = track.tpcNSigmaPr(); + } else if (nPDG == 321) { + nSigma = track.tpcNSigmaKa(); + } else if (nPDG == 111) { + nSigma = track.tpcNSigmaPi(); + + } else { + return nSigma = 100; //arbitarily large value + } + if (nSigma < nSigmaCut) { + return true; + } else { + return false; + } + } + + /// 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 = 0.0; + nPDG = TMath::Abs(nPDG); + if (nPDG == 2212) { + nSigma = track.tofNSigmaPr(); + } else if (nPDG == 321) { + nSigma = track.tofNSigmaKa(); + } else if (nPDG == 321) { + nSigma = track.tofNSigmaPi(); + } else { + return nSigma = 100; //arbitarily large value + } + if (nSigma < nSigmaCut) { + return true; + } else { + return false; + } + } + + /// 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; + } + } + + void process(aod::HfCandProng3 const& hfCandProng3s, aod::BigTracksPID const& tracks) + { + int statusLc; // final selection flag : 0-rejected 1-accepted + bool topolLc; + int pidLc, proton, kMinus, piPlus; + + for (auto& hfCandProng3 : hfCandProng3s) { //looping over 3 prong candidates + + auto trackPos1 = hfCandProng3.index0_as(); //positive daughter + auto trackNeg1 = hfCandProng3.index1_as(); //negative daughter + auto trackPos2 = hfCandProng3.index2_as(); //positive daughter + + statusLc = 0; + topolLc = true; + pidLc = -1; + proton = -1; + kMinus = -1; + piPlus = -1; + + // daughter track validity selection + if (!daughterSelection(trackPos1) || !daughterSelection(trackNeg1) || !daughterSelection(trackPos2)) { + hfSelLcCandidate(statusLc); + continue; + } + + //implement filter bit 4 cut - should be done before this task at the track selection level + + //conjugate independent topological selection + if (!selectionTopol(hfCandProng3)) { + hfSelLcCandidate(statusLc); + continue; + } + + //conjugate dependent toplogical selection for Lc + topolLc = selectionTopolConjugate(hfCandProng3, trackPos1, trackNeg1, trackPos2); + + if (!topolLc) { + hfSelLcCandidate(statusLc); + continue; + } + + proton = selectionPID(trackPos1, 2212); + kMinus = selectionPID(trackNeg1, 321); + piPlus = selectionPID(trackPos2, 211); + + if (proton == 0 || kMinus == 0 || piPlus == 0) { + pidLc = 0; //exclude Lc + } + if (proton == 1 && kMinus == 1 && piPlus == 1) { + pidLc = 1; //accept Lc + } + + if (pidLc == 0) { + hfSelLcCandidate(statusLc); + continue; + } + + if ((pidLc == -1 || pidLc == 1) && topolLc) { + statusLc = 1; //identified as Lc + } + + hfSelLcCandidate(statusLc); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const&) +{ + return WorkflowSpec{ + adaptAnalysisTask("hf-lc-candidate-selector")}; +} diff --git a/Analysis/Tasks/PWGHF/taskLc.cxx b/Analysis/Tasks/PWGHF/taskLc.cxx new file mode 100644 index 0000000000000..ca721dbc97e9e --- /dev/null +++ b/Analysis/Tasks/PWGHF/taskLc.cxx @@ -0,0 +1,84 @@ +// 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 taskDPlus.cxx +/// \brief D± analysis task +/// \note Extended from taskD0 +/// +/// \author Vít Kučera , CERN + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Analysis/HFSecondaryVertex.h" +#include "Analysis/HFCandidateSelectionTables.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::aod::hf_cand_prong3; +using namespace o2::framework::expressions; + +/// Lc+ analysis task +struct TaskLc { + OutputObj hmass{TH1F("hmass", "3-prong candidates;inv. mass (#pi K #pi) (GeV/#it{c}^{2});entries", 500, 1.6, 3.1)}; + OutputObj hptcand{TH1F("hptcand", "3-prong candidates;candidate #it{p}_{T} (GeV/#it{c});entries", 100, 0., 10.)}; + OutputObj hptprong0{TH1F("hptprong0", "3-prong candidates;prong 0 #it{p}_{T} (GeV/#it{c});entries", 100, 0., 10.)}; + OutputObj hptprong1{TH1F("hptprong1", "3-prong candidates;prong 1 #it{p}_{T} (GeV/#it{c});entries", 100, 0., 10.)}; + OutputObj hptprong2{TH1F("hptprong2", "3-prong candidates;prong 2 #it{p}_{T} (GeV/#it{c});entries", 100, 0., 10.)}; + OutputObj hdeclength{TH1F("declength", "3-prong candidates;decay length (cm);entries", 200, 0., 2.)}; + OutputObj hd0Prong0{TH1F("hd0Prong0", "3-prong candidates;prong 0 DCAxy to prim. vertex (cm);entries", 100, -1., 1.)}; + OutputObj hd0Prong1{TH1F("hd0Prong1", "3-prong candidates;prong 1 DCAxy to prim. vertex (cm);entries", 100, -1., 1.)}; + OutputObj hd0Prong2{TH1F("hd0Prong2", "3-prong candidates;prong 2 DCAxy to prim. vertex (cm);entries", 100, -1., 1.)}; + OutputObj hCt{TH1F("hCt", "3-prong candidates;proper lifetime (D^{#pm}) * #it{c} (cm);entries", 120, -20., 100.)}; + OutputObj hCPA{TH1F("hCPA", "3-prong candidates;cosine of pointing angle;entries", 110, -1.1, 1.1)}; + OutputObj hEta{TH1F("hEta", "3-prong candidates;candidate #it{#eta};entries", 100, -2., 2.)}; + //OutputObj hselectionstatus{TH1F("selectionstatus", "3-prong candidates;selection status;entries", 5, -0.5, 4.5)}; + OutputObj hImpParErr{TH1F("hImpParErr", "3-prong candidates;impact parameter error (cm);entries", 100, -1., 1.)}; + OutputObj hDecLenErr{TH1F("hDecLenErr", "3-prong candidates;decay length error (cm);entries", 100, 0., 1.)}; + //OutputObj hdca{TH1F("hdca", "3-prong candidates;prongs DCA to prim. vertex (cm);entries", 100, -1., 1.)}; + OutputObj hdca2{TH1F("hdca2", "3-prong candidates;prongs DCA to sec. vertex (cm);entries", 100, 0., 1.)}; + + Configurable d_selectionFlagLc{"d_selectionFlagLc", 1, "Selection Flag for Lc"}; + + Filter filterSelectCandidates = (aod::hf_selcandidate_lc::isSelLc >= d_selectionFlagLc); + + //void process(aod::HfCandProng3 const& candidates) + void process(soa::Filtered> const& candidates) + { + for (auto& candidate : candidates) { + if (candidate.isSelLc() >= d_selectionFlagLc) { + hmass->Fill(InvMassLc(candidate)); + } + hptcand->Fill(candidate.pt()); + hptprong0->Fill(candidate.ptProng0()); + hptprong1->Fill(candidate.ptProng1()); + hptprong2->Fill(candidate.ptProng2()); + hdeclength->Fill(candidate.decayLength()); + hd0Prong0->Fill(candidate.impactParameter0()); + hd0Prong1->Fill(candidate.impactParameter1()); + hd0Prong2->Fill(candidate.impactParameter2()); + hCt->Fill(CtLc(candidate)); + hCPA->Fill(candidate.cpa()); + hEta->Fill(candidate.eta()); + //hselectionstatus->Fill(candidate.isSelLc()); + hImpParErr->Fill(candidate.errorImpactParameter0()); + hImpParErr->Fill(candidate.errorImpactParameter1()); + hImpParErr->Fill(candidate.errorImpactParameter2()); + hDecLenErr->Fill(candidate.errorDecayLength()); + // hdca->Fill(candidate.dca()); + hdca2->Fill(candidate.chi2PCA()); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const&) +{ + return WorkflowSpec{ + adaptAnalysisTask("hf-task-lc")}; +}