Skip to content

PWGCF / FemtoUniverse : Fixing bug in two-particle pair kinematics #7628

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
saganatt merged 8 commits into from
Sep 9, 2024
Merged

PWGCF / FemtoUniverse : Fixing bug in two-particle pair kinematics #7628

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
bc00b07
Fixing bug in two-particle pair kinematics
prchakra Sep 9, 2024
32d3b66
Fixing bug in two-particle pair kinematics
prchakra Sep 9, 2024
5af7e31
Fixing bug in two-particle pair kinematics
prchakra Sep 9, 2024
9a090c7
Fixing bug in two-particle pair kinematics
prchakra Sep 9, 2024
41af238
Fixing bug in two-particle pair kinematics
prchakra Sep 9, 2024
467952a
Fixing bug in two-particle pair kinematics
prchakra Sep 9, 2024
65fb506
Fixing error in clang-format
prchakra Sep 9, 2024
51112bf
Fixing error in clang-format
prchakra Sep 9, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 2 additions & 3 deletions PWGCF/FemtoUniverse/Core/FemtoUniverse3DContainer.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,7 +74,6 @@ class FemtoUniverse3DContainer
/// \param mTAxis axis object for the mT axis
/// \param use3dplots Flag to fill 3D plots
/// \param isiden Identical or non-identical particle pair
/// \param islcms LCMS or PRF
template <typename T>
void init_base(std::string folderName, std::string femtoObs1D, std::string femtoObsKout, std::string femtoObsKside, std::string femtoObsKlong, T femtoObsAxis1D, T femtoObsAxisOut, T femtoObsAxisSide, T femtoObsAxisLong, T multAxis, T kTAxis, T mTAxis, T multAxis3D, T mTAxis3D, bool use3dplots, bool isiden)
{
Expand Down Expand Up @@ -214,13 +213,13 @@ class FemtoUniverse3DContainer
/// \param isiden Choosing identical or non-identical pairs
/// \param islcm Choosing LCMS or PRF
template <bool isMC, typename T>
void setPair(T const& part1, T const& part2, const int mult, bool use3dplots, bool isiden, bool islcms)
void setPair(T const& part1, T const& part2, const int mult, bool use3dplots, bool isiden)
{
std::vector<double> f3d;
const float kT = FemtoUniverseMath::getkT(part1, mMassOne, part2, mMassTwo);
const float mT = FemtoUniverseMath::getmT(part1, mMassOne, part2, mMassTwo);

f3d = FemtoUniverseMath::getpairmom3d(part1, mMassOne, part2, mMassTwo, isiden, islcms);
f3d = FemtoUniverseMath::newpairfunc(part1, mMassOne, part2, mMassTwo, isiden);

const float femtoObs1D = f3d[0];
const float femtoObsKout = f3d[1];
Expand Down
97 changes: 39 additions & 58 deletions PWGCF/FemtoUniverse/Core/FemtoUniverseMath.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -149,9 +149,8 @@ class FemtoUniverseMath
/// \param part2 Particle 2
/// \param mass2 Mass of particle 2
/// \param isiden Identical or non-identical particle pair
/// \param islcms LCMS or PRF
template <typename T>
static std::vector<double> getpairmom3d(const T& part1, const float mass1, const T& part2, const float mass2, bool isiden, bool islcms)
static std::vector<double> newpairfunc(const T& part1, const float mass1, const T& part2, const float mass2, bool isiden)
{
const double E1 = sqrt(pow(part1.px(), 2) + pow(part1.py(), 2) + pow(part1.pz(), 2) + pow(mass1, 2));
const double E2 = sqrt(pow(part2.px(), 2) + pow(part2.py(), 2) + pow(part2.pz(), 2) + pow(mass2, 2));
Expand All @@ -165,80 +164,62 @@ class FemtoUniverseMath
const double tPx = trackSum.px();
const double tPy = trackSum.py();
const double tPz = trackSum.pz();
const double tPE = trackSum.E();
const double tE = trackSum.E();

const double tPt = trackSum.pt();
const double tMt = trackSum.mt();
const double tPinv = std::sqrt((tMt * tMt) - (tPt * tPt));
double tPt = (tPx * tPx + tPy * tPy);
double tMt = (tE * tE - tPz * tPz);
double tM = sqrt(tMt - tPt);
tMt = sqrt(tMt);
tPt = sqrt(tPt);

float nullmass = 0.0;
const double m1 = std::max(nullmass, mass1);
const double m2 = std::max(nullmass, mass2);

const double tQinvL = std::pow((E1 - E2), 2) - std::pow((part1.px() - part2.px()), 2) -
std::pow((part1.py() - part2.py()), 2) - std::pow((part1.pz() - part2.pz()), 2);

double tQ = (m1 - m2) / tPinv;
tQ = ::sqrt(tQ * tQ - tQinvL);

const double fKStarCalc = tQ / 2.0;
vect.push_back(fKStarCalc);
double fDKOutLCMS, fDKSideLCMS, fDKLongLCMS;
double fDKOut, fDKSide, fDKLong, fDE;
double px1LCMS, py1LCMS, pz1LCMS;
double px2LCMS, py2LCMS, pz2LCMS;
double kstar;

// Boost to LCMS

const double beta = tPz / tPE;
const double gamma = tPE / tMt;
const double beta = tPz / tE;
const double gamma = tE / tMt;

const double px1L = (part1.px() * tPx + part1.py() * tPy) / tPt;
const double py1L = (-part1.px() * tPy + part1.py() * tPx) / tPt;
const double pz1L = gamma * (part1.pz() - beta * E1);
const double pE1L = gamma * (E1 - beta * part1.pz());
fDKOut = (part1.px() * tPx + part1.py() * tPy) / tPt;
fDKSide = (-part1.px() * tPy + part1.py() * tPx) / tPt;
fDKLong = gamma * (part1.pz() - beta * E1);
fDE = gamma * (E1 - beta * part1.pz());

const double px2L = (part2.px() * tPx + part2.py() * tPy) / tPt;
const double py2L = (-part2.px() * tPy + part2.py() * tPx) / tPt;
const double pz2L = gamma * (part2.pz() - beta * E2);
const double pE2L = gamma * (E2 - beta * part2.pz());
px1LCMS = fDKOut;
py1LCMS = fDKSide;
pz1LCMS = fDKLong;
// pE1LCMS = fDE;

double fDKOutLCMS;
double fDKSideLCMS;
double fDKLongLCMS;
px2LCMS = (part2.px() * tPx + part2.py() * tPy) / tPt;
py2LCMS = (part2.py() * tPx - part2.px() * tPy) / tPt;
pz2LCMS = gamma * (part2.pz() - beta * E2);
// pE2LCMS = gamma * (E2 - beta * part2.pz());

double fDKOutPRF;
double fDKSidePRF;
double fDKLongPRF;

if (!isiden) {
fDKOutLCMS = px1L;
fDKSideLCMS = py1L;
fDKLongLCMS = pz1L;
} else {
fDKOutLCMS = px1L - px2L;
fDKSideLCMS = py1L - py2L;
fDKLongLCMS = pz1L - pz2L;
}
fDKOutLCMS = px1LCMS - px2LCMS;
fDKSideLCMS = py1LCMS - py2LCMS;
fDKLongLCMS = pz1LCMS - pz2LCMS;

// Boost to PRF

const double betaOut = tPt / tMt;
const double gammaOut = tMt / tPinv;
const double gammaOut = tMt / tM;

if (!isiden) {
fDKOutPRF = gammaOut * (fDKOutLCMS - betaOut * pE1L);
fDKSidePRF = fDKSideLCMS;
fDKLongPRF = fDKLongLCMS;
} else {
fDKOutPRF = gammaOut * (fDKOutLCMS - betaOut * (pE1L - pE2L));
fDKSidePRF = fDKSideLCMS;
fDKLongPRF = fDKLongLCMS;
}
fDKOut = gammaOut * (fDKOut - betaOut * fDE);
kstar = sqrt(fDKOut * fDKOut + fDKSide * fDKSide + fDKLong * fDKLong);

if (islcms) {
if (isiden) {
vect.push_back(2.0 * (kstar));
vect.push_back(fDKOutLCMS);
vect.push_back(fDKSideLCMS);
vect.push_back(fDKLongLCMS);
} else {
vect.push_back(fDKOutPRF);
vect.push_back(fDKSidePRF);
vect.push_back(fDKLongPRF);
vect.push_back(kstar);
vect.push_back(fDKOut);
vect.push_back(fDKSide);
vect.push_back(fDKLong);
}
return vect;
}
Expand Down
14 changes: 7 additions & 7 deletions PWGCF/FemtoUniverse/Core/FemtoUniversePairSHCentMultKt.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -184,7 +184,7 @@ class PairSHCentMultKt
/// \param ktval kT value
template <typename T>
void fill_mult_NumDen(T const& part1, T const& part2, uint8_t ChosenEventType,
int maxl, int multval, float ktval)
int maxl, int multval, float ktval, bool isiden)
{
int multbinval;
int absmultval = multval;
Expand All @@ -201,7 +201,7 @@ class PairSHCentMultKt
return;
}
// std::cout<<"multbinval "<<multbinval<<std::endl;
fill_kT_NumDen(part1, part2, ChosenEventType, maxl, multbinval, ktval);
fill_kT_NumDen(part1, part2, ChosenEventType, maxl, multbinval, ktval, isiden);
}

/// Templated function to access different kT directory and call AddEventPair
Expand All @@ -213,7 +213,7 @@ class PairSHCentMultKt
/// \param ktval kT value
template <typename T>
void fill_kT_NumDen(T const& part1, T const& part2, uint8_t ChosenEventType,
int maxl, int multval, float ktval)
int maxl, int multval, float ktval, bool isiden)
{
int ktbinval = -1;
if ((ktval >= KtBins[0]) && (ktval < KtBins[1])) {
Expand All @@ -227,7 +227,7 @@ class PairSHCentMultKt
} else {
return;
}
AddEventPair(part1, part2, ChosenEventType, maxl, multval, ktbinval);
AddEventPair(part1, part2, ChosenEventType, maxl, multval, ktbinval, isiden);
}

/// Set the PDG codes of the two particles involved
Expand Down Expand Up @@ -262,14 +262,14 @@ class PairSHCentMultKt
/// \param ktval kT value
template <typename T>
void AddEventPair(T const& part1, T const& part2, uint8_t ChosenEventType,
int /*maxl*/, int multval, int ktval)
int /*maxl*/, int multval, int ktval, bool isiden)
{
int fMultBin = multval;
int fKtBin = ktval;
std::vector<std::complex<double>> fYlmBuffer(fMaxJM);
std::vector<double> f3d;
f3d = FemtoUniverseMath::getpairmom3d(part1, mMassOne, part2, mMassTwo,
true, true);
f3d = FemtoUniverseMath::newpairfunc(part1, mMassOne, part2, mMassTwo,
isiden);

const float qout = f3d[1];
const float qside = f3d[2];
Expand Down
9 changes: 3 additions & 6 deletions PWGCF/FemtoUniverse/Core/FemtoUniverseSHContainer.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -159,20 +159,17 @@ class FemtoUniverseSHContainer
/// \param ChosenEventType same or mixed event
/// \param maxl Maximum valie of L component of the spherical harmonics
template <bool isMC, typename T>
void AddEventPair(T const& part1, T const& part2, uint8_t ChosenEventType, int /*maxl*/)
void AddEventPair(T const& part1, T const& part2, uint8_t ChosenEventType, int /*maxl*/, bool isiden)
{
// int fMaxL = 2;
// int fMaxJM = (2+1)*(2+1);
std::vector<std::complex<double>> fYlmBuffer(fMaxJM);
std::vector<double> f3d;
f3d = FemtoUniverseMath::getpairmom3d(part1, mMassOne, part2, mMassTwo, true, true);
f3d = FemtoUniverseMath::newpairfunc(part1, mMassOne, part2, mMassTwo, isiden);

// const float qstar = f3d[0];
const float kv = f3d[0];
const float qout = f3d[1];
const float qside = f3d[2];
const float qlong = f3d[3];

double kv = sqrt(qout * qout + qside * qside + qlong * qlong);
int nqbin = fbinctn->GetXaxis()->FindFixBin(kv) - 1;

FemtoUniverseSpherHarMath Ylm;
Expand Down
32 changes: 16 additions & 16 deletions PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrack3DMultKtExtended.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -476,9 +476,9 @@ struct femtoUniversePairTaskTrackTrack3DMultKtExtended {
float kT = FemtoUniverseMath::getkT(p1, mass1, p2, mass2);

if (!cfgProcessMultBins) {
sameEventCont.setPair<isMC>(p1, p2, multCol, twotracksconfigs.ConfUse3D, ConfIsIden, ConfIsLCMS);
sameEventCont.setPair<isMC>(p1, p2, multCol, twotracksconfigs.ConfUse3D, ConfIsIden);
} else {
std::vector<double> k3d = FemtoUniverseMath::getpairmom3d(p1, mass1, p2, mass2, ConfIsIden, ConfIsLCMS);
std::vector<double> k3d = FemtoUniverseMath::newpairfunc(p1, mass1, p2, mass2, ConfIsIden);
sameEventMultCont.fill_3D<float>(k3d[1], k3d[2], k3d[3], multCol, kT);
}
}
Expand Down Expand Up @@ -515,16 +515,16 @@ struct femtoUniversePairTaskTrackTrack3DMultKtExtended {

if (rand > 0.5) {
if (!cfgProcessMultBins) {
sameEventContPP.setPair<isMC>(p1, p2, multCol, twotracksconfigs.ConfUse3D, ConfIsIden, ConfIsLCMS);
sameEventContPP.setPair<isMC>(p1, p2, multCol, twotracksconfigs.ConfUse3D, ConfIsIden);
} else {
k3d = FemtoUniverseMath::getpairmom3d(p1, mass1, p2, mass2, ConfIsIden, ConfIsLCMS);
k3d = FemtoUniverseMath::newpairfunc(p1, mass1, p2, mass2, ConfIsIden);
sameEventMultContPP.fill_3D<float>(k3d[1], k3d[2], k3d[3], multCol, kT);
}
} else {
if (!cfgProcessMultBins) {
sameEventContPP.setPair<isMC>(p2, p1, multCol, twotracksconfigs.ConfUse3D, ConfIsIden, ConfIsLCMS);
sameEventContPP.setPair<isMC>(p2, p1, multCol, twotracksconfigs.ConfUse3D, ConfIsIden);
} else {
k3d = FemtoUniverseMath::getpairmom3d(p2, mass2, p1, mass1, ConfIsIden, ConfIsLCMS);
k3d = FemtoUniverseMath::newpairfunc(p2, mass2, p1, mass1, ConfIsIden);
sameEventMultContPP.fill_3D<float>(k3d[1], k3d[2], k3d[3], multCol, kT);
}
}
Expand All @@ -538,16 +538,16 @@ struct femtoUniversePairTaskTrackTrack3DMultKtExtended {

if (rand > 0.5) {
if (!cfgProcessMultBins) {
sameEventContMM.setPair<isMC>(p1, p2, multCol, twotracksconfigs.ConfUse3D, ConfIsIden, ConfIsLCMS);
sameEventContMM.setPair<isMC>(p1, p2, multCol, twotracksconfigs.ConfUse3D, ConfIsIden);
} else {
k3d = FemtoUniverseMath::getpairmom3d(p1, mass1, p2, mass2, ConfIsIden, ConfIsLCMS);
k3d = FemtoUniverseMath::newpairfunc(p1, mass1, p2, mass2, ConfIsIden);
sameEventMultContMM.fill_3D<float>(k3d[1], k3d[2], k3d[3], multCol, kT);
}
} else {
if (!cfgProcessMultBins) {
sameEventContMM.setPair<isMC>(p2, p1, multCol, twotracksconfigs.ConfUse3D, ConfIsIden, ConfIsLCMS);
sameEventContMM.setPair<isMC>(p2, p1, multCol, twotracksconfigs.ConfUse3D, ConfIsIden);
} else {
k3d = FemtoUniverseMath::getpairmom3d(p2, mass2, p1, mass1, ConfIsIden, ConfIsLCMS);
k3d = FemtoUniverseMath::newpairfunc(p2, mass2, p1, mass1, ConfIsIden);
sameEventMultContMM.fill_3D<float>(k3d[1], k3d[2], k3d[3], multCol, kT);
}
}
Expand Down Expand Up @@ -650,19 +650,19 @@ struct femtoUniversePairTaskTrackTrack3DMultKtExtended {
case 1: {
float kT = FemtoUniverseMath::getkT(p1, mass1, p2, mass2);
if (!cfgProcessMultBins) {
mixedEventCont.setPair<isMC>(p1, p2, multCol, twotracksconfigs.ConfUse3D, ConfIsIden, ConfIsLCMS);
mixedEventCont.setPair<isMC>(p1, p2, multCol, twotracksconfigs.ConfUse3D, ConfIsIden);
} else {
std::vector<double> k3d = FemtoUniverseMath::getpairmom3d(p1, mass1, p2, mass2, ConfIsIden, ConfIsLCMS);
std::vector<double> k3d = FemtoUniverseMath::newpairfunc(p1, mass1, p2, mass2, ConfIsIden);
mixedEventMultCont.fill_3D<float>(k3d[1], k3d[2], k3d[3], multCol, kT);
}
break;
}
case 2: {
float kT = FemtoUniverseMath::getkT(p1, mass1, p2, mass1);
if (!cfgProcessMultBins) {
mixedEventContPP.setPair<isMC>(p1, p2, multCol, twotracksconfigs.ConfUse3D, ConfIsIden, ConfIsLCMS);
mixedEventContPP.setPair<isMC>(p1, p2, multCol, twotracksconfigs.ConfUse3D, ConfIsIden);
} else {
std::vector<double> k3d = FemtoUniverseMath::getpairmom3d(p1, mass1, p2, mass2, ConfIsIden, ConfIsLCMS);
std::vector<double> k3d = FemtoUniverseMath::newpairfunc(p1, mass1, p2, mass2, ConfIsIden);
mixedEventMultContPP.fill_3D<float>(k3d[1], k3d[2], k3d[3], multCol, kT);
}
break;
Expand All @@ -671,9 +671,9 @@ struct femtoUniversePairTaskTrackTrack3DMultKtExtended {
case 3: {
float kT = FemtoUniverseMath::getkT(p1, mass2, p2, mass2);
if (!cfgProcessMultBins) {
mixedEventContMM.setPair<isMC>(p1, p2, multCol, twotracksconfigs.ConfUse3D, ConfIsIden, ConfIsLCMS);
mixedEventContMM.setPair<isMC>(p1, p2, multCol, twotracksconfigs.ConfUse3D, ConfIsIden);
} else {
std::vector<double> k3d = FemtoUniverseMath::getpairmom3d(p1, mass1, p2, mass2, ConfIsIden, ConfIsLCMS);
std::vector<double> k3d = FemtoUniverseMath::newpairfunc(p1, mass1, p2, mass2, ConfIsIden);
mixedEventMultContMM.fill_3D<float>(k3d[1], k3d[2], k3d[3], multCol, kT);
}
break;
Expand Down
16 changes: 8 additions & 8 deletions PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackSpherHarMultKtExtended.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -475,7 +475,7 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended {
continue;
}
float kT = FemtoUniverseMath::getkT(p1, mass1, p2, mass2);
sameEventMultCont.fill_mult_NumDen(p1, p2, femtoUniverseSHContainer::EventType::same, 2, multCol, kT);
sameEventMultCont.fill_mult_NumDen(p1, p2, femtoUniverseSHContainer::EventType::same, 2, multCol, kT, ConfIsIden);
}
} else {
for (auto& [p1, p2] : combinations(CombinationsStrictlyUpperIndexPolicy(groupPartsOne, groupPartsOne))) {
Expand Down Expand Up @@ -507,19 +507,19 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended {
case 2: {
rand = randgen->Rndm();
if (rand > 0.5) {
sameEventMultContPP.fill_mult_NumDen(p1, p2, femtoUniverseSHContainer::EventType::same, 2, multCol, kT);
sameEventMultContPP.fill_mult_NumDen(p1, p2, femtoUniverseSHContainer::EventType::same, 2, multCol, kT, ConfIsIden);
} else if (rand <= 0.5) {
sameEventMultContPP.fill_mult_NumDen(p2, p1, femtoUniverseSHContainer::EventType::same, 2, multCol, kT);
sameEventMultContPP.fill_mult_NumDen(p2, p1, femtoUniverseSHContainer::EventType::same, 2, multCol, kT, ConfIsIden);
}
break;
}

case 3: {
rand = randgen->Rndm();
if (rand > 0.5) {
sameEventMultContMM.fill_mult_NumDen(p1, p2, femtoUniverseSHContainer::EventType::same, 2, multCol, kT);
sameEventMultContMM.fill_mult_NumDen(p1, p2, femtoUniverseSHContainer::EventType::same, 2, multCol, kT, ConfIsIden);
} else if (rand <= 0.5) {
sameEventMultContMM.fill_mult_NumDen(p2, p1, femtoUniverseSHContainer::EventType::same, 2, multCol, kT);
sameEventMultContMM.fill_mult_NumDen(p2, p1, femtoUniverseSHContainer::EventType::same, 2, multCol, kT, ConfIsIden);
}
break;
}
Expand Down Expand Up @@ -606,15 +606,15 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended {
float kT = FemtoUniverseMath::getkT(p1, mass1, p2, mass2);
switch (ContType) {
case 1: {
mixedEventMultCont.fill_mult_NumDen(p1, p2, femtoUniverseSHContainer::EventType::mixed, 2, multCol, kT);
mixedEventMultCont.fill_mult_NumDen(p1, p2, femtoUniverseSHContainer::EventType::mixed, 2, multCol, kT, ConfIsIden);
break;
}
case 2: {
mixedEventMultContPP.fill_mult_NumDen(p1, p2, femtoUniverseSHContainer::EventType::mixed, 2, multCol, kT);
mixedEventMultContPP.fill_mult_NumDen(p1, p2, femtoUniverseSHContainer::EventType::mixed, 2, multCol, kT, ConfIsIden);
break;
}
case 3: {
mixedEventMultContMM.fill_mult_NumDen(p1, p2, femtoUniverseSHContainer::EventType::mixed, 2, multCol, kT);
mixedEventMultContMM.fill_mult_NumDen(p1, p2, femtoUniverseSHContainer::EventType::mixed, 2, multCol, kT, ConfIsIden);
break;
}
default:
Expand Down