[PWGHF] ppbar decay channel for hidden charm

PWGHF/D2H/DataModel/ReducedDataModel.h

@@ -196,6 +196,7 @@ DECLARE_SOA_COLUMN(TpcNClsCrossedRowsProngMin, tpcNClsCrossedRowsProngMin, int);
 DECLARE_SOA_COLUMN(TpcChi2NClProngMax, tpcChi2NClProngMax, float);               //! maximum value of TPC chi2 for the decay daughter tracks
 DECLARE_SOA_COLUMN(PtProngMin, ptProngMin, float);                               //! minimum value of transverse momentum for the decay daughter tracks
 DECLARE_SOA_COLUMN(AbsEtaProngMin, absEtaProngMin, float);                       //! minimum value of absolute pseudorapidity for the decay daughter tracks
+DECLARE_SOA_COLUMN(TrackType, trackType, uint8_t);                               //! particle type according to PID selection (pion, kaon, proton)
 
 // dynamic columns
 DECLARE_SOA_DYNAMIC_COLUMN(Pt, pt, //! transverse momentum
@@ -1653,6 +1654,19 @@ DECLARE_SOA_TABLE(HfMcRecRedResos, "AOD", "HFMCRECREDRESO", //! Reconstruction-l
                   hf_reso_cand_reduced::InvMassGen,
                   hf_cand_mc_flag::NTracksDecayed,
                   o2::soa::Marker<1>);
+
+//! Table with selected tracks for Hc analysis
+DECLARE_SOA_TABLE(HcSelTracks, "AOD", "HCSELTRACKS",
+                  o2::soa::Index<>,
+                  // Indices
+                  hf_track_index_reduced::TrackId,
+                  hf_track_index_reduced::HfRedCollisionId,
+                  // Static
+                  hf_track_vars_reduced::Px,
+                  hf_track_vars_reduced::Py,
+                  hf_track_vars_reduced::Pz,
+                  hf_track_vars_reduced::Sign,
+                  hf_track_vars_reduced::TrackType);
 } // namespace aod
 
 namespace soa

PWGHF/D2H/TableProducer/CMakeLists.txt

@@ -90,6 +90,11 @@ o2physics_add_dpl_workflow(data-creator-jpsi-had-reduced
                     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::EventFilteringUtils
                     COMPONENT_NAME Analysis)
 
+o2physics_add_dpl_workflow(data-creator-hidden-charm-reduced
+                    SOURCES dataCreatorHiddenCharmReduced.cxx
+                    PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::EventFilteringUtils
+                    COMPONENT_NAME Analysis)
+
 # Converters
 
 o2physics_add_dpl_workflow(converter-reduced-3-prongs-ml

PWGHF/D2H/TableProducer/dataCreatorHiddenCharmReduced.cxx

@@ -0,0 +1,318 @@
+// 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 dataCreatorHiddenCharmReduced.cxx
+/// \brief Reduced data creator for hidden-charm analyses at midrapidity
+///
+/// \author A. Palasciano, <antonio.palasciano@cern.ch>, INFN Bari
+/// \author S. Politanò <stefano.politano@cern.ch>, CERN
+
+#include "PWGHF/Core/CentralityEstimation.h"
+#include "PWGHF/Core/SelectorCuts.h"
+#include "PWGHF/D2H/DataModel/ReducedDataModel.h"
+#include "PWGHF/D2H/Utils/utilsRedDataFormat.h"
+#include "PWGHF/DataModel/AliasTables.h"
+#include "PWGHF/DataModel/CandidateSelectionTables.h"
+#include "PWGHF/DataModel/TrackIndexSkimmingTables.h"
+#include "PWGHF/Utils/utilsAnalysis.h"
+#include "PWGHF/Utils/utilsBfieldCCDB.h"
+#include "PWGHF/Utils/utilsEvSelHf.h"
+
+#include "Common/Core/TPCVDriftManager.h"
+#include "Common/Core/ZorroSummary.h"
+#include "Common/DataModel/Centrality.h"
+#include "Common/DataModel/CollisionAssociationTables.h"
+#include "Common/DataModel/EventSelection.h"
+#include "Common/DataModel/PIDResponseTOF.h"
+#include "Common/DataModel/PIDResponseTPC.h"
+#include "Common/DataModel/TrackSelectionTables.h"
+
+#include <CCDB/BasicCCDBManager.h>
+#include <CCDB/CcdbApi.h>
+#include <DetectorsBase/MatLayerCylSet.h>
+#include <DetectorsBase/Propagator.h>
+#include <Framework/AnalysisDataModel.h>
+#include <Framework/AnalysisTask.h>
+#include <Framework/Configurable.h>
+#include <Framework/DeviceSpec.h>
+#include <Framework/HistogramRegistry.h>
+#include <Framework/HistogramSpec.h>
+#include <Framework/InitContext.h>
+#include <Framework/Logger.h>
+#include <Framework/runDataProcessing.h>
+
+#include <array>
+#include <cmath>
+#include <numeric>
+#include <string>
+#include <vector>
+
+using namespace o2;
+using namespace o2::analysis;
+using namespace o2::aod;
+using namespace o2::framework;
+using namespace o2::framework::expressions;
+using namespace o2::hf_centrality;
+
+enum TrackType : uint8_t {
+  Pion = 0,
+  Kaon,
+  Proton
+};
+
+struct HfDataCreatorHiddenCharmReduced {
+
+  Produces<aod::HfRedCollisions> hfReducedCollision;
+  Produces<aod::HfOrigColCounts> hfCollisionCounter;
+  Produces<aod::HcSelTracks> hfTrackLite;
+
+  struct : ConfigurableGroup {
+    // track quality
+    Configurable<bool> fillHistograms{"fillHistograms", true, "Fill proton QA histograms"};
+    Configurable<bool> selectProtons{"selectProtons", true, "Select protons"};
+    Configurable<int> itsNClsMin{"itsNClsMin", 5, "Minimum number of ITS clusters"};
+    Configurable<int> tpcNClsFoundMin{"tpcNClsFoundMin", 50, "Minimum number of found TPC clusters"};
+    Configurable<int> tpcNClsCrossedRowsMin{"tpcNClsCrossedRowsMin", 80, "Minimum number of crossed TPC rows"};
+    Configurable<double> ptMinTrack{"ptMinTrack", 0.5, "Minimum proton-track pT"};
+    Configurable<double> etaMaxTrack{"etaMaxTrack", 0.8, "Maximum proton-track |eta|"};
+    Configurable<float> momForCombinedPid{"momForCombinedPid", 0.75f, "Momentum threshold above which combined TPC+TOF proton PID is used"};
+    Configurable<std::vector<double>> binsPtTrack{"binsPtTrack", std::vector<double>{hf_cuts_single_track::vecBinsPtTrack}, "Track pT bin limits for DCA cuts"};
+    Configurable<LabeledArray<double>> cutsTrack{"cutsTrack", {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 DCA selections per pT bin"};
+    Configurable<float> maxNsigmaTofPi{"maxNsigmaTofPi", 2.f, "Maximum pion n-sigma in TOF for proton rejection"};
+    Configurable<float> maxNsigmaTofKa{"maxNsigmaTofKa", 2.f, "Maximum kaon n-sigma in TOF for proton rejection"};
+    Configurable<float> maxNsigmaCombinedPr{"maxNsigmaCombinedPr", 3.f, "Maximum combined proton n-sigma from TPC and TOF"};
+    Configurable<float> maxNsigmaTpcPi{"maxNsigmaTpcPi", 2.f, "Maximum pion n-sigma in TPC for proton rejection"};
+    Configurable<float> maxNsigmaTpcKa{"maxNsigmaTpcKa", 2.f, "Maximum kaon n-sigma in TPC for proton rejection"};
+    Configurable<float> maxNsigmaTpcPr{"maxNsigmaTpcPr", 3.f, "Maximum proton n-sigma in TPC"};
+    Configurable<bool> forceTOF{"forceTOF", false, "Require TOF PID information for proton selection"};
+  } config;
+
+  // Configurable (others)
+  Configurable<std::string> ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
+  Configurable<std::string> ccdbPathGrpMag{"ccdbPathGrpMag", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object (Run 3)"};
+  Configurable<bool> doMcRecQa{"doMcRecQa", true, "Fill QA histograms for Mc matching"};
+  Configurable<bool> rejectPairsWithCommonDaughter{"rejectPairsWithCommonDaughter", true, "flag to reject already at this stage the pairs that share a daughter track"};
+  Configurable<bool> rejectCollisionsWithBadEvSel{"rejectCollisionsWithBadEvSel", true, "flag to reject collisions with bad event selection"};
+
+  o2::hf_evsel::HfEventSelection hfEvSel;
+  o2::hf_evsel::HfEventSelectionMc hfEvSelMc;
+
+  double bz{0.};
+  int runNumber{0}; // needed to detect if the run changed and trigger update of calibrations etc.
+
+  // material correction for track propagation
+  o2::base::MatLayerCylSet* lut{};
+  o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrLUT;
+  o2::aod::common::TPCVDriftManager vDriftMgr;
+  o2::ccdb::CcdbApi ccdbApi;
+  Service<o2::ccdb::BasicCCDBManager> ccdb{};
+
+  using BCsInfo = soa::Join<aod::BCs, aod::Timestamps, aod::BcSels>;
+  using TracksWithPID = soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksDCA, aod::TrackSelection, aod::pidTPCFullPi, aod::pidTOFFullPi, aod::pidTPCFullKa, aod::pidTOFFullKa, aod::pidTPCFullPr, aod::pidTOFFullPr>;
+
+  Preslice<aod::TrackAssoc> trackIndicesPerCollision = aod::track_association::collisionId;
+
+  HistogramRegistry registry{"registry"};
+  OutputObj<ZorroSummary> zorroSummary{"zorroSummary"};
+
+  void init(InitContext& initContext)
+  {
+    // Configure CCDB access
+    ccdb->setURL(ccdbUrl.value);
+    ccdb->setCaching(true);
+    ccdb->setLocalObjectValidityChecking();
+    ccdb->setCreatedNotAfter(std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count());
+    ccdbApi.init(ccdbUrl);
+    runNumber = 0;
+    lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(ccdb->get<o2::base::MatLayerCylSet>("GLO/Param/MatLUT"));
+
+    if (config.fillHistograms) {
+      const AxisSpec axisPt{360, 0., 36., "#it{p}_{T}^{proton} (GeV/#it{c})"};
+      const AxisSpec axisEta{100, -1., 1., "#eta"};
+      const AxisSpec axisDca{400, -2., 2., "DCA_{xy} to primary vertex (cm)"};
+      const AxisSpec axisNSigma{100, -5., 5., "n#sigma"};
+
+      registry.add("hPzVtx", "Z position of primary vertex for selected tracks;z_{vtx} (cm);entries", {HistType::kTH1D, {AxisSpec{200, -20., 20., "z_{vtx} (cm)"}}});
+      registry.add("hPtNoCuts", "All associated tracks;#it{p}_{T}^{track} (GeV/#it{c});entries", {HistType::kTH1D, {axisPt}});
+      registry.add("hPtCutsProton", "Selected proton tracks;#it{p}_{T}^{track} (GeV/#it{c});entries", {HistType::kTH1D, {axisPt}});
+      registry.add("hEtaCutsProton", "Selected proton tracks;#eta;entries", {HistType::kTH1D, {axisEta}});
+      registry.add("hDCAToPrimXYVsPtCutsProton", "Selected proton tracks;#it{p}_{T}^{track} (GeV/#it{c});DCA_{xy} to primary vertex (cm)", {HistType::kTH2D, {axisPt, axisDca}});
+      registry.add("hNSigmaTPCProton", "Selected proton tracks;#it{p}_{T}^{track} (GeV/#it{c});n#sigma_{TPC}", {HistType::kTH2D, {axisPt, axisNSigma}});
+      registry.add("hNSigmaTOFProton", "Selected proton tracks;#it{p}_{T}^{track} (GeV/#it{c});n#sigma_{TOF}", {HistType::kTH2D, {axisPt, axisNSigma}});
+      registry.add("hInvMass", "Invariant mass of selected proton with all other tracks in the event;#it{p}_{T}^{proton} (GeV/#it{c});invariant mass with other tracks (GeV/#it{c}^{2})", {HistType::kTH2D, {axisPt, AxisSpec{100, 2.85, 3.25, "invariant mass with other tracks (GeV/#it{c}^{2})"}}});
+    }
+
+    // init HF event selection helper
+    hfEvSel.init(registry, &zorroSummary);
+
+    const auto& workflows = initContext.services().get<RunningWorkflowInfo const>();
+    for (const DeviceSpec& device : workflows.devices) {
+      if (device.name == "hf-data-creator-charm-reso-to-d0-reduced") {
+        // init HF event selection helper
+        hfEvSelMc.init(device, registry);
+        break;
+      }
+    }
+  }
+
+  template <typename TTrack>
+  bool isSelectedPid(TTrack const& track) const
+  {
+    const float momForCombinedPid = config.momForCombinedPid.value;
+    const float maxNsigmaTpcPr = config.maxNsigmaTpcPr.value;
+    const float maxNsigmaTpcPi = config.maxNsigmaTpcPi.value;
+    const float maxNsigmaTpcKa = config.maxNsigmaTpcKa.value;
+    const float maxNsigmaCombinedPr = config.maxNsigmaCombinedPr.value;
+    const float maxNsigmaTofPi = config.maxNsigmaTofPi.value;
+    const float maxNsigmaTofKa = config.maxNsigmaTofKa.value;
+    const bool forceTOF = config.forceTOF.value;
+
+    const float mom = std::hypot(track.px(), track.py(), track.pz());
+    const float nSigmaTPCPr = track.tpcNSigmaPr();
+    const float nSigmaTPCPi = track.tpcNSigmaPi();
+    const float nSigmaTPCKa = track.tpcNSigmaKa();
+    const bool hasTOF = track.hasTOF();
+    if (!hasTOF && forceTOF) {
+      return false;
+    }
+
+    bool isProton = false;
+    bool rejectAsPion = false;
+    bool rejectAsKaon = false;
+
+    if (mom < momForCombinedPid || !hasTOF) {
+      isProton = std::abs(nSigmaTPCPr) < maxNsigmaTpcPr;
+      rejectAsPion = std::abs(nSigmaTPCPi) < maxNsigmaTpcPi;
+      rejectAsKaon = std::abs(nSigmaTPCKa) < maxNsigmaTpcKa;
+    } else {
+      const float nSigmaTOFPr = track.tofNSigmaPr();
+      const float nSigmaTOFPi = track.tofNSigmaPi();
+      const float nSigmaTOFKa = track.tofNSigmaKa();
+      isProton = std::hypot(nSigmaTPCPr, nSigmaTOFPr) < maxNsigmaCombinedPr;
+      rejectAsPion = std::hypot(nSigmaTPCPi, nSigmaTOFPi) < maxNsigmaTofPi;
+      rejectAsKaon = std::hypot(nSigmaTPCKa, nSigmaTOFKa) < maxNsigmaTofKa;
+    }
+    return isProton && !rejectAsPion && !rejectAsKaon;
+  }
+
+  template <typename TTrack>
+  bool isSelectedTrack(TTrack const& track) const
+  {
+    const int tpcNClsFoundMin = config.tpcNClsFoundMin.value;
+    const int tpcNClsCrossedRowsMin = config.tpcNClsCrossedRowsMin.value;
+    const int itsNClsMin = config.itsNClsMin.value;
+    const double etaMaxTrack = config.etaMaxTrack.value;
+    const double ptMinTrack = config.ptMinTrack.value;
+
+    if (!track.isGlobalTrackWoDCA()) {
+      return false;
+    }
+    if (track.tpcNClsFound() < tpcNClsFoundMin) {
+      return false;
+    }
+    if (track.tpcNClsCrossedRows() < tpcNClsCrossedRowsMin) {
+      return false;
+    }
+    if (track.itsNCls() < itsNClsMin) {
+      return false;
+    }
+    if (std::abs(track.eta()) > etaMaxTrack) {
+      return false;
+    }
+    if (track.pt() < ptMinTrack) {
+      return false;
+    }
+    if (!isSelectedTrackDca(config.binsPtTrack, config.cutsTrack, track.pt(), track.dcaXY(), track.dcaZ())) {
+      return false;
+    }
+    return isSelectedPid(track);
+  }
+
+  void processEtaCTrack(soa::Join<aod::Collisions, aod::EvSels> const& collisions,
+                        aod::TrackAssoc const& trackIndices,
+                        TracksWithPID const& tracks,
+                        aod::BCsWithTimestamps const&)
+  {
+    hfTrackLite.reserve(tracks.size());
+
+    int zvtxColl{0};
+    int sel8Coll{0};
+    int zvtxAndSel8Coll{0};
+    int zvtxAndSel8CollAndSoftTrig{0};
+    int allSelColl{0};
+    constexpr int NDaughtersCharmMeson = 2;
+    for (const auto& collision : collisions) {
+      const auto hfRejMap = o2::hf_evsel::getEvSel<true, o2::hf_centrality::CentralityEstimator::None, aod::BCsWithTimestamps>(collision, hfEvSel, zvtxColl, sel8Coll, zvtxAndSel8Coll, zvtxAndSel8CollAndSoftTrig, allSelColl, ccdb, registry);
+      if (rejectCollisionsWithBadEvSel && hfRejMap != 0) {
+        continue;
+      }
+      const auto thisCollId = collision.globalIndex();
+      const auto trackIds = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId);
+      if (config.fillHistograms) {
+        registry.fill(HIST("hPzVtx"), collision.posZ());
+      }
+      std::vector<int> selectedTrackIds;
+      for (const auto& trkId : trackIds) {
+        auto trk = trkId.track_as<TracksWithPID>();
+        if (config.fillHistograms) {
+          registry.fill(HIST("hPtNoCuts"), trk.pt());
+        }
+        if (!isSelectedTrack(trk)) {
+          continue;
+        }
+        std::array pVecProton{trk.pVector()};
+        hfTrackLite(trk.globalIndex(), collision.globalIndex(), pVecProton[0], pVecProton[1], pVecProton[2], trk.sign(), static_cast<uint8_t>(TrackType::Proton));
+        selectedTrackIds.push_back(trk.globalIndex());
+        if (config.fillHistograms) {
+          registry.fill(HIST("hPtCutsProton"), trk.pt());
+          registry.fill(HIST("hEtaCutsProton"), trk.eta());
+          registry.fill(HIST("hDCAToPrimXYVsPtCutsProton"), trk.pt(), trk.dcaXY());
+          registry.fill(HIST("hNSigmaTPCProton"), trk.pt(), trk.tpcNSigmaPr());
+          if (trk.hasTOF()) {
+            registry.fill(HIST("hNSigmaTOFProton"), trk.pt(), trk.tofNSigmaPr());
+          }
+        }
+      }
+      if (selectedTrackIds.size() < NDaughtersCharmMeson) {
+        continue;
+      }
+      hfReducedCollision(collision.posX(), collision.posY(), collision.posZ(), collision.numContrib(), hfRejMap, bz);
+      for (size_t i = 0; i < selectedTrackIds.size(); ++i) {
+        auto t1 = tracks.rawIteratorAt(selectedTrackIds[i]);
+        std::array pVec1{t1.pVector()};
+        for (size_t j = i + 1; j < selectedTrackIds.size(); ++j) {
+          auto t2 = tracks.rawIteratorAt(selectedTrackIds[j]);
+          if (t1.sign() * t2.sign() > 0) {
+            continue;
+          }
+          std::array pVec2{t2.pVector()};
+          float invMass = RecoDecay::m(std::array{pVec1, pVec2}, std::array{o2::constants::physics::MassProton, o2::constants::physics::MassProton});
+          float ptEtac = RecoDecay::pt(RecoDecay::sumOfVec(pVec1, pVec2));
+          if (config.fillHistograms) {
+            registry.fill(HIST("hInvMass"), ptEtac, invMass);
+          }
+        }
+      }
+    }
+    hfCollisionCounter(collisions.tableSize(), zvtxColl, sel8Coll, zvtxAndSel8Coll, zvtxAndSel8CollAndSoftTrig, allSelColl);
+  }
+
+  PROCESS_SWITCH(HfDataCreatorHiddenCharmReduced, processEtaCTrack, "EtaC -> p pbar reconstruction", true);
+};
+
+// struct
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  WorkflowSpec workflow{};
+  workflow.push_back(adaptAnalysisTask<HfDataCreatorHiddenCharmReduced>(cfgc));
+  return workflow;
+}

PWGHF/D2H/Tasks/CMakeLists.txt

@@ -153,3 +153,8 @@ o2physics_add_dpl_workflow(task-xic0-to-xi-pi
                     SOURCES taskXic0ToXiPi.cxx
                     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
+
+o2physics_add_dpl_workflow(task-hidden-charm
+                    SOURCES taskHiddenCharm.cxx
+                    PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
+                    COMPONENT_NAME Analysis)