adding track selection to jet trigger - committing on behalf of Filip Krizek

Author: Nima Zardoshti

EventFiltering/CMakeLists.txt

@@ -56,7 +56,7 @@ o2physics_add_dpl_workflow(cf-threebodyfemto-filter
 
 o2physics_add_dpl_workflow(je-filter
                            SOURCES PWGJE/jetFilter.cxx
-                           PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore FastJet::FastJet FastJet::Contrib
+                           PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore FastJet::FastJet FastJet::Contrib
                            COMPONENT_NAME Analysis)
 
 o2physics_add_dpl_workflow(lf-strangeness-filter

EventFiltering/PWGJE/jetFilter.cxx

@@ -1,6 +1,6 @@
 // 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.
+// 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".
@@ -12,27 +12,29 @@
 
 // Author: Filip Krizek
 
-#include "ReconstructionDataFormats/Track.h"
-#include "Framework/runDataProcessing.h"
-#include "Framework/AnalysisTask.h"
-#include "Framework/AnalysisDataModel.h"
+#include "Framework/ASoA.h"
 #include "Framework/ASoAHelpers.h"
-//FK #include "Common/DataModel/PIDResponse.h"
+#include "Framework/AnalysisDataModel.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/runDataProcessing.h"
+#include "ReconstructionDataFormats/Track.h"
 
+#include "Common/Core/TrackSelection.h"
+#include "Common/Core/TrackSelectionDefaults.h"
 #include "Common/DataModel/EventSelection.h"
 #include "Common/DataModel/TrackSelectionTables.h"
 
-#include "PWGJE/DataModel/Jet.h"
-#include "PWGJE/DataModel/EMCALClusters.h"
 #include "PWGJE/Core/JetFinder.h"
+#include "PWGJE/DataModel/EMCALClusters.h"
+#include "PWGJE/DataModel/Jet.h"
 
 #include "../filterTables.h"
 
 #include "Framework/HistogramRegistry.h"
 
+#include <TMath.h>
 #include <cmath>
 #include <string>
-#include <TMath.h>
 
 using namespace o2;
 using namespace o2::framework;
@@ -44,53 +46,135 @@ struct jetFilter {
   enum { kJetChHighPt = 0,
          kHighPtObjects };
 
-  //event selection cuts
-  Configurable<float> selectionJetChHighPt{"selectionJetChHighPt", 33., "Minimum charged jet pT trigger threshold"}; //we want to keep all events having a charged jet with pT above this
+  // event selection cuts
+  Configurable<float> selectionJetChHighPt{
+    "selectionJetChHighPt", 33.,
+    "Minimum charged jet pT trigger threshold"}; // we want to keep all events
+                                                 // having a charged jet with
+                                                 // pT above this
 
   Produces<aod::JetFilters> tags;
 
-  //acceptance cuts
-  Configurable<float> cfgVertexCut{"cfgVertexCut", 10.0, "Accepted z-vertex range"};
-
-  Configurable<float> cfgJetChEtaCut{"cfgJetChEtaCut", 0.9, "Eta range for charged jets"};
-
-  HistogramRegistry spectra{"spectra", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
+  // acceptance cuts
+  Configurable<float> cfgVertexCut{"cfgVertexCut", 10.0,
+                                   "Accepted z-vertex range"};
+
+  Configurable<float> cfgTPCVolume{"cfgTPCVolume", 0.9,
+                                   "Full eta range"}; // eta range of TPC
+  Configurable<float> cfgJetR{"cfgJetR", 0.6,
+                              "jet resolution parameter"}; // jet cone radius
+  Configurable<float> cfgJetPtMin{
+    "cfgJetPtMin", 0.05,
+    "minimum jet pT constituent cut"}; // minimum jet constituent pT
+
+  HistogramRegistry spectra{
+    "spectra",
+    {},
+    OutputObjHandlingPolicy::AnalysisObject,
+    true,
+    true};
+
+  std::vector<fastjet::PseudoJet> jetConstituents;
+  std::vector<fastjet::PseudoJet> jetReclustered;
+  JetFinder jetReclusterer;
 
   void init(o2::framework::InitContext&)
   {
 
-    spectra.add("fCollZpos", "collision z position", HistType::kTH1F, {{200, -20., +20., "#it{z}_{vtx} position (cm)"}});
-    spectra.add("fJetChPtSelected", "pT of selected high pT charged jets", HistType::kTH1F, {{150, 0., +150., "charged jet #it{p}_{T} (GeV/#it{c})"}});
-
-    auto scalers{std::get<std::shared_ptr<TH1>>(spectra.add("fProcessedEvents", ";;Number of filtered events", HistType::kTH1F, {{kHighPtObjects, -0.5, kHighPtObjects - 0.5}}))};
+    spectra.add("fCollZpos", "collision z position", HistType::kTH1F,
+                {{200, -20., +20., "#it{z}_{vtx} position (cm)"}});
+    spectra.add(
+      "ptphiGoodTracks", "ptphiGoodTracks", HistType::kTH2F,
+      {{100, 0., +100., "#it{p}_{T} (GeV)"}, {60, 0, TMath::TwoPi()}});
+    spectra.add(
+      "ptphiRejectedTracks", "ptphiRejectedTracks", HistType::kTH2F,
+      {{100, 0., +100., "#it{p}_{T} (GeV)"}, {60, 0, TMath::TwoPi()}});
+    spectra.add("ptetaGoodTracks", "ptetaGoodTracks", HistType::kTH2F,
+                {{100, 0., +100., "#it{p}_{T} (GeV)"}, {36, -0.9, 0.9}});
+    spectra.add("ptetaRejectedTracks", "ptetaRejectedTracks", HistType::kTH2F,
+                {{100, 0., +100., "#it{p}_{T} (GeV)"}, {36, -0.9, 0.9}});
+
+    spectra.add("ptphiJetChSelected",
+                "pT of selected high pT charged jets vs phi", HistType::kTH2F,
+                {{150, 0., +150., "charged jet #it{p}_{T} (GeV/#it{c})"},
+                 {60, 0, TMath::TwoPi()}});
+    spectra.add("ptetaJetChSelected",
+                "pT of selected high pT charged jets vs eta", HistType::kTH2F,
+                {{150, 0., +150., "charged jet #it{p}_{T} (GeV/#it{c})"},
+                 {36, -0.9, 0.9}});
+
+    auto scalers{std::get<std::shared_ptr<TH1>>(spectra.add(
+      "fProcessedEvents", ";;Number of filtered events", HistType::kTH1F,
+      {{kHighPtObjects, -0.5, kHighPtObjects - 0.5}}))};
     for (uint32_t iS{1}; iS <= highPtObjectsNames.size(); ++iS) {
       scalers->GetXaxis()->SetBinLabel(iS, highPtObjectsNames[iS - 1].data());
     }
+
+    jetReclusterer.isReclustering = true;
+    jetReclusterer.algorithm = fastjet::JetAlgorithm::antikt_algorithm;
+    jetReclusterer.jetR = cfgJetR;
   }
 
-  //declare filters on tracks and charged jets
-  Filter collisionFilter = nabs(aod::collision::posZ) < cfgVertexCut;
-  Filter jetChFilter = (nabs(aod::jet::eta) < cfgJetChEtaCut);
+  // declare filters on tracks
+  // Filter collisionFilter = nabs(aod::collision::posZ) < cfgVertexCut;
 
-  using TrackCandidates = soa::Filtered<soa::Join<aod::Tracks, aod::TrackSelection>>;
+  using TrackCandidates = soa::Join<aod::Tracks, aod::TracksExtra,
+                                    aod::TracksDCA, aod::TrackSelection>;
 
-  void process(soa::Filtered<soa::Join<aod::Collisions, aod::EvSels>>::iterator const& collision, TrackCandidates const& tracks, aod::Jets const& jets)
+  void
+    process(soa::Join<aod::Collisions, aod::EvSels>::iterator const& collision,
+            TrackCandidates const& tracks)
   {
     // collision process loop
     bool keepEvent[kHighPtObjects]{false};
     //
     spectra.fill(HIST("fCollZpos"), collision.posZ());
 
-    //Check whether there is a high pT charged jet
-    for (auto& jet : jets) { // start loop over charged jets
-      if (jet.pt() >= selectionJetChHighPt) {
-        spectra.fill(HIST("fJetChPtSelected"), jet.pt()); //charged jet pT
-        keepEvent[kJetChHighPt] = true;
-        break;
+    jetConstituents.clear();
+    jetReclustered.clear();
+
+    for (auto& trk : tracks) {
+      if (!trk.isQualityTrack()) {
+        spectra.fill(HIST("ptphiRejectedTracks"), trk.pt(), trk.phi());
+        spectra.fill(HIST("ptetaRejectedTracks"), trk.pt(), trk.eta());
+        continue;
+      }
+      if (fabs(trk.eta()) < cfgTPCVolume) {
+        spectra.fill(HIST("ptphiGoodTracks"), trk.pt(), trk.phi());
+        spectra.fill(HIST("ptetaGoodTracks"), trk.pt(), trk.eta());
+
+        if (trk.pt() > cfgJetPtMin) { // jet constituents
+          fillConstituents(trk,
+                           jetConstituents); // ./PWGJE/Core/JetFinder.h
+                                             // recombination scheme is assumed
+                                             // to be Escheme with pion mass
+        }
+      }
+    }
+
+    // Reconstruct jet from tracks
+    fastjet::ClusterSequenceArea clusterSeq(
+      jetReclusterer.findJets(jetConstituents, jetReclustered));
+    jetReclustered = sorted_by_pt(jetReclustered);
+
+    // Check whether there is a high pT charged jet
+    for (auto& jet : jetReclustered) { // start loop over charged jets
+      if (fabs(jet.eta()) < cfgTPCVolume) {
+        if (jet.pt() >= selectionJetChHighPt) {
+          spectra.fill(HIST("ptphiJetChSelected"), jet.pt(),
+                       jet.phi()); // charged jet pT vs phi
+          spectra.fill(HIST("ptetaJetChSelected"), jet.pt(),
+                       jet.eta()); // charged jet pT vs eta
+          keepEvent[kJetChHighPt] = true;
+          break;
+        }
       }
     }
 
-    //count events which passed the selections
+    // count events which passed the selections
+    if (fabs(collision.posZ()) > cfgVertexCut)
+      keepEvent[kJetChHighPt] = false;
+
     for (int iDecision{0}; iDecision < kHighPtObjects; ++iDecision) {
       if (keepEvent[iDecision]) {
         spectra.fill(HIST("fProcessedEvents"), iDecision);
@@ -103,6 +187,5 @@ struct jetFilter {
 WorkflowSpec defineDataProcessing(ConfigContext const& cfg)
 {
 
-  return WorkflowSpec{
-    adaptAnalysisTask<jetFilter>(cfg)};
+  return WorkflowSpec{adaptAnalysisTask<jetFilter>(cfg)};
 }