reemoving global variables from jet finder

Author: Nima Zardoshti

PWGJE/TableProducer/jetfinder.cxx

@@ -35,9 +35,6 @@ struct JetFinderTask {
   OutputObj<TH1F> hJetEta{"h_jet_eta"};
   OutputObj<TH1F> hJetNTracks{"h_jet_ntracks"};
 
-  Service<O2DatabasePDG> pdg;
-  std::string trackSelection;
-
   // event level configurables
   Configurable<float> vertexZCut{"vertexZCut", 10.0f, "Accepted z-vertex range"};
 
@@ -74,16 +71,15 @@ struct JetFinderTask {
   Configurable<bool> DoRhoAreaSub{"DoRhoAreaSub", false, "do rho area subtraction"};
   Configurable<bool> DoConstSub{"DoConstSub", false, "do constituent subtraction"};
 
+  Service<O2DatabasePDG> pdg;
+  std::string trackSelection;
+
+  JetFinder jetFinder;
+  std::vector<fastjet::PseudoJet> inputParticles;
+
   void init(InitContext const&)
   {
-    // variables passed to the common header
-    trackEtaMin_ = static_cast<float>(trackEtaMin);
-    trackEtaMax_ = static_cast<float>(trackEtaMax);
-    jetRadius_ = static_cast<std::vector<double>>(jetRadius);
-    DoConstSub_ = static_cast<bool>(DoConstSub);
-
     trackSelection = static_cast<std::string>(trackSelections);
-    jetTypeParticleLevelCheck = jetTypeParticleLevel;
 
     h2JetPt.setObject(new TH2F("h2_jet_pt", "jet p_{T};p_{T} (GeV/#it{c})",
                                100, 0., 100., 10, 0.05, 1.05));
@@ -144,9 +140,8 @@ struct JetFinderTask {
 
     LOG(debug) << "Process data charged!";
     inputParticles.clear();
-    using ArgType = std::decay_t<decltype(tracks)>;
-    analyseTracks<ArgType, typename ArgType::iterator>(tracks, trackSelection);
-    findJets(collision, jetsTable, constituentsTable, constituentsSubTable);
+    analyseTracks<JetTracks, JetTracks::iterator>(inputParticles, tracks, trackSelection);
+    findJets(jetFinder, inputParticles, jetRadius, collision, jetsTable, constituentsTable, constituentsSubTable,DoConstSub);
   }
 
   PROCESS_SWITCH(JetFinderTask, processChargedJets, "Data jet finding for charged jets", false);
@@ -159,8 +154,8 @@ struct JetFinderTask {
     }
     LOG(debug) << "Process data neutral!";
     inputParticles.clear();
-    analyseClusters(&clusters);
-    findJets(collision, jetsTable, constituentsTable, constituentsSubTable);
+    analyseClusters(inputParticles, &clusters);
+    findJets(jetFinder, inputParticles, jetRadius, collision, jetsTable, constituentsTable, constituentsSubTable,DoConstSub);
   }
   PROCESS_SWITCH(JetFinderTask, processNeutralJets, "Data jet finding for neutral jets", false);
 
@@ -173,20 +168,18 @@ struct JetFinderTask {
     }
     LOG(debug) << "Process data full!";
     inputParticles.clear();
-    using ArgType = std::decay_t<decltype(tracks)>;
-    analyseTracks<ArgType, typename ArgType::iterator>(tracks, trackSelection);
-    analyseClusters(&clusters);
-    findJets(collision, jetsTable, constituentsTable, constituentsSubTable);
+    analyseTracks<JetTracks, JetTracks::iterator>(inputParticles,tracks, trackSelection);
+    analyseClusters(inputParticles,&clusters);
+    findJets(jetFinder, inputParticles, jetRadius, collision, jetsTable, constituentsTable, constituentsSubTable,DoConstSub);
   }
 
   PROCESS_SWITCH(JetFinderTask, processFullJets, "Data jet finding for full and neutral jets", false);
 
   void processParticleLevelJets(aod::McCollision const& collision, aod::McParticles const& particles)
   {
     // TODO: MC event selection?
-    using ArgType = std::decay_t<decltype(particles)>;
-    analyseParticles<ArgType, typename ArgType::iterator>(particles, pdg->Instance());
-    findJets(collision, jetsTable, constituentsTable, constituentsSubTable);
+    analyseParticles<aod::McParticles, aod::McParticles::iterator>(inputParticles, trackEtaMin, trackEtaMax, jetTypeParticleLevel, particles, pdg->Instance());
+    findJets(jetFinder, inputParticles, jetRadius, collision, jetsTable, constituentsTable, constituentsSubTable,DoConstSub);
   }
 
   PROCESS_SWITCH(JetFinderTask, processParticleLevelJets, "Particle level jet finding", false);

PWGJE/TableProducer/jetfinder.h

@@ -40,24 +40,6 @@
 #include "PWGJE/Core/JetFinder.h"
 #include "PWGJE/Core/FastJetUtilities.h"
 
-JetFinder jetFinder;
-std::vector<fastjet::PseudoJet> jets;
-std::vector<fastjet::PseudoJet> inputParticles;
-int jetTypeParticleLevelCheck = -1;
-
-bool doHFJetFinding = false;
-int candPDG;
-int candDecay;
-
-float trackEtaMin_;
-float trackEtaMax_;
-std::vector<double> jetRadius_;
-float candYMin_;
-float candYMax_;
-float candPtMin_;
-float candPtMax_;
-bool DoConstSub_;
-
 using JetTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TracksDCA, aod::TrackSelection>>;
 using JetClusters = o2::soa::Filtered<o2::aod::EMCALClusters>;
 
@@ -91,7 +73,7 @@ bool selectTrack(T const& track, std::string trackSelection)
 
 // function that adds tracks to the fastjet list, removing daughters of 2Prong candidates
 template <typename T, typename U>
-void analyseTracks(T const& tracks, std::string trackSelection, std::optional<U> const& candidate = std::nullopt)
+void analyseTracks(std::vector<fastjet::PseudoJet> &inputParticles, T const& tracks, std::string trackSelection, std::optional<U> const& candidate = std::nullopt)
 {
   for (auto& track : tracks) {
     if (!selectTrack(track, trackSelection)) {
@@ -124,7 +106,7 @@ void analyseTracks(T const& tracks, std::string trackSelection, std::optional<U>
 
 // function that adds clusters to the fastjet list
 template <typename T>
-void analyseClusters(T const& clusters)
+void analyseClusters(std::vector<fastjet::PseudoJet> &inputParticles, T const& clusters)
 {
   for (auto& cluster : *clusters) {
     // add cluster selections
@@ -134,12 +116,12 @@ void analyseClusters(T const& clusters)
 
 // function that takes any generic candidate, performs selections and adds the candidate to the fastjet list
 template <typename T>
-bool analyseCandidate(T const& candidate)
+bool analyseCandidate(std::vector<fastjet::PseudoJet> &inputParticles, int candPDG, float candPtMin, float candPtMax, float candYMin, float candYMax, T const& candidate)
 {
-  if (candidate.y(RecoDecay::getMassPDG(candPDG)) < candYMin_ || candidate.y(RecoDecay::getMassPDG(candPDG)) > candYMax_) {
+  if (candidate.y(RecoDecay::getMassPDG(candPDG)) < candYMin || candidate.y(RecoDecay::getMassPDG(candPDG)) > candYMax) {
     return false;
   }
-  if (candidate.pt() < candPtMin_ || candidate.pt() >= candPtMax_) {
+  if (candidate.pt() < candPtMin || candidate.pt() >= candPtMax) {
     return false;
   }
   FastJetUtilities::fillTracks(candidate, inputParticles, candidate.globalIndex(), static_cast<int>(JetConstituentStatus::candidateHF), RecoDecay::getMassPDG(candPDG));
@@ -148,23 +130,23 @@ bool analyseCandidate(T const& candidate)
 
 // function that checks the MC status of a candidate and then calls the function to analyseCandidates
 template <typename T>
-bool analyseCandidateMC(T const& candidate, bool rejectBackgroundMCCandidates)
+bool analyseCandidateMC(std::vector<fastjet::PseudoJet> &inputParticles, int candPDG, int candDecay, float candPtMin, float candPtMax, float candYMin, float candYMax, T const& candidate, bool rejectBackgroundMCCandidates)
 {
   if (rejectBackgroundMCCandidates && !(std::abs(candidate.flagMcMatchRec()) == 1 << candDecay)) {
     return false;
   }
-  return analyseCandidate(candidate);
+  return analyseCandidate(inputParticles, candPDG, candPtMin, candPtMax, candYMin, candYMax, candidate);
 }
 
 // function that calls the jet finding and fills the relevant tables
 template <typename T, typename U, typename V, typename W>
-void findJets(T const& collision, U& jetsTable, V& constituentsTable, W& constituentsSubTable)
+void findJets(JetFinder &jetFinder, std::vector<fastjet::PseudoJet> &inputParticles, std::vector<double> jetRadius, T const& collision, U& jetsTable, V& constituentsTable, W& constituentsSubTable, bool DoConstSub, bool doHFJetFinding = false)
 {
   // auto candidatepT = 0.0;
-  auto jetRValues = static_cast<std::vector<double>>(jetRadius_);
+  auto jetRValues = static_cast<std::vector<double>>(jetRadius);
   for (auto R : jetRValues) {
     jetFinder.jetR = R;
-    jets.clear();
+    std::vector<fastjet::PseudoJet> jets;
     fastjet::ClusterSequenceArea clusterSeq(jetFinder.findJets(inputParticles, jets));
     for (const auto& jet : jets) {
       bool isHFJet = false;
@@ -188,7 +170,7 @@ void findJets(T const& collision, U& jetsTable, V& constituentsTable, W& constit
                 jet.E(), jet.m(), jet.area(), std::round(R * 100));
       for (const auto& constituent : sorted_by_pt(jet.constituents())) {
         // need to add seperate thing for constituent subtraction
-        if (DoConstSub_) { // FIXME: needs to be addressed in Haadi's PR
+        if (DoConstSub) { // FIXME: needs to be addressed in Haadi's PR
           constituentsSubTable(jetsTable.lastIndex(), constituent.pt(), constituent.eta(), constituent.phi(),
                                constituent.E(), constituent.m(), constituent.user_index());
         }
@@ -234,23 +216,23 @@ bool checkDaughters(T const& particle, int globalIndex)
 }
 
 template <typename T, typename U>
-void analyseParticles(T const& particles, TDatabasePDG* pdg, std::optional<U> const& candidate = std::nullopt)
+void analyseParticles(std::vector<fastjet::PseudoJet> &inputParticles, float particleEtaMin, float particleEtaMax, int jetTypeParticleLevel, T const& particles, TDatabasePDG* pdg, std::optional<U> const& candidate = std::nullopt)
 {
   inputParticles.clear();
   for (auto& particle : particles) {
     // TODO: can we do this through the filter?
-    if (particle.eta() < trackEtaMin_ || particle.eta() > trackEtaMax_) {
+    if (particle.eta() < particleEtaMin || particle.eta() > particleEtaMax) {
       continue;
     }
     if (particle.getGenStatusCode() != 1) { // CHECK : Does this include HF hadrons that decay?
       continue;
     }
     auto pdgParticle = pdg->GetParticle(particle.pdgCode());
     auto pdgCharge = pdgParticle ? std::abs(pdgParticle->Charge()) : -1.0;
-    if (jetTypeParticleLevelCheck == static_cast<int>(JetType::charged) && pdgCharge < 3.0) {
+    if (jetTypeParticleLevel == static_cast<int>(JetType::charged) && pdgCharge < 3.0) {
       continue;
     }
-    if (jetTypeParticleLevelCheck == static_cast<int>(JetType::neutral) && pdgCharge != 0.0) {
+    if (jetTypeParticleLevel == static_cast<int>(JetType::neutral) && pdgCharge != 0.0) {
       continue;
     }
     if (candidate != std::nullopt) {

PWGJE/TableProducer/jetfinderhf.cxx

@@ -9,7 +9,7 @@
 // granted to it by virtue of its status as an Intergovernmental Organization
 // or submit itself to any jurisdiction.
 
-// jet finder task
+// jet finder hf task
 //
 // Authors: Nima Zardoshti, Jochen Klein
 
@@ -86,20 +86,15 @@ struct JetFinderHFTask {
   Service<O2DatabasePDG> pdg;
   std::string trackSelection;
 
+  JetFinder jetFinder;
+  std::vector<fastjet::PseudoJet> inputParticles;
+
+  int candPDG;
+  int candDecay;
+
   void init(InitContext const&)
   {
-    // variables passed to the common header
-    trackEtaMin_ = static_cast<float>(trackEtaMin);
-    trackEtaMax_ = static_cast<float>(trackEtaMax);
-    candPtMin_ = static_cast<float>(candPtMin);
-    candPtMax_ = static_cast<float>(candPtMax);
-    candYMin_ = static_cast<float>(candYMin);
-    candYMax_ = static_cast<float>(candYMax);
-    jetRadius_ = static_cast<std::vector<double>>(jetRadius);
-    DoConstSub_ = static_cast<bool>(DoConstSub);
-
     trackSelection = static_cast<std::string>(trackSelections);
-    jetTypeParticleLevelCheck = jetTypeParticleLevel;
 
     h2JetPt.setObject(new TH2F("h2_jet_pt", "jet p_{T};p_{T} (GeV/#it{c})",
                                100, 0., 100., 10, 0.05, 1.05));
@@ -160,19 +155,14 @@ struct JetFinderHFTask {
         candDecay = static_cast<int>(aod::hf_cand_2prong::DecayType::JpsiToMuMu);
       }
     }
-    doHFJetFinding = true;
   }
 
-  using JetParticles2Prong = soa::Filtered<soa::Join<aod::McParticles, aod::HfCand2ProngMcGen>>;
-  using JetParticles3Prong = soa::Filtered<soa::Join<aod::McParticles, aod::HfCand3ProngMcGen>>;
-  using JetParticlesBPlus = soa::Filtered<soa::Join<aod::McParticles, aod::HfCandBplusMcGen>>;
-
   o2::aod::EMCALClusterDefinition clusterDefinition = o2::aod::emcalcluster::getClusterDefinitionFromString(clusterDefinitionS.value);
   Filter collisionFilter = (nabs(aod::collision::posZ) < vertexZCut);
   Filter trackCuts = (aod::track::pt >= trackPtMin && aod::track::pt < trackPtMax && aod::track::eta > trackEtaMin && aod::track::eta < trackEtaMax && aod::track::phi >= trackPhiMin && aod::track::phi <= trackPhiMax);
   Filter partCuts = (aod::mcparticle::pt >= trackPtMin && aod::mcparticle::pt < trackPtMax);
   Filter clusterFilter = (o2::aod::emcalcluster::definition == static_cast<int>(clusterDefinition) && aod::emcalcluster::eta > clusterEtaMin && aod::emcalcluster::eta < clusterEtaMax && aod::emcalcluster::phi >= clusterPhiMin && aod::emcalcluster::phi <= clusterPhiMax);
-  // Filter candidateCuts = (aod::hf_cand::Pt >= candPtMin && aod::hf_cand::Pt < candPtMax); FIXME: why wont this work?
+  //Filter candidateCuts = (aod::hfcand::pt >= candPtMin && aod::hfcand::pt < candPtMax && aod::hfcand::y >= candYMin && aod::hfcand::y < candYMax);
   Filter candidateCutsD0 = (aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar);
   Filter candidateCutsLc = (aod::hf_sel_candidate_lc::isSelLcToPKPi >= selectionFlagLcToPKPi || aod::hf_sel_candidate_lc::isSelLcToPiKP >= selectionFlagLcToPiPK);
   Filter candidateCutsBPlus = (aod::hf_sel_candidate_bplus::isSelBplusToD0Pi >= selectionFlagBPlus);
@@ -187,11 +177,11 @@ struct JetFinderHFTask {
 
     for (auto& candidate : candidates) {
       inputParticles.clear();
-      if (!analyseCandidate(candidate)) {
+      if (!analyseCandidate(inputParticles, candPDG, candPtMin, candPtMax, candYMin, candYMax, candidate)) {
         continue;
       }
-      analyseTracks(tracks, trackSelection, std::optional{candidate});
-      findJets(collision, jetsTable, constituentsTable, constituentsSubTable);
+      analyseTracks(inputParticles, tracks, trackSelection, std::optional{candidate});
+      findJets(jetFinder, inputParticles, jetRadius, collision, jetsTable, constituentsTable, constituentsSubTable , DoConstSub, true);
     }
   }
 
@@ -205,11 +195,11 @@ struct JetFinderHFTask {
 
     for (auto& candidate : candidates) {
       inputParticles.clear();
-      if (!analyseCandidateMC(candidate, rejectBackgroundMCCandidates)) {
+      if (!analyseCandidateMC(inputParticles, candPDG, candDecay, candPtMin, candPtMax, candYMin, candYMax, candidate, rejectBackgroundMCCandidates)) {
         continue;
       }
-      analyseTracks(tracks, trackSelection, std::optional{candidate});
-      findJets(collision, jetsTable, constituentsTable, constituentsSubTable);
+      analyseTracks(inputParticles, tracks, trackSelection, std::optional{candidate});
+      findJets(jetFinder, inputParticles, jetRadius, collision, jetsTable, constituentsTable, constituentsSubTable, DoConstSub, true);
     }
   }
 
@@ -230,9 +220,9 @@ struct JetFinderHFTask {
       }
     }
     for (auto& candidate : candidates) {
-      analyseParticles(particles, pdg->Instance(), std::optional{candidate});
+      analyseParticles(inputParticles, trackEtaMin, trackEtaMax, jetTypeParticleLevel, particles, pdg->Instance(), std::optional{candidate});
       FastJetUtilities::fillTracks(candidate, inputParticles, candidate.globalIndex(), static_cast<int>(JetConstituentStatus::candidateHF), RecoDecay::getMassPDG(candidate.pdgCode()));
-      findJets(collision, jetsTable, constituentsTable, constituentsSubTable);
+      findJets(jetFinder, inputParticles, jetRadius, collision, jetsTable, constituentsTable, constituentsSubTable, DoConstSub, true);
     }
   }
 
@@ -241,7 +231,6 @@ struct JetFinderHFTask {
   template <typename T, typename U>
   void analyseMCGen2Prong(T const& collision, U const& particles)
   {
-    jets.clear();
     inputParticles.clear();
     LOG(debug) << "Per Event MCP";
     std::vector<JetParticles2Prong::iterator> candidates;
@@ -252,7 +241,6 @@ struct JetFinderHFTask {
   template <typename T, typename U>
   void analyseMCGen3Prong(T const& collision, U const& particles)
   {
-    jets.clear();
     inputParticles.clear();
     LOG(debug) << "Per Event MCP";
     std::vector<JetParticles3Prong::iterator> candidates;
@@ -262,7 +250,6 @@ struct JetFinderHFTask {
   template <typename T, typename U>
   void analyseMCGenBPlus(T const& collision, U const& particles)
   {
-    jets.clear();
     inputParticles.clear();
     LOG(debug) << "Per Event MCP";
     std::vector<JetParticlesBPlus::iterator> candidates;