Skip to content

Commit 9554c7f

Browse files
rahulverma012rahulverma012
authored
Merge branch 'AliceO2Group:master' into hParticleCorrelation
2 parents 061dc53 + d3c92c5 commit 9554c7f

29 files changed

Lines changed: 3575 additions & 564 deletions

‎Common/Core/RecoDecay.h‎

Lines changed: 59 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@
2727
#include <cmath> // std::abs, std::sqrt
2828
#include <cstddef> // std::size_t
2929
#include <cstdint> // intX_t
30+
#include <numeric> // std::iota
3031
#include <tuple> // std::apply
3132
#include <type_traits> // std::decay_t
3233
#include <utility> // std::move
@@ -49,6 +50,10 @@ struct RecoDecay {
4950
};
5051

5152
static constexpr int8_t StatusCodeAfterFlavourOscillation{92}; // decay products after B0(s) flavour oscillation
53+
static constexpr int8_t StatusCodePrimaryMin{81}; // primary particles (after hadronization) have status code >= 81
54+
static constexpr int8_t StatusCodePrimaryMax{89}; // primary particles (after hadronization) have status code <= 89
55+
static constexpr int8_t StatusCodeRhadronMin{101}; // R-hadrons (after hadronization) have status code >= 101
56+
static constexpr int8_t StatusCodeRhadronMax{109}; // R-hadrons (after hadronization) have status code <= 109
5257
static constexpr int PdgQuarkMax{8}; // largest quark PDG code; o2-linter: disable=pdg/explicit-code (t' does not have a named constant.)
5358
static constexpr int PdgBosonMin{PDG_t::kGluon}; // smallest boson (gauge or H) PDG code
5459
static constexpr int PdgBosonMax{37}; // largest boson (gauge or H) PDG code; o2-linter: disable=pdg/explicit-code (H+ does not have a named constant.)
@@ -561,6 +566,46 @@ struct RecoDecay {
561566
return maxNormDeltaIP;
562567
}
563568

569+
template <typename PPart>
570+
static void getOrderedMotherIds(const PPart& particle,
571+
std::vector<int64_t>& motherIdsOrdered,
572+
const bool searchUpToQuark = false)
573+
{
574+
if (!particle.has_mothers()) {
575+
return;
576+
}
577+
578+
if (particle.mothersIds().front() == particle.globalIndex()) {
579+
return;
580+
}
581+
582+
if (!searchUpToQuark) {
583+
motherIdsOrdered.push_back(particle.mothersIds().front());
584+
} else {
585+
if (particle.mothersIds().size() == 2) {
586+
if (particle.mothersIds().front() == particle.mothersIds().back()) {
587+
motherIdsOrdered.push_back(particle.mothersIds().front());
588+
} else if (particle.mothersIds().front() < particle.mothersIds().back()) {
589+
auto absStatusCode = std::abs(particle.getGenStatusCode());
590+
if ((absStatusCode >= StatusCodePrimaryMin && absStatusCode <= StatusCodePrimaryMax) || (absStatusCode >= StatusCodeRhadronMin && absStatusCode <= StatusCodeRhadronMax)) {
591+
motherIdsOrdered.resize(particle.mothersIds().back() - particle.mothersIds().front() + 1);
592+
std::iota(motherIdsOrdered.begin(), motherIdsOrdered.end(), particle.mothersIds().front());
593+
} else {
594+
motherIdsOrdered.push_back(particle.mothersIds().front());
595+
motherIdsOrdered.push_back(particle.mothersIds().back());
596+
}
597+
} else if (particle.mothersIds().front() > particle.mothersIds().back()) {
598+
if (particle.mothersIds().back() != 0) {
599+
motherIdsOrdered.push_back(particle.mothersIds().back());
600+
}
601+
motherIdsOrdered.push_back(particle.mothersIds().front());
602+
}
603+
} else if (particle.mothersIds().size() == 1) {
604+
motherIdsOrdered.push_back(particle.mothersIds().front());
605+
}
606+
}
607+
}
608+
564609
/// Finds the mother of an MC particle by looking for the expected PDG code in the mother chain.
565610
/// \tparam acceptFlavourOscillation switch to accept decays where the mother oscillated (e.g. B0 -> B0bar)
566611
/// \param particlesMC table with MC particles
@@ -599,10 +644,12 @@ struct RecoDecay {
599644
for (auto iPart : arrayIds[-stage]) { // check all the particles that were the mothers at the previous stage, o2-linter: disable=const-ref-in-for-loop (int elements)
600645
auto particleMother = particlesMC.rawIteratorAt(iPart - particlesMC.offset());
601646
if (particleMother.has_mothers()) {
602-
// If searchUpToQuark is false we only take the first mother (since decay products only have one mother)
603-
auto lastMotherIdxToCheck = searchUpToQuark ? particleMother.mothersIds().back() : particleMother.mothersIds().front();
604-
for (auto iMother = particleMother.mothersIds().front(); iMother <= lastMotherIdxToCheck; ++iMother) { // loop over the mother particles of the analysed particle
605-
if (std::find(arrayIdsStage.begin(), arrayIdsStage.end(), iMother) != arrayIdsStage.end()) { // if a mother is still present in the vector, do not check it again
647+
648+
std::vector<int64_t> motherIdsOrdered{};
649+
getOrderedMotherIds(particleMother, motherIdsOrdered, searchUpToQuark);
650+
651+
for (auto const& iMother : motherIdsOrdered) { // loop over the mother particles of the analysed particle
652+
if (std::find(arrayIdsStage.begin(), arrayIdsStage.end(), iMother) != arrayIdsStage.end()) { // if a mother is still present in the vector, do not check it again
606653
continue;
607654
}
608655
auto mother = particlesMC.rawIteratorAt(iMother - particlesMC.offset());
@@ -1059,6 +1106,7 @@ struct RecoDecay {
10591106
while (arrayIds[-stage].size() > 0) {
10601107
// vector of mother indices for the current stage
10611108
std::vector<int64_t> arrayIdsStage{};
1109+
std::vector<int64_t> motherIdsOrdered{};
10621110
for (auto iPart : arrayIds[-stage]) { // check all the particles that were the mothers at the previous stage, o2-linter: disable=const-ref-in-for-loop (int elements)
10631111
auto particleMother = particlesMC.rawIteratorAt(iPart - particlesMC.offset());
10641112
if (particleMother.has_mothers()) {
@@ -1072,8 +1120,9 @@ struct RecoDecay {
10721120
}
10731121
}
10741122

1075-
for (auto iMother = particleMother.mothersIds().front(); iMother <= particleMother.mothersIds().back(); ++iMother) { // loop over the mother particles of the analysed particle
1076-
if (std::find(arrayIdsStage.begin(), arrayIdsStage.end(), iMother) != arrayIdsStage.end()) { // if a mother is still present in the vector, do not check it again
1123+
getOrderedMotherIds(particleMother, motherIdsOrdered, searchUpToQuark);
1124+
for (auto const& iMother : motherIdsOrdered) { // loop over the mother particles of the analysed particle
1125+
if (std::find(arrayIdsStage.begin(), arrayIdsStage.end(), iMother) != arrayIdsStage.end()) { // if a mother is still present in the vector, do not check it again
10771126
continue;
10781127
}
10791128
auto mother = particlesMC.rawIteratorAt(iMother - particlesMC.offset());
@@ -1161,30 +1210,11 @@ struct RecoDecay {
11611210
auto particleMother = particlesMC.rawIteratorAt(iPart - particlesMC.offset());
11621211
if (particleMother.has_mothers()) {
11631212

1164-
// we break immediately if searchUpToQuark is false and the first mother is a quark or a boson (a hadron should never be the mother of a parton)
1165-
if (!searchUpToQuark) {
1166-
auto mother = particlesMC.rawIteratorAt(particleMother.mothersIds().front() - particlesMC.offset());
1167-
auto pdgParticleIMother = std::abs(mother.pdgCode()); // PDG code of the mother
1168-
if (pdgParticleIMother <= PdgQuarkMax || (pdgParticleIMother >= PdgBosonMin && pdgParticleIMother <= PdgBosonMax)) {
1169-
// auto PDGPaticle = std::abs(particleMother.pdgCode());
1170-
if (
1171-
(pdgParticle / PdgDivisorMeson == PDG_t::kBottom || // b mesons
1172-
pdgParticle / PdgDivisorBaryon == PDG_t::kBottom) // b baryons
1173-
) {
1174-
return OriginType::NonPrompt; // beauty
1175-
}
1176-
if (
1177-
(pdgParticle / PdgDivisorMeson == PDG_t::kCharm || // c mesons
1178-
pdgParticle / PdgDivisorBaryon == PDG_t::kCharm) // c baryons
1179-
) {
1180-
return OriginType::Prompt; // charm
1181-
}
1182-
break;
1183-
}
1184-
}
1213+
std::vector<int64_t> motherIdsOrdered{};
1214+
getOrderedMotherIds(particleMother, motherIdsOrdered, searchUpToQuark);
11851215

1186-
for (auto iMother = particleMother.mothersIds().front(); iMother <= particleMother.mothersIds().back(); ++iMother) { // loop over the mother particles of the analysed particle
1187-
if (std::find(arrayIdsStage.begin(), arrayIdsStage.end(), iMother) != arrayIdsStage.end()) { // if a mother is still present in the vector, do not check it again
1216+
for (auto const& iMother : motherIdsOrdered) { // loop over the mother particles of the analysed particle
1217+
if (std::find(arrayIdsStage.begin(), arrayIdsStage.end(), iMother) != arrayIdsStage.end()) { // if a mother is still present in the vector, do not check it again
11881218
continue;
11891219
}
11901220
auto mother = particlesMC.rawIteratorAt(iMother - particlesMC.offset());

‎PWGCF/GenericFramework/Core/FlowContainer.cxx‎

Lines changed: 23 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -415,6 +415,10 @@ TProfile* FlowContainer::GetCorrXXVsMulti(const char* order, int l_pti)
415415
}
416416
delete rethist;
417417
}
418+
if (!retSubset) {
419+
LOGF(error, "Correlation profile is null");
420+
return nullptr;
421+
}
418422
if (fMultiRebin > 0) {
419423
TString temp_name(retSubset->GetName());
420424
TProfile* tempprof = dynamic_cast<TProfile*>(retSubset->Clone("tempProfile"));
@@ -525,9 +529,9 @@ TH1D* FlowContainer::GetCN2VsX(int n, bool onPt, double arg1, double arg2)
525529
corrN2 = GetHistCorrXXVsMulti(Form("%i2", n), static_cast<int>(arg1));
526530
corrN2->SetName(Form("Corr_%s", corrN2->GetName()));
527531
TH1D* rethist = GetCN2(corrN2);
528-
TString* nam = new TString(corrN2->GetName());
532+
const TString nam(corrN2->GetName());
529533
delete corrN2;
530-
rethist->SetName(nam->Data());
534+
rethist->SetName(nam.Data());
531535
if (onPt) {
532536
int bins = fProf->GetXaxis()->FindBin(arg1);
533537
int bins2 = fProf->GetXaxis()->FindBin(arg2);
@@ -543,10 +547,10 @@ TH1D* FlowContainer::GetCN2VsX(int n, bool onPt, double arg1, double arg2)
543547
TH1D* FlowContainer::GetVN2VsX(int n, bool onPt, double arg1, double arg2)
544548
{
545549
TH1D* corrh = GetCN2VsX(n, onPt, arg1, arg2);
546-
TString* nam = new TString(corrh->GetName());
550+
const TString nam(corrh->GetName());
547551
TH1D* rethist = GetVN2(corrh);
548552
delete corrh;
549-
rethist->SetName(nam->Data());
553+
rethist->SetName(nam.Data());
550554
if (onPt) {
551555
int bins = fProf->GetXaxis()->FindBin(arg1);
552556
int bins2 = fProf->GetXaxis()->FindBin(arg2);
@@ -556,7 +560,6 @@ TH1D* FlowContainer::GetVN2VsX(int n, bool onPt, double arg1, double arg2)
556560
} else {
557561
rethist->SetTitle(Form(";#it{N}_{tr};v_{%i}{2}", n));
558562
}
559-
delete nam;
560563
return rethist;
561564
}
562565

@@ -706,10 +709,10 @@ TH1D* FlowContainer::GetCN4VsX(int n, bool onPt, double arg1, double arg2)
706709
corrN4->SetName(Form("Corr_%s", corrN4->GetName()));
707710
}
708711
TH1D* rethist = GetCN4(corrN4, corrN2);
709-
TString* nam = new TString(corrN4->GetName());
712+
const TString nam(corrN4->GetName());
710713
delete corrN2;
711714
delete corrN4;
712-
rethist->SetName(nam->Data());
715+
rethist->SetName(nam.Data());
713716
if (onPt) {
714717
int bins = fProf->GetXaxis()->FindBin(arg1);
715718
int bins2 = fProf->GetXaxis()->FindBin(arg2);
@@ -780,9 +783,9 @@ TH1D* FlowContainer::GetVN4VsX(int n, bool onPt, double arg1, double arg2)
780783
{
781784
TH1D* temph = GetCN4VsX(n, onPt, arg1, arg2);
782785
TH1D* rethist = GetVN4(temph);
783-
TString* nam = new TString(temph->GetName());
786+
const TString nam(temph->GetName());
784787
delete temph;
785-
rethist->SetName(nam->Data());
788+
rethist->SetName(nam.Data());
786789
if (onPt) {
787790
int bins = fProf->GetXaxis()->FindBin(arg1);
788791
int bins2 = fProf->GetXaxis()->FindBin(arg2);
@@ -812,9 +815,9 @@ TH1D* FlowContainer::GetCN6VsX(int n, bool onPt, double arg1, double arg2)
812815
TH1D* rethist = GetCN6(corrN6, corrN4, corrN2);
813816
delete corrN2;
814817
delete corrN4;
815-
TString* nam = new TString(corrN6->GetName());
818+
const TString nam(corrN6->GetName());
816819
delete corrN6;
817-
rethist->SetName(nam->Data());
820+
rethist->SetName(nam.Data());
818821
if (onPt) {
819822
int bins = fProf->GetXaxis()->FindBin(arg1);
820823
int bins2 = fProf->GetXaxis()->FindBin(arg2);
@@ -830,9 +833,9 @@ TH1D* FlowContainer::GetVN6VsX(int n, bool onPt, double arg1, double arg2)
830833
{
831834
TH1D* temph = GetCN6VsX(n, onPt, arg1, arg2);
832835
TH1D* rethist = GetVN6(temph);
833-
TString* nam = new TString(temph->GetName());
836+
const TString nam(temph->GetName());
834837
delete temph;
835-
rethist->SetName(nam->Data());
838+
rethist->SetName(nam.Data());
836839
if (onPt) {
837840
int bins = fProf->GetXaxis()->FindBin(arg1);
838841
int bins2 = fProf->GetXaxis()->FindBin(arg2);
@@ -866,9 +869,9 @@ TH1D* FlowContainer::GetCN8VsX(int n, bool onPt, double arg1, double arg2)
866869
delete corrN2;
867870
delete corrN4;
868871
delete corrN6;
869-
TString* nam = new TString(corrN8->GetName());
872+
const TString nam(corrN8->GetName());
870873
delete corrN8;
871-
rethist->SetName(nam->Data());
874+
rethist->SetName(nam.Data());
872875
if (onPt) {
873876
int bins = fProf->GetXaxis()->FindBin(arg1);
874877
int bins2 = fProf->GetXaxis()->FindBin(arg2);
@@ -884,9 +887,9 @@ TH1D* FlowContainer::GetVN8VsX(int n, bool onPt, double arg1, double arg2)
884887
{
885888
TH1D* temph = GetCN8VsX(n, onPt, arg1, arg2);
886889
TH1D* rethist = GetVN8(temph);
887-
TString* nam = new TString(temph->GetName());
890+
const TString nam(temph->GetName());
888891
delete temph;
889-
rethist->SetName(nam->Data());
892+
rethist->SetName(nam.Data());
890893
if (onPt) {
891894
int bins = fProf->GetXaxis()->FindBin(arg1);
892895
int bins2 = fProf->GetXaxis()->FindBin(arg2);
@@ -928,7 +931,7 @@ TProfile* FlowContainer::GetRefFlowProfile(const char* order, double m1, double
928931
if (nStopBin < nStartBin)
929932
nStopBin = fProf->GetXaxis()->GetNbins();
930933
int nBins = nStopBin - nStartBin + 1;
931-
double* l_bins = new double[nBins + 1];
934+
std::vector<double> l_bins(nBins + 1);
932935
for (int i = 0; i <= nBins; i++)
933936
l_bins[i] = i;
934937
TProfile* retpf = 0;
@@ -939,7 +942,7 @@ TProfile* FlowContainer::GetRefFlowProfile(const char* order, double m1, double
939942
while (fIDName.Tokenize(l_name, l_pos)) {
940943
l_name.Append(order);
941944
int ybin = fProf->GetYaxis()->FindBin(l_name.Data());
942-
TProfile* tempprof = rhSubset->GetSubset(kTRUE, "tempprof", ybin, ybin, nBins, l_bins);
945+
TProfile* tempprof = rhSubset->GetSubset(kTRUE, "tempprof", ybin, ybin, nBins, l_bins.data());
943946
if (!retpf)
944947
retpf = dynamic_cast<TProfile*>(tempprof->Clone("RefFlowProf"));
945948
else
@@ -1237,5 +1240,5 @@ double* FlowContainer::GetMultiRebin(int& nbins)
12371240
double* retBins = new double[fMultiRebin + 1];
12381241
for (int i = 0; i <= nbins; i++)
12391242
retBins[i] = fMultiRebinEdges[i];
1240-
return fMultiRebinEdges;
1243+
return retBins;
12411244
}

0 commit comments

Comments
 (0)