Skip to content

Commit 3e8d361

Browse files
siragonisiragoni
authored
[PWGUD] Exclusive phi->ee (DG) tree producer (#8401)
1 parent 25b99a9 commit 3e8d361

2 files changed

Lines changed: 315 additions & 0 deletions

File tree

‎PWGUD/Tasks/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -163,6 +163,11 @@ o2physics_add_dpl_workflow(exclusive-phi-leptons
163163
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::DGPIDSelector
164164
COMPONENT_NAME Analysis)
165165

166+
o2physics_add_dpl_workflow(exclusive-phi-leptons-trees
167+
SOURCES exclusivePhiLeptonsTrees.cxx
168+
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::DGPIDSelector
169+
COMPONENT_NAME Analysis)
170+
166171
o2physics_add_dpl_workflow(exclusive-pentaquark
167172
SOURCES exclusivePentaquark.cxx
168173
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::DGPIDSelector

‎PWGUD/Tasks/exclusivePhiLeptonsTrees.cxx‎

Lines changed: 310 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,310 @@
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+
#include <vector>
12+
#include "Framework/runDataProcessing.h"
13+
#include "Framework/AnalysisTask.h"
14+
#include "Framework/AnalysisDataModel.h"
15+
#include "iostream"
16+
#include "PWGUD/DataModel/UDTables.h"
17+
#include <TString.h>
18+
#include "TLorentzVector.h"
19+
#include "Common/DataModel/PIDResponse.h"
20+
#include "PWGUD/Core/SGSelector.h"
21+
using std::array;
22+
using namespace std;
23+
using namespace o2;
24+
using namespace o2::aod;
25+
using namespace o2::framework;
26+
using namespace o2::framework::expressions;
27+
28+
/// \brief Exclusive phi->ee tree producer, for ML applications, DG-based
29+
/// \author Simone Ragoni, Creighton
30+
/// \date 11/11/2024
31+
32+
namespace o2::aod
33+
{
34+
namespace tree
35+
{
36+
// track tables
37+
DECLARE_SOA_COLUMN(PX1, px1, float);
38+
DECLARE_SOA_COLUMN(PY1, py1, float);
39+
DECLARE_SOA_COLUMN(PZ1, pz1, float);
40+
DECLARE_SOA_COLUMN(PE1, pE1, float);
41+
DECLARE_SOA_COLUMN(PX2, px2, float);
42+
DECLARE_SOA_COLUMN(PY2, py2, float);
43+
DECLARE_SOA_COLUMN(PZ2, pz2, float);
44+
DECLARE_SOA_COLUMN(PE2, pE2, float);
45+
DECLARE_SOA_COLUMN(NCOUNTERPV, nCounterPV, int);
46+
DECLARE_SOA_COLUMN(NELECTRONSTOF, nElectronsTOF, int);
47+
} // namespace tree
48+
49+
DECLARE_SOA_TABLE(TREE, "AOD", "Tree", tree::PX1, tree::PY1, tree::PZ1, tree::PE1, tree::PX2, tree::PY2, tree::PZ2, tree::PE2, tree::NCOUNTERPV, tree::NELECTRONSTOF);
50+
} // namespace o2::aod
51+
52+
struct ExclusivePhiLeptonsTrees {
53+
Produces<aod::TREE> tree;
54+
Configurable<float> gap_Side{"gap", 2, "gap selection"};
55+
Configurable<float> pid2d_cut{"PID2D", 2., "PID cut in 2D"};
56+
Configurable<float> pid_cut{"PID", 2., "PID cut in 1D"};
57+
Configurable<std::size_t> electronsInTOF{"eTOF", 2, "electrons in TOF"};
58+
// defining histograms using histogram registry
59+
HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject};
60+
61+
//-----------------------------------------------------------------------------------------------------------------------
62+
void init(o2::framework::InitContext&)
63+
{
64+
registry.add("posx", "Vertex position in x", kTH1F, {{100, -0.5, 0.5}});
65+
registry.add("posy", "Vertex position in y", kTH1F, {{100, -0.5, 0.5}});
66+
registry.add("posz", "Vertex position in z", kTH1F, {{1000, -100., 100.}});
67+
registry.add("hITSCluster", "N_{cluster}", kTH1F, {{100, -0.5, 99.5}});
68+
registry.add("hChi2ITSTrkSegment", "N_{cluster}", kTH1F, {{100, -0.5, 99.5}});
69+
registry.add("hTPCCluster", "N_{cluster}", kTH1F, {{200, -0.5, 199.5}});
70+
registry.add("hdEdx", "p vs dE/dx Signal", kTH2F, {{100, 0.0, 3.0}, {1000, 0.0, 2000.0}});
71+
registry.add("hdEdx2", "p vs dE/dx Signal", kTH2F, {{100, 0.0, 3.0}, {1000, 0.0, 2000.0}});
72+
registry.add("hdSigmaElectron", "p vs dE/dx sigma electron", kTH2F, {{100, 0.0, 3.0}, {1000, -500.0, 500.0}});
73+
registry.add("hdSigmaElectron2", "p vs dE/dx sigma electron", kTH2F, {{100, 0.0, 3.0}, {1000, -500.0, 500.0}});
74+
registry.add("hdSigmaElectron3", "p vs dE/dx sigma electron", kTH2F, {{100, 0.0, 3.0}, {1000, -500.0, 500.0}});
75+
registry.add("hNsigEvsKa1", "NSigma(t1) vs NSigma (t2);n#sigma_{1};n#sigma_{2}", kTH2F, {{100, -15., 15.}, {100, -15., 15}});
76+
registry.add("hNsigEvsKa2", "NSigma(t1) vs NSigma (t2);n#sigma_{1};n#sigma_{2}", kTH2F, {{100, -15., 15.}, {100, -15., 15}});
77+
registry.add("hMomentum", "p_{#ka};#it{p_{trk}}, GeV/c;", kTH1F, {{100, 0., 3.}});
78+
registry.add("hEta1", "#eta_{#ka};#it{#eta_{trk}}, GeV/c;", kTH1F, {{100, -2., 2.}});
79+
registry.add("hPtLikeSignElectron", "Pt;#it{p_{t}}, GeV/c;", kTH1F, {{500, 0., 5.}});
80+
registry.add("hMassLikeSignElectron", "Raw Inv.M;#it{m_{ee}}, GeV/c^{2};", kTH1F, {{1000, 0., 10.}});
81+
registry.add("hMassPtLikeSignElectron", "Raw Inv.M;#it{m_{ee}}, GeV/c^{2};Pt;#it{p_{t}}, GeV/c;", kTH2F, {{1000, 0., 10.}, {400, 0., 4.}});
82+
83+
auto hSelectionCounter = registry.add<TH1>("hSelectionCounter", "hSelectionCounter;;NEvents", HistType::kTH1I, {{10, 0., 10.}});
84+
85+
TString SelectionCuts[9] = {"NoSelection", "GAPcondition", "PVtracks", "Good TPC-ITS track", "TPC/TOF PID track", "End trk loop", "Exactly 2e", "Like-sign ev", "Unlike-sign ev"};
86+
// now we can set BinLabel in histogram Registry
87+
88+
for (int i = 0; i < 9; i++) {
89+
hSelectionCounter->GetXaxis()->SetBinLabel(i + 1, SelectionCuts[i].Data());
90+
}
91+
92+
// Unlike sign pp
93+
registry.add("ee/hRapidity", "Rapidity;#it{y_{ee}};", kTH1F, {{100, -2., 2.}});
94+
registry.add("ee/hPtElectronVsElectron", "Pt1 vs Pt2;p_{T};p_{T};", kTH2F, {{100, 0., 3.}, {100, 0., 3.}});
95+
registry.add("ee/hMassPtUnlikeSignElectron", "Raw Inv.M;#it{m_{ee}}, GeV/c^{2};Pt;#it{p_{t}}, GeV/c;", kTH2F, {{400, 0., 4.}, {400, 0., 4.}});
96+
registry.add("ee/hMassUnlike", "m_{ee} [GeV/#it{c}^{2}]", kTH1F, {{1000, 0., 10.}});
97+
registry.add("ee/hUnlikePt", "Pt;#it{p_{t}}, GeV/c;", kTH1F, {{500, 0., 5.}});
98+
registry.add("ee/hCoherentMass", "Raw Inv.M;#it{m_{ee}}, GeV/c^{2};", kTH1F, {{1000, 0., 10.}});
99+
registry.add("ee/hIncoherentMass", "Raw Inv.M;#it{m_{ee}}, GeV/c^{2};", kTH1F, {{1000, 0., 10.}});
100+
}
101+
102+
using udtracks = soa::Join<aod::UDTracks, aod::UDTracksExtra, aod::UDTracksPID>;
103+
using udtracksfull = soa::Join<aod::UDTracks, aod::UDTracksPID, aod::UDTracksExtra, aod::UDTracksFlags, aod::UDTracksDCA>;
104+
// using UDCollisions = soa::Join<aod::UDCollisions, aod::UDCollisionsSels, aod::UDZdcsReduced>;
105+
//__________________________________________________________________________
106+
// Main process
107+
void process(UDCollisions::iterator const& collision, udtracksfull const& tracks)
108+
{
109+
registry.fill(HIST("hSelectionCounter"), 0);
110+
registry.fill(HIST("posx"), collision.posX());
111+
registry.fill(HIST("posy"), collision.posY());
112+
registry.fill(HIST("posz"), collision.posZ());
113+
TLorentzVector resonance; // lorentz vectors of tracks and the mother
114+
// ===================================
115+
// Task for ee pairs with PID
116+
// Topology:
117+
// - 2 TOF ee
118+
// - 1 TOF e + 1 TPC e
119+
// ===================================
120+
std::vector<TLorentzVector> onlyElectronTracks;
121+
std::vector<TLorentzVector> onlyElectronTracksTOF;
122+
std::vector<float> onlyElectronSigma;
123+
std::vector<float> onlyElectronSigmaTOF;
124+
std::vector<decltype(tracks.begin())> rawElectronTracks;
125+
std::vector<decltype(tracks.begin())> rawElectronTracksTOF;
126+
127+
// -------------------------------------------
128+
// TO BE SAVED:
129+
// - counterPV
130+
// - electronsTOF (0,1,2) = 2 - electronsTPC
131+
// - (px,py,pz,E)1
132+
// - (px,py,pz,E)2
133+
int counterPV = 0;
134+
for (auto trk : tracks) {
135+
// ----------------------------------------
136+
// SELECTIONS:
137+
// - PV track
138+
// - at least 70 TPC clusters
139+
// - at least 6 ITS clusters
140+
// - 0.3 < pT < 0.65 GeV from STARlight
141+
// - DCAxy, DCAz
142+
// - Nsigma^2 < 2^2
143+
// - |track eta| < 0.8
144+
if (!trk.isPVContributor()) {
145+
continue;
146+
}
147+
counterPV += 1;
148+
registry.fill(HIST("hSelectionCounter"), 2);
149+
150+
int NFindable = trk.tpcNClsFindable();
151+
int NMinusFound = trk.tpcNClsFindableMinusFound();
152+
int NCluster = NFindable - NMinusFound;
153+
registry.fill(HIST("hTPCCluster"), NCluster);
154+
registry.fill(HIST("hChi2ITSTrkSegment"), trk.itsChi2NCl());
155+
if (NCluster < 70) {
156+
continue;
157+
}
158+
if (trk.itsNCls() < 6) {
159+
continue;
160+
}
161+
if (trk.pt() < 0.300) {
162+
continue;
163+
}
164+
if (trk.pt() > 0.650) {
165+
continue;
166+
}
167+
if (!(std::abs(trk.dcaZ()) < 2.)) {
168+
continue;
169+
}
170+
double dcaLimit = 0.0105 + 0.035 / pow(trk.pt(), 1.1);
171+
if (!(std::abs(trk.dcaXY()) < dcaLimit)) {
172+
continue;
173+
}
174+
registry.fill(HIST("hSelectionCounter"), 3);
175+
registry.fill(HIST("hITSCluster"), trk.itsNCls());
176+
177+
double momentum = TMath::Sqrt(trk.px() * trk.px() + trk.py() * trk.py() + trk.pz() * trk.pz());
178+
double dEdx = trk.tpcSignal();
179+
registry.fill(HIST("hdEdx"), momentum, dEdx);
180+
181+
TLorentzVector electron;
182+
electron.SetXYZM(trk.px(), trk.py(), trk.pz(), o2::constants::physics::MassElectron);
183+
if (fabs(electron.Eta()) > 0.8) {
184+
return;
185+
}
186+
auto nSigmaEl = trk.tpcNSigmaEl();
187+
auto nSigmaElTOF = trk.tofNSigmaEl();
188+
189+
if (trk.hasTOF()) {
190+
registry.fill(HIST("hdSigmaElectron"), momentum, nSigmaElTOF);
191+
}
192+
if (fabs(nSigmaEl) < pid_cut) {
193+
registry.fill(HIST("hdEdx2"), momentum, dEdx);
194+
registry.fill(HIST("hdSigmaElectron2"), momentum, nSigmaEl);
195+
registry.fill(HIST("hMomentum"), momentum);
196+
if (trk.hasTOF() && fabs(nSigmaElTOF) < pid_cut) {
197+
registry.fill(HIST("hdSigmaElectron3"), momentum, nSigmaElTOF);
198+
onlyElectronTracksTOF.push_back(electron);
199+
onlyElectronSigmaTOF.push_back(nSigmaElTOF);
200+
rawElectronTracksTOF.push_back(trk);
201+
} else if (!trk.hasTOF()) {
202+
onlyElectronTracks.push_back(electron);
203+
onlyElectronSigma.push_back(nSigmaEl);
204+
rawElectronTracks.push_back(trk);
205+
}
206+
registry.fill(HIST("hSelectionCounter"), 4);
207+
}
208+
209+
} // trk loop
210+
211+
registry.fill(HIST("hSelectionCounter"), 5);
212+
if ((onlyElectronTracksTOF.size() >= electronsInTOF) && (onlyElectronTracks.size() + onlyElectronTracksTOF.size()) == 2) {
213+
registry.fill(HIST("hSelectionCounter"), 6);
214+
215+
int signSum = -999.;
216+
double sigmaTotal = -999.;
217+
TLorentzVector a, b;
218+
// two electrons in the TPC
219+
if (onlyElectronTracksTOF.size() == 0) {
220+
221+
registry.fill(HIST("hEta1"), onlyElectronTracks[0].Eta());
222+
registry.fill(HIST("hEta1"), onlyElectronTracks[1].Eta());
223+
resonance += onlyElectronTracks[0];
224+
resonance += onlyElectronTracks[1];
225+
a += onlyElectronTracks[0];
226+
b += onlyElectronTracks[1];
227+
sigmaTotal = 0;
228+
sigmaTotal = onlyElectronSigma[0] * onlyElectronSigma[0] + onlyElectronSigma[1] * onlyElectronSigma[1];
229+
;
230+
registry.fill(HIST("hNsigEvsKa1"), onlyElectronSigma[0], onlyElectronSigma[1]);
231+
signSum = rawElectronTracks[0].sign() + rawElectronTracks[1].sign();
232+
if (signSum == 0) {
233+
registry.fill(HIST("ee/hPtElectronVsElectron"), onlyElectronTracks[0].Pt(), onlyElectronTracks[1].Pt());
234+
}
235+
236+
} else if (onlyElectronTracksTOF.size() == 1) {
237+
238+
registry.fill(HIST("hEta1"), onlyElectronTracks[0].Eta());
239+
registry.fill(HIST("hEta1"), onlyElectronTracksTOF[0].Eta());
240+
resonance += onlyElectronTracks[0];
241+
resonance += onlyElectronTracksTOF[0];
242+
a += onlyElectronTracks[0];
243+
b += onlyElectronTracksTOF[0];
244+
sigmaTotal = 0;
245+
sigmaTotal = onlyElectronSigma[0] * onlyElectronSigma[0] + onlyElectronSigmaTOF[0] * onlyElectronSigmaTOF[0];
246+
;
247+
registry.fill(HIST("hNsigEvsKa1"), onlyElectronSigma[0], onlyElectronSigmaTOF[0]);
248+
signSum = rawElectronTracks[0].sign() + rawElectronTracksTOF[0].sign();
249+
if (signSum == 0) {
250+
registry.fill(HIST("ee/hPtElectronVsElectron"), onlyElectronTracks[0].Pt(), onlyElectronTracksTOF[0].Pt());
251+
}
252+
253+
} else if (onlyElectronTracksTOF.size() == 2) {
254+
255+
registry.fill(HIST("hEta1"), onlyElectronTracksTOF[0].Eta());
256+
registry.fill(HIST("hEta1"), onlyElectronTracksTOF[1].Eta());
257+
resonance += onlyElectronTracksTOF[0];
258+
resonance += onlyElectronTracksTOF[1];
259+
a += onlyElectronTracksTOF[0];
260+
b += onlyElectronTracksTOF[1];
261+
sigmaTotal = 0;
262+
sigmaTotal = onlyElectronSigmaTOF[0] * onlyElectronSigmaTOF[0] + onlyElectronSigmaTOF[1] * onlyElectronSigmaTOF[1];
263+
;
264+
registry.fill(HIST("hNsigEvsKa1"), onlyElectronSigmaTOF[0], onlyElectronSigmaTOF[1]);
265+
signSum = rawElectronTracksTOF[0].sign() + rawElectronTracksTOF[1].sign();
266+
if (signSum == 0) {
267+
registry.fill(HIST("ee/hPtElectronVsElectron"), onlyElectronTracksTOF[0].Pt(), onlyElectronTracksTOF[1].Pt());
268+
}
269+
}
270+
271+
if (sigmaTotal > pid2d_cut * pid2d_cut) {
272+
return;
273+
}
274+
if (onlyElectronTracksTOF.size() == 1) {
275+
registry.fill(HIST("hNsigEvsKa2"), onlyElectronSigma[0], onlyElectronSigmaTOF[0]);
276+
} else if (onlyElectronTracksTOF.size() == 2) {
277+
registry.fill(HIST("hNsigEvsKa2"), onlyElectronSigmaTOF[0], onlyElectronSigmaTOF[1]);
278+
}
279+
280+
if (signSum != 0) {
281+
registry.fill(HIST("hMassPtLikeSignElectron"), resonance.M(), resonance.Pt());
282+
registry.fill(HIST("hSelectionCounter"), 7);
283+
registry.fill(HIST("hPtLikeSignElectron"), resonance.Pt());
284+
registry.fill(HIST("hMassLikeSignElectron"), resonance.M());
285+
} else {
286+
registry.fill(HIST("ee/hMassPtUnlikeSignElectron"), resonance.M(), resonance.Pt());
287+
registry.fill(HIST("hSelectionCounter"), 8);
288+
registry.fill(HIST("ee/hMassUnlike"), resonance.M());
289+
registry.fill(HIST("ee/hRapidity"), resonance.Rapidity());
290+
if (resonance.Pt() > 0.1) {
291+
registry.fill(HIST("ee/hIncoherentMass"), resonance.M());
292+
} else {
293+
registry.fill(HIST("ee/hCoherentMass"), resonance.M());
294+
}
295+
if (resonance.M() > 1.01 && resonance.M() < 1.03) {
296+
registry.fill(HIST("ee/hUnlikePt"), resonance.Pt());
297+
}
298+
}
299+
// Filling tree, make to be consistent with the declared tables
300+
tree(a.Px(), a.Py(), a.Pz(), a.E(), b.Px(), b.Py(), b.Pz(), b.E(), counterPV, onlyElectronTracksTOF.size());
301+
}
302+
} // end of process
303+
304+
}; // end of struct
305+
306+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
307+
{
308+
return WorkflowSpec{
309+
adaptAnalysisTask<ExclusivePhiLeptonsTrees>(cfgc)};
310+
}

0 commit comments

Comments
 (0)