Skip to content

Commit 11246e6

Browse files
fixed syntax and bugs
1 parent 811dbd6 commit 11246e6

File tree

1 file changed

+87
-54
lines changed

1 file changed

+87
-54
lines changed

PWGLF/TableProducer/Nuspex/he3HadronFemto.cxx

Lines changed: 87 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -382,13 +382,14 @@ struct he3HadronFemto {
382382
mQaRegistry.fill(HIST("hEvents"), 0);
383383
mQaRegistry.fill(HIST("hVtxZBefore"), collision.posZ());
384384

385+
auto bc = collision.template bc_as<aod::BCsWithTimestamps>();
386+
initCCDB(bc);
387+
385388
if constexpr (isMC) {
386389
if (/*!collision.sel8() ||*/ std::abs(collision.posZ()) > settingCutVertex) {
387390
return false;
388391
}
389392
} else {
390-
auto bc = collision.template bc_as<aod::BCsWithTimestamps>();
391-
initCCDB(bc);
392393

393394
if (!collision.sel8() || std::abs(collision.posZ()) > settingCutVertex) {
394395
return false;
@@ -528,53 +529,66 @@ struct he3HadronFemto {
528529

529530
// ==================================================================================================================
530531

531-
template <typename Ttrack, typename Tcollisions, typename Ttracks>
532-
bool fillCandidateInfo(const Ttrack& trackHe3, const Ttrack& trackHad, const CollBracket& collBracket, const Tcollisions& collisions, He3HadCandidate& he3Hadcand, const Ttracks& /*trackTable*/, bool isMixedEvent)
532+
template <typename Tcollisions>
533+
std::array<float, 3> getCollisionVertex(const Tcollisions& collisions, int32_t collisionID)
533534
{
534-
const int numCoordinates = 3;
535-
if (!isMixedEvent) {
536-
auto trackCovHe3 = getTrackParCov(trackHe3);
537-
auto trackCovHad = getTrackParCov(trackHad);
538-
int nCand = CommonInite;
539-
try {
540-
nCand = mFitter.process(trackCovHe3, trackCovHad);
541-
} catch (...) {
542-
LOG(error) << "Exception caught in DCA fitter process call!";
543-
return false;
544-
}
545-
if (nCand == 0) {
546-
return false;
547-
}
535+
auto collision = collisions.rawIteratorAt(collisionID);
536+
std::array<float, 3> collisionVertex = {collision.posX(), collision.posY(), collision.posZ()};
537+
return collisionVertex;
538+
}
548539

549-
// associate collision id as the one that minimises the distance between the vertex and the PCAs of the daughters
550-
double distanceMin = -1;
551-
unsigned int collIdxMin = 0;
552-
const float defaultTodistance = 0.0f;
553-
for (int collIdx = collBracket.getMin(); collIdx <= collBracket.getMax(); collIdx++) {
554-
auto collision = collisions.rawIteratorAt(collIdx);
555-
std::array<float, 3> collVtx = {collision.posX(), collision.posY(), collision.posZ()};
556-
const auto& pca = mFitter.getPCACandidate();
557-
float distance = defaultTodistance;
558-
for (int i = 0; i < numCoordinates; i++) {
559-
distance += (pca[i] - collVtx[i]) * (pca[i] - collVtx[i]);
560-
}
561-
if (distanceMin < 0 || distance < distanceMin) {
562-
distanceMin = distance;
563-
collIdxMin = collIdx;
564-
}
565-
}
540+
template <typename Ttrack, typename Tcollisions>
541+
int32_t getCollisionID(const Ttrack& trackHe3, const Ttrack& trackHad, const CollBracket& collBracket, const Tcollisions& collisions, bool isMixedEvent)
542+
{
543+
if (isMixedEvent) {
544+
return collBracket.getMin();
545+
}
566546

567-
if (!mGoodCollisions[collIdxMin]) {
568-
return false;
547+
auto trackCovHe3 = getTrackParCov(trackHe3);
548+
auto trackCovHad = getTrackParCov(trackHad);
549+
int nCand = CommonInite;
550+
try {
551+
nCand = mFitter.process(trackCovHe3, trackCovHad);
552+
} catch (...) {
553+
LOG(error) << "Exception caught in DCA fitter process call!";
554+
return false;
555+
}
556+
if (nCand == 0) {
557+
return false;
558+
}
559+
560+
// associate collision id as the one that minimises the distance between the vertex and the PCAs of the daughters
561+
double distanceMin = -1;
562+
unsigned int collIdxMin = 0;
563+
const float defaultTodistance = 0.0f;
564+
for (int collIdx = collBracket.getMin(); collIdx <= collBracket.getMax(); collIdx++) {
565+
std::array<float, 3> collisionVertex = getCollisionVertex(collisions, collIdx);
566+
const auto& pca = mFitter.getPCACandidate();
567+
float distance = defaultTodistance;
568+
for (int i = 0; i < 3; i++) {
569+
distance += (pca[i] - collisionVertex[i]) * (pca[i] - collisionVertex[i]);
570+
}
571+
if (distanceMin < 0 || distance < distanceMin) {
572+
distanceMin = distance;
573+
collIdxMin = collIdx;
569574
}
570-
he3Hadcand.collisionID = collIdxMin;
575+
}
571576

572-
} else {
573-
he3Hadcand.collisionID = collBracket.getMin();
577+
if (!mGoodCollisions[collIdxMin]) {
578+
return false;
574579
}
580+
return collIdxMin;
581+
582+
}
583+
584+
template <typename Ttrack, typename Tcollisions, typename Ttracks>
585+
bool fillCandidateInfo(const Ttrack& trackHe3, const Ttrack& trackHad, const CollBracket& collBracket, const Tcollisions& collisions, He3HadCandidate& he3Hadcand, const Ttracks& /*trackTable*/, bool isMixedEvent)
586+
{
587+
he3Hadcand.collisionID = getCollisionID(trackHe3, trackHad, collBracket, collisions, isMixedEvent);
588+
std::array<float, 3> collisionVertex = getCollisionVertex(collisions, he3Hadcand.collisionID);
575589

576590
he3Hadcand.momHe3 = std::array{trackHe3.px(), trackHe3.py(), trackHe3.pz()};
577-
for (int i = 0; i < numCoordinates; i++)
591+
for (int i = 0; i < 3; i++)
578592
he3Hadcand.momHe3[i] = he3Hadcand.momHe3[i] * 2;
579593
he3Hadcand.momHad = std::array{trackHad.px(), trackHad.py(), trackHad.pz()};
580594
float invMass = CommonInite;
@@ -668,7 +682,6 @@ struct he3HadronFemto {
668682
template <typename Mc>
669683
void fillCandidateInfoMC(const Mc& mctrackHe3, const Mc& mctrackHad, He3HadCandidate& he3Hadcand)
670684
{
671-
LOG(info) << "--------------------------Filling candidate info MC";
672685
he3Hadcand.momHe3MC = mctrackHe3.pt() * (mctrackHe3.pdgCode() > 0 ? 1 : -1);
673686
he3Hadcand.etaHe3MC = mctrackHe3.eta();
674687
he3Hadcand.phiHe3MC = mctrackHe3.phi();
@@ -680,6 +693,7 @@ struct he3HadronFemto {
680693
template <typename Mc>
681694
void fillMotherInfoMC(const Mc& mctrackHe3, const Mc& mctrackHad, const Mc& mctrackMother, He3HadCandidate& he3Hadcand)
682695
{
696+
LOG(info) << "Mother track: " << mctrackMother.pdgCode() << " " << mctrackMother.pt() << " " << mctrackMother.eta() << " " << mctrackMother.phi();
683697
he3Hadcand.l4PtMC = mctrackMother.pt() * (mctrackMother.pdgCode() > 0 ? 1 : -1);
684698
const double eLit = mctrackHe3.e() + mctrackHad.e();
685699
he3Hadcand.l4MassMC = std::sqrt(eLit * eLit - mctrackMother.p() * mctrackMother.p());
@@ -848,15 +862,16 @@ struct he3HadronFemto {
848862
template <typename TmcParticle>
849863
void setMcParticleFlag(const TmcParticle& mcParticle, std::vector<unsigned int>& mothers, uint8_t& flag)
850864
{
851-
if (mcParticle.isPhysicalPrimary() == true) {
865+
if (mcParticle.isPhysicalPrimary()) {
852866

853867
flag |= ParticleFlags::kPhysicalPrimary;
854868
if (!mcParticle.has_mothers()) {
855869
return;
856870
}
857871

858-
for (const auto& mother : mcParticle.mothers_as<aod::McParticles>()) {
872+
for (const auto& mother : mcParticle.template mothers_as<aod::McParticles>()) {
859873
mothers.push_back(mother.globalIndex());
874+
LOG(info) << "Mother PDG code: " << mother.pdgCode() << ", global index: " << mother.globalIndex();
860875
if (std::abs(mother.pdgCode()) == Li4PDG) {
861876
flag |= ParticleFlags::kFromLi4;
862877
} else if (std::abs(mother.pdgCode()) == o2::constants::physics::Pdg::kHyperTriton) {
@@ -873,8 +888,9 @@ struct he3HadronFemto {
873888
return;
874889
}
875890

876-
for (const auto& mother : mothers) {
891+
for (const auto& mother : mcParticle.template mothers_as<aod::McParticles>()) {
877892
mothers.push_back(mother.globalIndex());
893+
LOG(info) << "Mother PDG code: " << mother.pdgCode() << ", global index: " << mother.globalIndex();
878894
if (std::abs(mother.pdgCode()) == Li4PDG) {
879895
flag |= ParticleFlags::kFromLi4;
880896
} else if (std::abs(mother.pdgCode()) == o2::constants::physics::Pdg::kHyperTriton) {
@@ -1014,44 +1030,61 @@ struct he3HadronFemto {
10141030

10151031
He3HadCandidate he3Hadcand;
10161032
McIter motherParticle;
1033+
unsigned int motherIdx;
10171034
std::vector<unsigned int> motherHe3Idxs, motherHadIdxs;
10181035
setMcParticleFlag(mctrackHe3, motherHe3Idxs, he3Hadcand.flagsHe3);
10191036
setMcParticleFlag(mctrackHad, motherHadIdxs, he3Hadcand.flagsHad);
10201037

1038+
bool isMixedPair = true;
1039+
10211040
if ((he3Hadcand.flagsHe3 == ParticleFlags::kPhysicalPrimary && he3Hadcand.flagsHad == ParticleFlags::kPhysicalPrimary)) {
10221041
he3Hadcand.flags |= Flags::kBothPrimaries;
1042+
isMixedPair = false;
10231043

10241044
} else if ((he3Hadcand.flagsHe3 & ParticleFlags::kFromLi4) && (he3Hadcand.flagsHad & ParticleFlags::kFromLi4)) {
1025-
he3Hadcand.flags |= Flags::kFromLi4;
10261045

10271046
std::unordered_set<unsigned int> motherHe3SetIdxs(motherHe3Idxs.begin(), motherHe3Idxs.end());
10281047
for (const auto& motherHadIdx : motherHadIdxs) {
1029-
if (motherHe3SetIdxs.contains(motherHadIdx)) {
1048+
if (!motherHe3SetIdxs.contains(motherHadIdx)) {
10301049
continue;
10311050
}
10321051

10331052
motherParticle = mcParticles.rawIteratorAt(motherHadIdx);
1034-
if (!motherParticle.pdgCode() == Li4PDG || std::abs(motherParticle.y()) > 1) {
1053+
motherIdx = motherHadIdx;
1054+
if (std::abs(motherParticle.pdgCode()) != Li4PDG || std::abs(motherParticle.y()) > 1) {
10351055
continue;
10361056
}
1057+
isMixedPair = false;
1058+
break;
1059+
}
1060+
if (!isMixedPair) {
1061+
he3Hadcand.flags |= Flags::kBothFromLi4;
10371062
}
10381063

10391064
} else if ((he3Hadcand.flagsHe3 & ParticleFlags::kFromHypertriton) && (he3Hadcand.flagsHad & ParticleFlags::kFromHypertriton)) {
1040-
he3Hadcand.flags |= Flags::kFromHypertriton;
10411065

10421066
std::unordered_set<unsigned int> motherHe3SetIdxs(motherHe3Idxs.begin(), motherHe3Idxs.end());
10431067
for (const auto& motherHadIdx : motherHadIdxs) {
1044-
if (motherHe3SetIdxs.contains(motherHadIdx)) {
1068+
if (!motherHe3SetIdxs.contains(motherHadIdx)) {
10451069
continue;
10461070
}
10471071

10481072
motherParticle = mcParticles.rawIteratorAt(motherHadIdx);
1049-
if (!motherParticle.pdgCode() == o2::constants::physics::Pdg::kHyperTriton || std::abs(motherParticle.y()) > 1) {
1073+
motherIdx = motherHadIdx;
1074+
if (std::abs(motherParticle.pdgCode()) != o2::constants::physics::Pdg::kHyperTriton || std::abs(motherParticle.y()) > 1) {
10501075
continue;
10511076
}
1077+
isMixedPair = false;
1078+
break;
10521079
}
10531080

1054-
} else {
1081+
if (!isMixedPair) {
1082+
he3Hadcand.flags |= Flags::kBothFromHypertriton;
1083+
}
1084+
1085+
}
1086+
1087+
if (isMixedPair) {
10551088
he3Hadcand.flags |= Flags::kMixedPair;
10561089
}
10571090

@@ -1061,8 +1094,8 @@ struct he3HadronFemto {
10611094
fillCandidateInfoMC(mctrackHe3, mctrackHad, he3Hadcand);
10621095

10631096
if ((he3Hadcand.flags == Flags::kBothFromLi4) || (he3Hadcand.flags == Flags::kBothFromHypertriton)) {
1064-
fillMotherInfoMC(mctrackHe3, mctrackHad, mothertrack, he3Hadcand);
1065-
filledMothers.push_back(motherParticle.globalIndex());
1097+
fillMotherInfoMC(mctrackHe3, mctrackHad, motherParticle, he3Hadcand);
1098+
filledMothers.push_back(motherIdx);
10661099
}
10671100

10681101
fillHistograms(he3Hadcand);

0 commit comments

Comments
 (0)