[PWGJE] V0 Jet Framework update

* Added option to either save all jets or only those containing V0 candidates

* Added option to select all MCP V0s or only those that decay to charged particles

Author: GijsvWeelden

PWGJE/Core/JetFindingUtilities.h

@@ -271,7 +271,7 @@ bool analyseV0s(std::vector<fastjet::PseudoJet>& inputParticles, T const& v0s, f
  * @param doHFJetFinding set whether only jets containing a HF candidate are saved
  */
 template <typename T, typename U, typename V>
-void findJets(JetFinder& jetFinder, std::vector<fastjet::PseudoJet>& inputParticles, float jetPtMin, float jetPtMax, std::vector<double> jetRadius, float jetAreaFractionMin, T const& collision, U& jetsTable, V& constituentsTable, std::shared_ptr<THn> thnSparseJet, bool fillThnSparse, bool doCandidateJetFinding = false)
+void findJets(JetFinder& jetFinder, std::vector<fastjet::PseudoJet>& inputParticles, float jetPtMin, float jetPtMax, std::vector<double> jetRadius, float jetAreaFractionMin, T const& collision, U& jetsTable, V& constituentsTable, std::shared_ptr<THn> thnSparseJet, bool fillThnSparse, bool doCandidateJetFinding = false, bool saveJetsWithoutCandidates = false)
 {
   auto jetRValues = static_cast<std::vector<double>>(jetRadius);
   jetFinder.jetPtMin = jetPtMin;
@@ -296,7 +296,7 @@ void findJets(JetFinder& jetFinder, std::vector<fastjet::PseudoJet>& inputPartic
             break;
           }
         }
-        if (!isCandidateJet) {
+        if (!isCandidateJet && !saveJetsWithoutCandidates) {
           continue;
         }
       }

PWGJE/Core/JetV0Utilities.h

@@ -138,9 +138,13 @@ auto slicedPerV0Candidate(T const& table, U const& candidate, V const& perV0Cand
 }
 
 template <typename T, typename U>
-bool isV0Particle(T const& particles, U const& particle)
+bool isV0Particle(T const& particles, U const& particle, bool onlyChargedDecays)
 {
-  return RecoDecay::isMatchedMCGen(particles, particle, +kK0Short, std::array{+kPiPlus, -kPiPlus}, true) || RecoDecay::isMatchedMCGen(particles, particle, +kLambda0, std::array{+kProton, -kPiPlus}, true) || RecoDecay::isMatchedMCGen(particles, particle, -kLambda0, std::array{-kProton, +kPiPlus}, true);
+  if (onlyChargedDecays) {
+    return RecoDecay::isMatchedMCGen(particles, particle, +kK0Short, std::array{+kPiPlus, -kPiPlus}, true) || RecoDecay::isMatchedMCGen(particles, particle, +kLambda0, std::array{+kProton, -kPiPlus}, true) || RecoDecay::isMatchedMCGen(particles, particle, -kLambda0, std::array{-kProton, +kPiPlus}, true);
+  } else {
+    return RecoDecay::isMatchedMCGen(particles, particle, +kK0Short, std::array{+kPiPlus, -kPiPlus}, true) || RecoDecay::isMatchedMCGen(particles, particle, +kK0Short, std::array{+kPi0, +kPi0}, true) || RecoDecay::isMatchedMCGen(particles, particle, +kLambda0, std::array{+kProton, -kPiPlus}, true) || RecoDecay::isMatchedMCGen(particles, particle, +kLambda0, std::array{+kNeutron, +kPi0}, true) || RecoDecay::isMatchedMCGen(particles, particle, -kLambda0, std::array{-kProton, +kPiPlus}, true) || RecoDecay::isMatchedMCGen(particles, particle, -kLambda0, std::array{+kNeutron, +kPi0}, true);
+  }
 }
 
 enum JV0ParticleDecays {

PWGJE/JetFinders/jetFinderV0.cxx

@@ -93,6 +93,7 @@ struct JetFinderV0Task {
   Configurable<bool> fillTHnSparse{"fillTHnSparse", true, "switch to fill the THnSparse"};
   Configurable<double> jetExtraParam{"jetExtraParam", -99.0, "sets the _extra_param in fastjet"};
   Configurable<bool> useV0SignalFlags{"useV0SignalFlags", true, "use V0 signal flags table"};
+  Configurable<bool> saveJetsWithoutCandidates{"saveJetsWithoutCandidates", false, "save all jets, even those without V0s"};
 
   Service<o2::framework::O2DatabasePDG> pdgDatabase;
   int trackSelection = -1;
@@ -171,7 +172,9 @@ struct JetFinderV0Task {
     }
     inputParticles.clear();
     if (!jetfindingutilities::analyseV0s(inputParticles, candidates, candPtMin, candPtMax, candYMin, candYMax, candIndex, useV0SignalFlags)) {
-      return;
+      if (!saveJetsWithoutCandidates) {
+        return;
+      }
     }
 
     /*
@@ -183,7 +186,7 @@ struct JetFinderV0Task {
         */
     jetfindingutilities::analyseTracksMultipleCandidates(inputParticles, tracks, trackSelection, trackingEfficiency, candidates);
 
-    jetfindingutilities::findJets(jetFinder, inputParticles, minJetPt, maxJetPt, jetRadius, jetAreaFractionMin, collision, jetsTableInput, constituentsTableInput, registry.get<THn>(HIST("hJet")), fillTHnSparse, true);
+    jetfindingutilities::findJets(jetFinder, inputParticles, minJetPt, maxJetPt, jetRadius, jetAreaFractionMin, collision, jetsTableInput, constituentsTableInput, registry.get<THn>(HIST("hJet")), fillTHnSparse, true, saveJetsWithoutCandidates);
   }
 
   template <typename T, typename U, typename V>
@@ -192,10 +195,12 @@ struct JetFinderV0Task {
 
     inputParticles.clear();
     if (!jetfindingutilities::analyseV0s(inputParticles, candidates, candPtMin, candPtMax, candYMin, candYMax, candIndex, useV0SignalFlags)) {
-      return;
+      if (!saveJetsWithoutCandidates) {
+        return;
+      }
     }
     jetfindingutilities::analyseParticles<true>(inputParticles, particleSelection, jetTypeParticleLevel, particles, pdgDatabase, &candidates);
-    jetfindingutilities::findJets(jetFinder, inputParticles, minJetPt, maxJetPt, jetRadius, jetAreaFractionMin, collision, jetsTable, constituentsTable, registry.get<THn>(HIST("hJetMCP")), fillTHnSparse, true);
+    jetfindingutilities::findJets(jetFinder, inputParticles, minJetPt, maxJetPt, jetRadius, jetAreaFractionMin, collision, jetsTable, constituentsTable, registry.get<THn>(HIST("hJetMCP")), fillTHnSparse, true, saveJetsWithoutCandidates);
   }
 
   void processDummy(aod::JetCollisions const&)

PWGJE/TableProducer/derivedDataProducer.cxx

@@ -763,8 +763,9 @@ struct JetDerivedDataProducerTask {
   void processV0MC(aod::McCollision const& mcCollision, aod::McParticles const& particles)
   { // can loop over McV0Labels tables if we want to only store matched V0Particles
     bool filledV0McCollisionTable = false;
+    bool onlyChargedDecays = false;
     for (auto const& particle : particles) {
-      if (jetv0utilities::isV0Particle(particles, particle)) {
+      if (jetv0utilities::isV0Particle(particles, particle, onlyChargedDecays)) {
         if (!filledV0McCollisionTable) {
           products.jV0McCollisionsTable(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ());
           products.jV0McCollisionIdsTable(mcCollision.globalIndex());
@@ -796,6 +797,43 @@ struct JetDerivedDataProducerTask {
   }
   PROCESS_SWITCH(JetDerivedDataProducerTask, processV0MC, "produces V0 particles", false);
 
+  void processV0MCChargedDecays(aod::McCollision const& mcCollision, aod::McParticles const& particles)
+  { // can loop over McV0Labels tables if we want to only store matched V0Particles
+    bool filledV0McCollisionTable = false;
+    bool onlyChargedDecays = true;
+    for (auto const& particle : particles) {
+      if (jetv0utilities::isV0Particle(particles, particle, onlyChargedDecays)) {
+        if (!filledV0McCollisionTable) {
+          products.jV0McCollisionsTable(mcCollision.posX(), mcCollision.posY(), mcCollision.posZ());
+          products.jV0McCollisionIdsTable(mcCollision.globalIndex());
+          filledV0McCollisionTable = true;
+        }
+        std::vector<int32_t> mothersId;
+        if (particle.has_mothers()) {
+          auto mothersIdTemps = particle.mothersIds();
+          for (auto mothersIdTemp : mothersIdTemps) {
+            mothersId.push_back(mothersIdTemp);
+          }
+        }
+        int daughtersId[2] = {-1, -1};
+        auto i = 0;
+        if (particle.has_daughters()) {
+          for (auto daughterId : particle.daughtersIds()) {
+            if (i > 1) {
+              break;
+            }
+            daughtersId[i] = daughterId;
+            i++;
+          }
+        }
+        auto pdgParticle = pdgDatabase->GetParticle(particle.pdgCode());
+        products.jV0McsTable(products.jV0McCollisionsTable.lastIndex(), particle.pt(), particle.eta(), particle.phi(), particle.y(), particle.e(), pdgParticle->Mass(), particle.pdgCode(), particle.getGenStatusCode(), particle.getHepMCStatusCode(), particle.isPhysicalPrimary(), jetv0utilities::setV0ParticleDecayBit(particles, particle));
+        products.jV0McIdsTable(mcCollision.globalIndex(), particle.globalIndex(), mothersId, daughtersId);
+      }
+    }
+  }
+  PROCESS_SWITCH(JetDerivedDataProducerTask, processV0MCChargedDecays, "produces V0 particles that decay to charged particles", false);
+
   void processDielectronCollisions(aod::ReducedEventsInfo::iterator const& DielectronCollision)
   {
     products.jDielectronCollisionIdsTable(DielectronCollision.collisionId());