Skip to content

Commit 99fc22f

Browse files
sawankumawatsawankumawat
authored
[PWGLF] Added Armenteros-Podolanski cut and QA histograms (#17351)
1 parent 45438cf commit 99fc22f

2 files changed

Lines changed: 253 additions & 135 deletions

File tree

‎PWGLF/Tasks/Resonances/higherMassResonances.cxx‎

Lines changed: 99 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -157,6 +157,7 @@ struct HigherMassResonances {
157157
Configurable<float> cfgETAcut{"cfgETAcut", 0.8f, "Track ETA cut"};
158158
Configurable<float> deltaRDaugherCut{"deltaRDaugherCut", 0.001f, "DeltaR cut on V0 daughters"};
159159
Configurable<float> deltaRK0sCut{"deltaRK0sCut", 0.1f, "Apply deltaR cut between two K0s"};
160+
Configurable<float> cfgArmenterosCut{"cfgArmenterosCut", 0.2f, "Armenteros-Podolanski cut"};
160161

161162
// Configurable for track selection and multiplicity
162163
Configurable<float> cfgPTcut{"cfgPTcut", 0.2f, "Track PT cut"};
@@ -174,7 +175,7 @@ struct HigherMassResonances {
174175
Configurable<bool> activateHelicityFrame{"activateHelicityFrame", false, "Activate the THnSparse with cosThStar w.r.t. helicity axis"};
175176
Configurable<bool> activateCollinsSoperFrame{"activateCollinsSoperFrame", false, "Activate the THnSparse with cosThStar w.r.t. Collins soper axis"};
176177
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)"};
178+
Configurable<bool> activateGJFrame{"activateGJFrame", true, "Activate the THnSparse with cosThStar w.r.t. beam axis (Gottified jackson frame)"};
178179
Configurable<bool> activateRandomFrame{"activateRandomFrame", false, "Activate the THnSparse with cosThStar w.r.t. random axis"};
179180
Configurable<int> cRotations{"cRotations", 3, "Number of random rotations in the rotational background"};
180181

@@ -212,6 +213,7 @@ struct HigherMassResonances {
212213
int refAId = 0;
213214
int refBId = 0;
214215
float minQvecAmp = 1e-5;
216+
double tolerance = 1e-12;
215217
} config;
216218

217219
// Service<o2::framework::O2DatabasePDG> PDGdatabase;
@@ -255,7 +257,7 @@ struct HigherMassResonances {
255257
AxisSpec axisEvtResPlQA = {102, -1.02, 1.02, ""};
256258

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

260262
if (std::accumulate(sparses.begin(), sparses.end(), 0) == 0) {
261263
LOGP(fatal, "No output THnSparses enabled");
@@ -269,7 +271,7 @@ struct HigherMassResonances {
269271
if (config.activateProductionFrame) {
270272
LOGP(info, "THnSparse with cosThStar w.r.t. production axis active.");
271273
}
272-
if (config.activateBeamAxisFrame) {
274+
if (config.activateGJFrame) {
273275
LOGP(info, "THnSparse with cosThStar w.r.t. beam axis active. (Gottified jackson frame)");
274276
}
275277
if (config.activateRandomFrame) {
@@ -311,10 +313,11 @@ struct HigherMassResonances {
311313
hv0label->GetXaxis()->SetBinLabel(5, "Daughter DCA");
312314
hv0label->GetXaxis()->SetBinLabel(6, "CosPA");
313315
hv0label->GetXaxis()->SetBinLabel(7, "Decay Radius");
314-
hv0label->GetXaxis()->SetBinLabel(8, "Lifetime");
315-
hv0label->GetXaxis()->SetBinLabel(9, "CompetingCascade");
316-
hv0label->GetXaxis()->SetBinLabel(10, "Standard V0");
317-
hv0label->GetXaxis()->SetBinLabel(11, "Mass Tolerance");
316+
hv0label->GetXaxis()->SetBinLabel(8, "Armenteros-Podolanski");
317+
hv0label->GetXaxis()->SetBinLabel(9, "Lifetime");
318+
hv0label->GetXaxis()->SetBinLabel(10, "CompetingCascade");
319+
hv0label->GetXaxis()->SetBinLabel(11, "Standard V0");
320+
hv0label->GetXaxis()->SetBinLabel(12, "Mass Tolerance");
318321

319322
std::shared_ptr<TH1> hv0DauLabel = rEventSelection.get<TH1>(HIST("htrackscheck_v0_daughters"));
320323
hv0DauLabel->GetXaxis()->SetBinLabel(1, "AllDau Tracks");
@@ -385,6 +388,7 @@ struct HigherMassResonances {
385388
rKzeroShort.add("hLT", "hLT", {HistType::kTH1F, {{100, 0.0f, 50.0f}}});
386389
rKzeroShort.add("angularSeparation", "Angular distribution between two K0s vs pT", {HistType::kTH1F, {{200, 0.0f, 4.0f}}});
387390
rKzeroShort.add("hDauDeltaR", "Delta R of positive and negative daughers", {HistType::kTHnSparseF, {angleSepAxis, angleSepAxis}});
391+
rKzeroShort.add("hArmenterosPodolanski", "Armenteros-Podolanski plot", HistType::kTH2D, {AxisSpec{100, -1, 1, "#alpha"}, {200, 0, 0.5, "qtArm"}});
388392
}
389393
rKzeroShort.add("NksProduced", "Number of K0s produced", kTH1I, {{15, -0.5, 14.5}});
390394

@@ -405,13 +409,17 @@ struct HigherMassResonances {
405409
hMChists.add("Genf17102", "Gen f_{0}(1710)", kTHnSparseF, {multiplicityAxis, ptAxis, thnAxisPOL});
406410
hMChists.add("Genf1710Calib2", "Calibrated Gen f_{0}(1710)", kTHnSparseF, {multiplicityAxis, ptAxis, thnAxisPOL});
407411
hMChists.add("Gen1710PWA2", "Gen f_{0}(1710) PWA", kTHnSparseF, {glueballMassAxis, thnAxisPOL, thnAxisPhi});
412+
hMChists.add("GenThetavsPhi1", "GenThetavsPhi1", kTH2F, {{thnAxisPOL}, {thnAxisPhi}});
413+
hMChists.add("GenThetavsPhi2", "GenThetavsPhi2", kTH2F, {{thnAxisPOL}, {thnAxisPhi}});
408414

409415
hMChists.add("Recf1710_pt1", "Rec f_{0}(1710) p_{T}", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL});
410416
hMChists.add("Recf1710Calib_pt1", "Calibrated Rec f_{0}(1710) p_{T}", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL});
411417
hMChists.add("Recf1710PWA_pt1", "Rec f_{0}(1710) PWA", kTHnSparseF, {glueballMassAxis, thnAxisPOL, thnAxisPhi});
412418
hMChists.add("Recf1710_pt2", "Rec f_{0}(1710) p_{T}", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL});
413419
hMChists.add("Recf1710Calib_pt2", "Calibrated Rec f_{0}(1710) p_{T}", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL});
414420
hMChists.add("Recf1710PWA_pt2", "Rec f_{0}(1710) PWA", kTHnSparseF, {glueballMassAxis, thnAxisPOL, thnAxisPhi});
421+
hMChists.add("RecThetavsPhi1", "RecTheta vs Phi1", kTH2F, {{thnAxisPOL}, {thnAxisPhi}});
422+
hMChists.add("RecThetavsPhi2", "RecTheta vs Phi2", kTH2F, {{thnAxisPOL}, {thnAxisPhi}});
415423

416424
hMChists.add("h1Recsplit", "Rec p_{T}2", kTH1F, {ptAxis});
417425
hMChists.add("Genf1710_mass", "Gen f_{0}(1710) mass", kTH1F, {glueballMassAxis});
@@ -630,9 +638,6 @@ struct HigherMassResonances {
630638
template <typename Collision, typename V0>
631639
bool selectionV0(Collision const& collision, V0 const& candidate, float /*multiplicity*/)
632640
{
633-
// const float qtarm = candidate.qtarm();
634-
// const float alph = candidate.alpha();
635-
// float arm = qtarm / alph;
636641
const float pT = candidate.pt();
637642
const float tranRad = candidate.v0radius();
638643
const float dcaDaughv0 = candidate.dcaV0daughters();
@@ -647,6 +652,7 @@ struct HigherMassResonances {
647652
rKzeroShort.fill(HIST("hLT"), ctauK0s);
648653
rKzeroShort.fill(HIST("hDCAV0Daughters"), candidate.dcaV0daughters());
649654
rKzeroShort.fill(HIST("hV0CosPA"), candidate.v0cosPA());
655+
rKzeroShort.fill(HIST("hArmenterosPodolanski"), candidate.alpha(), candidate.qtarm());
650656
}
651657
if (config.qAcorrelation2Dhist)
652658
rKzeroShort.fill(HIST("mass_lambda_kshort_before"), candidate.mK0Short(), candidate.mLambda());
@@ -683,9 +689,9 @@ struct HigherMassResonances {
683689
}
684690
rEventSelection.fill(HIST("htrackscheck_v0"), 6.5);
685691

686-
// if (tranRad > config.confV0TranRadV0Max) {
687-
// return false;
688-
// }
692+
if (candidate.qtarm() < config.cfgArmenterosCut * std::fabs(candidate.alpha())) {
693+
return false;
694+
}
689695
rEventSelection.fill(HIST("htrackscheck_v0"), 7.5);
690696

691697
if (std::fabs(ctauK0s) > config.cMaxV0LifeTime) {
@@ -799,36 +805,70 @@ struct HigherMassResonances {
799805
enum FrameType { kHelicity = 0,
800806
kCollinsSoper = 1,
801807
kProduction = 2,
802-
kBeamAxis = 3,
808+
kGottfriedJackson = 3,
803809
kRandom = 4 };
804810

805811
std::pair<double, double> getCosThetaPhi(const ROOT::Math::PxPyPzMVector& motherVec, const ROOT::Math::PxPyPzMVector& daughterVec, int frame)
806812
{
813+
//----------------------------------------------------------------
814+
// Boost daughter into mother rest frame
815+
//----------------------------------------------------------------
807816
ROOT::Math::Boost boost{motherVec.BoostToCM()};
817+
808818
auto dauCM = boost(daughterVec);
809-
ROOT::Math::XYZVectorF v1cm = ROOT::Math::XYZVectorF(dauCM.Vect()).Unit();
810819

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

814830
double cosTheta = 0.;
815831
double phi = 0.;
816832

817833
if (frame == kHelicity) {
834+
//==============================================================
835+
// HELICITY FRAME
836+
// z-axis : Mother momentum direction in laboratory frame
837+
// y-axis : Normal to beam plane
838+
// x-axis : Completes right-handed coordinate system
839+
//==============================================================
840+
818841
auto zaxisHELocal = ROOT::Math::XYZVectorF(motherVec.Vect()).Unit();
819842
auto yaxisHELocal = ROOT::Math::XYZVectorF(beam1cm.Cross(beam2cm)).Unit();
820843
auto xaxisHELocal = ROOT::Math::XYZVectorF(yaxisHELocal.Cross(zaxisHELocal)).Unit();
821844
cosTheta = motherVec.Vect().Dot(dauCM.Vect()) / (std::sqrt(dauCM.Vect().Mag2()) * std::sqrt(motherVec.Vect().Mag2()));
822845
phi = std::atan2(yaxisHELocal.Dot(v1cm), xaxisHELocal.Dot(v1cm));
823846
phi = RecoDecay::constrainAngle(phi, 0.0);
824847
} else if (frame == kCollinsSoper) {
848+
//==============================================================
849+
// COLLINS-SOPER FRAME
850+
// z-axis : Bisector of the two beam directions
851+
// y-axis : Normal to beam plane
852+
// x-axis : Completes right-handed system
853+
//==============================================================
854+
825855
auto zAxisCSLocal = ROOT::Math::XYZVectorF((beam1cm.Unit() - beam2cm.Unit())).Unit();
826856
auto yAxisCSLocal = ROOT::Math::XYZVectorF(beam1cm.Cross(beam2cm)).Unit();
827857
auto xAxisCSLocal = ROOT::Math::XYZVectorF(yAxisCSLocal.Cross(zAxisCSLocal)).Unit();
828858
cosTheta = zAxisCSLocal.Dot(v1cm);
829859
phi = std::atan2(yAxisCSLocal.Dot(v1cm), xAxisCSLocal.Dot(v1cm));
830860
phi = RecoDecay::constrainAngle(phi, 0.0);
831861
} else if (frame == kProduction) {
862+
//==============================================================
863+
// PRODUCTION FRAME
864+
// z-axis : Normal to production plane
865+
//
866+
// Production plane is defined by:
867+
// beam direction
868+
// mother momentum
869+
// Used occasionally in spin-alignment analyses.
870+
//==============================================================
871+
832872
ROOT::Math::XYZVector normalVecLocal = ROOT::Math::XYZVector(motherVec.Py(), -motherVec.Px(), 0.f);
833873
cosTheta = normalVecLocal.Dot(dauCM.Vect()) / (std::sqrt(dauCM.Vect().Mag2()) * std::sqrt(normalVecLocal.Mag2()));
834874
// for production frame reuse helicity-style azimuthal angle
@@ -837,13 +877,37 @@ struct HigherMassResonances {
837877
auto xaxisHELocal = ROOT::Math::XYZVectorF(yaxisHELocal.Cross(zaxisHELocal)).Unit();
838878
phi = std::atan2(yaxisHELocal.Dot(v1cm), xaxisHELocal.Dot(v1cm));
839879
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));
880+
} else if (frame == kGottfriedJackson) {
881+
// ------------------------------------------------------------------
882+
// Gottfried-Jackson (GJ') frame
883+
//
884+
// z-axis : beam direction in the mother rest frame
885+
// y-axis : normal to the production plane
886+
// x-axis : right-handed coordinate system
887+
// ------------------------------------------------------------------
888+
889+
// z-axis: beam direction boosted into mother rest frame
890+
auto zAxisGJ = beam1cm.Unit();
891+
892+
// Mother momentum in LAB (defines production plane)
893+
auto motherLab = ROOT::Math::XYZVectorF(motherVec.Vect()).Unit();
894+
895+
// Normal to production plane
896+
auto yAxisGJ = ROOT::Math::XYZVectorF(zAxisGJ.Cross(motherLab));
897+
898+
if (yAxisGJ.Mag2() > config.tolerance)
899+
yAxisGJ = yAxisGJ.Unit();
900+
else
901+
yAxisGJ = ROOT::Math::XYZVectorF(0.f, 1.f, 0.f);
902+
903+
// Complete right-handed system
904+
ROOT::Math::XYZVectorF xAxisGJ = ROOT::Math::XYZVectorF(yAxisGJ.Cross(zAxisGJ)).Unit();
905+
906+
// Daughter direction in mother rest frame
907+
cosTheta = zAxisGJ.Dot(v1cm);
908+
909+
phi = std::atan2(yAxisGJ.Dot(v1cm), xAxisGJ.Dot(v1cm));
910+
847911
phi = RecoDecay::constrainAngle(phi, 0.0);
848912
} else { // kRandom or fallback
849913
auto phiRandom = gRandom->Uniform(0.f, constants::math::TwoPI);
@@ -1008,16 +1072,16 @@ struct HigherMassResonances {
10081072
hglue.fill(HIST("h3glueInvMassME"), eventMultiplicity, motherVec.Pt(), motherVec.M(), pr.first, pr.second);
10091073
}
10101074
}
1011-
} else if (config.activateBeamAxisFrame) {
1012-
auto pr = getCosThetaPhi(motherVec, daughterVec1, kBeamAxis);
1075+
} else if (config.activateGJFrame) {
1076+
auto pr = getCosThetaPhi(motherVec, daughterVec1, kGottfriedJackson);
10131077
if (!isMixed) {
10141078
if (std::abs(motherVec.Rapidity()) < config.rapidityMotherData) {
10151079
hglue.fill(HIST("h3glueInvMassDS"), eventMultiplicity, motherVec.Pt(), motherVec.M(), pr.first, pr.second);
10161080
}
10171081
for (int i = 0; i < config.cRotations; i++) {
10181082
config.theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut);
10191083
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);
1084+
auto prrot = getCosThetaPhi(motherRot, daughterVec1, kGottfriedJackson);
10211085
if (std::abs(motherRot.Rapidity()) < config.rapidityMotherData) {
10221086
hglue.fill(HIST("h3glueInvMassRot"), eventMultiplicity, motherRot.Pt(), motherRot.M(), prrot.first, prrot.second);
10231087
}
@@ -1612,8 +1676,8 @@ struct HigherMassResonances {
16121676
frame = kCollinsSoper;
16131677
} else if (config.activateProductionFrame) {
16141678
frame = kProduction;
1615-
} else if (config.activateBeamAxisFrame) {
1616-
frame = kBeamAxis;
1679+
} else if (config.activateGJFrame) {
1680+
frame = kGottfriedJackson;
16171681
} else if (config.activateRandomFrame) {
16181682
frame = kRandom;
16191683
} else {
@@ -1630,6 +1694,7 @@ struct HigherMassResonances {
16301694
hMChists.fill(HIST("GenEta"), mcParticle.eta());
16311695
hMChists.fill(HIST("GenPhi"), mcParticle.phi());
16321696
hMChists.fill(HIST("Gen1710PWA"), lResonanceGen.M(), angularVar.first, angularVar.second);
1697+
hMChists.fill(HIST("GenThetavsPhi1"), angularVar.first, angularVar.second);
16331698

16341699
if (config.isapplyPairRapidityMC && std::abs(lResonanceGen1.Rapidity()) >= config.rapidityMotherData) {
16351700
continue;
@@ -1642,6 +1707,7 @@ struct HigherMassResonances {
16421707
hMChists.fill(HIST("GenEta2"), lResonanceGen1.Eta());
16431708
hMChists.fill(HIST("GenPhi2"), lResonanceGen1.Phi());
16441709
hMChists.fill(HIST("Gen1710PWA2"), lResonanceGen1.M(), angularVar1.first, angularVar1.second);
1710+
hMChists.fill(HIST("GenThetavsPhi2"), angularVar1.first, angularVar1.second);
16451711
}
16461712
passKs.clear(); // clear the vector for the next iteration
16471713
}
@@ -1903,8 +1969,8 @@ struct HigherMassResonances {
19031969
frame = kCollinsSoper;
19041970
} else if (config.activateProductionFrame) {
19051971
frame = kProduction;
1906-
} else if (config.activateBeamAxisFrame) {
1907-
frame = kBeamAxis;
1972+
} else if (config.activateGJFrame) {
1973+
frame = kGottfriedJackson;
19081974
} else if (config.activateRandomFrame) {
19091975
frame = kRandom;
19101976
}
@@ -1918,6 +1984,7 @@ struct HigherMassResonances {
19181984
hMChists.fill(HIST("RecPhi"), mothertrack1.phi());
19191985
hMChists.fill(HIST("RecEta"), mothertrack1.eta());
19201986
hMChists.fill(HIST("Recf1710PWA_pt1"), mother1.M(), angularVar1.first, angularVar1.second);
1987+
hMChists.fill(HIST("RecThetavsPhi1"), angularVar1.first, angularVar1.second);
19211988

19221989
if (config.isapplyPairRapidityMC && std::abs(mother.Rapidity()) >= config.rapidityMotherData) {
19231990
continue;
@@ -1929,6 +1996,7 @@ struct HigherMassResonances {
19291996
hMChists.fill(HIST("RecPhi2"), mother.Phi());
19301997
hMChists.fill(HIST("RecEta2"), mother.Eta());
19311998
hMChists.fill(HIST("Recf1710PWA_pt2"), mother.M(), angularVar.first, angularVar.second);
1999+
hMChists.fill(HIST("RecThetavsPhi2"), angularVar.first, angularVar.second);
19322000
}
19332001
gindex2.clear();
19342002
}

0 commit comments

Comments
 (0)