[PWGJE] Making hadronic corrections a part of the JE derived data and adding … (#9671)

Author: nzardosh

PWGJE/Core/FastJetUtilities.h

@@ -103,11 +103,27 @@ void fillTracks(const T& constituent, std::vector<fastjet::PseudoJet>& constitue
  */
 
 template <typename T>
-void fillClusters(const T& constituent, std::vector<fastjet::PseudoJet>& constituents, int index = -99999999, int status = static_cast<int>(JetConstituentStatus::cluster))
+void fillClusters(const T& constituent, std::vector<fastjet::PseudoJet>& constituents, int index = -99999999, int hadronicCorrectionType = 0, int status = static_cast<int>(JetConstituentStatus::cluster))
 {
   if (status == static_cast<int>(JetConstituentStatus::cluster)) {
-    double clusterpt = constituent.energy() / std::cosh(constituent.eta());
-    constituents.emplace_back(clusterpt * std::cos(constituent.phi()), clusterpt * std::sin(constituent.phi()), clusterpt * std::sinh(constituent.eta()), constituent.energy());
+    float constituentEnergy = 0.0;
+    if (hadronicCorrectionType == 0) {
+      constituentEnergy = constituent.energy();
+    }
+    if (hadronicCorrectionType == 1) {
+      constituentEnergy = constituent.energyCorrectedOneTrack1();
+    }
+    if (hadronicCorrectionType == 2) {
+      constituentEnergy = constituent.energyCorrectedOneTrack2();
+    }
+    if (hadronicCorrectionType == 3) {
+      constituentEnergy = constituent.energyCorrectedAllTracks1();
+    }
+    if (hadronicCorrectionType == 4) {
+      constituentEnergy = constituent.energyCorrectedAllTracks2();
+    }
+    float constituentPt = constituentEnergy / std::cosh(constituent.eta());
+    constituents.emplace_back(constituentPt * std::cos(constituent.phi()), constituentPt * std::sin(constituent.phi()), constituentPt * std::sinh(constituent.eta()), constituentEnergy);
   }
   setFastJetUserInfo(constituents, index, status);
 }

PWGJE/Core/JetFindingUtilities.h

@@ -159,11 +159,11 @@ 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)
+void analyseClusters(std::vector<fastjet::PseudoJet>& inputParticles, T const& clusters, int hadronicCorrectionType = 0)
 {
   for (auto& cluster : *clusters) {
     // add cluster selections
-    fastjetutilities::fillClusters(cluster, inputParticles, cluster.globalIndex());
+    fastjetutilities::fillClusters(cluster, inputParticles, cluster.globalIndex(), hadronicCorrectionType);
   }
 }
 

PWGJE/Core/JetMatchingUtilities.h

@@ -427,7 +427,7 @@ float getPtSum(T const& tracksBase, U const& clustersBase, V const& tracksTag, O
   if constexpr (isEMCAL) {
     if constexpr (jetsTagIsMc) {
       for (const auto& clusterBase : clustersBase) {
-        for (const auto& clusterBaseParticleId : clusterBase.mcParticleIds()) {
+        for (const auto& clusterBaseParticleId : clusterBase.mcParticlesIds()) {
           bool isClusterMatched = false;
           for (const auto& trackTag : tracksTag) {
             if (clusterBaseParticleId != -1 && clusterBaseParticleId == trackTag.globalIndex()) {
@@ -450,7 +450,7 @@ float getPtSum(T const& tracksBase, U const& clustersBase, V const& tracksTag, O
         auto trackBaseId = trackBase.globalIndex();
         for (const auto& clusterTag : clustersTag) {
           bool isClusterMatched = false;
-          for (const auto& clusterTagParticleId : clusterTag.mcParticleIds()) {
+          for (const auto& clusterTagParticleId : clusterTag.mcParticlesIds()) {
             if (trackBaseId != -1 && trackBaseId == clusterTagParticleId) {
               ptSum += trackBase.pt();
               isClusterMatched = true;

PWGJE/Core/JetSubstructureUtilities.h

@@ -53,15 +53,15 @@ namespace jetsubstructureutilities
  * @param pseudoJet converted pseudoJet object which is passed by reference
  */
 template <typename T, typename U, typename V, typename O>
-fastjet::ClusterSequenceArea jetToPseudoJet(T const& jet, U const& /*tracks*/, V const& /*clusters*/, O const& /*candidates*/, fastjet::PseudoJet& pseudoJet)
+fastjet::ClusterSequenceArea jetToPseudoJet(T const& jet, U const& /*tracks*/, V const& /*clusters*/, O const& /*candidates*/, fastjet::PseudoJet& pseudoJet, int hadronicCorrectionType = 0)
 {
   std::vector<fastjet::PseudoJet> jetConstituents;
   for (auto& jetConstituent : jet.template tracks_as<U>()) {
     fastjetutilities::fillTracks(jetConstituent, jetConstituents, jetConstituent.globalIndex());
   }
   if constexpr (std::is_same_v<std::decay_t<typename V::iterator>, o2::aod::JetClusters::iterator> || std::is_same_v<std::decay_t<typename V::iterator>, o2::aod::JetClusters::filtered_iterator>) {
     for (auto& jetClusterConstituent : jet.template clusters_as<V>()) {
-      fastjetutilities::fillClusters(jetClusterConstituent, jetConstituents, jetClusterConstituent.globalIndex());
+      fastjetutilities::fillClusters(jetClusterConstituent, jetConstituents, jetClusterConstituent.globalIndex(), hadronicCorrectionType);
     }
   }
   if constexpr (jetcandidateutilities::isCandidateTable<O>() || jetcandidateutilities::isCandidateMcTable<O>()) {
@@ -96,11 +96,11 @@ fastjet::ClusterSequenceArea jetToPseudoJet(T const& jet, U const& /*tracks*/, V
 
 // function that returns the N-subjettiness ratio and the distance betewwen the two axes considered for tau2, in the form of a vector
 template <typename T, typename U, typename V, typename O, typename M>
-std::vector<float> getNSubjettiness(T const& jet, U const& tracks, V const& clusters, O const& candidates, std::vector<fastjet::PseudoJet>::size_type nMax, M const& reclusteringAlgorithm, bool doSoftDrop = false, float zCut = 0.1, float beta = 0.0)
+std::vector<float> getNSubjettiness(T const& jet, U const& tracks, V const& clusters, O const& candidates, std::vector<fastjet::PseudoJet>::size_type nMax, M const& reclusteringAlgorithm, bool doSoftDrop = false, float zCut = 0.1, float beta = 0.0, int hadronicCorrectionType = 0)
 {
   std::vector<float> result;
   fastjet::PseudoJet pseudoJet;
-  fastjet::ClusterSequenceArea clusterSeq(jetToPseudoJet(jet, tracks, clusters, candidates, pseudoJet));
+  fastjet::ClusterSequenceArea clusterSeq(jetToPseudoJet(jet, tracks, clusters, candidates, pseudoJet, hadronicCorrectionType));
   if (doSoftDrop) {
     fastjet::contrib::SoftDrop softDrop(beta, zCut);
     pseudoJet = softDrop(pseudoJet);

PWGJE/DataModel/Jet.h

@@ -203,8 +203,8 @@ using JetCollisionMCD = o2::soa::Join<JetCollisions, JMcCollisionLbs>::iterator;
 using JetTracks = JTracks;
 using JetTracksMCD = o2::soa::Join<JetTracks, JMcTrackLbs>;
 using JetTracksSub = JTrackSubs;
-using JetClusters = JClusters;
-using JetClustersMCD = o2::soa::Join<JClusters, JMcClusterLbs>;
+using JetClusters = o2::soa::Join<JClusters, JClustersCorrectedEnergies, JClusterTracks>;
+using JetClustersMCD = o2::soa::Join<JClusters, JClustersCorrectedEnergies, JClusterTracks, JMcClusterLbs>;
 
 using JetMcCollisions = JMcCollisions;
 using JetMcCollision = JetMcCollisions::iterator;

PWGJE/DataModel/JetReducedData.h

@@ -358,9 +358,23 @@ DECLARE_SOA_TABLE_STAGED(JClusterPIs, "JCLUSTERPI",
 DECLARE_SOA_TABLE_STAGED(JClusterTracks, "JCLUSTERTRACK", //!
                          jcluster::JTrackIds);
 
+namespace jclusterhadroniccorrection
+{
+DECLARE_SOA_COLUMN(EnergyCorrectedOneTrack1, energyCorrectedOneTrack1, float);   //! with hadronic correction fraction (100%) for one matched track
+DECLARE_SOA_COLUMN(EnergyCorrectedOneTrack2, energyCorrectedOneTrack2, float);   //! with hadronic correction fraction (70%) for one matched track - systematic studies
+DECLARE_SOA_COLUMN(EnergyCorrectedAllTracks1, energyCorrectedAllTracks1, float); //! with hadronic correction fraction (100%) for all matched tracks
+DECLARE_SOA_COLUMN(EnergyCorrectedAllTracks2, energyCorrectedAllTracks2, float); //! with hadronic correction fraction (70%) for all matched tracks - for systematic studies
+} // namespace jclusterhadroniccorrection
+
+DECLARE_SOA_TABLE_STAGED(JClustersCorrectedEnergies, "JCLUSTERCORRE",            //! if this table changes it needs to be reflected in FastJetUtilities.h!!
+                         jclusterhadroniccorrection::EnergyCorrectedOneTrack1,   // corrected cluster energy for 1 matched track (f = 100%)
+                         jclusterhadroniccorrection::EnergyCorrectedOneTrack2,   // corrected cluster energy for 1 matched track (f = 70%)
+                         jclusterhadroniccorrection::EnergyCorrectedAllTracks1,  // corrected cluster energy for all matched tracks (f = 100%)
+                         jclusterhadroniccorrection::EnergyCorrectedAllTracks2); // corrected cluster energy for all matched tracks (f = 70%)
+
 namespace jmcclusterlb
 {
-DECLARE_SOA_ARRAY_INDEX_COLUMN(JMcParticle, mcParticle);
+DECLARE_SOA_ARRAY_INDEX_COLUMN(JMcParticle, mcParticles);
 DECLARE_SOA_COLUMN(AmplitudeA, amplitudeA, std::vector<float>);
 } // namespace jmcclusterlb
 

PWGJE/DataModel/emcalClusterHadronicCorrectionTask.h

@@ -65,6 +65,7 @@ struct JetFinderTask {
   Configurable<float> clusterTimeMin{"clusterTimeMin", -25., "minimum Cluster time (ns)"};
   Configurable<float> clusterTimeMax{"clusterTimeMax", 25., "maximum Cluster time (ns)"};
   Configurable<bool> clusterRejectExotics{"clusterRejectExotics", true, "Reject exotic clusters"};
+  Configurable<int> hadronicCorrectionType{"hadronicCorrectionType", 0, "0 = no correction, 1 = CorrectedOneTrack1, 2 = CorrectedOneTrack2, 3 = CorrectedAllTracks1, 4 = CorrectedAllTracks2"};
   Configurable<bool> doEMCALEventSelection{"doEMCALEventSelection", true, "apply the selection to the event alias_bit for full and neutral jets"};
   Configurable<bool> doEMCALEventSelectionChargedJets{"doEMCALEventSelectionChargedJets", false, "apply the selection to the event alias_bit for charged jets"};
 
@@ -195,7 +196,7 @@ struct JetFinderTask {
     }
     inputParticles.clear();
     jetfindingutilities::analyseTracks<soa::Filtered<aod::JetTracks>, soa::Filtered<aod::JetTracks>::iterator>(inputParticles, tracks, trackSelection, trackingEfficiency);
-    jetfindingutilities::analyseClusters(inputParticles, &clusters);
+    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);
   }
   PROCESS_SWITCH(JetFinderTask, processFullJets, "Data and reco level jet finding for full and neutral jets", false);

PWGJE/JetFinders/jetFinder.cxx

@@ -28,6 +28,11 @@ o2physics_add_dpl_workflow(jet-deriveddata-producer-dummy
                     PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)
 
+o2physics_add_dpl_workflow(jet-deriveddata-producer-dummy-charged
+                    SOURCES derivedDataProducerDummyCharged.cxx
+                    PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore
+                    COMPONENT_NAME Analysis)
+
 o2physics_add_dpl_workflow(jet-deriveddata-producer-dummy-d0
                     SOURCES derivedDataProducerDummyD0.cxx
                     PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore