diff --git a/PWGHF/TableProducer/CMakeLists.txt b/PWGHF/TableProducer/CMakeLists.txt index a9d71f86da1..9738128a300 100644 --- a/PWGHF/TableProducer/CMakeLists.txt +++ b/PWGHF/TableProducer/CMakeLists.txt @@ -284,6 +284,11 @@ o2physics_add_dpl_workflow(tree-creator-dstar-to-d0-pi PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(tree-creator-tcc-to-d0-d0-pi + SOURCES treeCreatorTccToD0D0Pi.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + # Derived-data creators o2physics_add_dpl_workflow(derived-data-creator-bplus-to-d0-pi diff --git a/PWGHF/TableProducer/treeCreatorTccToD0D0Pi.cxx b/PWGHF/TableProducer/treeCreatorTccToD0D0Pi.cxx new file mode 100644 index 00000000000..4af2614a32c --- /dev/null +++ b/PWGHF/TableProducer/treeCreatorTccToD0D0Pi.cxx @@ -0,0 +1,378 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// 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 treeCreatorTccToD0D0Pi.cxx +/// \brief tree creator for studying the charm exotic state Tcc to D0D0pi +/// \author Biao Zhang , Heidelberg University +/// \author Fabrizio Grosa , CERN + +#include "CommonConstants/PhysicsConstants.h" +#include "DCAFitter/DCAFitterN.h" +#include "Framework/AnalysisTask.h" +#include "Framework/O2DatabasePDGPlugin.h" +#include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/DCA.h" + +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/CollisionAssociationTables.h" + +#include "PWGHF/Core/HfHelper.h" +#include "PWGHF/Core/CentralityEstimation.h" +#include "PWGHF/DataModel/CandidateReconstructionTables.h" +#include "PWGHF/DataModel/CandidateSelectionTables.h" +#include "PWGHF/Utils/utilsBfieldCCDB.h" +#include "PWGHF/Utils/utilsEvSelHf.h" +#include "PWGHF/D2H/DataModel/ReducedDataModel.h" +#include "PWGHF/Utils/utilsTrkCandHf.h" +#include "PWGHF/D2H/Utils/utilsRedDataFormat.h" + +using namespace o2; +using namespace o2::analysis; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::hf_trkcandsel; + +namespace o2::aod +{ +namespace full +{ +DECLARE_SOA_INDEX_COLUMN(Collision, collision); +DECLARE_SOA_COLUMN(PtD1, ptD1, float); +DECLARE_SOA_COLUMN(PtD2, ptD2, float); +DECLARE_SOA_COLUMN(PtPi, ptPi, float); +DECLARE_SOA_COLUMN(PxProng0D1, pxProng0D1, float); +DECLARE_SOA_COLUMN(PxProng1D1, pxProng1D1, float); +DECLARE_SOA_COLUMN(PyProng0D1, pyProng0D1, float); +DECLARE_SOA_COLUMN(PyProng1D1, pyProng1D1, float); +DECLARE_SOA_COLUMN(PzProng0D1, pzProng0D1, float); +DECLARE_SOA_COLUMN(PzProng1D1, pzProng1D1, float); +DECLARE_SOA_COLUMN(PxProng0D2, pxProng0D2, float); +DECLARE_SOA_COLUMN(PxProng1D2, pxProng1D2, float); +DECLARE_SOA_COLUMN(PyProng0D2, pyProng0D2, float); +DECLARE_SOA_COLUMN(PyProng1D2, pyProng1D2, float); +DECLARE_SOA_COLUMN(PzProng0D2, pzProng0D2, float); +DECLARE_SOA_COLUMN(PzProng1D2, pzProng1D2, float); +DECLARE_SOA_COLUMN(PxSoftPi, pxSoftPi, float); +DECLARE_SOA_COLUMN(PySoftPi, pySoftPi, float); +DECLARE_SOA_COLUMN(PzSoftPi, pzSoftPi, float); +DECLARE_SOA_COLUMN(SelFlagD1, selFlagD1, int8_t); +DECLARE_SOA_COLUMN(SelFlagD2, selFlagD2, int8_t); +DECLARE_SOA_COLUMN(MD1, mD1, float); +DECLARE_SOA_COLUMN(MD2, mD2, float); +DECLARE_SOA_COLUMN(EtaD1, etaD1, float); +DECLARE_SOA_COLUMN(EtaD2, etaD2, float); +DECLARE_SOA_COLUMN(EtaSoftPi, etaSoftPi, float); +DECLARE_SOA_COLUMN(PhiD1, phiD1, float); +DECLARE_SOA_COLUMN(PhiD2, phiD2, float); +DECLARE_SOA_COLUMN(PhiSoftPi, phiSoftPi, float); +DECLARE_SOA_COLUMN(YD1, yD1, float); +DECLARE_SOA_COLUMN(YD2, yD2, float); +DECLARE_SOA_COLUMN(YSoftPi, ySoftPi, float); +DECLARE_SOA_COLUMN(NSigTpcSoftPi, nSigTpcSoftPi, float); +DECLARE_SOA_COLUMN(NSigTofSoftPi, nSigTofSoftPi, float); +DECLARE_SOA_COLUMN(MLScoreD1, mlScoreD1, float); +DECLARE_SOA_COLUMN(MLScoreD2, mlScoreD2, float); +DECLARE_SOA_COLUMN(DecayLengthD1, decayLengthD1, float); +DECLARE_SOA_COLUMN(DecayLengthD2, decayLengthD2, float); +DECLARE_SOA_COLUMN(CpaD1, cpaD1, float); +DECLARE_SOA_COLUMN(CpaD2, cpaD2, float); +DECLARE_SOA_COLUMN(SignSoftPi, signSoftPi, float); +DECLARE_SOA_COLUMN(DcaXYSoftPi, dcaXYSoftPi, float); +DECLARE_SOA_COLUMN(DcaZSoftPi, dcaZSoftPi, float); +DECLARE_SOA_COLUMN(NITSClsSoftPi, nitsClsSoftPi, float); +DECLARE_SOA_COLUMN(NTPCClsCrossedRowsSoftPi, ntpcClsCrossedRows, float); +DECLARE_SOA_COLUMN(NTPCChi2NClSoftPi, ntpcChi2NClSoftPi, float); +DECLARE_SOA_COLUMN(CentOfCand, centOfCand, float); +// Events +DECLARE_SOA_COLUMN(IsEventReject, isEventReject, int); +DECLARE_SOA_COLUMN(RunNumber, runNumber, int); +} // namespace full + +DECLARE_SOA_TABLE(HfCandTccFulls, "AOD", "HFCANDTCCFULL", + full::CollisionId, + full::PtD1, + full::PtD2, + full::PtPi, + full::PxProng0D1, + full::PxProng1D1, + full::PyProng0D1, + full::PyProng1D1, + full::PzProng0D1, + full::PzProng1D1, + full::PxProng0D2, + full::PxProng1D2, + full::PyProng0D2, + full::PyProng1D2, + full::PzProng0D2, + full::PzProng1D2, + full::PxSoftPi, + full::PySoftPi, + full::PzSoftPi, + full::SelFlagD1, + full::SelFlagD2, + full::MD1, + full::MD2, + full::EtaD1, + full::EtaD2, + full::EtaSoftPi, + full::PhiD1, + full::PhiD2, + full::PhiSoftPi, + full::YD1, + full::YD2, + full::YSoftPi, + full::NSigTpcSoftPi, + full::NSigTofSoftPi, + full::MLScoreD1, + full::MLScoreD2, + full::DecayLengthD1, + full::DecayLengthD2, + full::CpaD1, + full::CpaD2, + full::SignSoftPi, + full::DcaXYSoftPi, + full::DcaZSoftPi, + full::NITSClsSoftPi, + full::NTPCClsCrossedRowsSoftPi, + full::NTPCChi2NClSoftPi, + full::CentOfCand); + +DECLARE_SOA_TABLE(HfCandD0FullEvs, "AOD", "HFCANDD0FULLEV", + full::CollisionId, + collision::NumContrib, + collision::PosX, + collision::PosY, + collision::PosZ, + full::IsEventReject, + full::RunNumber); +} // namespace o2::aod + +/// Writes the full information in an output TTree +struct HfTreeCreatorTccToD0D0Pi { + Produces rowCandidateFull; + Produces rowCandidateFullEvents; + + Configurable ptMinSoftPion{"ptMinSoftPion", 0.4, "Min pt for the softpion"}; + Configurable usePionIsGlobalTrackWoDCA{"usePionIsGlobalTrackWoDCA", true, "check isGlobalTrackWoDCA status for pions, for Run3 studies"}; + Configurable> binsPtPion{"binsPtPion", std::vector{hf_cuts_single_track::vecBinsPtTrack}, "track pT bin limits for pion DCA XY pT-dependent cut"}; + Configurable> cutsTrackPionDCA{"cutsTrackPionDCA", {hf_cuts_single_track::cutsTrack[0], hf_cuts_single_track::nBinsPtTrack, hf_cuts_single_track::nCutVarsTrack, hf_cuts_single_track::labelsPtTrack, hf_cuts_single_track::labelsCutVarTrack}, "Single-track selections per pT bin for pions"}; + + HfHelper hfHelper; + + using TracksPid = soa::Join; + using TracksWPid = soa::Join; + + using CollisionsWCentMult = soa::Join; + using SelectedCandidatesMl = soa::Filtered>; + + Filter filterSelectCandidates = aod::hf_sel_candidate_d0::isSelD0 >= 1 || aod::hf_sel_candidate_d0::isSelD0bar >= 1; + + Preslice candsD0PerCollisionWithMl = aod::track_association::collisionId; + Preslice trackIndicesPerCollision = aod::track_association::collisionId; + // Partition candidatesMlAll = aod::hf_sel_candidate_d0::isSelD0 >= 0; + + void init(InitContext const&) + { + std::array doprocess{doprocessDataWithML}; + if (std::accumulate(doprocess.begin(), doprocess.end(), 0) != 1) { + LOGP(fatal, "Only one process function can be enabled at a time."); + } + } + + template + void fillEvent(const T& collision, int isEventReject, int runNumber) + { + rowCandidateFullEvents( + collision.globalIndex(), + collision.numContrib(), + collision.posX(), + collision.posY(), + collision.posZ(), + isEventReject, + runNumber); + } + + /// \param trackPion is a track with the pion hypothesis + /// \param dcaPion is the 2-D array with track DCAs of the pion + /// \return true if trackPion passes all cuts + template + bool isPionSelected(const T1& trackPion, const T2& dcaPion) + { + // check isGlobalTrackWoDCA status for pions if wanted + if (usePionIsGlobalTrackWoDCA && !trackPion.isGlobalTrackWoDCA()) { + return false; + } + float pTBinTrack = trackPion.pt(); + if (pTBinTrack < ptMinSoftPion) { + return false; + } + if (std::abs(dcaPion[0]) < cutsTrackPionDCA->get(pTBinTrack, "min_dcaxytoprimary")) { + return false; // minimum DCAxy + } + if (std::abs(dcaPion[0]) > cutsTrackPionDCA->get(pTBinTrack, "max_dcaxytoprimary")) { + return false; // maximum DCAxy + } + + return true; + } + + /// Evaluate centrality/multiplicity percentile (centrality estimator is automatically selected based on the used table) + /// \param candidate is candidate + /// \return centrality/multiplicity percentile of the collision + template + float evaluateCentralityColl(const Coll& collision) + { + return o2::hf_centrality::getCentralityColl(collision); + } + + template + void runCandCreatorData(CollType const& collision, + CandType const& candidates, + aod::TrackAssoc const& trackIndices, + TrkType const&, aod::BCs const&) + { + // Filling event properties + std::map selectedTracksPion; + + for (const auto& candidateD1 : candidates) { + for (auto candidateD2 = candidateD1 + 1; candidateD2 != candidates.end(); ++candidateD2) { + for (const auto& trackId : trackIndices) { + auto trackPion = trackId.template track_as(); + + auto trackParCovPion = getTrackParCov(trackPion); + o2::gpu::gpustd::array dcaPion{trackPion.dcaXY(), trackPion.dcaZ()}; + // apply selections on pion tracks + if (!isPionSelected(trackPion, dcaPion)) { + continue; + } + // avoid shared tracks + if ( + (candidateD1.prong0Id() == candidateD2.prong0Id()) || + (candidateD1.prong0Id() == candidateD2.prong1Id()) || + (candidateD1.prong1Id() == candidateD2.prong0Id()) || + (candidateD1.prong1Id() == candidateD2.prong1Id()) || + (candidateD1.prong0Id() == trackPion.globalIndex()) || + (candidateD1.prong1Id() == trackPion.globalIndex()) || + (candidateD2.prong0Id() == trackPion.globalIndex()) || + (candidateD2.prong1Id() == trackPion.globalIndex())) { + continue; + } + // Retrieve properties of the two D0 candidates + float yD1 = hfHelper.yD0(candidateD1); + float yD2 = hfHelper.yD0(candidateD2); + float massD01 = -999; + float massD02 = -999; + int candFlagD1 = -999; + int candFlagD2 = -999; + float cent = evaluateCentralityColl(collision); + + std::vector mlScoresD1; + std::vector mlScoresD2; + + if (candidateD1.isSelD0()) { + candFlagD1 = (candidateD1.isSelD0bar()) ? 3 : 1; + std::copy(candidateD1.mlProbD0().begin(), candidateD1.mlProbD0().end(), std::back_inserter(mlScoresD1)); + massD01 = hfHelper.invMassD0ToPiK(candidateD1); + } else if (candidateD1.isSelD0bar()) { + candFlagD1 = 2; + std::copy(candidateD1.mlProbD0bar().begin(), candidateD1.mlProbD0bar().end(), std::back_inserter(mlScoresD1)); + massD01 = hfHelper.invMassD0barToKPi(candidateD1); + } + + if (candidateD2.isSelD0()) { + candFlagD2 = (candidateD2.isSelD0bar()) ? 3 : 1; + std::copy(candidateD2.mlProbD0().begin(), candidateD2.mlProbD0().end(), std::back_inserter(mlScoresD2)); + massD02 = hfHelper.invMassD0ToPiK(candidateD2); + + } else if (candidateD2.isSelD0bar()) { + candFlagD2 = 2; + std::copy(candidateD2.mlProbD0bar().begin(), candidateD2.mlProbD0bar().end(), std::back_inserter(mlScoresD2)); + massD02 = hfHelper.invMassD0barToKPi(candidateD2); + } + + rowCandidateFull( + candidateD1.collisionId(), + candidateD1.pt(), + candidateD2.pt(), + trackPion.pt(), + candidateD1.pxProng0(), + candidateD1.pxProng1(), + candidateD1.pyProng0(), + candidateD1.pyProng1(), + candidateD1.pzProng0(), + candidateD1.pzProng1(), + candidateD2.pxProng0(), + candidateD2.pxProng1(), + candidateD2.pyProng0(), + candidateD2.pyProng1(), + candidateD2.pzProng0(), + candidateD2.pzProng1(), + trackPion.px(), + trackPion.py(), + trackPion.pz(), + candFlagD1, + candFlagD2, + massD01, + massD02, + candidateD1.eta(), + candidateD2.eta(), + trackPion.eta(), + candidateD1.phi(), + candidateD2.phi(), + trackPion.phi(), + yD1, + yD2, + trackPion.y(), + trackPion.tpcNSigmaPi(), + trackPion.tofNSigmaPi(), + mlScoresD1[0], + mlScoresD2[0], + candidateD1.decayLength(), + candidateD2.decayLength(), + candidateD1.cpa(), + candidateD2.cpa(), + trackPion.sign(), + dcaPion[0], + dcaPion[1], + trackPion.itsNCls(), + trackPion.tpcNClsCrossedRows(), + trackPion.tpcChi2NCl()); + cent); + } // end of loop track + } // end of loop second D0 + } // end of loop first D0 + } + + void processDataWithML(CollisionsWCentMult const& collisions, + SelectedCandidatesMl const& candidates, + aod::TrackAssoc const& trackIndices, + TracksWPid const& tracks, + aod::BCs const& bcs) + { + for (const auto& collision : collisions) { + rowCandidateFullEvents.reserve(collisions.size()); + fillEvent(collision, 0, collision.bc().runNumber()); + auto thisCollId = collision.globalIndex(); + auto candwD0ThisColl = candidates.sliceBy(candsD0PerCollisionWithMl, thisCollId); + if (candwD0ThisColl.size() <= 1) + continue; // only loop the collision that include at least 2 D candidates + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); + runCandCreatorData(collision, candwD0ThisColl, trackIdsThisCollision, tracks, bcs); + } + } + PROCESS_SWITCH(HfTreeCreatorTccToD0D0Pi, processDataWithML, "Process data with DCAFitterN", true); +}; +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow; + workflow.push_back(adaptAnalysisTask(cfgc)); + return workflow; +}