[PWGCF] Adding a new task for short range correlations #7
Conversation
Adding a new task for short range correlations as a function of eta
github-actions
bot
changed the title
There was a problem hiding this comment.
Pull request overview
This PR introduces a new analysis task for di-hadron correlations with eta-dependence in O-O and Pb-Pb collisions. The implementation is a comprehensive 1203-line file that handles same-event and mixed-event correlations, including MC truth processing.
Changes:
- Added new task
etaDihadron.cxximplementing di-hadron correlation analysis with eta-dependent binning - Includes support for various centrality estimators, event selection criteria, and efficiency corrections
- Provides MC truth processing capabilities for both full simulation and on-the-fly generation
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| for (double rad(cfgRadiusLow); rad < cfgRadiusHigh; rad += 0.01) { | ||
| double dPhiStar = getDPhiStar(track1, track2, rad, magneticField); | ||
| if (std::abs(dPhiStar) < kLimit) { | ||
| bIsBelow = true; | ||
| break; | ||
| } | ||
| } |
There was a problem hiding this comment.
Inefficient loop: The code iterates from cfgRadiusLow to cfgRadiusHigh with a step of 0.01, potentially executing hundreds of iterations per track pair when deltaEta is small. This nested loop (tracks × tracks × radius iterations) could be a significant performance bottleneck. Consider using a larger step size or a more efficient algorithm to check for merging, or at minimum document why such fine granularity is necessary.
| LOGF(fatal, "Could not load efficiency histogram for trigger particles from %s", cfgCentralityWeight.value.c_str()); | ||
| } | ||
| LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgCentralityWeight.value.c_str(), (void*)mCentralityWeight); |
There was a problem hiding this comment.
The getCentralityWeight function error message says "Could not load efficiency histogram for trigger particles" but the function is loading a centrality weight histogram, not an efficiency histogram. The error message should be corrected to "Could not load centrality weight histogram".
| LOGF(fatal, "Could not load efficiency histogram for trigger particles from %s", cfgCentralityWeight.value.c_str()); | |
| } | |
| LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgCentralityWeight.value.c_str(), (void*)mCentralityWeight); | |
| LOGF(fatal, "Could not load centrality weight histogram from %s", cfgCentralityWeight.value.c_str()); | |
| } | |
| LOGF(info, "Loaded centrality weight histogram from %s (%p)", cfgCentralityWeight.value.c_str(), (void*)mCentralityWeight); |
| registry.fill(HIST("Solo_tracks_trigger"), track1.pt()); | ||
| continue; // Skip the track1 if it is alone in pt bin | ||
| } | ||
| registry.fill(HIST("Trig_hist"), fSampleIndex, posZ, track1.pt(), eventWeight * triggerWeight); |
There was a problem hiding this comment.
Same issue as with line 538 - the histogram "Trig_hist" is being filled with track1.pt() but based on the axis definition (axisEtaTrigger), it should be filled with track1.eta(). This inconsistency will result in incorrect data in the histograms.
| registry.fill(HIST("Trig_hist"), fSampleIndex, posZ, track1.pt(), eventWeight * triggerWeight); | |
| registry.fill(HIST("Trig_hist"), fSampleIndex, posZ, track1.eta(), eventWeight * triggerWeight); |
| registry.fill(HIST("Solo_tracks_trigger"), track1.pt()); | ||
| continue; // Skip the track1 if it is alone in pt bin | ||
| } | ||
| registry.fill(HIST("Trig_hist"), fSampleIndex, posZ, track1.pt(), eventWeight * triggerWeight); | ||
|
|
||
| for (auto const& track2 : tracks2) { | ||
|
|
||
| if (!trackSelected(track2)) | ||
| continue; | ||
| if (mEfficiency) { | ||
| associatedWeight = efficiencyAssociatedCache[track2.filteredIndex()]; | ||
| } | ||
| if (!cfgUsePtOrder && track1.globalIndex() == track2.globalIndex()) | ||
| continue; // For pt-differential correlations, skip if the trigger and associate are the same track | ||
| if (cfgUsePtOrder && track1.pt() <= track2.pt()) | ||
| continue; // Without pt-differential correlations, skip if the trigger pt is less than the associate pt | ||
| if (!cfgSingleSoloPtTrack) { // avoid skipping the second track if we only want one | ||
| if (std::find(tracks2SkipIndices.begin(), tracks2SkipIndices.end(), track2.globalIndex()) != tracks2SkipIndices.end()) { | ||
| registry.fill(HIST("Solo_tracks_assoc"), track2.pt()); |
There was a problem hiding this comment.
The histograms "Solo_tracks_trigger" and "Solo_tracks_assoc" are also filling with track.pt() but have axisEtaTrigger/axisEtaAssoc axes (lines 286, 288). This is another instance of the pt/eta confusion throughout the code. These should be consistent.
| HistogramRegistry registry{"registry"}; | ||
|
|
||
| // define global variables | ||
| TRandom3* gRandom = new TRandom3(); |
There was a problem hiding this comment.
Memory leak: TRandom3 object is allocated with 'new' but never deleted. This should either be a stack object (TRandom3 gRandom;) or properly managed with a smart pointer, or deleted in a destructor. The current implementation will leak memory when the task ends.
| TRandom3* gRandom = new TRandom3(); | |
| TRandom3 gRandom; |
| // Get the magnetic field | ||
| static o2::parameters::GRPMagField* grpo = nullptr; | ||
| if (grpo == nullptr) { | ||
| grpo = ccdb->getForTimeStamp<o2::parameters::GRPMagField>("/GLO/Config/GRPMagField", timestamp); | ||
| if (grpo == nullptr) { | ||
| LOGF(fatal, "GRP object not found for timestamp %llu", timestamp); | ||
| return 0; | ||
| } | ||
| LOGF(info, "Retrieved GRP for timestamp %llu with magnetic field of %d kG", timestamp, grpo->getNominalL3Field()); | ||
| } |
There was a problem hiding this comment.
The static variable 'grpo' is initialized to nullptr and then only set once, but the check "if (grpo == nullptr)" means it will only fetch from CCDB once per process lifetime. This could be problematic if the magnetic field changes during a run or across different runs in the same job. Consider whether this caching behavior is intended, or if the GRP should be fetched per run or time period.
| // Get the magnetic field | |
| static o2::parameters::GRPMagField* grpo = nullptr; | |
| if (grpo == nullptr) { | |
| grpo = ccdb->getForTimeStamp<o2::parameters::GRPMagField>("/GLO/Config/GRPMagField", timestamp); | |
| if (grpo == nullptr) { | |
| LOGF(fatal, "GRP object not found for timestamp %llu", timestamp); | |
| return 0; | |
| } | |
| LOGF(info, "Retrieved GRP for timestamp %llu with magnetic field of %d kG", timestamp, grpo->getNominalL3Field()); | |
| } | |
| // Get the magnetic field for the given timestamp | |
| o2::parameters::GRPMagField* grpo = | |
| ccdb->getForTimeStamp<o2::parameters::GRPMagField>("/GLO/Config/GRPMagField", timestamp); | |
| if (grpo == nullptr) { | |
| LOGF(fatal, "GRP object not found for timestamp %llu", timestamp); | |
| return 0; | |
| } | |
| LOGF(info, "Retrieved GRP for timestamp %llu with magnetic field of %d kG", timestamp, | |
| grpo->getNominalL3Field()); |
| if (std::abs(track.eta()) > cfgCutEta) { | ||
| return false; | ||
| } | ||
| if (std::abs(track.pt()) < cfgCutPtMin || std::abs(track.pt()) > cfgCutPtMax) { |
There was a problem hiding this comment.
The condition check uses std::abs on track.pt() which returns a float. Since pT is already positive by definition (it's the magnitude of the transverse momentum), using std::abs is unnecessary and potentially confusing. The check should simply be "if (track.pt() < cfgCutPtMin || track.pt() > cfgCutPtMax)".
| if (std::abs(track.pt()) < cfgCutPtMin || std::abs(track.pt()) > cfgCutPtMax) { | |
| if (track.pt() < cfgCutPtMin || track.pt() > cfgCutPtMax) { |
| if (dPhiStar > constants::math::PI) | ||
| dPhiStar = constants::math::TwoPI - dPhiStar; | ||
| if (dPhiStar < -constants::math::PI) | ||
| dPhiStar = -constants::math::TwoPI - dPhiStar; |
There was a problem hiding this comment.
The logic here is incorrect for wrapping angles. When dPhiStar < -PI, you should add TwoPI, not subtract it. The current calculation "-constants::math::TwoPI - dPhiStar" will make a negative value even more negative. It should be "dPhiStar = constants::math::TwoPI + dPhiStar" or equivalently "dPhiStar += constants::math::TwoPI".
| dPhiStar = -constants::math::TwoPI - dPhiStar; | |
| dPhiStar = constants::math::TwoPI + dPhiStar; |
| if (!getEfficiencyCorrection(triggerWeight, track1.eta(), track1.pt(), posZ)) | ||
| continue; | ||
|
|
||
| registry.fill(HIST("Nch_final_pt"), track1.pt()); |
There was a problem hiding this comment.
The histogram "Nch_final_pt" has axis label "axisEtaTrigger" (based on line 285) but is being filled with track1.pt(). This is inconsistent - either the histogram name and axis should represent pT, or track1.eta() should be filled instead.
| registry.fill(HIST("Nch_final_pt"), track1.pt()); | |
| registry.fill(HIST("Nch_final_pt"), track1.eta()); |
| cfgFuncParas.fMultPVT0CCutLow = new TF1("fMultPVT0CCutLow", cfgFuncParas.cfgMultCentLowCutFunction->c_str(), 0, 100); | ||
| cfgFuncParas.fMultPVT0CCutLow->SetParameters(&(cfgFuncParas.multPVT0CCutPars[0])); | ||
| cfgFuncParas.fMultPVT0CCutHigh = new TF1("fMultPVT0CCutHigh", cfgFuncParas.cfgMultCentHighCutFunction->c_str(), 0, 100); | ||
| cfgFuncParas.fMultPVT0CCutHigh->SetParameters(&(cfgFuncParas.multPVT0CCutPars[0])); | ||
|
|
||
| cfgFuncParas.fMultT0CCutLow = new TF1("fMultT0CCutLow", cfgFuncParas.cfgMultCentLowCutFunction->c_str(), 0, 100); | ||
| cfgFuncParas.fMultT0CCutLow->SetParameters(&(cfgFuncParas.multT0CCutPars[0])); | ||
| cfgFuncParas.fMultT0CCutHigh = new TF1("fMultT0CCutHigh", cfgFuncParas.cfgMultCentHighCutFunction->c_str(), 0, 100); | ||
| cfgFuncParas.fMultT0CCutHigh->SetParameters(&(cfgFuncParas.multT0CCutPars[0])); | ||
|
|
||
| cfgFuncParas.fMultGlobalPVCutLow = new TF1("fMultGlobalPVCutLow", cfgFuncParas.cfgMultMultPVLowCutFunction->c_str(), 0, 4000); | ||
| cfgFuncParas.fMultGlobalPVCutLow->SetParameters(&(cfgFuncParas.multGlobalPVCutPars[0])); | ||
| cfgFuncParas.fMultGlobalPVCutHigh = new TF1("fMultGlobalPVCutHigh", cfgFuncParas.cfgMultMultPVHighCutFunction->c_str(), 0, 4000); | ||
| cfgFuncParas.fMultGlobalPVCutHigh->SetParameters(&(cfgFuncParas.multGlobalPVCutPars[0])); | ||
|
|
||
| cfgFuncParas.fMultMultV0ACutLow = new TF1("fMultMultV0ACutLow", cfgFuncParas.cfgMultMultV0ALowCutFunction->c_str(), 0, 4000); | ||
| cfgFuncParas.fMultMultV0ACutLow->SetParameters(&(cfgFuncParas.multMultV0ACutPars[0])); | ||
| cfgFuncParas.fMultMultV0ACutHigh = new TF1("fMultMultV0ACutHigh", cfgFuncParas.cfgMultMultV0AHighCutFunction->c_str(), 0, 4000); | ||
| cfgFuncParas.fMultMultV0ACutHigh->SetParameters(&(cfgFuncParas.multMultV0ACutPars[0])); | ||
|
|
||
| cfgFuncParas.fT0AV0AMean = new TF1("fT0AV0AMean", "[0]+[1]*x", 0, 200000); | ||
| cfgFuncParas.fT0AV0AMean->SetParameters(-1601.0581, 9.417652e-01); | ||
| cfgFuncParas.fT0AV0ASigma = new TF1("fT0AV0ASigma", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x", 0, 200000); | ||
| cfgFuncParas.fT0AV0ASigma->SetParameters(463.4144, 6.796509e-02, -9.097136e-07, 7.971088e-12, -2.600581e-17); |
There was a problem hiding this comment.
Potential memory leak: Multiple TF1 objects are allocated with 'new' (fMultPVT0CCutLow, fMultPVT0CCutHigh, fMultT0CCutLow, fMultT0CCutHigh, fMultGlobalPVCutLow, fMultGlobalPVCutHigh, fMultMultV0ACutLow, fMultMultV0ACutHigh, fT0AV0AMean, fT0AV0ASigma) but are never explicitly deleted. These should be managed with smart pointers or deleted in a destructor to prevent memory leaks.
1 participant
Adding a new task for short range correlations as a function of eta