Skip to content

Commit 5be7543

Browse files
MaolinZHMaolinZH
authored
Update taskSingleMuonSource.cxx
1. Fix the source separation bug: consider also the diquarks and gluons 2. Add also the source check function for correlated like-sign muon pairs decay from B hadrons
1 parent 64cd5f7 commit 5be7543

File tree

1 file changed

+184
-18
lines changed

1 file changed

+184
-18
lines changed

‎PWGHF/HFL/Tasks/taskSingleMuonSource.cxx‎

Lines changed: 184 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,6 @@
1414
// \author Maolin Zhang <maolin.zhang@cern.ch>, CCNU
1515

1616
#include "Common/Core/RecoDecay.h"
17-
#include "Common/DataModel/EventSelection.h"
1817
#include "Common/DataModel/TrackSelectionTables.h"
1918

2019
#include <Framework/ASoA.h>
@@ -31,6 +30,7 @@
3130
#include <TDatabasePDG.h>
3231
#include <TPDGCode.h>
3332
#include <TString.h>
33+
#include <Math/Vector4D.h>
3434

3535
#include <Rtypes.h>
3636

@@ -41,7 +41,7 @@
4141
using namespace o2;
4242
using namespace o2::aod;
4343
using namespace o2::framework;
44-
using MyCollisions = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels>;
44+
using MyCollisions = soa::Join<aod::Collisions, aod::McCollisionLabels>;
4545
using McMuons = soa::Join<aod::FwdTracks, aod::McFwdTrackLabels, aod::FwdTracksDCA>;
4646

4747
namespace
@@ -74,15 +74,19 @@ struct HfTaskSingleMuonSource {
7474
Produces<aod::HfMuonSource> singleMuonSource;
7575

7676
Configurable<int> mcMaskSelection{"mcMaskSelection", 0, "McMask for correct match, valid values are 0 and 128"};
77-
Configurable<int> trackType{"trackType", 0, "Muon track type, validated values are 0, 1, 2, 3 and 4"};
78-
Configurable<int> charge{"charge", -1, "Muon track charge, validated values are 0, 1 and -1, 0 represents both 1 and -1"};
79-
80-
double pDcaMax = 594.0; // p*DCA maximum value for large Rabs
81-
double rAbsMin = 26.5; // R at absorber end minimum value
82-
double rAbsMax = 89.5; // R at absorber end maximum value
83-
double etaLow = -3.6; // low edge of eta acceptance
84-
double etaUp = -2.5; // up edge of eta acceptance
85-
double edgeZ = 10.0; // edge of event position Z
77+
Configurable<int> trackType{"trackType", 3, "Muon track type, validated values are 0, 1, 2, 3 and 4"};
78+
Configurable<int> charge{"charge", 0, "Muon track charge, validated values are 0, 1 and -1, 0 represents both 1 and -1"};
79+
Configurable<bool> pairSource{"pairSource", true, "check also the source of like-sign muon pairs"};
80+
81+
double pDcaMax = 594.0; // p*DCA maximum value for small Rab
82+
double pDcaMax2 = 324.0; // p*DCA maximum value for large Rabs
83+
double rAbsMid = 26.5; // R at absorber end minimum value
84+
double rAbsMax = 89.5; // R at absorber end maximum value
85+
double rAbsMin = 17.6; // R at absorber end maximum value
86+
double etaLow = -4.0; // low edge of eta acceptance
87+
double etaUp = -2.5; // up edge of eta acceptance
88+
double edgeZ = 10.0; // edge of event position Z
89+
double ptLow = 1.0; // low edge of pT for muon pairs
8690

8791
HistogramRegistry registry{
8892
"registry",
@@ -108,14 +112,23 @@ struct HfTaskSingleMuonSource {
108112
AxisSpec const axisChi2{500, 0., 100., "#chi^{2} of MCH-MFT matching"};
109113
AxisSpec const axisPt{200, 0., 100., "#it{p}_{T,reco} (GeV/#it{c})"};
110114
AxisSpec const axisDeltaPt{1000, -50., 50., "#Delta #it{p}_{T} (GeV/#it{c})"};
115+
AxisSpec const axisMass{200, 0., 20., "Inv.Mass (GeV/#it{c}^2)"};
111116

112117
HistogramConfigSpec const h1ColNumber{HistType::kTH1F, {axisColNumber}};
113118
HistogramConfigSpec const h1Pt{HistType::kTH1F, {axisPt}};
119+
HistogramConfigSpec h1Mass{HistType::kTH1F, {axisMass}};
114120
HistogramConfigSpec const h2PtDCA{HistType::kTH2F, {axisPt, axisDCA}};
115121
HistogramConfigSpec const h2PtChi2{HistType::kTH2F, {axisPt, axisChi2}};
116122
HistogramConfigSpec const h2PtDeltaPt{HistType::kTH2F, {axisPt, axisDeltaPt}};
117123

118124
registry.add("h1ColNumber", "", h1ColNumber);
125+
registry.add("h1MuBeforeCuts", "", h1Pt);
126+
registry.add("h1MuonMass", "", h1Mass);
127+
registry.add("h1BeautyMass", "", h1Mass);
128+
registry.add("h1OtherMass", "", h1Mass);
129+
registry.add("h1MuonMassGen", "", h1Mass);
130+
registry.add("h1BeautyMassGen", "", h1Mass);
131+
registry.add("h1OtherMassGen", "", h1Mass);
119132
for (const auto& src : muonSources) {
120133
registry.add(Form("h1%sPt", src.Data()), "", h1Pt);
121134
registry.add(Form("h2%sPtDCA", src.Data()), "", h2PtDCA);
@@ -146,8 +159,8 @@ struct HfTaskSingleMuonSource {
146159
mcPart = *(mcPart.mothers_first_as<aod::McParticles>());
147160

148161
const auto pdgAbs(std::abs(mcPart.pdgCode()));
149-
if (pdgAbs < 10) {
150-
break; // Quark
162+
if (pdgAbs < 10 || pdgAbs == 21) {
163+
break; // Quark and gluon
151164
}
152165

153166
if (!mcPart.producedByGenerator()) { // Produced in transport code
@@ -170,6 +183,9 @@ struct HfTaskSingleMuonSource {
170183
if ((pdgRem < 100) || (pdgRem >= 10000)) {
171184
continue;
172185
}
186+
if ((pdgRem % 100 == 1 || pdgRem % 100 == 3) && pdgRem > 1000) { // diquarks
187+
continue;
188+
}
173189
// compute the flavor of constituent quark
174190
const int flv(pdgRem / std::pow(10, static_cast<int>(std::log10(pdgRem))));
175191
if (flv > 6) {
@@ -325,9 +341,113 @@ struct HfTaskSingleMuonSource {
325341
}
326342
}
327343

344+
int traceAncestor(const McMuons::iterator& muon, aod::McParticles const& mctracks)
345+
{
346+
int mcNum = 0;
347+
if (!muon.has_mcParticle()) {
348+
return 0;
349+
}
350+
auto mcPart(muon.mcParticle());
351+
if (std::abs(mcPart.pdgCode()) != kMuonMinus) {
352+
return 0;
353+
}
354+
while (mcPart.has_mothers()) { // the first hadron after hadronization
355+
auto mother = mcPart.mothers_first_as<aod::McParticles>();
356+
if (std::abs(mother.getGenStatusCode()) < 80) {
357+
break;
358+
}
359+
mcPart = mother;
360+
}
361+
int flv = mcPart.pdgCode() / std::pow(10, static_cast<int>(std::log10(std::abs(mcPart.pdgCode()))));
362+
if (abs(flv) == 5 && mcPart.pdgCode() < 1000)
363+
flv = -flv;
364+
for (int i = (mcPart.mothers_first_as<aod::McParticles>()).globalIndex(); i <= (mcPart.mothers_last_as<aod::McParticles>()).globalIndex(); i++) { // loop over the lund string
365+
for (auto mctrack : mctracks) {
366+
if (mctrack.globalIndex() != i) {
367+
continue;
368+
}
369+
if ((mctrack.pdgCode() != flv) && (abs(mctrack.pdgCode()) < abs(flv) * 1000)) {
370+
continue;
371+
}
372+
while (mctrack.has_mothers()) {
373+
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
374+
auto mother = (abs(motherflv) == abs(flv)) ? (mctrack.mothers_first_as<aod::McParticles>()) : (mctrack.mothers_last_as<aod::McParticles>());
375+
if ((mother.pdgCode() != mctrack.pdgCode()) && (abs(mctrack.pdgCode()) < 10)) { // both mother is not the the quark with same flavor
376+
mcNum = mctrack.globalIndex();
377+
return mcNum;
378+
}
379+
mctrack = mother;
380+
}
381+
}
382+
}
383+
return 0;
384+
}
385+
bool Corr(const McMuons::iterator& muon1, const McMuons::iterator& muon2, aod::McParticles const& mcParts)
386+
{
387+
388+
int moth11(0), moth12(0), moth21(1), moth22(1);
389+
int anc1 = traceAncestor(muon1, mcParts);
390+
int anc2 = traceAncestor(muon2, mcParts);
391+
if (anc1 == 0 || anc2 == 0) {
392+
return false;
393+
}
394+
for (auto mcPart : mcParts) {
395+
if (mcPart.globalIndex() == anc1) {
396+
moth11 = (mcPart.mothers_first_as<aod::McParticles>()).globalIndex();
397+
moth12 = (mcPart.mothers_last_as<aod::McParticles>()).globalIndex();
398+
}
399+
if (mcPart.globalIndex() == anc2) {
400+
moth21 = (mcPart.mothers_first_as<aod::McParticles>()).globalIndex();
401+
moth22 = (mcPart.mothers_last_as<aod::McParticles>()).globalIndex();
402+
}
403+
}
404+
if ((moth11 == moth21) && (moth12 == moth22)) {
405+
return true;
406+
}
407+
return false; // uncorrelated
408+
}
409+
void fillPairs(const McMuons::iterator& muon, const McMuons::iterator& muon2, aod::McParticles const& mcParts)
410+
{
411+
if (trackType != 3) {
412+
return;
413+
}
414+
float mm = o2::constants::physics::MassMuon;
415+
416+
const auto mask1(getMask(muon));
417+
const auto mask2(getMask(muon2));
418+
419+
ROOT::Math::PtEtaPhiMVector mu1Vec(muon.pt(), muon.eta(), muon.phi(), mm);
420+
ROOT::Math::PtEtaPhiMVector mu2Vec(muon2.pt(), muon2.eta(), muon2.phi(), mm);
421+
ROOT::Math::PtEtaPhiMVector dimuVec = mu1Vec + mu2Vec;
422+
auto InvM = dimuVec.M();
423+
424+
if (!muon.has_mcParticle() || !muon2.has_mcParticle()) {
425+
return;
426+
}
427+
auto mcPart1(muon.mcParticle());
428+
auto mcPart2(muon2.mcParticle());
429+
430+
ROOT::Math::PtEtaPhiMVector mu1VecGen(mcPart1.pt(), mcPart1.eta(), mcPart1.phi(), mm);
431+
ROOT::Math::PtEtaPhiMVector mu2VecGen(mcPart2.pt(), mcPart2.eta(), mcPart2.phi(), mm);
432+
ROOT::Math::PtEtaPhiMVector dimuVecGen = mu1VecGen + mu2VecGen;
433+
auto InvMGen = dimuVecGen.M();
434+
435+
if (isMuon(mask1) && isMuon(mask2)) {
436+
registry.fill(HIST("h1MuonMass"), InvM);
437+
registry.fill(HIST("h1MuonMassGen"), InvMGen);
438+
}
439+
if (Corr(muon, muon2, mcParts) && isBeautyMu(mask1) && isBeautyMu(mask2)) {
440+
registry.fill(HIST("h1BeautyMass"), InvM);
441+
registry.fill(HIST("h1BeautyMassGen"), InvMGen);
442+
} else {
443+
registry.fill(HIST("h1OtherMass"), InvM);
444+
registry.fill(HIST("h1OtherMassGen"), InvMGen);
445+
}
446+
}
447+
328448
void process(MyCollisions::iterator const& collision,
329449
McMuons const& muons,
330-
aod::McParticles const&)
450+
aod::McParticles const& mcParts)
331451
{
332452
// event selections
333453
if (std::abs(collision.posZ()) > edgeZ) {
@@ -347,11 +467,15 @@ struct HfTaskSingleMuonSource {
347467
if ((eta >= etaUp) || (eta < etaLow)) {
348468
continue;
349469
}
350-
if ((rAbs >= rAbsMax) || (rAbs < rAbsMin)) {
351-
continue;
470+
if ((rAbs >= rAbsMid) || (rAbs < rAbsMin)) {
471+
if (pDca >= pDcaMax || pDca < 0) {
472+
continue;
473+
}
352474
}
353-
if (pDca >= pDcaMax) {
354-
continue;
475+
if ((rAbs >= rAbsMax) || (rAbs < rAbsMid)) {
476+
if (pDca >= pDcaMax2 || pDca < 0) {
477+
continue;
478+
}
355479
}
356480
if ((muon.chi2() >= 1e6) || (muon.chi2() < 0)) {
357481
continue;
@@ -360,6 +484,48 @@ struct HfTaskSingleMuonSource {
360484
continue;
361485
}
362486
fillHistograms(muon);
487+
if (pairSource) {
488+
if (muon.pt() < ptLow) {
489+
continue;
490+
}
491+
for (const auto& muon2 : muons) {
492+
if (muon2.sign() != muon.sign()) {
493+
continue;
494+
}
495+
if (muon2.globalIndex() <= muon.globalIndex()) {
496+
continue;
497+
}
498+
// muon selections
499+
if (muon2.trackType() != trackType) {
500+
continue;
501+
}
502+
if (muon2.pt() < ptLow) {
503+
continue;
504+
}
505+
const auto eta2(muon2.eta()), pDca2(muon2.pDca()), rAbs2(muon2.rAtAbsorberEnd());
506+
if ((eta2 >= etaUp) || (eta2 < etaLow)) {
507+
continue;
508+
}
509+
if ((rAbs2 >= rAbsMid) || (rAbs2 < rAbsMin)) {
510+
if (pDca2 >= pDcaMax || pDca2 < 0) {
511+
continue;
512+
}
513+
}
514+
if ((rAbs2 >= rAbsMax) || (rAbs2 < rAbsMid)) {
515+
if (pDca2 >= pDcaMax2 || pDca2 < 0) {
516+
continue;
517+
}
518+
}
519+
520+
if ((muon2.chi2() >= 1e6) || (muon2.chi2() < 0)) {
521+
continue;
522+
}
523+
if ((muon2.chi2MatchMCHMID() >= 1e6) || (muon2.chi2MatchMCHMID() < 0)) {
524+
continue;
525+
}
526+
fillPairs(muon, muon2, mcParts);
527+
}
528+
}
363529
} // loop over muons
364530
}
365531
};

0 commit comments

Comments
 (0)