Skip to content

[PWGLF] Added Collins Soper frame and renamed some configurables #13005

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
dmallick2 merged 2 commits into from
Sep 17, 2025
Merged

[PWGLF] Added Collins Soper frame and renamed some configurables #13005

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d26158f
added Collins Soper frame and renamed some configurables
Sep 17, 2025
03fd420
fixed clang errors
Sep 17, 2025
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
93 changes: 68 additions & 25 deletions PWGLF/Tasks/Resonances/higherMassResonances.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -154,10 +154,11 @@ struct HigherMassResonances {
std::vector<int> pdgCodes = {10331, 335, 115, 10221, 9030221};

// output THnSparses
Configurable<bool> activateTHnSparseCosThStarHelicity{"activateTHnSparseCosThStarHelicity", false, "Activate the THnSparse with cosThStar w.r.t. helicity axis"};
Configurable<bool> activateTHnSparseCosThStarProduction{"activateTHnSparseCosThStarProduction", false, "Activate the THnSparse with cosThStar w.r.t. production axis"};
Configurable<bool> activateTHnSparseCosThStarBeam{"activateTHnSparseCosThStarBeam", true, "Activate the THnSparse with cosThStar w.r.t. beam axis (Gottified jackson frame)"};
Configurable<bool> activateTHnSparseCosThStarRandom{"activateTHnSparseCosThStarRandom", false, "Activate the THnSparse with cosThStar w.r.t. random axis"};
Configurable<bool> activateHelicityFrame{"activateHelicityFrame", false, "Activate the THnSparse with cosThStar w.r.t. helicity axis"};
Configurable<bool> activateCollinsSoperFrame{"activateCollinsSoperFrame", false, "Activate the THnSparse with cosThStar w.r.t. Collins soper axis"};
Configurable<bool> activateProductionFrame{"activateProductionFrame", false, "Activate the THnSparse with cosThStar w.r.t. production axis"};
Configurable<bool> activateBeamAxisFrame{"activateBeamAxisFrame", true, "Activate the THnSparse with cosThStar w.r.t. beam axis (Gottified jackson frame)"};
Configurable<bool> activateRandomFrame{"activateRandomFrame", false, "Activate the THnSparse with cosThStar w.r.t. random axis"};
Configurable<int> cRotations{"cRotations", 3, "Number of random rotations in the rotational background"};

// Other cuts on Ks and glueball
Expand Down Expand Up @@ -195,7 +196,7 @@ struct HigherMassResonances {
ROOT::Math::XYZVector zBeam; // ẑ: beam direction in lab frame
ROOT::Math::PxPyPzEVector beam1{0., 0., -config.beamMomentum, 13600. / 2.};
ROOT::Math::PxPyPzEVector beam2{0., 0., config.beamMomentum, 13600. / 2.};
ROOT::Math::XYZVectorF beam1CM, beam2CM;
ROOT::Math::XYZVectorF beam1CM, beam2CM, zAxisCS, yAxisCS, xAxisCS;

// const double massK0s = o2::constants::physics::MassK0Short;
bool isMix = false;
Expand All @@ -214,21 +215,24 @@ struct HigherMassResonances {
AxisSpec thnAxisPhi = {config.configThnAxisPhi, "Configurabel phi axis"}; // 0 to 2pi

// THnSparses
std::array<int, 4> sparses = {config.activateTHnSparseCosThStarHelicity, config.activateTHnSparseCosThStarProduction, config.activateTHnSparseCosThStarBeam, config.activateTHnSparseCosThStarRandom};
std::array<int, 5> sparses = {config.activateHelicityFrame, config.activateCollinsSoperFrame, config.activateProductionFrame, config.activateBeamAxisFrame, config.activateRandomFrame};

if (std::accumulate(sparses.begin(), sparses.end(), 0) == 0) {
LOGP(fatal, "No output THnSparses enabled");
} else {
if (config.activateTHnSparseCosThStarHelicity) {
if (config.activateHelicityFrame) {
LOGP(info, "THnSparse with cosThStar w.r.t. helicity axis active.");
}
if (config.activateTHnSparseCosThStarProduction) {
if (config.activateCollinsSoperFrame) {
LOGP(info, "THnSparse with cosThStar w.r.t. Collins Soper axis active.");
}
if (config.activateProductionFrame) {
LOGP(info, "THnSparse with cosThStar w.r.t. production axis active.");
}
if (config.activateTHnSparseCosThStarBeam) {
if (config.activateBeamAxisFrame) {
LOGP(info, "THnSparse with cosThStar w.r.t. beam axis active. (Gottified jackson frame)");
}
if (config.activateTHnSparseCosThStarRandom) {
if (config.activateRandomFrame) {
LOGP(info, "THnSparse with cosThStar w.r.t. random axis active.");
}
}
Expand Down Expand Up @@ -679,6 +683,7 @@ struct HigherMassResonances {
beam1CM = ROOT::Math::XYZVectorF((boost(beam1).Vect()).Unit());
beam2CM = ROOT::Math::XYZVectorF((boost(beam2).Vect()).Unit());

//========================Helicity and Production frame calculation==========================
// define y = zBeam x z: Normal to the production plane
// ẑ: mother direction in lab, boosted into mother's rest frame

Expand All @@ -696,12 +701,13 @@ struct HigherMassResonances {

// // Calculate φ in [-π, π]
// auto anglePhi = std::atan2(p_proj_y, p_proj_x); // φ in radians
//=============================================================================================

v1CM = ROOT::Math::XYZVectorF(boost(daughter1).Vect()).Unit();
// ROOT::Math::XYZVectorF v2_CM{(boost(daughter1).Vect()).Unit()};
// using positive sign convention for the first track
// ROOT::Math::XYZVectorF v_CM = (t1.sign() > 0 ? v1CM : v2_CM); // here selected decay daughter momentum is intested. here you can choose one decay daughter no need to check both case as it is neutral particle for our case
// Helicity frame
// Helicity Frame
zaxisHE = ROOT::Math::XYZVectorF(mother.Vect()).Unit();
yaxisHE = ROOT::Math::XYZVectorF(beam1CM.Cross(beam2CM)).Unit();
xaxisHE = ROOT::Math::XYZVectorF(yaxisHE.Cross(zaxisHE)).Unit();
Expand All @@ -714,8 +720,16 @@ struct HigherMassResonances {
// anglePhi += o2::constants::math::TwoPI; // ensure phi is in [0, 2pi]
// }

// CS Frame
zAxisCS = ROOT::Math::XYZVectorF((beam1CM.Unit() - beam2CM.Unit())).Unit();
yAxisCS = ROOT::Math::XYZVectorF(beam1CM.Cross(beam2CM)).Unit();
xAxisCS = ROOT::Math::XYZVectorF(yAxisCS.Cross(zAxisCS)).Unit();
double cosThetaStarCS = zAxisCS.Dot(v1CM);
auto phiCS = std::atan2(yAxisCS.Dot(v1CM), xAxisCS.Dot(v1CM));
phiCS = RecoDecay::constrainAngle(phiCS, 0.0);

// if (std::abs(mother.Rapidity()) < config.rapidityMotherData) {
if (config.activateTHnSparseCosThStarHelicity) {
if (config.activateHelicityFrame) {
// helicityVec = mother.Vect(); // 3 vector of mother in COM frame
// auto cosThetaStarHelicity = helicityVec.Dot(threeVecDauCM) / (std::sqrt(threeVecDauCM.Mag2()) * std::sqrt(helicityVec.Mag2()));
auto cosThetaStarHelicity = mother.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(mother.Vect().Mag2()));
Expand All @@ -735,34 +749,64 @@ struct HigherMassResonances {
daughterRotCM = boost2(daughterRot);

auto cosThetaStarHelicityRot = motherRot.Vect().Dot(daughterRotCM.Vect()) / (std::sqrt(daughterRotCM.Vect().Mag2()) * std::sqrt(motherRot.Vect().Mag2()));
auto phiHelicityRot = std::atan2(yaxisHE.Dot(daughterRotCM.Vect().Unit()), xaxisHE.Dot(daughterRotCM.Vect().Unit()));
phiHelicityRot = RecoDecay::constrainAngle(phiHelicityRot, 0.0);
if (motherRot.Rapidity() < config.rapidityMotherData)
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarHelicityRot, anglePhi);
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarHelicityRot, phiHelicityRot);
}
} else {
if (std::abs(mother.Rapidity()) < config.rapidityMotherData) {
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity, anglePhi);
}
}
} else if (config.activateTHnSparseCosThStarProduction) {
} else if (config.activateCollinsSoperFrame) {
if (!isMix) {
if (std::abs(mother.Rapidity()) < config.rapidityMotherData) {
hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarCS, phiCS);
}

for (int i = 0; i < config.cRotations; i++) {
theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut);

daughterRot = ROOT::Math::PxPyPzMVector(daughter1.Px() * std::cos(theta2) - daughter1.Py() * std::sin(theta2), daughter1.Px() * std::sin(theta2) + daughter1.Py() * std::cos(theta2), daughter1.Pz(), daughter1.M());

motherRot = daughterRot + daughter2;

ROOT::Math::Boost boost2{motherRot.BoostToCM()};
daughterRotCM = boost2(daughterRot);

auto cosThetaStarCSrot = zAxisCS.Dot(daughterRotCM.Vect()) / std::sqrt(daughterRotCM.Vect().Mag2());
auto phiCSrot = std::atan2(yAxisCS.Dot(daughterRotCM.Vect().Unit()), xAxisCS.Dot(daughterRotCM.Vect().Unit()));
phiCSrot = RecoDecay::constrainAngle(phiCSrot, 0.0);

if (motherRot.Rapidity() < config.rapidityMotherData)
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarCSrot, phiCSrot);
}
} else {
if (std::abs(mother.Rapidity()) < config.rapidityMotherData) {
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarCS, phiCS);
}
}
} else if (config.activateProductionFrame) {
normalVec = ROOT::Math::XYZVector(mother.Py(), -mother.Px(), 0.f);
auto cosThetaStarProduction = normalVec.Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(normalVec.Mag2()));
auto cosThetaProduction = normalVec.Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(normalVec.Mag2()));
if (!isMix) {
if (std::abs(mother.Rapidity()) < config.rapidityMotherData) {
hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarProduction, anglePhi);
hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaProduction, anglePhi);
}
for (int i = 0; i < config.cRotations; i++) {
theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut);
motherRot = ROOT::Math::PxPyPzMVector(mother.Px() * std::cos(theta2) - mother.Py() * std::sin(theta2), mother.Px() * std::sin(theta2) + mother.Py() * std::cos(theta2), mother.Pz(), mother.M());
if (std::abs(motherRot.Rapidity()) < config.rapidityMotherData) {
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarProduction, anglePhi);
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaProduction, anglePhi);
}
}
} else {
if (std::abs(mother.Rapidity()) < config.rapidityMotherData) {
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarProduction, anglePhi);
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaProduction, anglePhi);
}
}
} else if (config.activateTHnSparseCosThStarBeam) {
} else if (config.activateBeamAxisFrame) {
beamVec = ROOT::Math::XYZVector(0.f, 0.f, 1.f);
auto cosThetaStarBeam = beamVec.Dot(fourVecDauCM.Vect()) / std::sqrt(fourVecDauCM.Vect().Mag2());
if (!isMix) {
Expand All @@ -781,26 +825,26 @@ struct HigherMassResonances {
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarBeam, anglePhi);
}
}
} else if (config.activateTHnSparseCosThStarRandom) {
} else if (config.activateRandomFrame) {
auto phiRandom = gRandom->Uniform(0.f, constants::math::TwoPI);
auto thetaRandom = gRandom->Uniform(0.f, constants::math::PI);

randomVec = ROOT::Math::XYZVector(std::sin(thetaRandom) * std::cos(phiRandom), std::sin(thetaRandom) * std::sin(phiRandom), std::cos(thetaRandom));
auto cosThetaStarRandom = randomVec.Dot(fourVecDauCM.Vect()) / std::sqrt(fourVecDauCM.Vect().Mag2());
if (!isMix) {
if (std::abs(mother.Rapidity()) < config.rapidityMotherData) {
hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, anglePhi);
hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, phiRandom);
}
for (int i = 0; i < config.cRotations; i++) {
theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut);
motherRot = ROOT::Math::PxPyPzMVector(mother.Px() * std::cos(theta2) - mother.Py() * std::sin(theta2), mother.Px() * std::sin(theta2) + mother.Py() * std::cos(theta2), mother.Pz(), mother.M());
if (std::abs(motherRot.Rapidity()) < config.rapidityMotherData) {
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarRandom, anglePhi);
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarRandom, phiRandom);
}
}
} else {
if (std::abs(mother.Rapidity()) < config.rapidityMotherData) {
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, anglePhi);
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, phiRandom);
}
}
}
Expand Down Expand Up @@ -1444,6 +1488,5 @@ struct HigherMassResonances {

WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
{
return WorkflowSpec{
adaptAnalysisTask<HigherMassResonances>(cfgc)};
return WorkflowSpec{adaptAnalysisTask<HigherMassResonances>(cfgc)};
}
Loading