Skip to content

Commit cfe6e9b

Browse files
zhiwanxuzhiwanxu
authored
Issue 417 (#420)
1 parent 41ccb11 commit cfe6e9b

2 files changed

Lines changed: 165 additions & 26 deletions

File tree

‎src/main/java/org/jlab/clas/timeline/analysis/alert/alert_atof_tdc_minus_start_time.groovy‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -41,7 +41,7 @@ int max_index;
4141
if(h1!=null) {
4242
if (h1.getBinContent(h1.getMaximumBin()) > 30 && h1.getEntries()>300){
4343
data[run].put(String.format('atof_tdc_minus_start_time_%s', file_index), h1)
44-
def f1 = ALERTFitter.tdc_minus_start_time_fitter(h1)
44+
def f1 = ALERTFitter.tdc_minus_start_time_fitter(h1,component)
4545
data[run].put(String.format('fit_atof_tdc_minus_start_time_%s', file_index), f1)
4646
data[run].put(String.format('peak_location_atof_tdc_minus_start_time_%s', file_index), f1.getParameter(1).abs())
4747
data[run].put(String.format('sigma_atof_tdc_minus_start_time_%s', file_index), f1.getParameter(2).abs())

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

Lines changed: 164 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,8 @@
11
/**
22
*
3-
* Fitter package for CND
3+
* Fitter package for ALERT
44
*
5-
* Writer: Sangbaek Lee
5+
* Writer: Sangbaek Lee, Zhiwan Xu
66
*
77
**/
88
package org.jlab.clas.timeline.fitter
@@ -40,31 +40,170 @@ class ALERTFitter{
4040
return f1
4141
}
4242

43-
static F1D tdc_minus_start_time_fitter(H1F h1){
44-
def f1 =new F1D("fit:"+h1.getName(),"[amp]*gaus(x,[mean],[sigma])+[cst]", -5.0, 5.0);
45-
f1.setLineColor(33);
46-
f1.setLineWidth(10);
47-
f1.setOptStat("1111");
48-
double maxz = h1.getBinContent(h1.getMaximumBin());
49-
double peak_location = h1.getAxis().getBinCenter(h1.getMaximumBin());
50-
f1.setRange(peak_location - 5, peak_location + 5);
51-
f1.setParameter(0,maxz-h1.getBinContent(0));
52-
f1.setParameter(1, peak_location);
53-
f1.setParameter(2, 1.0);
54-
f1.setParameter(3, h1.getBinContent(0));
55-
if (maxz>0) f1.setParLimits(0, maxz*0.9,maxz*1.1);
56-
f1.setParLimits(3, 0.0, 0.1*maxz);
43+
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);
5758

58-
double hMean, hRMS
59-
def originalOut = System.out
60-
System.setOut(new PrintStream(OutputStream.nullOutputStream())) // Java 11+
61-
62-
// Code that prints to System.out
63-
DataFitter.fit(f1, h1, "");
59+
double hMean, hRMS
60+
def originalOut = System.out
61+
System.setOut(new PrintStream(OutputStream.nullOutputStream())) // Java 11+
6462

65-
System.setOut(originalOut) // Restore the original output
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+
if (hcut.integral() < 3) return null
121+
122+
int mb = hcut.getMaximumBin()
123+
if (mb >= cutBin - 2) return null // ensure left-side peak
124+
125+
double pk = hcut.getXaxis().getBinCenter(mb)
126+
double amp0 = hcut.getBinContent(mb)
127+
128+
F1D ftmp = new F1D("fgaus_left",
129+
"[amp]*gaus(x,[mean],[sigma])",
130+
pk - step, pk + step)
131+
132+
ftmp.setParameter(0, amp0)
133+
ftmp.setParameter(1, pk)
134+
ftmp.setParameter(2, 0.8)
135+
136+
ftmp.setParLimits(0, 0, 1.2 * amp0)
137+
ftmp.setParLimits(1, pk - step, pk + step)
138+
ftmp.setParLimits(2, 0, step)
139+
140+
System.setOut(new PrintStream(OutputStream.nullOutputStream()))
141+
DataFitter.fit(ftmp, hcut, "RQ")
142+
System.setOut(original)
143+
144+
double A = ftmp.getParameter(0)
145+
double M = ftmp.getParameter(1)
146+
double S = ftmp.getParameter(2)
147+
148+
// validation conditions
149+
if (A > maxY * 0.3 &&
150+
S < step && S > 0.1 &&
151+
M < prevMean - prevSigma &&
152+
hcut.integral() > 0.05 * entriesTotal)
153+
{
154+
return [A, M, S]
155+
}
156+
return null
157+
}
158+
159+
// -------------------------------------------------------------------------
160+
// Peak 1 (primary) already known: height, mean, sigma
161+
// Try Peak 2 (first left peak)
162+
// -------------------------------------------------------------------------
163+
def p2 = fitLeftPeak(hpy, mean, sigma)
164+
165+
// If Peak 2 exists, try Peak 3 (second left peak)
166+
def p3 = p2 ? fitLeftPeak(hpy, p2[1], p2[2]) : null
167+
168+
// If Peak 3 exists, try Peak 4 (third left peak)
169+
def p4 = p3 ? fitLeftPeak(hpy, p3[1], p3[2]) : null
170+
171+
// -------------------------------------------------------------------------
172+
// Choose deepest detected peak on the left: p4 > p3 > p2 > primary
173+
// -------------------------------------------------------------------------
174+
if (p4) {
175+
height_fit_set = p4[0]
176+
mean_fit_set = p4[1]
177+
sigma_fit_set = p4[2]
178+
} else if (p3) {
179+
height_fit_set = p3[0]
180+
mean_fit_set = p3[1]
181+
sigma_fit_set = p3[2]
182+
} else if (p2) {
183+
height_fit_set = p2[0]
184+
mean_fit_set = p2[1]
185+
sigma_fit_set = p2[2]
186+
} else {
187+
height_fit_set = height
188+
mean_fit_set = mean
189+
sigma_fit_set = sigma
190+
}
191+
192+
// Return final selected Gaussian
193+
F1D fout = new F1D("fit:" + h1.getName(),
194+
"[amp]*gaus(x,[mean],[sigma])",
195+
mean_fit_set - step * 2, mean_fit_set + step * 2)
196+
197+
fout.setParameter(0, height_fit_set)
198+
fout.setParameter(1, mean_fit_set)
199+
fout.setParameter(2, sigma_fit_set)
200+
201+
fout.setLineColor(33)
202+
fout.setLineWidth(10)
203+
204+
return fout
205+
}
66206

67-
return f1
68207
}
69208

70209

@@ -128,4 +267,4 @@ class ALERTFitter{
128267

129268
return f1
130269
}
131-
}
270+
}

0 commit comments

Comments
 (0)
⚡