[PWGJE,PWGLF,Trigger] Improvements to the JE framework

EventFiltering/PWGJE/jetFilter.cxx

@@ -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.
 //
@@ -11,32 +11,33 @@
 
 // Author: Filip Krizek
 
-#include <TMath.h>
-#include <cmath>
-#include <string>
-#include <vector>
+#include "../filterTables.h"
 
-#include "Framework/ASoA.h"
-#include "Framework/ASoAHelpers.h"
-#include "Framework/AnalysisDataModel.h"
-#include "Framework/AnalysisTask.h"
-#include "Framework/runDataProcessing.h"
-#include "ReconstructionDataFormats/Track.h"
+#include "PWGJE/Core/FastJetUtilities.h"
+#include "PWGJE/Core/JetBkgSubUtils.h"
+#include "PWGJE/Core/JetDerivedDataUtilities.h"
+#include "PWGJE/Core/JetFinder.h"
+#include "PWGJE/DataModel/EMCALClusters.h"
+#include "PWGJE/DataModel/Jet.h"
 
 #include "Common/Core/TrackSelection.h"
 #include "Common/Core/TrackSelectionDefaults.h"
 #include "Common/DataModel/EventSelection.h"
 #include "Common/DataModel/TrackSelectionTables.h"
 
-#include "PWGJE/Core/JetFinder.h"
-#include "PWGJE/Core/FastJetUtilities.h"
-#include "PWGJE/Core/JetBkgSubUtils.h"
-#include "PWGJE/DataModel/EMCALClusters.h"
-#include "PWGJE/DataModel/Jet.h"
+#include "Framework/ASoA.h"
+#include "Framework/ASoAHelpers.h"
+#include "Framework/AnalysisDataModel.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/HistogramRegistry.h"
+#include "Framework/runDataProcessing.h"
+#include "ReconstructionDataFormats/Track.h"
 
-#include "../filterTables.h"
+#include <TMath.h>
 
-#include "Framework/HistogramRegistry.h"
+#include <cmath>
+#include <string>
+#include <vector>
 
 using namespace o2;
 using namespace o2::framework;
@@ -66,8 +67,8 @@
 
   Configurable<std::string> trackSelections{"trackSelections", "globalTracks", "set track selections"};
   Configurable<float> trackPtMin{"trackPtMin", 0.15, "min track pT in filter"};
   Configurable<float> trackLowPtTriggerThreshold{"trackLowPtThreshold", 10.0, "low pT track trigger threshold"};
   Configurable<float> trackHighPtTriggerThreshold{"trackHighPtThreshold", 25.0, "high pT track trigger threshold"};
   Configurable<float> jetPtLowThreshold{"jetPtLowThreshold", 30.0, "threshold for low pT jet trigger"};
   Configurable<float> jetPtHighThreshold{"jetPtHighThreshold", 50.0, "threshold for high pT jet trigger"};
 
@@ -189,7 +190,7 @@
 
     // collision process loop
     bool keepEvent[kTriggerObjects]{false, false, false, false};
-    if (!collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) {
+    if (!jetderiveddatautilities::selectCollision(collision, jetderiveddatautilities::initialiseEventSelectionBits(static_cast<std::string>("NoTimeFrameBorder")))) {
       tags(keepEvent[kJetChLowPt], keepEvent[kJetChHighPt], keepEvent[kTrackLowPt], keepEvent[kTrackHighPt]);
       return;
     }
@@ -218,7 +219,7 @@
                 spectra.fill(HIST("hPtAKTJetsInclusiveBgSubtr"), jetr, jet.pt() - (collision.rho() * jet.area()));
               }
 
               if (jet.pt() > 10.) {
                 spectra.fill(HIST("hEtaAKTJetsInclusive"), jetr, jet.eta());
                 spectra.fill(HIST("hPhiAKTJetsInclusive"), jetr, jet.phi());
               }
@@ -287,13 +288,13 @@
     tags(keepEvent[kJetChLowPt], keepEvent[kJetChHighPt], keepEvent[kTrackLowPt], keepEvent[kTrackHighPt]);
   }
 
-  void processWithoutRho(soa::Join<aod::JetCollisions, aod::EvSels>::iterator const& collision, o2::aod::ChargedJets const& jets, soa::Filtered<aod::JetTracks> const& tracks)
+  void processWithoutRho(aod::JetCollision const& collision, o2::aod::ChargedJets const& jets, soa::Filtered<aod::JetTracks> const& tracks)
   {
     doTriggering<false>(collision, jets, tracks);
   }
   PROCESS_SWITCH(jetFilter, processWithoutRho, "Do charged jet triggering without background estimation for filling histograms", true);
 
-  void processWithRho(soa::Join<aod::JetCollisions, aod::BkgChargedRhos, aod::EvSels>::iterator const& collision, o2::aod::ChargedJets const& jets, soa::Filtered<aod::JetTracks> const& tracks)
+  void processWithRho(soa::Join<aod::JetCollisions, aod::BkgChargedRhos>::iterator const& collision, o2::aod::ChargedJets const& jets, soa::Filtered<aod::JetTracks> const& tracks)
   {
     doTriggering<true>(collision, jets, tracks);
   }

PWGJE/Core/JetBkgSubUtils.cxx

@@ -48,7 +48,7 @@
 
 void JetBkgSubUtils::initialise()
 {
-  // Note: if you are using the PerpCone method you should jetBkgR to be the same as the anit_kt jets R, otherwise use R=0.2
+  // Note: recommended to use R=0.2
   jetDefBkg = fastjet::JetDefinition(algorithmBkg, jetBkgR, recombSchemeBkg, fastjet::Best);
   areaDefBkg = fastjet::AreaDefinition(fastjet::active_area_explicit_ghosts, ghostAreaSpec);
   selRho = fastjet::SelectorRapRange(bkgEtaMin, bkgEtaMax) && fastjet::SelectorPhiRange(bkgPhiMin, bkgPhiMax) && !fastjet::SelectorNHardest(nHardReject); // here we have to put rap range, to be checked!
@@ -73,7 +73,7 @@
   std::vector<double> rhoMdvector;
 
   // Fill a vector for pT/area to be used for the median
   for (auto& ijet : alljets) {
 
     // Physical area/ Physical jets (no ghost)
     if (!clusterSeq.is_pure_ghost(ijet)) {
@@ -104,62 +104,6 @@
   return std::make_tuple(rho, rhoM);
 }
 
-std::tuple<double, double> JetBkgSubUtils::estimateRhoPerpCone(const std::vector<fastjet::PseudoJet>& inputParticles, const std::vector<fastjet::PseudoJet>& jets)
-{
-
-  JetBkgSubUtils::initialise();
-  if (inputParticles.size() == 0 || jets.size() == 0) {
-    return std::make_tuple(0.0, 0.0);
-  }
-
-  double perpPtDensity1 = 0;
-  double perpPtDensity2 = 0;
-  double perpMdDensity1 = 0;
-  double perpMdDensity2 = 0;
-
-  fastjet::Selector selectJet = fastjet::SelectorEtaRange(bkgEtaMin, bkgEtaMax) && fastjet::SelectorPhiRange(bkgPhiMin, bkgPhiMax);
-
-  std::vector<fastjet::PseudoJet> selectedJets = fastjet::sorted_by_pt(selectJet(jets));
-
-  if (selectedJets.size() == 0) {
-    return std::make_tuple(0.0, 0.0);
-  }
-
-  fastjet::PseudoJet leadingJet = selectedJets[0];
-
-  double dPhi1 = 999.;
-  double dPhi2 = 999.;
-  double dEta = 999.;
-  double PerpendicularConeAxisPhi1 = 999., PerpendicularConeAxisPhi2 = 999.;
-  // build 2 perp cones in phi around the leading jet (right and left of the jet)
-  PerpendicularConeAxisPhi1 = RecoDecay::constrainAngle<double, double>(leadingJet.phi() + (M_PI / 2.)); // This will contrain the angel between 0-2Pi
-  PerpendicularConeAxisPhi2 = RecoDecay::constrainAngle<double, double>(leadingJet.phi() - (M_PI / 2.)); // This will contrain the angel between 0-2Pi
-
-  for (auto& particle : inputParticles) {
-    // sum the momentum of all paricles that fill the two cones
-    dPhi1 = particle.phi() - PerpendicularConeAxisPhi1;
-    dPhi1 = RecoDecay::constrainAngle<double, double>(dPhi1, -M_PI); // This will contrain the angel between -pi & Pi
-    dPhi2 = particle.phi() - PerpendicularConeAxisPhi2;
-    dPhi2 = RecoDecay::constrainAngle<double, double>(dPhi2, -M_PI); // This will contrain the angel between -pi & Pi
-    dEta = leadingJet.eta() - particle.eta();                        // The perp cone eta is the same as the leading jet since the cones are perpendicular only in phi
-    if (TMath::Sqrt(dPhi1 * dPhi1 + dEta * dEta) <= jetBkgR) {
-      perpPtDensity1 += particle.perp();
-      perpMdDensity1 += TMath::Sqrt(particle.m() * particle.m() + particle.pt() * particle.pt()) - particle.pt();
-    }
-
-    if (TMath::Sqrt(dPhi2 * dPhi2 + dEta * dEta) <= jetBkgR) {
-      perpPtDensity2 += particle.perp();
-      perpMdDensity2 += TMath::Sqrt(particle.m() * particle.m() + particle.pt() * particle.pt()) - particle.pt();
-    }
-  }
-
-  // Caculate rho as the ratio of average pT of the two cones / the cone area
-  double perpPtDensity = (perpPtDensity1 + perpPtDensity2) / (2 * M_PI * jetBkgR * jetBkgR);
-  double perpMdDensity = (perpMdDensity1 + perpMdDensity2) / (2 * M_PI * jetBkgR * jetBkgR);
-
-  return std::make_tuple(perpPtDensity, perpMdDensity);
-}
-
 fastjet::PseudoJet JetBkgSubUtils::doRhoAreaSub(const fastjet::PseudoJet& jet, double rhoParam, double rhoMParam)
 {
 
@@ -219,7 +163,7 @@
 {
   // Refere to https://arxiv.org/abs/1211.2811 for the rhoM caclulation
   double sum(0);
   for (auto constituent : jet.constituents()) {
     sum += TMath::Sqrt(constituent.m() * constituent.m() + constituent.pt() * constituent.pt()) - constituent.pt();
   }
 

PWGJE/Core/JetBkgSubUtils.h

@@ -30,8 +30,8 @@
 
 enum class BkgSubEstimator { none = 0,
                              medianRho = 1,
-                             medianRhoSparse = 2,
-                             perpCone = 3
+                             medianRhoSparse = 2
+                             // perpendicular cone method is in JetUtilities
 };
 
 enum class BkgSubMode { none = 0,
@@ -68,12 +68,6 @@ class JetBkgSubUtils
   /// @return Rho, RhoM the underlying event density
   std::tuple<double, double> estimateRhoAreaMedian(const std::vector<fastjet::PseudoJet>& inputParticles, bool doSparseSub);
 
-  /// @brief Background estimator using the perpendicular cone method
-  /// @param inputParticles
-  /// @param jets (all jets in the event)
-  /// @return Rho, RhoM the underlying event density
-  std::tuple<double, double> estimateRhoPerpCone(const std::vector<fastjet::PseudoJet>& inputParticles, const std::vector<fastjet::PseudoJet>& jets);
-
   /// @brief method that subtracts the background from jets using the area method
   /// @param jet input jet to be background subtracted
   /// @param rhoParam the underlying evvent density vs pT (to be set)

PWGJE/Core/JetCandidateUtilities.h

@@ -118,17 +118,16 @@
  *
  * @param track track that is being checked
  * @param candidate candidate that is being checked
- * @param tracks the track table
  */
-template <typename T, typename U, typename V>
-bool isDaughterTrack(T& track, U& candidate, V const& tracks)
+template <typename T, typename U>
+bool isDaughterTrack(T& track, U& candidate)
 {
   if constexpr (jethfutilities::isHFCandidate<U>()) {
-    return jethfutilities::isHFDaughterTrack(track, candidate, tracks);
+    return jethfutilities::isHFDaughterTrack(track, candidate);
   } else if constexpr (jetv0utilities::isV0Candidate<U>()) {
-    return jetv0utilities::isV0DaughterTrack(track, candidate, tracks);
+    return jetv0utilities::isV0DaughterTrack(track, candidate);
   } else if constexpr (jetdqutilities::isDielectronCandidate<U>()) {
-    return jetdqutilities::isDielectronDaughterTrack(track, candidate, tracks);
+    return jetdqutilities::isDielectronDaughterTrack(track, candidate);
   } else {
     return false;
   }
@@ -146,7 +145,7 @@
   if (!particle.has_daughters()) {
     return false;
   }
   for (auto daughter : particle.template daughters_as<typename std::decay_t<T>::parent_t>()) {
     if (daughter.globalIndex() == globalIndex) {
       return true;
     }

PWGJE/Core/JetDQUtilities.h

@@ -94,10 +94,9 @@ constexpr bool isMatchedDielectronCandidate(T const& /*candidate*/)
  *
  * @param track track that is being checked
  * @param candidate Dielectron candidate that is being checked
- * @param tracks the track table
  */
-template <typename T, typename U, typename V>
-bool isDielectronDaughterTrack(T& track, U& candidate, V const& /*tracks*/)
+template <typename T, typename U>
+bool isDielectronDaughterTrack(T& track, U& candidate)
 {
   if constexpr (isDielectronCandidate<U>()) {
     if (candidate.prong0Id() == track.globalIndex() || candidate.prong1Id() == track.globalIndex()) {
@@ -283,7 +282,7 @@ bool selectDielectronParticleDecay(T const& dielectronParticle, int dielectronPa
   return (dielectronParticle.decayFlag() & (1 << dielectronParticleDecaySelection));
 }
 
-int initialiseDielectronParticleDecaySelection(std::string dielectronParticleDecaySelection)
+int initialiseDielectronParticleDecaySelection(const std::string& dielectronParticleDecaySelection)
 {
   if (dielectronParticleDecaySelection == "JPsiToEE") {
     return JDielectronParticleDecays::JPsiToEE;
@@ -353,7 +352,7 @@ void fillDielectronCandidateTable(T const& candidate, int32_t collisionIndex, U&
 template <typename T, typename U>
 void fillDielectronCandidateMcTable(T const& candidate, int32_t mcCollisionIndex, U& DielectronMcTable)
 {
-  DielectronMcTable(mcCollisionIndex, candidate.pt(), candidate.eta(), candidate.phi(), candidate.y(), candidate.e(), candidate.m(), candidate.pdgCode(), candidate.getGenStatusCode(), candidate.getHepMCStatusCode(), candidate.isPhysicalPrimary(), candidate.decayFlag(), candidate.origin());
+  DielectronMcTable(mcCollisionIndex, candidate.pt(), candidate.eta(), candidate.phi(), candidate.y(), candidate.e(), candidate.m(), candidate.pdgCode(), candidate.statusCode(), candidate.flags(), candidate.decayFlag(), candidate.origin());
 }
 
 }; // namespace jetdqutilities

PWGJE/Core/JetDerivedDataUtilities.h

@@ -18,6 +18,7 @@
 #define PWGJE_CORE_JETDERIVEDDATAUTILITIES_H_
 
 #include "Common/CCDB/EventSelectionParams.h"
+#include "Common/CCDB/RCTSelectionFlags.h"
 #include "Common/CCDB/TriggerAliases.h"
 
 #include <Rtypes.h>
@@ -32,7 +33,7 @@
 namespace jetderiveddatautilities
 {
 
 static constexpr float mPion = 0.139; // TDatabasePDG::Instance()->GetParticle(211)->Mass(); //can be removed when pion mass becomes default for unidentified tracks
 
 enum JCollisionSel {
   sel8 = 0,
@@ -53,15 +54,20 @@
 };
 
 template <typename T>
-bool selectCollision(T const& collision, std::vector<int> eventSelectionMaskBits, bool skipMBGapEvents = true)
+bool selectCollision(T const& collision, const std::vector<int>& eventSelectionMaskBits, bool skipMBGapEvents = true, bool rctSelection = true, std::string rctLabel = "CBT_hadronPID", bool rejectLimitedAcceptanceRct = false, bool requireZDCRct = false)
 {
   if (skipMBGapEvents && collision.subGeneratorId() == JCollisionSubGeneratorId::mbGap) {
     return false;
   }
+  o2::aod::rctsel::RCTFlagsChecker rctChecker;
+  rctChecker.init(rctLabel, requireZDCRct, rejectLimitedAcceptanceRct);
+  if (rctSelection && !rctChecker.checkTable(collision)) { // CBT_hadronPID given as default so that TOF is included in RCT selection to benefit from better timing for tracks. Impact of this for inclusive jets should be studied
+    return false;
+  }
   if (eventSelectionMaskBits.size() == 0) {
     return true;
   }
   for (auto eventSelectionMaskBit : eventSelectionMaskBits) {
     if (!(collision.eventSel() & (1 << eventSelectionMaskBit))) {
       return false;
     }
@@ -69,7 +75,7 @@
   return true;
 }
 
-bool eventSelectionMasksContainSelection(std::string eventSelectionMasks, std::string selection)
+bool eventSelectionMasksContainSelection(const std::string& eventSelectionMasks, std::string selection)
 {
   size_t position = 0;
   while ((position = eventSelectionMasks.find(selection, position)) != std::string::npos) {
@@ -83,7 +89,7 @@
   return false;
 }
 
-std::vector<int> initialiseEventSelectionBits(std::string eventSelectionMasks)
+std::vector<int> initialiseEventSelectionBits(const std::string& eventSelectionMasks)
 {
   std::vector<int> eventSelectionMaskBits;
   if (eventSelectionMasksContainSelection(eventSelectionMasks, "sel8")) {
@@ -230,7 +236,7 @@
 };
 
 template <typename T>
-bool selectTrigger(T const& collision, std::vector<int> triggerMaskBits)
+bool selectTrigger(T const& collision, const std::vector<int>& triggerMaskBits)
 {
   if (triggerMaskBits.size() == 0) {
     return true;
@@ -252,7 +258,7 @@
   return collision.triggerSel() & (1 << triggerMaskBit);
 }
 
-bool triggerMasksContainTrigger(std::string triggerMasks, std::string trigger)
+bool triggerMasksContainTrigger(const std::string& triggerMasks, std::string trigger)
 {
   size_t position = 0;
   while ((position = triggerMasks.find(trigger, position)) != std::string::npos) {
@@ -266,7 +272,7 @@
   return false;
 }
 
-std::vector<int> initialiseTriggerMaskBits(std::string triggerMasks)
+std::vector<int> initialiseTriggerMaskBits(const std::string& triggerMasks)
 {
   std::vector<int> triggerMaskBits;
   if (triggerMasksContainTrigger(triggerMasks, "fJetChLowPt")) {
@@ -335,7 +341,7 @@
   return triggerMaskBits;
 }
 
-uint64_t setTriggerSelectionBit(std::vector<bool> triggerDecisions)
+uint64_t setTriggerSelectionBit(const std::vector<bool>& triggerDecisions)
 {
   uint64_t bit = 0;
   for (std::vector<bool>::size_type i = 0; i < triggerDecisions.size(); i++) {
@@ -363,7 +369,7 @@
   return (collision.chargedTriggerSel() & (1 << triggerSelection));
 }
 
-int initialiseChargedTriggerSelection(std::string triggerSelection)
+int initialiseChargedTriggerSelection(const std::string& triggerSelection)
 {
   if (triggerSelection == "jetChLowPt") {
     return JTrigSelCh::jetChLowPt;
@@ -428,7 +434,7 @@
   return (collision.fullTriggerSel() & (1 << triggerSelection));
 }
 
-int initialiseFullTriggerSelection(std::string triggerSelection)
+int initialiseFullTriggerSelection(const std::string& triggerSelection)
 {
   if (triggerSelection == "emcalReadout") {
     return JTrigSelFull::emcalReadout;
@@ -523,7 +529,7 @@
   return (collision.chargedHFTriggerSel() & (1 << triggerSelection));
 }
 
-int initialiseChargedHFTriggerSelection(std::string triggerSelection)
+int initialiseChargedHFTriggerSelection(const std::string& triggerSelection)
 {
   if (triggerSelection == "jetD0ChLowPt") {
     return JTrigSelChHF::jetD0ChLowPt;
@@ -586,7 +592,7 @@
   return (track.trackSel() & (1 << trackSelection));
 }
 
-int initialiseTrackSelection(std::string trackSelection)
+int initialiseTrackSelection(const std::string& trackSelection)
 {
   if (trackSelection == "globalTracks") {
     return JTrackSel::globalTrack;
@@ -642,7 +648,7 @@
 template <typename T>
 bool selectTrackDcaZ(T const& track, double dcaZmax = 99.)
 {
   return abs(track.dcaZ()) < dcaZmax;
 }
 
 } // namespace jetderiveddatautilities