From 9f65830cc15f8ecb324ba790f4087185ba58dab6 Mon Sep 17 00:00:00 2001 From: Subhadeep Mandal Date: Sat, 22 Nov 2025 17:25:40 +0530 Subject: [PATCH] Added separate task file for Kstar analysis in Light ion collisions --- PWGLF/Tasks/Resonances/CMakeLists.txt | 5 + PWGLF/Tasks/Resonances/kstar892LightIon.cxx | 1539 +++++++++++++++++++ 2 files changed, 1544 insertions(+) create mode 100644 PWGLF/Tasks/Resonances/kstar892LightIon.cxx diff --git a/PWGLF/Tasks/Resonances/CMakeLists.txt b/PWGLF/Tasks/Resonances/CMakeLists.txt index f466c75e7e3..660eb015a1e 100644 --- a/PWGLF/Tasks/Resonances/CMakeLists.txt +++ b/PWGLF/Tasks/Resonances/CMakeLists.txt @@ -259,6 +259,11 @@ o2physics_add_dpl_workflow(phioo PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(kstar892-light-ion + SOURCES kstar892LightIon.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(phispectrapbpbqa SOURCES phispectrapbpbqa.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore diff --git a/PWGLF/Tasks/Resonances/kstar892LightIon.cxx b/PWGLF/Tasks/Resonances/kstar892LightIon.cxx new file mode 100644 index 00000000000..8c4da69d91f --- /dev/null +++ b/PWGLF/Tasks/Resonances/kstar892LightIon.cxx @@ -0,0 +1,1539 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file kstar892LightIon.cxx +/// \brief Code for K*0(892) resonance without resonance initializer in Light Ion collisions +/// \author Subhadeep Mandal +/// \since 22/11/2025 + +#include "PWGLF/DataModel/mcCentrality.h" + +#include "Common/Core/TrackSelection.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponseTOF.h" +#include "Common/DataModel/PIDResponseTPC.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "CCDB/BasicCCDBManager.h" +#include "CCDB/CcdbApi.h" +#include "CommonConstants/PhysicsConstants.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/StepTHn.h" +#include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/Track.h" + +#include "Math/GenVector/Boost.h" +#include "Math/Vector3D.h" +#include "Math/Vector4D.h" +#include "TRandom3.h" +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::soa; +using std::array; +using namespace o2::aod::rctsel; + +struct Kstar892LightIon { + SliceCache cache; + + struct : ConfigurableGroup { + Configurable requireRCTFlagChecker{"requireRCTFlagChecker", true, "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"}; + } rctCut; + RCTFlagsChecker rctChecker; + + struct : ConfigurableGroup { + // Configurables for event selections + Configurable cutVrtxZ{"cutVrtxZ", 10.0f, "Accepted z-vertex range (cm)"}; + Configurable isTriggerTVX{"isTriggerTVX", true, "TriggerTVX"}; + Configurable isGoodZvtxFT0vsPV{"isGoodZvtxFT0vsPV", true, "IsGoodZvtxFT0vsPV"}; + Configurable isApplyOccCut{"isApplyOccCut", false, "Apply occupancy cut"}; + Configurable cfgOccCut{"cfgOccCut", 1000., "Occupancy cut"}; + Configurable isNoSameBunchPileup{"isNoSameBunchPileup", true, "kNoSameBunchPileup"}; + Configurable isGoodITSLayersAll{"isGoodITSLayersAll", false, "Require all ITS layers to be good"}; + Configurable isNoTimeFrameBorder{"isNoTimeFrameBorder", true, "kNoTimeFrameBorder"}; + Configurable isNoITSROFrameBorder{"isNoITSROFrameBorder", true, "kNoITSROFrameBorder"}; + Configurable isApplyDeepAngle{"isApplyDeepAngle", false, "Deep Angle cut"}; + Configurable isNoCollInTimeRangeStandard{"isNoCollInTimeRangeStandard", false, "No collision in time range standard"}; + Configurable isVertexITSTPC{"isVertexITSTPC", false, "Vertex ITS TPC"}; + Configurable isVertexTOFMatched{"isVertexTOFMatched", false, "Vertex TOF Matched"}; + + // Configurables for track selections + Configurable cfgPVContributor{"cfgPVContributor", false, "PV contributor track selection"}; // PV Contriuibutor + Configurable cfgPrimaryTrack{"cfgPrimaryTrack", false, "Primary track selection"}; // kGoldenChi2 | kDCAxy | kDCAz + Configurable isGlobalTracks{"isGlobalTracks", true, "isGlobalTracks"}; + + Configurable cfgCutPT{"cfgCutPT", 0.1f, "PT cut on daughter track"}; + Configurable cfgCutEtaMax{"cfgCutEtaMax", 0.8f, "Eta cut on daughter track"}; + Configurable cfgCutDCAxyMax{"cfgCutDCAxyMax", 0.1f, "DCAxy range for tracks"}; + Configurable cfgCutDCAz{"cfgCutDCAz", 0.1f, "DCAz range for tracks"}; + Configurable cfgNoMixedEvents{"cfgNoMixedEvents", 15, "Number of mixed events per event"}; + Configurable cfgITScluster{"cfgITScluster", 0, "Number of ITS cluster"}; + Configurable cfgTPCcluster{"cfgTPCcluster", 70, "Number of TPC cluster"}; + Configurable cfgRCRFC{"cfgRCRFC", 0.8f, "Crossed Rows to Findable Clusters"}; + Configurable cfgITSChi2NCl{"cfgITSChi2NCl", 36.0, "ITS Chi2/NCl"}; + Configurable cfgTPCChi2NClMax{"cfgTPCChi2NClMax", 4.0, "TPC Chi2/NCl"}; + Configurable cfgTPCChi2NClMin{"cfgTPCChi2NClMin", 0.0, "TPC Chi2/NCl"}; + Configurable cfgUseITSTPCRefit{"cfgUseITSTPCRefit", false, "Require ITS Refit"}; + Configurable isApplyPtDepDCAxyCut{"isApplyPtDepDCAxyCut", false, "Apply pT dependent DCAxy cut"}; + Configurable isGoldenChi2{"isGoldenChi2", false, "Apply golden chi2 cut"}; + Configurable cfgDeepAngle{"cfgDeepAngle", 0.04, "Deep Angle cut value"}; + + // cuts on mother + // Configurable isApplyCutsOnMother{"isApplyCutsOnMother", false, "Enable additional cuts on Kstar mother"}; + // Configurable cMaxPtMotherCut{"cMaxPtMotherCut", 15.0, "Maximum pt of mother cut"}; + // Configurable cMaxMinvMotherCut{"cMaxMinvMotherCut", 1.5, "Maximum mass of mother cut"}; + Configurable motherRapidityCut{"motherRapidityCut", 0.5, "Maximum rapidity of mother"}; + // PID selections + Configurable nsigmaCutTPCPi{"nsigmaCutTPCPi", 3.0, "TPC Nsigma cut for pions"}; + Configurable nsigmaCutTPCKa{"nsigmaCutTPCKa", 3.0, "TPC Nsigma cut for kaons"}; + Configurable nsigmaCutTOFPi{"nsigmaCutTOFPi", 3.0, "TOF Nsigma cut for pions"}; + Configurable nsigmaCutTOFKa{"nsigmaCutTOFKa", 3.0, "TOF Nsigma cut for kaons"}; + Configurable nsigmaCutCombinedKa{"nsigmaCutCombinedKa", 3.0, "Combined Nsigma cut for kaon"}; + Configurable nsigmaCutCombinedPi{"nsigmaCutCombinedPi", 3.0, "Combined Nsigma cut for pion"}; + + // Fixed variables + float lowPtCutPID = 0.5; + } selectionConfig; + + // Histograms are defined with HistogramRegistry + HistogramRegistry hEventSelection{"eventSelection", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + HistogramRegistry hInvMass{"InvMass", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + HistogramRegistry hMC{"MC", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + HistogramRegistry hPID{"PID", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + HistogramRegistry hOthers{"Others", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + + Configurable calcLikeSign{"calcLikeSign", true, "Calculate Like Sign"}; + Configurable calcRotational{"calcRotational", true, "Calculate Rotational"}; + Configurable cRotations{"cRotations", 3, "Number of random rotations in the rotational background"}; + Configurable rotationalCut{"rotationalCut", 10, "Cut value (Rotation angle pi - pi/cut and pi + pi/cut)"}; + + // Confugrable for QA histograms + Configurable cQAplots{"cQAplots", true, "cQAplots"}; + Configurable cQAevents{"cQAevents", true, "centrality dist, DCAxy, DCAz"}; + Configurable onlyTOF{"onlyTOF", false, "only TOF tracks"}; + Configurable onlyTOFHIT{"onlyTOFHIT", false, "accept only TOF hit tracks at high pt"}; + Configurable onlyTPC{"onlyTPC", false, "only TPC tracks"}; + Configurable cSelectMultEstimator{"cSelectMultEstimator", 0, "Select centrality estimator: 0 - FT0M, 1 - FT0A, 2 - FT0C, 3 - FV0A"}; + Configurable applypTdepPID{"applypTdepPID", false, "Apply pT dependent PID"}; + + // Configurable cfgGlobalWoDCATrack{"cfgGlobalWoDCATrack", false, "Global track selection without DCA"}; // kQualityTracks (kTrackType | kTPCNCls | kTPCCrossedRows | kTPCCrossedRowsOverNCls | kTPCChi2NDF | kTPCRefit | kITSNCls | kITSChi2NDF | kITSRefit | kITSHits) | kInAcceptanceTracks (kPtRange | kEtaRange) + Configurable cBetaCutTOF{"cBetaCutTOF", 0.0, "cut TOF beta"}; + + // Configurable for histograms + Configurable avoidsplitrackMC{"avoidsplitrackMC", true, "avoid split track in MC"}; + ConfigurableAxis binsCentPlot{"binsCentPlot", {110, 0.0, 110}, "THnSpare centrality axis"}; + ConfigurableAxis axisdEdx{"axisdEdx", {1, 0.0f, 200.0f}, "dE/dx (a.u.)"}; + ConfigurableAxis axisPtfordEbydx{"axisPtfordEbydx", {1, 0, 20}, "pT (GeV/c)"}; + // ConfigurableAxis axisMultdist{"axisMultdist", {1, 0, 70000}, "centrality distribution"}; + ConfigurableAxis invMassKstarAxis{"invMassKstarAxis", {300, 0.7f, 1.3f}, "Kstar invariant mass axis"}; + ConfigurableAxis ptAxisKstar{"ptAxisKstar", {200, 0.0f, 20.0f}, "Kstar pT axis"}; + ConfigurableAxis binsImpactPar{"binsImpactPar", {100, 0, 25}, "Binning of the impact parameter axis"}; + ConfigurableAxis axisNch{"axisNch", {100, 0.0f, 100.0f}, "Number of charged particles in |y| < 0.5"}; + + Configurable indexCheck{"indexCheck", true, "Check if track2id < track1id matters"}; // check and remove + + enum MultEstimator { + kFT0M, + kFT0A, + kFT0C, + kFV0A, + kFV0C, + kFV0M, + kNEstimators // useful if you want to iterate or size things + }; + + enum PIDParticle { + kPion, + kKaon, + kProton + }; + + int noOfDaughters = 2; + + TRandom* rn = new TRandom(); + + void init(InitContext const&) + { + rctChecker.init(rctCut.cfgEvtRCTFlagCheckerLabel, rctCut.cfgEvtRCTFlagCheckerZDCCheck, rctCut.cfgEvtRCTFlagCheckerLimitAcceptAsBad); + // Axes + AxisSpec vertexZAxis = {60, -15., 15., "vrtx_{Z} [cm] for plots"}; + AxisSpec ptAxis = {ptAxisKstar, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec invmassAxis = {invMassKstarAxis, "Invariant mass (GeV/#it{c}^{2})"}; + AxisSpec centralityAxis = {binsCentPlot, "centrality Axis"}; + AxisSpec impactParAxis = {binsImpactPar, "Impact Parameter (cm)"}; + + // Histograms + // Event selection + hEventSelection.add("hVertexZ", "hVertexZ", kTH1F, {vertexZAxis}); + hEventSelection.add("hCentrality", "Centrality percentile", kTH1F, {{110, 0, 110}}); + + hEventSelection.add("hEventCut", "No. of event after cuts", kTH1D, {{20, 0, 20}}); + std::shared_ptr hCutFlow = hEventSelection.get(HIST("hEventCut")); + + auto check = [](bool enabled) { return enabled ? "" : " #otimes"; }; // check if a cut is enabled and put #otimes beside that label if not enabled + + std::vector eveCutLabels = { + "All Events", + Form("|Vz| < %.1f", selectionConfig.cutVrtxZ.value), + "sel8", + std::string("kNoTimeFrameBorder") + check(selectionConfig.isNoTimeFrameBorder.value), + std::string("kNoITSROFrameBorder") + check(selectionConfig.isNoITSROFrameBorder.value), + std::string("kNoSameBunchPileup") + check(selectionConfig.isNoSameBunchPileup.value), + std::string("kIsGoodITSLayersAll") + check(selectionConfig.isGoodITSLayersAll.value), + std::string("kNoCollInTimeRangeStandard") + check(selectionConfig.isNoCollInTimeRangeStandard.value), + Form("Occupancy < %.0f%s", selectionConfig.cfgOccCut.value, check(selectionConfig.isApplyOccCut.value)), + std::string("rctChecker") + check(rctCut.requireRCTFlagChecker.value), + std::string("kIsTriggerTVX") + check(selectionConfig.isTriggerTVX.value), + std::string("kIsGoodZvtxFT0vsPV") + check(selectionConfig.isGoodZvtxFT0vsPV.value), + std::string("isVertexITSTPC") + check(selectionConfig.isVertexITSTPC.value), + std::string("isVertexTOFMatched") + check(selectionConfig.isVertexTOFMatched.value)}; + // assign labels + for (size_t i = 0; i < eveCutLabels.size(); ++i) { + hCutFlow->GetXaxis()->SetBinLabel(i + 1, eveCutLabels[i].c_str()); + } + + // for primary tracksbinsCentPlot + if (cQAplots) { + hOthers.add("dE_by_dx_TPC", "dE/dx signal in the TPC as a function of pT", kTH2F, {axisPtfordEbydx, axisdEdx}); + hOthers.add("hEta_after", "Eta distribution", kTH1F, {{200, -1.0f, 1.0f}}); + hOthers.add("hCRFC_after", "CRFC after distribution", kTH1F, {{100, 0.0f, 10.0f}}); + hOthers.add("hCRFC_before", "CRFC before distribution", kTH1F, {{100, 0.0f, 10.0f}}); + + hOthers.add("hKstar_rap_pt", "Pair rapidity distribution; y; p_{T}; Counts", kTH2F, {{400, -2.0f, 2.0f}, ptAxis}); + hOthers.add("hKstar_eta_pt", "Pair eta distribution; #eta; p_{T}; Counts", kTH2F, {{400, -2.0f, 2.0f}, ptAxis}); + + hOthers.add("hDcaxyPi", "Dcaxy distribution of selected Pions", kTH1F, {{200, -1.0f, 1.0f}}); + hOthers.add("hDcaxyKa", "Dcaxy distribution of selected Kaons", kTH1F, {{200, -1.0f, 1.0f}}); + hOthers.add("hDcazPi", "Dcaz distribution of selected Pions", kTH1F, {{200, -1.0f, 1.0f}}); + hOthers.add("hDcazKa", "Dcaz distribution of selected Kaons", kTH1F, {{200, -1.0f, 1.0f}}); + + hOthers.add("hDcaxy_cent_pt", "Dcaxy distribution before PID", kTH3F, {{200, -1.0f, 1.0f}, centralityAxis, ptAxis}); + hOthers.add("hDcaz_cent_pt", "Dcaz distribution before PID", kTH3F, {{200, -1.0f, 1.0f}, centralityAxis, ptAxis}); + + hPID.add("Before/hNsigma_TPC_TOF_Ka_pt", "N #sigma Kaon TPC TOF before", kTH3F, {{50, -5.0f, 5.0f}, {50, -5.0f, 5.0f}, ptAxis}); + hPID.add("Before/hNsigma_TPC_TOF_Pi_pt", "N #sigma Pion TPC TOF before", kTH3F, {{50, -5.0f, 5.0f}, {50, -5.0f, 5.0f}, ptAxis}); + + hPID.add("Before/hTPCnsigKa_Neg_mult_pt", "TPC nsigma of K^{-} before PID with pt and centrality", kTH3F, {{100, -10.0f, 10.0f}, centralityAxis, ptAxis}); + hPID.add("Before/hTPCnsigPi_Neg_mult_pt", "TPC nsigma of #pi^{-} before PID with pt and centrality", kTH3F, {{100, -10.0f, 10.0f}, centralityAxis, ptAxis}); + hPID.add("Before/hTOFnsigKa_Neg_mult_pt", "TOF nsigma of K^{-} before PID with pt and centrality", kTH3F, {{100, -10.0f, 10.0f}, centralityAxis, ptAxis}); + hPID.add("Before/hTOFnsigPi_Neg_mult_pt", "TOF nsigma of #pi^{-} before PID with pt and centrality", kTH3F, {{100, -10.0f, 10.0f}, centralityAxis, ptAxis}); + + hPID.add("Before/hTPCnsigKa_Pos_mult_pt", "TPC nsigma of K^{+} before PID with pt and centrality", kTH3F, {{100, -10.0f, 10.0f}, centralityAxis, ptAxis}); + hPID.add("Before/hTPCnsigPi_Pos_mult_pt", "TPC nsigma of #pi^{+} before PID with pt and centrality", kTH3F, {{100, -10.0f, 10.0f}, centralityAxis, ptAxis}); + hPID.add("Before/hTOFnsigKa_Pos_mult_pt", "TOF nsigma of K^{+} before PID with pt and centrality", kTH3F, {{100, -10.0f, 10.0f}, centralityAxis, ptAxis}); + hPID.add("Before/hTOFnsigPi_Pos_mult_pt", "TOF nsigma of #pi^{+} before PID with pt and centrality", kTH3F, {{100, -10.0f, 10.0f}, centralityAxis, ptAxis}); + + hPID.add("After/hNsigma_TPC_TOF_Ka_pt", "N #sigma Kaon TPC TOF after", kTH3F, {{50, -5.0f, 5.0f}, {50, -5.0f, 5.0f}, ptAxis}); + hPID.add("After/hNsigma_TPC_TOF_Pi_pt", "N #sigma Pion TPC TOF after", kTH3F, {{50, -5.0f, 5.0f}, {50, -5.0f, 5.0f}, ptAxis}); + + hPID.add("After/hTPCnsigKa_Neg_mult_pt", "TPC nsigma of K^{-} after PID with pt and centrality", kTH3F, {{100, -10.0f, 10.0f}, centralityAxis, ptAxis}); + hPID.add("After/hTPCnsigPi_Neg_mult_pt", "TPC nsigma of #pi^{-} after PID with pt and centrality", kTH3F, {{100, -10.0f, 10.0f}, centralityAxis, ptAxis}); + hPID.add("After/hTOFnsigKa_Neg_mult_pt", "TOF nsigma of K^{-} after PID with pt and centrality", kTH3F, {{100, -10.0f, 10.0f}, centralityAxis, ptAxis}); + hPID.add("After/hTOFnsigPi_Neg_mult_pt", "TOF nsigma of #pi^{-} after PID with pt and centrality", kTH3F, {{100, -10.0f, 10.0f}, centralityAxis, ptAxis}); + + hPID.add("After/hTPCnsigKa_Pos_mult_pt", "TPC nsigma of K^{+} after PID with pt and centrality", kTH3F, {{100, -10.0f, 10.0f}, centralityAxis, ptAxis}); + hPID.add("After/hTPCnsigPi_Pos_mult_pt", "TPC nsigma of #pi^{+} after PID with pt and centrality", kTH3F, {{100, -10.0f, 10.0f}, centralityAxis, ptAxis}); + hPID.add("After/hTOFnsigKa_Pos_mult_pt", "TOF nsigma of K^{+} after PID with pt and centrality", kTH3F, {{100, -10.0f, 10.0f}, centralityAxis, ptAxis}); + hPID.add("After/hTOFnsigPi_Pos_mult_pt", "TOF nsigma of #pi^{+} after PID with pt and centrality", kTH3F, {{100, -10.0f, 10.0f}, centralityAxis, ptAxis}); + } + + // KStar histograms + hInvMass.add("h3KstarInvMassUnlikeSign", "kstar Unlike Sign", kTHnSparseF, {centralityAxis, ptAxis, invmassAxis}); + hInvMass.add("h3KstarInvMassMixed", "kstar Mixed", kTHnSparseF, {centralityAxis, ptAxis, invmassAxis}); + if (calcLikeSign) { + hInvMass.add("h3KstarInvMasslikeSignPP", "kstar like Sign", kTHnSparseF, {centralityAxis, ptAxis, invmassAxis}); + hInvMass.add("h3KstarInvMasslikeSignMM", "kstar like Sign", kTHnSparseF, {centralityAxis, ptAxis, invmassAxis}); + } + if (calcRotational) + hInvMass.add("h3KstarInvMassRotated", "kstar rotated", kTHnSparseF, {centralityAxis, ptAxis, invmassAxis}); + + // MC histograms + if (doprocessGen) { + hMC.add("hk892GenpT", "pT distribution of True MC K(892)0", kTHnSparseF, {ptAxis, centralityAxis}); + hMC.add("hk892GenpT2", "pT distribution of True MC K(892)0", kTHnSparseF, {ptAxis, centralityAxis}); + hMC.add("h1genmass", "Invariant mass of generated kstar meson", kTH1F, {invmassAxis}); + hMC.add("h1GenCent", "centrality generated", kTH1F, {centralityAxis}); + hMC.add("hAllGenCollisions", "All generated events", kTH1F, {centralityAxis}); + hMC.add("hAllGenCollisions1Rec", "All gen events with at least one rec event", kTH1F, {centralityAxis}); + hMC.add("hAllKstarGenCollisisons", "All generated Kstar in events with rapidity in 0.5", kTH2F, {ptAxis, centralityAxis}); + hMC.add("hAllKstarGenCollisisons1Rec", "All generated Kstar in events with at least one rec event in rapidity in 0.5", kTH2F, {ptAxis, centralityAxis}); + } + + if (doprocessRec) { + hMC.add("hAllRecCollisions", "All reconstructed events", kTH1F, {centralityAxis}); + hMC.add("h1KstarRecMass", "Invariant mass of kstar meson", kTH1F, {invmassAxis}); + hMC.add("h2KstarRecpt1", "pT of kstar meson", kTHnSparseF, {ptAxis, centralityAxis, invmassAxis}); + hMC.add("h2KstarRecpt2", "pT of kstar meson", kTHnSparseF, {ptAxis, centralityAxis, invmassAxis}); + hMC.add("h1RecCent", "centrality reconstructed", kTH1F, {centralityAxis}); + hMC.add("h1KSRecsplit", "KS meson Rec split", kTH1F, {{100, 0.0f, 10.0f}}); + } + + // Signal Loss & Event Loss + if (doprocessEvtLossSigLossMC) { + hMC.add("ImpactCorr/hImpactParameterGen", "Impact parameter of generated MC events", kTH1F, {impactParAxis}); + hMC.add("ImpactCorr/hImpactParameterRec", "Impact parameter of selected MC events", kTH1F, {impactParAxis}); + hMC.add("ImpactCorr/hImpactParvsCentrRec", "Impact parameter of selected MC events vs centrality", kTH2F, {{centralityAxis}, impactParAxis}); + hMC.add("ImpactCorr/hKstarGenBeforeEvtSel", "K*0 before event selections", kTH2F, {ptAxis, impactParAxis}); + hMC.add("ImpactCorr/hKstarGenAfterEvtSel", "K*0 after event selections", kTH2F, {ptAxis, impactParAxis}); + } + + if (doprocessCorrFactors) { + hMC.add("CorrFactors/hCentralityVsMultMC", "Event centrality vs MC centrality", kTH2F, {{101, 0.0f, 101.0f}, axisNch}); + hMC.add("CorrFactors/hEventCentrality", "Event centrality", kTH1F, {{101, 0, 101}}); + hMC.add("CorrFactors/hNrecInGen", "Number of collisions in MC", kTH1F, {{4, -0.5, 3.5}}); + hMC.add("CorrFactors/hGenEvents", "Generated events", kTH2F, {{axisNch}, {4, 0, 4}}); + auto hGenEvents = hMC.get(HIST("CorrFactors/hGenEvents")); + hGenEvents->GetYaxis()->SetBinLabel(1, "All generated events"); + hGenEvents->GetYaxis()->SetBinLabel(2, "Generated events with Mc collision V_{z} cut"); + hGenEvents->GetYaxis()->SetBinLabel(3, "Generated events with at least one reconstructed event"); + hMC.add("CorrFactors/h2dGenKstar", "Centrality vs p_{T}", kTH2D, {{101, 0.0f, 101.0f}, ptAxis}); + hMC.add("CorrFactors/h3dGenKstarVsMultMCVsCentrality", "MC centrality vs centrality vs p_{T}", kTH3D, {axisNch, {101, 0.0f, 101.0f}, ptAxis}); + } + + hEventSelection.add("tracksCheckData", "No. of events in the data", kTH1I, {{10, 0, 10}}); + hEventSelection.add("eventsCheckGen", "No. of events in the generated MC", kTH1I, {{10, 0, 10}}); + hEventSelection.add("recMCparticles", "No. of events in the reconstructed MC", kTH1I, {{20, 0, 20}}); + hEventSelection.add("hOccupancy", "Occupancy distribution", kTH1F, {{1000, 0, 15000}}); + + std::shared_ptr hrecLabel = hEventSelection.get(HIST("recMCparticles")); + hrecLabel->GetXaxis()->SetBinLabel(1, "All tracks"); + hrecLabel->GetXaxis()->SetBinLabel(2, "has_MC"); + hrecLabel->GetXaxis()->SetBinLabel(3, "Track selection"); + hrecLabel->GetXaxis()->SetBinLabel(4, "StrictlyUpperIndex"); + hrecLabel->GetXaxis()->SetBinLabel(5, "Unlike Sign"); + hrecLabel->GetXaxis()->SetBinLabel(6, "Physical Primary"); + hrecLabel->GetXaxis()->SetBinLabel(7, "TrackPDG Check1"); + hrecLabel->GetXaxis()->SetBinLabel(8, "TrackPDG Check2"); + hrecLabel->GetXaxis()->SetBinLabel(9, "Global Index"); + hrecLabel->GetXaxis()->SetBinLabel(10, "Generator"); + hrecLabel->GetXaxis()->SetBinLabel(11, "Mother y"); + hrecLabel->GetXaxis()->SetBinLabel(12, "Mother PDG"); + hrecLabel->GetXaxis()->SetBinLabel(13, "Track PID"); + hrecLabel->GetXaxis()->SetBinLabel(14, "Track MID"); + hrecLabel->GetXaxis()->SetBinLabel(15, "Track y"); + hrecLabel->GetXaxis()->SetBinLabel(16, "Split tracks"); + hrecLabel->GetXaxis()->SetBinLabel(17, "DeepAngle Cut"); + + std::shared_ptr hDataTracks = hEventSelection.get(HIST("tracksCheckData")); + hDataTracks->GetXaxis()->SetBinLabel(1, "All tracks"); + hDataTracks->GetXaxis()->SetBinLabel(2, "Track selection"); + hDataTracks->GetXaxis()->SetBinLabel(3, "PID Cut"); + hDataTracks->GetXaxis()->SetBinLabel(4, "Remove Fake Tracks"); + hDataTracks->GetXaxis()->SetBinLabel(5, "Rapidity Cut"); + hDataTracks->GetXaxis()->SetBinLabel(6, "MID"); + hDataTracks->GetXaxis()->SetBinLabel(7, "DeepAngle Cut"); + + std::shared_ptr hGenTracks = hEventSelection.get(HIST("eventsCheckGen")); + hGenTracks->GetXaxis()->SetBinLabel(1, "All events"); + hGenTracks->GetXaxis()->SetBinLabel(4, "Event Reconstructed"); + } + + double massPi = o2::constants::physics::MassPiPlus; + double massKa = o2::constants::physics::MassKPlus; + + template + bool selectionEvent(const Coll& collision, bool fillHist = false) // default to false + { + if (fillHist) + hEventSelection.fill(HIST("hEventCut"), 0); + + if (std::abs(collision.posZ()) > selectionConfig.cutVrtxZ) + return false; + if (fillHist) + hEventSelection.fill(HIST("hEventCut"), 1); + + if (!collision.sel8()) + return false; + if (fillHist) + hEventSelection.fill(HIST("hEventCut"), 2); + + if (selectionConfig.isNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) + return false; + if (fillHist) + hEventSelection.fill(HIST("hEventCut"), 3); + + if (selectionConfig.isNoITSROFrameBorder && !collision.selection_bit(aod::evsel::kNoITSROFrameBorder)) + return false; + if (fillHist) + hEventSelection.fill(HIST("hEventCut"), 4); + + if (selectionConfig.isNoSameBunchPileup && (!collision.selection_bit(aod::evsel::kNoSameBunchPileup))) + return false; + if (fillHist) + hEventSelection.fill(HIST("hEventCut"), 5); + + if (selectionConfig.isGoodITSLayersAll && !collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) + return false; + if (fillHist) + hEventSelection.fill(HIST("hEventCut"), 6); + + if (selectionConfig.isNoCollInTimeRangeStandard && (!collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard))) + return false; + if (fillHist) + hEventSelection.fill(HIST("hEventCut"), 7); + + if (selectionConfig.isApplyOccCut && (std::abs(collision.trackOccupancyInTimeRange()) > selectionConfig.cfgOccCut)) + return false; + if (fillHist) + hEventSelection.fill(HIST("hEventCut"), 8); + + if (rctCut.requireRCTFlagChecker && !rctChecker(collision)) + return false; + if (fillHist) + hEventSelection.fill(HIST("hEventCut"), 9); + + if (selectionConfig.isTriggerTVX && !collision.selection_bit(aod::evsel::kIsTriggerTVX)) + return false; + if (fillHist) + hEventSelection.fill(HIST("hEventCut"), 10); + + if (selectionConfig.isGoodZvtxFT0vsPV && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) + return false; + if (fillHist) + hEventSelection.fill(HIST("hEventCut"), 11); + + if (selectionConfig.isVertexITSTPC && !collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) { + return false; + } + if (fillHist) + hEventSelection.fill(HIST("hEventCut"), 12); + + if (selectionConfig.isVertexTOFMatched && !collision.selection_bit(aod::evsel::kIsVertexTOFmatched)) { + return false; + } + if (fillHist) + hEventSelection.fill(HIST("hEventCut"), 13); + + return true; + } + + template + bool selectionTrack(const T& candidate) + { + if (selectionConfig.isGlobalTracks) { + if (!candidate.isGlobalTrackWoDCA()) + return false; + if (std::abs(candidate.pt()) < selectionConfig.cfgCutPT) + return false; + + if (std::abs(candidate.eta()) > selectionConfig.cfgCutEtaMax) + return false; + if (!selectionConfig.isApplyPtDepDCAxyCut) { + if (std::abs(candidate.dcaXY()) > selectionConfig.cfgCutDCAxyMax) + return false; + } else { + if (std::abs(candidate.dcaXY()) > (0.0105 + 0.035 / std::pow(candidate.pt(), 1.1))) + return false; + } + if (selectionConfig.isGoldenChi2 && !candidate.passedGoldenChi2()) + return false; + if (std::abs(candidate.dcaZ()) > selectionConfig.cfgCutDCAz) + return false; + if (candidate.tpcCrossedRowsOverFindableCls() < selectionConfig.cfgRCRFC) + return false; + if (candidate.itsNCls() < selectionConfig.cfgITScluster) + return false; + if (candidate.tpcNClsFound() < selectionConfig.cfgTPCcluster) + return false; + if (candidate.itsChi2NCl() >= selectionConfig.cfgITSChi2NCl) + return false; + if (candidate.tpcChi2NCl() >= selectionConfig.cfgTPCChi2NClMax || candidate.tpcChi2NCl() < selectionConfig.cfgTPCChi2NClMin) + return false; + if (selectionConfig.cfgPVContributor && !candidate.isPVContributor()) + return false; + if (selectionConfig.cfgUseITSTPCRefit && (!(o2::aod::track::ITSrefit) || !(o2::aod::track::TPCrefit))) + return false; + } else if (!selectionConfig.isGlobalTracks) { + if (std::abs(candidate.pt()) < selectionConfig.cfgCutPT) + return false; + // if (std::abs(candidate.eta()) > selectionConfig.cfgCutEtaMax || std::abs(candidate.eta()) < selectionConfig.cfgCutEtaMin) + if (std::abs(candidate.eta()) > selectionConfig.cfgCutEtaMax) + return false; + // if (std::abs(candidate.dcaXY()) > selectionConfig.cfgCutDCAxyMax || std::abs(candidate.dcaXY()) < selectionConfig.cfgCutDCAxyMin) + if (std::abs(candidate.dcaXY()) > selectionConfig.cfgCutDCAxyMax) + return false; + if (std::abs(candidate.dcaZ()) > selectionConfig.cfgCutDCAz) + return false; + if (candidate.tpcCrossedRowsOverFindableCls() < selectionConfig.cfgRCRFC) + return false; + if (candidate.itsNCls() < selectionConfig.cfgITScluster) + return false; + if (candidate.tpcNClsFound() < selectionConfig.cfgTPCcluster) + return false; + if (candidate.itsChi2NCl() >= selectionConfig.cfgITSChi2NCl) + return false; + if (candidate.tpcChi2NCl() >= selectionConfig.cfgTPCChi2NClMax || candidate.tpcChi2NCl() < selectionConfig.cfgTPCChi2NClMin) + return false; + if (selectionConfig.cfgPVContributor && !candidate.isPVContributor()) + return false; + if (selectionConfig.cfgPrimaryTrack && !candidate.isPrimaryTrack()) + return false; + } + + return true; + } + + // 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 (selectionConfig.isApplyDeepAngle && angle < selectionConfig.cfgDeepAngle) { + return false; + } + return true; + } + + template + bool selectionPID(const T& candidate, int PID) + { + if (PID == PIDParticle::kPion) { + if (onlyTOF) { + if (candidate.hasTOF() && std::abs(candidate.tofNSigmaPi()) < selectionConfig.nsigmaCutTOFPi && candidate.beta() > cBetaCutTOF) { + return true; + } + } else if (onlyTOFHIT) { + if (candidate.hasTOF() && std::abs(candidate.tofNSigmaPi()) < selectionConfig.nsigmaCutTOFPi && candidate.beta() > cBetaCutTOF) { + return true; + } + if (!candidate.hasTOF() && std::abs(candidate.tpcNSigmaPi()) < selectionConfig.nsigmaCutTPCPi) { + return true; + } + } else if (onlyTPC) { + if (std::abs(candidate.tpcNSigmaPi()) < selectionConfig.nsigmaCutTPCPi) { + return true; + } + } else { + if (candidate.hasTOF() && (candidate.tofNSigmaPi() * candidate.tofNSigmaPi() + candidate.tpcNSigmaPi() * candidate.tpcNSigmaPi()) < (selectionConfig.nsigmaCutCombinedPi * selectionConfig.nsigmaCutCombinedPi) && candidate.beta() > cBetaCutTOF) { + return true; + } + if (!candidate.hasTOF() && std::abs(candidate.tpcNSigmaPi()) < selectionConfig.nsigmaCutTPCPi) { + return true; + } + } + } else if (PID == PIDParticle::kKaon) { + if (onlyTOF) { + if (candidate.hasTOF() && std::abs(candidate.tofNSigmaKa()) < selectionConfig.nsigmaCutTOFKa && candidate.beta() > cBetaCutTOF) { + return true; + } + } else if (onlyTOFHIT) { + if (candidate.hasTOF() && std::abs(candidate.tofNSigmaKa()) < selectionConfig.nsigmaCutTOFKa && candidate.beta() > cBetaCutTOF) { + return true; + } + if (!candidate.hasTOF() && std::abs(candidate.tpcNSigmaKa()) < selectionConfig.nsigmaCutTPCKa) { + return true; + } + } else if (onlyTPC) { + if (std::abs(candidate.tpcNSigmaKa()) < selectionConfig.nsigmaCutTPCKa) { + return true; + } + } else { + if (candidate.hasTOF() && (candidate.tofNSigmaKa() * candidate.tofNSigmaKa() + candidate.tpcNSigmaKa() * candidate.tpcNSigmaKa()) < (selectionConfig.nsigmaCutCombinedKa * selectionConfig.nsigmaCutCombinedKa) && candidate.beta() > cBetaCutTOF) { + return true; + } + if (!candidate.hasTOF() && std::abs(candidate.tpcNSigmaKa()) < selectionConfig.nsigmaCutTPCKa) { + return true; + } + } + } + return false; + } + + template + bool selectionPIDNew(const T& candidate, int PID) + { + if (PID == PIDParticle::kPion) { + if (candidate.pt() < selectionConfig.lowPtCutPID && std::abs(candidate.tpcNSigmaPi()) < selectionConfig.nsigmaCutTPCPi) { + return true; + } + if (candidate.pt() >= selectionConfig.lowPtCutPID && std::abs(candidate.tpcNSigmaPi()) < selectionConfig.nsigmaCutTPCPi && candidate.hasTOF() && std::abs(candidate.tofNSigmaPi()) < selectionConfig.nsigmaCutTOFPi) { + return true; + } + if (candidate.pt() >= selectionConfig.lowPtCutPID && std::abs(candidate.tpcNSigmaPi()) < selectionConfig.nsigmaCutTPCPi && !candidate.hasTOF()) { + return true; + } + } else if (PID == PIDParticle::kKaon) { + if (candidate.pt() < selectionConfig.lowPtCutPID && std::abs(candidate.tpcNSigmaKa()) < selectionConfig.nsigmaCutTPCKa) { + return true; + } + if (candidate.pt() >= selectionConfig.lowPtCutPID && std::abs(candidate.tpcNSigmaKa()) < selectionConfig.nsigmaCutTPCKa && candidate.hasTOF() && std::abs(candidate.tofNSigmaKa()) < selectionConfig.nsigmaCutTOFKa) { + return true; + } + if (candidate.pt() >= selectionConfig.lowPtCutPID && std::abs(candidate.tpcNSigmaKa()) < selectionConfig.nsigmaCutTPCKa && !candidate.hasTOF()) { + return true; + } + } + return false; + } + + // Defining filters for events (event selection) + // Processed events will be already fulfilling the event selection + // requirements + // Filter eventFilter = (o2::aod::evsel::sel8 == true); + Filter posZFilter = (nabs(o2::aod::collision::posZ) < selectionConfig.cutVrtxZ); + + // Filter acceptanceFilter = (nabs(aod::track::eta) < selectionConfig.cfgCutEtaMax && nabs(aod::track::pt) > selectionConfig.cfgCutPT) && (nabs(aod::track::eta) > selectionConfig.cfgCutEtaMin); + Filter acceptanceFilter = (nabs(aod::track::eta) < selectionConfig.cfgCutEtaMax && nabs(aod::track::pt) > selectionConfig.cfgCutPT); + // Filter fDCAcutFilter = (nabs(aod::track::dcaXY) < selectionConfig.cfgCutDCAxyMax) && (nabs(aod::track::dcaZ) < selectionConfig.cfgCutDCAz) && (nabs(aod::track::dcaXY) > selectionConfig.cfgCutDCAxyMin); + Filter fDCAcutFilter = (nabs(aod::track::dcaXY) < selectionConfig.cfgCutDCAxyMax) && (nabs(aod::track::dcaZ) < selectionConfig.cfgCutDCAz); + + using EventCandidates = soa::Filtered>; // aod::CentNGlobals, aod::CentNTPVs, aod::CentMFTs + using EventCandidatesMix = soa::Filtered>; // aod::CentNGlobals, aod::CentNTPVs, aod::CentMFTs + using TrackCandidates = soa::Filtered>; + using EventCandidatesMC = soa::Join; + // using EventCandidatesMC = soa::Filtered>; + + using TrackCandidatesMC = soa::Filtered>; + // using EventMCGenerated = soa::Join; // aod::CentNGlobals, aod::CentNTPVs, aod::CentMFTs + using EventMCGenerated = soa::Join; + + //*********Varibles declaration*************** + float centrality{-1.0}, theta2; + ROOT::Math::PxPyPzMVector daughter1, daughter2, daughterRot, mother, motherRot; + bool isMix = false; + + template + void fillInvMass(const T1& daughter1, const T1& daughter2, const T1& mother, float centrality, bool isMix, const T2& track1, const T2& track2) + { + if (track1.sign() * track2.sign() < 0) { + if (!isMix) { + if (std::abs(mother.Rapidity()) < selectionConfig.motherRapidityCut) { + hInvMass.fill(HIST("h3KstarInvMassUnlikeSign"), centrality, mother.Pt(), mother.M()); + } + for (int i = 0; i < cRotations; i++) { + theta2 = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / rotationalCut, o2::constants::math::PI + o2::constants::math::PI / 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()); + motherRot = daughterRot + daughter2; + + if (calcRotational && std::abs(motherRot.Rapidity()) < selectionConfig.motherRapidityCut) + hInvMass.fill(HIST("h3KstarInvMassRotated"), centrality, motherRot.Pt(), motherRot.M()); + } + } else if (isMix && std::abs(mother.Rapidity()) < selectionConfig.motherRapidityCut) { + hInvMass.fill(HIST("h3KstarInvMassMixed"), centrality, mother.Pt(), mother.M()); + } + } else { + if (!isMix) { + if (calcLikeSign && std::abs(mother.Rapidity()) < selectionConfig.motherRapidityCut) { + if (track1.sign() > 0 && track2.sign() > 0) { + hInvMass.fill(HIST("h3KstarInvMasslikeSignPP"), centrality, mother.Pt(), mother.M()); + } else if (track1.sign() < 0 && track2.sign() < 0) { + hInvMass.fill(HIST("h3KstarInvMasslikeSignMM"), centrality, mother.Pt(), mother.M()); + } + } + } + } + } + + void processSE(EventCandidates::iterator const& collision, TrackCandidates const& tracks, aod::BCs const&) + { + int occupancy = collision.trackOccupancyInTimeRange(); + hEventSelection.fill(HIST("hOccupancy"), occupancy); + + if (!selectionEvent(collision, true)) { // fill data event cut histogram + return; + } + + centrality = -1; + + if (cSelectMultEstimator == kFT0M) { + centrality = collision.centFT0M(); + } else if (cSelectMultEstimator == kFT0A) { + centrality = collision.centFT0A(); + } else if (cSelectMultEstimator == kFT0C) { + centrality = collision.centFT0C(); + } else if (cSelectMultEstimator == kFV0A) { + centrality = collision.centFV0A(); + } else { + centrality = collision.centFT0M(); // default + } + + /* else if (cSelectMultEstimator == 4) { + centrality = collision.centMFT(); + } */ + /* else if (cSelectMultEstimator == 5) { + centrality = collision.centNGlobal(); + } */ + /* else if (cSelectMultEstimator == 6) { + centrality = collision.centNTPV(); + } */ + + // Fill the event counter + if (cQAevents) { + hEventSelection.fill(HIST("hVertexZ"), collision.posZ()); + hEventSelection.fill(HIST("hCentrality"), centrality); + } + + for (const auto& [track1, track2] : combinations(CombinationsFullIndexPolicy(tracks, tracks))) { + hEventSelection.fill(HIST("tracksCheckData"), 0.5); + if (!selectionTrack(track1) || !selectionTrack(track2)) { + continue; + } + hEventSelection.fill(HIST("tracksCheckData"), 1.5); + + if (track1.globalIndex() == track2.globalIndex()) + continue; + + if (!selectionPair(track1, track2)) { + continue; + } + hEventSelection.fill(HIST("tracksCheckData"), 6.5); + + if (cQAplots) { + hOthers.fill(HIST("hCRFC_before"), track1.tpcCrossedRowsOverFindableCls()); + hOthers.fill(HIST("dE_by_dx_TPC"), track1.p(), track1.tpcSignal()); + + if (track1.sign() < 0) { + hPID.fill(HIST("Before/hTPCnsigKa_Neg_mult_pt"), track1.tpcNSigmaKa(), centrality, track1.pt()); + hPID.fill(HIST("Before/hTPCnsigPi_Neg_mult_pt"), track1.tpcNSigmaPi(), centrality, track1.pt()); + hPID.fill(HIST("Before/hTOFnsigKa_Neg_mult_pt"), track1.tofNSigmaKa(), centrality, track1.pt()); + hPID.fill(HIST("Before/hTOFnsigPi_Neg_mult_pt"), track1.tofNSigmaPi(), centrality, track1.pt()); + } else if (track1.sign() > 0) { + hPID.fill(HIST("Before/hTPCnsigKa_Pos_mult_pt"), track1.tpcNSigmaKa(), centrality, track1.pt()); + hPID.fill(HIST("Before/hTPCnsigPi_Pos_mult_pt"), track1.tpcNSigmaPi(), centrality, track1.pt()); + hPID.fill(HIST("Before/hTOFnsigKa_Pos_mult_pt"), track1.tofNSigmaKa(), centrality, track1.pt()); + hPID.fill(HIST("Before/hTOFnsigPi_Pos_mult_pt"), track1.tofNSigmaPi(), centrality, track1.pt()); + } + + hPID.fill(HIST("Before/hNsigma_TPC_TOF_Ka_pt"), track1.tpcNSigmaKa(), track1.tofNSigmaKa(), track1.pt()); + hPID.fill(HIST("Before/hNsigma_TPC_TOF_Pi_pt"), track1.tpcNSigmaPi(), track1.tofNSigmaPi(), track1.pt()); + } + + if (cQAevents) { + hOthers.fill(HIST("hDcaxy_cent_pt"), track1.dcaXY(), centrality, track1.pt()); + hOthers.fill(HIST("hDcaz_cent_pt"), track1.dcaZ(), centrality, track1.pt()); + } + + // since we are using combinations full index policy, so repeated pairs are allowed, so we can check one with Kaon and other with pion + if (!applypTdepPID && !selectionPID(track1, 1)) // Track 1 is checked with Kaon + continue; + if (!applypTdepPID && !selectionPID(track2, 0)) // Track 2 is checked with Pion + continue; + + if (applypTdepPID && !selectionPIDNew(track1, 1)) // Track 1 is checked with Kaon + continue; + if (applypTdepPID && !selectionPIDNew(track2, 0)) // Track 2 is checked with Pion + continue; + + hEventSelection.fill(HIST("tracksCheckData"), 2.5); + + hEventSelection.fill(HIST("tracksCheckData"), 5.5); + + if (cQAplots) { + hOthers.fill(HIST("hEta_after"), track1.eta()); + hOthers.fill(HIST("hCRFC_after"), track1.tpcCrossedRowsOverFindableCls()); + hOthers.fill(HIST("hDcaxyPi"), track2.dcaXY()); + hOthers.fill(HIST("hDcaxyKa"), track1.dcaXY()); + hOthers.fill(HIST("hDcazPi"), track2.dcaZ()); + hOthers.fill(HIST("hDcazKa"), track1.dcaZ()); + + if (track1.sign() < 0) { + hPID.fill(HIST("After/hTPCnsigKa_Neg_mult_pt"), track1.tpcNSigmaKa(), centrality, track1.pt()); + hPID.fill(HIST("After/hTOFnsigKa_Neg_mult_pt"), track1.tofNSigmaKa(), centrality, track1.pt()); + } else if (track1.sign() > 0) { + hPID.fill(HIST("After/hTPCnsigKa_Pos_mult_pt"), track1.tpcNSigmaKa(), centrality, track1.pt()); + hPID.fill(HIST("After/hTOFnsigKa_Pos_mult_pt"), track1.tofNSigmaKa(), centrality, track1.pt()); + } + + if (track2.sign() < 0) { + hPID.fill(HIST("After/hTPCnsigPi_Neg_mult_pt"), track2.tpcNSigmaPi(), centrality, track2.pt()); + hPID.fill(HIST("After/hTOFnsigPi_Neg_mult_pt"), track2.tofNSigmaPi(), centrality, track2.pt()); + } else if (track2.sign() > 0) { + hPID.fill(HIST("After/hTPCnsigPi_Pos_mult_pt"), track2.tpcNSigmaPi(), centrality, track2.pt()); + hPID.fill(HIST("After/hTOFnsigPi_Pos_mult_pt"), track2.tofNSigmaPi(), centrality, track2.pt()); + } + + hPID.fill(HIST("After/hNsigma_TPC_TOF_Ka_pt"), track1.tpcNSigmaKa(), track1.tofNSigmaKa(), track1.pt()); + hPID.fill(HIST("After/hNsigma_TPC_TOF_Pi_pt"), track2.tpcNSigmaPi(), track2.tofNSigmaPi(), track2.pt()); + } + + daughter1 = ROOT::Math::PxPyPzMVector(track1.px(), track1.py(), track1.pz(), massKa); + daughter2 = ROOT::Math::PxPyPzMVector(track2.px(), track2.py(), track2.pz(), massPi); + mother = daughter1 + daughter2; // Kstar meson + + /* if (selectionConfig.isApplyCutsOnMother) { + if (mother.Pt() >= selectionConfig.cMaxPtMotherCut) // excluding candidates in overflow + continue; + if (mother.M() >= selectionConfig.cMaxMinvMotherCut) // excluding candidates in overflow + continue; + } */ + + hOthers.fill(HIST("hKstar_rap_pt"), mother.Rapidity(), mother.Pt()); + hOthers.fill(HIST("hKstar_eta_pt"), mother.Eta(), mother.Pt()); + + isMix = false; + fillInvMass(daughter1, daughter2, mother, centrality, isMix, track1, track2); + } + } + PROCESS_SWITCH(Kstar892LightIon, processSE, "Process Same event", true); + + ConfigurableAxis axisVertex{"axisVertex", {20, -10, 10}, "vertex axis for ME mixing"}; + // ConfigurableAxis axisCentralityClass{"axisCentralityClass", {10, 0, 100}, "centrality percentile for ME mixing"}; + ConfigurableAxis axisCentrality{"axisCentrality", {2000, 0, 10000}, "TPC centrality axis for ME mixing"}; + + // using BinningTypeTPCcentrality = ColumnBinningPolicy; + using BinningTypeFT0M = ColumnBinningPolicy; + using BinningTypeFT0A = ColumnBinningPolicy; + using BinningTypeFT0C = ColumnBinningPolicy; + using BinningTypeFV0A = ColumnBinningPolicy; + + BinningTypeFT0M binningOnFT0M{{axisVertex, axisCentrality}, true}; + BinningTypeFT0A binningOnFT0A{{axisVertex, axisCentrality}, true}; + BinningTypeFT0C binningOnFT0C{{axisVertex, axisCentrality}, true}; + BinningTypeFV0A binningOnFV0A{{axisVertex, axisCentrality}, true}; + + SameKindPair pair1{binningOnFT0M, selectionConfig.cfgNoMixedEvents, -1, &cache}; + SameKindPair pair2{binningOnFT0A, selectionConfig.cfgNoMixedEvents, -1, &cache}; + SameKindPair pair3{binningOnFT0C, selectionConfig.cfgNoMixedEvents, -1, &cache}; + SameKindPair pair4{binningOnFV0A, selectionConfig.cfgNoMixedEvents, -1, &cache}; + + void processME(EventCandidatesMix const&, TrackCandidates const&) + { + // Map estimator to pair and centrality accessor + auto runMixing = [&](auto& pair, auto centralityGetter) { + for (const auto& [c1, tracks1, c2, tracks2] : pair) { + + if (!selectionEvent(c1, false) || !selectionEvent(c2, false)) { // don't fill event cut histogram + continue; + } + + centrality = centralityGetter(c1); + + for (const auto& [t1, t2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { + if (!selectionTrack(t1) || !selectionTrack(t2)) + continue; + if (!selectionPID(t1, 1) || !selectionPID(t2, 0)) + continue; + + if (!selectionPair(t1, t2)) { + continue; + } + + daughter1 = ROOT::Math::PxPyPzMVector(t1.px(), t1.py(), t1.pz(), massKa); + daughter2 = ROOT::Math::PxPyPzMVector(t2.px(), t2.py(), t2.pz(), massPi); + mother = daughter1 + daughter2; + + isMix = true; + + if (std::abs(mother.Rapidity()) < selectionConfig.motherRapidityCut) { + fillInvMass(daughter1, daughter2, mother, centrality, isMix, t1, t2); + } + } + } + }; + + // Call mixing based on selected estimator + if (cSelectMultEstimator == kFT0M) { + runMixing(pair1, [](const auto& c) { return c.centFT0M(); }); + } else if (cSelectMultEstimator == kFT0A) { + runMixing(pair2, [](const auto& c) { return c.centFT0A(); }); + } else if (cSelectMultEstimator == kFT0C) { + runMixing(pair3, [](const auto& c) { return c.centFT0C(); }); + } else if (cSelectMultEstimator == kFV0A) { + runMixing(pair4, [](const auto& c) { return c.centFV0A(); }); + } + } + PROCESS_SWITCH(Kstar892LightIon, processME, "Process Mixed event", true); + + using BinningTypeMCFT0M = ColumnBinningPolicy; + using BinningTypeMCFT0A = ColumnBinningPolicy; + using BinningTypeMCFT0C = ColumnBinningPolicy; + using BinningTypeMCFV0A = ColumnBinningPolicy; + + BinningTypeMCFT0M binningOnMCFT0M{{axisVertex, axisCentrality}, true}; + BinningTypeMCFT0A binningOnMCFT0A{{axisVertex, axisCentrality}, true}; + BinningTypeMCFT0C binningOnMCFT0C{{axisVertex, axisCentrality}, true}; + BinningTypeMCFV0A binningOnMCFV0A{{axisVertex, axisCentrality}, true}; + + SameKindPair pairmc1{binningOnMCFT0M, selectionConfig.cfgNoMixedEvents, -1, &cache}; + SameKindPair pairmc2{binningOnMCFT0A, selectionConfig.cfgNoMixedEvents, -1, &cache}; + SameKindPair pairmc3{binningOnMCFT0C, selectionConfig.cfgNoMixedEvents, -1, &cache}; + SameKindPair pairmc4{binningOnMCFV0A, selectionConfig.cfgNoMixedEvents, -1, &cache}; + + void processMEMC(EventCandidatesMC const&, TrackCandidatesMC const&, aod::McParticles const&, aod::McCollisions const&) + { + auto runMixing = [&](auto& pair, auto centralityGetter) { + for (const auto& [c1, tracks1, c2, tracks2] : pair) { + + if (!selectionEvent(c1, false) || !selectionEvent(c2, false)) { // don't fill event cut histogram + continue; + } + + if (!c1.has_mcCollision() || !c2.has_mcCollision()) { + continue; // skip if no MC collision associated + } + + centrality = centralityGetter(c1); + + for (const auto& [t1, t2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { + if (!selectionTrack(t1) || !selectionTrack(t2)) + continue; + + if (!selectionPID(t1, 1) || !selectionPID(t2, 0)) + continue; + + if (!t1.has_mcParticle() || !t2.has_mcParticle()) { + continue; // skip if no MC particle associated + } + + const auto mctrack1 = t1.mcParticle(); + const auto mctrack2 = t2.mcParticle(); + + if (!mctrack1.isPhysicalPrimary() || !mctrack2.isPhysicalPrimary()) { + continue; + } + + daughter1 = ROOT::Math::PxPyPzMVector(t1.px(), t1.py(), t1.pz(), massKa); + daughter2 = ROOT::Math::PxPyPzMVector(t2.px(), t2.py(), t2.pz(), massPi); + mother = daughter1 + daughter2; + + isMix = true; + + if (std::abs(mother.Rapidity()) < selectionConfig.motherRapidityCut) { + fillInvMass(daughter1, daughter2, mother, centrality, isMix, t1, t2); + } + } + } + }; + // Call mixing based on selected estimator + if (cSelectMultEstimator == kFT0M) { + runMixing(pairmc1, [](const auto& c) { return c.centFT0M(); }); + } else if (cSelectMultEstimator == kFT0A) { + runMixing(pairmc2, [](const auto& c) { return c.centFT0A(); }); + } else if (cSelectMultEstimator == kFT0C) { + runMixing(pairmc3, [](const auto& c) { return c.centFT0C(); }); + } else if (cSelectMultEstimator == kFV0A) { + runMixing(pairmc4, [](const auto& c) { return c.centFV0A(); }); + } + } + PROCESS_SWITCH(Kstar892LightIon, processMEMC, "Process mixed-event in MC", false); + + void processSEMC(EventCandidatesMC::iterator const& collision, TrackCandidatesMC const& tracks, aod::McParticles const&, aod::McCollisions const& /*mcCollisions*/) + { + // auto oldindex = -999; + if (!collision.has_mcCollision()) { + return; + } + int occupancy = collision.trackOccupancyInTimeRange(); + hEventSelection.fill(HIST("hOccupancy"), occupancy); + + if (!selectionEvent(collision, false)) { // don't fill event cut histogram + return; + } + + centrality = -1; + + if (cSelectMultEstimator == kFT0M) { + centrality = collision.centFT0M(); + } else if (cSelectMultEstimator == kFT0A) { + centrality = collision.centFT0A(); + } else if (cSelectMultEstimator == kFT0C) { + centrality = collision.centFT0C(); + } else if (cSelectMultEstimator == kFV0A) { + centrality = collision.centFV0A(); + } else { + centrality = collision.centFT0M(); // default + } + + // Fill the event counter + if (cQAevents) { + hEventSelection.fill(HIST("hVertexZ"), collision.posZ()); + hEventSelection.fill(HIST("hCentrality"), centrality); + } + + for (const auto& [track1, track2] : combinations(CombinationsFullIndexPolicy(tracks, tracks))) { + hEventSelection.fill(HIST("tracksCheckData"), 0.5); + if (!selectionTrack(track1) || !selectionTrack(track2)) + continue; + + const auto mctrack1 = track1.mcParticle(); + const auto mctrack2 = track2.mcParticle(); + + if (!track1.has_mcParticle() || !track2.has_mcParticle()) + continue; // skip if no MC particle associated + + if (!mctrack1.isPhysicalPrimary() || !mctrack2.isPhysicalPrimary()) + continue; + + hEventSelection.fill(HIST("tracksCheckData"), 1.5); + + if (track1.globalIndex() == track2.globalIndex()) + continue; + hEventSelection.fill(HIST("tracksCheckData"), 6.5); + + if (cQAplots) { + hOthers.fill(HIST("hCRFC_before"), track1.tpcCrossedRowsOverFindableCls()); + hOthers.fill(HIST("dE_by_dx_TPC"), track1.p(), track1.tpcSignal()); + + if (track1.sign() < 0) { + hPID.fill(HIST("Before/hTPCnsigKa_Neg_mult_pt"), track1.tpcNSigmaKa(), centrality, track1.pt()); + hPID.fill(HIST("Before/hTPCnsigPi_Neg_mult_pt"), track1.tpcNSigmaPi(), centrality, track1.pt()); + hPID.fill(HIST("Before/hTOFnsigKa_Neg_mult_pt"), track1.tofNSigmaKa(), centrality, track1.pt()); + hPID.fill(HIST("Before/hTOFnsigPi_Neg_mult_pt"), track1.tofNSigmaPi(), centrality, track1.pt()); + } else { + hPID.fill(HIST("Before/hTPCnsigKa_Pos_mult_pt"), track1.tpcNSigmaKa(), centrality, track1.pt()); + hPID.fill(HIST("Before/hTPCnsigPi_Pos_mult_pt"), track1.tpcNSigmaPi(), centrality, track1.pt()); + hPID.fill(HIST("Before/hTOFnsigKa_Pos_mult_pt"), track1.tofNSigmaKa(), centrality, track1.pt()); + hPID.fill(HIST("Before/hTOFnsigPi_Pos_mult_pt"), track1.tofNSigmaPi(), centrality, track1.pt()); + } + + hPID.fill(HIST("Before/hNsigma_TPC_TOF_Ka_pt"), track1.tpcNSigmaKa(), track1.tofNSigmaKa(), track1.pt()); + hPID.fill(HIST("Before/hNsigma_TPC_TOF_Pi_pt"), track1.tpcNSigmaPi(), track1.tofNSigmaPi(), track1.pt()); + } + + if (cQAevents) { + hOthers.fill(HIST("hDcaxy_cent_pt"), track1.dcaXY(), centrality, track1.pt()); + hOthers.fill(HIST("hDcaz_cent_pt"), track1.dcaZ(), centrality, track1.pt()); + } + + // since we are using combinations full index policy, so repeated pairs are allowed, so we can check one with Kaon and other with pion + if (!applypTdepPID && (!selectionPID(track1, 1) || !selectionPID(track2, 0))) // Track 1 is checked with Kaon, track 2 is checked with Pion + continue; + hEventSelection.fill(HIST("tracksCheckData"), 2.5); + + if (applypTdepPID && (!selectionPIDNew(track1, 1) || !selectionPIDNew(track2, 0))) // Track 1 is checked with Kaon, track 2 is checked with Pion + continue; + hEventSelection.fill(HIST("tracksCheckData"), 3.5); + + hEventSelection.fill(HIST("tracksCheckData"), 4.5); + + if (cQAplots) { + hOthers.fill(HIST("hEta_after"), track1.eta()); + hOthers.fill(HIST("hCRFC_after"), track1.tpcCrossedRowsOverFindableCls()); + hOthers.fill(HIST("hDcaxyPi"), track2.dcaXY()); + hOthers.fill(HIST("hDcaxyKa"), track1.dcaXY()); + hOthers.fill(HIST("hDcazPi"), track2.dcaZ()); + hOthers.fill(HIST("hDcazKa"), track1.dcaZ()); + + if (track1.sign() < 0) { + hPID.fill(HIST("After/hTPCnsigKa_Neg_mult_pt"), track1.tpcNSigmaKa(), centrality, track1.pt()); + hPID.fill(HIST("After/hTOFnsigKa_Neg_mult_pt"), track1.tofNSigmaKa(), centrality, track1.pt()); + } else if (track1.sign() > 0) { + hPID.fill(HIST("After/hTPCnsigKa_Pos_mult_pt"), track1.tpcNSigmaKa(), centrality, track1.pt()); + hPID.fill(HIST("After/hTOFnsigKa_Pos_mult_pt"), track1.tofNSigmaKa(), centrality, track1.pt()); + } + + if (track2.sign() < 0) { + hPID.fill(HIST("After/hTPCnsigPi_Neg_mult_pt"), track2.tpcNSigmaPi(), centrality, track2.pt()); + hPID.fill(HIST("After/hTOFnsigPi_Neg_mult_pt"), track2.tofNSigmaPi(), centrality, track2.pt()); + } else if (track2.sign() > 0) { + hPID.fill(HIST("After/hTPCnsigPi_Pos_mult_pt"), track2.tpcNSigmaPi(), centrality, track2.pt()); + hPID.fill(HIST("After/hTOFnsigPi_Pos_mult_pt"), track2.tofNSigmaPi(), centrality, track2.pt()); + } + + hPID.fill(HIST("After/hNsigma_TPC_TOF_Ka_pt"), track1.tpcNSigmaKa(), track1.tofNSigmaKa(), track1.pt()); + hPID.fill(HIST("After/hNsigma_TPC_TOF_Pi_pt"), track2.tpcNSigmaPi(), track2.tofNSigmaPi(), track2.pt()); + } + daughter1 = ROOT::Math::PxPyPzMVector(track1.px(), track1.py(), track1.pz(), massKa); + daughter2 = ROOT::Math::PxPyPzMVector(track2.px(), track2.py(), track2.pz(), massPi); + mother = daughter1 + daughter2; // Kstar meson + + /* if (selectionConfig.isApplyCutsOnMother) { + if (mother.Pt() >= selectionConfig.cMaxPtMotherCut) // excluding candidates in overflow + continue; + if (mother.M() >= selectionConfig.cMaxMinvMotherCut) // excluding candidates in overflow + continue; + } */ + + hOthers.fill(HIST("hKstar_rap_pt"), mother.Rapidity(), mother.Pt()); + hOthers.fill(HIST("hKstar_eta_pt"), mother.Eta(), mother.Pt()); + + isMix = false; + fillInvMass(daughter1, daughter2, mother, centrality, isMix, track1, track2); + } + } + PROCESS_SWITCH(Kstar892LightIon, processSEMC, "Process same event in MC", false); + + void processGen(EventMCGenerated::iterator const& mcCollision, aod::McParticles const& mcParticles, const soa::SmallGroups& collisions) + { + hEventSelection.fill(HIST("eventsCheckGen"), 0.5); + + if (std::abs(mcCollision.posZ()) < selectionConfig.cutVrtxZ) + hEventSelection.fill(HIST("eventsCheckGen"), 1.5); + + std::vector selectedEvents(collisions.size()); + int nevts = 0; + centrality = -1.0; + + hEventSelection.fill(HIST("eventsCheckGen"), 2.5); + + for (const auto& collision : collisions) { + if (!selectionEvent(collision, false)) { // don't fill event cut histogram + continue; + } + centrality = collision.centFT0M(); + + if (cSelectMultEstimator == kFT0M) { + centrality = collision.centFT0M(); + } else if (cSelectMultEstimator == kFT0A) { + centrality = collision.centFT0A(); + } else if (cSelectMultEstimator == kFT0C) { + centrality = collision.centFT0C(); + } else if (cSelectMultEstimator == kFV0A) { + centrality = collision.centFV0A(); + } else { + centrality = collision.centFT0M(); // default + } + hMC.fill(HIST("h1GenCent"), centrality); + + int occupancy = collision.trackOccupancyInTimeRange(); + hEventSelection.fill(HIST("hOccupancy"), occupancy); + + selectedEvents[nevts++] = collision.mcCollision_as().globalIndex(); + } + selectedEvents.resize(nevts); + + for (const auto& mcParticle : mcParticles) { + if (std::abs(mcParticle.y()) < selectionConfig.motherRapidityCut && std::abs(mcParticle.pdgCode()) == o2::constants::physics::kK0Star892) { + hMC.fill(HIST("hAllKstarGenCollisisons"), mcParticle.pt(), centrality); + } + } + + const auto evtReconstructedAndSelected = std::find(selectedEvents.begin(), selectedEvents.end(), mcCollision.globalIndex()) != selectedEvents.end(); + hMC.fill(HIST("hAllGenCollisions"), centrality); + if (!evtReconstructedAndSelected) { // Check that the event is reconstructed and that the reconstructed events pass the selection + return; + } + + hMC.fill(HIST("hAllGenCollisions1Rec"), centrality); + hEventSelection.fill(HIST("eventsCheckGen"), 3.5); + + for (const auto& mcParticle : mcParticles) { + + if (std::abs(mcParticle.y()) >= selectionConfig.motherRapidityCut) { + continue; + } + + /* if (selectionConfig.isApplyCutsOnMother) { + if (mcParticle.pt() >= selectionConfig.cMaxPtMotherCut) // excluding candidates in overflow + continue; + if ((std::sqrt(mcParticle.e() * mcParticle.e() - mcParticle.p() * mcParticle.p())) >= selectionConfig.cMaxMinvMotherCut) // excluding candidates in overflow + continue; + } */ + + if (std::abs(mcParticle.pdgCode()) != o2::constants::physics::kK0Star892) { + continue; + } + hMC.fill(HIST("hAllKstarGenCollisisons1Rec"), mcParticle.pt(), centrality); + + auto kDaughters = mcParticle.daughters_as(); + if (kDaughters.size() != noOfDaughters) { + continue; + } + + bool hasPos = false; + bool hasNeg = false; + + auto passkaon = false; + auto passpion = false; + + for (const auto& kCurrentDaughter : kDaughters) { + if (!kCurrentDaughter.isPhysicalPrimary()) { + continue; + } + + int pdgDau = kCurrentDaughter.pdgCode(); + int charge = (pdgDau > 0) - (pdgDau < 0); + + if (charge > 0) + hasPos = true; + if (charge < 0) + hasNeg = true; + + if (std::abs(pdgDau) == PDG_t::kKPlus) { + passkaon = true; + } else if (std::abs(pdgDau) == PDG_t::kPiPlus) { + passpion = true; + } + } + + if ((passkaon && passpion) && (hasPos && hasNeg)) { + hMC.fill(HIST("hk892GenpT"), mcParticle.pt(), centrality); + } + } + } + PROCESS_SWITCH(Kstar892LightIon, processGen, "Process Generated", false); + + void processRec(EventCandidatesMC::iterator const& collision, TrackCandidatesMC const& tracks, aod::McParticles const&, EventMCGenerated const&) + { + if (!collision.has_mcCollision()) { + return; + } + + // centrality = collision.centFT0M(); + + if (cSelectMultEstimator == kFT0M) { + centrality = collision.centFT0M(); + } else if (cSelectMultEstimator == kFT0A) { + centrality = collision.centFT0A(); + } else if (cSelectMultEstimator == kFT0C) { + centrality = collision.centFT0C(); + } else if (cSelectMultEstimator == kFV0A) { + centrality = collision.centFV0A(); + } else { + centrality = collision.centFT0M(); // default + } + + hMC.fill(HIST("hAllRecCollisions"), centrality); + + if (!selectionEvent(collision, true)) { // fill MC event cut histogram + return; + } + + hMC.fill(HIST("h1RecCent"), centrality); + + if (cQAevents) { + hEventSelection.fill(HIST("hVertexZ"), collision.posZ()); + } + + auto oldindex = -999; + + for (const auto& [track1, track2] : combinations(CombinationsFullIndexPolicy(tracks, tracks))) { + if (!selectionTrack(track1) || !selectionTrack(track2)) { + continue; + } + + if (!track1.has_mcParticle() || !track2.has_mcParticle()) { + continue; + } + + if (track1.index() == track2.index()) + continue; + + if (indexCheck && (track2.index() < track1.index())) + continue; + + hEventSelection.fill(HIST("recMCparticles"), 6.5); + + if (!selectionPair(track1, track2)) { + continue; + } + hEventSelection.fill(HIST("recMCparticles"), 16.5); + + if (cQAevents) { + hOthers.fill(HIST("hDcaxy_cent_pt"), track1.dcaXY(), centrality, track1.pt()); + hOthers.fill(HIST("hDcaz_cent_pt"), track1.dcaZ(), centrality, track1.pt()); + } + + if (track1.sign() * track2.sign() >= 0) { + continue; + } + hEventSelection.fill(HIST("recMCparticles"), 4.5); + + const auto mctrack1 = track1.mcParticle(); + const auto mctrack2 = track2.mcParticle(); + if (!mctrack1.isPhysicalPrimary() || !mctrack2.isPhysicalPrimary()) { + continue; + } + + int track1PDG = std::abs(mctrack1.pdgCode()); + int track2PDG = std::abs(mctrack2.pdgCode()); + + if (cQAplots) { + if (mctrack2.pdgCode() == PDG_t::kPiPlus) { // pion + hPID.fill(HIST("Before/hTPCnsigPi_Pos_mult_pt"), track2.tpcNSigmaPi(), centrality, track2.pt()); + hPID.fill(HIST("Before/hTOFnsigPi_Pos_mult_pt"), track2.tofNSigmaPi(), centrality, track2.pt()); + } + if (mctrack2.pdgCode() == PDG_t::kKPlus) { // kaon + hPID.fill(HIST("Before/hTPCnsigKa_Pos_mult_pt"), track2.tpcNSigmaKa(), centrality, track2.pt()); + hPID.fill(HIST("Before/hTOFnsigKa_Pos_mult_pt"), track2.tofNSigmaKa(), centrality, track2.pt()); + } + if (mctrack2.pdgCode() == PDG_t::kPiMinus) { // negative track pion + hPID.fill(HIST("Before/hTPCnsigPi_Neg_mult_pt"), track2.tpcNSigmaPi(), centrality, track2.pt()); + hPID.fill(HIST("Before/hTOFnsigPi_Neg_mult_pt"), track2.tofNSigmaPi(), centrality, track2.pt()); + } + if (mctrack2.pdgCode() == PDG_t::kKMinus) { // negative track kaon + hPID.fill(HIST("Before/hTPCnsigKa_Neg_mult_pt"), track2.tpcNSigmaKa(), centrality, track2.pt()); + hPID.fill(HIST("Before/hTOFnsigKa_Neg_mult_pt"), track2.tofNSigmaKa(), centrality, track2.pt()); + } + if (std::abs(mctrack1.pdgCode()) == PDG_t::kKPlus && std::abs(mctrack2.pdgCode()) == PDG_t::kPiPlus) { + hPID.fill(HIST("Before/hNsigma_TPC_TOF_Ka_pt"), track1.tpcNSigmaKa(), track1.tofNSigmaKa(), track1.pt()); + hPID.fill(HIST("Before/hNsigma_TPC_TOF_Pi_pt"), track2.tpcNSigmaPi(), track2.tofNSigmaPi(), track2.pt()); + } + } + + if (!(track1PDG == PDG_t::kPiPlus && track2PDG == PDG_t::kKPlus) && !(track1PDG == PDG_t::kKPlus && track2PDG == PDG_t::kPiPlus)) { + continue; + } + hEventSelection.fill(HIST("recMCparticles"), 7.5); + + for (const auto& mothertrack1 : mctrack1.mothers_as()) { + for (const auto& mothertrack2 : mctrack2.mothers_as()) { + if (mothertrack1.pdgCode() != mothertrack2.pdgCode()) { + continue; + } + + if (mothertrack1.globalIndex() != mothertrack2.globalIndex()) { + continue; + } + hEventSelection.fill(HIST("recMCparticles"), 8.5); + + if (!mothertrack1.producedByGenerator()) { + continue; + } + hEventSelection.fill(HIST("recMCparticles"), 9.5); + + if (std::abs(mothertrack1.y()) >= selectionConfig.motherRapidityCut) { + continue; + } + hEventSelection.fill(HIST("recMCparticles"), 10.5); + + if (std::abs(mothertrack1.pdgCode()) != o2::constants::physics::kK0Star892) { + continue; + } + hEventSelection.fill(HIST("recMCparticles"), 11.5); + + if (track1PDG == PDG_t::kPiPlus) { + if (!applypTdepPID && !(selectionPID(track1, 0) && selectionPID(track2, 1))) { // pion and kaon + continue; + } else if (applypTdepPID && !(selectionPIDNew(track1, 0) && selectionPIDNew(track2, 1))) { // pion and kaon + continue; + } + hEventSelection.fill(HIST("recMCparticles"), 12.5); + + hEventSelection.fill(HIST("recMCparticles"), 13.5); + + if (cQAplots) { + if (track1.sign() < 0 && track2.sign() > 0) { + hPID.fill(HIST("After/hTPCnsigPi_Neg_mult_pt"), track1.tpcNSigmaPi(), centrality, track1.pt()); + hPID.fill(HIST("After/hTOFnsigPi_Neg_mult_pt"), track1.tofNSigmaPi(), centrality, track1.pt()); + hPID.fill(HIST("After/hTPCnsigKa_Pos_mult_pt"), track2.tpcNSigmaKa(), centrality, track2.pt()); + hPID.fill(HIST("After/hTOFnsigKa_Pos_mult_pt"), track2.tofNSigmaKa(), centrality, track2.pt()); + } else { + hPID.fill(HIST("After/hTPCnsigPi_Pos_mult_pt"), track1.tpcNSigmaPi(), centrality, track1.pt()); + hPID.fill(HIST("After/hTOFnsigPi_Pos_mult_pt"), track1.tofNSigmaPi(), centrality, track1.pt()); + hPID.fill(HIST("After/hTPCnsigKa_Neg_mult_pt"), track2.tpcNSigmaKa(), centrality, track2.pt()); + hPID.fill(HIST("After/hTOFnsigKa_Neg_mult_pt"), track2.tofNSigmaKa(), centrality, track2.pt()); + } + } + + } else if (track1PDG == PDG_t::kKPlus) { + if (!applypTdepPID && !(selectionPID(track1, 1) && selectionPID(track2, 0))) { // kaon and pion + continue; + } else if (applypTdepPID && !(selectionPIDNew(track1, 1) && selectionPIDNew(track2, 0))) { // kaon and pion + continue; + } + hEventSelection.fill(HIST("recMCparticles"), 12.5); + + hEventSelection.fill(HIST("recMCparticles"), 13.5); + + if (cQAplots) { + if (track1.sign() < 0 && track2.sign() > 0) { + hPID.fill(HIST("After/hTPCnsigKa_Neg_mult_pt"), track1.tpcNSigmaKa(), centrality, track1.pt()); + hPID.fill(HIST("After/hTOFnsigKa_Neg_mult_pt"), track1.tofNSigmaKa(), centrality, track1.pt()); + hPID.fill(HIST("After/hTPCnsigPi_Pos_mult_pt"), track2.tpcNSigmaPi(), centrality, track2.pt()); + hPID.fill(HIST("After/hTOFnsigPi_Pos_mult_pt"), track2.tofNSigmaPi(), centrality, track2.pt()); + } else { + hPID.fill(HIST("After/hTPCnsigKa_Pos_mult_pt"), track1.tpcNSigmaKa(), centrality, track1.pt()); + hPID.fill(HIST("After/hTOFnsigKa_Pos_mult_pt"), track1.tofNSigmaKa(), centrality, track1.pt()); + hPID.fill(HIST("After/hTPCnsigPi_Neg_mult_pt"), track2.tpcNSigmaPi(), centrality, track2.pt()); + hPID.fill(HIST("After/hTOFnsigPi_Neg_mult_pt"), track2.tofNSigmaPi(), centrality, track2.pt()); + } + } + } + + /* if (selectionConfig.isApplyCutsOnMother) { + if (mothertrack1.pt() >= selectionConfig.cMaxPtMotherCut) // excluding candidates in overflow + continue; + if ((std::sqrt(mothertrack1.e() * mothertrack1.e() - mothertrack1.p() * mothertrack1.p())) >= selectionConfig.cMaxMinvMotherCut) // excluding candidates in overflow + continue; + } */ + + if (avoidsplitrackMC && oldindex == mothertrack1.globalIndex()) { + hMC.fill(HIST("h1KSRecsplit"), mothertrack1.pt()); + continue; + } + hEventSelection.fill(HIST("recMCparticles"), 15.5); + + 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); + } 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); + } + + mother = daughter1 + daughter2; // Kstar meson + + hMC.fill(HIST("h2KstarRecpt2"), mothertrack1.pt(), centrality, std::sqrt(mothertrack1.e() * mothertrack1.e() - mothertrack1.p() * mothertrack1.p())); + + if (mother.Rapidity() >= selectionConfig.motherRapidityCut) { + continue; + } + + hMC.fill(HIST("h1KstarRecMass"), mother.M()); + hMC.fill(HIST("h2KstarRecpt1"), mother.Pt(), centrality, mother.M()); + } + } + } + } + PROCESS_SWITCH(Kstar892LightIon, processRec, "Process Reconstructed", false); + + void processEvtLossSigLossMC(aod::McCollisions::iterator const& mcCollision, aod::McParticles const& mcParticles, const soa::SmallGroups& recCollisions) + { + auto impactPar = mcCollision.impactParameter(); + hMC.fill(HIST("ImpactCorr/hImpactParameterGen"), impactPar); + + bool isSelectedEvent = false; + auto centrality = -999.; + for (const auto& RecCollision : recCollisions) { + if (!RecCollision.has_mcCollision()) + continue; + if (!selectionEvent(RecCollision, false)) // don't fill event cut histogram + continue; + + if (cSelectMultEstimator == kFT0M) { + centrality = RecCollision.centFT0M(); + } else if (cSelectMultEstimator == kFT0A) { + centrality = RecCollision.centFT0A(); + } else if (cSelectMultEstimator == kFT0C) { + centrality = RecCollision.centFT0C(); + } else if (cSelectMultEstimator == kFV0A) { + centrality = RecCollision.centFV0A(); + } else { + centrality = RecCollision.centFT0M(); // default + } + + isSelectedEvent = true; + } + + if (isSelectedEvent) { + hMC.fill(HIST("ImpactCorr/hImpactParameterRec"), impactPar); + hMC.fill(HIST("ImpactCorr/hImpactParvsCentrRec"), centrality, impactPar); + } + + // Generated MC + for (const auto& mcPart : mcParticles) { + + if (std::abs(mcPart.y()) >= selectionConfig.motherRapidityCut || std::abs(mcPart.pdgCode()) != o2::constants::physics::kK0Star892) + continue; + + // signal loss estimation + hMC.fill(HIST("ImpactCorr/hKstarGenBeforeEvtSel"), mcPart.pt(), impactPar); + if (isSelectedEvent) { + // signal loss estimation + hMC.fill(HIST("ImpactCorr/hKstarGenAfterEvtSel"), mcPart.pt(), impactPar); + } + } // end loop on gen particles + } + PROCESS_SWITCH(Kstar892LightIon, processEvtLossSigLossMC, "Process Signal Loss, Event Loss for Kstar in Light Ion", false); + + using McCollisionMults = soa::Join; + using LabeledTracks = soa::Join; + + void processCorrFactors(McCollisionMults::iterator const& mcCollision, soa::SmallGroups const& collisions, LabeledTracks const& /*particles*/, aod::McParticles const& mcParticles) + { + hMC.fill(HIST("CorrFactors/hGenEvents"), mcCollision.multMCNParticlesEta08(), 0.5); + + if (std::abs(mcCollision.posZ()) > selectionConfig.cutVrtxZ) + return; + + hMC.fill(HIST("CorrFactors/hGenEvents"), mcCollision.multMCNParticlesEta08(), 1.5); + + float centrality = 100.5f; + for (auto const& collision : collisions) { + centrality = collision.centFT0M(); + } + + hMC.fill(HIST("CorrFactors/hCentralityVsMultMC"), centrality, mcCollision.multMCNParticlesEta08()); + hMC.fill(HIST("CorrFactors/hNrecInGen"), collisions.size()); + + for (const auto& mcParticle : mcParticles) { + + if (std::abs(mcParticle.y()) >= selectionConfig.motherRapidityCut) + continue; + + if (std::abs(mcParticle.pdgCode()) == o2::constants::physics::kK0Star892) { + + auto kDaughters = mcParticle.daughters_as(); + if (kDaughters.size() != noOfDaughters) { + continue; + } + + bool hasPos = false; + bool hasNeg = false; + + auto passkaon = false; + auto passpion = false; + for (const auto& kCurrentDaughter : kDaughters) { + // if (!kCurrentDaughter.isPhysicalPrimary()) + // continue; + + int pdgDau = kCurrentDaughter.pdgCode(); + int sign = (pdgDau > 0) - (pdgDau < 0); + + if (sign > 0) + hasPos = true; + if (sign < 0) + hasNeg = true; + + if (std::abs(kCurrentDaughter.pdgCode()) == PDG_t::kKPlus) { + passkaon = true; + daughter1 = ROOT::Math::PxPyPzMVector(kCurrentDaughter.px(), kCurrentDaughter.py(), kCurrentDaughter.pz(), massKa); + + } else if (std::abs(kCurrentDaughter.pdgCode()) == PDG_t::kPiPlus) { + passpion = true; + daughter2 = ROOT::Math::PxPyPzMVector(kCurrentDaughter.px(), kCurrentDaughter.py(), kCurrentDaughter.pz(), massPi); + } + } + + if ((passkaon && passpion) && (hasPos && hasNeg)) { + mother = daughter1 + daughter2; // Kstar meson + + hMC.fill(HIST("CorrFactors/h2dGenKstar"), centrality, mother.Pt()); + hMC.fill(HIST("CorrFactors/h3dGenKstarVsMultMCVsCentrality"), mcCollision.multMCNParticlesEta08(), centrality, mother.Pt()); + } + } + } + + if (collisions.size() == 0) + return; + + hMC.fill(HIST("CorrFactors/hGenEvents"), mcCollision.multMCNParticlesEta08(), 2.5); + } + PROCESS_SWITCH(Kstar892LightIon, processCorrFactors, "Process Signal Loss, Event Loss", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +}