Skip to content
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
61 changes: 31 additions & 30 deletions PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx
Original file line number Diff line number Diff line change
@@ -1,3 +1,3 @@
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
Expand All @@ -8,7 +8,7 @@
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
//

Check warning on line 11 in PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[doc/file]

Documentation for \author is missing, incorrect or misplaced.

Check warning on line 11 in PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[doc/file]

Documentation for \brief is missing, incorrect or misplaced.

Check warning on line 11 in PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[doc/file]

Documentation for \file is missing, incorrect or misplaced.
// author: Lars Jörgensen

#include <vector>
Expand Down Expand Up @@ -58,7 +58,7 @@
};

struct AngularCorrelationsInJets {
// Switches
// Switches
Configurable<bool> useRejectionCut{"useRejectionCut", true, "use nsigmaRejection for correlations"};
Configurable<bool> deuteronAnalysis{"deuteronAnalysis", true, "true [false]: analyse (anti)deuterons [(anti)helium-3]"};
Configurable<bool> measureYields{"measureYields", false, "measure yields"};
Expand All @@ -80,7 +80,7 @@
Configurable<float> maxChi2TPC{"maxChi2TPC", 4.0, "max chi2 per cluster TPC"};
Configurable<float> maxDCAxy{"maxDCAxy", 0.05, "max DCA to vertex xy"};
Configurable<float> maxDCAz{"maxDCAz", 0.05, "max DCA to vertex z"};
Configurable<float> maxEta{"maxEta", 0.8, "max pseudorapidity"}; // consider jet cone
Configurable<float> maxEta{"maxEta", 0.8, "max pseudorapidity"}; // consider jet cone
Configurable<float> deltaEtaEdge{"deltaEtaEdge", 0.05, "min eta distance of jet from acceptance edge"}; // consider jet cone
Configurable<float> minTrackPt{"minTrackPt", 0.3, "minimum track pT"};
Configurable<float> requirePVContributor{"requirePVContributor", false, "require track to be PV contributor"};
Expand Down Expand Up @@ -142,7 +142,6 @@
Configurable<float> kaonTPCnsigmaHighPt{"kaonTPCnsigmaHighPt", 3.0, "[kaon] max TPC nsigma with high pT"};
Configurable<float> kaonTOFnsigma{"kaonTOFnsigma", 3.0, "[kaon] max TOF nsigma"};


Configurable<float> nsigmaRejection{"nsigmaRejection", 1.0, "reject tracks with nsigma < nsigmaRejection for >1 species"};
Configurable<int> trackBufferSize{"trackBufferSize", 200, "Number of mixed-event tracks being stored"};

Expand Down Expand Up @@ -177,7 +176,7 @@
nabs(aod::track::eta) < 0.8f &&
aod::track::pt > 0.1f); // add more preliminary cuts to filter if possible
Filter collisionFilter = (nabs(aod::jcollision::posZ) < zVtx);
Filter jetTrackCuts = (nabs(aod::jtrack::eta) > maxEta/* && aod::jtrack::pt > minJetParticlePt */);
Filter jetTrackCuts = (nabs(aod::jtrack::eta) > maxEta /* && aod::jtrack::pt > minJetParticlePt */);
Filter jetFilter = (aod::jet::pt >= minJetPt && nabs(aod::jet::eta) < nabs(maxEta - aod::jet::r / 100.f));

Preslice<FullTracksRun2> perCollisionFullTracksRun2 = o2::aod::track::collisionId;
Expand Down Expand Up @@ -419,7 +418,7 @@
return false;
if (track.itsChi2NCl() > maxChi2ITS)
return false;
if (fabs(track.eta()) > maxEta)

Check warning on line 421 in PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix]

Use std:: prefix for names from the std namespace.
return false;
if (track.pt() < minTrackPt)
return false;
Expand Down Expand Up @@ -723,7 +722,7 @@
}

double phiToAxis = RecoDecay::constrainAngle(track.phi() - jetAxis.Phi(), 0);
double etaToAxis = track.eta() - jetAxis.Eta();
double etaToAxis = track.eta() - jetAxis.Eta();
double deltaPhi = RecoDecay::constrainAngle(phiToAxis - buffer.at(i).first, -constants::math::PIHalf);
double deltaEta = etaToAxis - buffer.at(i).second;

Expand Down Expand Up @@ -914,7 +913,7 @@

if (subtractedJetPerp.pt() < minJetPt) // cut on jet w/o bkg
return jetCounter;
if ((fabs(jet.eta()) + jetR) > (maxEta - deltaEtaEdge))

Check warning on line 916 in PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix]

Use std:: prefix for names from the std namespace.
return jetCounter;
jetCounter++;
registryData.fill(HIST("ptTotalSubJetPerp"), subtractedJetPerp.pt());
Expand Down Expand Up @@ -1070,7 +1069,7 @@
registryData.fill(HIST("ptJetAntiproton"), jetParticle.pt());
registryQC.fill(HIST("ptJetAntiprotonVsTotalJet"), jetParticle.pt(), subtractedJetPerp.pt());
registryData.fill(HIST("trackProtocol"), 6); // # antiprotons
} else if (isNucleus(jetParticle)) { // collect nuclei in jet
} else if (isNucleus(jetParticle)) { // collect nuclei in jet
registryData.fill(HIST("ptJetNuclei"), jetParticle.pt());
registryQC.fill(HIST("ptJetNucleiVsTotalJet"), jetParticle.pt(), subtractedJetPerp.pt());
registryData.fill(HIST("trackProtocol"), 8); // # nuclei
Expand All @@ -1087,7 +1086,8 @@
registryData.fill(HIST("trackProtocol"), 5); // # high purity protons
jetProtons.emplace_back(jetParticle);
registryData.fill(HIST("dcaZJetProton"), jetParticle.pt(), jetParticle.dcaZ());
} if (isAntiproton(jetParticle, true)) {
}
if (isAntiproton(jetParticle, true)) {
registryData.fill(HIST("trackProtocol"), 7); // # high purity antiprotons
jetAntiprotons.emplace_back(jetParticle);
registryData.fill(HIST("dcaZJetAntiproton"), jetParticle.pt(), jetParticle.dcaZ());
Expand Down Expand Up @@ -1301,14 +1301,14 @@
auto [rho, rhoM] = bkgSub.estimateRhoAreaMedian(jetInput, doSparse);
auto [rhoPerp, rhoMPerp] = bkgSub.estimateRhoPerpCone(jetInput, jets);

for (auto& jet : jets) {

Check warning on line 1304 in PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop]

Use constant references for non-modified iterators in range-based for loops.
if (!jet.has_constituents())
continue;
fastjet::PseudoJet subtractedJetPerp(0., 0., 0., 0.);
subtractedJetPerp = bkgSub.doRhoAreaSub(jet, rhoPerp, rhoMPerp);
fastjet::PseudoJet subtractedJetArea(0., 0., 0., 0.);
subtractedJetArea = bkgSub.doRhoAreaSub(jet, rho, rhoM);

if (subtractedJetPerp.pt() < minJetPt) // cut on jet w/o bkg
continue;
registryData.fill(HIST("ptTotalSubJetPerp"), subtractedJetPerp.pt());
Expand All @@ -1319,25 +1319,25 @@
registryQC.fill(HIST("rhoMEstimatePerp"), jet.pt(), rhoMPerp);
double jetBkgDeltaPt = jet.pt() - subtractedJetPerp.pt();
registryQC.fill(HIST("jetBkgDeltaPt"), jetBkgDeltaPt);

registryData.fill(HIST("eventProtocol"), 4);
std::vector<fastjet::PseudoJet> constituents = jet.constituents();

registryData.fill(HIST("eventProtocol"), 5);
registryData.fill(HIST("numberOfJets"), 0);
registryData.fill(HIST("ptTotalJet"), jet.pt());
registryData.fill(HIST("jetRapidity"), jet.rap());
registryData.fill(HIST("numPartInJet"), jet.constituents().size());
registryQC.fill(HIST("jetPtVsNumPart"), jet.pt(), jet.constituents().size());

double maxRadius = 0;
for (const auto& constituent : constituents) {
registryData.fill(HIST("ptJetParticle"), constituent.pt());
registryQC.fill(HIST("phiJet"), constituent.phi());
registryQC.fill(HIST("phiPtJet"), constituent.pt(), constituent.phi());
registryQC.fill(HIST("etaJet"), constituent.eta());
registryQC.fill(HIST("etaPtJet"), constituent.pt(), constituent.eta());

if (std::isnan(constituent.phi()) || std::isnan(jet.phi())) // geometric jet cone
continue;
double deltaPhi = RecoDecay::constrainAngle(constituent.phi() - jet.phi(), -constants::math::PIHalf);
Expand All @@ -1348,22 +1348,22 @@
maxRadius = delta;
}
registryQC.fill(HIST("maxRadiusVsPt"), jet.pt(), maxRadius);

// QA for comparison with nuclei_in_jets
TVector3 pJet(0., 0., 0.);
pJet.SetXYZ(jet.px(), jet.py(), jet.pz());
TVector3 UEAxis1(0.0, 0.0, 0.0);
TVector3 UEAxis2(0.0, 0.0, 0.0);
getPerpendicularAxis(pJet, UEAxis1, +1.0);
getPerpendicularAxis(pJet, UEAxis2, -1.0);

double NchJetPlusUE(0);
double NchJet(0);
double NchUE(0);
double ptJetPlusUE(0);
double ptJet(0);
double ptUE(0);

for (const auto& [index, track] : particles) {
TVector3 particleDir(track.px(), track.py(), track.pz());
double deltaEtaJet = particleDir.Eta() - pJet.Eta();
Expand All @@ -1375,7 +1375,7 @@
double deltaEtaUE2 = particleDir.Eta() - UEAxis2.Eta();
double deltaPhiUE2 = getDeltaPhi(particleDir.Phi(), UEAxis2.Phi());
double deltaRUE2 = std::abs(deltaEtaUE2 * deltaEtaUE2 + deltaPhiUE2 * deltaPhiUE2);

if (deltaRJet < Rmax) {
if (deltaPhiJet != -999)
registryQC.fill(HIST("deltaEtadeltaPhiJet"), deltaEtaJet, deltaPhiJet);
Expand All @@ -1395,7 +1395,7 @@
ptUE = ptUE + track.pt();
}
} // for (const auto& [index, track] : particles)

NchJet = NchJetPlusUE - 0.5 * NchUE;
ptJet = ptJetPlusUE - 0.5 * ptUE;
registryQC.fill(HIST("multiplicityJetPlusUE"), NchJetPlusUE);
Expand All @@ -1405,27 +1405,27 @@
registryQC.fill(HIST("ptJet"), ptJet);
registryQC.fill(HIST("ptUE"), 0.5 * ptUE);
registryQC.fill(HIST("deltaJetPt"), jet.pt() - ptJetPlusUE);

int nPartClusteredJet = static_cast<int>(constituents.size());

// Fill QA Histograms
if (ptJetPlusUE < minJetPt) { // swap for sub pt?

registryQC.fill(HIST("nParticlesClusteredInJet"), nPartClusteredJet);

for (const auto& track : constituents) {
registryQC.fill(HIST("ptParticlesClusteredInJet"), track.pt());
}
}

for (int i = 0; i < static_cast<int>(constituents.size()); i++) { // analyse jet constituents - this is where the magic happens
registryData.fill(HIST("trackProtocol"), 3);
fastjet::PseudoJet pseudoParticle = constituents.at(i);
int id = pseudoParticle.user_index();
const auto& jetParticle = particles.at(id);
if (!selectTrack(jetParticle))
continue;

registryData.fill(HIST("dcaXYFullJet"), jetParticle.pt() * jetParticle.sign(), jetParticle.dcaXY());
registryData.fill(HIST("dcaZFullJet"), jetParticle.pt() * jetParticle.sign(), jetParticle.dcaZ());
registryData.fill(HIST("tpcSignal"), jetParticle.pt() * jetParticle.sign(), jetParticle.tpcSignal());
Expand All @@ -1436,7 +1436,7 @@
double ptDiff = pseudoParticle.pt() - jetParticle.pt();
registryQC.fill(HIST("ptDiff"), ptDiff);
}

// if (jetParticle.pt() < minJetParticlePt)
// continue;
if (!jetParticle.has_mcParticle())
Expand Down Expand Up @@ -1495,7 +1495,7 @@
registryData.fill(HIST("ptJetAntinuclei"), jetParticle.pt());
registryQC.fill(HIST("ptJetAntinucleiVsTotalJet"), jetParticle.pt(), subtractedJetPerp.pt());
registryData.fill(HIST("trackProtocol"), 10); // # antinuclei
}
}
break;
case 1000020030:
if (!deuteronAnalysis) {
Expand Down Expand Up @@ -1560,8 +1560,9 @@

void processRun3(soa::Filtered<soa::Join<aod::JetCollisions, aod::JCollisionPIs, aod::BkgChargedRhos>>::iterator const& collision,
soa::Filtered<soa::Join<aod::ChargedJets, aod::ChargedJetConstituents>> const& allJets,
soa::Filtered<FullTracksRun3> const&/* ,
soa::Filtered<JetTracksRun3> const& */)
soa::Filtered<FullTracksRun3> const& /* ,
soa::Filtered<JetTracksRun3> const& */
)
{
registryData.fill(HIST("eventProtocol"), 0);
if (!jetderiveddatautilities::selectCollision(collision, eventSelection))
Expand Down Expand Up @@ -1604,7 +1605,7 @@
registryQC.fill(HIST("jetPtVsNumPart"), jet.pt(), jet.tracksIds().size());
registryQC.fill(HIST("maxRadiusVsPt"), jet.pt(), jet.r());

for (const auto& track/* jtrack */ : jet.template tracks_as<FullTracksRun3/* JetTracksRun3 */>()) {
for (const auto& track /* jtrack */ : jet.template tracks_as<FullTracksRun3 /* JetTracksRun3 */>()) {
// const auto& track = jtrack.track_as<FullTracksRun3>();
if (!selectTrack(track))
continue;
Expand Down Expand Up @@ -1652,11 +1653,11 @@
registryData.fill(HIST("ptJetProton"), track.pt());
registryQC.fill(HIST("ptJetProtonVsTotalJet"), track.pt(), jet.pt());
registryData.fill(HIST("trackProtocol"), 4); // # protons
} else if (isAntiproton(track, false)) { // collect antiprotons in jet
} else if (isAntiproton(track, false)) { // collect antiprotons in jet
registryData.fill(HIST("ptJetAntiproton"), track.pt());
registryQC.fill(HIST("ptJetAntiprotonVsTotalJet"), track.pt(), jet.pt());
registryData.fill(HIST("trackProtocol"), 6); // # antiprotons
} else if (isNucleus(track)) { // collect nuclei in jet
} else if (isNucleus(track)) { // collect nuclei in jet
registryData.fill(HIST("ptJetNuclei"), track.pt());
registryQC.fill(HIST("ptJetNucleiVsTotalJet"), track.pt(), jet.pt());
registryData.fill(HIST("trackProtocol"), 8); // # nuclei
Expand Down Expand Up @@ -1837,7 +1838,7 @@
}
PROCESS_SWITCH(AngularCorrelationsInJets, processRun3MCReco, "process Run 3 MC, not currently usable", false);

void processMCRun2old(McCollisions const& collisions, soa::Filtered<McTracksRun2> const& tracks, BCsWithRun2Info const&, aod::McParticles&, aod::McCollisions const&)

Check warning on line 1841 in PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-process]

Argument aod::McParticles& is not const&.
{
for (const auto& collision : collisions) {
auto bc = collision.bc_as<BCsWithRun2Info>();
Expand All @@ -1856,7 +1857,7 @@
}
PROCESS_SWITCH(AngularCorrelationsInJets, processMCRun2old, "process Run 2 MC w/o jet tables, not currently usable", false);

void processMCRun3old(McCollisions const& collisions, soa::Filtered<McTracksRun3old> const& tracks, aod::McParticles&, aod::McCollisions const&)

Check warning on line 1860 in PWGLF/Tasks/Nuspex/AngularCorrelationsInJets.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-process]

Argument aod::McParticles& is not const&.
{
for (const auto& collision : collisions) {
registryData.fill(HIST("eventProtocol"), 0);
Expand Down
Loading