diff --git a/PWGLF/Tasks/Resonances/chargedkstaranalysis.cxx b/PWGLF/Tasks/Resonances/chargedkstaranalysis.cxx index 71eaafc4209..6081e136eec 100644 --- a/PWGLF/Tasks/Resonances/chargedkstaranalysis.cxx +++ b/PWGLF/Tasks/Resonances/chargedkstaranalysis.cxx @@ -17,6 +17,7 @@ /// \author Navneet Kumar , Protay #include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGLF/DataModel/mcCentrality.h" #include "PWGLF/Utils/collisionCuts.h" #include "PWGLF/Utils/inelGt.h" @@ -69,6 +70,8 @@ #include #include #include +#include +#include #include using namespace o2; @@ -84,6 +87,9 @@ struct Chargedkstaranalysis { FT0C = 1, FT0M = 2 }; + SliceCache cache; + Preslice perCollision = aod::track::collisionId; + struct : ConfigurableGroup { ConfigurableAxis cfgvtxbins{"cfgvtxbins", {VARIABLE_WIDTH, -10.0f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - z-vertex"}; ConfigurableAxis cfgmultbins{"cfgmultbins", {VARIABLE_WIDTH, 0., 1., 5., 10., 30., 50., 70., 100., 110.}, "Mixing bins - multiplicity"}; @@ -106,7 +112,7 @@ struct Chargedkstaranalysis { Configurable activateProductionFrame{"activateProductionFrame", false, "Activate the THnSparse with cosThStar w.r.t. production axis"}; Configurable activateBeamAxisFrame{"activateBeamAxisFrame", false, "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"}; + Configurable cRotations{"cRotations", 5, "Number of random rotations in the rotational background"}; Configurable cBoostKShot{"cBoostKShot", true, "Boost the Kshot in Charged Kstar frame of reference"}; // Other cuts on Ks Configurable rotationalCut{"rotationalCut", 10, "Cut value (Rotation angle pi - pi/cut and pi + pi/cut)"}; @@ -125,9 +131,15 @@ struct Chargedkstaranalysis { using TrackCandidates = soa::Join; using V0Candidates = aod::V0Datas; + // for MC reco using MCEventCandidates = soa::Join; using MCTrackCandidates = soa::Join; using MCV0Candidates = soa::Join; + // for MC truth + using MCTrueEventCandidates = aod::McCollisions; + using MCTrueTrackCandidates = aod::McParticles; + + using LorentzVectorSetXYZM = ROOT::Math::LorentzVector>; HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; HistogramRegistry hChaKstar{"hChaKstar", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; @@ -158,7 +170,7 @@ struct Chargedkstaranalysis { Configurable cNbinsDiv{"cNbinsDiv", 1, "Integer to divide the number of bins"}; struct RCTCut : ConfigurableGroup { - Configurable requireRCTFlagChecker{"requireRCTFlagChecker", true, "Check event quality in run condition table"}; + Configurable requireRCTFlagChecker{"requireRCTFlagChecker", false, "Check event quality in run condition table"}; Configurable cfgEvtRCTFlagCheckerLabel{"cfgEvtRCTFlagCheckerLabel", "CBT_hadronPID", "Evt sel: RCT flag checker label"}; Configurable cfgEvtRCTFlagCheckerZDCCheck{"cfgEvtRCTFlagCheckerZDCCheck", false, "Evt sel: RCT flag checker ZDC check"}; Configurable cfgEvtRCTFlagCheckerLimitAcceptAsBad{"cfgEvtRCTFlagCheckerLimitAcceptAsBad", true, "Evt sel: RCT flag checker treat Limited Acceptance As Bad"}; @@ -273,14 +285,7 @@ struct Chargedkstaranalysis { Configurable cKstarMinRap{"cKstarMinRap", -0.5, "Kstar minimum rapidity"}; } kstarCutCfgs; - struct : ConfigurableGroup { - // For rotational background - Configurable cFillRotBkg{"cFillRotBkg", true, "Fill rotated background"}; - Configurable confMinRot{"confMinRot", 5.0 * o2::constants::math::PI / 6.0, "Minimum of rotation"}; - Configurable confMaxRot{"confMaxRot", 7.0 * o2::constants::math::PI / 6.0, "Maximum of rotation"}; - Configurable nBkgRotations{"nBkgRotations", 9, "Number of rotated copies (background) per each original candidate"}; - } rotBkgEstCfgs; - + Configurable isQaRequired{"isQaRequired", false, "Fill QA plots"}; float centrality; // PDG code @@ -341,224 +346,235 @@ struct Chargedkstaranalysis { // THnSparse AxisSpec mcLabelAxis = {5, -0.5, 4.5, "MC Label"}; - histos.add("hEvtSelInfo", "hEvtSelInfo", kTH1F, {{5, 0, 5.0}}); - auto hCutFlow = histos.get(HIST("hEvtSelInfo")); - hCutFlow->GetXaxis()->SetBinLabel(1, "All Events"); - hCutFlow->GetXaxis()->SetBinLabel(2, "coll cuts"); - hCutFlow->GetXaxis()->SetBinLabel(3, "rctChecker"); - hCutFlow->GetXaxis()->SetBinLabel(4, "Multiplicity"); - hCutFlow->GetXaxis()->SetBinLabel(5, "IsINELgt0"); - - constexpr int kNTrackCuts = 22; - - histos.add("QA/hTrackCutFlow", "Track cut flow", kTH1I, {{kNTrackCuts, 0.5, kNTrackCuts + 0.5}}); - - auto hTrackCutFlow = histos.get(HIST("QA/hTrackCutFlow")); - - int bin = 1; - hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "All tracks"); - hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "pT min"); - hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "|eta| max"); - hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "ITS clusters"); - hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "TPC clusters"); - hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "TPC crossed rows ratio"); - hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "ITS chi2/Ncl"); - hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "TPC chi2/Ncl"); - hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "Has ITS"); - hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "Has TPC"); - hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "Has TOF"); - hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "ITS refit"); - hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "TPC refit"); - hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "PV contributor"); - hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "Global w/o DCA"); - hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "Global track"); - hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "Primary track"); - hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "DCAxy max"); - hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "DCAz max"); - hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "pT-dep DCAxy"); - hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "pT-dep DCAz"); - - constexpr int kNK0sCuts = 14; - int iK0sbin = 1; - histos.add("QA/K0sCutCheck", "K0s cut flow", kTH1I, {{kNK0sCuts, 0.5, kNK0sCuts + 0.5}}); - auto hK0sCut = histos.get(HIST("QA/K0sCutCheck")); - hK0sCut->GetXaxis()->SetBinLabel(iK0sbin++, "All PASS"); - hK0sCut->GetXaxis()->SetBinLabel(iK0sbin++, "DauDCA>max"); - hK0sCut->GetXaxis()->SetBinLabel(iK0sbin++, "PosDCAtoPVGetXaxis()->SetBinLabel(iK0sbin++, "NegDCAtoPVGetXaxis()->SetBinLabel(iK0sbin++, "pTGetXaxis()->SetBinLabel(iK0sbin++, "|y|>max"); - hK0sCut->GetXaxis()->SetBinLabel(iK0sbin++, "Rmax"); - hK0sCut->GetXaxis()->SetBinLabel(iK0sbin++, "DCAtoPV>max"); - hK0sCut->GetXaxis()->SetBinLabel(iK0sbin++, "cosPAGetXaxis()->SetBinLabel(iK0sbin++, "ctau>max"); - hK0sCut->GetXaxis()->SetBinLabel(iK0sbin++, "qtarmGetXaxis()->SetBinLabel(iK0sbin++, "|M(K0s)-m0|>win"); - hK0sCut->GetXaxis()->SetBinLabel(iK0sbin++, "cross-mass veto"); - - histos.add("QA/before/CentDist", "Centrality distribution", {HistType::kTH1D, {centAxis}}); - histos.add("QA/before/CentDist1", "Centrality distribution", o2::framework::kTH1F, {{110, 0, 110}}); - histos.add("QA/before/VtxZ", "Centrality distribution", {HistType::kTH1D, {vtxzAxis}}); - histos.add("QA/before/hEvent", "Number of Events", HistType::kTH1F, {{1, 0.5, 1.5}}); - - histos.add("QA/trkbpionTPCPIDME", "TPC PID of bachelor pion candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); - - // Bachelor pion - histos.add("QA/before/trkbpionDCAxy", "DCAxy distribution of bachelor pion candidates", HistType::kTH1D, {dcaxyAxis}); - histos.add("QA/before/trkbpionDCAz", "DCAz distribution of bachelor pion candidates", HistType::kTH1D, {dcazAxis}); - histos.add("QA/before/trkbpionpT", "pT distribution of bachelor pion candidates", HistType::kTH1D, {ptAxisQA}); - histos.add("QA/before/trkbpionTPCPID", "TPC PID of bachelor pion candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); - histos.add("QA/before/trkbpionTOFPID", "TOF PID of bachelor pion candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); - histos.add("QA/before/trkbpionTPCTOFPID", "TPC-TOF PID map of bachelor pion candidates", HistType::kTH2D, {pidQAAxis, pidQAAxis}); - - histos.add("QA/after/trkbpionDCAxy", "DCAxy distribution of bachelor pion candidates", HistType::kTH1D, {dcaxyAxis}); - histos.add("QA/after/trkbpionDCAz", "DCAz distribution of bachelor pion candidates", HistType::kTH1D, {dcazAxis}); - histos.add("QA/after/trkbpionpT", "pT distribution of bachelor pion candidates", HistType::kTH1D, {ptAxisQA}); - histos.add("QA/after/trkbpionTPCPID", "TPC PID of bachelor pion candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); - histos.add("QA/after/trkbpionTOFPID", "TOF PID of bachelor pion candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); - histos.add("QA/after/trkbpionTPCTOFPID", "TPC-TOF PID map of bachelor pion candidates", HistType::kTH2D, {pidQAAxis, pidQAAxis}); - - // Secondary pion 1 - histos.add("QA/before/trkppionTPCPID", "TPC PID of secondary pion 1 (positive) candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); - histos.add("QA/before/trkppionTOFPID", "TOF PID of secondary pion 1 (positive) candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); - histos.add("QA/before/trkppionTPCTOFPID", "TPC-TOF PID map of secondary pion 1 (positive) candidates", HistType::kTH2D, {pidQAAxis, pidQAAxis}); - histos.add("QA/before/trkppionpT", "pT distribution of secondary pion 1 (positive) candidates", HistType::kTH1D, {ptAxisQA}); - histos.add("QA/before/trkppionDCAxy", "DCAxy distribution of secondary pion 1 (positive) candidates", HistType::kTH1D, {dcaxyAxis}); - histos.add("QA/before/trkppionDCAz", "DCAz distribution of secondary pion 1 (positive) candidates", HistType::kTH1D, {dcazAxis}); - - histos.add("QA/after/trkppionTPCPID", "TPC PID of secondary pion 1 (positive) candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); - histos.add("QA/after/trkppionTOFPID", "TOF PID of secondary pion 1 (positive) candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); - histos.add("QA/after/trkppionTPCTOFPID", "TPC-TOF PID map of secondary pion 1 (positive) candidates", HistType::kTH2D, {pidQAAxis, pidQAAxis}); - histos.add("QA/after/trkppionpT", "pT distribution of secondary pion 1 (positive) candidates", HistType::kTH1D, {ptAxisQA}); - histos.add("QA/after/trkppionDCAxy", "DCAxy distribution of secondary pion 1 (positive) candidates", HistType::kTH1D, {dcaxyAxis}); - histos.add("QA/after/trkppionDCAz", "DCAz distribution of secondary pion 1 (positive) candidates", HistType::kTH1D, {dcazAxis}); - - // Secondary pion 2 - histos.add("QA/before/trknpionTPCPID", "TPC PID of secondary pion 2 (negative) candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); - histos.add("QA/before/trknpionTOFPID", "TOF PID of secondary pion 2 (negative) candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); - histos.add("QA/before/trknpionTPCTOFPID", "TPC-TOF PID map of secondary pion 2 (negative) candidates", HistType::kTH2D, {pidQAAxis, pidQAAxis}); - histos.add("QA/before/trknpionpT", "pT distribution of secondary pion 2 (negative) candidates", HistType::kTH1D, {ptAxisQA}); - histos.add("QA/before/trknpionDCAxy", "DCAxy distribution of secondary pion 2 (negative) candidates", HistType::kTH1D, {dcaxyAxis}); - histos.add("QA/before/trknpionDCAz", "DCAz distribution of secondary pion 2 (negative) candidates", HistType::kTH1D, {dcazAxis}); - - histos.add("QA/after/trknpionTPCPID", "TPC PID of secondary pion 2 (negative) candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); - histos.add("QA/after/trknpionTOFPID", "TOF PID of secondary pion 2 (negative) candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); - histos.add("QA/after/trknpionTPCTOFPID", "TPC-TOF PID map of secondary pion 2 (negative) candidates", HistType::kTH2D, {pidQAAxis, pidQAAxis}); - histos.add("QA/after/trknpionpT", "pT distribution of secondary pion 2 (negative) candidates", HistType::kTH1D, {ptAxisQA}); - histos.add("QA/after/trknpionDCAxy", "DCAxy distribution of secondary pion 2 (negative) candidates", HistType::kTH1D, {dcaxyAxis}); - histos.add("QA/after/trknpionDCAz", "DCAz distribution of secondary pion 2 (negative) candidates", HistType::kTH1D, {dcazAxis}); - - // K0s - histos.add("QA/before/hDauDCASecondary", "DCA of daughters of secondary resonance", HistType::kTH1D, {dcaAxis}); - histos.add("QA/before/hDauPosDCAtoPVSecondary", "Pos DCA to PV of daughters secondary resonance", HistType::kTH1D, {dcaAxis}); - histos.add("QA/before/hDauNegDCAtoPVSecondary", "Neg DCA to PV of daughters secondary resonance", HistType::kTH1D, {dcaAxis}); - histos.add("QA/before/hpT_Secondary", "pT distribution of secondary resonance", HistType::kTH1D, {ptAxisQA}); - histos.add("QA/before/hy_Secondary", "Rapidity distribution of secondary resonance", HistType::kTH1D, {yAxis}); - histos.add("QA/before/hRadiusSecondary", "Radius distribution of secondary resonance", HistType::kTH1D, {radiusAxis}); - histos.add("QA/before/hCPASecondary", "Cosine pointing angle distribution of secondary resonance", HistType::kTH1D, {cpaAxis}); - histos.add("QA/before/hDCAtoPVSecondary", "DCA to PV distribution of secondary resonance", HistType::kTH1D, {dcaAxis}); - histos.add("QA/before/hPropTauSecondary", "Proper Lifetime distribution of secondary resonance", HistType::kTH1D, {tauAxis}); - histos.add("QA/before/hPtAsymSecondary", "pT asymmetry distribution of secondary resonance", HistType::kTH1D, {AxisSpec{100, -1, 1, "Pair asymmetry"}}); - histos.add("QA/before/hInvmassSecondary", "Invariant mass of unlike-sign secondary resonance", HistType::kTH1D, {invMassAxisK0s}); - - histos.add("QA/after/hDauDCASecondary", "DCA of daughters of secondary resonance", HistType::kTH1D, {dcaAxis}); - histos.add("QA/after/hDauPosDCAtoPVSecondary", "Pos DCA to PV of daughters secondary resonance", HistType::kTH1D, {dcaAxis}); - histos.add("QA/after/hDauNegDCAtoPVSecondary", "Neg DCA to PV of daughters secondary resonance", HistType::kTH1D, {dcaAxis}); - histos.add("QA/after/hpT_Secondary", "pT distribution of secondary resonance", HistType::kTH1D, {ptAxisQA}); - histos.add("QA/after/hy_Secondary", "Rapidity distribution of secondary resonance", HistType::kTH1D, {yAxis}); - histos.add("QA/after/hRadiusSecondary", "Radius distribution of secondary resonance", HistType::kTH1D, {radiusAxis}); - histos.add("QA/after/hCPASecondary", "Cosine pointing angle distribution of secondary resonance", HistType::kTH1D, {cpaAxis}); - histos.add("QA/after/hDCAtoPVSecondary", "DCA to PV distribution of secondary resonance", HistType::kTH1D, {dcaAxis}); - histos.add("QA/after/hPropTauSecondary", "Proper Lifetime distribution of secondary resonance", HistType::kTH1D, {tauAxis}); - histos.add("QA/after/hPtAsymSecondary", "pT asymmetry distribution of secondary resonance", HistType::kTH1D, {AxisSpec{100, -1, 1, "Pair asymmetry"}}); - histos.add("QA/after/hInvmassSecondary", "Invariant mass of unlike-sign secondary resonance", HistType::kTH1D, {invMassAxisK0s}); - - // Kstar - // Invariant mass nSparse - histos.add("QA/before/KstarRapidity", "Rapidity distribution of chK(892)", HistType::kTH1D, {yAxis}); - histos.add("hInvmass_Kstar", "Invariant mass of unlike-sign chK(892)", HistType::kTHnSparseD, {centAxis, ptAxis, invMassAxisReso}); - histos.add("hInvmass_KstarME", "Invariant mass of unlike-sign chK(892)ME", HistType::kTHnSparseD, {centAxis, ptAxis, invMassAxisReso}); - histos.add("hInvmass_KstarRotated", "Invariant mass of unlike-sign chK(892)Rota", HistType::kTHnSparseD, {centAxis, ptAxis, invMassAxisReso}); - - // Mass QA (quick check) - histos.add("QA/before/kstarinvmass", "Invariant mass of unlike-sign chK(892)", HistType::kTH1D, {invMassAxisReso}); - histos.add("QA/before/kstarinvmass_Mix", "Invariant mass of unlike-sign chK(892) from mixed event", HistType::kTH1D, {invMassAxisReso}); - - histos.add("QA/after/KstarRapidity", "Rapidity distribution of chK(892)", HistType::kTH1D, {yAxis}); - histos.add("QA/after/kstarinvmass", "Invariant mass of unlike-sign chK(892)", HistType::kTH1D, {invMassAxisReso}); - histos.add("QA/after/kstarinvmass_Mix", "Invariant mass of unlike-sign chK(892) from mixed event", HistType::kTH1D, {invMassAxisReso}); - - if (!helicityCfgs.qAOptimisation) { - hChaKstar.add("h3ChaKstarInvMassDS", "h3ChaKstarInvMassDS", kTHnSparseF, {centAxis, ptAxis, invMassAxisReso, thnAxisPOL, thnAxisPhi}, true); - hChaKstar.add("h3ChaKstarInvMassME", "h3ChaKstarInvMassME", kTHnSparseF, {centAxis, ptAxis, invMassAxisReso, thnAxisPOL, thnAxisPhi}, true); - hChaKstar.add("h3ChaKstarInvMassRot", "h3ChaKstarInvMassRot", kTHnSparseF, {centAxis, ptAxis, invMassAxisReso, thnAxisPOL, thnAxisPhi}, true); - } + if (!doprocessMC) { + histos.add("hEvtSelInfo", "hEvtSelInfo", kTH1F, {{5, 0, 5.0}}); + auto hCutFlow = histos.get(HIST("hEvtSelInfo")); + hCutFlow->GetXaxis()->SetBinLabel(1, "All Events"); + hCutFlow->GetXaxis()->SetBinLabel(2, "coll cuts"); + hCutFlow->GetXaxis()->SetBinLabel(3, "rctChecker"); + hCutFlow->GetXaxis()->SetBinLabel(4, "Multiplicity"); + hCutFlow->GetXaxis()->SetBinLabel(5, "IsINELgt0"); + if (isQaRequired) { + constexpr int kNTrackCuts = 22; + + histos.add("QA/hTrackCutFlow", "Track cut flow", kTH1I, {{kNTrackCuts, 0.5, kNTrackCuts + 0.5}}); + + auto hTrackCutFlow = histos.get(HIST("QA/hTrackCutFlow")); + + int bin = 1; + hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "All tracks"); + hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "pT min"); + hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "|eta| max"); + hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "ITS clusters"); + hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "TPC clusters"); + hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "TPC crossed rows ratio"); + hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "ITS chi2/Ncl"); + hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "TPC chi2/Ncl"); + hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "Has ITS"); + hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "Has TPC"); + hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "Has TOF"); + hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "ITS refit"); + hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "TPC refit"); + hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "PV contributor"); + hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "Global w/o DCA"); + hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "Global track"); + hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "Primary track"); + hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "DCAxy max"); + hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "DCAz max"); + hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "pT-dep DCAxy"); + hTrackCutFlow->GetXaxis()->SetBinLabel(bin++, "pT-dep DCAz"); + + constexpr int kNK0sCuts = 14; + int iK0sbin = 1; + histos.add("QA/K0sCutCheck", "K0s cut flow", kTH1I, {{kNK0sCuts, 0.5, kNK0sCuts + 0.5}}); + auto hK0sCut = histos.get(HIST("QA/K0sCutCheck")); + hK0sCut->GetXaxis()->SetBinLabel(iK0sbin++, "All PASS"); + hK0sCut->GetXaxis()->SetBinLabel(iK0sbin++, "DauDCA>max"); + hK0sCut->GetXaxis()->SetBinLabel(iK0sbin++, "PosDCAtoPVGetXaxis()->SetBinLabel(iK0sbin++, "NegDCAtoPVGetXaxis()->SetBinLabel(iK0sbin++, "pTGetXaxis()->SetBinLabel(iK0sbin++, "|y|>max"); + hK0sCut->GetXaxis()->SetBinLabel(iK0sbin++, "Rmax"); + hK0sCut->GetXaxis()->SetBinLabel(iK0sbin++, "DCAtoPV>max"); + hK0sCut->GetXaxis()->SetBinLabel(iK0sbin++, "cosPAGetXaxis()->SetBinLabel(iK0sbin++, "ctau>max"); + hK0sCut->GetXaxis()->SetBinLabel(iK0sbin++, "qtarmGetXaxis()->SetBinLabel(iK0sbin++, "|M(K0s)-m0|>win"); + hK0sCut->GetXaxis()->SetBinLabel(iK0sbin++, "cross-mass veto"); + + histos.add("QA/before/CentDist", "Centrality distribution", {HistType::kTH1D, {centAxis}}); + histos.add("QA/before/CentDist1", "Centrality distribution", o2::framework::kTH1F, {{110, 0, 110}}); + histos.add("QA/before/VtxZ", "Centrality distribution", {HistType::kTH1D, {vtxzAxis}}); + histos.add("QA/before/hEvent", "Number of Events", HistType::kTH1F, {{1, 0.5, 1.5}}); + + histos.add("QA/trkbpionTPCPIDME", "TPC PID of bachelor pion candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); + + // Bachelor pion + histos.add("QA/before/trkbpionDCAxy", "DCAxy distribution of bachelor pion candidates", HistType::kTH1D, {dcaxyAxis}); + histos.add("QA/before/trkbpionDCAz", "DCAz distribution of bachelor pion candidates", HistType::kTH1D, {dcazAxis}); + histos.add("QA/before/trkbpionpT", "pT distribution of bachelor pion candidates", HistType::kTH1D, {ptAxisQA}); + histos.add("QA/before/trkbpionTPCPID", "TPC PID of bachelor pion candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); + histos.add("QA/before/trkbpionTOFPID", "TOF PID of bachelor pion candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); + histos.add("QA/before/trkbpionTPCTOFPID", "TPC-TOF PID map of bachelor pion candidates", HistType::kTH2D, {pidQAAxis, pidQAAxis}); + + histos.add("QA/after/trkbpionDCAxy", "DCAxy distribution of bachelor pion candidates", HistType::kTH1D, {dcaxyAxis}); + histos.add("QA/after/trkbpionDCAz", "DCAz distribution of bachelor pion candidates", HistType::kTH1D, {dcazAxis}); + histos.add("QA/after/trkbpionpT", "pT distribution of bachelor pion candidates", HistType::kTH1D, {ptAxisQA}); + histos.add("QA/after/trkbpionTPCPID", "TPC PID of bachelor pion candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); + histos.add("QA/after/trkbpionTOFPID", "TOF PID of bachelor pion candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); + histos.add("QA/after/trkbpionTPCTOFPID", "TPC-TOF PID map of bachelor pion candidates", HistType::kTH2D, {pidQAAxis, pidQAAxis}); + + // Secondary pion 1 + histos.add("QA/before/trkppionTPCPID", "TPC PID of secondary pion 1 (positive) candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); + histos.add("QA/before/trkppionTOFPID", "TOF PID of secondary pion 1 (positive) candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); + histos.add("QA/before/trkppionTPCTOFPID", "TPC-TOF PID map of secondary pion 1 (positive) candidates", HistType::kTH2D, {pidQAAxis, pidQAAxis}); + histos.add("QA/before/trkppionpT", "pT distribution of secondary pion 1 (positive) candidates", HistType::kTH1D, {ptAxisQA}); + histos.add("QA/before/trkppionDCAxy", "DCAxy distribution of secondary pion 1 (positive) candidates", HistType::kTH1D, {dcaxyAxis}); + histos.add("QA/before/trkppionDCAz", "DCAz distribution of secondary pion 1 (positive) candidates", HistType::kTH1D, {dcazAxis}); + + histos.add("QA/after/trkppionTPCPID", "TPC PID of secondary pion 1 (positive) candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); + histos.add("QA/after/trkppionTOFPID", "TOF PID of secondary pion 1 (positive) candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); + histos.add("QA/after/trkppionTPCTOFPID", "TPC-TOF PID map of secondary pion 1 (positive) candidates", HistType::kTH2D, {pidQAAxis, pidQAAxis}); + histos.add("QA/after/trkppionpT", "pT distribution of secondary pion 1 (positive) candidates", HistType::kTH1D, {ptAxisQA}); + histos.add("QA/after/trkppionDCAxy", "DCAxy distribution of secondary pion 1 (positive) candidates", HistType::kTH1D, {dcaxyAxis}); + histos.add("QA/after/trkppionDCAz", "DCAz distribution of secondary pion 1 (positive) candidates", HistType::kTH1D, {dcazAxis}); + + // Secondary pion 2 + histos.add("QA/before/trknpionTPCPID", "TPC PID of secondary pion 2 (negative) candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); + histos.add("QA/before/trknpionTOFPID", "TOF PID of secondary pion 2 (negative) candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); + histos.add("QA/before/trknpionTPCTOFPID", "TPC-TOF PID map of secondary pion 2 (negative) candidates", HistType::kTH2D, {pidQAAxis, pidQAAxis}); + histos.add("QA/before/trknpionpT", "pT distribution of secondary pion 2 (negative) candidates", HistType::kTH1D, {ptAxisQA}); + histos.add("QA/before/trknpionDCAxy", "DCAxy distribution of secondary pion 2 (negative) candidates", HistType::kTH1D, {dcaxyAxis}); + histos.add("QA/before/trknpionDCAz", "DCAz distribution of secondary pion 2 (negative) candidates", HistType::kTH1D, {dcazAxis}); + + histos.add("QA/after/trknpionTPCPID", "TPC PID of secondary pion 2 (negative) candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); + histos.add("QA/after/trknpionTOFPID", "TOF PID of secondary pion 2 (negative) candidates", HistType::kTH2D, {ptAxisQA, pidQAAxis}); + histos.add("QA/after/trknpionTPCTOFPID", "TPC-TOF PID map of secondary pion 2 (negative) candidates", HistType::kTH2D, {pidQAAxis, pidQAAxis}); + histos.add("QA/after/trknpionpT", "pT distribution of secondary pion 2 (negative) candidates", HistType::kTH1D, {ptAxisQA}); + histos.add("QA/after/trknpionDCAxy", "DCAxy distribution of secondary pion 2 (negative) candidates", HistType::kTH1D, {dcaxyAxis}); + histos.add("QA/after/trknpionDCAz", "DCAz distribution of secondary pion 2 (negative) candidates", HistType::kTH1D, {dcazAxis}); + + // K0s + histos.add("QA/before/hDauDCASecondary", "DCA of daughters of secondary resonance", HistType::kTH1D, {dcaAxis}); + histos.add("QA/before/hDauPosDCAtoPVSecondary", "Pos DCA to PV of daughters secondary resonance", HistType::kTH1D, {dcaAxis}); + histos.add("QA/before/hDauNegDCAtoPVSecondary", "Neg DCA to PV of daughters secondary resonance", HistType::kTH1D, {dcaAxis}); + histos.add("QA/before/hpT_Secondary", "pT distribution of secondary resonance", HistType::kTH1D, {ptAxisQA}); + histos.add("QA/before/hy_Secondary", "Rapidity distribution of secondary resonance", HistType::kTH1D, {yAxis}); + histos.add("QA/before/hRadiusSecondary", "Radius distribution of secondary resonance", HistType::kTH1D, {radiusAxis}); + histos.add("QA/before/hCPASecondary", "Cosine pointing angle distribution of secondary resonance", HistType::kTH1D, {cpaAxis}); + histos.add("QA/before/hDCAtoPVSecondary", "DCA to PV distribution of secondary resonance", HistType::kTH1D, {dcaAxis}); + histos.add("QA/before/hPropTauSecondary", "Proper Lifetime distribution of secondary resonance", HistType::kTH1D, {tauAxis}); + histos.add("QA/before/hPtAsymSecondary", "pT asymmetry distribution of secondary resonance", HistType::kTH1D, {AxisSpec{100, -1, 1, "Pair asymmetry"}}); + histos.add("QA/before/hInvmassSecondary", "Invariant mass of unlike-sign secondary resonance", HistType::kTH1D, {invMassAxisK0s}); + + histos.add("QA/after/hDauDCASecondary", "DCA of daughters of secondary resonance", HistType::kTH1D, {dcaAxis}); + histos.add("QA/after/hDauPosDCAtoPVSecondary", "Pos DCA to PV of daughters secondary resonance", HistType::kTH1D, {dcaAxis}); + histos.add("QA/after/hDauNegDCAtoPVSecondary", "Neg DCA to PV of daughters secondary resonance", HistType::kTH1D, {dcaAxis}); + histos.add("QA/after/hpT_Secondary", "pT distribution of secondary resonance", HistType::kTH1D, {ptAxisQA}); + histos.add("QA/after/hy_Secondary", "Rapidity distribution of secondary resonance", HistType::kTH1D, {yAxis}); + histos.add("QA/after/hRadiusSecondary", "Radius distribution of secondary resonance", HistType::kTH1D, {radiusAxis}); + histos.add("QA/after/hCPASecondary", "Cosine pointing angle distribution of secondary resonance", HistType::kTH1D, {cpaAxis}); + histos.add("QA/after/hDCAtoPVSecondary", "DCA to PV distribution of secondary resonance", HistType::kTH1D, {dcaAxis}); + histos.add("QA/after/hPropTauSecondary", "Proper Lifetime distribution of secondary resonance", HistType::kTH1D, {tauAxis}); + histos.add("QA/after/hPtAsymSecondary", "pT asymmetry distribution of secondary resonance", HistType::kTH1D, {AxisSpec{100, -1, 1, "Pair asymmetry"}}); + histos.add("QA/after/hInvmassSecondary", "Invariant mass of unlike-sign secondary resonance", HistType::kTH1D, {invMassAxisK0s}); + + // Kstar + // Invariant mass nSparse + histos.add("QA/before/KstarRapidity", "Rapidity distribution of chK(892)", HistType::kTH1D, {yAxis}); + + // Mass QA (quick check) + histos.add("QA/before/kstarinvmass", "Invariant mass of unlike-sign chK(892)", HistType::kTH1D, {invMassAxisReso}); + histos.add("QA/before/kstarinvmass_Mix", "Invariant mass of unlike-sign chK(892) from mixed event", HistType::kTH1D, {invMassAxisReso}); + + histos.add("QA/after/KstarRapidity", "Rapidity distribution of chK(892)", HistType::kTH1D, {yAxis}); + histos.add("QA/after/kstarinvmass", "Invariant mass of unlike-sign chK(892)", HistType::kTH1D, {invMassAxisReso}); + histos.add("QA/after/kstarinvmass_Mix", "Invariant mass of unlike-sign chK(892) from mixed event", HistType::kTH1D, {invMassAxisReso}); + } + if (!helicityCfgs.qAOptimisation) { + hChaKstar.add("h3ChaKstarInvMassDS", "h3ChaKstarInvMassDS", kTHnSparseF, {centAxis, ptAxis, invMassAxisReso, thnAxisPOL}, true); + hChaKstar.add("h3ChaKstarInvMassRot", "h3ChaKstarInvMassRot", kTHnSparseF, {centAxis, ptAxis, invMassAxisReso, thnAxisPOL}, true); - if (rotBkgEstCfgs.cFillRotBkg) { - histos.add("hRotation", "hRotation", kTH1F, {{360, 0.0, o2::constants::math::TwoPI}}); + // hChaKstar.add("h3ChaKstarInvMassDS", "h3ChaKstarInvMassDS", kTHnSparseF, {centAxis, ptAxis, invMassAxisReso, thnAxisPOL, thnAxisPhi}, true); + // hChaKstar.add("h3ChaKstarInvMassRot", "h3ChaKstarInvMassRot", kTHnSparseF, {centAxis, ptAxis, invMassAxisReso, thnAxisPOL, thnAxisPhi}, true); + } } // MC if (doprocessMC) { + if (isQaRequired) { + histos.add("QA/MC/QACent_woCut", "Centrality without cut", HistType::kTH1F, {centAxis}); + histos.add("QA/MC/QACent_woCentCut", "Centrality without cent cut", HistType::kTH1F, {centAxis}); + histos.add("QA/MC/QACent_wCentCut", "Centrality with cent cut", HistType::kTH1F, {centAxis}); + histos.add("QA/MC/QAvtxz_woCut", "z-vertex without cut", HistType::kTH1F, {vtxzAxis}); + histos.add("QA/MC/QAvtxz_wVtxzCut", "z-vertex with vtxz cut", HistType::kTH1F, {vtxzAxis}); + + histos.add("QAMC/hEvent", "Number of Events", HistType::kTH1F, {{1, 0.5, 1.5}}); + // Bachelor pion + histos.add("QAMC/trkbpionDCAxy", "DCAxy distribution of bachelor pion candidates", HistType::kTH1D, {dcaxyAxis}); + histos.add("QAMC/trkbpionDCAz", "DCAz distribution of bachelor pion candidates", HistType::kTH1D, {dcazAxis}); + histos.add("QAMC/trkbpionpT", "pT distribution of bachelor pion candidates", HistType::kTH1D, {ptAxis}); + histos.add("QAMC/trkbpionTPCPID", "TPC PID of bachelor pion candidates", HistType::kTH2D, {ptAxis, pidQAAxis}); + histos.add("QAMC/trkbpionTOFPID", "TOF PID of bachelor pion candidates", HistType::kTH2D, {ptAxis, pidQAAxis}); + histos.add("QAMC/trkbpionTPCTOFPID", "TPC-TOF PID map of bachelor pion candidates", HistType::kTH2D, {pidQAAxis, pidQAAxis}); + + // Secondary pion 1 + histos.add("QAMC/trkppionDCAxy", "DCAxy distribution of secondary pion 1 (positive) candidates", HistType::kTH1D, {dcaxyAxis}); + histos.add("QAMC/trkppionDCAz", "DCAz distribution of secondary pion 1 (positive) candidates", HistType::kTH1D, {dcazAxis}); + histos.add("QAMC/trkppionpT", "pT distribution of secondary pion 1 (positive) candidates", HistType::kTH1D, {ptAxis}); + histos.add("QAMC/trkppionTPCPID", "TPC PID of secondary pion 1 (positive) candidates", HistType::kTH2D, {ptAxis, pidQAAxis}); + histos.add("QAMC/trkppionTOFPID", "TOF PID of secondary pion 1 (positive) candidates", HistType::kTH2D, {ptAxis, pidQAAxis}); + histos.add("QAMC/trkppionTPCTOFPID", "TPC-TOF PID map of secondary pion 1 (positive) candidates", HistType::kTH2D, {pidQAAxis, pidQAAxis}); + + // Secondary pion 2 + histos.add("QAMC/trknpionTPCPID", "TPC PID of secondary pion 2 (negative) candidates", HistType::kTH2D, {ptAxis, pidQAAxis}); + histos.add("QAMC/trknpionTOFPID", "TOF PID of secondary pion 2 (negative) candidates", HistType::kTH2D, {ptAxis, pidQAAxis}); + histos.add("QAMC/trknpionTPCTOFPID", "TPC-TOF PID map of secondary pion 2 (negative) candidates", HistType::kTH2D, {pidQAAxis, pidQAAxis}); + histos.add("QAMC/trknpionpT", "pT distribution of secondary pion 2 (negative) candidates", HistType::kTH1D, {ptAxis}); + histos.add("QAMC/trknpionDCAxy", "DCAxy distribution of secondary pion 2 (negative) candidates", HistType::kTH1D, {dcaxyAxis}); + histos.add("QAMC/trknpionDCAz", "DCAz distribution of secondary pion 2 (negative) candidates", HistType::kTH1D, {dcazAxis}); + + // Secondary Resonance (K0s cand) + histos.add("QAMC/hDauDCASecondary", "DCA of daughters of secondary resonance", HistType::kTH1D, {dcaAxis}); + histos.add("QAMC/hDauPosDCAtoPVSecondary", "Pos DCA to PV of daughters secondary resonance", HistType::kTH1D, {dcaAxis}); + histos.add("QAMC/hDauNegDCAtoPVSecondary", "Neg DCA to PV of daughters secondary resonance", HistType::kTH1D, {dcaAxis}); + + histos.add("QAMC/hpT_Secondary", "pT distribution of secondary resonance", HistType::kTH1D, {ptAxis}); + histos.add("QAMC/hy_Secondary", "Rapidity distribution of secondary resonance", HistType::kTH1D, {yAxis}); + histos.add("QAMC/hRadiusSecondary", "Radius distribution of secondary resonance", HistType::kTH1D, {radiusAxis}); + histos.add("QAMC/hCPASecondary", "Cosine pointing angle distribution of secondary resonance", HistType::kTH1D, {cpaAxis}); + histos.add("QAMC/hDCAtoPVSecondary", "DCA to PV distribution of secondary resonance", HistType::kTH1D, {dcaAxis}); + histos.add("QAMC/hPropTauSecondary", "Proper Lifetime distribution of secondary resonance", HistType::kTH1D, {tauAxis}); + histos.add("QAMC/hPtAsymSecondary", "pT asymmetry distribution of secondary resonance", HistType::kTH1D, {AxisSpec{100, -1, 1, "Pair asymmetry"}}); + histos.add("QAMC/hInvmassSecondary", "Invariant mass of unlike-sign secondary resonance", HistType::kTH1D, {invMassAxisK0s}); + + // K892 + histos.add("QAMC/KstarOA", "Opening angle of chK(892)", HistType::kTH1D, {AxisSpec{100, 0, 3.14, "Opening angle"}}); + histos.add("QAMC/KstarRapidity", "Rapidity distribution of chK(892)", HistType::kTH1D, {yAxis}); + + histos.add("QAMC/kstarinvmass", "Invariant mass of unlike-sign chK(892)", HistType::kTH1D, {invMassAxisReso}); + histos.add("QAMC/kstarinvmass_noKstar", "Invariant mass of unlike-sign no chK(892)", HistType::kTH1D, {invMassAxisReso}); + } + histos.add("EffKstar/genKstar", "Gen Kstar (|y|<0.5)", HistType::kTH2F, {ptAxis, centAxis}); + histos.add("EffKstar/genKstar_pri", "Gen primary Kstar (|y|<0.5)", HistType::kTH2F, {ptAxis, centAxis}); - histos.add("QAMC/hEvent", "Number of Events", HistType::kTH1F, {{1, 0.5, 1.5}}); - // Bachelor pion - histos.add("QAMC/trkbpionDCAxy", "DCAxy distribution of bachelor pion candidates", HistType::kTH1D, {dcaxyAxis}); - histos.add("QAMC/trkbpionDCAz", "DCAz distribution of bachelor pion candidates", HistType::kTH1D, {dcazAxis}); - histos.add("QAMC/trkbpionpT", "pT distribution of bachelor pion candidates", HistType::kTH1D, {ptAxis}); - histos.add("QAMC/trkbpionTPCPID", "TPC PID of bachelor pion candidates", HistType::kTH2D, {ptAxis, pidQAAxis}); - histos.add("QAMC/trkbpionTOFPID", "TOF PID of bachelor pion candidates", HistType::kTH2D, {ptAxis, pidQAAxis}); - histos.add("QAMC/trkbpionTPCTOFPID", "TPC-TOF PID map of bachelor pion candidates", HistType::kTH2D, {pidQAAxis, pidQAAxis}); - - // Secondary pion 1 - histos.add("QAMC/trkppionDCAxy", "DCAxy distribution of secondary pion 1 (positive) candidates", HistType::kTH1D, {dcaxyAxis}); - histos.add("QAMC/trkppionDCAz", "DCAz distribution of secondary pion 1 (positive) candidates", HistType::kTH1D, {dcazAxis}); - histos.add("QAMC/trkppionpT", "pT distribution of secondary pion 1 (positive) candidates", HistType::kTH1D, {ptAxis}); - histos.add("QAMC/trkppionTPCPID", "TPC PID of secondary pion 1 (positive) candidates", HistType::kTH2D, {ptAxis, pidQAAxis}); - histos.add("QAMC/trkppionTOFPID", "TOF PID of secondary pion 1 (positive) candidates", HistType::kTH2D, {ptAxis, pidQAAxis}); - histos.add("QAMC/trkppionTPCTOFPID", "TPC-TOF PID map of secondary pion 1 (positive) candidates", HistType::kTH2D, {pidQAAxis, pidQAAxis}); - - // Secondary pion 2 - histos.add("QAMC/trknpionTPCPID", "TPC PID of secondary pion 2 (negative) candidates", HistType::kTH2D, {ptAxis, pidQAAxis}); - histos.add("QAMC/trknpionTOFPID", "TOF PID of secondary pion 2 (negative) candidates", HistType::kTH2D, {ptAxis, pidQAAxis}); - histos.add("QAMC/trknpionTPCTOFPID", "TPC-TOF PID map of secondary pion 2 (negative) candidates", HistType::kTH2D, {pidQAAxis, pidQAAxis}); - histos.add("QAMC/trknpionpT", "pT distribution of secondary pion 2 (negative) candidates", HistType::kTH1D, {ptAxis}); - histos.add("QAMC/trknpionDCAxy", "DCAxy distribution of secondary pion 2 (negative) candidates", HistType::kTH1D, {dcaxyAxis}); - histos.add("QAMC/trknpionDCAz", "DCAz distribution of secondary pion 2 (negative) candidates", HistType::kTH1D, {dcazAxis}); - - // Secondary Resonance (K0s cand) - histos.add("QAMC/hDauDCASecondary", "DCA of daughters of secondary resonance", HistType::kTH1D, {dcaAxis}); - histos.add("QAMC/hDauPosDCAtoPVSecondary", "Pos DCA to PV of daughters secondary resonance", HistType::kTH1D, {dcaAxis}); - histos.add("QAMC/hDauNegDCAtoPVSecondary", "Neg DCA to PV of daughters secondary resonance", HistType::kTH1D, {dcaAxis}); - - histos.add("QAMC/hpT_Secondary", "pT distribution of secondary resonance", HistType::kTH1D, {ptAxis}); - histos.add("QAMC/hy_Secondary", "Rapidity distribution of secondary resonance", HistType::kTH1D, {yAxis}); - histos.add("QAMC/hRadiusSecondary", "Radius distribution of secondary resonance", HistType::kTH1D, {radiusAxis}); - histos.add("QAMC/hCPASecondary", "Cosine pointing angle distribution of secondary resonance", HistType::kTH1D, {cpaAxis}); - histos.add("QAMC/hDCAtoPVSecondary", "DCA to PV distribution of secondary resonance", HistType::kTH1D, {dcaAxis}); - histos.add("QAMC/hPropTauSecondary", "Proper Lifetime distribution of secondary resonance", HistType::kTH1D, {tauAxis}); - histos.add("QAMC/hPtAsymSecondary", "pT asymmetry distribution of secondary resonance", HistType::kTH1D, {AxisSpec{100, -1, 1, "Pair asymmetry"}}); - histos.add("QAMC/hInvmassSecondary", "Invariant mass of unlike-sign secondary resonance", HistType::kTH1D, {invMassAxisK0s}); - - // K892 - histos.add("QAMC/KstarOA", "Opening angle of chK(892)", HistType::kTH1D, {AxisSpec{100, 0, 3.14, "Opening angle"}}); - histos.add("QAMC/KstarRapidity", "Rapidity distribution of chK(892)", HistType::kTH1D, {yAxis}); - - histos.add("QAMC/kstarinvmass", "Invariant mass of unlike-sign chK(892)", HistType::kTH1D, {invMassAxisReso}); - histos.add("QAMC/kstarinvmass_noKstar", "Invariant mass of unlike-sign no chK(892)", HistType::kTH1D, {invMassAxisReso}); - - histos.add("hInvmass_Kstar_MC", "Invariant mass of unlike chK(892)", HistType::kTHnSparseD, {centAxis, ptAxis, invMassAxisReso}); - - ccdb->setURL(cfgURL); - ccdbApi.init("http://alice-ccdb.cern.ch"); - ccdb->setCaching(true); - ccdb->setLocalObjectValidityChecking(); - ccdb->setCreatedNotAfter(std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count()); + histos.add("EffKstar/recoKstar", "Kstar Reco matched (final all)", HistType::kTH2F, {ptAxis, centAxis}); + histos.add("MCReco/hInvmass_Kstar_true", "MC-reco truth-tagged chK(892)", HistType::kTHnSparseD, {centAxis, ptAxis, invMassAxisReso}); } + ccdb->setURL(cfgURL); + ccdbApi.init("http://alice-ccdb.cern.ch"); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setCreatedNotAfter(std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count()); + // Print output histograms statistics LOG(info) << "Size of the histograms in chK(892) Analysis Task"; histos.print(); } + + std::unordered_set allowedMcIds; + std::unordered_map centTruthByAllowed; + float lMultiplicity; template float getCentrality(CollisionType const& collision) @@ -575,107 +591,86 @@ struct Chargedkstaranalysis { } } - // Track selection template bool trackCut(TrackType const& track) { int ibin = 1; - histos.fill(HIST("QA/hTrackCutFlow"), ibin++); - // basic track cuts - if (std::abs(track.pt()) < trackCutCfgs.cMinPtcut) - return false; - histos.fill(HIST("QA/hTrackCutFlow"), ibin++); + auto applyCut = [&](bool condition) -> bool { + if (!condition) + return false; + if (!doprocessMC && isQaRequired) + histos.fill(HIST("QA/hTrackCutFlow"), ibin); + ibin++; + return true; + }; - if (std::abs(track.eta()) > trackCutCfgs.cMaxEtacut) - return false; - histos.fill(HIST("QA/hTrackCutFlow"), ibin++); + // First bin (before any cuts) + if (!doprocessMC && isQaRequired) + histos.fill(HIST("QA/hTrackCutFlow"), ibin++); - if (track.itsNCls() < trackCutCfgs.cfgITScluster) + // Cuts + if (!applyCut(std::abs(track.pt()) >= trackCutCfgs.cMinPtcut)) return false; - histos.fill(HIST("QA/hTrackCutFlow"), ibin++); - - if (track.tpcNClsFound() < trackCutCfgs.cfgTPCcluster) + if (!applyCut(std::abs(track.eta()) <= trackCutCfgs.cMaxEtacut)) return false; - histos.fill(HIST("QA/hTrackCutFlow"), ibin++); - - if (track.tpcCrossedRowsOverFindableCls() < trackCutCfgs.cfgRatioTPCRowsOverFindableCls) + if (!applyCut(track.itsNCls() >= trackCutCfgs.cfgITScluster)) return false; - histos.fill(HIST("QA/hTrackCutFlow"), ibin++); - - if (track.itsChi2NCl() >= trackCutCfgs.cfgITSChi2NCl) + if (!applyCut(track.tpcNClsFound() >= trackCutCfgs.cfgTPCcluster)) return false; - histos.fill(HIST("QA/hTrackCutFlow"), ibin++); - - if (track.tpcChi2NCl() >= trackCutCfgs.cfgTPCChi2NCl) + if (!applyCut(track.tpcCrossedRowsOverFindableCls() >= trackCutCfgs.cfgRatioTPCRowsOverFindableCls)) return false; - histos.fill(HIST("QA/hTrackCutFlow"), ibin++); - - if (trackCutCfgs.cfgHasITS && !track.hasITS()) + if (!applyCut(track.itsChi2NCl() < trackCutCfgs.cfgITSChi2NCl)) return false; - histos.fill(HIST("QA/hTrackCutFlow"), ibin++); - - if (trackCutCfgs.cfgHasTPC && !track.hasTPC()) + if (!applyCut(track.tpcChi2NCl() < trackCutCfgs.cfgTPCChi2NCl)) return false; - histos.fill(HIST("QA/hTrackCutFlow"), ibin++); - if (trackCutCfgs.cfgHasTOF && !track.hasTOF()) + if (!applyCut(!trackCutCfgs.cfgHasITS || track.hasITS())) return false; - histos.fill(HIST("QA/hTrackCutFlow"), ibin++); - - if (trackCutCfgs.cfgUseITSRefit && !track.passedITSRefit()) + if (!applyCut(!trackCutCfgs.cfgHasTPC || track.hasTPC())) return false; - histos.fill(HIST("QA/hTrackCutFlow"), ibin++); - - if (trackCutCfgs.cfgUseTPCRefit && !track.passedTPCRefit()) + if (!applyCut(!trackCutCfgs.cfgHasTOF || track.hasTOF())) return false; - histos.fill(HIST("QA/hTrackCutFlow"), ibin++); - if (trackCutCfgs.cfgPVContributor && !track.isPVContributor()) + if (!applyCut(!trackCutCfgs.cfgUseITSRefit || track.passedITSRefit())) return false; - histos.fill(HIST("QA/hTrackCutFlow"), ibin++); - - if (trackCutCfgs.cfgGlobalWoDCATrack && !track.isGlobalTrackWoDCA()) + if (!applyCut(!trackCutCfgs.cfgUseTPCRefit || track.passedTPCRefit())) return false; - histos.fill(HIST("QA/hTrackCutFlow"), ibin++); - if (trackCutCfgs.cfgGlobalTrack && !track.isGlobalTrack()) + if (!applyCut(!trackCutCfgs.cfgPVContributor || track.isPVContributor())) return false; - histos.fill(HIST("QA/hTrackCutFlow"), ibin++); - - if (trackCutCfgs.cfgPrimaryTrack && !track.isPrimaryTrack()) + if (!applyCut(!trackCutCfgs.cfgGlobalWoDCATrack || track.isGlobalTrackWoDCA())) return false; - histos.fill(HIST("QA/hTrackCutFlow"), ibin++); - - if (std::abs(track.dcaXY()) > trackCutCfgs.cMaxbDCArToPVcut) + if (!applyCut(!trackCutCfgs.cfgGlobalTrack || track.isGlobalTrack())) + return false; + if (!applyCut(!trackCutCfgs.cfgPrimaryTrack || track.isPrimaryTrack())) return false; - histos.fill(HIST("QA/hTrackCutFlow"), ibin++); - if (std::abs(track.dcaZ()) > trackCutCfgs.cMaxbDCAzToPVcut) + if (!applyCut(std::abs(track.dcaXY()) <= trackCutCfgs.cMaxbDCArToPVcut)) + return false; + if (!applyCut(std::abs(track.dcaZ()) <= trackCutCfgs.cMaxbDCAzToPVcut)) return false; - histos.fill(HIST("QA/hTrackCutFlow"), ibin++); + // pT dependent DCA XY if (trackCutCfgs.cfgpTdepDCAxyCut) { - if (std::abs(track.dcaXY()) > (0.004 + (0.013 / track.pt()))) + if (!applyCut(std::abs(track.dcaXY()) <= (0.004 + (0.013 / track.pt())))) return false; } else { - if (std::abs(track.dcaXY()) > trackCutCfgs.cfgMaxbDCArToPVcut) + if (!applyCut(std::abs(track.dcaXY()) <= trackCutCfgs.cfgMaxbDCArToPVcut)) return false; } - histos.fill(HIST("QA/hTrackCutFlow"), ibin++); + // pT dependent DCA Z if (trackCutCfgs.cfgpTdepDCAzCut) { - // Tuned on the LHC22f anchored MC LHC23d1d on primary pions. 7 Sigmas of the resolution - if (std::abs(track.dcaZ()) > (0.004 + (0.013 / track.pt()))) + if (!applyCut(std::abs(track.dcaZ()) <= (0.004 + (0.013 / track.pt())))) return false; } else { - if (std::abs(track.dcaZ()) > trackCutCfgs.cfgMaxbDCAzToPVcut) + if (!applyCut(std::abs(track.dcaZ()) <= trackCutCfgs.cfgMaxbDCAzToPVcut)) return false; } - histos.fill(HIST("QA/hTrackCutFlow"), ibin++); + return true; } - // PID selection tools template bool selectionPIDPion(TrackType const& candidate) @@ -703,6 +698,19 @@ struct Chargedkstaranalysis { bool selectionK0s(CollisionType const& collision, K0sType const& candidate) { int ibin = 1; + bool returnFlag = true; + + auto applyCut = [&](bool condition) { + if (!condition) { + returnFlag = false; + } + if (returnFlag && (!doprocessMC && isQaRequired)) { + histos.fill(HIST("QA/K0sCutCheck"), ibin); + } + ibin++; + }; + + // Precompute variables auto dauDCA = std::fabs(candidate.dcaV0daughters()); auto dauPosDCAtoPV = std::fabs(candidate.dcapostopv()); auto dauNegDCAtoPV = std::fabs(candidate.dcanegtopv()); @@ -716,95 +724,25 @@ struct Chargedkstaranalysis { auto mLambda = candidate.mLambda(); auto mALambda = candidate.mAntiLambda(); - bool returnFlag = true; - histos.fill(HIST("QA/K0sCutCheck"), ibin); - if (std::fabs(dauDCA) > secondaryCutsCfgs.cSecondaryDauDCAMax) { - returnFlag = false; - } - ibin++; - if (returnFlag == true) - histos.fill(HIST("QA/K0sCutCheck"), ibin); - if (std::fabs(dauPosDCAtoPV) < secondaryCutsCfgs.cSecondaryDauPosDCAtoPVMin) { - returnFlag = false; - } - ibin++; - if (returnFlag == true) - histos.fill(HIST("QA/K0sCutCheck"), ibin); - - if (std::fabs(dauNegDCAtoPV) < secondaryCutsCfgs.cSecondaryDauNegDCAtoPVMin) { - returnFlag = false; - } - ibin++; - if (returnFlag == true) - histos.fill(HIST("QA/K0sCutCheck"), ibin); - - if (pT < secondaryCutsCfgs.cSecondaryPtMin) { - returnFlag = false; - } - ibin++; - if (returnFlag == true) - histos.fill(HIST("QA/K0sCutCheck"), ibin); - - if (std::fabs(rapidity) > secondaryCutsCfgs.cSecondaryRapidityMax) { - returnFlag = false; - } - ibin++; - if (returnFlag == true) - histos.fill(HIST("QA/K0sCutCheck"), ibin); - - if (v0Radius < secondaryCutsCfgs.cSecondaryRadiusMin || v0Radius > secondaryCutsCfgs.cSecondaryRadiusMax) { - returnFlag = false; - } - ibin++; - if (returnFlag == true) - histos.fill(HIST("QA/K0sCutCheck"), ibin); - - if (dcaToPV > secondaryCutsCfgs.cSecondaryDCAtoPVMax) { - returnFlag = false; - } - ibin++; - if (returnFlag == true) - histos.fill(HIST("QA/K0sCutCheck"), ibin); - - if (cosPA < secondaryCutsCfgs.cSecondaryCosPAMin) { - returnFlag = false; - } - ibin++; - if (returnFlag == true) - histos.fill(HIST("QA/K0sCutCheck"), ibin); - - if (propTauK0s > secondaryCutsCfgs.cSecondaryProperLifetimeMax) { - returnFlag = false; - } - ibin++; - if (returnFlag == true) + if (!doprocessMC && isQaRequired) histos.fill(HIST("QA/K0sCutCheck"), ibin); - if (candidate.qtarm() < secondaryCutsCfgs.cfgSecondaryparamArmenterosCut * std::fabs(candidate.alpha())) { - returnFlag = false; - } - ibin++; - if (returnFlag == true) - histos.fill(HIST("QA/K0sCutCheck"), ibin); + applyCut(dauDCA <= secondaryCutsCfgs.cSecondaryDauDCAMax); + applyCut(dauPosDCAtoPV >= secondaryCutsCfgs.cSecondaryDauPosDCAtoPVMin); + applyCut(dauNegDCAtoPV >= secondaryCutsCfgs.cSecondaryDauNegDCAtoPVMin); + applyCut(pT >= secondaryCutsCfgs.cSecondaryPtMin); + applyCut(std::fabs(rapidity) <= secondaryCutsCfgs.cSecondaryRapidityMax); + applyCut(v0Radius >= secondaryCutsCfgs.cSecondaryRadiusMin && v0Radius <= secondaryCutsCfgs.cSecondaryRadiusMax); + applyCut(dcaToPV <= secondaryCutsCfgs.cSecondaryDCAtoPVMax); + applyCut(cosPA >= secondaryCutsCfgs.cSecondaryCosPAMin); + applyCut(propTauK0s <= secondaryCutsCfgs.cSecondaryProperLifetimeMax); + applyCut(candidate.qtarm() >= secondaryCutsCfgs.cfgSecondaryparamArmenterosCut * std::fabs(candidate.alpha())); + applyCut(std::fabs(mK0s - MassK0Short) <= secondaryCutsCfgs.cSecondaryMassWindow); - if (std::fabs(mK0s - MassK0Short) > secondaryCutsCfgs.cSecondaryMassWindow) { - returnFlag = false; - } - ibin++; - if (returnFlag == true) - histos.fill(HIST("QA/K0sCutCheck"), ibin); + applyCut(!secondaryCutsCfgs.cfgSecondaryCrossMassHypothesisCut || ((std::fabs(mLambda - MassLambda0) >= secondaryCutsCfgs.cfgSecondaryCrossMassCutWindow) && (std::fabs(mALambda - MassLambda0Bar) >= secondaryCutsCfgs.cfgSecondaryCrossMassCutWindow))); - if (secondaryCutsCfgs.cfgSecondaryCrossMassHypothesisCut && - ((std::fabs(mLambda - MassLambda0) < secondaryCutsCfgs.cfgSecondaryCrossMassCutWindow) || (std::fabs(mALambda - MassLambda0Bar) < secondaryCutsCfgs.cfgSecondaryCrossMassCutWindow))) { - returnFlag = false; - } - ibin++; - if (returnFlag == true) - histos.fill(HIST("QA/K0sCutCheck"), ibin); return returnFlag; - - } // selectionK0s - + } template bool isTrueKstar(const TrackTemplate& bTrack, const V0Template& K0scand) { @@ -843,6 +781,53 @@ struct Chargedkstaranalysis { double massPi = o2::constants::physics::MassPionCharged; double massK0s = o2::constants::physics::MassK0Short; + template + bool matchRecoToTruthKstar(V0T const& v0, TrkT const& trk) + { + if (!v0.has_mcParticle() || !trk.has_mcParticle()) + return false; + + auto mcK0s = v0.template mcParticle_as(); + auto mcPi = trk.template mcParticle_as(); + + if (std::abs(mcK0s.pdgCode()) != kPDGK0s) + return false; + if (std::abs(mcPi.pdgCode()) != kPiPlus) + return false; + + MCTrueTrackCandidates::iterator kstarFromPi; + bool havePiKstar = false; + for (const auto& m1 : mcPi.template mothers_as()) { + if (std::abs(m1.pdgCode()) == kKstarPlus) { + kstarFromPi = m1; + havePiKstar = true; + break; + } + } + if (!havePiKstar) { + return false; + } + + bool shareSameKstar = false; + for (const auto& m1 : mcK0s.template mothers_as()) { + if (std::abs(m1.pdgCode()) == kPDGK0) { + for (const auto& m2 : m1.template mothers_as()) { + if (m2.globalIndex() == kstarFromPi.globalIndex()) { + shareSameKstar = true; + break; + } + } + if (shareSameKstar) + break; + } + } + if (!shareSameKstar) { + return false; + } + + return true; + } // matchRecoToTruthKstar + template void fillInvMass(const T& mother, float multiplicity, const T& daughter1, const T& daughter2, bool isMix) { @@ -896,7 +881,7 @@ struct Chargedkstaranalysis { auto cosThetaStarHelicity = mother.Vect().Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(mother.Vect().Mag2())); if (!isMix) { if (std::abs(mother.Rapidity()) < helicityCfgs.rapidityMotherData) { - hChaKstar.fill(HIST("h3ChaKstarInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity, anglePhi); + hChaKstar.fill(HIST("h3ChaKstarInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity); //, anglePhi); } for (int i = 0; i < helicityCfgs.cRotations; i++) { @@ -912,18 +897,14 @@ struct Chargedkstaranalysis { 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); - if (motherRot.Rapidity() < helicityCfgs.rapidityMotherData) - hChaKstar.fill(HIST("h3ChaKstarInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarHelicityRot, phiHelicityRot); - } - } else { - if (std::abs(mother.Rapidity()) < helicityCfgs.rapidityMotherData) { - hChaKstar.fill(HIST("h3ChaKstarInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarHelicity, anglePhi); + if (std::abs(motherRot.Rapidity()) < helicityCfgs.rapidityMotherData) + hChaKstar.fill(HIST("h3ChaKstarInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarHelicityRot); //, phiHelicityRot); } } } else if (helicityCfgs.activateCollinsSoperFrame) { if (!isMix) { if (std::abs(mother.Rapidity()) < helicityCfgs.rapidityMotherData) { - hChaKstar.fill(HIST("h3ChaKstarInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarCS, phiCS); + hChaKstar.fill(HIST("h3ChaKstarInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarCS); //, phiCS); } for (int i = 0; i < helicityCfgs.cRotations; i++) { @@ -940,12 +921,8 @@ struct Chargedkstaranalysis { auto phiCSrot = std::atan2(yAxisCS.Dot(daughterRotCM.Vect().Unit()), xAxisCS.Dot(daughterRotCM.Vect().Unit())); phiCSrot = RecoDecay::constrainAngle(phiCSrot, 0.0); - if (motherRot.Rapidity() < helicityCfgs.rapidityMotherData) - hChaKstar.fill(HIST("h3ChaKstarInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarCSrot, phiCSrot); - } - } else { - if (std::abs(mother.Rapidity()) < helicityCfgs.rapidityMotherData) { - hChaKstar.fill(HIST("h3ChaKstarInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarCS, phiCS); + if (std::abs(motherRot.Rapidity()) < helicityCfgs.rapidityMotherData) + hChaKstar.fill(HIST("h3ChaKstarInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarCSrot); //, phiCSrot); } } } else if (helicityCfgs.activateProductionFrame) { @@ -953,38 +930,30 @@ struct Chargedkstaranalysis { auto cosThetaProduction = normalVec.Dot(fourVecDauCM.Vect()) / (std::sqrt(fourVecDauCM.Vect().Mag2()) * std::sqrt(normalVec.Mag2())); if (!isMix) { if (std::abs(mother.Rapidity()) < helicityCfgs.rapidityMotherData) { - hChaKstar.fill(HIST("h3ChaKstarInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaProduction, anglePhi); + hChaKstar.fill(HIST("h3ChaKstarInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaProduction); //, anglePhi); } for (int i = 0; i < helicityCfgs.cRotations; i++) { theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / helicityCfgs.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / helicityCfgs.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()); if (std::abs(motherRot.Rapidity()) < helicityCfgs.rapidityMotherData) { - hChaKstar.fill(HIST("h3ChaKstarInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaProduction, anglePhi); + hChaKstar.fill(HIST("h3ChaKstarInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaProduction); //, anglePhi); } } - } else { - if (std::abs(mother.Rapidity()) < helicityCfgs.rapidityMotherData) { - hChaKstar.fill(HIST("h3ChaKstarInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaProduction, anglePhi); - } } } else if (helicityCfgs.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()) < helicityCfgs.rapidityMotherData) { - hChaKstar.fill(HIST("h3ChaKstarInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarBeam, anglePhi); + hChaKstar.fill(HIST("h3ChaKstarInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarBeam); //, anglePhi); } for (int i = 0; i < helicityCfgs.cRotations; i++) { theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / helicityCfgs.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / helicityCfgs.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()); if (std::abs(motherRot.Rapidity()) < helicityCfgs.rapidityMotherData) { - hChaKstar.fill(HIST("h3ChaKstarInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarBeam, anglePhi); + hChaKstar.fill(HIST("h3ChaKstarInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarBeam); //, anglePhi); } } - } else { - if (std::abs(mother.Rapidity()) < helicityCfgs.rapidityMotherData) { - hChaKstar.fill(HIST("h3ChaKstarInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarBeam, anglePhi); - } } } else if (helicityCfgs.activateRandomFrame) { auto phiRandom = gRandom->Uniform(0.f, constants::math::TwoPI); @@ -994,19 +963,15 @@ struct Chargedkstaranalysis { auto cosThetaStarRandom = randomVec.Dot(fourVecDauCM.Vect()) / std::sqrt(fourVecDauCM.Vect().Mag2()); if (!isMix) { if (std::abs(mother.Rapidity()) < helicityCfgs.rapidityMotherData) { - hChaKstar.fill(HIST("h3ChaKstarInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, phiRandom); + hChaKstar.fill(HIST("h3ChaKstarInvMassDS"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom); //, phiRandom); } for (int i = 0; i < helicityCfgs.cRotations; i++) { theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / helicityCfgs.rotationalCut, o2::constants::math::PI + o2::constants::math::PI / helicityCfgs.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()); if (std::abs(motherRot.Rapidity()) < helicityCfgs.rapidityMotherData) { - hChaKstar.fill(HIST("h3ChaKstarInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarRandom, phiRandom); + hChaKstar.fill(HIST("h3ChaKstarInvMassRot"), multiplicity, motherRot.Pt(), motherRot.M(), cosThetaStarRandom); //, phiRandom); } } - } else { - if (std::abs(mother.Rapidity()) < helicityCfgs.rapidityMotherData) { - hChaKstar.fill(HIST("h3ChaKstarInvMassME"), multiplicity, mother.Pt(), mother.M(), cosThetaStarRandom, phiRandom); - } } } // } @@ -1015,69 +980,82 @@ struct Chargedkstaranalysis { template void fillHistograms(const CollisionType& collision, const TracksType& dTracks1, const TracksTypeK0s& dTracks2) { - histos.fill(HIST("QA/before/CentDist"), collision.centFT0M()); - histos.fill(HIST("QA/before/CentDist1"), collision.centFT0M()); - ROOT::Math::PxPyPzMVector lDecayDaughter1, lDecayDaughter2, lResoSecondary, lDecayDaughter_bach, lResoKstar, chargeKstarrot; - std::vector trackIndicies = {}; - std::vector k0sIndicies = {}; + if (!doprocessMC && isQaRequired) { + histos.fill(HIST("QA/before/CentDist"), collision.centFT0M()); + histos.fill(HIST("QA/before/CentDist1"), collision.centFT0M()); + } + + ROOT::Math::PxPyPzMVector lResoSecondary, lDecayDaughter_bach, lResoKstar, chargeKstarrot; + std::vector trackIndicies; + std::vector k0sIndicies; + + // ========================= + // Bachelor tracks + // ========================= for (const auto& bTrack : dTracks1) { + auto trkbpt = bTrack.pt(); auto istrkbhasTOF = bTrack.hasTOF(); auto trkbNSigmaPiTPC = bTrack.tpcNSigmaPi(); - auto trkbNSigmaPiTOF = (istrkbhasTOF) ? bTrack.tofNSigmaPi() : -999.; + auto trkbNSigmaPiTOF = istrkbhasTOF ? bTrack.tofNSigmaPi() : -999.; if constexpr (!IsMix) { - // Bachelor pion QA plots - histos.fill(HIST("QA/before/trkbpionTPCPID"), trkbpt, trkbNSigmaPiTPC); - if (istrkbhasTOF) { - histos.fill(HIST("QA/before/trkbpionTOFPID"), trkbpt, trkbNSigmaPiTOF); - histos.fill(HIST("QA/before/trkbpionTPCTOFPID"), trkbNSigmaPiTPC, trkbNSigmaPiTOF); + if (!doprocessMC && isQaRequired) { + histos.fill(HIST("QA/before/trkbpionTPCPID"), trkbpt, trkbNSigmaPiTPC); + if (istrkbhasTOF) { + histos.fill(HIST("QA/before/trkbpionTOFPID"), trkbpt, trkbNSigmaPiTOF); + histos.fill(HIST("QA/before/trkbpionTPCTOFPID"), trkbNSigmaPiTPC, trkbNSigmaPiTOF); + } + histos.fill(HIST("QA/before/trkbpionpT"), trkbpt); + histos.fill(HIST("QA/before/trkbpionDCAxy"), bTrack.dcaXY()); + histos.fill(HIST("QA/before/trkbpionDCAz"), bTrack.dcaZ()); } - histos.fill(HIST("QA/before/trkbpionpT"), trkbpt); - histos.fill(HIST("QA/before/trkbpionDCAxy"), bTrack.dcaXY()); - histos.fill(HIST("QA/before/trkbpionDCAz"), bTrack.dcaZ()); } else { - - histos.fill(HIST("QA/trkbpionTPCPIDME"), trkbpt, trkbNSigmaPiTPC); + if (!doprocessMC && isQaRequired) { + histos.fill(HIST("QA/trkbpionTPCPIDME"), trkbpt, trkbNSigmaPiTPC); + } } - if (!trackCut(bTrack)) - continue; - if (!selectionPIDPion(bTrack)) + if (!trackCut(bTrack) || !selectionPIDPion(bTrack)) continue; if constexpr (!IsMix) { - // Bachelor pion QA plots after applying cuts - histos.fill(HIST("QA/after/trkbpionTPCPID"), trkbpt, trkbNSigmaPiTPC); - if (istrkbhasTOF) { - histos.fill(HIST("QA/after/trkbpionTOFPID"), trkbpt, trkbNSigmaPiTOF); - histos.fill(HIST("QA/after/trkbpionTPCTOFPID"), trkbNSigmaPiTPC, trkbNSigmaPiTOF); + if (!doprocessMC && isQaRequired) { + histos.fill(HIST("QA/after/trkbpionTPCPID"), trkbpt, trkbNSigmaPiTPC); + if (istrkbhasTOF) { + histos.fill(HIST("QA/after/trkbpionTOFPID"), trkbpt, trkbNSigmaPiTOF); + histos.fill(HIST("QA/after/trkbpionTPCTOFPID"), trkbNSigmaPiTPC, trkbNSigmaPiTOF); + } + histos.fill(HIST("QA/after/trkbpionpT"), trkbpt); + histos.fill(HIST("QA/after/trkbpionDCAxy"), bTrack.dcaXY()); + histos.fill(HIST("QA/after/trkbpionDCAz"), bTrack.dcaZ()); } - histos.fill(HIST("QA/after/trkbpionpT"), trkbpt); - histos.fill(HIST("QA/after/trkbpionDCAxy"), bTrack.dcaXY()); - histos.fill(HIST("QA/after/trkbpionDCAz"), bTrack.dcaZ()); } + trackIndicies.push_back(bTrack.index()); } + // ========================= + // K0s loop + // ========================= for (const auto& K0scand : dTracks2) { - auto posDauTrack = K0scand.template posTrack_as(); - auto negDauTrack = K0scand.template negTrack_as(); - - /// Daughters - // Positve pion - auto trkppt = posDauTrack.pt(); - auto istrkphasTOF = posDauTrack.hasTOF(); - auto trkpNSigmaPiTPC = posDauTrack.tpcNSigmaPi(); - auto trkpNSigmaPiTOF = (istrkphasTOF) ? posDauTrack.tofNSigmaPi() : -999.; - // Negative pion - auto trknpt = negDauTrack.pt(); - auto istrknhasTOF = negDauTrack.hasTOF(); - auto trknNSigmaPiTPC = negDauTrack.tpcNSigmaPi(); - auto trknNSigmaPiTOF = (istrknhasTOF) ? negDauTrack.tofNSigmaPi() : -999.; - - /// K0s + + auto pos = K0scand.template posTrack_as(); + auto neg = K0scand.template negTrack_as(); + + auto trkppt = pos.pt(); + auto trknpt = neg.pt(); + + auto istrkphasTOF = pos.hasTOF(); + auto istrknhasTOF = neg.hasTOF(); + + auto trkpNSigmaPiTPC = pos.tpcNSigmaPi(); + auto trknNSigmaPiTPC = neg.tpcNSigmaPi(); + + auto trkpNSigmaPiTOF = istrkphasTOF ? pos.tofNSigmaPi() : -999.; + auto trknNSigmaPiTOF = istrknhasTOF ? neg.tofNSigmaPi() : -999.; + auto trkkDauDCA = K0scand.dcaV0daughters(); auto trkkDauDCAPostoPV = K0scand.dcapostopv(); auto trkkDauDCANegtoPV = K0scand.dcanegtopv(); @@ -1089,139 +1067,99 @@ struct Chargedkstaranalysis { auto trkkMass = K0scand.mK0Short(); if constexpr (!IsMix) { - // Seconddary QA plots - histos.fill(HIST("QA/before/trkppionTPCPID"), trkppt, trkpNSigmaPiTPC); - if (istrkphasTOF) { - histos.fill(HIST("QA/before/trkppionTOFPID"), trkppt, trkpNSigmaPiTOF); - histos.fill(HIST("QA/before/trkppionTPCTOFPID"), trkpNSigmaPiTPC, trkpNSigmaPiTOF); - } - histos.fill(HIST("QA/before/trkppionpT"), trkppt); - histos.fill(HIST("QA/before/trkppionDCAxy"), posDauTrack.dcaXY()); - histos.fill(HIST("QA/before/trkppionDCAz"), posDauTrack.dcaZ()); - - histos.fill(HIST("QA/before/trknpionTPCPID"), trknpt, trknNSigmaPiTPC); - if (istrknhasTOF) { - histos.fill(HIST("QA/before/trknpionTOFPID"), trknpt, trknNSigmaPiTOF); - histos.fill(HIST("QA/before/trknpionTPCTOFPID"), trknNSigmaPiTPC, trknNSigmaPiTOF); + if (!doprocessMC && isQaRequired) { + // positive pion + histos.fill(HIST("QA/before/trkppionTPCPID"), trkppt, trkpNSigmaPiTPC); + if (istrkphasTOF) { + histos.fill(HIST("QA/before/trkppionTOFPID"), trkppt, trkpNSigmaPiTOF); + histos.fill(HIST("QA/before/trkppionTPCTOFPID"), trkpNSigmaPiTPC, trkpNSigmaPiTOF); + } + + // negative pion + histos.fill(HIST("QA/before/trknpionTPCPID"), trknpt, trknNSigmaPiTPC); + if (istrknhasTOF) { + histos.fill(HIST("QA/before/trknpionTOFPID"), trknpt, trknNSigmaPiTOF); + histos.fill(HIST("QA/before/trknpionTPCTOFPID"), trknNSigmaPiTPC, trknNSigmaPiTOF); + } + + // K0s + histos.fill(HIST("QA/before/hDauDCASecondary"), trkkDauDCA); + histos.fill(HIST("QA/before/hDauPosDCAtoPVSecondary"), trkkDauDCAPostoPV); + histos.fill(HIST("QA/before/hDauNegDCAtoPVSecondary"), trkkDauDCANegtoPV); + histos.fill(HIST("QA/before/hpT_Secondary"), trkkpt); + histos.fill(HIST("QA/before/hy_Secondary"), trkky); + histos.fill(HIST("QA/before/hRadiusSecondary"), trkkRadius); + histos.fill(HIST("QA/before/hDCAtoPVSecondary"), trkkDCAtoPV); + histos.fill(HIST("QA/before/hCPASecondary"), trkkCPA); + histos.fill(HIST("QA/before/hInvmassSecondary"), trkkMass); } - histos.fill(HIST("QA/before/trknpionpT"), trknpt); - histos.fill(HIST("QA/before/trknpionDCAxy"), negDauTrack.dcaXY()); - histos.fill(HIST("QA/before/trknpionDCAz"), negDauTrack.dcaZ()); - - histos.fill(HIST("QA/before/hDauDCASecondary"), trkkDauDCA); - histos.fill(HIST("QA/before/hDauPosDCAtoPVSecondary"), trkkDauDCAPostoPV); - histos.fill(HIST("QA/before/hDauNegDCAtoPVSecondary"), trkkDauDCANegtoPV); - - histos.fill(HIST("QA/before/hpT_Secondary"), trkkpt); - histos.fill(HIST("QA/before/hy_Secondary"), trkky); - histos.fill(HIST("QA/before/hRadiusSecondary"), trkkRadius); - histos.fill(HIST("QA/before/hDCAtoPVSecondary"), trkkDCAtoPV); - histos.fill(HIST("QA/before/hCPASecondary"), trkkCPA); - histos.fill(HIST("QA/before/hInvmassSecondary"), trkkMass); } - if (!secondaryCutsCfgs.cfgByPassDauPIDSelection && !selectionPIDPion(posDauTrack)) // Perhaps it's already applied in trackCut (need to check QA plots) - continue; - if (!secondaryCutsCfgs.cfgByPassDauPIDSelection && !selectionPIDPion(negDauTrack)) - continue; - if (!selectionK0s(collision, K0scand)) + if ((!secondaryCutsCfgs.cfgByPassDauPIDSelection && + (!selectionPIDPion(pos) || !selectionPIDPion(neg))) || + !selectionK0s(collision, K0scand)) continue; if constexpr (!IsMix) { - // Seconddary QA plots after applying cuts - - histos.fill(HIST("QA/after/trkppionTPCPID"), trkppt, trkpNSigmaPiTPC); - if (istrkphasTOF) { - histos.fill(HIST("QA/after/trkppionTOFPID"), trkppt, trkpNSigmaPiTOF); - histos.fill(HIST("QA/after/trkppionTPCTOFPID"), trkpNSigmaPiTPC, trkpNSigmaPiTOF); - } - histos.fill(HIST("QA/after/trkppionpT"), trkppt); - histos.fill(HIST("QA/after/trkppionDCAxy"), posDauTrack.dcaXY()); - histos.fill(HIST("QA/after/trkppionDCAz"), posDauTrack.dcaZ()); - - histos.fill(HIST("QA/after/trknpionTPCPID"), trknpt, trknNSigmaPiTPC); - if (istrknhasTOF) { - histos.fill(HIST("QA/after/trknpionTOFPID"), trknpt, trknNSigmaPiTOF); - histos.fill(HIST("QA/after/trknpionTPCTOFPID"), trknNSigmaPiTPC, trknNSigmaPiTOF); + if (!doprocessMC && isQaRequired) { + histos.fill(HIST("QA/after/hDauDCASecondary"), trkkDauDCA); + histos.fill(HIST("QA/after/hDauPosDCAtoPVSecondary"), trkkDauDCAPostoPV); + histos.fill(HIST("QA/after/hDauNegDCAtoPVSecondary"), trkkDauDCANegtoPV); + histos.fill(HIST("QA/after/hpT_Secondary"), trkkpt); + histos.fill(HIST("QA/after/hy_Secondary"), trkky); + histos.fill(HIST("QA/after/hRadiusSecondary"), trkkRadius); + histos.fill(HIST("QA/after/hDCAtoPVSecondary"), trkkDCAtoPV); + histos.fill(HIST("QA/after/hCPASecondary"), trkkCPA); + histos.fill(HIST("QA/after/hInvmassSecondary"), trkkMass); } - histos.fill(HIST("QA/after/trknpionpT"), trknpt); - histos.fill(HIST("QA/after/trknpionDCAxy"), negDauTrack.dcaXY()); - histos.fill(HIST("QA/after/trknpionDCAz"), negDauTrack.dcaZ()); - - histos.fill(HIST("QA/after/hDauDCASecondary"), trkkDauDCA); - histos.fill(HIST("QA/after/hDauPosDCAtoPVSecondary"), trkkDauDCAPostoPV); - histos.fill(HIST("QA/after/hDauNegDCAtoPVSecondary"), trkkDauDCANegtoPV); - - histos.fill(HIST("QA/after/hpT_Secondary"), trkkpt); - histos.fill(HIST("QA/after/hy_Secondary"), trkky); - histos.fill(HIST("QA/after/hRadiusSecondary"), trkkRadius); - histos.fill(HIST("QA/after/hDCAtoPVSecondary"), trkkDCAtoPV); - histos.fill(HIST("QA/after/hCPASecondary"), trkkCPA); - histos.fill(HIST("QA/after/hInvmassSecondary"), trkkMass); } + k0sIndicies.push_back(K0scand.index()); } - for (const auto& trackIndex : trackIndicies) { - for (const auto& k0sIndex : k0sIndicies) { - auto bTrack = dTracks1.rawIteratorAt(trackIndex); - auto k0Scand = dTracks2.rawIteratorAt(k0sIndex); - lDecayDaughter_bach = ROOT::Math::PxPyPzMVector(bTrack.px(), bTrack.py(), bTrack.pz(), massPi); - lResoSecondary = ROOT::Math::PxPyPzMVector(k0Scand.px(), k0Scand.py(), k0Scand.pz(), massK0s); + // ========================= + // Pairing + // ========================= + for (auto tIdx : trackIndicies) { + for (auto kIdx : k0sIndicies) { + + auto bTrack = dTracks1.rawIteratorAt(tIdx); + auto k0s = dTracks2.rawIteratorAt(kIdx); + + lDecayDaughter_bach = {bTrack.px(), bTrack.py(), bTrack.pz(), massPi}; + lResoSecondary = {k0s.px(), k0s.py(), k0s.pz(), massK0s}; lResoKstar = lResoSecondary + lDecayDaughter_bach; - // QA plots if constexpr (!IsMix) { - histos.fill(HIST("QA/before/KstarRapidity"), lResoKstar.Rapidity()); - histos.fill(HIST("QA/before/kstarinvmass"), lResoKstar.M()); + if (!doprocessMC && isQaRequired) { + histos.fill(HIST("QA/before/KstarRapidity"), lResoKstar.Rapidity()); + histos.fill(HIST("QA/before/kstarinvmass"), lResoKstar.M()); + } } - if (lResoKstar.Rapidity() > kstarCutCfgs.cKstarMaxRap || lResoKstar.Rapidity() < kstarCutCfgs.cKstarMinRap) + if (lResoKstar.Rapidity() > kstarCutCfgs.cKstarMaxRap || + lResoKstar.Rapidity() < kstarCutCfgs.cKstarMinRap) continue; if constexpr (!IsMix) { + if (!doprocessMC && isQaRequired) { + histos.fill(HIST("QA/after/KstarRapidity"), lResoKstar.Rapidity()); + histos.fill(HIST("QA/after/kstarinvmass"), lResoKstar.M()); + } - histos.fill(HIST("QA/after/KstarRapidity"), lResoKstar.Rapidity()); - histos.fill(HIST("QA/after/kstarinvmass"), lResoKstar.M()); - histos.fill(HIST("hInvmass_Kstar"), collision.centFT0M(), lResoKstar.Pt(), lResoKstar.M()); if (helicityCfgs.cBoostKShot) { fillInvMass(lResoKstar, collision.centFT0M(), lResoSecondary, lDecayDaughter_bach, IsMix); } else { fillInvMass(lResoKstar, collision.centFT0M(), lDecayDaughter_bach, lResoSecondary, IsMix); } } else { - - histos.fill(HIST("hInvmass_KstarME"), collision.centFT0M(), lResoKstar.Pt(), lResoKstar.M()); fillInvMass(lResoKstar, collision.centFT0M(), lResoSecondary, lDecayDaughter_bach, IsMix); } - if constexpr (!IsMix) { - if (rotBkgEstCfgs.cFillRotBkg) { - for (int nrotbkg = 0; nrotbkg < rotBkgEstCfgs.nBkgRotations; nrotbkg++) { - auto rotangle = o2::constants::math::PI; // If there is only one rotation then it should be pi ): - if (rotBkgEstCfgs.nBkgRotations > 1) { - auto anglestart = rotBkgEstCfgs.confMinRot; - auto angleend = rotBkgEstCfgs.confMaxRot; - auto anglestep = (angleend - anglestart) / (1.0 * (rotBkgEstCfgs.nBkgRotations - 1)); - rotangle = anglestart + nrotbkg * anglestep; - } - histos.fill(HIST("hRotation"), rotangle); - auto rotpionPx = lDecayDaughter_bach.Px() * std::cos(rotangle) - lDecayDaughter_bach.Py() * std::sin(rotangle); - auto rotpionPy = lDecayDaughter_bach.Px() * std::sin(rotangle) + lDecayDaughter_bach.Py() * std::cos(rotangle); - ROOT::Math::PtEtaPhiMVector pionrot; - pionrot = ROOT::Math::PxPyPzMVector(rotpionPx, rotpionPy, lDecayDaughter_bach.Pz(), massPi); - chargeKstarrot = pionrot + lResoSecondary; - if (chargeKstarrot.Rapidity() > kstarCutCfgs.cKstarMaxRap || chargeKstarrot.Rapidity() < kstarCutCfgs.cKstarMinRap) - continue; - histos.fill(HIST("hInvmass_KstarRotated"), collision.centFT0M(), chargeKstarrot.Pt(), chargeKstarrot.M()); - } - } - } - } // K0scand - } // bTrack + } + } count++; - - } // fillHistograms + } // process data void processDataSE(EventCandidates::iterator const& collision, @@ -1250,7 +1188,6 @@ struct Chargedkstaranalysis { } PROCESS_SWITCH(Chargedkstaranalysis, processDataSE, "Process Event for data without Partitioning", true); - SliceCache cache; using BinningTypeVertexContributor = ColumnBinningPolicy; BinningTypeVertexContributor binningOnPositions{{axisCfgs.cfgvtxbins, axisCfgs.cfgmultbins}, true}; Pair pair{binningOnPositions, nEvtMixing, -1, &cache}; @@ -1296,15 +1233,178 @@ struct Chargedkstaranalysis { } PROCESS_SWITCH(Chargedkstaranalysis, processDataME, "Process Event for data without Partitioning", true); - // process MC reconstructed level - void processMC(MCEventCandidates::iterator const& collision, - MCTrackCandidates const& tracks, - MCV0Candidates const& v0s) + void processMC(soa::Join const&, aod::McParticles& mcParticles, soa::Join const& events, MCV0Candidates const& v0s, MCTrackCandidates const& tracks) { + allowedMcIds.clear(); + centTruthByAllowed.clear(); + + // To apply event selection and store the collision IDs of reconstructed tracks that pass the selection criteria + for (const auto& coll : events) { - // histos.fill(HIST("QAMC/hEvent"), 1.0); + if (!coll.has_mcCollision()) + continue; + + const auto mcid = coll.mcCollisionId(); + + const auto mccoll = coll.template mcCollision_as>(); + const float lCentrality = mccoll.centFT0M(); + + if (doprocessMC && isQaRequired) { + histos.fill(HIST("QA/MC/QACent_woCut"), lCentrality); + histos.fill(HIST("QA/MC/QAvtxz_woCut"), coll.posZ()); + } + + if (!colCuts.isSelected(coll)) + continue; + if (rctCut.requireRCTFlagChecker && !rctCut.rctChecker(coll)) + continue; + if (!coll.isInelGt0()) + continue; + + if (doprocessMC && isQaRequired) { + histos.fill(HIST("QA/MC/QACent_woCentCut"), lCentrality); + histos.fill(HIST("QA/MC/QAvtxz_wVtxzCut"), coll.posZ()); + } + + if (lCentrality < eventCutCfgs.cfgEventCentralityMin || lCentrality > eventCutCfgs.cfgEventCentralityMax) + continue; - fillHistograms(collision, tracks, v0s); + if (doprocessMC && isQaRequired) { + histos.fill(HIST("QA/MC/QACent_wCentCut"), lCentrality); + } + allowedMcIds.insert(mcid); + centTruthByAllowed.emplace(mcid, lCentrality); + } + + // Calculating the generated Kstar + for (const auto& part : mcParticles) { + if (!part.has_mcCollision()) + continue; + if (std::abs(part.pdgCode()) != kKstarPlus) + continue; + if (std::abs(part.y()) > kstarCutCfgs.cKstarMaxRap) + continue; + + const int pionWanted = (part.pdgCode() > 0) ? +kPiPlus : -kPiPlus; + bool hasRightPion = false; + bool hasK0sToPipi = false; + + for (const auto& d1 : part.template daughters_as()) { + const int pdg1 = d1.pdgCode(); + if (pdg1 == pionWanted) { + hasRightPion = true; + } else if (std::abs(pdg1) == kPDGK0) { + for (const auto& d2 : d1.template daughters_as()) { + if (std::abs(d2.pdgCode()) == kPDGK0s) { + bool seenPip = false, seenPim = false; + for (const auto& d3 : d2.template daughters_as()) { + if (d3.pdgCode() == +kPiPlus) + seenPip = true; + else if (d3.pdgCode() == -kPiPlus) + seenPim = true; + } + if (seenPip && seenPim) { + hasK0sToPipi = true; + break; + } + } + } + } + if (hasRightPion && hasK0sToPipi) + break; + } + + if (!(hasRightPion && hasK0sToPipi)) + continue; + + const auto mcid = part.mcCollisionId(); + if (allowedMcIds.count(mcid) == 0) + continue; + + auto iter = centTruthByAllowed.find(mcid); + if (iter == centTruthByAllowed.end()) + continue; + + const float lCentrality = iter->second; + + histos.fill(HIST("EffKstar/genKstar"), part.pt(), lCentrality); + + if (part.vt() == 0) { + histos.fill(HIST("EffKstar/genKstar_pri"), part.pt(), lCentrality); // To check the primary particle + } + } + // To store the recoKstar + for (const auto& v0 : v0s) { + auto coll = v0.template collision_as(); + + if (!coll.has_mcCollision()) + continue; + + const auto mcid = coll.mcCollisionId(); + + if (allowedMcIds.count(mcid) == 0) + continue; // To check the event is allowed or not + + const auto mccoll = coll.template mcCollision_as>(); + const float lCentrality = mccoll.centFT0M(); + + if (!secondaryCutsCfgs.cfgByPassDauPIDSelection) { + auto posDauTrack = v0.template posTrack_as(); + auto negDauTrack = v0.template negTrack_as(); + if (!selectionPIDPion(posDauTrack)) + continue; + if (!selectionPIDPion(negDauTrack)) + continue; + } + if (!selectionK0s(coll, v0)) + continue; + + auto trks = tracks.sliceBy(perCollision, v0.collisionId()); // Grouping the tracks with the v0s, means only those tracks that belong to the same collision as v0 + for (const auto& bTrack : trks) { + if (bTrack.collisionId() != v0.collisionId()) + continue; + if (!trackCut(bTrack)) + continue; + if (!selectionPIDPion(bTrack)) + continue; + + LorentzVectorSetXYZM lResoSecondary, lDecayDaughter_bach, lResoKstar, lDaughterRot; + + lResoSecondary = LorentzVectorSetXYZM(v0.px(), v0.py(), v0.pz(), MassK0Short); + lDecayDaughter_bach = LorentzVectorSetXYZM(bTrack.px(), bTrack.py(), bTrack.pz(), MassPionCharged); + lResoKstar = lResoSecondary + lDecayDaughter_bach; + + const double ptreco = lResoKstar.Pt(); + const double yreco = lResoKstar.Rapidity(); + if (std::abs(yreco) > kstarCutCfgs.cKstarMaxRap) + continue; + + // Since we are doing the MC study and we know about the PDG code of each particle let's try to check the things which we have + if (!v0.has_mcParticle() || !bTrack.has_mcParticle()) + continue; + auto mcK0s = v0.template mcParticle_as(); // To get the MC truth particle corressponds to the V0 candidate + auto mcPi = bTrack.template mcParticle_as(); + if (std::abs(mcK0s.pdgCode()) != kPDGK0s) + continue; + if (std::abs(mcPi.pdgCode()) != kPiPlus) + continue; + MCTrueTrackCandidates::iterator kstarFromPi; + bool havePiKstar = false; + // Loops over all the mother's of pions and check if this pion comming from a kstar + for (const auto& m1 : mcPi.template mothers_as()) { + if (std::abs(m1.pdgCode()) == kKstarPlus) { + kstarFromPi = m1; + havePiKstar = true; + break; + } + } + if (!havePiKstar) { + continue; + } + histos.fill(HIST("EffKstar/recoKstar"), ptreco, lCentrality); + histos.fill(HIST("MCReco/hInvmass_Kstar_true"), lCentrality, ptreco, lResoKstar.M()); + } + } } PROCESS_SWITCH(Chargedkstaranalysis, processMC, "Process Event for MC", false); };