diff --git a/PWGJE/Core/JetUtilities.h b/PWGJE/Core/JetUtilities.h index a681d76bf89..dac55dfbaa1 100644 --- a/PWGJE/Core/JetUtilities.h +++ b/PWGJE/Core/JetUtilities.h @@ -141,12 +141,20 @@ std::tuple>, std::vector>> MatchCl template float deltaR(T const& A, U const& B) { - float dPhi = RecoDecay::constrainAngle(A.phi() - B.phi(), -M_PI); + float dPhi = RecoDecay::constrainAngle(RecoDecay::constrainAngle(A.phi(), -M_PI) - RecoDecay::constrainAngle(B.phi(), -M_PI), -M_PI); float dEta = A.eta() - B.eta(); - return TMath::Sqrt(dEta * dEta + dPhi * dPhi); + return std::sqrt(dEta * dEta + dPhi * dPhi); } +// same as deltaR but explicit specification of the eta and phi components +template +float deltaR(T const& eta1, U const& phi1, V const& eta2, W const& phi2) +{ + float dPhi = RecoDecay::constrainAngle(RecoDecay::constrainAngle(phi1, -M_PI) - RecoDecay::constrainAngle(phi2, -M_PI), -M_PI); + float dEta = eta1 - eta2; + return std::sqrt(dEta * dEta + dPhi * dPhi); +} }; // namespace jetutilities #endif // PWGJE_CORE_JETUTILITIES_H_ diff --git a/PWGJE/DataModel/GammaJetAnalysisTree.h b/PWGJE/DataModel/GammaJetAnalysisTree.h new file mode 100644 index 00000000000..beb927f381c --- /dev/null +++ b/PWGJE/DataModel/GammaJetAnalysisTree.h @@ -0,0 +1,75 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// +/// \brief Table definitions for gamma-jet analyses +/// +/// \author Florian Jonas + +#ifndef PWGJE_DATAMODEL_GAMMAJETANALYSISTREE_H_ +#define PWGJE_DATAMODEL_GAMMAJETANALYSISTREE_H_ + +#include "Framework/AnalysisDataModel.h" +#include "PWGJE/DataModel/EMCALClusters.h" +#include "PWGJE/Core/JetDerivedDataUtilities.h" +#include "PWGJE/DataModel/Jet.h" + +namespace o2::aod +{ + +namespace gjevent +{ // TODO add rho //! event index +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_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) + +using GjEvent = GjEvents::iterator; + +namespace gjgamma +{ +DECLARE_SOA_INDEX_COLUMN(GjEvent, gjevent); //! event index +DECLARE_SOA_COLUMN(Energy, energy, float); //! cluster energy (GeV) +DECLARE_SOA_COLUMN(Eta, eta, float); //! cluster pseudorapidity (calculated using vertex) +DECLARE_SOA_COLUMN(Phi, phi, float); //! cluster azimuthal angle (calculated using vertex) +DECLARE_SOA_COLUMN(M02, m02, float); //! shower shape long axis +DECLARE_SOA_COLUMN(M20, m20, float); //! shower shape short axis +DECLARE_SOA_COLUMN(NCells, nCells, ushort); //! number of cells in cluster +DECLARE_SOA_COLUMN(Time, time, float); //! cluster time (ns) +DECLARE_SOA_COLUMN(IsExotic, isExotic, bool); //! flag to mark cluster as exotic +DECLARE_SOA_COLUMN(DistanceToBadChannel, distanceToBadChannel, float); //! distance to bad channel +DECLARE_SOA_COLUMN(NLM, nlm, ushort); //! number of local maxima +DECLARE_SOA_COLUMN(IsoRaw, isoraw, ushort); //! isolation in cone not corrected for Rho +DECLARE_SOA_COLUMN(PerpConeRho, perpconerho, float); //! rho in perpendicular cone +DECLARE_SOA_COLUMN(TMdeltaPhi, tmdeltaphi, float); //! delta phi between cluster and closest match +DECLARE_SOA_COLUMN(TMdeltaEta, tmdeltaeta, float); //! delta eta between cluster and closest match +DECLARE_SOA_COLUMN(TMtrackP, tmtrackp, float); //! track momentum of closest match, -1 if no match found +} // namespace gjgamma +DECLARE_SOA_TABLE(GjGammas, "AOD", "GJGAMMA", + gjgamma::GjEventId, gjgamma::Energy, gjgamma::Eta, gjgamma::Phi, gjgamma::M02, gjgamma::M20, gjgamma::NCells, gjgamma::Time, gjgamma::IsExotic, gjgamma::DistanceToBadChannel, gjgamma::NLM, gjgamma::IsoRaw, gjgamma::PerpConeRho, gjgamma::TMdeltaPhi, gjgamma::TMdeltaEta, gjgamma::TMtrackP) +namespace gjchjet +{ +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(Energy, energy, float); +DECLARE_SOA_COLUMN(Mass, mass, float); +DECLARE_SOA_COLUMN(Area, area, 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) +} // namespace o2::aod + +#endif // PWGJE_DATAMODEL_GAMMAJETANALYSISTREE_H_ diff --git a/PWGJE/Tasks/CMakeLists.txt b/PWGJE/Tasks/CMakeLists.txt index adbd4e141a5..a9ac88c8f28 100644 --- a/PWGJE/Tasks/CMakeLists.txt +++ b/PWGJE/Tasks/CMakeLists.txt @@ -168,6 +168,10 @@ if(FastJet_FOUND) SOURCES fulljetspectrapp.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(gamma-jet-tree-producer + SOURCES gammajettreeproducer.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::PWGJECore O2Physics::AnalysisCore + COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(bjet-tagging-ml SOURCES bjetTaggingML.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore O2Physics::MLCore diff --git a/PWGJE/Tasks/gammajettreeproducer.cxx b/PWGJE/Tasks/gammajettreeproducer.cxx new file mode 100644 index 00000000000..d721ca725df --- /dev/null +++ b/PWGJE/Tasks/gammajettreeproducer.cxx @@ -0,0 +1,295 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "Framework/ASoA.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" + +#include "Common/Core/RecoDecay.h" +#include "Common/Core/TrackSelection.h" +#include "Common/Core/TrackSelectionDefaults.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "PWGJE/Core/FastJetUtilities.h" +#include "PWGJE/Core/JetDerivedDataUtilities.h" +#include "PWGJE/Core/JetUtilities.h" +#include "PWGJE/DataModel/Jet.h" +#include "PWGJE/DataModel/GammaJetAnalysisTree.h" + +#include "EMCALBase/Geometry.h" +#include "EMCALCalib/BadChannelMap.h" +#include "PWGJE/DataModel/EMCALClusters.h" +#include "DataFormatsEMCAL/Cell.h" +#include "DataFormatsEMCAL/Constants.h" +#include "DataFormatsEMCAL/AnalysisCluster.h" +#include "TVector2.h" + +#include "CommonDataFormat/InteractionRecord.h" + +#include "EventFiltering/filterTables.h" + +// \struct GammaJetTreeProducer +/// \brief Task to produce a tree for gamma-jet analysis, including photons (and information of isolation) and charged and full jets +/// \author Florian Jonas , UC Berkeley/LBNL +/// \since 02.08.2024 +/// +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using selectedClusters = o2::soa::Filtered; + +#include "Framework/runDataProcessing.h" + +struct GammaJetTreeProducer { + // analysis tree + // charged jets + // photon candidates + Produces chargedJetsTable; + Produces eventsTable; + Produces gammasTable; + + HistogramRegistry mHistograms{"GammaJetTreeProducerHisto"}; + + // --------------- + // Configureables + // --------------- + + // event cuts + Configurable mVertexCut{"vertexCut", 10.0, "apply z-vertex cut with value in cm"}; + Configurable eventSelections{"eventSelections", "sel8", "choose event selection"}; + Configurable triggerMasks{"triggerMasks", "", "possible JE Trigger masks: fJetChLowPt,fJetChHighPt,fTrackLowPt,fTrackHighPt,fJetD0ChLowPt,fJetD0ChHighPt,fJetLcChLowPt,fJetLcChHighPt,fEMCALReadout,fJetFullHighPt,fJetFullLowPt,fJetNeutralHighPt,fJetNeutralLowPt,fGammaVeryHighPtEMCAL,fGammaVeryHighPtDCAL,fGammaHighPtEMCAL,fGammaHighPtDCAL,fGammaLowPtEMCAL,fGammaLowPtDCAL,fGammaVeryLowPtEMCAL,fGammaVeryLowPtDCAL"}; + Configurable + trackSelections{"trackSelections", "globalTracks", "set track selections"}; + Configurable trackMinPt{"trackMinPt", 0.15, "minimum track pT cut"}; + Configurable jetPtMin{"jetPtMin", 5.0, "minimum jet pT cut"}; + Configurable jetR{"jetR", 0.4, "jet resolution parameter"}; + Configurable isoR{"isoR", 0.4, "isolation cone radius"}; + + // cluster cuts + Configurable mClusterDefinition{"clusterDefinition", 10, "cluster definition to be selected, e.g. 10=kV3Default"}; + // Preslice perClusterMatchedTracks = o2::aod::jcluster::clusterId; + + int mRunNumber = 0; + int eventSelection = -1; + int trackSelection = -1; + + std::unordered_map collisionMapping; + std::vector triggerMaskBits; + + void init(InitContext const&) + { + using o2HistType = HistType; + using o2Axis = AxisSpec; + + eventSelection = jetderiveddatautilities::initialiseEventSelection(static_cast(eventSelections)); + triggerMaskBits = jetderiveddatautilities::initialiseTriggerMaskBits(triggerMasks); + trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast(trackSelections)); + + // create histograms + LOG(info) << "Creating histograms"; + + const o2Axis ptAxis{100, 0, 100, "p_{T} (GeV/c)"}; + const o2Axis energyAxis{100, 0, 100, "E (GeV)"}; + const o2Axis m02Axis{100, 0, 3, "m02"}; + + 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("chjetpt_vs_constpt", "pT of charged jet vs pT of constituents", o2HistType::kTH2F, {ptAxis, ptAxis}); + } + + // --------------------- + // Helper functions + // --------------------- + bool isTrackSelected(const auto& track) + { + if (!jetderiveddatautilities::selectTrack(track, trackSelection)) { + return false; + } + if (track.pt() < trackMinPt) { + return false; + } + + return true; + } + + bool isEventAccepted(const auto& collision) + { + + if (collision.posZ() > mVertexCut) { + return false; + } + if (!jetderiveddatautilities::selectCollision(collision, eventSelection)) { + return false; + } + if (!jetderiveddatautilities::selectTrigger(collision, triggerMaskBits)) { + return false; + } + if (!jetderiveddatautilities::eventEMCAL(collision)) { + return false; + } + return true; + } + + double ch_iso_in_cone(const auto& cluster, JetTracks const& tracks, float radius = 0.4) + { + double iso = 0; + for (auto track : tracks) { + if (!isTrackSelected(track)) { + continue; + } + // make dR function live somwhere else + float dR = jetutilities::deltaR(cluster, track); + if (dR < radius) { + iso += track.pt(); + } + } + return iso; + } + double ch_perp_cone_rho(const auto& cluster, JetTracks const& tracks, float radius = 0.4) + { + double ptSumLeft = 0; + double ptSumRight = 0; + + double cPhi = TVector2::Phi_0_2pi(cluster.phi()); + + // rotate cone left by 90 degrees + float cPhiLeft = cPhi - TMath::Pi() / 2; + float cPhiRight = cPhi + TMath::Pi() / 2; + + // loop over tracks + float dRLeft, dRRight; + for (auto track : tracks) { + if (!isTrackSelected(track)) { + continue; + } + dRLeft = jetutilities::deltaR(cluster.eta(), cPhiLeft, track.eta(), track.phi()); + dRRight = jetutilities::deltaR(cluster.eta(), cPhiRight, track.eta(), track.phi()); + + if (dRLeft < radius) { + ptSumLeft += track.pt(); + } + if (dRRight < radius) { + ptSumRight += track.pt(); + } + } + + float rho = (ptSumLeft + ptSumRight) / (2 * TMath::Pi() * radius * radius); + return rho; + } + + // --------------------- + // Processing functions + // --------------------- + void processClearMaps(JetCollisions const&) + { + collisionMapping.clear(); + } + PROCESS_SWITCH(GammaJetTreeProducer, processClearMaps, "process function that clears all the maps in each dataframe", true); + + // define cluster filter. It selects only those clusters which are of the type + // 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); + // Process clusters + void processClusters(soa::Join::iterator const& collision, selectedClusters const& clusters, JetTracks const& tracks) + { + if (!isEventAccepted(collision)) { + return; + } + + eventsTable(collision.multiplicity(), collision.centrality(), collision.rho(), collision.eventSel(), collision.alias_raw()); + collisionMapping[collision.globalIndex()] = eventsTable.lastIndex(); + + // loop over clusters + for (auto cluster : clusters) { + + // fill histograms + mHistograms.fill(HIST("clusterE"), cluster.energy()); + + double isoraw = ch_iso_in_cone(cluster, tracks, isoR); + double perpconerho = ch_perp_cone_rho(cluster, tracks, isoR); + + // find closest matched track + double dEta = 0; + double dPhi = 0; + // 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().Eta(), match.tracks_as().Phi()); + // if (dR < dRMin) { + // dRMin = dR; + // dEta = cluster.eta() - match.tracks_as().eta(); + // dPhi = TVector2::Phi_0_2pi(cluster.phi()) - TVector2::Phi_0_2pi(match.tracks_as().phi()); + // if (abs(dPhi) > M_PI) { + // dPhi = 2 * M_PI - abs(dPhi); + // } + // p = match.tracks_as().p(); + // } + // } + + // // for compression reasons make dPhi and dEta 0 if no match is found + // if (p == -1) { + // dPhi = 0; + // dEta = 0; + // } + + 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); + } + + // dummy loop over tracks + for (auto track : tracks) { + mHistograms.fill(HIST("trackPt"), track.pt()); + } + } + PROCESS_SWITCH(GammaJetTreeProducer, processClusters, "Process EMCal clusters", true); + + Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * 100.0f); + // Process charged jets + void processChargedJets(soa::Join::iterator const& collision, soa::Filtered> const& chargedJets, JetTracks const&) + { + // event selection + if (!isEventAccepted(collision)) { + return; + } + + // loop over charged jets + for (auto jet : chargedJets) { + if (jet.pt() < jetPtMin) + continue; + ushort nconst = 0; + // loop over constituents + for (auto& constituent : jet.template tracks_as()) { + mHistograms.fill(HIST("chjetpt_vs_constpt"), jet.pt(), constituent.pt()); + nconst++; + } + 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); + // fill histograms + mHistograms.fill(HIST("chjetPt"), jet.pt()); + } + } + PROCESS_SWITCH(GammaJetTreeProducer, processChargedJets, "Process charged jets", true); +}; +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + WorkflowSpec workflow{ + adaptAnalysisTask(cfgc, TaskName{"gamma-jet-tree-producer"})}; + return workflow; +}