diff --git a/PWGHF/HFL/Tasks/taskSingleMuonSource.cxx b/PWGHF/HFL/Tasks/taskSingleMuonSource.cxx index aea4ee26bea..e2d209966dc 100644 --- a/PWGHF/HFL/Tasks/taskSingleMuonSource.cxx +++ b/PWGHF/HFL/Tasks/taskSingleMuonSource.cxx @@ -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 , CCNU +/// +/// \file taskSingleMuonSource.cxx +/// \brief Task used to seperate single muons source in Monte Carlo simulation. +/// \author Maolin Zhang , CCNU #include "Common/Core/RecoDecay.h" #include "Common/DataModel/TrackSelectionTables.h" @@ -24,11 +24,11 @@ #include #include #include +#include #include #include #include -#include #include #include @@ -78,6 +78,8 @@ struct HfTaskSingleMuonSource { Configurable charge{"charge", 0, "Muon track charge, validated values are 0, 1 and -1, 0 represents both 1 and -1"}; Configurable pairSource{"pairSource", true, "check also the source of like-sign muon pairs"}; + Service 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 @@ -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}}; @@ -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); @@ -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 { @@ -159,7 +165,7 @@ struct HfTaskSingleMuonSource { mcPart = *(mcPart.mothers_first_as()); const auto pdgAbs(std::abs(mcPart.pdgCode())); - if (pdgAbs < 10 || pdgAbs == 21) { + if (pdgAbs < kElectron || pdgAbs == kGluon) { break; // Quark and gluon } @@ -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(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); @@ -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; } @@ -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; } @@ -353,36 +362,38 @@ struct HfTaskSingleMuonSource { } while (mcPart.has_mothers()) { // the first hadron after hadronization auto mother = mcPart.mothers_first_as(); - if (std::abs(mother.getGenStatusCode()) < 80) { + if (std::abs(mother.getGenStatusCode()) < hadronStatus) { break; } mcPart = mother; } int flv = mcPart.pdgCode() / std::pow(10, static_cast(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()).globalIndex(); i <= (mcPart.mothers_last_as()).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()).pdgCode() / std::pow(10, static_cast(std::log10(abs((mctrack.mothers_first_as()).pdgCode())))); // find the mother with same flavor - auto mother = (abs(motherflv) == abs(flv)) ? (mctrack.mothers_first_as()) : (mctrack.mothers_last_as()); - 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()).pdgCode() / std::pow(10, static_cast(std::log10(std::abs((currentTrk.mothers_first_as()).pdgCode())))); // find the mother with same flavor + auto mother = (std::abs(motherflv) == std::abs(flv)) ? (currentTrk.mothers_first_as()) : (currentTrk.mothers_last_as()); + 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); @@ -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()).globalIndex(); moth12 = (mcPart.mothers_last_as()).globalIndex(); @@ -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; @@ -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; @@ -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); } } @@ -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) { @@ -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);