Skip to content

Commit 2a950fd

Browse files
SangbaekSangbaek
committed
Avoid using ternary operator inside groot
1 parent 63db340 commit 2a950fd

1 file changed

Lines changed: 163 additions & 115 deletions

File tree

‎src/main/java/org/jlab/clas/timeline/fitter/ALERTFitter.groovy‎

Lines changed: 163 additions & 115 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,8 @@
11
/**
22
*
3-
* Fitter package for ALERT
3+
* Fitter package for CND
44
*
5-
* Writer: Sangbaek Lee, Zhiwan Xu
5+
* Writer: Sangbaek Lee
66
*
77
**/
88
package org.jlab.clas.timeline.fitter
@@ -41,119 +41,167 @@ class ALERTFitter{
4141
}
4242

4343
static F1D tdc_minus_start_time_fitter(H1F h1, int component){
44-
if (component > 9) {//bars
45-
def f1 =new F1D("fit:"+h1.getName(),"[amp]*gaus(x,[mean],[sigma])+[cst]", -5.0, 5.0);
46-
f1.setLineColor(33);
47-
f1.setLineWidth(10);
48-
f1.setOptStat("1111");
49-
double maxz = h1.getBinContent(h1.getMaximumBin());
50-
double peak_location = h1.getAxis().getBinCenter(h1.getMaximumBin());
51-
f1.setRange(peak_location - 5, peak_location + 5);
52-
f1.setParameter(0,maxz-h1.getBinContent(0));
53-
f1.setParameter(1, peak_location);
54-
f1.setParameter(2, 1.0);
55-
f1.setParameter(3, h1.getBinContent(0));
56-
if (maxz>0) f1.setParLimits(0, maxz*0.9,maxz*1.1);
57-
f1.setParLimits(3, 0.0, 0.1*maxz);
44+
if(component>9){//bars
45+
def f1 =new F1D("fit:"+h1.getName(),"[amp]*gaus(x,[mean],[sigma])+[cst]", -5.0, 5.0);
46+
f1.setLineColor(33);
47+
f1.setLineWidth(10);
48+
f1.setOptStat("1111");
49+
double maxz = h1.getBinContent(h1.getMaximumBin());
50+
double peak_location = h1.getAxis().getBinCenter(h1.getMaximumBin());
51+
f1.setRange(peak_location - 5, peak_location + 5);
52+
f1.setParameter(0,maxz-h1.getBinContent(0));
53+
f1.setParameter(1, peak_location);
54+
f1.setParameter(2, 1.0);
55+
f1.setParameter(3, h1.getBinContent(0));
56+
if (maxz>0) f1.setParLimits(0, maxz*0.9,maxz*1.1);
57+
f1.setParLimits(3, 0.0, 0.1*maxz);
58+
59+
double hMean, hRMS
60+
def originalOut = System.out
61+
System.setOut(new PrintStream(OutputStream.nullOutputStream())) // Java 11+
62+
63+
// Code that prints to System.out
64+
DataFitter.fit(f1, h1, "");
65+
66+
System.setOut(originalOut) // Restore the original output
67+
68+
return f1
69+
}
70+
else {
71+
double height_fit_set = 0
72+
double mean_fit_set = 0
73+
double sigma_fit_set = 0
74+
double step = 2
75+
76+
// ----- Clone & analyze primary peak -----
77+
H1F hpy = h1.histClone("hpy_zoom")
78+
int maxBin = hpy.getMaximumBin()
79+
double maxY = hpy.getBinContent(maxBin)
80+
double peak = hpy.getXaxis().getBinCenter(maxBin)
81+
double sigma0 = Math.min(step, getRestrictedRMS(hpy, peak - step, peak + step))
82+
83+
// ----- Primary Gaussian Fit -----
84+
F1D fgaus = new F1D("fgaus", "[amp]*gaus(x,[mean],[sigma])",
85+
peak - step, peak + step)
86+
87+
fgaus.setParameter(0, maxY)
88+
fgaus.setParameter(1, peak)
89+
if (sigma0 != null && sigma0 > 0) {
90+
fgaus.setParameter(2, sigma0)
91+
} else {
92+
fgaus.setParameter(2, 0.8)
93+
}
94+
fgaus.setParLimits(0, 0, 1.2 * maxY)
95+
fgaus.setParLimits(1, peak - step, peak + step)
96+
fgaus.setParLimits(2, 0, step)
97+
98+
PrintStream original = System.out
99+
System.setOut(new PrintStream(OutputStream.nullOutputStream()))
100+
DataFitter.fit(fgaus, hpy, "RQ")
101+
System.setOut(original)
102+
103+
double height = fgaus.getParameter(0)
104+
double mean = fgaus.getParameter(1)
105+
double sigma = fgaus.getParameter(2)
106+
double entriesTotal = hpy.integral()
107+
108+
// -------------------------------------------------------------------------
109+
// Utility closure to fit a left-side peak after cutting the histogram
110+
// -------------------------------------------------------------------------
111+
def fitLeftPeak = { H1F h, double prevMean, double prevSigma ->
112+
// Cut histogram to the left of the previous peak
113+
H1F hcut = h.histClone("hcut")
114+
int cutBin = hcut.getXaxis().getBin(prevMean - prevSigma * 2)
115+
for (int b = cutBin; b <= hcut.getXaxis().getNBins(); b++) {
116+
hcut.setBinContent(b, 0)
117+
hcut.setBinError(b, 0)
118+
}
119+
120+
int mb = hcut.getMaximumBin()
121+
if (mb >= cutBin - 2) return null // ensure left-side peak
122+
123+
double pk = hcut.getXaxis().getBinCenter(mb)
124+
double amp0 = hcut.getBinContent(mb)
125+
126+
F1D ftmp = new F1D("fgaus_left",
127+
"[amp]*gaus(x,[mean],[sigma])",
128+
pk - step, pk + step)
129+
130+
ftmp.setParameter(0, amp0)
131+
ftmp.setParameter(1, pk)
132+
ftmp.setParameter(2, 0.8)
133+
134+
ftmp.setParLimits(0, 0, 1.2 * amp0)
135+
ftmp.setParLimits(1, pk - step, pk + step)
136+
ftmp.setParLimits(2, 0, step)
137+
138+
System.setOut(new PrintStream(OutputStream.nullOutputStream()))
139+
DataFitter.fit(ftmp, hcut, "RQ")
140+
System.setOut(original)
141+
142+
double A = ftmp.getParameter(0)
143+
double M = ftmp.getParameter(1)
144+
double S = ftmp.getParameter(2)
145+
146+
// validation conditions
147+
if (A > maxY * 0.3 &&
148+
S < step && S > 0.1 &&
149+
M < prevMean - prevSigma &&
150+
hcut.integral() > 0.05 * entriesTotal)
151+
{
152+
return [A, M, S]
153+
}
154+
return null
155+
}
156+
157+
// -------------------------------------------------------------------------
158+
// Peak 1 (primary) already known: height, mean, sigma
159+
// Try Peak 2 (first left peak)
160+
// -------------------------------------------------------------------------
161+
def p2 = fitLeftPeak(hpy, mean, sigma)
162+
163+
// If Peak 2 exists, try Peak 3 (second left peak)
164+
def p3 = p2 ? fitLeftPeak(hpy, p2[1], p2[2]) : null
165+
166+
// If Peak 3 exists, try Peak 4 (third left peak)
167+
def p4 = p3 ? fitLeftPeak(hpy, p3[1], p3[2]) : null
168+
169+
// -------------------------------------------------------------------------
170+
// Choose deepest detected peak on the left: p4 > p3 > p2 > primary
171+
// -------------------------------------------------------------------------
172+
if (p4) {
173+
height_fit_set = p4[0]
174+
mean_fit_set = p4[1]
175+
sigma_fit_set = p4[2]
176+
} else if (p3) {
177+
height_fit_set = p3[0]
178+
mean_fit_set = p3[1]
179+
sigma_fit_set = p3[2]
180+
} else if (p2) {
181+
height_fit_set = p2[0]
182+
mean_fit_set = p2[1]
183+
sigma_fit_set = p2[2]
184+
} else {
185+
height_fit_set = height
186+
mean_fit_set = mean
187+
sigma_fit_set = sigma
188+
}
189+
190+
// Return final selected Gaussian
191+
F1D fout = new F1D("fit:" + h1.getName(),
192+
"[amp]*gaus(x,[mean],[sigma])",
193+
mean_fit_set - step * 2, mean_fit_set + step * 2)
194+
195+
fout.setParameter(0, height_fit_set)
196+
fout.setParameter(1, mean_fit_set)
197+
fout.setParameter(2, sigma_fit_set)
198+
199+
fout.setLineColor(33)
200+
fout.setLineWidth(10)
201+
202+
return fout
203+
}
58204

59-
double hMean, hRMS
60-
def originalOut = System.out
61-
System.setOut(new PrintStream(OutputStream.nullOutputStream())) // Java 11+
62-
63-
// Code that prints to System.out
64-
DataFitter.fit(f1, h1, "");
65-
66-
System.setOut(originalOut) // Restore the original output
67-
68-
return f1
69-
}
70-
else {//wedges
71-
double height_fit_set = 0
72-
double mean_fit_set = 0
73-
double sigma_fit_set = 0
74-
double step = 2
75-
76-
// ----- Clone & analyze primary peak -----
77-
H1F hpy = h1.histClone("hpy_zoom")
78-
int maxBin = hpy.getMaximumBin()
79-
double maxY = hpy.getBinContent(maxBin)
80-
double peak = hpy.getXaxis().getBinCenter(maxBin)
81-
double sigma0 = Math.min(step, getRestrictedRMS(hpy, peak - step, peak + step))
82-
83-
// ----- Primary Gaussian Fit -----
84-
F1D fgaus = new F1D("fgaus", "[amp]*gaus(x,[mean],[sigma])",
85-
peak - step, peak + step)
86-
87-
fgaus.setParameter(0, maxY)
88-
fgaus.setParameter(1, peak)
89-
fgaus.setParameter(2, sigma0 > 0 ? sigma0 : 0.8)
90-
fgaus.setParLimits(0, 0, 1.2 * maxY)
91-
fgaus.setParLimits(1, peak - step, peak + step)
92-
fgaus.setParLimits(2, 0, step)
93-
94-
def original = System.out
95-
System.setOut(new PrintStream(OutputStream.nullOutputStream()))
96-
DataFitter.fit(fgaus, hpy, "RQ")
97-
System.setOut(original)
98-
99-
double height = fgaus.getParameter(0)
100-
double mean = fgaus.getParameter(1)
101-
double sigma = fgaus.getParameter(2)
102-
double entriesTotal = hpy.integral()
103-
104-
// -------------------------------------------------------------------------
105-
// Utility closure to fit a left-side peak after cutting the histogram
106-
// -------------------------------------------------------------------------
107-
def fitLeftPeak = { H1F h, double prevMean, double prevSigma ->
108-
// Cut histogram to the left of the previous peak
109-
H1F hcut = h.histClone("hcut")
110-
int cutBin = hcut.getXaxis().getBin(prevMean - prevSigma * 2)
111-
for (int b = cutBin; b <= hcut.getXaxis().getNBins(); b++) {
112-
hcut.setBinContent(b, 0)
113-
hcut.setBinError(b, 0)
114-
}
115-
116-
if (hcut.integral() < 3) return null
117-
118-
int mb = hcut.getMaximumBin()
119-
if (mb >= cutBin - 2) return null // ensure left-side peak
120-
121-
double pk = hcut.getXaxis().getBinCenter(mb)
122-
double amp0 = hcut.getBinContent(mb)
123-
124-
F1D ftmp = new F1D("fgaus_left",
125-
"[amp]*gaus(x,[mean],[sigma])",
126-
pk - step, pk + step)
127-
128-
ftmp.setParameter(0, amp0)
129-
ftmp.setParameter(1, pk)
130-
ftmp.setParameter(2, 0.8)
131-
132-
ftmp.setParLimits(0, 0, 1.2 * amp0)
133-
ftmp.setParLimits(1, pk - step, pk + step)
134-
ftmp.setParLimits(2, 0, step)
135-
136-
System.setOut(new PrintStream(OutputStream.nullOutputStream()))
137-
DataFitter.fit(ftmp, hcut, "RQ")
138-
System.setOut(original)
139-
140-
double A = ftmp.getParameter(0)
141-
double M = ftmp.getParameter(1)
142-
double S = ftmp.getParameter(2)
143-
144-
// validation conditions
145-
if (A > maxY * 0.3 &&
146-
S < step && S > 0.1 &&
147-
M < prevMean - prevSigma &&
148-
hcut.integral() > 0.05 * entriesTotal)
149-
{
150-
return [A, M, S]
151-
}
152-
return null
153-
}
154-
155-
return null
156-
}
157205
}
158206

159207

@@ -217,4 +265,4 @@ def fitLeftPeak = { H1F h, double prevMean, double prevSigma ->
217265

218266
return f1
219267
}
220-
}
268+
}

0 commit comments

Comments
 (0)
⚡