Skip to content

PWGJE: removing global variables from jet finder #2473

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
nzardosh merged 1 commit into from
Apr 24, 2023
+46 −84
Merged

PWGJE: removing global variables from jet finder #2473

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
37 changes: 15 additions & 22 deletions PWGJE/TableProducer/jetfinder.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"};

Expand Down Expand Up @@ -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));
Expand Down Expand Up @@ -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);
Expand All @@ -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);

Expand All @@ -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);
Expand Down
48 changes: 15 additions & 33 deletions PWGJE/TableProducer/jetfinder.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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>;

Expand Down Expand Up @@ -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)) {
Expand Down Expand Up @@ -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
Expand All @@ -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));
Expand All @@ -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;
Expand All @@ -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());
}
Expand Down Expand Up @@ -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) {
Expand Down
45 changes: 16 additions & 29 deletions PWGJE/TableProducer/jetfinderhf.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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));
Expand Down Expand Up @@ -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);
Expand All @@ -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);
}
}

Expand All @@ -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);
}
}

Expand All @@ -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);
}
}

Expand All @@ -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;
Expand All @@ -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;
Expand All @@ -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;
Expand Down