Skip to content

[PWGJE] additional features for GammaJetTrees #8565

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
ktf merged 7 commits into from
Nov 25, 2024
Merged

[PWGJE] additional features for GammaJetTrees #8565

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
242d740
added TM, occupancy, multiple Jet radii and more
fjonasALICE Nov 18, 2024
89dc3d4
changed occupancy to int and added QA
fjonasALICE Nov 18, 2024
7bb11ca
fixing merge conflicts
fjonasALICE Nov 18, 2024
8337d62
adding emcal eta phi on surface to tree
fjonasALICE Nov 20, 2024
d0a7e3f
adding correct namespace for JE tables
fjonasALICE Nov 22, 2024
6d4d72c
Please consider the following formatting changes
alibuild Nov 22, 2024
a3b7e2e
Merge pull request #4 from alibuild/alibot-cleanup-8565
fjonasALICE Nov 22, 2024
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
8 changes: 6 additions & 2 deletions PWGJE/DataModel/GammaJetAnalysisTree.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,9 +31,10 @@ DECLARE_SOA_COLUMN(Multiplicity, multiplicity, float);
DECLARE_SOA_COLUMN(Centrality, centrality, float);
DECLARE_SOA_COLUMN(Rho, rho, float);
DECLARE_SOA_COLUMN(EventSel, eventSel, uint8_t);
DECLARE_SOA_COLUMN(Occupancy, occupancy, int);
DECLARE_SOA_BITMAP_COLUMN(Alias, alias, 32);
} // namespace gjevent
DECLARE_SOA_TABLE(GjEvents, "AOD", "GJEVENT", o2::soa::Index<>, gjevent::Multiplicity, gjevent::Centrality, gjevent::Rho, gjevent::EventSel, gjevent::Alias)
DECLARE_SOA_TABLE(GjEvents, "AOD", "GJEVENT", o2::soa::Index<>, gjevent::Multiplicity, gjevent::Centrality, gjevent::Rho, gjevent::EventSel, gjevent::Occupancy, gjevent::Alias)

using GjEvent = GjEvents::iterator;

Expand Down Expand Up @@ -64,12 +65,15 @@ DECLARE_SOA_INDEX_COLUMN(GjEvent, gjevent);
DECLARE_SOA_COLUMN(Pt, pt, float);
DECLARE_SOA_COLUMN(Eta, eta, float);
DECLARE_SOA_COLUMN(Phi, phi, float);
DECLARE_SOA_COLUMN(Radius, radius, float);
DECLARE_SOA_COLUMN(Energy, energy, float);
DECLARE_SOA_COLUMN(Mass, mass, float);
DECLARE_SOA_COLUMN(Area, area, float);
DECLARE_SOA_COLUMN(LeadingTrackPt, leadingtrackpt, float);
DECLARE_SOA_COLUMN(PerpConeRho, perpconerho, float);
DECLARE_SOA_COLUMN(NConstituents, nConstituents, ushort);
} // namespace gjchjet
DECLARE_SOA_TABLE(GjChargedJets, "AOD", "GJCHJET", gjchjet::GjEventId, gjchjet::Pt, gjchjet::Eta, gjchjet::Phi, gjchjet::Energy, gjchjet::Mass, gjchjet::Area, gjchjet::NConstituents)
DECLARE_SOA_TABLE(GjChargedJets, "AOD", "GJCHJET", gjchjet::GjEventId, gjchjet::Pt, gjchjet::Eta, gjchjet::Phi, gjchjet::Radius, gjchjet::Energy, gjchjet::Mass, gjchjet::Area, gjchjet::LeadingTrackPt, gjchjet::PerpConeRho, gjchjet::NConstituents)
} // namespace o2::aod

#endif // PWGJE_DATAMODEL_GAMMAJETANALYSISTREE_H_
5 changes: 4 additions & 1 deletion PWGJE/TableProducer/emcalCorrectionTask.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -184,13 +184,15 @@ struct EmcalCorrectionTask {
o2Axis energyAxis{200, 0., 100., "E (GeV)"},
timeAxis{300, -100, 200., "t (ns)"},
etaAxis{160, -0.8, 0.8, "#eta"},
phiAxis{72, 0, 2 * 3.14159, "phi"};
phiAxis{72, 0, 2 * 3.14159, "phi"},
nlmAxis{50, -0.5, 49.5, "NLM"};
mHistManager.add("hCellE", "hCellE", o2HistType::kTH1F, {energyAxis});
mHistManager.add("hCellTowerID", "hCellTowerID", o2HistType::kTH1D, {{20000, 0, 20000}});
mHistManager.add("hCellEtaPhi", "hCellEtaPhi", o2HistType::kTH2F, {etaAxis, phiAxis});
// NOTE: Reversed column and row because it's more natural for presentation.
mHistManager.add("hCellRowCol", "hCellRowCol;Column;Row", o2HistType::kTH2D, {{97, 0, 97}, {600, 0, 600}});
mHistManager.add("hClusterE", "hClusterE", o2HistType::kTH1F, {energyAxis});
mHistManager.add("hClusterNLM", "hClusterNLM", o2HistType::kTH1F, {nlmAxis});
mHistManager.add("hClusterEtaPhi", "hClusterEtaPhi", o2HistType::kTH2F, {etaAxis, phiAxis});
mHistManager.add("hClusterTime", "hClusterTime", o2HistType::kTH1F, {timeAxis});
mHistManager.add("hGlobalTrackEtaPhi", "hGlobalTrackEtaPhi", o2HistType::kTH2F, {etaAxis, phiAxis});
Expand Down Expand Up @@ -627,6 +629,7 @@ struct EmcalCorrectionTask {
} // end of cells of cluser loop
// fill histograms
mHistManager.fill(HIST("hClusterE"), cluster.E());
mHistManager.fill(HIST("hClusterNLM"), cluster.getNExMax());
mHistManager.fill(HIST("hClusterTime"), cluster.getClusterTime());
mHistManager.fill(HIST("hClusterEtaPhi"), pos.Eta(), TVector2::Phi_0_2pi(pos.Phi()));
if (IndexMapPair && trackGlobalIndex) {
Expand Down
111 changes: 74 additions & 37 deletions PWGJE/Tasks/gammajettreeproducer.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,12 @@
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.

// C++ system headers first
#include <string>
#include <unordered_map>
#include <vector>

// Framework and other headers after
#include "Framework/ASoA.h"
#include "Framework/AnalysisDataModel.h"
#include "Framework/AnalysisTask.h"
Expand Down Expand Up @@ -44,9 +50,10 @@
/// \since 02.08.2024
///
using namespace o2;
using namespace o2::aod;
using namespace o2::framework;
using namespace o2::framework::expressions;
using selectedClusters = o2::soa::Filtered<o2::aod::JClusters>;
using selectedClusters = o2::soa::Filtered<o2::soa::Join<o2::aod::JClusters, o2::aod::JClusterTracks>>;

#include "Framework/runDataProcessing.h"

Expand All @@ -72,12 +79,11 @@ struct GammaJetTreeProducer {
trackSelections{"trackSelections", "globalTracks", "set track selections"};
Configurable<float> trackMinPt{"trackMinPt", 0.15, "minimum track pT cut"};
Configurable<float> jetPtMin{"jetPtMin", 5.0, "minimum jet pT cut"};
Configurable<float> jetR{"jetR", 0.4, "jet resolution parameter"};
Configurable<float> isoR{"isoR", 0.4, "isolation cone radius"};

Configurable<float> perpConeJetR{"perpConeJetR", 0.4, "perpendicular cone radius used to calculate perp cone rho for jet"};
Configurable<float> trackMatchingEoverP{"trackMatchingEoverP", 2.0, "closest track is required to have E/p < value"};
// cluster cuts
Configurable<int> mClusterDefinition{"clusterDefinition", 10, "cluster definition to be selected, e.g. 10=kV3Default"};
// Preslice<o2::aod::JClusterTracks> perClusterMatchedTracks = o2::aod::jcluster::clusterId;

int mRunNumber = 0;
int eventSelection = -1;
Expand All @@ -98,14 +104,21 @@ struct GammaJetTreeProducer {
// create histograms
LOG(info) << "Creating histograms";

const o2Axis ptAxis{100, 0, 100, "p_{T} (GeV/c)"};
const o2Axis ptAxis{100, 0, 200, "p_{T} (GeV/c)"};
const o2Axis energyAxis{100, 0, 100, "E (GeV)"};
const o2Axis m02Axis{100, 0, 3, "m02"};

const o2Axis etaAxis{100, -1, 1, "#eta"};
const o2Axis phiAxis{100, 0, 2 * TMath::Pi(), "#phi"};
const o2Axis occupancyAxis{300, 0, 30000, "occupancy"};
mHistograms.add("clusterE", "Energy of cluster", o2HistType::kTH1F, {energyAxis});
mHistograms.add("trackPt", "pT of track", o2HistType::kTH1F, {ptAxis});
mHistograms.add("chjetPt", "pT of charged jet", o2HistType::kTH1F, {ptAxis});
mHistograms.add("chjetPtEtaPhi", "pT of charged jet", o2HistType::kTHnSparseF, {ptAxis, etaAxis, phiAxis});
mHistograms.add("chjetpt_vs_constpt", "pT of charged jet vs pT of constituents", o2HistType::kTH2F, {ptAxis, ptAxis});

// track QA THnSparse
mHistograms.add("trackPtEtaPhi", "Track QA", o2HistType::kTHnSparseF, {ptAxis, etaAxis, phiAxis});
mHistograms.add("trackPtEtaOccupancy", "Track QA vs occupancy", o2HistType::kTHnSparseF, {ptAxis, etaAxis, occupancyAxis});
}

// ---------------------
Expand Down Expand Up @@ -156,12 +169,25 @@ struct GammaJetTreeProducer {
}
return iso;
}
double ch_perp_cone_rho(const auto& cluster, aod::JetTracks const& tracks, float radius = 0.4)

void runTrackQA(const auto& collision, aod::JetTracks const& tracks)
{
for (auto track : tracks) {
if (!isTrackSelected(track)) {
continue;
}
mHistograms.fill(HIST("trackPt"), track.pt());
mHistograms.fill(HIST("trackPtEtaPhi"), track.pt(), track.eta(), track.phi());
mHistograms.fill(HIST("trackPtEtaOccupancy"), track.pt(), track.eta(), collision.trackOccupancyInTimeRange());
}
}

double ch_perp_cone_rho(const auto& object, aod::JetTracks const& tracks, float radius = 0.4)
{
double ptSumLeft = 0;
double ptSumRight = 0;

double cPhi = TVector2::Phi_0_2pi(cluster.phi());
double cPhi = TVector2::Phi_0_2pi(object.phi());

// rotate cone left by 90 degrees
float cPhiLeft = cPhi - TMath::Pi() / 2;
Expand All @@ -173,8 +199,8 @@ struct GammaJetTreeProducer {
if (!isTrackSelected(track)) {
continue;
}
dRLeft = jetutilities::deltaR(cluster.eta(), cPhiLeft, track.eta(), track.phi());
dRRight = jetutilities::deltaR(cluster.eta(), cPhiRight, track.eta(), track.phi());
dRLeft = jetutilities::deltaR(object.eta(), cPhiLeft, track.eta(), track.phi());
dRRight = jetutilities::deltaR(object.eta(), cPhiRight, track.eta(), track.phi());

if (dRLeft < radius) {
ptSumLeft += track.pt();
Expand All @@ -201,16 +227,21 @@ struct GammaJetTreeProducer {
// sadly passing of the string at runtime is not possible for technical region so cluster definition is
// an integer instead
Filter clusterDefinitionSelection = (o2::aod::jcluster::definition == mClusterDefinition);
PresliceUnsorted<aod::JEMCTracks> EMCTrackPerTrack = aod::jemctrack::trackId;

// Process clusters
void processClusters(soa::Join<aod::JetCollisions, aod::BkgChargedRhos, aod::JCollisionBCs>::iterator const& collision, selectedClusters const& clusters, aod::JetTracks const& tracks)
void processClusters(soa::Join<aod::JetCollisions, aod::BkgChargedRhos, aod::JCollisionBCs>::iterator const& collision, selectedClusters const& clusters, aod::JetTracks const& tracks, aod::JEMCTracks const& emctracks)
{
if (!isEventAccepted(collision)) {
return;
}

eventsTable(collision.multiplicity(), collision.centrality(), collision.rho(), collision.eventSel(), collision.alias_raw());
eventsTable(collision.multiplicity(), collision.centrality(), collision.rho(), collision.eventSel(), collision.trackOccupancyInTimeRange(), collision.alias_raw());
collisionMapping[collision.globalIndex()] = eventsTable.lastIndex();

// loop over tracks one time for QA
runTrackQA(collision, tracks);

// loop over clusters
for (auto cluster : clusters) {

Expand All @@ -226,26 +257,24 @@ struct GammaJetTreeProducer {
// double dRMin = 100;
double p = -1;

// auto tracksofcluster = matchedtracks.sliceBy(perClusterMatchedTracks, cluster.globalIndex());
// for (const auto& match : tracksofcluster) {
// // ask the jtracks table for track with ID trackID
// double dR = deltaR(cluster.eta(), cluster.phi(), match.tracks_as<o2::aod::JTracks>().Eta(), match.tracks_as<o2::aod::JTracks>().Phi());
// if (dR < dRMin) {
// dRMin = dR;
// dEta = cluster.eta() - match.tracks_as<o2::aod::JTracks>().eta();
// dPhi = TVector2::Phi_0_2pi(cluster.phi()) - TVector2::Phi_0_2pi(match.tracks_as<o2::aod::JTracks>().phi());
// if (abs(dPhi) > M_PI) {
// dPhi = 2 * M_PI - abs(dPhi);
// }
// p = match.tracks_as<o2::aod::JTracks>().p();
// }
// }

// // for compression reasons make dPhi and dEta 0 if no match is found
// if (p == -1) {
// dPhi = 0;
// dEta = 0;
// }
// do track matching
auto tracksofcluster = cluster.matchedTracks_as<aod::JetTracks>();
for (auto track : tracksofcluster) {
if (!isTrackSelected(track)) {
continue;
}
auto emcTracksPerTrack = emctracks.sliceBy(EMCTrackPerTrack, track.globalIndex());
auto emcTrack = emcTracksPerTrack.iteratorAt(0);
// find closest track that still has E/p < trackMatchingEoverP
if (cluster.energy() / track.p() > trackMatchingEoverP) {
continue;
} else {
dEta = cluster.eta() - emcTrack.etaEmcal();
dPhi = RecoDecay::constrainAngle(RecoDecay::constrainAngle(emcTrack.phiEmcal(), -M_PI) - RecoDecay::constrainAngle(cluster.phi(), -M_PI), -M_PI);
p = track.p();
break;
}
}

gammasTable(eventsTable.lastIndex(), cluster.energy(), cluster.eta(), cluster.phi(), cluster.m02(), cluster.m20(), cluster.nCells(), cluster.time(), cluster.isExotic(), cluster.distanceToBadChannel(), cluster.nlm(), isoraw, perpconerho, dPhi, dEta, p);
}
Expand All @@ -257,30 +286,38 @@ struct GammaJetTreeProducer {
}
PROCESS_SWITCH(GammaJetTreeProducer, processClusters, "Process EMCal clusters", true);

Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * 100.0f);
Filter jetCuts = aod::jet::pt > jetPtMin;
// Process charged jets
void processChargedJets(soa::Join<aod::JetCollisions, aod::BkgChargedRhos, aod::JCollisionBCs>::iterator const& collision, soa::Filtered<soa::Join<aod::ChargedJets, aod::ChargedJetConstituents>> const& chargedJets, aod::JetTracks const&)
void processChargedJets(soa::Join<aod::JetCollisions, aod::BkgChargedRhos, aod::JCollisionBCs>::iterator const& collision, soa::Filtered<soa::Join<aod::ChargedJets, aod::ChargedJetConstituents>> const& chargedJets, aod::JetTracks const& tracks)
{
// event selection
if (!isEventAccepted(collision)) {
return;
}

float leadingTrackPt = 0;
ushort nconst = 0;
// loop over charged jets
for (auto jet : chargedJets) {
if (jet.pt() < jetPtMin)
continue;
ushort nconst = 0;
nconst = 0;
leadingTrackPt = 0;
// loop over constituents
for (auto& constituent : jet.template tracks_as<aod::JetTracks>()) {
mHistograms.fill(HIST("chjetpt_vs_constpt"), jet.pt(), constituent.pt());
nconst++;
if (constituent.pt() > leadingTrackPt) {
leadingTrackPt = constituent.pt();
}
}
int32_t storedColIndex = -1;
if (auto foundCol = collisionMapping.find(collision.globalIndex()); foundCol != collisionMapping.end()) {
storedColIndex = foundCol->second;
}
chargedJetsTable(storedColIndex, jet.pt(), jet.eta(), jet.phi(), jet.energy(), jet.mass(), jet.area(), nconst);
// calculate perp cone rho
double perpconerho = ch_perp_cone_rho(jet, tracks, perpConeJetR);
mHistograms.fill(HIST("chjetPtEtaPhi"), jet.pt(), jet.eta(), jet.phi());
chargedJetsTable(storedColIndex, jet.pt(), jet.eta(), jet.phi(), jet.r(), jet.energy(), jet.mass(), jet.area(), leadingTrackPt, perpconerho, nconst);
// fill histograms
mHistograms.fill(HIST("chjetPt"), jet.pt());
}
Expand Down
Loading