Skip to content

[PWGUD] fix efficiency and cuts #16001

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
miedema-11 wants to merge 3 commits into
base: master
Choose a base branch
from
+13 −85
Draft

[PWGUD] fix efficiency and cuts #16001

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
967f799
fix efficiency and cuts
miedema-11 Apr 28, 2026
05f4434
delete unused parameters
miedema-11 Apr 29, 2026
cc8be7a
fix conflits
miedema-11 Apr 29, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
98 changes: 13 additions & 85 deletions PWGUD/Tasks/flowCorrelationsUpc.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -153,28 +153,20 @@ 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(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(cfgCutOccupancy, 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")
O2_DEFINE_CONFIGURABLE(cfgCutTPCCrossedRows, float, 70.0f, "minimum number of crossed TPC Rows")
O2_DEFINE_CONFIGURABLE(cfgCutTPCclu, float, 50.0f, "minimum number of found TPC clusters")
O2_DEFINE_CONFIGURABLE(cfgCutITSclu, float, 5.0f, "minimum number of ITS clusters")
O2_DEFINE_CONFIGURABLE(cfgCutTPCChi2NCl, int, 4, "max chi2 per TPC clusters")
O2_DEFINE_CONFIGURABLE(cfgGlobalTrack, bool, true, "require TPC+ITS track")
O2_DEFINE_CONFIGURABLE(cfgUseNchCorrected, bool, true, "use corrected Nch for X axis")
O2_DEFINE_CONFIGURABLE(cfgOutputNUAWeights, bool, false, "Fill and output NUA weights")
O2_DEFINE_CONFIGURABLE(cfgOutputNUAWeightsRefPt, bool, false, "NUA weights are filled in ref pt bins")
O2_DEFINE_CONFIGURABLE(cfgOutputNUAWeightsRunbyRun, bool, false, "NUA weights are filled run-by-run")
O2_DEFINE_CONFIGURABLE(cfgEfficiency, std::string, "", "CCDB path to efficiency object")
O2_DEFINE_CONFIGURABLE(cfgAcceptance, std::string, "", "CCDB path to acceptance object")
O2_DEFINE_CONFIGURABLE(cfgUseEventWeights, bool, false, "Use event weights for mixed event")
O2_DEFINE_CONFIGURABLE(cfgLocalEfficiency, bool, false, "Use local efficiency object")

ConfigurableAxis axisVertex{"axisVertex", {10, -10, 10}, "vertex axis for histograms"};
ConfigurableAxis axisEta{"axisEta", {40, -1., 1.}, "eta axis for histograms"};
Expand Down Expand Up @@ -203,25 +195,22 @@ struct FlowCorrelationsUpc {
ConfigurableAxis axisNch{"axisNch", {300, 0, 300}, "N_{ch}"};

// Corrections
TH1D* mEfficiency = nullptr;
GFWWeights* mAcceptance = nullptr;
TH3D* mEfficiency = nullptr;
bool correctionsLoaded = false;

// make the filters and cuts.
Filter trackFilter = (aod::udtrack::isPVContributor == true);
Filter collisionFilter = ((aod::udcollision::gapSide == (uint8_t)1 || aod::udcollision::gapSide == (uint8_t)0) && (cfgIfVertex == false || aod::collision::posZ < cfgZVtxCut) && (aod::udcollision::occupancyInTime > 0 && aod::udcollision::occupancyInTime < cfgCutOccupancyHigh) && (aod::flowcorrupc::truegapside == 1 || aod::flowcorrupc::truegapside == 0));
Filter collisionFilter = ((aod::udcollision::gapSide == (uint8_t)1 || aod::udcollision::gapSide == (uint8_t)0) && (cfgIfVertex == false || aod::collision::posZ < cfgZVtxCut) && (!cfgCutOccupancy || (aod::udcollision::occupancyInTime > 0 && aod::udcollision::occupancyInTime < cfgCutOccupancyHigh)) && (aod::flowcorrupc::truegapside == 1 || aod::flowcorrupc::truegapside == 0));

// Connect to ccdb
Service<ccdb::BasicCCDBManager> ccdb;
Configurable<std::string> ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};

OutputObj<GFWWeights> fWeights{GFWWeights("weights")};
OutputObj<GFWWeights> fWeightsMc{GFWWeights("weightsMC")};

TAxis* fPtAxis = nullptr;
int lastRunNumber = -1;
std::vector<int> runNumbers;
std::map<int, std::shared_ptr<TH3>> th3sPerRun; // map of TH3 histograms for all runs
std::vector<int> runNumbers; // map of TH3 histograms for all runs
std::vector<float> efficiencyCache;

using UdTracks = soa::Filtered<soa::Join<aod::UDTracks, aod::UDTracksExtra, aod::UDTracksPID>>;
Expand All @@ -248,9 +237,6 @@ struct FlowCorrelationsUpc {
registry.add("pT", "pT", {HistType::kTH1D, {axisPtTrigger}});
registry.add("Nch", "N_{ch}", {HistType::kTH1D, {axisMultiplicity}});
registry.add("zVtx", "zVtx", {HistType::kTH1D, {axisVertex}});
registry.add("Nch_vs_zVtx", "Nch vs zVtx", {HistType::kTH2D, {axisVertex, axisMultiplicity}});
registry.add("Nch_same", "Nch same event", {HistType::kTH1D, {axisMultiplicity}});
registry.add("Nch_mixed", "Nch mixed event", {HistType::kTH1D, {axisMultiplicity}});
registry.add("EtaCorrected", "Eta corrected", {HistType::kTH1D, {axisEta}});
registry.add("pTCorrected", "pT corrected", {HistType::kTH1D, {axisPtTrigger}});

Expand All @@ -273,13 +259,6 @@ struct FlowCorrelationsUpc {
double* ptBins = &(axis.binEdges)[0];
fPtAxis = new TAxis(nPtBins, ptBins);

if (cfgOutputNUAWeights) {
fWeights->setPtBins(nPtBins, ptBins);
fWeights->init(true, false);
fWeightsMc->setPtBins(nPtBins, ptBins);
fWeightsMc->init(true, false);
}

std::vector<AxisSpec> corrAxis = {{axisSample, "Sample"},
{axisVertex, "z-vtx (cm)"},
{axisIndependent, "Independent (N_{ch} corrected)"},
Expand Down Expand Up @@ -333,6 +312,9 @@ struct FlowCorrelationsUpc {
if (track.pt() < cfgPtCutMin || track.pt() > cfgPtCutMax) {
return false;
}
if (cfgGlobalTrack && !(track.hasITS() && track.hasTPC())) {
return false;
}
// registry.fill(HIST("hTrackCount"), 1.5);
if (cfgDcaz && !(std::fabs(track.dcaZ()) < cfgDcazCut)) {
return false;
Expand Down Expand Up @@ -361,28 +343,13 @@ struct FlowCorrelationsUpc {
return true;
}

void createOutputObjectsForRun(int runNumber)
{
const AxisSpec axisPhi{60, 0.0, constants::math::TwoPI, "#varphi"};
std::shared_ptr<TH3> histPhiEtaVtxz = registry.add<TH3>(Form("%d/hPhiEtaVtxz", runNumber), ";#varphi;#eta;v_{z}", {HistType::kTH3D, {axisPhi, {64, -1.8, 1.8}, {40, -10, 10}}});
th3sPerRun.insert(std::make_pair(runNumber, histPhiEtaVtxz));
}

void loadCorrections(uint64_t timestamp)
{
if (correctionsLoaded) {
return;
}
if (cfgAcceptance.value.empty() == false) {
mAcceptance = ccdb->getForTimeStamp<GFWWeights>(cfgAcceptance, timestamp);
if (mAcceptance) {
LOGF(info, "Loaded acceptance weights from %s (%p)", cfgAcceptance.value.c_str(), (void*)mAcceptance);
} else {
LOGF(warning, "Could not load acceptance weights from %s (%p)", cfgAcceptance.value.c_str(), (void*)mAcceptance);
}
}
if (cfgEfficiency.value.empty() == false) {
mEfficiency = ccdb->getForTimeStamp<TH1D>(cfgEfficiency, timestamp);
mEfficiency = ccdb->getForTimeStamp<TH3D>(cfgEfficiency, timestamp);
if (mEfficiency == nullptr) {
LOGF(fatal, "Could not load efficiency histogram for trigger particles from %s", cfgEfficiency.value.c_str());
}
Expand All @@ -408,7 +375,7 @@ struct FlowCorrelationsUpc {
return true;
}

bool setCurrentParticleWeights(float& weight_nue, float& weight_nua, float phi, float eta, float pt, float vtxz)
bool setCurrentParticleWeights(float& weight_nue, float& weight_nua, float pt)
{
float eff = 1.;
if (mEfficiency)
Expand All @@ -418,17 +385,13 @@ struct FlowCorrelationsUpc {
if (eff == 0)
return false;
weight_nue = 1. / eff;

if (mAcceptance)
weight_nua = mAcceptance->getNUA(phi, eta, vtxz);
else
weight_nua = 1;
weight_nua = 1.; // Set to 1 as NUA weight is not being used
return true;
}

// fill multiple histograms
template <typename TCollision, typename TTracks>
void fillYield(TCollision collision, TTracks tracks, int runNumber, float vtxz) // function to fill the yield and etaphi histograms.
void fillYield(TCollision collision, TTracks tracks, float vtxz) // function to fill the yield and etaphi histograms.
{
registry.fill(HIST("Nch"), tracks.size());
registry.fill(HIST("zVtx"), collision.posZ());
Expand All @@ -449,20 +412,6 @@ struct FlowCorrelationsUpc {
registry.fill(HIST("pT"), pt);
registry.fill(HIST("EtaCorrected"), eta, weff);
registry.fill(HIST("pTCorrected"), pt, weff);

if (cfgOutputNUAWeights) {
if (cfgOutputNUAWeightsRefPt) {
if (pt > cfgPtCutMin && pt < cfgPtCutMax) {
fWeights->fill(phi, eta, vtxz, pt, 0., 0);
if (cfgOutputNUAWeightsRunbyRun)
th3sPerRun[runNumber]->Fill(phi, eta, vtxz);
}
} else {
fWeights->fill(phi, eta, vtxz, pt, 0., 0);
if (cfgOutputNUAWeightsRunbyRun)
th3sPerRun[runNumber]->Fill(phi, eta, vtxz);
}
}
}
}

Expand Down Expand Up @@ -497,7 +446,7 @@ struct FlowCorrelationsUpc {

// 计算track1的权重
float weff1 = 1., wacc1 = 1.;
if (!setCurrentParticleWeights(weff1, wacc1, phi1, eta1, pt1, vtxz))
if (!setCurrentParticleWeights(weff1, wacc1, pt1))
continue;

if (system == SameEvent) {
Expand All @@ -520,24 +469,16 @@ struct FlowCorrelationsUpc {
double phi2 = RecoDecay::phi(momentum);
double eta2 = RecoDecay::eta(momentum);

// 计算track2的权重
float weff2 = 1., wacc2 = 1.;
if (mEfficiency) {
weff2 = efficiencyCache[track2.filteredIndex()];
} else {
getEfficiencyCorrection(weff2, eta2, pt2, vtxz);
}

if (mAcceptance) {
wacc2 = mAcceptance->getNUA(phi2, eta2, vtxz);
} else {
wacc2 = 1;
}

float deltaPhi = RecoDecay::constrainAngle(phi1 - phi2, -PIHalf);
float deltaEta = eta1 - eta2;

// 计算组合权重
float weight = eventWeight * weff1 * weff2 * wacc1 * wacc2;

// Merging cut
Expand Down Expand Up @@ -585,19 +526,6 @@ struct FlowCorrelationsUpc {

loadCorrections(runDuration.first);

if (cfgOutputNUAWeightsRunbyRun && currentRunNumber != lastRunNumber) {
lastRunNumber = currentRunNumber;
if (std::find(runNumbers.begin(), runNumbers.end(), currentRunNumber) == runNumbers.end()) {
createOutputObjectsForRun(currentRunNumber);
runNumbers.push_back(currentRunNumber);
}

if (th3sPerRun.find(currentRunNumber) == th3sPerRun.end()) {
LOGF(fatal, "RunNumber %d not found in th3sPerRun", currentRunNumber);
return;
}
}

registry.fill(HIST("eventcont"), 3.5);

//-----------independent---------------
Expand Down Expand Up @@ -628,7 +556,7 @@ struct FlowCorrelationsUpc {
independent = nTracksCorrected;
}

fillYield(collision, tracks, currentRunNumber, vtxz);
fillYield(collision, tracks, vtxz);

fillCorrelations<CorrelationContainer::kCFStepReconstructed>(
tracks, tracks, collision.posZ(), SameEvent,
Expand Down
Loading