Skip to content

Commit 57b4bff

Browse files
sawankumawatSawan Sawan
sawankumawat
and
Sawan Sawan
authored
[PWGLF] Added Collins Soper frame and renamed some configurables (#13005)
Co-authored-by: Sawan Sawan <sawan.sawan@cern.ch>
1 parent dfc407d commit 57b4bff

File tree

1 file changed

+68
-25
lines changed

1 file changed

+68
-25
lines changed

‎PWGLF/Tasks/Resonances/higherMassResonances.cxx‎

Lines changed: 68 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -154,10 +154,11 @@ struct HigherMassResonances {
154154
std::vector<int> pdgCodes = {10331, 335, 115, 10221, 9030221};
155155

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

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

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

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

219220
if (std::accumulate(sparses.begin(), sparses.end(), 0) == 0) {
220221
LOGP(fatal, "No output THnSparses enabled");
221222
} else {
222-
if (config.activateTHnSparseCosThStarHelicity) {
223+
if (config.activateHelicityFrame) {
223224
LOGP(info, "THnSparse with cosThStar w.r.t. helicity axis active.");
224225
}
225-
if (config.activateTHnSparseCosThStarProduction) {
226+
if (config.activateCollinsSoperFrame) {
227+
LOGP(info, "THnSparse with cosThStar w.r.t. Collins Soper axis active.");
228+
}
229+
if (config.activateProductionFrame) {
226230
LOGP(info, "THnSparse with cosThStar w.r.t. production axis active.");
227231
}
228-
if (config.activateTHnSparseCosThStarBeam) {
232+
if (config.activateBeamAxisFrame) {
229233
LOGP(info, "THnSparse with cosThStar w.r.t. beam axis active. (Gottified jackson frame)");
230234
}
231-
if (config.activateTHnSparseCosThStarRandom) {
235+
if (config.activateRandomFrame) {
232236
LOGP(info, "THnSparse with cosThStar w.r.t. random axis active.");
233237
}
234238
}
@@ -679,6 +683,7 @@ struct HigherMassResonances {
679683
beam1CM = ROOT::Math::XYZVectorF((boost(beam1).Vect()).Unit());
680684
beam2CM = ROOT::Math::XYZVectorF((boost(beam2).Vect()).Unit());
681685

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

@@ -696,12 +701,13 @@ struct HigherMassResonances {
696701

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

700706
v1CM = ROOT::Math::XYZVectorF(boost(daughter1).Vect()).Unit();
701707
// ROOT::Math::XYZVectorF v2_CM{(boost(daughter1).Vect()).Unit()};
702708
// using positive sign convention for the first track
703709
// 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
704-
// Helicity frame
710+
// Helicity Frame
705711
zaxisHE = ROOT::Math::XYZVectorF(mother.Vect()).Unit();
706712
yaxisHE = ROOT::Math::XYZVectorF(beam1CM.Cross(beam2CM)).Unit();
707713
xaxisHE = ROOT::Math::XYZVectorF(yaxisHE.Cross(zaxisHE)).Unit();
@@ -714,8 +720,16 @@ struct HigherMassResonances {
714720
// anglePhi += o2::constants::math::TwoPI; // ensure phi is in [0, 2pi]
715721
// }
716722

723+
// CS Frame
724+
zAxisCS = ROOT::Math::XYZVectorF((beam1CM.Unit() - beam2CM.Unit())).Unit();
725+
yAxisCS = ROOT::Math::XYZVectorF(beam1CM.Cross(beam2CM)).Unit();
726+
xAxisCS = ROOT::Math::XYZVectorF(yAxisCS.Cross(zAxisCS)).Unit();
727+
double cosThetaStarCS = zAxisCS.Dot(v1CM);
728+
auto phiCS = std::atan2(yAxisCS.Dot(v1CM), xAxisCS.Dot(v1CM));
729+
phiCS = RecoDecay::constrainAngle(phiCS, 0.0);
730+
717731
// if (std::abs(mother.Rapidity()) < config.rapidityMotherData) {
718-
if (config.activateTHnSparseCosThStarHelicity) {
732+
if (config.activateHelicityFrame) {
719733
// helicityVec = mother.Vect(); // 3 vector of mother in COM frame
720734
// auto cosThetaStarHelicity = helicityVec.Dot(threeVecDauCM) / (std::sqrt(threeVecDauCM.Mag2()) * std::sqrt(helicityVec.Mag2()));
721735
auto cosThetaStarHelicity = mother.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(mother.Vect().Mag2()));
@@ -735,34 +749,64 @@ struct HigherMassResonances {
735749
daughterRotCM = boost2(daughterRot);
736750

737751
auto cosThetaStarHelicityRot = motherRot.Vect().Dot(daughterRotCM.Vect()) / (std::sqrt(daughterRotCM.Vect().Mag2()) * std::sqrt(motherRot.Vect().Mag2()));
752+
auto phiHelicityRot = std::atan2(yaxisHE.Dot(daughterRotCM.Vect().Unit()), xaxisHE.Dot(daughterRotCM.Vect().Unit()));
753+
phiHelicityRot = RecoDecay::constrainAngle(phiHelicityRot, 0.0);
738754
if (motherRot.Rapidity() < config.rapidityMotherData)
739-
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarHelicityRot, anglePhi);
755+
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarHelicityRot, phiHelicityRot);
740756
}
741757
} else {
742758
if (std::abs(mother.Rapidity()) < config.rapidityMotherData) {
743759
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity, anglePhi);
744760
}
745761
}
746-
} else if (config.activateTHnSparseCosThStarProduction) {
762+
} else if (config.activateCollinsSoperFrame) {
763+
if (!isMix) {
764+
if (std::abs(mother.Rapidity()) < config.rapidityMotherData) {
765+
hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarCS, phiCS);
766+
}
767+
768+
for (int i = 0; i < config.cRotations; i++) {
769+
theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut);
770+
771+
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());
772+
773+
motherRot = daughterRot + daughter2;
774+
775+
ROOT::Math::Boost boost2{motherRot.BoostToCM()};
776+
daughterRotCM = boost2(daughterRot);
777+
778+
auto cosThetaStarCSrot = zAxisCS.Dot(daughterRotCM.Vect()) / std::sqrt(daughterRotCM.Vect().Mag2());
779+
auto phiCSrot = std::atan2(yAxisCS.Dot(daughterRotCM.Vect().Unit()), xAxisCS.Dot(daughterRotCM.Vect().Unit()));
780+
phiCSrot = RecoDecay::constrainAngle(phiCSrot, 0.0);
781+
782+
if (motherRot.Rapidity() < config.rapidityMotherData)
783+
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarCSrot, phiCSrot);
784+
}
785+
} else {
786+
if (std::abs(mother.Rapidity()) < config.rapidityMotherData) {
787+
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarCS, phiCS);
788+
}
789+
}
790+
} else if (config.activateProductionFrame) {
747791
normalVec = ROOT::Math::XYZVector(mother.Py(), -mother.Px(), 0.f);
748-
auto cosThetaStarProduction = normalVec.Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(normalVec.Mag2()));
792+
auto cosThetaProduction = normalVec.Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(normalVec.Mag2()));
749793
if (!isMix) {
750794
if (std::abs(mother.Rapidity()) < config.rapidityMotherData) {
751-
hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarProduction, anglePhi);
795+
hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaProduction, anglePhi);
752796
}
753797
for (int i = 0; i < config.cRotations; i++) {
754798
theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut);
755799
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());
756800
if (std::abs(motherRot.Rapidity()) < config.rapidityMotherData) {
757-
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarProduction, anglePhi);
801+
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaProduction, anglePhi);
758802
}
759803
}
760804
} else {
761805
if (std::abs(mother.Rapidity()) < config.rapidityMotherData) {
762-
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarProduction, anglePhi);
806+
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaProduction, anglePhi);
763807
}
764808
}
765-
} else if (config.activateTHnSparseCosThStarBeam) {
809+
} else if (config.activateBeamAxisFrame) {
766810
beamVec = ROOT::Math::XYZVector(0.f, 0.f, 1.f);
767811
auto cosThetaStarBeam = beamVec.Dot(fourVecDauCM.Vect()) / std::sqrt(fourVecDauCM.Vect().Mag2());
768812
if (!isMix) {
@@ -781,26 +825,26 @@ struct HigherMassResonances {
781825
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarBeam, anglePhi);
782826
}
783827
}
784-
} else if (config.activateTHnSparseCosThStarRandom) {
828+
} else if (config.activateRandomFrame) {
785829
auto phiRandom = gRandom->Uniform(0.f, constants::math::TwoPI);
786830
auto thetaRandom = gRandom->Uniform(0.f, constants::math::PI);
787831

788832
randomVec = ROOT::Math::XYZVector(std::sin(thetaRandom) * std::cos(phiRandom), std::sin(thetaRandom) * std::sin(phiRandom), std::cos(thetaRandom));
789833
auto cosThetaStarRandom = randomVec.Dot(fourVecDauCM.Vect()) / std::sqrt(fourVecDauCM.Vect().Mag2());
790834
if (!isMix) {
791835
if (std::abs(mother.Rapidity()) < config.rapidityMotherData) {
792-
hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, anglePhi);
836+
hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, phiRandom);
793837
}
794838
for (int i = 0; i < config.cRotations; i++) {
795839
theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut);
796840
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());
797841
if (std::abs(motherRot.Rapidity()) < config.rapidityMotherData) {
798-
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarRandom, anglePhi);
842+
hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarRandom, phiRandom);
799843
}
800844
}
801845
} else {
802846
if (std::abs(mother.Rapidity()) < config.rapidityMotherData) {
803-
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, anglePhi);
847+
hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, phiRandom);
804848
}
805849
}
806850
}
@@ -1444,6 +1488,5 @@ struct HigherMassResonances {
14441488

14451489
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
14461490
{
1447-
return WorkflowSpec{
1448-
adaptAnalysisTask<HigherMassResonances>(cfgc)};
1491+
return WorkflowSpec{adaptAnalysisTask<HigherMassResonances>(cfgc)};
14491492
}

0 commit comments

Comments
 (0)