From 5a0026bf4274b2040282e90055a4c22961013979 Mon Sep 17 00:00:00 2001 From: miedema-11 <1260971129@qq.com> Date: Thu, 23 Oct 2025 14:52:08 +0800 Subject: [PATCH 1/8] resorted MC --- PWGUD/Tasks/flowCumulantsUpc.cxx | 356 ++++++++++++++++++------------- 1 file changed, 202 insertions(+), 154 deletions(-) diff --git a/PWGUD/Tasks/flowCumulantsUpc.cxx b/PWGUD/Tasks/flowCumulantsUpc.cxx index 5e2d64b6d94..9065dc8d958 100644 --- a/PWGUD/Tasks/flowCumulantsUpc.cxx +++ b/PWGUD/Tasks/flowCumulantsUpc.cxx @@ -21,6 +21,7 @@ #include "GFWWeights.h" #include "PWGUD/Core/SGSelector.h" +#include "PWGUD/Core/UPCHelpers.h" #include "PWGUD/DataModel/UDTables.h" #include "Common/CCDB/ctpRateFetcher.h" @@ -39,7 +40,10 @@ #include "Framework/runDataProcessing.h" #include +#include "Math/LorentzVector.h" +#include "Math/PxPyPzM4D.h" #include "TList.h" +#include "TMath.h" #include "TVector3.h" #include #include @@ -48,6 +52,7 @@ #include #include +#include #include #include #include @@ -61,6 +66,11 @@ using namespace o2::framework::expressions; struct FlowCumulantsUpc { + // resort + // Preslice perCollision = aod::track::collisionId; + PresliceUnsorted perMcCollision = o2::aod::udmcparticle::udMcCollisionId; + // Preslice perCol = aod::track::collisionId; + O2_DEFINE_CONFIGURABLE(cfgCutVertex, float, 10.0f, "Accepted z-vertex range") O2_DEFINE_CONFIGURABLE(cfgCentEstimator, int, 0, "0:FT0C; 1:FT0CVariant1; 2:FT0M; 3:FT0A") O2_DEFINE_CONFIGURABLE(cfgCentFT0CMin, float, 0.0f, "Minimum centrality (FT0C) to cut events in filter") @@ -117,12 +127,12 @@ struct FlowCumulantsUpc { ConfigurableAxis axisDCAxy{"axisDCAxy", {200, -1, 1}, "DCA_{xy} (cm)"}; // Added UPC Cuts - SGSelector sgSelector; - Configurable cfgCutFV0{"cfgCutFV0", 50., "FV0A threshold"}; - Configurable cfgCutFT0A{"cfgCutFT0A", 150., "FT0A threshold"}; - Configurable cfgCutFT0C{"cfgCutFT0C", 50., "FT0C threshold"}; - Configurable cfgCutZDC{"cfgCutZDC", 10., "ZDC threshold"}; - Configurable cfgGapSideSelection{"cfgGapSideSelection", 2, "gap selection"}; + // SGSelector sgSelector; + // Configurable cfgCutFV0{"cfgCutFV0", 50., "FV0A threshold"}; + // Configurable cfgCutFT0A{"cfgCutFT0A", 150., "FT0A threshold"}; + // Configurable cfgCutFT0C{"cfgCutFT0C", 50., "FT0C threshold"}; + // Configurable cfgCutZDC{"cfgCutZDC", 10., "ZDC threshold"}; + // Configurable cfgGapSideSelection{"cfgGapSideSelection", 2, "gap selection"}; // Corrections TH1D* mEfficiency = nullptr; @@ -166,9 +176,10 @@ struct FlowCumulantsUpc { TH2* gCurrentHadronicRate; // + // using MCparticles = aod::UDMcParticles::iterator; using UdTracks = soa::Join; using UdTracksFull = soa::Join; - using UDCollisionsFull = soa::Join; + // using UDCollisionsFull = soa::Join; // Track selection TrackSelection myTrackSel; @@ -231,6 +242,10 @@ struct FlowCumulantsUpc { registry.add("hVtxZ", "Vexter Z distribution", {HistType::kTH1D, {axisVertex}}); registry.add("hMult", "Multiplicity distribution", {HistType::kTH1D, {{3000, 0.5, 3000.5}}}); std::string hCentTitle = "Centrality distribution, Estimator " + std::to_string(cfgCentEstimator); + registry.add("hCentMC", hCentTitle.c_str(), {HistType::kTH1D, {{90, 0, 90}}}); + registry.add("hVtxZMC", "Vexter Z distribution", {HistType::kTH1D, {axisVertex}}); + registry.add("hMultMC", "Multiplicity distribution", {HistType::kTH1D, {{3000, 0.5, 3000.5}}}); + registry.add("numberOfTracksMC", "Number of tracks per event", {HistType::kTH1D, {{1000, 0, 1000}}}); registry.add("hCent", hCentTitle.c_str(), {HistType::kTH1D, {{90, 0, 90}}}); if (!cfgUseSmallMemory) { registry.add("BeforeSel8_globalTracks_centT0C", "before sel8;Centrality T0C;mulplicity global tracks", {HistType::kTH2D, {axisCentForQA, axisNch}}); @@ -267,6 +282,10 @@ struct FlowCumulantsUpc { registry.add("hDCAxy", "DCAxy after cuts; DCAxy (cm); Pt", {HistType::kTH2D, {{200, -0.5, 0.5}, {200, 0, 5}}}); registry.add("hTrackCorrection2d", "Correlation table for number of tracks table; uncorrected track; corrected track", {HistType::kTH2D, {axisNch, axisNch}}); + registry.add("hPxMc", "Px distribution of MC truth particles", {HistType::kTH1D, {{100, -10, 10}}}); + registry.add("hPyMc", "Px distribution of MC truth particles", {HistType::kTH1D, {{100, -10, 10}}}); + registry.add("hPzMc", "Px distribution of MC truth particles", {HistType::kTH1D, {{100, -10, 10}}}); + registry.add("hweightMc", "weight distribution of MC truth particles", {HistType::kTH1D, {{100, 0, 1}}}); registry.add("hPhiMC", "#phi distribution", {HistType::kTH1D, {axisPhi}}); registry.add("hPhiWeightedMC", "corrected #phi distribution", {HistType::kTH1D, {axisPhi}}); registry.add("hEtaMC", "#eta distribution", {HistType::kTH1D, {axisEta}}); @@ -279,6 +298,7 @@ struct FlowCumulantsUpc { registry.add("hnTPCCrossedRowMC", "Number of crossed TPC Rows", {HistType::kTH1D, {{100, 40, 180}}}); registry.add("hDCAzMC", "DCAz after cuts; DCAz (cm); Pt", {HistType::kTH2D, {{200, -0.5, 0.5}, {200, 0, 5}}}); registry.add("hDCAxyMC", "DCAxy after cuts; DCAxy (cm); Pt", {HistType::kTH2D, {{200, -0.5, 0.5}, {200, 0, 5}}}); + registry.add("eventCounterMC", "Number of MC Event;; Count", {HistType::kTH1D, {{5, 0, 5}}}); registry.add("hTrackCorrection2dMC", "Correlation table for number of tracks table; uncorrected track; corrected track", {HistType::kTH2D, {axisNch, axisNch}}); o2::framework::AxisSpec axis = axisPt; @@ -412,6 +432,7 @@ struct FlowCumulantsUpc { fGFWMC->AddRegion("olN", -0.8, -0.4, 1 + fPtAxis->GetNbins(), 4); fGFWMC->AddRegion("olN10", -0.8, -0.5, 1 + fPtAxis->GetNbins(), 4); fGFWMC->AddRegion("olfull", -0.8, 0.8, 1 + fPtAxis->GetNbins(), 4); + fGFWMC->CreateRegions(); corrconfigs.push_back(fGFW->GetCorrelatorConfig("full {2 -2}", "ChFull22", kFALSE)); corrconfigs.push_back(fGFW->GetCorrelatorConfig("full {3 -3}", "ChFull32", kFALSE)); @@ -833,171 +854,198 @@ struct FlowCumulantsUpc { gCurrentHadronicRate = gHadronicRate[mRunNumber]; } - void process(UDCollisionsFull::iterator const& collision, UdTracksFull const& tracks) - { - registry.fill(HIST("hEventCount"), 0.5); - int gapSide = collision.gapSide(); - constexpr int kGapSideSelection = 0; - constexpr int kGapSideOppositeSelection = 2; - if (gapSide < kGapSideSelection || gapSide > kGapSideOppositeSelection) { - return; - } + // void process(UDCollisionsFull::iterator const& collision, UdTracksFull const& tracks) + // { + // std::cout << "Processing collision============================================== " << std::endl; - int trueGapSide = sgSelector.trueGap(collision, cfgCutFV0, cfgCutFT0A, cfgCutFT0C, cfgCutZDC); - gapSide = trueGapSide; - if (gapSide == cfgGapSideSelection) { - return; - } - registry.fill(HIST("hEventCount"), 1.5); - float cent = 100; - float lRandom = fRndm->Rndm(); - float vtxz = collision.posZ(); - registry.fill(HIST("hVtxZ"), vtxz); - registry.fill(HIST("hMult"), tracks.size()); - registry.fill(HIST("hCent"), cent); - fGFW->Clear(); - - // // track weights - float weff = 1, wacc = 1; - double nTracksCorrected = 0; - float independent = cent; - if (cfgUseNch) { - independent = static_cast(tracks.size()); - } - - for (const auto& track : tracks) { - if (!trackSelected(track)) - continue; - auto momentum = std::array{track.px(), track.py(), track.pz()}; - double pt = RecoDecay::pt(momentum); - double phi = RecoDecay::phi(momentum); - double eta = RecoDecay::eta(momentum); - bool withinPtPOI = (cfgCutPtPOIMin < pt) && (pt < cfgCutPtPOIMax); // within POI pT range - bool withinPtRef = (cfgCutPtRefMin < pt) && (pt < cfgCutPtRefMax); // within RF pT range - if (cfgOutputNUAWeights) { - if (cfgOutputNUAWeightsRefPt) { - if (withinPtRef) { - fWeights->fill(phi, eta, vtxz, pt, cent, 0); - } - } else { - fWeights->fill(phi, eta, vtxz, pt, cent, 0); - } - } - if (!setCurrentParticleWeights(weff, wacc, phi, eta, pt, vtxz)) { - continue; - } - registry.fill(HIST("hPt"), track.pt()); - if (withinPtRef) { - registry.fill(HIST("hPhi"), phi); - registry.fill(HIST("hPhiWeighted"), phi, wacc); - registry.fill(HIST("hEta"), eta); - registry.fill(HIST("hPtRef"), pt); - registry.fill(HIST("hDCAz"), track.dcaZ(), track.pt()); - registry.fill(HIST("hDCAxy"), track.dcaXY(), track.pt()); - nTracksCorrected += weff; - } - if (withinPtRef) { - fGFW->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 1); - } - if (withinPtPOI) { - fGFW->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 2); - } - if (withinPtPOI && withinPtRef) { - fGFW->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 4); - } - } - registry.fill(HIST("hTrackCorrection2d"), tracks.size(), nTracksCorrected); + // registry.fill(HIST("hEventCount"), 0.5); + // int gapSide = collision.gapSide(); + // constexpr int kGapSideSelection = 0; + // constexpr int kGapSideOppositeSelection = 2; + // if (gapSide < kGapSideSelection || gapSide > kGapSideOppositeSelection) { + // return; + // } - // Filling Flow Container - for (uint l_ind = 0; l_ind < corrconfigs.size(); l_ind++) { - fillFC(corrconfigs.at(l_ind), independent, lRandom); - } - } - PROCESS_SWITCH(FlowCumulantsUpc, process, "process", true); + // int trueGapSide = sgSelector.trueGap(collision, cfgCutFV0, cfgCutFT0A, cfgCutFT0C, cfgCutZDC); + // gapSide = trueGapSide; + // if (gapSide == cfgGapSideSelection) { + // return; + // } + // registry.fill(HIST("hEventCount"), 1.5); + // float cent = 100; + // float lRandom = fRndm->Rndm(); + // float vtxz = collision.posZ(); + // registry.fill(HIST("hVtxZ"), vtxz); + // registry.fill(HIST("hMult"), tracks.size()); + // registry.fill(HIST("hCent"), cent); + // fGFW->Clear(); + + // // // track weights + // float weff = 1, wacc = 1; + // double nTracksCorrected = 0; + // float independent = cent; + // if (cfgUseNch) { + // independent = static_cast(tracks.size()); + // } + + // for (const auto& track : tracks) { + // if (!trackSelected(track)) + // continue; + // auto momentum = std::array{track.px(), track.py(), track.pz()}; + // double pt = RecoDecay::pt(momentum); + // double phi = RecoDecay::phi(momentum); + // double eta = RecoDecay::eta(momentum); + // bool withinPtPOI = (cfgCutPtPOIMin < pt) && (pt < cfgCutPtPOIMax); // within POI pT range + // bool withinPtRef = (cfgCutPtRefMin < pt) && (pt < cfgCutPtRefMax); // within RF pT range + // if (cfgOutputNUAWeights) { + // if (cfgOutputNUAWeightsRefPt) { + // if (withinPtRef) { + // fWeights->fill(phi, eta, vtxz, pt, cent, 0); + // } + // } else { + // fWeights->fill(phi, eta, vtxz, pt, cent, 0); + // } + // } + // if (!setCurrentParticleWeights(weff, wacc, phi, eta, pt, vtxz)) { + // continue; + // } + // registry.fill(HIST("hPt"), track.pt()); + // if (withinPtRef) { + // registry.fill(HIST("hPhi"), phi); + // registry.fill(HIST("hPhiWeighted"), phi, wacc); + // registry.fill(HIST("hEta"), eta); + // registry.fill(HIST("hPtRef"), pt); + // registry.fill(HIST("hDCAz"), track.dcaZ(), track.pt()); + // registry.fill(HIST("hDCAxy"), track.dcaXY(), track.pt()); + // nTracksCorrected += weff; + // } + // if (withinPtRef) { + // fGFW->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 1); + // } + // if (withinPtPOI) { + // fGFW->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 2); + // } + // if (withinPtPOI && withinPtRef) { + // fGFW->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 4); + // } + // } + // registry.fill(HIST("hTrackCorrection2d"), tracks.size(), nTracksCorrected); + + // // Filling Flow Container + // for (uint l_ind = 0; l_ind < corrconfigs.size(); l_ind++) { + // fillFC(corrconfigs.at(l_ind), independent, lRandom); + // } + // } + // PROCESS_SWITCH(FlowCumulantsUpc, process, "process", false); //----------------------------------------------------------------------------------------------------------------------- - void processSim(aod::UDMcCollision const& mcCollision, aod::UDMcParticles const& mcParticles) + void processSim(aod::UDMcCollisions const& mcCollisions, aod::UDMcParticles const& mcParticles) { - registry.fill(HIST("eventCounterMC"), 0.5); + std::cout << "Processing MC collision======================== " << std::endl; - registry.fill(HIST("hEventCount"), 1.5); - float cent = 100; - float vtxz = mcCollision.posZ(); - registry.fill(HIST("hVtxZMC"), vtxz); - registry.fill(HIST("hMultMC"), mcParticles.size()); - registry.fill(HIST("hCentMC"), cent); + // std::cout << mcParicles::IndexudmcCollisions.is_sorted(); - auto massPion = o2::constants::physics::MassPionCharged; - registry.fill(HIST("numberOfTracksMC"), mcParticles.size()); - // LOGF(info, "New event! mcParticles.size() = %d", mcParticles.size()); + for (auto& mcCollision : mcCollisions) { + auto groupedUDMcParticles = mcParticles.sliceBy(perMcCollision, mcCollision.globalIndex()); + registry.fill(HIST("eventCounterMC"), 0.5); - float lRandomMc = fRndmMc->Rndm(); - fGFWMC->Clear(); + // registry.fill(HIST("hEventCount"), 1.5); + float cent = 100; + float vtxz = 0; - // // track weights - float weff = 1, wacc = 1; - double nTracksCorrected = 0; - float independent = cent; - if (cfgUseNch) { - independent = static_cast(mcParticles.size()); - } + vtxz = mcCollision.posZ(); + registry.fill(HIST("hVtxZMC"), vtxz); + registry.fill(HIST("eventCounterMC"), mcCollision.size()); + registry.fill(HIST("hMultMC"), groupedUDMcParticles.size()); + registry.fill(HIST("hCentMC"), cent); - for (const auto& mcParticle : mcParticles) { - if (!mcParticle.isPhysicalPrimary()) - continue; - std::array momentum = {mcParticle.px(), mcParticle.py(), mcParticle.pz()}; - double energy = std::sqrt(momentum[0] * momentum[0] + momentum[1] * momentum[1] + momentum[2] * momentum[2] + massPion * massPion); - ROOT::Math::LorentzVector> protoMC(momentum[0], momentum[1], momentum[2], energy); - constexpr double kEtaCut = 0.8; - constexpr double kPtCut = 0.1; - if (!(std::fabs(protoMC.Eta()) < kEtaCut && protoMC.Pt() > kPtCut)) { - continue; + auto massPion = o2::constants::physics::MassPionCharged; + registry.fill(HIST("numberOfTracksMC"), groupedUDMcParticles.size()); + // LOGF(info, "New event! mcParticles.size() = %d", mcParticles.size()); + + float lRandomMc = fRndmMc->Rndm(); + fGFWMC->Clear(); + + // // track weights + float weff = 1, wacc = 1; + double nTracksCorrected = 0; + float independent = cent; + if (cfgUseNch) { + independent = static_cast(groupedUDMcParticles.size()); } - // auto momentum = std::array{mcParticle.px(), mcParticle.py(), mcParticle.pz()}; - double pt = RecoDecay::pt(momentum); - double phi = RecoDecay::phi(momentum); - double eta = RecoDecay::eta(momentum); - bool withinPtPOI = (cfgCutPtPOIMin < pt) && (pt < cfgCutPtPOIMax); // within POI pT range - bool withinPtRef = (cfgCutPtRefMin < pt) && (pt < cfgCutPtRefMax); // within RF pT range - if (cfgOutputNUAWeights) { - if (cfgOutputNUAWeightsRefPt) { - if (withinPtRef) { + + std::cout << "mcParticles.size() = " << groupedUDMcParticles.size() << std::endl; + + for (const auto& mcParticle : groupedUDMcParticles) { + + std::cout << "filling mc particle px: " << mcParticle.px() << ", py: " << mcParticle.py() << ", pz: " << mcParticle.pz() << std::endl; + + // output information from mcparticles + registry.fill(HIST("hPxMc"), mcParticle.px()); + registry.fill(HIST("hPyMc"), mcParticle.py()); + registry.fill(HIST("hPzMc"), mcParticle.pz()); + registry.fill(HIST("hweightMc"), mcParticle.weight()); + + if (!mcParticle.isPhysicalPrimary()) { + std::cout << "mcParticle.isPhysicalPrimary() = " << mcParticle.isPhysicalPrimary() << std::endl; + continue; + } + std::array momentum = {mcParticle.px(), mcParticle.py(), mcParticle.pz()}; + double energy = std::sqrt(momentum[0] * momentum[0] + momentum[1] * momentum[1] + momentum[2] * momentum[2] + massPion * massPion); + ROOT::Math::LorentzVector> protoMC(momentum[0], momentum[1], momentum[2], energy); + constexpr double kEtaCut = 0.8; + constexpr double kPtCut = 0.1; + if (!(std::fabs(protoMC.Eta()) < kEtaCut && protoMC.Pt() > kPtCut)) { + std::cout << "protoMC.Eta() = " << protoMC.Eta() << ", protoMC.Pt() = " << protoMC.Pt() << std::endl; + continue; + } + // auto momentum = std::array{mcParticle.px(), mcParticle.py(), mcParticle.pz()}; + double pt = RecoDecay::pt(momentum); + double phi = RecoDecay::phi(momentum); + double eta = RecoDecay::eta(momentum); + bool withinPtPOI = (cfgCutPtPOIMin < pt) && (pt < cfgCutPtPOIMax); // within POI pT range + bool withinPtRef = (cfgCutPtRefMin < pt) && (pt < cfgCutPtRefMax); // within RF pT range + if (cfgOutputNUAWeights) { + if (cfgOutputNUAWeightsRefPt) { + if (withinPtRef) { + fWeightsMc->fill(phi, eta, vtxz, pt, cent, 0); + } + } else { fWeightsMc->fill(phi, eta, vtxz, pt, cent, 0); } - } else { - fWeightsMc->fill(phi, eta, vtxz, pt, cent, 0); + } + if (!setCurrentParticleWeights(weff, wacc, phi, eta, pt, vtxz)) { + continue; + } + + if (withinPtRef) { + registry.fill(HIST("hPhiMC"), phi); + registry.fill(HIST("hPhiWeightedMC"), phi, wacc); + registry.fill(HIST("hEtaMC"), eta); + registry.fill(HIST("hPtRefMC"), pt); + // registry.fill(HIST("hDCAzMC"), track.dcaZ(), track.pt()); + // registry.fill(HIST("hDCAxyMC"), track.dcaXY(), track.pt()); + nTracksCorrected += weff; + } + if (withinPtRef) { + fGFWMC->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 1); + } + if (withinPtPOI) { + fGFWMC->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 2); + } + if (withinPtPOI && withinPtRef) { + fGFWMC->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 4); } } - if (!setCurrentParticleWeights(weff, wacc, phi, eta, pt, vtxz)) { - continue; - } - if (withinPtRef) { - registry.fill(HIST("hPhiMC"), phi); - registry.fill(HIST("hPhiWeightedMC"), phi, wacc); - registry.fill(HIST("hEtaMC"), eta); - registry.fill(HIST("hPtRefMC"), pt); - // registry.fill(HIST("hDCAzMC"), track.dcaZ(), track.pt()); - // registry.fill(HIST("hDCAxyMC"), track.dcaXY(), track.pt()); - nTracksCorrected += weff; - } - if (withinPtRef) { - fGFWMC->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 1); - } - if (withinPtPOI) { - fGFWMC->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 2); - } - if (withinPtPOI && withinPtRef) { - fGFWMC->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 4); + registry.fill(HIST("hTrackCorrection2dMC"), mcParticles.size(), nTracksCorrected); + // Filling Flow Container + for (uint l_ind = 0; l_ind < corrconfigsmc.size(); l_ind++) { + std::cout << "filling flow container for MC" << std::endl; + fillFCMC(corrconfigsmc.at(l_ind), independent, lRandomMc); } } - registry.fill(HIST("hTrackCorrection2dMC"), mcParticles.size(), nTracksCorrected); - - // Filling Flow Container - for (uint l_ind = 0; l_ind < corrconfigs.size(); l_ind++) { - fillFCMC(corrconfigs.at(l_ind), independent, lRandomMc); - } - PROCESS_SWITCH(FlowCumulantsUpc, processSim, "processSim", false); } + PROCESS_SWITCH(FlowCumulantsUpc, processSim, "processSim", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) From c88abbdb0cf2b3146570609311f23584d94df5e4 Mon Sep 17 00:00:00 2001 From: miedema-11 <1260971129@qq.com> Date: Tue, 2 Dec 2025 01:23:46 +0800 Subject: [PATCH 2/8] modify process switch --- PWGUD/Tasks/flowCumulantsUpc.cxx | 36 ++++++++++++++++++-------------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/PWGUD/Tasks/flowCumulantsUpc.cxx b/PWGUD/Tasks/flowCumulantsUpc.cxx index 9065dc8d958..147325ee9c5 100644 --- a/PWGUD/Tasks/flowCumulantsUpc.cxx +++ b/PWGUD/Tasks/flowCumulantsUpc.cxx @@ -52,7 +52,6 @@ #include #include -#include #include #include #include @@ -127,12 +126,12 @@ struct FlowCumulantsUpc { ConfigurableAxis axisDCAxy{"axisDCAxy", {200, -1, 1}, "DCA_{xy} (cm)"}; // Added UPC Cuts - // SGSelector sgSelector; - // Configurable cfgCutFV0{"cfgCutFV0", 50., "FV0A threshold"}; - // Configurable cfgCutFT0A{"cfgCutFT0A", 150., "FT0A threshold"}; - // Configurable cfgCutFT0C{"cfgCutFT0C", 50., "FT0C threshold"}; - // Configurable cfgCutZDC{"cfgCutZDC", 10., "ZDC threshold"}; - // Configurable cfgGapSideSelection{"cfgGapSideSelection", 2, "gap selection"}; + SGSelector sgSelector; + Configurable cfgCutFV0{"cfgCutFV0", 50., "FV0A threshold"}; + Configurable cfgCutFT0A{"cfgCutFT0A", 150., "FT0A threshold"}; + Configurable cfgCutFT0C{"cfgCutFT0C", 50., "FT0C threshold"}; + Configurable cfgCutZDC{"cfgCutZDC", 10., "ZDC threshold"}; + Configurable cfgGapSideSelection{"cfgGapSideSelection", 2, "gap selection"}; // Corrections TH1D* mEfficiency = nullptr; @@ -432,7 +431,6 @@ struct FlowCumulantsUpc { fGFWMC->AddRegion("olN", -0.8, -0.4, 1 + fPtAxis->GetNbins(), 4); fGFWMC->AddRegion("olN10", -0.8, -0.5, 1 + fPtAxis->GetNbins(), 4); fGFWMC->AddRegion("olfull", -0.8, 0.8, 1 + fPtAxis->GetNbins(), 4); - fGFWMC->CreateRegions(); corrconfigs.push_back(fGFW->GetCorrelatorConfig("full {2 -2}", "ChFull22", kFALSE)); corrconfigs.push_back(fGFW->GetCorrelatorConfig("full {3 -3}", "ChFull32", kFALSE)); @@ -526,9 +524,11 @@ struct FlowCumulantsUpc { } LOGF(info, "%d: %s %s", i, userDefineGFWCorr.at(i).c_str(), userDefineGFWName.at(i).c_str()); corrconfigs.push_back(fGFW->GetCorrelatorConfig(userDefineGFWCorr.at(i).c_str(), userDefineGFWName.at(i).c_str(), kFALSE)); + corrconfigsmc.push_back(fGFWMC->GetCorrelatorConfig(userDefineGFWCorr.at(i).c_str(), userDefineGFWName.at(i).c_str(), kFALSE)); } } fGFW->CreateRegions(); + fGFWMC->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); @@ -644,6 +644,8 @@ struct FlowCumulantsUpc { return; } + // ... 其余代码保持不变 ... + void loadCorrections(uint64_t timestamp, int runNumber) { if (correctionsLoaded) { @@ -941,16 +943,16 @@ struct FlowCumulantsUpc { //----------------------------------------------------------------------------------------------------------------------- void processSim(aod::UDMcCollisions const& mcCollisions, aod::UDMcParticles const& mcParticles) { - std::cout << "Processing MC collision======================== " << std::endl; + LOG(info) << "Processing MC collision======================== " << std::endl; // std::cout << mcParicles::IndexudmcCollisions.is_sorted(); - for (auto& mcCollision : mcCollisions) { + for (const auto& mcCollision : mcCollisions) { auto groupedUDMcParticles = mcParticles.sliceBy(perMcCollision, mcCollision.globalIndex()); registry.fill(HIST("eventCounterMC"), 0.5); // registry.fill(HIST("hEventCount"), 1.5); - float cent = 100; + float cent = 50; float vtxz = 0; vtxz = mcCollision.posZ(); @@ -974,11 +976,11 @@ struct FlowCumulantsUpc { independent = static_cast(groupedUDMcParticles.size()); } - std::cout << "mcParticles.size() = " << groupedUDMcParticles.size() << std::endl; + LOG(info) << "mcParticles.size() = " << groupedUDMcParticles.size() << std::endl; for (const auto& mcParticle : groupedUDMcParticles) { - std::cout << "filling mc particle px: " << mcParticle.px() << ", py: " << mcParticle.py() << ", pz: " << mcParticle.pz() << std::endl; + LOG(info) << "filling mc particle px: " << mcParticle.px() << ", py: " << mcParticle.py() << ", pz: " << mcParticle.pz() << std::endl; // output information from mcparticles registry.fill(HIST("hPxMc"), mcParticle.px()); @@ -987,16 +989,17 @@ struct FlowCumulantsUpc { registry.fill(HIST("hweightMc"), mcParticle.weight()); if (!mcParticle.isPhysicalPrimary()) { - std::cout << "mcParticle.isPhysicalPrimary() = " << mcParticle.isPhysicalPrimary() << std::endl; + LOG(info) << "mcParticle.isPhysicalPrimary() = " << mcParticle.isPhysicalPrimary() << std::endl; continue; } + std::array momentum = {mcParticle.px(), mcParticle.py(), mcParticle.pz()}; double energy = std::sqrt(momentum[0] * momentum[0] + momentum[1] * momentum[1] + momentum[2] * momentum[2] + massPion * massPion); ROOT::Math::LorentzVector> protoMC(momentum[0], momentum[1], momentum[2], energy); constexpr double kEtaCut = 0.8; constexpr double kPtCut = 0.1; if (!(std::fabs(protoMC.Eta()) < kEtaCut && protoMC.Pt() > kPtCut)) { - std::cout << "protoMC.Eta() = " << protoMC.Eta() << ", protoMC.Pt() = " << protoMC.Pt() << std::endl; + LOG(info) << "protoMC.Eta() = " << protoMC.Eta() << ", protoMC.Pt() = " << protoMC.Pt() << std::endl; continue; } // auto momentum = std::array{mcParticle.px(), mcParticle.py(), mcParticle.pz()}; @@ -1036,11 +1039,12 @@ struct FlowCumulantsUpc { if (withinPtPOI && withinPtRef) { fGFWMC->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 4); } + LOG(info) << "successfully filled" << std::endl; } registry.fill(HIST("hTrackCorrection2dMC"), mcParticles.size(), nTracksCorrected); // Filling Flow Container for (uint l_ind = 0; l_ind < corrconfigsmc.size(); l_ind++) { - std::cout << "filling flow container for MC" << std::endl; + LOG(info) << "filling flow container for MC" << std::endl; fillFCMC(corrconfigsmc.at(l_ind), independent, lRandomMc); } } From 43f7419b2c93374b55f9f6e18f74cd2f9727f24b Mon Sep 17 00:00:00 2001 From: miedema-11 <1260971129@qq.com> Date: Fri, 5 Dec 2025 11:23:30 +0800 Subject: [PATCH 3/8] decrease log output --- PWGUD/Tasks/flowCumulantsUpc.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGUD/Tasks/flowCumulantsUpc.cxx b/PWGUD/Tasks/flowCumulantsUpc.cxx index 147325ee9c5..8b0923a7e99 100644 --- a/PWGUD/Tasks/flowCumulantsUpc.cxx +++ b/PWGUD/Tasks/flowCumulantsUpc.cxx @@ -1039,12 +1039,12 @@ struct FlowCumulantsUpc { if (withinPtPOI && withinPtRef) { fGFWMC->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 4); } - LOG(info) << "successfully filled" << std::endl; + // LOG(info) << "successfully filled" << std::endl; } registry.fill(HIST("hTrackCorrection2dMC"), mcParticles.size(), nTracksCorrected); // Filling Flow Container for (uint l_ind = 0; l_ind < corrconfigsmc.size(); l_ind++) { - LOG(info) << "filling flow container for MC" << std::endl; + // LOG(info) << "filling flow container for MC" << std::endl; fillFCMC(corrconfigsmc.at(l_ind), independent, lRandomMc); } } From fba50dabfd4f840baa0e907f8b02c65b87f7078a Mon Sep 17 00:00:00 2001 From: miedema-11 <1260971129@qq.com> Date: Mon, 8 Dec 2025 14:07:00 +0800 Subject: [PATCH 4/8] decrease more log --- PWGUD/Tasks/flowCumulantsUpc.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/PWGUD/Tasks/flowCumulantsUpc.cxx b/PWGUD/Tasks/flowCumulantsUpc.cxx index 8b0923a7e99..b97b8e26edc 100644 --- a/PWGUD/Tasks/flowCumulantsUpc.cxx +++ b/PWGUD/Tasks/flowCumulantsUpc.cxx @@ -976,11 +976,11 @@ struct FlowCumulantsUpc { independent = static_cast(groupedUDMcParticles.size()); } - LOG(info) << "mcParticles.size() = " << groupedUDMcParticles.size() << std::endl; + // LOG(info) << "mcParticles.size() = " << groupedUDMcParticles.size() << std::endl; for (const auto& mcParticle : groupedUDMcParticles) { - LOG(info) << "filling mc particle px: " << mcParticle.px() << ", py: " << mcParticle.py() << ", pz: " << mcParticle.pz() << std::endl; + // LOG(info) << "filling mc particle px: " << mcParticle.px() << ", py: " << mcParticle.py() << ", pz: " << mcParticle.pz() << std::endl; // output information from mcparticles registry.fill(HIST("hPxMc"), mcParticle.px()); @@ -989,7 +989,7 @@ struct FlowCumulantsUpc { registry.fill(HIST("hweightMc"), mcParticle.weight()); if (!mcParticle.isPhysicalPrimary()) { - LOG(info) << "mcParticle.isPhysicalPrimary() = " << mcParticle.isPhysicalPrimary() << std::endl; + // LOG(info) << "mcParticle.isPhysicalPrimary() = " << mcParticle.isPhysicalPrimary() << std::endl; continue; } @@ -999,7 +999,7 @@ struct FlowCumulantsUpc { constexpr double kEtaCut = 0.8; constexpr double kPtCut = 0.1; if (!(std::fabs(protoMC.Eta()) < kEtaCut && protoMC.Pt() > kPtCut)) { - LOG(info) << "protoMC.Eta() = " << protoMC.Eta() << ", protoMC.Pt() = " << protoMC.Pt() << std::endl; + // LOG(info) << "protoMC.Eta() = " << protoMC.Eta() << ", protoMC.Pt() = " << protoMC.Pt() << std::endl; continue; } // auto momentum = std::array{mcParticle.px(), mcParticle.py(), mcParticle.pz()}; From a85669d5add796fed6c750e46c4c603799082314 Mon Sep 17 00:00:00 2001 From: miedema-11 <1260971129@qq.com> Date: Sat, 17 Jan 2026 22:27:36 +0800 Subject: [PATCH 5/8] add default cut --- PWGUD/Tasks/flowCumulantsUpc.cxx | 396 +++++++++++++++---------------- 1 file changed, 194 insertions(+), 202 deletions(-) diff --git a/PWGUD/Tasks/flowCumulantsUpc.cxx b/PWGUD/Tasks/flowCumulantsUpc.cxx index b97b8e26edc..84575df58f3 100644 --- a/PWGUD/Tasks/flowCumulantsUpc.cxx +++ b/PWGUD/Tasks/flowCumulantsUpc.cxx @@ -21,7 +21,6 @@ #include "GFWWeights.h" #include "PWGUD/Core/SGSelector.h" -#include "PWGUD/Core/UPCHelpers.h" #include "PWGUD/DataModel/UDTables.h" #include "Common/CCDB/ctpRateFetcher.h" @@ -40,10 +39,7 @@ #include "Framework/runDataProcessing.h" #include -#include "Math/LorentzVector.h" -#include "Math/PxPyPzM4D.h" #include "TList.h" -#include "TMath.h" #include "TVector3.h" #include #include @@ -65,11 +61,6 @@ using namespace o2::framework::expressions; struct FlowCumulantsUpc { - // resort - // Preslice perCollision = aod::track::collisionId; - PresliceUnsorted perMcCollision = o2::aod::udmcparticle::udMcCollisionId; - // Preslice perCol = aod::track::collisionId; - O2_DEFINE_CONFIGURABLE(cfgCutVertex, float, 10.0f, "Accepted z-vertex range") O2_DEFINE_CONFIGURABLE(cfgCentEstimator, int, 0, "0:FT0C; 1:FT0CVariant1; 2:FT0M; 3:FT0A") O2_DEFINE_CONFIGURABLE(cfgCentFT0CMin, float, 0.0f, "Minimum centrality (FT0C) to cut events in filter") @@ -175,10 +166,9 @@ struct FlowCumulantsUpc { TH2* gCurrentHadronicRate; // - // using MCparticles = aod::UDMcParticles::iterator; using UdTracks = soa::Join; using UdTracksFull = soa::Join; - // using UDCollisionsFull = soa::Join; + using UDCollisionsFull = soa::Join; // Track selection TrackSelection myTrackSel; @@ -241,10 +231,6 @@ struct FlowCumulantsUpc { registry.add("hVtxZ", "Vexter Z distribution", {HistType::kTH1D, {axisVertex}}); registry.add("hMult", "Multiplicity distribution", {HistType::kTH1D, {{3000, 0.5, 3000.5}}}); std::string hCentTitle = "Centrality distribution, Estimator " + std::to_string(cfgCentEstimator); - registry.add("hCentMC", hCentTitle.c_str(), {HistType::kTH1D, {{90, 0, 90}}}); - registry.add("hVtxZMC", "Vexter Z distribution", {HistType::kTH1D, {axisVertex}}); - registry.add("hMultMC", "Multiplicity distribution", {HistType::kTH1D, {{3000, 0.5, 3000.5}}}); - registry.add("numberOfTracksMC", "Number of tracks per event", {HistType::kTH1D, {{1000, 0, 1000}}}); registry.add("hCent", hCentTitle.c_str(), {HistType::kTH1D, {{90, 0, 90}}}); if (!cfgUseSmallMemory) { registry.add("BeforeSel8_globalTracks_centT0C", "before sel8;Centrality T0C;mulplicity global tracks", {HistType::kTH2D, {axisCentForQA, axisNch}}); @@ -281,10 +267,6 @@ struct FlowCumulantsUpc { registry.add("hDCAxy", "DCAxy after cuts; DCAxy (cm); Pt", {HistType::kTH2D, {{200, -0.5, 0.5}, {200, 0, 5}}}); registry.add("hTrackCorrection2d", "Correlation table for number of tracks table; uncorrected track; corrected track", {HistType::kTH2D, {axisNch, axisNch}}); - registry.add("hPxMc", "Px distribution of MC truth particles", {HistType::kTH1D, {{100, -10, 10}}}); - registry.add("hPyMc", "Px distribution of MC truth particles", {HistType::kTH1D, {{100, -10, 10}}}); - registry.add("hPzMc", "Px distribution of MC truth particles", {HistType::kTH1D, {{100, -10, 10}}}); - registry.add("hweightMc", "weight distribution of MC truth particles", {HistType::kTH1D, {{100, 0, 1}}}); registry.add("hPhiMC", "#phi distribution", {HistType::kTH1D, {axisPhi}}); registry.add("hPhiWeightedMC", "corrected #phi distribution", {HistType::kTH1D, {axisPhi}}); registry.add("hEtaMC", "#eta distribution", {HistType::kTH1D, {axisEta}}); @@ -297,7 +279,6 @@ struct FlowCumulantsUpc { registry.add("hnTPCCrossedRowMC", "Number of crossed TPC Rows", {HistType::kTH1D, {{100, 40, 180}}}); registry.add("hDCAzMC", "DCAz after cuts; DCAz (cm); Pt", {HistType::kTH2D, {{200, -0.5, 0.5}, {200, 0, 5}}}); registry.add("hDCAxyMC", "DCAxy after cuts; DCAxy (cm); Pt", {HistType::kTH2D, {{200, -0.5, 0.5}, {200, 0, 5}}}); - registry.add("eventCounterMC", "Number of MC Event;; Count", {HistType::kTH1D, {{5, 0, 5}}}); registry.add("hTrackCorrection2dMC", "Correlation table for number of tracks table; uncorrected track; corrected track", {HistType::kTH2D, {axisNch, axisNch}}); o2::framework::AxisSpec axis = axisPt; @@ -524,11 +505,9 @@ struct FlowCumulantsUpc { } LOGF(info, "%d: %s %s", i, userDefineGFWCorr.at(i).c_str(), userDefineGFWName.at(i).c_str()); corrconfigs.push_back(fGFW->GetCorrelatorConfig(userDefineGFWCorr.at(i).c_str(), userDefineGFWName.at(i).c_str(), kFALSE)); - corrconfigsmc.push_back(fGFWMC->GetCorrelatorConfig(userDefineGFWCorr.at(i).c_str(), userDefineGFWName.at(i).c_str(), kFALSE)); } } fGFW->CreateRegions(); - fGFWMC->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); @@ -644,8 +623,6 @@ struct FlowCumulantsUpc { return; } - // ... 其余代码保持不变 ... - void loadCorrections(uint64_t timestamp, int runNumber) { if (correctionsLoaded) { @@ -815,10 +792,10 @@ struct FlowCumulantsUpc { if (!((multNTracksPV < fMultPVCutLow->Eval(centrality)) || (multNTracksPV > fMultPVCutHigh->Eval(centrality)) || (multTrk < fMultCutLow->Eval(centrality)) || (multTrk > fMultCutHigh->Eval(centrality)))) { registry.fill(HIST("hEventCountTentative"), 8.5); } - constexpr int kSigmaCut = 5; - if (!(std::fabs(collision.multFV0A() - fT0AV0AMean->Eval(collision.multFT0A())) > kSigmaCut * fT0AV0ASigma->Eval(collision.multFT0A()))) { - registry.fill(HIST("hEventCountTentative"), 9.5); - } + // constexpr int kSigmaCut = 5; + // if (!(std::fabs(collision.multFV0A() - fT0AV0AMean->Eval(collision.multFT0A())) > kSigmaCut * fT0AV0ASigma->Eval(collision.multFT0A()))) { + // registry.fill(HIST("hEventCountTentative"), 9.5); + // } } template @@ -836,6 +813,15 @@ struct FlowCumulantsUpc { if (!(std::fabs(track.dcaXY()) < dcaLimit)) { return false; } + if (track.tpcNClsFindableMinusCrossedRows() <= 70) { + return false; + } + if (track.itsClusterSizes() <= 5) { + return false; + } + if (track.tpcChi2NCl() >= 4) { + return false; + } return true; } @@ -856,200 +842,206 @@ struct FlowCumulantsUpc { gCurrentHadronicRate = gHadronicRate[mRunNumber]; } - // void process(UDCollisionsFull::iterator const& collision, UdTracksFull const& tracks) - // { - // std::cout << "Processing collision============================================== " << std::endl; - - // registry.fill(HIST("hEventCount"), 0.5); - // int gapSide = collision.gapSide(); - // constexpr int kGapSideSelection = 0; - // constexpr int kGapSideOppositeSelection = 2; - // if (gapSide < kGapSideSelection || gapSide > kGapSideOppositeSelection) { - // return; - // } - - // int trueGapSide = sgSelector.trueGap(collision, cfgCutFV0, cfgCutFT0A, cfgCutFT0C, cfgCutZDC); - // gapSide = trueGapSide; - // if (gapSide == cfgGapSideSelection) { - // return; - // } - // registry.fill(HIST("hEventCount"), 1.5); - // float cent = 100; - // float lRandom = fRndm->Rndm(); - // float vtxz = collision.posZ(); - // registry.fill(HIST("hVtxZ"), vtxz); - // registry.fill(HIST("hMult"), tracks.size()); - // registry.fill(HIST("hCent"), cent); - // fGFW->Clear(); - - // // // track weights - // float weff = 1, wacc = 1; - // double nTracksCorrected = 0; - // float independent = cent; - // if (cfgUseNch) { - // independent = static_cast(tracks.size()); - // } - - // for (const auto& track : tracks) { - // if (!trackSelected(track)) - // continue; - // auto momentum = std::array{track.px(), track.py(), track.pz()}; - // double pt = RecoDecay::pt(momentum); - // double phi = RecoDecay::phi(momentum); - // double eta = RecoDecay::eta(momentum); - // bool withinPtPOI = (cfgCutPtPOIMin < pt) && (pt < cfgCutPtPOIMax); // within POI pT range - // bool withinPtRef = (cfgCutPtRefMin < pt) && (pt < cfgCutPtRefMax); // within RF pT range - // if (cfgOutputNUAWeights) { - // if (cfgOutputNUAWeightsRefPt) { - // if (withinPtRef) { - // fWeights->fill(phi, eta, vtxz, pt, cent, 0); - // } - // } else { - // fWeights->fill(phi, eta, vtxz, pt, cent, 0); - // } - // } - // if (!setCurrentParticleWeights(weff, wacc, phi, eta, pt, vtxz)) { - // continue; - // } - // registry.fill(HIST("hPt"), track.pt()); - // if (withinPtRef) { - // registry.fill(HIST("hPhi"), phi); - // registry.fill(HIST("hPhiWeighted"), phi, wacc); - // registry.fill(HIST("hEta"), eta); - // registry.fill(HIST("hPtRef"), pt); - // registry.fill(HIST("hDCAz"), track.dcaZ(), track.pt()); - // registry.fill(HIST("hDCAxy"), track.dcaXY(), track.pt()); - // nTracksCorrected += weff; - // } - // if (withinPtRef) { - // fGFW->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 1); - // } - // if (withinPtPOI) { - // fGFW->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 2); - // } - // if (withinPtPOI && withinPtRef) { - // fGFW->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 4); - // } - // } - // registry.fill(HIST("hTrackCorrection2d"), tracks.size(), nTracksCorrected); - - // // Filling Flow Container - // for (uint l_ind = 0; l_ind < corrconfigs.size(); l_ind++) { - // fillFC(corrconfigs.at(l_ind), independent, lRandom); - // } - // } - // PROCESS_SWITCH(FlowCumulantsUpc, process, "process", false); - - //----------------------------------------------------------------------------------------------------------------------- - void processSim(aod::UDMcCollisions const& mcCollisions, aod::UDMcParticles const& mcParticles) + template + float getDPhiStar(TTrack const& track1, TTrack const& track2, float radius, int runnum) { - LOG(info) << "Processing MC collision======================== " << std::endl; - - // std::cout << mcParicles::IndexudmcCollisions.is_sorted(); + float charge1 = track1.sign(); + float charge2 = track2.sign(); - for (const auto& mcCollision : mcCollisions) { - auto groupedUDMcParticles = mcParticles.sliceBy(perMcCollision, mcCollision.globalIndex()); - registry.fill(HIST("eventCounterMC"), 0.5); + float phi1 = track1.phi(); + float phi2 = track2.phi(); - // registry.fill(HIST("hEventCount"), 1.5); - float cent = 50; - float vtxz = 0; + float pt1 = track1.pt(); + float pt2 = track2.pt(); - vtxz = mcCollision.posZ(); - registry.fill(HIST("hVtxZMC"), vtxz); - registry.fill(HIST("eventCounterMC"), mcCollision.size()); - registry.fill(HIST("hMultMC"), groupedUDMcParticles.size()); - registry.fill(HIST("hCentMC"), cent); + int fbSign = 1; - auto massPion = o2::constants::physics::MassPionCharged; - registry.fill(HIST("numberOfTracksMC"), groupedUDMcParticles.size()); - // LOGF(info, "New event! mcParticles.size() = %d", mcParticles.size()); + if (runnum >= 544868) { + fbSign = -1; + } - float lRandomMc = fRndmMc->Rndm(); - fGFWMC->Clear(); + float dPhiStar = phi1 - phi2 - charge1 * fbSign * std::asin(0.075 * radius / pt1) + charge2 * fbSign * std::asin(0.075 * radius / pt2); - // // track weights - float weff = 1, wacc = 1; - double nTracksCorrected = 0; - float independent = cent; - if (cfgUseNch) { - independent = static_cast(groupedUDMcParticles.size()); - } + if (dPhiStar > constants::math::PI) + dPhiStar = constants::math::TwoPI - dPhiStar; + if (dPhiStar < -constants::math::PI) + dPhiStar = -constants::math::TwoPI - dPhiStar; - // LOG(info) << "mcParticles.size() = " << groupedUDMcParticles.size() << std::endl; + return dPhiStar; + } - for (const auto& mcParticle : groupedUDMcParticles) { + void process(UDCollisionsFull::iterator const& collision, UdTracksFull const& tracks) + { + registry.fill(HIST("hEventCount"), 0.5); + // if(!eventSelected(collision, tracks.size(), 100.0f)) { + // eventCounterQA(collision, tracks.size(), 100.0f); + // return; + // } + int gapSide = collision.gapSide(); + constexpr int kGapSideSelection = 0; + constexpr int kGapSideOppositeSelection = 2; + if (gapSide > kGapSideSelection && gapSide < kGapSideOppositeSelection) { + return; + } + if (collision.trs() == 0) { + return; + } + int trueGapSide = sgSelector.trueGap(collision, cfgCutFV0, cfgCutFT0A, cfgCutFT0C, cfgCutZDC); + gapSide = trueGapSide; + if (gapSide == cfgGapSideSelection) { + return; + } + registry.fill(HIST("hEventCount"), 1.5); + float cent = 100; + float lRandom = fRndm->Rndm(); + float vtxz = collision.posZ(); + registry.fill(HIST("hVtxZ"), vtxz); + registry.fill(HIST("hMult"), tracks.size()); + registry.fill(HIST("hCent"), cent); + fGFW->Clear(); + + // // track weights + float weff = 1, wacc = 1; + double nTracksCorrected = 0; + float independent = cent; + if (cfgUseNch) { + independent = static_cast(tracks.size()); + } - // LOG(info) << "filling mc particle px: " << mcParticle.px() << ", py: " << mcParticle.py() << ", pz: " << mcParticle.pz() << std::endl; + for (const auto& track : tracks) { + registry.fill(HIST("hChi2prTPCcls"), track.tpcChi2NCl()); + if (!trackSelected(track)) + continue; + auto momentum = std::array{track.px(), track.py(), track.pz()}; + double pt = RecoDecay::pt(momentum); + double phi = RecoDecay::phi(momentum); + double eta = RecoDecay::eta(momentum); + bool withinPtPOI = (cfgCutPtPOIMin < pt) && (pt < cfgCutPtPOIMax); // within POI pT range + bool withinPtRef = (cfgCutPtRefMin < pt) && (pt < cfgCutPtRefMax); // within RF pT range + if (cfgOutputNUAWeights) { + if (cfgOutputNUAWeightsRefPt) { + if (withinPtRef) { + fWeights->fill(phi, eta, vtxz, pt, cent, 0); + } + } else { + fWeights->fill(phi, eta, vtxz, pt, cent, 0); + } + } + if (!setCurrentParticleWeights(weff, wacc, phi, eta, pt, vtxz)) { + continue; + } + registry.fill(HIST("hPt"), track.pt()); + if (withinPtRef) { + registry.fill(HIST("hPhi"), phi); + registry.fill(HIST("hPhiWeighted"), phi, wacc); + registry.fill(HIST("hEta"), eta); + registry.fill(HIST("hPtRef"), pt); + registry.fill(HIST("hDCAz"), track.dcaZ(), track.pt()); + registry.fill(HIST("hDCAxy"), track.dcaXY(), track.pt()); + nTracksCorrected += weff; + } + if (withinPtRef) { + fGFW->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 1); + } + if (withinPtPOI) { + fGFW->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 2); + } + if (withinPtPOI && withinPtRef) { + fGFW->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 4); + } + } + registry.fill(HIST("hTrackCorrection2d"), tracks.size(), nTracksCorrected); - // output information from mcparticles - registry.fill(HIST("hPxMc"), mcParticle.px()); - registry.fill(HIST("hPyMc"), mcParticle.py()); - registry.fill(HIST("hPzMc"), mcParticle.pz()); - registry.fill(HIST("hweightMc"), mcParticle.weight()); + // Filling Flow Container + for (uint l_ind = 0; l_ind < corrconfigs.size(); l_ind++) { + fillFC(corrconfigs.at(l_ind), independent, lRandom); + } + } + PROCESS_SWITCH(FlowCumulantsUpc, process, "process", true); - if (!mcParticle.isPhysicalPrimary()) { - // LOG(info) << "mcParticle.isPhysicalPrimary() = " << mcParticle.isPhysicalPrimary() << std::endl; - continue; - } + //----------------------------------------------------------------------------------------------------------------------- + void processSim(aod::UDMcCollision const& mcCollision, aod::UDMcParticles const& mcParticles) + { + registry.fill(HIST("eventCounterMC"), 0.5); + + registry.fill(HIST("hEventCount"), 1.5); + float cent = 100; + float vtxz = mcCollision.posZ(); + registry.fill(HIST("hVtxZMC"), vtxz); + registry.fill(HIST("hMultMC"), mcParticles.size()); + registry.fill(HIST("hCentMC"), cent); + + auto massPion = o2::constants::physics::MassPionCharged; + registry.fill(HIST("numberOfTracksMC"), mcParticles.size()); + // LOGF(info, "New event! mcParticles.size() = %d", mcParticles.size()); + + float lRandomMc = fRndmMc->Rndm(); + fGFWMC->Clear(); + + // // track weights + float weff = 1, wacc = 1; + double nTracksCorrected = 0; + float independent = cent; + if (cfgUseNch) { + independent = static_cast(mcParticles.size()); + } - std::array momentum = {mcParticle.px(), mcParticle.py(), mcParticle.pz()}; - double energy = std::sqrt(momentum[0] * momentum[0] + momentum[1] * momentum[1] + momentum[2] * momentum[2] + massPion * massPion); - ROOT::Math::LorentzVector> protoMC(momentum[0], momentum[1], momentum[2], energy); - constexpr double kEtaCut = 0.8; - constexpr double kPtCut = 0.1; - if (!(std::fabs(protoMC.Eta()) < kEtaCut && protoMC.Pt() > kPtCut)) { - // LOG(info) << "protoMC.Eta() = " << protoMC.Eta() << ", protoMC.Pt() = " << protoMC.Pt() << std::endl; - continue; - } - // auto momentum = std::array{mcParticle.px(), mcParticle.py(), mcParticle.pz()}; - double pt = RecoDecay::pt(momentum); - double phi = RecoDecay::phi(momentum); - double eta = RecoDecay::eta(momentum); - bool withinPtPOI = (cfgCutPtPOIMin < pt) && (pt < cfgCutPtPOIMax); // within POI pT range - bool withinPtRef = (cfgCutPtRefMin < pt) && (pt < cfgCutPtRefMax); // within RF pT range - if (cfgOutputNUAWeights) { - if (cfgOutputNUAWeightsRefPt) { - if (withinPtRef) { - fWeightsMc->fill(phi, eta, vtxz, pt, cent, 0); - } - } else { + for (const auto& mcParticle : mcParticles) { + if (!mcParticle.isPhysicalPrimary()) + continue; + std::array momentum = {mcParticle.px(), mcParticle.py(), mcParticle.pz()}; + double energy = std::sqrt(momentum[0] * momentum[0] + momentum[1] * momentum[1] + momentum[2] * momentum[2] + massPion * massPion); + ROOT::Math::LorentzVector> protoMC(momentum[0], momentum[1], momentum[2], energy); + constexpr double kEtaCut = 0.8; + constexpr double kPtCut = 0.1; + if (!(std::fabs(protoMC.Eta()) < kEtaCut && protoMC.Pt() > kPtCut)) { + continue; + } + // auto momentum = std::array{mcParticle.px(), mcParticle.py(), mcParticle.pz()}; + double pt = RecoDecay::pt(momentum); + double phi = RecoDecay::phi(momentum); + double eta = RecoDecay::eta(momentum); + bool withinPtPOI = (cfgCutPtPOIMin < pt) && (pt < cfgCutPtPOIMax); // within POI pT range + bool withinPtRef = (cfgCutPtRefMin < pt) && (pt < cfgCutPtRefMax); // within RF pT range + if (cfgOutputNUAWeights) { + if (cfgOutputNUAWeightsRefPt) { + if (withinPtRef) { fWeightsMc->fill(phi, eta, vtxz, pt, cent, 0); } + } else { + fWeightsMc->fill(phi, eta, vtxz, pt, cent, 0); } - if (!setCurrentParticleWeights(weff, wacc, phi, eta, pt, vtxz)) { - continue; - } - - if (withinPtRef) { - registry.fill(HIST("hPhiMC"), phi); - registry.fill(HIST("hPhiWeightedMC"), phi, wacc); - registry.fill(HIST("hEtaMC"), eta); - registry.fill(HIST("hPtRefMC"), pt); - // registry.fill(HIST("hDCAzMC"), track.dcaZ(), track.pt()); - // registry.fill(HIST("hDCAxyMC"), track.dcaXY(), track.pt()); - nTracksCorrected += weff; - } - if (withinPtRef) { - fGFWMC->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 1); - } - if (withinPtPOI) { - fGFWMC->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 2); - } - if (withinPtPOI && withinPtRef) { - fGFWMC->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 4); - } - // LOG(info) << "successfully filled" << std::endl; } - registry.fill(HIST("hTrackCorrection2dMC"), mcParticles.size(), nTracksCorrected); - // Filling Flow Container - for (uint l_ind = 0; l_ind < corrconfigsmc.size(); l_ind++) { - // LOG(info) << "filling flow container for MC" << std::endl; - fillFCMC(corrconfigsmc.at(l_ind), independent, lRandomMc); + if (!setCurrentParticleWeights(weff, wacc, phi, eta, pt, vtxz)) { + continue; + } + if (withinPtRef) { + registry.fill(HIST("hPhiMC"), phi); + registry.fill(HIST("hPhiWeightedMC"), phi, wacc); + registry.fill(HIST("hEtaMC"), eta); + registry.fill(HIST("hPtRefMC"), pt); + // registry.fill(HIST("hDCAzMC"), track.dcaZ(), track.pt()); + // registry.fill(HIST("hDCAxyMC"), track.dcaXY(), track.pt()); + nTracksCorrected += weff; + } + if (withinPtRef) { + fGFWMC->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 1); + } + if (withinPtPOI) { + fGFWMC->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 2); } + if (withinPtPOI && withinPtRef) { + fGFWMC->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 4); + } + } + registry.fill(HIST("hTrackCorrection2dMC"), mcParticles.size(), nTracksCorrected); + + // Filling Flow Container + for (uint l_ind = 0; l_ind < corrconfigs.size(); l_ind++) { + fillFCMC(corrconfigs.at(l_ind), independent, lRandomMc); } + PROCESS_SWITCH(FlowCumulantsUpc, processSim, "processSim", false); } - PROCESS_SWITCH(FlowCumulantsUpc, processSim, "processSim", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) From 94b30406d7e80f3c05bc0046aa9d40753ac11f9c Mon Sep 17 00:00:00 2001 From: miedema-11 <1260971129@qq.com> Date: Sun, 25 Jan 2026 23:50:30 +0800 Subject: [PATCH 6/8] remove tpc cluster --- PWGUD/Tasks/flowCorrelationsUpc.cxx | 42 ++++++++++++++++++++++------- PWGUD/Tasks/flowCumulantsUpc.cxx | 4 --- 2 files changed, 33 insertions(+), 13 deletions(-) diff --git a/PWGUD/Tasks/flowCorrelationsUpc.cxx b/PWGUD/Tasks/flowCorrelationsUpc.cxx index 81235ce9c6f..190d20e2b1a 100644 --- a/PWGUD/Tasks/flowCorrelationsUpc.cxx +++ b/PWGUD/Tasks/flowCorrelationsUpc.cxx @@ -103,7 +103,7 @@ struct FlowCorrelationsUpc { O2_DEFINE_CONFIGURABLE(cfgSampleSize, double, 10, "Sample size 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(cfgCutMerging, float, 0.0, "Merging cut on track merge") + O2_DEFINE_CONFIGURABLE(cfgCutMerging, float, 0.02, "Merging cut on track merge") O2_DEFINE_CONFIGURABLE(cfgRadiusLow, float, 0.8, "Low radius for merging cut") O2_DEFINE_CONFIGURABLE(cfgRadiusHigh, float, 2.5, "High radius for merging cut") @@ -160,7 +160,7 @@ struct FlowCorrelationsUpc { registry.add("Trig_hist", "", {HistType::kTHnSparseF, {{axisSample, axisVertex, axisPtTrigger}}}); - registry.add("eventcount", "bin", {HistType::kTH1F, {{3, 0, 3, "bin"}}}); // histogram to see how many events are in the same and mixed event + registry.add("eventcount", "bin", {HistType::kTH1F, {{4, 0, 4, "bin"}}}); // histogram to see how many events are in the same and mixed event std::vector corrAxis = {{axisSample, "Sample"}, {axisVertex, "z-vtx (cm)"}, @@ -224,19 +224,17 @@ struct FlowCorrelationsUpc { if (!(std::fabs(track.dcaXY()) < dcaLimit)) { return false; } - constexpr int kMinTPCClusters = 70; constexpr int kMinITSClusters = 5; constexpr int kMaxTPCChi2NCl = 4; - if (track.tpcNClsFindableMinusCrossedRows() <= kMinTPCClusters) { - return false; - } if (track.itsClusterSizes() <= kMinITSClusters) { return false; } if (track.tpcChi2NCl() >= kMaxTPCChi2NCl) { return false; } + if (track.pt() < cfgPtCutMin || track.pt() > cfgPtCutMax) + return false; return true; } @@ -264,6 +262,9 @@ struct FlowCorrelationsUpc { // loop over all tracks for (auto const& track1 : tracks1) { + if (!trackSelected(track1)) + continue; + if (system == SameEvent) { registry.fill(HIST("Trig_hist"), fSampleIndex, posZ, track1.pt()); } @@ -285,9 +286,6 @@ struct FlowCorrelationsUpc { float deltaPhi = RecoDecay::constrainAngle(phi1 - phi2, -PIHalf); float deltaEta = RecoDecay::eta(momentum1) - RecoDecay::eta(momentum2); - if (pt2 < cfgPtCutMin || pt2 > cfgPtCutMax) - continue; - if (std::abs(deltaEta) < cfgCutMerging) { double dPhiStarHigh = getDPhiStar(track1, track2, cfgRadiusHigh, runnum, phi1, phi2); @@ -365,6 +363,32 @@ struct FlowCorrelationsUpc { SameKindPair pairs{binningOnVtxAndMult, cfgMinMixEventNum, -1, collisions, tracksTuple, &cache}; // -1 is the number of the bin to skip for (auto const& [collision1, tracks1, collision2, tracks2] : pairs) { + if (tracks1.size() < cfgMinMult || tracks1.size() > cfgMaxMult || tracks2.size() < cfgMinMult || tracks2.size() > cfgMaxMult) { + continue; + } + if (collision1.trs() == 0 || collision2.trs() == 0) { + continue; + } + + const int minGapSide = 0; + const int maxGapSide = 2; + if (collision1.gapSide() > minGapSide && collision1.gapSide() < maxGapSide) { + continue; + } + if (collision2.gapSide() > minGapSide && collision2.gapSide() < maxGapSide) { + continue; + } + + int trueGapSide = sgSelector.trueGap(collision1, cfgCutFV0, cfgCutFT0A, cfgCutFT0C, cfgCutZDC); + int gapSide = trueGapSide; + if (gapSide == cfgGapSideSelection) { + continue; + } + trueGapSide = sgSelector.trueGap(collision2, cfgCutFV0, cfgCutFT0A, cfgCutFT0C, cfgCutZDC); + gapSide = trueGapSide; + if (gapSide == cfgGapSideSelection) { + continue; + } registry.fill(HIST("eventcount"), MixedEvent); // fill the mixed event in the 3 bin fillCorrelations(tracks1, tracks2, collision1.posZ(), MixedEvent, collision1.runNumber()); // fill the ME histogram and Sparse } diff --git a/PWGUD/Tasks/flowCumulantsUpc.cxx b/PWGUD/Tasks/flowCumulantsUpc.cxx index 12332abbd6a..8ede8d2d136 100644 --- a/PWGUD/Tasks/flowCumulantsUpc.cxx +++ b/PWGUD/Tasks/flowCumulantsUpc.cxx @@ -813,13 +813,9 @@ struct FlowCumulantsUpc { if (!(std::fabs(track.dcaXY()) < dcaLimit)) { return false; } - constexpr int kMinTPCClusters = 70; constexpr int kMinITSClusters = 5; constexpr int kMaxTPCChi2NCl = 4; - if (track.tpcNClsFindableMinusCrossedRows() <= kMinTPCClusters) { - return false; - } if (track.itsClusterSizes() <= kMinITSClusters) { return false; } From de81db088f17db70b2212f6797dc46056b30f75a Mon Sep 17 00:00:00 2001 From: dyx-11 <1260971129@qq.com> Date: Mon, 26 Jan 2026 15:28:20 +0800 Subject: [PATCH 7/8] Remove unused variable pt2 from flowCorrelationsUpc.cxx --- PWGUD/Tasks/flowCorrelationsUpc.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/PWGUD/Tasks/flowCorrelationsUpc.cxx b/PWGUD/Tasks/flowCorrelationsUpc.cxx index 190d20e2b1a..4b0cc17e487 100644 --- a/PWGUD/Tasks/flowCorrelationsUpc.cxx +++ b/PWGUD/Tasks/flowCorrelationsUpc.cxx @@ -280,7 +280,6 @@ struct FlowCorrelationsUpc { auto momentum1 = std::array{track1.px(), track1.py(), track1.pz()}; auto momentum2 = std::array{track2.px(), track2.py(), track2.pz()}; - double pt2 = RecoDecay::pt(momentum2); double phi1 = RecoDecay::phi(momentum1); double phi2 = RecoDecay::phi(momentum2); float deltaPhi = RecoDecay::constrainAngle(phi1 - phi2, -PIHalf); From b5f637d0e3c257595c7573cd534b9df263f2f69e Mon Sep 17 00:00:00 2001 From: miedema-11 <1260971129@qq.com> Date: Sat, 14 Feb 2026 14:58:37 +0800 Subject: [PATCH 8/8] fix gapside selection --- PWGUD/Tasks/flowCorrelationsUpc.cxx | 141 ++++++--- PWGUD/Tasks/flowCumulantsUpc.cxx | 442 ++++++++++++++++------------ 2 files changed, 356 insertions(+), 227 deletions(-) diff --git a/PWGUD/Tasks/flowCorrelationsUpc.cxx b/PWGUD/Tasks/flowCorrelationsUpc.cxx index 4b0cc17e487..3de478cddbf 100644 --- a/PWGUD/Tasks/flowCorrelationsUpc.cxx +++ b/PWGUD/Tasks/flowCorrelationsUpc.cxx @@ -11,13 +11,14 @@ /// \file flowCorrelationsUpc.cxx /// \brief Provides a sparse with usefull two particle correlation info -/// \author Mingrui Zhao (mingrui.zhao@cern.ch, mingrui.zhao@mail.labz0.org) +/// \author Yongxi Du (yongxi.du@cern.ch), Mingrui Zhao (mingrui.zhao@cern.ch, mingrui.zhao@mail.labz0.org) /// copied from Thor Jensen (thor.kjaersgaard.jensen@cern.ch) and Debojit Sarkar (debojit.sarkar@cern.ch) #include "PWGCF/Core/CorrelationContainer.h" #include "PWGCF/Core/PairCuts.h" #include "PWGCF/DataModel/CorrelationsDerived.h" #include "PWGUD/Core/SGSelector.h" +#include "PWGUD/DataModel/SGTables.h" #include "PWGUD/DataModel/UDTables.h" #include "Common/Core/RecoDecay.h" @@ -94,7 +95,8 @@ struct CalcNchUpc { struct FlowCorrelationsUpc { O2_DEFINE_CONFIGURABLE(cfgZVtxCut, float, 10.0f, "Accepted z-vertex range") - O2_DEFINE_CONFIGURABLE(cfgPtCutMin, float, 0.2f, "minimum accepted track pT") + O2_DEFINE_CONFIGURABLE(cfgIfVertex, bool, false, "choose vertex or not") + O2_DEFINE_CONFIGURABLE(cfgPtCutMin, float, 0.1f, "minimum accepted track pT") O2_DEFINE_CONFIGURABLE(cfgPtCutMax, float, 10.0f, "maximum accepted track pT") O2_DEFINE_CONFIGURABLE(cfgEtaCut, float, 0.8f, "Eta cut") O2_DEFINE_CONFIGURABLE(cfgMinMixEventNum, int, 5, "Minimum number of events to mix") @@ -106,6 +108,17 @@ struct FlowCorrelationsUpc { O2_DEFINE_CONFIGURABLE(cfgCutMerging, float, 0.02, "Merging cut on track merge") O2_DEFINE_CONFIGURABLE(cfgRadiusLow, float, 0.8, "Low radius for merging cut") O2_DEFINE_CONFIGURABLE(cfgRadiusHigh, float, 2.5, "High radius for merging cut") + O2_DEFINE_CONFIGURABLE(cfgIsGoodItsLayers, bool, false, "whether choose itslayers") + O2_DEFINE_CONFIGURABLE(cfgGapSideA, bool, true, "choose gapside A") + O2_DEFINE_CONFIGURABLE(cfgGapSideC, bool, false, "choose gapside C") + O2_DEFINE_CONFIGURABLE(cfgDcaxy, bool, true, "choose dcaxy") + O2_DEFINE_CONFIGURABLE(cfgDcaz, bool, false, "choose dcaz") + O2_DEFINE_CONFIGURABLE(cfgDcazCut, float, 10.0, "dcaz cut") + O2_DEFINE_CONFIGURABLE(cfgItsClusterSize, unsigned int, 5, "ITS cluster size") + O2_DEFINE_CONFIGURABLE(cfgMaxTPCChi2NCl, int, 4, "tpcchi2") + O2_DEFINE_CONFIGURABLE(cfgEvSelOccupancy, bool, true, "Occupancy cut") + O2_DEFINE_CONFIGURABLE(cfgCutOccupancyHigh, int, 1000, "High cut on TPC occupancy") + O2_DEFINE_CONFIGURABLE(cfgCutOccupancyLow, int, 0, "Low cut on TPC occupancy") ConfigurableAxis axisVertex{"axisVertex", {10, -10, 10}, "vertex axis for histograms"}; ConfigurableAxis axisEta{"axisEta", {40, -1., 1.}, "eta axis for histograms"}; @@ -130,7 +143,6 @@ struct FlowCorrelationsUpc { Configurable cfgCutFT0A{"cfgCutFT0A", 150., "FT0A threshold"}; Configurable cfgCutFT0C{"cfgCutFT0C", 50., "FT0C threshold"}; Configurable cfgCutZDC{"cfgCutZDC", 10., "ZDC threshold"}; - Configurable cfgGapSideSelection{"cfgGapSideSelection", 2, "gap selection"}; // make the filters and cuts. // Filter collisionFilter = (nabs(aod::collision::posZ) < cfgZVtxCut) && (aod::flowcorrupc::multiplicity) > cfgMinMult && (aod::flowcorrupc::multiplicity) < cfgMaxMult && (aod::evsel::sel8) == true; @@ -212,29 +224,29 @@ struct FlowCorrelationsUpc { template bool trackSelected(TTrack track) { + registry.fill(HIST("hTrackCount"), 0.5); // UPC selection if (!track.isPVContributor()) { return false; } - constexpr float kDcazCut = 2.0; - if (!(std::fabs(track.dcaZ()) < kDcazCut)) { + registry.fill(HIST("hTrackCount"), 1.5); + if (cfgDcaz && !(std::fabs(track.dcaZ()) < cfgDcazCut)) { return false; } + registry.fill(HIST("hTrackCount"), 2.5); double dcaLimit = 0.0105 + 0.035 / std::pow(track.pt(), 1.1); - if (!(std::fabs(track.dcaXY()) < dcaLimit)) { + if (cfgDcaxy && !(std::fabs(track.dcaXY()) < dcaLimit)) { return false; } - constexpr int kMinITSClusters = 5; - constexpr int kMaxTPCChi2NCl = 4; - - if (track.itsClusterSizes() <= kMinITSClusters) { + registry.fill(HIST("hTrackCount"), 3.5); + if (track.itsClusterSizes() <= cfgItsClusterSize) { return false; } - if (track.tpcChi2NCl() >= kMaxTPCChi2NCl) { + registry.fill(HIST("hTrackCount"), 4.5); + if (track.tpcChi2NCl() >= cfgMaxTPCChi2NCl) { return false; } - if (track.pt() < cfgPtCutMin || track.pt() > cfgPtCutMax) - return false; + registry.fill(HIST("hTrackCount"), 5.5); return true; } @@ -324,23 +336,47 @@ struct FlowCorrelationsUpc { if (tracks.size() < cfgMinMult || tracks.size() > cfgMaxMult) { return; } - if (collision.trs() == 0) { + if (cfgIsGoodItsLayers && collision.trs() == 0) { return; } int gapSide = collision.gapSide(); - const int minGapSide = 0; - const int maxGapSide = 2; - if (gapSide > minGapSide && gapSide < maxGapSide) { + if (gapSide == 0) { + if (!cfgGapSideA) { + return; + } + } + if (gapSide == 1) { + if (!cfgGapSideC) { + return; + } + } + if (gapSide != 0 && gapSide != 1) { return; } - int trueGapSide = sgSelector.trueGap(collision, cfgCutFV0, cfgCutFT0A, cfgCutFT0C, cfgCutZDC); gapSide = trueGapSide; - if (gapSide == cfgGapSideSelection) { + if (gapSide == 0) { + if (!cfgGapSideA) { + return; + } + } + if (gapSide == 1) { + if (!cfgGapSideC) { + return; + } + } + if (gapSide != 0 && gapSide != 1) { + return; + } + float vtxz = collision.posZ(); + if (cfgIfVertex && abs(vtxz) > cfgZVtxCut) { + return; + } + int occupancy = collision.occupancyInTime(); + if (cfgEvSelOccupancy && (occupancy < cfgCutOccupancyLow || occupancy > cfgCutOccupancyHigh)) { return; } - int runIndex = collision.runNumber(); registry.fill(HIST("eventcount"), SameEvent); // because its same event i put it in the 1 bin @@ -365,29 +401,66 @@ struct FlowCorrelationsUpc { if (tracks1.size() < cfgMinMult || tracks1.size() > cfgMaxMult || tracks2.size() < cfgMinMult || tracks2.size() > cfgMaxMult) { continue; } - if (collision1.trs() == 0 || collision2.trs() == 0) { + if (cfgIsGoodItsLayers && (collision1.trs() == 0 || collision2.trs() == 0)) { continue; } - const int minGapSide = 0; - const int maxGapSide = 2; - if (collision1.gapSide() > minGapSide && collision1.gapSide() < maxGapSide) { - continue; + int gapSide1 = collision1.gapSide(); + if (gapSide1 == 0) { + if (!cfgGapSideA) { + continue; + } } - if (collision2.gapSide() > minGapSide && collision2.gapSide() < maxGapSide) { - continue; + if (gapSide1 == 1) { + if (!cfgGapSideC) { + continue; + } } - - int trueGapSide = sgSelector.trueGap(collision1, cfgCutFV0, cfgCutFT0A, cfgCutFT0C, cfgCutZDC); - int gapSide = trueGapSide; - if (gapSide == cfgGapSideSelection) { + int gapSide2 = collision2.gapSide(); + if (gapSide2 == 0) { + if (!cfgGapSideA) { + continue; + } + } + if (gapSide2 == 1) { + if (!cfgGapSideC) { + continue; + } + } + int trueGapSide1 = sgSelector.trueGap(collision1, cfgCutFV0, cfgCutFT0A, cfgCutFT0C, cfgCutZDC); + int trueGapSide2 = sgSelector.trueGap(collision2, cfgCutFV0, cfgCutFT0A, cfgCutFT0C, cfgCutZDC); + if (trueGapSide1 != trueGapSide2) { continue; } - trueGapSide = sgSelector.trueGap(collision2, cfgCutFV0, cfgCutFT0A, cfgCutFT0C, cfgCutZDC); - gapSide = trueGapSide; - if (gapSide == cfgGapSideSelection) { + if (trueGapSide1 == 0) { + if (!cfgGapSideA) { + continue; + } + } + if (trueGapSide2 == 1) { + if (!cfgGapSideC) { + continue; + } + } + if ((gapSide1 != 0 && gapSide1 != 1) || (gapSide2 != 0 && gapSide2 != 1)) { continue; } + float vtxz = collision1.posZ(); + if (cfgIfVertex && abs(vtxz) > cfgZVtxCut) { + return; + } + int occupancy = collision1.occupancyInTime(); + if (cfgEvSelOccupancy && (occupancy < cfgCutOccupancyLow || occupancy > cfgCutOccupancyHigh)) { + return; + } + vtxz = collision2.posZ(); + if (cfgIfVertex && abs(vtxz) > cfgZVtxCut) { + return; + } + occupancy = collision2.occupancyInTime(); + if (cfgEvSelOccupancy && (occupancy < cfgCutOccupancyLow || occupancy > cfgCutOccupancyHigh)) { + return; + } registry.fill(HIST("eventcount"), MixedEvent); // fill the mixed event in the 3 bin fillCorrelations(tracks1, tracks2, collision1.posZ(), MixedEvent, collision1.runNumber()); // fill the ME histogram and Sparse } diff --git a/PWGUD/Tasks/flowCumulantsUpc.cxx b/PWGUD/Tasks/flowCumulantsUpc.cxx index 8ede8d2d136..98006bbde6a 100644 --- a/PWGUD/Tasks/flowCumulantsUpc.cxx +++ b/PWGUD/Tasks/flowCumulantsUpc.cxx @@ -10,7 +10,7 @@ // or submit itself to any jurisdiction. /// \file flowCumulantsUpc.cxx -/// \author Mingrui Zhao (mingrui.zhao@mail.labz0.org, mingrui.zhao@cern.ch) +/// \author Yongxi Du (yongxi.du@cern.ch), Mingrui Zhao (mingrui.zhao@mail.labz0.org, mingrui.zhao@cern.ch) /// \since Mar/2025 /// \brief jira: , task to measure flow observables with cumulant method @@ -21,6 +21,7 @@ #include "GFWWeights.h" #include "PWGUD/Core/SGSelector.h" +#include "PWGUD/DataModel/SGTables.h" #include "PWGUD/DataModel/UDTables.h" #include "Common/CCDB/ctpRateFetcher.h" @@ -62,14 +63,15 @@ using namespace o2::framework::expressions; struct FlowCumulantsUpc { O2_DEFINE_CONFIGURABLE(cfgCutVertex, float, 10.0f, "Accepted z-vertex range") + O2_DEFINE_CONFIGURABLE(cfgIfVertex, bool, false, "choose vertex or not") O2_DEFINE_CONFIGURABLE(cfgCentEstimator, int, 0, "0:FT0C; 1:FT0CVariant1; 2:FT0M; 3:FT0A") O2_DEFINE_CONFIGURABLE(cfgCentFT0CMin, float, 0.0f, "Minimum centrality (FT0C) to cut events in filter") O2_DEFINE_CONFIGURABLE(cfgCentFT0CMax, float, 100.0f, "Maximum centrality (FT0C) to cut events in filter") - O2_DEFINE_CONFIGURABLE(cfgCutPtPOIMin, float, 0.2f, "Minimal pT for poi tracks") + O2_DEFINE_CONFIGURABLE(cfgCutPtPOIMin, float, 0.1f, "Minimal pT for poi tracks") O2_DEFINE_CONFIGURABLE(cfgCutPtPOIMax, float, 10.0f, "Maximal pT for poi tracks") - O2_DEFINE_CONFIGURABLE(cfgCutPtRefMin, float, 0.2f, "Minimal pT for ref tracks") + O2_DEFINE_CONFIGURABLE(cfgCutPtRefMin, float, 0.1f, "Minimal pT for ref tracks") O2_DEFINE_CONFIGURABLE(cfgCutPtRefMax, float, 3.0f, "Maximal pT for ref tracks") - O2_DEFINE_CONFIGURABLE(cfgCutPtMin, float, 0.2f, "Minimal pT for all tracks") + O2_DEFINE_CONFIGURABLE(cfgCutPtMin, float, 0.1f, "Minimal pT for all tracks") O2_DEFINE_CONFIGURABLE(cfgCutPtMax, float, 10.0f, "Maximal pT for all tracks") O2_DEFINE_CONFIGURABLE(cfgCutEta, float, 0.8f, "Eta range for tracks") O2_DEFINE_CONFIGURABLE(cfgCutChi2prTPCcls, float, 2.5f, "max chi2 per TPC clusters") @@ -77,22 +79,22 @@ struct FlowCumulantsUpc { O2_DEFINE_CONFIGURABLE(cfgCutITSclu, float, 5.0f, "minimum ITS clusters") O2_DEFINE_CONFIGURABLE(cfgCutDCAz, float, 2.0f, "max DCA to vertex z") O2_DEFINE_CONFIGURABLE(cfgCutDCAxyppPass3Enabled, bool, false, "switch of ppPass3 DCAxy pt dependent cut") - O2_DEFINE_CONFIGURABLE(cfgCutDCAzPtDepEnabled, bool, false, "switch of DCAz pt dependent cut") - O2_DEFINE_CONFIGURABLE(cfgTrkSelSwitch, bool, false, "switch for self-defined track selection") - O2_DEFINE_CONFIGURABLE(cfgUseAdditionalEventCut, bool, false, "Use additional event cut on mult correlations") - O2_DEFINE_CONFIGURABLE(cfgUseTentativeEventCounter, bool, false, "After sel8(), count events regardless of real event selection") - O2_DEFINE_CONFIGURABLE(cfgEvSelkNoSameBunchPileup, bool, false, "rejects collisions which are associated with the same found-by-T0 bunch crossing") - O2_DEFINE_CONFIGURABLE(cfgEvSelkIsGoodZvtxFT0vsPV, bool, false, "removes collisions with large differences between z of PV by tracks and z of PV from FT0 A-C time difference, use this cut at low multiplicities with caution") - O2_DEFINE_CONFIGURABLE(cfgEvSelkNoCollInTimeRangeStandard, bool, false, "no collisions in specified time range") - O2_DEFINE_CONFIGURABLE(cfgEvSelkIsGoodITSLayersAll, bool, true, "cut time intervals with dead ITS staves") - O2_DEFINE_CONFIGURABLE(cfgEvSelkNoCollInRofStandard, bool, false, "no other collisions in this Readout Frame with per-collision multiplicity above threshold") - O2_DEFINE_CONFIGURABLE(cfgEvSelkNoHighMultCollInPrevRof, bool, false, "veto an event if FT0C amplitude in previous ITS ROF is above threshold") - O2_DEFINE_CONFIGURABLE(cfgEvSelMultCorrelation, bool, true, "Multiplicity correlation cut") - O2_DEFINE_CONFIGURABLE(cfgEvSelV0AT0ACut, bool, true, "V0A T0A 5 sigma cut") - O2_DEFINE_CONFIGURABLE(cfgGetInteractionRate, bool, false, "Get interaction rate from CCDB") - O2_DEFINE_CONFIGURABLE(cfgUseInteractionRateCut, bool, false, "Use events with low interaction rate") - O2_DEFINE_CONFIGURABLE(cfgCutMaxIR, float, 50.0f, "maximum interaction rate (kHz)") - O2_DEFINE_CONFIGURABLE(cfgCutMinIR, float, 0.0f, "minimum interaction rate (kHz)") + // O2_DEFINE_CONFIGURABLE(cfgCutDCAzPtDepEnabled, bool, false, "switch of DCAz pt dependent cut") + // O2_DEFINE_CONFIGURABLE(cfgTrkSelSwitch, bool, false, "switch for self-defined track selection") + // O2_DEFINE_CONFIGURABLE(cfgUseAdditionalEventCut, bool, false, "Use additional event cut on mult correlations") + // O2_DEFINE_CONFIGURABLE(cfgUseTentativeEventCounter, bool, false, "After sel8(), count events regardless of real event selection") + // O2_DEFINE_CONFIGURABLE(cfgEvSelkNoSameBunchPileup, bool, false, "rejects collisions which are associated with the same found-by-T0 bunch crossing") + // O2_DEFINE_CONFIGURABLE(cfgEvSelkIsGoodZvtxFT0vsPV, bool, false, "removes collisions with large differences between z of PV by tracks and z of PV from FT0 A-C time difference, use this cut at low multiplicities with caution") + // O2_DEFINE_CONFIGURABLE(cfgEvSelkNoCollInTimeRangeStandard, bool, false, "no collisions in specified time range") + // O2_DEFINE_CONFIGURABLE(cfgEvSelkIsGoodITSLayersAll, bool, true, "cut time intervals with dead ITS staves") + // O2_DEFINE_CONFIGURABLE(cfgEvSelkNoCollInRofStandard, bool, false, "no other collisions in this Readout Frame with per-collision multiplicity above threshold") + // O2_DEFINE_CONFIGURABLE(cfgEvSelkNoHighMultCollInPrevRof, bool, false, "veto an event if FT0C amplitude in previous ITS ROF is above threshold") + // O2_DEFINE_CONFIGURABLE(cfgEvSelMultCorrelation, bool, true, "Multiplicity correlation cut") + // O2_DEFINE_CONFIGURABLE(cfgEvSelV0AT0ACut, bool, true, "V0A T0A 5 sigma cut") + // O2_DEFINE_CONFIGURABLE(cfgGetInteractionRate, bool, false, "Get interaction rate from CCDB") + // O2_DEFINE_CONFIGURABLE(cfgUseInteractionRateCut, bool, false, "Use events with low interaction rate") + // O2_DEFINE_CONFIGURABLE(cfgCutMaxIR, float, 50.0f, "maximum interaction rate (kHz)") + // O2_DEFINE_CONFIGURABLE(cfgCutMinIR, float, 0.0f, "minimum interaction rate (kHz)") O2_DEFINE_CONFIGURABLE(cfgUseNch, bool, false, "Use Nch for flow observables") O2_DEFINE_CONFIGURABLE(cfgNbootstrap, int, 30, "Number of subsamples") O2_DEFINE_CONFIGURABLE(cfgOutputNUAWeights, bool, false, "Fill and output NUA weights") @@ -102,9 +104,17 @@ struct FlowCumulantsUpc { O2_DEFINE_CONFIGURABLE(cfgAcceptanceList, std::string, "", "CCDB path to acceptance lsit object") O2_DEFINE_CONFIGURABLE(cfgAcceptanceListEnabled, bool, false, "switch of acceptance list") O2_DEFINE_CONFIGURABLE(cfgEvSelOccupancy, bool, true, "Occupancy cut") - O2_DEFINE_CONFIGURABLE(cfgCutOccupancyHigh, int, 500, "High cut on TPC occupancy") + O2_DEFINE_CONFIGURABLE(cfgCutOccupancyHigh, int, 1000, "High cut on TPC occupancy") O2_DEFINE_CONFIGURABLE(cfgCutOccupancyLow, int, 0, "Low cut on TPC occupancy") O2_DEFINE_CONFIGURABLE(cfgUseSmallMemory, bool, false, "Use small memory mode") + O2_DEFINE_CONFIGURABLE(cfgIsGoodItsLayers, bool, false, "whether choose itslayers") + O2_DEFINE_CONFIGURABLE(cfgGapSideA, bool, true, "choose gapside A") + O2_DEFINE_CONFIGURABLE(cfgGapSideC, bool, false, "choose gapside C") + O2_DEFINE_CONFIGURABLE(cfgDcaxy, bool, true, "choose dcaxy") + O2_DEFINE_CONFIGURABLE(cfgDcaz, bool, false, "choose dcaz") + O2_DEFINE_CONFIGURABLE(cfgDcazCut, float, 10.0, "dcaz cut") + O2_DEFINE_CONFIGURABLE(cfgItsClusterSize, unsigned int, 5, "ITS cluster size") + O2_DEFINE_CONFIGURABLE(cfgMaxTPCChi2NCl, int, 4, "tpcchi2") Configurable> cfgUserDefineGFWCorr{"cfgUserDefineGFWCorr", std::vector{"refN02 {2} refP02 {-2}", "refN12 {2} refP12 {-2}"}, "User defined GFW CorrelatorConfig"}; Configurable> cfgUserDefineGFWName{"cfgUserDefineGFWName", std::vector{"Ch02Gap22", "Ch12Gap22"}, "User defined GFW Name"}; Configurable> cfgRunRemoveList{"cfgRunRemoveList", std::vector{-1}, "excluded run numbers"}; @@ -122,7 +132,6 @@ struct FlowCumulantsUpc { Configurable cfgCutFT0A{"cfgCutFT0A", 150., "FT0A threshold"}; Configurable cfgCutFT0C{"cfgCutFT0C", 50., "FT0C threshold"}; Configurable cfgCutZDC{"cfgCutZDC", 10., "ZDC threshold"}; - Configurable cfgGapSideSelection{"cfgGapSideSelection", 2, "gap selection"}; // Corrections TH1D* mEfficiency = nullptr; @@ -200,10 +209,17 @@ struct FlowCumulantsUpc { // Event QA registry.add("hEventCount", "Number of Event;; Count", {HistType::kTH1D, {{5, 0, 5}}}); registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(1, "Filtered event"); - registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(2, "after sel8"); - registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(3, "after supicious Runs removal"); - registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(4, "after additional event cut"); - registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(5, "after correction loads"); + registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(2, "after gapside selection"); + registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(3, "after its selection"); + registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(4, "after pt selection"); + registry.get(HIST("hEventCount"))->GetXaxis()->SetBinLabel(5, "after occupancy"); + registry.add("hTrackCount", "Number of tracks;; Count", {HistType::kTH1D, {{5, 0, 5}}}); + registry.get(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(1, "after event selection"); + registry.get(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(2, "PVContributor"); + registry.get(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(3, "dcaz"); + registry.get(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(4, "dcaxy"); + registry.get(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(5, "its clusters"); + registry.get(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(6, "tpc chi2"); registry.add("hEventCountSpecific", "Number of Event;; Count", {HistType::kTH1D, {{10, 0, 10}}}); registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(1, "after sel8"); registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(2, "kNoSameBunchPileup"); @@ -215,19 +231,19 @@ struct FlowCumulantsUpc { registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(8, "occupancy"); registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(9, "MultCorrelation"); registry.get(HIST("hEventCountSpecific"))->GetXaxis()->SetBinLabel(10, "cfgEvSelV0AT0ACut"); - if (cfgUseTentativeEventCounter) { - registry.add("hEventCountTentative", "Number of Event;; Count", {HistType::kTH1D, {{10, 0, 10}}}); - registry.get(HIST("hEventCountTentative"))->GetXaxis()->SetBinLabel(1, "after sel8"); - registry.get(HIST("hEventCountTentative"))->GetXaxis()->SetBinLabel(2, "kNoSameBunchPileup"); - registry.get(HIST("hEventCountTentative"))->GetXaxis()->SetBinLabel(3, "kIsGoodZvtxFT0vsPV"); - registry.get(HIST("hEventCountTentative"))->GetXaxis()->SetBinLabel(4, "kNoCollInTimeRangeStandard"); - registry.get(HIST("hEventCountTentative"))->GetXaxis()->SetBinLabel(5, "kIsGoodITSLayersAll"); - registry.get(HIST("hEventCountTentative"))->GetXaxis()->SetBinLabel(6, "kNoCollInRofStandard"); - registry.get(HIST("hEventCountTentative"))->GetXaxis()->SetBinLabel(7, "kNoHighMultCollInPrevRof"); - registry.get(HIST("hEventCountTentative"))->GetXaxis()->SetBinLabel(8, "occupancy"); - registry.get(HIST("hEventCountTentative"))->GetXaxis()->SetBinLabel(9, "MultCorrelation"); - registry.get(HIST("hEventCountTentative"))->GetXaxis()->SetBinLabel(10, "cfgEvSelV0AT0ACut"); - } + // if (cfgUseTentativeEventCounter) { + // registry.add("hEventCountTentative", "Number of Event;; Count", {HistType::kTH1D, {{10, 0, 10}}}); + // registry.get(HIST("hEventCountTentative"))->GetXaxis()->SetBinLabel(1, "after sel8"); + // registry.get(HIST("hEventCountTentative"))->GetXaxis()->SetBinLabel(2, "kNoSameBunchPileup"); + // registry.get(HIST("hEventCountTentative"))->GetXaxis()->SetBinLabel(3, "kIsGoodZvtxFT0vsPV"); + // registry.get(HIST("hEventCountTentative"))->GetXaxis()->SetBinLabel(4, "kNoCollInTimeRangeStandard"); + // registry.get(HIST("hEventCountTentative"))->GetXaxis()->SetBinLabel(5, "kIsGoodITSLayersAll"); + // registry.get(HIST("hEventCountTentative"))->GetXaxis()->SetBinLabel(6, "kNoCollInRofStandard"); + // registry.get(HIST("hEventCountTentative"))->GetXaxis()->SetBinLabel(7, "kNoHighMultCollInPrevRof"); + // registry.get(HIST("hEventCountTentative"))->GetXaxis()->SetBinLabel(8, "occupancy"); + // registry.get(HIST("hEventCountTentative"))->GetXaxis()->SetBinLabel(9, "MultCorrelation"); + // registry.get(HIST("hEventCountTentative"))->GetXaxis()->SetBinLabel(10, "cfgEvSelV0AT0ACut"); + // } registry.add("hVtxZ", "Vexter Z distribution", {HistType::kTH1D, {axisVertex}}); registry.add("hMult", "Multiplicity distribution", {HistType::kTH1D, {{3000, 0.5, 3000.5}}}); std::string hCentTitle = "Centrality distribution, Estimator " + std::to_string(cfgCentEstimator); @@ -266,6 +282,7 @@ struct FlowCumulantsUpc { registry.add("hDCAz", "DCAz after cuts; DCAz (cm); Pt", {HistType::kTH2D, {{200, -0.5, 0.5}, {200, 0, 5}}}); registry.add("hDCAxy", "DCAxy after cuts; DCAxy (cm); Pt", {HistType::kTH2D, {{200, -0.5, 0.5}, {200, 0, 5}}}); registry.add("hTrackCorrection2d", "Correlation table for number of tracks table; uncorrected track; corrected track", {HistType::kTH2D, {axisNch, axisNch}}); + registry.add("hEtaNch2D", "Eta vs Nch; #eta; Nch", {HistType::kTH2D, {axisEta, axisNch}}); registry.add("hPhiMC", "#phi distribution", {HistType::kTH1D, {axisPhi}}); registry.add("hPhiWeightedMC", "corrected #phi distribution", {HistType::kTH1D, {axisPhi}}); @@ -281,6 +298,13 @@ struct FlowCumulantsUpc { registry.add("hDCAxyMC", "DCAxy after cuts; DCAxy (cm); Pt", {HistType::kTH2D, {{200, -0.5, 0.5}, {200, 0, 5}}}); registry.add("hTrackCorrection2dMC", "Correlation table for number of tracks table; uncorrected track; corrected track", {HistType::kTH2D, {axisNch, axisNch}}); + // // MC event QA histograms + // registry.add("eventCounterMC", "Number of MC Events;; Count", {HistType::kTH1D, {{5, 0, 5}}}); + // registry.add("hVtxZMC", "Vexter Z distribution (MC)", {HistType::kTH1D, {axisVertex}}); + // registry.add("hMultMC", "Multiplicity distribution (MC)", {HistType::kTH1D, {{3000, 0.5, 3000.5}}}); + // registry.add("hCentMC", "Centrality distribution (MC)", {HistType::kTH1D, {{90, 0, 90}}}); + // registry.add("numberOfTracksMC", "Number of MC tracks;; Count", {HistType::kTH1D, {{3000, 0.5, 3000.5}}}); + o2::framework::AxisSpec axis = axisPt; int nPtBins = axis.binEdges.size() - 1; double* ptBins = &(axis.binEdges)[0]; @@ -509,22 +533,22 @@ struct FlowCumulantsUpc { } 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); - - fT0AV0AMean = new TF1("fT0AV0AMean", "[0]+[1]*x", 0, 200000); - fT0AV0AMean->SetParameters(-1601.0581, 9.417652e-01); - fT0AV0ASigma = new TF1("fT0AV0ASigma", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x", 0, 200000); - fT0AV0ASigma->SetParameters(463.4144, 6.796509e-02, -9.097136e-07, 7.971088e-12, -2.600581e-17); - } + // 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); + + // fT0AV0AMean = new TF1("fT0AV0AMean", "[0]+[1]*x", 0, 200000); + // fT0AV0AMean->SetParameters(-1601.0581, 9.417652e-01); + // fT0AV0ASigma = new TF1("fT0AV0ASigma", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x", 0, 200000); + // fT0AV0ASigma->SetParameters(463.4144, 6.796509e-02, -9.097136e-07, 7.971088e-12, -2.600581e-17); + // } myTrackSel = getGlobalTrackSelectionRun3ITSMatch(TrackSelection::GlobalTrackRun3ITSMatching::Run3ITSibAny, TrackSelection::GlobalTrackRun3DCAxyCut::Default); myTrackSel.SetMinNClustersTPC(cfgCutTPCclu); @@ -677,151 +701,153 @@ struct FlowCumulantsUpc { return true; } - template - bool eventSelected(TCollision collision, const int multTrk, const float centrality) - { - registry.fill(HIST("hEventCountSpecific"), 0.5); - if (cfgEvSelkNoSameBunchPileup && !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 (cfgEvSelkNoSameBunchPileup) { - registry.fill(HIST("hEventCountSpecific"), 1.5); - } - if (cfgEvSelkIsGoodZvtxFT0vsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { - // removes collisions with large differences between z of PV by tracks and z of PV from FT0 A-C time difference - // use this cut at low multiplicities with caution - return 0; - } - if (cfgEvSelkIsGoodZvtxFT0vsPV) { - registry.fill(HIST("hEventCountSpecific"), 2.5); - } - if (cfgEvSelkNoCollInTimeRangeStandard && !collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) { - // no collisions in specified time range - return 0; - } - if (cfgEvSelkNoCollInTimeRangeStandard) { - registry.fill(HIST("hEventCountSpecific"), 3.5); - } - if (cfgEvSelkIsGoodITSLayersAll && !collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) { - // from Jan 9 2025 AOT meeting - // cut time intervals with dead ITS staves - return 0; - } - if (cfgEvSelkIsGoodITSLayersAll) { - registry.fill(HIST("hEventCountSpecific"), 4.5); - } - if (cfgEvSelkNoCollInRofStandard && !collision.selection_bit(o2::aod::evsel::kNoCollInRofStandard)) { - // no other collisions in this Readout Frame with per-collision multiplicity above threshold - return 0; - } - if (cfgEvSelkNoCollInRofStandard) { - registry.fill(HIST("hEventCountSpecific"), 5.5); - } - if (cfgEvSelkNoHighMultCollInPrevRof && !collision.selection_bit(o2::aod::evsel::kNoHighMultCollInPrevRof)) { - // veto an event if FT0C amplitude in previous ITS ROF is above threshold - return 0; - } - if (cfgEvSelkNoHighMultCollInPrevRof) { - registry.fill(HIST("hEventCountSpecific"), 6.5); - } - auto multNTracksPV = collision.multNTracksPV(); - auto occupancy = collision.trackOccupancyInTimeRange(); - if (cfgEvSelOccupancy && (occupancy < cfgCutOccupancyLow || occupancy > cfgCutOccupancyHigh)) { - return 0; - } - if (cfgEvSelOccupancy) { - registry.fill(HIST("hEventCountSpecific"), 7.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 (cfgEvSelMultCorrelation) { - registry.fill(HIST("hEventCountSpecific"), 8.5); - } - - // V0A T0A 5 sigma cut - constexpr int kSigmaCut = 5; - if (cfgEvSelV0AT0ACut && (std::fabs(collision.multFV0A() - fT0AV0AMean->Eval(collision.multFT0A())) > kSigmaCut * fT0AV0ASigma->Eval(collision.multFT0A()))) { - return 0; - } - if (cfgEvSelV0AT0ACut) { - registry.fill(HIST("hEventCountSpecific"), 9.5); - } - - return 1; - } - - template - void eventCounterQA(TCollision collision, const int multTrk, const float centrality) - { - registry.fill(HIST("hEventCountTentative"), 0.5); - // Regradless of the event selection, fill the event counter histograms - if (collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { - registry.fill(HIST("hEventCountTentative"), 1.5); - } - if (collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { - registry.fill(HIST("hEventCountTentative"), 2.5); - } - if (collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) { - registry.fill(HIST("hEventCountTentative"), 3.5); - } - if (collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) { - registry.fill(HIST("hEventCountTentative"), 4.5); - } - if (collision.selection_bit(o2::aod::evsel::kNoCollInRofStandard)) { - registry.fill(HIST("hEventCountTentative"), 5.5); - } - if (collision.selection_bit(o2::aod::evsel::kNoHighMultCollInPrevRof)) { - registry.fill(HIST("hEventCountTentative"), 6.5); - } - auto multNTracksPV = collision.multNTracksPV(); - auto occupancy = collision.trackOccupancyInTimeRange(); - if (!(occupancy < cfgCutOccupancyLow || occupancy > cfgCutOccupancyHigh)) { - registry.fill(HIST("hEventCountTentative"), 7.5); - } - if (!((multNTracksPV < fMultPVCutLow->Eval(centrality)) || (multNTracksPV > fMultPVCutHigh->Eval(centrality)) || (multTrk < fMultCutLow->Eval(centrality)) || (multTrk > fMultCutHigh->Eval(centrality)))) { - registry.fill(HIST("hEventCountTentative"), 8.5); - } - // constexpr int kSigmaCut = 5; - // if (!(std::fabs(collision.multFV0A() - fT0AV0AMean->Eval(collision.multFT0A())) > kSigmaCut * fT0AV0ASigma->Eval(collision.multFT0A()))) { - // registry.fill(HIST("hEventCountTentative"), 9.5); - // } - } + // template + // bool eventSelected(TCollision collision, const int multTrk, const float centrality) + // { + // registry.fill(HIST("hEventCountSpecific"), 0.5); + // if (cfgEvSelkNoSameBunchPileup && !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 (cfgEvSelkNoSameBunchPileup) { + // // registry.fill(HIST("hEventCountSpecific"), 1.5); + // // } + // if (cfgEvSelkIsGoodZvtxFT0vsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { + // // removes collisions with large differences between z of PV by tracks and z of PV from FT0 A-C time difference + // // use this cut at low multiplicities with caution + // return 0; + // } + // if (cfgEvSelkIsGoodZvtxFT0vsPV) { + // registry.fill(HIST("hEventCountSpecific"), 2.5); + // } + // if (cfgEvSelkNoCollInTimeRangeStandard && !collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) { + // // no collisions in specified time range + // return 0; + // } + // if (cfgEvSelkNoCollInTimeRangeStandard) { + // registry.fill(HIST("hEventCountSpecific"), 3.5); + // } + // if (cfgEvSelkIsGoodITSLayersAll && !collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) { + // // from Jan 9 2025 AOT meeting + // // cut time intervals with dead ITS staves + // return 0; + // } + // if (cfgEvSelkIsGoodITSLayersAll) { + // registry.fill(HIST("hEventCountSpecific"), 4.5); + // } + // if (cfgEvSelkNoCollInRofStandard && !collision.selection_bit(o2::aod::evsel::kNoCollInRofStandard)) { + // // no other collisions in this Readout Frame with per-collision multiplicity above threshold + // return 0; + // } + // if (cfgEvSelkNoCollInRofStandard) { + // registry.fill(HIST("hEventCountSpecific"), 5.5); + // } + // if (cfgEvSelkNoHighMultCollInPrevRof && !collision.selection_bit(o2::aod::evsel::kNoHighMultCollInPrevRof)) { + // // veto an event if FT0C amplitude in previous ITS ROF is above threshold + // return 0; + // } + // if (cfgEvSelkNoHighMultCollInPrevRof) { + // registry.fill(HIST("hEventCountSpecific"), 6.5); + // } + // auto multNTracksPV = collision.multNTracksPV(); + // auto occupancy = collision.occupancyInTime(); + // if (cfgEvSelOccupancy && (occupancy < cfgCutOccupancyLow || occupancy > cfgCutOccupancyHigh)) { + // return 0; + // } + // if (cfgEvSelOccupancy) { + // registry.fill(HIST("hEventCountSpecific"), 7.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 (cfgEvSelMultCorrelation) { + // registry.fill(HIST("hEventCountSpecific"), 8.5); + // } + + // // V0A T0A 5 sigma cut + // constexpr int kSigmaCut = 5; + // if (cfgEvSelV0AT0ACut && (std::fabs(collision.multFV0A() - fT0AV0AMean->Eval(collision.multFT0A())) > kSigmaCut * fT0AV0ASigma->Eval(collision.multFT0A()))) { + // return 0; + // } + // if (cfgEvSelV0AT0ACut) { + // registry.fill(HIST("hEventCountSpecific"), 9.5); + // } + + // return 1; + // } + + // // template + // // void eventCounterQA(TCollision collision, const int multTrk, const float centrality) + // // { + // // registry.fill(HIST("hEventCountTentative"), 0.5); + // // // Regradless of the event selection, fill the event counter histograms + // // if (collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { + // // registry.fill(HIST("hEventCountTentative"), 1.5); + // // } + // // if (collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { + // // registry.fill(HIST("hEventCountTentative"), 2.5); + // } + // if (collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) { + // registry.fill(HIST("hEventCountTentative"), 3.5); + // } + // if (collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) { + // registry.fill(HIST("hEventCountTentative"), 4.5); + // } + // if (collision.selection_bit(o2::aod::evsel::kNoCollInRofStandard)) { + // registry.fill(HIST("hEventCountTentative"), 5.5); + // } + // if (collision.selection_bit(o2::aod::evsel::kNoHighMultCollInPrevRof)) { + // registry.fill(HIST("hEventCountTentative"), 6.5); + // } + // auto multNTracksPV = collision.multNTracksPV(); + // auto occupancy = collision.trackOccupancyInTimeRange(); + // if (!(occupancy < cfgCutOccupancyLow || occupancy > cfgCutOccupancyHigh)) { + // registry.fill(HIST("hEventCountTentative"), 7.5); + // } + // if (!((multNTracksPV < fMultPVCutLow->Eval(centrality)) || (multNTracksPV > fMultPVCutHigh->Eval(centrality)) || (multTrk < fMultCutLow->Eval(centrality)) || (multTrk > fMultCutHigh->Eval(centrality)))) { + // registry.fill(HIST("hEventCountTentative"), 8.5); + // } + // // constexpr int kSigmaCut = 5; + // // if (!(std::fabs(collision.multFV0A() - fT0AV0AMean->Eval(collision.multFT0A())) > kSigmaCut * fT0AV0ASigma->Eval(collision.multFT0A()))) { + // // registry.fill(HIST("hEventCountTentative"), 9.5); + // // } + // } template bool trackSelected(TTrack track) { + registry.fill(HIST("hTrackCount"), 0.5); // UPC selection if (!track.isPVContributor()) { return false; } - constexpr float kDcazCut = 2.0; - if (!(std::fabs(track.dcaZ()) < kDcazCut)) { + registry.fill(HIST("hTrackCount"), 1.5); + if (cfgDcaz && !(std::fabs(track.dcaZ()) < cfgDcazCut)) { return false; } + registry.fill(HIST("hTrackCount"), 2.5); double dcaLimit = 0.0105 + 0.035 / std::pow(track.pt(), 1.1); - if (!(std::fabs(track.dcaXY()) < dcaLimit)) { + if (cfgDcaxy && !(std::fabs(track.dcaXY()) < dcaLimit)) { return false; } - constexpr int kMinITSClusters = 5; - constexpr int kMaxTPCChi2NCl = 4; - - if (track.itsClusterSizes() <= kMinITSClusters) { + registry.fill(HIST("hTrackCount"), 3.5); + if (track.itsClusterSizes() <= cfgItsClusterSize) { return false; } - if (track.tpcChi2NCl() >= kMaxTPCChi2NCl) { + registry.fill(HIST("hTrackCount"), 4.5); + if (track.tpcChi2NCl() >= cfgMaxTPCChi2NCl) { return false; } + registry.fill(HIST("hTrackCount"), 5.5); return true; } @@ -879,20 +905,39 @@ struct FlowCumulantsUpc { // return; // } int gapSide = collision.gapSide(); - constexpr int kGapSideSelection = 0; - constexpr int kGapSideOppositeSelection = 2; - if (gapSide > kGapSideSelection && gapSide < kGapSideOppositeSelection) { - return; + if (gapSide == 0) { + if (!cfgGapSideA) { + return; + } + } + if (gapSide == 1) { + if (!cfgGapSideC) { + return; + } } - if (collision.trs() == 0) { + if (gapSide != 0 && gapSide != 1) { return; } int trueGapSide = sgSelector.trueGap(collision, cfgCutFV0, cfgCutFT0A, cfgCutFT0C, cfgCutZDC); gapSide = trueGapSide; - if (gapSide == cfgGapSideSelection) { + if (gapSide == 0) { + if (!cfgGapSideA) { + return; + } + } + if (gapSide == 1) { + if (!cfgGapSideC) { + return; + } + } + if (gapSide != 0 && gapSide != 1) { return; } registry.fill(HIST("hEventCount"), 1.5); + if (cfgIsGoodItsLayers && collision.trs() == 0) { + return; + } + registry.fill(HIST("hEventCount"), 2.5); float cent = 100; float lRandom = fRndm->Rndm(); float vtxz = collision.posZ(); @@ -900,7 +945,15 @@ struct FlowCumulantsUpc { registry.fill(HIST("hMult"), tracks.size()); registry.fill(HIST("hCent"), cent); fGFW->Clear(); - + if (cfgIfVertex && abs(vtxz) > cfgCutVertex) { + return; + } + registry.fill(HIST("hEventCount"), 3.5); + int occupancy = collision.occupancyInTime(); + if (cfgEvSelOccupancy && (occupancy < cfgCutOccupancyLow || occupancy > cfgCutOccupancyHigh)) { + return; + } + registry.fill(HIST("hEventCount"), 4.5); // // track weights float weff = 1, wacc = 1; double nTracksCorrected = 0; @@ -911,8 +964,9 @@ struct FlowCumulantsUpc { for (const auto& track : tracks) { registry.fill(HIST("hChi2prTPCcls"), track.tpcChi2NCl()); - if (!trackSelected(track)) + if (!trackSelected(track)) { continue; + } auto momentum = std::array{track.px(), track.py(), track.pz()}; double pt = RecoDecay::pt(momentum); double phi = RecoDecay::phi(momentum); @@ -950,6 +1004,7 @@ struct FlowCumulantsUpc { if (withinPtPOI && withinPtRef) { fGFW->Fill(eta, fPtAxis->FindBin(pt) - 1, phi, wacc * weff, 4); } + registry.fill(HIST("hEtaNch2D"), eta, tracks.size()); } registry.fill(HIST("hTrackCorrection2d"), tracks.size(), nTracksCorrected); @@ -1041,12 +1096,13 @@ struct FlowCumulantsUpc { for (uint l_ind = 0; l_ind < corrconfigs.size(); l_ind++) { fillFCMC(corrconfigs.at(l_ind), independent, lRandomMc); } - PROCESS_SWITCH(FlowCumulantsUpc, processSim, "processSim", false); } + PROCESS_SWITCH(FlowCumulantsUpc, processSim, "processSim", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask(cfgc)}; + adaptAnalysisTask(cfgc), + }; }