CM Distortion Updates#4169
Conversation
📝 WalkthroughWalkthroughAdds graph-based region splitting and verbosity-controlled refinement to LaserClusterizer; loads CDB truth and adds total-distortion, outlier-skipping, and manual-interpolation modes in TpcCentralMembraneMatching; updates TpcLaminationFitting to internal interpolation, field-off handling, parameter scans, and a global R-matching routine. Changes
✨ Finishing touches
Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out. Comment |
![coderabbitai[bot]](https://avatars.githubusercontent.com/in/347564?s=60&v=4){kind=small}
| std::function<void(int)> dfs = [&](int u) | ||
| { | ||
| disc[u] = low[u] = ++time; | ||
| for(auto v : adj[u]) | ||
| { | ||
| if(disc[v] == -1) | ||
| { | ||
| parent[v] = u; | ||
| dfs(v); | ||
| low[u] = std::min(low[u], low[v]); | ||
| if(low[v] > disc[u]) | ||
| { | ||
| bridges.emplace_back(u,v); | ||
| } | ||
| } | ||
| else if(v != parent[u]) | ||
| { | ||
| low[u] = std::min(low[u], disc[v]); | ||
| } | ||
| } | ||
| }; | ||
|
|
||
| for(int i=0; i<N; i++) | ||
| { | ||
| if(disc[i] == -1) dfs(i); | ||
| } |
There was a problem hiding this comment.
Recursive DFS risks stack overflow on large proto-clusters.
The dfs lambda recurses to a depth equal to the connected-component size. TPC proto-clusters before splitting can easily reach hundreds to thousands of hits. With the default pthread stack (typically 2–8 MB), a deep recursion chain will overflow the stack and crash the process.
Consider converting to an iterative DFS using an explicit std::stack — the algorithm translates directly and eliminates the risk.
Sketch of iterative approach
- std::function<void(int)> dfs = [&](int u)
- {
- disc[u] = low[u] = ++time;
- for(auto v : adj[u])
- {
- if(disc[v] == -1)
- {
- parent[v] = u;
- dfs(v);
- low[u] = std::min(low[u], low[v]);
- if(low[v] > disc[u])
- {
- bridges.emplace_back(u,v);
- }
- }
- else if(v != parent[u])
- {
- low[u] = std::min(low[u], disc[v]);
- }
- }
- };
-
- for(int i=0; i<N; i++)
- {
- if(disc[i] == -1) dfs(i);
- }
+ // Iterative Tarjan bridge-finding using explicit stack
+ struct Frame { int u; size_t idx; };
+ std::stack<Frame> stk;
+ for(int i = 0; i < N; i++)
+ {
+ if(disc[i] != -1) continue;
+ stk.push({i, 0});
+ disc[i] = low[i] = ++timer;
+ while(!stk.empty())
+ {
+ auto &[u, idx] = stk.top();
+ if(idx < adj[u].size())
+ {
+ int v = adj[u][idx++];
+ if(disc[v] == -1)
+ {
+ parent[v] = u;
+ disc[v] = low[v] = ++timer;
+ stk.push({v, 0});
+ }
+ else if(v != parent[u])
+ {
+ low[u] = std::min(low[u], disc[v]);
+ }
+ }
+ else
+ {
+ stk.pop();
+ if(!stk.empty())
+ {
+ int p = stk.top().u;
+ low[p] = std::min(low[p], low[u]);
+ if(low[u] > disc[p])
+ {
+ bridges.emplace_back(p, u);
+ }
+ }
+ }
+ }
+ }| int TpcLaminationFitting::doGlobalRMatching(int side) | ||
| { | ||
|
|
||
| std::vector<double> distortedPhi; | ||
| TF1 *tmpLamFit = (TF1*)m_fLamination[0][side]->Clone(); | ||
|
|
||
| double meanB = 0.0; | ||
|
|
||
| if(m_fieldOff) | ||
| { | ||
| tmpLamFit->SetParameters(0.0, 0.0); | ||
| meanB = -999.99; | ||
| } | ||
| else | ||
| { | ||
| double meanA = 0.0; | ||
| double meanC = 0.0; | ||
| double meanOffset = 0.0; | ||
| int nGoodFits = 0; | ||
| for(int l = 0; l < 18; l++) | ||
| { | ||
| if(!m_laminationGoodFit[l][side]) | ||
| { | ||
| continue; | ||
| } | ||
| meanA += m_fLamination[l][side]->GetParameter(0); | ||
| meanB += m_fLamination[l][side]->GetParameter(1); | ||
| meanC += m_fLamination[l][side]->GetParameter(2); | ||
| meanOffset += m_laminationOffset[l][side]; | ||
| nGoodFits++; | ||
| } | ||
| if(nGoodFits == 0) | ||
| { | ||
| return Fun4AllReturnCodes::EVENT_OK; | ||
| } |
There was a problem hiding this comment.
tmpLamFit leaks, and early return leaves m_bestRMatch[side] null — later dereferenced in End().
Two issues in doGlobalRMatching:
-
Memory leak (line 826):
tmpLamFitis created viaClone()but never freed. Adddelete tmpLamFit;before eachreturn. -
Null-pointer crash (line 853-856): When
nGoodFits == 0, the function returnsEVENT_OKwithout allocatingm_bestRMatch[side]. Back inEnd()at line 1191,m_bestRMatch[s]->Write()is called unconditionally, resulting in a null dereference. This will crash whenever field is on and all lamination fits are bad (e.g., low statistics run).
Suggested fix
if(nGoodFits == 0)
{
+ delete tmpLamFit;
return Fun4AllReturnCodes::EVENT_OK;
}And at the end of the function, before the final return:
+ delete tmpLamFit;
return Fun4AllReturnCodes::EVENT_OK;
}And guard the write in End():
- m_bestRMatch[s]->Write();
+ if(m_bestRMatch[s]) { m_bestRMatch[s]->Write(); }
There was a problem hiding this comment.
Actionable comments posted: 4
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (2)
offline/packages/tpccalib/TpcLaminationFitting.cc (2)
988-993:⚠️ Potential issue | 🔴 CriticalDivision by zero when
scalers[0][2]is zero.
scalersis zero-initialized. If the GL1 database returns no row withindex == 0, or if the raw scaler count for that channel is genuinely zero, line 991 performsscalers[3][2] / scalers[0][2]→ IEEE 754 infinity (or NaN if both are zero). Thatm_ZDC_coincidencethen propagates intofitLaminationsas the ZDC rate used to seed fit parameters (Af[s]->Eval(m_ZDC_coincidence), etc.), producing garbage fits or downstream failures.Suggested guard
delete resultSet; + if(scalers[0][2] == 0) + { + std::cerr << "GL1 raw scaler for index 0 is zero for run " << m_runnumber + << ". Cannot compute ZDC coincidence rate, aborting." << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } + m_ZDC_coincidence = (1.0*scalers[3][2]/scalers[0][2])/(106e-9);
557-574:⚠️ Potential issue | 🟠 MajorLeaked
TGraphandTF1objects inside the 18 × 2 lamination loop.
gr,proj(line 557–558),fitSeed(line 571), andlocalFit(line 574) are heap-allocated every iteration but never freed. Over onefitLaminationscall that is 36 × 4 = 144 leaked objects.TGraphis not owned by ROOT's directory system, so these genuinely accumulate.Af[2]andBf[2](lines 479–480) also leak.Add cleanup at loop end and function exit
m_distanceToFit[l][s] = distToFunc / nBinsUsed; m_nBinsFit[l][s] = nBinsUsed; if(c>0) m_fitRMSE[l][s] = sqrt(wc / c); else m_fitRMSE[l][s] = -999; if (nBinsUsed < 10 || distToFunc / nBinsUsed > 1.0 || nBinsUsed_R_lt_45 < 5) { m_laminationGoodFit[l][s] = false; } else { m_laminationGoodFit[l][s] = true; } + + delete gr; + delete proj; + delete fitSeed; + delete localFit; } } + delete Af[0]; + delete Af[1]; + delete Bf[0]; + delete Bf[1]; + return Fun4AllReturnCodes::EVENT_OK; }
…oderabbit suggestions
There was a problem hiding this comment.
Actionable comments posted: 1
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (1)
offline/packages/tpccalib/TpcLaminationFitting.cc (1)
997-997:⚠️ Potential issue | 🟠 MajorDivision by zero when scaler data is missing or the clock scaler is zero.
scalersis zero-initialized. If the DB query returns no rows (e.g., run not yet in the scaler table), or ifscalers[0][2](raw clock scaler) happens to be zero, this line producesinforNaN. That value then seeds the lamination fit parameters viaAf[s]->Eval(m_ZDC_coincidence), silently corrupting all fits for the run.Suggested guard
+ if(scalers[0][2] == 0) + { + std::cerr << "TpcLaminationFitting::End - raw clock scaler is zero for run " + << m_runnumber << ", cannot compute ZDC rate. Aborting." << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } m_ZDC_coincidence = (1.0*scalers[3][2]/scalers[0][2])/(106e-9);
There was a problem hiding this comment.
Actionable comments posted: 1
Caution
Some comments are outside the diff and can’t be posted inline due to platform limitations.
⚠️ Outside diff range comments (2)
offline/packages/tpccalib/TpcLaminationFitting.cc (2)
557-636:⚠️ Potential issue | 🟠 MajorMemory leak: four ROOT objects allocated per lamination iteration are never freed.
gr,proj,fitSeed, andlocalFitarenew'd each iteration of the(s, l)double loop (36 iterations) but none are deleted. Over a long-running job this is 144 leaked ROOT objects per call tofitLaminations().Proposed fix — delete at end of inner loop body
m_distanceToFit[l][s] = distToFunc / nBinsUsed; m_nBinsFit[l][s] = nBinsUsed; ... + delete gr; + delete proj; + delete fitSeed; + delete localFit; } }Add these before the closing brace of the
for (int l …)loop (after the good-fit decision around line 708).
995-998:⚠️ Potential issue | 🔴 CriticalDivision by zero if
scalers[0][2]is zero.If the GL1 scaler index 0 has a raw count of zero (e.g., very short or empty run),
scalers[0][2]is 0 and line 997 divides by it, producinginfornanform_ZDC_coincidence. This value then seeds the fit parameters viaAf[s]->Eval(m_ZDC_coincidence)/Bf[s]->Eval(…), corrupting all lamination fits for the run.Proposed guard
delete resultSet; + if (scalers[0][2] == 0) + { + std::cerr << "TpcLaminationFitting: GL1 raw scaler[0] is zero for run " << m_runnumber << ", cannot compute ZDC rate" << std::endl; + return Fun4AllReturnCodes::ABORTRUN; + } m_ZDC_coincidence = (1.0*scalers[3][2]/scalers[0][2])/(106e-9);
| } | ||
|
|
||
| m_distanceToFit[l][s] = distToFunc / nBinsUsed; | ||
| m_nBinsFit[l][s] = nBinsUsed; | ||
| if (nBinsUsed < 10 || distToFunc / nBinsUsed > 1.0) | ||
| if(c>0) { m_fitRMSE[l][s] = sqrt(wc / c); | ||
| } else { m_fitRMSE[l][s] = -999; | ||
| } | ||
| if (nBinsUsed < 10 || distToFunc / nBinsUsed > 1.0 || nBinsUsed_R_lt_45 < 5) |
There was a problem hiding this comment.
Division by zero when nBinsUsed == 0.
If no bin along the fit passes the content threshold, nBinsUsed stays at 0 and line 696 divides by it. The same zero value feeds the comparison on line 701 (distToFunc / nBinsUsed > 1.0). This is undefined behavior and will produce nan/inf that propagates into m_distanceToFit.
Proposed fix
+ if (nBinsUsed == 0)
+ {
+ m_distanceToFit[l][s] = -999;
+ m_nBinsFit[l][s] = 0;
+ m_fitRMSE[l][s] = -999;
+ m_laminationGoodFit[l][s] = false;
+ continue;
+ }
m_distanceToFit[l][s] = distToFunc / nBinsUsed;
m_nBinsFit[l][s] = nBinsUsed;
Build & test reportReport for commit 54c8a9eda7c85ed830305b4448cb65efa73ff376:
Automatically generated by sPHENIX Jenkins continuous integration |
Build & test reportReport for commit 4dfbee8874259f8c65b43ada8fc5e678d96604a4:
Automatically generated by sPHENIX Jenkins continuous integration |
Build & test reportReport for commit 00ecb7f952ed535bec20cb1cdde6f8807bfd790d:
Automatically generated by sPHENIX Jenkins continuous integration |
|
The previous commits are clang-tidy fixes, merging this now so that we can get it into the new build and start working with it |
1 similar comment
|
The previous commits are clang-tidy fixes, merging this now so that we can get it into the new build and start working with it |
3 participants
Types of changes
What kind of change does this PR introduce? (Bug fix, feature, ...)
Updated laser clustering to look for dumbbell shaped proto-clusters and split them to keep only the portion with the highest ADC hit, which gets turned into a cluster. This has minimal impact in field on data, but is important to prevent cluster merging when the magnetic field is off.
Stripe matching now uses measured stripe pattern from field off data as truth (loaded from a CDBTTree) rather than ideal geometry. Also new methods to remove outliers from distortion calculation (flag false by default) as well as do a manual interpolation using weighted average where distance squared from the bin center is the weight (also has the flag as false by default).
Many updates to lamination fitting. Additional metric for goodness of fit (not fully implemented yet). New constant fit option added (with toggle) for use in field off data. Nominal lamination positions shifted by offset determined from field off data (currently hard-coded for each side).
New method for determining radial distortions to first order by assuming linear distortion as function of R. Scan parameter space of slope and y-intercept and use along with average lamination correction in phi to distort truth positions (again from measurements in field off data), calculate 1/chi2 by summing bin content in small window around each distorted truth position, highest 1/chi2 gives the best fit which is then put into distortion map for all phi values. Map of parameter space is stored in distortion file to keep track of scan.
TODOs (if applicable)
Links to other PRs in macros and calibration repositories (if applicable)
CM Distortion Updates
Motivation & Context
Reduce reconstruction artifacts (notably in field‑off laser runs) and improve TPC central membrane (CM) distortion estimation by moving to data‑driven stripe/lamination information, making clustering more robust against merged “dumbbell” proto‑clusters, and adding flexible, more-inspectable distortion interpolation and a first‑order radial distortion scan that produces QA and best‑fit parameters.
Key changes
Potential risk areas
Possible future improvements
Note: AI assistance can make mistakes—please inspect the code, run tests, and validate outputs before relying on this summary or merging.