diff --git a/PWGCF/Flow/Tasks/flowPbpbPikp.cxx b/PWGCF/Flow/Tasks/flowPbpbPikp.cxx index 311cccc5b84..0f3848d0098 100644 --- a/PWGCF/Flow/Tasks/flowPbpbPikp.cxx +++ b/PWGCF/Flow/Tasks/flowPbpbPikp.cxx @@ -52,6 +52,7 @@ #include #include +#include using namespace o2; using namespace o2::framework; @@ -88,6 +89,10 @@ struct FlowPbpbPikp { O2_DEFINE_CONFIGURABLE(cfgCutOccupancy, int, 3000, "Occupancy cut") O2_DEFINE_CONFIGURABLE(cfgUseGlobalTrack, bool, true, "use Global track") O2_DEFINE_CONFIGURABLE(cfgITScluster, int, 0, "Number of ITS cluster") + O2_DEFINE_CONFIGURABLE(cfgTrackDensityCorrUse, bool, true, "Use track density efficiency correction") + + 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"}; ConfigurableAxis axisVertex{"axisVertex", {20, -10, 10}, "vertex axis for histograms"}; ConfigurableAxis axisPhi{"axisPhi", {60, 0.0, constants::math::TwoPI}, "phi axis for histograms"}; @@ -127,6 +132,13 @@ struct FlowPbpbPikp { std::vector mAcceptance; bool correctionsLoaded = false; + // local track density correction + std::vector funcEff; + TH1D* hFindPtBin; + TF1* funcV2; + TF1* funcV3; + TF1* funcV4; + void init(InitContext const&) { ccdb->setURL(ccdbUrl.value); @@ -199,6 +211,10 @@ struct FlowPbpbPikp { for (int i = 0; i < fPtAxis->GetNbins(); i++) oba->Add(new TNamed(Form("Pr08Gap24_pt_%i", i + 1), "Pr08Gap24_pTDiff")); + oba->Add(new TNamed("Pi08BGap22", "Pi08BGap22")); + for (int i = 0; i < fPtAxis->GetNbins(); i++) + oba->Add(new TNamed(Form("Pi08BGap22_pt_%i", i + 1), "Pi08BGap22_pTDiff")); + fFC->SetName("FlowContainer"); fFC->SetXAxis(fPtAxis); fFC->Initialize(oba, axisMultiplicity, cfgNbootstrap); @@ -219,6 +235,9 @@ struct FlowPbpbPikp { fGFW->AddRegion("poiNpi", -0.8, -0.4, 1 + fPtAxis->GetNbins(), 2); fGFW->AddRegion("olNpi", -0.8, -0.4, 1 + fPtAxis->GetNbins(), 16); + fGFW->AddRegion("poiPpi", 0.4, 0.8, 1 + fPtAxis->GetNbins(), 2); + fGFW->AddRegion("olPpi", 0.4, 0.8, 1 + fPtAxis->GetNbins(), 16); + // kaon fGFW->AddRegion("poiNk", -0.8, -0.4, 1 + fPtAxis->GetNbins(), 4); fGFW->AddRegion("olNk", -0.8, -0.4, 1 + fPtAxis->GetNbins(), 32); @@ -251,7 +270,30 @@ struct FlowPbpbPikp { corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiNk refN08 | olNk {2 2} refP08 {-2 -2}", "Ka08Gap24", kTRUE)); corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiNpr refN08 | olNpr {2 2} refP08 {-2 -2}", "Pr08Gap24", kTRUE)); + // Backward correlations + corrconfigs.push_back(fGFW->GetCorrelatorConfig("refN08 {2} refP08 {-2}", "Pi08BGap22", kFALSE)); + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPpi refP08 | olPpi {2} refN08 {-2}", "Pi08BGap22", kTRUE)); + fGFW->CreateRegions(); + + if (cfgTrackDensityCorrUse) { + std::vector pTEffBins = {0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.4, 1.8, 2.2, 2.6, 3.0}; + 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; + 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]); + } + funcV2 = new TF1("funcV2", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x", 0, 100); + funcV2->SetParameters(0.0186111, 0.00351907, -4.38264e-05, 1.35383e-07, -3.96266e-10); + funcV3 = new TF1("funcV3", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x", 0, 100); + funcV3->SetParameters(0.0174056, 0.000703329, -1.45044e-05, 1.91991e-07, -1.62137e-09); + funcV4 = new TF1("funcV4", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x", 0, 100); + funcV4->SetParameters(0.008845, 0.000259668, -3.24435e-06, 4.54837e-08, -6.01825e-10); + } } enum Particles { @@ -523,6 +565,37 @@ struct FlowPbpbPikp { int pidIndex; loadCorrections(bc); // load corrections for the each event + // Track loop for calculating the Qn angles + double psi2Est = 0, psi3Est = 0, psi4Est = 0; + float wEPeff = 1; + double v2 = 0, v3 = 0, v4 = 0; + // be cautious, this only works for Pb-Pb + // esimate the Qn angles and vn for this event + if (cfgTrackDensityCorrUse) { + double q2x = 0, q2y = 0; + double q3x = 0, q3y = 0; + double q4x = 0, q4y = 0; + for (const auto& track : tracks) { + bool withinPtRef = (cfgCutPtMin < track.pt()) && (track.pt() < cfgCutPtMax); // within RF pT range + if (withinPtRef) { + q2x += std::cos(2 * track.phi()); + q2y += std::sin(2 * track.phi()); + q3x += std::cos(3 * track.phi()); + q3y += std::sin(3 * track.phi()); + q4x += std::cos(4 * track.phi()); + q4y += std::sin(4 * track.phi()); + } + } + psi2Est = std::atan2(q2y, q2x) / 2.; + psi3Est = std::atan2(q3y, q3x) / 3.; + psi4Est = std::atan2(q4y, q4x) / 4.; + + v2 = funcV2->Eval(cent); + v3 = funcV3->Eval(cent); + v4 = funcV4->Eval(cent); + } + + // Actual track loop for (auto const& track : tracks) { if (!selectionTrack(track)) continue; @@ -540,6 +613,9 @@ struct FlowPbpbPikp { // pidIndex = getBayesPIDIndex(track); pidIndex = getNsigmaPID(track); + + weff = 1; // Initializing weff for each track + // NUA weights if (cfgOutputNUAWeights) fillWeights(track, vtxz, pidIndex, runNumber); @@ -550,10 +626,24 @@ struct FlowPbpbPikp { if (withinPtRef && withinPtPOI && pidIndex) waccRef = waccPOI; // if particle is both (then it's overlap), override ref with POI + // Track density correction + if (cfgTrackDensityCorrUse && withinPtRef) { + double fphi = v2 * std::cos(2 * (track.phi() - psi2Est)) + v3 * std::cos(3 * (track.phi() - psi3Est)) + v4 * std::cos(4 * (track.phi() - psi4Est)); + fphi = (1 + 2 * fphi); + int pTBinForEff = hFindPtBin->FindBin(track.pt()); + if (pTBinForEff >= 1 && pTBinForEff <= hFindPtBin->GetNbinsX()) { + wEPeff = funcEff[pTBinForEff - 1]->Eval(fphi * tracks.size()); + if (wEPeff > 0.) { + wEPeff = 1. / wEPeff; + weff *= wEPeff; + } + } + } // end of track density correction loop + if (withinPtRef) { histos.fill(HIST("hPhiWeighted"), track.phi(), waccRef); fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), waccRef * weff, 1); - fGFW->Fill(track.eta(), 1, track.phi(), waccRef * weff, 512); + fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), waccRef * weff, 512); } if (withinPtPOI) { fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), waccPOI * weff, 128); diff --git a/PWGCF/Flow/Tasks/resonancesGfwFlow.cxx b/PWGCF/Flow/Tasks/resonancesGfwFlow.cxx index 6bde6c05c99..145c2f82e3a 100644 --- a/PWGCF/Flow/Tasks/resonancesGfwFlow.cxx +++ b/PWGCF/Flow/Tasks/resonancesGfwFlow.cxx @@ -38,6 +38,7 @@ #include "Common/DataModel/Multiplicity.h" #include "Common/Core/TrackSelection.h" #include "Common/Core/trackUtilities.h" +#include "Common/Core/RecoDecay.h" #include "CommonConstants/PhysicsConstants.h" #include "PWGLF/DataModel/EPCalibrationTables.h" @@ -135,10 +136,19 @@ struct ResonancesGfwFlow { O2_DEFINE_CONFIGURABLE(cfgUseWeightPhiEtaPt, bool, false, "Use Phi, Eta, Pt dependent NUA weights") O2_DEFINE_CONFIGURABLE(cfgNbootstrap, int, 10, "Number of subsamples") O2_DEFINE_CONFIGURABLE(cfgUseBootStrap, bool, true, "Use bootstrap for error estimation") - O2_DEFINE_CONFIGURABLE(cfgTrackDensityCorrUse, bool, false, "Use track density efficiency correction") + O2_DEFINE_CONFIGURABLE(cfgTrackDensityCorrUse, bool, true, "Use track density efficiency correction") + O2_DEFINE_CONFIGURABLE(cfgUsePhi, bool, true, "Analyze Phi") O2_DEFINE_CONFIGURABLE(cfgUseK0, bool, true, "Analyze K0") O2_DEFINE_CONFIGURABLE(cfgUseLambda, bool, true, "Analyze Lambda") - O2_DEFINE_CONFIGURABLE(cfgUsePhi, bool, true, "Analyze Phi") + + enum OutputSpecies { + Ref = 0, + K0 = 1, + Lambda = 2, + AnLambda = 3, + Phi = 4, + kCount_OutputSpecies + }; struct : ConfigurableGroup { Configurable> cfgCosPAs{"cfgCosPAs", std::vector{0.97f, 0.995f, 0.04f}, "Minimum Pointing angle for resonances [K0, Lambda, Phi]"}; @@ -146,7 +156,7 @@ struct ResonancesGfwFlow { Configurable> cfgMassMin{"cfgMassMin", std::vector{0.44f, 1.1f, 0.99f}, "Minimum mass for resonances [K0, Lambda, Phi]"}; Configurable> cfgMassMax{"cfgMassMax", std::vector{0.56f, 1.16f, 1.06f}, "Maximum mass for resonances [K0, Lambda, Phi]"}; Configurable> cfgNMassBins{"cfgNMassBins", std::vector{70, 70, 70}, "Invariant mass bins for resonances [K0, Lambda, Phi]"}; - Configurable> cfgMccCut{"cfgMccCut", std::vector{0.005f, 0.01f}, "MCC cut for resonances [K0, Lambda]"}; + Configurable> cfgMccCut{"cfgMccCut", std::vector{0.005f, 0.01f, 0.0f}, "MCC cut for resonances [K0, Lambda, Phi]"}; Configurable> cfgPosTrackPt{"cfgPosTrackPt", std::vector{0.15f, 0.15f, 0.15f}, "Pt cut for positive track of resonances [K0, Lambda, Phi]"}; Configurable> cfgNegTrackPt{"cfgNegTrackPt", std::vector{0.15f, 0.15f, 0.15f}, "Pt cut for negative track of resonances [K0, Lambda, Phi]"}; } resoCuts; @@ -200,15 +210,6 @@ struct ResonancesGfwFlow { TAxis* fLambdaMassAxis; TRandom3* fRndm = new TRandom3(0); - enum OutputSpecies { - Ref = 0, - K0 = 1, - Lambda = 2, - AnLambda = 3, - Phi = 4, - kCount_OutputSpecies - }; - std::vector mAcceptance; bool correctionsLoaded = false; @@ -243,13 +244,41 @@ struct ResonancesGfwFlow { histos.add("hPhiPhi", "", {HistType::kTH1D, {axisPhi}}); histos.add("hPhiEta", "", {HistType::kTH1D, {axisEta}}); histos.add("hPhiMass_sparse", "", {HistType::kTHnSparseD, {{axisPhiMass, axisPt, axisMultiplicity}}}); + histos.add("hPhimassSparse_RD", "", {HistType::kTHnSparseD, {{axisPhiMass, axisPt, axisMultiplicity}}}); histos.add("Phid22Fpt", "", {HistType::kTProfile3D, {axisPt, axisPhiMass, axisMultiplicity}}); histos.add("Phid24Fpt", "", {HistType::kTProfile3D, {axisPt, axisPhiMass, axisMultiplicity}}); histos.add("Phid22Bpt", "", {HistType::kTProfile3D, {axisPt, axisPhiMass, axisMultiplicity}}); histos.add("Phid24Bpt", "", {HistType::kTProfile3D, {axisPt, axisPhiMass, axisMultiplicity}}); } + if (cfgUseK0) { + histos.add("PlusTPC_K0", "", {HistType::kTH2D, {{axisPt, axisTPCsignal}}}); + histos.add("MinusTPC_K0", "", {HistType::kTH2D, {{axisPt, axisTPCsignal}}}); + histos.add("PlusTOF_K0", "", {HistType::kTH2D, {{axisPt, axisTOFsignal}}}); + histos.add("MinusTOF_K0", "", {HistType::kTH2D, {{axisPt, axisTOFsignal}}}); + histos.add("hK0Phi", "", {HistType::kTH1D, {axisPhi}}); + histos.add("hK0Eta", "", {HistType::kTH1D, {axisEta}}); + histos.add("hK0Mass_sparse", "", {HistType::kTHnSparseF, {{axisK0Mass, axisPt, axisMultiplicity}}}); + histos.add("hK0s", "", {HistType::kTH1D, {singleCount}}); + histos.add("K0d22Fpt", "", {HistType::kTProfile3D, {axisPt, axisK0Mass, axisMultiplicity}}); + histos.add("K0d24Fpt", "", {HistType::kTProfile3D, {axisPt, axisK0Mass, axisMultiplicity}}); + histos.add("K0d22Bpt", "", {HistType::kTProfile3D, {axisPt, axisK0Mass, axisMultiplicity}}); + histos.add("K0d24Bpt", "", {HistType::kTProfile3D, {axisPt, axisK0Mass, axisMultiplicity}}); + + histos.add("hK0Count", "Number of K0;; Count", {HistType::kTH1D, {{11, 0, 11}}}); + histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(1, "K0 candidates"); + histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(2, "Daughter pt"); + histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(3, "Mass cut"); + histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(4, "Rapidity cut"); + histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(5, "DCA to PV"); + histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(6, "DCA between daughters"); + histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(7, "V0radius"); + histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(8, "CosPA"); + histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(9, "Proper lifetime"); + histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(10, "Daughter track selection"); + histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(11, "Mass cross check"); + } if (cfgUseLambda) { histos.add("PlusTPC_L", "", {HistType::kTH2D, {{axisPt, axisTPCsignal}}}); histos.add("MinusTPC_L", "", {HistType::kTH2D, {{axisPt, axisTPCsignal}}}); @@ -291,35 +320,6 @@ struct ResonancesGfwFlow { histos.get(HIST("hLambdaCount"))->GetXaxis()->SetBinLabel(11, "Mass cross check"); } - if (cfgUseK0) { - histos.add("PlusTPC_K0", "", {HistType::kTH2D, {{axisPt, axisTPCsignal}}}); - histos.add("MinusTPC_K0", "", {HistType::kTH2D, {{axisPt, axisTPCsignal}}}); - histos.add("PlusTOF_K0", "", {HistType::kTH2D, {{axisPt, axisTOFsignal}}}); - histos.add("MinusTOF_K0", "", {HistType::kTH2D, {{axisPt, axisTOFsignal}}}); - histos.add("hK0Phi", "", {HistType::kTH1D, {axisPhi}}); - histos.add("hK0Eta", "", {HistType::kTH1D, {axisEta}}); - histos.add("hK0Mass_sparse", "", {HistType::kTHnSparseF, {{axisK0Mass, axisPt, axisMultiplicity}}}); - histos.add("hK0s", "", {HistType::kTH1D, {singleCount}}); - - histos.add("K0d22Fpt", "", {HistType::kTProfile3D, {axisPt, axisK0Mass, axisMultiplicity}}); - histos.add("K0d24Fpt", "", {HistType::kTProfile3D, {axisPt, axisK0Mass, axisMultiplicity}}); - histos.add("K0d22Bpt", "", {HistType::kTProfile3D, {axisPt, axisK0Mass, axisMultiplicity}}); - histos.add("K0d24Bpt", "", {HistType::kTProfile3D, {axisPt, axisK0Mass, axisMultiplicity}}); - - histos.add("hK0Count", "Number of K0;; Count", {HistType::kTH1D, {{11, 0, 11}}}); - histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(1, "K0 candidates"); - histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(2, "Daughter pt"); - histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(3, "Mass cut"); - histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(4, "Rapidity cut"); - histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(5, "DCA to PV"); - histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(6, "DCA between daughters"); - histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(7, "V0radius"); - histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(8, "CosPA"); - histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(9, "Proper lifetime"); - histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(10, "Daughter track selection"); - histos.get(HIST("hK0Count"))->GetXaxis()->SetBinLabel(11, "Mass cross check"); - } - histos.add("hEventCount", "Number of Event;; Count", {HistType::kTH1D, {{8, 0, 8}}}); histos.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(1, "Filtered event"); histos.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(2, "After sel8"); @@ -355,10 +355,10 @@ struct ResonancesGfwFlow { refC22Boot[i] = histos.add(Form("BootStrap/Refc22_bootstrap_%d", i), "", {HistType::kTProfile, {axisMultiplicity}}); refC24Boot[i] = histos.add(Form("BootStrap/Refc24_bootstrap_%d", i), "", {HistType::kTProfile, {axisMultiplicity}}); if (cfgUsePhi) { - phiD22FPtBoot[i] = histos.add(Form("BootStrap/Phid22Fpt_bootstrap_%d", i), "", {HistType::kTProfile3D, {axisPt, axisPhiMass, axisMultiplicity}}); - phiD24FPtBoot[i] = histos.add(Form("BootStrap/Phid24Fpt_bootstrap_%d", i), "", {HistType::kTProfile3D, {axisPt, axisPhiMass, axisMultiplicity}}); - phiD22BPtBoot[i] = histos.add(Form("BootStrap/Phid22Bpt_bootstrap_%d", i), "", {HistType::kTProfile3D, {axisPt, axisPhiMass, axisMultiplicity}}); - phiD24BPtBoot[i] = histos.add(Form("BootStrap/Phid24Bpt_bootstrap_%d", i), "", {HistType::kTProfile3D, {axisPt, axisPhiMass, axisMultiplicity}}); + phiD22FPtBoot[i] = histos.add(Form("BootStrap/phid22Fpt_bootstrap_%d", i), "", {HistType::kTProfile3D, {axisPt, axisPhiMass, axisMultiplicity}}); + phiD24FPtBoot[i] = histos.add(Form("BootStrap/phid24Fpt_bootstrap_%d", i), "", {HistType::kTProfile3D, {axisPt, axisPhiMass, axisMultiplicity}}); + phiD22BPtBoot[i] = histos.add(Form("BootStrap/phid22Bpt_bootstrap_%d", i), "", {HistType::kTProfile3D, {axisPt, axisPhiMass, axisMultiplicity}}); + phiD24BPtBoot[i] = histos.add(Form("BootStrap/phid24Bpt_bootstrap_%d", i), "", {HistType::kTProfile3D, {axisPt, axisPhiMass, axisMultiplicity}}); } if (cfgUseK0) { k0D22FPtBoot[i] = histos.add(Form("BootStrap/k0d22Fpt_bootstrap_%d", i), "", {HistType::kTProfile3D, {axisPt, axisK0Mass, axisMultiplicity}}); @@ -385,7 +385,7 @@ struct ResonancesGfwFlow { double* ptBins = &(axis.binEdges)[0]; fPtAxis = new TAxis(nPtBins, ptBins); - fPhiMassAxis = new TAxis(vMassBins[Phi - 2], 0.99, 1.06); + fPhiMassAxis = new TAxis(vMassBins[Phi - 2], vMassMin[Phi - 2], vMassMax[Phi - 2]); fK0MassAxis = new TAxis(vMassBins[K0 - 1], vMassMin[K0 - 1], vMassMax[K0 - 1]); fLambdaMassAxis = new TAxis(vMassBins[Lambda - 1], vMassMin[Lambda - 1], vMassMax[Lambda - 1]); @@ -433,29 +433,29 @@ struct ResonancesGfwFlow { corrconfigs.push_back(fGFW->GetCorrelatorConfig("refN08 {2 2} refP08 {-2 -2}", "Ref08Gap24", kFALSE)); //--------- pt differential pois - if (cfgUsePhi) { - corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiNphi refN08 | olNphi {2} refP08 {-2}", "PhiF08Gap22", kTRUE)); - corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiNphi refN08 | olNphi {2 2} refP08 {-2 -2}", "PhiF08Gap24", kTRUE)); - corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPphi refP08 | olPphi {2} refN08 {-2}", "PhiB08Gap22", kTRUE)); - corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPphi refP08 | olPphi {2 2} refN08 {-2 -2}", "PhiB08Gap24", kTRUE)); - } - if (cfgUseK0) { - corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiNk0 refN08 | olNk0 {2} refP08 {-2}", "KsF08Gap22", kTRUE)); - corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiNk0 refN08 | olNk0 {2 2} refP08 {-2 -2}", "KsF08Gap24", kTRUE)); - corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPk0 refP08 | olPk0 {2} refN08 {-2}", "KsB08Gap22", kTRUE)); - corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPk0 refP08 | olPk0 {2 2} refN08 {-2 -2}", "KsB08Gap24", kTRUE)); - } - if (cfgUseLambda) { - corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiNlam refN08 | olNlam {2} refP08 {-2}", "LamF08Gap22", kTRUE)); - corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiNlam refN08 | olNlam {2 2} refP08 {-2 -2}", "LamF08Gap24", kTRUE)); - corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPlam refP08 | olPlam {2} refN08 {-2}", "LamB08Gap22", kTRUE)); - corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPlam refP08 | olPlam {2 2} refN08 {-2 -2}", "LamB08Gap24", kTRUE)); - - corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiNantilam refN08 | olNantilam {2} refP08 {-2}", "AnLamF08Gap22", kTRUE)); - corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiNantilam refN08 | olNantilam {2 2} refP08 {-2 -2}", "AnLamF08Gap24", kTRUE)); - corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPantilam refP08 | olPantilam {2} refN08 {-2}", "AnLamB08Gap22", kTRUE)); - corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPantilam refP08 | olPantilam {2 2} refN08 {-2 -2}", "AnLamB08Gap24", kTRUE)); - } + // Phi + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiNphi refN08 | olNphi {2} refP08 {-2}", "PhiF08Gap22", kTRUE)); + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiNphi refN08 | olNphi {2 2} refP08 {-2 -2}", "PhiF08Gap24", kTRUE)); + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPphi refP08 | olPphi {2} refN08 {-2}", "PhiB08Gap22", kTRUE)); + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPphi refP08 | olPphi {2 2} refN08 {-2 -2}", "PhiB08Gap24", kTRUE)); + + // K0 + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiNk0 refN08 | olNk0 {2} refP08 {-2}", "KsF08Gap22", kTRUE)); + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiNk0 refN08 | olNk0 {2 2} refP08 {-2 -2}", "KsF08Gap24", kTRUE)); + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPk0 refP08 | olPk0 {2} refN08 {-2}", "KsB08Gap22", kTRUE)); + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPk0 refP08 | olPk0 {2 2} refN08 {-2 -2}", "KsB08Gap24", kTRUE)); + + // Lambda + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiNlam refN08 | olNlam {2} refP08 {-2}", "LamF08Gap22", kTRUE)); + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiNlam refN08 | olNlam {2 2} refP08 {-2 -2}", "LamF08Gap24", kTRUE)); + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPlam refP08 | olPlam {2} refN08 {-2}", "LamB08Gap22", kTRUE)); + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPlam refP08 | olPlam {2 2} refN08 {-2 -2}", "LamB08Gap24", kTRUE)); + + // Antilambda + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiNantilam refN08 | olNantilam {2} refP08 {-2}", "AnLamF08Gap22", kTRUE)); + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiNantilam refN08 | olNantilam {2 2} refP08 {-2 -2}", "AnLamF08Gap24", kTRUE)); + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPantilam refP08 | olPantilam {2} refN08 {-2}", "AnLamB08Gap22", kTRUE)); + corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiPantilam refP08 | olPantilam {2 2} refN08 {-2 -2}", "AnLamB08Gap24", kTRUE)); fGFW->CreateRegions(); @@ -480,7 +480,7 @@ struct ResonancesGfwFlow { } template - void fillResoProfile(const GFW::CorrConfig& corrconf, const ConstStr& tarName, const double& cent, TAxis* partaxis) + void fillResoProfile(const GFW::CorrConfig& corrconf, const ConstStr& tarName, const double cent, TAxis* partaxis) { double dnx, val; if (!corrconf.pTDif) { @@ -684,6 +684,8 @@ struct ResonancesGfwFlow { double wacc = 1; double cent = collision.centFT0C(); double vtxz = collision.posZ(); + double phi = mom.Phi(); + phi = RecoDecay::constrainAngle(phi, 0.0, 1); // constrain azimuthal angle to [0,2pi] if ((cfgUseWeightPhiEtaVtxz && cfgUseWeightPhiPtCent) || (cfgUseWeightPhiEtaPt && cfgUseWeightPhiPtCent) || (cfgUseWeightPhiEtaVtxz && cfgUseWeightPhiEtaPt)) { LOGF(fatal, "Only one of the three weight options can be used at a time"); @@ -694,11 +696,11 @@ struct ResonancesGfwFlow { return 1; } if (cfgUseWeightPhiEtaVtxz) - wacc = mAcceptance[pid_index_reso]->getNUA(mom.Phi(), mom.Eta(), vtxz); + wacc = mAcceptance[pid_index_reso]->getNUA(phi, mom.Eta(), vtxz); if (cfgUseWeightPhiPtCent) - wacc = mAcceptance[pid_index_reso]->getNUA(mom.Phi(), mom.Pt(), cent); + wacc = mAcceptance[pid_index_reso]->getNUA(phi, mom.Pt(), cent); if (cfgUseWeightPhiEtaPt) - wacc = mAcceptance[pid_index_reso]->getNUA(mom.Phi(), mom.Eta(), mom.Pt()); + wacc = mAcceptance[pid_index_reso]->getNUA(phi, mom.Eta(), mom.Pt()); } return wacc; } @@ -761,19 +763,41 @@ struct ResonancesGfwFlow { histos.fill(HIST("KaminusTPC"), partminus.pt(), partminus.tpcNSigmaKa()); histos.fill(HIST("KaminusTOF"), partminus.pt(), partminus.tofNSigmaKa()); + std::array, 2> ptarr = {{{partplus.px(), partplus.py(), partplus.pz()}, {partminus.px(), partminus.py(), partminus.pz()}}}; + std::array massarr = {plusmass, plusmass}; + + // Calculation using RecoDecay + double invMassRD = RecoDecay::m2(ptarr, massarr); + double ptRD = std::sqrt(RecoDecay::sumOfSquares(partplus.pt(), partminus.pt())); + + // Calculation using ROOT vectors plusdaug = ROOT::Math::PxPyPzMVector(partplus.px(), partplus.py(), partplus.pz(), plusmass); minusdaug = ROOT::Math::PxPyPzMVector(partminus.px(), partminus.py(), partminus.pz(), plusmass); mom = plusdaug + minusdaug; double pt = mom.Pt(); double invMass = mom.M(); + double phi = mom.Phi(); bool withinPtPOI = (cfgCutPtPOIMin < pt) && (pt < cfgCutPtPOIMax); // within POI pT range bool withinPtRef = (cfgCutPtMin < pt) && (pt < cfgCutPtMax); + phi = RecoDecay::constrainAngle(phi, 0.0, 1); // constrain azimuthal angle to [0,2pi] + if (std::abs(mom.Rapidity()) < cfgRapidityCut) { histos.fill(hist, invMass, pt, collision.centFT0C()); - histos.fill(HIST("hPhiPhi"), mom.Phi()); + histos.fill(HIST("hPhiPhi"), phi); histos.fill(HIST("hPhiEta"), mom.Eta()); + histos.fill(HIST("hPhimassSparse_RD"), invMassRD, ptRD, collision.centFT0C()); // fill RecoDecay mass and pt + + // Fill Phi weights + if (cfgOutputNUAWeights && withinPtPOI) { + histos.fill(HIST("NUA/hPhiEtaVtxz_phi"), phi, mom.Eta(), collision.posZ()); + histos.fill(HIST("NUA/hPhiPtCent_phi"), phi, pt, collision.centFT0C()); + histos.fill(HIST("NUA/hPhiEtaPt_phi"), phi, mom.Eta(), pt); + } + + double weff = 1; + double waccPOI = getAcceptancePhi(mom, collision, Phi); // Fill Phi weights if (cfgOutputNUAWeights && withinPtPOI) { @@ -786,9 +810,9 @@ struct ResonancesGfwFlow { double waccPOI = getAcceptancePhi(mom, collision, Phi); if (withinPtPOI) - fGFW->Fill(mom.Eta(), ((fPtAxis->FindBin(pt) - 1) * fPhiMassAxis->GetNbins()) + (fPhiMassAxis->FindBin(invMass) - 1), mom.Phi(), weff * waccPOI, 2); + fGFW->Fill(mom.Eta(), ((fPtAxis->FindBin(pt) - 1) * fPhiMassAxis->GetNbins()) + (fPhiMassAxis->FindBin(invMass) - 1), phi, weff * waccPOI, 2); if (withinPtPOI && withinPtRef) - fGFW->Fill(mom.Eta(), ((fPtAxis->FindBin(pt) - 1) * fPhiMassAxis->GetNbins()) + (fPhiMassAxis->FindBin(invMass) - 1), mom.Phi(), weff * waccPOI, 32); + fGFW->Fill(mom.Eta(), ((fPtAxis->FindBin(pt) - 1) * fPhiMassAxis->GetNbins()) + (fPhiMassAxis->FindBin(invMass) - 1), phi, weff * waccPOI, 32); } } return; @@ -1087,6 +1111,8 @@ struct ResonancesGfwFlow { double pt = track.pt(); bool withinPtRef = (cfgCutPtMin < pt) && (pt < cfgCutPtMax); + weff = 1; // Initializing weff for each track + if (withinPtRef) if (cfgOutputNUAWeights) fillWeights(track, collision, Ref); @@ -1117,19 +1143,18 @@ struct ResonancesGfwFlow { } for (auto const& v0s : V0s) { - if (cfgUseLambda) { - if (selectionLambda(collision, v0s) == true) - histos.fill(HIST("hLambdas"), 1); - } if (cfgUseK0) { if (selectionK0(collision, v0s) == true) histos.fill(HIST("hK0s"), 1); } + if (cfgUseLambda) { + if (selectionLambda(collision, v0s) == true) + histos.fill(HIST("hLambdas"), 1); + } } // End of v0 loop fillResoProfile(corrconfigs.at(0), HIST("Refc22"), cent, fPhiMassAxis); fillResoProfile(corrconfigs.at(1), HIST("Refc24"), cent, fPhiMassAxis); - if (cfgUsePhi) { fillResoProfile(corrconfigs.at(2), HIST("Phid22Fpt"), cent, fPhiMassAxis); fillResoProfile(corrconfigs.at(3), HIST("Phid24Fpt"), cent, fPhiMassAxis); @@ -1185,7 +1210,6 @@ struct ResonancesGfwFlow { fillProfileBoot3D(corrconfigs.at(16), anLambdaD22BPtBoot[bootId], cent, fLambdaMassAxis); fillProfileBoot3D(corrconfigs.at(17), anLambdaD24BPtBoot[bootId], cent, fLambdaMassAxis); } - } // end of bootstrap condition } // end of process };