diff --git a/PWGHF/HFC/DataModel/DerivedDataCorrelationTables.h b/PWGHF/HFC/DataModel/DerivedDataCorrelationTables.h index 39a3231977c..d48cd0f9d10 100644 --- a/PWGHF/HFC/DataModel/DerivedDataCorrelationTables.h +++ b/PWGHF/HFC/DataModel/DerivedDataCorrelationTables.h @@ -24,6 +24,7 @@ namespace hf_collisions_reduced { DECLARE_SOA_COLUMN(NumPvContrib, numPvContrib, int); //! Event multiplicity from PV contributors DECLARE_SOA_COLUMN(Multiplicity, multiplicity, float); //! Event multiplicity +DECLARE_SOA_COLUMN(Centrality, centrality, float); //! Event centrality DECLARE_SOA_COLUMN(PosZ, posZ, float); //! Primary vertex z position } // namespace hf_collisions_reduced @@ -34,6 +35,13 @@ DECLARE_SOA_TABLE(HfcRedCollisions, "AOD", "HFCREDCOLLISION", //! Table with col aod::hf_collisions_reduced::NumPvContrib, aod::hf_collisions_reduced::PosZ); +DECLARE_SOA_TABLE(HfcRedFlowColls, "AOD", "HFCREDFLOWCOLL", //! Table with collision info + soa::Index<>, + aod::hf_collisions_reduced::Multiplicity, + aod::hf_collisions_reduced::NumPvContrib, + aod::hf_collisions_reduced::Centrality, + aod::hf_collisions_reduced::PosZ); + using HfcRedCollision = HfcRedCollisions::iterator; // DECLARE_SOA_TABLE(HfCandColCounts, "AOD", "HFCANDCOLCOUNT", //! Table with number of collisions which contain at least one candidate @@ -41,16 +49,20 @@ using HfcRedCollision = HfcRedCollisions::iterator; namespace hf_candidate_reduced { -DECLARE_SOA_INDEX_COLUMN(HfcRedCollision, hfcRedCollision); //! ReducedCollision index -DECLARE_SOA_COLUMN(Prong0Id, prong0Id, int); //! Prong 0 index -DECLARE_SOA_COLUMN(Prong1Id, prong1Id, int); //! Prong 1 index -DECLARE_SOA_COLUMN(Prong2Id, prong2Id, int); //! Prong2 index -DECLARE_SOA_COLUMN(PhiCand, phiCand, float); //! Phi of the candidate -DECLARE_SOA_COLUMN(EtaCand, etaCand, float); //! Eta of the candidate -DECLARE_SOA_COLUMN(PtCand, ptCand, float); //! Pt of the candidate -DECLARE_SOA_COLUMN(InvMassDs, invMassDs, float); //! Invariant mass of Ds candidate -DECLARE_SOA_COLUMN(BdtScorePrompt, bdtScorePrompt, float); //! BDT output score for prompt hypothesis -DECLARE_SOA_COLUMN(BdtScoreBkg, bdtScoreBkg, float); //! BDT output score for backgronud hypothesis +DECLARE_SOA_INDEX_COLUMN(HfcRedCollision, hfcRedCollision); //! ReducedCollision index +DECLARE_SOA_INDEX_COLUMN(HfcRedFlowColl, hfcRedFlowColl); //! ReducedCollision index +DECLARE_SOA_COLUMN(Prong0Id, prong0Id, int); //! Prong 0 index +DECLARE_SOA_COLUMN(Prong1Id, prong1Id, int); //! Prong 1 index +DECLARE_SOA_COLUMN(Prong2Id, prong2Id, int); //! Prong2 index +DECLARE_SOA_COLUMN(PhiCand, phiCand, float); //! Phi of the candidate +DECLARE_SOA_COLUMN(EtaCand, etaCand, float); //! Eta of the candidate +DECLARE_SOA_COLUMN(PtCand, ptCand, float); //! Pt of the candidate +DECLARE_SOA_COLUMN(InvMassDs, invMassDs, float); //! Invariant mass of Ds candidate +DECLARE_SOA_COLUMN(InvMassCharmHad, invMassCharmHad, float); //! Invariant mass of CharmHad candidate +DECLARE_SOA_COLUMN(BdtScorePrompt, bdtScorePrompt, float); //! BDT output score for prompt hypothesis +DECLARE_SOA_COLUMN(BdtScoreBkg, bdtScoreBkg, float); //! BDT output score for backgronud hypothesis +DECLARE_SOA_COLUMN(BdtScore0, bdtScore0, float); //! First BDT output score +DECLARE_SOA_COLUMN(BdtScore1, bdtScore1, float); //! Second BDT output score } // namespace hf_candidate_reduced DECLARE_SOA_TABLE(DsCandReduceds, "AOD", "DSCANDREDUCED", //! Table with Ds candidate info soa::Index<>, @@ -69,6 +81,23 @@ DECLARE_SOA_TABLE(DsCandSelInfos, "AOD", "DSCANDSELINFO", //! Table with Ds cand aod::hf_candidate_reduced::BdtScorePrompt, aod::hf_candidate_reduced::BdtScoreBkg); +DECLARE_SOA_TABLE(HfcRedCharmHads, "AOD", "HFCREDCHARMHAD", //! Table with charm hadron candidate info + soa::Index<>, + aod::hf_candidate_reduced::HfcRedFlowCollId, + aod::hf_candidate_reduced::PhiCand, + aod::hf_candidate_reduced::EtaCand, + aod::hf_candidate_reduced::PtCand, + aod::hf_candidate_reduced::InvMassCharmHad, + aod::hf_candidate_reduced::Prong0Id, + aod::hf_candidate_reduced::Prong1Id, + aod::hf_candidate_reduced::Prong2Id); + +DECLARE_SOA_TABLE(HfcRedCharmMls, "AOD", "HFCREDCHARMML", //! Table with charm hadron candidate selection info + soa::Index<>, + aod::hf_candidate_reduced::HfcRedFlowCollId, + aod::hf_candidate_reduced::BdtScore0, + aod::hf_candidate_reduced::BdtScore1); + namespace hf_assoc_track_reduced { DECLARE_SOA_COLUMN(OriginTrackId, originTrackId, int); //! Original track index @@ -97,6 +126,23 @@ DECLARE_SOA_TABLE(AssocTrackSels, "AOD", "ASSOCTRACKSEL", //! Table with associa aod::hf_assoc_track_reduced::ItsNCls, aod::hf_assoc_track_reduced::DcaXY, aod::hf_assoc_track_reduced::DcaZ) + +DECLARE_SOA_TABLE(HfcRedTrkAssoc, "AOD", "HFCREDTRKASSOC", //! Table with associated track info + soa::Index<>, + aod::hf_candidate_reduced::HfcRedFlowCollId, + aod::hf_assoc_track_reduced::OriginTrackId, + aod::hf_assoc_track_reduced::PhiAssocTrack, + aod::hf_assoc_track_reduced::EtaAssocTrack, + aod::hf_assoc_track_reduced::PtAssocTrack); + +DECLARE_SOA_TABLE(HfcRedTrkSels, "AOD", "HFCREDTRKSELS", //! Table with associated track info + soa::Index<>, + aod::hf_candidate_reduced::HfcRedFlowCollId, + aod::hf_assoc_track_reduced::NTpcCrossedRows, + aod::hf_assoc_track_reduced::ItsClusterMap, + aod::hf_assoc_track_reduced::ItsNCls, + aod::hf_assoc_track_reduced::DcaXY, + aod::hf_assoc_track_reduced::DcaZ) } // namespace o2::aod #endif // PWGHF_HFC_DATAMODEL_DERIVEDDATACORRELATIONTABLES_H_ diff --git a/PWGHF/HFC/TableProducer/CMakeLists.txt b/PWGHF/HFC/TableProducer/CMakeLists.txt index 99f2c4fb152..eff0f6557b3 100644 --- a/PWGHF/HFC/TableProducer/CMakeLists.txt +++ b/PWGHF/HFC/TableProducer/CMakeLists.txt @@ -44,6 +44,11 @@ o2physics_add_dpl_workflow(correlator-ds-hadrons PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(correlator-flow-charm-hadrons + SOURCES correlatorFlowCharmHadrons.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::EventFilteringUtils + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(correlator-ds-hadrons-reduced SOURCES correlatorDsHadronsReduced.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore diff --git a/PWGHF/HFC/TableProducer/correlatorFlowCharmHadrons.cxx b/PWGHF/HFC/TableProducer/correlatorFlowCharmHadrons.cxx new file mode 100644 index 00000000000..ec5f2a40c19 --- /dev/null +++ b/PWGHF/HFC/TableProducer/correlatorFlowCharmHadrons.cxx @@ -0,0 +1,272 @@ +// 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 correlatorFlowCharmHadrons.cxx +/// \brief CharmHadrons-Hadrons correlator tree creator for data and MC-reco analyses +/// \author Marcello Di Costanzo , Politecnico and INFN Torino +/// \author Stefano Politanò , CERN + +#include "PWGHF/Core/HfHelper.h" +#include "PWGHF/DataModel/CandidateReconstructionTables.h" +#include "PWGHF/DataModel/CandidateSelectionTables.h" +#include "PWGHF/HFC/DataModel/DerivedDataCorrelationTables.h" +#include "PWGHF/Utils/utilsEvSelHf.h" + +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +using namespace o2; +using namespace o2::hf_centrality; +using namespace o2::hf_evsel; + +enum DecayChannel { + DplusToPiKPi = 0, + DsToKKPi, + DsToPiKK +}; + +/// Code to select collisions with at least one Ds meson +struct HfCorrelatorFlowCharmHadrons { + Produces rowCollisions; + Produces rowCharmCandidates; + Produces rowCharmCandidatesMl; + Produces rowAssocTrackReduced; + Produces rowAssocTrackSelInfo; + + Configurable centEstimator{"centEstimator", 2, "Centrality estimation (FT0A: 1, FT0C: 2, FT0M: 3, FV0A: 4)"}; + Configurable selectionFlag{"selectionFlag", 1, "Selection Flag for hadron (e.g. 1 for skimming, 3 for topo. and kine., 7 for PID)"}; + Configurable forceCharmInCollision{"forceCharmInCollision", false, "Flag to force charm in collision"}; + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable> classMl{"classMl", {0, 2}, "Indexes of BDT scores to be stored. Two indexes max."}; + Configurable yCandMax{"yCandMax", 0.8, "max. cand. rapidity"}; + Configurable ptCandMin{"ptCandMin", 1., "min. cand. pT"}; + Configurable ptCandMax{"ptCandMax", 24., "max. cand. pT"}; + Configurable etaTrackMax{"etaTrackMax", 1., "max. track eta"}; + Configurable ptTrackMin{"ptTrackMin", 0.15, "min. track pT"}; + Configurable ptTrackMax{"ptTrackMax", 5., "max. track pT"}; + Configurable dcaXYTrackMax{"dcaXYTrackMax", 1., "max. track DCA XY"}; + Configurable dcaZTrackMax{"dcaZTrackMax", 1., "max. track DCA Z"}; + + HfHelper hfHelper; + HfEventSelection hfEvSel; // event selection and monitoring + o2::framework::Service ccdb; + SliceCache cache; + + double massCharm{0.}; + + using CollsWithCentMult = soa::Join; + using CandDsDataWMl = soa::Filtered>; + using CandDplusDataWMl = soa::Filtered>; + using TracksData = soa::Filtered>; + + Filter filterSelectDsCandidates = aod::hf_sel_candidate_ds::isSelDsToKKPi >= selectionFlag || aod::hf_sel_candidate_ds::isSelDsToPiKK >= selectionFlag; + Filter filterSelectDplusCandidates = aod::hf_sel_candidate_dplus::isSelDplusToPiKPi >= selectionFlag; + Filter filterSelectTrackData = (nabs(aod::track::eta) < etaTrackMax) && (aod::track::pt > ptTrackMin) && (aod::track::pt < ptTrackMax) && (nabs(aod::track::dcaXY) < dcaXYTrackMax) && (nabs(aod::track::dcaZ) < dcaZTrackMax); + + Preslice candsDsPerCollWMl = aod::hf_cand::collisionId; + Preslice candsDplusPerCollWMl = aod::hf_cand::collisionId; + Preslice trackIndicesPerColl = aod::track::collisionId; + + Partition selectedDsToKKPiWMl = aod::hf_sel_candidate_ds::isSelDsToKKPi >= selectionFlag; + Partition selectedDsToPiKKWMl = aod::hf_sel_candidate_ds::isSelDsToPiKK >= selectionFlag; + + HistogramRegistry registry{"registry", {}}; + + void init(InitContext&) + { + if (doprocessDplusWithMl) { + massCharm = o2::constants::physics::MassDPlus; + } else if (doprocessDsWithMl) { + massCharm = o2::constants::physics::MassDS; + } + + hfEvSel.addHistograms(registry); // collision monitoring + ccdb->setURL(ccdbUrl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + }; // end init + + /// Check event selections for collision and fill the collision table + /// \param collision is the collision + template + bool checkAndFillCollision(Coll const& collision) + { + float cent{-1.f}; + float mult{-1.f}; + o2::hf_evsel::HfCollisionRejectionMask collRejMask{}; + if (centEstimator == CentralityEstimator::FT0A) { + collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); + mult = collision.multFT0A(); + } else if (centEstimator == CentralityEstimator::FT0C) { + collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); + mult = collision.multFT0C(); + } else if (centEstimator == CentralityEstimator::FT0M) { + collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); + mult = collision.multFT0M(); + } else if (centEstimator == CentralityEstimator::FV0A) { + collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, ccdb, registry); + mult = collision.multFV0A(); + } else { + LOG(fatal) << "Centrality estimator not recognized for collision selection"; + std::abort(); + } + hfEvSel.fillHistograms(collision, collRejMask, cent); + if (collRejMask != 0) { + return false; + } + rowCollisions(mult, collision.numContrib(), cent, collision.posZ()); + return true; + } + + /// Get charm hadron candidate mass + /// \param candidate is the charm hadron candidate + template + double getCandMass(const TCand& candidate) + { + if constexpr (channel == DecayChannel::DsToKKPi) { + return hfHelper.invMassDsToKKPi(candidate); + } + if constexpr (channel == DecayChannel::DsToPiKK) { + return hfHelper.invMassDsToPiKK(candidate); + } + if constexpr (channel == DecayChannel::DplusToPiKPi) { + return hfHelper.invMassDplusToPiKPi(candidate); + } + return -1.; + } + + /// Get charm hadron bdt scores + /// \param candidate is the charm hadron candidate + template + std::vector getCandMlScores(const TCand& candidate) + { + std::vector outputMl{-999., -999.}; + if constexpr (channel == DecayChannel::DsToKKPi) { + for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { + outputMl[iclass] = candidate.mlProbDsToKKPi()[classMl->at(iclass)]; + } + } + if constexpr (channel == DecayChannel::DsToPiKK) { + for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { + outputMl[iclass] = candidate.mlProbDsToPiKK()[classMl->at(iclass)]; + } + } + if constexpr (channel == DecayChannel::DplusToPiKPi) { + for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { + outputMl[iclass] = candidate.mlProbDplusToPiKPi()[classMl->at(iclass)]; + } + } + return outputMl; + } + + /// Fill charm hadron tables + /// \param candidates are the selected charm hadron candidates + template + void fillCharmHadronTables(TCand const& candidates) + { + int indexRedColl = rowCollisions.lastIndex(); + for (const auto& candidate : candidates) { + if (std::abs(candidate.y(massCharm)) > yCandMax || candidate.pt() < ptCandMin || candidate.pt() > ptCandMax) { + continue; + } + double massCand = getCandMass(candidate); + rowCharmCandidates(indexRedColl, candidate.phi(), candidate.eta(), candidate.pt(), massCand, candidate.prong0Id(), candidate.prong1Id(), candidate.prong2Id()); + + std::vector outputMl = getCandMlScores(candidate); + rowCharmCandidatesMl(indexRedColl, outputMl[0], outputMl[1]); + } + } + + /// Fill tracks tables + /// \param tracks are the selected tracks + template + void fillTracksTables(TTrack const& tracks) + { + int indexRedColl = rowCollisions.lastIndex(); + for (const auto& track : tracks) { + if (!track.isGlobalTrackWoDCA()) { + continue; + } + rowAssocTrackReduced(indexRedColl, track.globalIndex(), track.phi(), track.eta(), track.pt()); + rowAssocTrackSelInfo(indexRedColl, track.tpcNClsCrossedRows(), track.itsClusterMap(), track.itsNCls(), track.dcaXY(), track.dcaZ()); + } + } + + // Dplus with ML selections + void processDplusWithMl(CollsWithCentMult const& colls, + CandDplusDataWMl const& candsDplus, + TracksData const& tracks) + { + for (const auto& coll : colls) { + auto thisCollId = coll.globalIndex(); + auto candsCThisColl = candsDplus.sliceBy(candsDplusPerCollWMl, thisCollId); + if (forceCharmInCollision && candsCThisColl.size() < 1) { + continue; + } + if (!checkAndFillCollision(coll)) { + continue; + } + auto trackIdsThisColl = tracks.sliceBy(trackIndicesPerColl, thisCollId); + fillCharmHadronTables(candsCThisColl); + fillTracksTables(trackIdsThisColl); + } + } + PROCESS_SWITCH(HfCorrelatorFlowCharmHadrons, processDplusWithMl, "Process Dplus candidates with ML info", false); + + // Ds with ML selections + void processDsWithMl(CollsWithCentMult const& colls, + TracksData const& tracks, + CandDsDataWMl const&) + { + for (const auto& coll : colls) { + auto thisCollId = coll.globalIndex(); + auto candsDsToKKPiWMl = selectedDsToKKPiWMl->sliceByCached(aod::hf_cand::collisionId, thisCollId, cache); + auto candsDsToPiKKWMl = selectedDsToPiKKWMl->sliceByCached(aod::hf_cand::collisionId, thisCollId, cache); + if (forceCharmInCollision && candsDsToKKPiWMl.size() < 1 && candsDsToPiKKWMl.size() < 1) { + continue; + } + if (!checkAndFillCollision(coll)) { + continue; + } + auto trackIdsThisColl = tracks.sliceBy(trackIndicesPerColl, thisCollId); + fillCharmHadronTables(candsDsToPiKKWMl); + fillCharmHadronTables(candsDsToKKPiWMl); + fillTracksTables(trackIdsThisColl); + } + } + PROCESS_SWITCH(HfCorrelatorFlowCharmHadrons, processDsWithMl, "Process Ds candidates with ML info", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +}