ATLAS Offline Software
TFCS1DFunction.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 
7 #include "TH1.h"
8 #include "TCanvas.h"
9 #include "TRandom.h"
10 #include <string>
11 #include <iostream>
12 #include <cmath>
13 
14 //=============================================
15 //======= TFCS1DFunction =========
16 //=============================================
17 
18 void TFCS1DFunction::rnd_to_fct(float value[], const float rnd[]) const {
19  value[0] = rnd_to_fct(rnd[0]);
20 }
21 
22 double TFCS1DFunction::get_maxdev(TH1 *h_input1, TH1 *h_approx1) {
23  TH1D *h_input = (TH1D *)h_input1->Clone("h_input");
24  TH1D *h_approx = (TH1D *)h_approx1->Clone("h_approx");
25 
26  double maxdev = 0.0;
27 
28  // normalize the histos to the same area:
29  double integral_input = h_input->Integral();
30  double integral_approx = 0.0;
31  for (int b = 1; b <= h_input->GetNbinsX(); b++)
32  integral_approx +=
33  h_approx->GetBinContent(h_approx->FindBin(h_input->GetBinCenter(b)));
34  h_approx->Scale(integral_input / integral_approx);
35 
36  double ymax = h_approx->GetBinContent(h_approx->GetNbinsX()) -
37  h_approx->GetBinContent(h_approx->GetMinimumBin());
38  for (int i = 1; i <= h_input->GetNbinsX(); i++) {
39  double val = std::abs(h_approx->GetBinContent(
40  h_approx->FindBin(h_input->GetBinCenter(i))) -
41  h_input->GetBinContent(i)) /
42  ymax;
43  if (val > maxdev)
44  maxdev = val;
45  }
46 
47  delete h_input;
48  delete h_approx;
49 
50  return maxdev * 100.0;
51 }
52 
54  const TH1 *hist, std::vector<double> &integral_vec, int &first, int &last) {
56  Int_t nbins = hist->GetNbinsX();
57 
58  float integral = 0;
59  integral_vec.resize(nbins);
60  for (int ix = 1; ix <= nbins; ix++) {
61  float binval = hist->GetBinContent(ix);
62  if (binval < 0) {
63  // Can't work if a bin is negative, forcing bins to 0 in this case
64  double fraction = binval / hist->Integral();
65  if (std::abs(fraction) > 1e-5) {
66  ATH_MSG_NOCLASS(logger, "Warning: bin content is negative in histogram "
67  << hist->GetName() << " : "
68  << hist->GetTitle() << " binval=" << binval
69  << " " << fraction * 100
70  << "% of integral=" << hist->Integral()
71  << ". Forcing bin to 0.");
72  }
73  binval = 0;
74  }
75  integral += binval;
76  integral_vec[ix - 1] = integral;
77  }
78 
79  for (first = 0; first < nbins; first++)
80  if (integral_vec[first] != 0)
81  break;
82  for (last = nbins - 1; last > 0; last--)
83  if (integral_vec[last] != integral)
84  break;
85  last++;
86 
87  if (integral <= 0) {
88  ATH_MSG_NOCLASS(logger, "Error: histogram "
89  << hist->GetName() << " : " << hist->GetTitle()
90  << " integral=" << integral << " is <=0");
91  }
92  return integral;
93 }
94 
95 TH1 *TFCS1DFunction::generate_histogram_random_slope
96 ATLAS_NOT_THREAD_SAFE(int nbinsx, double xmin, double xmax,
97  double zerothreshold) {
98  TH1 *hist = new TH1D("test_slope1D", "test_slope1D", nbinsx, xmin, xmax);
99  hist->Sumw2();
100  for (int ix = 1; ix <= nbinsx; ++ix) {
101  double val = (0.5 + gRandom->Rndm()) * (nbinsx + ix);
102  if (gRandom->Rndm() < zerothreshold)
103  val = 0;
104  hist->SetBinContent(ix, val);
105  hist->SetBinError(ix, 0);
106  }
107  return hist;
108 }
109 
110 TH1 *TFCS1DFunction::generate_histogram_random_gauss
111 ATLAS_NOT_THREAD_SAFE(int nbinsx, int ntoy, double xmin, double xmax,
112  double xpeak, double sigma) {
113  TH1 *hist = new TH1D("test_gauss1D", "test_gauss1D", nbinsx, xmin, xmax);
114  hist->Sumw2();
115  for (int i = 1; i <= ntoy; ++i) {
116  double x = gRandom->Gaus(xpeak, sigma);
117  if (x >= xmin && x < xmax)
118  hist->Fill(x);
119  }
120  return hist;
121 }
122 
123 void TFCS1DFunction::unit_test ATLAS_NOT_THREAD_SAFE(TH1 *hist,
124  TFCS1DFunction *rtof,
125  int nrnd, TH1 *histfine) {
127  ATH_MSG_NOCLASS(logger, "========= " << hist->GetName() << " funcsize="
128  << rtof->MemorySize() << " ========");
129  int nbinsx = hist->GetNbinsX();
130  double integral = hist->Integral();
131 
132  float value[2];
133  float rnd[2];
134  //cppcheck-suppress uninitvar
135  for (rnd[0] = 0; rnd[0] < 0.9999; rnd[0] += 0.25) {
136  rtof->rnd_to_fct(value, rnd);
137  ATH_MSG_NOCLASS(logger, "rnd0=" << rnd[0] << " -> x=" << value[0]);
138  }
139 
140  TH1 *hist_val = nullptr;
141  if (histfine)
142  hist_val = (TH1 *)histfine->Clone(TString(hist->GetName()) + "hist_val");
143  else
144  hist_val = (TH1 *)hist->Clone(TString(hist->GetName()) + "hist_val");
145  double weightfine = hist_val->Integral() / nrnd;
146  hist_val->SetTitle("toy simulation");
147  hist_val->Reset();
148  hist_val->SetLineColor(2);
149  hist_val->Sumw2();
150 
151  TH1 *hist_diff = (TH1 *)hist->Clone(TString(hist->GetName()) + "_difference");
152  hist_diff->SetTitle("cut efficiency difference");
153  hist_diff->Reset();
154  hist_diff->Sumw2();
155 
156  double weight = integral / nrnd;
157  for (int i = 0; i < nrnd; ++i) {
158  rnd[0] = gRandom->Rndm();
159  rtof->rnd_to_fct(value, rnd);
160  hist_val->Fill(value[0], weightfine);
161  hist_diff->Fill(value[0], weight);
162  }
163  hist_diff->Add(hist, -1);
164  hist_diff->Scale(1.0 / integral);
165 
166  TH1F *hist_pull =
167  new TH1F(TString(hist->GetName()) + "_pull", "pull", 200, -10, 10);
168  for (int ix = 1; ix <= nbinsx; ++ix) {
169  float val = hist_diff->GetBinContent(ix);
170  float err = hist_diff->GetBinError(ix);
171  if (err > 0)
172  hist_pull->Fill(val / err);
173  }
174 
175 // Screen output in athena won't make sense and would require linking of
176 // additional libraries
177 #if defined(__FastCaloSimStandAlone__)
178  TCanvas *c = new TCanvas(hist->GetName(), hist->GetName());
179  c->Divide(2, 2);
180 
181  c->cd(1);
182  if (histfine) {
183  histfine->SetLineColor(kGray);
184  histfine->DrawClone("hist");
185  hist->DrawClone("same");
186  } else {
187  hist->DrawClone();
188  }
189  hist_val->DrawClone("sameshist");
190 
191  c->cd(2);
192  if (histfine) {
193  histfine->SetLineColor(kGray);
194  histfine->DrawClone("hist");
195  hist->DrawClone("same");
196  } else {
197  hist->DrawClone();
198  }
199  hist_val->DrawClone("sameshist");
200  gPad->SetLogy();
201 
202  c->cd(3);
203  hist_diff->Draw();
204 
205  c->cd(4);
206  hist_pull->Draw();
207 
208  c->SaveAs(".png");
209 #endif
210 }
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
ISF_FCS::MLogging
Cut down AthMessaging.
Definition: MLogging.h:176
pdg_comparison.sigma
sigma
Definition: pdg_comparison.py:324
integral
double integral(TH1 *h)
Definition: computils.cxx:59
plotmaker.hist
hist
Definition: plotmaker.py:148
TFCS1DFunction::rnd_to_fct
virtual void rnd_to_fct(float value[], const float rnd[]) const
Function gets array of random numbers rnd[] in the range [0,1) as arguments and returns function valu...
Definition: TFCS1DFunction.cxx:18
athena.value
value
Definition: athena.py:124
x
#define x
dqt_zlumi_pandas.weight
int weight
Definition: dqt_zlumi_pandas.py:190
ATLAS_NOT_THREAD_SAFE
TH1 *TFCS1DFunction::generate_histogram_random_slope ATLAS_NOT_THREAD_SAFE(int nbinsx, double xmin, double xmax, double zerothreshold)
Definition: TFCS1DFunction.cxx:96
TFCS1DFunction::get_maxdev
static double get_maxdev(TH1 *, TH1 *)
Definition: TFCS1DFunction.cxx:22
dqt_zlumi_pandas.err
err
Definition: dqt_zlumi_pandas.py:183
lumiFormat.i
int i
Definition: lumiFormat.py:85
xmin
double xmin
Definition: listroot.cxx:60
TFCS1DFunction::CheckAndIntegrate1DHistogram
static double CheckAndIntegrate1DHistogram(const TH1 *hist, std::vector< double > &integral_vec, int &first, int &last)
Definition: TFCS1DFunction.cxx:53
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:76
SCT_CalibAlgs::nbins
@ nbins
Definition: SCT_CalibNumbers.h:10
ATH_MSG_NOCLASS
#define ATH_MSG_NOCLASS(logger_name, x)
Definition: MLogging.h:52
Pythia8_RapidityOrderMPI.val
val
Definition: Pythia8_RapidityOrderMPI.py:14
DeMoScan.first
bool first
Definition: DeMoScan.py:534
xmax
double xmax
Definition: listroot.cxx:61
TFCS1DFunction.h
python.TrigEgammaMonitorHelper.TH1F
def TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:24
TFCS1DFunction
Definition: TFCS1DFunction.h:17
python.compressB64.c
def c
Definition: compressB64.py:93
python.iconfTool.gui.pad.logger
logger
Definition: pad.py:14
TFCSFunction::MemorySize
virtual std::size_t MemorySize() const
Gives the total memory size, including the size of additional memory allocated inside the class.
Definition: TFCSFunction.h:19
ymax
double ymax
Definition: listroot.cxx:64