Belle II Software  release-08-01-10
DQMCommonUtils.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 #include <dqm/utilities/DQMCommonUtils.h>
10 
11 #include <framework/database/DBImportObjPtr.h>
12 #include <framework/database/IntervalOfValidity.h>
13 #include <framework/database/DBObjPtr.h>
14 
15 #include <TVectorT.h>
16 
17 using namespace std;
18 using namespace Belle2;
19 
20 
21 int DQMCommonUtils::SetFlag(int Type, int bin, const double* pars, double ratio, TH1F* hist, TH1F* refhist,
22  TH1I* flaghist)
23 {
24  int iret = 0;
25  float WarningLevel = 6.0;
26  float ErrorLevel = 10.0;
27  auto temp = std::unique_ptr<TH1F>(new TH1F("temp", "temp", hist->GetNbinsX(), hist->GetXaxis()->GetXmin(),
28  hist->GetXaxis()->GetXmax()));
29  double NEvents = 0;
30  double flagInt = 0;
31  double flagrInt = 0;
32  for (int j = 0; j < hist->GetNbinsX(); j++) {
33  double val = hist->GetBinContent(j + 1);
34  NEvents += val;
35  val = val / ratio;
36  temp->SetBinContent(j + 1, val);
37  flagInt += temp->GetBinContent(j + 1);
38  flagrInt += refhist->GetBinContent(j + 1);
39  }
40  if (NEvents < 100) { // not enough information for comparition
41  iret = -1;
42  flaghist->SetBinContent(bin + 1, -1);
43  return iret;
44  }
45  double flag = temp->GetMean();
46  double flagErr = temp->GetMeanError();
47  double flagRMS = temp->GetRMS();
48  double flagRMSErr = temp->GetRMSError();
49  double flagr = refhist->GetMean();
50  double flagrErr = refhist->GetMeanError();
51  double flagrRMS = refhist->GetRMS();
52  double flagrRMSErr = refhist->GetRMSError();
53  TString strDebugInfo = Form("Conditions for Flag--->\n source %f %f+-%f %f+-%f\n referen %f %f+-%f %f+-%f\n",
54  flagInt, flag, flagErr, flagRMS, flagRMSErr,
55  flagrInt, flagr, flagrErr, flagrRMS, flagrRMSErr
56  );
57  B2DEBUG(130, strDebugInfo.Data());
58  if (Type == 1) { // counts+mean+RMS use
59  if ((fabs(flag - flagr) > ErrorLevel * (flagErr + flagrErr)) ||
60  (fabs(flagRMS - flagrRMS) > ErrorLevel * (flagRMSErr + flagrRMSErr)) ||
61  (fabs(flagInt - flagrInt) > ErrorLevel * (sqrt(flagInt) + sqrt(flagrInt)))
62  ) {
63  flaghist->SetBinContent(bin + 1, 2);
64  } else if ((fabs(flag - flagr) > WarningLevel * (flagErr + flagrErr)) ||
65  (fabs(flagRMS - flagrRMS) > WarningLevel * (flagRMSErr + flagrRMSErr)) ||
66  (fabs(flagInt - flagrInt) > WarningLevel * (sqrt(flagInt) + sqrt(flagrInt)))
67  ) {
68  flaghist->SetBinContent(bin + 1, 1);
69  } else {
70  flaghist->SetBinContent(bin + 1, 0);
71  }
72  iret = 1;
73  } else if (Type == 2) { // counts use
74  if (fabs(flagInt - flagrInt) > ErrorLevel * (sqrt(flagInt) + sqrt(flagrInt))) {
75  flaghist->SetBinContent(bin + 1, 2);
76  } else if (fabs(flagInt - flagrInt) > WarningLevel * (sqrt(flagInt) + sqrt(flagrInt))) {
77  flaghist->SetBinContent(bin + 1, 1);
78  } else {
79  flaghist->SetBinContent(bin + 1, 0);
80  }
81  iret = 1;
82  } else if (Type == 3) { // mean use
83  if (fabs(flag - flagr) > ErrorLevel * (flagErr + flagrErr)) {
84  flaghist->SetBinContent(bin + 1, 2);
85  } else if (fabs(flag - flagr) > WarningLevel * (flagErr + flagrErr)) {
86  flaghist->SetBinContent(bin + 1, 1);
87  } else {
88  flaghist->SetBinContent(bin + 1, 0);
89  }
90  iret = 1;
91  } else if (Type == 4) { // RMS use
92  if (fabs(flagRMS - flagrRMS) > ErrorLevel * (flagRMSErr + flagrRMSErr)) {
93  flaghist->SetBinContent(bin + 1, 2);
94  } else if (fabs(flagRMS - flagrRMS) > WarningLevel * (flagRMSErr + flagrRMSErr)) {
95  flaghist->SetBinContent(bin + 1, 1);
96  } else {
97  flaghist->SetBinContent(bin + 1, 0);
98  }
99  iret = 1;
100  } else if (Type == 5) { // counts+mean use
101  if ((fabs(flag - flagr) > ErrorLevel * (flagErr + flagrErr)) ||
102  (fabs(flagInt - flagrInt) > ErrorLevel * (sqrt(flagInt) + sqrt(flagrInt)))
103  ) {
104  flaghist->SetBinContent(bin + 1, 2);
105  } else if ((fabs(flag - flagr) > WarningLevel * (flagErr + flagrErr)) ||
106  (fabs(flagInt - flagrInt) > WarningLevel * (sqrt(flagInt) + sqrt(flagrInt)))
107  ) {
108  flaghist->SetBinContent(bin + 1, 1);
109  } else {
110  flaghist->SetBinContent(bin + 1, 0);
111  }
112  iret = 1;
113  } else if (Type == 9) { // bin content use
114  flagInt = temp->GetBinContent(bin + 1);
115  flagrInt = refhist->GetBinContent(bin + 1);
116  if (fabs(flagInt - flagrInt) > ErrorLevel * (sqrt(flagInt) + sqrt(flagrInt))) {
117  flaghist->SetBinContent(bin + 1, 2);
118  } else if (fabs(flagInt - flagrInt) > WarningLevel * (sqrt(flagInt) + sqrt(flagrInt))) {
119  flaghist->SetBinContent(bin + 1, 1);
120  } else {
121  flaghist->SetBinContent(bin + 1, 0);
122  }
123  iret = 1;
124  } else if (Type == 10) {
125  float flag2 = refhist->Chi2Test(temp.get());
126  flaghist->SetBinContent(bin + 1, 0);
127  if (flag2 > pars[1])
128  flaghist->SetBinContent(bin + 1, 2);
129  if (flag2 > pars[0])
130  flaghist->SetBinContent(bin + 1, 1);
131  iret = 1;
132  } else if (Type == 100) {
133  flaghist->SetBinContent(bin + 1, 0);
134  iret = 1;
135  } else {
136  flaghist->SetBinContent(bin + 1, -3);
137  iret = -1;
138  }
139  strDebugInfo = Form("SetFlag---> %f, type %i\n", flaghist->GetBinContent(bin + 1), Type);
140  B2DEBUG(130, strDebugInfo.Data());
141  return iret;
142 }
143 
144 
145 int DQMCommonUtils::SetFlag(int Type, int bin, const double* pars, double ratio, TH1I* hist, TH1I* refhist,
146  TH1I* flaghist)
147 {
148  auto histF = std::unique_ptr<TH1F>(new TH1F("histF", "histF", hist->GetNbinsX(), hist->GetXaxis()->GetXmin(),
149  hist->GetXaxis()->GetXmax()));
150  auto refhistF = std::unique_ptr<TH1F>(new TH1F("refhistF", "refhistF", refhist->GetNbinsX(), refhist->GetXaxis()->GetXmin(),
151  refhist->GetXaxis()->GetXmax()));
152  for (int j = 0; j < hist->GetNbinsX(); j++) {
153  histF->SetBinContent(j + 1, hist->GetBinContent(j + 1));
154  refhistF->SetBinContent(j + 1, refhist->GetBinContent(j + 1));
155  }
156  int ret = SetFlag(Type, bin, pars, ratio, histF.get(), refhistF.get(), flaghist);
157  return ret;
158 }
159 
160 
161 void DQMCommonUtils::CreateDBHisto(TH1F* HistoDB)
162 {
163  IntervalOfValidity iov(0, 0, -1, -1);
164  TString Name = Form("%s_Ref", HistoDB->GetName());
165  DBImportObjPtr<TVectorD> DBHisto(Name.Data());
166  DBHisto.construct(HistoDB->GetNbinsX() + 3);
167  double* Content = new double[HistoDB->GetNbinsX() + 3];
168  Content[0] = HistoDB->GetNbinsX();
169  Content[1] = HistoDB->GetXaxis()->GetXmin();
170  Content[2] = HistoDB->GetXaxis()->GetXmax();
171  for (int i = 0; i < HistoDB->GetNbinsX(); i++) {
172  Content[i + 3] = HistoDB->GetBinContent(i + 1);
173  }
174  DBHisto->SetElements(Content);
175  DBHisto.import(iov);
176  delete [] Content;
177 }
178 
179 
180 void DQMCommonUtils::CreateDBHisto(TH1I* HistoDB)
181 {
182  IntervalOfValidity iov(0, 0, -1, -1);
183  TString Name = Form("%s_Ref", HistoDB->GetName());
184  DBImportObjPtr<TVectorD> DBHisto(Name.Data());
185  DBHisto.construct(HistoDB->GetNbinsX() + 3);
186  double* Content = new double[HistoDB->GetNbinsX() + 3];
187  Content[0] = HistoDB->GetNbinsX();
188  Content[1] = HistoDB->GetXaxis()->GetXmin();
189  Content[2] = HistoDB->GetXaxis()->GetXmax();
190  for (int i = 0; i < HistoDB->GetNbinsX(); i++) {
191  Content[i + 3] = HistoDB->GetBinContent(i + 1);
192  }
193  DBHisto->SetElements(Content);
194  DBHisto.import(iov);
195  delete [] Content;
196 }
197 
198 
199 void DQMCommonUtils::CreateDBHistoGroup(TH1F** HistoDB, int number)
200 {
201  IntervalOfValidity iov(0, 0, -1, -1);
202  TString Name = Form("%s_Ref", HistoDB[0]->GetName());
203  DBImportObjPtr<TVectorD> DBHisto(Name.Data());
204  DBHisto.construct(number * HistoDB[0]->GetNbinsX() + 3);
205  double* Content = new double[number * HistoDB[0]->GetNbinsX() + 3];
206  Content[0] = HistoDB[0]->GetNbinsX();
207  Content[1] = HistoDB[0]->GetXaxis()->GetXmin();
208  Content[2] = HistoDB[0]->GetXaxis()->GetXmax();
209  for (int j = 0; j < number; j++) {
210  for (int i = 0; i < HistoDB[j]->GetNbinsX(); i++) {
211  Content[j * HistoDB[j]->GetNbinsX() + i + 3] = HistoDB[j]->GetBinContent(i + 1);
212  }
213  }
214  DBHisto->SetElements(Content);
215  DBHisto.import(iov);
216  delete [] Content;
217 }
218 
219 
220 void DQMCommonUtils::CreateDBHistoGroup(TH1I** HistoDB, int number)
221 {
222  IntervalOfValidity iov(0, 0, -1, -1);
223  TString Name = Form("%s_Ref", HistoDB[0]->GetName());
224  DBImportObjPtr<TVectorD> DBHisto(Name.Data());
225  DBHisto.construct(number * HistoDB[0]->GetNbinsX() + 3);
226  double* Content = new double[number * HistoDB[0]->GetNbinsX() + 3];
227  Content[0] = HistoDB[0]->GetNbinsX();
228  Content[1] = HistoDB[0]->GetXaxis()->GetXmin();
229  Content[2] = HistoDB[0]->GetXaxis()->GetXmax();
230  for (int j = 0; j < number; j++) {
231  for (int i = 0; i < HistoDB[j]->GetNbinsX(); i++) {
232  Content[j * HistoDB[j]->GetNbinsX() + i + 3] = HistoDB[j]->GetBinContent(i + 1);
233  }
234  }
235  DBHisto->SetElements(Content);
236  DBHisto.import(iov);
237  delete [] Content;
238 }
239 
240 
241 int DQMCommonUtils::LoadDBHisto(TH1F* HistoDB)
242 {
243  TString Name = Form("%s_Ref", HistoDB->GetName());
244  DBObjPtr<TVectorD> DBHisto(Name.Data());
245  int ret = 0;
246  if (DBHisto.isValid()) {
247  ret = 1;
248  if (HistoDB->GetNbinsX() != (int)DBHisto->GetMatrixArray()[0]) ret = 0;
249  if (HistoDB->GetXaxis()->GetXmin() != DBHisto->GetMatrixArray()[1]) ret = 0;
250  if (HistoDB->GetXaxis()->GetXmax() != DBHisto->GetMatrixArray()[2]) ret = 0;
251  if (ret) {
252  for (int i = 0; i < HistoDB->GetNbinsX(); i++) {
253  HistoDB->SetBinContent(i + 1, (int)DBHisto->GetMatrixArray()[i + 3]);
254  }
255  }
256  }
257  if (!ret) {
258  B2INFO("ERROR to open reference histogram: " << Name.Data());
259  }
260  return ret;
261 }
262 
263 
264 int DQMCommonUtils::LoadDBHisto(TH1I* HistoDB)
265 {
266  TString Name = Form("%s_Ref", HistoDB->GetName());
267  DBObjPtr<TVectorD> DBHisto(Name.Data());
268  int ret = 0;
269  if (DBHisto.isValid()) {
270  ret = 1;
271  if (HistoDB->GetNbinsX() != (int)DBHisto->GetMatrixArray()[0]) ret = 0;
272  if (HistoDB->GetXaxis()->GetXmin() != DBHisto->GetMatrixArray()[1]) ret = 0;
273  if (HistoDB->GetXaxis()->GetXmax() != DBHisto->GetMatrixArray()[2]) ret = 0;
274  if (ret) {
275  for (int i = 0; i < HistoDB->GetNbinsX(); i++) {
276  HistoDB->SetBinContent(i + 1, (int)DBHisto->GetMatrixArray()[i + 3]);
277  }
278  }
279  }
280  if (!ret) {
281  B2INFO("ERROR to open reference histogram: " << Name.Data());
282  }
283  return ret;
284 }
285 
286 
287 int DQMCommonUtils::LoadDBHistoGroup(TH1F** HistoDB, int number)
288 {
289  TString Name = Form("%s_Ref", HistoDB[0]->GetName());
290  DBObjPtr<TVectorD> DBHisto(Name.Data());
291  int ret = 0;
292  if (DBHisto.isValid()) {
293  ret = 1;
294  if (HistoDB[0]->GetNbinsX() != (int)DBHisto->GetMatrixArray()[0]) ret = 0;
295  if (HistoDB[0]->GetXaxis()->GetXmin() != DBHisto->GetMatrixArray()[1]) ret = 0;
296  if (HistoDB[0]->GetXaxis()->GetXmax() != DBHisto->GetMatrixArray()[2]) ret = 0;
297  for (int j = 0; j < number; j++) {
298  for (int i = 0; i < HistoDB[j]->GetNbinsX(); i++) {
299  HistoDB[j]->SetBinContent(i + 1, DBHisto->GetMatrixArray()[j * HistoDB[j]->GetNbinsX() + i + 3]);
300  }
301  }
302  }
303  if (!ret) {
304  B2INFO("ERROR to open reference histogram: " << Name.Data());
305  }
306  return ret;
307 }
308 
309 
310 int DQMCommonUtils::LoadDBHistoGroup(TH1I** HistoDB, int number)
311 {
312  TString Name = Form("%s_Ref", HistoDB[0]->GetName());
313  DBObjPtr<TVectorD> DBHisto(Name.Data());
314  int ret = 0;
315  if (DBHisto.isValid()) {
316  ret = 1;
317  if (HistoDB[0]->GetNbinsX() != (int)DBHisto->GetMatrixArray()[0]) ret = 0;
318  if (HistoDB[0]->GetXaxis()->GetXmin() != DBHisto->GetMatrixArray()[1]) ret = 0;
319  if (HistoDB[0]->GetXaxis()->GetXmax() != DBHisto->GetMatrixArray()[2]) ret = 0;
320  for (int j = 0; j < number; j++) {
321  for (int i = 0; i < HistoDB[j]->GetNbinsX(); i++) {
322  HistoDB[j]->SetBinContent(i + 1, DBHisto->GetMatrixArray()[j * HistoDB[j]->GetNbinsX() + i + 3]);
323  }
324  }
325  }
326  if (!ret) {
327  B2INFO("ERROR to open reference histogram: " << Name.Data());
328  }
329  return ret;
330 }
331 
Belle2::DBAccessorBase::isValid
bool isValid() const
Check whether a valid object was obtained from the database.
Definition: DBAccessorBase.h:72
Belle2::DBImportBase::import
bool import(const IntervalOfValidity &iov)
Import the object to database.
Definition: DBImportBase.cc:36
Belle2::DBImportObjPtr
Class for importing a single object to the database.
Definition: DBImportObjPtr.h:23
Belle2::DBImportObjPtr::construct
void construct(Args &&... params)
Construct an object of type T in this DBImportObjPtr using the provided constructor arguments.
Definition: DBImportObjPtr.h:47
Belle2::DBObjPtr
Class for accessing objects in the database.
Definition: DBObjPtr.h:21
Belle2::IntervalOfValidity
A class that describes the interval of experiments/runs for which an object in the database is valid.
Definition: IntervalOfValidity.h:25
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17