diff --git a/Analysis/Core/include/AnalysisCore/RecoDecay.h b/Analysis/Core/include/AnalysisCore/RecoDecay.h index 740d4ec43247c..369f2dae6ca24 100644 --- a/Analysis/Core/include/AnalysisCore/RecoDecay.h +++ b/Analysis/Core/include/AnalysisCore/RecoDecay.h @@ -717,7 +717,7 @@ class RecoDecay } // Check that the number of direct daughters is not larger than the number of expected final daughters. if (indexDaughterFirst > -1 && indexDaughterLast > -1 && indexDaughterLast - indexDaughterFirst + 1 > N) { - //Printf("MC Rec: Rejected: too many direct daughters: %d (expected %d final)", indexDaughterLast - indexDaughterFirst + 1, N); + //Printf("MC Rec: Rejected: too many direct daughters: %d (expected %ld final)", indexDaughterLast - indexDaughterFirst + 1, N); return -1; } // Get the list of actual final daughters. @@ -729,7 +729,7 @@ class RecoDecay //printf("\n"); // Check whether the number of actual final daughters is equal to the number of expected final daughters (i.e. the number of provided prongs). if (arrAllDaughtersIndex.size() != N) { - //Printf("MC Rec: Rejected: incorrect number of final daughters: %d (expected %d)", arrAllDaughtersIndex.size(), N); + //Printf("MC Rec: Rejected: incorrect number of final daughters: %ld (expected %ld)", arrAllDaughtersIndex.size(), N); return -1; } } @@ -837,19 +837,19 @@ class RecoDecay } // Check that the number of direct daughters is not larger than the number of expected final daughters. if (indexDaughterFirst > -1 && indexDaughterLast > -1 && indexDaughterLast - indexDaughterFirst + 1 > N) { - //Printf("MC Gen: Rejected: too many direct daughters: %d (expected %d final)", indexDaughterLast - indexDaughterFirst + 1, N); + //Printf("MC Gen: Rejected: too many direct daughters: %d (expected %ld final)", indexDaughterLast - indexDaughterFirst + 1, N); return false; } // Get the list of actual final daughters. getDaughters(particlesMC, candidate, &arrAllDaughtersIndex, arrPDGDaughters, depthMax); - //printf("MC Gen: Mother %d has %d final daughters:", candidate.globalIndex(), arrAllDaughtersIndex.size()); + //printf("MC Gen: Mother %ld has %ld final daughters:", candidate.globalIndex(), arrAllDaughtersIndex.size()); //for (auto i : arrAllDaughtersIndex) { // printf(" %d", i); //} //printf("\n"); // Check whether the number of final daughters is equal to the required number. if (arrAllDaughtersIndex.size() != N) { - //Printf("MC Gen: Rejected: incorrect number of final daughters: %d (expected %d)", arrAllDaughtersIndex.size(), N); + //Printf("MC Gen: Rejected: incorrect number of final daughters: %ld (expected %ld)", arrAllDaughtersIndex.size(), N); return false; } // Check daughters' PDG codes. diff --git a/Analysis/DataModel/include/AnalysisDataModel/HFCandidateSelectionTables.h b/Analysis/DataModel/include/AnalysisDataModel/HFCandidateSelectionTables.h index 8a664269bb901..50c5057e30ca8 100644 --- a/Analysis/DataModel/include/AnalysisDataModel/HFCandidateSelectionTables.h +++ b/Analysis/DataModel/include/AnalysisDataModel/HFCandidateSelectionTables.h @@ -42,6 +42,13 @@ DECLARE_SOA_COLUMN(IsSelJpsiToEE, isSelJpsiToEE, int); //! } // namespace hf_selcandidate_jpsi DECLARE_SOA_TABLE(HFSelJpsiToEECandidate, "AOD", "HFSELJPSICAND", //! hf_selcandidate_jpsi::IsSelJpsiToEE); +namespace hf_selcandidate_lc_k0sp +{ +DECLARE_SOA_COLUMN(IsSelLcK0sP, isSelLcK0sP, int); +} // namespace hf_selcandidate_lc_k0sp +DECLARE_SOA_TABLE(HFSelLcK0sPCandidate, "AOD", "HFSELLCK0SPCAND", //! + hf_selcandidate_lc_k0sp::IsSelLcK0sP); + } // namespace o2::aod namespace o2::aod diff --git a/Analysis/DataModel/include/AnalysisDataModel/HFSecondaryVertex.h b/Analysis/DataModel/include/AnalysisDataModel/HFSecondaryVertex.h index 8d4a61c89cf4a..e8b6566021563 100644 --- a/Analysis/DataModel/include/AnalysisDataModel/HFSecondaryVertex.h +++ b/Analysis/DataModel/include/AnalysisDataModel/HFSecondaryVertex.h @@ -21,6 +21,7 @@ #include "AnalysisCore/RecoDecay.h" #include "AnalysisCore/HFSelectorCuts.h" #include "AnalysisDataModel/PID/PIDResponse.h" +#include "AnalysisDataModel/StrangenessTables.h" using namespace o2::analysis; @@ -52,6 +53,7 @@ DECLARE_SOA_INDEX_COLUMN_FULL(Index0, index0, int, Tracks, "_0"); //! DECLARE_SOA_INDEX_COLUMN_FULL(Index1, index1, int, Tracks, "_1"); //! DECLARE_SOA_INDEX_COLUMN_FULL(Index2, index2, int, Tracks, "_2"); //! DECLARE_SOA_INDEX_COLUMN_FULL(Index3, index3, int, Tracks, "_3"); //! +DECLARE_SOA_INDEX_COLUMN_FULL(IndexV0, indexV0, int, aod::V0Datas, "_V0"); //! DECLARE_SOA_COLUMN(HFflag, hfflag, uint8_t); //! DECLARE_SOA_COLUMN(D0ToKPiFlag, d0ToKPiFlag, uint8_t); //! @@ -68,6 +70,11 @@ DECLARE_SOA_TABLE(HfTrackIndexProng2, "AOD", "HFTRACKIDXP2", //! hf_track_index::Index1Id, hf_track_index::HFflag); +DECLARE_SOA_TABLE(HfTrackIndexCasc, "AOD", "HFTRACKIDXCASC", //! + hf_track_index::Index0Id, + hf_track_index::IndexV0Id, + hf_track_index::HFflag); + DECLARE_SOA_TABLE(HfCutStatusProng2, "AOD", "HFCUTSTATUSP2", //! hf_track_index::D0ToKPiFlag, hf_track_index::JpsiToEEFlag); @@ -357,6 +364,105 @@ DECLARE_SOA_TABLE(HfCandProng2MCGen, "AOD", "HFCANDP2MCGEN", //! hf_cand_prong2::FlagMCMatchGen, hf_cand_prong2::OriginMCGen); +// cascade decay candidate table + +namespace hf_cand_casc +{ +DECLARE_SOA_EXPRESSION_COLUMN(Px, px, //! + float, 1.f * aod::hf_cand::pxProng0 + 1.f * aod::hf_cand::pxProng1); +DECLARE_SOA_EXPRESSION_COLUMN(Py, py, //! + float, 1.f * aod::hf_cand::pyProng0 + 1.f * aod::hf_cand::pyProng1); +DECLARE_SOA_EXPRESSION_COLUMN(Pz, pz, //! + float, 1.f * aod::hf_cand::pzProng0 + 1.f * aod::hf_cand::pzProng1); +//DECLARE_SOA_DYNAMIC_COLUMN(M, m, [](float px0, float py0, float pz0, float px1, float py1, float pz1, const array& m) { return RecoDecay::M(array{array{px0, py0, pz0}, array{px1, py1, pz1}}, m); }); +DECLARE_SOA_DYNAMIC_COLUMN(PtV0Pos, ptV0Pos, //! + [](float px, float py) { return RecoDecay::Pt(px, py); }); +DECLARE_SOA_DYNAMIC_COLUMN(PtV0Neg, ptV0Neg, //! + [](float px, float py) { return RecoDecay::Pt(px, py); }); +DECLARE_SOA_COLUMN(FlagMCMatchRec, flagMCMatchRec, int8_t); //! reconstruction level +DECLARE_SOA_COLUMN(FlagMCMatchGen, flagMCMatchGen, int8_t); //! generator level + +template +auto InvMassLcToK0sP(const T& candidate) +{ + return candidate.m(array{RecoDecay::getMassPDG(kK0Short), RecoDecay::getMassPDG(kProton)}); // first daughter is K0s +} + +template +auto InvMassGamma(const T& candidate) +{ + return candidate.m(array{RecoDecay::getMassPDG(kElectron), RecoDecay::getMassPDG(kElectron)}); +} + +} // namespace hf_cand_casc + +DECLARE_SOA_TABLE(HfCandCascBase, "AOD", "HFCANDCASCBASE", //! + // general columns + HFCAND_COLUMNS, + // cascade specific columns + hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, + hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1, + hf_cand::ImpactParameter0, hf_cand::ImpactParameter1, + hf_cand::ErrorImpactParameter0, hf_cand::ErrorImpactParameter1, + hf_track_index::Index0Id, + hf_track_index::IndexV0Id, // V0 index + hf_track_index::HFflag, + // V0 + v0data::X, v0data::Y, v0data::Z, + v0data::PosTrackId, v0data::NegTrackId, // indices of V0 tracks in FullTracks table + v0data::PxPos, v0data::PyPos, v0data::PzPos, v0data::PxNeg, v0data::PyNeg, v0data::PzNeg, + v0data::DCAV0Daughters, + v0data::DCAPosToPV, // this is the impact param wrt prim vtx in xy! + v0data::DCANegToPV, // this is the impact param wrt prim vtx in xy! + /* dynamic columns */ + hf_cand_prong2::M, + hf_cand_prong2::M2, + hf_cand_prong2::ImpactParameterProduct, + hf_cand_prong2::CosThetaStar, + hf_cand_prong2::ImpactParameterProngSqSum, + /* dynamic columns that use candidate momentum components */ + hf_cand::Pt, + hf_cand::Pt2, + hf_cand::P, + hf_cand::P2, + hf_cand::PVector, + hf_cand::CPA, + hf_cand::CPAXY, + hf_cand::Ct, + hf_cand::ImpactParameterXY, + hf_cand_prong2::MaxNormalisedDeltaIP, + hf_cand::Eta, + hf_cand::Phi, + hf_cand::Y, + hf_cand::E, + hf_cand::E2, + // dynamic columns from V0 + hf_cand_casc::PtV0Pos, // pT of positive V0 daughter + hf_cand_casc::PtV0Neg, // pT of negative V0 daughter + v0data::V0Radius, + v0data::V0CosPA, + v0data::MLambda, + v0data::MAntiLambda, + v0data::MK0Short); +/*, + v0data::MLambda, + v0data::MAntiLambda, + v0data::MK0Short); +*/ +// extended table with expression columns that can be used as arguments of dynamic columns +DECLARE_SOA_EXTENDED_TABLE_USER(HfCandCascExt, HfCandCascBase, "HFCANDCASCEXT", //! + hf_cand_casc::Px, hf_cand_casc::Py, hf_cand_casc::Pz); + +using HfCandCascade = HfCandCascExt; + +// table with results of reconstruction level MC matching for Cascade +DECLARE_SOA_TABLE(HfCandCascadeMCRec, "AOD", "HFCANDCASCMCREC", //! + hf_cand_casc::FlagMCMatchRec); + +// table with results of generator level MC matching +DECLARE_SOA_TABLE(HfCandCascadeMCGen, "AOD", "HFCANDCASCMCGEN", //! + hf_cand_casc::FlagMCMatchGen); + // specific 3-prong decay properties namespace hf_cand_prong3 { diff --git a/Analysis/Tasks/CMakeLists.txt b/Analysis/Tasks/CMakeLists.txt index 970e02bdb2dd0..39b3f394b20ef 100644 --- a/Analysis/Tasks/CMakeLists.txt +++ b/Analysis/Tasks/CMakeLists.txt @@ -14,6 +14,7 @@ add_subdirectory(PWGHF) add_subdirectory(PWGJE) add_subdirectory(PWGLF) add_subdirectory(PWGUD) +add_subdirectory(Utils) add_subdirectory(SkimmingTutorials) add_subdirectory(PID) diff --git a/Analysis/Tasks/PWGHF/CMakeLists.txt b/Analysis/Tasks/PWGHF/CMakeLists.txt index 6fb285562c45d..37517c72640bf 100644 --- a/Analysis/Tasks/PWGHF/CMakeLists.txt +++ b/Analysis/Tasks/PWGHF/CMakeLists.txt @@ -25,7 +25,7 @@ o2_add_dpl_workflow(qa-simple o2_add_dpl_workflow(hf-track-index-skims-creator SOURCES HFTrackIndexSkimsCreator.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing ROOT::EG + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing O2::AnalysisTasksUtils ROOT::EG COMPONENT_NAME Analysis) o2_add_dpl_workflow(hf-candidate-creator-2prong @@ -38,6 +38,11 @@ o2_add_dpl_workflow(hf-tree-creator-d0-tokpi PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing ROOT::EG COMPONENT_NAME Analysis) +o2_add_dpl_workflow(hf-candidate-creator-cascade + SOURCES HFCandidateCreatorCascade.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing O2::AnalysisTasksUtils ROOT::EG + COMPONENT_NAME Analysis) + o2_add_dpl_workflow(hf-candidate-creator-3prong SOURCES HFCandidateCreator3Prong.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing ROOT::EG @@ -73,6 +78,11 @@ o2_add_dpl_workflow(hf-xic-topkpi-candidate-selector PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing COMPONENT_NAME Analysis) +o2_add_dpl_workflow(hf-lc-tok0sp-candidate-selector + SOURCES HFLcK0sPCandidateSelector.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing O2::AnalysisTasksUtils + COMPONENT_NAME Analysis) + o2_add_dpl_workflow(hf-task-d0 SOURCES taskD0.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing @@ -112,3 +122,9 @@ o2_add_dpl_workflow(hf-mc-validation SOURCES HFMCValidation.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing COMPONENT_NAME Analysis) + +o2_add_dpl_workflow(hf-task-lc-tok0sp + SOURCES taskLcK0sP.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing O2::AnalysisTasksUtils + COMPONENT_NAME Analysis) + diff --git a/Analysis/Tasks/PWGHF/HFCandidateCreatorCascade.cxx b/Analysis/Tasks/PWGHF/HFCandidateCreatorCascade.cxx new file mode 100644 index 0000000000000..951de55a845f2 --- /dev/null +++ b/Analysis/Tasks/PWGHF/HFCandidateCreatorCascade.cxx @@ -0,0 +1,306 @@ +// 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 HFCandidateCreatorCascade.cxx +/// \brief Reconstruction of heavy-flavour cascade decay candidates +/// + +#include "Framework/AnalysisTask.h" +#include "DetectorsVertexing/DCAFitterN.h" +#include "AnalysisDataModel/HFSecondaryVertex.h" +#include "AnalysisCore/trackUtilities.h" +#include "ReconstructionDataFormats/DCA.h" +#include "ReconstructionDataFormats/V0.h" +#include "AnalysisTasksUtils/UtilsDebugLcK0Sp.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::aod::hf_cand_prong2; + +void customize(std::vector& workflowOptions) +{ + ConfigParamSpec optionDoMC{"doMC", VariantType::Bool, false, {"Perform MC matching."}}; + workflowOptions.push_back(optionDoMC); +} + +#include "Framework/runDataProcessing.h" + +//#define MY_DEBUG + +#ifdef MY_DEBUG +using MyBigTracks = aod::BigTracksMC; +#define MY_DEBUG_MSG(condition, cmd) \ + if (condition) { \ + cmd; \ + } +#else +using MyBigTracks = aod::BigTracks; +#define MY_DEBUG_MSG(condition, cmd) +#endif + +/// Reconstruction of heavy-flavour cascade decay candidates +struct HFCandidateCreatorCascade { + + Produces rowCandidateBase; + + Configurable bZ{"bZ", 5., "magnetic field"}; + Configurable propDCA{"propDCA", true, "create tracks version propagated to PCA"}; + Configurable maxR{"maxR", 200., "reject PCA's above this radius"}; + Configurable maxDZIni{"maxDZIni", 4., "reject (if>0) PCA candidate if tracks DZ exceeds threshold"}; + Configurable minParamChange{"minParamChange", 1.e-3, "stop iterations if largest change of any X is smaller than this"}; + Configurable minRelChi2Change{"minRelChi2Change", 0.9, "stop iterations is chi2/chi2old > this"}; + Configurable doValPlots{"doValPlots", true, "do validation plots"}; + + // for debugging +#ifdef MY_DEBUG + Configurable> indexK0Spos{"indexK0Spos", {729, 2866, 4754, 5457, 6891, 7824, 9243, 9810}, "indices of K0S positive daughters, for debug"}; + Configurable> indexK0Sneg{"indexK0Sneg", {730, 2867, 4755, 5458, 6892, 7825, 9244, 9811}, "indices of K0S negative daughters, for debug"}; + Configurable> indexProton{"indexProton", {717, 2810, 4393, 5442, 6769, 7793, 9002, 9789}, "indices of protons, for debug"}; +#endif + + OutputObj hmass2{TH1F("hmass2", "2-prong candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", 500, 0., 5.)}; + OutputObj hCovPVXX{TH1F("hCovPVXX", "2-prong candidates;XX element of cov. matrix of prim. vtx. position (cm^{2});entries", 100, 0., 1.e-4)}; + OutputObj hCovSVXX{TH1F("hCovSVXX", "2-prong candidates;XX element of cov. matrix of sec. vtx. position (cm^{2});entries", 100, 0., 0.2)}; + + double massP = RecoDecay::getMassPDG(kProton); + double massK0s = RecoDecay::getMassPDG(kK0Short); + double massPi = RecoDecay::getMassPDG(kPiPlus); + double massLc = RecoDecay::getMassPDG(pdg::Code::kLambdaCPlus); + double mass2K0sP{0.}; + + void process(aod::Collisions const& collisions, + aod::HfTrackIndexCasc const& rowsTrackIndexCasc, + MyBigTracks const& tracks, + aod::V0Datas const& V0s +#ifdef MY_DEBUG + , + aod::McParticles& mcParticles +#endif + ) + { + // 2-prong vertex fitter + o2::vertexing::DCAFitterN<2> df; + df.setBz(bZ); + df.setPropagateToPCA(propDCA); + df.setMaxR(maxR); + df.setMaxDZIni(maxDZIni); + df.setMinParamChange(minParamChange); + df.setMinRelChi2Change(minRelChi2Change); + df.setUseAbsDCA(true); + + // loop over pairs of track indeces + for (const auto& casc : rowsTrackIndexCasc) { + + const auto& bach = casc.index0_as(); + const auto& v0 = casc.indexV0_as(); + const auto& trackV0DaughPos = v0.posTrack_as(); + const auto& trackV0DaughNeg = v0.negTrack_as(); + +#ifdef MY_DEBUG + auto indexBach = bach.mcParticleId(); + auto indexV0DaughPos = trackV0DaughPos.mcParticleId(); + auto indexV0DaughNeg = trackV0DaughNeg.mcParticleId(); + bool isLc = isLcK0SpFunc(indexBach, indexV0DaughPos, indexV0DaughNeg, indexProton, indexK0Spos, indexK0Sneg); +#endif + + MY_DEBUG_MSG(isLc, LOG(INFO) << "Processing the Lc with proton " << indexBach << " trackV0DaughPos " << indexV0DaughPos << " trackV0DaughNeg " << indexV0DaughNeg); + + auto trackParCovBach = getTrackParCov(bach); + auto trackParCovV0DaughPos = getTrackParCov(trackV0DaughPos); // check that MyBigTracks does not need TracksExtended! + auto trackParCovV0DaughNeg = getTrackParCov(trackV0DaughNeg); // check that MyBigTracks does not need TracksExtended! + trackParCovV0DaughPos.propagateTo(v0.posX(), bZ); // propagate the track to the X closest to the V0 vertex + trackParCovV0DaughNeg.propagateTo(v0.negX(), bZ); // propagate the track to the X closest to the V0 vertex + const std::array vertexV0 = {v0.x(), v0.y(), v0.z()}; + const std::array momentumV0 = {v0.px(), v0.py(), v0.pz()}; + // we build the neutral track to then build the cascade + auto trackV0 = o2::dataformats::V0(vertexV0, momentumV0, {0, 0, 0, 0, 0, 0}, trackParCovV0DaughPos, trackParCovV0DaughNeg, {0, 0}, {0, 0}); // build the V0 track (indices for v0 daughters set to 0 for now) + + auto collision = bach.collision(); + + // reconstruct the cascade secondary vertex + if (df.process(trackV0, trackParCovBach) == 0) { + MY_DEBUG_MSG(isLc, LOG(INFO) << "Vertexing failed for Lc candidate"); + // if (isLc) { + // LOG(INFO) << "Vertexing failed for Lc with proton " << indexBach << " trackV0DaughPos " << indexV0DaughPos << " trackV0DaughNeg " << indexV0DaughNeg; + //} + continue; + } else { + //LOG(INFO) << "Vertexing succeeded for Lc candidate"; + } + + const auto& secondaryVertex = df.getPCACandidate(); + auto chi2PCA = df.getChi2AtPCACandidate(); + auto covMatrixPCA = df.calcPCACovMatrix().Array(); + hCovSVXX->Fill(covMatrixPCA[0]); // FIXME: Calculation of errorDecayLength(XY) gives wrong values without this line. + // do I have to call "df.propagateTracksToVertex();"? + auto trackParVarV0 = df.getTrack(0); + auto trackParVarBach = df.getTrack(1); + + // get track momenta + array pVecV0; + array pVecBach; + trackParVarV0.getPxPyPzGlo(pVecV0); + trackParVarBach.getPxPyPzGlo(pVecBach); + + // get track impact parameters + // This modifies track momenta! + auto primaryVertex = getPrimaryVertex(collision); + auto covMatrixPV = primaryVertex.getCov(); + hCovPVXX->Fill(covMatrixPV[0]); + o2::dataformats::DCA impactParameterV0; + o2::dataformats::DCA impactParameterBach; + trackParVarV0.propagateToDCA(primaryVertex, bZ, &impactParameterV0); // we do this wrt the primary vtx + trackParVarBach.propagateToDCA(primaryVertex, bZ, &impactParameterBach); + + // get uncertainty of the decay length + double phi, theta; + getPointDirection(array{collision.posX(), collision.posY(), collision.posZ()}, secondaryVertex, phi, theta); + auto errorDecayLength = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta)); + auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.)); + + // fill candidate table rows + MY_DEBUG_MSG(isLc, LOG(INFO) << "IT IS A Lc! Filling for Lc with proton " << indexBach << " trackV0DaughPos " << indexV0DaughPos << " trackV0DaughNeg " << indexV0DaughNeg); + rowCandidateBase(collision.globalIndex(), + collision.posX(), collision.posY(), collision.posZ(), + secondaryVertex[0], secondaryVertex[1], secondaryVertex[2], + errorDecayLength, errorDecayLengthXY, + chi2PCA, + pVecBach[0], pVecBach[1], pVecBach[2], + pVecV0[0], pVecV0[1], pVecV0[2], + impactParameterBach.getY(), impactParameterV0.getY(), + std::sqrt(impactParameterBach.getSigmaY2()), std::sqrt(impactParameterV0.getSigmaY2()), + casc.index0Id(), casc.indexV0Id(), + casc.hfflag(), + v0.x(), v0.y(), v0.z(), + //v0.posTrack(), v0.negTrack(), // why this was not fine? + trackV0DaughPos.globalIndex(), trackV0DaughNeg.globalIndex(), + v0.pxpos(), v0.pypos(), v0.pzpos(), + v0.pxneg(), v0.pyneg(), v0.pzneg(), + v0.dcaV0daughters(), + v0.dcapostopv(), + v0.dcanegtopv()); + + // fill histograms + if (doValPlots) { + // calculate invariant masses + mass2K0sP = RecoDecay::M(array{pVecBach, pVecV0}, array{massP, massK0s}); + hmass2->Fill(mass2K0sP); + } + } + } +}; + +/// Extends the base table with expression columns. +struct HFCandidateCreatorCascadeExpressions { + Spawns rowCandidateCasc; + void init(InitContext const&) {} +}; + +//___________________________________________________________________________________________ + +/// Performs MC matching. +struct HFCandidateCreatorCascadeMC { + Produces rowMCMatchRec; + Produces rowMCMatchGen; + +#ifdef MY_DEBUG + Configurable> indexK0Spos{"indexK0Spos", {729, 2866, 4754, 5457, 6891, 7824, 9243, 9810}, "indices of K0S positive daughters, for debug"}; + Configurable> indexK0Sneg{"indexK0Sneg", {730, 2867, 4755, 5458, 6892, 7825, 9244, 9811}, "indices of K0S negative daughters, for debug"}; + Configurable> indexProton{"indexProton", {717, 2810, 4393, 5442, 6769, 7793, 9002, 9789}, "indices of protons, for debug"}; +#endif + + void process(aod::HfCandCascade const& candidates, + aod::BigTracksMC const& tracks, + aod::McParticles const& particlesMC) + { + int8_t sign = 0; + std::vector arrDaughLcIndex; + std::array arrDaughLcPDG; + std::array arrDaughLcPDGRef = {2212, 211, -211}; + + // Match reconstructed candidates. + for (auto& candidate : candidates) { + + const auto& bach = candidate.index0_as(); + const auto& trackV0DaughPos = candidate.posTrack_as(); + const auto& trackV0DaughNeg = candidate.negTrack_as(); + + auto arrayDaughtersV0 = array{trackV0DaughPos, trackV0DaughNeg}; + auto arrayDaughtersLc = array{bach, trackV0DaughPos, trackV0DaughNeg}; + + // First we check the K0s + LOG(DEBUG) << "\n"; + LOG(DEBUG) << "Checking MC for candidate!"; + LOG(DEBUG) << "Looking for K0s"; + auto indexV0DaughPos = trackV0DaughPos.mcParticleId(); + auto indexV0DaughNeg = trackV0DaughNeg.mcParticleId(); + auto indexBach = bach.mcParticleId(); + +#ifdef MY_DEBUG + bool isLc = isLcK0SpFunc(indexBach, indexV0DaughPos, indexV0DaughNeg, indexProton, indexK0Spos, indexK0Sneg); + bool isK0SfromLc = isK0SfromLcFunc(indexV0DaughPos, indexV0DaughNeg, indexK0Spos, indexK0Sneg); +#endif + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "correct K0S in the Lc daughters: posTrack --> " << indexV0DaughPos << ", negTrack --> " << indexV0DaughNeg); + + //if (isLc) { + RecoDecay::getMatchedMCRec(particlesMC, arrayDaughtersV0, kK0Short, array{+kPiPlus, -kPiPlus}, true, &sign, 1); // does it matter the "acceptAntiParticle" in the K0s case? In principle, there is no anti-K0s + + if (sign != 0) { // we have already positively checked the K0s + // then we check the Lc + MY_DEBUG_MSG(sign, LOG(INFO) << "K0S was correct! now we check the Lc"); + MY_DEBUG_MSG(sign, LOG(INFO) << "index proton = " << indexBach); + RecoDecay::getMatchedMCRec(particlesMC, arrayDaughtersLc, pdg::Code::kLambdaCPlus, array{+kProton, +kPiPlus, -kPiPlus}, true, &sign, 3); // 3-levels Lc --> p + K0 --> p + K0s --> p + pi+ pi- + MY_DEBUG_MSG(sign, LOG(INFO) << "Lc found with sign " << sign; printf("\n")); + } + + rowMCMatchRec(sign); + } + //} + + // Match generated particles. + for (auto& particle : particlesMC) { + // checking if I have a Lc --> K0S + p + RecoDecay::isMatchedMCGen(particlesMC, particle, pdg::Code::kLambdaCPlus, array{+kProton, +kK0Short}, true, &sign, 2); + if (sign != 0) { + MY_DEBUG_MSG(sign, LOG(INFO) << "Lc in K0S p"); + arrDaughLcIndex.clear(); + // checking that the final daughters (decay depth = 3) are p, pi+, pi- + RecoDecay::getDaughters(particlesMC, particle, &arrDaughLcIndex, arrDaughLcPDGRef, 3); // best would be to check the K0S daughters + if (arrDaughLcIndex.size() == 3) { + for (auto iProng = 0; iProng < arrDaughLcIndex.size(); ++iProng) { + auto daughI = particlesMC.iteratorAt(arrDaughLcIndex[iProng]); + arrDaughLcPDG[iProng] = daughI.pdgCode(); + } + if (!(arrDaughLcPDG[0] == arrDaughLcPDGRef[0] && arrDaughLcPDG[1] == arrDaughLcPDGRef[1] && arrDaughLcPDG[2] == arrDaughLcPDGRef[2])) { // this should be the condition, first bach, then v0 + sign = 0; + } else { + LOG(INFO) << "Lc --> K0S+p found in MC table"; + } + MY_DEBUG_MSG(sign == 0, LOG(INFO) << "Pity, the three final daughters are not p, pi+, pi-, but " << arrDaughLcPDG[0] << ", " << arrDaughLcPDG[1] << ", " << arrDaughLcPDG[2]); + } + } + rowMCMatchGen(sign); + } + } +}; + +//____________________________________________________________________ + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow{ + adaptAnalysisTask(cfgc, TaskName{"hf-cand-creator-cascade"}), + adaptAnalysisTask(cfgc, TaskName{"hf-cand-creator-cascade-expressions"})}; + const bool doMC = cfgc.options().get("doMC"); + if (doMC) { + workflow.push_back(adaptAnalysisTask(cfgc, TaskName{"hf-cand-creator-cascade-mc"})); + } + return workflow; +} diff --git a/Analysis/Tasks/PWGHF/HFLcK0sPCandidateSelector.cxx b/Analysis/Tasks/PWGHF/HFLcK0sPCandidateSelector.cxx new file mode 100644 index 0000000000000..b7b80669065e0 --- /dev/null +++ b/Analysis/Tasks/PWGHF/HFLcK0sPCandidateSelector.cxx @@ -0,0 +1,394 @@ +// 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 HFLcK0spCandidateSelector.cxx +/// \brief Lc --> K0s+p selection task. +/// +/// \author Chiara Zampolli , CERN + +/// based on HFD0CandidateSelector.cxx + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "AnalysisDataModel/HFSecondaryVertex.h" +#include "AnalysisDataModel/HFCandidateSelectionTables.h" +#include "AnalysisTasksUtils/UtilsDebugLcK0Sp.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::aod::hf_cand_casc; + +static const int nPtBins = 8; +static const int nCutVars = 8; +//temporary until 2D array in configurable is solved - then move to json +//mK0s(GeV) mLambdas(GeV) mGammas(GeV) ptp ptK0sdau pTLc d0p d0K0 +constexpr double cuts[nPtBins][nCutVars] = {{0.008, 0.005, 0.1, 0.5, 0.3, 0.6, 0.05, 999999.}, // 1 < pt < 2 + {0.008, 0.005, 0.1, 0.5, 0.4, 1.3, 0.05, 999999.}, // 2 < pt < 3 + {0.009, 0.005, 0.1, 0.6, 0.4, 1.3, 0.05, 999999.}, // 3 < pt < 4 + {0.011, 0.005, 0.1, 0.6, 0.4, 1.4, 0.05, 999999.}, // 4 < pt < 5 + {0.013, 0.005, 0.1, 0.6, 0.4, 1.4, 0.06, 999999.}, // 5 < pt < 6 + {0.013, 0.005, 0.1, 0.9, 0.4, 1.6, 0.09, 999999.}, // 6 < pt < 8 + {0.016, 0.005, 0.1, 0.9, 0.4, 1.7, 0.10, 999999.}, // 8 < pt < 12 + {0.019, 0.005, 0.1, 1.0, 0.4, 1.9, 0.20, 999999.}}; // 12 < pt < 24 +/// Struct for applying D0 selection cuts + +//#define MY_DEBUG + +#ifdef MY_DEBUG +#define MY_DEBUG_MSG(condition, cmd) \ + if (condition) { \ + cmd; \ + } +using MyBigTracks = soa::Join; +#else +#define MY_DEBUG_MSG(condition, cmd) +using MyBigTracks = aod::BigTracksPID; +#endif + +struct HFLcK0sPCandidateSelector { + + Produces hfSelLcK0sPCandidate; + + Configurable pTCandMin{"pTCandMin", 0., "Lower bound of candidate pT"}; + Configurable pTCandMax{"pTCandMax", 50., "Upper bound of candidate pT"}; + + // PID + Configurable applyPidTPCMinPt{"applyPidTPCMinPt", 4., "Lower bound of track pT to apply TPC PID"}; + Configurable pidTPCMinPt{"pidTPCMinPt", 0., "Lower bound of track pT for TPC PID"}; + Configurable pidTPCMaxPt{"pidTPCMaxPt", 100., "Upper bound of track pT for TPC PID"}; + Configurable pidCombMaxP{"pidCombMaxP", 4., "Upper bound of track p to use TOF + TPC Bayes PID"}; + Configurable nSigmaTPC{"nSigmaTPC", 3., "Nsigma cut on TPC only"}; + + // track quality + Configurable TPCNClsFindablePIDCut{"TPCNClsFindablePIDCut", 50., "Lower bound of TPC findable clusters for good PID"}; + Configurable requireTPC{"requireTPC", true, "Flag to require a positive Number of found clusters in TPC"}; + + // for debugging +#ifdef MY_DEBUG + Configurable> indexK0Spos{"indexK0Spos", {729, 2866, 4754, 5457, 6891, 7824, 9243, 9810}, "indices of K0S positive daughters, for debug"}; + Configurable> indexK0Sneg{"indexK0Sneg", {730, 2867, 4755, 5458, 6892, 7825, 9244, 9811}, "indices of K0S negative daughters, for debug"}; + Configurable> indexProton{"indexProton", {717, 2810, 4393, 5442, 6769, 7793, 9002, 9789}, "indices of protons, for debug"}; +#endif + + /// 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) // This should be taken out of the selector, since it is something in common to everyone; + // it should become parameterized with the pt intervals, and also the pt intervals + // should be configurable from outside + { + double ptBins[nPtBins + 1] = {1., 2., 3., 4., 5., 6., 8., 12., 24.}; + 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) // there could be more checks done on the bachelor here + { + // this is for now just a placeholder, in case we want to add extra checks + return true; + } + + /// Conjugate independent toplogical cuts + /// \param hfCandCascade is candidate + /// \return true if candidate passes all cuts + template + bool selectionTopol(const T& hfCandCascade) + { + auto candPt = hfCandCascade.pt(); + int ptBin = getPtBin(candPt); + if (ptBin == -1) { + return false; + } + + if (candPt < pTCandMin || candPt >= pTCandMax) { + LOG(DEBUG) << "cand pt (first check) cut failed: from cascade --> " << candPt << ", cut --> " << pTCandMax; + return false; //check that the candidate pT is within the analysis range + } + + if (std::abs(hfCandCascade.mK0Short() - RecoDecay::getMassPDG(kK0Short)) > cuts[ptBin][0]) { + LOG(DEBUG) << "massK0s cut failed: from v0 in cascade, K0s --> " << hfCandCascade.mK0Short() << ", in PDG K0s --> " << RecoDecay::getMassPDG(kK0Short) << ", cut --> " << cuts[ptBin][0]; + return false; // mass of the K0s + } + + if ((std::abs(hfCandCascade.mLambda() - RecoDecay::getMassPDG(kLambda0)) < cuts[ptBin][1]) || (std::abs(hfCandCascade.mAntiLambda() - RecoDecay::getMassPDG(kLambda0)) < cuts[ptBin][1])) { + LOG(DEBUG) << "mass L cut failed: from v0 in cascade, Lambda --> " << hfCandCascade.mLambda() << ", AntiLambda --> " << hfCandCascade.mAntiLambda() << ", in PDG, Lambda --> " << RecoDecay::getMassPDG(kLambda0) << ", cut --> " << cuts[ptBin][1]; + return false; // mass of the Lambda + } + + if (std::abs(InvMassGamma(hfCandCascade) - RecoDecay::getMassPDG(kGamma)) < cuts[ptBin][2]) { + LOG(DEBUG) << "mass gamma cut failed: from v0 in cascade, gamma --> " << InvMassGamma(hfCandCascade) << ", cut --> " << cuts[ptBin][2]; + return false; // mass of the Gamma + } + + if (hfCandCascade.ptProng0() < cuts[ptBin][3]) { + LOG(DEBUG) << "bach pt cut failed, from cascade --> " << hfCandCascade.ptProng0() << " , cut --> " << cuts[ptBin][3]; + return false; // pt of the p + } + + if (hfCandCascade.ptV0Pos() < cuts[ptBin][4]) { + LOG(DEBUG) << "v0 pos daugh pt cut failed, from cascade --> " << hfCandCascade.ptV0Pos() << ", cut --> " << cuts[ptBin][4]; + return false; // pt of the K0 + } + + if (hfCandCascade.ptV0Neg() < cuts[ptBin][4]) { + LOG(DEBUG) << "v0 neg daugh pt cut failed, from cascade --> " << hfCandCascade.ptV0Neg() << ", cut --> " << cuts[ptBin][4]; + return false; // pt of the K0 + } + + if (hfCandCascade.pt() < cuts[ptBin][5]) { + LOG(DEBUG) << "cand pt cut failed, from cascade --> " << hfCandCascade.pt() << ", cut --> " << cuts[ptBin][5]; + return false; // pt of the Lc + } + + if (std::abs(hfCandCascade.impactParameter0()) > cuts[ptBin][6]) { + LOG(DEBUG) << "d0 bach cut failed, in cascade --> " << hfCandCascade.impactParameter0() << ", cut --> " << cuts[ptBin][6]; + return false; // d0 of the bachelor + } + + /* + if ((std::abs(hfCandCascade.dcapostopv()) > cuts[ptBin][7]) || (std::abs(hfCandCascade.dcanegtopv()) > cuts[ptBin][7])) { + LOG(DEBUG) << "v0 daugh cut failed, positive v0 daugh --> " << hfCandCascade.dcapostopv() << ", negative v0 daugh --> " << hfCandCascade.dcanegtopv() << " , cut --> " << cuts[ptBin][7]; + return false; // d0 of the K0s daughters + } + */ + + if (std::abs(hfCandCascade.impactParameter1()) > cuts[ptBin][7]) { + LOG(DEBUG) << "d0 v0 cut failed, in cascade --> " << hfCandCascade.impactParameter1() << ", cut --> " << cuts[ptBin][7]; + return false; // d0 of the v0 + } + + 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 (track.pt() < pidTPCMinPt || track.pt() >= pidTPCMaxPt) { + LOG(DEBUG) << "Bachelor pt is " << track.pt() << ", we trust TPC PID in [" << pidTPCMinPt << ", " << pidTPCMaxPt << "]"; + return false; + } + return true; + } + + /// Check if we will use TPC PID + /// \param track is the track + /// \note function to be expanded + /// \return true if track is ok for TPC PID + template + bool applyTPCPID(const T& track) + { + if (track.pt() < applyPidTPCMinPt) { + LOG(DEBUG) << "Bachelor pt is " << track.pt() << ", we apply TPC PID from " << applyPidTPCMinPt; + return false; + } + LOG(DEBUG) << "Bachelor pt is " << track.pt() << ", we apply TPC PID from " << applyPidTPCMinPt; + 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 validCombPID(const T& track) + { + if (track.pt() > pidCombMaxP) { // is the pt sign used for the charge? If it is always positive, we should remove the abs + 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=211 pion nPDG=321 kaon + /// \return true if track satisfies TPC PID hypothesis for given Nsigma cut + template + bool selectionPIDTPC(const T& track, double nSigmaCut) + { + double nSigma = 100.0; //arbitarily large value + nSigma = track.tpcNSigmaPr(); + LOG(DEBUG) << "nSigma for bachelor = " << nSigma << ", cut at " << nSigmaCut; + return std::abs(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=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, double nSigmaCut) + { + double nSigma = 100.0; //arbitarily large value + nPDG = TMath::Abs(nPDG); + if (nPDG == 111) { + nSigma = track.tofNSigmaPi(); + } else if (nPDG == 321) { + nSigma = track.tofNSigmaKa(); + } 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=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 statusTPC = -1; + // int statusTOF = -1; + + if (!applyTPCPID(track)) { + // we do not apply TPC PID in this range + return 1; + } + + if (validTPCPID(track)) { + LOG(DEBUG) << "We check the TPC PID now"; + if (!selectionPIDTPC(track, nSigmaTPC)) { + statusTPC = 0; + /* + if (!selectionPIDTPC(track, nPDG, 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; + } + } + + return statusTPC; + /* + if (validTOFPID(track)) { + if (!selectionPIDTOF(track, nPDG, nSigmaTOF)) { + if (!selectionPIDTOF(track, nPDG, 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::HfCandCascade const& candidates, MyBigTracks const& tracks) + { + int statusLc = 0; // final selection flag : 0-rejected 1-accepted + bool topolLc = 0; + int pidProton = -1; + int pidLc = -1; + + for (auto& candidate : candidates) { //looping over cascade candidates + + const auto& bach = candidate.index0_as(); //bachelor track +#ifdef MY_DEBUG + auto indexV0DaughPos = candidate.posTrack_as().mcParticleId(); + auto indexV0DaughNeg = candidate.negTrack_as().mcParticleId(); + auto indexBach = bach.mcParticleId(); + bool isLc = isLcK0SpFunc(indexBach, indexV0DaughPos, indexV0DaughNeg, indexProton, indexK0Spos, indexK0Sneg); + bool isK0SfromLc = isK0SfromLcFunc(indexV0DaughPos, indexV0DaughNeg, indexK0Spos, indexK0Sneg); +#endif + MY_DEBUG_MSG(isLc, printf("\n"); LOG(INFO) << "In selector: correct Lc found: proton --> " << indexBach << ", posTrack --> " << indexV0DaughPos << ", negTrack --> " << indexV0DaughNeg); + //MY_DEBUG_MSG(isLc != 1, printf("\n"); LOG(INFO) << "In selector: wrong Lc found: proton --> " << indexBach << ", posTrack --> " << indexV0DaughPos << ", negTrack --> " << indexV0DaughNeg); + + statusLc = 0; + /* // not needed for the Lc + if (!(candidate.hfflag() & 1 << D0ToPiK)) { + hfSelD0Candidate(statusLc); + continue; + } + */ + + topolLc = true; + pidProton = -1; + + // daughter track validity selection + LOG(DEBUG) << "daughterSelection(bach) = " << daughterSelection(bach); + if (!daughterSelection(bach)) { + MY_DEBUG_MSG(isLc, LOG(INFO) << "In selector: Lc rejected due to selections on bachelor"); + hfSelLcK0sPCandidate(statusLc); + continue; + } + + //implement filter bit 4 cut - should be done before this task at the track selection level + //need to add special cuts (additional cuts on decay length and d0 norm) + LOG(DEBUG) << "selectionTopol(candidate) = " << selectionTopol(candidate); + if (!selectionTopol(candidate)) { + MY_DEBUG_MSG(isLc, LOG(INFO) << "In selector: Lc rejected due to topological selections"); + hfSelLcK0sPCandidate(statusLc); + continue; + } + + pidProton = selectionPID(bach); + + LOG(DEBUG) << "pidProton = " << pidProton; + + if (pidProton == 1) { + statusLc = 1; + } + + MY_DEBUG_MSG(isLc && pidProton != 1, LOG(INFO) << "In selector: Lc rejected due to PID selections on bachelor"); + MY_DEBUG_MSG(isLc && pidProton == 1, LOG(INFO) << "In selector: Lc ACCEPTED"); + + hfSelLcK0sPCandidate(statusLc); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfcg) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfcg, TaskName{"hf-lc-tok0sp-candidate-selector"})}; +} diff --git a/Analysis/Tasks/PWGHF/HFTrackIndexSkimsCreator.cxx b/Analysis/Tasks/PWGHF/HFTrackIndexSkimsCreator.cxx index 4de9237723504..a2684c9e08cf3 100644 --- a/Analysis/Tasks/PWGHF/HFTrackIndexSkimsCreator.cxx +++ b/Analysis/Tasks/PWGHF/HFTrackIndexSkimsCreator.cxx @@ -15,7 +15,6 @@ /// \author Vít Kučera , CERN /// \author Nima Zardoshti , CERN -#include "Framework/runDataProcessing.h" #include "Framework/AnalysisTask.h" #include "Framework/HistogramRegistry.h" #include "DetectorsVertexing/DCAFitterN.h" @@ -23,15 +22,42 @@ #include "AnalysisCore/trackUtilities.h" #include "AnalysisCore/HFConfigurables.h" //#include "AnalysisDataModel/Centrality.h" +#include "AnalysisDataModel/StrangenessTables.h" +#include "AnalysisDataModel/TrackSelectionTables.h" +#include "ReconstructionDataFormats/V0.h" +#include "AnalysisTasksUtils/UtilsDebugLcK0Sp.h" + #include using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::aod; -using namespace o2::analysis; using namespace o2::analysis::hf_cuts_single_track; +void customize(std::vector& workflowOptions) +{ + ConfigParamSpec optionDoMC{"do-LcK0Sp", VariantType::Bool, false, {"Skim also Lc --> K0S+p"}}; + workflowOptions.push_back(optionDoMC); +} + +#include "Framework/runDataProcessing.h" + +//#define MY_DEBUG + +#ifdef MY_DEBUG +using MY_TYPE1 = soa::Join; +using MyTracks = soa::Join; +#define MY_DEBUG_MSG(condition, cmd) \ + if (condition) { \ + cmd; \ + } +#else +using MY_TYPE1 = soa::Join; +using MyTracks = soa::Join; +#define MY_DEBUG_MSG(condition, cmd) +#endif + /// Track selection struct SelectTracks { @@ -43,11 +69,11 @@ struct SelectTracks { Produces rowSelectedTrack; - Configurable b_dovalplots{"b_dovalplots", true, "fill histograms"}; - Configurable d_bz{"d_bz", 5., "bz field"}; + Configurable dovalplots{"dovalplots", true, "fill histograms"}; + Configurable bz{"bz", 5., "bz field"}; // quality cut Configurable doCutQuality{"doCutQuality", true, "apply quality cuts"}; - Configurable d_tpcnclsfound{"d_tpcnclsfound", 70, ">= min. number of TPC clusters needed"}; + Configurable tpcnclsfound{"tpcnclsfound", 70, ">= min. number of TPC clusters needed"}; // pT bins for single-track cuts Configurable> pTBinsTrack{"ptbins_singletrack", std::vector{hf_cuts_single_track::pTBinsTrack_v}, "track pT bin limits for 2-prong DCAXY pT-depentend cut"}; // 2-prong cuts @@ -58,6 +84,19 @@ struct SelectTracks { Configurable ptmintrack_3prong{"ptmintrack_3prong", -1., "min. track pT for 3 prong candidate"}; Configurable> cutsTrack3Prong{"cuts_singletrack_3prong", {hf_cuts_single_track::cutsTrack[0], npTBinsTrack, nCutVarsTrack, pTBinLabelsTrack, cutVarLabelsTrack}, "Single-track selections per pT bin for 3-prong candidates"}; Configurable etamax_3prong{"etamax_3prong", 4., "max. pseudorapidity for 3 prong candidate"}; + // bachelor cuts (when using cascades) + Configurable ptMinTrackBach{"ptMinTrackBach", 0.3, "min. track pT for bachelor in cascade candidate"}; // 0.5 for PbPb 2015? + Configurable dcaToPrimXYMaxPtBach{"dcaToPrimXYMaxPtBach", 2., "max pt cut for min. DCAXY to prim. vtx. for bachelor in cascade candidate"}; + Configurable dcaToPrimXYMinBach{"dcaToPrimXYMinBach", 0., "min. DCAXY to prim. vtx. for bachelor in cascade candidate"}; // for PbPb 2018, the cut should be 0.0025 + Configurable dcaToPrimXYMaxBach{"dcaToPrimXYMaxBach", 1.0, "max. DCAXY to prim. vtx. for bachelor in cascade candidate"}; + Configurable etaMaxBach{"etaMaxBach", 0.8, "max. pseudorapidity for bachelor in cascade candidate"}; + + // for debugging +#ifdef MY_DEBUG + Configurable> indexK0Spos{"indexK0Spos", {729, 2866, 4754, 5457, 6891, 7824, 9243, 9810}, "indices of K0S positive daughters, for debug"}; + Configurable> indexK0Sneg{"indexK0Sneg", {730, 2867, 4755, 5458, 6892, 7825, 9244, 9811}, "indices of K0S negative daughters, for debug"}; + Configurable> indexProton{"indexProton", {717, 2810, 4393, 5442, 6769, 7793, 9002, 9789}, "indices of protons, for debug"}; +#endif HistogramRegistry registry{ "registry", @@ -69,7 +108,13 @@ struct SelectTracks { // 3-prong histograms {"hpt_cuts_3prong", "tracks selected for 3-prong vertexing;#it{p}_{T}^{track} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, {"hdcatoprimxy_cuts_3prong", "tracks selected for 3-prong vertexing;DCAxy to prim. vtx. (cm);entries", {HistType::kTH1F, {{400, -2., 2.}}}}, - {"heta_cuts_3prong", "tracks selected for 3-prong vertexing;#it{#eta};entries", {HistType::kTH1F, {{static_cast(1.2 * etamax_3prong * 100), -1.2 * etamax_3prong, 1.2 * etamax_3prong}}}}}}; + {"heta_cuts_3prong", "tracks selected for 3-prong vertexing;#it{#eta};entries", {HistType::kTH1F, {{static_cast(1.2 * etamax_3prong * 100), -1.2 * etamax_3prong, 1.2 * etamax_3prong}}}}, + // bachelor (for cascades) histograms + {"hpt_cuts_bach", "bachelor tracks selected for cascade vertexing;#it{p}_{T}^{track} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, + {"hdcatoprimxy_cuts_bach", "bachelor tracks selected for cascade vertexing;DCAxy to prim. vtx. (cm);entries", {HistType::kTH1F, {{100, -1., 1.}}}}, + {"heta_cuts_bach", "bachelortracks selected for cascade vertexing;#it{#eta};entries", {HistType::kTH1F, {{100, -1., 1.}}}} + + }}; // array of 2-prong and 3-prong single-track cuts std::array, 2> cutsSingleTrack; @@ -91,25 +136,39 @@ struct SelectTracks { return false; } - if (abs(dca[0]) < cutsSingleTrack[candType].get(pTBinTrack, "min_dcaxytoprimary")) { + if (std::abs(dca[0]) < cutsSingleTrack[candType].get(pTBinTrack, "min_dcaxytoprimary")) { return false; //minimum DCAxy } - if (abs(dca[0]) > cutsSingleTrack[candType].get(pTBinTrack, "max_dcaxytoprimary")) { + if (std::abs(dca[0]) > cutsSingleTrack[candType].get(pTBinTrack, "max_dcaxytoprimary")) { return false; //maximum DCAxy } return true; } void process(aod::Collision const& collision, - soa::Join const& tracks) + MY_TYPE1 const& tracks +#ifdef MY_DEBUG + , + aod::McParticles& mcParticles +#endif + ) { math_utils::Point3D vtxXYZ(collision.posX(), collision.posY(), collision.posZ()); for (auto& track : tracks) { - int status_prong = 3; // selection flag , 2 bits on +#ifdef MY_DEBUG + auto indexBach = track.mcParticleId(); + // LOG(INFO) << "Checking label " << indexBach; + bool isProtonFromLc = isProtonFromLcFunc(indexBach, indexProton); + +#endif + + MY_DEBUG_MSG(isProtonFromLc, LOG(INFO) << "\nWe found the proton " << indexBach); + + int status_prong = 7; // selection flag , 3 bits on auto trackPt = track.pt(); - if (b_dovalplots.value) { + if (dovalplots.value) { registry.get(HIST("hpt_nocuts"))->Fill(trackPt); } @@ -120,31 +179,45 @@ struct SelectTracks { if (trackPt < ptmintrack_3prong) { status_prong = status_prong & ~(1 << 1); } + MY_DEBUG_MSG(isProtonFromLc, LOG(INFO) << "proton " << indexBach << " pt = " << track.pt() << " (cut " << ptMinTrackBach << ")"); + + if (track.pt() < ptMinTrackBach) { + status_prong = status_prong & ~(1 << 2); + } auto trackEta = track.eta(); // eta cut - if ((status_prong & (1 << 0)) && abs(trackEta) > etamax_2prong) { + if ((status_prong & (1 << 0)) && std::abs(trackEta) > etamax_2prong) { status_prong = status_prong & ~(1 << 0); } - if ((status_prong & (1 << 1)) && abs(trackEta) > etamax_3prong) { + if ((status_prong & (1 << 1)) && std::abs(trackEta) > etamax_3prong) { status_prong = status_prong & ~(1 << 1); } + MY_DEBUG_MSG(isProtonFromLc, LOG(INFO) << "proton " << indexBach << " eta = " << trackEta << " (cut " << etaMaxBach << ")"); + + if ((status_prong & (1 << 2)) && std::abs(trackEta) > etaMaxBach) { + status_prong = status_prong & ~(1 << 2); + } // quality cut + MY_DEBUG_MSG(isProtonFromLc, LOG(INFO) << "proton " << indexBach << " tpcNClsFound = " << track.tpcNClsFound() << " (cut " << tpcnclsfound.value << ")"); + if (doCutQuality.value && status_prong > 0) { // FIXME to make a more complete selection e.g track.flags() & o2::aod::track::TPCrefit && track.flags() & o2::aod::track::GoldenChi2 && UChar_t clustermap = track.itsClusterMap(); - if (!(track.tpcNClsFound() >= d_tpcnclsfound.value && + if (!(track.tpcNClsFound() >= tpcnclsfound.value && // is this the number of TPC clusters? It should not be used track.flags() & o2::aod::track::ITSrefit && (TESTBIT(clustermap, 0) || TESTBIT(clustermap, 1)))) { status_prong = 0; + MY_DEBUG_MSG(isProtonFromLc, LOG(INFO) << "proton " << indexBach << " did not pass clusters cut"); } } // DCA cut array dca; if (status_prong > 0) { + double dcaToPrimXYMinBach_ptDep = dcaToPrimXYMinBach * TMath::Max(0., (1 - TMath::Floor(trackPt / dcaToPrimXYMaxPtBach))); auto trackparvar0 = getTrackParCov(track); - if (!trackparvar0.propagateParamToDCA(vtxXYZ, d_bz, &dca, 100.)) { // get impact parameters + if (!trackparvar0.propagateParamToDCA(vtxXYZ, bz, &dca, 100.)) { // get impact parameters status_prong = 0; } if ((status_prong & (1 << 0)) && !isSelectedTrack(track, dca, Cand2Prong)) { @@ -153,10 +226,17 @@ struct SelectTracks { if ((status_prong & (1 << 1)) && !isSelectedTrack(track, dca, Cand3Prong)) { status_prong = status_prong & ~(1 << 1); } + if ((status_prong & (1 << 2)) && (std::abs(dca[0]) < dcaToPrimXYMinBach_ptDep || std::abs(dca[0]) > dcaToPrimXYMaxBach)) { + MY_DEBUG_MSG(isProtonFromLc, + LOG(INFO) << "proton " << indexBach << " did not pass DCA cut"; + LOG(INFO) << "dca[0] = " << dca[0] << " (lower cut " << dcaToPrimXYMinBach_ptDep << ", upper cut " << dcaToPrimXYMaxBach << ")";); + status_prong = status_prong & ~(1 << 2); + } } + MY_DEBUG_MSG(isProtonFromLc, LOG(INFO) << "status_prong = " << status_prong; printf("\n")); // fill histograms - if (b_dovalplots) { + if (dovalplots) { if (status_prong & (1 << 0)) { registry.get(HIST("hpt_cuts_2prong"))->Fill(trackPt); registry.get(HIST("hdcatoprimxy_cuts_2prong"))->Fill(dca[0]); @@ -167,6 +247,12 @@ struct SelectTracks { registry.get(HIST("hdcatoprimxy_cuts_3prong"))->Fill(dca[0]); registry.get(HIST("heta_cuts_3prong"))->Fill(trackEta); } + if (status_prong & (1 << 2)) { + MY_DEBUG_MSG(isProtonFromLc, LOG(INFO) << "Will be kept: Proton from Lc " << indexBach); + registry.get(HIST("hpt_cuts_bach"))->Fill(trackPt); + registry.get(HIST("hdcatoprimxy_cuts_bach"))->Fill(dca[0]); + registry.get(HIST("heta_cuts_bach"))->Fill(trackEta); + } } // fill table row @@ -175,6 +261,8 @@ struct SelectTracks { } }; +//____________________________________________________________________________________________________________________________________________ + /// Pre-selection of 2-prong and 3-prong secondary vertices struct HFTrackIndexSkimsCreator { Produces rowTrackIndexProng2; @@ -183,20 +271,20 @@ struct HFTrackIndexSkimsCreator { Produces rowProng3CutStatus; //Configurable nCollsMax{"nCollsMax", -1, "Max collisions per file"}; //can be added to run over limited collisions per file - for tesing purposes - Configurable b_dovalplots{"b_dovalplots", true, "fill histograms"}; + Configurable dovalplots{"dovalplots", true, "fill histograms"}; Configurable do3prong{"do3prong", 0, "do 3 prong"}; // event selection Configurable triggerindex{"triggerindex", -1, "trigger index"}; // vertexing parameters - Configurable d_bz{"d_bz", 5., "magnetic field kG"}; - Configurable b_propdca{"b_propdca", true, "create tracks version propagated to PCA"}; + Configurable bz{"bz", 5., "magnetic field kG"}; + Configurable propdca{"propdca", true, "create tracks version propagated to PCA"}; Configurable useAbsDCA{"useAbsDCA", true, "Minimise abs. distance rather than chi2"}; - Configurable d_maxr{"d_maxr", 200., "reject PCA's above this radius"}; - Configurable d_maxdzini{"d_maxdzini", 4., "reject (if>0) PCA candidate if tracks DZ exceeds threshold"}; - Configurable d_minparamchange{"d_minparamchange", 1.e-3, "stop iterations if largest change of any X is smaller than this"}; - Configurable d_minrelchi2change{"d_minrelchi2change", 0.9, "stop iterations if chi2/chi2old > this"}; + Configurable maxr{"maxr", 200., "reject PCA's above this radius"}; + Configurable maxdzini{"maxdzini", 4., "reject (if>0) PCA candidate if tracks DZ exceeds threshold"}; + Configurable minparamchange{"minparamchange", 1.e-3, "stop iterations if largest change of any X is smaller than this"}; + Configurable minrelchi2change{"minrelchi2change", 0.9, "stop iterations if chi2/chi2old > this"}; Configurable configs{"configs", {}, "configurables"}; - Configurable b_debug{"b_debug", false, "debug mode"}; + Configurable debug{"debug", false, "debug mode"}; HistogramRegistry registry{ "registry", @@ -220,8 +308,7 @@ struct HFTrackIndexSkimsCreator { {"hmassDsToPiKK", "D_{s} candidates;inv. mass (K K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}}, {"hmassXicToPKPi", "#Xi_{c} candidates;inv. mass (p K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}}}}; - Filter filterSelectTracks = (aod::hf_seltrack::isSelProng > 0); - + Filter filterSelectTracks = (aod::hf_seltrack::isSelProng > 0 && aod::hf_seltrack::isSelProng < 4); using SelectedTracks = soa::Filtered>; // FIXME @@ -353,22 +440,22 @@ struct HFTrackIndexSkimsCreator { // 2-prong vertex fitter o2::vertexing::DCAFitterN<2> df2; - df2.setBz(d_bz); - df2.setPropagateToPCA(b_propdca); - df2.setMaxR(d_maxr); - df2.setMaxDZIni(d_maxdzini); - df2.setMinParamChange(d_minparamchange); - df2.setMinRelChi2Change(d_minrelchi2change); + df2.setBz(bz); + df2.setPropagateToPCA(propdca); + df2.setMaxR(maxr); + df2.setMaxDZIni(maxdzini); + df2.setMinParamChange(minparamchange); + df2.setMinRelChi2Change(minrelchi2change); df2.setUseAbsDCA(useAbsDCA); // 3-prong vertex fitter o2::vertexing::DCAFitterN<3> df3; - df3.setBz(d_bz); - df3.setPropagateToPCA(b_propdca); - df3.setMaxR(d_maxr); - df3.setMaxDZIni(d_maxdzini); - df3.setMinParamChange(d_minparamchange); - df3.setMinRelChi2Change(d_minrelchi2change); + df3.setBz(bz); + df3.setPropagateToPCA(propdca); + df3.setMaxR(maxr); + df3.setMaxDZIni(maxdzini); + df3.setMinParamChange(minparamchange); + df3.setMinRelChi2Change(minrelchi2change); df3.setUseAbsDCA(useAbsDCA); // used to calculate number of candidiates per event @@ -410,7 +497,7 @@ struct HFTrackIndexSkimsCreator { int isSelected2ProngCand = n2ProngBit; //bitmap for checking status of two-prong candidates (1 is true, 0 is rejected) - if (b_debug) { + if (debug) { for (int n2 = 0; n2 < n2ProngDecays; n2++) { for (int n2cut = 0; n2cut < nCuts2Prong; n2cut++) { cutStatus2Prong[n2][n2cut] = true; @@ -427,11 +514,11 @@ struct HFTrackIndexSkimsCreator { for (int n2 = 0; n2 < n2ProngDecays; n2++) { mass2ProngHypo1[n2] = RecoDecay::M(arrMom, arr2Mass1[n2]); mass2ProngHypo2[n2] = RecoDecay::M(arrMom, arr2Mass2[n2]); - if ((b_debug || (isSelected2ProngCand & 1 << n2)) && cut2ProngInvMassCandMin[n2] >= 0. && cut2ProngInvMassCandMax[n2] > 0.) { //no need to check isSelected2Prong but to avoid mistakes + if ((debug || (isSelected2ProngCand & 1 << n2)) && cut2ProngInvMassCandMin[n2] >= 0. && cut2ProngInvMassCandMax[n2] > 0.) { //no need to check isSelected2Prong but to avoid mistakes if ((mass2ProngHypo1[n2] < cut2ProngInvMassCandMin[n2] || mass2ProngHypo1[n2] >= cut2ProngInvMassCandMax[n2]) && (mass2ProngHypo2[n2] < cut2ProngInvMassCandMin[n2] || mass2ProngHypo2[n2] >= cut2ProngInvMassCandMax[n2])) { isSelected2ProngCand = isSelected2ProngCand & ~(1 << n2); - if (b_debug) { + if (debug) { cutStatus2Prong[n2][iDebugCut] = false; } } @@ -452,12 +539,12 @@ struct HFTrackIndexSkimsCreator { auto pVecCandProng2 = RecoDecay::PVec(pvec0, pvec1); // candidate pT cut - if ((b_debug || isSelected2ProngCand > 0) && (std::count_if(std::begin(cut2ProngPtCandMin), std::end(cut2ProngPtCandMin), [](double d) { return d >= 0.; }) > 0)) { + if ((debug || isSelected2ProngCand > 0) && (std::count_if(std::begin(cut2ProngPtCandMin), std::end(cut2ProngPtCandMin), [](double d) { return d >= 0.; }) > 0)) { double cand2ProngPt = RecoDecay::Pt(pVecCandProng2); for (int n2 = 0; n2 < n2ProngDecays; n2++) { - if ((b_debug || (isSelected2ProngCand & 1 << n2)) && cand2ProngPt < cut2ProngPtCandMin[n2]) { + if ((debug || (isSelected2ProngCand & 1 << n2)) && cand2ProngPt < cut2ProngPtCandMin[n2]) { isSelected2ProngCand = isSelected2ProngCand & ~(1 << n2); - if (b_debug) { + if (debug) { cutStatus2Prong[n2][iDebugCut] = false; } } @@ -466,12 +553,12 @@ struct HFTrackIndexSkimsCreator { iDebugCut++; // imp. par. product cut - if ((b_debug || isSelected2ProngCand > 0) && (std::count_if(std::begin(cut2ProngImpParProductCandMax), std::end(cut2ProngImpParProductCandMax), [](double d) { return d < 100.; }) > 0)) { + if ((debug || isSelected2ProngCand > 0) && (std::count_if(std::begin(cut2ProngImpParProductCandMax), std::end(cut2ProngImpParProductCandMax), [](double d) { return d < 100.; }) > 0)) { auto impParProduct = trackPos1.dcaPrim0() * trackNeg1.dcaPrim0(); for (int n2 = 0; n2 < n2ProngDecays; n2++) { - if ((b_debug || (isSelected2ProngCand & 1 << n2)) && impParProduct > cut2ProngImpParProductCandMax[n2]) { + if ((debug || (isSelected2ProngCand & 1 << n2)) && impParProduct > cut2ProngImpParProductCandMax[n2]) { isSelected2ProngCand = isSelected2ProngCand & ~(1 << n2); - if (b_debug) { + if (debug) { cutStatus2Prong[n2][iDebugCut] = false; } } @@ -480,12 +567,12 @@ struct HFTrackIndexSkimsCreator { iDebugCut++; // CPA cut - if ((b_debug || isSelected2ProngCand > 0) && (std::count_if(std::begin(cut2ProngCPACandMin), std::end(cut2ProngCPACandMin), [](double d) { return d > -2.; }) > 0)) { + if ((debug || isSelected2ProngCand > 0) && (std::count_if(std::begin(cut2ProngCPACandMin), std::end(cut2ProngCPACandMin), [](double d) { return d > -2.; }) > 0)) { auto cpa = RecoDecay::CPA(array{collision.posX(), collision.posY(), collision.posZ()}, secondaryVertex2, pVecCandProng2); for (int n2 = 0; n2 < n2ProngDecays; n2++) { - if ((b_debug || (isSelected2ProngCand & 1 << n2)) && cpa < cut2ProngCPACandMin[n2]) { + if ((debug || (isSelected2ProngCand & 1 << n2)) && cpa < cut2ProngCPACandMin[n2]) { isSelected2ProngCand = isSelected2ProngCand & ~(1 << n2); - if (b_debug) { + if (debug) { cutStatus2Prong[n2][iDebugCut] = false; } } @@ -497,7 +584,7 @@ struct HFTrackIndexSkimsCreator { // fill table row rowTrackIndexProng2(trackPos1.globalIndex(), trackNeg1.globalIndex(), isSelected2ProngCand); - if (b_debug) { + if (debug) { int Prong2CutStatus[n2ProngDecays]; for (int n2 = 0; n2 < n2ProngDecays; n2++) { Prong2CutStatus[n2] = nCutStatus2ProngBit; @@ -511,7 +598,7 @@ struct HFTrackIndexSkimsCreator { } // fill histograms - if (b_dovalplots) { + if (dovalplots) { registry.get(HIST("hvtx2_x"))->Fill(secondaryVertex2[0]); registry.get(HIST("hvtx2_y"))->Fill(secondaryVertex2[1]); @@ -559,7 +646,7 @@ struct HFTrackIndexSkimsCreator { int isSelected3ProngCand = n3ProngBit; - if (b_debug) { + if (debug) { for (int n3 = 0; n3 < n3ProngDecays; n3++) { for (int n3cut = 0; n3cut < nCuts3Prong; n3cut++) { cutStatus3Prong[n3][n3cut] = true; @@ -581,13 +668,13 @@ struct HFTrackIndexSkimsCreator { if ((mass3ProngHypo1[n3] < cut3ProngInvMassCandMin[n3] || mass3ProngHypo1[n3] >= cut3ProngInvMassCandMax[n3]) && (mass3ProngHypo2[n3] < cut3ProngInvMassCandMin[n3] || mass3ProngHypo2[n3] >= cut3ProngInvMassCandMax[n3])) { isSelected3ProngCand = isSelected3ProngCand & ~(1 << n3); - if (b_debug) { + if (debug) { cutStatus3Prong[n3][iDebugCut] = false; } } } } - if (!b_debug && isSelected3ProngCand == 0) { + if (!debug && isSelected3ProngCand == 0) { continue; } iDebugCut++; @@ -616,11 +703,11 @@ struct HFTrackIndexSkimsCreator { if (cand3ProngPt < cut3ProngPtCandMin[n3]) { isSelected3ProngCand = isSelected3ProngCand & ~(1 << n3); } - if (b_debug) { + if (debug) { cutStatus3Prong[n3][iDebugCut] = false; } } - if (!b_debug && isSelected3ProngCand == 0) { + if (!debug && isSelected3ProngCand == 0) { continue; //this and all further instances should be changed if 4 track loop is added } } @@ -633,11 +720,11 @@ struct HFTrackIndexSkimsCreator { if ((isSelected3ProngCand & 1 << n3) && cpa < cut3ProngCPACandMin[n3]) { isSelected3ProngCand = isSelected3ProngCand & ~(1 << n3); } - if (b_debug) { + if (debug) { cutStatus3Prong[n3][iDebugCut] = false; } } - if (!b_debug && isSelected3ProngCand == 0) { + if (!debug && isSelected3ProngCand == 0) { continue; } } @@ -649,12 +736,12 @@ struct HFTrackIndexSkimsCreator { for (int n3 = 0; n3 < n3ProngDecays; n3++) { if ((isSelected3ProngCand & 1 << n3) && decayLength < cut3ProngDecLenCandMin[n3]) { isSelected3ProngCand = isSelected3ProngCand & ~(1 << n3); - if (b_debug) { + if (debug) { cutStatus3Prong[n3][iDebugCut] = false; } } } - if (!b_debug && isSelected3ProngCand == 0) { + if (!debug && isSelected3ProngCand == 0) { continue; } } @@ -665,7 +752,7 @@ struct HFTrackIndexSkimsCreator { trackNeg1.globalIndex(), trackPos2.globalIndex(), isSelected3ProngCand); - if (b_debug) { + if (debug) { int Prong3CutStatus[n3ProngDecays]; for (int n3 = 0; n3 < n3ProngDecays; n3++) { Prong3CutStatus[n3] = nCutStatus3ProngBit; @@ -679,7 +766,7 @@ struct HFTrackIndexSkimsCreator { } // fill histograms - if (b_dovalplots) { + if (dovalplots) { registry.get(HIST("hvtx3_x"))->Fill(secondaryVertex3[0]); registry.get(HIST("hvtx3_y"))->Fill(secondaryVertex3[1]); @@ -731,7 +818,7 @@ struct HFTrackIndexSkimsCreator { int isSelected3ProngCand = n3ProngBit; - if (b_debug) { + if (debug) { for (int n3 = 0; n3 < n3ProngDecays; n3++) { for (int n3cut = 0; n3cut < nCuts3Prong; n3cut++) { cutStatus3Prong[n3][n3cut] = true; @@ -753,13 +840,13 @@ struct HFTrackIndexSkimsCreator { if ((mass3ProngHypo1[n3] < cut3ProngInvMassCandMin[n3] || mass3ProngHypo1[n3] >= cut3ProngInvMassCandMax[n3]) && (mass3ProngHypo2[n3] < cut3ProngInvMassCandMin[n3] || mass3ProngHypo2[n3] >= cut3ProngInvMassCandMax[n3])) { isSelected3ProngCand = isSelected3ProngCand & ~(1 << n3); - if (b_debug) { + if (debug) { cutStatus3Prong[n3][iDebugCut] = false; } } } } - if (!b_debug && isSelected3ProngCand == 0) { + if (!debug && isSelected3ProngCand == 0) { continue; } iDebugCut++; @@ -788,12 +875,12 @@ struct HFTrackIndexSkimsCreator { for (int n3 = 0; n3 < n3ProngDecays; n3++) { if (cand3ProngPt < cut3ProngPtCandMin[n3]) { isSelected3ProngCand = isSelected3ProngCand & ~(1 << n3); - if (b_debug) { + if (debug) { cutStatus3Prong[n3][iDebugCut] = false; } } } - if (!b_debug && isSelected3ProngCand == 0) { + if (!debug && isSelected3ProngCand == 0) { continue; } } @@ -805,12 +892,12 @@ struct HFTrackIndexSkimsCreator { for (int n3 = 0; n3 < n3ProngDecays; n3++) { if ((isSelected3ProngCand & 1 << n3) && cpa < cut3ProngCPACandMin[n3]) { isSelected3ProngCand = isSelected3ProngCand & ~(1 << n3); - if (b_debug) { + if (debug) { cutStatus3Prong[n3][iDebugCut] = false; } } } - if (!b_debug && isSelected3ProngCand == 0) { + if (!debug && isSelected3ProngCand == 0) { continue; } } @@ -822,12 +909,12 @@ struct HFTrackIndexSkimsCreator { for (int n3 = 0; n3 < n3ProngDecays; n3++) { if ((isSelected3ProngCand & 1 << n3) && decayLength < cut3ProngDecLenCandMin[n3]) { isSelected3ProngCand = isSelected3ProngCand & ~(1 << n3); - if (b_debug) { + if (debug) { cutStatus3Prong[n3][iDebugCut] = false; } } } - if (!b_debug && isSelected3ProngCand == 0) { + if (!debug && isSelected3ProngCand == 0) { continue; } } @@ -838,7 +925,7 @@ struct HFTrackIndexSkimsCreator { trackPos1.globalIndex(), trackNeg2.globalIndex(), isSelected3ProngCand); - if (b_debug) { + if (debug) { int Prong3CutStatus[n3ProngDecays]; for (int n3 = 0; n3 < n3ProngDecays; n3++) { Prong3CutStatus[n3] = nCutStatus3ProngBit; @@ -852,7 +939,7 @@ struct HFTrackIndexSkimsCreator { } // fill histograms - if (b_dovalplots) { + if (dovalplots) { registry.get(HIST("hvtx3_x"))->Fill(secondaryVertex3[0]); registry.get(HIST("hvtx3_y"))->Fill(secondaryVertex3[1]); @@ -907,9 +994,284 @@ struct HFTrackIndexSkimsCreator { } }; +//________________________________________________________________________________________________________________________ + +/// Pre-selection of cascade secondary vertices +/// It will produce in any case a HfTrackIndexProng2 object, but mixing a V0 +/// with a track, instead of 2 tracks + +/// to run: o2-analysis-weak-decay-indices --aod-file AO2D.root -b | o2-analysis-lambdakzerobuilder -b | +/// o2-analysis-trackextension -b | o2-analysis-hf-track-index-skims-creator -b + +struct HFTrackIndexSkimsCreatorCascades { + Produces rowTrackIndexCasc; + // Produces rowTrackIndexCasc; + + // whether to do or not validation plots + Configurable doValPlots{"doValPlots", true, "fill histograms"}; + + // event selection + //Configurable triggerindex{"triggerindex", -1, "trigger index"}; + + // vertexing parameters + Configurable bZ{"bZ", 5., "magnetic field"}; + Configurable propDCA{"propDCA", true, "create tracks version propagated to PCA"}; + Configurable maxR{"maxR", 200., "reject PCA's above this radius"}; + Configurable maxDZIni{"maxDZIni", 4., "reject (if>0) PCA candidate if tracks DZ exceeds threshold"}; + Configurable minParamChange{"minParamChange", 1.e-3, "stop iterations if largest change of any X is smaller than this"}; + Configurable minRelChi2Change{"minRelChi2Change", 0.9, "stop iterations if chi2/chi2old > this"}; + Configurable UseAbsDCA{"UseAbsDCA", true, "Use Abs DCAs"}; + + // quality cut + Configurable doCutQuality{"doCutQuality", true, "apply quality cuts"}; + + // track cuts for V0 daughters + Configurable TPCRefit{"TPCRefit", true, "request TPC refit V0 daughters"}; + Configurable i_minCrossedRows{"i_minCrossedRows", 50, "min crossed rows V0 daughters"}; + Configurable etaMax{"etaMax", 1.1, "max. pseudorapidity V0 daughters"}; + Configurable ptMin{"ptMin", 0.05, "min. pT V0 daughters"}; + + // bachelor cuts + // Configurable dcabachtopv{"dcabachtopv", .1, "DCA Bach To PV"}; + // Configurable ptminbach{"ptminbach", -1., "min. track pT bachelor"}; + + // v0 cuts + Configurable cosPAV0{"cosPAV0", .995, "CosPA V0"}; // as in the task that create the V0s + Configurable dcaXYNegToPV{"dcaXYNegToPV", .1, "DCA_XY Neg To PV"}; // check: in HF Run 2, it was 0 at filtering + Configurable dcaXYPosToPV{"dcaXYPosToPVS", .1, "DCA_XY Pos To PV"}; // check: in HF Run 2, it was 0 at filtering + Configurable cutInvMassV0{"cutInvMassV0", 0.05, "V0 candidate invariant mass difference wrt PDG"}; + + // cascade cuts + Configurable cutCascPtCandMin{"cutCascPtCandMin", -1., "min. pT of the 2-prong candidate"}; // PbPb 2018: use 1 + Configurable cutCascInvMassLc{"cutCascInvMassLc", 1., "Lc candidate invariant mass difference wrt PDG"}; // for PbPb 2018: use 0.2 + //Configurable cutCascDCADaughters{"cutCascDCADaughters", .1, "DCA between V0 and bachelor in cascade"}; + + // for debugging +#ifdef MY_DEBUG + Configurable> indexK0Spos{"indexK0Spos", {729, 2866, 4754, 5457, 6891, 7824, 9243, 9810}, "indices of K0S positive daughters, for debug"}; + Configurable> indexK0Sneg{"indexK0Sneg", {730, 2867, 4755, 5458, 6892, 7825, 9244, 9811}, "indices of K0S negative daughters, for debug"}; + Configurable> indexProton{"indexProton", {717, 2810, 4393, 5442, 6769, 7793, 9002, 9789}, "indices of protons, for debug"}; +#endif + + // histograms + HistogramRegistry registry{ + "registry", + {{"hvtx2_x", "2-prong candidates;#it{x}_{sec. vtx.} (cm);entries", {HistType::kTH1F, {{1000, -2., 2.}}}}, + {"hvtx2_y", "2-prong candidates;#it{y}_{sec. vtx.} (cm);entries", {HistType::kTH1F, {{1000, -2., 2.}}}}, + {"hvtx2_z", "2-prong candidates;#it{z}_{sec. vtx.} (cm);entries", {HistType::kTH1F, {{1000, -20., 20.}}}}, + {"hmass2", "2-prong candidates;inv. mass (K0s p) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0., 5.}}}}}}; + + // NB: using FullTracks = soa::Join; defined in Framework/Core/include/Framework/AnalysisDataModel.h + //using MyTracks = aod::BigTracksMC; + //Partition selectedTracks = aod::hf_seltrack::isSelProng >= 4; + // using SelectedV0s = soa::Filtered; + + double massP = RecoDecay::getMassPDG(kProton); + double massK0s = RecoDecay::getMassPDG(kK0Short); + double massPi = RecoDecay::getMassPDG(kPiPlus); + double massLc = RecoDecay::getMassPDG(pdg::Code::kLambdaCPlus); + double mass2K0sP{0.}; // WHY HERE? + + using FullTracksExt = soa::Join; + + void process(aod::Collision const& collision, + aod::BCs const& bcs, + //soa::Filtered const& V0s, + aod::V0Datas const& V0s, + MyTracks const& tracks +#ifdef MY_DEBUG + , + aod::McParticles& mcParticles +#endif + ) // TODO: I am now assuming that the V0s are already filtered with my cuts (David's work to come) + { + + //Define o2 fitter, 2-prong + o2::vertexing::DCAFitterN<2> fitter; + fitter.setBz(bZ); + fitter.setPropagateToPCA(propDCA); + fitter.setMaxR(maxR); + fitter.setMinParamChange(minParamChange); + fitter.setMinRelChi2Change(minRelChi2Change); + //fitter.setMaxDZIni(1e9); // used in cascadeproducer.cxx, but not for the 2 prongs + //fitter.setMaxChi2(1e9); // used in cascadeproducer.cxx, but not for the 2 prongs + fitter.setUseAbsDCA(UseAbsDCA); + + // fist we loop over the bachelor candidate + + //for (const auto& bach : selectedTracks) { + for (const auto& bach : tracks) { + + MY_DEBUG_MSG(1, printf("\n"); LOG(INFO) << "Bachelor loop"); +#ifdef MY_DEBUG + auto indexBach = bach.mcParticleId(); + bool isProtonFromLc = isProtonFromLcFunc(indexBach, indexProton); +#endif + // selections on the bachelor + // pT cut + if (bach.isSelProng() < 4) { + MY_DEBUG_MSG(isProtonFromLc, LOG(INFO) << "proton " << indexBach << ": rejected due to HFsel"); + continue; + } + + if (TPCRefit) { + if (!(bach.trackType() & o2::aod::track::TPCrefit)) { + MY_DEBUG_MSG(isProtonFromLc, LOG(INFO) << "proton " << indexBach << ": rejected due to TPCrefit"); + continue; + } + } + if (bach.tpcNClsCrossedRows() < i_minCrossedRows) { + MY_DEBUG_MSG(isProtonFromLc, LOG(INFO) << "proton " << indexBach << ": rejected due to minNUmberOfCrossedRows " << bach.tpcNClsCrossedRows() << " (cut " << i_minCrossedRows << ")"); + continue; + } + MY_DEBUG_MSG(isProtonFromLc, LOG(INFO) << "KEPT! proton from Lc with daughters " << indexBach); + + auto trackBach = getTrackParCov(bach); + // now we loop over the V0s + for (const auto& v0 : V0s) { + MY_DEBUG_MSG(1, LOG(INFO) << "*** Checking next K0S"); + // selections on the V0 daughters + const auto& trackV0DaughPos = v0.posTrack_as(); + const auto& trackV0DaughNeg = v0.negTrack_as(); +#ifdef MY_DEBUG + auto indexV0DaughPos = trackV0DaughPos.mcParticleId(); + auto indexV0DaughNeg = trackV0DaughNeg.mcParticleId(); + bool isK0SfromLc = isK0SfromLcFunc(indexV0DaughPos, indexV0DaughNeg, indexK0Spos, indexK0Sneg); + + bool isLc = isLcK0SpFunc(indexBach, indexV0DaughPos, indexV0DaughNeg, indexProton, indexK0Spos, indexK0Sneg); +#endif + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "K0S from Lc found, trackV0DaughPos --> " << indexV0DaughPos << ", trackV0DaughNeg --> " << indexV0DaughNeg); + + MY_DEBUG_MSG(isK0SfromLc && isProtonFromLc, + LOG(INFO) << "ACCEPTED!!!"; + LOG(INFO) << "proton belonging to a Lc found: label --> " << indexBach; + LOG(INFO) << "K0S belonging to a Lc found: trackV0DaughPos --> " << indexV0DaughPos << ", trackV0DaughNeg --> " << indexV0DaughNeg); + + MY_DEBUG_MSG(isLc, LOG(INFO) << "Combination of K0S and p which correspond to a Lc found!"); + + if (TPCRefit) { + if (!(trackV0DaughPos.trackType() & o2::aod::track::TPCrefit) || + !(trackV0DaughNeg.trackType() & o2::aod::track::TPCrefit)) { + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "K0S with daughters " << indexV0DaughPos << " and " << indexV0DaughNeg << ": rejected due to TPCrefit"); + continue; + } + } + if (trackV0DaughPos.tpcNClsCrossedRows() < i_minCrossedRows || + trackV0DaughNeg.tpcNClsCrossedRows() < i_minCrossedRows) { + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "K0S with daughters " << indexV0DaughPos << " and " << indexV0DaughNeg << ": rejected due to minCrossedRows"); + continue; + } + // + // if (trackV0DaughPos.dcaXY() < dcaXYPosToPV || // to the filters? + // trackV0DaughNeg.dcaXY() < dcaXYNegToPV) { + // continue; + // } + // + if (trackV0DaughPos.pt() < ptMin || // to the filters? I can't for now, it is not in the tables + trackV0DaughNeg.pt() < ptMin) { + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "K0S with daughters " << indexV0DaughPos << " and " << indexV0DaughNeg << ": rejected due to minPt --> pos " << trackV0DaughPos.pt() << ", neg " << trackV0DaughNeg.pt() << " (cut " << ptMin << ")"); + continue; + } + if (std::abs(trackV0DaughPos.eta()) > etaMax || // to the filters? I can't for now, it is not in the tables + std::abs(trackV0DaughNeg.eta()) > etaMax) { + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "K0S with daughters " << indexV0DaughPos << " and " << indexV0DaughNeg << ": rejected due to eta --> pos " << trackV0DaughPos.eta() << ", neg " << trackV0DaughNeg.eta() << " (cut " << etaMax << ")"); + continue; + } + + // V0 invariant mass selection + if (std::abs(v0.mK0Short() - massK0s) > cutInvMassV0) { + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "K0S with daughters " << indexV0DaughPos << " and " << indexV0DaughNeg << ": rejected due to invMass --> " << v0.mK0Short() - massK0s << " (cut " << cutInvMassV0 << ")"); + continue; // should go to the filter, but since it is a dynamic column, I cannot use it there + } + + // V0 cosPointingAngle selection + if (v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()) < cosPAV0) { + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "K0S with daughters " << indexV0DaughPos << " and " << indexV0DaughNeg << ": rejected due to cosPA --> " << v0.v0cosPA(collision.posX(), collision.posY(), collision.posZ()) << " (cut " << cosPAV0 << ")"); + continue; + } + + const std::array momentumV0 = {v0.px(), v0.py(), v0.pz()}; + + // invariant-mass cut: we do it here, before updating the momenta of bach and V0 during the fitting to save CPU + // TODO: but one should better check that the value here and after the fitter do not change significantly!!! + mass2K0sP = RecoDecay::M(array{array{bach.px(), bach.py(), bach.pz()}, momentumV0}, array{massP, massK0s}); + if ((cutCascInvMassLc >= 0.) && (std::abs(mass2K0sP - massLc) > cutCascInvMassLc)) { + MY_DEBUG_MSG(isK0SfromLc && isProtonFromLc, LOG(INFO) << "True Lc from proton " << indexBach << " and K0S pos " << indexV0DaughPos << " and neg " << indexV0DaughNeg << " rejected due to invMass cut: " << mass2K0sP << ", mass Lc " << massLc << " (cut " << cutCascInvMassLc << ")"); + continue; + } + + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "KEPT! K0S from Lc with daughters " << indexV0DaughPos << " and " << indexV0DaughNeg); + + auto trackParCovV0DaughPos = getTrackParCov(trackV0DaughPos); + trackParCovV0DaughPos.propagateTo(v0.posX(), bZ); // propagate the track to the X closest to the V0 vertex + auto trackParCovV0DaughNeg = getTrackParCov(trackV0DaughNeg); + trackParCovV0DaughNeg.propagateTo(v0.negX(), bZ); // propagate the track to the X closest to the V0 vertex + std::array pVecV0 = {0., 0., 0.}; + std::array pVecBach = {0., 0., 0.}; + + const std::array vertexV0 = {v0.x(), v0.y(), v0.z()}; + // we build the neutral track to then build the cascade + auto trackV0 = o2::dataformats::V0(vertexV0, momentumV0, {0, 0, 0, 0, 0, 0}, trackParCovV0DaughPos, trackParCovV0DaughNeg, {0, 0}, {0, 0}); // build the V0 track + + // now we find the DCA between the V0 and the bachelor, for the cascade + int nCand2 = fitter.process(trackV0, trackBach); + MY_DEBUG_MSG(isK0SfromLc && isProtonFromLc, LOG(INFO) << "Fitter result = " << nCand2 << " proton = " << indexBach << " and K0S pos " << indexV0DaughPos << " and neg " << indexV0DaughNeg); + MY_DEBUG_MSG(isLc, LOG(INFO) << "Fitter result for true Lc = " << nCand2); + if (nCand2 == 0) { + continue; + } + fitter.propagateTracksToVertex(); // propagate the bach and V0 to the Lc vertex + fitter.getTrack(0).getPxPyPzGlo(pVecV0); // take the momentum at the Lc vertex + fitter.getTrack(1).getPxPyPzGlo(pVecBach); + + // cascade candidate pT cut + auto ptCascCand = RecoDecay::Pt(pVecBach, pVecV0); + if (ptCascCand < cutCascPtCandMin) { + MY_DEBUG_MSG(isK0SfromLc && isProtonFromLc, LOG(INFO) << "True Lc from proton " << indexBach << " and K0S pos " << indexV0DaughPos << " and neg " << indexV0DaughNeg << " rejected due to pt cut: " << ptCascCand << " (cut " << cutCascPtCandMin << ")"); + continue; + } + + // invariant mass + // re-calculate invariant masses with updated momenta, to fill the histogram + mass2K0sP = RecoDecay::M(array{pVecBach, pVecV0}, array{massP, massK0s}); + + std::array posCasc = {0., 0., 0.}; + const auto& cascVtx = fitter.getPCACandidate(); + for (int i = 0; i < 3; i++) { + posCasc[i] = cascVtx[i]; + } + + // fill table row + rowTrackIndexCasc(bach.globalIndex(), + v0.globalIndex(), + 1); // 1 should be the value for the Lc + // fill histograms + if (doValPlots) { + MY_DEBUG_MSG(isK0SfromLc && isProtonFromLc && isLc, LOG(INFO) << "KEPT! True Lc from proton " << indexBach << " and K0S pos " << indexV0DaughPos << " and neg " << indexV0DaughNeg); + registry.get(HIST("hvtx2_x"))->Fill(posCasc[0]); + registry.get(HIST("hvtx2_y"))->Fill(posCasc[1]); + registry.get(HIST("hvtx2_z"))->Fill(posCasc[2]); + registry.get(HIST("hmass2"))->Fill(mass2K0sP); + } + + } // loop over V0s + + } // loop over tracks + } // process +}; + +//________________________________________________________________________________________________________________________ + WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - return WorkflowSpec{ + WorkflowSpec workflow{ adaptAnalysisTask(cfgc, TaskName{"hf-produce-sel-track"}), adaptAnalysisTask(cfgc, TaskName{"hf-track-index-skims-creator"})}; + + const bool doLcK0Sp = cfgc.options().get("do-LcK0Sp"); + if (doLcK0Sp) { + workflow.push_back(adaptAnalysisTask(cfgc, TaskName{"hf-track-index-skims-cascades-creator"})); + } + + return workflow; } diff --git a/Analysis/Tasks/PWGHF/taskLcK0sP.cxx b/Analysis/Tasks/PWGHF/taskLcK0sP.cxx new file mode 100644 index 0000000000000..237461d281c4e --- /dev/null +++ b/Analysis/Tasks/PWGHF/taskLcK0sP.cxx @@ -0,0 +1,153 @@ +// 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 taskLcK0sp.cxx +/// \brief Lc -> K0S+p analysis task +/// +/// \author Chiara Zampolli, , CERN +/// +/// based on taskD0.cxx, taskLc.cxx + +#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_casc; +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" + +/// LcK0sp analysis task +struct TaskLcK0sP { + HistogramRegistry registry{ + "registry", + {{"hMass", "cascade candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 0.0f, 5.0f}}}}, + {"hPtCand", "cascade candidates;candidate #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0.0f, 10.0f}}}}, + {"hPtBach", "cascade candidates;bachelor #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0.0f, 10.0f}}}}, + {"hPtV0", "cascade candidates;v0 #it{p}_{T} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0.0f, 10.0f}}}}, + {"hd0Bach", "cascade candidates;bachelor DCAxy to prim. vertex (cm);entries", {HistType::kTH1F, {{100, -1.0f, 1.0f}}}}, + {"hd0V0pos", "cascade candidates;pos daugh v0 DCAxy to prim. vertex (cm);entries", {HistType::kTH1F, {{100, -5.0f, 5.0f}}}}, + {"hd0V0neg", "cascade candidates;neg daugh v0 DCAxy to prim. vertex (cm);entries", {HistType::kTH1F, {{100, -5.0f, 5.0f}}}}, + {"hV0CPA", "cascade candidates;v0 cosine of pointing angle;entries", {HistType::kTH1F, {{110, -0.98f, 1.1f}}}}, + {"hEta", "cascade candidates;candidate #it{#eta};entries", {HistType::kTH1F, {{100, -2.0f, 2.0f}}}}, + {"hSelectionStatus", "cascade candidates;selection status;entries", {HistType::kTH1F, {{5, -0.5f, 4.5f}}}}}}; + + Configurable selectionFlagLcK0sp{"selectionFlagLcK0sp", 1, "Selection Flag for LcK0sp"}; + Configurable cutEtaCandMax{"cutEtaCandMax", -1., "max. cand. pseudorapidity"}; + + Filter filterSelectCandidates = (aod::hf_selcandidate_lc_k0sp::isSelLcK0sP >= selectionFlagLcK0sp); + + void process(soa::Filtered> const& candidates) + { + //Printf("Candidates: %d", candidates.size()); + for (auto& candidate : candidates) { + /* + // no such selection for LcK0sp for now - it is the only cascade + if (!(candidate.hfflag() & 1 << D0ToPiK)) { + continue; + } + */ + if (cutEtaCandMax >= 0. && std::abs(candidate.eta()) > cutEtaCandMax) { + //Printf("Candidate: eta rejection: %g", candidate.eta()); + continue; + } + + registry.fill(HIST("hMass"), InvMassLcToK0sP(candidate)); + registry.fill(HIST("hPtCand"), candidate.pt()); + registry.fill(HIST("hPtBach"), candidate.ptProng0()); + registry.fill(HIST("hPtV0"), candidate.ptProng1()); + registry.fill(HIST("hd0Bach"), candidate.impactParameter0()); + registry.fill(HIST("hd0V0pos"), candidate.dcapostopv()); + registry.fill(HIST("hd0V0neg"), candidate.dcanegtopv()); + registry.fill(HIST("hV0CPA"), candidate.v0cosPA()); + registry.fill(HIST("hEta"), candidate.eta()); + registry.fill(HIST("hSelectionStatus"), candidate.isSelLcK0sP()); + } + } +}; + +/// Fills MC histograms. +struct TaskLcK0SpMC { + HistogramRegistry registry{ + "registry", + {{"hPtRecSig", "cascade candidates (MC);#it{p}_{T}^{rec.} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, + {"hPtRecBg", "cascade candidates (unmatched);#it{p}_{T}^{rec.} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, + {"hPtGen", "cascade (MC);#it{p}_{T}^{gen.} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, + {"hPtGenSig", "cascade candidates (MC);#it{p}_{T}^{gen.} (GeV/#it{c});entries", {HistType::kTH1F, {{100, 0., 10.}}}}, + {"hCPARecSig", "cascade candidates (matched);cosine of pointing angle;entries", {HistType::kTH1F, {{110, -1.1, 1.1}}}}, + {"hCPARecBg", "cascade candidates (unmatched);cosine of pointing angle;entries", {HistType::kTH1F, {{110, -1.1, 1.1}}}}, + {"hEtaRecSig", "cascade candidates (matched);#it{#eta};entries", {HistType::kTH1F, {{100, -2., 2.}}}}, + {"hEtaRecBg", "cascade candidates (unmatched);#it{#eta};entries", {HistType::kTH1F, {{100, -2., 2.}}}}, + {"hEtaGen", "MC particles (MC);#it{#eta};entries", {HistType::kTH1F, {{100, -2., 2.}}}}}}; + + Configurable selectionFlagLc{"selectionFlagLc", 1, "Selection Flag for Lc"}; + Configurable selectionFlagLcbar{"selectionFlagLcbar", 1, "Selection Flag for Lcbar"}; + Configurable cutEtaCandMax{"cutEtaCandMax", -1., "max. cand. pseudorapidity"}; + + Filter filterSelectCandidates = (aod::hf_selcandidate_lc_k0sp::isSelLcK0sP >= selectionFlagLc || aod::hf_selcandidate_lc_k0sp::isSelLcK0sP >= selectionFlagLcbar); + + void process(soa::Filtered> const& candidates, + soa::Join const& particlesMC, aod::BigTracksMC const& tracks) + { + // MC rec. + //Printf("MC Candidates: %d", candidates.size()); + for (auto& candidate : candidates) { + if (cutEtaCandMax >= 0. && std::abs(candidate.eta()) > cutEtaCandMax) { + //Printf("MC Rec.: eta rejection: %g", candidate.eta()); + continue; + } + if (std::abs(candidate.flagMCMatchRec()) == 1) { + // Get the corresponding MC particle. + auto indexMother = RecoDecay::getMother(particlesMC, candidate.index0_as().mcParticle_as>(), pdg::Code::kLambdaCPlus, true); + auto particleMother = particlesMC.iteratorAt(indexMother); + registry.fill(HIST("hPtGenSig"), particleMother.pt()); // gen. level pT + registry.fill(HIST("hPtRecSig"), candidate.pt()); // rec. level 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()) == 1) { + registry.fill(HIST("hPtGen"), particle.pt()); + registry.fill(HIST("hEtaGen"), particle.eta()); + } + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow{ + adaptAnalysisTask(cfgc, TaskName{"hf-task-lc-tok0sP"})}; + const bool doMC = cfgc.options().get("doMC"); + if (doMC) { + workflow.push_back(adaptAnalysisTask(cfgc, TaskName{"hf-task-lc-tok0sP-mc"})); + } + return workflow; +} diff --git a/Analysis/Tasks/PWGLF/CMakeLists.txt b/Analysis/Tasks/PWGLF/CMakeLists.txt index 3103556034c03..5e865fba760e9 100644 --- a/Analysis/Tasks/PWGLF/CMakeLists.txt +++ b/Analysis/Tasks/PWGLF/CMakeLists.txt @@ -60,7 +60,7 @@ o2_add_dpl_workflow(nuclei-spectra o2_add_dpl_workflow(lambdakzerobuilder SOURCES lambdakzerobuilder.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2::AnalysisDataModel O2::AnalysisCore O2::DetectorsVertexing O2::AnalysisTasksUtils COMPONENT_NAME Analysis) o2_add_dpl_workflow(lambdakzeroanalysis diff --git a/Analysis/Tasks/PWGLF/lambdakzerobuilder.cxx b/Analysis/Tasks/PWGLF/lambdakzerobuilder.cxx index d830b40b9c3c9..c71af8da32fda 100644 --- a/Analysis/Tasks/PWGLF/lambdakzerobuilder.cxx +++ b/Analysis/Tasks/PWGLF/lambdakzerobuilder.cxx @@ -53,6 +53,7 @@ #include #include #include "Framework/ASoAHelpers.h" +#include "AnalysisTasksUtils/UtilsDebugLcK0Sp.h" using namespace o2; using namespace o2::framework; @@ -73,6 +74,19 @@ DECLARE_SOA_TABLE(V0GoodIndices, "AOD", "V0GOODINDICES", o2::soa::Index<>, } // namespace o2::aod using FullTracksExt = soa::Join; +using FullTracksExtMC = soa::Join; + +//#define MY_DEBUG +#ifdef MY_DEBUG +using MyTracks = FullTracksExtMC; +#define MY_DEBUG_MSG(condition, cmd) \ + if (condition) { \ + cmd; \ + } +#else +using MyTracks = FullTracksExt; +#define MY_DEBUG_MSG(condition, cmd) +#endif //This prefilter creates a skimmed list of good V0s to re-reconstruct so that //CPU is saved in case there are specific selections that are to be done @@ -83,6 +97,12 @@ struct lambdakzeroprefilterpairs { Configurable mincrossedrows{"mincrossedrows", 70, "min crossed rows"}; Configurable tpcrefit{"tpcrefit", 1, "demand TPC refit"}; + // for debugging +#ifdef MY_DEBUG + Configurable> v_labelK0Spos{"v_labelK0Spos", {729, 2866, 4754}, "labels of K0S positive daughters, for debug"}; + Configurable> v_labelK0Sneg{"v_labelK0Sneg", {730, 2867, 4755}, "labels of K0S negative daughters, for debug"}; +#endif + HistogramRegistry registry{ "registry", { @@ -93,36 +113,56 @@ struct lambdakzeroprefilterpairs { Produces v0goodindices; void process(aod::Collision const& collision, aod::V0s const& V0s, - soa::Join const& tracks) + MyTracks const& tracks +#ifdef MY_DEBUG + , + aod::McParticles const& particlesMC +#endif + ) { for (auto& V0 : V0s) { + +#ifdef MY_DEBUG + auto labelPos = V0.posTrack_as().mcParticleId(); + auto labelNeg = V0.negTrack_as().mcParticleId(); + bool isK0SfromLc = isK0SfromLcFunc(labelPos, labelNeg, v_labelK0Spos, v_labelK0Sneg); +#endif + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "V0 builder: found K0S from Lc, posTrack --> " << labelPos << ", negTrack --> " << labelNeg); + registry.fill(HIST("hGoodIndices"), 0.5); if (tpcrefit) { - if (!(V0.posTrack_as().trackType() & o2::aod::track::TPCrefit)) { + if (!(V0.posTrack_as().trackType() & o2::aod::track::TPCrefit)) { + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "posTrack " << labelPos << " has no TPC refit"); continue; //TPC refit } - if (!(V0.negTrack_as().trackType() & o2::aod::track::TPCrefit)) { + if (!(V0.negTrack_as().trackType() & o2::aod::track::TPCrefit)) { + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "negTrack " << labelNeg << " has no TPC refit"); continue; //TPC refit } } registry.fill(HIST("hGoodIndices"), 1.5); - if (V0.posTrack_as().tpcNClsCrossedRows() < mincrossedrows) { + if (V0.posTrack_as().tpcNClsCrossedRows() < mincrossedrows) { + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "posTrack " << labelPos << " has " << V0.posTrack_as().tpcNClsCrossedRows() << " crossed rows, cut at " << mincrossedrows); continue; } - if (V0.negTrack_as().tpcNClsCrossedRows() < mincrossedrows) { + if (V0.negTrack_as().tpcNClsCrossedRows() < mincrossedrows) { + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "negTrack " << labelNeg << " has " << V0.negTrack_as().tpcNClsCrossedRows() << " crossed rows, cut at " << mincrossedrows); continue; } registry.fill(HIST("hGoodIndices"), 2.5); - if (fabs(V0.posTrack_as().dcaXY()) < dcapostopv) { + if (fabs(V0.posTrack_as().dcaXY()) < dcapostopv) { + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "posTrack " << labelPos << " has dcaXY " << V0.posTrack_as().dcaXY() << " , cut at " << dcanegtopv); continue; } - if (fabs(V0.negTrack_as().dcaXY()) < dcanegtopv) { + if (fabs(V0.negTrack_as().dcaXY()) < dcanegtopv) { + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "negTrack " << labelNeg << " has dcaXY " << V0.negTrack_as().dcaXY() << " , cut at " << dcanegtopv); continue; } + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "Filling good indices: posTrack --> " << labelPos << ", negTrack --> " << labelNeg); registry.fill(HIST("hGoodIndices"), 3.5); - v0goodindices(V0.posTrack_as().globalIndex(), - V0.negTrack_as().globalIndex(), - V0.posTrack_as().collisionId()); + v0goodindices(V0.posTrack_as().globalIndex(), + V0.negTrack_as().globalIndex(), + V0.posTrack_as().collisionId()); } } }; @@ -149,13 +189,22 @@ struct lambdakzerobuilder { Configurable dcav0dau{"dcav0dau", 1.0, "DCA V0 Daughters"}; Configurable v0radius{"v0radius", 5.0, "v0radius"}; + // for debugging +#ifdef MY_DEBUG + Configurable> v_labelK0Spos{"v_labelK0Spos", {729, 2866, 4754}, "labels of K0S positive daughters, for debug"}; + Configurable> v_labelK0Sneg{"v_labelK0Sneg", {730, 2867, 4755}, "labels of K0S positive daughters, for debug"}; +#endif + double massPi = TDatabasePDG::Instance()->GetParticle(kPiPlus)->Mass(); double massKa = TDatabasePDG::Instance()->GetParticle(kKPlus)->Mass(); double massPr = TDatabasePDG::Instance()->GetParticle(kProton)->Mass(); - using FullTracksExt = soa::Join; - - void process(aod::Collision const& collision, aod::V0GoodIndices const& V0s, soa::Join const& tracks) + void process(aod::Collision const& collision, aod::V0GoodIndices const& V0s, MyTracks const& tracks +#ifdef MY_DEBUG + , + aod::McParticles const& particlesMC +#endif + ) { //Define o2 fitter, 2-prong o2::vertexing::DCAFitterN<2> fitter; @@ -176,11 +225,20 @@ struct lambdakzerobuilder { std::array pvec0 = {0.}; std::array pvec1 = {0.}; +#ifdef MY_DEBUG + auto labelPos = V0.posTrack_as().mcParticleId(); + auto labelNeg = V0.negTrack_as().mcParticleId(); + bool isK0SfromLc = isK0SfromLcFunc(labelPos, labelNeg, v_labelK0Spos, v_labelK0Sneg); +#endif + + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "labelPos = " << labelPos << ", labelNeg = " << labelNeg); + registry.fill(HIST("hV0Candidate"), 0.5); - auto pTrack = getTrackParCov(V0.posTrack_as()); - auto nTrack = getTrackParCov(V0.negTrack_as()); + auto pTrack = getTrackParCov(V0.posTrack_as()); + auto nTrack = getTrackParCov(V0.negTrack_as()); int nCand = fitter.process(pTrack, nTrack); + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "fitter: nCand = " << nCand << " for posTrack --> " << labelPos << ", negTrack --> " << labelNeg); if (nCand != 0) { fitter.propagateTracksToVertex(); const auto& vtx = fitter.getPCACandidate(); @@ -193,33 +251,40 @@ struct lambdakzerobuilder { continue; } + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "in builder 0: posTrack --> " << labelPos << ", negTrack --> " << labelNeg); + //Apply selections so a skimmed table is created only if (fitter.getChi2AtPCACandidate() > dcav0dau) { + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "posTrack --> " << labelPos << ", negTrack --> " << labelNeg << " will be skipped due to dca cut"); continue; } auto V0CosinePA = RecoDecay::CPA(array{collision.posX(), collision.posY(), collision.posZ()}, array{pos[0], pos[1], pos[2]}, array{pvec0[0] + pvec1[0], pvec0[1] + pvec1[1], pvec0[2] + pvec1[2]}); if (V0CosinePA < v0cospa) { + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "posTrack --> " << labelPos << ", negTrack --> " << labelNeg << " will be skipped due to CPA cut"); continue; } auto V0radius = RecoDecay::sqrtSumOfSquares(pos[0], pos[1]); //probably find better name to differentiate the cut from the variable if (V0radius < v0radius) { + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "posTrack --> " << labelPos << ", negTrack --> " << labelNeg << " will be skipped due to radius cut"); continue; } + MY_DEBUG_MSG(isK0SfromLc, LOG(INFO) << "in builder 1, keeping K0S candidate: posTrack --> " << labelPos << ", negTrack --> " << labelNeg); + registry.fill(HIST("hV0Candidate"), 1.5); v0data( - V0.posTrack_as().globalIndex(), - V0.negTrack_as().globalIndex(), - V0.negTrack_as().collisionId(), + V0.posTrack_as().globalIndex(), + V0.negTrack_as().globalIndex(), + V0.negTrack_as().collisionId(), fitter.getTrack(0).getX(), fitter.getTrack(1).getX(), pos[0], pos[1], pos[2], pvec0[0], pvec0[1], pvec0[2], pvec1[0], pvec1[1], pvec1[2], fitter.getChi2AtPCACandidate(), - V0.posTrack_as().dcaXY(), - V0.negTrack_as().dcaXY()); + V0.posTrack_as().dcaXY(), + V0.negTrack_as().dcaXY()); } } }; diff --git a/Analysis/Tasks/Utils/CMakeLists.txt b/Analysis/Tasks/Utils/CMakeLists.txt new file mode 100644 index 0000000000000..c8e0adfb43e32 --- /dev/null +++ b/Analysis/Tasks/Utils/CMakeLists.txt @@ -0,0 +1,13 @@ +# 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. + +o2_add_header_only_library(AnalysisTasksUtils) + + diff --git a/Analysis/Tasks/Utils/include/AnalysisTasksUtils/UtilsDebugLcK0Sp.h b/Analysis/Tasks/Utils/include/AnalysisTasksUtils/UtilsDebugLcK0Sp.h new file mode 100644 index 0000000000000..03758fd2d59bf --- /dev/null +++ b/Analysis/Tasks/Utils/include/AnalysisTasksUtils/UtilsDebugLcK0Sp.h @@ -0,0 +1,82 @@ +// 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 UtilsDebugLcK0Sp +/// \brief Some utilities to do debugging for the LcK0Sp task +/// +/// For example, for: /alice/sim/2020/LHC20l3a/286350/PWGZZ/Run3_Conversion/156_20210308-1000/0001/AO2D.root +/// +/// listLabelsProton = {717, 2810, 4393, 5442, 6769, 7793, 9002, 9789, 10385, 12601, 13819, 17443, 19936, 21715, 23285, 23999, 25676, 28660, 30684, 32049, 34846, 37042, 39278, 40278, 41567, 43696, 44819, 45232, 46578, 47831, 48981, 53527, 57152, 58155, 59296, 61010, 61932, 63577, 65683, 68047, 69657, 70721, 73141, 75320, 76581, 78165, 79504, 80662, 81552, 84208, 86342, 86829, 87347, 89922, 95032, 96087, 97410, 100154, 101743, 103368, 104635, 106089, 108676, 110308, 111335, 111342, 112606, 113469, 114529, 114533, 116388, 117754, 118765, 118784, 119733, 120691, 121773, 123970, 125316, 127253, 129008, 129883, 131408, 132349, 132395, 133263, 134063, 135417, 136194, 137441, 138569, 139167, 141741, 143422, 144322, 145013, 145695, 147134, 148503, 149283, 149728, 153214, 154481, 155193, 158187, 159039, 160009, 161126, 162191, 163774, 165792, 166934, 168768, 174013, 174017, 175467, 177005, 177258, 178380, 179566, 179577, 181145, 187316, 188512, 189094, 191582, 193018, 194159, 195821, 197459, 201840, 202646, 203119, 203763, 205553, 207745, 208851, 211636, 216231, 217125, 217516, 218522, 219477, 219960, 223246, 224677, 224702, 225438, 227194, 230507, 231304, 232070, 232772, 234765, 235877, 236893, 237989, 239575, 241469, 243404, 244872, 245511, 246688, 249625, 250580, 251879, 253031, 254465, 254511, 255917, 256782, 258734, 261436, 262878, 264465, 264467, 268907, 269974, 271856, 274044, 276071, 276915, 278461, 279559, 280441, 281783, 281976, 283405, 284722, 286324, 287929, 289681, 291005, 292324, 293478, 296484, 300536} +/// listLabelsK0SPos = {729, 2866, 4754, 5457, 6891, 7824, 9243, 9810, 10388, 12665, 13830, 17460, 19955, 21786, 23319, 24010, 26234, 28713, 30699, 32056, 34866, 37075, 39314, 40287, 41617, 43790, 44868, 45252, 46587, 47866, 48992, 53534, 57211, 58160, 59355, 61019, 62003, 63580, 65691, 68164, 69721, 70852, 73167, 75331, 76641, 78188, 79595, 80678, 81613, 84311, 86376, 86840, 87374, 89998, 95078, 96094, 97425, 100157, 101788, 103468, 104646, 106155, 108696, 110388, 111355, 111348, 112712, 113578, 114556, 114550, 116432, 117766, 118796, 118790, 119766, 120728, 121872, 124027, 125407, 127284, 129068, 129924, 131453, 132456, 132436, 133266, 134084, 135470, 136207, 137446, 138616, 139217, 141832, 143513, 144416, 145018, 145702, 147139, 148508, 149290, 149737, 153271, 154515, 155206, 158200, 159059, 160082, 161140, 162208, 163777, 165874, 167001, 168852, 174033, 174027, 175501, 177396, 177265, 178383, 179595, 179580, 181168, 187326, 188608, 189128, 191591, 193025, 194170, 195824, 197472, 201905, 202664, 203163, 203835, 205556, 207838, 208866, 211654, 216312, 217132, 217554, 218616, 219503, 219996, 223354, 224713, 224707, 225482, 227207, 230562, 231324, 232083, 232809, 234774, 235892, 236900, 238038, 239629, 241517, 243460, 244881, 245553, 246721, 249672, 250589, 251882, 253034, 254543, 254539, 255941, 256808, 258807, 261612, 262887, 264473, 264471, 268912, 270023, 271998, 274081, 276125, 276952, 278530, 279613, 280446, 282152, 282034, 283408, 284753, 286615, 288037, 289736, 291009, 292336, 293565, 296529, 300564} +/// listLabelsK0SNeg = {730, 2867, 4755, 5458, 6892, 7825, 9244, 9811, 10389, 12666, 13831, 17461, 19956, 21787, 23320, 24011, 26235, 28714, 30700, 32057, 34867, 37076, 39315, 40288, 41618, 43791, 44869, 45253, 46588, 47867, 48993, 53535, 57212, 58161, 59356, 61020, 62004, 63581, 65692, 68165, 69722, 70853, 73168, 75332, 76642, 78189, 79596, 80679, 81614, 84312, 86377, 86841, 87375, 89999, 95079, 96095, 97426, 100158, 101789, 103469, 104647, 106156, 108697, 110389, 111356, 111349, 112713, 113579, 114557, 114551, 116433, 117767, 118797, 118791, 119767, 120729, 121873, 124028, 125408, 127285, 129069, 129925, 131454, 132457, 132437, 133267, 134085, 135471, 136208, 137447, 138617, 139218, 141833, 143514, 144417, 145019, 145703, 147140, 148509, 149291, 149738, 153272, 154516, 155207, 158201, 159060, 160083, 161141, 162209, 163778, 165875, 167002, 168853, 174034, 174028, 175502, 177397, 177266, 178384, 179596, 179581, 181169, 187327, 188609, 189129, 191592, 193026, 194171, 195825, 197473, 201906, 202665, 203164, 203836, 205557, 207839, 208867, 211655, 216313, 217133, 217555, 218617, 219504, 219997, 223355, 224714, 224708, 225483, 227208, 230563, 231325, 232084, 232810, 234775, 235893, 236901, 238039, 239630, 241518, 243461, 244882, 245554, 246722, 249673, 250590, 251883, 253035, 254544, 254540, 255942, 256809, 258808, 261613, 262888, 264474, 264472, 268913, 270024, 271999, 274082, 276126, 276953, 278531, 279614, 280447, 282153, 282035, 283409, 284754, 286616, 288038, 289737, 291010, 292337, 293566, 296530, 300565} +/// + +#include +#include + +inline bool isK0SfromLcFunc(int labelK0SPos, int labelK0SNeg, std::vector listLabelsK0SPos, std::vector listLabelsK0SNeg) +{ + + auto nPositiveDau = listLabelsK0SPos.size(); + auto nNegativeDau = listLabelsK0SNeg.size(); + + // checking sizes of vectors: they should be identical + if (nPositiveDau != nNegativeDau) { + LOG(ERROR) << "Number of elements in vector of positive daughters of K0S different from the number of elements in the vector for the negative ones: " << nPositiveDau << " : " << nNegativeDau; + throw std::runtime_error("sizes of configurables for debug do not match"); + } + + // checking if we find the candidate + bool matchesK0Spositive = std::any_of(listLabelsK0SPos.begin(), listLabelsK0SPos.end(), [&labelK0SPos](const int& label) { return label == labelK0SPos; }); + bool matchesK0Snegative = std::any_of(listLabelsK0SNeg.begin(), listLabelsK0SNeg.end(), [&labelK0SNeg](const int& label) { return label == labelK0SNeg; }); + if (matchesK0Spositive && matchesK0Snegative) { + return true; + } + + return false; +} + +//------------------------------- + +inline bool isProtonFromLcFunc(int labelProton, std::vector listLabelsProton) +{ + // checking if we find the candidate` + bool matchesProton = std::any_of(listLabelsProton.begin(), listLabelsProton.end(), [&labelProton](const int& label) { return label == labelProton; }); + if (matchesProton) { + return true; + } + return false; +} + +//--------------------------------- + +inline bool isLcK0SpFunc(int labelProton, int labelK0SPos, int labelK0SNeg, std::vector listLabelsProton, std::vector listLabelsK0SPos, std::vector listLabelsK0SNeg) +{ + + auto nPositiveDau = listLabelsK0SPos.size(); + auto nNegativeDau = listLabelsK0SNeg.size(); + auto nProtons = listLabelsProton.size(); + + // checking sizes of vectors: they should be identical + if (nPositiveDau != nNegativeDau || nPositiveDau != nProtons) { + LOG(ERROR) << "Number of elements in vector of positive daughters of K0S, in vector of negative daughters of K0S, and in vector of protons differ: " << nPositiveDau << " : " << nNegativeDau << " : " << nProtons; + throw std::runtime_error("sizes of configurables for debug do not match"); + } + + // checking if we find the candidate + bool matchesK0Spositive = std::any_of(listLabelsK0SPos.begin(), listLabelsK0SPos.end(), [&labelK0SPos](const int& label) { return label == labelK0SPos; }); + bool matchesK0Snegative = std::any_of(listLabelsK0SNeg.begin(), listLabelsK0SNeg.end(), [&labelK0SNeg](const int& label) { return label == labelK0SNeg; }); + bool matchesProton = std::any_of(listLabelsProton.begin(), listLabelsProton.end(), [&labelProton](const int& label) { return label == labelProton; }); + if (matchesK0Spositive && matchesK0Snegative && matchesProton) { + return true; + } + + return false; +}