Skip to content

Commit 4e5fda7

Browse files
alcalivaalcaliva
authored
[PWGLF] updated jet selection (#9778)
1 parent ce050cb commit 4e5fda7

File tree

1 file changed

+65
-52
lines changed

1 file changed

+65
-52
lines changed

‎PWGLF/Tasks/Nuspex/nucleiInJets.cxx‎

Lines changed: 65 additions & 52 deletions
Original file line numberDiff line numberDiff line change
@@ -98,6 +98,8 @@ struct NucleiInJets {
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"};
100100
Configurable<int> nGhosts{"nGhosts", 1000, "number of ghost particles"};
101+
Configurable<double> alpha{"alpha", 1.0, "Alpha"};
102+
Configurable<double> averagePtUE{"averagePtUE", 0.1, "Average pt of UE"};
101103

102104
// Track Parameters
103105
Configurable<double> par0{"par0", 0.00164, "par 0"};
@@ -162,10 +164,12 @@ struct NucleiInJets {
162164
registryQC.add("sumPtUE", "sumPtUE", HistType::kTH1F, {{500, 0, 50, "#it{p}_{T} (GeV/#it{c})"}});
163165
registryQC.add("nJets_found", "nJets_found", HistType::kTH1F, {{10, 0, 10, "#it{n}_{Jet}"}});
164166
registryQC.add("nJets_selected", "nJets_selected", HistType::kTH1F, {{10, 0, 10, "#it{n}_{Jet}"}});
165-
registryQC.add("event_selection_jets", "event_selection_jets", HistType::kTH1F, {{10, 0, 10, "counter"}});
166167
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)"}});
167168
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)"}});
169+
registryQC.add("jet_jet_overlaps", "jet_jet_overlaps", HistType::kTH2F, {{20, 0.0, 20.0, "#it{n}_{jet}"}, {200, 0.0, 200.0, "#it{n}_{overlaps}"}});
168170
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}"}});
171+
registryQC.add("ue_ue_overlaps", "ue_ue_overlaps", HistType::kTH2F, {{20, 0.0, 20.0, "#it{n}_{jet}"}, {200, 0.0, 200.0, "#it{n}_{overlaps}"}});
172+
registryQC.add("tot_overlaps", "tot_overlaps", HistType::kTH2F, {{20, 0.0, 20.0, "#it{n}_{jet}"}, {200, 0.0, 200.0, "#it{n}_{overlaps}"}});
169173
registryQC.add("hJetArea", "hJetArea", HistType::kTH1F, {{450, 0, 15, "Area"}});
170174

171175
// Event Counters
@@ -243,7 +247,9 @@ struct NucleiInJets {
243247
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}"}});
244248

245249
// 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})"}});
250+
if (doprocessDetResponseMatrix) {
251+
registryMC.add("detectorResponseMatrix", "detectorResponseMatrix", HistType::kTH2F, {{5000, 0.0, 50.0, "#it{p}_{T}^{gen} (GeV/#it{c})"}, {5000, 0.0, 50.0, "#it{p}_{T}^{rec} (GeV/#it{c})"}});
252+
}
247253
}
248254

249255
// ITS Hit
@@ -443,6 +449,12 @@ struct NucleiInJets {
443449
return false;
444450
}
445451

452+
double getCorrectedPt(double ptRec)
453+
{
454+
// to be developed
455+
return ptRec;
456+
}
457+
446458
void getReweightingHistograms(o2::framework::Service<o2::ccdb::BasicCCDBManager> const& ccdbObj, TString filepath, TString histname_antip_jet, TString histname_antip_ue)
447459
{
448460
TList* l = ccdbObj->get<TList>(filepath.Data());
@@ -469,7 +481,6 @@ struct NucleiInJets {
469481
{
470482
// Event Counter: before event selection
471483
registryData.fill(HIST("number_of_events_data"), 0.5);
472-
registryQC.fill(HIST("event_selection_jets"), 0.5); // all events before jet selection
473484

474485
// Event Selection
475486
if (!collision.sel8())
@@ -490,6 +501,7 @@ struct NucleiInJets {
490501

491502
// List of Tracks
492503
std::vector<TVector3> trk;
504+
std::vector<int> ntrk;
493505

494506
for (auto track : tracks) { // o2-linter: disable=[const-ref-in-for-loop]
495507

@@ -500,13 +512,20 @@ struct NucleiInJets {
500512

501513
TVector3 momentum(track.px(), track.py(), track.pz());
502514
trk.push_back(momentum);
515+
ntrk.push_back(1);
503516
}
504517

518+
// Reject Empty Events
519+
if (static_cast<int>(trk.size()) < 1)
520+
return;
521+
registryData.fill(HIST("number_of_events_data"), 3.5);
522+
505523
// Anti-kt Jet Finder
506524
int nParticlesRemoved(0);
507525
std::vector<TVector3> jet;
508526
std::vector<TVector3> ue1;
509527
std::vector<TVector3> ue2;
528+
std::vector<int> nParticlesInjet;
510529

511530
do {
512531
double dijMin(1e+06), diBmin(1e+06);
@@ -532,11 +551,14 @@ struct NucleiInJets {
532551
}
533552
if (dijMin < diBmin) {
534553
trk[iMin] = trk[iMin] + trk[jMin];
554+
ntrk[iMin] = ntrk[iMin] + ntrk[jMin];
535555
trk[jMin].SetXYZ(0, 0, 0);
556+
ntrk[jMin] = 0;
536557
nParticlesRemoved++;
537558
}
538559
if (dijMin > diBmin) {
539560
jet.push_back(trk[iBmin]);
561+
nParticlesInjet.push_back(ntrk[iBmin]);
540562
trk[iBmin].SetXYZ(0, 0, 0);
541563
nParticlesRemoved++;
542564
}
@@ -546,14 +568,16 @@ struct NucleiInJets {
546568

547569
// Jet Selection
548570
std::vector<int> isSelected;
549-
for (int i = 0; i < static_cast<int>(jet.size()); i++) { // o2-linter: disable=[const-ref-in-for-loop]
550-
isSelected.push_back(0);
551-
}
552-
553571
int nJetsSelected(0);
554572
for (int i = 0; i < static_cast<int>(jet.size()); i++) { // o2-linter: disable=[const-ref-in-for-loop]
555573

556-
if ((std::fabs(jet[i].Eta()) + rJet) > maxEta)
574+
// Initialization
575+
isSelected.push_back(0);
576+
577+
// Jet fully contained inside acceptance
578+
if ((std::fabs(jet[i].Eta()) + rJet) > (maxEta - 0.5))
579+
continue;
580+
if (nParticlesInjet[i] < minNparticlesInJet)
557581
continue;
558582

559583
// Perpendicular cones
@@ -564,93 +588,82 @@ struct NucleiInJets {
564588
ue1.push_back(ueAxis1);
565589
ue2.push_back(ueAxis2);
566590

567-
double nPartJetPlusUE(0);
568-
double nPartJet(0);
569-
double nPartUE(0);
570-
double ptJetPlusUE(0);
571-
double ptJet(0);
572591
double ptUE(0);
573-
574592
for (auto track : tracks) { // o2-linter: disable=[const-ref-in-for-loop]
575593

576594
if (!passedTrackSelectionForJetReconstruction(track))
577595
continue;
578596
TVector3 selectedTrack(track.px(), track.py(), track.pz());
579597

580-
double deltaEtaJet = selectedTrack.Eta() - jet[i].Eta();
581-
double deltaPhiJet = getDeltaPhi(selectedTrack.Phi(), jet[i].Phi());
582-
double deltaRjet = std::sqrt(deltaEtaJet * deltaEtaJet + deltaPhiJet * deltaPhiJet);
583598
double deltaEtaUe1 = selectedTrack.Eta() - ueAxis1.Eta();
584599
double deltaPhiUe1 = getDeltaPhi(selectedTrack.Phi(), ueAxis1.Phi());
585600
double deltaRUe1 = std::sqrt(deltaEtaUe1 * deltaEtaUe1 + deltaPhiUe1 * deltaPhiUe1);
586601
double deltaEtaUe2 = selectedTrack.Eta() - ueAxis2.Eta();
587602
double deltaPhiUe2 = getDeltaPhi(selectedTrack.Phi(), ueAxis2.Phi());
588603
double deltaRUe2 = std::sqrt(deltaEtaUe2 * deltaEtaUe2 + deltaPhiUe2 * deltaPhiUe2);
589604

590-
if (deltaRjet < rJet) {
591-
registryQC.fill(HIST("deltaEtadeltaPhiJet"), deltaEtaJet, deltaPhiJet);
592-
nPartJetPlusUE++;
593-
ptJetPlusUE = ptJetPlusUE + selectedTrack.Pt();
594-
}
595-
if (deltaRUe1 < rJet) {
605+
if (deltaRUe1 < alpha * rJet) {
596606
registryQC.fill(HIST("deltaEtadeltaPhi_ue"), deltaEtaUe1, deltaPhiUe1);
597-
nPartUE++;
598607
ptUE = ptUE + selectedTrack.Pt();
599608
}
600-
if (deltaRUe2 < rJet) {
609+
if (deltaRUe2 < alpha * rJet) {
601610
registryQC.fill(HIST("deltaEtadeltaPhi_ue"), deltaEtaUe2, deltaPhiUe2);
602-
nPartUE++;
603611
ptUE = ptUE + selectedTrack.Pt();
604612
}
605613
}
606-
nPartJet = nPartJetPlusUE - 0.5 * nPartUE;
607-
ptJet = ptJetPlusUE - 0.5 * ptUE;
608-
registryQC.fill(HIST("NchJetPlusUE"), nPartJetPlusUE);
609-
registryQC.fill(HIST("NchJet"), nPartJet);
610-
registryQC.fill(HIST("NchUE"), nPartUE);
611-
registryQC.fill(HIST("sumPtJetPlusUE"), ptJetPlusUE);
612-
registryQC.fill(HIST("sumPtJet"), ptJet);
613-
registryQC.fill(HIST("sumPtUE"), ptUE);
614-
615-
if (ptJet < minJetPt)
616-
continue;
617-
if (nPartJetPlusUE < minNparticlesInJet)
614+
registryQC.fill(HIST("sumPtUE"), 0.5 * ptUE);
615+
registryQC.fill(HIST("NchJetPlusUE"), nParticlesInjet[i]);
616+
617+
double ptJetRec = jet[i].Pt() - averagePtUE;
618+
double ptJetCorr = getCorrectedPt(ptJetRec);
619+
620+
if (ptJetCorr < minJetPt)
618621
continue;
622+
619623
nJetsSelected++;
620624
isSelected[i] = 1;
621625
}
622626
registryQC.fill(HIST("nJets_selected"), nJetsSelected);
623627

624628
if (nJetsSelected == 0)
625629
return;
626-
registryData.fill(HIST("number_of_events_data"), 3.5);
627-
registryQC.fill(HIST("event_selection_jets"), 1.5); // events with pTjet>10 GeV/c selected
628-
//************************************************************************************************************************************
630+
registryData.fill(HIST("number_of_events_data"), 4.5);
629631

630632
// Overlaps
631-
int nOverlaps(0);
632-
for (int i = 0; i < static_cast<int>(jet.size()); i++) { // o2-linter: disable=[const-ref-in-for-loop]
633+
int nOverlapsJetJet(0);
634+
int nOverlapsJetUe(0);
635+
int nOverlapsUeUe(0);
636+
int nOverlapsTot(0);
637+
for (int i = 0; i < static_cast<int>(jet.size()); i++) {
633638
if (isSelected[i] == 0)
634639
continue;
635640

636-
for (int j = 0; j < static_cast<int>(jet.size()); j++) { // o2-linter: disable=[const-ref-in-for-loop]
637-
if (isSelected[j] == 0 || i == j)
641+
for (int j = (i + 1); j < static_cast<int>(jet.size()); j++) {
642+
if (isSelected[j] == 0)
638643
continue;
639-
if (overlap(jet[i], ue1[j], rJet) || overlap(jet[i], ue2[j], rJet) || overlap(jet[i], jet[j], rJet))
640-
nOverlaps++;
644+
if (overlap(jet[i], jet[j], rJet))
645+
nOverlapsJetJet++;
646+
if (overlap(jet[i], ue1[j], rJet) || overlap(jet[i], ue2[j], rJet))
647+
nOverlapsJetUe++;
648+
if (overlap(ue1[i], ue1[j], rJet) || overlap(ue1[i], ue2[j], rJet) || overlap(ue2[i], ue2[j], rJet))
649+
nOverlapsUeUe++;
641650
}
642651
}
643-
registryQC.fill(HIST("jet_ue_overlaps"), nJetsSelected, nOverlaps);
652+
nOverlapsTot = nOverlapsJetJet + nOverlapsJetUe + nOverlapsUeUe;
653+
registryQC.fill(HIST("jet_jet_overlaps"), nJetsSelected, nOverlapsJetJet);
654+
registryQC.fill(HIST("jet_ue_overlaps"), nJetsSelected, nOverlapsJetUe);
655+
registryQC.fill(HIST("ue_ue_overlaps"), nJetsSelected, nOverlapsUeUe);
656+
registryQC.fill(HIST("tot_overlaps"), nJetsSelected, nOverlapsTot);
644657

645658
if (nJetsSelected > nJetsPerEventMax)
646659
return;
647-
registryData.fill(HIST("number_of_events_data"), 4.5);
660+
registryData.fill(HIST("number_of_events_data"), 5.5);
648661

649-
if (requireNoOverlap && nOverlaps > 0)
662+
if (requireNoOverlap && nOverlapsTot > 0)
650663
return;
651-
registryData.fill(HIST("number_of_events_data"), 5.5);
664+
registryData.fill(HIST("number_of_events_data"), 6.5);
652665

653-
//************************************************************************************************************************************
666+
//**************************************************************************************************************
654667

655668
for (int i = 0; i < static_cast<int>(jet.size()); i++) { // o2-linter: disable=[const-ref-in-for-loop]
656669

0 commit comments

Comments
 (0)