diff --git a/PWGLF/Tasks/Resonances/higherMassResonances.cxx b/PWGLF/Tasks/Resonances/higherMassResonances.cxx index 792916cb6fd..1a8b6d003f7 100644 --- a/PWGLF/Tasks/Resonances/higherMassResonances.cxx +++ b/PWGLF/Tasks/Resonances/higherMassResonances.cxx @@ -213,8 +213,8 @@ struct HigherMassResonances { // variables declaration float multiplicity = 0.0f; - float theta2; - ROOT::Math::PxPyPzMVector daughter1, daughter2, daughterRot, daughterRotCM, mother, motherRot, fourVecDauCM, fourVecDauCM1; + float theta2 -= 0.0f; + ROOT::Math::PxPyPzMVector daughter1, daughter2, daughterRot, daughterRotCM, mother, motherRot, fourVecDauCM, fourVecDauCM1, mother2, mother3; ROOT::Math::PxPyPzEVector mother1; ROOT::Math::XYZVector randomVec, beamVec, normalVec; ROOT::Math::XYZVectorF v1CM, zaxisHE, yaxisHE, xaxisHE; @@ -228,9 +228,9 @@ struct HigherMassResonances { bool isMix = false; EventPlaneHelper helperEP; - int detId; - int refAId; - int refBId; + int detId = 0; + int refAId = 0; + int refBId = 0; float minQvecAmp = 1e-5; void init(InitContext const&) @@ -399,12 +399,18 @@ struct HigherMassResonances { hMChists.add("MCcorrections/hGenNo", "Generated collisions before and after event selection", kTH1F, {{5, 0.0f, 5.0f}}); hMChists.add("Genf1710", "Gen f_{0}(1710)", kTHnSparseF, {multiplicityAxis, ptAxis, thnAxisPOL}); hMChists.add("Genf1710Calib", "Calibrated Gen f_{0}(1710)", kTHnSparseF, {multiplicityAxis, ptAxis, thnAxisPOL}); + hMChists.add("Gen1710PWA", "Gen f_{0}(1710) PWA", kTHnSparseF, {glueballMassAxis, thnAxisPOL, thnAxisPhi}); hMChists.add("Genf17102", "Gen f_{0}(1710)", kTHnSparseF, {multiplicityAxis, ptAxis, thnAxisPOL}); hMChists.add("Genf1710Calib2", "Calibrated Gen f_{0}(1710)", kTHnSparseF, {multiplicityAxis, ptAxis, thnAxisPOL}); + hMChists.add("Gen1710PWA2", "Gen f_{0}(1710) PWA", kTHnSparseF, {glueballMassAxis, thnAxisPOL, thnAxisPhi}); + hMChists.add("Recf1710_pt1", "Rec f_{0}(1710) p_{T}", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL}); hMChists.add("Recf1710Calib_pt1", "Calibrated Rec f_{0}(1710) p_{T}", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL}); + hMChists.add("Recf1710PWA_pt1", "Rec f_{0}(1710) PWA", kTHnSparseF, {glueballMassAxis, thnAxisPOL, thnAxisPhi}); hMChists.add("Recf1710_pt2", "Rec f_{0}(1710) p_{T}", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL}); hMChists.add("Recf1710Calib_pt2", "Calibrated Rec f_{0}(1710) p_{T}", kTHnSparseF, {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL}); + hMChists.add("Recf1710PWA_pt2", "Rec f_{0}(1710) PWA", kTHnSparseF, {glueballMassAxis, thnAxisPOL, thnAxisPhi}); + hMChists.add("h1Recsplit", "Rec p_{T}2", kTH1F, {ptAxis}); hMChists.add("Genf1710_mass", "Gen f_{0}(1710) mass", kTH1F, {glueballMassAxis}); hMChists.add("Genf1710_mass2", "Gen f_{0}(1710) mass", kTH1F, {glueballMassAxis}); @@ -783,6 +789,66 @@ struct HigherMassResonances { return std::sqrt(d1 * d1 + d2 * d2); } + enum FrameType { kHelicity = 0, + kCollinsSoper = 1, + kProduction = 2, + kBeamAxis = 3, + kRandom = 4 }; + + std::pair getCosThetaPhi(const ROOT::Math::PxPyPzMVector& motherVec, const ROOT::Math::PxPyPzMVector& daughterVec, int frame) + { + ROOT::Math::Boost boost{motherVec.BoostToCM()}; + auto dauCM = boost(daughterVec); + ROOT::Math::XYZVectorF v1cm = ROOT::Math::XYZVectorF(dauCM.Vect()).Unit(); + + ROOT::Math::XYZVectorF beam1cm = ROOT::Math::XYZVectorF((boost(beam1).Vect()).Unit()); + ROOT::Math::XYZVectorF beam2cm = ROOT::Math::XYZVectorF((boost(beam2).Vect()).Unit()); + + double cosTheta = 0.; + double phi = 0.; + + if (frame == kHelicity) { + auto zaxisHE = ROOT::Math::XYZVectorF(motherVec.Vect()).Unit(); + auto yaxisHE = ROOT::Math::XYZVectorF(beam1cm.Cross(beam2cm)).Unit(); + auto xaxisHE = ROOT::Math::XYZVectorF(yaxisHE.Cross(zaxisHE)).Unit(); + cosTheta = motherVec.Vect().Dot(dauCM.Vect()) / (std::sqrt(dauCM.Vect().Mag2()) * std::sqrt(motherVec.Vect().Mag2())); + phi = std::atan2(yaxisHE.Dot(v1cm), xaxisHE.Dot(v1cm)); + phi = RecoDecay::constrainAngle(phi, 0.0); + } else if (frame == kCollinsSoper) { + auto zAxisCS = ROOT::Math::XYZVectorF((beam1cm.Unit() - beam2cm.Unit())).Unit(); + auto yAxisCS = ROOT::Math::XYZVectorF(beam1cm.Cross(beam2cm)).Unit(); + auto xAxisCS = ROOT::Math::XYZVectorF(yAxisCS.Cross(zAxisCS)).Unit(); + cosTheta = zAxisCS.Dot(v1cm); + phi = std::atan2(yAxisCS.Dot(v1cm), xAxisCS.Dot(v1cm)); + phi = RecoDecay::constrainAngle(phi, 0.0); + } else if (frame == kProduction) { + ROOT::Math::XYZVector normalVecLocal = ROOT::Math::XYZVector(motherVec.Py(), -motherVec.Px(), 0.f); + cosTheta = normalVecLocal.Dot(dauCM.Vect()) / (std::sqrt(dauCM.Vect().Mag2()) * std::sqrt(normalVecLocal.Mag2())); + // for production frame reuse helicity-style azimuthal angle + auto zaxisHE = ROOT::Math::XYZVectorF(motherVec.Vect()).Unit(); + auto yaxisHE = ROOT::Math::XYZVectorF(beam1cm.Cross(beam2cm)).Unit(); + auto xaxisHE = ROOT::Math::XYZVectorF(yaxisHE.Cross(zaxisHE)).Unit(); + phi = std::atan2(yaxisHE.Dot(v1cm), xaxisHE.Dot(v1cm)); + phi = RecoDecay::constrainAngle(phi, 0.0); + } else if (frame == kBeamAxis) { + ROOT::Math::XYZVector beamVecLocal = ROOT::Math::XYZVector(0.f, 0.f, 1.f); + cosTheta = beamVecLocal.Dot(dauCM.Vect()) / std::sqrt(dauCM.Vect().Mag2()); + auto zaxisHE = ROOT::Math::XYZVectorF(motherVec.Vect()).Unit(); + auto yaxisHE = ROOT::Math::XYZVectorF(beam1cm.Cross(beam2cm)).Unit(); + auto xaxisHE = ROOT::Math::XYZVectorF(yaxisHE.Cross(zaxisHE)).Unit(); + phi = std::atan2(yaxisHE.Dot(v1cm), xaxisHE.Dot(v1cm)); + phi = RecoDecay::constrainAngle(phi, 0.0); + } else { // kRandom or fallback + auto phiRandom = gRandom->Uniform(0.f, constants::math::TwoPI); + auto thetaRandom = gRandom->Uniform(0.f, constants::math::PI); + ROOT::Math::XYZVector randomVecLocal = ROOT::Math::XYZVector(std::sin(thetaRandom) * std::cos(phiRandom), std::sin(thetaRandom) * std::sin(phiRandom), std::cos(thetaRandom)); + cosTheta = randomVecLocal.Dot(dauCM.Vect()) / std::sqrt(dauCM.Vect().Mag2()); + phi = phiRandom; + } + + return {cosTheta, phi}; + } + using EventCandidatesDerivedData = soa::Join; using V0CandidatesDerivedData = soa::Join; // using DauTracks = soa::Join; @@ -869,69 +935,11 @@ struct HigherMassResonances { void fillInvMass(const T& mother, float multiplicity, const T& daughter1, const T& daughter2, bool isMix) { - // //polarization calculations - // zBeam = ROOT::Math::XYZVector(0.f, 0.f, 1.f); // ẑ: beam direction in lab frame - - ROOT::Math::Boost boost{mother.BoostToCM()}; // define the boost to the center of mass frame - fourVecDauCM = boost(daughter1); // boost the frame of daughter to the center of mass frame - // threeVecDauCM = fourVecDauCM.Vect(); // get the 3 vector of daughter in the frame of mother - - 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 - - // auto motherLabDirection = ROOT::Math::XYZVector(0, 0, mother.Vect().Z()); // ẑ axis in lab frame - - // // ŷ = zBeam × ẑ - // auto y_axis = zBeam.Cross(motherLabDirection).Unit(); - - // // x̂ = ŷ × ẑ - // auto x_axis = y_axis.Cross(motherLabDirection).Unit(); - - // // Project daughter momentum onto x–y plane - // auto p_proj_x = threeVecDauCM.Dot(x_axis); - // auto p_proj_y = threeVecDauCM.Dot(y_axis); - - // // 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 - zaxisHE = ROOT::Math::XYZVectorF(mother.Vect()).Unit(); - yaxisHE = ROOT::Math::XYZVectorF(beam1CM.Cross(beam2CM)).Unit(); - xaxisHE = ROOT::Math::XYZVectorF(yaxisHE.Cross(zaxisHE)).Unit(); - - // CosThetaHE = zaxisHE.Dot(v_CM); - - auto anglePhi = std::atan2(yaxisHE.Dot(v1CM), xaxisHE.Dot(v1CM)); - anglePhi = RecoDecay::constrainAngle(anglePhi, 0.0); - // if (anglePhi < 0) { - // 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.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())); + auto pr = getCosThetaPhi(mother, daughter1, kHelicity); if (!isMix) { if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity, anglePhi); + hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), pr.first, pr.second); } for (int i = 0; i < config.cRotations; i++) { @@ -941,24 +949,20 @@ struct HigherMassResonances { motherRot = daughterRot + daughter2; - ROOT::Math::Boost boost2{motherRot.BoostToCM()}; - 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); + auto prrot = getCosThetaPhi(motherRot, daughterRot, kHelicity); if (motherRot.Rapidity() < config.rapidityMotherData) - hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarHelicityRot, phiHelicityRot); + hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), prrot.first, prrot.second); } } else { if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity, anglePhi); + hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), pr.first, pr.second); } } } else if (config.activateCollinsSoperFrame) { + auto pr = getCosThetaPhi(mother, daughter1, kCollinsSoper); if (!isMix) { if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarCS, phiCS); + hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), pr.first, pr.second); } for (int i = 0; i < config.cRotations; i++) { @@ -968,83 +972,74 @@ struct HigherMassResonances { 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); + auto prrot = getCosThetaPhi(motherRot, daughterRot, kCollinsSoper); if (motherRot.Rapidity() < config.rapidityMotherData) - hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarCSrot, phiCSrot); + hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), prrot.first, prrot.second); } } else { if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarCS, phiCS); + hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), pr.first, pr.second); } } } else if (config.activateProductionFrame) { - normalVec = ROOT::Math::XYZVector(mother.Py(), -mother.Px(), 0.f); - auto cosThetaProduction = normalVec.Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(normalVec.Mag2())); + auto pr = getCosThetaPhi(mother, daughter1, kProduction); if (!isMix) { if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaProduction, anglePhi); + hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), pr.first, pr.second); } 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()); + auto prrot = getCosThetaPhi(motherRot, daughter1, kProduction); if (std::abs(motherRot.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaProduction, anglePhi); + hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), prrot.first, prrot.second); } } } else { if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaProduction, anglePhi); + hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), pr.first, pr.second); } } } 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()); + auto pr = getCosThetaPhi(mother, daughter1, kBeamAxis); if (!isMix) { if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarBeam, anglePhi); + hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), pr.first, pr.second); } 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()); + auto prrot = getCosThetaPhi(motherRot, daughter1, kBeamAxis); if (std::abs(motherRot.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarBeam, anglePhi); + hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), prrot.first, prrot.second); } } } else { if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarBeam, anglePhi); + hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), pr.first, pr.second); } } } 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()); + auto pr = getCosThetaPhi(mother, daughter1, kRandom); if (!isMix) { if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, phiRandom); + hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), pr.first, pr.second); } 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()); + auto prrot = getCosThetaPhi(motherRot, daughter1, kRandom); if (std::abs(motherRot.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarRandom, phiRandom); + hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), prrot.first, prrot.second); } } } else { if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, phiRandom); + hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), pr.first, pr.second); } } } - // } } void processSE(EventCandidates::iterator const& collision, TrackCandidates const& /*tracks*/, V0TrackCandidate const& V0s) @@ -1605,33 +1600,49 @@ struct HigherMassResonances { if (static_cast(passKs.size()) == config.noOfDaughters) { lResonanceGen = ROOT::Math::PxPyPzEVector(mcParticle.pt(), mcParticle.eta(), mcParticle.phi(), mcParticle.e()); lResonanceGen1 = daughter1 + daughter2; + mother2 = ROOT::Math::PxPyPzMVector(lResonanceGen1.Px(), lResonanceGen1.Py(), lResonanceGen1.Pz(), lResonanceGen1.M()); + mother3 = ROOT::Math::PxPyPzMVector(lResonanceGen.Px(), lResonanceGen.Py(), lResonanceGen.Pz(), lResonanceGen.M()); + + std::pair angularVar, angularVar1; + if (config.activateHelicityFrame) { + angularVar = getCosThetaPhi(mother3, daughter1, kHelicity); + angularVar1 = getCosThetaPhi(mother2, daughter1, kHelicity); + } else if (config.activateCollinsSoperFrame) { + angularVar = getCosThetaPhi(mother3, daughter1, kCollinsSoper); + angularVar1 = getCosThetaPhi(mother2, daughter1, kCollinsSoper); + } else if (config.activateProductionFrame) { + angularVar = getCosThetaPhi(mother3, daughter1, kProduction); + angularVar1 = getCosThetaPhi(mother2, daughter1, kProduction); + } else if (config.activateBeamAxisFrame) { + angularVar = getCosThetaPhi(mother3, daughter1, kBeamAxis); + angularVar1 = getCosThetaPhi(mother2, daughter1, kBeamAxis); + } else if (config.activateRandomFrame) { + angularVar = getCosThetaPhi(mother3, daughter1, kRandom); + angularVar1 = getCosThetaPhi(mother2, daughter1, kRandom); + } else { + angularVar = getCosThetaPhi(mother3, daughter1, kHelicity); + angularVar1 = getCosThetaPhi(mother2, daughter1, kHelicity); + } - ROOT::Math::Boost boost{lResonanceGen.BoostToCM()}; - ROOT::Math::Boost boost1{lResonanceGen1.BoostToCM()}; - - fourVecDauCM = boost(daughter1); - fourVecDauCM1 = boost1(daughter1); - - auto helicityGen = lResonanceGen.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(lResonanceGen.Vect().Mag2())); - auto helicityGen1 = lResonanceGen1.Vect().Dot(fourVecDauCM1.Vect()) / (std::sqrt(fourVecDauCM1.Vect().Mag2()) * std::sqrt(lResonanceGen1.Vect().Mag2())); - - hMChists.fill(HIST("Genf1710"), multiplicityGen, lResonanceGen.pt(), helicityGen); - hMChists.fill(HIST("Genf1710Calib"), genMultiplicity, lResonanceGen.pt(), helicityGen); + hMChists.fill(HIST("Genf1710"), multiplicityGen, lResonanceGen.pt(), angularVar.first); + hMChists.fill(HIST("Genf1710Calib"), genMultiplicity, lResonanceGen.pt(), angularVar.first); hMChists.fill(HIST("Genf1710_mass"), lResonanceGen.M()); hMChists.fill(HIST("GenRapidity"), mcParticle.y()); hMChists.fill(HIST("GenEta"), mcParticle.eta()); hMChists.fill(HIST("GenPhi"), mcParticle.phi()); + hMChists.fill(HIST("Gen1710PWA"), lResonanceGen.M(), angularVar.first, angularVar.second); if (config.isapplyPairRapidityMC && std::abs(lResonanceGen1.Rapidity()) >= config.rapidityMotherData) { continue; } - hMChists.fill(HIST("Genf17102"), multiplicityGen, lResonanceGen1.pt(), helicityGen1); - hMChists.fill(HIST("Genf1710Calib2"), genMultiplicity, lResonanceGen1.pt(), helicityGen1); + hMChists.fill(HIST("Genf17102"), multiplicityGen, lResonanceGen1.pt(), angularVar1.first); + hMChists.fill(HIST("Genf1710Calib2"), genMultiplicity, lResonanceGen1.pt(), angularVar1.first); hMChists.fill(HIST("Genf1710_mass2"), lResonanceGen1.M()); hMChists.fill(HIST("GenRapidity2"), lResonanceGen1.Rapidity()); hMChists.fill(HIST("GenEta2"), lResonanceGen1.Eta()); hMChists.fill(HIST("GenPhi2"), lResonanceGen1.Phi()); + hMChists.fill(HIST("Gen1710PWA2"), lResonanceGen1.M(), angularVar1.first, angularVar1.second); } passKs.clear(); // clear the vector for the next iteration } @@ -1887,37 +1898,46 @@ struct HigherMassResonances { daughter2 = ROOT::Math::PxPyPzMVector(v02.px(), v02.py(), v02.pz(), o2::constants::physics::MassK0Short); mother = daughter1 + daughter2; mother1 = ROOT::Math::PxPyPzEVector(mothertrack1.px(), mothertrack1.py(), mothertrack1.pz(), mothertrack1.e()); + mother2 = ROOT::Math::PxPyPzMVector(mother1.Px(), mother1.Py(), mother1.Pz(), mother1.M()); + + std::pair angularVar, angularVar1; + if (config.activateHelicityFrame) { + angularVar = getCosThetaPhi(mother, daughter1, kHelicity); + angularVar1 = getCosThetaPhi(mother2, daughter1, kHelicity); + } else if (config.activateCollinsSoperFrame) { + angularVar = getCosThetaPhi(mother, daughter1, kCollinsSoper); + angularVar1 = getCosThetaPhi(mother2, daughter1, kCollinsSoper); + } else if (config.activateProductionFrame) { + angularVar = getCosThetaPhi(mother, daughter1, kProduction); + angularVar1 = getCosThetaPhi(mother2, daughter1, kProduction); + } else if (config.activateBeamAxisFrame) { + angularVar = getCosThetaPhi(mother, daughter1, kBeamAxis); + angularVar1 = getCosThetaPhi(mother2, daughter1, kBeamAxis); + } else if (config.activateRandomFrame) { + angularVar = getCosThetaPhi(mother, daughter1, kRandom); + angularVar1 = getCosThetaPhi(mother2, daughter1, kRandom); + } else { + angularVar = getCosThetaPhi(mother, daughter1, kHelicity); + angularVar1 = getCosThetaPhi(mother2, daughter1, kHelicity); + } - ROOT::Math::Boost boost{mother.BoostToCM()}; - ROOT::Math::Boost boost1{mother1.BoostToCM()}; - - fourVecDauCM = boost(daughter1); - fourVecDauCM1 = boost1(daughter1); - - auto helicityRec = mother.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(mother.Vect().Mag2())); - - auto helicityRec2 = mother1.Vect().Dot(fourVecDauCM1.Vect()) / (std::sqrt(fourVecDauCM1.Vect().Mag2()) * std::sqrt(mother1.Vect().Mag2())); - - // const double deltaMassRec = deltaM(mctrackv01.mK0Short(), mctrackv02.mK0Short()); - // if (deltaMassRec > config.cMaxDeltaM) { - // continue; - // } - - hMChists.fill(HIST("Recf1710_pt1"), multiplicity, mothertrack1.pt(), mother1.M(), helicityRec2); - hMChists.fill(HIST("Recf1710Calib_pt1"), multiplicityRec, mothertrack1.pt(), mother1.M(), helicityRec2); + hMChists.fill(HIST("Recf1710_pt1"), multiplicity, mothertrack1.pt(), mother1.M(), angularVar1.first); + hMChists.fill(HIST("Recf1710Calib_pt1"), multiplicityRec, mothertrack1.pt(), mother1.M(), angularVar1.first); hMChists.fill(HIST("RecRapidity"), mothertrack1.y()); hMChists.fill(HIST("RecPhi"), mothertrack1.phi()); hMChists.fill(HIST("RecEta"), mothertrack1.eta()); + hMChists.fill(HIST("Recf1710PWA_pt1"), mother1.M(), angularVar1.first, angularVar1.second); if (config.isapplyPairRapidityMC && std::abs(mother.Rapidity()) >= config.rapidityMotherData) { continue; } - hMChists.fill(HIST("Recf1710_pt2"), multiplicity, mother.Pt(), mother.M(), helicityRec); - hMChists.fill(HIST("Recf1710Calib_pt2"), multiplicityRec, mother.Pt(), mother.M(), helicityRec); + hMChists.fill(HIST("Recf1710_pt2"), multiplicity, mother.Pt(), mother.M(), angularVar.first); + hMChists.fill(HIST("Recf1710Calib_pt2"), multiplicityRec, mother.Pt(), mother.M(), angularVar.first); hMChists.fill(HIST("RecRapidity2"), mother.Rapidity()); hMChists.fill(HIST("RecPhi2"), mother.Phi()); hMChists.fill(HIST("RecEta2"), mother.Eta()); + hMChists.fill(HIST("Recf1710PWA_pt2"), mother.M(), angularVar.first, angularVar.second); } gindex2.clear(); }