26 #include "TPaveText.h"
61 double fitMax,
int lwb,
int upb,
unsigned int chi2Pars,
64 TH1D*
h =
new TH1D(
name,
"#chi^{2} distribution",
nBins, 0,
max);
65 unsigned int nDofEff = 0, nDofLast = 0;
66 std::unique_ptr<UniformShapeErrorGetter> kFactorGetter;
67 std::unique_ptr<CombinedShapeErrorGetter> shapeErrorGetter;
69 kFactorGetter = std::make_unique<UniformShapeErrorGetter>(
kFactor);
70 shapeErrorGetter = std::make_unique<CombinedShapeErrorGetter>();
72 shapeErrorGetter->
add(*kFactorGetter);
75 if ((
i+1) % 10000 == 0) cout <<
"Processing entry # " <<
i+1 << endl;
77 if (!history)
continue;
79 for (
unsigned int j = 0; j < history->
nData(); j++) {
80 h->Fill(history->
chi2(j, lwb, upb, chi2Pars, shapeErrorType, &nDofEff));
81 if (nDofLast > 0 && nDofLast != nDofEff) cout <<
"WARNING @ hash = " <<
i <<
", index = " << j <<
" : nDof varied from " << nDofLast <<
" to " << nDofEff << endl;
86 if (nDof > 0 && nDof != nDofEff) cout <<
"WARNING : fixing nDof = " << nDof <<
", but computer nDofEff = " << nDofEff << endl;
88 if (fitMax < 0)
return h;
89 TF1* fChi2 =
fitChi2(*
h,
"chi2", 0,
max, nDof, 0, fitMax);
90 TPaveText*
p =
new TPaveText(0.6, 0.85, 0.85, 0.65,
"NDC");
91 p->AddText(Form(
"k = %.1f%%",
kFactor*100));
92 p->AddText(Form(
"Fit #chi^{2} = %.1f", fChi2->GetChisquare()));
93 p->AddText(Form(
"Fit nDOF = %.1f", fChi2->GetParameter(1)));
101 TF1* fChi2 =
new TF1(
name,
"[0]*TMath::Power(x/2, [1]/2)/x*TMath::Exp(-x/2)/TMath::Gamma([1]/2)", xMin, xMax);
102 fChi2->SetParameter(0, 100);
103 fChi2->SetParameter(1, 3);
104 fChi2->SetParLimits(1, 1, 7);
110 double nDof,
double fitMin,
double fitMax)
116 fChi2->SetParameter(0, 4*
h.GetMaximum());
117 if (nDof > 0) fChi2->FixParameter(1, nDof);
119 h.Fit(fChi2,
"",
"", fitMin, fitMax);
120 fChi2->SetRange(xMin, xMax);
126 double refErrMin,
double refErrMax,
unsigned int refErrNBins,
127 int lwb,
int upb,
unsigned int chi2Pars,
unsigned int nDof)
const
130 if (nDof > 0) fChi2->FixParameter(1, nDof);
132 double delta = (refErrMax - refErrMin)/refErrNBins;
133 TH1D* fitResults =
new TH1D(
"fitResults",
"", refErrNBins, refErrMin - delta/2, refErrMax - delta/2);
135 for (
unsigned int k = 0;
k < refErrNBins;
k++) {
136 double refErr = refErrMin + delta*
k;
137 TH1D*
h =
new TH1D(
name,
"#chi^{2} distribution",
nBins, 0,
max);
141 shapeErrorGetter.
add(kFactorGetter);
144 if (!history)
continue;
146 for (
unsigned int j = 0; j < history->
nData(); j++) {
147 h->Fill(history->
chi2(j, lwb, upb, chi2Pars));
151 fitResults->SetBinContent(
k+1, fChi2->GetChisquare());
155 TF1* fPol2 =
new TF1(
"fPol2",
"[0]*(x - [1])*(x - [1]) + [2]", refErrMin, refErrMax);
157 double center = fitResults->GetBinCenter(fitResults->GetMinimumBin());
158 fPol2->SetParameter(0, (fitResults->GetMaximum() - fitResults->GetMinimum())/center/center);
159 fPol2->SetParameter(1, center);
160 fPol2->SetParameter(2, fitResults->GetMinimum());
163 fitResults->Fit(fPol2);
165 TPaveText*
p =
new TPaveText(0.6, 0.85, 0.85, 0.65,
"NDC");
166 p->AddText(Form(
"Fit k_{best} = %.1f%%", fPol2->GetParameter(1)*100));
167 p->AddText(Form(
"Fit #delta k_{best} = %.2f%%", 1/TMath::Sqrt(fPol2->GetParameter(0))*100));
168 p->AddText(Form(
"Fit #chi^{2}_{min} = %.2f", fPol2->GetParameter(2)));
178 TH1D*
h =
new TH1D(
name,
"gain distribution", 3, -0.5, 2.5);
179 h->GetXaxis()->SetBinLabel(1,
"Low");
180 h->GetXaxis()->SetBinLabel(2,
"Medium");
181 h->GetXaxis()->SetBinLabel(3,
"High");
184 if (!history)
continue;
185 for (
unsigned int j = 0; j < history->
nData(); j++)
h->Fill(history->
data(j)->
gain());
193 TH1D*
h =
new TH1D(
name,
"layer distribution", 4, -0.5, 3.5);
194 h->GetXaxis()->SetBinLabel(1,
"PS");
195 h->GetXaxis()->SetBinLabel(2,
"Strips");
196 h->GetXaxis()->SetBinLabel(3,
"Middle");
197 h->GetXaxis()->SetBinLabel(4,
"Back");
200 if (!history)
continue;
209 TH1D*
h =
new TH1D(
name,
"Residual distribution",
nBins, rMin, rMax);
210 h->GetXaxis()->SetTitle(
norm ?
"Residual/ADCMax (%)" :
"Residual (ADC counts)");
214 if (!history)
continue;
215 for (
unsigned int j = 0; j < history->
nData(); j++)
223 double lo,
double hi,
const TString&
fileName)
227 TString globalTitle = TString(xip ?
"Derivative" :
"Offset") +
" correction";
228 for (
unsigned short k = 0;
k < 5;
k++) {
229 TString
title = globalTitle +
" for sampling " + Form(
"%d",
k);
230 h1[
k] =
new TH1D(Form(
"h%d",
k),
title, 100, lo, hi);
231 h1[
k]->GetXaxis()->SetTitle(
"correction");
233 hEtaPhi[
k]->SetMinimum(lo);
234 hEtaPhi[
k]->SetMaximum(hi);
239 TH2D* hAll =
new TH2D(
"hAll", globalTitle +
" vs. sampling", 5, 0, 5, 100, lo, hi);
240 hAll->GetXaxis()->SetTitle(
"sampling");
241 hAll->GetYaxis()->SetTitle(
"correction");
243 if ((
i+1) % 10000 == 0) cout <<
"Cell # " <<
i+1 <<
"/" <<
nChannels() << endl;
245 if (!history || history->
nData() == 0)
continue;
246 for (
unsigned short k = 0;
k < 5;
k++) {
255 for (
unsigned int i = 1;
i <= 100;
i++)
256 for (
unsigned short k = 0;
k < 5;
k++)
257 hAll->SetBinContent(
k+1,
i,
h1[
k]->GetBinContent(
i));
258 TFile*
ff = TFile::Open(
fileName +
".root",
"RECREATE");
259 for (
unsigned short k = 0;
k < 5;
k++) {
265 TCanvas*
c1 =
new TCanvas(
"c1",
Id::str(calo) +
" layer " + Form(
"%d",
layer), 1600, 2000);
267 for (
unsigned short k = 0;
k < 5;
k++) {
269 hEtaPhi[
k]->Draw(
"COLZ");
280 double lo,
double hi,
const TString&
fileName)
282 unsigned int nBins = 100;
284 TString globalTitle = TString(xip ?
"Derivative" :
"Offset") +
" correction";
285 for (
unsigned short l = 0;
l < 5;
l++)
286 for (
unsigned short k = 0;
k < 5;
k++)
287 h[
l][
k] =
new TH1D(Form(
"h%d%d",
l,
k),
"",
nBins,
l < 2 ? -15 : lo,
l < 2 ? 15 : hi);
292 hAll[0] =
new TH2D(
"h0All",
"#xi_{cell}/#delta#xi vs. sampling", 5, 0, 5,
nBins, -15, 15);
293 hAll[1] =
new TH2D(
"hDAll",
"(#xi_{cell} - #xi_{ring})/#delta#xi vs. sampling", 5, 0, 5,
nBins, -15, 15);
294 hAll[2] =
new TH2D(
"hEAll", globalTitle +
" per-cell value vs. sampling", 5, 0, 5,
nBins, 0, hi);
295 hAll[3] =
new TH2D(
"hRAll", globalTitle +
" per-ring value vs. sampling", 5, 0, 5,
nBins, lo, hi);
296 hAll[4] =
new TH2D(
"hCAll", globalTitle +
" per-cell error vs. sampling", 5, 0, 5,
nBins, lo, hi);
297 TH2D* hCor =
new TH2D(
"hCor", globalTitle +
" per-ring vs. per-cell values",
nBins, lo, hi,
nBins, lo, hi);
298 hCor->GetXaxis()->SetTitle(
"#xi_{cell}");
299 hCor->GetYaxis()->SetTitle(
"#xi_{ring}");
300 for (
unsigned short l = 0;
l < 4;
l++) {
301 hAll[
l]->GetXaxis()->SetTitle(
"sampling");
302 hAll[
l]->GetYaxis()->SetTitle(
"significance");
305 if ((
i+1) % 10000 == 0) cout <<
"Cell # " <<
i+1 <<
"/" <<
nChannels() << endl;
307 if (!history || history->
nData() == 0)
continue;
308 for (
unsigned short k = 0;
k < 5;
k++) {
311 if (!cellSED || !ringSED) {
312 if (cellSED)
delete cellSED;
313 if (ringSED)
delete ringSED;
316 double cellVal = (xip ? cellSED->
xi()(
k) : cellSED->
xip()(
k));
317 double ringVal = (xip ? ringSED->
xi()(
k) : ringSED->
xip()(
k));
318 double cellErr = TMath::Sqrt(xip ? cellSED->
xiErr()(
k,
k) : cellSED->
xipErr()(
k,
k));
321 h[0][
k]->Fill(cellVal/cellErr);
322 h[1][
k]->Fill((cellVal - ringVal)/cellErr);
323 h[2][
k]->Fill(cellVal);
324 h[3][
k]->Fill(ringVal);
325 h[4][
k]->Fill(cellErr);
326 hCor->Fill(cellVal, ringVal);
329 for (
unsigned int i = 1;
i <=
nBins;
i++)
330 for (
unsigned short l = 0;
l < 5;
l++)
331 for (
unsigned short k = 0;
k < 5;
k++)
332 hAll[
l]->SetBinContent(
k+1,
i,
h[
l][
k]->GetBinContent(
i));
333 TFile*
ff = TFile::Open(
fileName +
".root",
"RECREATE");
334 for (
unsigned short l = 0;
l < 5;
l++) {
335 for (
unsigned short k = 0;
k < 5;
k++)
340 TCanvas*
c1 =
new TCanvas(
"c1",
Id::str(calo) +
" layer " + Form(
"%d",
layer), 1600, 500);
342 c1->cd(1); hCor->Draw();
343 TLatex
l;
l.SetNDC();
344 TString corrStr = Form(
"Correlation : %.1f%%", 100*hCor->GetCorrelationFactor());
345 l.DrawLatex(0.2, 0.8, corrStr);
346 c1->cd(2); hAll[0]->Draw(
"COL");
347 c1->cd(3); hAll[1]->Draw(
"COL");
355 double lo,
double hi,
const TString&
fileName)
357 unsigned int nBins = 100;
359 TString globalTitle = TString(xip ?
"Derivative" :
"Offset") +
" correction";
360 for (
unsigned short l = 0;
l < 5;
l++)
361 for (
unsigned short k = 0;
k < 5;
k++)
362 h[
l][
k] =
new TH1D(Form(
"h%d%d",
l,
k),
"",
nBins,
l < 2 ? -15 : lo,
l < 2 ? 15 : hi);
367 hAll[0] =
new TH2D(
"h0All",
"#xi_{high}/#delta#xi vs. sampling", 5, 0, 5,
nBins, -15, 15);
368 hAll[1] =
new TH2D(
"hDAll",
"(#xi_{high} - #xi_{low})/#delta#xi vs. sampling", 5, 0, 5,
nBins, -15, 15);
369 hAll[2] =
new TH2D(
"hEAll", globalTitle +
" higher-gain value vs. sampling", 5, 0, 5,
nBins, 0, hi);
370 hAll[3] =
new TH2D(
"hRAll", globalTitle +
" lower-gain value vs. sampling", 5, 0, 5,
nBins, lo, hi);
371 hAll[4] =
new TH2D(
"hCAll", globalTitle +
" higher-gain error vs. sampling", 5, 0, 5,
nBins, lo, hi);
372 TH2D* hCor =
new TH2D(
"hCor", globalTitle +
" lower-gain vs. higher-gain values",
nBins, lo, hi,
nBins, lo, hi);
373 hCor->GetXaxis()->SetTitle(
"#xi_{high}");
374 hCor->GetYaxis()->SetTitle(
"#xi_{low}");
375 for (
unsigned short l = 0;
l < 4;
l++) {
376 hAll[
l]->GetXaxis()->SetTitle(
"sampling");
377 hAll[
l]->GetYaxis()->SetTitle(
"significance");
384 if ((
i+1) % 10000 == 0) cout <<
"Cell # " <<
i+1 <<
"/" <<
nChannels() << endl;
386 if (!history || history->
nData() == 0)
continue;
387 for (
unsigned short k = 0;
k < 5;
k++) {
390 if (!hgSED || !lgSED) {
391 if (hgSED)
delete hgSED;
392 if (lgSED)
delete lgSED;
395 double hgVal = (xip ? hgSED->
xi()(
k) : hgSED->
xip()(
k));
396 double lgVal = (xip ? lgSED->
xi()(
k) : lgSED->
xip()(
k));
397 double hgErr = TMath::Sqrt(xip ? hgSED->
xiErr()(
k,
k) : hgSED->
xipErr()(
k,
k));
401 h[0][
k]->Fill(hgVal/hgErr);
402 h[1][
k]->Fill((hgVal - lgVal)/
error);
403 h[2][
k]->Fill(hgVal);
404 h[3][
k]->Fill(lgVal);
405 h[4][
k]->Fill(hgErr);
406 hCor->Fill(hgVal, lgVal);
409 for (
unsigned int i = 1;
i <=
nBins;
i++)
410 for (
unsigned short l = 0;
l < 5;
l++)
411 for (
unsigned short k = 0;
k < 5;
k++)
412 hAll[
l]->SetBinContent(
k+1,
i,
h[
l][
k]->GetBinContent(
i));
413 TFile*
ff = TFile::Open(
fileName +
".root",
"RECREATE");
414 for (
unsigned short l = 0;
l < 5;
l++) {
415 for (
unsigned short k = 0;
k < 5;
k++)
420 TCanvas*
c1 =
new TCanvas(
"c1",
Id::str(calo) +
" layer " + Form(
"%d",
layer), 1600, 500);
422 c1->cd(1); hCor->Draw();
423 TLatex
l;
l.SetNDC();
424 TString corrStr = Form(
"Correlation : %.1f%%", 100*hCor->GetCorrelationFactor());
425 l.DrawLatex(0.2, 0.8, corrStr);
426 c1->cd(2); hAll[0]->Draw(
"COL");
427 c1->cd(3); hAll[1]->Draw(
"COL");
436 TH1D*
h =
new TH1D(
name,
"Shape error",
nBins, rMin, rMax);
437 h->GetXaxis()->SetTitle(
"#epsilon^{2}");
441 if (!history)
continue;
442 for (
unsigned int j = 0; j < history->
nData(); j++) {
446 h->Fill((delta*delta -
ofc->Gamma()(
k,
k))/(adcMax*adcMax));
449 cout << sqrt(
h->GetMean())*100 <<
" +/- " << sqrt(
h->GetMeanError())*100 << endl;
461 if (!history)
continue;
462 for (
unsigned int j = 0; j < history->
nData(); j++) {
463 double residual1 = history->
data(j)->
delta(k1);
464 double residual2 = history->
data(j)->
delta(k2);
465 sum += residual1*residual2;
476 CovMatrix sum2N(lwb, upb), sumCN(lwb, upb);
477 TVectorD
sum(lwb, upb), sumN(lwb, upb);
478 double sumA = 0,
n = 0;
481 if (!history)
continue;
482 for (
unsigned int j = 0; j < history->
nData(); j++) {
485 if (!
ofc ||
ofc->hasSameRange(lwb, upb)) { cout <<
"Invalid index bounds!" << endl;
return false; }
486 for (
int i1 = lwb; i1 <= upb; i1++) {
488 sumN(i1) = sumN(i1) + history->
data(j)->
delta(i1)/TMath::Power(history->
data(j)->
adcMax(), 2);
489 for (
int i2 = lwb; i2 <= upb; i2++) {
490 sumCN(i1, i2) = sumCN(i1, i2) +
ofc->Gamma()(i1, i2)/TMath::Power(history->
data(j)->
adcMax(), 2);
494 sumA += 1/TMath::Power(history->
data(j)->
adcMax(), 2);
500 k.ResizeTo(lwb, upb);
501 means.ResizeTo(lwb, upb);
503 for (
int i1 = lwb; i1 <= upb; i1++) {
505 sumN(i1) = sumN(i1)/
n;
509 for (
int i1 = lwb; i1 <= upb; i1++) {
510 for (
int i2 = lwb; i2 <= upb; i2++) {
511 sumCN(i1, i2) = sumCN(i1, i2)/
n;
512 sum2N(i1, i2) = sum2N(i1, i2)/
n;
513 k(i1, i2) = sum2N(i1, i2) - sumN(i1)*
sum(i2) -
sum(i1)*sumN(i2)
514 +
sum(i1)*
sum(i2)*sumA - sumCN(i1,i2);
515 k(i1, i2) = (
k(i1, i2) >= 0 ? sqrt(
k(i1, i2)) : -sqrt(-
k(i1, i2)));
527 pshape = pref =
nullptr;
529 std::unique_ptr<SimpleShape> shape;
530 std::unique_ptr<SimpleShape>
ref;
534 if ((
i+1) % 50000 == 0) cout <<
"Cell # " <<
i+1 <<
"/" <<
nChannels() << endl;
536 if (!history)
continue;
537 for (
unsigned int j = 0; j < history->
nData(); j++) {
538 if (!
f.passEvent(*history->
data(j)))
continue;
540 cout <<
"Adding pulse (" <<
n+1 <<
") at hash " <<
i <<
", index " << j <<
", max = " << maxSum/(
n+1) << endl;
542 auto thisData = std::make_unique<SimpleShape>(*history->
data(j));
543 auto thisRef = std::unique_ptr<SimpleShape>(history->
referenceShape(j));
558 pshape = shape.release();
559 pref =
ref.release();
566 if (!history)
return nullptr;
569 if (
adjust)
delete history;
575 unsigned int minSize,
bool weigh,
bool adjust,
bool zeroTime)
const
578 std::vector< std::vector<ResidualCalculator> > ringCalcs;
582 cout <<
"Processing cells" << endl;
585 if ((
i+1) % 10000 == 0) cout <<
"Processing hash = " << (
i+1) << endl;
588 if (residuals && (resTrunc > 0 || timeTrunc > 0)) {
591 residuals = truncated;
594 if (residuals && residuals->
size() < minSize) {
618 cout <<
"Processing phi-symmetric data" << endl;
622 if (ringCalcs[
g][
i].
size() > 0)
623 cout <<
"Saving data for gain " <<
g <<
", phi ring " <<
i <<
", n = " << ringCalcs[
g][
i].
size() <<
", range = " << ringCalcs[
g][
i].rangeStr() << endl;
628 cout <<
"Done!" << endl;