diff --git a/PWGLF/Tasks/Nuspex/CMakeLists.txt b/PWGLF/Tasks/Nuspex/CMakeLists.txt index b34ef39cf73..ffab28c2aa9 100644 --- a/PWGLF/Tasks/Nuspex/CMakeLists.txt +++ b/PWGLF/Tasks/Nuspex/CMakeLists.txt @@ -54,6 +54,12 @@ o2physics_add_dpl_workflow(spectra-tof PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) +#New +o2physics_add_dpl_workflow(spectra-tof-light + SOURCES spectraTOFLight.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(spectra-tof-run2 SOURCES spectraTOFRun2.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore diff --git a/PWGLF/Tasks/Nuspex/spectraTOFLight.cxx b/PWGLF/Tasks/Nuspex/spectraTOFLight.cxx new file mode 100644 index 00000000000..b137f8c4089 --- /dev/null +++ b/PWGLF/Tasks/Nuspex/spectraTOFLight.cxx @@ -0,0 +1,1212 @@ +// 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 spectraTOFLight.cxx +/// \author Nicolò Jacazio nicolo.jacazio@cern.ch, Flavio Pio Volpe +/// +/// \brief Task for the analysis of the spectra with the TPC and TOF detector. +/// + +#include "PWGLF/DataModel/LFParticleIdentification.h" // IWYU pragma: keep +#include "PWGLF/DataModel/mcCentrality.h" +#include "PWGLF/DataModel/spectraTOF.h" +#include "PWGLF/Utils/inelGt.h" + +#include "Common/CCDB/EventSelectionParams.h" +#include "Common/Core/RecoDecay.h" +#include "Common/Core/TrackSelection.h" +#include "Common/Core/TrackSelectionDefaults.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponseTOF.h" +#include "Common/DataModel/PIDResponseTPC.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::track; +using namespace o2::framework; +using namespace o2::framework::expressions; + +// Histograms +std::array, NpCharge> hDcaXYZ; +std::array, NpCharge> hDcaXYZPrm; +std::array, NpCharge> hDcaXYZStr; +std::array, NpCharge> hDcaXYZMat; + +// inel cut constants +static constexpr int EvSelInelGt0Cut = 1; + +// mc dca filling constants +static constexpr float DcaMaxCombinedSigma = 2.f; +static constexpr float DcaMaxTPCSigma = 1.f; +static constexpr float DcaTrkPtCut = 0.4f; + +// Spectra task +struct tofSpectra { + struct : ConfigurableGroup { + Configurable cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted z-vertex range"}; + Configurable cfgINELCut{"cfgINELCut", 0, "INEL event selection: 0 sel, 1 INEL>0"}; + Configurable requireSel8{"requireSel8", false, "Ask for sel8"}; + Configurable removeITSROFrameBorder{"removeITSROFrameBorder", false, "Remove TF border"}; + Configurable removeNoTimeFrameBorder{"removeNoTimeFrameBorder", false, "Remove TF border"}; + } evselOptions; + + struct : ConfigurableGroup { + Configurable cfgCutEtaMax{"cfgCutEtaMax", 0.8f, "Max eta range for tracks"}; + Configurable cfgCutEtaMin{"cfgCutEtaMin", -0.8f, "Min eta range for tracks"}; + Configurable cfgCutNsigma{"cfgCutNsigma", 100.0f, "nsigma cut range for tracks"}; + Configurable cfgCutY{"cfgCutY", 0.5f, "Y range for tracks"}; + } trkselOptions; + + TrackSelection customTrackCuts; + Configurable kaonIsPvContrib{"kaonIsPvContrib", false, "Flag to check if kaon tracks are from pv"}; + Configurable useCustomTrackCuts{"useCustomTrackCuts", false, "Flag to use custom track cuts"}; + Configurable itsPattern{"itsPattern", 0, "0 = Run3ITSibAny, 1 = Run3ITSallAny, 2 = Run3ITSall7Layers, 3 = Run3ITSibTwo"}; + Configurable requireITS{"requireITS", true, "Additional cut on the ITS requirement"}; + Configurable requireTPC{"requireTPC", true, "Additional cut on the TPC requirement"}; + Configurable requireGoldenChi2{"requireGoldenChi2", true, "Additional cut on the GoldenChi2"}; + Configurable minNCrossedRowsTPC{"minNCrossedRowsTPC", 70.f, "Additional cut on the minimum number of crossed rows in the TPC"}; + Configurable minNCrossedRowsOverFindableClustersTPC{"minNCrossedRowsOverFindableClustersTPC", 0.8f, "Additional cut on the minimum value of the ratio between crossed rows and findable clusters in the TPC"}; + Configurable maxChi2PerClusterTPC{"maxChi2PerClusterTPC", 4.f, "Additional cut on the maximum value of the chi2 per cluster in the TPC"}; + Configurable minChi2PerClusterTPC{"minChi2PerClusterTPC", 0.5f, "Additional cut on the minimum value of the chi2 per cluster in the TPC"}; + Configurable maxChi2PerClusterITS{"maxChi2PerClusterITS", 36.f, "Additional cut on the maximum value of the chi2 per cluster in the ITS"}; + Configurable maxDcaXYFactor{"maxDcaXYFactor", 1.f, "Additional cut on the maximum value of the DCA xy (multiplicative factor)"}; + Configurable maxDcaZ{"maxDcaZ", 2.f, "Additional cut on the maximum value of the DCA z"}; + Configurable minTPCNClsFound{"minTPCNClsFound", 100.f, "Additional cut on the minimum value of the number of found clusters in the TPC"}; + Configurable makeTHnSparseChoice{"makeTHnSparseChoice", false, "choose if produce thnsparse"}; // RD + Configurable enableTPCTOFvsEtaHistograms{"enableTPCTOFvsEtaHistograms", false, "choose if produce TPC tof vs Eta"}; + Configurable includeCentralityMC{"includeCentralityMC", false, "choose if include Centrality to MC"}; + Configurable enableTPCTOFVsMult{"enableTPCTOFVsMult", false, "Produce TPC-TOF plots vs multiplicity"}; + Configurable minITSnClusters{"minITSnClusters", 5, "minimum number of found ITS clusters"}; + + struct : ConfigurableGroup { + ConfigurableAxis binsPt{"binsPt", {VARIABLE_WIDTH, 0.0, 0.1, 0.12, 0.14, 0.16, 0.18, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 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}, "Binning of the pT axis"}; + ConfigurableAxis binsEta{"binsEta", {100, -1, 1}, "Binning of the eta axis"}; + ConfigurableAxis binsY{"binsY", {100, -1, 1}, "Binning of the y axis"}; + ConfigurableAxis binsnsigmaTPC{"binsnsigmaTPC", {200, -10, 10}, "Binning of the nsigmaTPC axis"}; + ConfigurableAxis binsnsigmaTOF{"binsnsigmaTOF", {200, -10, 10}, "Binning of the nsigmaTOF axis"}; + ConfigurableAxis binsDca{"binsDca", {VARIABLE_WIDTH, -3.0, -2.95, -2.9, -2.85, -2.8, -2.75, -2.7, -2.65, -2.6, -2.55, -2.5, -2.45, -2.4, -2.35, -2.3, -2.25, -2.2, -2.15, -2.1, -2.05, -2.0, -1.975, -1.95, -1.925, -1.9, -1.875, -1.85, -1.825, -1.8, -1.775, -1.75, -1.725, -1.7, -1.675, -1.65, -1.625, -1.6, -1.575, -1.55, -1.525, -1.5, -1.475, -1.45, -1.425, -1.4, -1.375, -1.35, -1.325, -1.3, -1.275, -1.25, -1.225, -1.2, -1.175, -1.15, -1.125, -1.1, -1.075, -1.05, -1.025, -1.0, -0.99, -0.98, -0.97, -0.96, -0.95, -0.94, -0.93, -0.92, -0.91, -0.9, -0.89, -0.88, -0.87, -0.86, -0.85, -0.84, -0.83, -0.82, -0.81, -0.8, -0.79, -0.78, -0.77, -0.76, -0.75, -0.74, -0.73, -0.72, -0.71, -0.7, -0.69, -0.68, -0.67, -0.66, -0.65, -0.64, -0.63, -0.62, -0.61, -0.6, -0.59, -0.58, -0.57, -0.56, -0.55, -0.54, -0.53, -0.52, -0.51, -0.5, -0.49, -0.48, -0.47, -0.46, -0.45, -0.44, -0.43, -0.42, -0.41, -0.4, -0.396, -0.392, -0.388, -0.384, -0.38, -0.376, -0.372, -0.368, -0.364, -0.36, -0.356, -0.352, -0.348, -0.344, -0.34, -0.336, -0.332, -0.328, -0.324, -0.32, -0.316, -0.312, -0.308, -0.304, -0.3, -0.296, -0.292, -0.288, -0.284, -0.28, -0.276, -0.272, -0.268, -0.264, -0.26, -0.256, -0.252, -0.248, -0.244, -0.24, -0.236, -0.232, -0.228, -0.224, -0.22, -0.216, -0.212, -0.208, -0.204, -0.2, -0.198, -0.196, -0.194, -0.192, -0.19, -0.188, -0.186, -0.184, -0.182, -0.18, -0.178, -0.176, -0.174, -0.172, -0.17, -0.168, -0.166, -0.164, -0.162, -0.16, -0.158, -0.156, -0.154, -0.152, -0.15, -0.148, -0.146, -0.144, -0.142, -0.14, -0.138, -0.136, -0.134, -0.132, -0.13, -0.128, -0.126, -0.124, -0.122, -0.12, -0.118, -0.116, -0.114, -0.112, -0.11, -0.108, -0.106, -0.104, -0.102, -0.1, -0.099, -0.098, -0.097, -0.096, -0.095, -0.094, -0.093, -0.092, -0.091, -0.09, -0.089, -0.088, -0.087, -0.086, -0.085, -0.084, -0.083, -0.082, -0.081, -0.08, -0.079, -0.078, -0.077, -0.076, -0.075, -0.074, -0.073, -0.072, -0.071, -0.07, -0.069, -0.068, -0.067, -0.066, -0.065, -0.064, -0.063, -0.062, -0.061, -0.06, -0.059, -0.058, -0.057, -0.056, -0.055, -0.054, -0.053, -0.052, -0.051, -0.05, -0.049, -0.048, -0.047, -0.046, -0.045, -0.044, -0.043, -0.042, -0.041, -0.04, -0.039, -0.038, -0.037, -0.036, -0.035, -0.034, -0.033, -0.032, -0.031, -0.03, -0.029, -0.028, -0.027, -0.026, -0.025, -0.024, -0.023, -0.022, -0.021, -0.02, -0.019, -0.018, -0.017, -0.016, -0.015, -0.014, -0.013, -0.012, -0.011, -0.01, -0.009, -0.008, -0.007, -0.006, -0.005, -0.004, -0.003, -0.002, -0.001, -0.0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01, 0.011, 0.012, 0.013, 0.014, 0.015, 0.016, 0.017, 0.018, 0.019, 0.02, 0.021, 0.022, 0.023, 0.024, 0.025, 0.026, 0.027, 0.028, 0.029, 0.03, 0.031, 0.032, 0.033, 0.034, 0.035, 0.036, 0.037, 0.038, 0.039, 0.04, 0.041, 0.042, 0.043, 0.044, 0.045, 0.046, 0.047, 0.048, 0.049, 0.05, 0.051, 0.052, 0.053, 0.054, 0.055, 0.056, 0.057, 0.058, 0.059, 0.06, 0.061, 0.062, 0.063, 0.064, 0.065, 0.066, 0.067, 0.068, 0.069, 0.07, 0.071, 0.072, 0.073, 0.074, 0.075, 0.076, 0.077, 0.078, 0.079, 0.08, 0.081, 0.082, 0.083, 0.084, 0.085, 0.086, 0.087, 0.088, 0.089, 0.09, 0.091, 0.092, 0.093, 0.094, 0.095, 0.096, 0.097, 0.098, 0.099, 0.1, 0.102, 0.104, 0.106, 0.108, 0.11, 0.112, 0.114, 0.116, 0.118, 0.12, 0.122, 0.124, 0.126, 0.128, 0.13, 0.132, 0.134, 0.136, 0.138, 0.14, 0.142, 0.144, 0.146, 0.148, 0.15, 0.152, 0.154, 0.156, 0.158, 0.16, 0.162, 0.164, 0.166, 0.168, 0.17, 0.172, 0.174, 0.176, 0.178, 0.18, 0.182, 0.184, 0.186, 0.188, 0.19, 0.192, 0.194, 0.196, 0.198, 0.2, 0.204, 0.208, 0.212, 0.216, 0.22, 0.224, 0.228, 0.232, 0.236, 0.24, 0.244, 0.248, 0.252, 0.256, 0.26, 0.264, 0.268, 0.272, 0.276, 0.28, 0.284, 0.288, 0.292, 0.296, 0.3, 0.304, 0.308, 0.312, 0.316, 0.32, 0.324, 0.328, 0.332, 0.336, 0.34, 0.344, 0.348, 0.352, 0.356, 0.36, 0.364, 0.368, 0.372, 0.376, 0.38, 0.384, 0.388, 0.392, 0.396, 0.4, 0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.6, 0.61, 0.62, 0.63, 0.64, 0.65, 0.66, 0.67, 0.68, 0.69, 0.7, 0.71, 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.8, 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88, 0.89, 0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.0, 1.025, 1.05, 1.075, 1.1, 1.125, 1.15, 1.175, 1.2, 1.225, 1.25, 1.275, 1.3, 1.325, 1.35, 1.375, 1.4, 1.425, 1.45, 1.475, 1.5, 1.525, 1.55, 1.575, 1.6, 1.625, 1.65, 1.675, 1.7, 1.725, 1.75, 1.775, 1.8, 1.825, 1.85, 1.875, 1.9, 1.925, 1.95, 1.975, 2.0, 2.05, 2.1, 2.15, 2.2, 2.25, 2.3, 2.35, 2.4, 2.45, 2.5, 2.55, 2.6, 2.65, 2.7, 2.75, 2.8, 2.85, 2.9, 2.95, 3.0}, "Binning of DCA xy and z axis"}; + ConfigurableAxis binsMultiplicity{"binsMultiplicity", {100, 0, 100}, "Binning for multiplicity"}; + ConfigurableAxis binsPercentile{"binsPercentile", {100, 0, 100}, "Binning for percentiles"}; + } binsOptions; + + Configurable enablePureDCAHistogram{"enablePureDCAHistogram", false, "Enables the pure DCA histograms"}; + Configurable enableTPCTOFHistograms{"enableTPCTOFHistograms", true, "Enables TPC TOF histograms"}; + Configurable enableDCAxyzHistograms{"enableDCAxyzHistograms", false, "Enables DCAxyz correlation histograms"}; + Configurable multiplicityEstimator{"multiplicityEstimator", 9, "Flag to use a multiplicity estimator: 0 no multiplicity, 1 MultFV0M, 2 MultFT0M, 3 MultFDDM, 4 MultTracklets, 5 MultTPC, 6 MultNTracksPV, 7 MultNTracksPVeta1, 8 CentralityFT0C, 9 CentralityFT0M, 10 CentralityFV0A"}; + Configurable enableTrackCutHistograms{"enableTrackCutHistograms", true, "Enables track cut histograms, before and after the cut"}; + Configurable doDCAxyCut{"doDCAxyCut", false, "do DCAxy Cut"}; + Configurable doDCAzCut{"doDCAzCut", false, "do DCAz Cut"}; + + // Histograms + HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + void init(o2::framework::InitContext&) + { + LOG(info) << "doprocessMC = " << doprocessMC; + LOG(info) << "doprocessStandard = " << doprocessStandard; + // Custom track cuts + if (useCustomTrackCuts.value) { + LOG(info) << "Using custom track cuts from values:"; + LOG(info) << "\trequireITS=" << requireITS.value; + LOG(info) << "\trequireTPC=" << requireTPC.value; + LOG(info) << "\trequireGoldenChi2=" << requireGoldenChi2.value; + LOG(info) << "\tmaxChi2PerClusterTPC=" << maxChi2PerClusterTPC.value; + LOG(info) << "\tminChi2PerClusterTPC=" << minChi2PerClusterTPC.value; + LOG(info) << "\tminNCrossedRowsTPC=" << minNCrossedRowsTPC.value; + LOG(info) << "\tminITSnClusters=" << minITSnClusters.value; + LOG(info) << "\tminTPCNClsFound=" << minTPCNClsFound.value; + LOG(info) << "\tmaxChi2PerClusterITS=" << maxChi2PerClusterITS.value; + LOG(info) << "\tmaxDcaZ=" << maxDcaZ.value; + LOG(info) << "\tmakeTHnSparseChoice=" << makeTHnSparseChoice.value; + + customTrackCuts = getGlobalTrackSelectionRun3ITSMatch(itsPattern.value); + LOG(info) << "Customizing track cuts:"; + customTrackCuts.SetRequireITSRefit(requireITS.value); + customTrackCuts.SetRequireTPCRefit(requireTPC.value); + customTrackCuts.SetMinNClustersITS(minITSnClusters.value); + customTrackCuts.SetRequireGoldenChi2(requireGoldenChi2.value); + customTrackCuts.SetMaxChi2PerClusterTPC(maxChi2PerClusterTPC.value); + customTrackCuts.SetMaxChi2PerClusterITS(maxChi2PerClusterITS.value); + customTrackCuts.SetMinNCrossedRowsTPC(minNCrossedRowsTPC.value); + customTrackCuts.SetMinNClustersTPC(minTPCNClsFound.value); + customTrackCuts.SetMinNCrossedRowsOverFindableClustersTPC(minNCrossedRowsOverFindableClustersTPC.value); + customTrackCuts.SetMaxDcaXYPtDep([](float /*pt*/) { return 10000.f; }); // No DCAxy cut will be used, this is done via the member function of the task + customTrackCuts.SetMaxDcaZ(maxDcaZ.value); + customTrackCuts.print(); + } + // Histograms + const AxisSpec vtxZAxis{100, -20, 20, "Vtx_{z} (cm)"}; + const AxisSpec pAxis{binsOptions.binsPt, "#it{p} (GeV/#it{c})"}; + const AxisSpec ptAxis{binsOptions.binsPt, "#it{p}_{T} (GeV/#it{c})"}; + const AxisSpec etaAxis{binsOptions.binsEta, "#eta"}; + const AxisSpec yAxis{binsOptions.binsY, "y"}; + AxisSpec multAxis{binsOptions.binsMultiplicity, "Undefined multiplicity estimator"}; + switch (multiplicityEstimator) { + case MultCodes::kNoMultiplicity: // No multiplicity + break; + case MultCodes::kMultFV0M: // MultFV0M + multAxis.name = "MultFV0M"; + break; + case MultCodes::kMultFT0M: // MultFT0M + multAxis.name = "MultFT0M"; + break; + case MultCodes::kMultFDDM: // MultFDDM + multAxis.name = "MultFDDM"; + break; + case MultCodes::kMultTracklets: // MultTracklets + multAxis.name = "MultTracklets"; + break; + case MultCodes::kMultTPC: // MultTPC + multAxis.name = "MultTPC"; + break; + case MultCodes::kMultNTracksPV: // MultNTracksPV + multAxis.name = "MultNTracksPV"; + break; + case MultCodes::kMultNTracksPVeta1: // MultNTracksPVeta1 + multAxis.name = "MultNTracksPVeta1"; + break; + case MultCodes::kCentralityFT0C: // Centrality FT0C + multAxis = {binsOptions.binsPercentile, "Centrality FT0C"}; + break; + case MultCodes::kCentralityFT0M: // Centrality FT0M + multAxis = {binsOptions.binsPercentile, "Centrality FT0M"}; + break; + case MultCodes::kCentralityFV0A: // Centrality FV0A + multAxis = {binsOptions.binsPercentile, "Centrality FV0A"}; + break; + default: + LOG(fatal) << "Unrecognized option for multiplicity " << multiplicityEstimator; + } + histos.add("event/vertexz", "", HistType::kTH1D, {vtxZAxis}); + + enum EvSelBin { + kAllCollisions = 1, + kTriggerTVXCut, + kNoTFBorderCut, + kNoITSROFBorderCut, + kPosZCut, + kINELGt0Cut, + kXimXimb + }; + auto h = histos.add("Evsel", "Evsel", HistType::kTH1D, {{7, 0.5, 7.5}}); + h->GetXaxis()->SetBinLabel(kAllCollisions, "All collisions"); + h->GetXaxis()->SetBinLabel(kTriggerTVXCut, "Trigger TVX cut"); + h->GetXaxis()->SetBinLabel(kNoTFBorderCut, "kNoTFBorder cut"); + h->GetXaxis()->SetBinLabel(kNoITSROFBorderCut, "kNoITSROFBorder cut"); + h->GetXaxis()->SetBinLabel(kPosZCut, "posZ cut"); + h->GetXaxis()->SetBinLabel(kINELGt0Cut, "INEL>0 cut"); + h->GetXaxis()->SetBinLabel(kXimXimb, "With at least a #Xi^{-}-#bar{#Xi^{-}} pair"); + + h = histos.add("tracksel", "tracksel", HistType::kTH1D, {{10, 0.5, 10.5}}); + h->GetXaxis()->SetBinLabel(1, "Tracks read"); + h->GetXaxis()->SetBinLabel(2, Form(" %.2f < #eta < %.2f ", trkselOptions.cfgCutEtaMin.value, trkselOptions.cfgCutEtaMax.value)); + h->GetXaxis()->SetBinLabel(3, "Quality passed"); + h->GetXaxis()->SetBinLabel(4, "TOF passed (partial)"); + + histos.add("Centrality/FT0M", "FT0M", HistType::kTH1D, {{binsOptions.binsPercentile, "Centrality FT0M"}}); + histos.add("Mult/FT0M", "MultFT0M", HistType::kTH1D, {{binsOptions.binsMultiplicity, "MultFT0M"}}); + + histos.add("Mult/TPC", "MultTPC", HistType::kTH1D, {{binsOptions.binsMultiplicity, "MultTPC"}}); + histos.add("Mult/NTracksPV", "MultNTracksPV", HistType::kTH1D, {{binsOptions.binsMultiplicity, "MultNTracksPV"}}); + histos.add("Mult/NTracksPVeta1", "MultNTracksPVeta1", HistType::kTH1D, {{binsOptions.binsMultiplicity, "MultNTracksPVeta1"}}); + + const AxisSpec dcaXyAxis{binsOptions.binsDca, "DCA_{xy} (cm)"}; + const AxisSpec dcaZAxis{binsOptions.binsDca, "DCA_{z} (cm)"}; + + if (enableTrackCutHistograms) { + const AxisSpec chargeAxis{2, -2.f, 2.f, "Charge"}; + histos.add("track/pos/Eta", "Eta Positive tracks", HistType::kTH1D, {{binsOptions.binsEta, "#eta tracks"}}); + histos.add("track/neg/Eta", "Eta Negative tracks", HistType::kTH1D, {{binsOptions.binsEta, "#eta tracks"}}); + // its histograms + histos.add("track/ITS/itsNCls", "number of found ITS clusters;# clusters ITS", kTH2D, {{8, -0.5, 7.5}, chargeAxis}); + histos.add("track/ITS/itsChi2NCl", "chi2 per ITS cluster;chi2 / cluster ITS", kTH2D, {{100, 0, 40}, chargeAxis}); + // tpc histograms + histos.add("track/TPC/tpcNClsFindable", "number of findable TPC clusters;# findable clusters TPC", kTH2D, {{165, -0.5, 164.5}, chargeAxis}); + histos.add("track/TPC/tpcNClsFound", "number of found TPC clusters;# clusters TPC", kTH2D, {{165, -0.5, 164.5}, chargeAxis}); + histos.add("track/TPC/tpcNClsShared", "number of shared TPC clusters;# shared clusters TPC", kTH2D, {{165, -0.5, 164.5}, chargeAxis}); + histos.add("track/TPC/tpcCrossedRows", "number of crossed TPC rows;# crossed rows TPC", kTH2D, {{165, -0.5, 164.5}, chargeAxis}); + histos.add("track/TPC/tpcFractionSharedCls", "fraction of shared TPC clusters;fraction shared clusters TPC", kTH2D, {{100, 0., 1.}, chargeAxis}); + histos.add("track/TPC/tpcCrossedRowsOverFindableCls", "crossed TPC rows over findable clusters;crossed rows / findable clusters TPC", kTH2D, {{60, 0.7, 1.3}, chargeAxis}); + histos.add("track/TPC/tpcChi2NCl", "chi2 per cluster in TPC;chi2 / cluster TPC", kTH2D, {{100, 0, 10}, chargeAxis}); + histos.add("Vertex/histGenVtxMC", "MC generated vertex z position", HistType::kTH1F, {{400, -40., +40., "z position (cm)"}}); + + histos.addClone("track/ITS/itsNCls", "track/selected/ITS/itsNCls"); + histos.addClone("track/ITS/itsChi2NCl", "track/selected/ITS/itsChi2NCl"); + histos.addClone("track/TPC/tpcNClsFindable", "track/selected/TPC/tpcNClsFindable"); + histos.addClone("track/TPC/tpcNClsFound", "track/selected/TPC/tpcNClsFound"); + histos.addClone("track/TPC/tpcNClsShared", "track/selected/TPC/tpcNClsShared"); + histos.addClone("track/TPC/tpcCrossedRows", "track/selected/TPC/tpcCrossedRows"); + histos.addClone("track/TPC/tpcFractionSharedCls", "track/selected/TPC/tpcFractionSharedCls"); + histos.addClone("track/TPC/tpcCrossedRowsOverFindableCls", "track/selected/TPC/tpcCrossedRowsOverFindableCls"); + histos.addClone("track/TPC/tpcChi2NCl", "track/selected/TPC/tpcChi2NCl"); + } + + // 2 detectors + histos.add("Data/pos/pt/its_tpc", "pos ITS-TPC", kTH1D, {ptAxis}); + histos.add("Data/neg/pt/its_tpc", "neg ITS-TPC", kTH1D, {ptAxis}); + + histos.add("Data/pos/pt/tpc_tof", "pos TPC-TOF", kTH1D, {ptAxis}); + histos.add("Data/neg/pt/tpc_tof", "neg TPC-TOF", kTH1D, {ptAxis}); + + histos.add("Data/pos/pt/its_tof", "pos ITS-TOF", kTH1D, {ptAxis}); + histos.add("Data/neg/pt/its_tof", "neg ITS-TOF", kTH1D, {ptAxis}); + + // 1 detectors + histos.add("Data/pos/pt/its", "pos ITS", kTH1D, {ptAxis}); + histos.add("Data/neg/pt/its", "neg ITS", kTH1D, {ptAxis}); + histos.add("Data/pos/pt/tpc", "pos TPC", kTH1D, {ptAxis}); + histos.add("Data/neg/pt/tpc", "neg TPC", kTH1D, {ptAxis}); + if (doprocessMC) { + histos.add("MC/fake/pos", "Fake positive tracks", kTH1D, {ptAxis}); + histos.add("MC/fake/neg", "Fake negative tracks", kTH1D, {ptAxis}); + histos.add("MC/no_collision/pos", "No collision pos track", kTH1D, {ptAxis}); + histos.add("MC/no_collision/neg", "No collision neg track", kTH1D, {ptAxis}); + + auto hh = histos.add("MC/GenRecoCollisions", "Generated and Reconstructed MC Collisions", kTH1D, {{10, 0.5, 10.5}}); + hh->GetXaxis()->SetBinLabel(1, "Collisions generated"); + hh->GetXaxis()->SetBinLabel(2, "Collisions reconstructed"); + hh->GetXaxis()->SetBinLabel(3, "INEL>0"); + hh->GetXaxis()->SetBinLabel(4, "INEL>1"); + hh->GetXaxis()->SetBinLabel(5, "hasParticleInFT0C && hasParticleInFT0A"); + histos.add("MC/MultiplicityRecoEv", "MC multiplicity", kTH1D, {multAxis}); + histos.add("MC/Multiplicity", "MC multiplicity", kTH1D, {multAxis}); + histos.add("MC/MultiplicityMCINELgt0", "MC multiplicity", kTH1D, {multAxis}); + } + for (int i = 0; i < NpCharge; i++) { + const AxisSpec nsigmaTPCAxis{binsOptions.binsnsigmaTPC, Form("N_{#sigma}^{TPC}(%s)", pTCharge[i])}; + const AxisSpec nsigmaTOFAxis{binsOptions.binsnsigmaTOF, Form("N_{#sigma}^{TOF}(%s)", pTCharge[i])}; + + if (multiplicityEstimator == MultCodes::kNoMultiplicity) { + histos.add(hnsigmatof[i].data(), pTCharge[i], kTH2D, {ptAxis, nsigmaTOFAxis}); + histos.add(hnsigmatpc[i].data(), pTCharge[i], kTH2D, {ptAxis, nsigmaTPCAxis}); + if (enableTPCTOFHistograms) { + if (enableTPCTOFvsEtaHistograms) { + histos.add(hnsigmatpctof[i].data(), pTCharge[i], kTHnSparseD, {ptAxis, etaAxis, nsigmaTPCAxis, nsigmaTOFAxis}); + } else { + histos.add(hnsigmatpctof[i].data(), pTCharge[i], kTH3D, {ptAxis, nsigmaTPCAxis, nsigmaTOFAxis}); + } + } + } else { + if (makeTHnSparseChoice) { // RD + histos.add(hnsigmatof[i].data(), pTCharge[i], kTHnSparseD, {ptAxis, nsigmaTOFAxis, multAxis, dcaXyAxis}); // RD + histos.add(hnsigmatpc[i].data(), pTCharge[i], kTHnSparseD, {ptAxis, nsigmaTPCAxis, multAxis, dcaXyAxis}); // RD + } else { + histos.add(hnsigmatof[i].data(), pTCharge[i], kTH3D, {ptAxis, nsigmaTOFAxis, multAxis}); + histos.add(hnsigmatpc[i].data(), pTCharge[i], kTH3D, {ptAxis, nsigmaTPCAxis, multAxis}); + } + if (enableTPCTOFHistograms) { + if (enableTPCTOFVsMult) { + if (enableTPCTOFvsEtaHistograms) { + histos.add(hnsigmatpctof[i].data(), pTCharge[i], kTHnSparseD, {ptAxis, etaAxis, nsigmaTPCAxis, nsigmaTOFAxis, multAxis}); + } else { + histos.add(hnsigmatpctof[i].data(), pTCharge[i], kTHnSparseD, {ptAxis, nsigmaTPCAxis, nsigmaTOFAxis, multAxis}); + } + } else { + if (enableTPCTOFvsEtaHistograms) { + histos.add(hnsigmatpctof[i].data(), pTCharge[i], kTHnSparseD, {ptAxis, etaAxis, nsigmaTPCAxis, nsigmaTOFAxis}); + } else { + histos.add(hnsigmatpctof[i].data(), pTCharge[i], kTH3D, {ptAxis, nsigmaTPCAxis, nsigmaTOFAxis}); + } + } + } + } + + if (enableDCAxyzHistograms) { + hDcaXYZ[i] = histos.add(Form("dca/%s/%s", (i < Np) ? "pos" : "neg", pN[i % Np]), pTCharge[i], kTH3D, {ptAxis, dcaXyAxis, dcaZAxis}); + } else { + histos.add(hdcaxy[i].data(), pTCharge[i], kTH2D, {ptAxis, dcaXyAxis}); + histos.add(hdcaz[i].data(), pTCharge[i], kTH2D, {ptAxis, dcaZAxis}); + } + // --- MC --- + if (doprocessMC) { + const std::string cpName = Form("/%s/%s", (i < Np) ? "pos" : "neg", pN[i % Np]); + if (includeCentralityMC) { + histos.add(hpt_num_prm[i].data(), pTCharge[i], kTHnSparseD, {ptAxis, multAxis, dcaXyAxis}); + histos.add(hpt_num_str[i].data(), pTCharge[i], kTHnSparseD, {ptAxis, multAxis, dcaXyAxis}); + histos.add(hpt_num_mat[i].data(), pTCharge[i], kTHnSparseD, {ptAxis, multAxis, dcaXyAxis}); + + histos.add(hpt_den_prm[i].data(), pTCharge[i], kTHnSparseD, {ptAxis, multAxis}); + histos.add(hpt_den_str[i].data(), pTCharge[i], kTHnSparseD, {ptAxis, multAxis}); + histos.add(hpt_den_mat[i].data(), pTCharge[i], kTHnSparseD, {ptAxis, multAxis}); + + histos.add(hpt_numtof_prm[i].data(), pTCharge[i], kTHnSparseD, {ptAxis, multAxis, dcaXyAxis}); + histos.add(hpt_numtof_str[i].data(), pTCharge[i], kTHnSparseD, {ptAxis, multAxis, dcaXyAxis}); + histos.add(hpt_numtof_mat[i].data(), pTCharge[i], kTHnSparseD, {ptAxis, multAxis, dcaXyAxis}); + + histos.add(hpt_den_prm_goodev[i].data(), pTCharge[i], kTH3D, {ptAxis, multAxis, etaAxis}); + + } else { + histos.add(hpt_num_prm[i].data(), pTCharge[i], kTH2D, {ptAxis, dcaXyAxis}); + histos.add(hpt_num_str[i].data(), pTCharge[i], kTH2D, {ptAxis, dcaXyAxis}); + histos.add(hpt_num_mat[i].data(), pTCharge[i], kTH2D, {ptAxis, dcaXyAxis}); + + histos.add(hpt_numtof_prm[i].data(), pTCharge[i], kTH2D, {ptAxis, dcaXyAxis}); + histos.add(hpt_numtof_str[i].data(), pTCharge[i], kTH2D, {ptAxis, dcaXyAxis}); + histos.add(hpt_numtof_mat[i].data(), pTCharge[i], kTH2D, {ptAxis, dcaXyAxis}); + + histos.add(hpt_den_prm[i].data(), pTCharge[i], kTH1D, {ptAxis}); + histos.add(hpt_den_str[i].data(), pTCharge[i], kTH1D, {ptAxis}); + histos.add(hpt_den_mat[i].data(), pTCharge[i], kTH1D, {ptAxis}); + + histos.add(hpt_den_prm_goodev[i].data(), pTCharge[i], kTH2D, {ptAxis, etaAxis}); + } + histos.add(hpt_den_prm_mcgoodev[i].data(), pTCharge[i], kTH2D, {ptAxis, multAxis}); + if (enableDCAxyzHistograms) { + hDcaXYZPrm[i] = histos.add("dcaprm" + cpName, pTCharge[i], kTH3D, {ptAxis, dcaXyAxis, dcaZAxis}); + hDcaXYZStr[i] = histos.add("dcastr" + cpName, pTCharge[i], kTH3D, {ptAxis, dcaXyAxis, dcaZAxis}); + hDcaXYZMat[i] = histos.add("dcamat" + cpName, pTCharge[i], kTH3D, {ptAxis, dcaXyAxis, dcaZAxis}); + } else { + histos.add(hdcaxyprm[i].data(), pTCharge[i], kTH2D, {ptAxis, dcaXyAxis}); + histos.add(hdcazprm[i].data(), pTCharge[i], kTH2D, {ptAxis, dcaZAxis}); + histos.add(hdcaxystr[i].data(), pTCharge[i], kTH2D, {ptAxis, dcaXyAxis}); + histos.add(hdcazstr[i].data(), pTCharge[i], kTH2D, {ptAxis, dcaZAxis}); + histos.add(hdcaxymat[i].data(), pTCharge[i], kTH2D, {ptAxis, dcaXyAxis}); + histos.add(hdcazmat[i].data(), pTCharge[i], kTH2D, {ptAxis, dcaZAxis}); + } + } + } + // Print output histograms statistics + LOG(info) << "Size of the histograms in spectraTOF"; + histos.print(); + } + + template + void fillParticleHistos(const T& track, const C& collision) + { + if (std::abs(track.rapidity(PID::getMass(id))) > trkselOptions.cfgCutY) { + return; + } + if constexpr (id == PID::Kaon) { + if (kaonIsPvContrib && !track.isPVContributor()) { + return; + } + } + + const auto& nsigmaTOF = o2::aod::pidutils::tofNSigma(track); + const auto& nsigmaTPC = o2::aod::pidutils::tpcNSigma(track); + + // Filling DCA info with the TPC+TOF PID + bool isDCAPureSample = (std::sqrt(nsigmaTOF * nsigmaTOF + nsigmaTPC * nsigmaTPC) < DcaMaxCombinedSigma); + if (track.pt() <= DcaTrkPtCut) { + isDCAPureSample = (fabs(nsigmaTPC) < DcaMaxTPCSigma); + } + if (isDCAPureSample) { + if (enableDCAxyzHistograms) { + if (track.sign() > 0) { + hDcaXYZ[id]->Fill(track.pt(), track.dcaXY(), track.dcaZ()); + } else { + hDcaXYZ[id + Np]->Fill(track.pt(), track.dcaXY(), track.dcaZ()); + } + } else { + if (track.sign() > 0) { + histos.fill(HIST(hdcaxy[id]), track.pt(), track.dcaXY()); + histos.fill(HIST(hdcaz[id]), track.pt(), track.dcaZ()); + } else { + histos.fill(HIST(hdcaxy[id + Np]), track.pt(), track.dcaXY()); + histos.fill(HIST(hdcaz[id + Np]), track.pt(), track.dcaZ()); + } + } + } + + // -- nsigma histograms -- + if (doDCAxyCut && !passesDCAxyCut(track)) { + return; + } + if (doDCAzCut && !passesDCAzCut(track)) { + return; + } + const float multiplicity = getMultiplicity(collision); + if (makeTHnSparseChoice) { // RD + if (track.sign() > 0) { // RD + histos.fill(HIST(hnsigmatpc[id]), track.pt(), nsigmaTPC, multiplicity, track.dcaXY()); // RD + } else { // RD + histos.fill(HIST(hnsigmatpc[id + Np]), track.pt(), nsigmaTPC, multiplicity, track.dcaXY()); // RD + } + } else { + if (track.sign() > 0) { + histos.fill(HIST(hnsigmatpc[id]), track.pt(), nsigmaTPC, multiplicity); + } else { + histos.fill(HIST(hnsigmatpc[id + Np]), track.pt(), nsigmaTPC, multiplicity); + } + } + // TOF part + if (!track.hasTOF()) { + return; + } + + if (makeTHnSparseChoice) { // RD + if (track.sign() > 0) { // RD + histos.fill(HIST(hnsigmatof[id]), track.pt(), nsigmaTOF, multiplicity, track.dcaXY()); // RD + } else { // RD + histos.fill(HIST(hnsigmatof[id + Np]), track.pt(), nsigmaTOF, multiplicity, track.dcaXY()); // RD + } + } else { + if (track.sign() > 0) { + histos.fill(HIST(hnsigmatof[id]), track.pt(), nsigmaTOF, multiplicity); + } else { + histos.fill(HIST(hnsigmatof[id + Np]), track.pt(), nsigmaTOF, multiplicity); + } + } + + if (enableTPCTOFHistograms) { + if (enableTPCTOFVsMult) { + if (enableTPCTOFvsEtaHistograms) { + if (track.sign() > 0) { + histos.fill(HIST(hnsigmatpctof[id]), track.pt(), track.eta(), nsigmaTPC, nsigmaTOF, multiplicity); + } else { + histos.fill(HIST(hnsigmatpctof[id + Np]), track.pt(), track.eta(), nsigmaTPC, nsigmaTOF, multiplicity); + } + } else { + if (track.sign() > 0) { + histos.fill(HIST(hnsigmatpctof[id]), track.pt(), nsigmaTPC, nsigmaTOF, multiplicity); + } else { + histos.fill(HIST(hnsigmatpctof[id + Np]), track.pt(), nsigmaTPC, nsigmaTOF, multiplicity); + } + } + } else { + if (enableTPCTOFvsEtaHistograms) { + if (track.sign() > 0) { + histos.fill(HIST(hnsigmatpctof[id]), track.pt(), track.eta(), nsigmaTPC, nsigmaTOF); + } else { + histos.fill(HIST(hnsigmatpctof[id + Np]), track.pt(), track.eta(), nsigmaTPC, nsigmaTOF); + } + } else { + if (track.sign() > 0) { + histos.fill(HIST(hnsigmatpctof[id]), track.pt(), nsigmaTPC, nsigmaTOF); + } else { + histos.fill(HIST(hnsigmatpctof[id + Np]), track.pt(), nsigmaTPC, nsigmaTOF); + } + } + } + } + + } // RD + + template + bool isEventSelected(CollisionType const& collision) + { + if constexpr (fillHistograms) { + histos.fill(HIST("Evsel"), 1.f); + } + + if (evselOptions.requireSel8) { + if (!collision.sel8()) { + return false; + } + histos.fill(HIST("Evsel"), 2.f); + histos.fill(HIST("Evsel"), 3.f); + histos.fill(HIST("Evsel"), 4.f); + } else { + if (!collision.selection_bit(aod::evsel::kIsTriggerTVX)) { + return false; + } + histos.fill(HIST("Evsel"), 2.f); + if (evselOptions.removeNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) { + return false; + } + histos.fill(HIST("Evsel"), 3.f); + if (evselOptions.removeITSROFrameBorder && !collision.selection_bit(aod::evsel::kNoITSROFrameBorder)) { + return false; + } + histos.fill(HIST("Evsel"), 4.f); + } + if (std::abs(collision.posZ()) > evselOptions.cfgCutVertex) { + return false; + } + if constexpr (fillHistograms) { + histos.fill(HIST("Evsel"), 5.f); + } + if (evselOptions.cfgINELCut == EvSelInelGt0Cut && !collision.isInelGt0()) { + return false; + } + if constexpr (fillHistograms) { + histos.fill(HIST("Evsel"), 6.f); + histos.fill(HIST("event/vertexz"), collision.posZ()); + } + + if constexpr (fillMultiplicity) { + histos.fill(HIST("Centrality/FT0M"), collision.centFT0M()); + histos.fill(HIST("Mult/FT0M"), collision.multZeqFT0A() + collision.multZeqFT0C()); + histos.fill(HIST("Mult/TPC"), collision.multTPC()); + histos.fill(HIST("Mult/NTracksPV"), collision.multZeqNTracksPV()); + histos.fill(HIST("Mult/NTracksPVeta1"), collision.multNTracksPVeta1()); + } + + return true; + } + template + bool passesDCAxyCut(TrackType const& track) const + { + if (useCustomTrackCuts.value) { + for (int i = 0; i < static_cast(TrackSelection::TrackCuts::kNCuts); i++) { + if (i == static_cast(TrackSelection::TrackCuts::kDCAxy)) { + continue; + } + if (!customTrackCuts.IsSelected(track, static_cast(i))) { + return false; + } + } + return (std::abs(track.dcaXY()) <= (maxDcaXYFactor.value * (0.004f + 0.013f / track.pt()))); + } + + return track.isGlobalTrack(); + } + + template + bool passesDCAzCut(TrackType const& track) const + { + if (useCustomTrackCuts.value) { + for (int i = 0; i < static_cast(TrackSelection::TrackCuts::kNCuts); i++) { + if (i == static_cast(TrackSelection::TrackCuts::kDCAz)) { + continue; + } + if (!customTrackCuts.IsSelected(track, static_cast(i))) { + return false; + } + } + // return (std::abs(track.dcaZ()) <= maxDcaZ.value); + return true; + } + + return track.isGlobalTrack(); + } + + template + bool passesCutWoDCA(TrackType const& track) const + { + if (useCustomTrackCuts.value) { + for (int i = 0; i < static_cast(TrackSelection::TrackCuts::kNCuts); i++) { + if (i == static_cast(TrackSelection::TrackCuts::kDCAxy)) { + continue; + } + if (i == static_cast(TrackSelection::TrackCuts::kDCAz)) { + continue; + } + if (!customTrackCuts.IsSelected(track, static_cast(i))) { + return false; + } + } + return true; + } + return track.isGlobalTrackWoDCA(); + } + template + bool isTrackSelected(TrackType const& track, CollisionType const& /*collision*/) + { + if constexpr (fillHistograms) { + histos.fill(HIST("tracksel"), 1); + } + if (track.eta() < trkselOptions.cfgCutEtaMin || track.eta() > trkselOptions.cfgCutEtaMax) { + return false; + } + + if constexpr (fillHistograms) { + histos.fill(HIST("tracksel"), 2); + if (enableTrackCutHistograms) { + if (track.sign() > 0) { + histos.fill(HIST("track/pos/Eta"), track.eta()); + } else { + histos.fill(HIST("track/neg/Eta"), track.eta()); + } + + if (track.hasITS() && track.hasTPC()) { + histos.fill(HIST("track/ITS/itsNCls"), track.itsNCls(), track.sign()); + histos.fill(HIST("track/ITS/itsChi2NCl"), track.itsChi2NCl(), track.sign()); + + histos.fill(HIST("track/TPC/tpcNClsFindable"), track.tpcNClsFindable(), track.sign()); + histos.fill(HIST("track/TPC/tpcNClsFound"), track.tpcNClsFound(), track.sign()); + histos.fill(HIST("track/TPC/tpcNClsShared"), track.tpcNClsShared(), track.sign()); + histos.fill(HIST("track/TPC/tpcCrossedRows"), track.tpcNClsCrossedRows(), track.sign()); + histos.fill(HIST("track/TPC/tpcCrossedRowsOverFindableCls"), track.tpcCrossedRowsOverFindableCls(), track.sign()); + histos.fill(HIST("track/TPC/tpcFractionSharedCls"), track.tpcFractionSharedCls(), track.sign()); + histos.fill(HIST("track/TPC/tpcChi2NCl"), track.tpcChi2NCl(), track.sign()); + } + } + + if (track.hasITS() && track.isQualityTrackITS() && track.isInAcceptanceTrack()) { + if (track.sign() > 0) { + histos.fill(HIST("Data/pos/pt/its"), track.pt()); + } else { + histos.fill(HIST("Data/neg/pt/its"), track.pt()); + } + } + if (track.hasTPC() && track.isQualityTrackTPC() && track.isInAcceptanceTrack()) { + if (track.sign() > 0) { + histos.fill(HIST("Data/pos/pt/tpc"), track.pt()); + } else { + histos.fill(HIST("Data/neg/pt/tpc"), track.pt()); + } + } + } + + if (track.tpcChi2NCl() < minChi2PerClusterTPC || track.tpcChi2NCl() > maxChi2PerClusterTPC) { + return false; + } + + if (!passesCutWoDCA(track)) { + return false; + } + + if constexpr (fillHistograms) { + histos.fill(HIST("tracksel"), 3); + if (track.hasTOF()) { + histos.fill(HIST("tracksel"), 4); + } + if (enableTrackCutHistograms) { + histos.fill(HIST("track/selected/ITS/itsNCls"), track.itsNCls(), track.sign()); + histos.fill(HIST("track/selected/ITS/itsChi2NCl"), track.itsChi2NCl(), track.sign()); + + histos.fill(HIST("track/selected/TPC/tpcNClsFindable"), track.tpcNClsFindable(), track.sign()); + histos.fill(HIST("track/selected/TPC/tpcNClsFound"), track.tpcNClsFound(), track.sign()); + histos.fill(HIST("track/selected/TPC/tpcNClsShared"), track.tpcNClsShared(), track.sign()); + histos.fill(HIST("track/selected/TPC/tpcCrossedRows"), track.tpcNClsCrossedRows(), track.sign()); + histos.fill(HIST("track/selected/TPC/tpcCrossedRowsOverFindableCls"), track.tpcCrossedRowsOverFindableCls(), track.sign()); + histos.fill(HIST("track/selected/TPC/tpcFractionSharedCls"), track.tpcFractionSharedCls(), track.sign()); + histos.fill(HIST("track/selected/TPC/tpcChi2NCl"), track.tpcChi2NCl(), track.sign()); + } + } + + if constexpr (fillHistograms) { + + if (track.hasITS() && track.hasTPC()) { + if (track.sign() > 0) { + histos.fill(HIST("Data/pos/pt/its_tpc"), track.pt()); + } else { + histos.fill(HIST("Data/neg/pt/its_tpc"), track.pt()); + } + } + if (track.hasTPC() && track.hasTOF()) { + if (track.sign() > 0) { + histos.fill(HIST("Data/pos/pt/tpc_tof"), track.pt()); + } else { + histos.fill(HIST("Data/neg/pt/tpc_tof"), track.pt()); + } + } + if (track.hasITS() && track.hasTOF()) { + if (track.sign() > 0) { + histos.fill(HIST("Data/pos/pt/its_tof"), track.pt()); + } else { + histos.fill(HIST("Data/neg/pt/its_tof"), track.pt()); + } + } + } + + return true; + } + + using CollisionCandidate = soa::Join; + using TrackCandidates = soa::Join; + + void processStandard(CollisionCandidate::iterator const& collision, + soa::Join const& tracks, + aod::BCs const&) + { + if (!isEventSelected(collision)) { + return; + } + + for (const auto& track : tracks) { + if (!isTrackSelected(track, collision)) { + continue; + } + fillParticleHistos(track, collision); + fillParticleHistos(track, collision); + // fillParticleHistos(track, collision); + } + + } // end of the process function + PROCESS_SWITCH(tofSpectra, processStandard, "Standard data processor", true); + + template + float getMultiplicity(const CollisionType& collision) + { + switch (multiplicityEstimator) { + case MultCodes::kNoMultiplicity: // No multiplicity + return 50.f; // to check if its filled + break; + case MultCodes::kMultFV0M: // MultFV0M + if constexpr (!isMC) { + return collision.multZeqFV0A(); + } else { + return 50.f; // Not implemented yet + } + break; + case MultCodes::kMultFT0M: + if constexpr (!isMC) { + return collision.multZeqFT0A() + collision.multZeqFT0C(); + } else { + return 50.f; // Not implemented yet + } + break; + case MultCodes::kMultFDDM: // MultFDDM + if constexpr (!isMC) { + return collision.multZeqFDDA() + collision.multZeqFDDC(); + } else { + return 50.f; // Not implemented yet + } + break; + case MultCodes::kMultTracklets: // MultTracklets + if constexpr (!isMC) { + // return collision.multTracklets(); + } else { + return 50.f; // Not implemented yet + } + return 0.f; // Undefined in Run3 + break; + case MultCodes::kMultTPC: // MultTPC + if constexpr (!isMC) { + return collision.multTPC(); + } else { + return 50.f; // Not implemented yet + } + break; + case MultCodes::kMultNTracksPV: // MultNTracksPV + if constexpr (!isMC) { + // return collision.multNTracksPV(); + return collision.multZeqNTracksPV(); + } else { + return 50.f; // Not implemented yet + } + break; + case MultCodes::kMultNTracksPVeta1: // MultNTracksPVeta1 + if constexpr (!isMC) { + return collision.multNTracksPVeta1(); + } else { + return 50.f; // Not implemented yet + } + break; + case MultCodes::kCentralityFT0C: // Centrality FT0C + if constexpr (!isMC) { + return collision.centFT0C(); + } else { + return 50.f; // Not implemented yet + } + break; + case MultCodes::kCentralityFT0M: // Centrality FT0M + return collision.centFT0M(); + break; + default: + LOG(fatal) << "Unknown multiplicity estimator: " << multiplicityEstimator; + return 0.f; + } + } + using GenMCCollisions = soa::Join; + float getMultiplicityMC(const GenMCCollisions::iterator& collision) { return getMultiplicity(collision); } + + template + bool isParticleEnabled() + { + if constexpr (id == PID::Pion || id == Np + PID::Pion) { + return true; + } else if constexpr (id == PID::Kaon || id == Np + PID::Kaon) { + return true; + } + return false; + } + + using RecoMCCollisions = soa::Join; // RD + template + void fillTrackHistogramsMC(TrackType const& track, + ParticleType::iterator const& mcParticle, + RecoMCCollisions::iterator const& collision, + ParticleType const& /*mcParticles*/) + { + if (!isParticleEnabled()) { // Check if the particle is enabled + return; + } + if (!collision.has_mcCollision()) { + return; // Skips processing if no corresponding MC collision is found (rare case!) + } + + const float multiplicity = getMultiplicity(collision); + + if (mcParticle.pdgCode() != PDGs[i]) { + return; + } + + if (track.eta() < trkselOptions.cfgCutEtaMin || track.eta() > trkselOptions.cfgCutEtaMax) { + return; + } + + if (std::abs(mcParticle.y()) > trkselOptions.cfgCutY) { + return; + } + if (enablePureDCAHistogram) { + float nsigmaTPC = 999.f; + float nsigmaTOF = 999.f; + + if constexpr (i == PID::Pion || i == Np + PID::Pion) { + nsigmaTPC = o2::aod::pidutils::tpcNSigma<2>(track); + nsigmaTOF = o2::aod::pidutils::tofNSigma<2>(track); + } else if constexpr (i == PID::Kaon || i == Np + PID::Kaon) { + nsigmaTPC = o2::aod::pidutils::tpcNSigma<3>(track); + nsigmaTOF = o2::aod::pidutils::tofNSigma<3>(track); + } + + // Filling DCA info with the TPC+TOF PID + // -- FEED DOWN : DCAxy distribution cor prm, str e mat per DCA Pure Sample. NO MULT per ora -- + bool isDCAPureSample = (std::sqrt(nsigmaTOF * nsigmaTOF + nsigmaTPC * nsigmaTPC) < DcaMaxCombinedSigma); + if (track.pt() <= DcaTrkPtCut) { + isDCAPureSample = (nsigmaTPC < DcaMaxTPCSigma); + } + + if (isDCAPureSample) { + if (!mcParticle.isPhysicalPrimary()) { // Secondaries (weak decays and material) + if (mcParticle.getProcess() == kPDecay) { // Particles from decay + if (enableDCAxyzHistograms) { + hDcaXYZStr[i]->Fill(track.pt(), track.dcaXY(), track.dcaZ()); + } else { + histos.fill(HIST(hdcaxystr[i]), track.pt(), track.dcaXY()); + histos.fill(HIST(hdcazstr[i]), track.pt(), track.dcaZ()); + } + } else { // Particles from the material + if (enableDCAxyzHistograms) { + hDcaXYZMat[i]->Fill(track.pt(), track.dcaXY(), track.dcaZ()); + } else { + histos.fill(HIST(hdcaxymat[i]), track.pt(), track.dcaXY()); + histos.fill(HIST(hdcazmat[i]), track.pt(), track.dcaZ()); + } + } + } else { // Primaries + if (enableDCAxyzHistograms) { + hDcaXYZPrm[i]->Fill(track.pt(), track.dcaXY(), track.dcaZ()); + } else { + histos.fill(HIST(hdcaxyprm[i]), track.pt(), track.dcaXY()); + histos.fill(HIST(hdcazprm[i]), track.pt(), track.dcaZ()); + } + } + } + } else { + if (!mcParticle.isPhysicalPrimary()) { + if (mcParticle.getProcess() == kPDecay) { // particles from decay + histos.fill(HIST(hdcaxystr[i]), track.pt(), track.dcaXY()); + } else { + histos.fill(HIST(hdcaxymat[i]), track.pt(), track.dcaXY()); + } + } else { + histos.fill(HIST(hdcaxyprm[i]), track.pt(), track.dcaXY()); + } + } + + if (!passesDCAxyCut(track)) { // Skipping tracks that don't pass the standard cuts + return; + } + if (!passesDCAzCut(track)) { + return; + } + + // -- Number of MC recostructed tracks after DCAxy cut -- + if (!mcParticle.isPhysicalPrimary()) { // Is not physical primary + if (mcParticle.getProcess() == kPDecay) { // Is from decay + if (includeCentralityMC) { + histos.fill(HIST(hpt_num_str[i]), track.pt(), multiplicity, track.dcaXY()); + } else { + histos.fill(HIST(hpt_num_str[i]), track.pt(), track.dcaXY()); + } + if (track.hasTOF()) { + if (includeCentralityMC) { + histos.fill(HIST(hpt_numtof_str[i]), track.pt(), multiplicity, track.dcaXY()); + } else { + histos.fill(HIST(hpt_numtof_str[i]), track.pt(), track.dcaXY()); + } + } + } else { + if (includeCentralityMC) { + histos.fill(HIST(hpt_num_mat[i]), track.pt(), multiplicity, track.dcaXY()); + if (track.hasTOF()) { + histos.fill(HIST(hpt_numtof_mat[i]), track.pt(), multiplicity, track.dcaXY()); + } + } else { + histos.fill(HIST(hpt_num_mat[i]), track.pt(), track.dcaXY()); + if (track.hasTOF()) { + histos.fill(HIST(hpt_numtof_mat[i]), track.pt(), track.dcaXY()); + } + } + } + } else { // Is physical primary + if (includeCentralityMC) { + histos.fill(HIST(hpt_num_prm[i]), track.pt(), multiplicity, track.dcaXY()); + if (track.hasTOF()) { + histos.fill(HIST(hpt_numtof_prm[i]), track.pt(), multiplicity, track.dcaXY()); + } + } else { + histos.fill(HIST(hpt_num_prm[i]), track.pt(), track.dcaXY()); + if (track.hasTOF()) { + histos.fill(HIST(hpt_numtof_prm[i]), track.pt(), track.dcaXY()); + } + } + } + } + + template + void fillParticleHistogramsMC(const float multiplicity, ParticleType const& mcParticle) + { + if (!isParticleEnabled()) { // Check if the particle is enabled + return; + } + + if (mcParticle.pdgCode() != PDGs[i]) { + return; + } + + if (!mcParticle.isPhysicalPrimary()) { + if (mcParticle.getProcess() == kPDecay) { // Particles from decay + if (includeCentralityMC) { + histos.fill(HIST(hpt_den_str[i]), mcParticle.pt(), multiplicity); + } else { + histos.fill(HIST(hpt_den_str[i]), mcParticle.pt()); + } + } else { + if (includeCentralityMC) { + histos.fill(HIST(hpt_den_mat[i]), mcParticle.pt(), multiplicity); + } else { + histos.fill(HIST(hpt_den_mat[i]), mcParticle.pt()); + } + } + } else { + if (includeCentralityMC) { + histos.fill(HIST(hpt_den_prm[i]), mcParticle.pt(), multiplicity); + } else { + histos.fill(HIST(hpt_den_prm[i]), mcParticle.pt()); + } + } + } + + template + void fillParticleHistogramsMCRecoEvs(ParticleType const& mcParticle, RecoMCCollisions::iterator const& collision) + { + if (!isParticleEnabled()) { // Check if the particle is enabled + return; + } + + if (mcParticle.pdgCode() != PDGs[i]) { + return; + } + + const auto& mcCollision = collision.mcCollision_as(); + const float multiplicity = getMultiplicityMC(mcCollision); + + // Number of particles that pass the final selection + if (mcParticle.isPhysicalPrimary()) { + if (isEventSelected(collision)) { + if (includeCentralityMC) { + histos.fill(HIST(hpt_den_prm_goodev[i]), mcParticle.pt(), multiplicity, mcParticle.eta()); + } else { + histos.fill(HIST(hpt_den_prm_goodev[i]), mcParticle.pt(), mcParticle.eta()); + } + } + } + } + + template + void fillParticleHistogramsMCGenEvs(ParticleType const& mcParticle, GenMCCollisions::iterator const& mcCollision) + { + + if (!isParticleEnabled()) { // Check if the particle is enabled + return; + } + + if (mcParticle.pdgCode() != PDGs[i]) { + return; + } + + const float multiplicity = getMultiplicityMC(mcCollision); + if (mcParticle.isPhysicalPrimary()) { + // Number of primary particles generated from good events (z<10cm) + histos.fill(HIST(hpt_den_prm_mcgoodev[i]), mcParticle.pt(), multiplicity); + } + } + + //-- Process MC .--- + + Service pdgDB; + + Preslice perMCCol = aod::mcparticle::mcCollisionId; + SliceCache cache; + void processMC(soa::Join const& tracks, + aod::McParticles const& mcParticles, + GenMCCollisions const& mcCollisions, + RecoMCCollisions const& collisions) + { + // Fill number of generated and reconstructed collisions for normalization + histos.fill(HIST("MC/GenRecoCollisions"), 1.f, mcCollisions.size()); + histos.fill(HIST("MC/GenRecoCollisions"), 2.f, collisions.size()); + + for (const auto& track : tracks) { + if (!track.has_collision()) { + if (track.sign() > 0) { + histos.fill(HIST("MC/no_collision/pos"), track.pt()); + } else { + histos.fill(HIST("MC/no_collision/neg"), track.pt()); + } + continue; + } + if (!isEventSelected(track.collision_as())) { + continue; + } + if (!passesCutWoDCA(track)) { + continue; + } + if (!track.has_mcParticle()) { + if (track.sign() > 0) { + histos.fill(HIST("MC/fake/pos"), track.pt()); + } else { + histos.fill(HIST("MC/fake/neg"), track.pt()); + } + continue; + } + const auto& mcParticle = track.mcParticle(); + + static_for<0, 17>([&](auto i) { + fillTrackHistogramsMC(track, mcParticle, track.collision_as(), mcParticles); + }); + } + if (includeCentralityMC) { + for (const auto& collision : collisions) { + if (!collision.has_mcCollision()) { + continue; + } + const auto& mcCollision = collision.mcCollision_as(); + const auto& particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cache); + const float multiplicity = getMultiplicity(collision); + for (const auto& mcParticle : particlesInCollision) { + + if (std::abs(mcParticle.y()) > trkselOptions.cfgCutY) { + continue; + } + static_for<0, 17>([&](auto i) { + fillParticleHistogramsMC(multiplicity, mcParticle); + }); + } + } + } else { + for (const auto& mcParticle : mcParticles) { + if (std::abs(mcParticle.y()) > trkselOptions.cfgCutY) { + continue; + } + + const auto& mcCollision = mcParticle.mcCollision_as(); + const float multiplicity = getMultiplicityMC(mcCollision); + + static_for<0, 17>([&](auto i) { + fillParticleHistogramsMC(multiplicity, mcParticle); + }); + } + } + // Loop on reconstructed collisions + for (const auto& collision : collisions) { + if (!collision.has_mcCollision()) { + continue; + } + const auto& mcCollision = collision.mcCollision_as(); + const auto& particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cache); + if (evselOptions.cfgINELCut.value == 1) { + if (!o2::pwglf::isINELgt0mc(particlesInCollision, pdgDB)) { + continue; + } + } + if (evselOptions.cfgINELCut.value == 2) { + if (!o2::pwglf::isINELgt1mc(particlesInCollision, pdgDB)) { + continue; + } + } + + if (isEventSelected(collision)) { + histos.fill(HIST("MC/MultiplicityRecoEv"), getMultiplicityMC(mcCollision)); + } + for (const auto& mcParticle : particlesInCollision) { + if (std::abs(mcParticle.y()) > trkselOptions.cfgCutY) { + continue; + } + static_for<0, 17>([&](auto i) { + fillParticleHistogramsMCRecoEvs(mcParticle, collision); + }); + } + } + + // Loop on generated collisions + for (const auto& mcCollision : mcCollisions) { + if (std::abs(mcCollision.posZ()) > evselOptions.cfgCutVertex) { + continue; + } + histos.fill(HIST("MC/Multiplicity"), getMultiplicityMC(mcCollision)); + const auto& particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cache); + bool hasParticleInFT0C = false; + bool hasParticleInFT0A = false; + if (evselOptions.cfgINELCut.value == EvSelInelGt0Cut) { + if (!o2::pwglf::isINELgt0mc(particlesInCollision, pdgDB)) { + continue; + } + } + histos.fill(HIST("MC/MultiplicityMCINELgt0"), getMultiplicityMC(mcCollision)); + + for (const auto& mcParticle : particlesInCollision) { + if (std::abs(mcParticle.y()) > trkselOptions.cfgCutY) { + continue; + } + static_for<0, 17>([&](auto i) { + fillParticleHistogramsMCGenEvs(mcParticle, mcCollision); + }); + } /* + if (mcCollision.isInelGt0()) { + histos.fill(HIST("MC/GenRecoCollisions"), 3.f); + } + if (mcCollision.isInelGt1()) { + histos.fill(HIST("MC/GenRecoCollisions"), 4.f); + }*/ + if (hasParticleInFT0C && hasParticleInFT0A) { + histos.fill(HIST("MC/GenRecoCollisions"), 5.f); + } + } + } + PROCESS_SWITCH(tofSpectra, processMC, "Process MC", false); + +}; // end of spectra task + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc)}; }