Skip to content

Commit 23920fa

Browse files
sawankumawatsawankumawat
committed
added correlated histogram
1 parent 225de2c commit 23920fa

1 file changed

Lines changed: 95 additions & 21 deletions

File tree

‎PWGLF/Tasks/Resonances/higherMassResonances.cxx‎

Lines changed: 95 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -174,7 +174,7 @@ struct HigherMassResonances {
174174
Configurable<bool> activateHelicityFrame{"activateHelicityFrame", false, "Activate the THnSparse with cosThStar w.r.t. helicity axis"};
175175
Configurable<bool> activateCollinsSoperFrame{"activateCollinsSoperFrame", false, "Activate the THnSparse with cosThStar w.r.t. Collins soper axis"};
176176
Configurable<bool> activateProductionFrame{"activateProductionFrame", false, "Activate the THnSparse with cosThStar w.r.t. production axis"};
177-
Configurable<bool> activateBeamAxisFrame{"activateBeamAxisFrame", true, "Activate the THnSparse with cosThStar w.r.t. beam axis (Gottified jackson frame)"};
177+
Configurable<bool> activateGJFrame{"activateGJFrame", true, "Activate the THnSparse with cosThStar w.r.t. beam axis (Gottified jackson frame)"};
178178
Configurable<bool> activateRandomFrame{"activateRandomFrame", false, "Activate the THnSparse with cosThStar w.r.t. random axis"};
179179
Configurable<int> cRotations{"cRotations", 3, "Number of random rotations in the rotational background"};
180180

@@ -255,7 +255,7 @@ struct HigherMassResonances {
255255
AxisSpec axisEvtResPlQA = {102, -1.02, 1.02, ""};
256256

257257
// THnSparses
258-
std::array<bool, 5> sparses = {config.activateHelicityFrame, config.activateCollinsSoperFrame, config.activateProductionFrame, config.activateBeamAxisFrame, config.activateRandomFrame};
258+
std::array<bool, 5> sparses = {config.activateHelicityFrame, config.activateCollinsSoperFrame, config.activateProductionFrame, config.activateGJFrame, config.activateRandomFrame};
259259

260260
if (std::accumulate(sparses.begin(), sparses.end(), 0) == 0) {
261261
LOGP(fatal, "No output THnSparses enabled");
@@ -269,7 +269,7 @@ struct HigherMassResonances {
269269
if (config.activateProductionFrame) {
270270
LOGP(info, "THnSparse with cosThStar w.r.t. production axis active.");
271271
}
272-
if (config.activateBeamAxisFrame) {
272+
if (config.activateGJFrame) {
273273
LOGP(info, "THnSparse with cosThStar w.r.t. beam axis active. (Gottified jackson frame)");
274274
}
275275
if (config.activateRandomFrame) {
@@ -405,13 +405,17 @@ struct HigherMassResonances {
405405
hMChists.add("Genf17102", "Gen f_{0}(1710)", kTHnSparseF, {multiplicityAxis, ptAxis, thnAxisPOL});
406406
hMChists.add("Genf1710Calib2", "Calibrated Gen f_{0}(1710)", kTHnSparseF, {multiplicityAxis, ptAxis, thnAxisPOL});
407407
hMChists.add("Gen1710PWA2", "Gen f_{0}(1710) PWA", kTHnSparseF, {glueballMassAxis, thnAxisPOL, thnAxisPhi});
408+
hMChists.add("GenThetavsPhi1", "GenThetavsPhi1", kTH2F, {{thnAxisPOL}, {thnAxisPhi}});
409+
hMChists.add("GenThetavsPhi2", "GenThetavsPhi2", kTH2F, {{thnAxisPOL}, {thnAxisPhi}});
408410

409411
hMChists.add("Recf1710_pt1", "Rec f_{0}(1710) p_{T}", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL});
410412
hMChists.add("Recf1710Calib_pt1", "Calibrated Rec f_{0}(1710) p_{T}", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL});
411413
hMChists.add("Recf1710PWA_pt1", "Rec f_{0}(1710) PWA", kTHnSparseF, {glueballMassAxis, thnAxisPOL, thnAxisPhi});
412414
hMChists.add("Recf1710_pt2", "Rec f_{0}(1710) p_{T}", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL});
413415
hMChists.add("Recf1710Calib_pt2", "Calibrated Rec f_{0}(1710) p_{T}", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL});
414416
hMChists.add("Recf1710PWA_pt2", "Rec f_{0}(1710) PWA", kTHnSparseF, {glueballMassAxis, thnAxisPOL, thnAxisPhi});
417+
hMChists.add("RecThetavsPhi1", "RecTheta vs Phi1", kTH2F, {{thnAxisPOL}, {thnAxisPhi}});
418+
hMChists.add("RecThetavsPhi2", "RecTheta vs Phi2", kTH2F, {{thnAxisPOL}, {thnAxisPhi}});
415419

416420
hMChists.add("h1Recsplit", "Rec p_{T}2", kTH1F, {ptAxis});
417421
hMChists.add("Genf1710_mass", "Gen f_{0}(1710) mass", kTH1F, {glueballMassAxis});
@@ -799,36 +803,73 @@ struct HigherMassResonances {
799803
enum FrameType { kHelicity = 0,
800804
kCollinsSoper = 1,
801805
kProduction = 2,
802-
kBeamAxis = 3,
806+
kGottfriedJackson = 3,
803807
kRandom = 4 };
804808

805809
std::pair<double, double> getCosThetaPhi(const ROOT::Math::PxPyPzMVector& motherVec, const ROOT::Math::PxPyPzMVector& daughterVec, int frame)
806810
{
811+
//----------------------------------------------------------------
812+
// Boost daughter into mother rest frame
813+
//----------------------------------------------------------------
807814
ROOT::Math::Boost boost{motherVec.BoostToCM()};
815+
808816
auto dauCM = boost(daughterVec);
809-
ROOT::Math::XYZVectorF v1cm = ROOT::Math::XYZVectorF(dauCM.Vect()).Unit();
810817

811-
ROOT::Math::XYZVectorF beam1cm = ROOT::Math::XYZVectorF((boost(beam1).Vect()).Unit());
812-
ROOT::Math::XYZVectorF beam2cm = ROOT::Math::XYZVectorF((boost(beam2).Vect()).Unit());
818+
// Unit vector of daughter momentum in mother rest frame
819+
ROOT::Math::XYZVectorF v1cm =
820+
ROOT::Math::XYZVectorF(dauCM.Vect()).Unit();
821+
822+
//----------------------------------------------------------------
823+
// Beam directions boosted into mother rest frame
824+
//----------------------------------------------------------------
825+
ROOT::Math::XYZVectorF beam1cm =
826+
ROOT::Math::XYZVectorF((boost(beam1).Vect()).Unit());
827+
828+
ROOT::Math::XYZVectorF beam2cm =
829+
ROOT::Math::XYZVectorF((boost(beam2).Vect()).Unit());
813830

814831
double cosTheta = 0.;
815832
double phi = 0.;
816833

817834
if (frame == kHelicity) {
835+
//==============================================================
836+
// HELICITY FRAME
837+
// z-axis : Mother momentum direction in laboratory frame
838+
// y-axis : Normal to beam plane
839+
// x-axis : Completes right-handed coordinate system
840+
//==============================================================
841+
818842
auto zaxisHELocal = ROOT::Math::XYZVectorF(motherVec.Vect()).Unit();
819843
auto yaxisHELocal = ROOT::Math::XYZVectorF(beam1cm.Cross(beam2cm)).Unit();
820844
auto xaxisHELocal = ROOT::Math::XYZVectorF(yaxisHELocal.Cross(zaxisHELocal)).Unit();
821845
cosTheta = motherVec.Vect().Dot(dauCM.Vect()) / (std::sqrt(dauCM.Vect().Mag2()) * std::sqrt(motherVec.Vect().Mag2()));
822846
phi = std::atan2(yaxisHELocal.Dot(v1cm), xaxisHELocal.Dot(v1cm));
823847
phi = RecoDecay::constrainAngle(phi, 0.0);
824848
} else if (frame == kCollinsSoper) {
849+
//==============================================================
850+
// COLLINS-SOPER FRAME
851+
// z-axis : Bisector of the two beam directions
852+
// y-axis : Normal to beam plane
853+
// x-axis : Completes right-handed system
854+
//==============================================================
855+
825856
auto zAxisCSLocal = ROOT::Math::XYZVectorF((beam1cm.Unit() - beam2cm.Unit())).Unit();
826857
auto yAxisCSLocal = ROOT::Math::XYZVectorF(beam1cm.Cross(beam2cm)).Unit();
827858
auto xAxisCSLocal = ROOT::Math::XYZVectorF(yAxisCSLocal.Cross(zAxisCSLocal)).Unit();
828859
cosTheta = zAxisCSLocal.Dot(v1cm);
829860
phi = std::atan2(yAxisCSLocal.Dot(v1cm), xAxisCSLocal.Dot(v1cm));
830861
phi = RecoDecay::constrainAngle(phi, 0.0);
831862
} else if (frame == kProduction) {
863+
//==============================================================
864+
// PRODUCTION FRAME
865+
// z-axis : Normal to production plane
866+
//
867+
// Production plane is defined by:
868+
// beam direction
869+
// mother momentum
870+
// Used occasionally in spin-alignment analyses.
871+
//==============================================================
872+
832873
ROOT::Math::XYZVector normalVecLocal = ROOT::Math::XYZVector(motherVec.Py(), -motherVec.Px(), 0.f);
833874
cosTheta = normalVecLocal.Dot(dauCM.Vect()) / (std::sqrt(dauCM.Vect().Mag2()) * std::sqrt(normalVecLocal.Mag2()));
834875
// for production frame reuse helicity-style azimuthal angle
@@ -837,13 +878,42 @@ struct HigherMassResonances {
837878
auto xaxisHELocal = ROOT::Math::XYZVectorF(yaxisHELocal.Cross(zaxisHELocal)).Unit();
838879
phi = std::atan2(yaxisHELocal.Dot(v1cm), xaxisHELocal.Dot(v1cm));
839880
phi = RecoDecay::constrainAngle(phi, 0.0);
840-
} else if (frame == kBeamAxis) {
841-
ROOT::Math::XYZVector beamVecLocal = ROOT::Math::XYZVector(0.f, 0.f, 1.f);
842-
cosTheta = beamVecLocal.Dot(dauCM.Vect()) / std::sqrt(dauCM.Vect().Mag2());
843-
auto zaxisHELocal = ROOT::Math::XYZVectorF(motherVec.Vect()).Unit();
844-
auto yaxisHELocal = ROOT::Math::XYZVectorF(beam1cm.Cross(beam2cm)).Unit();
845-
auto xaxisHELocal = ROOT::Math::XYZVectorF(yaxisHELocal.Cross(zaxisHELocal)).Unit();
846-
phi = std::atan2(yaxisHELocal.Dot(v1cm), xaxisHELocal.Dot(v1cm));
881+
} else if (frame == kGottfriedJackson) {
882+
// ------------------------------------------------------------------
883+
// Gottfried-Jackson (GJ') frame
884+
//
885+
// z-axis : beam direction in the mother rest frame
886+
// y-axis : normal to the production plane
887+
// x-axis : right-handed coordinate system
888+
// ------------------------------------------------------------------
889+
890+
// z-axis: beam direction boosted into mother rest frame
891+
ROOT::Math::XYZVectorF zAxisGJ = beam1cm.Unit();
892+
893+
// Mother momentum in LAB (defines production plane)
894+
ROOT::Math::XYZVectorF motherLab =
895+
ROOT::Math::XYZVectorF(motherVec.Vect()).Unit();
896+
897+
// Normal to production plane
898+
ROOT::Math::XYZVectorF yAxisGJ =
899+
ROOT::Math::XYZVectorF(zAxisGJ.Cross(motherLab));
900+
901+
if (yAxisGJ.Mag2() > 1e-12)
902+
yAxisGJ = yAxisGJ.Unit();
903+
else
904+
yAxisGJ = ROOT::Math::XYZVectorF(0.f, 1.f, 0.f);
905+
906+
// Complete right-handed system
907+
ROOT::Math::XYZVectorF xAxisGJ =
908+
ROOT::Math::XYZVectorF(yAxisGJ.Cross(zAxisGJ)).Unit();
909+
910+
// Daughter direction in mother rest frame
911+
cosTheta = zAxisGJ.Dot(v1cm);
912+
913+
phi = std::atan2(
914+
yAxisGJ.Dot(v1cm),
915+
xAxisGJ.Dot(v1cm));
916+
847917
phi = RecoDecay::constrainAngle(phi, 0.0);
848918
} else { // kRandom or fallback
849919
auto phiRandom = gRandom->Uniform(0.f, constants::math::TwoPI);
@@ -1008,16 +1078,16 @@ struct HigherMassResonances {
10081078
hglue.fill(HIST("h3glueInvMassME"), eventMultiplicity, motherVec.Pt(), motherVec.M(), pr.first, pr.second);
10091079
}
10101080
}
1011-
} else if (config.activateBeamAxisFrame) {
1012-
auto pr = getCosThetaPhi(motherVec, daughterVec1, kBeamAxis);
1081+
} else if (config.activateGJFrame) {
1082+
auto pr = getCosThetaPhi(motherVec, daughterVec1, kGottfriedJackson);
10131083
if (!isMixed) {
10141084
if (std::abs(motherVec.Rapidity()) < config.rapidityMotherData) {
10151085
hglue.fill(HIST("h3glueInvMassDS"), eventMultiplicity, motherVec.Pt(), motherVec.M(), pr.first, pr.second);
10161086
}
10171087
for (int i = 0; i < config.cRotations; i++) {
10181088
config.theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut);
10191089
motherRot = ROOT::Math::PxPyPzMVector(motherVec.Px() * std::cos(config.theta2) - motherVec.Py() * std::sin(config.theta2), motherVec.Px() * std::sin(config.theta2) + motherVec.Py() * std::cos(config.theta2), motherVec.Pz(), motherVec.M());
1020-
auto prrot = getCosThetaPhi(motherRot, daughterVec1, kBeamAxis);
1090+
auto prrot = getCosThetaPhi(motherRot, daughterVec1, kGottfriedJackson);
10211091
if (std::abs(motherRot.Rapidity()) < config.rapidityMotherData) {
10221092
hglue.fill(HIST("h3glueInvMassRot"), eventMultiplicity, motherRot.Pt(), motherRot.M(), prrot.first, prrot.second);
10231093
}
@@ -1612,8 +1682,8 @@ struct HigherMassResonances {
16121682
frame = kCollinsSoper;
16131683
} else if (config.activateProductionFrame) {
16141684
frame = kProduction;
1615-
} else if (config.activateBeamAxisFrame) {
1616-
frame = kBeamAxis;
1685+
} else if (config.activateGJFrame) {
1686+
frame = kGottfriedJackson;
16171687
} else if (config.activateRandomFrame) {
16181688
frame = kRandom;
16191689
} else {
@@ -1630,6 +1700,7 @@ struct HigherMassResonances {
16301700
hMChists.fill(HIST("GenEta"), mcParticle.eta());
16311701
hMChists.fill(HIST("GenPhi"), mcParticle.phi());
16321702
hMChists.fill(HIST("Gen1710PWA"), lResonanceGen.M(), angularVar.first, angularVar.second);
1703+
hMChists.fill(HIST("GenThetavsPhi1"), angularVar.first, angularVar.second);
16331704

16341705
if (config.isapplyPairRapidityMC && std::abs(lResonanceGen1.Rapidity()) >= config.rapidityMotherData) {
16351706
continue;
@@ -1642,6 +1713,7 @@ struct HigherMassResonances {
16421713
hMChists.fill(HIST("GenEta2"), lResonanceGen1.Eta());
16431714
hMChists.fill(HIST("GenPhi2"), lResonanceGen1.Phi());
16441715
hMChists.fill(HIST("Gen1710PWA2"), lResonanceGen1.M(), angularVar1.first, angularVar1.second);
1716+
hMChists.fill(HIST("GenThetavsPhi2"), angularVar1.first, angularVar1.second);
16451717
}
16461718
passKs.clear(); // clear the vector for the next iteration
16471719
}
@@ -1903,8 +1975,8 @@ struct HigherMassResonances {
19031975
frame = kCollinsSoper;
19041976
} else if (config.activateProductionFrame) {
19051977
frame = kProduction;
1906-
} else if (config.activateBeamAxisFrame) {
1907-
frame = kBeamAxis;
1978+
} else if (config.activateGJFrame) {
1979+
frame = kGottfriedJackson;
19081980
} else if (config.activateRandomFrame) {
19091981
frame = kRandom;
19101982
}
@@ -1918,6 +1990,7 @@ struct HigherMassResonances {
19181990
hMChists.fill(HIST("RecPhi"), mothertrack1.phi());
19191991
hMChists.fill(HIST("RecEta"), mothertrack1.eta());
19201992
hMChists.fill(HIST("Recf1710PWA_pt1"), mother1.M(), angularVar1.first, angularVar1.second);
1993+
hMChists.fill(HIST("RecThetavsPhi1"), angularVar1.first, angularVar1.second);
19211994

19221995
if (config.isapplyPairRapidityMC && std::abs(mother.Rapidity()) >= config.rapidityMotherData) {
19231996
continue;
@@ -1929,6 +2002,7 @@ struct HigherMassResonances {
19292002
hMChists.fill(HIST("RecPhi2"), mother.Phi());
19302003
hMChists.fill(HIST("RecEta2"), mother.Eta());
19312004
hMChists.fill(HIST("Recf1710PWA_pt2"), mother.M(), angularVar.first, angularVar.second);
2005+
hMChists.fill(HIST("RecThetavsPhi2"), angularVar.first, angularVar.second);
19322006
}
19332007
gindex2.clear();
19342008
}

0 commit comments

Comments
 (0)