Skip to content

Commit 8a35660

Browse files
skundu692skundu692
authored
[PWGLF] Task for high mas reconstruction in 2 phi decay mode (#9090)
1 parent fa29514 commit 8a35660

File tree

2 files changed

+324
-0
lines changed

2 files changed

+324
-0
lines changed

‎PWGLF/Tasks/Resonances/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -173,3 +173,8 @@ o2physics_add_dpl_workflow(chk892pp
173173
SOURCES chk892pp.cxx
174174
PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
175175
COMPONENT_NAME Analysis)
176+
177+
o2physics_add_dpl_workflow(doublephimeson
178+
SOURCES doublephimeson.cxx
179+
PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
180+
COMPONENT_NAME Analysis)

‎PWGLF/Tasks/Resonances/doublephimeson.cxx‎

Lines changed: 319 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,319 @@
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+
/// \brief this is a starting point for the Resonances tutorial
13+
/// \author sourav kundu
14+
/// \since 02/11/2023
15+
16+
#include <Framework/Configurable.h>
17+
#include <TLorentzVector.h>
18+
#include <Math/GenVector/Boost.h>
19+
#include <Math/Vector4D.h>
20+
#include <TMath.h>
21+
#include <fairlogger/Logger.h>
22+
#include <iostream>
23+
#include <iterator>
24+
#include <string>
25+
26+
#include "Framework/AnalysisTask.h"
27+
#include "Framework/ASoAHelpers.h"
28+
#include "Framework/runDataProcessing.h"
29+
#include "Framework/AnalysisDataModel.h"
30+
#include "Framework/StepTHn.h"
31+
#include "Common/Core/trackUtilities.h"
32+
#include "PWGLF/DataModel/ReducedDoublePhiTables.h"
33+
#include "CommonConstants/PhysicsConstants.h"
34+
35+
using namespace o2;
36+
using namespace o2::framework;
37+
using namespace o2::framework::expressions;
38+
using namespace o2::soa;
39+
40+
struct doublephimeson {
41+
HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
42+
43+
Configurable<bool> fillRotation{"fillRotation", 1, "Fill rotation"};
44+
Configurable<int> strategyPID{"strategyPID", 0, "PID strategy"};
45+
Configurable<float> minPhiMass{"minPhiMass", 1.01, "Minimum phi mass"};
46+
Configurable<float> maxPhiMass{"maxPhiMass", 1.03, "Maximum phi mass"};
47+
Configurable<bool> additionalEvsel{"additionalEvsel", false, "Additional event selection"};
48+
Configurable<float> cutNsigmaTPC{"cutNsigmaTPC", 2.5, "nsigma cut TPC"};
49+
Configurable<float> cutNsigmaTOF{"cutNsigmaTOF", 3.0, "nsigma cut TOF"};
50+
// Event Mixing
51+
Configurable<int> nEvtMixing{"nEvtMixing", 10, "Number of events to mix"};
52+
ConfigurableAxis CfgVtxBins{"CfgVtxBins", {10, -10, 10}, "Mixing bins - z-vertex"};
53+
ConfigurableAxis CfgMultBins{"CfgMultBins", {VARIABLE_WIDTH, 0.0, 20.0, 40.0, 60.0, 80.0, 500.0}, "Mixing bins - number of contributor"};
54+
55+
// THnsparse bining
56+
ConfigurableAxis configThnAxisInvMass{"configThnAxisInvMass", {1500, 2.0, 3.5}, "#it{M} (GeV/#it{c}^{2})"};
57+
ConfigurableAxis configThnAxisDaugherPt{"configThnAxisDaugherPt", {25, 0.0, 50.}, "#it{p}_{T} (GeV/#it{c})"};
58+
ConfigurableAxis configThnAxisPt{"configThnAxisPt", {40, 0.0, 20.}, "#it{p}_{T} (GeV/#it{c})"};
59+
ConfigurableAxis configThnAxisKstar{"configThnAxisKstar", {50, 0.0, 0.5}, "#it{k}^{*} (GeV/#it{c})"};
60+
61+
// Initialize the ananlysis task
62+
void init(o2::framework::InitContext&)
63+
{
64+
// register histograms
65+
histos.add("hnsigmaTPCKaonPlus", "hnsigmaTPCKaonPlus", kTH2F, {{1000, -3.0, 3.0f}, {100, 0.0f, 10.0f}});
66+
histos.add("hnsigmaTPCKaonMinus", "hnsigmaTPCKaonMinus", kTH2F, {{1000, -3.0, 3.0f}, {100, 0.0f, 10.0f}});
67+
histos.add("hPhid1Mass", "hPhid1Mass", kTH2F, {{40, 1.0, 1.04f}, {100, 0.0f, 10.0f}});
68+
histos.add("hPhid2Mass", "hPhid2Mass", kTH2F, {{40, 1.0, 1.04f}, {100, 0.0f, 10.0f}});
69+
70+
const AxisSpec thnAxisInvMass{configThnAxisInvMass, "#it{M} (GeV/#it{c}^{2})"};
71+
const AxisSpec thnAxisDaughterPt{configThnAxisDaugherPt, "#it{p}_{T} (GeV/#it{c})"};
72+
const AxisSpec thnAxisPt{configThnAxisPt, "#it{p}_{T} (GeV/#it{c})"};
73+
const AxisSpec thnAxisKstar{configThnAxisKstar, "#it{k}^{*} (GeV/#it{c})"};
74+
75+
histos.add("SEMassUnlike", "SEMassUnlike", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisDaughterPt, thnAxisDaughterPt, thnAxisKstar});
76+
histos.add("SEMassRot", "SEMassRot", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisDaughterPt, thnAxisDaughterPt, thnAxisKstar});
77+
78+
histos.add("MEMassUnlike", "MEMassUnlike", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisDaughterPt, thnAxisDaughterPt, thnAxisKstar});
79+
histos.add("MEMassRot", "MEMassRot", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisDaughterPt, thnAxisDaughterPt, thnAxisKstar});
80+
}
81+
82+
// get kstar
83+
TLorentzVector trackSum, PartOneCMS, PartTwoCMS, trackRelK;
84+
float getkstar(const TLorentzVector part1,
85+
const TLorentzVector part2)
86+
{
87+
// const TLorentzVector trackSum = part1 + part2;
88+
trackSum = part1 + part2;
89+
const float beta = trackSum.Beta();
90+
const float betax = beta * std::cos(trackSum.Phi()) * std::sin(trackSum.Theta());
91+
const float betay = beta * std::sin(trackSum.Phi()) * std::sin(trackSum.Theta());
92+
const float betaz = beta * std::cos(trackSum.Theta());
93+
// TLorentzVector PartOneCMS(part1);
94+
// TLorentzVector PartTwoCMS(part2);
95+
PartOneCMS.SetXYZM(part1.Px(), part1.Py(), part1.Pz(), part1.M());
96+
PartTwoCMS.SetXYZM(part2.Px(), part2.Py(), part2.Pz(), part2.M());
97+
const ROOT::Math::Boost boostPRF = ROOT::Math::Boost(-betax, -betay, -betaz);
98+
PartOneCMS = boostPRF(PartOneCMS);
99+
PartTwoCMS = boostPRF(PartTwoCMS);
100+
// const TLorentzVector trackRelK = PartOneCMS - PartTwoCMS;
101+
trackRelK = PartOneCMS - PartTwoCMS;
102+
return 0.5 * trackRelK.P();
103+
}
104+
bool selectionPID(float nsigmaTPC, float nsigmaTOF, int TOFHit, int PIDStrategy, float ptcand)
105+
{
106+
if (PIDStrategy == 0) {
107+
if (TOFHit != 1) {
108+
if (TMath::Abs(nsigmaTPC) < cutNsigmaTPC) {
109+
return true;
110+
}
111+
}
112+
if (TOFHit == 1) {
113+
if (TMath::Abs(nsigmaTOF) < cutNsigmaTOF) {
114+
return true;
115+
}
116+
}
117+
}
118+
if (PIDStrategy == 1) {
119+
if (ptcand < 0.6) {
120+
if (TOFHit != 1 && TMath::Abs(nsigmaTPC) < cutNsigmaTPC) {
121+
return true;
122+
}
123+
if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cutNsigmaTOF) {
124+
return true;
125+
}
126+
}
127+
if (ptcand >= 0.6) {
128+
if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cutNsigmaTOF) {
129+
return true;
130+
}
131+
}
132+
}
133+
if (PIDStrategy == 2) {
134+
if (ptcand < 0.6) {
135+
if (TOFHit != 1 && TMath::Abs(nsigmaTPC) < cutNsigmaTPC) {
136+
return true;
137+
}
138+
if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cutNsigmaTOF) {
139+
return true;
140+
}
141+
}
142+
if (ptcand >= 0.6 && ptcand < 1.2) {
143+
if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cutNsigmaTOF) {
144+
return true;
145+
}
146+
if (TOFHit != 1 && nsigmaTPC > -1.5 && nsigmaTPC < cutNsigmaTPC) {
147+
return true;
148+
}
149+
}
150+
if (ptcand >= 1.2) {
151+
if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cutNsigmaTOF) {
152+
return true;
153+
}
154+
if (TOFHit != 1 && TMath::Abs(nsigmaTPC) < cutNsigmaTPC) {
155+
return true;
156+
}
157+
}
158+
}
159+
if (PIDStrategy == 3) {
160+
if (ptcand < 0.6) {
161+
if (TOFHit != 1 && TMath::Abs(nsigmaTPC) < cutNsigmaTPC) {
162+
return true;
163+
}
164+
if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cutNsigmaTOF) {
165+
return true;
166+
}
167+
}
168+
if (ptcand >= 0.6 && ptcand < 1.2) {
169+
if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cutNsigmaTOF) {
170+
return true;
171+
}
172+
}
173+
if (ptcand >= 1.2) {
174+
if (TOFHit == 1 && TMath::Abs(nsigmaTOF) < cutNsigmaTOF) {
175+
return true;
176+
}
177+
if (TOFHit != 1 && TMath::Abs(nsigmaTPC) < cutNsigmaTPC) {
178+
return true;
179+
}
180+
}
181+
}
182+
return false;
183+
}
184+
185+
TLorentzVector exotic, Phid1, Phid2;
186+
TLorentzVector exoticRot, Phid1Rot;
187+
void process(aod::RedPhiEvents::iterator const& collision, aod::PhiTracks const& phitracks)
188+
{
189+
if (additionalEvsel && (collision.numPos() < 2 || collision.numNeg() < 2)) {
190+
return;
191+
}
192+
for (auto phitrackd1 : phitracks) {
193+
if (phitrackd1.phiMass() < minPhiMass || phitrackd1.phiMass() > maxPhiMass) {
194+
continue;
195+
}
196+
197+
auto kaonplusd1pt = TMath::Sqrt(phitrackd1.phid1Px() * phitrackd1.phid1Px() + phitrackd1.phid1Py() * phitrackd1.phid1Py());
198+
auto kaonminusd1pt = TMath::Sqrt(phitrackd1.phid2Px() * phitrackd1.phid2Px() + phitrackd1.phid2Py() * phitrackd1.phid2Py());
199+
200+
if (!selectionPID(phitrackd1.phid1TPC(), phitrackd1.phid1TOF(), phitrackd1.phid1TOFHit(), strategyPID, kaonplusd1pt)) {
201+
continue;
202+
}
203+
if (!selectionPID(phitrackd1.phid2TPC(), phitrackd1.phid2TOF(), phitrackd1.phid2TOFHit(), strategyPID, kaonminusd1pt)) {
204+
continue;
205+
}
206+
207+
histos.fill(HIST("hnsigmaTPCKaonPlus"), phitrackd1.phid1TPC(), kaonplusd1pt);
208+
histos.fill(HIST("hnsigmaTPCKaonMinus"), phitrackd1.phid2TPC(), kaonminusd1pt);
209+
Phid1.SetXYZM(phitrackd1.phiPx(), phitrackd1.phiPy(), phitrackd1.phiPz(), phitrackd1.phiMass());
210+
histos.fill(HIST("hPhid1Mass"), Phid1.M(), Phid1.Pt());
211+
auto phid1id = phitrackd1.index();
212+
for (auto phitrackd2 : phitracks) {
213+
auto phid2id = phitrackd2.index();
214+
if (phid2id <= phid1id) {
215+
continue;
216+
}
217+
if (phitrackd2.phiMass() < minPhiMass || phitrackd2.phiMass() > maxPhiMass) {
218+
continue;
219+
}
220+
auto kaonplusd2pt = TMath::Sqrt(phitrackd2.phid1Px() * phitrackd2.phid1Px() + phitrackd2.phid1Py() * phitrackd2.phid1Py());
221+
auto kaonminusd2pt = TMath::Sqrt(phitrackd2.phid2Px() * phitrackd2.phid2Px() + phitrackd2.phid2Py() * phitrackd2.phid2Py());
222+
223+
if (!selectionPID(phitrackd2.phid1TPC(), phitrackd2.phid1TOF(), phitrackd2.phid1TOFHit(), strategyPID, kaonplusd2pt)) {
224+
continue;
225+
}
226+
if (!selectionPID(phitrackd2.phid2TPC(), phitrackd2.phid2TOF(), phitrackd2.phid2TOFHit(), strategyPID, kaonminusd2pt)) {
227+
continue;
228+
}
229+
if (phitrackd1.phid1Index() == phitrackd2.phid1Index()) {
230+
continue;
231+
}
232+
if (phitrackd1.phid2Index() == phitrackd2.phid2Index()) {
233+
continue;
234+
}
235+
Phid2.SetXYZM(phitrackd2.phiPx(), phitrackd2.phiPy(), phitrackd2.phiPz(), phitrackd2.phiMass());
236+
exotic = Phid1 + Phid2;
237+
auto kstar = getkstar(Phid1, Phid2);
238+
histos.fill(HIST("SEMassUnlike"), exotic.M(), exotic.Pt(), Phid1.Pt(), Phid2.Pt(), kstar);
239+
if (fillRotation) {
240+
for (int nrotbkg = 0; nrotbkg < 9; nrotbkg++) {
241+
auto anglestart = 5.0 * TMath::Pi() / 6.0;
242+
auto angleend = 7.0 * TMath::Pi() / 6.0;
243+
auto anglestep = (angleend - anglestart) / (1.0 * (9.0 - 1.0));
244+
auto rotangle = anglestart + nrotbkg * anglestep;
245+
auto rotd1px = Phid1.Px() * std::cos(rotangle) - Phid1.Py() * std::sin(rotangle);
246+
auto rotd1py = Phid1.Px() * std::sin(rotangle) + Phid1.Py() * std::cos(rotangle);
247+
Phid1Rot.SetXYZM(rotd1px, rotd1py, Phid1.Pz(), Phid1.M());
248+
exoticRot = Phid1Rot + Phid2;
249+
auto kstar_rot = getkstar(Phid1Rot, Phid2);
250+
histos.fill(HIST("SEMassRot"), exoticRot.M(), exoticRot.Pt(), Phid1Rot.Pt(), Phid2.Pt(), kstar_rot);
251+
}
252+
}
253+
}
254+
}
255+
}
256+
257+
SliceCache cache;
258+
using BinningTypeVertexContributor = ColumnBinningPolicy<aod::collision::PosZ, aod::collision::NumContrib>;
259+
void processMixedEvent(aod::RedPhiEvents& collisions, aod::PhiTracks& phitracks)
260+
{
261+
auto tracksTuple = std::make_tuple(phitracks);
262+
BinningTypeVertexContributor binningOnPositions{{CfgVtxBins, CfgMultBins}, true};
263+
SameKindPair<aod::RedPhiEvents, aod::PhiTracks, BinningTypeVertexContributor> pair{binningOnPositions, nEvtMixing, -1, collisions, tracksTuple, &cache};
264+
265+
for (auto& [collision1, tracks1, collision2, tracks2] : pair) {
266+
if (collision1.index() == collision2.index()) {
267+
continue;
268+
}
269+
for (auto& [phitrackd1, phitrackd2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) {
270+
if (phitrackd1.phiMass() < minPhiMass || phitrackd1.phiMass() > maxPhiMass) {
271+
continue;
272+
}
273+
274+
auto kaonplusd1pt = TMath::Sqrt(phitrackd1.phid1Px() * phitrackd1.phid1Px() + phitrackd1.phid1Py() * phitrackd1.phid1Py());
275+
auto kaonminusd1pt = TMath::Sqrt(phitrackd1.phid2Px() * phitrackd1.phid2Px() + phitrackd1.phid2Py() * phitrackd1.phid2Py());
276+
auto kaonplusd2pt = TMath::Sqrt(phitrackd2.phid1Px() * phitrackd2.phid1Px() + phitrackd2.phid1Py() * phitrackd2.phid1Py());
277+
auto kaonminusd2pt = TMath::Sqrt(phitrackd2.phid2Px() * phitrackd2.phid2Px() + phitrackd2.phid2Py() * phitrackd2.phid2Py());
278+
279+
if (!selectionPID(phitrackd1.phid1TPC(), phitrackd1.phid1TOF(), phitrackd1.phid1TOFHit(), strategyPID, kaonplusd1pt)) {
280+
continue;
281+
}
282+
if (!selectionPID(phitrackd1.phid2TPC(), phitrackd1.phid2TOF(), phitrackd1.phid2TOFHit(), strategyPID, kaonminusd1pt)) {
283+
continue;
284+
}
285+
Phid1.SetXYZM(phitrackd1.phiPx(), phitrackd1.phiPy(), phitrackd1.phiPz(), phitrackd1.phiMass());
286+
if (phitrackd2.phiMass() < minPhiMass || phitrackd2.phiMass() > maxPhiMass) {
287+
continue;
288+
}
289+
if (!selectionPID(phitrackd2.phid1TPC(), phitrackd2.phid1TOF(), phitrackd2.phid1TOFHit(), strategyPID, kaonplusd2pt)) {
290+
continue;
291+
}
292+
if (!selectionPID(phitrackd2.phid2TPC(), phitrackd2.phid2TOF(), phitrackd2.phid2TOFHit(), strategyPID, kaonminusd2pt)) {
293+
continue;
294+
}
295+
Phid2.SetXYZM(phitrackd2.phiPx(), phitrackd2.phiPy(), phitrackd2.phiPz(), phitrackd2.phiMass());
296+
exotic = Phid1 + Phid2;
297+
auto kstar = getkstar(Phid1, Phid2);
298+
histos.fill(HIST("MEMassUnlike"), exotic.M(), exotic.Pt(), Phid1.Pt(), Phid2.Pt(), kstar);
299+
if (fillRotation) {
300+
for (int nrotbkg = 0; nrotbkg < 9; nrotbkg++) {
301+
auto anglestart = 5.0 * TMath::Pi() / 6.0;
302+
auto angleend = 7.0 * TMath::Pi() / 6.0;
303+
auto anglestep = (angleend - anglestart) / (1.0 * (9.0 - 1.0));
304+
auto rotangle = anglestart + nrotbkg * anglestep;
305+
auto rotd1px = Phid1.Px() * std::cos(rotangle) - Phid1.Py() * std::sin(rotangle);
306+
auto rotd1py = Phid1.Px() * std::sin(rotangle) + Phid1.Py() * std::cos(rotangle);
307+
Phid1Rot.SetXYZM(rotd1px, rotd1py, Phid1.Pz(), Phid1.M());
308+
exoticRot = Phid1Rot + Phid2;
309+
auto kstar_rot = getkstar(Phid1Rot, Phid2);
310+
histos.fill(HIST("MEMassRot"), exoticRot.M(), exoticRot.Pt(), Phid1Rot.Pt(), Phid2.Pt(), kstar_rot);
311+
}
312+
}
313+
}
314+
}
315+
}
316+
PROCESS_SWITCH(doublephimeson, processMixedEvent, "Process EventMixing for combinatorial background", true);
317+
};
318+
319+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask<doublephimeson>(cfgc)}; }

0 commit comments

Comments
 (0)