diff --git a/PWGLF/Tasks/Resonances/higherMassResonances.cxx b/PWGLF/Tasks/Resonances/higherMassResonances.cxx index 792916cb6fd..1e8c7358f73 100644 --- a/PWGLF/Tasks/Resonances/higherMassResonances.cxx +++ b/PWGLF/Tasks/Resonances/higherMassResonances.cxx @@ -64,6 +64,7 @@ #include #include #include +#include #include using namespace o2; @@ -206,6 +207,12 @@ struct HigherMassResonances { float beamEnergy = 13600.0; double beamMomentum = std::sqrt(beamEnergy * beamEnergy / 4 - o2::constants::physics::MassProton * o2::constants::physics::MassProton); // GeV int noOfDaughters = 2; + float theta2 = 0.0f; + int detId = 0; + int refAId = 0; + int refBId = 0; + float minQvecAmp = 1e-5; + Service pdgDB{}; } config; // Service PDGdatabase; @@ -213,8 +220,7 @@ struct HigherMassResonances { // variables declaration float multiplicity = 0.0f; - float theta2; - ROOT::Math::PxPyPzMVector daughter1, daughter2, daughterRot, daughterRotCM, mother, motherRot, fourVecDauCM, fourVecDauCM1; + 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,10 +234,6 @@ struct HigherMassResonances { bool isMix = false; EventPlaneHelper helperEP; - int detId; - int refAId; - int refBId; - float minQvecAmp = 1e-5; void init(InitContext const&) { @@ -253,7 +255,7 @@ struct HigherMassResonances { AxisSpec axisEvtResPlQA = {102, -1.02, 1.02, ""}; // THnSparses - std::array sparses = {config.activateHelicityFrame, config.activateCollinsSoperFrame, config.activateProductionFrame, config.activateBeamAxisFrame, config.activateRandomFrame}; + std::array 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"); @@ -399,12 +401,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}); @@ -441,9 +449,9 @@ struct HigherMassResonances { } if (doprocessSEEP) { - detId = getdetId(config.cfgDetName); - refAId = getdetId(config.cfgRefAName); - refBId = getdetId(config.cfgRefBName); + config.detId = getdetId(config.cfgDetName); + config.refAId = getdetId(config.cfgRefAName); + config.refBId = getdetId(config.cfgRefBName); hglue.add("EpDet", "", {HistType::kTH2F, {axisCentQA, axisEvtPlQA}}); hglue.add("EpRefA", "", {HistType::kTH2F, {axisCentQA, axisEvtPlQA}}); @@ -460,21 +468,26 @@ struct HigherMassResonances { { if (name.value == "FT0C") { return 0; - } else if (name.value == "FT0A") { + } + if (name.value == "FT0A") { return 1; - } else if (name.value == "FT0M") { + } + if (name.value == "FT0M") { return 2; - } else if (name.value == "FV0A") { + } + if (name.value == "FV0A") { return 3; - } else if (name.value == "TPCPos") { + } + if (name.value == "TPCPos") { return 4; - } else if (name.value == "TPCNeg") { + } + if (name.value == "TPCNeg") { return 5; - } else if (name.value == "TPCTot") { + } + if (name.value == "TPCTot") { return 6; - } else { - return 0; } + return 0; } template @@ -490,8 +503,8 @@ struct HigherMassResonances { // if (config.isSel8 && !collision.sel8()) // return false; - if (fillHist) - rEventSelection.fill(HIST("hEventCut"), 2); + // if (fillHist) + // rEventSelection.fill(HIST("hEventCut"), 2); if (config.isNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) return false; @@ -505,8 +518,8 @@ struct HigherMassResonances { // if (config.isNoSameBunchPileup && (!collision.selection_bit(aod::evsel::kNoSameBunchPileup))) // return false; - if (fillHist) - rEventSelection.fill(HIST("hEventCut"), 5); + // if (fillHist) + // rEventSelection.fill(HIST("hEventCut"), 5); if (config.isAllLayersGoodITS && !collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) return false; @@ -518,8 +531,8 @@ struct HigherMassResonances { // if (config.isApplyOccCut && (std::abs(collision.trackOccupancyInTimeRange()) > config.configOccCut)) // return false; - if (fillHist) - rEventSelection.fill(HIST("hEventCut"), 7); + // if (fillHist) + // rEventSelection.fill(HIST("hEventCut"), 7); if (rctCut.requireRCTFlagChecker && !rctChecker(collision)) return false; @@ -570,8 +583,8 @@ struct HigherMassResonances { // if (config.isSel8 && !collision.sel8()) // return false; - if (fillHist) - rEventSelection.fill(HIST("hEventCut"), 2); + // if (fillHist) + // rEventSelection.fill(HIST("hEventCut"), 2); if (config.isNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) return false; @@ -585,8 +598,8 @@ struct HigherMassResonances { // if (config.isNoSameBunchPileup && (!collision.selection_bit(aod::evsel::kNoSameBunchPileup))) // return false; - if (fillHist) - rEventSelection.fill(HIST("hEventCut"), 5); + // if (fillHist) + // rEventSelection.fill(HIST("hEventCut"), 5); if (config.isAllLayersGoodITS && !collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) return false; @@ -598,8 +611,8 @@ struct HigherMassResonances { // if (config.isApplyOccCut && (std::abs(collision.trackOccupancyInTimeRange()) > config.configOccCut)) // return false; - if (fillHist) - rEventSelection.fill(HIST("hEventCut"), 7); + // if (fillHist) + // rEventSelection.fill(HIST("hEventCut"), 7); if (rctCut.requireRCTFlagChecker && !rctChecker(collision)) return false; @@ -783,6 +796,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 zaxisHELocal = ROOT::Math::XYZVectorF(motherVec.Vect()).Unit(); + auto yaxisHELocal = ROOT::Math::XYZVectorF(beam1cm.Cross(beam2cm)).Unit(); + auto xaxisHELocal = ROOT::Math::XYZVectorF(yaxisHELocal.Cross(zaxisHELocal)).Unit(); + cosTheta = motherVec.Vect().Dot(dauCM.Vect()) / (std::sqrt(dauCM.Vect().Mag2()) * std::sqrt(motherVec.Vect().Mag2())); + phi = std::atan2(yaxisHELocal.Dot(v1cm), xaxisHELocal.Dot(v1cm)); + phi = RecoDecay::constrainAngle(phi, 0.0); + } else if (frame == kCollinsSoper) { + auto zAxisCSLocal = ROOT::Math::XYZVectorF((beam1cm.Unit() - beam2cm.Unit())).Unit(); + auto yAxisCSLocal = ROOT::Math::XYZVectorF(beam1cm.Cross(beam2cm)).Unit(); + auto xAxisCSLocal = ROOT::Math::XYZVectorF(yAxisCSLocal.Cross(zAxisCSLocal)).Unit(); + cosTheta = zAxisCSLocal.Dot(v1cm); + phi = std::atan2(yAxisCSLocal.Dot(v1cm), xAxisCSLocal.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 zaxisHELocal = ROOT::Math::XYZVectorF(motherVec.Vect()).Unit(); + auto yaxisHELocal = ROOT::Math::XYZVectorF(beam1cm.Cross(beam2cm)).Unit(); + auto xaxisHELocal = ROOT::Math::XYZVectorF(yaxisHELocal.Cross(zaxisHELocal)).Unit(); + phi = std::atan2(yaxisHELocal.Dot(v1cm), xaxisHELocal.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 zaxisHELocal = ROOT::Math::XYZVectorF(motherVec.Vect()).Unit(); + auto yaxisHELocal = ROOT::Math::XYZVectorF(beam1cm.Cross(beam2cm)).Unit(); + auto xaxisHELocal = ROOT::Math::XYZVectorF(yaxisHELocal.Cross(zaxisHELocal)).Unit(); + phi = std::atan2(yaxisHELocal.Dot(v1cm), xaxisHELocal.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; @@ -866,185 +939,114 @@ struct HigherMassResonances { // zBeam direction in lab frame template - void fillInvMass(const T& mother, float multiplicity, const T& daughter1, const T& daughter2, bool isMix) + void fillInvMass(const T& motherVec, float eventMultiplicity, const T& daughterVec1, const T& daughterVec2, bool isMixed) { - // //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())); - if (!isMix) { - if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity, anglePhi); + auto pr = getCosThetaPhi(motherVec, daughterVec1, kHelicity); + if (!isMixed) { + if (std::abs(motherVec.Rapidity()) < config.rapidityMotherData) { + hglue.fill(HIST("h3glueInvMassDS"), eventMultiplicity, motherVec.Pt(), motherVec.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); + config.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()); + daughterRot = ROOT::Math::PxPyPzMVector(daughterVec1.Px() * std::cos(config.theta2) - daughterVec1.Py() * std::sin(config.theta2), daughterVec1.Px() * std::sin(config.theta2) + daughterVec1.Py() * std::cos(config.theta2), daughterVec1.Pz(), daughterVec1.M()); - motherRot = daughterRot + daughter2; - - ROOT::Math::Boost boost2{motherRot.BoostToCM()}; - daughterRotCM = boost2(daughterRot); + motherRot = daughterRot + daughterVec2; - 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"), eventMultiplicity, 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); + if (std::abs(motherVec.Rapidity()) < config.rapidityMotherData) { + hglue.fill(HIST("h3glueInvMassME"), eventMultiplicity, motherVec.Pt(), motherVec.M(), pr.first, pr.second); } } } else if (config.activateCollinsSoperFrame) { - if (!isMix) { - if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarCS, phiCS); + auto pr = getCosThetaPhi(motherVec, daughterVec1, kCollinsSoper); + if (!isMixed) { + if (std::abs(motherVec.Rapidity()) < config.rapidityMotherData) { + hglue.fill(HIST("h3glueInvMassDS"), eventMultiplicity, motherVec.Pt(), motherVec.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); - - 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()); + config.theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut); - motherRot = daughterRot + daughter2; + daughterRot = ROOT::Math::PxPyPzMVector(daughterVec1.Px() * std::cos(config.theta2) - daughterVec1.Py() * std::sin(config.theta2), daughterVec1.Px() * std::sin(config.theta2) + daughterVec1.Py() * std::cos(config.theta2), daughterVec1.Pz(), daughterVec1.M()); - ROOT::Math::Boost boost2{motherRot.BoostToCM()}; - daughterRotCM = boost2(daughterRot); + motherRot = daughterRot + daughterVec2; - 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"), eventMultiplicity, 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); + if (std::abs(motherVec.Rapidity()) < config.rapidityMotherData) { + hglue.fill(HIST("h3glueInvMassME"), eventMultiplicity, motherVec.Pt(), motherVec.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())); - if (!isMix) { - if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaProduction, anglePhi); + auto pr = getCosThetaPhi(motherVec, daughterVec1, kProduction); + if (!isMixed) { + if (std::abs(motherVec.Rapidity()) < config.rapidityMotherData) { + hglue.fill(HIST("h3glueInvMassDS"), eventMultiplicity, motherVec.Pt(), motherVec.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()); + config.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(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()); + auto prrot = getCosThetaPhi(motherRot, daughterVec1, kProduction); if (std::abs(motherRot.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaProduction, anglePhi); + hglue.fill(HIST("h3glueInvMassRot"), eventMultiplicity, 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); + if (std::abs(motherVec.Rapidity()) < config.rapidityMotherData) { + hglue.fill(HIST("h3glueInvMassME"), eventMultiplicity, motherVec.Pt(), motherVec.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()); - if (!isMix) { - if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarBeam, anglePhi); + auto pr = getCosThetaPhi(motherVec, daughterVec1, kBeamAxis); + if (!isMixed) { + if (std::abs(motherVec.Rapidity()) < config.rapidityMotherData) { + hglue.fill(HIST("h3glueInvMassDS"), eventMultiplicity, motherVec.Pt(), motherVec.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()); + config.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(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()); + auto prrot = getCosThetaPhi(motherRot, daughterVec1, kBeamAxis); if (std::abs(motherRot.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarBeam, anglePhi); + hglue.fill(HIST("h3glueInvMassRot"), eventMultiplicity, 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); + if (std::abs(motherVec.Rapidity()) < config.rapidityMotherData) { + hglue.fill(HIST("h3glueInvMassME"), eventMultiplicity, motherVec.Pt(), motherVec.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()); - if (!isMix) { - if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, phiRandom); + auto pr = getCosThetaPhi(motherVec, daughterVec1, kRandom); + if (!isMixed) { + if (std::abs(motherVec.Rapidity()) < config.rapidityMotherData) { + hglue.fill(HIST("h3glueInvMassDS"), eventMultiplicity, motherVec.Pt(), motherVec.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()); + config.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(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()); + auto prrot = getCosThetaPhi(motherRot, daughterVec1, kRandom); if (std::abs(motherRot.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarRandom, phiRandom); + hglue.fill(HIST("h3glueInvMassRot"), eventMultiplicity, 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); + if (std::abs(motherVec.Rapidity()) < config.rapidityMotherData) { + hglue.fill(HIST("h3glueInvMassME"), eventMultiplicity, motherVec.Pt(), motherVec.M(), pr.first, pr.second); } } } - // } } void processSE(EventCandidates::iterator const& collision, TrackCandidates const& /*tracks*/, V0TrackCandidate const& V0s) @@ -1059,8 +1061,6 @@ struct HigherMassResonances { multiplicity = collision.centFT0C(); } else if (config.cSelectMultEstimator == kFV0A) { multiplicity = collision.centFV0A(); - } else { - multiplicity = collision.centFT0M(); // default } if (!selectionEvent(collision, true)) { @@ -1085,7 +1085,7 @@ struct HigherMassResonances { } std::vector v0indexes; - bool allConditionsMet = 0; + bool allConditionsMet = false; for (const auto& [v1, v2] : combinations(CombinationsFullIndexPolicy(V0s, V0s))) { @@ -1147,7 +1147,7 @@ struct HigherMassResonances { continue; } - allConditionsMet = 1; + allConditionsMet = true; daughter1 = ROOT::Math::PxPyPzMVector(v1.px(), v1.py(), v1.pz(), o2::constants::physics::MassK0Short); // Kshort daughter2 = ROOT::Math::PxPyPzMVector(v2.px(), v2.py(), v2.pz(), o2::constants::physics::MassK0Short); // Kshort @@ -1193,9 +1193,9 @@ struct HigherMassResonances { } for (int i = 0; i < config.cRotations; i++) { - double theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut); + config.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()); + daughterRot = ROOT::Math::PxPyPzMVector(daughter1.Px() * std::cos(config.theta2) - daughter1.Py() * std::sin(config.theta2), daughter1.Px() * std::sin(config.theta2) + daughter1.Py() * std::cos(config.theta2), daughter1.Pz(), daughter1.M()); motherRot = daughterRot + daughter2; @@ -1228,15 +1228,13 @@ struct HigherMassResonances { multiplicity = collision.centFT0C(); } else if (config.cSelectMultEstimator == kFV0A) { multiplicity = collision.centFV0A(); - } else { - multiplicity = collision.centFT0C(); // default } if (!selectionEvent(collision, true)) { return; } - if (collision.qvecAmp()[detId] < minQvecAmp || collision.qvecAmp()[refAId] < minQvecAmp || collision.qvecAmp()[refBId] < minQvecAmp) { + if (collision.qvecAmp()[config.detId] < config.minQvecAmp || collision.qvecAmp()[config.refAId] < config.minQvecAmp || collision.qvecAmp()[config.refBId] < config.minQvecAmp) { return; } @@ -1245,10 +1243,11 @@ struct HigherMassResonances { rEventSelection.fill(HIST("hmultiplicity"), multiplicity); } - float eps[3] = {0.}; - eps[0] = helperEP.GetEventPlane(collision.qvecRe()[4 * detId + 3], collision.qvecIm()[4 * detId + 3], 2); - eps[1] = helperEP.GetEventPlane(collision.qvecRe()[4 * refAId + 3], collision.qvecIm()[4 * refAId + 3], 2); - eps[2] = helperEP.GetEventPlane(collision.qvecRe()[4 * refBId + 3], collision.qvecIm()[4 * refBId + 3], 2); + std::array eps{}; + + eps[0] = helperEP.GetEventPlane(collision.qvecRe()[4 * config.detId + 3], collision.qvecIm()[4 * config.detId + 3], 2); + eps[1] = helperEP.GetEventPlane(collision.qvecRe()[4 * config.refAId + 3], collision.qvecIm()[4 * config.refAId + 3], 2); + eps[2] = helperEP.GetEventPlane(collision.qvecRe()[4 * config.refBId + 3], collision.qvecIm()[4 * config.refBId + 3], 2); float resNumA = helperEP.GetResolution(eps[0], eps[1], 2); float resNumB = helperEP.GetResolution(eps[0], eps[2], 2); @@ -1361,9 +1360,9 @@ struct HigherMassResonances { } for (int i = 0; i < config.cRotations; i++) { - double theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut); + config.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()); + daughterRot = ROOT::Math::PxPyPzMVector(daughter1.Px() * std::cos(config.theta2) - daughter1.Py() * std::sin(config.theta2), daughter1.Px() * std::sin(config.theta2) + daughter1.Py() * std::cos(config.theta2), daughter1.Pz(), daughter1.M()); motherRot = daughterRot + daughter2; @@ -1404,7 +1403,7 @@ struct HigherMassResonances { SameKindPair pair3{binningOnFT0C, config.cfgNmixedEvents, -1, collisions, tracksTuple, &cache}; SameKindPair pair4{binningOnFV0A, config.cfgNmixedEvents, -1, collisions, tracksTuple, &cache}; - auto runMixing = [&](auto& pair, auto multiplicityGetter) { + auto runMixing = [&](const auto& pair, auto multiplicityGetter) { for (const auto& [c1, tracks1, c2, tracks2] : pair) { multiplicity = multiplicityGetter(c1); @@ -1481,11 +1480,11 @@ struct HigherMassResonances { if (config.qAOptimisation) { double deltaRvalue = std::sqrt(TVector2::Phi_mpi_pi(daughter1.phi() - daughter2.phi()) * TVector2::Phi_mpi_pi(daughter1.phi() - daughter2.phi()) + (daughter1.eta() - daughter2.eta()) * (daughter1.eta() - daughter2.eta())); - const double deltaMass = deltaM(t1.mK0Short(), t2.mK0Short()); + const double deltaMassOpt = deltaM(t1.mK0Short(), t2.mK0Short()); // const double ptCorr = std::abs(daughter1.Pt() + daughter2.Pt()) / mother.Pt(); const double ptCorr = (mother.Pt() - daughter1.Pt() != 0.) ? daughter1.Pt() / (mother.Pt() - daughter1.Pt()) : 0.; if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), deltaMass, deltaRvalue, ptCorr); + hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), deltaMassOpt, deltaRvalue, ptCorr); } } } @@ -1509,7 +1508,6 @@ struct HigherMassResonances { std::vector passKs; ROOT::Math::PxPyPzMVector lResonanceGen1; ROOT::Math::PxPyPzEVector lResonanceGen; - Service pdgDB; // void processGen(aod::McCollision const& mcCollision, aod::McParticles const& mcParticles, const soa::SmallGroups& collisions) void processGen(EventMCGenerated::iterator const& mcCollision, aod::McParticles const& mcParticles, const soa::SmallGroups& collisions) @@ -1521,7 +1519,7 @@ struct HigherMassResonances { bool isINELgt0true = false; - if (pwglf::isINELgtNmc(mcParticles, 0, pdgDB)) { + if (pwglf::isINELgtNmc(mcParticles, 0, config.pdgDB)) { isINELgt0true = true; } if (config.isINELgt0 && !isINELgt0true) { @@ -1544,8 +1542,6 @@ struct HigherMassResonances { multiplicityGen = collision.centFT0C(); } else if (config.cSelectMultEstimator == kFV0A) { multiplicityGen = collision.centFV0A(); - } else { - multiplicityGen = collision.centFT0M(); // default } if (!selectionEvent(collision, true)) { @@ -1605,33 +1601,47 @@ 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()); + + int frame = -1; + + if (config.activateHelicityFrame) { + frame = kHelicity; + } else if (config.activateCollinsSoperFrame) { + frame = kCollinsSoper; + } else if (config.activateProductionFrame) { + frame = kProduction; + } else if (config.activateBeamAxisFrame) { + frame = kBeamAxis; + } else if (config.activateRandomFrame) { + frame = kRandom; + } else { + LOGP(fatal, "No reference frame selected"); + } - 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())); + auto angularVar = getCosThetaPhi(mother3, daughter1, frame); + auto angularVar1 = getCosThetaPhi(mother2, daughter1, frame); - 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 } @@ -1662,8 +1672,6 @@ struct HigherMassResonances { multiplicity1 = RecCollision.centFT0C(); } else if (config.cSelectMultEstimator == kFV0A) { multiplicity1 = RecCollision.centFV0A(); - } else { - multiplicity1 = RecCollision.centFT0M(); // default } isSelectedEvent = true; @@ -1671,7 +1679,7 @@ struct HigherMassResonances { bool isINELgt0true = false; - if (pwglf::isINELgtNmc(mcParticles, 0, pdgDB)) { + if (pwglf::isINELgtNmc(mcParticles, 0, config.pdgDB)) { isINELgt0true = true; } if (config.isINELgt0 && !isINELgt0true) { @@ -1682,15 +1690,15 @@ struct HigherMassResonances { return; } - auto multiplicityGen = -1; - multiplicityGen = mcCollision.centFT0M(); - hMChists.fill(HIST("MCcorrections/MultiplicityGen"), multiplicityGen); + auto multiplicityGenEvt = -1.f; + multiplicityGenEvt = mcCollision.centFT0M(); + hMChists.fill(HIST("MCcorrections/MultiplicityGen"), multiplicityGenEvt); hMChists.fill(HIST("MCcorrections/MultiplicityGen2"), multiplicity1); hMChists.fill(HIST("MCcorrections/hMultvsCent"), multiplicity1, multMC); // Event loss if (isSelectedEvent) { - hMChists.fill(HIST("MCcorrections/MultiplicityRec"), multiplicityGen); + hMChists.fill(HIST("MCcorrections/MultiplicityRec"), multiplicityGenEvt); hMChists.fill(HIST("MCcorrections/MultiplicityRec2"), multiplicity1); } @@ -1699,11 +1707,11 @@ struct HigherMassResonances { if (std::abs(mcPart.y()) >= config.rapidityMotherData || std::abs(mcPart.pdgCode()) != config.pdgCodes[config.selectMCparticles]) continue; - hMChists.fill(HIST("MCcorrections/hSignalLossDenominator"), mcPart.pt(), multiplicityGen); + hMChists.fill(HIST("MCcorrections/hSignalLossDenominator"), mcPart.pt(), multiplicityGenEvt); hMChists.fill(HIST("MCcorrections/hSignalLossDenominator2"), mcPart.pt(), multiplicity1); hMChists.fill(HIST("MCcorrections/hSignalLossDenominator3"), mcPart.pt(), multMC); if (isSelectedEvent) { - hMChists.fill(HIST("MCcorrections/hSignalLossNumerator"), mcPart.pt(), multiplicityGen); + hMChists.fill(HIST("MCcorrections/hSignalLossNumerator"), mcPart.pt(), multiplicityGenEvt); hMChists.fill(HIST("MCcorrections/hSignalLossNumerator2"), mcPart.pt(), multiplicity1); hMChists.fill(HIST("MCcorrections/hSignalLossNumerator3"), mcPart.pt(), multMC); } @@ -1753,26 +1761,24 @@ struct HigherMassResonances { const auto& mcCollisionRec = collision.mcCollision_as(); multiplicityRec = mcCollisionRec.centFT0M(); - auto multiplicity = -999.0; + auto eventMultiplicity = -999.0; if (config.cSelectMultEstimator == kFT0M) { - multiplicity = collision.centFT0M(); + eventMultiplicity = collision.centFT0M(); } else if (config.cSelectMultEstimator == kFT0A) { - multiplicity = collision.centFT0A(); + eventMultiplicity = collision.centFT0A(); } else if (config.cSelectMultEstimator == kFT0C) { - multiplicity = collision.centFT0C(); + eventMultiplicity = collision.centFT0C(); } else if (config.cSelectMultEstimator == kFV0A) { - multiplicity = collision.centFV0A(); - } else { - multiplicity = collision.centFT0M(); // default + eventMultiplicity = collision.centFV0A(); } if (!selectionEvent(collision, false)) { return; } - hMChists.fill(HIST("Rec_Multiplicity"), multiplicity); + hMChists.fill(HIST("Rec_Multiplicity"), eventMultiplicity); rEventSelection.fill(HIST("hVertexZRec"), collision.posZ()); - hMChists.fill(HIST("MC_mult_after_event_sel"), multiplicity); + hMChists.fill(HIST("MC_mult_after_event_sel"), eventMultiplicity); eventCounter++; for (const auto& v01 : V0s) { @@ -1800,21 +1806,23 @@ struct HigherMassResonances { if (!negtrack1.has_mcParticle() || !negtrack2.has_mcParticle()) continue; - double nTPCSigmaPos1[1]{postrack1.tpcNSigmaPi()}; - double nTPCSigmaNeg1[1]{negtrack1.tpcNSigmaPi()}; - double nTPCSigmaPos2[1]{postrack2.tpcNSigmaPi()}; - double nTPCSigmaNeg2[1]{negtrack2.tpcNSigmaPi()}; + double nTPCSigmaPos1 = postrack1.tpcNSigmaPi(); + double nTPCSigmaNeg1 = negtrack1.tpcNSigmaPi(); + double nTPCSigmaPos2 = postrack2.tpcNSigmaPi(); + double nTPCSigmaNeg2 = negtrack2.tpcNSigmaPi(); - if (!isSelectedV0Daughter(postrack1, 1, nTPCSigmaPos1[0], v01) || !isSelectedV0Daughter(postrack2, 1, nTPCSigmaPos2[0], v02)) { + if (!isSelectedV0Daughter(postrack1, 1, nTPCSigmaPos1, v01) || + !isSelectedV0Daughter(postrack2, 1, nTPCSigmaPos2, v02)) { continue; } - if (!isSelectedV0Daughter(negtrack1, -1, nTPCSigmaNeg1[0], v01) || !isSelectedV0Daughter(negtrack2, -1, nTPCSigmaNeg2[0], v02)) { + if (!isSelectedV0Daughter(negtrack1, -1, nTPCSigmaNeg1, v01) || + !isSelectedV0Daughter(negtrack2, -1, nTPCSigmaNeg2, v02)) { continue; } hMChists.fill(HIST("events_checkrec"), 1.5); - if (!selectionV0(collision, v01, multiplicity) || !selectionV0(collision, v02, multiplicity)) { + if (!selectionV0(collision, v01, eventMultiplicity) || !selectionV0(collision, v02, eventMultiplicity)) { continue; } hMChists.fill(HIST("events_checkrec"), 2.5); @@ -1887,37 +1895,40 @@ 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()); + + int frame = kHelicity; + + if (config.activateCollinsSoperFrame) { + frame = kCollinsSoper; + } else if (config.activateProductionFrame) { + frame = kProduction; + } else if (config.activateBeamAxisFrame) { + frame = kBeamAxis; + } else if (config.activateRandomFrame) { + frame = kRandom; + } - 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; - // } + auto angularVar = getCosThetaPhi(mother, daughter1, frame); + auto angularVar1 = getCosThetaPhi(mother2, daughter1, frame); - 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"), eventMultiplicity, 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"), eventMultiplicity, 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(); } @@ -1938,8 +1949,6 @@ struct HigherMassResonances { multiplicity = collision.centFT0C(); } else if (config.cSelectMultEstimator == kFV0A) { multiplicity = collision.centFV0A(); - } else { - multiplicity = collision.centFT0M(); // default } if (!selectionEventDerived(collision, true)) { @@ -1952,7 +1961,7 @@ struct HigherMassResonances { } std::vector v0indexes; - bool allConditionsMet = 0; + bool allConditionsMet = false; for (const auto& [v1, v2] : combinations(CombinationsFullIndexPolicy(V0s, V0s))) { @@ -1979,7 +1988,7 @@ struct HigherMassResonances { continue; } - allConditionsMet = 1; + allConditionsMet = true; daughter1 = ROOT::Math::PxPyPzMVector(v1.px(), v1.py(), v1.pz(), o2::constants::physics::MassK0Short); // Kshort daughter2 = ROOT::Math::PxPyPzMVector(v2.px(), v2.py(), v2.pz(), o2::constants::physics::MassK0Short); // Kshort @@ -2022,9 +2031,9 @@ struct HigherMassResonances { } for (int i = 0; i < config.cRotations; i++) { - double theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / config.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / config.rotationalCut); + config.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()); + daughterRot = ROOT::Math::PxPyPzMVector(daughter1.Px() * std::cos(config.theta2) - daughter1.Py() * std::sin(config.theta2), daughter1.Px() * std::sin(config.theta2) + daughter1.Py() * std::cos(config.theta2), daughter1.Pz(), daughter1.M()); motherRot = daughterRot + daughter2; // double pTcorrRot = std::abs(daughterRot.Pt() + daughter2.Pt()) / motherRot.Pt(); @@ -2095,11 +2104,11 @@ struct HigherMassResonances { if (config.qAOptimisation) { double deltaRvalue = std::sqrt(TVector2::Phi_mpi_pi(daughter1.phi() - daughter2.phi()) * TVector2::Phi_mpi_pi(daughter1.phi() - daughter2.phi()) + (daughter1.eta() - daughter2.eta()) * (daughter1.eta() - daughter2.eta())); - const double deltaMass = deltaM(t1.mK0Short(), t2.mK0Short()); + const double deltaMassOpt = deltaM(t1.mK0Short(), t2.mK0Short()); // const double ptCorr = std::abs(daughter1.Pt() + daughter2.Pt()) / mother.Pt(); const double ptCorr = (mother.Pt() - daughter1.Pt() != 0.) ? daughter1.Pt() / (mother.Pt() - daughter1.Pt()) : 0.; if (std::abs(mother.Rapidity()) < config.rapidityMotherData) { - hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), deltaMass, deltaRvalue, ptCorr); + hglue.fill(HIST("h3glueInvMassME"), multiplicity, mother.Pt(), mother.M(), deltaMassOpt, deltaRvalue, ptCorr); } } }