Skip to content

Commit a707602

Browse files
MaolinZHfjonasALICE
MaolinZH
authored and
fjonasALICE
committed
[PWGHF] add new classifications for like-sign muon pairs (#15683)
1 parent 249446f commit a707602

File tree

1 file changed

+61
-42
lines changed

1 file changed

+61
-42
lines changed

‎PWGHF/HFL/Tasks/taskSingleMuonSource.cxx‎

Lines changed: 61 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -8,10 +8,10 @@
88
// In applying this license CERN does not waive the privileges and immunities
99
// granted to it by virtue of its status as an Intergovernmental Organization
1010
// or submit itself to any jurisdiction.
11-
//
12-
// \file taskSingleMuonSource.cxx
13-
// \brief Task used to seperate single muons source in Monte Carlo simulation.
14-
// \author Maolin Zhang <maolin.zhang@cern.ch>, CCNU
11+
///
12+
/// \file taskSingleMuonSource.cxx
13+
/// \brief Task used to seperate single muons source in Monte Carlo simulation.
14+
/// \author Maolin Zhang <maolin.zhang@cern.ch>, CCNU
1515

1616
#include "Common/Core/RecoDecay.h"
1717
#include "Common/DataModel/TrackSelectionTables.h"
@@ -24,11 +24,11 @@
2424
#include <Framework/HistogramRegistry.h>
2525
#include <Framework/HistogramSpec.h>
2626
#include <Framework/InitContext.h>
27+
#include <Framework/O2DatabasePDGPlugin.h>
2728
#include <Framework/OutputObjHeader.h>
2829
#include <Framework/runDataProcessing.h>
2930

3031
#include <Math/Vector4D.h>
31-
#include <TDatabasePDG.h>
3232
#include <TPDGCode.h>
3333
#include <TString.h>
3434

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

81+
Service<o2::framework::O2DatabasePDG> pdgDB;
82+
8183
double pDcaMax = 594.0; // p*DCA maximum value for small Rab
8284
double pDcaMax2 = 324.0; // p*DCA maximum value for large Rabs
8385
double rAbsMid = 26.5; // R at absorber end minimum value
@@ -116,7 +118,7 @@ struct HfTaskSingleMuonSource {
116118

117119
HistogramConfigSpec const h1ColNumber{HistType::kTH1F, {axisColNumber}};
118120
HistogramConfigSpec const h1Pt{HistType::kTH1F, {axisPt}};
119-
HistogramConfigSpec h1Mass{HistType::kTH1F, {axisMass}};
121+
HistogramConfigSpec const h1Mass{HistType::kTH1F, {axisMass}};
120122
HistogramConfigSpec const h2PtDCA{HistType::kTH2F, {axisPt, axisDCA}};
121123
HistogramConfigSpec const h2PtChi2{HistType::kTH2F, {axisPt, axisChi2}};
122124
HistogramConfigSpec const h2PtDeltaPt{HistType::kTH2F, {axisPt, axisDeltaPt}};
@@ -125,9 +127,11 @@ struct HfTaskSingleMuonSource {
125127
registry.add("h1MuBeforeCuts", "", h1Pt);
126128
registry.add("h1MuonMass", "", h1Mass);
127129
registry.add("h1BeautyMass", "", h1Mass);
130+
registry.add("h1CorrBeautyMass", "", h1Mass);
128131
registry.add("h1OtherMass", "", h1Mass);
129132
registry.add("h1MuonMassGen", "", h1Mass);
130133
registry.add("h1BeautyMassGen", "", h1Mass);
134+
registry.add("h1CorrBeautyMassGen", "", h1Mass);
131135
registry.add("h1OtherMassGen", "", h1Mass);
132136
for (const auto& src : muonSources) {
133137
registry.add(Form("h1%sPt", src.Data()), "", h1Pt);
@@ -141,6 +145,8 @@ struct HfTaskSingleMuonSource {
141145
uint8_t getMask(const McMuons::iterator& muon)
142146
{
143147
uint8_t mask(0);
148+
const int diquarkEdge = 1000;
149+
const int hadronEdge = 10000;
144150
if (muon.has_mcParticle()) {
145151
SETBIT(mask, IsIdentified);
146152
} else {
@@ -159,7 +165,7 @@ struct HfTaskSingleMuonSource {
159165
mcPart = *(mcPart.mothers_first_as<aod::McParticles>());
160166

161167
const auto pdgAbs(std::abs(mcPart.pdgCode()));
162-
if (pdgAbs < 10 || pdgAbs == 21) {
168+
if (pdgAbs < kElectron || pdgAbs == kGluon) {
163169
break; // Quark and gluon
164170
}
165171

@@ -168,40 +174,39 @@ struct HfTaskSingleMuonSource {
168174
continue;
169175
}
170176

171-
if (pdgAbs == kTauMinus) {
172-
// Tau
177+
if (pdgAbs == kTauMinus) { // Tau
173178
SETBIT(mask, HasTauParent);
174179
continue;
175180
}
176181

177182
const int pdgRem(pdgAbs % 100000);
183+
const int pdgRemRem(pdgRem % 100);
178184

179185
if (pdgRem == kProton) {
180186
continue;
181187
} // Beam particle
182188

183-
if ((pdgRem < 100) || (pdgRem >= 10000)) {
189+
if ((pdgRem < kPi0) || (pdgRem >= hadronEdge)) {
184190
continue;
185191
}
186-
if ((pdgRem % 100 == 1 || pdgRem % 100 == 3) && pdgRem > 1000) { // diquarks
192+
if ((pdgRemRem == kDown || pdgRemRem == kStrange) && pdgRem > diquarkEdge) { // diquarks
187193
continue;
188194
}
189195
// compute the flavor of constituent quark
190196
const int flv(pdgRem / std::pow(10, static_cast<int>(std::log10(pdgRem))));
191-
if (flv > 6) {
197+
if (flv > kTop) {
192198
// no more than 6 flavors
193199
continue;
194200
}
195-
if (flv < 4) {
201+
if (flv < kCharm) {
196202
// light flavor
197203
SETBIT(mask, HasLightParent);
198204
continue;
199205
}
200-
201-
auto* pdgData(TDatabasePDG::Instance()->GetParticle(mcPart.pdgCode()));
206+
auto* pdgData = pdgDB->GetParticle(mcPart.pdgCode());
202207
if ((pdgData != nullptr) && (pdgData->AntiParticle() == nullptr)) {
203208
SETBIT(mask, HasQuarkoniumParent);
204-
} else if (flv == 4) {
209+
} else if (flv == kCharm) {
205210
SETBIT(mask, HasCharmParent);
206211
} else {
207212
SETBIT(mask, HasBeautyParent);
@@ -274,11 +279,13 @@ struct HfTaskSingleMuonSource {
274279
// fill the histograms of each particle types
275280
void fillHistograms(const McMuons::iterator& muon)
276281
{
282+
const int type0 = 0;
283+
const int type2 = 2;
277284
const auto mask(getMask(muon));
278285
const auto pt(muon.pt()), chi2(muon.chi2MatchMCHMFT());
279286
const auto dca(RecoDecay::sqrtSumOfSquares(muon.fwdDcaX(), muon.fwdDcaY()));
280287

281-
if (trackType == 0 || trackType == 2) {
288+
if (trackType == type0 || trackType == type2) {
282289
if (!muon.has_matchMCHTrack()) {
283290
return;
284291
}
@@ -344,6 +351,8 @@ struct HfTaskSingleMuonSource {
344351
int traceAncestor(const McMuons::iterator& muon, aod::McParticles const& mctracks)
345352
{
346353
int mcNum = 0;
354+
const int hadronStatus = 80;
355+
const int diquarkEdge = 1000;
347356
if (!muon.has_mcParticle()) {
348357
return 0;
349358
}
@@ -353,36 +362,38 @@ struct HfTaskSingleMuonSource {
353362
}
354363
while (mcPart.has_mothers()) { // the first hadron after hadronization
355364
auto mother = mcPart.mothers_first_as<aod::McParticles>();
356-
if (std::abs(mother.getGenStatusCode()) < 80) {
365+
if (std::abs(mother.getGenStatusCode()) < hadronStatus) {
357366
break;
358367
}
359368
mcPart = mother;
360369
}
361370
int flv = mcPart.pdgCode() / std::pow(10, static_cast<int>(std::log10(std::abs(mcPart.pdgCode()))));
362-
if (abs(flv) == 5 && mcPart.pdgCode() < 1000)
371+
if (std::abs(flv) == kBottom && mcPart.pdgCode() < diquarkEdge) {
363372
flv = -flv;
373+
}
364374
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) {
375+
for (const auto& mctrack : mctracks) {
366376
if (mctrack.globalIndex() != i) {
367377
continue;
368378
}
369-
if ((mctrack.pdgCode() != flv) && (abs(mctrack.pdgCode()) < abs(flv) * 1000)) {
379+
if ((mctrack.pdgCode() != flv) && (std::abs(mctrack.pdgCode()) < std::abs(flv) * 1000)) {
370380
continue;
371381
}
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();
382+
auto currentTrk = mctrack;
383+
while (currentTrk.has_mothers()) {
384+
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
385+
auto mother = (std::abs(motherflv) == std::abs(flv)) ? (currentTrk.mothers_first_as<aod::McParticles>()) : (currentTrk.mothers_last_as<aod::McParticles>());
386+
if ((mother.pdgCode() != currentTrk.pdgCode()) && (std::abs(currentTrk.pdgCode()) < kElectron)) { // both mother is not the the quark with same flavor
387+
mcNum = currentTrk.globalIndex();
377388
return mcNum;
378389
}
379-
mctrack = mother;
390+
currentTrk = mother;
380391
}
381392
}
382393
}
383394
return 0;
384395
}
385-
bool Corr(const McMuons::iterator& muon1, const McMuons::iterator& muon2, aod::McParticles const& mcParts)
396+
bool isCorr(const McMuons::iterator& muon1, const McMuons::iterator& muon2, aod::McParticles const& mcParts)
386397
{
387398

388399
int moth11(0), moth12(0), moth21(1), moth22(1);
@@ -391,7 +402,7 @@ struct HfTaskSingleMuonSource {
391402
if (anc1 == 0 || anc2 == 0) {
392403
return false;
393404
}
394-
for (auto mcPart : mcParts) {
405+
for (const auto& mcPart : mcParts) {
395406
if (mcPart.globalIndex() == anc1) {
396407
moth11 = (mcPart.mothers_first_as<aod::McParticles>()).globalIndex();
397408
moth12 = (mcPart.mothers_last_as<aod::McParticles>()).globalIndex();
@@ -408,7 +419,8 @@ struct HfTaskSingleMuonSource {
408419
}
409420
void fillPairs(const McMuons::iterator& muon, const McMuons::iterator& muon2, aod::McParticles const& mcParts)
410421
{
411-
if (trackType != 3) {
422+
const int type3 = 3;
423+
if (trackType != type3) {
412424
return;
413425
}
414426
float mm = o2::constants::physics::MassMuon;
@@ -419,7 +431,7 @@ struct HfTaskSingleMuonSource {
419431
ROOT::Math::PtEtaPhiMVector mu1Vec(muon.pt(), muon.eta(), muon.phi(), mm);
420432
ROOT::Math::PtEtaPhiMVector mu2Vec(muon2.pt(), muon2.eta(), muon2.phi(), mm);
421433
ROOT::Math::PtEtaPhiMVector dimuVec = mu1Vec + mu2Vec;
422-
auto InvM = dimuVec.M();
434+
auto invMass = dimuVec.M();
423435

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

435447
if (isMuon(mask1) && isMuon(mask2)) {
436-
registry.fill(HIST("h1MuonMass"), InvM);
437-
registry.fill(HIST("h1MuonMassGen"), InvMGen);
448+
registry.fill(HIST("h1MuonMass"), invMass);
449+
registry.fill(HIST("h1MuonMassGen"), invMassGen);
438450
}
439-
if (Corr(muon, muon2, mcParts) && isBeautyMu(mask1) && isBeautyMu(mask2)) {
440-
registry.fill(HIST("h1BeautyMass"), InvM);
441-
registry.fill(HIST("h1BeautyMassGen"), InvMGen);
451+
if (isBeautyMu(mask1) && isBeautyMu(mask2)) {
452+
registry.fill(HIST("h1BeautyMass"), invMass);
453+
registry.fill(HIST("h1BeautyMassGen"), invMassGen);
454+
if (isCorr(muon, muon2, mcParts)) {
455+
registry.fill(HIST("h1CorrBeautyMass"), invMass);
456+
registry.fill(HIST("h1CorrBeautyMassGen"), invMassGen);
457+
}
442458
} else {
443-
registry.fill(HIST("h1OtherMass"), InvM);
444-
registry.fill(HIST("h1OtherMassGen"), InvMGen);
459+
registry.fill(HIST("h1OtherMass"), invMass);
460+
registry.fill(HIST("h1OtherMassGen"), invMassGen);
445461
}
446462
}
447463

@@ -477,7 +493,10 @@ struct HfTaskSingleMuonSource {
477493
continue;
478494
}
479495
}
480-
if ((muon.chi2() >= 1e6) || (muon.chi2() < 0)) {
496+
if (muon.chi2() < 0) {
497+
continue;
498+
}
499+
if (muon.chi2MatchMCHMID() < 0) {
481500
continue;
482501
}
483502
if (charge != 0 && muon.sign() != charge) {
@@ -517,10 +536,10 @@ struct HfTaskSingleMuonSource {
517536
}
518537
}
519538

520-
if ((muon2.chi2() >= 1e6) || (muon2.chi2() < 0)) {
539+
if (muon2.chi2() < 0) {
521540
continue;
522541
}
523-
if ((muon2.chi2MatchMCHMID() >= 1e6) || (muon2.chi2MatchMCHMID() < 0)) {
542+
if (muon2.chi2MatchMCHMID() < 0) {
524543
continue;
525544
}
526545
fillPairs(muon, muon2, mcParts);

0 commit comments

Comments
 (0)