Skip to content

PWGCF: Added a Monte-Carlo process function for generated and reconstructed data #7940

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
victor-gonzalez merged 5 commits into from
Oct 10, 2024
Merged

PWGCF: Added a Monte-Carlo process function for generated and reconstructed data #7940

Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
b13c6c0
Added a configurable for enabling the FakeV0Filter function, as well …
JStaa Oct 9, 2024
a124700
Added isPhysicalPrimary() check
JStaa Oct 10, 2024
ff12754
Removed Counter variable and changed a comment
JStaa Oct 10, 2024
2a05ddc
Please consider the following formatting changes
alibuild Oct 10, 2024
b9af4a8
Merge pull request #11 from alibuild/alibot-cleanup-7940
JStaa Oct 10, 2024
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
176 changes: 134 additions & 42 deletions PWGCF/MultiparticleCorrelations/Tasks/ThreeParticleCorrelations.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,11 +26,13 @@ struct ThreePartCorr {
// Histogram registry
HistogramRegistry MECorrRegistry{"MECorrRegistry", {}, OutputObjHandlingPolicy::AnalysisObject, false, true};
HistogramRegistry SECorrRegistry{"SECorrRegistry", {}, OutputObjHandlingPolicy::AnalysisObject, false, true};
HistogramRegistry MCRegistry{"MCRegistry", {}, OutputObjHandlingPolicy::AnalysisObject, false, true};
HistogramRegistry QARegistry{"QARegistry", {}, OutputObjHandlingPolicy::AnalysisObject, false, true};

// Collision filters
Filter CollCent = aod::cent::centFT0C > 0.0f && aod::cent::centFT0C < 90.0f;
Filter CollZvtx = nabs(aod::collision::posZ) < 7.0f;
Filter MCCollZvtx = nabs(aod::mccollision::posZ) < 7.0f;

// V0 filters
Filter V0Pt = aod::v0data::pt > 0.6f && aod::v0data::pt < 12.0f;
Expand All @@ -40,14 +42,24 @@ struct ThreePartCorr {
Filter TrackPt = aod::track::pt > 0.2f && aod::track::pt < 3.0f;
Filter TrackEta = nabs(aod::track::eta) < 0.8f;

// Table aliases
// Particle filters
Filter ParticlePt = aod::mcparticle::pt > 0.2f && aod::mcparticle::pt < 3.0f;
Filter ParticleEta = nabs(aod::mcparticle::eta) < 0.8f;

// Table aliases - Data
using MyFilteredCollisions = soa::Filtered<soa::Join<aod::Collisions, aod::CentFT0Cs>>;
using MyFilteredCollision = MyFilteredCollisions::iterator;
using MyFilteredV0s = soa::Filtered<aod::V0Datas>;
using MyFilteredTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra,
aod::pidTPCPi, aod::pidTPCKa, aod::pidTPCPr,
aod::pidTOFPi, aod::pidTOFKa, aod::pidTOFPr, aod::pidTOFbeta>>;

// Table aliases - MC
using MyFilteredMCGenCollision = soa::Filtered<aod::McCollisions>::iterator;
using MyFilteredMCParticles = soa::Filtered<aod::McParticles>;
using MyFilteredMCRecCollision = soa::Filtered<soa::Join<aod::Collisions, aod::McCollisionLabels>>::iterator;
using MyFilteredMCTracks = soa::Filtered<soa::Join<aod::Tracks, aod::McTrackLabels>>;

// Mixed-events binning policy
SliceCache cache;
ConfigurableAxis ConfCentBins{"ConfCentBins", {VARIABLE_WIDTH, 0.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f}, "ME Centrality binning"};
Expand All @@ -57,6 +69,9 @@ struct ThreePartCorr {
BinningType CollBinning{{ConfCentBins, ConfZvtxBins}, true};
Pair<MyFilteredCollisions, MyFilteredV0s, MyFilteredTracks, BinningType> pair{CollBinning, 5, -1, &cache};

// Process configurables
Configurable<bool> FilterSwitch{"FilterSwitch", false, "Switch for the FakeV0Filter function"};

// Particle masses
Double_t massLambda = 1.115683;
Double_t DGaussSigma = 0.0021;
Expand All @@ -77,7 +92,8 @@ struct ThreePartCorr {
const AxisSpec ZvtxAxis{ConfZvtxBins};
const AxisSpec PhiAxis{36, (-1. / 2) * M_PI, (3. / 2) * M_PI};
const AxisSpec EtaAxis{32, -1.52, 1.52};
const AxisSpec PtAxis{120, 0, 12};
const AxisSpec V0PtAxis{114, 0.6, 12};
const AxisSpec TrackPtAxis{28, 0.2, 3};
const AxisSpec LambdaInvMassAxis{100, 1.08, 1.16};

QARegistry.add("hTrackPt", "hTrackPt", {HistType::kTH1D, {{100, 0, 4}}});
Expand All @@ -98,8 +114,21 @@ struct ThreePartCorr {
QARegistry.add("hNSigmaKaon", "hNSigmaKaon", {HistType::kTH2D, {{201, -5.025, 5.025}, {201, -5.025, 5.025}}});
QARegistry.add("hNSigmaProton", "hNSigmaProton", {HistType::kTH2D, {{201, -5.025, 5.025}, {201, -5.025, 5.025}}});

QARegistry.add("hInvMassLambda", "hInvMassLambda", {HistType::kTH3D, {{LambdaInvMassAxis}, {PtAxis}, {CentralityAxis}}});
QARegistry.add("hInvMassAntiLambda", "hInvMassAntiLambda", {HistType::kTH3D, {{LambdaInvMassAxis}, {PtAxis}, {CentralityAxis}}});
QARegistry.add("hInvMassLambda", "hInvMassLambda", {HistType::kTH3D, {{LambdaInvMassAxis}, {V0PtAxis}, {CentralityAxis}}});
QARegistry.add("hInvMassAntiLambda", "hInvMassAntiLambda", {HistType::kTH3D, {{LambdaInvMassAxis}, {V0PtAxis}, {CentralityAxis}}});

MCRegistry.add("hGenPionP", "hGenMomPionP", {HistType::kTH1D, {TrackPtAxis}});
MCRegistry.add("hGenPionN", "hGenMomPionN", {HistType::kTH1D, {TrackPtAxis}});
MCRegistry.add("hGenKaonP", "hGenMomKaonP", {HistType::kTH1D, {TrackPtAxis}});
MCRegistry.add("hGenKaonN", "hGenMomKaonN", {HistType::kTH1D, {TrackPtAxis}});
MCRegistry.add("hGenProtonP", "hGenMomProtonP", {HistType::kTH1D, {TrackPtAxis}});
MCRegistry.add("hGenProtonN", "hGenMomProtonN", {HistType::kTH1D, {TrackPtAxis}});
MCRegistry.add("hRecPionP", "hRecMomPionP", {HistType::kTH1D, {TrackPtAxis}});
MCRegistry.add("hRecPionN", "hRecMomPionN", {HistType::kTH1D, {TrackPtAxis}});
MCRegistry.add("hRecKaonP", "hRecMomKaonP", {HistType::kTH1D, {TrackPtAxis}});
MCRegistry.add("hRecKaonN", "hRecMomKaonN", {HistType::kTH1D, {TrackPtAxis}});
MCRegistry.add("hRecProtonP", "hRecMomProtonP", {HistType::kTH1D, {TrackPtAxis}});
MCRegistry.add("hRecProtonN", "hRecMomProtonN", {HistType::kTH1D, {TrackPtAxis}});

SECorrRegistry.add("hSameLambdaPion_SGNL", "Same-event #Lambda - #pi correlator (SGNL region)", {HistType::kTHnSparseD, {{PhiAxis}, {EtaAxis}, {CentralityAxis}, {ZvtxAxis}, {2, -2, 2}, {2, -2, 2}}});
SECorrRegistry.add("hSameLambdaPion_SB", "Same-event #Lambda - #pi correlator (SB region)", {HistType::kTHnSparseD, {{PhiAxis}, {EtaAxis}, {CentralityAxis}, {ZvtxAxis}, {2, -2, 2}, {2, -2, 2}}});
Expand Down Expand Up @@ -172,19 +201,19 @@ struct ThreePartCorr {
DeltaEta = trigger.eta() - associate.eta();

if (CandMass >= massLambda - 4 * DGaussSigma && CandMass <= massLambda + 4 * DGaussSigma) {
if (A_PID[0] == 0) { // Pions
if (A_PID[0] == 0.0) { // Pions
SECorrRegistry.fill(HIST("hSameLambdaPion_SGNL"), DeltaPhi, DeltaEta, collision.centFT0C(), collision.posZ(), T_Sign, associate.sign());
} else if (A_PID[0] == 1) { // Kaons
} else if (A_PID[0] == 1.0) { // Kaons
SECorrRegistry.fill(HIST("hSameLambdaKaon_SGNL"), DeltaPhi, DeltaEta, collision.centFT0C(), collision.posZ(), T_Sign, associate.sign());
} else if (A_PID[0] == 2) { // Protons
} else if (A_PID[0] == 2.0) { // Protons
SECorrRegistry.fill(HIST("hSameLambdaProton_SGNL"), DeltaPhi, DeltaEta, collision.centFT0C(), collision.posZ(), T_Sign, associate.sign());
}
} else if (CandMass >= massLambda - 8 * DGaussSigma && CandMass <= massLambda + 8 * DGaussSigma) {
if (A_PID[0] == 0) { // Pions
if (A_PID[0] == 0.0) { // Pions
SECorrRegistry.fill(HIST("hSameLambdaPion_SB"), DeltaPhi, DeltaEta, collision.centFT0C(), collision.posZ(), T_Sign, associate.sign());
} else if (A_PID[0] == 1) { // Kaons
} else if (A_PID[0] == 1.0) { // Kaons
SECorrRegistry.fill(HIST("hSameLambdaKaon_SB"), DeltaPhi, DeltaEta, collision.centFT0C(), collision.posZ(), T_Sign, associate.sign());
} else if (A_PID[0] == 2) { // Protons
} else if (A_PID[0] == 2.0) { // Protons
SECorrRegistry.fill(HIST("hSameLambdaProton_SB"), DeltaPhi, DeltaEta, collision.centFT0C(), collision.posZ(), T_Sign, associate.sign());
}
}
Expand All @@ -195,13 +224,10 @@ struct ThreePartCorr {
}
// End of the V0-Track Correlations
}
PROCESS_SWITCH(ThreePartCorr, processSame, "Process same-event correlations", true);

void processMixed(MyFilteredCollisions const& collisions, MyFilteredV0s const& v0s, MyFilteredTracks const& tracks)
void processMixed(MyFilteredCollisions const&, MyFilteredV0s const&, MyFilteredTracks const&)
{

LOGF(info, "Input data Collisions %d, V0s %d, Tracks %d ", collisions.size(), v0s.size(), tracks.size());

// Start of the Mixed-events Correlations
for (const auto& [coll_1, v0_1, coll_2, track_2] : pair) {
for (const auto& [trigger, associate] : soa::combinations(soa::CombinationsFullIndexPolicy(v0_1, track_2))) {
Expand All @@ -220,19 +246,19 @@ struct ThreePartCorr {
DeltaEta = trigger.eta() - associate.eta();

if (CandMass >= massLambda - 4 * DGaussSigma && CandMass <= massLambda + 4 * DGaussSigma) {
if (A_PID[0] == 0) { // Pions
if (A_PID[0] == 0.0) { // Pions
MECorrRegistry.fill(HIST("hMixLambdaPion_SGNL"), DeltaPhi, DeltaEta, coll_1.centFT0C(), coll_1.posZ(), T_Sign, associate.sign());
} else if (A_PID[0] == 1) { // Kaons
} else if (A_PID[0] == 1.0) { // Kaons
MECorrRegistry.fill(HIST("hMixLambdaKaon_SGNL"), DeltaPhi, DeltaEta, coll_1.centFT0C(), coll_1.posZ(), T_Sign, associate.sign());
} else if (A_PID[0] == 2) { // Protons
} else if (A_PID[0] == 2.0) { // Protons
MECorrRegistry.fill(HIST("hMixLambdaProton_SGNL"), DeltaPhi, DeltaEta, coll_1.centFT0C(), coll_1.posZ(), T_Sign, associate.sign());
}
} else if (CandMass >= massLambda - 8 * DGaussSigma && CandMass <= massLambda + 8 * DGaussSigma) {
if (A_PID[0] == 0) { // Pions
if (A_PID[0] == 0.0) { // Pions
MECorrRegistry.fill(HIST("hMixLambdaPion_SB"), DeltaPhi, DeltaEta, coll_1.centFT0C(), coll_1.posZ(), T_Sign, associate.sign());
} else if (A_PID[0] == 1) { // Kaons
} else if (A_PID[0] == 1.0) { // Kaons
MECorrRegistry.fill(HIST("hMixLambdaKaon_SB"), DeltaPhi, DeltaEta, coll_1.centFT0C(), coll_1.posZ(), T_Sign, associate.sign());
} else if (A_PID[0] == 2) { // Protons
} else if (A_PID[0] == 2.0) { // Protons
MECorrRegistry.fill(HIST("hMixLambdaProton_SB"), DeltaPhi, DeltaEta, coll_1.centFT0C(), coll_1.posZ(), T_Sign, associate.sign());
}
}
Expand All @@ -242,7 +268,70 @@ struct ThreePartCorr {
}
// End of the Mixed-events Correlations
}

void processMCGen(MyFilteredMCGenCollision const&, MyFilteredMCParticles const& particles)
{

// Start of the Monte-Carlo generated QA
for (const auto& particle : particles) {
if (particle.isPhysicalPrimary()) {

if (particle.pdgCode() == 211) { // Pos pions
MCRegistry.fill(HIST("hGenPionP"), particle.pt());
} else if (particle.pdgCode() == -211) { // Neg pions
MCRegistry.fill(HIST("hGenPionN"), particle.pt());
} else if (particle.pdgCode() == 310) { // Pos kaons
MCRegistry.fill(HIST("hGenKaonP"), particle.pt());
} else if (particle.pdgCode() == -310) { // Neg kaons
MCRegistry.fill(HIST("hGenKaonN"), particle.pt());
} else if (particle.pdgCode() == 2212) { // Pos protons
MCRegistry.fill(HIST("hGenProtonP"), particle.pt());
} else if (particle.pdgCode() == -2212) { // Neg protons
MCRegistry.fill(HIST("hGenProtonN"), particle.pt());
}
}
}
// End of the Monte-Carlo generated QA
}

void processMCRec(MyFilteredMCRecCollision const& collision, MyFilteredMCTracks const& tracks, aod::McCollisions const&, aod::McParticles const&)
{

if (!collision.has_mcCollision()) {
return;
}

// Start of the Monte-Carlo reconstructed QA
for (const auto& track : tracks) {
if (!track.has_mcParticle()) {
continue;
}

auto particle = track.mcParticle();
if (particle.isPhysicalPrimary()) {

if (particle.pdgCode() == 211) { // Pos pions
MCRegistry.fill(HIST("hRecPionP"), particle.pt());
} else if (particle.pdgCode() == -211) { // Neg pions
MCRegistry.fill(HIST("hRecPionN"), particle.pt());
} else if (particle.pdgCode() == 310) { // Pos kaons
MCRegistry.fill(HIST("hRecKaonP"), particle.pt());
} else if (particle.pdgCode() == -310) { // Neg kaons
MCRegistry.fill(HIST("hRecKaonN"), particle.pt());
} else if (particle.pdgCode() == 2212) { // Pos protons
MCRegistry.fill(HIST("hRecProtonP"), particle.pt());
} else if (particle.pdgCode() == -2212) { // Neg protons
MCRegistry.fill(HIST("hRecProtonN"), particle.pt());
}
Comment on lines +310 to +325
Copy link
Collaborator

@victor-gonzalez victor-gonzalez Oct 10, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would suggest the same processing as for data
For not replicating code, a single function with the common functionality called from both process when appropriate will do the magic

}
}
// End of the Monte-Carlo reconstructed QA
}

PROCESS_SWITCH(ThreePartCorr, processSame, "Process same-event correlations", true);
PROCESS_SWITCH(ThreePartCorr, processMixed, "Process mixed-event correlations", true);
PROCESS_SWITCH(ThreePartCorr, processMCGen, "Process Monte-Carlo, generator level", false);
PROCESS_SWITCH(ThreePartCorr, processMCRec, "Process Monte-Carlo, reconstructed level", false);

//================================================================================================================================================================================================================

Expand Down Expand Up @@ -345,29 +434,32 @@ struct ThreePartCorr {
Bool_t FakeV0Filter(const V0Cand& V0, const TrackCand& Track)
{

TLorentzVector Daughter, Associate;
if (TrackPID(Track)[0] == 1.0) { // Kaons
return kTRUE;
} else if (V0Sign(V0) == 1 && TrackPID(Track)[0] == 0 && Track.sign() == -1) { // Lambda - Pi_min
const auto& dTrack = V0.template posTrack_as<MyFilteredTracks>();
Daughter.SetPtEtaPhiM(dTrack.pt(), dTrack.eta(), dTrack.phi(), o2::constants::physics::MassProton);
Associate.SetPtEtaPhiM(Track.pt(), Track.eta(), Track.phi(), o2::constants::physics::MassPionCharged);
} else if (V0Sign(V0) == -1 && TrackPID(Track)[0] == 0 && Track.sign() == 1) { // Antilambda - Pi_plus
const auto& dTrack = V0.template negTrack_as<MyFilteredTracks>();
Daughter.SetPtEtaPhiM(dTrack.pt(), dTrack.eta(), dTrack.phi(), o2::constants::physics::MassProton);
Associate.SetPtEtaPhiM(Track.pt(), Track.eta(), Track.phi(), o2::constants::physics::MassPionCharged);
} else if (V0Sign(V0) == 1 && TrackPID(Track)[0] == 2 && Track.sign() == 1) { // Lambda - Proton
const auto& dTrack = V0.template negTrack_as<MyFilteredTracks>();
Daughter.SetPtEtaPhiM(dTrack.pt(), dTrack.eta(), dTrack.phi(), o2::constants::physics::MassPionCharged);
Associate.SetPtEtaPhiM(Track.pt(), Track.eta(), Track.phi(), o2::constants::physics::MassProton);
} else if (V0Sign(V0) == -1 && TrackPID(Track)[0] == 2 && Track.sign() == -1) { // Antilambda - Antiproton
const auto& dTrack = V0.template posTrack_as<MyFilteredTracks>();
Daughter.SetPtEtaPhiM(dTrack.pt(), dTrack.eta(), dTrack.phi(), o2::constants::physics::MassPionCharged);
Associate.SetPtEtaPhiM(Track.pt(), Track.eta(), Track.phi(), o2::constants::physics::MassProton);
}
if (FilterSwitch) {

TLorentzVector Daughter, Associate;
if (TrackPID(Track)[0] == 1.0) { // Kaons
return kTRUE;
} else if (V0Sign(V0) == 1 && TrackPID(Track)[0] == 0.0 && Track.sign() == -1) { // Lambda - Pi_min
const auto& dTrack = V0.template posTrack_as<MyFilteredTracks>();
Daughter.SetPtEtaPhiM(dTrack.pt(), dTrack.eta(), dTrack.phi(), o2::constants::physics::MassProton);
Associate.SetPtEtaPhiM(Track.pt(), Track.eta(), Track.phi(), o2::constants::physics::MassPionCharged);
} else if (V0Sign(V0) == -1 && TrackPID(Track)[0] == 0.0 && Track.sign() == 1) { // Antilambda - Pi_plus
const auto& dTrack = V0.template negTrack_as<MyFilteredTracks>();
Daughter.SetPtEtaPhiM(dTrack.pt(), dTrack.eta(), dTrack.phi(), o2::constants::physics::MassProton);
Associate.SetPtEtaPhiM(Track.pt(), Track.eta(), Track.phi(), o2::constants::physics::MassPionCharged);
} else if (V0Sign(V0) == 1 && TrackPID(Track)[0] == 2.0 && Track.sign() == 1) { // Lambda - Proton
const auto& dTrack = V0.template negTrack_as<MyFilteredTracks>();
Daughter.SetPtEtaPhiM(dTrack.pt(), dTrack.eta(), dTrack.phi(), o2::constants::physics::MassPionCharged);
Associate.SetPtEtaPhiM(Track.pt(), Track.eta(), Track.phi(), o2::constants::physics::MassProton);
} else if (V0Sign(V0) == -1 && TrackPID(Track)[0] == 2.0 && Track.sign() == -1) { // Antilambda - Antiproton
const auto& dTrack = V0.template posTrack_as<MyFilteredTracks>();
Daughter.SetPtEtaPhiM(dTrack.pt(), dTrack.eta(), dTrack.phi(), o2::constants::physics::MassPionCharged);
Associate.SetPtEtaPhiM(Track.pt(), Track.eta(), Track.phi(), o2::constants::physics::MassProton);
}

if ((Daughter + Associate).M() >= massLambda - 4 * DGaussSigma && (Daughter + Associate).M() <= massLambda + 4 * DGaussSigma) {
return kFALSE;
if ((Daughter + Associate).M() >= massLambda - 4 * DGaussSigma && (Daughter + Associate).M() <= massLambda + 4 * DGaussSigma) {
return kFALSE;
}
Comment on lines +439 to +462
Copy link
Collaborator

@victor-gonzalez victor-gonzalez Oct 10, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Consider using the functionality on https://github.com/AliceO2Group/O2Physics/blob/master/Common/Core/RecoDecay.h
for extracting the invariant mass instead of using TLorentzVector. It will be much more efficient

}

return kTRUE;
Expand Down