|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
|
#include <cmath>
#include "TH1D.h"
#include "TProfile2D.h"
TH1D* projectProfile2DAlongX(TProfile2D* prof2d) {
TH1D* res = nullptr;
if (prof2d) {
char name[200];
sprintf(name, "%s_proj", prof2d->GetName());
res = new TH1D(
name, prof2d->GetTitle(), prof2d->GetNbinsY(), prof2d->GetYaxis()->GetXmin(), prof2d->GetYaxis()->GetXmax());
res->SetDirectory(nullptr);
res->Sumw2();
for (int iy = 1; iy < prof2d->GetNbinsY() + 1; ++iy) {
double sum = 0.;
double sumsq = 0.;
double nevt = 0.;
for (int ix = 1; ix < prof2d->GetNbinsX() + 1; ++ix) {
const int ibin = prof2d->GetBin(ix, iy);
sum += prof2d->GetBinContent(ibin) * prof2d->GetBinEntries(ibin);
sumsq += prof2d->GetBinError(ibin) * prof2d->GetBinError(ibin) * prof2d->GetBinEntries(ibin) *
prof2d->GetBinEntries(ibin) +
prof2d->GetBinContent(ibin) * prof2d->GetBinContent(ibin) * prof2d->GetBinEntries(ibin);
nevt += prof2d->GetBinEntries(ibin);
}
double mean = nevt == 0 ? 0 : sum / nevt;
double meansq = nevt == 0 ? 0 : sumsq / nevt;
double err = meansq >= mean * mean ? sqrt(meansq - mean * mean) : 0;
err = nevt == 0 ? 0 : err / sqrt(nevt);
res->SetBinContent(iy, mean);
res->SetBinError(iy, err);
}
}
return res;
}
TH1D* projectProfile2DAlongY(TProfile2D* prof2d) {
TH1D* res = nullptr;
if (prof2d) {
char name[200];
sprintf(name, "%s_proj", prof2d->GetName());
res = new TH1D(
name, prof2d->GetTitle(), prof2d->GetNbinsX(), prof2d->GetXaxis()->GetXmin(), prof2d->GetXaxis()->GetXmax());
res->SetDirectory(nullptr);
res->Sumw2();
for (int ix = 1; ix < prof2d->GetNbinsX() + 1; ++ix) {
double sum = 0.;
double sumsq = 0.;
double nevt = 0.;
for (int iy = 1; iy < prof2d->GetNbinsY() + 1; ++iy) {
const int ibin = prof2d->GetBin(ix, iy);
sum += prof2d->GetBinContent(ibin) * prof2d->GetBinEntries(ibin);
sumsq += prof2d->GetBinError(ibin) * prof2d->GetBinError(ibin) * prof2d->GetBinEntries(ibin) *
prof2d->GetBinEntries(ibin) +
prof2d->GetBinContent(ibin) * prof2d->GetBinContent(ibin) * prof2d->GetBinEntries(ibin);
nevt += prof2d->GetBinEntries(ibin);
}
double mean = nevt == 0 ? 0 : sum / nevt;
double meansq = nevt == 0 ? 0 : sumsq / nevt;
double err = meansq >= mean * mean ? sqrt(meansq - mean * mean) : 0;
err = nevt == 0 ? 0 : err / sqrt(nevt);
res->SetBinContent(ix, mean);
res->SetBinError(ix, err);
}
}
return res;
}
|