81 const string myname =
"StickyCodeMetrics::evaluateMetrics: ";
83 if ( dbg ) cout << myname <<
"Debugging level: " << dbg <<
endl;
86 Index badadc = 999999;
87 Index maxadc = badadc;
94 double countSum = 0.0;
95 for ( BinCounter::value_type ibco : counts ) {
96 Index iadc = ibco.first;
97 double count = ibco.second;
99 if ( count > maxCount ) {
104 if ( adcmod == 0 ) nmod0 +=
count;
105 if ( adcmod == 1 ) nmod1 +=
count;
106 if ( adcmod == 63 ) nmod63 +=
count;
108 adcSum += count*iadc;
110 double countSum2 = 0;
112 Index maxadc2 = badadc;
113 for ( BinCounter::value_type ibco : counts ) {
114 Index iadc = ibco.first;
115 double count = ibco.second;
116 if ( iadc != maxadc ) {
118 adcSum2 += count*iadc;
119 if ( count > maxCount2 ) {
127 m_meanAdc = countSum>0 ? adcSum/countSum : -1.0;
128 m_meanAdc2 = countSum2>0 ? adcSum2/countSum2 : -1.0;
134 if ( counts.size() == 0 )
return 0;
135 if (
m_hnam.size() == 0 )
return 0;
136 if (
m_nbin == 0 )
return 1;
144 Index idadc1 = counts.begin()->first;
145 Index idadc2 = counts.rbegin()->first;
150 Index halfHist = nbinHist/2;
151 Index iadc1 = idadc1 < halfHist ? 0 : idadc1 - halfHist;
152 iadc1 = nrebin*(iadc1/nrebin);
153 Index binOffset = iadc1;
154 Index iadc2= idadc2 + halfHist;
155 if ( iadc2 > 4096 && idadc2 < 4096 ) iadc2 = 4096;
157 for (
Index iadc=iadc1; iadc<=iadc2; ++iadc ) {
159 double count = icnt != counts.end() ? fabs(icnt->second) : 0;
161 binCounts.push_back(count);
164 Index nreb = (nbin+nrebin-1)/nrebin;
167 for (
Index ibin=0; ibin<nbin; ++ibin ) {
168 double count = binCounts[ibin];
169 Index iadc = binOffset + ibin;
171 Index ireb = ibin/nrebin;
172 rebinCounts[ireb] +=
count;
173 rebinAdcCounts[ireb] += iadc*
count;
177 cout << myname <<
" Data ADC range: [" << idadc1 <<
", " << idadc2 <<
"]" <<
endl;
178 cout << myname <<
" Bin offset: " << binOffset <<
endl;
179 cout << myname <<
" Data bin count: " << nbin <<
endl;
180 cout << myname <<
" Data rebin count: " << rebinCounts.size() <<
endl;
181 cout << myname <<
" Hist bin count: " << nbinHist <<
endl;
184 double countUnder = 0.0;
185 double countOver = 0.0;
188 Index ihadc1 = binOffset;
189 while ( ihadc1 + nrebin <= idadc1 ) ihadc1 += nrebin;
190 Index ihadc2 = ihadc1 + nbinHist;
192 if ( ihadc2 > idadc2 ) {
193 Index nbinLeft = idadc1 - ihadc1;
194 Index nbinRight = ihadc2 - idadc2 - 1;
195 if ( dbg >= 3 ) cout << myname <<
"All data are in histogram range." <<
endl;
196 if ( dbg >= 4 ) cout << myname <<
"nl, nr, hleft: " << nbinLeft <<
", " 197 << nbinRight <<
", " << ihadc1 <<
endl;
198 while ( nbinRight > nbinLeft + nrebin ) {
203 if ( dbg >= 4 ) cout << myname <<
"nl, nr, hleft: " << nbinLeft <<
", " 204 << nbinRight <<
", " << ihadc1 <<
endl;
205 if ( ihadc1 < nrebin )
break;
210 ihadc2 = ihadc1 + nbinHist;
211 if ( dbg >= 3 ) cout << myname <<
"Some data are outside histogram range." <<
endl;
213 double adcHistCenter =
maxAdc();
218 cout << myname <<
" max ADC: " <<
maxAdc() <<
endl;
219 cout << myname <<
" max ADC2: " <<
maxAdc() <<
endl;
220 cout << myname <<
" mean ADC: " <<
meanAdc() <<
endl;
221 cout << myname <<
" Preferred hist center: " << adcHistCenter <<
endl;
224 double countSumHist = 0.0;
225 double adcSumHist = 0.0;
226 for (
Index iadc=ihadc1; iadc<ihadc2; ++iadc ) {
227 Index ibin = iadc - binOffset;
228 countSumHist += binCounts[ibin];
229 adcSumHist += binCounts[ibin]*iadc;
233 Index ihadc1Selected = ihadc1;
234 double countSumHistSelected = countSumHist;
235 double hmean = countSumHist == 0.0 ? 0.0 : adcSumHist/countSumHist;
236 double hmeanSelected = hmean;
237 if ( dbg >= 4 ) cout <<
" hleft, hmean, area: " << ihadc1 <<
", " << hmean
238 <<
", " << countSumHist <<
endl;
239 while ( ihadc1 <= idadc2 ) {
241 for (
Index iadc=ihadc2; iadc<ihadc2+nrebin; ++iadc ) {
242 if ( iadc > idadc2 )
break;
243 Index ibin = iadc - binOffset;
244 countSumHist += binCounts[ibin];
245 adcSumHist += iadc*binCounts[ibin];
248 Index ireb = (ihadc1-binOffset)/nrebin;
249 if ( dbg >= 4 ) cout <<
" ireb/nreb: " << ireb <<
"/" << nrebin <<
endl;
250 countSumHist -= rebinCounts[ireb];
251 adcSumHist -= rebinAdcCounts[ireb];
252 hmean = countSumHist == 0.0 ? 0.0 : adcSumHist/countSumHist;
256 if ( dbg >= 4 ) cout <<
" hleft, hmean, area: " << ihadc1 <<
", " << hmean
257 <<
", " << countSumHist <<
endl;
260 double diffCount = countSumHist - countSumHistSelected;
261 bool sameCount = diffCount == 0.0;
263 double rat = diffCount/(countSumHist + countSumHistSelected);
264 sameCount = fabs(rat) < 0.001;
266 if ( ihadc1Selected > ihadc1 ) {
267 cout << myname <<
"Unexpected value for ihadc1Selected: " << ihadc1Selected <<
endl;
268 cout << myname <<
" with ihadc1: " << ihadc1 <<
endl;
271 bool haveOverlap = ihadc1Selected + nbinHist > ihadc1;
272 bool updateSel =
false;
273 if ( sameCount && haveOverlap ) {
274 if ( dbg >= 4 ) cout <<
" Same count." <<
endl;
276 double diffOld = fabs(ihadc1Selected + nbinHist/2 - hmeanSelected);
277 double diffNew = fabs(ihadc1 + nbinHist/2 - hmean);
278 updateSel = diffNew < diffOld;
280 double diffOld = fabs(ihadc1Selected + nbinHist/2 - adcHistCenter);
281 double diffNew = fabs(ihadc1 + nbinHist/2 - adcHistCenter);
282 updateSel = diffNew < diffOld;
285 updateSel = diffCount > 0.0;
288 countSumHistSelected = countSumHist;
289 hmeanSelected = hmean;
290 ihadc1Selected = ihadc1;
291 if ( dbg >= 4 ) cout <<
" Updated selected hleft, hmean to " << ihadc1Selected
292 <<
", " << hmeanSelected <<
endl;
295 ihadc1 = ihadc1Selected;
296 ihadc2 = ihadc1 + nbinHist;
298 for ( BinCounter::value_type icnt : counts ) {
299 Index iadc = icnt.first;
300 double count = icnt.second;
301 if ( iadc < ihadc1 ) countUnder +=
count;
302 if ( iadc >= ihadc2 ) countOver +=
count;
306 TH1* ph =
new TH1F(
m_hnam.c_str(),
m_httl.c_str(), nbinHist, ihadc1, ihadc2);
307 string::size_type iposx =
m_httl.find(
";");
308 bool haveXlab = iposx != string::npos;
309 bool haveYlab = haveXlab &&
m_httl.find(
";", iposx+1) != string::npos;
310 if ( ! haveXlab ) ph->GetXaxis()->SetTitle(
"ADC count");
311 if ( ! haveYlab ) ph->GetYaxis()->SetTitle(
"# samples");
316 for (
Index iadc=binOffset; iadc<=idadc2; ++iadc ) {
318 Index count = icnt != counts.end() ? fabs(icnt->second) : 0;
319 ph->SetBinContent(iadc+1-ihadc1, count);
321 ph->SetBinContent(0, countUnder);
322 ph->SetBinContent(nbinHist+1, countOver);
325 TF1 fitter(
"pedgaus",
"gaus", ihadc1, ihadc2, TF1::EAddToList::kNo);
326 fitter.SetParameters(0.1*ph->Integral(), ph->GetMean(), 5.0);
332 fitter.SetParLimits(2, 0.1, 100.0);
334 TF1* pfinit =
dynamic_cast<TF1*
>(fitter.Clone(
"pedgaus0"));
335 pfinit->SetLineColor(3);
336 pfinit->SetLineStyle(2);
346 ErrorHandlerFunc_t pehSave =
nullptr;
347 ErrorHandlerFunc_t pehDefault = DefaultErrorHandler;
348 if ( GetErrorHandler() != pehDefault ) {
349 pehSave = SetErrorHandler(pehDefault);
354 if ( pehSave !=
nullptr ) SetErrorHandler(pehSave);
356 ph->GetListOfFunctions()->AddLast(pfinit,
"0");
357 ph->GetListOfFunctions()->Last()->SetBit(TF1::kNotDraw,
true);
counts_as<> counts
Number of ADC counts, represented by signed short int.
StickyCodeMetrics::BinCounter BinCounter
QTextStream & endl(QTextStream &s)