diff --git a/PWGLF/Tasks/Resonances/higherMassResonances.cxx b/PWGLF/Tasks/Resonances/higherMassResonances.cxx index 4c18869d78d..36491e50864 100644 --- a/PWGLF/Tasks/Resonances/higherMassResonances.cxx +++ b/PWGLF/Tasks/Resonances/higherMassResonances.cxx @@ -157,6 +157,7 @@ struct HigherMassResonances { Configurable cfgETAcut{"cfgETAcut", 0.8f, "Track ETA cut"}; Configurable deltaRDaugherCut{"deltaRDaugherCut", 0.001f, "DeltaR cut on V0 daughters"}; Configurable deltaRK0sCut{"deltaRK0sCut", 0.1f, "Apply deltaR cut between two K0s"}; + Configurable cfgArmenterosCut{"cfgArmenterosCut", 0.2f, "Armenteros-Podolanski cut"}; // Configurable for track selection and multiplicity Configurable cfgPTcut{"cfgPTcut", 0.2f, "Track PT cut"}; @@ -174,7 +175,7 @@ struct HigherMassResonances { Configurable activateHelicityFrame{"activateHelicityFrame", false, "Activate the THnSparse with cosThStar w.r.t. helicity axis"}; Configurable activateCollinsSoperFrame{"activateCollinsSoperFrame", false, "Activate the THnSparse with cosThStar w.r.t. Collins soper axis"}; Configurable activateProductionFrame{"activateProductionFrame", false, "Activate the THnSparse with cosThStar w.r.t. production axis"}; - Configurable activateBeamAxisFrame{"activateBeamAxisFrame", true, "Activate the THnSparse with cosThStar w.r.t. beam axis (Gottified jackson frame)"}; + Configurable activateGJFrame{"activateGJFrame", true, "Activate the THnSparse with cosThStar w.r.t. beam axis (Gottified jackson frame)"}; Configurable activateRandomFrame{"activateRandomFrame", false, "Activate the THnSparse with cosThStar w.r.t. random axis"}; Configurable cRotations{"cRotations", 3, "Number of random rotations in the rotational background"}; @@ -212,6 +213,7 @@ struct HigherMassResonances { int refAId = 0; int refBId = 0; float minQvecAmp = 1e-5; + double tolerance = 1e-12; } config; // Service PDGdatabase; @@ -255,7 +257,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.activateGJFrame, config.activateRandomFrame}; if (std::accumulate(sparses.begin(), sparses.end(), 0) == 0) { LOGP(fatal, "No output THnSparses enabled"); @@ -269,7 +271,7 @@ struct HigherMassResonances { if (config.activateProductionFrame) { LOGP(info, "THnSparse with cosThStar w.r.t. production axis active."); } - if (config.activateBeamAxisFrame) { + if (config.activateGJFrame) { LOGP(info, "THnSparse with cosThStar w.r.t. beam axis active. (Gottified jackson frame)"); } if (config.activateRandomFrame) { @@ -311,10 +313,11 @@ struct HigherMassResonances { hv0label->GetXaxis()->SetBinLabel(5, "Daughter DCA"); hv0label->GetXaxis()->SetBinLabel(6, "CosPA"); hv0label->GetXaxis()->SetBinLabel(7, "Decay Radius"); - hv0label->GetXaxis()->SetBinLabel(8, "Lifetime"); - hv0label->GetXaxis()->SetBinLabel(9, "CompetingCascade"); - hv0label->GetXaxis()->SetBinLabel(10, "Standard V0"); - hv0label->GetXaxis()->SetBinLabel(11, "Mass Tolerance"); + hv0label->GetXaxis()->SetBinLabel(8, "Armenteros-Podolanski"); + hv0label->GetXaxis()->SetBinLabel(9, "Lifetime"); + hv0label->GetXaxis()->SetBinLabel(10, "CompetingCascade"); + hv0label->GetXaxis()->SetBinLabel(11, "Standard V0"); + hv0label->GetXaxis()->SetBinLabel(12, "Mass Tolerance"); std::shared_ptr hv0DauLabel = rEventSelection.get(HIST("htrackscheck_v0_daughters")); hv0DauLabel->GetXaxis()->SetBinLabel(1, "AllDau Tracks"); @@ -385,6 +388,7 @@ struct HigherMassResonances { rKzeroShort.add("hLT", "hLT", {HistType::kTH1F, {{100, 0.0f, 50.0f}}}); rKzeroShort.add("angularSeparation", "Angular distribution between two K0s vs pT", {HistType::kTH1F, {{200, 0.0f, 4.0f}}}); rKzeroShort.add("hDauDeltaR", "Delta R of positive and negative daughers", {HistType::kTHnSparseF, {angleSepAxis, angleSepAxis}}); + rKzeroShort.add("hArmenterosPodolanski", "Armenteros-Podolanski plot", HistType::kTH2D, {AxisSpec{100, -1, 1, "#alpha"}, {200, 0, 0.5, "qtArm"}}); } rKzeroShort.add("NksProduced", "Number of K0s produced", kTH1I, {{15, -0.5, 14.5}}); @@ -405,6 +409,8 @@ struct HigherMassResonances { 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("GenThetavsPhi1", "GenThetavsPhi1", kTH2F, {{thnAxisPOL}, {thnAxisPhi}}); + hMChists.add("GenThetavsPhi2", "GenThetavsPhi2", kTH2F, {{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}); @@ -412,6 +418,8 @@ struct HigherMassResonances { 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("RecThetavsPhi1", "RecTheta vs Phi1", kTH2F, {{thnAxisPOL}, {thnAxisPhi}}); + hMChists.add("RecThetavsPhi2", "RecTheta vs Phi2", kTH2F, {{thnAxisPOL}, {thnAxisPhi}}); hMChists.add("h1Recsplit", "Rec p_{T}2", kTH1F, {ptAxis}); hMChists.add("Genf1710_mass", "Gen f_{0}(1710) mass", kTH1F, {glueballMassAxis}); @@ -630,9 +638,6 @@ struct HigherMassResonances { template bool selectionV0(Collision const& collision, V0 const& candidate, float /*multiplicity*/) { - // const float qtarm = candidate.qtarm(); - // const float alph = candidate.alpha(); - // float arm = qtarm / alph; const float pT = candidate.pt(); const float tranRad = candidate.v0radius(); const float dcaDaughv0 = candidate.dcaV0daughters(); @@ -647,6 +652,7 @@ struct HigherMassResonances { rKzeroShort.fill(HIST("hLT"), ctauK0s); rKzeroShort.fill(HIST("hDCAV0Daughters"), candidate.dcaV0daughters()); rKzeroShort.fill(HIST("hV0CosPA"), candidate.v0cosPA()); + rKzeroShort.fill(HIST("hArmenterosPodolanski"), candidate.alpha(), candidate.qtarm()); } if (config.qAcorrelation2Dhist) rKzeroShort.fill(HIST("mass_lambda_kshort_before"), candidate.mK0Short(), candidate.mLambda()); @@ -683,9 +689,9 @@ struct HigherMassResonances { } rEventSelection.fill(HIST("htrackscheck_v0"), 6.5); - // if (tranRad > config.confV0TranRadV0Max) { - // return false; - // } + if (candidate.qtarm() < config.cfgArmenterosCut * std::fabs(candidate.alpha())) { + return false; + } rEventSelection.fill(HIST("htrackscheck_v0"), 7.5); if (std::fabs(ctauK0s) > config.cMaxV0LifeTime) { @@ -799,22 +805,39 @@ struct HigherMassResonances { enum FrameType { kHelicity = 0, kCollinsSoper = 1, kProduction = 2, - kBeamAxis = 3, + kGottfriedJackson = 3, kRandom = 4 }; std::pair getCosThetaPhi(const ROOT::Math::PxPyPzMVector& motherVec, const ROOT::Math::PxPyPzMVector& daughterVec, int frame) { + //---------------------------------------------------------------- + // Boost daughter into mother rest 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()); + // Unit vector of daughter momentum in mother rest frame + auto v1cm = ROOT::Math::XYZVectorF(dauCM.Vect()).Unit(); + + //---------------------------------------------------------------- + // Beam directions boosted into mother rest frame + //---------------------------------------------------------------- + auto beam1cm = ROOT::Math::XYZVectorF((boost(beam1).Vect()).Unit()); + + auto beam2cm = ROOT::Math::XYZVectorF((boost(beam2).Vect()).Unit()); double cosTheta = 0.; double phi = 0.; if (frame == kHelicity) { + //============================================================== + // HELICITY FRAME + // z-axis : Mother momentum direction in laboratory frame + // y-axis : Normal to beam plane + // x-axis : Completes right-handed coordinate system + //============================================================== + 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(); @@ -822,6 +845,13 @@ struct HigherMassResonances { phi = std::atan2(yaxisHELocal.Dot(v1cm), xaxisHELocal.Dot(v1cm)); phi = RecoDecay::constrainAngle(phi, 0.0); } else if (frame == kCollinsSoper) { + //============================================================== + // COLLINS-SOPER FRAME + // z-axis : Bisector of the two beam directions + // y-axis : Normal to beam plane + // x-axis : Completes right-handed system + //============================================================== + 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(); @@ -829,6 +859,16 @@ struct HigherMassResonances { phi = std::atan2(yAxisCSLocal.Dot(v1cm), xAxisCSLocal.Dot(v1cm)); phi = RecoDecay::constrainAngle(phi, 0.0); } else if (frame == kProduction) { + //============================================================== + // PRODUCTION FRAME + // z-axis : Normal to production plane + // + // Production plane is defined by: + // beam direction + // mother momentum + // Used occasionally in spin-alignment analyses. + //============================================================== + 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 @@ -837,13 +877,37 @@ struct HigherMassResonances { 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)); + } else if (frame == kGottfriedJackson) { + // ------------------------------------------------------------------ + // Gottfried-Jackson (GJ') frame + // + // z-axis : beam direction in the mother rest frame + // y-axis : normal to the production plane + // x-axis : right-handed coordinate system + // ------------------------------------------------------------------ + + // z-axis: beam direction boosted into mother rest frame + auto zAxisGJ = beam1cm.Unit(); + + // Mother momentum in LAB (defines production plane) + auto motherLab = ROOT::Math::XYZVectorF(motherVec.Vect()).Unit(); + + // Normal to production plane + auto yAxisGJ = ROOT::Math::XYZVectorF(zAxisGJ.Cross(motherLab)); + + if (yAxisGJ.Mag2() > config.tolerance) + yAxisGJ = yAxisGJ.Unit(); + else + yAxisGJ = ROOT::Math::XYZVectorF(0.f, 1.f, 0.f); + + // Complete right-handed system + ROOT::Math::XYZVectorF xAxisGJ = ROOT::Math::XYZVectorF(yAxisGJ.Cross(zAxisGJ)).Unit(); + + // Daughter direction in mother rest frame + cosTheta = zAxisGJ.Dot(v1cm); + + phi = std::atan2(yAxisGJ.Dot(v1cm), xAxisGJ.Dot(v1cm)); + phi = RecoDecay::constrainAngle(phi, 0.0); } else { // kRandom or fallback auto phiRandom = gRandom->Uniform(0.f, constants::math::TwoPI); @@ -1008,8 +1072,8 @@ struct HigherMassResonances { hglue.fill(HIST("h3glueInvMassME"), eventMultiplicity, motherVec.Pt(), motherVec.M(), pr.first, pr.second); } } - } else if (config.activateBeamAxisFrame) { - auto pr = getCosThetaPhi(motherVec, daughterVec1, kBeamAxis); + } else if (config.activateGJFrame) { + auto pr = getCosThetaPhi(motherVec, daughterVec1, kGottfriedJackson); if (!isMixed) { if (std::abs(motherVec.Rapidity()) < config.rapidityMotherData) { hglue.fill(HIST("h3glueInvMassDS"), eventMultiplicity, motherVec.Pt(), motherVec.M(), pr.first, pr.second); @@ -1017,7 +1081,7 @@ struct HigherMassResonances { for (int i = 0; i < config.cRotations; i++) { 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); + auto prrot = getCosThetaPhi(motherRot, daughterVec1, kGottfriedJackson); if (std::abs(motherRot.Rapidity()) < config.rapidityMotherData) { hglue.fill(HIST("h3glueInvMassRot"), eventMultiplicity, motherRot.Pt(), motherRot.M(), prrot.first, prrot.second); } @@ -1612,8 +1676,8 @@ struct HigherMassResonances { frame = kCollinsSoper; } else if (config.activateProductionFrame) { frame = kProduction; - } else if (config.activateBeamAxisFrame) { - frame = kBeamAxis; + } else if (config.activateGJFrame) { + frame = kGottfriedJackson; } else if (config.activateRandomFrame) { frame = kRandom; } else { @@ -1630,6 +1694,7 @@ struct HigherMassResonances { hMChists.fill(HIST("GenEta"), mcParticle.eta()); hMChists.fill(HIST("GenPhi"), mcParticle.phi()); hMChists.fill(HIST("Gen1710PWA"), lResonanceGen.M(), angularVar.first, angularVar.second); + hMChists.fill(HIST("GenThetavsPhi1"), angularVar.first, angularVar.second); if (config.isapplyPairRapidityMC && std::abs(lResonanceGen1.Rapidity()) >= config.rapidityMotherData) { continue; @@ -1642,6 +1707,7 @@ struct HigherMassResonances { hMChists.fill(HIST("GenEta2"), lResonanceGen1.Eta()); hMChists.fill(HIST("GenPhi2"), lResonanceGen1.Phi()); hMChists.fill(HIST("Gen1710PWA2"), lResonanceGen1.M(), angularVar1.first, angularVar1.second); + hMChists.fill(HIST("GenThetavsPhi2"), angularVar1.first, angularVar1.second); } passKs.clear(); // clear the vector for the next iteration } @@ -1903,8 +1969,8 @@ struct HigherMassResonances { frame = kCollinsSoper; } else if (config.activateProductionFrame) { frame = kProduction; - } else if (config.activateBeamAxisFrame) { - frame = kBeamAxis; + } else if (config.activateGJFrame) { + frame = kGottfriedJackson; } else if (config.activateRandomFrame) { frame = kRandom; } @@ -1918,6 +1984,7 @@ struct HigherMassResonances { hMChists.fill(HIST("RecPhi"), mothertrack1.phi()); hMChists.fill(HIST("RecEta"), mothertrack1.eta()); hMChists.fill(HIST("Recf1710PWA_pt1"), mother1.M(), angularVar1.first, angularVar1.second); + hMChists.fill(HIST("RecThetavsPhi1"), angularVar1.first, angularVar1.second); if (config.isapplyPairRapidityMC && std::abs(mother.Rapidity()) >= config.rapidityMotherData) { continue; @@ -1929,6 +1996,7 @@ struct HigherMassResonances { hMChists.fill(HIST("RecPhi2"), mother.Phi()); hMChists.fill(HIST("RecEta2"), mother.Eta()); hMChists.fill(HIST("Recf1710PWA_pt2"), mother.M(), angularVar.first, angularVar.second); + hMChists.fill(HIST("RecThetavsPhi2"), angularVar.first, angularVar.second); } gindex2.clear(); } diff --git a/PWGLF/Tasks/Resonances/kstarqa.cxx b/PWGLF/Tasks/Resonances/kstarqa.cxx index 8dd87cb002a..91cc416aca4 100644 --- a/PWGLF/Tasks/Resonances/kstarqa.cxx +++ b/PWGLF/Tasks/Resonances/kstarqa.cxx @@ -179,6 +179,8 @@ struct Kstarqa { } configGp; + Service pdgDB{}; + enum MultEstimator { kFT0M, kFT0A, @@ -361,6 +363,12 @@ struct Kstarqa { hInvMass.add("h1RecMult", "Multiplicity reconstructed", kTH1F, {multiplicityAxis}); hInvMass.add("h1RecMult2", "Multiplicity reconstructed", kTH1F, {multiplicityAxis}); hInvMass.add("h1KSRecsplit", "KS meson Rec split", kTH1F, {{100, 0.0f, 10.0f}}); + + hInvMass.add("h3KstarMassGen", "Gen pt, multiplicity, mass", kTHnSparseF, {ptAxis, multiplicityAxis, invmassAxis}); + hInvMass.add("h3KstarMassGenCalib", "Gen pt, multiplicity, mass", kTHnSparseF, {ptAxis, multiplicityAxis, invmassAxis}); + hInvMass.add("h3KstarMassRec", "Rec pt, multiplicity, mass", kTHnSparseF, {ptAxis, multiplicityAxis, invmassAxis}); + hInvMass.add("h3KstarMassRecCalib", "Rec pt, multiplicity, mass", kTHnSparseF, {ptAxis, multiplicityAxis, invmassAxis}); + // hInvMass.add("hAllGenCollisionsImpact", "All generated collisions vs impact parameter", kTH1F, {multiplicityAxis}); hInvMass.add("hAllGenCollisions", "All generated events", kTH1F, {multiplicityAxis}); hInvMass.add("hAllGenCollisions1Rec", "All gen events with at least one rec event", kTH1F, {multiplicityAxis}); @@ -561,52 +569,73 @@ struct Kstarqa { if (configGp.isGlobalTracks) { if (!candidate.isGlobalTrack()) return false; + if (std::abs(candidate.pt()) < configGp.cfgCutPT) return false; + if (std::abs(candidate.eta()) > configGp.cfgCutEtaMax) return false; + if (!configGp.isApplyPtDepDCAxyCut) { if (std::abs(candidate.dcaXY()) > configGp.cfgCutDCAxyMax) return false; + } else { if (std::abs(candidate.dcaXY()) > (0.0105 + 0.035 / std::pow(candidate.pt(), 1.1))) return false; } if (std::abs(candidate.dcaZ()) > configGp.cfgCutDCAz) return false; + if (candidate.itsNCls() < configGp.cfgITScluster) return false; + if (candidate.tpcNClsFound() < configGp.cfgTPCcluster) return false; + if (configGp.hasITS && !candidate.hasITS()) return false; - if (configGp.isITSTPCRefit && (!(o2::aod::track::ITSrefit) || !(o2::aod::track::TPCrefit))) - return false; + + if (configGp.isITSTPCRefit) { + if (!(candidate.flags() & o2::aod::track::ITSrefit) || + !(candidate.flags() & o2::aod::track::TPCrefit)) { + return false; + } + } if (configGp.cfgPVContributor && !candidate.isPVContributor()) return false; + } else if (!configGp.isGlobalTracks) { if (std::abs(candidate.pt()) < configGp.cfgCutPT) return false; - // if (std::abs(candidate.eta()) > configGp.cfgCutEtaMax || std::abs(candidate.eta()) < configGp.cfgCutEtaMin) + if (std::abs(candidate.eta()) > configGp.cfgCutEtaMax) return false; - // if (std::abs(candidate.dcaXY()) > configGp.cfgCutDCAxyMax || std::abs(candidate.dcaXY()) < configGp.cfgCutDCAxyMin) + if (std::abs(candidate.dcaXY()) > configGp.cfgCutDCAxyMax) return false; + if (std::abs(candidate.dcaZ()) > configGp.cfgCutDCAz) return false; + // if (candidate.tpcCrossedRowsOverFindableCls() < configGp.cfgRCRFC) - return false; + // return false; + if (candidate.itsNCls() < configGp.cfgITScluster) return false; + if (candidate.tpcNClsFound() < configGp.cfgTPCcluster) return false; + // if (candidate.itsChi2NCl() >= configGp.cfgITSChi2NCl) // return false; + // if (candidate.tpcChi2NCl() >= configGp.cfgTPCChi2NClMax || candidate.tpcChi2NCl() < configGp.cfgTPCChi2NClMin) // return false; + if (configGp.cfgPVContributor && !candidate.isPVContributor()) return false; + if (!candidate.isPrimaryTrack()) return false; } @@ -619,32 +648,22 @@ struct Kstarqa { { const auto pglobal = track.p(); const auto ptpc = track.tpcInnerParam(); - if (std::abs(pglobal - ptpc) > configGp.cFakeTrackCutKa) { - return true; - } - return false; + return std::abs(pglobal - ptpc) > configGp.cFakeTrackCutKa; } // deep angle cut on pair to remove photon conversion template bool selectionPair(const T1& candidate1, const T2& candidate2) { - double pt1, pt2, pz1, pz2, p1, p2, angle; - pt1 = candidate1.pt(); - pt2 = candidate2.pt(); - pz1 = candidate1.pz(); - pz2 = candidate2.pz(); - p1 = candidate1.p(); - p2 = candidate2.p(); - angle = std::acos((pt1 * pt2 + pz1 * pz2) / (p1 * p2)); - if (configGp.isApplyDeepAngle && angle < configGp.cfgDeepAngle) { - return false; - } - // double deltaRvalue = std::sqrt(TVector2::Phi_mpi_pi(candidate1.phi() - candidate2.phi()) * TVector2::Phi_mpi_pi(candidate1.phi() - candidate2.phi()) + (candidate1.eta() - candidate2.eta()) * (candidate1.eta() - candidate2.eta())); - // if (deltaRvalue < configGp.deltaRCut) { - // return false; - // } - return true; + const double pt1 = candidate1.pt(); + const double pt2 = candidate2.pt(); + const double pz1 = candidate1.pz(); + const double pz2 = candidate2.pz(); + const double p1 = candidate1.p(); + const double p2 = candidate2.p(); + const double angle = std::acos((pt1 * pt2 + pz1 * pz2) / (p1 * p2)); + + return !configGp.isApplyDeepAngle || angle >= configGp.cfgDeepAngle; } template @@ -1021,8 +1040,8 @@ struct Kstarqa { return false; } - std::array pvec0; - std::array pvec1; + std::array pvec0 = {0.0, 0.0, 0.0}; + std::array pvec1 = {0.0, 0.0, 0.0}; // Defining filters for events (event selection) // Processed events will be already fulfilling the event selection @@ -1049,53 +1068,75 @@ struct Kstarqa { using LabeledTracks = soa::Join; //*********Varibles declaration*************** - float multiplicity{-1.0}, theta2; - ROOT::Math::PxPyPzMVector daughter1, daughter2, daughterRot, mother, motherRot, daughterSelected, fourVecDauCM, daughterRotCM; + float multiplicity{-1.0}, theta2{0.0}; + ROOT::Math::PxPyPzMVector daughter1, daughter2, daughterRot, mother, motherRot, daughterSelected, fourVecDauCM, daughterRotCM, genDaughter1, genDaughter2, genMother; ROOT::Math::XYZVector randomVec, beamVec, normalVec; bool isMix = false; template - void fillInvMass(const T1& daughter1, const T1& daughter2, const T1& mother, float multiplicity, bool isMix, const T2& track1, const T2& track2) + void fillInvMass(const T1& dau1, const T1& dau2, const T1& motherVec, float multiplicityFill, bool isMixInput, const T2& track1, const T2& track2) { - daughterSelected = (boostDaugter1) ? daughter1 : daughter2; // polarization calculations - ROOT::Math::Boost boost{mother.BoostToCM()}; // boost mother to center of mass frame - fourVecDauCM = boost(daughterSelected); // boost the frame of daughter same as mother + daughterSelected = (boostDaugter1) ? dau1 : dau2; // polarization calculations + + // Boost mother to the center-of-mass frame + ROOT::Math::Boost boost{motherVec.BoostToCM()}; + + // Boost selected daughter to the mother rest frame + fourVecDauCM = boost(daughterSelected); - // if (std::abs(mother.Rapidity()) < configGp.rapidityMotherData) { if (activateTHnSparseCosThStarHelicity) { - auto cosThetaStarHelicity = mother.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(mother.Vect().Mag2())); + auto cosThetaStarHelicity = motherVec.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(motherVec.Vect().Mag2())); + + // Unlike-sign pairs if (track1.sign() * track2.sign() < 0) { - if (!isMix) { - if (std::abs(mother.Rapidity()) < configGp.rapidityMotherData) { - hInvMass.fill(HIST("h3KstarInvMassUnlikeSign"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity); + + // Same-event + if (!isMixInput) { + + if (std::abs(motherVec.Rapidity()) < configGp.rapidityMotherData) { + hInvMass.fill(HIST("h3KstarInvMassUnlikeSign"), multiplicityFill, motherVec.Pt(), motherVec.M(), cosThetaStarHelicity); } + // Rotational background for (int i = 0; i < cRotations; i++) { + theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / configGp.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / configGp.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(dau1.Px() * std::cos(theta2) - dau1.Py() * std::sin(theta2), dau1.Px() * std::sin(theta2) + dau1.Py() * std::cos(theta2), dau1.Pz(), dau1.M()); - motherRot = daughterRot + daughter2; + motherRot = daughterRot + dau2; 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())); - if (calcRotational && std::abs(motherRot.Rapidity()) < configGp.rapidityMotherData) - hInvMass.fill(HIST("h3KstarInvMassRotated"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarHelicityRot); + if (calcRotational && std::abs(motherRot.Rapidity()) < configGp.rapidityMotherData) { + hInvMass.fill(HIST("h3KstarInvMassRotated"), multiplicityFill, motherRot.Pt(), motherRot.M(), cosThetaStarHelicityRot); + } } - } else if (isMix && std::abs(mother.Rapidity()) < configGp.rapidityMotherData) { - hInvMass.fill(HIST("h3KstarInvMassMixed"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity); + + // Mixed-event + } else if (std::abs(motherVec.Rapidity()) < configGp.rapidityMotherData) { + + hInvMass.fill(HIST("h3KstarInvMassMixed"), multiplicityFill, motherVec.Pt(), motherVec.M(), cosThetaStarHelicity); } + + // Like-sign pairs } else { - if (!isMix) { - if (calcLikeSign && std::abs(mother.Rapidity()) < configGp.rapidityMotherData) { + + if (!isMixInput) { + + if (calcLikeSign && std::abs(motherVec.Rapidity()) < configGp.rapidityMotherData) { + if (track1.sign() > 0 && track2.sign() > 0) { - hInvMass.fill(HIST("h3KstarInvMasslikeSignPP"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity); + + hInvMass.fill(HIST("h3KstarInvMasslikeSignPP"), multiplicityFill, motherVec.Pt(), motherVec.M(), cosThetaStarHelicity); + } else if (track1.sign() < 0 && track2.sign() < 0) { - hInvMass.fill(HIST("h3KstarInvMasslikeSignMM"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity); + + hInvMass.fill(HIST("h3KstarInvMasslikeSignMM"), multiplicityFill, motherVec.Pt(), motherVec.M(), cosThetaStarHelicity); } } } @@ -1124,8 +1165,6 @@ struct Kstarqa { multiplicity = collision.centFT0C(); } else if (cSelectMultEstimator == kFV0A) { multiplicity = collision.centFV0A(); - } else { - multiplicity = collision.centFT0M(); // default } // Fill the event counter @@ -1317,7 +1356,7 @@ struct Kstarqa { void processME(EventCandidatesMix const&, TrackCandidates const&) { // Map estimator to pair and multiplicity accessor - auto runMixing = [&](auto& pair, auto multiplicityGetter) { + auto runMixing = [&](const auto& pair, auto multiplicityGetter) { for (const auto& [c1, tracks1, c2, tracks2] : pair) { // if (!c1.sel8() || !c2.sel8()) // continue; @@ -1419,7 +1458,7 @@ struct Kstarqa { void processMEMC(EventCandidatesMC const&, TrackCandidatesMC const&, aod::McParticles const&, aod::McCollisions const&) { - auto runMixing = [&](auto& pair, auto multiplicityGetter) { + auto runMixing = [&](const auto& pair, auto multiplicityGetter) { for (const auto& [c1, tracks1, c2, tracks2] : pair) { if (!selectionEvent(c1, false) || !selectionEvent(c2, false)) { // don't fill event cut histogram @@ -1539,8 +1578,6 @@ struct Kstarqa { multiplicity = collision.centFT0C(); } else if (cSelectMultEstimator == kFV0A) { multiplicity = collision.centFV0A(); - } else { - multiplicity = collision.centFT0M(); // default } // Fill the event counter @@ -1713,8 +1750,6 @@ struct Kstarqa { } PROCESS_SWITCH(Kstarqa, processSEMC, "Process same event in MC", false); - Service pdgDB; - void processGen(EventMCGenerated::iterator const& mcCollision, aod::McParticles const& mcParticles, const soa::SmallGroups& collisions) // void processGen(aod::McCollision const& mcCollision, aod::McParticles const& mcParticles, const soa::SmallGroups& collisions) { @@ -1760,8 +1795,6 @@ struct Kstarqa { multiplicity = collision.centFT0C(); } else if (cSelectMultEstimator == kFV0A) { multiplicity = collision.centFV0A(); - } else { - multiplicity = collision.centFT0M(); // default } hInvMass.fill(HIST("h1GenMult"), multiplicity); @@ -1879,8 +1912,6 @@ struct Kstarqa { multiplicity1 = RecCollision.centFT0C(); } else if (cSelectMultEstimator == kFV0A) { multiplicity1 = RecCollision.centFV0A(); - } else { - multiplicity1 = RecCollision.centFT0M(); // default } isSelectedEvent = true; @@ -1896,7 +1927,7 @@ struct Kstarqa { } auto impactPar = mcCollision.impactParameter(); - auto multiplicityGen = -1; + auto multiplicityGen = -1.0f; multiplicityGen = mcCollision.centFT0M(); hInvMass.fill(HIST("MCcorrections/hImpactParameterGen"), impactPar); hInvMass.fill(HIST("MCcorrections/MultiplicityGen"), multiplicityGen); @@ -1941,7 +1972,7 @@ struct Kstarqa { } hInvMass.fill(HIST("CorrFactors/hGenEvents"), multiplicityNch, 2.5); - float multiplicity = -1.0; + float multiplicitySigLoss = -1.0; bool isSelectedEvent = false; for (auto const& collision : collisions) { @@ -1951,15 +1982,13 @@ struct Kstarqa { continue; if (cSelectMultEstimator == kFT0M) { - multiplicity = collision.centFT0M(); + multiplicitySigLoss = collision.centFT0M(); } else if (cSelectMultEstimator == kFT0A) { - multiplicity = collision.centFT0A(); + multiplicitySigLoss = collision.centFT0A(); } else if (cSelectMultEstimator == kFT0C) { - multiplicity = collision.centFT0C(); + multiplicitySigLoss = collision.centFT0C(); } else if (cSelectMultEstimator == kFV0A) { - multiplicity = collision.centFV0A(); - } else { - multiplicity = collision.centFT0M(); // default + multiplicitySigLoss = collision.centFV0A(); } isSelectedEvent = true; } @@ -1967,7 +1996,7 @@ struct Kstarqa { // auto multiplicityGen = -1; // multiplicityGen = mcCollision.centFT0M(); - hInvMass.fill(HIST("CorrFactors/hMultiplicityVsMultMC"), multiplicity, multiplicityNch); + hInvMass.fill(HIST("CorrFactors/hMultiplicityVsMultMC"), multiplicitySigLoss, multiplicityNch); hInvMass.fill(HIST("CorrFactors/hNrecInGen"), collisions.size()); hInvMass.fill(HIST("CorrFactors/MultiplicityGen"), multiplicityNch); if (isSelectedEvent) { @@ -1996,7 +2025,12 @@ struct Kstarqa { // continue; int pdgDau = kCurrentDaughter.pdgCode(); - int sign = (pdgDau > 0) - (pdgDau < 0); + int sign = 0; + if (pdgDau > 0) { + sign = 1; + } else if (pdgDau < 0) { + sign = -1; + } if (sign > 0) hasPos = true; @@ -2016,8 +2050,8 @@ struct Kstarqa { if ((passkaon && passpion) && (hasPos && hasNeg)) { mother = daughter1 + daughter2; // Kstar meson - hInvMass.fill(HIST("CorrFactors/h2dGenKstar"), multiplicity, mother.Pt()); - hInvMass.fill(HIST("CorrFactors/h3dGenKstarVsMultMCVsMultiplicity"), multiplicityNch, multiplicity, mother.Pt()); + hInvMass.fill(HIST("CorrFactors/h2dGenKstar"), multiplicitySigLoss, mother.Pt()); + hInvMass.fill(HIST("CorrFactors/h3dGenKstarVsMultMCVsMultiplicity"), multiplicityNch, multiplicitySigLoss, mother.Pt()); hInvMass.fill(HIST("CorrFactors/hSignalLoss1"), mother.pt(), multiplicityNch); if (isSelectedEvent) { hInvMass.fill(HIST("CorrFactors/hSignalLoss2"), mother.pt(), multiplicityNch); @@ -2033,6 +2067,8 @@ struct Kstarqa { } PROCESS_SWITCH(Kstarqa, processEvSigLossFactors, "Process Event and Signal loss", false); + double genMass = 0.0, recMass = 0.0, recPt = 0.0, genPt = 0.0; + void processRec(EventCandidatesMC::iterator const& collision, TrackCandidatesMC const& tracks, aod::McParticles const&, EventMCGenerated const&) { rEventSelection.fill(HIST("eventsCheckRec"), 0.5); @@ -2061,8 +2097,6 @@ struct Kstarqa { multiplicity = collision.centFT0C(); } else if (cSelectMultEstimator == kFV0A) { multiplicity = collision.centFV0A(); - } else { - multiplicity = collision.centFT0M(); // default } hInvMass.fill(HIST("hAllRecCollisions"), multiplicity); @@ -2144,7 +2178,8 @@ struct Kstarqa { } rEventSelection.fill(HIST("recMCparticles"), 6.5); - if (!(track1PDG == PDG_t::kPiPlus && track2PDG == PDG_t::kKPlus) && !(track1PDG == PDG_t::kKPlus && track2PDG == PDG_t::kPiPlus)) { + if ((track1PDG != PDG_t::kPiPlus || track2PDG != PDG_t::kKPlus) && + (track1PDG != PDG_t::kKPlus || track2PDG != PDG_t::kPiPlus)) { continue; } rEventSelection.fill(HIST("recMCparticles"), 7.5); @@ -2300,13 +2335,35 @@ struct Kstarqa { oldindex = mothertrack1.globalIndex(); if (track1PDG == PDG_t::kPiPlus) { + daughter1 = ROOT::Math::PxPyPzMVector(track1.px(), track1.py(), track1.pz(), massPi); daughter2 = ROOT::Math::PxPyPzMVector(track2.px(), track2.py(), track2.pz(), massKa); + + genDaughter1 = ROOT::Math::PxPyPzMVector(mctrack1.px(), mctrack1.py(), mctrack1.pz(), massPi); + genDaughter2 = ROOT::Math::PxPyPzMVector(mctrack2.px(), mctrack2.py(), mctrack2.pz(), massKa); + } else if (track1PDG == PDG_t::kKPlus) { + daughter1 = ROOT::Math::PxPyPzMVector(track1.px(), track1.py(), track1.pz(), massKa); daughter2 = ROOT::Math::PxPyPzMVector(track2.px(), track2.py(), track2.pz(), massPi); + + genDaughter1 = ROOT::Math::PxPyPzMVector(mctrack1.px(), mctrack1.py(), mctrack1.pz(), massKa); + genDaughter2 = ROOT::Math::PxPyPzMVector(mctrack2.px(), mctrack2.py(), mctrack2.pz(), massPi); } - mother = daughter1 + daughter2; // Kstar meson + mother = daughter1 + daughter2; // Kstar meson + genMother = genDaughter1 + genDaughter2; // Gen Kstar from MC daughters + + recMass = mother.M(); + recPt = mother.Pt(); + + genMass = genMother.M(); + genPt = mothertrack1.pt(); + + hInvMass.fill(HIST("h3KstarMassGen"), genPt, multiplicity, genMass); + hInvMass.fill(HIST("h3KstarMassGenCalib"), genPt, multiplicityRec, genMass); + + hInvMass.fill(HIST("h3KstarMassRec"), recPt, multiplicity, recMass); + hInvMass.fill(HIST("h3KstarMassRecCalib"), recPt, multiplicityRec, recMass); hInvMass.fill(HIST("h2KstarRecpt2"), mothertrack1.pt(), multiplicity, std::sqrt(mothertrack1.e() * mothertrack1.e() - mothertrack1.p() * mothertrack1.p())); hInvMass.fill(HIST("h2KstarRecptCalib2"), mothertrack1.pt(), multiplicityRec, std::sqrt(mothertrack1.e() * mothertrack1.e() - mothertrack1.p() * mothertrack1.p())); @@ -2318,6 +2375,7 @@ struct Kstarqa { hInvMass.fill(HIST("h1KstarRecMass"), mother.M()); hInvMass.fill(HIST("h2KstarRecpt1"), mother.Pt(), multiplicity, mother.M()); hInvMass.fill(HIST("h2KstarRecptCalib1"), mother.Pt(), multiplicityRec, mother.M()); + // h3KstarMassGen } } } @@ -2349,8 +2407,6 @@ struct Kstarqa { multiplicity = collision.centFT0C(); } else if (cSelectMultEstimator == kFV0A) { multiplicity = collision.centFV0A(); - } else { - multiplicity = collision.centFT0M(); // default } hInvMass.fill(HIST("hAllRecCollisions"), multiplicity); @@ -2561,8 +2617,6 @@ struct Kstarqa { multiplicity = collision.centFT0C(); } else if (cSelectMultEstimator == kFV0A) { multiplicity = collision.centFV0A(); - } else { - multiplicity = collision.centFT0M(); // default } // Fill the event counter @@ -2658,7 +2712,7 @@ struct Kstarqa { void processMEPhi(EventCandidatesMix const&, TrackCandidates const&) { // Map estimator to pair and multiplicity accessor - auto runMixing = [&](auto& pair, auto multiplicityGetter) { + auto runMixing = [&](const auto& pair, auto multiplicityGetter) { for (const auto& [c1, tracks1, c2, tracks2] : pair) { // if (!c1.sel8() || !c2.sel8()) // continue; @@ -2758,8 +2812,6 @@ struct Kstarqa { multiplicity = collision.centFT0C(); } else if (cSelectMultEstimator == kFV0A) { multiplicity = collision.centFV0A(); - } else { - multiplicity = collision.centFT0M(); // default } hInvMass.fill(HIST("h1GenMult"), multiplicity); @@ -2846,7 +2898,7 @@ struct Kstarqa { } hInvMass.fill(HIST("CorrFactors/hGenEvents"), multiplicityNch, 2.5); - float multiplicity = -1.0; + float multiplicityPhi = -1.0; bool isSelectedEvent = false; for (auto const& collision : collisions) { @@ -2856,15 +2908,13 @@ struct Kstarqa { continue; if (cSelectMultEstimator == kFT0M) { - multiplicity = collision.centFT0M(); + multiplicityPhi = collision.centFT0M(); } else if (cSelectMultEstimator == kFT0A) { - multiplicity = collision.centFT0A(); + multiplicityPhi = collision.centFT0A(); } else if (cSelectMultEstimator == kFT0C) { - multiplicity = collision.centFT0C(); + multiplicityPhi = collision.centFT0C(); } else if (cSelectMultEstimator == kFV0A) { - multiplicity = collision.centFV0A(); - } else { - multiplicity = collision.centFT0M(); // default + multiplicityPhi = collision.centFV0A(); } isSelectedEvent = true; } @@ -2872,7 +2922,7 @@ struct Kstarqa { // auto multiplicityGen = -1; // multiplicityGen = mcCollision.centFT0M(); - hInvMass.fill(HIST("CorrFactors/hMultiplicityVsMultMC"), multiplicity, multiplicityNch); + hInvMass.fill(HIST("CorrFactors/hMultiplicityVsMultMC"), multiplicityPhi, multiplicityNch); hInvMass.fill(HIST("CorrFactors/hNrecInGen"), collisions.size()); hInvMass.fill(HIST("CorrFactors/MultiplicityGen"), multiplicityNch); if (isSelectedEvent) { @@ -2901,7 +2951,12 @@ struct Kstarqa { // continue; int pdgDau = kCurrentDaughter.pdgCode(); - int sign = (pdgDau > 0) - (pdgDau < 0); + int sign = 0; + if (pdgDau > 0) { + sign = 1; + } else if (pdgDau < 0) { + sign = -1; + } if (sign > 0) hasPos = true; @@ -2921,8 +2976,8 @@ struct Kstarqa { if ((passkaon && passpion) && (hasPos && hasNeg)) { mother = daughter1 + daughter2; // Kstar meson - hInvMass.fill(HIST("CorrFactors/h2dGenKstar"), multiplicity, mother.Pt()); - hInvMass.fill(HIST("CorrFactors/h3dGenKstarVsMultMCVsMultiplicity"), multiplicityNch, multiplicity, mother.Pt()); + hInvMass.fill(HIST("CorrFactors/h2dGenKstar"), multiplicityPhi, mother.Pt()); + hInvMass.fill(HIST("CorrFactors/h3dGenKstarVsMultMCVsMultiplicity"), multiplicityNch, multiplicityPhi, mother.Pt()); hInvMass.fill(HIST("CorrFactors/hSignalLoss1"), mother.pt(), multiplicityNch); if (isSelectedEvent) { hInvMass.fill(HIST("CorrFactors/hSignalLoss2"), mother.pt(), multiplicityNch); @@ -2963,8 +3018,6 @@ struct Kstarqa { multiplicity = collision.centFT0C(); } else if (cSelectMultEstimator == kFV0A) { multiplicity = collision.centFV0A(); - } else { - multiplicity = collision.centFT0M(); // default } hInvMass.fill(HIST("hAllRecCollisions"), multiplicity); @@ -3044,7 +3097,7 @@ struct Kstarqa { } rEventSelection.fill(HIST("recMCparticles"), 5.5); - if ((track1PDG != PDG_t::kKPlus) || (track2PDG != PDG_t::kKPlus)) { + if (track1PDG != PDG_t::kKPlus || track2PDG != PDG_t::kKPlus) { continue; } rEventSelection.fill(HIST("recMCparticles"), 6.5); @@ -3076,9 +3129,10 @@ struct Kstarqa { } rEventSelection.fill(HIST("recMCparticles"), 11.5); - if (!configGp.isapplypTdepPID && !(selectionPID(track1, 1) && selectionPID(track2, 1))) { // kaon and kaon + if (!configGp.isapplypTdepPID && !(selectionPID(track1, 1) && selectionPID(track2, 1))) { continue; - } else if (configGp.isapplypTdepPID && !(selectionPIDPtDep(track1, 1) && selectionPIDPtDep(track2, 1))) { // kaon and kaon + } + if (configGp.isapplypTdepPID && !(selectionPIDPtDep(track1, 1) && selectionPIDPtDep(track2, 1))) { continue; } rEventSelection.fill(HIST("recMCparticles"), 12.5); @@ -3147,8 +3201,6 @@ struct Kstarqa { centrality = RecCollision.centFT0C(); } else if (cSelectMultEstimator == kFV0A) { centrality = RecCollision.centFV0A(); - } else { - centrality = RecCollision.centFT0M(); // default } isSelectedEvent = true; @@ -3197,7 +3249,7 @@ struct Kstarqa { auto impactPar = mcCollision.impactParameter(); // auto multiplicityRec = -1; - auto multiplicityGen = -1; + auto multiplicityGen = -1.0f; multiplicityGen = mcCollision.centFT0M(); hInvMass.fill(HIST("MCcorrections/hImpactParameterGen"), impactPar); hInvMass.fill(HIST("MCcorrections/MultiplicityGen"), multiplicityGen); @@ -3221,8 +3273,6 @@ struct Kstarqa { multiplicity1 = RecCollision.centFT0C(); } else if (cSelectMultEstimator == kFV0A) { multiplicity1 = RecCollision.centFV0A(); - } else { - multiplicity1 = RecCollision.centFT0M(); // default } isSelectedEvent = true; }