From 3b31d32da14e31e10c2b03b5b6ba5af5b87a1891 Mon Sep 17 00:00:00 2001 From: SuJeong Ji <120470463+SuJeong-Ji@users.noreply.github.com> Date: Tue, 10 Dec 2024 19:05:42 +0900 Subject: [PATCH] [PWGLF] Charged K(892) flow analysis fill histograms (#8877) Co-authored-by: ALICE Action Bot --- PWGLF/Tasks/Resonances/CMakeLists.txt | 5 + PWGLF/Tasks/Resonances/chk892flow.cxx | 975 ++++++++++++++++++++++++++ PWGLF/Utils/collisionCuts.h | 66 +- 3 files changed, 1011 insertions(+), 35 deletions(-) create mode 100644 PWGLF/Tasks/Resonances/chk892flow.cxx diff --git a/PWGLF/Tasks/Resonances/CMakeLists.txt b/PWGLF/Tasks/Resonances/CMakeLists.txt index e6f7b573f2e..d49deced8f3 100644 --- a/PWGLF/Tasks/Resonances/CMakeLists.txt +++ b/PWGLF/Tasks/Resonances/CMakeLists.txt @@ -163,3 +163,8 @@ o2physics_add_dpl_workflow(highmasslambdasvx SOURCES highmasslambdasvx.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(chk892flow + SOURCES chk892flow.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) diff --git a/PWGLF/Tasks/Resonances/chk892flow.cxx b/PWGLF/Tasks/Resonances/chk892flow.cxx new file mode 100644 index 00000000000..7bc523eb9e0 --- /dev/null +++ b/PWGLF/Tasks/Resonances/chk892flow.cxx @@ -0,0 +1,975 @@ +// 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 chk892flow.cxx +/// \brief Reconstruction of track-track decay resonance candidates +/// +/// +/// \author Su-Jeong Ji + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include // FIXME +#include // FIXME + +#include +#include +#include +#include +#include +#include + +#include "TRandom3.h" +#include "TF1.h" +#include "TVector2.h" +#include "Math/Vector3D.h" +#include "Math/Vector4D.h" +#include "Math/GenVector/Boost.h" +#include + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/StepTHn.h" +#include "Framework/O2DatabasePDGPlugin.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/StaticFor.h" +#include "DCAFitter/DCAFitterN.h" + +#include "Common/DataModel/PIDResponse.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Qvectors.h" + +#include "Common/Core/trackUtilities.h" +#include "Common/Core/TrackSelection.h" +#include "Common/Core/RecoDecay.h" + +#include "CommonConstants/PhysicsConstants.h" + +#include "ReconstructionDataFormats/Track.h" + +#include "DataFormatsParameters/GRPObject.h" +#include "DataFormatsParameters/GRPMagField.h" + +#include "CCDB/CcdbApi.h" +#include "CCDB/BasicCCDBManager.h" + +#include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGLF/Utils/collisionCuts.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::soa; +using namespace o2::constants::physics; + +struct chk892flow { + enum binType : unsigned int { + kKstarP = 0, + kKstarN, + kKstarP_Mix, + kKstarN_Mix, + kKstarP_GenINEL10, + kKstarN_GenINEL10, + kKstarP_GenINELgt10, + kKstarN_GenINELgt10, + kKstarP_GenTrig10, + kKstarN_GenTrig10, + kKstarP_GenEvtSel, + kKstarN_GenEvtSel, + kKstarP_Rec, + kKstarN_Rec, + kTYEnd + }; + + SliceCache cache; + Preslice perCollision = aod::track::collisionId; + + using EventCandidates = soa::Join; + using TrackCandidates = soa::Join; + using V0Candidates = aod::V0Datas; + + using MCEventCandidates = soa::Join; + using MCTrackCandidates = soa::Join; + using MCV0Candidates = soa::Join; + + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + Service ccdb; + Service pdg; + o2::ccdb::CcdbApi ccdbApi; + + Configurable cfgURL{"cfgURL", "http://alice-ccdb.cern.ch", "Address of the CCDB to browse"}; + Configurable nolaterthan{"ccdb-no-later-than", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "Latest acceptable timestamp of creation for the object"}; + + // Configurables + ConfigurableAxis cfgBinsPt{"cfgBinsPt", {VARIABLE_WIDTH, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5.0, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6.0, 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6.9, 7.0, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8.0, 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7, 8.8, 8.9, 9.0, 9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8, 9.9, 10.0, 10.1, 10.2, 10.3, 10.4, 10.5, 10.6, 10.7, 10.8, 10.9, 11.0, 11.1, 11.2, 11.3, 11.4, 11.5, 11.6, 11.7, 11.8, 11.9, 12.0, 12.1, 12.2, 12.3, 12.4, 12.5, 12.6, 12.7, 12.8, 12.9, 13.0, 13.1, 13.2, 13.3, 13.4, 13.5, 13.6, 13.7, 13.8, 13.9, 14.0, 14.1, 14.2, 14.3, 14.4, 14.5, 14.6, 14.7, 14.8, 14.9, 15.0}, "Binning of the pT axis"}; + ConfigurableAxis cfgBinsPtQA{"cfgBinsPtQA", {VARIABLE_WIDTH, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.2, 4.4, 4.6, 4.8, 5.0, 5.2, 5.4, 5.6, 5.8, 6.0, 6.2, 6.4, 6.6, 6.8, 7.0, 7.2, 7.4, 7.6, 7.8, 8.0, 8.2, 8.4, 8.6, 8.8, 9.0, 9.2, 9.4, 9.6, 9.8, 10.0}, "Binning of the pT axis"}; + ConfigurableAxis cfgBinsCent{"cfgBinsCent", {VARIABLE_WIDTH, 0.0, 1.0, 5.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 110.0}, "Binning of the centrality axis"}; + ConfigurableAxis cfgBinsVtxZ{"cfgBinsVtxZ", {VARIABLE_WIDTH, -10.0, -9.0, -8.0, -7.0, -6.0, -5.0, -4.0, -3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0}, "Binning of the z-vertex axis"}; + Configurable cNbinsDiv{"cNbinsDiv", 1, "Integer to divide the number of bins"}; + + /// Event cuts + o2::analysis::CollisonCuts colCuts; + Configurable ConfEvtZvtx{"ConfEvtZvtx", 10.f, "Evt sel: Max. z-Vertex (cm)"}; + Configurable ConfEvtOccupancyInTimeRangeMax{"ConfEvtOccupancyInTimeRangeMax", -1, "Evt sel: maximum track occupancy"}; + Configurable ConfEvtOccupancyInTimeRangeMin{"ConfEvtOccupancyInTimeRangeMin", -1, "Evt sel: minimum track occupancy"}; + Configurable ConfEvtTriggerCheck{"ConfEvtTriggerCheck", false, "Evt sel: check for trigger"}; + Configurable ConfEvtTriggerSel{"ConfEvtTriggerSel", 8, "Evt sel: trigger"}; + Configurable ConfEvtOfflineCheck{"ConfEvtOfflineCheck", true, "Evt sel: check for offline selection"}; + Configurable ConfEvtTriggerTVXSel{"ConfEvtTriggerTVXSel", false, "Evt sel: triggerTVX selection (MB)"}; + Configurable ConfEvtTFBorderCut{"ConfEvtTFBorderCut", false, "Evt sel: apply TF border cut"}; + Configurable ConfEvtUseITSTPCvertex{"ConfEvtUseITSTPCvertex", false, "Evt sel: use at lease on ITS-TPC track for vertexing"}; + Configurable ConfEvtZvertexTimedifference{"ConfEvtZvertexTimedifference", true, "Evt sel: apply Z-vertex time difference"}; + Configurable ConfEvtPileupRejection{"ConfEvtPileupRejection", true, "Evt sel: apply pileup rejection"}; + Configurable ConfEvtNoITSROBorderCut{"ConfEvtNoITSROBorderCut", false, "Evt sel: apply NoITSRO border cut"}; + Configurable ConfincludeCentralityMC{"ConfincludeCentralityMC", false, "Include centrality in MC"}; + Configurable ConfEvtCollInTimeRangeStandard{"ConfEvtCollInTimeRangeStandard", true, "Evt sel: apply NoCollInTimeRangeStandard"}; + + /// Track selections + Configurable cMinPtcut{"cMinPtcut", 0.15, "Track minium pt cut"}; + Configurable cMaxEtacut{"cMaxEtacut", 0.8, "Track maximum eta cut"}; + + // Cuts from polarization analysis + Configurable cfgQvecSel{"cfgQvecSel", true, "Reject events when no QVector"}; + Configurable cfgOccupancySel{"cfgOccupancySel", false, "Occupancy selection"}; + Configurable cfgMaxOccupancy{"cfgMaxOccupancy", 999999, "maximum occupancy of tracks in neighbouring collisions in a given time range"}; + Configurable cfgMinOccupancy{"cfgMinOccupancy", -100, "maximum occupancy of tracks in neighbouring collisions in a given time range"}; + Configurable cfgNCollinTR{"cfgNCollinTR", false, "Additional selection for the number of coll in time range"}; + + Configurable cfgCentEst{"cfgCentEst", 1, "Centrality estimator, 1: FT0C, 2: FT0M"}; + + // DCAr to PV + Configurable cMaxbDCArToPVcut{"cMaxbDCArToPVcut", 0.1, "Track DCAr cut to PV Maximum"}; + // DCAz to PV + Configurable cMaxbDCAzToPVcut{"cMaxbDCAzToPVcut", 0.1, "Track DCAz cut to PV Maximum"}; + + /// PID Selections, pion + Configurable cMaxTPCnSigmaPion{"cMaxTPCnSigmaPion", 3.0, "TPC nSigma cut for Pion"}; // TPC + Configurable cMaxTOFnSigmaPion{"cMaxTOFnSigmaPion", 3.0, "TOF nSigma cut for Pion"}; // TOF + Configurable nsigmaCutCombinedPion{"nsigmaCutCombinedPion", -999, "Combined nSigma cut for Pion"}; // Combined + Configurable cTOFVeto{"cTOFVeto", true, "TOF Veto, if false, TOF is nessessary for PID selection"}; // TOF Veto + + // Track selections + Configurable cfgPrimaryTrack{"cfgPrimaryTrack", true, "Primary track selection"}; // kGoldenChi2 | kDCAxy | kDCAz + Configurable cfgGlobalWoDCATrack{"cfgGlobalWoDCATrack", true, "Global track selection without DCA"}; // kQualityTracks (kTrackType | kTPCNCls | kTPCCrossedRows | kTPCCrossedRowsOverNCls | kTPCChi2NDF | kTPCRefit | kITSNCls | kITSChi2NDF | kITSRefit | kITSHits) | kInAcceptanceTracks (kPtRange | kEtaRange) + Configurable cfgGlobalTrack{"cfgGlobalTrack", false, "Global track selection"}; // kGoldenChi2 | kDCAxy | kDCAz + Configurable cfgPVContributor{"cfgPVContributor", false, "PV contributor track selection"}; // PV Contriuibutor + + Configurable cfgITScluster{"cfgITScluster", 0, "Number of ITS cluster"}; + Configurable cfgTPCcluster{"cfgTPCcluster", 0, "Number of TPC cluster"}; + Configurable cfgRatioTPCRowsOverFindableCls{"cfgRatioTPCRowsOverFindableCls", 0.0f, "TPC Crossed Rows to Findable Clusters"}; + Configurable cfgITSChi2NCl{"cfgITSChi2NCl", 999.0, "ITS Chi2/NCl"}; + Configurable cfgTPCChi2NCl{"cfgTPCChi2NCl", 999.0, "TPC Chi2/NCl"}; + Configurable cfgUseTPCRefit{"cfgUseTPCRefit", false, "Require TPC Refit"}; + Configurable cfgUseITSRefit{"cfgUseITSRefit", false, "Require ITS Refit"}; + Configurable cfgHasITS{"cfgHasITS", false, "Require ITS"}; + Configurable cfgHasTPC{"cfgHasTPC", false, "Require TPC"}; + Configurable cfgHasTOF{"cfgHasTOF", false, "Require TOF"}; + + // Secondary Selection + Configurable cfgReturnFlag{"boolReturnFlag", false, "Return Flag for debugging"}; + Configurable cSecondaryArmenterosCut{"boolArmenterosCut", true, "cut on Armenteros-Podolanski graph"}; + + Configurable cfgByPassDauPIDSelection{"cfgByPassDauPIDSelection", true, "Bypass Daughters PID selection"}; + Configurable cSecondaryDauDCAMax{"cSecondaryDauDCAMax", 1., "Maximum DCA Secondary daughters to PV"}; + Configurable cSecondaryDauPosDCAtoPVMin{"cSecondaryDauPosDCAtoPVMin", 0.0, "Minimum DCA Secondary positive daughters to PV"}; + Configurable cSecondaryDauNegDCAtoPVMin{"cSecondaryDauNegDCAtoPVMin", 0.0, "Minimum DCA Secondary negative daughters to PV"}; + + Configurable cSecondaryPtMin{"cSecondaryPtMin", 0.f, "Minimum transverse momentum of Secondary"}; + Configurable cSecondaryRapidityMax{"cSecondaryRapidityMax", 0.5, "Maximum rapidity of Secondary"}; + Configurable cSecondaryRadiusMin{"cSecondaryRadiusMin", 1.2, "Minimum transverse radius of Secondary"}; + Configurable cSecondaryCosPAMin{"cSecondaryCosPAMin", 0.995, "Mininum cosine pointing angle of Secondary"}; + Configurable cSecondaryDCAtoPVMax{"cSecondaryDCAtoPVMax", 0.3, "Maximum DCA Secondary to PV"}; + Configurable cSecondaryProperLifetimeMax{"cSecondaryProperLifetimeMax", 20, "Maximum Secondary Lifetime"}; + Configurable cSecondaryparamArmenterosCut{"paramArmenterosCut", 0.2, "parameter for Armenteros Cut"}; + Configurable cSecondaryMassWindow{"cSecondaryMassWindow", 0.075, "Secondary inv mass selciton window"}; + + // K* selection + Configurable cKstarMaxRap{"cKstarMaxRap", 0.5, "Kstar maximum rapidity"}; + Configurable cKstarMinRap{"cKstarMinRap", -0.5, "Kstar minimum rapidity"}; + + // Confs from flow analysis + Configurable cfgnMods{"cfgnMods", 1, "The number of modulations of interest starting from 2"}; + Configurable cfgNQvec{"cfgNQvec", 7, "The number of total Qvectors for looping over the task"}; + + Configurable cfgQvecDetName{"cfgQvecDetName", "FT0C", "The name of detector to be analyzed"}; + Configurable cfgQvecRefAName{"cfgQvecRefAName", "TPCpos", "The name of detector for reference A"}; + Configurable cfgQvecRefBName{"cfgQvecRefBName", "TPCneg", "The name of detector for reference B"}; + + int DetId; + int RefAId; + int RefBId; + + int QvecDetInd; + int QvecRefAInd; + int QvecRefBInd; + + float centrality; + + // PDG code + int kPDGK0s = 310; + int kPDGK0 = 311; + int kKstarPlus = 323; + int kPiPlus = 211; + + void init(o2::framework::InitContext&) + { + centrality = -999; + + colCuts.setCuts(ConfEvtZvtx, ConfEvtTriggerCheck, ConfEvtTriggerSel, ConfEvtOfflineCheck, /*checkRun3*/ true, /*triggerTVXsel*/ false, ConfEvtOccupancyInTimeRangeMax, ConfEvtOccupancyInTimeRangeMin); + colCuts.init(&histos); + colCuts.setTriggerTVX(ConfEvtTriggerTVXSel); + colCuts.setApplyTFBorderCut(ConfEvtTFBorderCut); + colCuts.setApplyITSTPCvertex(ConfEvtUseITSTPCvertex); + colCuts.setApplyZvertexTimedifference(ConfEvtZvertexTimedifference); + colCuts.setApplyPileupRejection(ConfEvtPileupRejection); + colCuts.setApplyNoITSROBorderCut(ConfEvtNoITSROBorderCut); + colCuts.setApplyCollInTimeRangeStandard(ConfEvtCollInTimeRangeStandard); + + AxisSpec centAxis = {cfgBinsCent, "T0M (%)"}; + AxisSpec vtxzAxis = {cfgBinsVtxZ, "Z Vertex (cm)"}; + AxisSpec epAxis = {100, -1.0 * constants::math::PI, constants::math::PI}; + AxisSpec epresAxis = {100, -1.02, 1.02}; + AxisSpec ptAxis = {cfgBinsPt, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec ptAxisQA = {cfgBinsPtQA, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec v2Axis = {200, -1, 1, "#v_{2}"}; + AxisSpec radiusAxis = {50, 0, 5, "Radius (cm)"}; + AxisSpec cpaAxis = {50, 0.95, 1.0, "CPA"}; + AxisSpec tauAxis = {250, 0, 25, "Lifetime (cm)"}; + AxisSpec dcaAxis = {200, 0, 2, "DCA (cm)"}; + AxisSpec dcaxyAxis = {200, 0, 2, "DCA_{#it{xy}} (cm)"}; + AxisSpec dcazAxis = {200, 0, 2, "DCA_{#it{z}} (cm)"}; + AxisSpec yAxis = {100, -1, 1, "Rapidity"}; + AxisSpec invMassAxisK0s = {400 / cNbinsDiv, 0.3, 0.7, "Invariant Mass (GeV/#it{c}^2)"}; // K0s ~497.611 + AxisSpec invMassAxisReso = {900 / cNbinsDiv, 0.5f, 1.4f, "Invariant Mass (GeV/#it{c}^2)"}; // chK(892) ~892 + AxisSpec invMassAxisScan = {150, 0, 1.5, "Invariant Mass (GeV/#it{c}^2)"}; // For selection + AxisSpec pidQAAxis = {130, -6.5, 6.5}; + AxisSpec dataTypeAxis = {9, 0, 9, "Histogram types"}; + AxisSpec mcTypeAxis = {4, 0, 4, "Histogram types"}; + + // THnSparse + AxisSpec axisType = {binType::kTYEnd, 0, binType::kTYEnd, "Type of bin with charge and mix"}; + AxisSpec mcLabelAxis = {5, -0.5, 4.5, "MC Label"}; + + histos.add("QA/K0sCutCheck", "Check K0s cut", HistType::kTH1D, {AxisSpec{12, -0.5, 11.5, "Check"}}); + + histos.add("QA/before/CentDist", "Centrality distribution", {HistType::kTH1D, {centAxis}}); + 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}}); + + // EventPlane + histos.add("QA/EP/EPhEPDet", "Event plane distribution of FT0C (Det = A)", {HistType::kTH2D, {centAxis, epAxis}}); + histos.add("QA/EP/EPhEPB", "Event plane distribution of TPCpos (B)", {HistType::kTH2D, {centAxis, epAxis}}); + histos.add("QA/EP/EPhEPC", "Event plane distribution of TPCneg (C)", {HistType::kTH2D, {centAxis, epAxis}}); + + histos.add("QA/EP/hEPDet", "Event plane distribution of FT0C (Det = A)", {HistType::kTH2D, {centAxis, epAxis}}); + histos.add("QA/EP/hEPB", "Event plane distribution of TPCpos (B)", {HistType::kTH2D, {centAxis, epAxis}}); + histos.add("QA/EP/hEPC", "Event plane distribution of TPCneg (C)", {HistType::kTH2D, {centAxis, epAxis}}); + histos.add("QA/EP/hEPResAB", "cos(n(A-B))", {HistType::kTH2D, {centAxis, epAxis}}); + histos.add("QA/EP/hEPResAC", "cos(n(A-C))", {HistType::kTH2D, {centAxis, epAxis}}); + histos.add("QA/EP/hEPResBC", "cos(n(B-C))", {HistType::kTH2D, {centAxis, epAxis}}); + + // 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, {axisType, centAxis, ptAxis, invMassAxisReso, v2Axis}); + histos.add("hInvmass_Kstar_Mix", "Invariant mass of unlike-sign chK(892) from mixed event", HistType::kTHnSparseD, {axisType, centAxis, ptAxis, invMassAxisReso, v2Axis}); + + // Mass QA (quick check) + histos.add("QA/before/kstarinvmass", "Invariant mass of unlike-sign chK(892)", HistType::kTH1D, {invMassAxisReso}); + histos.add("QA/before/k0sv2vsinvmass", "Invariant mass vs v2 of unlike-sign K0s", HistType::kTH2D, {invMassAxisK0s, v2Axis}); + histos.add("QA/before/kstarv2vsinvmass", "Invariant mass vs v2 of unlike-sign chK(892)", HistType::kTH2D, {invMassAxisReso, v2Axis}); + histos.add("QA/before/kstarinvmass_Mix", "Invariant mass of unlike-sign chK(892) from mixed event", HistType::kTH1D, {invMassAxisReso}); + histos.add("QA/before/kstarv2vsinvmass_Mix", "Invariant mass vs v2 of unlike-sign chK(892) from mixed event", HistType::kTH2D, {invMassAxisReso, v2Axis}); + + 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/k0sv2vsinvmass", "Invariant mass vs v2 of unlike-sign K0s", HistType::kTH2D, {invMassAxisK0s, v2Axis}); + histos.add("QA/after/kstarv2vsinvmass", "Invariant mass vs v2 of unlike-sign chK(892)", HistType::kTH2D, {invMassAxisReso, v2Axis}); + histos.add("QA/after/kstarinvmass_Mix", "Invariant mass of unlike-sign chK(892) from mixed event", HistType::kTH1D, {invMassAxisReso}); + histos.add("QA/after/kstarv2vsinvmass_Mix", "Invariant mass vs v2 of unlike-sign chK(892) from mixed event", HistType::kTH2D, {invMassAxisReso, v2Axis}); + + DetId = GetDetId(cfgQvecDetName); + RefAId = GetDetId(cfgQvecRefAName); + RefBId = GetDetId(cfgQvecRefBName); + + if (DetId == RefAId || DetId == RefBId || RefAId == RefBId) { + LOGF(info, "Wrong detector configuration \n The FT0C will be used to get Q-Vector \n The TPCpos and TPCneg will be used as reference systems"); + // LOGF(info) << "Wrong detector configuration \n The FT0C will be used to get Q-Vector \n The TPCpos and TPCneg will be used as reference systems"; + DetId = 0; + RefAId = 4; + RefBId = 5; + } + + // MC + if (doprocessMC) { + + 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/KstarPairAsym", "Pair asymmetry of chK(892)", HistType::kTH1D, {AxisSpec{100, -1, 1, "Pair asymmetry"}}); + 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/k0sv2vsinvmass", "Invariant mass vs v2 of unlike-sign K0s", HistType::kTH2D, {invMassAxisK0s, v2Axis}); + histos.add("QAMC/kstarv2vsinvmass", "Invariant mass vs v2 of unlike-sign chK(892)", HistType::kTH2D, {invMassAxisReso, v2Axis}); + histos.add("QAMC/kstarinvmass_noKstar", "Invariant mass of unlike-sign no chK(892)", HistType::kTH1D, {invMassAxisReso}); + histos.add("QAMC/kstarv2vsinvmass_noKstar", "Invariant mass vs v2 of unlike-sign no chK(892)", HistType::kTH2D, {invMassAxisReso, v2Axis}); + + histos.add("hInvmass_Kstar_MC", "Invariant mass of unlike chK(892)", HistType::kTHnSparseD, {axisType, centAxis, ptAxis, invMassAxisReso, v2Axis}); + + 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(); + } + + template + float GetCentrality(CollisionType const& collision) + { + if (cfgCentEst == 1) { + return collision.centFT0C(); + } else if (cfgCentEst == 2) { + return collision.centFT0M(); + } else { + return -999; + } + } + + template + int GetDetId(DetNameType const& name) + { + LOGF(info, "GetDetID running"); + if (name.value == "FT0C") { + return 0; + } else if (name.value == "FT0A") { + return 1; + } else if (name.value == "FT0M") { + return 2; + } else if (name.value == "FV0A") { + return 3; + } else if (name.value == "TPCpos") { + return 4; + } else if (name.value == "TPCneg") { + return 5; + } else { + return false; + } + } + + // Track selection + template + bool trackCut(TrackType const& track) + { + // basic track cuts + if (std::abs(track.pt()) < cMinPtcut) + return false; + if (std::abs(track.eta()) > cMaxEtacut) + return false; + if (track.itsNCls() < cfgITScluster) + return false; + if (track.tpcNClsFound() < cfgTPCcluster) + return false; + if (track.tpcCrossedRowsOverFindableCls() < cfgRatioTPCRowsOverFindableCls) + return false; + if (track.itsChi2NCl() >= cfgITSChi2NCl) + return false; + if (track.tpcChi2NCl() >= cfgTPCChi2NCl) + return false; + if (cfgHasITS && !track.hasITS()) + return false; + if (cfgHasTPC && !track.hasTPC()) + return false; + if (cfgHasTOF && !track.hasTOF()) + return false; + if (cfgUseITSRefit && !track.passedITSRefit()) + return false; + if (cfgUseTPCRefit && !track.passedTPCRefit()) + return false; + if (cfgPVContributor && !track.isPVContributor()) + return false; + if (cfgGlobalWoDCATrack && !track.isGlobalTrackWoDCA()) + return false; + if (cfgGlobalTrack && !track.isGlobalTrack()) + return false; + if (cfgPrimaryTrack && !track.isPrimaryTrack()) + return false; + if (std::abs(track.dcaXY()) > cMaxbDCArToPVcut) + return false; + if (std::abs(track.dcaZ()) > cMaxbDCAzToPVcut) + return false; + return true; + } + + // PID selection tools + template + bool selectionPIDPion(TrackType const& candidate) + { + bool tpcPIDPassed{false}, tofPIDPassed{false}; + if (std::abs(candidate.tpcNSigmaPi()) < cMaxTPCnSigmaPion) { + tpcPIDPassed = true; + } else { + return false; + } + if (candidate.hasTOF()) { + if (std::abs(candidate.tofNSigmaPi()) < cMaxTOFnSigmaPion) { + tofPIDPassed = true; + } + if ((nsigmaCutCombinedPion > 0) && (candidate.tpcNSigmaPi() * candidate.tpcNSigmaPi() + candidate.tofNSigmaPi() * candidate.tofNSigmaPi() < nsigmaCutCombinedPion * nsigmaCutCombinedPion)) { + tofPIDPassed = true; + } + } else { + if (!cTOFVeto) { + return false; + } + tofPIDPassed = true; + } + if (tpcPIDPassed && tofPIDPassed) { + return true; + } + return false; + } + + template + bool selectionK0s(CollisionType const& collision, K0sType const& candidate) + { + auto DauDCA = candidate.dcaV0daughters(); + auto DauPosDCAtoPV = candidate.dcapostopv(); + auto DauNegDCAtoPV = candidate.dcanegtopv(); + auto pT = candidate.pt(); + auto Rapidity = candidate.yK0Short(); + auto Radius = candidate.v0radius(); + auto DCAtoPV = candidate.dcav0topv(); + auto CPA = candidate.v0cosPA(); + auto PropTauK0s = candidate.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * MassK0Short; + auto mK0s = candidate.mK0Short(); + + if (cfgReturnFlag) { + bool returnFlag = true; + + histos.fill(HIST("QA/K0sCutCheck"), 0); + if (DauDCA > cSecondaryDauDCAMax) { + histos.fill(HIST("QA/K0sCutCheck"), 1); + returnFlag = false; + } + if (DauPosDCAtoPV < cSecondaryDauPosDCAtoPVMin) { + histos.fill(HIST("QA/K0sCutCheck"), 2); + returnFlag = false; + } + if (DauNegDCAtoPV < cSecondaryDauNegDCAtoPVMin) { + histos.fill(HIST("QA/K0sCutCheck"), 3); + returnFlag = false; + } + if (pT < cSecondaryPtMin) { + histos.fill(HIST("QA/K0sCutCheck"), 4); + returnFlag = false; + } + if (Rapidity > cSecondaryRapidityMax) { + histos.fill(HIST("QA/K0sCutCheck"), 5); + returnFlag = false; + } + if (Radius < cSecondaryRadiusMin) { + histos.fill(HIST("QA/K0sCutCheck"), 6); + returnFlag = false; + } + if (DCAtoPV > cSecondaryDCAtoPVMax) { + histos.fill(HIST("QA/K0sCutCheck"), 7); + returnFlag = false; + } + if (CPA < cSecondaryCosPAMin) { + histos.fill(HIST("QA/K0sCutCheck"), 8); + returnFlag = false; + } + if (PropTauK0s > cSecondaryProperLifetimeMax) { + histos.fill(HIST("QA/K0sCutCheck"), 9); + returnFlag = false; + } + if (fabs(mK0s - MassK0Short) > cSecondaryMassWindow) { + histos.fill(HIST("QA/K0sCutCheck"), 10); + returnFlag = false; + } + if (candidate.qtarm() < cSecondaryparamArmenterosCut * TMath::Abs(candidate.alpha())) { + histos.fill(HIST("QA/K0sCutCheck"), 11); + returnFlag = false; + } + return returnFlag; + + } else { + histos.fill(HIST("QA/K0sCutCheck"), 0); + if (DauDCA > cSecondaryDauDCAMax) { + histos.fill(HIST("QA/K0sCutCheck"), 1); + return false; + } + if (DauPosDCAtoPV < cSecondaryDauPosDCAtoPVMin) { + histos.fill(HIST("QA/K0sCutCheck"), 2); + return false; + } + if (DauNegDCAtoPV < cSecondaryDauNegDCAtoPVMin) { + histos.fill(HIST("QA/K0sCutCheck"), 3); + return false; + } + if (pT < cSecondaryPtMin) { + histos.fill(HIST("QA/K0sCutCheck"), 4); + return false; + } + if (Rapidity > cSecondaryRapidityMax) { + histos.fill(HIST("QA/K0sCutCheck"), 5); + return false; + } + if (Radius < cSecondaryRadiusMin) { + histos.fill(HIST("QA/K0sCutCheck"), 6); + return false; + } + if (DCAtoPV > cSecondaryDCAtoPVMax) { + histos.fill(HIST("QA/K0sCutCheck"), 7); + return false; + } + if (CPA < cSecondaryCosPAMin) { + histos.fill(HIST("QA/K0sCutCheck"), 8); + return false; + } + if (PropTauK0s > cSecondaryProperLifetimeMax) { + histos.fill(HIST("QA/K0sCutCheck"), 9); + return false; + } + if (fabs(mK0s - MassK0Short) > cSecondaryMassWindow) { + histos.fill(HIST("QA/K0sCutCheck"), 10); + return false; + } + if (candidate.qtarm() < cSecondaryparamArmenterosCut * TMath::Abs(candidate.alpha())) { + histos.fill(HIST("QA/K0sCutCheck"), 11); + return false; + } + return true; + } + + } // selectionK0s + + double GetPhiInRange(double phi) + { + double result = phi; + while (result < 0) { + result = result + 2. * TMath::Pi() / 2; + } + while (result > 2. * TMath::Pi() / 2) { + result = result - 2. * TMath::Pi() / 2; + } + return result; + } + + template + bool isTrueKstar(const TrackTemplate& bTrack, const V0Template& K0scand) + { + if (abs(bTrack.PDGCode()) != kPiPlus) // Are you pion? + return false; + if (abs(K0scand.PDGCode()) != kPDGK0s) // Are you K0s? + return false; + + auto motherbTrack = bTrack.template mothers_as(); + auto motherkV0 = K0scand.template mothers_as(); + + // Check bTrack first + if (abs(motherbTrack.pdgCode()) != kKstarPlus) // Are you charged Kstar's daughter? + return false; // Apply first since it's more restrictive + + if (abs(motherkV0.pdgCode()) != 310) // Is it K0s? + return false; + // Check if K0s's mother is K0 (311) + auto motherK0 = motherkV0.template mothers_as(); + if (abs(motherK0.pdgCode()) != 311) + return false; + + // Check if K0's mother is Kstar (323) + auto motherKstar = motherK0.template mothers_as(); + if (abs(motherKstar.pdgCode()) != 323) + return false; + + // Check if bTrack and K0 have the same mother (global index) + if (motherbTrack.globalIndex() != motherK0.globalIndex()) + return false; + + return true; + } + + int count = 0; + + template + void fillHistograms(const CollisionType& collision, const TracksType& dTracks1, const TracksTypeK0s& dTracks2, int nmode) + { + histos.fill(HIST("QA/before/CentDist"), centrality); + + QvecDetInd = DetId * 4 + 3 + (nmode - 2) * cfgNQvec * 4; + QvecRefAInd = RefAId * 4 + 3 + (nmode - 2) * cfgNQvec * 4; + QvecRefBInd = RefBId * 4 + 3 + (nmode - 2) * cfgNQvec * 4; + + double EPDet = TMath::ATan2(collision.qvecIm()[QvecDetInd], collision.qvecRe()[QvecDetInd]) / static_cast(nmode); + double EPRefB = TMath::ATan2(collision.qvecIm()[QvecRefAInd], collision.qvecRe()[QvecRefAInd]) / static_cast(nmode); + double EPRefC = TMath::ATan2(collision.qvecIm()[QvecRefBInd], collision.qvecRe()[QvecRefBInd]) / static_cast(nmode); + + double EPResAB = TMath::Cos(static_cast(nmode) * (EPDet - EPRefB)); + double EPResAC = TMath::Cos(static_cast(nmode) * (EPDet - EPRefC)); + double EPResBC = TMath::Cos(static_cast(nmode) * (EPRefB - EPRefC)); + + histos.fill(HIST("QA/EP/hEPDet"), centrality, EPDet); + histos.fill(HIST("QA/EP/hEPB"), centrality, EPRefB); + histos.fill(HIST("QA/EP/hEPC"), centrality, EPRefC); + histos.fill(HIST("QA/EP/hEPResAB"), centrality, EPResAB); + histos.fill(HIST("QA/EP/hEPResAC"), centrality, EPResAC); + histos.fill(HIST("QA/EP/hEPResBC"), centrality, EPResBC); + + TLorentzVector lDecayDaughter1, lDecayDaughter2, lResoSecondary, lDecayDaughter_bach, lResoKstar; + std::vector trackIndicies = {}; + std::vector k0sIndicies = {}; + + for (auto& bTrack : dTracks1) { + auto trkbpt = bTrack.pt(); + auto istrkbhasTOF = bTrack.hasTOF(); + auto trkbNSigmaPiTPC = bTrack.tpcNSigmaPi(); + 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); + } + histos.fill(HIST("QA/before/trkbpionpT"), trkbpt); + histos.fill(HIST("QA/before/trkbpionDCAxy"), bTrack.dcaXY()); + histos.fill(HIST("QA/before/trkbpionDCAz"), bTrack.dcaZ()); + } + + if (!trackCut(bTrack)) + continue; + if (!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); + } + 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()); + } + + for (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 trkkDauDCA = K0scand.dcaV0daughters(); + auto trkkDauDCAPostoPV = K0scand.dcapostopv(); + auto trkkDauDCANegtoPV = K0scand.dcanegtopv(); + auto trkkpt = K0scand.pt(); + auto trkky = K0scand.yK0Short(); + auto trkkRadius = K0scand.v0radius(); + auto trkkDCAtoPV = K0scand.dcav0topv(); + auto trkkCPA = K0scand.v0cosPA(); + auto trkkPropTau = K0scand.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * MassK0Short; + 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); + } + 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/hPropTauSecondary"), trkkPropTau); + histos.fill(HIST("QA/before/hInvmassSecondary"), trkkMass); + } + + // if (!trackCut(posDauTrack) || !trackCut(negDauTrack)) // Too tight cut for K0s daugthers + // continue; + if (!cfgByPassDauPIDSelection && !selectionPIDPion(posDauTrack)) // Perhaps it's already applied in trackCut (need to check QA plots) + continue; + if (!cfgByPassDauPIDSelection && !selectionPIDPion(negDauTrack)) + continue; + if (!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); + } + 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/hPropTauSecondary"), trkkPropTau); + histos.fill(HIST("QA/after/hInvmassSecondary"), trkkMass); + } + k0sIndicies.push_back(K0scand.index()); + } + + for (auto& trackIndex : trackIndicies) { + for (auto& k0sIndex : k0sIndicies) { + auto bTrack = dTracks1.rawIteratorAt(trackIndex); + auto K0scand = dTracks2.rawIteratorAt(k0sIndex); + + lDecayDaughter_bach.SetXYZM(bTrack.px(), bTrack.py(), bTrack.pz(), MassPionCharged); + lResoSecondary.SetXYZM(K0scand.px(), K0scand.py(), K0scand.pz(), MassK0Short); + lResoKstar = lResoSecondary + lDecayDaughter_bach; + + auto phiminuspsi_k0s = GetPhiInRange(lResoSecondary.Phi() - EPDet); + auto phiminuspsi_kstar = GetPhiInRange(lResoKstar.Phi() - EPDet); + auto v2_k0s = TMath::Cos(static_cast(nmode) * phiminuspsi_k0s); + auto v2_kstar = TMath::Cos(static_cast(nmode) * phiminuspsi_kstar); + + // QA plots + if constexpr (!IsMix) { + histos.fill(HIST("QA/before/KstarRapidity"), lResoKstar.Rapidity()); + histos.fill(HIST("QA/before/kstarinvmass"), lResoKstar.M()); + histos.fill(HIST("QA/before/k0sv2vsinvmass"), lResoSecondary.M(), v2_k0s); + histos.fill(HIST("QA/before/kstarv2vsinvmass"), lResoKstar.M(), v2_kstar); + } + + if (lResoKstar.Rapidity() > cKstarMaxRap || lResoKstar.Rapidity() < cKstarMinRap) + continue; + + if constexpr (!IsMix) { + unsigned int typeKstar = bTrack.sign() > 0 ? binType::kKstarP : binType::kKstarN; + + histos.fill(HIST("QA/after/KstarRapidity"), lResoKstar.Rapidity()); + histos.fill(HIST("QA/after/kstarinvmass"), lResoKstar.M()); + histos.fill(HIST("QA/after/k0sv2vsinvmass"), lResoSecondary.M(), v2_k0s); + histos.fill(HIST("QA/after/kstarv2vsinvmass"), lResoKstar.M(), v2_kstar); + histos.fill(HIST("hInvmass_Kstar"), typeKstar, centrality, lResoKstar.Pt(), lResoKstar.M(), v2_kstar); + + } // IsMix + } // K0scand + } // bTrack + + count++; + + } // fillHistograms + + // process data + void processData(EventCandidates::iterator const& collision, + TrackCandidates const& tracks, + V0Candidates const& v0s, + aod::BCsWithTimestamps const&) + { + if (!colCuts.isSelected(collision)) // Default event selection + return; + if (cfgQvecSel && (collision.qvecAmp()[DetId] < 1e-4 || collision.qvecAmp()[RefAId] < 1e-4 || collision.qvecAmp()[RefBId] < 1e-4)) + return; // If we don't have a Q-vector + colCuts.fillQA(collision); + centrality = GetCentrality(collision); + + fillHistograms(collision, tracks, v0s, 2); // second order + } + PROCESS_SWITCH(chk892flow, processData, "Process Event for data without Partitioning", true); + + // process MC reconstructed level + void processMC(EventCandidates::iterator const& collision, + MCTrackCandidates const& tracks, + MCV0Candidates const& v0s) + { + + histos.fill(HIST("QAMC/hEvent"), 1.0); + + fillHistograms(collision, tracks, v0s, 2); + } + PROCESS_SWITCH(chk892flow, processMC, "Process Event for MC", false); +}; +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc, TaskName{"lf-chk892flow"})}; +} diff --git a/PWGLF/Utils/collisionCuts.h b/PWGLF/Utils/collisionCuts.h index 6322e272083..a272cb32670 100644 --- a/PWGLF/Utils/collisionCuts.h +++ b/PWGLF/Utils/collisionCuts.h @@ -24,6 +24,8 @@ #include "Framework/Logger.h" #include "Common/DataModel/EventSelection.h" +using namespace o2::framework; + namespace o2::analysis { @@ -77,26 +79,26 @@ class CollisonCuts /// Initializes histograms for the task /// \param registry Histogram registry to be passed - void init(o2::framework::HistogramRegistry* registry) + void init(HistogramRegistry* registry) { if (!mCutsSet) { LOGF(error, "Event selection not set - quitting!"); } mHistogramRegistry = registry; - mHistogramRegistry->add("Event/posZ", "; vtx_{z} (cm); Entries", o2::framework::kTH1F, {{250, -12.5, 12.5}}); // z-vertex histogram after event selections - mHistogramRegistry->add("Event/posZ_noCut", "; vtx_{z} (cm); Entries", o2::framework::kTH1F, {{250, -12.5, 12.5}}); // z-vertex histogram before all selections + mHistogramRegistry->add("Event/posZ", "; vtx_{z} (cm); Entries", kTH1F, {{250, -12.5, 12.5}}); // z-vertex histogram after event selections + mHistogramRegistry->add("Event/posZ_noCut", "; vtx_{z} (cm); Entries", kTH1F, {{250, -12.5, 12.5}}); // z-vertex histogram before all selections if (mCheckIsRun3) { - mHistogramRegistry->add("Event/CentFV0A", "; vCentV0A; Entries", o2::framework::kTH1F, {{110, 0, 110}}); - mHistogramRegistry->add("Event/CentFT0M", "; vCentT0M; Entries", o2::framework::kTH1F, {{110, 0, 110}}); - mHistogramRegistry->add("Event/CentFT0C", "; vCentT0C; Entries", o2::framework::kTH1F, {{110, 0, 110}}); - mHistogramRegistry->add("Event/CentFT0A", "; vCentT0A; Entries", o2::framework::kTH1F, {{110, 0, 110}}); - mHistogramRegistry->add("Event/posZ_ITSOnly", "; vtx_{z} (cm); Entries", o2::framework::kTH1F, {{250, -12.5, 12.5}}); - mHistogramRegistry->add("Event/posZ_ITSTPC", "; vtx_{z} (cm); Entries", o2::framework::kTH1F, {{250, -12.5, 12.5}}); - mHistogramRegistry->add("Event/trackOccupancyInTimeRange_noCut", "; Occupancy; Entries", o2::framework::kTH1F, {{500, 0., 20000.}}); + mHistogramRegistry->add("Event/CentFV0A", "; vCentV0A; Entries", kTH1F, {{110, 0, 110}}); + mHistogramRegistry->add("Event/CentFT0M", "; vCentT0M; Entries", kTH1F, {{110, 0, 110}}); + mHistogramRegistry->add("Event/CentFT0C", "; vCentT0C; Entries", kTH1F, {{110, 0, 110}}); + mHistogramRegistry->add("Event/CentFT0A", "; vCentT0A; Entries", kTH1F, {{110, 0, 110}}); + mHistogramRegistry->add("Event/posZ_ITSOnly", "; vtx_{z} (cm); Entries", kTH1F, {{250, -12.5, 12.5}}); + mHistogramRegistry->add("Event/posZ_ITSTPC", "; vtx_{z} (cm); Entries", kTH1F, {{250, -12.5, 12.5}}); + mHistogramRegistry->add("Event/trackOccupancyInTimeRange_noCut", "; Occupancy; Entries", kTH1F, {{500, 0., 20000.}}); } else { - mHistogramRegistry->add("Event/CentRun2V0M", "; vCentV0M; Entries", o2::framework::kTH1F, {{110, 0, 110}}); + mHistogramRegistry->add("Event/CentRun2V0M", "; vCentV0M; Entries", kTH1F, {{110, 0, 110}}); } - mHistogramRegistry->add("CollCutCounts", "; ; Entries", o2::framework::kTH1F, {{11, 0., 11.}}); + mHistogramRegistry->add("CollCutCounts", "; ; Entries", kTH1F, {{11, 0., 11.}}); mHistogramRegistry->get(HIST("CollCutCounts"))->GetXaxis()->SetBinLabel(binLabel(EvtSel::kAllEvent), "all"); mHistogramRegistry->get(HIST("CollCutCounts"))->GetXaxis()->SetBinLabel(binLabel(EvtSel::kFlagZvertex), "Zvtx"); mHistogramRegistry->get(HIST("CollCutCounts"))->GetXaxis()->SetBinLabel(binLabel(EvtSel::kFlagTrigerTVX), "IsTriggerTVX"); @@ -219,11 +221,6 @@ class CollisonCuts LOGF(debug, "Offline selection failed (sel7)"); return false; } - auto bc = col.template bc_as(); - if (!(bc.eventCuts() & BIT(aod::Run2EventCuts::kAliEventCutsAccepted))) { - LOGF(debug, "Offline selection failed (AliEventCuts)"); - return false; - } mHistogramRegistry->fill(HIST("CollCutCounts"), EvtSel::kFlagTrigerTVX); } if (mCheckTrigger && !col.alias_bit(mTrigger)) { @@ -276,24 +273,23 @@ class CollisonCuts } private: - using BCsWithRun2Info = soa::Join; - o2::framework::HistogramRegistry* mHistogramRegistry = nullptr; ///< For QA output - bool mCutsSet = false; ///< Protection against running without cuts - bool mCheckTrigger = false; ///< Check for trigger - bool mTriggerTVXselection = false; ///< Check for trigger TVX selection - bool mCheckOffline = false; ///< Check for offline criteria (might change) - bool mCheckIsRun3 = false; ///< Check if running on Pilot Beam - bool mInitialTriggerScan = false; ///< Check trigger when the event is first selected - bool mApplyTFBorderCut = false; ///< Apply time frame border cut - bool mApplyITSTPCvertex = false; ///< Apply at least one ITS-TPC track for vertexing - bool mApplyZvertexTimedifference = false; ///< removes collisions with large differences between z of PV by tracks and z of PV from FT0 A-C time difference. - bool mApplyPileupRejection = false; ///< Pileup rejection - bool mApplyNoITSROBorderCut = false; ///< Apply NoITSRO frame border cut - bool mApplyCollInTimeRangeStandard = false; ///< Apply NoCollInTimeRangeStandard selection - int mTrigger = kINT7; ///< Trigger to check for - float mZvtxMax = 999.f; ///< Maximal deviation from nominal z-vertex (cm) - int mtrackOccupancyInTimeRangeMax = -1; ///< Maximum trackOccupancyInTimeRange cut (-1 no cut) - int mtrackOccupancyInTimeRangeMin = -1; ///< Minimum trackOccupancyInTimeRange cut (-1 no cut) + HistogramRegistry* mHistogramRegistry = nullptr; ///< For QA output + bool mCutsSet = false; ///< Protection against running without cuts + bool mCheckTrigger = false; ///< Check for trigger + bool mTriggerTVXselection = false; ///< Check for trigger TVX selection + bool mCheckOffline = false; ///< Check for offline criteria (might change) + bool mCheckIsRun3 = false; ///< Check if running on Pilot Beam + bool mInitialTriggerScan = false; ///< Check trigger when the event is first selected + bool mApplyTFBorderCut = false; ///< Apply time frame border cut + bool mApplyITSTPCvertex = false; ///< Apply at least one ITS-TPC track for vertexing + bool mApplyZvertexTimedifference = false; ///< removes collisions with large differences between z of PV by tracks and z of PV from FT0 A-C time difference. + bool mApplyPileupRejection = false; ///< Pileup rejection + bool mApplyNoITSROBorderCut = false; ///< Apply NoITSRO frame border cut + bool mApplyCollInTimeRangeStandard = false; ///< Apply NoCollInTimeRangeStandard selection + int mTrigger = kINT7; ///< Trigger to check for + float mZvtxMax = 999.f; ///< Maximal deviation from nominal z-vertex (cm) + int mtrackOccupancyInTimeRangeMax = -1; ///< Maximum trackOccupancyInTimeRange cut (-1 no cut) + int mtrackOccupancyInTimeRangeMin = -1; ///< Minimum trackOccupancyInTimeRange cut (-1 no cut) }; } // namespace o2::analysis