Skip to content

Commit eead5eb

Browse files
alcalivaalcaliva
authored
[PWGLF] calculation of overlap between jet-jet, ue-ue, jet-ue (#9669)
1 parent b4cf386 commit eead5eb

File tree

1 file changed

+198
-13
lines changed

1 file changed

+198
-13
lines changed

‎PWGLF/Tasks/Nuspex/nucleiInJets.cxx‎

Lines changed: 198 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -97,6 +97,7 @@ struct NucleiInJets {
9797
Configurable<int> minNparticlesInJet{"minNparticlesInJet", 2, "Minimum number of particles inside jet"};
9898
Configurable<int> nJetsPerEventMax{"nJetsPerEventMax", 1000, "Maximum number of jets per event"};
9999
Configurable<bool> requireNoOverlap{"requireNoOverlap", true, "require no overlap between jets and UE cones"};
100+
Configurable<int> nGhosts{"nGhosts", 1000, "number of ghost particles"};
100101

101102
// Track Parameters
102103
Configurable<double> par0{"par0", 0.00164, "par 0"};
@@ -165,6 +166,7 @@ struct NucleiInJets {
165166
registryQC.add("dcaxy_vs_pt", "dcaxy_vs_pt", HistType::kTH2F, {{100, 0.0, 5.0, "#it{p}_{T} (GeV/#it{c})"}, {2000, -0.05, 0.05, "DCA_{xy} (cm)"}});
166167
registryQC.add("dcaz_vs_pt", "dcaz_vs_pt", HistType::kTH2F, {{100, 0.0, 5.0, "#it{p}_{T} (GeV/#it{c})"}, {2000, -0.05, 0.05, "DCA_{z} (cm)"}});
167168
registryQC.add("jet_ue_overlaps", "jet_ue_overlaps", HistType::kTH2F, {{20, 0.0, 20.0, "#it{n}_{jet}"}, {200, 0.0, 200.0, "#it{n}_{overlaps}"}});
169+
registryQC.add("hJetArea", "hJetArea", HistType::kTH1F, {{450, 0, 15, "Area"}});
168170

169171
// Event Counters
170172
registryData.add("number_of_events_data", "number of events in data", HistType::kTH1F, {{10, 0, 10, "counter"}});
@@ -239,6 +241,17 @@ struct NucleiInJets {
239241
registryMC.add("antiproton_eta_pt_pythia", "antiproton_eta_pt_pythia", HistType::kTH2F, {{200, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {20, -1.0, 1.0, "#it{#eta}"}});
240242
registryMC.add("antiproton_eta_pt_jet", "antiproton_eta_pt_jet", HistType::kTH2F, {{200, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {20, -1.0, 1.0, "#it{#eta}"}});
241243
registryMC.add("antiproton_eta_pt_ue", "antiproton_eta_pt_ue", HistType::kTH2F, {{200, 0.0, 10.0, "#it{p}_{T} (GeV/#it{c})"}, {20, -1.0, 1.0, "#it{#eta}"}});
244+
245+
// Detector Response Matrix
246+
registryMC.add("detectorResponseMatrix", "detectorResponseMatrix", HistType::kTH2F, {{500, 0.0, 50.0, "#it{p}_{T}^{gen} (GeV/#it{c})"}, {500, 0.0, 50.0, "#it{p}_{T}^{rec} (GeV/#it{c})"}});
247+
}
248+
249+
// ITS Hit
250+
template <typename T>
251+
bool hasITSHit(T const& track, int layer)
252+
{
253+
int ibit = layer - 1;
254+
return (track.itsClusterMap() & (1 << ibit));
242255
}
243256

244257
// Single-Track Selection for Particles inside Jets
@@ -247,12 +260,14 @@ struct NucleiInJets {
247260
{
248261
if (!track.hasITS())
249262
return false;
250-
if (track.itsNCls() < 3)
263+
if ((!hasITSHit(track, 1)) && (!hasITSHit(track, 2)) && (!hasITSHit(track, 3)))
251264
return false;
252265
if (!track.hasTPC())
253266
return false;
254267
if (track.tpcNClsCrossedRows() < 70)
255268
return false;
269+
if (track.tpcNClsCrossedRows() / track.tpcNClsFindable() < 0.8)
270+
return false;
256271
if (track.tpcChi2NCl() > 4)
257272
return false;
258273
if (track.itsChi2NCl() > 36)
@@ -261,7 +276,12 @@ struct NucleiInJets {
261276
return false;
262277
if (track.pt() < 0.15)
263278
return false;
279+
if (std::fabs(track.dcaXY()) > 0.25)
280+
return false;
281+
if (std::fabs(track.dcaZ()) > 2.0)
282+
return false;
264283

284+
/*
265285
// pt-dependent selection
266286
if (setDCAselectionPtDep) {
267287
if (std::fabs(track.dcaXY()) > (par0 + par1 / track.pt()))
@@ -277,6 +297,7 @@ struct NucleiInJets {
277297
if (std::fabs(track.dcaZ()) > maxDcaz)
278298
return false;
279299
}
300+
*/
280301

281302
return true;
282303
}
@@ -489,14 +510,14 @@ struct NucleiInJets {
489510

490511
do {
491512
double dijMin(1e+06), diBmin(1e+06);
492-
int iMin(0), jMin(0), iB_min(0);
513+
int iMin(0), jMin(0), iBmin(0);
493514
for (int i = 0; i < static_cast<int>(trk.size()); i++) { // o2-linter: disable=[const-ref-in-for-loop]
494515
if (trk[i].Mag() == 0)
495516
continue;
496517
double diB = 1.0 / (trk[i].Pt() * trk[i].Pt());
497518
if (diB < diBmin) {
498519
diBmin = diB;
499-
iB_min = i;
520+
iBmin = i;
500521
}
501522
for (int j = (i + 1); j < static_cast<int>(trk.size()); j++) { // o2-linter: disable=[const-ref-in-for-loop]
502523
if (trk[j].Mag() == 0)
@@ -515,8 +536,8 @@ struct NucleiInJets {
515536
nParticlesRemoved++;
516537
}
517538
if (dijMin > diBmin) {
518-
jet.push_back(trk[iB_min]);
519-
trk[iB_min].SetXYZ(0, 0, 0);
539+
jet.push_back(trk[iBmin]);
540+
trk[iBmin].SetXYZ(0, 0, 0);
520541
nParticlesRemoved++;
521542
}
522543
} while (nParticlesRemoved < static_cast<int>(trk.size()));
@@ -969,14 +990,14 @@ struct NucleiInJets {
969990

970991
do {
971992
double dijMin(1e+06), diBmin(1e+06);
972-
int iMin(0), jMin(0), iB_min(0);
993+
int iMin(0), jMin(0), iBmin(0);
973994
for (int i = 0; i < static_cast<int>(trk.size()); i++) { // o2-linter: disable=[const-ref-in-for-loop]
974995
if (trk[i].Mag() == 0)
975996
continue;
976997
double diB = 1.0 / (trk[i].Pt() * trk[i].Pt());
977998
if (diB < diBmin) {
978999
diBmin = diB;
979-
iB_min = i;
1000+
iBmin = i;
9801001
}
9811002
for (int j = (i + 1); j < static_cast<int>(trk.size()); j++) { // o2-linter: disable=[const-ref-in-for-loop]
9821003
if (trk[j].Mag() == 0)
@@ -995,8 +1016,8 @@ struct NucleiInJets {
9951016
nParticlesRemoved++;
9961017
}
9971018
if (dijMin > diBmin) {
998-
jet.push_back(trk[iB_min]);
999-
trk[iB_min].SetXYZ(0, 0, 0);
1019+
jet.push_back(trk[iBmin]);
1020+
trk[iBmin].SetXYZ(0, 0, 0);
10001021
nParticlesRemoved++;
10011022
}
10021023
} while (nParticlesRemoved < static_cast<int>(trk.size()));
@@ -1181,14 +1202,14 @@ struct NucleiInJets {
11811202

11821203
do {
11831204
double dijMin(1e+06), diBmin(1e+06);
1184-
int iMin(0), jMin(0), iB_min(0);
1205+
int iMin(0), jMin(0), iBmin(0);
11851206
for (int i = 0; i < static_cast<int>(trk.size()); i++) { // o2-linter: disable=[const-ref-in-for-loop]
11861207
if (trk[i].Mag() == 0)
11871208
continue;
11881209
double diB = 1.0 / (trk[i].Pt() * trk[i].Pt());
11891210
if (diB < diBmin) {
11901211
diBmin = diB;
1191-
iB_min = i;
1212+
iBmin = i;
11921213
}
11931214
for (int j = (i + 1); j < static_cast<int>(trk.size()); j++) { // o2-linter: disable=[const-ref-in-for-loop]
11941215
if (trk[j].Mag() == 0)
@@ -1207,8 +1228,8 @@ struct NucleiInJets {
12071228
nParticlesRemoved++;
12081229
}
12091230
if (dijMin > diBmin) {
1210-
jet.push_back(trk[iB_min]);
1211-
trk[iB_min].SetXYZ(0, 0, 0);
1231+
jet.push_back(trk[iBmin]);
1232+
trk[iBmin].SetXYZ(0, 0, 0);
12121233
nParticlesRemoved++;
12131234
}
12141235
} while (nParticlesRemoved < static_cast<int>(trk.size()));
@@ -1340,6 +1361,170 @@ struct NucleiInJets {
13401361
}
13411362
}
13421363
PROCESS_SWITCH(NucleiInJets, processAntiprotonReweighting, "Process antiproton reweighting", false);
1364+
1365+
void processGhosts(SelectedCollisions::iterator const& collision, FullNucleiTracks const& tracks)
1366+
{
1367+
// Event Selection
1368+
if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx)
1369+
return;
1370+
1371+
// Track Selection for Jet Finding
1372+
std::vector<TVector3> trk;
1373+
for (auto track : tracks) { // o2-linter: disable=[const-ref-in-for-loop]
1374+
1375+
if (!passedTrackSelectionForJetReconstruction(track))
1376+
continue;
1377+
TVector3 momentum(track.px(), track.py(), track.pz());
1378+
trk.push_back(momentum);
1379+
}
1380+
// int nTracks = static_cast<int>(trk.size());
1381+
1382+
// Generate Ghosts
1383+
for (int i = 0; i < nGhosts; i++) { // o2-linter: disable=[const-ref-in-for-loop]
1384+
1385+
double eta = gRandom->Uniform(-0.8, 0.8);
1386+
double phi = gRandom->Uniform(0.0, TwoPI);
1387+
double pt = 1e-100;
1388+
TVector3 ghost;
1389+
ghost.SetPtEtaPhi(pt, eta, phi);
1390+
trk.push_back(ghost);
1391+
}
1392+
1393+
// Anti-kt Jet Finder
1394+
int nParticlesRemoved(0);
1395+
std::vector<TVector3> jet;
1396+
std::vector<double> jetArea;
1397+
1398+
do {
1399+
double dijMin(1e+06), diBmin(1e+06);
1400+
int iMin(0), jMin(0), iBmin(0);
1401+
int nGhostsInJet(0);
1402+
for (int i = 0; i < static_cast<int>(trk.size()); i++) { // o2-linter: disable=[const-ref-in-for-loop]
1403+
if (trk[i].Mag() == 0)
1404+
continue;
1405+
double diB = 1.0 / (trk[i].Pt() * trk[i].Pt());
1406+
if (diB < diBmin) {
1407+
diBmin = diB;
1408+
iBmin = i;
1409+
}
1410+
for (int j = (i + 1); j < static_cast<int>(trk.size()); j++) { // o2-linter: disable=[const-ref-in-for-loop]
1411+
if (trk[j].Mag() == 0)
1412+
continue;
1413+
double dij = calculateDij(trk[i], trk[j], rJet);
1414+
if (dij < dijMin) {
1415+
dijMin = dij;
1416+
iMin = i;
1417+
jMin = j;
1418+
}
1419+
}
1420+
}
1421+
if (dijMin < diBmin) {
1422+
if (trk[iMin].Pt() == 1e-100)
1423+
nGhostsInJet++;
1424+
if (trk[jMin].Pt() == 1e-100)
1425+
nGhostsInJet++;
1426+
trk[iMin] = trk[iMin] + trk[jMin];
1427+
trk[jMin].SetXYZ(0, 0, 0);
1428+
nParticlesRemoved++;
1429+
}
1430+
if (dijMin > diBmin) {
1431+
double area = (static_cast<double>(nGhostsInJet) / static_cast<double>(nGhosts)) * TwoPI * 1.6;
1432+
jetArea.push_back(area);
1433+
jet.push_back(trk[iBmin]);
1434+
trk[iBmin].SetXYZ(0, 0, 0);
1435+
nParticlesRemoved++;
1436+
}
1437+
} while (nParticlesRemoved < static_cast<int>(trk.size()));
1438+
1439+
for (int i = 0; i < static_cast<int>(jet.size()); i++) { // o2-linter: disable=[const-ref-in-for-loop]
1440+
1441+
if ((std::fabs(jet[i].Eta()) + rJet) > maxEta)
1442+
continue;
1443+
registryQC.fill(HIST("hJetArea"), jetArea[i]);
1444+
}
1445+
}
1446+
PROCESS_SWITCH(NucleiInJets, processGhosts, "Process Ghosts", false);
1447+
1448+
void processDetResponseMatrix(SimCollisions const& collisions, MCTracks const& mcTracks, aod::McCollisions const&, const aod::McParticles&)
1449+
{
1450+
for (const auto& collision : collisions) { // o2-linter: disable=[const-ref-in-for-loop]
1451+
1452+
// Event Selection
1453+
if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx)
1454+
continue;
1455+
1456+
// List of Tracks and Particles
1457+
std::vector<TVector3> trk;
1458+
std::vector<TVector3> part;
1459+
auto tracksPerColl = mcTracks.sliceBy(perCollision, collision.globalIndex());
1460+
1461+
for (auto track : tracksPerColl) { // o2-linter: disable=[const-ref-in-for-loop]
1462+
1463+
if (!passedTrackSelectionForJetReconstruction(track))
1464+
continue;
1465+
if (!track.has_mcParticle())
1466+
continue;
1467+
const auto particle = track.mcParticle();
1468+
1469+
TVector3 recMomentum(track.px(), track.py(), track.pz());
1470+
TVector3 genMomentum(particle.px(), particle.py(), particle.pz());
1471+
trk.push_back(recMomentum);
1472+
part.push_back(genMomentum);
1473+
}
1474+
1475+
// Anti-kt Jet Finder
1476+
int nParticlesRemoved(0);
1477+
std::vector<TVector3> jetRecMomentum;
1478+
std::vector<TVector3> jetGenMomentum;
1479+
1480+
do {
1481+
double dijMin(1e+06), diBmin(1e+06);
1482+
int iMin(0), jMin(0), iBmin(0);
1483+
for (int i = 0; i < static_cast<int>(trk.size()); i++) { // o2-linter: disable=[const-ref-in-for-loop]
1484+
if (trk[i].Mag() == 0)
1485+
continue;
1486+
double diB = 1.0 / (trk[i].Pt() * trk[i].Pt());
1487+
if (diB < diBmin) {
1488+
diBmin = diB;
1489+
iBmin = i;
1490+
}
1491+
for (int j = (i + 1); j < static_cast<int>(trk.size()); j++) { // o2-linter: disable=[const-ref-in-for-loop]
1492+
if (trk[j].Mag() == 0)
1493+
continue;
1494+
double dij = calculateDij(trk[i], trk[j], rJet);
1495+
if (dij < dijMin) {
1496+
dijMin = dij;
1497+
iMin = i;
1498+
jMin = j;
1499+
}
1500+
}
1501+
}
1502+
if (dijMin < diBmin) {
1503+
trk[iMin] = trk[iMin] + trk[jMin];
1504+
part[iMin] = part[iMin] + part[jMin];
1505+
trk[jMin].SetXYZ(0, 0, 0);
1506+
nParticlesRemoved++;
1507+
}
1508+
if (dijMin > diBmin) {
1509+
jetRecMomentum.push_back(trk[iBmin]);
1510+
jetGenMomentum.push_back(part[iBmin]);
1511+
trk[iBmin].SetXYZ(0, 0, 0);
1512+
nParticlesRemoved++;
1513+
}
1514+
} while (nParticlesRemoved < static_cast<int>(trk.size()));
1515+
1516+
for (int i = 0; i < static_cast<int>(jetRecMomentum.size()); i++) { // o2-linter: disable=[const-ref-in-for-loop]
1517+
1518+
if ((std::fabs(jetRecMomentum[i].Eta()) + rJet) > maxEta)
1519+
continue;
1520+
1521+
double ptGen = jetGenMomentum[i].Pt();
1522+
double ptRec = jetRecMomentum[i].Pt();
1523+
registryMC.fill(HIST("detectorResponseMatrix"), ptGen, ptRec);
1524+
}
1525+
}
1526+
}
1527+
PROCESS_SWITCH(NucleiInJets, processDetResponseMatrix, "process detector response matrix", false);
13431528
};
13441529

13451530
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)