CF 2-prong candidate filter and tables (AliceO2Group#4362)

PWGCF/DataModel/CorrelationsDerived.h

@@ -8,8 +8,8 @@
 // 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.
-#ifndef O2_ANALYSIS_CFDERIVED_H
-#define O2_ANALYSIS_CFDERIVED_H
+#ifndef PWGCF_DATAMODEL_CORRELATIONSDERIVED_H_
+#define PWGCF_DATAMODEL_CORRELATIONSDERIVED_H_
 
 #include "Framework/ASoA.h"
 #include "Framework/AnalysisDataModel.h"
@@ -79,6 +79,42 @@ using CFTrack = CFTracks::iterator;
 using CFTrackLabel = CFTrackLabels::iterator;
 using CFTracksWithLabel = soa::Join<CFTracks, CFTrackLabels>;
 using CFTrackWithLabel = CFTracksWithLabel::iterator;
+
+//------transient CF-filter to CF-2prong-filter
+DECLARE_SOA_TABLE(CFCollRefs, "AOD", "CFCOLLREF", o2::soa::Index<>, track::CollisionId); //! Transient cf collision index table
+
+using CFCollRef = CFCollRefs::iterator;
+
+namespace cftrackref
+{
+DECLARE_SOA_INDEX_COLUMN(Track, track);
+} // namespace cftrackref
+DECLARE_SOA_TABLE(CFTrackRefs, "AOD", "CFTRACKREF", o2::soa::Index<>, track::CollisionId, cftrackref::TrackId); //! Transient cf track index table
+
+using CFTrackRef = CFTrackRefs::iterator;
+//------
+
+namespace cf2prongtrack
+{
+DECLARE_SOA_INDEX_COLUMN_FULL(CFTrackProng0, cfTrackProng0, int, CFTracks, "_0"); //! Index to prong 1 CFTrack
+DECLARE_SOA_INDEX_COLUMN_FULL(CFTrackProng1, cfTrackProng1, int, CFTracks, "_1"); //! Index to prong 2 CFTrack
+DECLARE_SOA_COLUMN(Pt, pt, float);                                                //! pT (GeV/c)
+DECLARE_SOA_COLUMN(Eta, eta, float);                                              //! Pseudorapidity
+DECLARE_SOA_COLUMN(Phi, phi, float);                                              //! Phi angle
+DECLARE_SOA_COLUMN(Decay, decay, uint8_t);                                        //! Particle decay
+enum ParticleDecay {
+  D0ToPiK,
+  JPsiToEE,
+  JPsiToMuMu
+};
+} // namespace cf2prongtrack
+DECLARE_SOA_TABLE(CF2ProngTracks, "AOD", "CF2PRONGTRACK", //! Reduced track table
+                  o2::soa::Index<>,
+                  cftrack::CFCollisionId,
+                  cf2prongtrack::CFTrackProng0Id,
+                  cf2prongtrack::CFTrackProng1Id,
+                  cf2prongtrack::Pt, cf2prongtrack::Eta, cf2prongtrack::Phi, cf2prongtrack::Decay);
+using CF2ProngTrack = CF2ProngTracks::iterator;
 } // namespace o2::aod
 
-#endif // O2_ANALYSIS_CFDERIVED_H
+#endif // PWGCF_DATAMODEL_CORRELATIONSDERIVED_H_

PWGCF/TableProducer/CMakeLists.txt

@@ -14,7 +14,12 @@ o2physics_add_dpl_workflow(filter-correlations
                     PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGCFCore
                     COMPONENT_NAME Analysis)
 
+o2physics_add_dpl_workflow(filter-correlations-2prong
+                    SOURCES filter2Prong.cxx
+                    PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGCFCore
+                    COMPONENT_NAME Analysis)
+
 o2physics_add_dpl_workflow(dptdpt-filter
                            SOURCES dptdptfilter.cxx
                            PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGCFCore
-                           COMPONENT_NAME Analysis)
+                           COMPONENT_NAME Analysis)

PWGCF/TableProducer/filter2Prong.cxx

@@ -0,0 +1,73 @@
+// 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.
+#include "Framework/runDataProcessing.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/AnalysisDataModel.h"
+#include "Framework/ASoAHelpers.h"
+
+#include "MathUtils/detail/TypeTruncation.h"
+
+#include "PWGCF/DataModel/CorrelationsDerived.h"
+
+#include "PWGHF/Core/HfHelper.h"
+#include "PWGHF/DataModel/CandidateReconstructionTables.h"
+#include "PWGHF/DataModel/CandidateSelectionTables.h"
+
+using namespace o2;
+using namespace o2::framework;
+using namespace o2::framework::expressions;
+using namespace o2::math_utils::detail;
+
+#define FLOAT_PRECISION 0xFFFFFFF0
+#define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable<TYPE> NAME{#NAME, DEFAULT, HELP};
+
+struct Filter2Prong {
+  O2_DEFINE_CONFIGURABLE(cfgVerbosity, int, 0, "Verbosity level (0 = major, 1 = per collision)")
+
+  Produces<aod::CF2ProngTracks> output2ProngTracks;
+
+  using HFCandidates = soa::Join<aod::HfCand2Prong, aod::HfSelD0>;
+  void processData(aod::Collisions::iterator const& collision, aod::BCsWithTimestamps const&, aod::CFCollRefs const& cfcollisions, aod::CFTrackRefs const& cftracks, HFCandidates const& candidates)
+  {
+    if (cfcollisions.size() <= 0 || cftracks.size() <= 0)
+      return; // rejected collision
+    if (cfgVerbosity > 0 && candidates.size() > 0)
+      LOGF(info, "Candidates for collision: %lu, cfcollisions: %lu, CFTracks: %lu", candidates.size(), cfcollisions.size(), cftracks.size());
+    for (auto& c : candidates) {
+      int prongCFId[2] = {-1, -1};
+      for (auto& cftrack : cftracks) {
+        if (c.prong0Id() == cftrack.trackId()) {
+          prongCFId[0] = cftrack.globalIndex();
+          break;
+        }
+      }
+      for (auto& cftrack : cftracks) {
+        if (c.prong1Id() == cftrack.trackId()) {
+          prongCFId[1] = cftrack.globalIndex();
+          break;
+        }
+      }
+      // look-up the collision id
+      auto collisionId = cfcollisions.begin().globalIndex();
+      if ((c.hfflag() & (1 << aod::hf_cand_2prong::DecayType::D0ToPiK)) == 0)
+        continue;
+      output2ProngTracks(collisionId,
+                         prongCFId[0], prongCFId[1], c.pt(), c.eta(), c.phi(), aod::cf2prongtrack::D0ToPiK);
+    }
+  }
+  PROCESS_SWITCH(Filter2Prong, processData, "Process data D0 candidates", true);
+}; // struct
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  return WorkflowSpec{
+    adaptAnalysisTask<Filter2Prong>(cfgc, TaskName{"filter-2prong"})};
+}

PWGCF/TableProducer/filterCorrelations.cxx

@@ -55,6 +55,7 @@ struct FilterCF {
   O2_DEFINE_CONFIGURABLE(cfgVerbosity, int, 1, "Verbosity level (0 = major, 1 = per collision)")
   O2_DEFINE_CONFIGURABLE(cfgTrigger, int, 7, "Trigger choice: (0 = none, 7 = sel7, 8 = sel8)")
   O2_DEFINE_CONFIGURABLE(cfgCollisionFlags, uint16_t, aod::collision::CollisionFlagsRun2::Run2VertexerTracks, "Request collision flags if non-zero (0 = off, 1 = Run2VertexerTracks)")
+  O2_DEFINE_CONFIGURABLE(cfgTransientTables, bool, false, "Output transient tables for collision and track IDs")
 
   // Filters and input definitions
   Filter collisionZVtxFilter = nabs(aod::collision::posZ) < cfgCutVertex;
@@ -78,6 +79,9 @@ struct FilterCF {
   Produces<aod::CFMcCollisions> outputMcCollisions;
   Produces<aod::CFMcParticles> outputMcParticles;
 
+  Produces<aod::CFCollRefs> outputCollRefs;
+  Produces<aod::CFTrackRefs> outputTrackRefs;
+
   template <typename TCollision>
   bool keepCollision(TCollision& collision)
   {
@@ -104,6 +108,9 @@ struct FilterCF {
     auto bc = collision.bc_as<aod::BCsWithTimestamps>();
     outputCollisions(bc.runNumber(), collision.posZ(), collision.multiplicity(), bc.timestamp());
 
+    if (cfgTransientTables)
+      outputCollRefs(collision.globalIndex());
+
     for (auto& track : tracks) {
       uint8_t trackType = 0;
       if (track.isGlobalTrack()) {
@@ -113,6 +120,8 @@ struct FilterCF {
       }
 
       outputTracks(outputCollisions.lastIndex(), track.pt(), track.eta(), track.phi(), track.sign(), trackType);
+      if (cfgTransientTables)
+        outputTrackRefs(collision.globalIndex(), track.globalIndex());
 
       yields->Fill(collision.multiplicity(), track.pt(), track.eta());
       etaphi->Fill(collision.multiplicity(), track.eta(), track.phi());

PWGCF/Tasks/correlations.cxx

@@ -8,11 +8,18 @@
 // 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.
+#include <experimental/type_traits>
+
+#include <TH1F.h>
+#include <cmath>
+#include <TDirectory.h>
+#include <THn.h>
+
 #include "Framework/runDataProcessing.h"
 #include "Framework/AnalysisTask.h"
 #include "Framework/AnalysisDataModel.h"
 #include "Framework/ASoAHelpers.h"
-#include <CCDB/BasicCCDBManager.h>
+#include "CCDB/BasicCCDBManager.h"
 #include "Framework/StepTHn.h"
 #include "Framework/HistogramRegistry.h"
 #include "Framework/RunningWorkflowInfo.h"
@@ -27,11 +34,6 @@
 #include "DataFormatsParameters/GRPObject.h"
 #include "DataFormatsParameters/GRPMagField.h"
 
-#include <TH1F.h>
-#include <cmath>
-#include <TDirectory.h>
-#include <THn.h>
-
 using namespace o2;
 using namespace o2::framework;
 using namespace o2::framework::expressions;
@@ -78,6 +80,8 @@ struct CorrelationTask {
 
   O2_DEFINE_CONFIGURABLE(cfgVerbosity, int, 1, "Verbosity level (0 = major, 1 = per collision)")
 
+  O2_DEFINE_CONFIGURABLE(cfgDecayParticleMask, int, 0, "Selection bitmask for the decay particles: 0 = no selection")
+
   ConfigurableAxis axisVertex{"axisVertex", {7, -7, 7}, "vertex axis for histograms"};
   ConfigurableAxis axisDeltaPhi{"axisDeltaPhi", {72, -PIHalf, PIHalf * 3}, "delta phi axis for histograms"};
   ConfigurableAxis axisDeltaEta{"axisDeltaEta", {40, -2, 2}, "delta eta axis for histograms"};
@@ -102,6 +106,9 @@ struct CorrelationTask {
   Filter cfMCCollisionFilter = nabs(aod::mccollision::posZ) < cfgCutVertex;
   Filter cfMCParticleFilter = (nabs(aod::cfmcparticle::eta) < cfgCutEta) && (aod::cfmcparticle::pt > cfgCutPt) && (aod::cfmcparticle::sign != 0);
 
+  // HF filters
+  Filter track2pFilter = (aod::cf2prongtrack::eta < cfgCutEta) && (aod::cf2prongtrack::pt > cfgCutPt);
+
   // Output definitions
   OutputObj<CorrelationContainer> same{"sameEvent"};
   OutputObj<CorrelationContainer> mixed{"mixedEvent"};
@@ -128,6 +135,10 @@ struct CorrelationTask {
   {
     registry.add("yields", "multiplicity/centrality vs pT vs eta", {HistType::kTH3F, {{100, 0, 100, "/multiplicity/centrality"}, {40, 0, 20, "p_{T}"}, {100, -2, 2, "#eta"}}});
     registry.add("etaphi", "multiplicity/centrality vs eta vs phi", {HistType::kTH3F, {{100, 0, 100, "multiplicity/centrality"}, {100, -2, 2, "#eta"}, {200, 0, 2 * M_PI, "#varphi"}}});
+    if (doprocessSame2ProngDerived) {
+      registry.add("yieldsTrack2", "multiplicity/centrality vs pT vs eta (track2)", {HistType::kTH3F, {{100, 0, 100, "/multiplicity/centrality"}, {40, 0, 20, "p_{T}"}, {100, -2, 2, "#eta"}}});
+      registry.add("etaphiTrack2", "multiplicity/centrality vs eta vs phi (track2)", {HistType::kTH3F, {{100, 0, 100, "multiplicity/centrality"}, {100, -2, 2, "#eta"}, {200, 0, 2 * M_PI, "#varphi"}}});
+    }
 
     const int maxMixBin = AxisSpec(axisMultiplicity).getNbins() * AxisSpec(axisVertex).getNbins();
     registry.add("eventcount_same", "bin", {HistType::kTH1F, {{maxMixBin + 2, -2.5, -0.5 + maxMixBin, "bin"}}});
@@ -171,7 +182,7 @@ struct CorrelationTask {
     ccdb->setCaching(true);
     ccdb->setLocalObjectValidityChecking();
 
-    long now = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
+    auto now = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
     ccdb->setCreatedNotAfter(now); // TODO must become global parameter from the train creation time
   }
 
@@ -193,14 +204,28 @@ struct CorrelationTask {
   }
 
   template <typename TCollision, typename TTracks>
-  void fillQA(TCollision collision, float multiplicity, TTracks tracks)
+  void fillQA(const TCollision& collision, float multiplicity, const TTracks& tracks)
   {
     for (auto& track1 : tracks) {
       registry.fill(HIST("yields"), multiplicity, track1.pt(), track1.eta());
       registry.fill(HIST("etaphi"), multiplicity, track1.eta(), track1.phi());
     }
   }
 
+  template <typename TCollision, typename TTracks1, typename TTracks2>
+  void fillQA(const TCollision& collision, float multiplicity, const TTracks1& tracks1, const TTracks2& tracks2)
+  {
+    fillQA(collision, multiplicity, tracks1);
+    for (auto& track2 : tracks2) {
+      if constexpr (std::is_same<TTracks2, aod::CF2ProngTracks>::value) {
+        if (cfgDecayParticleMask != 0 && (cfgDecayParticleMask & (1u << (uint32_t)track2.decay())) == 0u)
+          continue;
+      }
+      registry.fill(HIST("yieldsTrack2"), multiplicity, track2.pt(), track2.eta());
+      registry.fill(HIST("etaphiTrack2"), multiplicity, track2.eta(), track2.phi());
+    }
+  }
+
   template <typename TTarget, typename TCollision>
   bool fillCollisionAOD(TTarget target, TCollision collision, float multiplicity)
   {
@@ -229,8 +254,11 @@ struct CorrelationTask {
     return true;
   }
 
-  template <CorrelationContainer::CFStep step, typename TTarget, typename TTracks>
-  void fillCorrelations(TTarget target, TTracks& tracks1, TTracks& tracks2, float multiplicity, float posZ, int magField, float eventWeight)
+  template <class T>
+  using hasSign = decltype(std::declval<T&>().sign());
+
+  template <CorrelationContainer::CFStep step, typename TTarget, typename TTracks1, typename TTracks2>
+  void fillCorrelations(TTarget target, TTracks1& tracks1, TTracks2& tracks2, float multiplicity, float posZ, int magField, float eventWeight)
   {
     // Cache efficiency for particles (too many FindBin lookups)
     float* efficiencyAssociated = nullptr;
@@ -267,9 +295,17 @@ struct CorrelationTask {
       target->getTriggerHist()->Fill(step, track1.pt(), multiplicity, posZ, triggerWeight);
 
       for (auto& track2 : tracks2) {
-        if (track1.globalIndex() == track2.globalIndex()) {
-          // LOGF(info, "Track identical: %f | %f | %f || %f | %f | %f", track1.eta(), track1.phi(), track1.pt(),  track2.eta(), track2.phi(), track2.pt());
-          continue;
+        if constexpr (std::is_same<TTracks1, TTracks2>::value) {
+          if (track1.globalIndex() == track2.globalIndex()) {
+            // LOGF(info, "Track identical: %f | %f | %f || %f | %f | %f", track1.eta(), track1.phi(), track1.pt(),  track2.eta(), track2.phi(), track2.pt());
+            continue;
+          }
+        }
+        if constexpr (std::is_same<TTracks2, aod::CF2ProngTracks>::value) {
+          if (track1.globalIndex() == track2.cfTrackProng0Id() || track1.globalIndex() == track2.cfTrackProng1Id()) // do not correlate daughter tracks of the same event
+            continue;
+          if (cfgDecayParticleMask != 0 && (cfgDecayParticleMask & (1u << (uint32_t)track2.decay())) == 0u)
+            continue;
         }
 
         if constexpr (step <= CorrelationContainer::kCFStepTracked) {
@@ -282,23 +318,27 @@ struct CorrelationTask {
           continue;
         }
 
-        if (cfgAssociatedCharge != 0 && cfgAssociatedCharge * track2.sign() < 0) {
-          continue;
-        }
-        if (cfgPairCharge != 0 && cfgPairCharge * track1.sign() * track2.sign() < 0) {
-          continue;
-        }
-
-        if constexpr (step >= CorrelationContainer::kCFStepReconstructed) {
-          if (cfg.mPairCuts && mPairCuts.conversionCuts(track1, track2)) {
+        if constexpr (std::experimental::is_detected<hasSign, TTracks2>::value) {
+          if (cfgAssociatedCharge != 0 && cfgAssociatedCharge * track2.sign() < 0) {
             continue;
           }
-
-          if (cfgTwoTrackCut > 0 && mPairCuts.twoTrackCut(track1, track2, magField)) {
+          if (cfgPairCharge != 0 && cfgPairCharge * track1.sign() * track2.sign() < 0) {
             continue;
           }
         }
 
+        if constexpr (std::is_same<TTracks1, TTracks2>::value) {
+          if constexpr (step >= CorrelationContainer::kCFStepReconstructed) {
+            if (cfg.mPairCuts && mPairCuts.conversionCuts(track1, track2)) {
+              continue;
+            }
+
+            if (cfgTwoTrackCut > 0 && mPairCuts.twoTrackCut(track1, track2, magField)) {
+              continue;
+            }
+          }
+        }
+
         float associatedWeight = triggerWeight;
         if constexpr (step == CorrelationContainer::kCFStepCorrected) {
           if (cfg.mEfficiencyAssociated) {
@@ -405,6 +445,29 @@ struct CorrelationTask {
   }
   PROCESS_SWITCH(CorrelationTask, processSameDerived, "Process same event on derived data", false);
 
+  void processSame2ProngDerived(derivedCollisions::iterator const& collision, soa::Filtered<aod::CFTracks> const& tracks, soa::Filtered<aod::CF2ProngTracks> const& p2tracks)
+  {
+    if (cfgVerbosity > 0) {
+      LOGF(info, "processSame2ProngDerived: Tracks for collision: %d | 2-prong candidates: %d | Vertex: %.1f | Multiplicity/Centrality: %.1f", tracks.size(), p2tracks.size(), collision.posZ(), collision.multiplicity());
+    }
+    loadEfficiency(collision.timestamp());
+
+    const auto multiplicity = collision.multiplicity();
+
+    int bin = configurableBinningDerived.getBin({collision.posZ(), collision.multiplicity()});
+    registry.fill(HIST("eventcount_same"), bin);
+    fillQA(collision, multiplicity, tracks, p2tracks);
+
+    same->fillEvent(multiplicity, CorrelationContainer::kCFStepReconstructed);
+    fillCorrelations<CorrelationContainer::kCFStepReconstructed>(same, tracks, p2tracks, multiplicity, collision.posZ(), 0, 1.0f);
+
+    if (cfg.mEfficiencyAssociated || cfg.mEfficiencyTrigger) {
+      same->fillEvent(multiplicity, CorrelationContainer::kCFStepCorrected);
+      fillCorrelations<CorrelationContainer::kCFStepCorrected>(same, tracks, p2tracks, multiplicity, collision.posZ(), 0, 1.0f);
+    }
+  }
+  PROCESS_SWITCH(CorrelationTask, processSame2ProngDerived, "Process same event on derived data", false);
+
   using BinningTypeAOD = ColumnBinningPolicy<aod::collision::PosZ, aod::cent::CentRun2V0M>;
   void processMixedAOD(aodCollisions& collisions, aodTracks const& tracks, aod::BCsWithTimestamps const&)
   {