Skip to content

Commit a941324

Browse files
author
sandeep dudi
committed
new analysis for pion and kaon using kink topology
1 parent 602db8e commit a941324

File tree

2 files changed

+313
-0
lines changed

2 files changed

+313
-0
lines changed

‎PWGLF/Tasks/Nuspex/CMakeLists.txt‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -139,4 +139,10 @@ o2physics_add_dpl_workflow(antinuclei-in-jets
139139
SOURCES antinucleiInJets.cxx
140140
PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::PWGJECore FastJet::FastJet FastJet::Contrib
141141
COMPONENT_NAME Analysis)
142+
143+
o2physics_add_dpl_workflow(kink-pika
144+
SOURCES spectraKinkPiKa.cxx
145+
PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
146+
COMPONENT_NAME Analysis)
147+
142148
endif()

‎PWGLF/Tasks/Nuspex/spectraKinkPiKa.cxx‎

Lines changed: 307 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,307 @@
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+
/// \file spectraKinkPiKa.cxx
13+
/// \brief Example of a simple task for the analysis of the muon from Kaon pion using kink topology
14+
/// \author
15+
16+
#include "PWGLF/DataModel/LFKinkDecayTables.h"
17+
18+
#include "Common/DataModel/EventSelection.h"
19+
#include "Common/DataModel/PIDResponse.h"
20+
21+
#include "Framework/AnalysisTask.h"
22+
#include "Framework/O2DatabasePDGPlugin.h"
23+
#include "Framework/runDataProcessing.h"
24+
#include "ReconstructionDataFormats/PID.h"
25+
26+
#include "Math/GenVector/Boost.h"
27+
#include "Math/Vector3D.h"
28+
#include "Math/Vector4D.h"
29+
#include "TPDGCode.h"
30+
#include "TVector3.h"
31+
#include <TMath.h>
32+
#include <TString.h>
33+
34+
#include <cstdlib>
35+
#include <vector>
36+
37+
using namespace std;
38+
using namespace o2;
39+
using namespace o2::aod;
40+
using namespace o2::framework;
41+
using namespace o2::framework::expressions;
42+
using namespace o2::constants::physics;
43+
44+
using TracksFull = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksCovIU, aod::pidTPCPi, aod::pidTPCKa, aod::pidTPCMu>;
45+
using CollisionsFull = soa::Join<aod::Collisions, aod::EvSel>;
46+
using CollisionsFullMC = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSel>;
47+
struct spectraKinkPiKa {
48+
Service<o2::framework::O2DatabasePDG> pdg;
49+
// Histograms are defined with HistogramRegistry
50+
HistogramRegistry rEventSelection{"eventSelection", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
51+
HistogramRegistry rpiKkink{"rpiKkink", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
52+
53+
// Configurable for event selection
54+
Configurable<float> cutzvertex{"cutzvertex", 10.0f, "Accepted z-vertex range (cm)"};
55+
Configurable<float> cutNSigmaPi{"cutNSigmaPi", 4, "NSigmaTPCPion"};
56+
Configurable<float> cutNSigmaKa{"cutNSigmaKa", 4, "NSigmaTPCKaon"};
57+
Configurable<int> pid{"pidMother", 321, ""};
58+
Configurable<bool> d0pid{"dopid", 0, ""};
59+
60+
Preslice<aod::KinkCands> mPerCol = aod::track::collisionId;
61+
62+
void init(InitContext const&)
63+
{
64+
// Axes
65+
const AxisSpec ptAxis{100, 0, 10, "#it{p}_{T} (GeV/#it{c})"};
66+
const AxisSpec qtAxis{2000, 0, 2, "#it{q}_{T} (GeV/#it{c})"};
67+
const AxisSpec kinkAxis{200, 0, 4, "#theta"};
68+
const AxisSpec etaAxis{200, -5.0, 5.0, "#eta"};
69+
const AxisSpec vertexZAxis{100, -15., 15., "vrtx_{Z} [cm]"};
70+
71+
// Event selection
72+
rEventSelection.add("hVertexZRec", "hVertexZRec", {HistType::kTH1F, {vertexZAxis}});
73+
74+
rpiKkink.add("h2_dau_pt_vs_eta_rec", "pt_vs_eta_dau", {HistType::kTH2F, {ptAxis, etaAxis}});
75+
rpiKkink.add("h2_moth_pt_vs_eta_rec", "pt_vs_eta_moth", {HistType::kTH2F, {ptAxis, etaAxis}});
76+
rpiKkink.add("h2_pt_moth_vs_dau_rec", "pt_moth_vs_dau", {HistType::kTH2F, {ptAxis, ptAxis}});
77+
78+
rpiKkink.add("h2_qt", "qt", {HistType::kTH1F, {qtAxis}});
79+
rpiKkink.add("h2_qt_vs_pt", "qt_pt", {HistType::kTH2F, {qtAxis, ptAxis}});
80+
81+
rpiKkink.add("h2_kink_angle", "kink angle", {HistType::kTH1F, {kinkAxis}});
82+
83+
// pion
84+
rpiKkink.add("h2_dau_pt_vs_eta_rec_pion", "pt_vs_eta_dau", {HistType::kTH2F, {ptAxis, etaAxis}});
85+
rpiKkink.add("h2_moth_pt_vs_eta_rec_pion", "pt_vs_eta_moth", {HistType::kTH2F, {ptAxis, etaAxis}});
86+
rpiKkink.add("h2_pt_moth_vs_dau_rec_pion", "pt_moth_vs_dau", {HistType::kTH2F, {ptAxis, ptAxis}});
87+
88+
rpiKkink.add("h2_qt_pion", "qt", {HistType::kTH1F, {qtAxis}});
89+
rpiKkink.add("h2_qt_vs_ptpion", "qt_pt", {HistType::kTH2F, {qtAxis, ptAxis}});
90+
rpiKkink.add("h2_kink_angle_pion", "kink angle", {HistType::kTH1F, {kinkAxis}});
91+
92+
if (doprocessMC) {
93+
rpiKkink.add("h2_dau_pt_vs_eta_gen", "pt_vs_eta_dau", {HistType::kTH2F, {ptAxis, etaAxis}});
94+
rpiKkink.add("h2_moth_pt_vs_eta_gen", "pt_vs_eta_moth", {HistType::kTH2F, {ptAxis, etaAxis}});
95+
rpiKkink.add("h2_pt_moth_vs_dau_gen", "pt_moth_vs_dau", {HistType::kTH2F, {ptAxis, ptAxis}});
96+
97+
rpiKkink.add("h2_qt_gen", "qt", {HistType::kTH1F, {qtAxis}});
98+
rpiKkink.add("h2_qt_rec", "qt", {HistType::kTH1F, {qtAxis}});
99+
rpiKkink.add("h2_kink_angle_gen", "kink angle", {HistType::kTH1F, {kinkAxis}});
100+
}
101+
}
102+
103+
void processData(CollisionsFull::iterator const& collision, aod::KinkCands const& KinkCands, TracksFull const&)
104+
{
105+
ROOT::Math::PxPyPzMVector v0;
106+
ROOT::Math::PxPyPzMVector v1;
107+
108+
if (std::abs(collision.posZ()) > cutzvertex || !collision.sel8() || !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) {
109+
return;
110+
}
111+
if (!collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder) || !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) {
112+
return;
113+
}
114+
115+
rEventSelection.fill(HIST("hVertexZRec"), collision.posZ());
116+
for (const auto& kinkCand : KinkCands) {
117+
auto dauTrack = kinkCand.trackDaug_as<TracksFull>();
118+
auto mothTrack = kinkCand.trackMoth_as<TracksFull>();
119+
bool kaon = false;
120+
bool pion = false;
121+
if (abs(mothTrack.tpcNSigmaKa()) < cutNSigmaKa) {
122+
kaon = true;
123+
}
124+
if (abs(mothTrack.tpcNSigmaPi()) < cutNSigmaPi) {
125+
pion = true;
126+
}
127+
if (!kaon && !pion)
128+
continue;
129+
v0.SetCoordinates(mothTrack.px(), mothTrack.py(), mothTrack.pz(), o2::constants::physics::MassPionCharged);
130+
v1.SetCoordinates(dauTrack.px(), dauTrack.py(), dauTrack.pz(), o2::constants::physics::MassMuon);
131+
132+
float pMoth = v0.P();
133+
float pDaug = v1.P();
134+
float spKink = mothTrack.px() * dauTrack.px() + mothTrack.py() * dauTrack.py() + mothTrack.pz() * dauTrack.pz();
135+
float kink_angle = std::acos(spKink / (pMoth * pDaug));
136+
if (kaon) {
137+
rpiKkink.fill(HIST("h2_moth_pt_vs_eta_rec"), v0.Pt(), v0.Eta());
138+
rpiKkink.fill(HIST("h2_dau_pt_vs_eta_rec"), v1.Pt(), v1.Eta());
139+
rpiKkink.fill(HIST("h2_pt_moth_vs_dau_rec"), v0.Pt(), v1.Pt());
140+
rpiKkink.fill(HIST("h2_kink_angle"), kink_angle);
141+
}
142+
if (pion) {
143+
rpiKkink.fill(HIST("h2_moth_pt_vs_eta_rec_pion"), v0.Pt(), v0.Eta());
144+
rpiKkink.fill(HIST("h2_dau_pt_vs_eta_rec_pion"), v1.Pt(), v1.Eta());
145+
rpiKkink.fill(HIST("h2_pt_moth_vs_dau_rec_pion"), v0.Pt(), v1.Pt());
146+
rpiKkink.fill(HIST("h2_kink_angle_pion"), kink_angle);
147+
}
148+
ROOT::Math::Boost boost(-v0.BoostToCM());
149+
ROOT::Math::PxPyPzMVector dBoosted = boost(v1);
150+
TVector3 p_d_rest(dBoosted.Px(), dBoosted.Py(), dBoosted.Pz()); // Daughter in mother rest frame
151+
// Compute transverse component
152+
TVector3 boostDir(v0.Px(), v0.Py(), v0.Pz());
153+
boostDir = boostDir.Unit();
154+
double pt_d = p_d_rest.Perp(boostDir); // or p_d_rest.Mag() * sin(theta)
155+
156+
if (kaon) {
157+
rpiKkink.fill(HIST("h2_qt"), pt_d);
158+
rpiKkink.fill(HIST("h2_qt_vs_pt"), pt_d, v1.Pt());
159+
}
160+
if (pion) {
161+
rpiKkink.fill(HIST("h2_qt_pion"), pt_d);
162+
rpiKkink.fill(HIST("h2_qt_vs_ptpion"), pt_d, v1.Pt());
163+
}
164+
}
165+
}
166+
PROCESS_SWITCH(spectraKinkPiKa, processData, "Data processing", true);
167+
168+
void processMC(CollisionsFullMC const& collisions, aod::KinkCands const& KinkCands, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const& particlesMC, TracksFull const&)
169+
{
170+
for (const auto& collision : collisions) {
171+
ROOT::Math::PxPyPzMVector v0;
172+
ROOT::Math::PxPyPzMVector v1;
173+
if (std::abs(collision.posZ()) > cutzvertex || !collision.sel8() || !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) {
174+
continue;
175+
}
176+
if (!collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder) || !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) {
177+
continue;
178+
}
179+
rEventSelection.fill(HIST("hVertexZRec"), collision.posZ());
180+
auto kinkCandPerColl = KinkCands.sliceBy(mPerCol, collision.globalIndex());
181+
for (const auto& kinkCand : kinkCandPerColl) {
182+
auto dauTrack = kinkCand.trackDaug_as<TracksFull>();
183+
auto mothTrack = kinkCand.trackMoth_as<TracksFull>();
184+
if (dauTrack.sign() != mothTrack.sign()) {
185+
LOG(info) << "Skipping kink candidate with opposite sign daughter and mother: " << kinkCand.globalIndex();
186+
continue; // Skip if the daughter has the opposite sign as the mother
187+
}
188+
bool kaon = false;
189+
bool pion = false;
190+
if (abs(mothTrack.tpcNSigmaKa()) < cutNSigmaKa) {
191+
kaon = true;
192+
}
193+
if (abs(mothTrack.tpcNSigmaPi()) < cutNSigmaPi) {
194+
pion = true;
195+
}
196+
if (!kaon && !pion)
197+
continue;
198+
199+
v0.SetCoordinates(mothTrack.px(), mothTrack.py(), mothTrack.pz(), o2::constants::physics::MassPionCharged);
200+
v1.SetCoordinates(dauTrack.px(), dauTrack.py(), dauTrack.pz(), o2::constants::physics::MassMuon);
201+
202+
float pMoth = v0.P();
203+
float pDaug = v1.P();
204+
float spKink = mothTrack.px() * dauTrack.px() + mothTrack.py() * dauTrack.py() + mothTrack.pz() * dauTrack.pz();
205+
float kink_angle = std::acos(spKink / (pMoth * pDaug));
206+
207+
rpiKkink.fill(HIST("h2_moth_pt_vs_eta_rec"), v0.Pt(), v0.Eta());
208+
rpiKkink.fill(HIST("h2_dau_pt_vs_eta_rec"), v1.Pt(), v1.Eta());
209+
rpiKkink.fill(HIST("h2_pt_moth_vs_dau_rec"), v0.Pt(), v1.Pt());
210+
rpiKkink.fill(HIST("h2_kink_angle"), kink_angle);
211+
212+
ROOT::Math::Boost boost(-v0.BoostToCM());
213+
ROOT::Math::PxPyPzMVector dBoosted = boost(v1);
214+
TVector3 p_d_rest(dBoosted.Px(), dBoosted.Py(), dBoosted.Pz()); // Daughter in mother rest frame
215+
// Compute transverse component
216+
TVector3 boostDir(v0.Px(), v0.Py(), v0.Pz());
217+
boostDir = boostDir.Unit();
218+
double pt_d = p_d_rest.Perp(boostDir); // or p_d_rest.Mag() * sin(theta)
219+
rpiKkink.fill(HIST("h2_qt"), pt_d);
220+
221+
// do MC association
222+
auto mcLabMoth = trackLabelsMC.rawIteratorAt(mothTrack.globalIndex());
223+
auto mcLabDau = trackLabelsMC.rawIteratorAt(dauTrack.globalIndex());
224+
if (mcLabMoth.has_mcParticle() && mcLabDau.has_mcParticle()) {
225+
auto mcTrackMoth = mcLabMoth.mcParticle_as<aod::McParticles>();
226+
auto mcTrackDau = mcLabDau.mcParticle_as<aod::McParticles>();
227+
if (!mcTrackDau.has_mothers()) {
228+
continue;
229+
}
230+
for (auto& piMother : mcTrackDau.mothers_as<aod::McParticles>()) {
231+
if (piMother.globalIndex() != mcTrackMoth.globalIndex()) {
232+
continue;
233+
}
234+
if (std::abs(mcTrackMoth.pdgCode()) != pid || std::abs(mcTrackDau.pdgCode()) != kMuonPlus) {
235+
continue;
236+
}
237+
// rpiKkink.fill(HIST("h2MassPtMCRec"), kinkCand.ptMoth(), v1.Pt());
238+
rpiKkink.fill(HIST("h2_qt_rec"), pt_d);
239+
}
240+
}
241+
}
242+
}
243+
for (const auto& mcPart : particlesMC) {
244+
ROOT::Math::PxPyPzMVector v0;
245+
ROOT::Math::PxPyPzMVector v1;
246+
247+
if (!d0pid && (std::abs(mcPart.pdgCode()) != pid || std::abs(mcPart.y()) > 0.8)) {
248+
continue;
249+
}
250+
if (d0pid && (std::abs(mcPart.pdgCode()) != kD0 || std::abs(mcPart.pdgCode()) != kDPlus || std::abs(mcPart.pdgCode()) != kDStar || std::abs(mcPart.y()) > 0.8)) {
251+
continue;
252+
}
253+
254+
if (!mcPart.has_daughters()) {
255+
continue; // Skip if no daughters
256+
}
257+
bool hasKaonpionDaughter = false;
258+
for (const auto& daughter : mcPart.daughters_as<aod::McParticles>()) {
259+
if (std::abs(daughter.pdgCode()) == kMuonPlus) { // muon PDG code
260+
hasKaonpionDaughter = true;
261+
v1.SetCoordinates(daughter.px(), daughter.py(), daughter.pz(), o2::constants::physics::MassMuon);
262+
break; // Found a muon daughter, exit loop
263+
}
264+
}
265+
if (!hasKaonpionDaughter) {
266+
continue; // Skip if no muon daughter found
267+
}
268+
if (pid == kKPlus) {
269+
v0.SetCoordinates(mcPart.px(), mcPart.py(), mcPart.pz(), o2::constants::physics::MassKaonCharged);
270+
}
271+
272+
if (pid == kPiPlus) {
273+
v0.SetCoordinates(mcPart.px(), mcPart.py(), mcPart.pz(), o2::constants::physics::MassPionCharged);
274+
}
275+
if (d0pid) {
276+
v0.SetCoordinates(mcPart.px(), mcPart.py(), mcPart.pz(), o2::constants::physics::MassD0);
277+
}
278+
279+
float pMoth = v0.P();
280+
float pDaug = v1.P();
281+
float spKink = v0.Px() * v1.Px() + v0.Py() * v1.Py() + v0.Pz() * v0.Pz();
282+
float kink_angle = std::acos(spKink / (pMoth * pDaug));
283+
284+
// std::cout<< kinkCand.ptMoth()<<" check "<<v0.Pt()<<std::endl;
285+
rpiKkink.fill(HIST("h2_moth_pt_vs_eta_gen"), v0.Pt(), v0.Eta());
286+
rpiKkink.fill(HIST("h2_dau_pt_vs_eta_gen"), v1.Pt(), v1.Eta());
287+
rpiKkink.fill(HIST("h2_pt_moth_vs_dau_gen"), v0.Pt(), v1.Pt());
288+
rpiKkink.fill(HIST("h2_kink_angle_gen"), kink_angle);
289+
290+
ROOT::Math::Boost boost(-v0.BoostToCM());
291+
ROOT::Math::PxPyPzMVector dBoosted = boost(v1);
292+
TVector3 p_d_rest(dBoosted.Px(), dBoosted.Py(), dBoosted.Pz()); // Daughter in mother rest frame
293+
// Compute transverse component
294+
TVector3 boostDir(v0.Px(), v0.Py(), v0.Pz());
295+
boostDir = boostDir.Unit();
296+
double pt_d = p_d_rest.Perp(boostDir); // or p_d_rest.Mag() * sin(theta)
297+
rpiKkink.fill(HIST("h2_qt_gen"), pt_d);
298+
}
299+
}
300+
PROCESS_SWITCH(spectraKinkPiKa, processMC, "MC processing", false);
301+
};
302+
303+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
304+
{
305+
return WorkflowSpec{
306+
adaptAnalysisTask<spectraKinkPiKa>(cfgc)};
307+
}

0 commit comments

Comments
 (0)