diff --git a/PWGCF/Flow/Tasks/flowMc.cxx b/PWGCF/Flow/Tasks/flowMc.cxx index b310096ef40..682eeece738 100644 --- a/PWGCF/Flow/Tasks/flowMc.cxx +++ b/PWGCF/Flow/Tasks/flowMc.cxx @@ -54,13 +54,19 @@ struct FlowMc { Configurable minB{"minB", 0.0f, "min impact parameter"}; Configurable maxB{"maxB", 20.0f, "max impact parameter"}; + O2_DEFINE_CONFIGURABLE(cfgCutVertex, float, 10.0f, "Accepted z-vertex range") + O2_DEFINE_CONFIGURABLE(cfgCutPtMin, float, 0.2f, "Minimal pT for tracks") + O2_DEFINE_CONFIGURABLE(cfgCutPtMax, float, 1000.0f, "Maximal pT for tracks") + O2_DEFINE_CONFIGURABLE(cfgCutEta, float, 0.8f, "Eta range for tracks") O2_DEFINE_CONFIGURABLE(cfgOutputNUAWeights, bool, false, "Fill and output NUA weights") O2_DEFINE_CONFIGURABLE(cfgCutPtRefMin, float, 0.2f, "Minimal pT for ref tracks") O2_DEFINE_CONFIGURABLE(cfgCutPtRefMax, float, 3.0f, "Maximal pT for ref tracks") O2_DEFINE_CONFIGURABLE(cfgCutPtPOIMin, float, 0.2f, "Minimal pT for poi tracks") O2_DEFINE_CONFIGURABLE(cfgCutPtPOIMax, float, 10.0f, "Maximal pT for poi tracks") - O2_DEFINE_CONFIGURABLE(cfgCutTPCclu, float, 70.0f, "minimum TPC clusters") - O2_DEFINE_CONFIGURABLE(cfgCutITSclu, float, 6.0f, "minimum ITS clusters") + O2_DEFINE_CONFIGURABLE(cfgCutTPCclu, float, 50.0f, "minimum TPC clusters") + O2_DEFINE_CONFIGURABLE(cfgCutITSclu, float, 5.0f, "minimum ITS clusters") + O2_DEFINE_CONFIGURABLE(cfgCutChi2prTPCcls, float, 2.5f, "max chi2 per TPC clusters") + O2_DEFINE_CONFIGURABLE(cfgCutTPCcrossedrows, float, 70.0f, "minimum TPC crossed rows") O2_DEFINE_CONFIGURABLE(cfgFlowAcceptance, std::string, "", "CCDB path to acceptance object") O2_DEFINE_CONFIGURABLE(cfgFlowEfficiency, std::string, "", "CCDB path to efficiency object") O2_DEFINE_CONFIGURABLE(cfgCentVsIPTruth, std::string, "", "CCDB path to centrality vs IP truth") @@ -70,6 +76,10 @@ struct FlowMc { O2_DEFINE_CONFIGURABLE(cfgFlowCumulantNbootstrap, int, 30, "Number of subsamples") O2_DEFINE_CONFIGURABLE(cfgTrackDensityCorrUse, bool, false, "Use track density efficiency correction") O2_DEFINE_CONFIGURABLE(cfgTrackDensityCorrSlopeFactor, float, 1.0f, "A factor to scale the track density efficiency slope") + O2_DEFINE_CONFIGURABLE(cfgRecoEvRejectMC, bool, false, "reject both MC and Reco events when reco do not pass") + O2_DEFINE_CONFIGURABLE(cfgRecoEvSel8, bool, false, "require sel8 for reconstruction events") + O2_DEFINE_CONFIGURABLE(cfgRecoEvkIsGoodITSLayersAll, bool, false, "require kIsGoodITSLayersAll for reconstruction events") + O2_DEFINE_CONFIGURABLE(cfgRecoEvkNoSameBunchPileup, bool, false, "require kNoSameBunchPileup for reconstruction events") Configurable> cfgTrackDensityP0{"cfgTrackDensityP0", std::vector{0.6003720411, 0.6152630970, 0.6288860646, 0.6360694031, 0.6409494798, 0.6450540203, 0.6482117301, 0.6512592056, 0.6640008690, 0.6862631416, 0.7005738691, 0.7106567432, 0.7170728333}, "parameter 0 for track density efficiency correction"}; Configurable> cfgTrackDensityP1{"cfgTrackDensityP1", std::vector{-1.007592e-05, -8.932635e-06, -9.114538e-06, -1.054818e-05, -1.220212e-05, -1.312304e-05, -1.376433e-05, -1.412813e-05, -1.289562e-05, -1.050065e-05, -8.635725e-06, -7.380821e-06, -6.201250e-06}, "parameter 1 for track density efficiency correction"}; float maxEta = 0.8; @@ -81,6 +91,18 @@ struct FlowMc { ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f}, "pt axis"}; + // Filter for MCcollisions + Filter mccollisionFilter = nabs(aod::mccollision::posZ) < cfgCutVertex; + using FilteredMcCollisions = soa::Filtered; + // Filter for MCParticle + Filter particleFilter = (nabs(aod::mcparticle::eta) < cfgCutEta) && (aod::mcparticle::pt > cfgCutPtMin) && (aod::mcparticle::pt < cfgCutPtMax); + using FilteredMcParticles = soa::Filtered>; + // Filter for reco tracks + Filter trackFilter = (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPtMin) && (aod::track::pt < cfgCutPtMax) && (aod::track::tpcChi2NCl < cfgCutChi2prTPCcls); + using FilteredTracks = soa::Filtered>; + + // using FilteredTracks = soa::Join; + // Cent vs IP TH1D* mCentVsIPTruth = nullptr; bool centVsIPTruthLoaded = false; @@ -98,7 +120,6 @@ struct FlowMc { // Connect to ccdb Service ccdb; - Configurable ccdbNoLaterThan{"ccdbNoLaterThan", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object"}; Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; OutputObj fWeights{GFWWeights("weights")}; @@ -110,10 +131,22 @@ struct FlowMc { std::vector corrconfigsTruth; std::vector corrconfigsReco; TRandom3* fRndm = new TRandom3(0); + double epsilon = 1e-6; void init(InitContext&) { + ccdb->setURL(ccdbUrl.value); + ccdb->setCaching(true); + auto now = std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(); + ccdb->setCreatedNotAfter(now); + + const AxisSpec axisVertex{20, -10, 10, "Vtxz (cm)"}; + const AxisSpec axisEta{20, -1., 1., "#eta"}; + const AxisSpec axisCounter{1, 0, +1, ""}; // QA histograms + histos.add("mcEventCounter", "Monte Carlo Truth EventCounter", HistType::kTH1F, {axisCounter}); + histos.add("numberOfRecoCollisions", "numberOfRecoCollisions", HistType::kTH1F, {{10, -0.5f, 9.5f}}); + histos.add("RecoEventCounter", "Reconstruction EventCounter", HistType::kTH1F, {axisCounter}); histos.add("hnTPCClu", "Number of found TPC clusters", HistType::kTH1D, {{100, 40, 180}}); histos.add("hnITSClu", "Number of found ITS clusters", HistType::kTH1D, {{100, 0, 20}}); // pT histograms @@ -156,7 +189,9 @@ struct FlowMc { histos.add("hPtNchGeneratedLambda", "Reco production; pT (GeV/c); multiplicity", HistType::kTH2D, {axisPt, axisNch}); histos.add("hPtNchGlobalLambda", "Global production; pT (GeV/c); multiplicity", HistType::kTH2D, {axisPt, axisNch}); histos.add("hPtMCGen", "Monte Carlo Truth; pT (GeV/c);", {HistType::kTH1D, {axisPt}}); + histos.add("hEtaPtVtxzMCGen", "Monte Carlo Truth; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}}); histos.add("hPtMCGlobal", "Monte Carlo Global; pT (GeV/c);", {HistType::kTH1D, {axisPt}}); + histos.add("hEtaPtVtxzMCGlobal", "Monte Carlo Global; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}}); histos.add("hPhiWeightedTrDen", "corrected #phi distribution, considering track density", {HistType::kTH1D, {axisPhi}}); o2::framework::AxisSpec axis = axisPt; @@ -318,9 +353,35 @@ struct FlowMc { } } - using RecoTracks = soa::Join; + template + bool eventSelected(TCollision collision) + { + if (std::fabs(collision.posZ()) > cfgCutVertex) { + return 0; + } + if (cfgRecoEvSel8 && !collision.sel8()) { + return 0; + } + if (cfgRecoEvkNoSameBunchPileup && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { + // rejects collisions which are associated with the same "found-by-T0" bunch crossing + // https://indico.cern.ch/event/1396220/#1-event-selection-with-its-rof + return 0; + } + if (cfgRecoEvkIsGoodITSLayersAll && !collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) { + // from Jan 9 2025 AOT meeting + // cut time intervals with dead ITS staves + return 0; + } + return 1; + } + + template + bool trackSelected(TTrack track) + { + return ((track.tpcNClsFound() >= cfgCutTPCclu) && (track.tpcNClsCrossedRows() >= cfgCutTPCcrossedrows) && (track.itsNCls() >= cfgCutITSclu)); + } - void process(aod::McCollision const& mcCollision, aod::BCsWithTimestamps const&, soa::Join const& mcParticles, RecoTracks const&) + void process(FilteredMcCollisions::iterator const& mcCollision, aod::BCsWithTimestamps const&, soa::SmallGroups> const& collisions, FilteredMcParticles const& mcParticles, FilteredTracks const&) { float imp = mcCollision.impactParameter(); @@ -337,6 +398,21 @@ struct FlowMc { auto bc = mcCollision.bc_as(); loadCorrections(bc.timestamp()); + if (collisions.size() > -1) { + histos.fill(HIST("mcEventCounter"), 0.5); + histos.fill(HIST("numberOfRecoCollisions"), collisions.size()); // number of times coll was reco-ed + if (cfgRecoEvRejectMC) { + if (collisions.size() != 1) { // only pass those have one reconstruction event + return; + } + for (auto const& collision : collisions) { + if (!eventSelected(collision)) + return; + } + } + histos.fill(HIST("RecoEventCounter"), 0.5); + } + if (imp > minB && imp < maxB) { // event within range histos.fill(HIST("hImpactParameter"), imp); @@ -363,9 +439,9 @@ struct FlowMc { if (std::fabs(mcParticle.eta()) > maxEta) // main acceptance continue; if (mcParticle.has_tracks()) { - auto const& tracks = mcParticle.tracks_as(); + auto const& tracks = mcParticle.tracks_as(); for (auto const& track : tracks) { - if (!((track.tpcNClsFound() >= cfgCutTPCclu) && (track.itsNCls() >= cfgCutITSclu))) { + if (!trackSelected(track)) { continue; } if (cfgIsGlobalTrack && track.isGlobalTrack()) { @@ -418,6 +494,7 @@ struct FlowMc { histos.fill(HIST("hBVsPtVsPhiGenerated"), imp, deltaPhi, mcParticle.pt()); histos.fill(HIST("hPtNchGenerated"), mcParticle.pt(), nChGlobal); histos.fill(HIST("hPtMCGen"), mcParticle.pt()); + histos.fill(HIST("hEtaPtVtxzMCGen"), mcParticle.eta(), mcParticle.pt(), vtxz); if (pdgCode == PDG_t::kPiPlus) histos.fill(HIST("hPtNchGeneratedPion"), mcParticle.pt(), nChGlobal); if (pdgCode == PDG_t::kKPlus) @@ -437,9 +514,9 @@ struct FlowMc { bool validITSTrack = false; bool validITSABTrack = false; if (mcParticle.has_tracks()) { - auto const& tracks = mcParticle.tracks_as(); + auto const& tracks = mcParticle.tracks_as(); for (auto const& track : tracks) { - if (!((track.tpcNClsFound() >= cfgCutTPCclu) && (track.itsNCls() >= cfgCutITSclu))) { + if (!trackSelected(track)) { continue; } histos.fill(HIST("hnTPCClu"), track.tpcNClsFound()); @@ -456,10 +533,10 @@ struct FlowMc { if (track.hasTPC()) { validTPCTrack = true; } - if (track.hasITS() && track.itsChi2NCl() > -1e-6) { + if (track.hasITS() && track.itsChi2NCl() > -1. * epsilon) { validITSTrack = true; } - if (track.hasITS() && track.itsChi2NCl() < -1e-6) { + if (track.hasITS() && track.itsChi2NCl() < -1. * epsilon) { validITSABTrack = true; } } @@ -519,6 +596,7 @@ struct FlowMc { histos.fill(HIST("hBVsPtVsPhiGlobal"), imp, deltaPhi, mcParticle.pt(), wacc * weff); histos.fill(HIST("hPtNchGlobal"), mcParticle.pt(), nChGlobal); histos.fill(HIST("hPtMCGlobal"), mcParticle.pt()); + histos.fill(HIST("hEtaPtVtxzMCGlobal"), mcParticle.eta(), mcParticle.pt(), vtxz); if (pdgCode == PDG_t::kPiPlus) histos.fill(HIST("hPtNchGlobalPion"), mcParticle.pt(), nChGlobal); if (pdgCode == PDG_t::kKPlus) diff --git a/PWGCF/Flow/Tasks/flowTask.cxx b/PWGCF/Flow/Tasks/flowTask.cxx index b552add5223..158cb8820d9 100644 --- a/PWGCF/Flow/Tasks/flowTask.cxx +++ b/PWGCF/Flow/Tasks/flowTask.cxx @@ -110,8 +110,23 @@ struct FlowTask { Configurable> cfgUserDefineGFWName{"cfgUserDefineGFWName", std::vector{"Ch02Gap22", "Ch12Gap22"}, "User defined GFW Name"}; Configurable cfgUserPtVnCorrConfig{"cfgUserPtVnCorrConfig", {{"refP {2} refN {-2}", "refP {3} refN {-3}"}, {"ChGap22", "ChGap32"}, {0, 0}, {3, 3}}, "Configurations for vn-pt correlations"}; Configurable> cfgRunRemoveList{"cfgRunRemoveList", std::vector{-1}, "excluded run numbers"}; - Configurable> cfgTrackDensityP0{"cfgTrackDensityP0", std::vector{0.7217476707, 0.7384792571, 0.7542625668, 0.7640680200, 0.7701951667, 0.7755299053, 0.7805901710, 0.7849446786, 0.7957356586, 0.8113039262, 0.8211968966, 0.8280558878, 0.8329342135}, "parameter 0 for track density efficiency correction"}; - Configurable> cfgTrackDensityP1{"cfgTrackDensityP1", std::vector{-2.169488e-05, -2.191913e-05, -2.295484e-05, -2.556538e-05, -2.754463e-05, -2.816832e-05, -2.846502e-05, -2.843857e-05, -2.705974e-05, -2.477018e-05, -2.321730e-05, -2.203315e-05, -2.109474e-05}, "parameter 1 for track density efficiency correction"}; + struct : ConfigurableGroup { + Configurable> cfgTrackDensityP0{"cfgTrackDensityP0", std::vector{0.7217476707, 0.7384792571, 0.7542625668, 0.7640680200, 0.7701951667, 0.7755299053, 0.7805901710, 0.7849446786, 0.7957356586, 0.8113039262, 0.8211968966, 0.8280558878, 0.8329342135}, "parameter 0 for track density efficiency correction"}; + Configurable> cfgTrackDensityP1{"cfgTrackDensityP1", std::vector{-2.169488e-05, -2.191913e-05, -2.295484e-05, -2.556538e-05, -2.754463e-05, -2.816832e-05, -2.846502e-05, -2.843857e-05, -2.705974e-05, -2.477018e-05, -2.321730e-05, -2.203315e-05, -2.109474e-05}, "parameter 1 for track density efficiency correction"}; + O2_DEFINE_CONFIGURABLE(cfgMultCentHighCutFunction, std::string, "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x + 10.*([5] + [6]*x + [7]*x*x + [8]*x*x*x + [9]*x*x*x*x)", "Functional for multiplicity correlation cut"); + O2_DEFINE_CONFIGURABLE(cfgMultCentLowCutFunction, std::string, "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x - 3.*([5] + [6]*x + [7]*x*x + [8]*x*x*x + [9]*x*x*x*x)", "Functional for multiplicity correlation cut"); + O2_DEFINE_CONFIGURABLE(cfgMultT0CCutEnabled, bool, false, "Enable Global multiplicity vs T0C centrality cut") + Configurable> cfgMultT0CCutPars{"cfgMultT0CCutPars", std::vector{143.04, -4.58368, 0.0766055, -0.000727796, 2.86153e-06, 23.3108, -0.36304, 0.00437706, -4.717e-05, 1.98332e-07}, "Global multiplicity vs T0C centrality cut parameter values"}; + O2_DEFINE_CONFIGURABLE(cfgMultPVT0CCutEnabled, bool, false, "Enable PV multiplicity vs T0C centrality cut") + Configurable> cfgMultPVT0CCutPars{"cfgMultPVT0CCutPars", std::vector{195.357, -6.15194, 0.101313, -0.000955828, 3.74793e-06, 30.0326, -0.43322, 0.00476265, -5.11206e-05, 2.13613e-07}, "PV multiplicity vs T0C centrality cut parameter values"}; + O2_DEFINE_CONFIGURABLE(cfgMultMultHighCutFunction, std::string, "[0]+[1]*x + 5.*([2]+[3]*x)", "Functional for multiplicity correlation cut"); + O2_DEFINE_CONFIGURABLE(cfgMultMultLowCutFunction, std::string, "[0]+[1]*x - 5.*([2]+[3]*x)", "Functional for multiplicity correlation cut"); + O2_DEFINE_CONFIGURABLE(cfgMultGlobalPVCutEnabled, bool, false, "Enable global multiplicity vs PV multiplicity cut") + Configurable> cfgMultGlobalPVCutPars{"cfgMultGlobalPVCutPars", std::vector{-0.140809, 0.734344, 2.77495, 0.0165935}, "PV multiplicity vs T0C centrality cut parameter values"}; + std::vector multT0CCutPars; + std::vector multPVT0CCutPars; + std::vector multGlobalPVCutPars; + } cfgFuncParas; ConfigurableAxis axisPtHist{"axisPtHist", {100, 0., 10.}, "pt axis for histograms"}; ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 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, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.2, 2.4, 2.6, 2.8, 3, 3.5, 4, 5, 6, 8, 10}, "pt axis for histograms"}; @@ -184,10 +199,12 @@ struct FlowTask { TF1* funcV4; // Additional Event selection cuts - Copy from flowGenericFramework.cxx - TF1* fMultPVCutLow = nullptr; - TF1* fMultPVCutHigh = nullptr; - TF1* fMultCutLow = nullptr; - TF1* fMultCutHigh = nullptr; + TF1* fMultPVT0CCutLow = nullptr; + TF1* fMultPVT0CCutHigh = nullptr; + TF1* fMultT0CCutLow = nullptr; + TF1* fMultT0CCutHigh = nullptr; + TF1* fMultGlobalPVCutLow = nullptr; + TF1* fMultGlobalPVCutHigh = nullptr; TF1* fT0AV0AMean = nullptr; TF1* fT0AV0ASigma = nullptr; @@ -448,16 +465,24 @@ struct FlowTask { } fGFW->CreateRegions(); - if (cfgUseAdditionalEventCut) { - fMultPVCutLow = new TF1("fMultPVCutLow", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x - 3.5*([5]+[6]*x+[7]*x*x+[8]*x*x*x+[9]*x*x*x*x)", 0, 100); - fMultPVCutLow->SetParameters(3257.29, -121.848, 1.98492, -0.0172128, 6.47528e-05, 154.756, -1.86072, -0.0274713, 0.000633499, -3.37757e-06); - fMultPVCutHigh = new TF1("fMultPVCutHigh", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x + 3.5*([5]+[6]*x+[7]*x*x+[8]*x*x*x+[9]*x*x*x*x)", 0, 100); - fMultPVCutHigh->SetParameters(3257.29, -121.848, 1.98492, -0.0172128, 6.47528e-05, 154.756, -1.86072, -0.0274713, 0.000633499, -3.37757e-06); - - fMultCutLow = new TF1("fMultCutLow", "[0]+[1]*x+[2]*x*x+[3]*x*x*x - 2.*([4]+[5]*x+[6]*x*x+[7]*x*x*x+[8]*x*x*x*x)", 0, 100); - fMultCutLow->SetParameters(1654.46, -47.2379, 0.449833, -0.0014125, 150.773, -3.67334, 0.0530503, -0.000614061, 3.15956e-06); - fMultCutHigh = new TF1("fMultCutHigh", "[0]+[1]*x+[2]*x*x+[3]*x*x*x + 3.*([4]+[5]*x+[6]*x*x+[7]*x*x*x+[8]*x*x*x*x)", 0, 100); - fMultCutHigh->SetParameters(1654.46, -47.2379, 0.449833, -0.0014125, 150.773, -3.67334, 0.0530503, -0.000614061, 3.15956e-06); + if (cfgEvSelMultCorrelation) { + cfgFuncParas.multT0CCutPars = cfgFuncParas.cfgMultT0CCutPars; + cfgFuncParas.multPVT0CCutPars = cfgFuncParas.cfgMultPVT0CCutPars; + cfgFuncParas.multGlobalPVCutPars = cfgFuncParas.cfgMultGlobalPVCutPars; + fMultPVT0CCutLow = new TF1("fMultPVT0CCutLow", cfgFuncParas.cfgMultCentLowCutFunction->c_str(), 0, 100); + fMultPVT0CCutLow->SetParameters(&(cfgFuncParas.multPVT0CCutPars[0])); + fMultPVT0CCutHigh = new TF1("fMultPVT0CCutHigh", cfgFuncParas.cfgMultCentHighCutFunction->c_str(), 0, 100); + fMultPVT0CCutHigh->SetParameters(&(cfgFuncParas.multPVT0CCutPars[0])); + + fMultT0CCutLow = new TF1("fMultT0CCutLow", cfgFuncParas.cfgMultCentLowCutFunction->c_str(), 0, 100); + fMultT0CCutLow->SetParameters(&(cfgFuncParas.multT0CCutPars[0])); + fMultT0CCutHigh = new TF1("fMultT0CCutHigh", cfgFuncParas.cfgMultCentHighCutFunction->c_str(), 0, 100); + fMultT0CCutHigh->SetParameters(&(cfgFuncParas.multT0CCutPars[0])); + + fMultGlobalPVCutLow = new TF1("fMultGlobalPVCutLow", cfgFuncParas.cfgMultMultLowCutFunction->c_str(), 0, 4000); + fMultGlobalPVCutLow->SetParameters(&(cfgFuncParas.multGlobalPVCutPars[0])); + fMultGlobalPVCutHigh = new TF1("fMultGlobalPVCutHigh", cfgFuncParas.cfgMultMultHighCutFunction->c_str(), 0, 4000); + fMultGlobalPVCutHigh->SetParameters(&(cfgFuncParas.multGlobalPVCutPars[0])); fT0AV0AMean = new TF1("fT0AV0AMean", "[0]+[1]*x", 0, 200000); fT0AV0AMean->SetParameters(-1601.0581, 9.417652e-01); @@ -470,8 +495,8 @@ struct FlowTask { hFindPtBin = new TH1D("hFindPtBin", "hFindPtBin", pTEffBins.size() - 1, &pTEffBins[0]); funcEff.resize(pTEffBins.size() - 1); // LHC24g3 Eff - std::vector f1p0 = cfgTrackDensityP0; - std::vector f1p1 = cfgTrackDensityP1; + std::vector f1p0 = cfgFuncParas.cfgTrackDensityP0; + std::vector f1p1 = cfgFuncParas.cfgTrackDensityP1; for (uint ifunc = 0; ifunc < pTEffBins.size() - 1; ifunc++) { funcEff[ifunc] = new TF1(Form("funcEff%i", ifunc), "[0]+[1]*x", 0, 3000); funcEff[ifunc]->SetParameters(f1p0[ifunc], f1p1[ifunc]); @@ -656,14 +681,24 @@ struct FlowTask { registry.fill(HIST("hEventCountSpecific"), 9.5); if (cfgEvSelMultCorrelation) { - if (multNTracksPV < fMultPVCutLow->Eval(centrality)) - return 0; - if (multNTracksPV > fMultPVCutHigh->Eval(centrality)) - return 0; - if (multTrk < fMultCutLow->Eval(centrality)) - return 0; - if (multTrk > fMultCutHigh->Eval(centrality)) - return 0; + if (cfgFuncParas.cfgMultPVT0CCutEnabled) { + if (multNTracksPV < fMultPVT0CCutLow->Eval(centrality)) + return 0; + if (multNTracksPV > fMultPVT0CCutHigh->Eval(centrality)) + return 0; + } + if (cfgFuncParas.cfgMultT0CCutEnabled) { + if (multTrk < fMultT0CCutLow->Eval(centrality)) + return 0; + if (multTrk > fMultT0CCutHigh->Eval(centrality)) + return 0; + } + if (cfgFuncParas.cfgMultGlobalPVCutEnabled) { + if (multTrk < fMultGlobalPVCutLow->Eval(multNTracksPV)) + return 0; + if (multTrk > fMultGlobalPVCutHigh->Eval(multNTracksPV)) + return 0; + } } if (cfgEvSelMultCorrelation) registry.fill(HIST("hEventCountSpecific"), 10.5); diff --git a/PWGCF/TwoParticleCorrelations/Tasks/diHadronCor.cxx b/PWGCF/TwoParticleCorrelations/Tasks/diHadronCor.cxx index 6a1b61718f2..d47aff3c749 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/diHadronCor.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/diHadronCor.cxx @@ -28,7 +28,6 @@ #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/PIDResponse.h" -#include "Common/DataModel/PIDResponseITS.h" #include "Common/DataModel/TrackSelectionTables.h" #include "CommonConstants/MathConstants.h" @@ -41,8 +40,6 @@ #include "Framework/RunningWorkflowInfo.h" #include "Framework/StepTHn.h" #include "Framework/runDataProcessing.h" -#include "ReconstructionDataFormats/PID.h" -#include "ReconstructionDataFormats/Track.h" #include #include "TF1.h" @@ -98,14 +95,27 @@ struct DiHadronCor { O2_DEFINE_CONFIGURABLE(cfgCutOccupancyHigh, int, 2000, "High cut on TPC occupancy") O2_DEFINE_CONFIGURABLE(cfgCutOccupancyLow, int, 0, "Low cut on TPC occupancy") O2_DEFINE_CONFIGURABLE(cfgEfficiency, std::string, "", "CCDB path to efficiency object") + O2_DEFINE_CONFIGURABLE(cfgCentralityWeight, std::string, "", "CCDB path to centrality weight object") O2_DEFINE_CONFIGURABLE(cfgLocalEfficiency, bool, false, "Use local efficiency object") O2_DEFINE_CONFIGURABLE(cfgVerbosity, bool, false, "Verbose output") O2_DEFINE_CONFIGURABLE(cfgUseEventWeights, bool, false, "Use event weights for mixed event") O2_DEFINE_CONFIGURABLE(cfgUsePtOrder, bool, true, "enable trigger pT < associated pT cut") O2_DEFINE_CONFIGURABLE(cfgUsePtOrderInMixEvent, bool, true, "enable trigger pT < associated pT cut in mixed event") - O2_DEFINE_CONFIGURABLE(cfgPIDUseITSPID, bool, true, "Use ITS PID for particle identification") - O2_DEFINE_CONFIGURABLE(cfgPIDTofPtCut, float, 0.5f, "Minimum pt to use TOF N-sigma") - O2_DEFINE_CONFIGURABLE(cfgPIDParticle, int, 0, "1 = pion, 2 = kaon, 3 = proton, 0 for no PID") + struct : ConfigurableGroup { + O2_DEFINE_CONFIGURABLE(cfgMultCentHighCutFunction, std::string, "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x + 10.*([5] + [6]*x + [7]*x*x + [8]*x*x*x + [9]*x*x*x*x)", "Functional for multiplicity correlation cut"); + O2_DEFINE_CONFIGURABLE(cfgMultCentLowCutFunction, std::string, "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x - 3.*([5] + [6]*x + [7]*x*x + [8]*x*x*x + [9]*x*x*x*x)", "Functional for multiplicity correlation cut"); + O2_DEFINE_CONFIGURABLE(cfgMultT0CCutEnabled, bool, false, "Enable Global multiplicity vs T0C centrality cut") + Configurable> cfgMultT0CCutPars{"cfgMultT0CCutPars", std::vector{143.04, -4.58368, 0.0766055, -0.000727796, 2.86153e-06, 23.3108, -0.36304, 0.00437706, -4.717e-05, 1.98332e-07}, "Global multiplicity vs T0C centrality cut parameter values"}; + O2_DEFINE_CONFIGURABLE(cfgMultPVT0CCutEnabled, bool, false, "Enable PV multiplicity vs T0C centrality cut") + Configurable> cfgMultPVT0CCutPars{"cfgMultPVT0CCutPars", std::vector{195.357, -6.15194, 0.101313, -0.000955828, 3.74793e-06, 30.0326, -0.43322, 0.00476265, -5.11206e-05, 2.13613e-07}, "PV multiplicity vs T0C centrality cut parameter values"}; + O2_DEFINE_CONFIGURABLE(cfgMultMultHighCutFunction, std::string, "[0]+[1]*x + 5.*([2]+[3]*x)", "Functional for multiplicity correlation cut"); + O2_DEFINE_CONFIGURABLE(cfgMultMultLowCutFunction, std::string, "[0]+[1]*x - 5.*([2]+[3]*x)", "Functional for multiplicity correlation cut"); + O2_DEFINE_CONFIGURABLE(cfgMultGlobalPVCutEnabled, bool, false, "Enable global multiplicity vs PV multiplicity cut") + Configurable> cfgMultGlobalPVCutPars{"cfgMultGlobalPVCutPars", std::vector{-0.140809, 0.734344, 2.77495, 0.0165935}, "PV multiplicity vs T0C centrality cut parameter values"}; + std::vector multT0CCutPars; + std::vector multPVT0CCutPars; + std::vector multGlobalPVCutPars; + } cfgFuncParas; SliceCache cache; @@ -120,9 +130,6 @@ struct DiHadronCor { ConfigurableAxis axisVtxMix{"axisVtxMix", {VARIABLE_WIDTH, -10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, "vertex axis for mixed event histograms"}; ConfigurableAxis axisMultMix{"axisMultMix", {VARIABLE_WIDTH, 0, 10, 20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220, 240, 260}, "multiplicity / centrality axis for mixed event histograms"}; ConfigurableAxis axisSample{"axisSample", {cfgSampleSize, 0, cfgSampleSize}, "sample axis for histograms"}; - Configurable> pidTofNsigmaCut{"pidTofNsigmaCut", std::vector{1.5, 1.5, 1.5, -1.5, -1.5, -1.5}, "TOF n-sigma cut for pions_posNsigma, kaons_posNsigma, protons_posNsigma, pions_negNsigma, kaons_negNsigma, protons_negNsigma"}; - Configurable> pidItsNsigmaCut{"pidItsNsigmaCut", std::vector{3, 3, 3, -3, -3, -3}, "ITS n-sigma cut for pions_posNsigma, kaons_posNsigma, protons_posNsigma, pions_negNsigma, kaons_negNsigma, protons_negNsigma"}; - Configurable> pidTpcNsigmaCut{"pidTpcNsigmaCut", std::vector{10, 10, 10, -10, -10, -10}, "TOF n-sigma cut for pions_posNsigma, kaons_posNsigma, protons_posNsigma, pions_negNsigma, kaons_negNsigma, protons_negNsigma"}; ConfigurableAxis axisVertexEfficiency{"axisVertexEfficiency", {10, -10, 10}, "vertex axis for efficiency histograms"}; ConfigurableAxis axisEtaEfficiency{"axisEtaEfficiency", {20, -1.0, 1.0}, "eta axis for efficiency histograms"}; @@ -132,7 +139,7 @@ struct DiHadronCor { Filter collisionFilter = (nabs(aod::collision::posZ) < cfgCutVtxZ); Filter trackFilter = (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPtMin) && (aod::track::pt < cfgCutPtMax) && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true)) && (aod::track::tpcChi2NCl < cfgCutChi2prTPCcls) && (nabs(aod::track::dcaZ) < cfgCutDCAz); using FilteredCollisions = soa::Filtered>; - using FilteredTracks = soa::Filtered>; + using FilteredTracks = soa::Filtered>; using FilteredTracksWithMCLabels = soa::Filtered>; // Filter for MCParticle @@ -150,6 +157,7 @@ struct DiHadronCor { // Corrections TH3D* mEfficiency = nullptr; + TH1D* mCentralityWeight = nullptr; bool correctionsLoaded = false; // Define the outputs @@ -171,25 +179,17 @@ struct DiHadronCor { SameEvent = 1, MixedEvent = 3 }; - std::vector tofNsigmaCut; - std::vector itsNsigmaCut; - std::vector tpcNsigmaCut; - o2::aod::ITSResponse itsResponse; - enum Particles { - PIONS, - KAONS, - PROTONS - }; // persistent caches std::vector efficiencyAssociatedCache; // Additional Event selection cuts - Copy from flowGenericFramework.cxx - TF1* fMultPVCutLow = nullptr; - TF1* fMultPVCutHigh = nullptr; - TF1* fMultCutLow = nullptr; - TF1* fMultCutHigh = nullptr; - TF1* fMultMultPVCut = nullptr; + TF1* fMultPVT0CCutLow = nullptr; + TF1* fMultPVT0CCutHigh = nullptr; + TF1* fMultT0CCutLow = nullptr; + TF1* fMultT0CCutHigh = nullptr; + TF1* fMultGlobalPVCutLow = nullptr; + TF1* fMultGlobalPVCutHigh = nullptr; TF1* fT0AV0AMean = nullptr; TF1* fT0AV0ASigma = nullptr; @@ -225,16 +225,24 @@ struct DiHadronCor { registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(12, "cfgEvSelV0AT0ACut"); } - if (cfgUseAdditionalEventCut) { - fMultPVCutLow = new TF1("fMultPVCutLow", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x - 3.5*([5]+[6]*x+[7]*x*x+[8]*x*x*x+[9]*x*x*x*x)", 0, 100); - fMultPVCutLow->SetParameters(3257.29, -121.848, 1.98492, -0.0172128, 6.47528e-05, 154.756, -1.86072, -0.0274713, 0.000633499, -3.37757e-06); - fMultPVCutHigh = new TF1("fMultPVCutHigh", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x + 3.5*([5]+[6]*x+[7]*x*x+[8]*x*x*x+[9]*x*x*x*x)", 0, 100); - fMultPVCutHigh->SetParameters(3257.29, -121.848, 1.98492, -0.0172128, 6.47528e-05, 154.756, -1.86072, -0.0274713, 0.000633499, -3.37757e-06); - - fMultCutLow = new TF1("fMultCutLow", "[0]+[1]*x+[2]*x*x+[3]*x*x*x - 2.*([4]+[5]*x+[6]*x*x+[7]*x*x*x+[8]*x*x*x*x)", 0, 100); - fMultCutLow->SetParameters(1654.46, -47.2379, 0.449833, -0.0014125, 150.773, -3.67334, 0.0530503, -0.000614061, 3.15956e-06); - fMultCutHigh = new TF1("fMultCutHigh", "[0]+[1]*x+[2]*x*x+[3]*x*x*x + 3.*([4]+[5]*x+[6]*x*x+[7]*x*x*x+[8]*x*x*x*x)", 0, 100); - fMultCutHigh->SetParameters(1654.46, -47.2379, 0.449833, -0.0014125, 150.773, -3.67334, 0.0530503, -0.000614061, 3.15956e-06); + if (cfgEvSelMultCorrelation) { + cfgFuncParas.multT0CCutPars = cfgFuncParas.cfgMultT0CCutPars; + cfgFuncParas.multPVT0CCutPars = cfgFuncParas.cfgMultPVT0CCutPars; + cfgFuncParas.multGlobalPVCutPars = cfgFuncParas.cfgMultGlobalPVCutPars; + fMultPVT0CCutLow = new TF1("fMultPVT0CCutLow", cfgFuncParas.cfgMultCentLowCutFunction->c_str(), 0, 100); + fMultPVT0CCutLow->SetParameters(&(cfgFuncParas.multPVT0CCutPars[0])); + fMultPVT0CCutHigh = new TF1("fMultPVT0CCutHigh", cfgFuncParas.cfgMultCentHighCutFunction->c_str(), 0, 100); + fMultPVT0CCutHigh->SetParameters(&(cfgFuncParas.multPVT0CCutPars[0])); + + fMultT0CCutLow = new TF1("fMultT0CCutLow", cfgFuncParas.cfgMultCentLowCutFunction->c_str(), 0, 100); + fMultT0CCutLow->SetParameters(&(cfgFuncParas.multT0CCutPars[0])); + fMultT0CCutHigh = new TF1("fMultT0CCutHigh", cfgFuncParas.cfgMultCentHighCutFunction->c_str(), 0, 100); + fMultT0CCutHigh->SetParameters(&(cfgFuncParas.multT0CCutPars[0])); + + fMultGlobalPVCutLow = new TF1("fMultGlobalPVCutLow", cfgFuncParas.cfgMultMultLowCutFunction->c_str(), 0, 4000); + fMultGlobalPVCutLow->SetParameters(&(cfgFuncParas.multGlobalPVCutPars[0])); + fMultGlobalPVCutHigh = new TF1("fMultGlobalPVCutHigh", cfgFuncParas.cfgMultMultHighCutFunction->c_str(), 0, 4000); + fMultGlobalPVCutHigh->SetParameters(&(cfgFuncParas.multGlobalPVCutPars[0])); fT0AV0AMean = new TF1("fT0AV0AMean", "[0]+[1]*x", 0, 200000); fT0AV0AMean->SetParameters(-1601.0581, 9.417652e-01); @@ -310,10 +318,6 @@ struct DiHadronCor { same.setObject(new CorrelationContainer("sameEvent", "sameEvent", corrAxis, effAxis, userAxis)); mixed.setObject(new CorrelationContainer("mixedEvent", "mixedEvent", corrAxis, effAxis, userAxis)); - tofNsigmaCut = pidTofNsigmaCut; - itsNsigmaCut = pidItsNsigmaCut; - tpcNsigmaCut = pidTpcNsigmaCut; - LOGF(info, "End of init"); } @@ -358,9 +362,6 @@ struct DiHadronCor { template bool trackSelected(TTrack track) { - if (cfgPIDParticle && getNsigmaPID(track) != cfgPIDParticle) { - return false; - } return ((track.tpcNClsFound() >= cfgCutTPCclu) && (track.tpcNClsCrossedRows() >= cfgCutTPCCrossedRows) && (track.itsNCls() >= cfgCutITSclu)); } @@ -382,7 +383,7 @@ struct DiHadronCor { return true; } - void loadEfficiency(uint64_t timestamp) + void loadCorrection(uint64_t timestamp) { if (correctionsLoaded) { return; @@ -399,6 +400,13 @@ struct DiHadronCor { } LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgEfficiency.value.c_str(), (void*)mEfficiency); } + if (cfgCentralityWeight.value.empty() == false) { + mCentralityWeight = ccdb->getForTimeStamp(cfgCentralityWeight, timestamp); + if (mCentralityWeight == nullptr) { + LOGF(fatal, "Could not load efficiency histogram for trigger particles from %s", cfgCentralityWeight.value.c_str()); + } + LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgCentralityWeight.value.c_str(), (void*)mCentralityWeight); + } correctionsLoaded = true; } @@ -419,6 +427,19 @@ struct DiHadronCor { return true; } + bool getCentralityWeight(float& weightCent, const float centrality) + { + float weight = 1.; + if (mCentralityWeight) + weight = mCentralityWeight->GetBinContent(mCentralityWeight->FindBin(centrality)); + else + weight = 1.0; + if (weight == 0) + return false; + weightCent = weight; + return true; + } + // fill multiple histograms template void fillYield(TCollision collision, TTracks tracks) // function to fill the yield and etaphi histograms. @@ -659,14 +680,24 @@ struct DiHadronCor { auto multNTracksPV = collision.multNTracksPV(); if (cfgEvSelMultCorrelation) { - if (multNTracksPV < fMultPVCutLow->Eval(centrality)) - return 0; - if (multNTracksPV > fMultPVCutHigh->Eval(centrality)) - return 0; - if (multTrk < fMultCutLow->Eval(centrality)) - return 0; - if (multTrk > fMultCutHigh->Eval(centrality)) - return 0; + if (cfgFuncParas.cfgMultPVT0CCutEnabled) { + if (multNTracksPV < fMultPVT0CCutLow->Eval(centrality)) + return 0; + if (multNTracksPV > fMultPVT0CCutHigh->Eval(centrality)) + return 0; + } + if (cfgFuncParas.cfgMultT0CCutEnabled) { + if (multTrk < fMultT0CCutLow->Eval(centrality)) + return 0; + if (multTrk > fMultT0CCutHigh->Eval(centrality)) + return 0; + } + if (cfgFuncParas.cfgMultGlobalPVCutEnabled) { + if (multTrk < fMultGlobalPVCutLow->Eval(multNTracksPV)) + return 0; + if (multTrk > fMultGlobalPVCutHigh->Eval(multNTracksPV)) + return 0; + } } if (fillCounter && cfgEvSelMultCorrelation) registry.fill(HIST("hEventCountSpecific"), 10.5); @@ -704,12 +735,15 @@ struct DiHadronCor { return; } - loadEfficiency(bc.timestamp()); + loadCorrection(bc.timestamp()); registry.fill(HIST("eventcount"), SameEvent); // because its same event i put it in the 1 bin fillYield(collision, tracks); + float weightCent = 1.0f; + if (!cfgCentTableUnavailable) + getCentralityWeight(weightCent, cent); same->fillEvent(tracks.size(), CorrelationContainer::kCFStepReconstructed); - fillCorrelations(tracks, tracks, collision.posZ(), SameEvent, getMagneticField(bc.timestamp()), cent, 1.0f); + fillCorrelations(tracks, tracks, collision.posZ(), SameEvent, getMagneticField(bc.timestamp()), cent, weightCent); } PROCESS_SWITCH(DiHadronCor, processSame, "Process same event", true); @@ -759,13 +793,16 @@ struct DiHadronCor { registry.fill(HIST("eventcount"), MixedEvent); // fill the mixed event in the 3 bin auto bc = collision1.bc_as(); - loadEfficiency(bc.timestamp()); + loadCorrection(bc.timestamp()); float eventWeight = 1.0f; if (cfgUseEventWeights) { eventWeight = 1.0f / it.currentWindowNeighbours(); } + float weightCent = 1.0f; + if (!cfgCentTableUnavailable) + getCentralityWeight(weightCent, cent1); - fillCorrelations(tracks1, tracks2, collision1.posZ(), MixedEvent, getMagneticField(bc.timestamp()), cent1, eventWeight); + fillCorrelations(tracks1, tracks2, collision1.posZ(), MixedEvent, getMagneticField(bc.timestamp()), cent1, eventWeight * weightCent); } } @@ -989,56 +1026,6 @@ struct DiHadronCor { } } PROCESS_SWITCH(DiHadronCor, processOntheflyMixed, "Process on-the-fly mixed events", false); - - template - int getNsigmaPID(TTrack track) - { - // Computing Nsigma arrays for pion, kaon, and protons - std::array nSigmaTPC = {track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr()}; - std::array nSigmaTOF = {track.tofNSigmaPi(), track.tofNSigmaKa(), track.tofNSigmaPr()}; - std::array nSigmaITS = {itsResponse.nSigmaITS(track), itsResponse.nSigmaITS(track), itsResponse.nSigmaITS(track)}; - int pid = -1; - - std::array nSigmaToUse = cfgPIDUseITSPID ? nSigmaITS : nSigmaTPC; // Choose which nSigma to use: TPC or ITS - std::vector detectorNsigmaCut = cfgPIDUseITSPID ? itsNsigmaCut : tpcNsigmaCut; // Choose which nSigma to use: TPC or ITS - - bool isPion, isKaon, isProton; - bool isDetectedPion = nSigmaToUse[0] < detectorNsigmaCut[0] && nSigmaToUse[0] > detectorNsigmaCut[0 + 3]; - bool isDetectedKaon = nSigmaToUse[1] < detectorNsigmaCut[1] && nSigmaToUse[1] > detectorNsigmaCut[1 + 3]; - bool isDetectedProton = nSigmaToUse[2] < detectorNsigmaCut[2] && nSigmaToUse[2] > detectorNsigmaCut[2 + 3]; - - bool isTofPion = nSigmaTOF[0] < tofNsigmaCut[0] && nSigmaTOF[0] > tofNsigmaCut[0 + 3]; - bool isTofKaon = nSigmaTOF[1] < tofNsigmaCut[1] && nSigmaTOF[1] > tofNsigmaCut[1 + 3]; - bool isTofProton = nSigmaTOF[2] < tofNsigmaCut[2] && nSigmaTOF[2] > tofNsigmaCut[2 + 3]; - - if (track.pt() > cfgPIDTofPtCut && !track.hasTOF()) { - return 0; - } else if (track.pt() > cfgPIDTofPtCut && track.hasTOF()) { - isPion = isTofPion && isDetectedPion; - isKaon = isTofKaon && isDetectedKaon; - isProton = isTofProton && isDetectedProton; - } else { - isPion = isDetectedPion; - isKaon = isDetectedKaon; - isProton = isDetectedProton; - } - - if ((isPion && isKaon) || (isPion && isProton) || (isKaon && isProton)) { - return 0; // more than one particle satisfy the criteria - } - - if (isPion) { - pid = PIONS; - } else if (isKaon) { - pid = KAONS; - } else if (isProton) { - pid = PROTONS; - } else { - return 0; // no particle satisfies the criteria - } - - return pid + 1; // shift the pid by 1, 1 = pion, 2 = kaon, 3 = proton - } }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)