[PWGHF] Tcc via DDPi

PWGHF/TableProducer/CMakeLists.txt

@@ -1,4 +1,4 @@
 # 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.
 #
@@ -244,12 +244,12 @@
                     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter O2Physics::EventFilteringUtils
                     COMPONENT_NAME Analysis)
 
 o2physics_add_dpl_workflow(tree-creator-omegac0-to-omega-ka
                     SOURCES treeCreatorOmegacToOmegaKa.cxx
                     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
 
 o2physics_add_dpl_workflow(tree-creator-omegac0-to-omega-pi
                     SOURCES treeCreatorOmegacToOmegaPi.cxx
                     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
@@ -284,6 +284,11 @@
                     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

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 <biao.zhang@cern.ch>, Heidelberg University
+/// \author Fabrizio Grosa <fabrizio.grosa@cern.ch>, 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<o2::aod::HfCandTccFulls> rowCandidateFull;
+  Produces<o2::aod::HfCandD0FullEvs> rowCandidateFullEvents;
+
+  Configurable<float> ptMinSoftPion{"ptMinSoftPion", 0.4, "Min pt for the softpion"};
+  Configurable<bool> usePionIsGlobalTrackWoDCA{"usePionIsGlobalTrackWoDCA", true, "check isGlobalTrackWoDCA status for pions, for Run3 studies"};
+  Configurable<std::vector<double>> binsPtPion{"binsPtPion", std::vector<double>{hf_cuts_single_track::vecBinsPtTrack}, "track pT bin limits for pion DCA XY pT-dependent cut"};
+  Configurable<LabeledArray<double>> 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<aod::pidTPCFullPi, aod::pidTOFFullPi, aod::pidTPCFullKa, aod::pidTOFFullKa>;
+  using TracksWPid = soa::Join<aod::TracksWCovDcaExtra, TracksPid, aod::TrackSelection>;
+
+  using CollisionsWCentMult = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs>;
+  using SelectedCandidatesMl = soa::Filtered<soa::Join<aod::HfCand2ProngWPid, aod::HfSelD0, aod::HfMlD0>>;
+
+  Filter filterSelectCandidates = aod::hf_sel_candidate_d0::isSelD0 >= 1 || aod::hf_sel_candidate_d0::isSelD0bar >= 1;
+
+  Preslice<SelectedCandidatesMl> candsD0PerCollisionWithMl = aod::track_association::collisionId;
+  Preslice<aod::TrackAssoc> trackIndicesPerCollision = aod::track_association::collisionId;
+  // Partition<SelectedCandidatesMl> candidatesMlAll = aod::hf_sel_candidate_d0::isSelD0 >= 0;
+
+  void init(InitContext const&)
+  {
+    std::array<bool, 1> doprocess{doprocessDataWithML};
+    if (std::accumulate(doprocess.begin(), doprocess.end(), 0) != 1) {
+      LOGP(fatal, "Only one process function can be enabled at a time.");
+    }
+  }
+
+  template <typename T>
+  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 <typename T1, typename T2>
+  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 <typename Coll>
+  float evaluateCentralityColl(const Coll& collision)
+  {
+    return o2::hf_centrality::getCentralityColl<Coll>(collision);
+  }
+
+  template <int reconstructionType, typename CollType, typename CandType, typename TrkType>
+  void runCandCreatorData(CollType const& collision,
+                          CandType const& candidates,
+                          aod::TrackAssoc const& trackIndices,
+                          TrkType const&, aod::BCs const&)
+  {
+    // Filling event properties
+    std::map<int64_t, int64_t> 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<TrkType>();
+
+          auto trackParCovPion = getTrackParCov(trackPion);
+          o2::gpu::gpustd::array<float, 2> 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<float> mlScoresD1;
+          std::vector<float> 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<aod::hf_cand::VertexerType::DCAFitter>(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<HfTreeCreatorTccToD0D0Pi>(cfgc));
+  return workflow;
+}