Skip to content

Commit e5c480b

Browse files
monamelopmonamelop
committed
[PWGJE] Add jet-ds-spectrum-subs task
1 parent 37278bf commit e5c480b

File tree

2 files changed

+233
-1
lines changed

2 files changed

+233
-1
lines changed

‎PWGJE/Tasks/CMakeLists.txt‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -358,5 +358,9 @@ if(FastJet_FOUND)
358358
o2physics_add_dpl_workflow(substructure-debug
359359
SOURCES substructureDebug.cxx
360360
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore
361+
COMPONENT_NAME Analysis)
362+
o2physics_add_dpl_workflow(jet-ds-spectrum-subs
363+
SOURCES jetDsSpectrumAndSubstructure.cxx
364+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore
361365
COMPONENT_NAME Analysis)
362-
endif()
366+
endif()

‎PWGJE/Tasks/jetDsSpectrumAndSubstructure.cxx‎

Lines changed: 228 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,228 @@
1+
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
//
5+
// This software is distributed under the terms of the GNU General Public
6+
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
//
8+
// In applying this license CERN does not waive the privileges and immunities
9+
// granted to it by virtue of its status as an Intergovernmental Organization
10+
// or submit itself to any jurisdiction.
11+
12+
// Jet substructure and spectrum task for D_s mesons
13+
//
14+
// This task is used to reconstruct and analyse jets containing charged D_s
15+
// mesons
16+
//
17+
/// \author Monalisa Melo <monalisa.melo@cern.ch>
18+
//
19+
20+
21+
#include "PWGJE/DataModel/Jet.h"
22+
#include "PWGJE/DataModel/JetSubtraction.h"
23+
#include "PWGJE/DataModel/JetReducedData.h"
24+
25+
#include "PWGJE/Core/JetDerivedDataUtilities.h"
26+
#include "PWGJE/Core/JetUtilities.h"
27+
#include "PWGHF/Core/DecayChannels.h"
28+
29+
#include "Common/Core/RecoDecay.h"
30+
#include "Common/DataModel/TrackSelectionTables.h"
31+
32+
#include "Framework/ASoA.h"
33+
#include "Framework/AnalysisDataModel.h"
34+
#include <Framework/AnalysisTask.h>
35+
#include "Framework/HistogramRegistry.h"
36+
#include <CommonConstants/MathConstants.h>
37+
#include <Framework/Configurable.h>
38+
#include <Framework/HistogramSpec.h>
39+
#include <Framework/InitContext.h>
40+
#include <Framework/ConfigContext.h>
41+
#include <Framework/DataProcessorSpec.h>
42+
#include <Framework/runDataProcessing.h>
43+
44+
45+
46+
#include "TVector3.h"
47+
#include <cmath>
48+
#include <cstdlib>
49+
#include <string>
50+
#include <vector>
51+
52+
using namespace o2;
53+
using namespace o2::framework;
54+
using namespace o2::framework::expressions;
55+
56+
namespace o2::aod
57+
{
58+
namespace jet_distance
59+
{
60+
DECLARE_SOA_COLUMN(JetHfDist, jetHfDist, float);
61+
DECLARE_SOA_COLUMN(JetPt, jetPt, float);
62+
DECLARE_SOA_COLUMN(JetEta, jetEta, float);
63+
DECLARE_SOA_COLUMN(JetPhi, jetPhi, float);
64+
DECLARE_SOA_COLUMN(JetNConst, jetNConst, int);
65+
DECLARE_SOA_COLUMN(HfPt, hfPt, float);
66+
DECLARE_SOA_COLUMN(HfEta, hfEta, float);
67+
DECLARE_SOA_COLUMN(HfPhi, hfPhi, float);
68+
DECLARE_SOA_COLUMN(HfMass, hfMass, float);
69+
DECLARE_SOA_COLUMN(HfY, hfY, float);
70+
DECLARE_SOA_COLUMN(HfMlScore0, hfMlScore0, float);
71+
DECLARE_SOA_COLUMN(HfMlScore1, hfMlScore1, float);
72+
DECLARE_SOA_COLUMN(HfMlScore2, hfMlScore2, float);
73+
} // namespace jet_distance
74+
75+
DECLARE_SOA_TABLE(JetDistanceTable, "AOD", "JETDISTTABLE",
76+
jet_distance::JetHfDist,
77+
jet_distance::JetPt,
78+
jet_distance::JetEta,
79+
jet_distance::JetPhi,
80+
jet_distance::JetNConst,
81+
jet_distance::HfPt,
82+
jet_distance::HfEta,
83+
jet_distance::HfPhi,
84+
jet_distance::HfMass,
85+
jet_distance::HfY,
86+
jet_distance::HfMlScore0,
87+
jet_distance::HfMlScore1,
88+
jet_distance::HfMlScore2);
89+
}
90+
91+
struct JetDsSpecSubs {
92+
HistogramRegistry registry{"registry",
93+
{{"h_collisions", "event status;event status;entries", {HistType::kTH1F, {{4, 0.0, 4.0}}}},
94+
{"h_track_pt", "track pT;#it{p}_{T,track} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}},
95+
{"h_track_eta", "track #eta;#eta_{track};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}},
96+
{"h_track_phi", "track #varphi;#varphi_{track};entries", {HistType::kTH1F, {{80, -1.0, 7.}}}},
97+
{"h_jet_pt", "jet pT;#it{p}_{T,jet} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}},
98+
{"h_jet_eta", "jet #eta;#eta_{jet};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}},
99+
{"h_jet_phi", "jet #phi;#phi_{jet};entries", {HistType::kTH1F, {{80, -1.0, 7.}}}},
100+
{"h_collision_counter", "# of collisions;", {HistType::kTH1F, {{200, 0., 200.}}}},
101+
{"h_jet_counter", ";# of D_{S} jets;", {HistType::kTH1F, {{6, 0., 3.0}}}},
102+
{"h_ds_jet_projection", ";z^{D_{S},jet}_{||};dN/dz^{D_{S},jet}_{||}", {HistType::kTH1F, {{1000, 0., 10.}}}},
103+
{"h_ds_jet_distance_vs_projection", ";#DeltaR_{D_{S},jet};z^{D_{S},jet}_{||}", {HistType::kTH2F, {{1000, 0., 10.}, {1000, 0., 10.}}}},
104+
{"h_ds_jet_distance", ";#DeltaR_{D_{S},jet};dN/d(#DeltaR)", {HistType::kTH1F, {{1000, 0., 10.}}}},
105+
{"h_ds_jet_pt", ";p_{T,D_{S} jet};dN/dp_{T,D_{S} jet}", {HistType::kTH1F, {{200, 0., 10.}}}},
106+
{"h_ds_jet_eta", ";#eta_{T,D_{S} jet};dN/d#eta_{D_{S} jet}", {HistType::kTH1F, {{250, -5., 5.}}}},
107+
{"h_ds_jet_phi", ";#phi_{T,D_{S} jet};dN/d#phi_{D_{S} jet}", {HistType::kTH1F, {{250, -10., 10.}}}},
108+
{"h_ds_mass", ";m_{D_{S}} (GeV/c^{2});dN/dm_{D_{S}}", {HistType::kTH1F, {{1000, 0., 10.}}}},
109+
{"h_ds_eta", ";#eta_{D_{S}} (GeV/c^{2});dN/d#eta_{D_{S}}", {HistType::kTH1F, {{250, -5., 5.}}}},
110+
{"h_ds_phi", ";#phi_{D_{S}} (GeV/c^{2});dN/d#phi_{D_{S}}", {HistType::kTH1F, {{250, -10., 10.}}}}}
111+
};
112+
113+
Configurable<float> vertexZCut{"vertexZCut", 10.0f, "Accepted z-vertex range"};
114+
115+
Configurable<float> jetPtMin{"jetPtMin", 5.0, "minimum jet pT cut"};
116+
Configurable<float> jetR{"jetR", 0.4, "jet resolution parameter"};
117+
118+
Configurable<std::string> eventSelections{"eventSelections", "sel8", "choose event selection"};
119+
Configurable<std::string> trackSelections{"trackSelections", "globalTracks", "set track selections"};
120+
121+
std::vector<int> eventSelectionBits;
122+
int trackSelection = -1;
123+
124+
Produces<aod::JetDistanceTable> distJetTable;
125+
126+
127+
void init(o2::framework::InitContext&)
128+
{
129+
eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits(static_cast<std::string>(eventSelections));
130+
trackSelection = jetderiveddatautilities::initialiseTrackSelection(static_cast<std::string>(trackSelections));
131+
}
132+
133+
134+
Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * 100.0f);
135+
Filter collisionFilter = nabs(aod::jcollision::posZ) < vertexZCut;
136+
137+
138+
void processCollisions(aod::JetCollision const& collision, aod::JetTracks const& tracks)
139+
{
140+
141+
registry.fill(HIST("h_collisions"), 0.5);
142+
if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) {
143+
return;
144+
}
145+
registry.fill(HIST("h_collisions"), 1.5);
146+
for (auto const& track : tracks) {
147+
if (!jetderiveddatautilities::selectTrack(track, trackSelection)) {
148+
continue;
149+
}
150+
registry.fill(HIST("h_track_pt"), track.pt());
151+
registry.fill(HIST("h_track_eta"), track.eta());
152+
registry.fill(HIST("h_track_phi"), track.phi());
153+
}
154+
}
155+
PROCESS_SWITCH(JetDsSpecSubs, processCollisions, "process JE collisions", false);
156+
157+
void processDataCharged(soa::Filtered<aod::JetCollisions>::iterator const& collision, soa::Filtered<aod::ChargedJets> const& jets)
158+
{
159+
if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) {
160+
return;
161+
}
162+
//jets -> charged jet
163+
for (auto& jet : jets) {
164+
registry.fill(HIST("h_jet_pt"), jet.pt());
165+
registry.fill(HIST("h_jet_eta"), jet.eta());
166+
registry.fill(HIST("h_jet_phi"), jet.phi());
167+
}
168+
}
169+
PROCESS_SWITCH(JetDsSpecSubs, processDataCharged, "charged jets in data", false);
170+
171+
void processDataChargedSubstructure(aod::JetCollision const& collision,
172+
soa::Join<aod::DsChargedJets, aod::DsChargedJetConstituents> const& jets,
173+
aod::CandidatesDsData const&,
174+
aod::JetTracks const&)
175+
{
176+
// apply event selection and fill histograms for sanity check
177+
registry.fill(HIST("h_collision_counter"), 2.0);
178+
if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) || !(std::abs(collision.posZ()) < vertexZCut)) {
179+
return;
180+
}
181+
registry.fill(HIST("h_collision_counter"), 3.0);
182+
// jets -> charged jet with Ds
183+
for (const auto& jet : jets) {
184+
//number of charged jets with Ds
185+
registry.fill(HIST("h_jet_counter"), 0.5);
186+
// obtaining jet 3-vector
187+
TVector3 jetVector(jet.px(), jet.py(), jet.pz());
188+
189+
for (const auto& dsCandidate : jet.candidates_as<aod::CandidatesDsData>()) {
190+
191+
// obtaining jet 3-vector
192+
TVector3 dsVector(dsCandidate.px(), dsCandidate.py(), dsCandidate.pz());
193+
194+
// calculating fraction of the jet momentum carried by the Ds along the direction of the jet axis
195+
double zParallel = (jetVector * dsVector) / (jetVector * jetVector);
196+
197+
// calculating angular distance in eta-phi plane
198+
double axisDistance = jetutilities::deltaR(jet, dsCandidate);
199+
200+
// filling histograms
201+
registry.fill(HIST("h_ds_jet_projection"), zParallel);
202+
registry.fill(HIST("h_ds_jet_distance_vs_projection"), axisDistance, zParallel);
203+
registry.fill(HIST("h_ds_jet_distance"), axisDistance);
204+
registry.fill(HIST("h_ds_jet_pt"), jet.pt());
205+
registry.fill(HIST("h_ds_jet_eta"), jet.eta());
206+
registry.fill(HIST("h_ds_jet_phi"), jet.phi());
207+
registry.fill(HIST("h_ds_mass"), dsCandidate.m());
208+
registry.fill(HIST("h_ds_eta"), dsCandidate.eta());
209+
registry.fill(HIST("h_ds_phi"), dsCandidate.phi());
210+
211+
// filling table
212+
distJetTable(axisDistance,
213+
jet.pt(), jet.eta(), jet.phi(), jet.tracks_as<aod::JetTracks>().size(),
214+
dsCandidate.pt(), dsCandidate.eta(), dsCandidate.phi(), dsCandidate.m(), dsCandidate.y(), dsCandidate.mlScores()[0], dsCandidate.mlScores()[1], dsCandidate.mlScores()[2]);
215+
216+
break; // get out of candidates' loop after first HF particle is found in jet
217+
} // end of DS candidates loop
218+
219+
} // end of jets loop
220+
221+
} // end of process function
222+
PROCESS_SWITCH(JetDsSpecSubs, processDataChargedSubstructure, "charged HF jet substructure", false);
223+
224+
225+
226+
};
227+
228+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask<JetDsSpecSubs>(cfgc, TaskName{"jet-ds-spectrum-subs"})}; }

0 commit comments

Comments
 (0)