Skip to content

[PWGHF] add new classifications for like-sign muon pairs #15683

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
NicoleBastid merged 5 commits into from
Apr 10, 2026
+61 −42
Merged

[PWGHF] add new classifications for like-sign muon pairs #15683

Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
bddb6d0
Update taskSingleMuonSource.cxx
MaolinZH Apr 8, 2026
97a764a
Fix the O2linter errors
MaolinZH Apr 9, 2026
857c66c
Remove magic numbers to fix the O2Linter
MaolinZH Apr 9, 2026
aa94fd9
Fix the headerfile format and last magic number
MaolinZH Apr 9, 2026
69edf1e
Fix the O2linter warnings
MaolinZH Apr 9, 2026
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
103 changes: 61 additions & 42 deletions PWGHF/HFL/Tasks/taskSingleMuonSource.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,10 +8,10 @@
// 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.
//
// \file taskSingleMuonSource.cxx
// \brief Task used to seperate single muons source in Monte Carlo simulation.
// \author Maolin Zhang <maolin.zhang@cern.ch>, CCNU
///
/// \file taskSingleMuonSource.cxx
/// \brief Task used to seperate single muons source in Monte Carlo simulation.
/// \author Maolin Zhang <maolin.zhang@cern.ch>, CCNU

#include "Common/Core/RecoDecay.h"
#include "Common/DataModel/TrackSelectionTables.h"
Expand All @@ -24,11 +24,11 @@
#include <Framework/HistogramRegistry.h>
#include <Framework/HistogramSpec.h>
#include <Framework/InitContext.h>
#include <Framework/O2DatabasePDGPlugin.h>
#include <Framework/OutputObjHeader.h>
#include <Framework/runDataProcessing.h>

#include <Math/Vector4D.h>
#include <TDatabasePDG.h>
#include <TPDGCode.h>
#include <TString.h>

Expand Down Expand Up @@ -78,6 +78,8 @@ struct HfTaskSingleMuonSource {
Configurable<int> charge{"charge", 0, "Muon track charge, validated values are 0, 1 and -1, 0 represents both 1 and -1"};
Configurable<bool> pairSource{"pairSource", true, "check also the source of like-sign muon pairs"};

Service<o2::framework::O2DatabasePDG> pdgDB;

double pDcaMax = 594.0; // p*DCA maximum value for small Rab
double pDcaMax2 = 324.0; // p*DCA maximum value for large Rabs
double rAbsMid = 26.5; // R at absorber end minimum value
Expand Down Expand Up @@ -116,7 +118,7 @@ struct HfTaskSingleMuonSource {

HistogramConfigSpec const h1ColNumber{HistType::kTH1F, {axisColNumber}};
HistogramConfigSpec const h1Pt{HistType::kTH1F, {axisPt}};
HistogramConfigSpec h1Mass{HistType::kTH1F, {axisMass}};
HistogramConfigSpec const h1Mass{HistType::kTH1F, {axisMass}};
HistogramConfigSpec const h2PtDCA{HistType::kTH2F, {axisPt, axisDCA}};
HistogramConfigSpec const h2PtChi2{HistType::kTH2F, {axisPt, axisChi2}};
HistogramConfigSpec const h2PtDeltaPt{HistType::kTH2F, {axisPt, axisDeltaPt}};
Expand All @@ -125,9 +127,11 @@ struct HfTaskSingleMuonSource {
registry.add("h1MuBeforeCuts", "", h1Pt);
registry.add("h1MuonMass", "", h1Mass);
registry.add("h1BeautyMass", "", h1Mass);
registry.add("h1CorrBeautyMass", "", h1Mass);
registry.add("h1OtherMass", "", h1Mass);
registry.add("h1MuonMassGen", "", h1Mass);
registry.add("h1BeautyMassGen", "", h1Mass);
registry.add("h1CorrBeautyMassGen", "", h1Mass);
registry.add("h1OtherMassGen", "", h1Mass);
for (const auto& src : muonSources) {
registry.add(Form("h1%sPt", src.Data()), "", h1Pt);
Expand All @@ -141,6 +145,8 @@ struct HfTaskSingleMuonSource {
uint8_t getMask(const McMuons::iterator& muon)
{
uint8_t mask(0);
const int diquarkEdge = 1000;
const int hadronEdge = 10000;
if (muon.has_mcParticle()) {
SETBIT(mask, IsIdentified);
} else {
Expand All @@ -159,7 +165,7 @@ struct HfTaskSingleMuonSource {
mcPart = *(mcPart.mothers_first_as<aod::McParticles>());

const auto pdgAbs(std::abs(mcPart.pdgCode()));
if (pdgAbs < 10 || pdgAbs == 21) {
if (pdgAbs < kElectron || pdgAbs == kGluon) {
break; // Quark and gluon
}

Expand All @@ -168,40 +174,39 @@ struct HfTaskSingleMuonSource {
continue;
}

if (pdgAbs == kTauMinus) {
// Tau
if (pdgAbs == kTauMinus) { // Tau
SETBIT(mask, HasTauParent);
continue;
}

const int pdgRem(pdgAbs % 100000);
const int pdgRemRem(pdgRem % 100);

if (pdgRem == kProton) {
continue;
} // Beam particle

if ((pdgRem < 100) || (pdgRem >= 10000)) {
if ((pdgRem < kPi0) || (pdgRem >= hadronEdge)) {
continue;
}
if ((pdgRem % 100 == 1 || pdgRem % 100 == 3) && pdgRem > 1000) { // diquarks
if ((pdgRemRem == kDown || pdgRemRem == kStrange) && pdgRem > diquarkEdge) { // diquarks
continue;
}
// compute the flavor of constituent quark
const int flv(pdgRem / std::pow(10, static_cast<int>(std::log10(pdgRem))));
if (flv > 6) {
if (flv > kTop) {
// no more than 6 flavors
continue;
}
if (flv < 4) {
if (flv < kCharm) {
// light flavor
SETBIT(mask, HasLightParent);
continue;
}

auto* pdgData(TDatabasePDG::Instance()->GetParticle(mcPart.pdgCode()));
auto* pdgData = pdgDB->GetParticle(mcPart.pdgCode());
if ((pdgData != nullptr) && (pdgData->AntiParticle() == nullptr)) {
SETBIT(mask, HasQuarkoniumParent);
} else if (flv == 4) {
} else if (flv == kCharm) {
SETBIT(mask, HasCharmParent);
} else {
SETBIT(mask, HasBeautyParent);
Expand Down Expand Up @@ -274,11 +279,13 @@ struct HfTaskSingleMuonSource {
// fill the histograms of each particle types
void fillHistograms(const McMuons::iterator& muon)
{
const int type0 = 0;
const int type2 = 2;
const auto mask(getMask(muon));
const auto pt(muon.pt()), chi2(muon.chi2MatchMCHMFT());
const auto dca(RecoDecay::sqrtSumOfSquares(muon.fwdDcaX(), muon.fwdDcaY()));

if (trackType == 0 || trackType == 2) {
if (trackType == type0 || trackType == type2) {
if (!muon.has_matchMCHTrack()) {
return;
}
Expand Down Expand Up @@ -344,6 +351,8 @@ struct HfTaskSingleMuonSource {
int traceAncestor(const McMuons::iterator& muon, aod::McParticles const& mctracks)
{
int mcNum = 0;
const int hadronStatus = 80;
const int diquarkEdge = 1000;
if (!muon.has_mcParticle()) {
return 0;
}
Expand All @@ -353,36 +362,38 @@ struct HfTaskSingleMuonSource {
}
while (mcPart.has_mothers()) { // the first hadron after hadronization
auto mother = mcPart.mothers_first_as<aod::McParticles>();
if (std::abs(mother.getGenStatusCode()) < 80) {
if (std::abs(mother.getGenStatusCode()) < hadronStatus) {
break;
}
mcPart = mother;
}
int flv = mcPart.pdgCode() / std::pow(10, static_cast<int>(std::log10(std::abs(mcPart.pdgCode()))));
if (abs(flv) == 5 && mcPart.pdgCode() < 1000)
if (std::abs(flv) == kBottom && mcPart.pdgCode() < diquarkEdge) {
flv = -flv;
}
for (int i = (mcPart.mothers_first_as<aod::McParticles>()).globalIndex(); i <= (mcPart.mothers_last_as<aod::McParticles>()).globalIndex(); i++) { // loop over the lund string
for (auto mctrack : mctracks) {
for (const auto& mctrack : mctracks) {
if (mctrack.globalIndex() != i) {
continue;
}
if ((mctrack.pdgCode() != flv) && (abs(mctrack.pdgCode()) < abs(flv) * 1000)) {
if ((mctrack.pdgCode() != flv) && (std::abs(mctrack.pdgCode()) < std::abs(flv) * 1000)) {
continue;
}
while (mctrack.has_mothers()) {
int motherflv = (mctrack.mothers_first_as<aod::McParticles>()).pdgCode() / std::pow(10, static_cast<int>(std::log10(abs((mctrack.mothers_first_as<aod::McParticles>()).pdgCode())))); // find the mother with same flavor
auto mother = (abs(motherflv) == abs(flv)) ? (mctrack.mothers_first_as<aod::McParticles>()) : (mctrack.mothers_last_as<aod::McParticles>());
if ((mother.pdgCode() != mctrack.pdgCode()) && (abs(mctrack.pdgCode()) < 10)) { // both mother is not the the quark with same flavor
mcNum = mctrack.globalIndex();
auto currentTrk = mctrack;
while (currentTrk.has_mothers()) {
int motherflv = (currentTrk.mothers_first_as<aod::McParticles>()).pdgCode() / std::pow(10, static_cast<int>(std::log10(std::abs((currentTrk.mothers_first_as<aod::McParticles>()).pdgCode())))); // find the mother with same flavor
auto mother = (std::abs(motherflv) == std::abs(flv)) ? (currentTrk.mothers_first_as<aod::McParticles>()) : (currentTrk.mothers_last_as<aod::McParticles>());
if ((mother.pdgCode() != currentTrk.pdgCode()) && (std::abs(currentTrk.pdgCode()) < kElectron)) { // both mother is not the the quark with same flavor
mcNum = currentTrk.globalIndex();
return mcNum;
}
mctrack = mother;
currentTrk = mother;
}
}
}
return 0;
}
bool Corr(const McMuons::iterator& muon1, const McMuons::iterator& muon2, aod::McParticles const& mcParts)
bool isCorr(const McMuons::iterator& muon1, const McMuons::iterator& muon2, aod::McParticles const& mcParts)
{

int moth11(0), moth12(0), moth21(1), moth22(1);
Expand All @@ -391,7 +402,7 @@ struct HfTaskSingleMuonSource {
if (anc1 == 0 || anc2 == 0) {
return false;
}
for (auto mcPart : mcParts) {
for (const auto& mcPart : mcParts) {
if (mcPart.globalIndex() == anc1) {
moth11 = (mcPart.mothers_first_as<aod::McParticles>()).globalIndex();
moth12 = (mcPart.mothers_last_as<aod::McParticles>()).globalIndex();
Expand All @@ -408,7 +419,8 @@ struct HfTaskSingleMuonSource {
}
void fillPairs(const McMuons::iterator& muon, const McMuons::iterator& muon2, aod::McParticles const& mcParts)
{
if (trackType != 3) {
const int type3 = 3;
if (trackType != type3) {
return;
}
float mm = o2::constants::physics::MassMuon;
Expand All @@ -419,7 +431,7 @@ struct HfTaskSingleMuonSource {
ROOT::Math::PtEtaPhiMVector mu1Vec(muon.pt(), muon.eta(), muon.phi(), mm);
ROOT::Math::PtEtaPhiMVector mu2Vec(muon2.pt(), muon2.eta(), muon2.phi(), mm);
ROOT::Math::PtEtaPhiMVector dimuVec = mu1Vec + mu2Vec;
auto InvM = dimuVec.M();
auto invMass = dimuVec.M();

if (!muon.has_mcParticle() || !muon2.has_mcParticle()) {
return;
Expand All @@ -430,18 +442,22 @@ struct HfTaskSingleMuonSource {
ROOT::Math::PtEtaPhiMVector mu1VecGen(mcPart1.pt(), mcPart1.eta(), mcPart1.phi(), mm);
ROOT::Math::PtEtaPhiMVector mu2VecGen(mcPart2.pt(), mcPart2.eta(), mcPart2.phi(), mm);
ROOT::Math::PtEtaPhiMVector dimuVecGen = mu1VecGen + mu2VecGen;
auto InvMGen = dimuVecGen.M();
auto invMassGen = dimuVecGen.M();

if (isMuon(mask1) && isMuon(mask2)) {
registry.fill(HIST("h1MuonMass"), InvM);
registry.fill(HIST("h1MuonMassGen"), InvMGen);
registry.fill(HIST("h1MuonMass"), invMass);
registry.fill(HIST("h1MuonMassGen"), invMassGen);
}
if (Corr(muon, muon2, mcParts) && isBeautyMu(mask1) && isBeautyMu(mask2)) {
registry.fill(HIST("h1BeautyMass"), InvM);
registry.fill(HIST("h1BeautyMassGen"), InvMGen);
if (isBeautyMu(mask1) && isBeautyMu(mask2)) {
registry.fill(HIST("h1BeautyMass"), invMass);
registry.fill(HIST("h1BeautyMassGen"), invMassGen);
if (isCorr(muon, muon2, mcParts)) {
registry.fill(HIST("h1CorrBeautyMass"), invMass);
registry.fill(HIST("h1CorrBeautyMassGen"), invMassGen);
}
} else {
registry.fill(HIST("h1OtherMass"), InvM);
registry.fill(HIST("h1OtherMassGen"), InvMGen);
registry.fill(HIST("h1OtherMass"), invMass);
registry.fill(HIST("h1OtherMassGen"), invMassGen);
}
}

Expand Down Expand Up @@ -477,7 +493,10 @@ struct HfTaskSingleMuonSource {
continue;
}
}
if ((muon.chi2() >= 1e6) || (muon.chi2() < 0)) {
if (muon.chi2() < 0) {
continue;
}
if (muon.chi2MatchMCHMID() < 0) {
continue;
}
if (charge != 0 && muon.sign() != charge) {
Expand Down Expand Up @@ -517,10 +536,10 @@ struct HfTaskSingleMuonSource {
}
}

if ((muon2.chi2() >= 1e6) || (muon2.chi2() < 0)) {
if (muon2.chi2() < 0) {
continue;
}
if ((muon2.chi2MatchMCHMID() >= 1e6) || (muon2.chi2MatchMCHMID() < 0)) {
if (muon2.chi2MatchMCHMID() < 0) {
continue;
}
fillPairs(muon, muon2, mcParts);
Expand Down
Loading