Skip to content

Commit 4a0e2d7

Browse files
prchakraprchakra
authored
[PWGCF] FemtoUniverse: Add function to pair MCTruth particles
1 parent 894709e commit 4a0e2d7

File tree

1 file changed

+232
-20
lines changed

1 file changed

+232
-20
lines changed

‎PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackSpherHarMultKtExtended.cxx‎

Lines changed: 232 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -61,6 +61,7 @@ static const float cutsTable[nPart][nCuts]{
6161
struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended {
6262

6363
Service<o2::framework::O2DatabasePDG> pdg;
64+
Service<o2::framework::O2DatabasePDG> pdgMC;
6465

6566
/// Particle selection part
6667

@@ -101,6 +102,9 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended {
101102
SliceCache cache;
102103
Preslice<FilteredFemtoFullParticles> perCol = aod::femtouniverseparticle::fdCollisionId;
103104

105+
using FemtoTruthParticles = soa::Join<aod::FDParticles, aod::FDExtParticles, aod::FDMCLabels>;
106+
Preslice<FemtoTruthParticles> perColMCTruth = aod::femtouniverseparticle::fdCollisionId;
107+
104108
/// Particle 1
105109
struct : o2::framework::ConfigurableGroup {
106110
Configurable<int> ConfPDGCodePartOne{"ConfPDGCodePartOne", 211, "Particle 1 - PDG code"};
@@ -114,6 +118,8 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended {
114118
/// Partition for particle 1
115119
Partition<FilteredFemtoFullParticles> partsOne = (aod::femtouniverseparticle::partType == uint8_t(aod::femtouniverseparticle::ParticleType::kTrack)) && aod::femtouniverseparticle::sign == as<int8_t>(trackonefilter.ConfChargePart1) && aod::femtouniverseparticle::pt < trackonefilter.ConfPtHighPart1 && aod::femtouniverseparticle::pt > trackonefilter.ConfPtLowPart1;
116120
Partition<FilteredFemtoRecoParticles> partsOneMC = (aod::femtouniverseparticle::partType == uint8_t(aod::femtouniverseparticle::ParticleType::kTrack)) && aod::femtouniverseparticle::sign == as<int8_t>(trackonefilter.ConfChargePart1) && aod::femtouniverseparticle::pt < trackonefilter.ConfPtHighPart1 && aod::femtouniverseparticle::pt > trackonefilter.ConfPtLowPart1;
121+
Partition<FemtoTruthParticles> partsOneMCTruth = aod::femtouniverseparticle::partType == uint8_t(aod::femtouniverseparticle::ParticleType::kMCTruthTrack);
122+
117123
//
118124

119125
/// Histogramming for particle 1
@@ -133,6 +139,7 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended {
133139
/// Partition for particle 2
134140
Partition<FilteredFemtoFullParticles> partsTwo = (aod::femtouniverseparticle::partType == uint8_t(aod::femtouniverseparticle::ParticleType::kTrack)) && (aod::femtouniverseparticle::sign == as<int8_t>(tracktwofilter.ConfChargePart2)) && aod::femtouniverseparticle::pt < tracktwofilter.ConfPtHighPart2 && aod::femtouniverseparticle::pt > tracktwofilter.ConfPtLowPart2;
135141
Partition<FilteredFemtoRecoParticles> partsTwoMC = aod::femtouniverseparticle::partType == uint8_t(aod::femtouniverseparticle::ParticleType::kTrack) && (aod::femtouniverseparticle::sign == as<int8_t>(tracktwofilter.ConfChargePart2)) && aod::femtouniverseparticle::pt < tracktwofilter.ConfPtHighPart2 && aod::femtouniverseparticle::pt > tracktwofilter.ConfPtLowPart2;
142+
Partition<FemtoTruthParticles> partsTwoMCTruth = aod::femtouniverseparticle::partType == uint8_t(aod::femtouniverseparticle::ParticleType::kMCTruthTrack);
136143

137144
/// Histogramming for particle 2
138145
FemtoUniverseParticleHisto<aod::femtouniverseparticle::ParticleType::kTrack, 2> trackHistoPartTwo;
@@ -695,6 +702,138 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended {
695702
}
696703
PROCESS_SWITCH(femtoUniversePairTaskTrackTrackSpherHarMultKtExtended, processSameEventMC, "Enable processing same event for Monte Carlo", false);
697704

705+
/// This function processes the same event and takes care of all the histogramming
706+
/// \todo the trivial loops over the tracks should be factored out since they will be common to all combinations of T-T, T-V0, V0-V0, ...
707+
/// @tparam PartitionType
708+
/// @tparam PartType
709+
/// @tparam isMC: enables Monte Carlo truth specific histograms
710+
/// @param groupPartsOne partition for the first particle passed by the process function
711+
/// @param groupPartsTwo partition for the second particle passed by the process function
712+
/// @param parts femtoUniverseParticles table (in case of Monte Carlo joined with FemtoUniverseMCLabels)
713+
/// @param magFieldTesla magnetic field of the collision
714+
/// @param multCol multiplicity of the collision
715+
template <bool isMC, typename PartitionType>
716+
void doSameEventMCTruth(PartitionType groupPartsOne, PartitionType groupPartsTwo, int multCol, int ContType, bool fillQA)
717+
{
718+
719+
/// Histogramming same event
720+
if ((ContType == 1 || ContType == 2) && fillQA) {
721+
for (const auto& part : groupPartsOne) {
722+
if (part.partType() == uint8_t(aod::femtouniverseparticle::ParticleType::kMCTruthTrack)) {
723+
int pdgCode = static_cast<int>(part.tempFitVar());
724+
const auto& pdgParticle = pdgMC->GetParticle(pdgCode);
725+
if (pdgParticle) {
726+
trackHistoPartOne.fillQA<isMC, false>(part);
727+
}
728+
}
729+
}
730+
}
731+
732+
if ((ContType == 1 || ContType == 3) && fillQA) {
733+
for (const auto& part : groupPartsTwo) {
734+
if (part.partType() == uint8_t(aod::femtouniverseparticle::ParticleType::kMCTruthTrack)) {
735+
int pdgCode = static_cast<int>(part.tempFitVar());
736+
const auto& pdgParticle = pdgMC->GetParticle(pdgCode);
737+
if (pdgParticle) {
738+
trackHistoPartTwo.fillQA<isMC, false>(part);
739+
}
740+
}
741+
}
742+
}
743+
744+
if (ContType == 1) {
745+
746+
/// Now build the combinations for non-identical particle pairs
747+
for (const auto& [p1, p2] : combinations(CombinationsFullIndexPolicy(groupPartsOne, groupPartsTwo))) {
748+
749+
int pdgCodePartOne = static_cast<int>(p1.tempFitVar());
750+
const auto& pdgParticleOne = pdgMC->GetParticle(pdgCodePartOne);
751+
int pdgCodePartTwo = static_cast<int>(p2.tempFitVar());
752+
const auto& pdgParticleTwo = pdgMC->GetParticle(pdgCodePartTwo);
753+
if (pdgParticleOne && pdgParticleTwo && (pdgCodePartOne == trackonefilter.ConfPDGCodePartOne) && (pdgCodePartTwo == tracktwofilter.ConfPDGCodePartTwo)) {
754+
float kT = FemtoUniverseMath::getkT(p1, mass1, p2, mass2);
755+
sameEventMultCont.fillMultNumDen(p1, p2, femto_universe_sh_container::EventType::same, 2, multCol, kT, ConfIsIden);
756+
}
757+
}
758+
} else {
759+
for (const auto& [p1, p2] : combinations(CombinationsStrictlyUpperIndexPolicy(groupPartsOne, groupPartsOne))) {
760+
761+
int pdgCodePartOne = static_cast<int>(p1.tempFitVar());
762+
const auto& pdgParticleOne = pdgMC->GetParticle(pdgCodePartOne);
763+
int pdgCodePartTwo = static_cast<int>(p2.tempFitVar());
764+
const auto& pdgParticleTwo = pdgMC->GetParticle(pdgCodePartTwo);
765+
766+
float kT = FemtoUniverseMath::getkT(p1, mass1, p2, mass2);
767+
double rand;
768+
rand = randgen->Rndm();
769+
std::vector<double> f3d;
770+
if (pdgParticleOne && pdgParticleTwo && (pdgCodePartOne == trackonefilter.ConfPDGCodePartOne) && (pdgCodePartTwo == tracktwofilter.ConfPDGCodePartTwo)) {
771+
772+
switch (ContType) {
773+
case 2: {
774+
if (rand > 0.5) {
775+
sameEventMultContPP.fillMultNumDen(p1, p2, femto_universe_sh_container::EventType::same, 2, multCol, kT, ConfIsIden);
776+
} else if (rand <= 0.5) {
777+
sameEventMultContPP.fillMultNumDen(p2, p1, femto_universe_sh_container::EventType::same, 2, multCol, kT, ConfIsIden);
778+
}
779+
break;
780+
}
781+
782+
case 3: {
783+
if (rand > 0.5) {
784+
sameEventMultContMM.fillMultNumDen(p1, p2, femto_universe_sh_container::EventType::same, 2, multCol, kT, ConfIsIden);
785+
} else if (rand <= 0.5) {
786+
sameEventMultContMM.fillMultNumDen(p2, p1, femto_universe_sh_container::EventType::same, 2, multCol, kT, ConfIsIden);
787+
}
788+
break;
789+
}
790+
default:
791+
break;
792+
}
793+
}
794+
}
795+
}
796+
}
797+
798+
/// process function for to call doSameEvent with Monte Carlo
799+
/// \param col subscribe to the collision table (Monte Carlo Reconstructed reconstructed)
800+
/// \param parts subscribe to joined table FemtoUniverseParticles and FemtoUniverseMCLables to access Monte Carlo truth
801+
/// \param FemtoUniverseMCParticles subscribe to the Monte Carlo truth table
802+
void processSameEventMCTruth(o2::aod::FdCollision const& col,
803+
FemtoTruthParticles const&)
804+
{
805+
fillCollision(col, ConfIsCent);
806+
807+
auto thegroupPartsOne = partsOneMCTruth->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache);
808+
auto thegroupPartsTwo = partsTwoMCTruth->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache);
809+
bool fillQA = true;
810+
randgen = new TRandom2(0);
811+
812+
if (ConfIsCent) {
813+
if (cfgProcessPM) {
814+
doSameEventMCTruth<false>(thegroupPartsOne, thegroupPartsTwo, col.multV0M(), 1, fillQA);
815+
}
816+
if (cfgProcessPP) {
817+
doSameEventMCTruth<false>(thegroupPartsOne, thegroupPartsOne, col.multV0M(), 2, fillQA);
818+
}
819+
if (cfgProcessMM) {
820+
doSameEventMCTruth<false>(thegroupPartsTwo, thegroupPartsTwo, col.multV0M(), 3, fillQA);
821+
}
822+
} else {
823+
if (cfgProcessPM) {
824+
doSameEventMCTruth<false>(thegroupPartsOne, thegroupPartsTwo, col.multNtr(), 1, fillQA);
825+
}
826+
if (cfgProcessPP) {
827+
doSameEventMCTruth<false>(thegroupPartsOne, thegroupPartsOne, col.multNtr(), 2, fillQA);
828+
}
829+
if (cfgProcessMM) {
830+
doSameEventMCTruth<false>(thegroupPartsTwo, thegroupPartsTwo, col.multNtr(), 3, fillQA);
831+
}
832+
}
833+
delete randgen;
834+
}
835+
PROCESS_SWITCH(femtoUniversePairTaskTrackTrackSpherHarMultKtExtended, processSameEventMCTruth, "Enable processing same event for MC truth", false);
836+
698837
/// This function processes the mixed event
699838
/// \todo the trivial loops over the collisions and tracks should be factored out since they will be common to all combinations of T-T, T-V0, V0-V0, ...
700839
/// \tparam PartitionType
@@ -864,26 +1003,6 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended {
8641003
}
8651004
PROCESS_SWITCH(femtoUniversePairTaskTrackTrackSpherHarMultKtExtended, processMixedEventNtr, "Enable processing mixed events for centrality", false);
8661005

867-
/// process function for to fill covariance histograms
868-
/// \param col subscribe to the collision table (Data)
869-
/// \param parts subscribe to the femtoUniverseParticleTable
870-
void processCov(soa::Filtered<o2::aod::FdCollisions>::iterator const& /*col*/,
871-
FilteredFemtoFullParticles const& /*parts*/)
872-
{
873-
int JMax = (ConfLMax + 1) * (ConfLMax + 1);
874-
if (cfgProcessMM) {
875-
sameEventMultContMM.fillMultkTCov(femto_universe_sh_container::EventType::same, JMax);
876-
mixedEventMultContMM.fillMultkTCov(femto_universe_sh_container::EventType::mixed, JMax);
877-
} else if (cfgProcessPP) {
878-
sameEventMultContPP.fillMultkTCov(femto_universe_sh_container::EventType::same, JMax);
879-
mixedEventMultContPP.fillMultkTCov(femto_universe_sh_container::EventType::mixed, JMax);
880-
} else if (cfgProcessPM) {
881-
sameEventMultCont.fillMultkTCov(femto_universe_sh_container::EventType::same, JMax);
882-
mixedEventMultCont.fillMultkTCov(femto_universe_sh_container::EventType::mixed, JMax);
883-
}
884-
}
885-
PROCESS_SWITCH(femtoUniversePairTaskTrackTrackSpherHarMultKtExtended, processCov, "Enable processing same event covariance", false);
886-
8871006
/// brief process function for to call doMixedEvent with Monte Carlo
8881007
/// @param cols subscribe to the collisions table (Monte Carlo Reconstructed reconstructed)
8891008
/// @param parts subscribe to joined table FemtoUniverseParticles and FemtoUniverseMCLables to access Monte Carlo truth
@@ -971,6 +1090,99 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended {
9711090
delete randgen;
9721091
}
9731092
PROCESS_SWITCH(femtoUniversePairTaskTrackTrackSpherHarMultKtExtended, processMixedEventMCNtr, "Enable processing mixed events MC", false);
1093+
1094+
/// This function processes the mixed event
1095+
/// \todo the trivial loops over the collisions and tracks should be factored out since they will be common to all combinations of T-T, T-V0, V0-V0, ...
1096+
/// \tparam PartitionType
1097+
/// \tparam PartType
1098+
/// \tparam isMC: enables Monte Carlo truth specific histograms
1099+
/// \param groupPartsOne partition for the first particle passed by the process function
1100+
/// \param groupPartsTwo partition for the second particle passed by the process function
1101+
/// \param parts femtoUniverseParticles table (in case of Monte Carlo joined with FemtoUniverseMCLabels)
1102+
/// \param magFieldTesla magnetic field of the collision
1103+
/// \param multCol multiplicity of the collision
1104+
template <bool isMC, typename PartitionType>
1105+
void doMixedEventMCTruth(PartitionType groupPartsOne, PartitionType groupPartsTwo, int multCol, int ContType)
1106+
{
1107+
1108+
for (const auto& [p1, p2] : combinations(CombinationsFullIndexPolicy(groupPartsOne, groupPartsTwo))) {
1109+
1110+
int pdgCodePartOne = static_cast<int>(p1.tempFitVar());
1111+
const auto& pdgParticleOne = pdgMC->GetParticle(pdgCodePartOne);
1112+
int pdgCodePartTwo = static_cast<int>(p2.tempFitVar());
1113+
const auto& pdgParticleTwo = pdgMC->GetParticle(pdgCodePartTwo);
1114+
1115+
float kT = FemtoUniverseMath::getkT(p1, mass1, p2, mass2);
1116+
double rand;
1117+
rand = randgen->Rndm();
1118+
if (pdgParticleOne && pdgParticleTwo && (pdgCodePartOne == trackonefilter.ConfPDGCodePartOne) && (pdgCodePartTwo == tracktwofilter.ConfPDGCodePartTwo)) {
1119+
1120+
switch (ContType) {
1121+
case 1: {
1122+
if (rand > 0.5) {
1123+
mixedEventMultCont.fillMultNumDen(p1, p2, femto_universe_sh_container::EventType::mixed, 2, multCol, kT, ConfIsIden);
1124+
} else {
1125+
mixedEventMultCont.fillMultNumDen(p2, p1, femto_universe_sh_container::EventType::mixed, 2, multCol, kT, ConfIsIden);
1126+
}
1127+
break;
1128+
}
1129+
1130+
case 2: {
1131+
if (rand > 0.5) {
1132+
mixedEventMultContPP.fillMultNumDen(p1, p2, femto_universe_sh_container::EventType::mixed, 2, multCol, kT, ConfIsIden);
1133+
} else {
1134+
mixedEventMultContPP.fillMultNumDen(p2, p1, femto_universe_sh_container::EventType::mixed, 2, multCol, kT, ConfIsIden);
1135+
}
1136+
break;
1137+
}
1138+
1139+
case 3: {
1140+
if (rand > 0.5) {
1141+
mixedEventMultContMM.fillMultNumDen(p1, p2, femto_universe_sh_container::EventType::mixed, 2, multCol, kT, ConfIsIden);
1142+
} else {
1143+
mixedEventMultContMM.fillMultNumDen(p2, p1, femto_universe_sh_container::EventType::mixed, 2, multCol, kT, ConfIsIden);
1144+
}
1145+
break;
1146+
}
1147+
default:
1148+
break;
1149+
}
1150+
}
1151+
}
1152+
}
1153+
1154+
/// process function for to call doMixedEvent with Data
1155+
/// @param cols subscribe to the collisions table (Data)
1156+
/// @param parts subscribe to the femtoUniverseParticleTable
1157+
void processMixedEventNtrMCTruth(o2::aod::FdCollisions const& cols,
1158+
FemtoTruthParticles const&)
1159+
{
1160+
randgen = new TRandom2(0);
1161+
1162+
for (const auto& [collision1, collision2] : soa::selfCombinations(colBinningNtr, ConfNEventsMix, -1, cols, cols)) {
1163+
1164+
const int multiplicityCol = collision1.multNtr();
1165+
MixQaRegistry.fill(HIST("MixingQA/hMECollisionBins"), colBinningNtr.getBin({collision1.posZ(), multiplicityCol}));
1166+
1167+
if (cfgProcessPM) {
1168+
auto groupPartsOne = partsOneMCTruth->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision1.globalIndex(), cache);
1169+
auto groupPartsTwo = partsTwoMCTruth->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision2.globalIndex(), cache);
1170+
doMixedEventMCTruth<false>(groupPartsOne, groupPartsTwo, multiplicityCol, 1);
1171+
}
1172+
if (cfgProcessPP) {
1173+
auto groupPartsOne = partsOneMCTruth->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision1.globalIndex(), cache);
1174+
auto groupPartsTwo = partsOneMCTruth->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision2.globalIndex(), cache);
1175+
doMixedEventMCTruth<false>(groupPartsOne, groupPartsTwo, multiplicityCol, 2);
1176+
}
1177+
if (cfgProcessMM) {
1178+
auto groupPartsOne = partsTwoMCTruth->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision1.globalIndex(), cache);
1179+
auto groupPartsTwo = partsTwoMCTruth->sliceByCached(aod::femtouniverseparticle::fdCollisionId, collision2.globalIndex(), cache);
1180+
doMixedEventMCTruth<false>(groupPartsOne, groupPartsTwo, multiplicityCol, 3);
1181+
}
1182+
}
1183+
delete randgen;
1184+
}
1185+
PROCESS_SWITCH(femtoUniversePairTaskTrackTrackSpherHarMultKtExtended, processMixedEventNtrMCTruth, "Enable processing MC Truth mixed events for multiplicity", false);
9741186
};
9751187

9761188
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)