ATLAS Offline Software
TFCS1DFunctionInt16Histogram.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 #include <algorithm>
7 #include <iostream>
8 #include "TMath.h"
9 #include "TCanvas.h"
10 #include "TH2F.h"
11 #include "TRandom.h"
12 #include "TFile.h"
13 
14 //=============================================
15 //======= TFCS1DFunctionInt16Histogram =========
16 //=============================================
17 
20 
22  Int_t nbinsx = hist->GetNbinsX();
23  Int_t nbins = nbinsx;
24 
25  float integral = 0;
26  m_HistoBorders.resize(nbinsx + 1);
27  m_HistoContents.resize(nbins);
28  std::vector<double> temp_HistoContents(nbins);
29  int ibin = 0;
30  for (int ix = 1; ix <= nbinsx; ix++) {
31  float binval = hist->GetBinContent(ix);
32  if (binval < 0) {
33  // Can't work if a bin is negative, forcing bins to 0 in this case
34  double fraction = binval / hist->Integral();
35  if (TMath::Abs(fraction) > 1e-5) {
36  ATH_MSG_WARNING("bin content is negative in histogram "
37  << hist->GetName() << " : " << hist->GetTitle()
38  << " binval=" << binval << " " << fraction * 100
39  << "% of integral=" << hist->Integral()
40  << ". Forcing bin to 0.");
41  }
42  binval = 0;
43  }
44  integral += binval;
45  temp_HistoContents[ibin] = integral;
46  ++ibin;
47  }
48  if (integral <= 0) {
49  ATH_MSG_ERROR("histogram " << hist->GetName() << " : " << hist->GetTitle()
50  << " integral=" << integral << " is <=0");
51  m_HistoBorders.resize(0);
52  m_HistoContents.resize(0);
53  return;
54  }
55 
56  for (int ix = 1; ix <= nbinsx; ix++)
57  m_HistoBorders[ix - 1] = hist->GetXaxis()->GetBinLowEdge(ix);
58  m_HistoBorders[nbinsx] = hist->GetXaxis()->GetXmax();
59 
60  for (ibin = 0; ibin < nbins; ++ibin) {
61  m_HistoContents[ibin] = s_MaxValue * (temp_HistoContents[ibin] / integral);
62  // ATH_MSG_INFO("bin="<<ibin<<" val="<<m_HistoContents[ibin]);
63  }
64 }
65 
66 double TFCS1DFunctionInt16Histogram::rnd_to_fct(double rnd) const {
67  if (m_HistoContents.empty()) {
68  return 0;
69  }
70  HistoContent_t int_rnd = s_MaxValue * rnd;
71  auto it =
72  std::upper_bound(m_HistoContents.begin(), m_HistoContents.end(), int_rnd);
73  int ibin = std::distance(m_HistoContents.begin(), it);
74  if (ibin >= (int)m_HistoContents.size())
75  ibin = m_HistoContents.size() - 1;
76  Int_t binx = ibin;
77 
78  HistoContent_t basecont = 0;
79  if (ibin > 0)
80  basecont = m_HistoContents[ibin - 1];
81 
82  HistoContent_t dcont = m_HistoContents[ibin] - basecont;
83  if (dcont > 0) {
84  return m_HistoBorders[binx] +
85  ((m_HistoBorders[binx + 1] - m_HistoBorders[binx]) *
86  (int_rnd - basecont)) /
87  dcont;
88  } else {
89  return m_HistoBorders[binx] +
90  (m_HistoBorders[binx + 1] - m_HistoBorders[binx]) / 2;
91  }
92 }
93 
94 void TFCS1DFunctionInt16Histogram::unit_test ATLAS_NOT_THREAD_SAFE(TH1 *hist) {
96  int nbinsx;
97  if (hist == nullptr) {
98  nbinsx = 50;
99  hist = new TH1D("test1D", "test1D", nbinsx, 0, 1);
100  hist->Sumw2();
101  for (int ix = 1; ix <= nbinsx; ++ix) {
102  double val = (0.5 + gRandom->Rndm()) * (nbinsx + ix);
103  if (gRandom->Rndm() < 0.1)
104  val = 0;
105  hist->SetBinContent(ix, val);
106  hist->SetBinError(ix, 0);
107  }
108  }
110  nbinsx = hist->GetNbinsX();
111 
112  float value[2];
113  float rnd[2];
114  //cppcheck-suppress uninitvar
115  for (rnd[0] = 0; rnd[0] < 0.9999; rnd[0] += 0.25) {
116  rtof.rnd_to_fct(value, rnd);
117  ATH_MSG_NOCLASS(logger, "rnd0=" << rnd[0] << " -> x=" << value[0]);
118  }
119 
120  TH1 *hist_val = (TH1 *)hist->Clone("hist_val");
121  hist_val->SetTitle("difference");
122  hist_val->Reset();
123  int nrnd = 100000;
124  double weight = hist->Integral() / nrnd;
125  hist_val->Sumw2();
126  for (int i = 0; i < nrnd; ++i) {
127  rnd[0] = gRandom->Rndm();
128  rtof.rnd_to_fct(value, rnd);
129  hist_val->Fill(value[0], weight);
130  }
131  hist_val->Add(hist, -1);
132 
133  TH1F *hist_pull = new TH1F("pull", "pull", 200, -10, 10);
134  for (int ix = 1; ix <= nbinsx; ++ix) {
135  float val = hist_val->GetBinContent(ix);
136  float err = hist_val->GetBinError(ix);
137  if (err > 0)
138  hist_pull->Fill(val / err);
139  ATH_MSG_NOCLASS(logger, "val=" << val << " err=" << err);
140  }
141 
142 // Screen output in athena won't make sense and would require linking of
143 // additional libraries
144 #if defined(__FastCaloSimStandAlone__)
145  new TCanvas("input", "Input");
146  hist->Draw();
147 
148  new TCanvas("validation", "Validation");
149  hist_val->Draw();
150 
151  new TCanvas("pull", "Pull");
152  hist_pull->Draw();
153 #endif
154 }
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
TFCS1DFunctionInt16Histogram.h
ISF_FCS::MLogging
Cut down AthMessaging.
Definition: MLogging.h:176
integral
double integral(TH1 *h)
Definition: computils.cxx:58
plotmaker.hist
hist
Definition: plotmaker.py:148
skel.it
it
Definition: skel.GENtoEVGEN.py:396
athena.value
value
Definition: athena.py:124
TFCS1DFunctionInt16Histogram::Initialize
void Initialize(const TH1 *hist)
Definition: TFCS1DFunctionInt16Histogram.cxx:21
TFCS1DFunctionInt16Histogram::m_HistoContents
std::vector< HistoContent_t > m_HistoContents
Definition: TFCS1DFunctionInt16Histogram.h:47
dqt_zlumi_pandas.weight
int weight
Definition: dqt_zlumi_pandas.py:189
TFCS1DFunctionInt16Histogram::m_HistoBorders
std::vector< float > m_HistoBorders
Definition: TFCS1DFunctionInt16Histogram.h:46
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
dqt_zlumi_pandas.err
err
Definition: dqt_zlumi_pandas.py:182
lumiFormat.i
int i
Definition: lumiFormat.py:85
TFCS1DFunctionInt16Histogram::s_MaxValue
static const HistoContent_t s_MaxValue
Definition: TFCS1DFunctionInt16Histogram.h:28
TFCS1DFunctionInt16Histogram
Definition: TFCS1DFunctionInt16Histogram.h:15
TFCS1DFunctionInt16Histogram::rnd_to_fct
virtual double rnd_to_fct(double rnd) const
Function gets random number rnd in the range [0,1) as argument and returns function value according t...
Definition: TFCS1DFunctionInt16Histogram.cxx:66
ATLAS_NOT_THREAD_SAFE
void TFCS1DFunctionInt16Histogram::unit_test ATLAS_NOT_THREAD_SAFE(TH1 *hist)
Definition: TFCS1DFunctionInt16Histogram.cxx:94
SCT_CalibAlgs::nbins
@ nbins
Definition: SCT_CalibNumbers.h:10
ATH_MSG_NOCLASS
#define ATH_MSG_NOCLASS(logger_name, x)
Definition: MLogging.h:52
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
Pythia8_RapidityOrderMPI.val
val
Definition: Pythia8_RapidityOrderMPI.py:14
python.TrigEgammaMonitorHelper.TH1F
def TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:24
TFCS1DFunctionInt16Histogram::HistoContent_t
uint16_t HistoContent_t
Definition: TFCS1DFunctionInt16Histogram.h:27
Amg::distance
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
Definition: GeoPrimitivesHelpers.h:54
python.iconfTool.gui.pad.logger
logger
Definition: pad.py:14