-
Notifications
You must be signed in to change notification settings - Fork 632
[PWGJE] Adding D0+JET correlation task #14795
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.
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
Open
MattOckleton
wants to merge
10
commits into
AliceO2Group:master
Choose a base branch
from
MattOckleton:master
base: master
Could not load branches
Branch not found: {{ refName }}
Loading
Could not load tags
Nothing to show
Loading
Are you sure you want to change the base?
Some commits from the old base branch may be removed from the timeline,
and old review comments may become outdated.
Open
Changes from all commits
Commits
Show all changes
10 commits
Select commit
Hold shift + click to select a range
507a196
Adding jetCorrelationD0 task to O2Physics
MattOckleton 08fc8d7
Formatting changes
MattOckleton e349acd
Formatting changes
MattOckleton bd2e3c3
Collision index dev
MattOckleton 9c73715
Table idx fixes in jetCorrelationD0 task
MattOckleton 4178945
jetCorrelationD0 indexing updates
MattOckleton 84d4b7d
Formatting changes jetCorrelationD0
MattOckleton 985f17a
Removing processMCMatched in jetCorrelationD0
MattOckleton 7d1f09d
Tidying up table definitions
MattOckleton d53cf9c
Formatting
MattOckleton File filter
Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,327 @@ | ||
| // 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. | ||
| // | ||
| /// \file jetCorrelationD0.cxx | ||
| /// \brief Task for analysing D0 triggered jet events. | ||
| /// \author Matthew Ockleton <matthew.ockleton@cern.ch>, University of Liverpool | ||
|
|
||
| #include "PWGHF/Core/DecayChannels.h" | ||
| #include "PWGHF/DataModel/AliasTables.h" | ||
| #include "PWGJE/Core/JetDerivedDataUtilities.h" | ||
| #include "PWGJE/Core/JetUtilities.h" | ||
| #include "PWGJE/DataModel/Jet.h" | ||
| #include "PWGJE/DataModel/JetReducedData.h" | ||
|
|
||
| #include "Common/Core/RecoDecay.h" | ||
|
|
||
| #include "Framework/AnalysisDataModel.h" | ||
| #include "Framework/AnalysisTask.h" | ||
| #include "Framework/HistogramRegistry.h" | ||
| #include "Framework/Logger.h" | ||
| #include "Framework/runDataProcessing.h" | ||
| #include <Framework/ASoA.h> | ||
| #include <Framework/HistogramSpec.h> | ||
|
|
||
| #include <fairlogger/Logger.h> | ||
|
|
||
| #include <string> | ||
| #include <vector> | ||
|
|
||
| using namespace o2; | ||
| using namespace o2::framework; | ||
| using namespace o2::framework::expressions; | ||
|
|
||
| namespace o2::aod | ||
| { | ||
| DECLARE_SOA_TABLE(CollisionTables, "AOD", "COLLISIONINFOTABLE", | ||
| o2::soa::Index<>, | ||
| collision::PosZ); | ||
|
|
||
| namespace collisionInfo | ||
| { | ||
| DECLARE_SOA_INDEX_COLUMN(CollisionTable, collisionTable); | ||
| } // namespace collisionInfo | ||
| namespace d0Info | ||
| { | ||
| // D0 | ||
| DECLARE_SOA_COLUMN(D0PromptBDT, d0PromptBDT, float); | ||
| DECLARE_SOA_COLUMN(D0NonPromptBDT, d0NonPromptBDT, float); | ||
| DECLARE_SOA_COLUMN(D0BkgBDT, d0BkgBDT, float); | ||
| DECLARE_SOA_COLUMN(D0M, d0M, float); | ||
| DECLARE_SOA_COLUMN(D0Pt, d0Pt, float); | ||
| DECLARE_SOA_COLUMN(D0Eta, d0Eta, float); | ||
| DECLARE_SOA_COLUMN(D0Phi, d0Phi, float); | ||
| DECLARE_SOA_COLUMN(D0McOrigin, d0McOrigin, float); | ||
| DECLARE_SOA_COLUMN(D0MD, d0MD, float); | ||
| DECLARE_SOA_COLUMN(D0PtD, d0PtD, float); | ||
| DECLARE_SOA_COLUMN(D0EtaD, d0EtaD, float); | ||
| DECLARE_SOA_COLUMN(D0PhiD, d0PhiD, float); | ||
| DECLARE_SOA_COLUMN(D0Reflection, d0Reflection, int); | ||
| } // namespace d0Info | ||
|
|
||
| DECLARE_SOA_TABLE(D0DataTables, "AOD", "D0DATATABLE", | ||
| o2::soa::Index<>, | ||
| collisionInfo::CollisionTableId, | ||
| d0Info::D0PromptBDT, | ||
| d0Info::D0NonPromptBDT, | ||
| d0Info::D0BkgBDT, | ||
| d0Info::D0M, | ||
| d0Info::D0Pt, | ||
| d0Info::D0Eta, | ||
| d0Info::D0Phi); | ||
|
|
||
| DECLARE_SOA_TABLE(D0McPTables, "AOD", "D0MCPARTICLELEVELTABLE", | ||
| o2::soa::Index<>, | ||
| collisionInfo::CollisionTableId, | ||
| d0Info::D0McOrigin, | ||
| d0Info::D0Pt, | ||
| d0Info::D0Eta, | ||
| d0Info::D0Phi); | ||
|
|
||
| namespace jetInfo | ||
| { | ||
| // D0 tables | ||
| DECLARE_SOA_INDEX_COLUMN(D0DataTable, d0Data); | ||
| DECLARE_SOA_INDEX_COLUMN(D0McPTable, d0MCP); | ||
| // Jet | ||
| DECLARE_SOA_COLUMN(JetPt, jetPt, float); | ||
| DECLARE_SOA_COLUMN(JetEta, jetEta, float); | ||
| DECLARE_SOA_COLUMN(JetPhi, jetPhi, float); | ||
| // D0-jet | ||
| DECLARE_SOA_COLUMN(D0JetDeltaPhi, d0JetDeltaPhi, float); | ||
| } // namespace jetInfo | ||
|
|
||
| DECLARE_SOA_TABLE_STAGED(JetDataTables, "JETDATATABLE", | ||
| o2::soa::Index<>, | ||
| collisionInfo::CollisionTableId, | ||
| jetInfo::D0DataTableId, | ||
| jetInfo::JetPt, | ||
| jetInfo::JetEta, | ||
| jetInfo::JetPhi, | ||
| jetInfo::D0JetDeltaPhi); | ||
|
|
||
| DECLARE_SOA_TABLE_STAGED(JetMCPTables, "JETMCPARTICLELEVELTABLE", | ||
| o2::soa::Index<>, | ||
| collisionInfo::CollisionTableId, | ||
| jetInfo::D0McPTableId, | ||
| jetInfo::JetPt, | ||
| jetInfo::JetEta, | ||
| jetInfo::JetPhi, | ||
| jetInfo::D0JetDeltaPhi); | ||
|
|
||
| } // namespace o2::aod | ||
|
|
||
| struct JetCorrelationD0 { | ||
| // Define new table | ||
| Produces<aod::CollisionTables> tableCollision; | ||
| Produces<aod::D0DataTables> tableD0; | ||
| Produces<aod::D0McPTables> tableD0MCParticle; | ||
| Produces<aod::JetDataTables> tableJet; | ||
| Produces<aod::JetMCPTables> tableJetMCParticle; | ||
|
|
||
| using TracksSelQuality = soa::Join<aod::TracksExtra, aod::TracksWMc>; | ||
|
|
||
| // Configurables | ||
| Configurable<std::string> eventSelections{"eventSelections", "sel8", "choose event selection"}; | ||
| // Configurable<std::string> 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<float> jetPtCutMin{"jetPtCutMin", 5.0, "minimum value of jet pt"}; | ||
| Configurable<float> d0PtCutMin{"d0PtCutMin", 3.0, "minimum value of d0 pt"}; | ||
| Configurable<bool> doSumw{"doSumw", false, "enable sumw2 for weighted histograms"}; | ||
| Configurable<float> vertexZCut{"vertexZCut", 10.0f, "Accepted z-vertex range"}; | ||
| Configurable<float> pTHatExponent{"pTHatExponent", 6.0, "exponent of the event weight for the calculation of pTHat"}; | ||
| Configurable<float> pTHatMaxMCD{"pTHatMaxMCD", 999.0, "maximum fraction of hard scattering for jet acceptance in detector MC"}; | ||
| Configurable<float> pTHatMaxMCP{"pTHatMaxMCP", 999.0, "maximum fraction of hard scattering for jet acceptance in particle MC"}; | ||
| Configurable<float> pTHatAbsoluteMin{"pTHatAbsoluteMin", -99.0, "minimum value of pTHat"}; | ||
|
|
||
| // Filters | ||
| Filter eventCuts = (nabs(aod::jcollision::posZ) < vertexZCut); | ||
| std::vector<int> eventSelectionBits; | ||
|
|
||
| // Histograms | ||
| HistogramRegistry registry{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; | ||
| template <typename T, typename U> | ||
| void fillD0Histograms(T const& d0, U const& scores) | ||
| { | ||
| registry.fill(HIST("hD0MlBkg"), scores[0]); | ||
| registry.fill(HIST("hD0MlNonPrompt"), scores[1]); | ||
| registry.fill(HIST("hD0MlPrompt"), scores[2]); | ||
|
|
||
| registry.fill(HIST("hD0Pt"), d0.pt()); | ||
| registry.fill(HIST("hD0M"), d0.m()); | ||
| registry.fill(HIST("hD0Eta"), d0.eta()); | ||
| registry.fill(HIST("hD0Phi"), d0.phi()); | ||
| } | ||
| template <typename T, typename U> | ||
| void fillJetHistograms(T const& jet, U const& dphi) | ||
| { | ||
| registry.fill(HIST("hJetPt"), jet.pt()); | ||
| registry.fill(HIST("hJetEta"), jet.eta()); | ||
| registry.fill(HIST("hJetPhi"), jet.phi()); | ||
| registry.fill(HIST("hJet3D"), jet.pt(), jet.eta(), jet.phi()); | ||
| registry.fill(HIST("h_Jet_pT_D0_Jet_dPhi"), jet.pt(), dphi); | ||
| } | ||
|
|
||
| template <typename T, typename U> | ||
| void applyCollisionSelections(T const& collision, U const& eventSelectionBits) | ||
| { | ||
| registry.fill(HIST("hCollisions"), 0.5); // All collisions | ||
| if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { | ||
| return; | ||
| } | ||
| registry.fill(HIST("hCollisions"), 1.5); // Selected collisions | ||
| registry.fill(HIST("hZvtxSelected"), collision.posZ()); | ||
| } | ||
|
|
||
| template <typename T, typename U> | ||
| // Jetbase is an MCD jet. We then loop through jettagv(MCP jets) to test if they match | ||
| // void fillMatchedHistograms(T const& jetBase, float weight = 1.0) // float leadingTrackPtBase, | ||
| void fillMatchedHistograms(T const& jetsBase, U const&, float weight = 1.0, float rho = 0.0) | ||
| { | ||
| for (const auto& jetBase : jetsBase) { | ||
| if (jetBase.has_matchedJetGeo()) { // geometric matching | ||
| for (auto& jetTag : jetBase.template matchedJetGeo_as<std::decay_t<U>>()) { | ||
| registry.fill(HIST("hPtMatched"), jetBase.pt() - (rho * jetBase.area()), jetTag.pt(), weight); | ||
| registry.fill(HIST("hPtMatched1d"), jetTag.pt(), weight); | ||
| registry.fill(HIST("hPhiMatched"), jetBase.phi(), jetTag.phi(), weight); | ||
| registry.fill(HIST("hEtaMatched"), jetBase.eta(), jetTag.eta(), weight); | ||
| registry.fill(HIST("hPtResolution"), jetTag.pt(), (jetTag.pt() - (jetBase.pt() - (rho * jetBase.area()))) / jetTag.pt(), weight); | ||
| registry.fill(HIST("hPhiResolution"), jetTag.pt(), jetTag.phi() - jetBase.phi(), weight); | ||
| registry.fill(HIST("hEtaResolution"), jetTag.pt(), jetTag.eta() - jetBase.eta(), weight); | ||
| } | ||
| } | ||
| } | ||
| } | ||
|
|
||
| void init(InitContext const&) | ||
| { | ||
| eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast<std::string>(eventSelections)); | ||
|
|
||
| // General Axes | ||
| AxisSpec axisEta = {100, -1.0, 1.0, "#eta"}; | ||
| AxisSpec axisPhi = {100, 0.0, o2::constants::math::TwoPI, "#phi"}; | ||
| AxisSpec axisInvMass = {500, 0, 10, "M (GeV/c)"}; | ||
|
|
||
| // General Histograms | ||
| registry.add("hCollisions", "event status;event status;entries", {HistType::kTH1F, {{4, 0.0, 4.0}}}); | ||
| registry.add("hZvtxSelected", "Z vertex position;Z_{vtx};entries", {HistType::kTH1F, {{80, -20, 20}}}); | ||
|
|
||
| // D0 Histograms | ||
| registry.add("hD0MlPrompt", "D0 ML Prompt Scores", {HistType::kTH1F, {{100, -1.0, 2.0}}}); | ||
| registry.add("hD0MlNonPrompt", "D0 ML NonPrompt Scores", {HistType::kTH1F, {{100, -1.0, 2.0}}}); | ||
| registry.add("hD0MlBkg", "D0 ML Background Scores", {HistType::kTH1F, {{100, -1.0, 2.0}}}); | ||
|
|
||
| registry.add("hD0Pt", "D^{0} p_{T};p_{T}^{D^{0}} (GeV/c);entries", {HistType::kTH1F, {{500, -100, 400, "p_{T}^{D^{0}} (GeV/c)"}}}); | ||
| registry.add("hD0M", "D^{0} Mass;M_{#pi K} (GeV/c);entries", HistType::kTH1F, {axisInvMass}); | ||
| registry.add("hD0Eta", "D^{0} #eta ;#eta_{D^{0}};entries", HistType::kTH1F, {axisEta}); | ||
| registry.add("hD0Phi", "D^{0} #phi ;#phi_{D^{0}};entries", HistType::kTH1F, {axisPhi}); | ||
|
|
||
| // Jet Histograms | ||
| registry.add("hJetPt", "jet p_{T};p_{T,jet};entries", {HistType::kTH1F, {{500, -100, 400}}}); | ||
| registry.add("hJetEta", "jet #eta;#eta_{jet};entries", HistType::kTH1F, {axisEta}); | ||
| registry.add("hJetPhi", "jet #phi;#phi_{jet};entries", HistType::kTH1F, {axisPhi}); | ||
| registry.add("hJet3D", "3D jet distribution;p_{T};#eta;#phi", {HistType::kTH3F, {{500, -100, 400}, {100, -1.0, 1.0}, {100, 0.0, o2::constants::math::TwoPI}}}); | ||
| registry.add("h_Jet_pT_D0_Jet_dPhi", "p_{T, jet} vs #Delta #phi _{D^{0}, jet}", kTH2F, {{100, 0, 100}, {100, 0, o2::constants::math::TwoPI}}); | ||
|
|
||
| // Matching histograms | ||
| registry.add("hPtMatched", "p_{T} matching;p_{T,det};p_{T,part}", {HistType::kTH2F, {{500, -100, 400}, {400, 0, 400}}}, doSumw); | ||
| registry.add("hPtMatched1d", "p_{T} matching 1d;p_{T,part}", {HistType::kTH1F, {{400, 0, 400}}}, doSumw); | ||
| registry.add("hPhiMatched", "#phi matching;#phi_{det};#phi_{part}", {HistType::kTH2F, {{100, 0.0, o2::constants::math::TwoPI}, {100, 0.0, o2::constants::math::TwoPI}}}, doSumw); | ||
| registry.add("hEtaMatched", "#eta matching;#eta_{det};#eta_{part}", {HistType::kTH2F, {{100, -1, 1}, {100, -1, 1}}}, doSumw); | ||
| registry.add("hPtResolution", "p_{T} resolution;p_{T,part};Relative Resolution", {HistType::kTH2F, {{400, 0, 400}, {1000, -5.0, 5.0}}}, doSumw); | ||
| registry.add("hPhiResolution", "#phi resolution;#p_{T,part};Resolution", {HistType::kTH2F, {{400, 0, 400}, {1000, -7.0, 7.0}}}, doSumw); | ||
| registry.add("hEtaResolution", "#eta resolution;#p_{T,part};Resolution", {HistType::kTH2F, {{400, 0, 400}, {1000, -1.0, 1.0}}}, doSumw); | ||
| } | ||
|
|
||
| void processData(soa::Filtered<aod::JetCollisions>::iterator const& collision, | ||
| aod::CandidatesD0Data const& d0Candidates, | ||
| soa::Join<aod::ChargedJets, aod::ChargedJetConstituents> const& jets) | ||
| { | ||
| applyCollisionSelections(collision, eventSelectionBits); | ||
| tableCollision(collision.posZ()); | ||
|
|
||
| for (const auto& d0Candidate : d0Candidates) { | ||
| const auto scores = d0Candidate.mlScores(); | ||
| if (d0Candidate.pt() < d0PtCutMin) { | ||
| continue; | ||
| } | ||
| fillD0Histograms(d0Candidate, scores); | ||
| tableD0(tableCollision.lastIndex(), | ||
| scores[2], | ||
| scores[1], | ||
| scores[0], | ||
| d0Candidate.m(), | ||
| d0Candidate.pt(), | ||
| d0Candidate.eta(), | ||
| d0Candidate.phi()); | ||
| for (const auto& jet : jets) { | ||
| if (jet.pt() < jetPtCutMin) { | ||
| continue; | ||
| } | ||
| float dphi = RecoDecay::constrainAngle(jet.phi() - d0Candidate.phi()); | ||
| if (abs(dphi - M_PI) > (M_PI / 2)) { | ||
| continue; | ||
| } | ||
| fillJetHistograms(jet, dphi); | ||
| tableJet(tableCollision.lastIndex(), | ||
| tableD0.lastIndex(), | ||
| jet.pt(), | ||
| jet.eta(), | ||
| jet.phi(), | ||
| dphi); | ||
| } | ||
| } | ||
| } | ||
| PROCESS_SWITCH(JetCorrelationD0, processData, "charged particle level jet analysis", true); | ||
|
|
||
| void processMCParticle(aod::JetMcCollision const& collision, | ||
| aod::CandidatesD0MCP const& d0MCPCandidates, | ||
| soa::Filtered<soa::Join<aod::ChargedMCParticleLevelJets, aod::ChargedMCParticleLevelJetConstituents>> const& jets) | ||
| { | ||
| // applyCollisionSelections(collision, eventSelectionBits); | ||
| // Need to select on pz | ||
| tableCollision(collision.posZ()); | ||
|
|
||
| for (const auto& d0MCPCandidate : d0MCPCandidates) { | ||
| if (d0MCPCandidate.pt() < d0PtCutMin) { | ||
| continue; | ||
| } | ||
| tableD0MCParticle(tableCollision.lastIndex(), | ||
| d0MCPCandidate.originMcGen(), | ||
| d0MCPCandidate.pt(), | ||
| d0MCPCandidate.eta(), | ||
| d0MCPCandidate.phi()); | ||
|
|
||
| for (const auto& jet : jets) { | ||
| if (jet.pt() < jetPtCutMin) { | ||
| continue; | ||
| } | ||
| float dphi = RecoDecay::constrainAngle(jet.phi() - d0MCPCandidate.phi()); | ||
| if (abs(dphi - M_PI) > (M_PI / 2)) { | ||
| continue; | ||
| } | ||
| fillJetHistograms(jet, dphi); | ||
| tableJetMCParticle(tableCollision.lastIndex(), | ||
| tableD0MCParticle.lastIndex(), | ||
| jet.pt(), | ||
| jet.eta(), | ||
| jet.phi(), | ||
| dphi); | ||
| } | ||
| } | ||
| } | ||
| PROCESS_SWITCH(JetCorrelationD0, processMCParticle, "process MC Particle jets", false); | ||
| }; | ||
|
|
||
| WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) | ||
| { | ||
| return WorkflowSpec{ | ||
| adaptAnalysisTask<JetCorrelationD0>(cfgc)}; | ||
| } | ||
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
does this need to be constrained in [-pi,pi] or is [0-2pi] as it is now ok?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Constrained now with abs(dphi - pi) > pi /2