diff --git a/.gitignore b/.gitignore deleted file mode 100644 index 73ebcda4683..00000000000 --- a/.gitignore +++ /dev/null @@ -1,31 +0,0 @@ -CMakeCache.txt -CMakeFiles/ -.vscode -build -compile_commands.json -.DS_Store -megalinter-reports/ -.mypy_cache/ -github_conf/ -.cache - -#VScode -.vscode/* -!.vscode/settings.json -!.vscode/tasks.json -!.vscode/launch.json -!.vscode/extensions.json -!.vscode/*.code-snippets - -# Local History for Visual Studio Code -.history/ - -# Built Visual Studio Code Extensions -*.vsix - -#Xcode -## User settings -xcuserdata/ - -#QtCreator -CMakeLists.txt.user diff --git a/PWGJE/.DS_Store b/PWGJE/.DS_Store new file mode 100644 index 00000000000..ea48dab3283 Binary files /dev/null and b/PWGJE/.DS_Store differ diff --git a/PWGJE/Tasks/CMakeLists.txt b/PWGJE/Tasks/CMakeLists.txt index 050decc2ba7..5cb34d756f7 100644 --- a/PWGJE/Tasks/CMakeLists.txt +++ b/PWGJE/Tasks/CMakeLists.txt @@ -132,6 +132,10 @@ if(FastJet_FOUND) SOURCES phiInJets.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(kstar-in-jets + SOURCES kstarInJets.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore + COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(jet-taggerhf-qa SOURCES jettaggerhfQA.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore diff --git a/PWGJE/Tasks/kstarInJets.cxx b/PWGJE/Tasks/kstarInJets.cxx new file mode 100644 index 00000000000..d5c8cb1246f --- /dev/null +++ b/PWGJE/Tasks/kstarInJets.cxx @@ -0,0 +1,1160 @@ +// 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 kstarInJets.cxx +/// \brief Reconstruction of Phi yield through track-track Minv correlations for resonance hadrochemistry analysis. +/// +/// +/// \author Jimun Lee +/// \author Adrian Fereydon Nassirpour + +#include +#include + +#include "Framework/ASoA.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/Track.h" + +#include "Common/Core/RecoDecay.h" +#include "Common/Core/TrackSelection.h" +#include "Common/Core/TrackSelectionDefaults.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponse.h" +#include "CommonConstants/PhysicsConstants.h" + +#include "PWGJE/Core/FastJetUtilities.h" +#include "PWGJE/Core/JetDerivedDataUtilities.h" +#include "PWGJE/DataModel/Jet.h" + +#include "PWGLF/DataModel/LFResonanceTables.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +struct kstarInJets { + SliceCache cache; + HistogramRegistry JEhistos{"JEhistos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + HistogramRegistry registry{"registry", + {{"h_jet_pt", "jet pT;#it{p}_{T,jet} (GeV/#it{c});entries", {HistType::kTH1F, {{4000, 0., 200.}}}}, + {"h_jet_eta", "jet #eta;#eta_{jet};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}}, + {"h_jet_phi", "jet #phi;#phi_{jet};entries", {HistType::kTH1F, {{80, -1.0, 7.}}}}, + {"h_matched_REC_jet_pt", "matched_REC level jet pT;#it{p}_{T,jet part} (GeV/#it{c});Delta", {HistType::kTH2F, {{200, 0., 200.}, {400, -20., 20.}}}}, + {"h_matched_REC_jet_eta", "matched_REC level jet #eta;#eta_{jet part};Delta", {HistType::kTH2F, {{100, -1.0, 1.0}, {400, -20., 20.}}}}, + {"h_matched_REC_jet_phi", "matched_REC level jet #phi;#phi_{jet part};Delta", {HistType::kTH2F, {{80, -1.0, 7.}, {400, -20., 20.}}}}, + {"h_matched_GEN_jet_pt", "matched_GEN level jet pT;#it{p}_{T,jet part} (GeV/#it{c});Delta", {HistType::kTH2F, {{200, 0., 200.}, {400, -20., 20.}}}}, + {"h_matched_GEN_jet_eta", "matched_GEN level jet #eta;#eta_{jet part};Delta", {HistType::kTH2F, {{100, -1.0, 1.0}, {400, -20., 20.}}}}, + {"h_matched_GEN_jet_phi", "matched_GEN level jet #phi;#phi_{jet part};Delta", {HistType::kTH2F, {{80, -1.0, 7.}, {400, -20., 20.}}}}, + {"h_part_jet_pt", "particle level jet pT;#it{p}_{T,jet part} (GeV/#it{c});entries", {HistType::kTH1F, {{4000, 0., 200.}}}}, + {"h_part_jet_eta", "particle level jet #eta;#eta_{jet part};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}}, + {"h_part_jet_phi", "particle level jet #phi;#phi_{jet part};entries", {HistType::kTH1F, {{80, -1.0, 7.}}}}}}; + + Configurable cfgeventSelections{"cfgeventSelections", "sel8", "choose event selection"}; + Configurable cfgtrackSelections{"cfgtrackSelections", "globalTracks", "set track selections"}; + + Configurable cfgtrkMinPt{"cfgtrkMinPt", 0.15, "set track min pT"}; + Configurable cfgtrkMaxEta{"cfgtrkMaxEta", 0.9, "set track max Eta"}; + Configurable cfgMaxDCArToPVcut{"cfgMaxDCArToPVcut", 0.5, "Track DCAr cut to PV Maximum"}; + Configurable cfgMaxDCAzToPVcut{"cfgMaxDCAzToPVcut", 2.0, "Track DCAz cut to PV Maximum"}; + Configurable cfgPrimaryTrack{"cfgPrimaryTrack", true, "Primary track selection"}; // kGoldenChi2 | kDCAxy | kDCAz + Configurable cfgConnectedToPV{"cfgConnectedToPV", true, "PV contributor track selection"}; // PV Contriuibutor + Configurable cfgGlobalWoDCATrack{"cfgGlobalWoDCATrack", true, "Global track selection without DCA"}; // kQualityTracks (kTrackType | kTPCNCls | kTPCCrossedRows | kTPCCrossedRowsOverNCls | kTPCChi2NDF | kTPCRefit | kITSNCls | kITSChi2NDF | kITSRefit | kITSHits) | kInAcceptanceTracks (kPtRange | kEtaRange) + Configurable cfgnFindableTPCClusters{"cfgnFindableTPCClusters", 50, "nFindable TPC Clusters"}; + Configurable cfgnTPCCrossedRows{"cfgnTPCCrossedRows", 70, "nCrossed TPC Rows"}; + Configurable cfgnRowsOverFindable{"cfgnRowsOverFindable", 1.2, "nRowsOverFindable TPC CLusters"}; + Configurable cfgnTPCChi2{"cfgnTPChi2", 4.0, "nTPC Chi2 per Cluster"}; + Configurable cfgnITSChi2{"cfgnITShi2", 36.0, "nITS Chi2 per Cluster"}; + Configurable cfgnTPCPID{"cfgnTPCPID", 4, "nTPC PID"}; + Configurable cfgnTOFPID{"cfgnTOFPID", 4, "nTOF PID"}; + Configurable cfgjetPtMin{"cfgjetPtMin", 5.0, "minimum jet pT cut"}; + Configurable cfgjetR{"cfgjetR", 0.4, "jet resolution parameter"}; + Configurable cfgVtxCut{"cfgVtxCut", 10.0, "V_z cut selection"}; + Configurable cDebugLevel{"cDebugLevel", 0, "Resolution of Debug"}; + Configurable cfgBR{"cfgBR", false, "Forces Gen. Charged BR Only"}; + Configurable cfgSimPID{"cfgSimPID", false, "Enforces PID on the Gen. Rec level"}; + // CONFIG DONE + ///////////////////////////////////////// //INIT + + int eventSelection = -1; + + void init(o2::framework::InitContext&) + { + // HISTOGRAMS + const AxisSpec axisEta{30, -1.5, +1.5, "#eta"}; + const AxisSpec axisPhi{200, -1, +7, "#phi"}; + const AxisSpec axisPt{200, 0, +200, "#pt"}; + const AxisSpec MinvAxis = {400, 0.65, 1.15}; + const AxisSpec PtAxis = {200, 0, 20.0}; + const AxisSpec MultAxis = {100, 0, 100}; + const AxisSpec dRAxis = {100, 0, 3.6}; + + JEhistos.add("ptJEHistogramPion", "ptJEHistogramPion", kTH1F, {PtAxis}); + JEhistos.add("ptJEHistogramKaon", "ptJEHistogramKaon", kTH1F, {PtAxis}); + JEhistos.add("ptJEHistogramProton", "ptJEHistogramProton", kTH1F, {PtAxis}); + JEhistos.add("ptJEHistogramPhi", "ptJEHistogramPhi", kTH1F, {PtAxis}); + JEhistos.add("ptJEHistogramPhi_JetTrigger", "ptJEHistogramPhi_JetTrigger", kTH1F, {PtAxis}); + JEhistos.add("minvJEHistogramPhi", "minvJEHistogramPhi", kTH1F, {MinvAxis}); + + JEhistos.add("Resp_Matrix", "Resp_Matrix", HistType::kTHnSparseD, {PtAxis, axisPt, PtAxis, axisPt}); // REC(Phi,Jet), GEN(Phi,Jet) + JEhistos.add("Resp_Matrix_MATCHED", "Resp_Matrix_MATCHED", HistType::kTHnSparseD, {PtAxis, axisPt, PtAxis, axisPt}); // REC(Phi,Jet), GEN(Phi,Jet) + + JEhistos.add("ptGeneratedPion", "ptGeneratedPion", kTH1F, {PtAxis}); + JEhistos.add("ptGeneratedKaon", "ptGeneratedKaon", kTH1F, {PtAxis}); + JEhistos.add("ptGeneratedProton", "ptGeneratedProton", kTH1F, {PtAxis}); + JEhistos.add("ptGeneratedPhi", "ptGeneratedPhi", kTH1F, {PtAxis}); + JEhistos.add("ptGeneratedPhi_ALLBR", "ptGeneratedPhi_ALLBR", kTH1F, {PtAxis}); + + JEhistos.add("ptGeneratedPhi_JetTrigger", "ptGeneratedPhi_JetTrigger", kTH1F, {PtAxis}); + + JEhistos.add("mGeneratedPhi", "mGeneratedPhi", kTH1F, {MinvAxis}); + + JEhistos.add("etaHistogram", "etaHistogram", kTH1F, {axisEta}); + JEhistos.add("phiHistogram", "phiHistogram", kTH1F, {axisPhi}); + JEhistos.add("FJetaHistogram", "FJetaHistogram", kTH1F, {axisEta}); + JEhistos.add("FJphiHistogram", "FJphiHistogram", kTH1F, {axisPhi}); + JEhistos.add("FJptHistogram", "FJptHistogram", kTH1F, {axisPt}); + JEhistos.add("FJetaHistogram_MCRec", "FJetaHistogram_MCRec", kTH1F, {axisEta}); + JEhistos.add("FJphiHistogram_MCRec", "FJphiHistogram_MCRec", kTH1F, {axisPhi}); + JEhistos.add("FJptHistogram_MCRec", "FJptHistogram_MCRec", kTH1F, {axisPt}); + JEhistos.add("FJetaHistogram_MCTrue", "FJetaHistogram_MCTrue", kTH1F, {axisEta}); + JEhistos.add("FJphiHistogram_MCTrue", "FJphiHistogram_MCTrue", kTH1F, {axisPhi}); + JEhistos.add("FJptHistogram_MCTrue", "FJptHistogram_MCTrue", kTH1F, {axisPt}); + + JEhistos.add("FJnchHistogram", "FJnchHistogram", kTH1F, {MultAxis}); + JEhistos.add("nEvents", "nEvents", kTH1F, {{4, 0.0, 4.0}}); + JEhistos.add("nEvents_MCRec", "nEvents_MCRec", kTH1F, {{4, 0.0, 4.0}}); + JEhistos.add("nEvents_MCGen", "nEvents_MCGen", kTH1F, {{4, 0.0, 4.0}}); + JEhistos.add("nEvents_MCRec_MATCHED", "nEvents_MCRec_MATCHED", kTH1F, {{4, 0.0, 4.0}}); + JEhistos.add("nEvents_MCGen_MATCHED", "nEvents_MCGen_MATCHED", kTH1F, {{4, 0.0, 4.0}}); + + JEhistos.add("hMCTrue_nonmatch_hUSS_Kangle_v_pt", "hMCTrue_nonmatch_hUSS_Kangle_v_pt", kTH2F, {axisEta, PtAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_Kangle_v_pt", "hMCRec_nonmatch_hUSS_Kangle_v_pt", kTH2F, {axisEta, PtAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_INSIDE_pt_v_eta", "hMCRec_nonmatch_hUSS_INSIDE_pt_v_eta", kTH2F, {PtAxis, axisEta}); + JEhistos.add("hMCTrue_nonmatch_hUSS_INSIDE_pt_v_eta", "hMCTrue_nonmatch_hUSS_INSIDE_pt_v_eta", kTH2F, {PtAxis, axisEta}); + JEhistos.add("JetVsPhi_GEN", "JetVsPhi_GEN", kTH2F, {{4000, 0., 200.}, {200, 0, 20.0}}); + JEhistos.add("JetVsPhi_REC", "JetVsPhi_REC", kTH2F, {{4000, 0., 200.}, {200, 0, 20.0}}); + JEhistos.add("nJetsPerEvent", "nJetsPerEvent", kTH1F, {{10, 0.0, 10.0}}); + + JEhistos.add("hDCArToPv", "DCArToPv", kTH1F, {{300, 0.0, 3.0}}); + JEhistos.add("hDCAzToPv", "DCAzToPv", kTH1F, {{300, 0.0, 3.0}}); + JEhistos.add("rawpT", "rawpT", kTH1F, {{1000, 0.0, 10.0}}); + JEhistos.add("rawDpT", "rawDpT", kTH2F, {{1000, 0.0, 10.0}, {300, -1.5, 1.5}}); + JEhistos.add("hIsPrim", "hIsPrim", kTH1F, {{2, -0.5, +1.5}}); + JEhistos.add("hIsGood", "hIsGood", kTH1F, {{2, -0.5, +1.5}}); + JEhistos.add("hIsPrimCont", "hIsPrimCont", kTH1F, {{2, -0.5, +1.5}}); + JEhistos.add("hFindableTPCClusters", "hFindableTPCClusters", kTH1F, {{200, 0, 200}}); + JEhistos.add("hFindableTPCRows", "hFindableTPCRows", kTH1F, {{200, 0, 200}}); + JEhistos.add("hClustersVsRows", "hClustersVsRows", kTH1F, {{200, 0, 2}}); + JEhistos.add("hTPCChi2", "hTPCChi2", kTH1F, {{200, 0, 100}}); + JEhistos.add("hITSChi2", "hITSChi2", kTH1F, {{200, 0, 100}}); + + JEhistos.add("hUSS", "hUSS", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hUSS_1D", "hUSS_1D", kTH1F, {MinvAxis}); + JEhistos.add("hUSS_1D_2_3", "hUSS_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hLSS", "hLSS", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hLSS_1D", "hLSS_1D", kTH1F, {MinvAxis}); + JEhistos.add("hLSS_1D_2_3", "hLSS_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hUSS_INSIDE", "hUSS_INSIDE", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hUSS_INSIDE_1D", "hUSS_INSIDE_1D", kTH1F, {MinvAxis}); + JEhistos.add("hUSS_INSIDE_1D_2_3", "hUSS_INSIDE_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hLSS_INSIDE", "hLSS_INSIDE", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hLSS_INSIDE_1D", "hLSS_INSIDE_1D", kTH1F, {MinvAxis}); + JEhistos.add("hLSS_INSIDE_1D_2_3", "hLSS_INSIDE_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hUSS_OUTSIDE", "hUSS_OUTSIDE", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hUSS_OUTSIDE_1D", "hUSS_OUTSIDE_1D", kTH1F, {MinvAxis}); + JEhistos.add("hUSS_OUTSIDE_1D_2_3", "hUSS_OUTSIDE_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hLSS_OUTSIDE", "hLSS_OUTSIDE", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hLSS_OUTSIDE_1D", "hLSS_OUTSIDE_1D", kTH1F, {MinvAxis}); + JEhistos.add("hLSS_OUTSIDE_1D_2_3", "hLSS_OUTSIDE_1D_2_3", kTH1F, {MinvAxis}); + + JEhistos.add("hMCTrue_hUSS_INSIDE", "hMCTrue_hUSS_INSIDE", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCTrue_hUSS_INSIDE_1D", "hMCTrue_hUSS_INSIDE_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCTrue_hUSS_INSIDE_1D_2_3", "hMCTrue_hUSS_INSIDE_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hMCTrue_hUSS_OUTSIDE", "hMCTrue_hUSS_OUTSIDE", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCTrue_hUSS_OUTSIDE_1D", "hMCTrue_hUSS_OUTSIDE_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCTrue_hUSS_OUTSIDE_1D_2_3", "hMCTrue_hUSS_OUTSIDE_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hMCTrue_hUSS_OUTSIDE_TRIG", "hMCTrue_hUSS_OUTSIDE_TRIG", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCTrue_hUSS_OUTSIDE_TRIG_1D", "hMCTrue_hUSS_OUTSIDE_TRIG_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCTrue_hUSS_OUTSIDE_TRIG_1D_2_3", "hMCTrue_hUSS_OUTSIDE_TRIG_1D_2_3", kTH1F, {MinvAxis}); + + JEhistos.add("hMCTrue_nonmatch_hUSS_INSIDE", "hMCTrue_nonmatch_hUSS_INSIDE", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCTrue_nonmatch_hUSS_INSIDE_1D", "hMCTrue_nonmatch_hUSS_INSIDE_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCTrue_nonmatch_hUSS_INSIDE_1D_2_3", "hMCTrue_nonmatch_hUSS_INSIDE_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hMCTrue_nonmatch_hUSS_OUTSIDE", "hMCTrue_nonmatch_hUSS_OUTSIDE", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCTrue_nonmatch_hUSS_OUTSIDE_1D", "hMCTrue_nonmatch_hUSS_OUTSIDE_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCTrue_nonmatch_hUSS_OUTSIDE_1D_2_3", "hMCTrue_nonmatch_hUSS_OUTSIDE_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG", "hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG_1D", "hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG_1D_2_3", "hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG_1D_2_3", kTH1F, {MinvAxis}); + + JEhistos.add("hMCRec_hUSS", "hMCRec_hUSS", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCRec_hUSS_1D", "hMCRec_hUSS_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_hUSS_1D_2_3", "hMCRec_hUSS_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_hUSS_INSIDE", "hMCRec_hUSS_INSIDE", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCRec_hUSS_INSIDE_1D", "hMCRec_hUSS_INSIDE_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_hUSS_INSIDE_1D_2_3", "hMCRec_hUSS_INSIDE_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_hUSS_OUTSIDE", "hMCRec_hUSS_OUTSIDE", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCRec_hUSS_OUTSIDE_1D", "hMCRec_hUSS_OUTSIDE_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_hUSS_OUTSIDE_1D_2_3", "hMCRec_hUSS_OUTSIDE_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_hUSS_OUTSIDE_TRIG", "hMCRec_hUSS_OUTSIDE_TRIG", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCRec_hUSS_OUTSIDE_TRIG_1D", "hMCRec_hUSS_OUTSIDE_TRIG_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_hUSS_OUTSIDE_TRIG_1D_2_3", "hMCRec_hUSS_OUTSIDE_TRIG_1D_2_3", kTH1F, {MinvAxis}); + + JEhistos.add("hMCRec_nonmatch_hUSS", "hMCRec_nonmatch_hUSS", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_1D", "hMCRec_nonmatch_hUSS_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_1D_2_3", "hMCRec_nonmatch_hUSS_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_INSIDE", "hMCRec_nonmatch_hUSS_INSIDE", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_INSIDE_1D", "hMCRec_nonmatch_hUSS_INSIDE_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_INSIDE_1D_2_3", "hMCRec_nonmatch_hUSS_INSIDE_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_OUTSIDE", "hMCRec_nonmatch_hUSS_OUTSIDE", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_OUTSIDE_1D", "hMCRec_nonmatch_hUSS_OUTSIDE_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_OUTSIDE_1D_2_3", "hMCRec_nonmatch_hUSS_OUTSIDE_1D_2_3", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_OUTSIDE_TRIG", "hMCRec_nonmatch_hUSS_OUTSIDE_TRIG", kTH3F, {dRAxis, PtAxis, MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_OUTSIDE_TRIG_1D", "hMCRec_nonmatch_hUSS_OUTSIDE_TRIG_1D", kTH1F, {MinvAxis}); + JEhistos.add("hMCRec_nonmatch_hUSS_OUTSIDE_TRIG_1D_2_3", "hMCRec_nonmatch_hUSS_OUTSIDE_TRIG_1D_2_3", kTH1F, {MinvAxis}); + + // EVENT SELECTION + eventSelection = jetderiveddatautilities::initialiseEventSelection(static_cast(cfgeventSelections)); + + } // end of init + + double massKa = o2::constants::physics::MassKPlus; + double massPi = o2::constants::physics::MassPiMinus; + + using EventCandidates = soa::Join; // , aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs + using TrackCandidates = soa::Join; + Filter jetCuts = aod::jet::pt > cfgjetPtMin&& aod::jet::r == nround(cfgjetR.node() * 100.0f); + + // Function for track quality cuts + ///////////////////////////////////////////////////////////////////////////// + ///////////////////////////////////////////////////////////////////////////// + template + bool trackSelection(const TrackType track) + { + // basic track cuts + if (track.pt() < cfgtrkMinPt) + return false; + + if (std::abs(track.eta()) > cfgtrkMaxEta) + return false; + + if (std::abs(track.dcaXY()) > cfgMaxDCArToPVcut) + return false; + + if (std::abs(track.dcaZ()) > cfgMaxDCAzToPVcut) + return false; + + if (cfgPrimaryTrack && !track.isPrimaryTrack()) + return false; + + if (track.tpcNClsFindable() < cfgnFindableTPCClusters) + return false; + + if (track.tpcNClsCrossedRows() < cfgnTPCCrossedRows) + return false; + + if (track.tpcCrossedRowsOverFindableCls() > cfgnRowsOverFindable) + return false; + + if (track.tpcChi2NCl() > cfgnTPCChi2) + return false; + + if (track.itsChi2NCl() > cfgnITSChi2) + return false; + + if (cfgConnectedToPV && !track.isPVContributor()) + return false; + + return true; + }; + ///////////////////////////////////////////////////////////////////////////// + ///////////////////////////////////////////////////////////////////////////// + + template + bool trackPIDKaon(const T& candidate) + { + bool tpcPIDPassed{false}, tofPIDPassed{false}; + if (std::abs(candidate.tpcNSigmaKa()) < cfgnTPCPID) + tpcPIDPassed = true; + + if (candidate.hasTOF()) { + if (std::abs(candidate.tofNSigmaKa()) < cfgnTOFPID) { + tofPIDPassed = true; + } + } else { + tofPIDPassed = true; + } + if (tpcPIDPassed && tofPIDPassed) { + return true; + } + return false; + } + + template + bool trackPIDPion(const T& candidate) + { + bool tpcPIDPassed{false}, tofPIDPassed{false}; + if (std::abs(candidate.tpcNSigmaPi()) < cfgnTPCPID) + tpcPIDPassed = true; + + if (candidate.hasTOF()) { + if (std::abs(candidate.tofNSigmaPi()) < cfgnTOFPID) { + tofPIDPassed = true; + } + } else { + tofPIDPassed = true; + } + if (tpcPIDPassed && tofPIDPassed) { + return true; + } + return false; + } + ///////////////////////////////////////////////////////////////////////////// + ///////////////////////////////////////////////////////////////////////////// + + template + void minvReconstruction(double mult, const TracksType& trk1, const TracksType& trk2, const JetType& jets) + { + TLorentzVector lDecayDaughter1, lDecayDaughter2, lResonance; + + if (!trackSelection(trk1) || !trackSelection(trk2)) + return; + + if (!trackPIDKaon(trk1) || !trackPIDPion(trk2)) + return; + + if (trk1.index() == trk2.index()) + return; // We need to run (0,1), (1,0) pairs as well. but same id pairs are not needed. + + lDecayDaughter1.SetXYZM(trk1.px(), trk1.py(), trk1.pz(), massKa); + lDecayDaughter2.SetXYZM(trk2.px(), trk2.py(), trk2.pz(), massPi); + lResonance = lDecayDaughter1 + lDecayDaughter2; + + if (std::abs(lResonance.Eta()) > cfgtrkMaxEta) + return; + + ///////////////////////////////////////////////////////////////////////////// + // Fill Global Event Minv + if (trk1.sign() * trk2.sign() < 0) { + if (!IsMC) { + JEhistos.fill(HIST("hUSS_1D"), lResonance.M()); + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hUSS_1D_2_3"), lResonance.M()); + JEhistos.fill(HIST("hUSS"), mult, lResonance.Pt(), lResonance.M()); + } + + if (IsMC) { + JEhistos.fill(HIST("hMCRec_hUSS_1D"), lResonance.M()); + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCRec_hUSS_1D_2_3"), lResonance.M()); + JEhistos.fill(HIST("hMCRec_hUSS"), mult, lResonance.Pt(), lResonance.M()); + } + + } else if (trk1.sign() * trk2.sign() > 0) { + + JEhistos.fill(HIST("hLSS_1D"), lResonance.M()); + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hLSS_1D_2_3"), lResonance.M()); + JEhistos.fill(HIST("hLSS"), mult, lResonance.Pt(), lResonance.M()); + } + ///////////////////////////////////////////////////////////////////////////// + + bool jetFlag = false; + for (auto const& jet : jets) { + + double phidiff = TVector2::Phi_mpi_pi(jet.phi() - lResonance.Phi()); + double etadiff = jet.eta() - lResonance.Eta(); + double R = TMath::Sqrt((etadiff * etadiff) + (phidiff * phidiff)); + if (R < cfgjetR) + jetFlag = true; + } + + ///////////////////////////////////////////////////////////////////////////// + // Fill inside Jet + if (jetFlag) { + if (trk1.sign() * trk2.sign() < 0) { + if (!IsMC) { + JEhistos.fill(HIST("hUSS_INSIDE_1D"), lResonance.M()); + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hUSS_INSIDE_1D_2_3"), lResonance.M()); + JEhistos.fill(HIST("hUSS_INSIDE"), mult, lResonance.Pt(), lResonance.M()); + } + + if (IsMC) { + JEhistos.fill(HIST("hMCRec_hUSS_INSIDE_1D"), lResonance.M()); + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCRec_hUSS_INSIDE_1D_2_3"), lResonance.M()); + JEhistos.fill(HIST("hMCRec_hUSS_INSIDE"), mult, lResonance.Pt(), lResonance.M()); + } + + } else if (trk1.sign() * trk2.sign() > 0) { + + JEhistos.fill(HIST("hLSS_INSIDE_1D"), lResonance.M()); + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hLSS_INSIDE_1D_2_3"), lResonance.M()); + JEhistos.fill(HIST("hLSS_INSIDE"), mult, lResonance.Pt(), lResonance.M()); + } + } // jetflag + ///////////////////////////////////////////////////////////////////////////// + + ///////////////////////////////////////////////////////////////////////////// + // Fill outside Jet + if (!jetFlag) { + if (trk1.sign() * trk2.sign() < 0) { + if (!IsMC) { + JEhistos.fill(HIST("hUSS_OUTSIDE_1D"), lResonance.M()); + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hUSS_OUTSIDE_1D_2_3"), lResonance.M()); + JEhistos.fill(HIST("hUSS_OUTSIDE"), mult, lResonance.Pt(), lResonance.M()); + } + + if (IsMC) { + JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_1D"), lResonance.M()); + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_1D_2_3"), lResonance.M()); + JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE"), mult, lResonance.Pt(), lResonance.M()); + } + + } else if (trk1.sign() * trk2.sign() > 0) { + + JEhistos.fill(HIST("hLSS_OUTSIDE_1D"), lResonance.M()); + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hLSS_OUTSIDE_1D_2_3"), lResonance.M()); + JEhistos.fill(HIST("hLSS_OUTSIDE"), mult, lResonance.Pt(), lResonance.M()); + } + } //! jetflag + ///////////////////////////////////////////////////////////////////////////// + } // MinvReconstruction + + ///////////////////////////////////////////////////////////////////////////// + ///////////////////////////////////////////////////////////////////////////// + int nEvents = 0; + void processJetTracks(aod::JCollision const& collision, soa::Filtered> const& chargedjets, soa::Join const& tracks, TrackCandidates const&) + { + if (cDebugLevel > 0) { + nEvents++; + if ((nEvents + 1) % 10000 == 0) + std::cout << "Processed Data Events: " << nEvents << std::endl; + } + JEhistos.fill(HIST("nEvents"), 0.5); + + if (fabs(collision.posZ()) > cfgVtxCut) + return; + if (!jetderiveddatautilities::selectCollision(collision, jetderiveddatautilities::JCollisionSel::sel8)) + return; + + for (auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(tracks, tracks))) { + auto trk1 = track1.track_as>(); + auto trk2 = track2.track_as>(); + minvReconstruction(1.0, trk1, trk2, chargedjets); + } + int nJets = 0; + for (auto chargedjet : chargedjets) { + JEhistos.fill(HIST("FJetaHistogram"), chargedjet.eta()); + JEhistos.fill(HIST("FJphiHistogram"), chargedjet.phi()); + JEhistos.fill(HIST("FJptHistogram"), chargedjet.pt()); + // JEhistos.fill(HIST("FJnchHistogram"), chargedjet.tracks().size()); + nJets++; + } + + JEhistos.fill(HIST("nJetsPerEvent"), nJets); + + JEhistos.fill(HIST("nEvents"), 1.5); + // return; + + for (auto& trackC : tracks) { + auto originalTrack = trackC.track_as>(); + JEhistos.fill(HIST("hDCArToPv"), originalTrack.dcaXY()); + JEhistos.fill(HIST("hDCAzToPv"), originalTrack.dcaZ()); + JEhistos.fill(HIST("rawpT"), originalTrack.pt()); + JEhistos.fill(HIST("rawDpT"), trackC.pt(), trackC.pt() - originalTrack.pt()); + JEhistos.fill(HIST("hIsPrim"), originalTrack.isPrimaryTrack()); + JEhistos.fill(HIST("hIsGood"), originalTrack.isGlobalTrackWoDCA()); + JEhistos.fill(HIST("hIsPrimCont"), originalTrack.isPVContributor()); + JEhistos.fill(HIST("hFindableTPCClusters"), originalTrack.tpcNClsFindable()); + JEhistos.fill(HIST("hFindableTPCRows"), originalTrack.tpcNClsCrossedRows()); + JEhistos.fill(HIST("hClustersVsRows"), originalTrack.tpcCrossedRowsOverFindableCls()); + JEhistos.fill(HIST("hTPCChi2"), originalTrack.tpcChi2NCl()); + JEhistos.fill(HIST("hITSChi2"), originalTrack.itsChi2NCl()); + + if (!trackSelection(originalTrack)) + continue; + + JEhistos.fill(HIST("etaHistogram"), trackC.eta()); + JEhistos.fill(HIST("phiHistogram"), trackC.phi()); + } // JTrack Loop + + }; // Process Switch + PROCESS_SWITCH(kstarInJets, processJetTracks, "process JE Framework", true); + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + ////////////////////////////////////////////////////////////////MC STUFF//////////////////////////////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + using myCompleteTracks = soa::Join; + using myCompleteJetTracks = soa::Join; + int nJEEvents = 0; + int nprocessRecEvents = 0; + void processRec(o2::aod::JCollision const& collision, myCompleteJetTracks const& tracks, soa::Filtered const& mcdjets, aod::McParticles const&, myCompleteTracks const& /*originalTracks*/) + { + if (cDebugLevel > 0) { + nprocessRecEvents++; + if ((nprocessRecEvents + 1) % 10000 == 0) + std::cout << "processRec: " << nprocessRecEvents << std::endl; + } + + JEhistos.fill(HIST("nEvents_MCRec"), 0.5); + if (fabs(collision.posZ()) > cfgVtxCut) + return; + if (!jetderiveddatautilities::selectCollision(collision, jetderiveddatautilities::JCollisionSel::sel8)) + return; + + bool INELgt0 = false; + for (const auto& track : tracks) { + if (fabs(track.eta()) < cfgtrkMaxEta) { + INELgt0 = true; + break; + } + } + if (!INELgt0) + return; + + JEhistos.fill(HIST("nEvents_MCRec"), 1.5); + + std::vector mcd_pt{}; + std::vector mcd_phi{}; + std::vector mcd_eta{}; + + bool hasJets = kFALSE; + for (auto& mcdjet : mcdjets) { + hasJets = kTRUE; + mcd_pt.push_back(mcdjet.pt()); + mcd_eta.push_back(mcdjet.eta()); + mcd_phi.push_back(mcdjet.phi()); + registry.fill(HIST("h_jet_pt"), mcdjet.pt()); + registry.fill(HIST("h_jet_eta"), mcdjet.eta()); + registry.fill(HIST("h_jet_phi"), mcdjet.phi()); + } + if (hasJets) { + JEhistos.fill(HIST("nEvents_MCRec"), 2.5); + } + + // Track Eff + for (const auto& track : tracks) { + auto originalTrack = track.track_as(); + if (!trackSelection(originalTrack)) + continue; + if (cfgSimPID) { + if (!trackPIDKaon(originalTrack)) + continue; + if (!trackPIDPion(originalTrack)) + continue; + } + if (track.has_mcParticle()) { + auto mcParticle = track.mcParticle(); + + if (mcParticle.isPhysicalPrimary() && fabs(mcParticle.eta()) <= cfgtrkMaxEta) { + if (abs(mcParticle.pdgCode()) == 211) + JEhistos.fill(HIST("ptJEHistogramPion"), mcParticle.pt()); + if (abs(mcParticle.pdgCode()) == 321) + JEhistos.fill(HIST("ptJEHistogramKaon"), mcParticle.pt()); + if (abs(mcParticle.pdgCode()) == 2212) + JEhistos.fill(HIST("ptJEHistogramProton"), mcParticle.pt()); + } + } + for (const auto& track2 : tracks) { + auto originalTrack2 = track2.track_as(); + if (!trackSelection(originalTrack2)) + continue; + if (cfgSimPID) { + if (!trackPIDKaon(originalTrack)) + continue; + if (!trackPIDPion(originalTrack)) + continue; + } + + if (originalTrack.index() >= originalTrack2.index()) + continue; + + if (fabs(originalTrack.eta()) > cfgtrkMaxEta || fabs(originalTrack2.eta()) > cfgtrkMaxEta) + continue; + + // check PID + + if (track.has_mcParticle() && track2.has_mcParticle()) { + auto part1 = track.mcParticle(); + auto part2 = track2.mcParticle(); + if (fabs(part1.pdgCode()) != 321) + continue; // Not Kaon + if (fabs(part2.pdgCode()) != 321) + continue; // Not Kaon + if (!part1.has_mothers()) + continue; // Not decaying Kaon + if (!part2.has_mothers()) + continue; // Not decaying Kaon + + std::vector mothers1{}; + std::vector mothers1PDG{}; + for (auto& part1_mom : part1.mothers_as()) { + mothers1.push_back(part1_mom.globalIndex()); + mothers1PDG.push_back(part1_mom.pdgCode()); + } + + std::vector mothers2{}; + std::vector mothers2PDG{}; + for (auto& part2_mom : part2.mothers_as()) { + mothers2.push_back(part2_mom.globalIndex()); + mothers2PDG.push_back(part2_mom.pdgCode()); + } + + if (mothers1PDG[0] != 333) + continue; // mother not phi + if (mothers2PDG[0] != 333) + continue; // mother not phi + if (mothers1[0] != mothers2[0]) + continue; // Kaons not from the same phi + + TLorentzVector lDecayDaughter1, lDecayDaughter2, lResonance; + lDecayDaughter1.SetXYZM(originalTrack.px(), originalTrack.py(), originalTrack.pz(), massKa); + lDecayDaughter2.SetXYZM(originalTrack2.px(), originalTrack2.py(), originalTrack2.pz(), massPi); + lResonance = lDecayDaughter1 + lDecayDaughter2; + if (fabs(lResonance.Eta()) > cfgtrkMaxEta) + continue; + JEhistos.fill(HIST("ptJEHistogramPhi"), lResonance.Pt()); + + bool jetFlag = false; + for (int i = 0; i < mcd_pt.size(); i++) { + double phidiff = TVector2::Phi_mpi_pi(mcd_phi[i] - lResonance.Phi()); + double etadiff = mcd_eta[i] - lResonance.Eta(); + double R = TMath::Sqrt((etadiff * etadiff) + (phidiff * phidiff)); + + double phidiff_K1 = TVector2::Phi_mpi_pi(mcd_phi[i] - lDecayDaughter1.Phi()); + double etadiff_K1 = mcd_eta[i] - lDecayDaughter1.Eta(); + double R_K1 = TMath::Sqrt((etadiff_K1 * etadiff_K1) + (phidiff_K1 * phidiff_K1)); + double phidiff_K2 = TVector2::Phi_mpi_pi(mcd_phi[i] - lDecayDaughter2.Phi()); + double etadiff_K2 = mcd_eta[i] - lDecayDaughter2.Eta(); + double R_K2 = TMath::Sqrt((etadiff_K2 * etadiff_K2) + (phidiff_K2 * phidiff_K2)); + if (R < cfgjetR) { + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_Kangle_v_pt"), R_K1, lResonance.Pt()); + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_Kangle_v_pt"), R_K2, lResonance.Pt()); + } + if (R < cfgjetR) + jetFlag = true; + } + + if (jetFlag) { + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_INSIDE_pt_v_eta"), lResonance.Pt(), lResonance.Eta()); + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_INSIDE_1D"), lResonance.M()); + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_INSIDE_1D_2_3"), lResonance.M()); + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_INSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + + } else if (!jetFlag && mcd_pt.size() > 0) { + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE_TRIG_1D"), lResonance.M()); + + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE_TRIG_1D_2_3"), lResonance.M()); + + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE_TRIG"), 1.0, lResonance.Pt(), lResonance.M()); + + } else if (!jetFlag) { + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE_1D"), lResonance.M()); + + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE_1D_2_3"), lResonance.M()); + + JEhistos.fill(HIST("hMCRec_nonmatch_hUSS_OUTSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + } //! jetflag + + if (hasJets) { + JEhistos.fill(HIST("ptJEHistogramPhi_JetTrigger"), lResonance.Pt()); + auto triggerjet = std::min_element(mcd_pt.begin(), mcd_pt.end()); + double triggerjet_pt = *triggerjet; + JEhistos.fill(HIST("JetVsPhi_REC"), triggerjet_pt, lResonance.Pt()); + } + JEhistos.fill(HIST("minvJEHistogramPhi"), lResonance.M()); + } // mcpart check + } // tracks2 + } // tracks1 + // Jet Eff + } + PROCESS_SWITCH(kstarInJets, processRec, "pikp detector level MC JE", true); + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + int nprocessSimEvents = 0; + // Preslice slice = o2::aod::JCollision::collisionId; + void processSim(o2::aod::JMcCollision const& collision, soa::SmallGroups> const& recocolls, aod::JMcParticles const& mcParticles, soa::Filtered const& mcpjets) + { + if (cDebugLevel > 0) { + nprocessSimEvents++; + if ((nprocessSimEvents + 1) % 10000 == 0) + std::cout << "processSim: " << nprocessSimEvents << std::endl; + } + + JEhistos.fill(HIST("nEvents_MCGen"), 0.5); + + if (recocolls.size() <= 0) // not reconstructed + return; + + for (auto& recocoll : recocolls) { // poorly reconstructed + if (!jetderiveddatautilities::selectCollision(recocoll, jetderiveddatautilities::JCollisionSel::sel8)) + return; + } + + if (fabs(collision.posZ()) > cfgVtxCut) // bad vertex + return; + bool INELgt0 = false; + for (const auto& mcParticle : mcParticles) { + if (fabs(mcParticle.eta()) < cfgtrkMaxEta) { + INELgt0 = true; + break; + } + } + + if (!INELgt0) // not INEL + return; + std::vector mcp_pt{}; + std::vector mcp_phi{}; + std::vector mcp_eta{}; + JEhistos.fill(HIST("nEvents_MCGen"), 1.5); + + bool hasJets = kFALSE; + // Check jets + for (auto& mcpjet : mcpjets) { + hasJets = kTRUE; + mcp_pt.push_back(mcpjet.pt()); + mcp_eta.push_back(mcpjet.eta()); + mcp_phi.push_back(mcpjet.phi()); + + registry.fill(HIST("h_part_jet_pt"), mcpjet.pt()); + registry.fill(HIST("h_part_jet_eta"), mcpjet.eta()); + registry.fill(HIST("h_part_jet_phi"), mcpjet.phi()); + } + if (hasJets) + JEhistos.fill(HIST("nEvents_MCGen"), 2.5); + + // Check pikp and phi + for (const auto& mcParticle : mcParticles) { + if (mcParticle.isPhysicalPrimary() && fabs(mcParticle.eta()) <= cfgtrkMaxEta) { // watch out for context!!! + if (abs(mcParticle.pdgCode()) == 211) + JEhistos.fill(HIST("ptGeneratedPion"), mcParticle.pt()); + if (abs(mcParticle.pdgCode()) == 321) + JEhistos.fill(HIST("ptGeneratedKaon"), mcParticle.pt()); + if (abs(mcParticle.pdgCode()) == 2212) + JEhistos.fill(HIST("ptGeneratedProton"), mcParticle.pt()); + } + if (fabs(mcParticle.eta()) <= cfgtrkMaxEta) { // watch out for context!!! + TLorentzVector lResonance; + lResonance.SetPxPyPzE(mcParticle.px(), mcParticle.py(), mcParticle.pz(), mcParticle.e()); + if (abs(mcParticle.pdgCode()) == 333) { + JEhistos.fill(HIST("ptGeneratedPhi_ALLBR"), mcParticle.pt()); + + bool skip = false; + // First we check for Forced BR + if (mcParticle.has_daughters()) + for (auto& dgth : mcParticle.daughters_as()) + if (fabs(dgth.pdgCode()) != 321) + skip = true; + + if (skip && cfgBR) + continue; + JEhistos.fill(HIST("ptGeneratedPhi"), mcParticle.pt()); + JEhistos.fill(HIST("mGeneratedPhi"), lResonance.M()); + + ////////////////////////////Implementation of phi finding + + TLorentzVector lResonance; + lResonance.SetPxPyPzE(mcParticle.px(), mcParticle.py(), mcParticle.pz(), mcParticle.e()); + bool jetFlag = false; + for (int i = 0; i < mcp_pt.size(); i++) { + double phidiff = TVector2::Phi_mpi_pi(mcp_phi[i] - lResonance.Phi()); + double etadiff = mcp_eta[i] - lResonance.Eta(); + double R = TMath::Sqrt((etadiff * etadiff) + (phidiff * phidiff)); + if (mcParticle.has_daughters()) { + for (auto& dgth : mcParticle.daughters_as()) { + double phidiff_K = TVector2::Phi_mpi_pi(mcp_phi[i] - dgth.phi()); + double etadiff_K = mcp_eta[i] - dgth.eta(); + double R_K = TMath::Sqrt((etadiff_K * etadiff_K) + (phidiff_K * phidiff_K)); + if (R < cfgjetR) + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_Kangle_v_pt"), R_K, lResonance.Pt()); + } + } + if (R < cfgjetR) + jetFlag = true; + } + + if (jetFlag) { + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_INSIDE_pt_v_eta"), lResonance.Pt(), lResonance.Eta()); + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_INSIDE_1D"), lResonance.M()); + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_INSIDE_1D_2_3"), lResonance.M()); + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_INSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + + } else if (!jetFlag && mcp_pt.size() > 0) { + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG_1D"), lResonance.M()); + + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG_1D_2_3"), lResonance.M()); + + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE_TRIG"), 1.0, lResonance.Pt(), lResonance.M()); + + } else if (!jetFlag) { + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE_1D"), lResonance.M()); + + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE_1D_2_3"), lResonance.M()); + + JEhistos.fill(HIST("hMCTrue_nonmatch_hUSS_OUTSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + + } //! jetflag + + ////////////////////////////Phi found + if (hasJets) { + JEhistos.fill(HIST("ptGeneratedPhi_JetTrigger"), mcParticle.pt()); + auto triggerjet = std::min_element(mcp_pt.begin(), mcp_pt.end()); + double triggerjet_pt = *triggerjet; + JEhistos.fill(HIST("JetVsPhi_GEN"), triggerjet_pt, mcParticle.pt()); + } // check for jets + + } // check for phi + } // check for rapidity + } // loop over particles + } // process switch + PROCESS_SWITCH(kstarInJets, processSim, "pikp particle level MC", true); + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + using JetMCPTable = soa::Filtered>; + using JetMCDTable = soa::Filtered>; + + // void processMatchedGen(o2::aod::JMcCollision const& collision, aod::JMcParticles const& mcParticles, soa::Filtered const& mcpjets) + int nprocessSimJEEvents = 0; + void processMatchedGen(aod::JMcCollision const& collision, + soa::SmallGroups> const& recocolls, + JetMCDTable const& /*mcdjets*/, + JetMCPTable const& mcpjets, + aod::JMcParticles const& mcParticles) + + { + if (cDebugLevel > 0) { + nprocessSimJEEvents++; + if ((nprocessSimJEEvents + 1) % 10000 == 0) + std::cout << "JEvents: " << nprocessSimJEEvents << std::endl; + } + + JEhistos.fill(HIST("nEvents_MCGen_MATCHED"), 0.5); + + if (fabs(collision.posZ()) > cfgVtxCut) + return; + + if (recocolls.size() <= 0) // not reconstructed + return; + for (auto& recocoll : recocolls) { // poorly reconstructed + if (!jetderiveddatautilities::selectCollision(recocoll, jetderiveddatautilities::JCollisionSel::sel8)) + return; + } + + bool INELgt0 = false; + for (const auto& mcParticle : mcParticles) { + if (TMath::Abs(mcParticle.eta()) < cfgtrkMaxEta) { + INELgt0 = true; + break; + } + } + if (!INELgt0) + return; + + JEhistos.fill(HIST("nEvents_MCGen_MATCHED"), 1.5); + + std::vector mcd_pt{}; + std::vector mcd_phi{}; + std::vector mcd_eta{}; + std::vector mcp_pt{}; + std::vector mcp_phi{}; + std::vector mcp_eta{}; + + bool hasJets = kFALSE; + + for (auto& mcpjet : mcpjets) { + if (!mcpjet.has_matchedJetGeo()) + continue; + + for (auto& mcdjet : mcpjet.template matchedJetGeo_as()) { + if (!mcdjet.has_matchedJetGeo()) + continue; + hasJets = kTRUE; + registry.fill(HIST("h_matched_GEN_jet_pt"), mcpjet.pt(), mcdjet.pt() - mcpjet.pt()); + registry.fill(HIST("h_matched_GEN_jet_phi"), mcpjet.phi(), mcdjet.phi() - mcpjet.phi()); + registry.fill(HIST("h_matched_GEN_jet_eta"), mcpjet.eta(), mcdjet.eta() - mcpjet.eta()); + + mcd_pt.push_back(mcdjet.pt()); + mcd_eta.push_back(mcdjet.eta()); + mcd_phi.push_back(mcdjet.phi()); + mcp_pt.push_back(mcpjet.pt()); + mcp_eta.push_back(mcpjet.eta()); + mcp_phi.push_back(mcpjet.phi()); + } // mcpjets + } // mcdjets + + if (hasJets) + JEhistos.fill(HIST("nEvents_MCGen_MATCHED"), 2.5); + + // First we do GEN part + for (const auto& mcParticle : mcParticles) { + if (fabs(mcParticle.eta()) > cfgtrkMaxEta) + continue; + + if (fabs(mcParticle.pdgCode()) == 333) { + bool skip = false; + // First we check for Forced BR + if (mcParticle.has_daughters()) + for (auto& dgth : mcParticle.daughters_as()) + if (fabs(dgth.pdgCode()) != 321) + skip = true; + + if (skip && cfgBR) + continue; + TLorentzVector lResonance; + lResonance.SetPxPyPzE(mcParticle.px(), mcParticle.py(), mcParticle.pz(), mcParticle.e()); + bool jetFlag = false; + for (int i = 0; i < mcp_pt.size(); i++) { + double phidiff = TVector2::Phi_mpi_pi(mcp_phi[i] - lResonance.Phi()); + double etadiff = mcp_eta[i] - lResonance.Eta(); + double R = TMath::Sqrt((etadiff * etadiff) + (phidiff * phidiff)); + if (R < cfgjetR) + jetFlag = true; + } + + if (jetFlag) { + JEhistos.fill(HIST("hMCTrue_hUSS_INSIDE_1D"), lResonance.M()); + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCTrue_hUSS_INSIDE_1D_2_3"), lResonance.M()); + JEhistos.fill(HIST("hMCTrue_hUSS_INSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + + } else if (!jetFlag && mcp_pt.size() > 0) { + JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE_TRIG_1D"), lResonance.M()); + + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE_TRIG_1D_2_3"), lResonance.M()); + + JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE_TRIG"), 1.0, lResonance.Pt(), lResonance.M()); + + } else if (!jetFlag) { + JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE_1D"), lResonance.M()); + + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE_1D_2_3"), lResonance.M()); + + JEhistos.fill(HIST("hMCTrue_hUSS_OUTSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + + } //! jetflag + } // chech for phi + } // MC Particles + } // main fcn + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + int nprocessRecJEEvents = 0; + // void processMatchedRec(o2::aod::JCollision const& collision, myCompleteJetTracks const& tracks, soa::Filtered const& mcdjets, aod::McParticles const&, myCompleteTracks const& originalTracks) + void processMatchedRec(aod::JCollision const& collision, + JetMCDTable const& mcdjets, + JetMCPTable const&, + myCompleteJetTracks const& tracks, + myCompleteTracks const&, + aod::McParticles const&) + { + + if (cDebugLevel > 0) { + nprocessRecJEEvents++; + if ((nprocessRecJEEvents + 1) % 10000 == 0) + std::cout << "JEvents: " << nprocessRecJEEvents << std::endl; + } + JEhistos.fill(HIST("nEvents_MCRec_MATCHED"), 0.5); + + if (fabs(collision.posZ()) > cfgVtxCut) + return; + if (!jetderiveddatautilities::selectCollision(collision, jetderiveddatautilities::JCollisionSel::sel8)) + return; + + bool INELgt0 = false; + for (const auto& track : tracks) { + if (fabs(track.eta()) < cfgtrkMaxEta) { + INELgt0 = true; + break; + } + } + if (!INELgt0) + return; + + JEhistos.fill(HIST("nEvents_MCRec_MATCHED"), 1.5); + + std::vector mcd_pt{}; + std::vector mcd_phi{}; + std::vector mcd_eta{}; + std::vector mcp_pt{}; + std::vector mcp_phi{}; + std::vector mcp_eta{}; + + bool hasJets = kFALSE; + for (auto& mcdjet : mcdjets) { + if (!mcdjet.has_matchedJetGeo()) + continue; + for (auto& mcpjet : mcdjet.template matchedJetGeo_as()) { + if (!mcpjet.has_matchedJetGeo()) + continue; + hasJets = kTRUE; + registry.fill(HIST("h_matched_REC_jet_pt"), mcpjet.pt(), mcdjet.pt() - mcpjet.pt()); + registry.fill(HIST("h_matched_REC_jet_phi"), mcpjet.phi(), mcdjet.phi() - mcpjet.phi()); + registry.fill(HIST("h_matched_REC_jet_eta"), mcpjet.eta(), mcdjet.eta() - mcpjet.eta()); + + mcd_pt.push_back(mcdjet.pt()); + mcd_eta.push_back(mcdjet.eta()); + mcd_phi.push_back(mcdjet.phi()); + mcp_pt.push_back(mcpjet.pt()); + mcp_eta.push_back(mcpjet.eta()); + mcp_phi.push_back(mcpjet.phi()); + } // mcpjets + } // mcdjets + // Now we do REC part + if (hasJets) + JEhistos.fill(HIST("nEvents_MCRec_MATCHED"), 2.5); + + for (const auto& track1 : tracks) { + auto trk1 = track1.track_as(); + for (const auto& track2 : tracks) { + auto trk2 = track2.track_as(); + if (trk1.index() >= trk2.index()) + continue; + if (fabs(trk1.eta()) > cfgtrkMaxEta || fabs(trk2.eta()) > cfgtrkMaxEta) + continue; + if ((trk1.sign() * trk2.sign()) > 0) + continue; // Not K+K- + if (trackSelection(trk1) && trackSelection(trk2)) { + if (cfgSimPID) { + if (!trackPIDKaon(trk1)) + continue; + if (!trackPIDPion(trk2)) + continue; + } + if (track1.has_mcParticle() && track2.has_mcParticle()) { + auto part1 = track1.mcParticle(); + auto part2 = track2.mcParticle(); + if (fabs(part1.pdgCode()) != 321) + continue; // Not Kaon + if (fabs(part2.pdgCode()) != 321) + continue; // Not Kaon + if (!part1.has_mothers()) + continue; // Not decaying Kaon + if (!part2.has_mothers()) + continue; // Not decaying Kaon + + std::vector mothers1{}; + std::vector mothers1PDG{}; + std::vector mothers1Pt{}; + for (auto& part1_mom : part1.mothers_as()) { + mothers1.push_back(part1_mom.globalIndex()); + mothers1PDG.push_back(part1_mom.pdgCode()); + mothers1Pt.push_back(part1_mom.pt()); + } + + std::vector mothers2{}; + std::vector mothers2PDG{}; + std::vector mothers2Pt{}; + for (auto& part2_mom : part2.mothers_as()) { + mothers2.push_back(part2_mom.globalIndex()); + mothers2PDG.push_back(part2_mom.pdgCode()); + mothers2Pt.push_back(part2_mom.pt()); + } + if (mothers1PDG[0] != 333) + continue; // mother not phi + if (mothers2PDG[0] != 333) + continue; // mother not phi + if (mothers1[0] != mothers2[0]) + continue; // Kaons not from the same phi + + TLorentzVector lDecayDaughter1, lDecayDaughter2, lResonance; + lDecayDaughter1.SetXYZM(trk1.px(), trk1.py(), trk1.pz(), massKa); + lDecayDaughter2.SetXYZM(trk2.px(), trk2.py(), trk2.pz(), massPi); + lResonance = lDecayDaughter1 + lDecayDaughter2; + + if (fabs(lResonance.Eta()) > cfgtrkMaxEta) + continue; + + bool jetFlag = false; + for (int i = 0; i < mcd_pt.size(); i++) { + double phidiff = TVector2::Phi_mpi_pi(mcd_phi[i] - lResonance.Phi()); + double etadiff = mcd_eta[i] - lResonance.Eta(); + double R = TMath::Sqrt((etadiff * etadiff) + (phidiff * phidiff)); + if (R < cfgjetR) + jetFlag = true; + + if (jetFlag) { // Fill Resp. Matrix + if (cDebugLevel > 0) { + std::cout << "******************************************" << std::endl; + std::cout << "Rec. Phi Pt: " << lResonance.Pt() << std::endl; + std::cout << "Rec. Jet Pt: " << mcd_pt[i] << std::endl; + std::cout << "Gen. Phi Pt: " << mothers1Pt[0] << std::endl; + std::cout << "Gen. Jet Pt: " << mcp_pt[i] << std::endl; + std::cout << "******************************************" << std::endl; + } + JEhistos.fill(HIST("Resp_Matrix_MATCHED"), lResonance.Pt(), mcd_pt[i], mothers1Pt[0], mcp_pt[i]); + } + } + + if (jetFlag) { + JEhistos.fill(HIST("hMCRec_hUSS_INSIDE_1D"), lResonance.M()); + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCRec_hUSS_INSIDE_1D_2_3"), lResonance.M()); + JEhistos.fill(HIST("hMCRec_hUSS_INSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + + } else if (!jetFlag && mcd_pt.size() > 0) { + JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_TRIG_1D"), lResonance.M()); + + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_TRIG_1D_2_3"), lResonance.M()); + + JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_TRIG"), 1.0, lResonance.Pt(), lResonance.M()); + + } else if (!jetFlag) { + JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_1D"), lResonance.M()); + + if (lResonance.Pt() > 2.0 && lResonance.Pt() < 3) + JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE_1D_2_3"), lResonance.M()); + + JEhistos.fill(HIST("hMCRec_hUSS_OUTSIDE"), 1.0, lResonance.Pt(), lResonance.M()); + + } //! jetflag + + } // pass track cut + } // has mc particle + + } // tracks 2 + } // tracks 1 + // tracks + } // main fcn + PROCESS_SWITCH(kstarInJets, processMatchedRec, "phi matched Rec level MC", true); +}; // end of main struct + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +}; diff --git a/compile_commands.json b/compile_commands.json new file mode 120000 index 00000000000..95d26ca7216 --- /dev/null +++ b/compile_commands.json @@ -0,0 +1 @@ +/Users/jimun/alice/sw/BUILD/6124bd93563403cdd582bf495e22b311782aa20d/O2Physics/compile_commands.json \ No newline at end of file