Skip to content

Commit 9052156

Browse files
alcalivaalcaliva
authored
modified UE calculation + small improvements (#4106)
1 parent e5e20c3 commit 9052156

File tree

1 file changed

+67
-44
lines changed

1 file changed

+67
-44
lines changed

‎PWGLF/Tasks/nuclei_in_jets.cxx‎

Lines changed: 67 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -104,13 +104,14 @@ struct nuclei_in_jets {
104104
void init(InitContext const&)
105105
{
106106
// Global Properties and QC
107-
registryQC.add("number_of_events_data", "number of events in data", HistType::kTH1F, {{5, 0, 5, "1 = all events, 2 = selected events, 3 = events with pt>pt_threshold, 4 = events with pt>pt_threshold and particle of interest"}});
108-
registryQC.add("jet_plus_ue_multiplicity", "jet + underlying-event multiplicity", HistType::kTH1F, {{300, 0, 300, "#it{N}_{ch}"}});
109-
registryQC.add("jet_multiplicity", "jet multiplicity", HistType::kTH1F, {{300, 0, 300, "#it{N}_{ch}"}});
110-
registryQC.add("ue_multiplicity", "underlying-event multiplicity", HistType::kTH1F, {{300, 0, 300, "#it{N}_{ch}"}});
107+
registryQC.add("number_of_events_data", "number of events in data", HistType::kTH1F, {{10, 0, 10, "counter"}});
108+
registryQC.add("jet_plus_ue_multiplicity", "jet + underlying-event multiplicity", HistType::kTH1F, {{100, 0, 100, "#it{N}_{ch}"}});
109+
registryQC.add("jet_multiplicity", "jet multiplicity", HistType::kTH1F, {{100, 0, 100, "#it{N}_{ch}"}});
110+
registryQC.add("ue_multiplicity", "underlying-event multiplicity", HistType::kTH1F, {{100, 0, 100, "#it{N}_{ch}"}});
111111
registryQC.add("pt_leading", "pt leading", HistType::kTH1F, {{500, 0, 50, "#it{p}_{T} (GeV/#it{c})"}});
112112
registryQC.add("eta_phi_jet", "DeltaEta DeltaPhi jet", HistType::kTH2F, {{500, -0.5, 0.5, "#Delta#eta"}, {500, 0.0, TMath::Pi(), "#Delta#phi"}});
113113
registryQC.add("eta_phi_ue", "DeltaEta DeltaPhi UE", HistType::kTH2F, {{500, -0.5, 0.5, "#Delta#eta"}, {500, 0.0, TMath::Pi(), "#Delta#phi"}});
114+
registryQC.add("r_max_jet", "R Max jet", HistType::kTH1F, {{400, 0.0, 0.8, "#it{R}_{max}"}});
114115
registryQC.add("r_jet", "R jet", HistType::kTH1F, {{400, 0.0, 0.8, "#it{R}"}});
115116
registryQC.add("r_ue", "R ue", HistType::kTH1F, {{400, 0.0, 0.8, "#it{R}"}});
116117

@@ -249,6 +250,10 @@ struct nuclei_in_jets {
249250
// Process Data
250251
void processData(SelectedCollisions::iterator const& collision, FullTracks const& tracks)
251252
{
253+
254+
// Seed
255+
gRandom->SetSeed(0);
256+
252257
// Event Counter (before event selection)
253258
registryQC.fill(HIST("number_of_events_data"), 0.5);
254259

@@ -299,6 +304,7 @@ struct nuclei_in_jets {
299304
// Skip Events with no trigger Particle
300305
if (pt_max == 0)
301306
return;
307+
registryQC.fill(HIST("number_of_events_data"), 2.5);
302308

303309
// Histogram with pt_leading
304310
registryQC.fill(HIST("pt_leading"), pt_max);
@@ -309,26 +315,25 @@ struct nuclei_in_jets {
309315
// Selection of Events with pt > pt_leading
310316
if (nParticles < 2)
311317
return;
318+
registryQC.fill(HIST("number_of_events_data"), 3.5);
319+
312320
if (pt_max < min_pt_leading)
313321
return;
314322

315323
// Event Counter (after pt > pt_max selection)
316-
registryQC.fill(HIST("number_of_events_data"), 2.5);
324+
registryQC.fill(HIST("number_of_events_data"), 4.5);
317325

318326
// Skip Events with no Particle of Interest
319327
if (!containsParticleOfInterest)
320328
return;
321329

322330
// Event Counter (events with pt > pt_max that contain particle of interest)
323-
registryQC.fill(HIST("number_of_events_data"), 3.5);
331+
registryQC.fill(HIST("number_of_events_data"), 5.5);
324332

325333
// Momentum of the Leading Particle
326334
auto const& leading_track = tracks.iteratorAt(leading_ID);
327335
TVector3 p_leading(leading_track.px(), leading_track.py(), leading_track.pz());
328336

329-
// Instruction to be removed
330-
registryQC.fill(HIST("number_of_events_data"), 4.5);
331-
332337
// Array of Particles inside Jet
333338
std::vector<int> jet_particle_ID;
334339
jet_particle_ID.push_back(leading_ID);
@@ -405,16 +410,49 @@ struct nuclei_in_jets {
405410
// Multiplicity inside Jet + UE
406411
int nParticlesJetUE = static_cast<int>(jet_particle_ID.size());
407412

413+
// Find Maximum Distance from Jet Axis
414+
float Rmax(0);
415+
416+
for (int i = 0; i < nParticlesJetUE; i++) {
417+
418+
const auto& jet_track = tracks.iteratorAt(jet_particle_ID[i]);
419+
TVector3 p_i(jet_track.px(), jet_track.py(), jet_track.pz());
420+
421+
// Track Selection
422+
if (!passedTrackSelection(jet_track))
423+
continue;
424+
if (require_primVtx_contributor && !(jet_track.isPVContributor()))
425+
continue;
426+
427+
float deltaEta = p_i.Eta() - p_leading.Eta();
428+
float deltaPhi = TVector2::Phi_0_2pi(p_i.Phi() - p_leading.Phi());
429+
float R = TMath::Sqrt(deltaEta * deltaEta + deltaPhi * deltaPhi);
430+
if (R > Rmax)
431+
Rmax = R;
432+
}
433+
434+
registryQC.fill(HIST("r_max_jet"), Rmax);
435+
436+
// Cut on eta jet
437+
float eta_jet_axis = p_leading.Eta();
438+
if ((TMath::Abs(eta_jet_axis) + Rmax) > max_eta)
439+
return;
440+
registryQC.fill(HIST("number_of_events_data"), 6.5);
441+
408442
// Fill Jet Multiplicity
409-
registryQC.fill(HIST("jet_plus_ue_multiplicity"), static_cast<float>(jet_particle_ID.size()));
443+
registryQC.fill(HIST("jet_plus_ue_multiplicity"), nParticlesJetUE);
410444

411445
// Perpendicular Cones for UE Estimate
412-
TVector3 z_positive(0.0, 0.0, 1.0);
413-
TVector3 z_negative(0.0, 0.0, -1.0);
414-
TVector3 v1 = (z_positive.Cross(p_leading)).Unit();
415-
TVector3 v2 = (z_negative.Cross(p_leading)).Unit();
416-
TVector3 v3 = (p_leading.Cross(v1)).Unit();
417-
TVector3 v4 = (p_leading.Cross(v2)).Unit();
446+
TVector3 z(0.0, 0.0, p_leading.Mag());
447+
TVector3 ue_axis = z.Cross(p_leading);
448+
449+
double dEta(200);
450+
do {
451+
double angle = gRandom->Uniform(0.0, TMath::TwoPi());
452+
ue_axis.Rotate(angle, p_leading);
453+
double eta_ue_axis = ue_axis.Eta();
454+
dEta = (TMath::Abs(eta_ue_axis) + Rmax);
455+
} while (dEta > max_eta);
418456

419457
// Store UE
420458
std::vector<int> ue_particle_ID;
@@ -429,43 +467,26 @@ struct nuclei_in_jets {
429467
const auto& ue_track = tracks.iteratorAt(particle_ID[i]);
430468

431469
// Variables
432-
float deltaEta1 = ue_track.eta() - v1.Eta();
433-
float deltaEta2 = ue_track.eta() - v2.Eta();
434-
float deltaEta3 = ue_track.eta() - v3.Eta();
435-
float deltaEta4 = ue_track.eta() - v4.Eta();
436-
437-
float deltaPhi1 = TVector2::Phi_0_2pi(ue_track.phi() - v1.Phi());
438-
float deltaPhi2 = TVector2::Phi_0_2pi(ue_track.phi() - v2.Phi());
439-
float deltaPhi3 = TVector2::Phi_0_2pi(ue_track.phi() - v3.Phi());
440-
float deltaPhi4 = TVector2::Phi_0_2pi(ue_track.phi() - v4.Phi());
441-
442-
float dr1 = TMath::Sqrt(deltaEta1 * deltaEta1 + deltaPhi1 * deltaPhi1);
443-
float dr2 = TMath::Sqrt(deltaEta2 * deltaEta2 + deltaPhi2 * deltaPhi2);
444-
float dr3 = TMath::Sqrt(deltaEta3 * deltaEta3 + deltaPhi3 * deltaPhi3);
445-
float dr4 = TMath::Sqrt(deltaEta4 * deltaEta4 + deltaPhi4 * deltaPhi4);
446-
447-
registryQC.fill(HIST("eta_phi_ue"), deltaEta1, deltaPhi1);
448-
registryQC.fill(HIST("eta_phi_ue"), deltaEta2, deltaPhi2);
449-
registryQC.fill(HIST("eta_phi_ue"), deltaEta3, deltaPhi3);
450-
registryQC.fill(HIST("eta_phi_ue"), deltaEta4, deltaPhi4);
451-
registryQC.fill(HIST("r_ue"), TMath::Sqrt(deltaEta1 * deltaEta1 + deltaPhi1 * deltaPhi1));
452-
registryQC.fill(HIST("r_ue"), TMath::Sqrt(deltaEta2 * deltaEta2 + deltaPhi2 * deltaPhi2));
453-
registryQC.fill(HIST("r_ue"), TMath::Sqrt(deltaEta3 * deltaEta3 + deltaPhi3 * deltaPhi3));
454-
registryQC.fill(HIST("r_ue"), TMath::Sqrt(deltaEta4 * deltaEta4 + deltaPhi4 * deltaPhi4));
470+
float deltaEta = ue_track.eta() - ue_axis.Eta();
471+
float deltaPhi = TVector2::Phi_0_2pi(ue_track.phi() - ue_axis.Phi());
472+
float dr = TMath::Sqrt(deltaEta * deltaEta + deltaPhi * deltaPhi);
473+
474+
registryQC.fill(HIST("eta_phi_ue"), deltaEta, deltaPhi);
475+
registryQC.fill(HIST("r_ue"), dr);
455476

456477
// Store Particles in the UE
457-
if (dr1 < max_jet_radius || dr2 < max_jet_radius || dr3 < max_jet_radius || dr4 < max_jet_radius) {
478+
if (dr < Rmax) {
458479
ue_particle_ID.push_back(particle_ID[i]);
459480
}
460481
}
461482

462483
// UE Multiplicity
463484
int nParticlesUE = static_cast<int>(ue_particle_ID.size());
464485

465-
registryQC.fill(HIST("ue_multiplicity"), static_cast<float>(ue_particle_ID.size()) / 4.0);
486+
registryQC.fill(HIST("ue_multiplicity"), nParticlesUE);
466487

467488
// Jet Multiplicity
468-
float jet_Nch = static_cast<float>(jet_particle_ID.size()) - static_cast<float>(ue_particle_ID.size()) / 4.0;
489+
int jet_Nch = nParticlesJetUE - nParticlesUE;
469490
registryQC.fill(HIST("jet_multiplicity"), jet_Nch);
470491

471492
// Loop over particles inside Jet
@@ -476,8 +497,10 @@ struct nuclei_in_jets {
476497

477498
float deltaEta = p_i.Eta() - p_leading.Eta();
478499
float deltaPhi = TVector2::Phi_0_2pi(p_i.Phi() - p_leading.Phi());
479-
registryQC.fill(HIST("eta_phi_jet"), deltaEta, deltaPhi);
480-
registryQC.fill(HIST("r_jet"), TMath::Sqrt(deltaEta * deltaEta + deltaPhi * deltaPhi));
500+
if (deltaEta != 0 && deltaPhi != 0) {
501+
registryQC.fill(HIST("eta_phi_jet"), deltaEta, deltaPhi);
502+
registryQC.fill(HIST("r_jet"), TMath::Sqrt(deltaEta * deltaEta + deltaPhi * deltaPhi));
503+
}
481504

482505
// Track Selection
483506
if (!passedTrackSelection(jet_track))

0 commit comments

Comments
 (0)