[PWGJE] adding ability to do jet finding over multiple cluster definitions an… (#14961)

Author: nzardosh

PWGJE/Core/JetDerivedDataUtilities.h

@@ -17,6 +17,7 @@
 #ifndef PWGJE_CORE_JETDERIVEDDATAUTILITIES_H_
 #define PWGJE_CORE_JETDERIVEDDATAUTILITIES_H_
 
+#include "PWGJE/DataModel/EMCALClusters.h"
 #include "PWGUD/Core/SGSelector.h"
 
 #include "Common/CCDB/EventSelectionParams.h"
@@ -708,6 +709,25 @@ bool selectTrackDcaZ(T const& track, double dcaZmax = 99.)
   return std::abs(track.dcaZ()) < dcaZmax;
 }
 
+std::vector<int> initialiseClusterDefinitions(const std::string clusterDefinitions)
+{
+  std::vector<int> clusterDefinitionsVec;
+  if (clusterDefinitions.empty()) {
+    return clusterDefinitionsVec;
+  }
+  size_t start = 0;
+  size_t end;
+  while ((end = clusterDefinitions.find(',', start)) != std::string::npos) {
+    clusterDefinitionsVec.push_back(static_cast<int>(o2::aod::emcalcluster::getClusterDefinitionFromString(clusterDefinitions.substr(start, end - start))));
+    start = end + 1;
+  }
+  // Process the last element (after the final comma or if no comma exists)
+  if (start < clusterDefinitions.length()) {
+    clusterDefinitionsVec.push_back(static_cast<int>(o2::aod::emcalcluster::getClusterDefinitionFromString(clusterDefinitions.substr(start))));
+  }
+  return clusterDefinitionsVec;
+}
+
 } // namespace jetderiveddatautilities
 
 #endif // PWGJE_CORE_JETDERIVEDDATAUTILITIES_H_

PWGJE/Core/JetFindingUtilities.h

@@ -156,10 +156,13 @@ void analyseTracksMultipleCandidates(std::vector<fastjet::PseudoJet>& inputParti
  * @param clusters track table to be added
  */
 template <typename T>
-void analyseClusters(std::vector<fastjet::PseudoJet>& inputParticles, T const& clusters, int hadronicCorrectionType = 0)
+void analyseClusters(std::vector<fastjet::PseudoJet>& inputParticles, T const& clusters, int clusterDefinition, int hadronicCorrectionType)
 {
-  for (auto& cluster : *clusters) {
+  for (auto const& cluster : clusters) {
     // add cluster selections
+    if (cluster.definition() != clusterDefinition) {
+      continue;
+    }
     fastjetutilities::fillClusters(cluster, inputParticles, cluster.globalIndex(), hadronicCorrectionType);
   }
 }

PWGJE/JetFinders/jetFinder.h

@@ -76,7 +76,7 @@ struct JetFinderTask {
   o2::framework::Configurable<std::string> particleSelections{"particleSelections", "PhysicalPrimary", "set particle selections"};
 
   // cluster level configurables
-  o2::framework::Configurable<std::string> clusterDefinitionS{"clusterDefinition", "kV3Default", "cluster definition to be selected, e.g. V3Default"};
+  o2::framework::Configurable<std::string> clusterDefinitions{"clusterDefinitions", "kV3Default", "cluster definition to be selected, e.g. V3Default"};
   o2::framework::Configurable<float> clusterEtaMin{"clusterEtaMin", -0.71, "minimum cluster eta"}; // For ECMAL: |eta| < 0.7, phi = 1.40 - 3.26
   o2::framework::Configurable<float> clusterEtaMax{"clusterEtaMax", 0.71, "maximum cluster eta"};  // For ECMAL: |eta| < 0.7, phi = 1.40 - 3.26
   o2::framework::Configurable<float> clusterPhiMin{"clusterPhiMin", 1.39, "minimum cluster phi"};
@@ -85,7 +85,7 @@ struct JetFinderTask {
   o2::framework::Configurable<float> clusterTimeMin{"clusterTimeMin", -25., "minimum Cluster time (ns)"};
   o2::framework::Configurable<float> clusterTimeMax{"clusterTimeMax", 25., "maximum Cluster time (ns)"};
   o2::framework::Configurable<bool> clusterRejectExotics{"clusterRejectExotics", true, "Reject exotic clusters"};
-  o2::framework::Configurable<int> hadronicCorrectionType{"hadronicCorrectionType", 0, "0 = no correction, 1 = CorrectedOneTrack1, 2 = CorrectedOneTrack2, 3 = CorrectedAllTracks1, 4 = CorrectedAllTracks2"};
+  o2::framework::Configurable<std::vector<int>> hadronicCorrectionTypes{"hadronicCorrectionTypes", {0}, "0 = no correction, 1 = CorrectedOneTrack1, 2 = CorrectedOneTrack2, 3 = CorrectedAllTracks1, 4 = CorrectedAllTracks2 note :analyses should only use one"};
   o2::framework::Configurable<bool> doEMCALEventSelection{"doEMCALEventSelection", true, "apply the selection to the event alias_bit for full and neutral jets"};
   o2::framework::Configurable<bool> doEMCALEventSelectionChargedJets{"doEMCALEventSelectionChargedJets", false, "apply the selection to the event alias_bit for charged jets"};
 
@@ -119,16 +119,18 @@ struct JetFinderTask {
 
   std::vector<int> triggerMaskBits;
 
-  o2::aod::EMCALClusterDefinition clusterDefinition;
+  std::vector<int> clusterDefinitionsVec;
+  std::vector<int> hadronicCorrectionTypesVec;
 
   void init(o2::framework::InitContext const&)
   {
+    hadronicCorrectionTypesVec = (std::vector<int>)hadronicCorrectionTypes;
     eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast<std::string>(eventSelections));
     triggerMaskBits = jetderiveddatautilities::initialiseTriggerMaskBits(triggerMasks);
     trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast<std::string>(trackSelections));
     particleSelection = static_cast<std::string>(particleSelections);
 
-    clusterDefinition = o2::aod::emcalcluster::getClusterDefinitionFromString(clusterDefinitionS.value);
+    clusterDefinitionsVec = jetderiveddatautilities::initialiseClusterDefinitions(clusterDefinitions.value);
 
     jetFinder.etaMin = trackEtaMin;
     jetFinder.etaMax = trackEtaMax;
@@ -185,7 +187,7 @@ struct JetFinderTask {
   o2::framework::expressions::Filter mcCollisionFilter = (nabs(o2::aod::jmccollision::posZ) < vertexZCut);
   o2::framework::expressions::Filter trackCuts = (o2::aod::jtrack::pt >= trackPtMin && o2::aod::jtrack::pt < trackPtMax && o2::aod::jtrack::eta >= trackEtaMin && o2::aod::jtrack::eta <= trackEtaMax && o2::aod::jtrack::phi >= trackPhiMin && o2::aod::jtrack::phi <= trackPhiMax); // do we need eta cut both here and in globalselection?
   o2::framework::expressions::Filter partCuts = (o2::aod::jmcparticle::pt >= trackPtMin && o2::aod::jmcparticle::pt < trackPtMax && o2::aod::jmcparticle::eta >= trackEtaMin && o2::aod::jmcparticle::eta <= trackEtaMax && o2::aod::jmcparticle::phi >= trackPhiMin && o2::aod::jmcparticle::phi <= trackPhiMax);
-  o2::framework::expressions::Filter clusterFilter = (o2::aod::jcluster::definition == static_cast<int>(clusterDefinition) && o2::aod::jcluster::eta >= clusterEtaMin && o2::aod::jcluster::eta <= clusterEtaMax && o2::aod::jcluster::phi >= clusterPhiMin && o2::aod::jcluster::phi <= clusterPhiMax && o2::aod::jcluster::energy >= clusterEnergyMin && o2::aod::jcluster::time > clusterTimeMin && o2::aod::jcluster::time < clusterTimeMax && (!clusterRejectExotics || o2::aod::jcluster::isExotic != true));
+  o2::framework::expressions::Filter clusterFilter = (o2::aod::jcluster::eta >= clusterEtaMin && o2::aod::jcluster::eta <= clusterEtaMax && o2::aod::jcluster::phi >= clusterPhiMin && o2::aod::jcluster::phi <= clusterPhiMax && o2::aod::jcluster::energy >= clusterEnergyMin && o2::aod::jcluster::time > clusterTimeMin && o2::aod::jcluster::time < clusterTimeMax && (!clusterRejectExotics || o2::aod::jcluster::isExotic != true));
 
   void processChargedJets(o2::soa::Filtered<o2::aod::JetCollisions>::iterator const& collision,
                           o2::soa::Filtered<o2::aod::JetTracks> const& tracks)
@@ -219,9 +221,13 @@ struct JetFinderTask {
     if ((doEMCALEventSelection && !jetderiveddatautilities::eventEMCAL(collision)) || !jetderiveddatautilities::selectCollision(collision, eventSelectionBits, skipMBGapEvents, applyRCTSelections, "CBT_calo") || !jetderiveddatautilities::selectTrigger(collision, triggerMaskBits)) {
       return;
     }
-    inputParticles.clear();
-    jetfindingutilities::analyseClusters(inputParticles, &clusters);
-    jetfindingutilities::findJets(jetFinder, inputParticles, jetPtMin, jetPtMax, jetRadius, jetAreaFractionMin, collision, jetsTable, constituentsTable, fillTHnSparse ? registry.get<THn>(HIST("hJet")) : std::shared_ptr<THn>(nullptr), fillTHnSparse);
+    for (auto const& clusterDefinition : clusterDefinitionsVec) {
+      for (auto const& hadronicCorrectionType : hadronicCorrectionTypesVec) {
+        inputParticles.clear();
+        jetfindingutilities::analyseClusters(inputParticles, clusters, clusterDefinition, hadronicCorrectionType);
+        jetfindingutilities::findJets(jetFinder, inputParticles, jetPtMin, jetPtMax, jetRadius, jetAreaFractionMin, collision, jetsTable, constituentsTable, fillTHnSparse ? registry.get<THn>(HIST("hJet")) : std::shared_ptr<THn>(nullptr), fillTHnSparse);
+      }
+    }
   }
   PROCESS_SWITCH(JetFinderTask, processNeutralJets, "Data and reco level jet finding for neutral jets", false);
 
@@ -232,10 +238,14 @@ struct JetFinderTask {
     if ((doEMCALEventSelection && !jetderiveddatautilities::eventEMCAL(collision)) || !jetderiveddatautilities::selectCollision(collision, eventSelectionBits, skipMBGapEvents, applyRCTSelections, "CBT_calo") || !jetderiveddatautilities::selectTrigger(collision, triggerMaskBits)) {
       return;
     }
-    inputParticles.clear();
-    jetfindingutilities::analyseTracks<o2::soa::Filtered<o2::aod::JetTracks>, o2::soa::Filtered<o2::aod::JetTracks>::iterator>(inputParticles, tracks, trackSelection);
-    jetfindingutilities::analyseClusters(inputParticles, &clusters, hadronicCorrectionType);
-    jetfindingutilities::findJets(jetFinder, inputParticles, jetPtMin, jetPtMax, jetRadius, jetAreaFractionMin, collision, jetsTable, constituentsTable, fillTHnSparse ? registry.get<THn>(HIST("hJet")) : std::shared_ptr<THn>(nullptr), fillTHnSparse);
+    for (auto const& clusterDefinition : clusterDefinitionsVec) {
+      for (auto const& hadronicCorrectionType : hadronicCorrectionTypesVec) {
+        inputParticles.clear();
+        jetfindingutilities::analyseTracks<o2::soa::Filtered<o2::aod::JetTracks>, o2::soa::Filtered<o2::aod::JetTracks>::iterator>(inputParticles, tracks, trackSelection);
+        jetfindingutilities::analyseClusters(inputParticles, clusters, clusterDefinition, hadronicCorrectionType);
+        jetfindingutilities::findJets(jetFinder, inputParticles, jetPtMin, jetPtMax, jetRadius, jetAreaFractionMin, collision, jetsTable, constituentsTable, fillTHnSparse ? registry.get<THn>(HIST("hJet")) : std::shared_ptr<THn>(nullptr), fillTHnSparse);
+      }
+    }
   }
   PROCESS_SWITCH(JetFinderTask, processFullJets, "Data and reco level jet finding for full and neutral jets", false);
 

PWGJE/TableProducer/derivedDataWriter.cxx

@@ -55,7 +55,6 @@ struct JetDerivedDataWriter {
     Configurable<float> trackEtaSelectionMax{"trackEtaSelectionMax", 0.9, "only save tracks that have an eta smaller than this eta"};
     Configurable<bool> savePartonLevelInfo{"savePartonLevelInfo", true, "save parton level info at MCP level"};
 
-    Configurable<std::string> triggerMasks{"triggerMasks", "", "possible JE Trigger masks: fJetChLowPt,fJetChHighPt,fTrackLowPt,fTrackHighPt,fJetD0ChLowPt,fJetD0ChHighPt,fJetLcChLowPt,fJetLcChHighPt,fEMCALReadout,fJetFullHighPt,fJetFullLowPt,fJetNeutralHighPt,fJetNeutralLowPt,fGammaVeryHighPtEMCAL,fGammaVeryHighPtDCAL,fGammaHighPtEMCAL,fGammaHighPtDCAL,fGammaLowPtEMCAL,fGammaLowPtDCAL,fGammaVeryLowPtEMCAL,fGammaVeryLowPtDCAL"};
   } config;
 
   struct : ProducesGroup {
@@ -79,7 +78,6 @@ struct JetDerivedDataWriter {
     Produces<aod::StoredJMcParticles> storedJMcParticlesTable;
     Produces<aod::StoredJMcParticlePIs> storedJParticlesParentIndexTable;
     Produces<aod::StoredJClusters> storedJClustersTable;
-    Produces<aod::StoredJClustersCorrectedEnergies> storedJClustersCorrectedEnergiesTable;
     Produces<aod::StoredJClusterPIs> storedJClustersParentIndexTable;
     Produces<aod::StoredJClusterTracks> storedJClustersMatchedTracksTable;
     Produces<aod::StoredJMcClusterLbs> storedJMcClustersLabelTable;
@@ -530,15 +528,14 @@ struct JetDerivedDataWriter {
   }
   PROCESS_SWITCH(JetDerivedDataWriter, processTracks, "write out output tables for tracks", true);
 
-  void processClusters(soa::Join<aod::JCollisions, aod::JEMCCollisionLbs, aod::JCollisionSelections>::iterator const& collision, aod::JTracks const&, aod::JEMCTracks const& emcTracks, soa::Join<aod::JClusters, aod::JClustersCorrectedEnergies, aod::JClusterPIs, aod::JClusterTracks> const& clusters)
+  void processClusters(soa::Join<aod::JCollisions, aod::JEMCCollisionLbs, aod::JCollisionSelections>::iterator const& collision, aod::JTracks const&, aod::JEMCTracks const& emcTracks, soa::Join<aod::JClusters, aod::JClusterPIs, aod::JClusterTracks> const& clusters)
   {
     if (collision.isCollisionSelected()) {
       products.storedJCollisionsEMCalLabelTable(collision.isAmbiguous(), collision.isEmcalReadout());
       for (const auto& cluster : clusters) {
         products.storedJClustersTable(collisionMapping[collision.globalIndex()], cluster.id(), cluster.energy(), cluster.coreEnergy(), cluster.rawEnergy(),
                                       cluster.eta(), cluster.phi(), cluster.m02(), cluster.m20(), cluster.nCells(), cluster.time(), cluster.isExotic(), cluster.distanceToBadChannel(),
                                       cluster.nlm(), cluster.definition(), cluster.leadingCellEnergy(), cluster.subleadingCellEnergy(), cluster.leadingCellNumber(), cluster.subleadingCellNumber());
-        products.storedJClustersCorrectedEnergiesTable(cluster.energyCorrectedOneTrack1(), cluster.energyCorrectedOneTrack2(), cluster.energyCorrectedAllTracks1(), cluster.energyCorrectedAllTracks2());
         products.storedJClustersParentIndexTable(cluster.clusterId());
 
         std::vector<int32_t> clusterStoredJTrackIDs;