Heavy ion Analysis Libriares
Loading...
Searching...
No Matches
DividedHisto.cxx
1/*
2 * HalDividedHisto1D.cxx
3 *
4 * Created on: 13-10-2014
5 * Author: Daniel Wielanek
6 * E-mail: daniel.wielanek@gmail.com
7 * Warsaw University of Technology, Faculty of Physics
8 */
9
10#include "DividedHisto.h"
11
12#include "Cout.h"
13#include "HistoStyle.h"
14#include "HistogramManager.h"
15#include "HtmlCore.h"
16#include "HtmlFile.h"
17#include "HtmlTable.h"
18#include "Std.h"
19#include "StdHist.h"
20#include "StdString.h"
21#include "Style.h"
22
23#include <TAttFill.h>
24#include <TAttLine.h>
25#include <TAxis.h>
26#include <TBrowser.h>
27#include <TCanvas.h>
28#include <TCollection.h>
29#include <TH1.h>
30#include <TLegend.h>
31#include <TList.h>
32#include <TMath.h>
33#include <TMathBase.h>
34#include <TNamed.h>
35#include <TObjString.h>
36#include <TROOT.h>
37#include <TSystem.h>
38#include <TVirtualPad.h>
39#include <fstream>
40#include <vector>
41
42namespace Hal {
44 fNumIsCloned(kFALSE),
45 fDenIsCloned(kFALSE),
46 fNum(nullptr),
47 fDen(nullptr),
48 fBinomial(kFALSE),
49 fDim(0),
50 fNormMin(nullptr),
51 fNormMax(nullptr),
52 fScale(0.0),
53 fAxisName(""),
54 fComment(""),
55 fLabels(nullptr) {}
56
57 DividedHisto1D::DividedHisto1D(TString name, const Int_t dim) :
58 fNumIsCloned(kFALSE),
59 fDenIsCloned(kFALSE),
60 fNum(NULL),
61 fDen(NULL),
62 fBinomial(kFALSE),
63 fDim(dim),
64 fNormMin(NULL),
65 fNormMax(NULL),
66 fScale(0),
67 fLabels(NULL) {
68 fNormMin = new Double_t[fDim];
69 fNormMax = new Double_t[fDim];
70 fName = name;
71 for (int i = 0; i < fDim; i++) {
72 fNormMin[i] = fNormMax[i] = 0;
73 }
74 fComment = " ";
75 fLabels = new TList();
76 fLabels->SetOwner(kTRUE);
77 fLabels->SetName("Labels");
78 }
79
81 fNumIsCloned(kFALSE),
82 fDenIsCloned(kFALSE),
83 fNum(NULL),
84 fDen(NULL),
85 fBinomial(other.fBinomial),
86 fDim(other.fDim),
87 fScale(other.fScale),
88 fAxisName(other.fAxisName),
89 fComment(other.fComment),
90 fLabels(NULL) {
91 if (other.fLabels) { fLabels = (TList*) other.fLabels->Clone(); }
92 if (other.fNum) {
93 fNum = (TH1*) other.fNum->Clone();
94 fNumIsCloned = kTRUE;
95 }
96 if (other.fDen) {
97 fDen = (TH1*) other.fDen->Clone();
98 fDenIsCloned = kTRUE;
99 }
100 fNormMin = new Double_t[fDim];
101 fNormMax = new Double_t[fDim];
102 for (int i = 0; i < fDim; i++) {
103 fNormMin[i] = other.fNormMin[i];
104 fNormMax[i] = other.fNormMax[i];
105 }
106 }
107
108 void DividedHisto1D::Normalize(TH1* h) const {
109 if (fScale == 0) {
110 for (int i = 0; i < fDim; i++) {
111 Int_t binMinX = fNum->GetXaxis()->FindBin(fNormMin[0]);
112 Int_t binMaxX = fNum->GetXaxis()->FindBin(fNormMax[0]);
113 Double_t num = fNum->Integral(binMinX, binMaxX);
114 binMinX = fDen->GetXaxis()->FindBin(fNormMin[0]);
115 binMaxX = fDen->GetXaxis()->FindBin(fNormMax[0]);
116 Double_t den = fDen->Integral(binMinX, binMaxX);
117 Double_t scale = 1.0;
118 if (num == 0) {
119 Cout::PrintInfo("Numerator integral is zero, set scale to one", Hal::EInfo::kLowWarning);
120 } else {
121 scale = den / num;
122 }
123 if (fNormMin[i] == fNormMax[i]) {
124 Cout::PrintInfo(Form("Same norm min and norm max on %i axis", i), Hal::EInfo::kLowWarning);
125 scale = 1.0;
126 }
127 if (scale == 0.0) scale = 1.0;
128 h->Scale(scale);
129 }
130 } else {
131 h->Scale(fScale);
132 }
133 }
134
135 void DividedHisto1D::AddNum(TH1* num, Option_t* opt) {
136 if (fNumIsCloned && fNum != NULL) delete fNum;
137 if (num == NULL) {
138 Cout::PrintInfo("HalDividedHisto - trying to assign NULL Numerator", Hal::EInfo::kLowWarning);
139 fNum = NULL;
140 return;
141 }
142 TString option = opt;
143 if (option.EqualTo("clone")) {
144 fNumIsCloned = kTRUE;
145 fNum = (TH1*) num->Clone();
146 } else {
147 fNumIsCloned = kFALSE;
148 fNum = num;
149 }
150 }
151
152 void DividedHisto1D::AddDen(TH1* den, Option_t* opt) {
153 if (fDenIsCloned && fDen != NULL) delete fDen;
154 if (den == NULL) {
155 Cout::PrintInfo("HalDividedHisto - trying to assign NULL Denominator", Hal::EInfo::kLowWarning);
156 fDen = NULL;
157 return;
158 }
159 TString option = opt;
160 if (option.EqualTo("clone")) {
161 fDenIsCloned = kTRUE;
162 fDen = (TH1*) den->Clone();
163 } else {
164 fDenIsCloned = kFALSE;
165 fDen = den;
166 }
167 }
168
169 void DividedHisto1D::AddNumDen(TH1* num, TH1* den, Option_t* opt) {
170 AddNum(num, opt);
171 AddDen(den, opt);
172 }
173
174 TH1* DividedHisto1D::GetNum() const { return fNum; }
175
176 TH1* DividedHisto1D::GetDen() const { return fDen; }
177
178 TH1* DividedHisto1D::GetHist(Bool_t normalized) const {
179 TH1* res = (TH1*) fNum->Clone();
180 TString name = res->GetTitle();
181 name = name + "Divided";
182 res->SetTitle(name);
183 if (fBinomial) {
184 res->Divide(res, fDen, 1, 1, "B");
185 } else {
186 res->Divide(fDen);
187 }
188 if (normalized) Normalize(res);
189 SetHistoName(res);
190 return res;
191 }
192
193 void DividedHisto1D::SetNorm(Double_t min, Double_t max, Int_t axis) {
194 if (axis < 0) {
195 Cout::PrintInfo("HalDividedHisto too small norm axis", Hal::EInfo::kLowWarning);
196 return;
197 } else if (axis > fDim) {
198 Cout::PrintInfo(("HalDividedHist too large norm axis"), Hal::EInfo::kLowWarning);
199 return;
200 }
201 if (min > max) {
202 Cout::PrintInfo("HalDividedHisto min> max", Hal::EInfo::kLowWarning);
203 return;
204 } else if (min == max) {
205 Cout::PrintInfo("HalDividedHisto min == max", Hal::EInfo::kLowWarning);
206 return;
207 }
208 fNormMin[axis] = min;
209 fNormMax[axis] = max;
210 }
211
212 void DividedHisto1D::SetScale(Double_t scale) { fScale = scale; }
213
215 DividedHisto1D* div = (DividedHisto1D*) h;
216 if (div->fDim != this->fDim) {
217 Cout::PrintInfo(Form("Invalid dimensions in %s and %s", this->ClassName(), h->ClassName()), Hal::EInfo::kCriticalError);
218 }
219 for (int i = 0; i < fDim; i++) {
220 if (this->fNormMin[i] != div->fNormMin[i]) {
221 Cout::PrintInfo("Invalid minimum values in HalDividedHisto", Hal::EInfo::kError);
222 } else if (this->fNormMax[i] != div->fNormMax[i]) {
223 Cout::PrintInfo("Invalid maximum values in HalDividedHisto", Hal::EInfo::kError);
224 }
225 }
226 if ((this->fScale == 0 && div->fScale != 0) || (this->fScale > 0 && div->fScale == 0)) {
227 Cout::PrintInfo(("Incompatibile scales !"), Hal::EInfo::kError);
228 fScale = 1.0;
229 }
230 if (fScale == 0) { // use norms
231 this->fNum->Add(div->fNum);
232 this->fDen->Add(div->fDen);
233 } else {
234 this->fScale = (div->fScale + this->fScale) / 2.0;
235 this->fNum->Add(div->fNum);
236 this->fDen->Add(div->fDen);
237 }
238 }
239
240 void DividedHisto1D::SetHistoName(TH1* h) const {
241 switch (fDim) {
242 case 1: h->GetYaxis()->SetTitle(fAxisName); break;
243 case 2: h->GetZaxis()->SetTitle(fAxisName); break;
244 default: break;
245 }
246 }
247
248 void DividedHisto1D::SetAxisName(TString name) { fAxisName = name; }
249
250 Double_t DividedHisto1D::GetNormMin(Int_t no) const {
251 if (no >= 0 && no < fDim) { return fNormMin[no]; }
252 Cout::PrintInfo(Form("%s::GetNormMin invalid no", ClassName()), Hal::EInfo::kLowWarning);
253 return 0;
254 }
255
256 Double_t DividedHisto1D::GetNormMax(Int_t no) const {
257 if (no >= 0 && no < fDim) { return fNormMax[no]; }
258 Cout::PrintInfo(Form("%s::GetNormMax invalid no", ClassName()), Hal::EInfo::kLowWarning);
259 return 0;
260 }
261
262 Double_t DividedHisto1D::GetScale() const { return fScale; }
263
264 TString DividedHisto1D::HTMLExtract(Int_t counter, TString dir) const {
265 Bool_t batch = gROOT->IsBatch();
266 TString path = Form("%s/divided_%i", dir.Data(), counter);
267 gSystem->MakeDirectory(path);
268 gROOT->SetBatch(kTRUE);
269 TCanvas* c1 = new TCanvas("canvas", "canvas", 0, 0, 1000, 1500);
270 c1->Divide(2, 1);
271 c1->cd(1);
272 TH1* cf = GetHist(kTRUE);
273 cf->Draw();
274 c1->cd(2);
275 TH1* num = (TH1*) GetNum()->Clone();
276 TH1* den = (TH1*) GetDen()->Clone();
277 den->SetLineColor(kRed);
278 num->SetLineColor(kGreen);
279 den->Draw();
280 num->Draw("SAME");
281 TLegend* leg = new TLegend(0.4, 0.8, 0.9, 0.9);
282 leg->AddEntry(num, "Numerator", "LPM");
283 leg->AddEntry(den, "Denominator", "LPM");
284 leg->Draw("SAME");
285 gSystem->mkdir(path);
286 c1->SaveAs(Form("%s/divided.root", path.Data()));
287 delete c1;
288 delete leg;
289 delete num;
290 delete den;
291 delete cf;
292 gROOT->SetBatch(batch);
293 TString page = CommonExtract(counter, dir);
294 return page;
295 }
296
297 TString DividedHisto1D::GetExtractType() const { return "Page"; }
298
299 TString DividedHisto1D::CommonExtract(Int_t counter, TString dir) const {
300 TString path = Form("%s/divided_%i", dir.Data(), counter);
301 gSystem->MakeDirectory(path);
302 TString filename = Form("%s/divided.html", path.Data());
303 HtmlFile file(filename, kFALSE);
304 Int_t dimN = GetNDim();
305 Double_t* min = new Double_t[dimN];
306 Double_t* max = new Double_t[dimN];
307 for (int i = 0; i < dimN; i++) {
308 min[i] = GetNormMin(i);
309 max[i] = GetNormMax(i);
310 }
311 HtmlTable table("", "haltable", "");
312 HtmlRow row1("", "dark_blue", "");
313 row1.AddContent(HtmlCell("label"));
314 row1.AddContent(HtmlCellCol("value", 2));
315 table.AddContent(row1);
316 HtmlRow row2("", "dark_blue", "");
317 row2.AddContent(HtmlCell("Name"));
318 row2.AddContent(HtmlCellCol(GetName(), 2));
319 table.AddContent(row2);
320 HtmlRow row3("", "dark_blue", "");
321 row3.AddContent(HtmlCell("Comment"));
323 table.AddContent(row3);
324 if (GetScale() == 0) {
325 HtmlRow rowS("", "dark_blue", "");
326 rowS.AddContent(HtmlCell("Scale"));
327 HtmlCellCol cellscale("0(disabled)", 2);
328 rowS.AddContent(cellscale);
329 table.AddContent(rowS);
330 } else {
331 }
332 for (int i = 0; i < dimN; i++) {
333 TString label;
334 switch (i) {
335 case 0: label = "x"; break;
336 case 1: label = "y"; break;
337 default: label = "z"; break;
338 }
339 HtmlRow row4("", "dark_blue", "");
340 row4.AddContent(HtmlCell(Form("Norm on %s axis", label.Data())));
341 row4.AddContent(HtmlCell(Form("%4.2f", min[i])));
342 row4.AddContent(HtmlCell(Form("%4.2f", max[i])));
343 table.AddContent(row4);
344 }
345 if (fBinomial) {
346 HtmlRow row5("", "dark_blue", "");
347 row5.AddContent(HtmlCell("Binomial"));
348 row5.AddContent(HtmlCellCol("kTRUE", 2));
349 table.AddContent(row5);
350 } else {
351 HtmlRow row5("", "dark_blue", "");
352 row5.AddContent(HtmlCell("Binomial"));
353 row5.AddContent(HtmlCellCol("kFALSE", 2));
354 table.AddContent(row5);
355 }
356 file.AddContent(table);
357 HtmlTable table2("", "haltable", "");
358 HtmlRow row2_1("", "dark_blue", "");
359 row2_1.AddContent(HtmlCell("No."));
360 row2_1.AddContent(HtmlCell("Name/Value"));
361 table2.AddContent(row2_1);
362 if (fLabels) // for backward compatilibiyt
363 for (int i = 0; i < fLabels->GetEntries(); i++) {
364 TObjString* obj = (TObjString*) fLabels->At(i);
365 HtmlRow singleRow("", "light_blue", "");
366 singleRow.AddContent(HtmlCell(Form("%i", i)));
367 singleRow.AddContent(HtmlCell(obj->GetString()));
368 table2.AddContent(singleRow);
369 }
370 file.AddContent(table2);
371 TString file_name = Form("%s/divided.root", path.Data());
372 std::ifstream f(file_name.Data());
373 if (f.good()) {
374 file.AddStringContent(HtmlCore::GetJsDiv("divided.root", "canvas;1"));
375 } else {
376 file.AddStringContent(GetPic());
377 }
378 file.Save();
379 return HtmlCore::GetUrl(Form("divided_%i/divided.html", counter), this->ClassName());
380 }
381
382 TString DividedHisto1D::GetLabel(Int_t i) const {
383 TObjString* obj_str = (TObjString*) fLabels->At(i);
384 if (obj_str == NULL) return "";
385 return obj_str->GetString();
386 }
387
388 Long64_t DividedHisto1D::Merge(TCollection* collection) {
389 if (collection) {
390 DividedHisto1D* pack = NULL;
391 TIter iterator(collection);
392 while ((pack = (DividedHisto1D*) iterator())) {
393 Add(pack);
394 }
395 }
396 return 1;
397 }
398
399 void DividedHisto1D::Browse(TBrowser* b) {
400 gPad->Clear();
401 TVirtualPad* c1 = gPad;
402 Draw();
403 gPad = c1;
404 b->Add(fNum);
405 b->Add(fDen);
406 }
407
408 void DividedHisto1D::SetBinomial(Bool_t binomial) { fBinomial = binomial; }
409
410 void DividedHisto1D::Fold1D(Double_t val, Option_t* opt) {
411 TString option = opt;
412 if (fNum == NULL || fDen == NULL) {
413 Cout::PrintInfo("lack of numerator and/or denominator folding cannot be done", Hal::EInfo::kLowWarning);
414 return;
415 }
416 TAxis* ax = NULL;
417 if (Hal::Std::FindParam(option, "x")) {
418 ax = fNum->GetXaxis();
419 } else if (Hal::Std::FindParam(option, "y")) {
420 ax = fNum->GetYaxis();
421 } else if (Hal::Std::FindParam(option, "z")) {
422 ax = fNum->GetZaxis();
423 }
424 Int_t bin = ax->FindBin(val);
425 Double_t cent = ax->GetBinCenter(bin);
426 Double_t low = ax->GetBinLowEdge(bin);
427 if (TMath::Abs(cent - val) > 1E-10 && TMath::Abs(low - val) > 1E-10) {
428 Cout::PrintInfo("value of folding is not equal a bin center of bin edge, folding "
429 "cannot be done",
430 Hal::EInfo::kLowWarning);
431 return;
432 }
433 Folding1D(val, opt);
434 }
435
436 void DividedHisto1D::Folding1D(Double_t val, Option_t* /*axis*/) {
437 // TODO improve this
438 TH1* num_c = (TH1*) fNum->Clone();
439 TH1* den_c = (TH1*) fDen->Clone();
440 Double_t min = fNum->GetXaxis()->GetBinLowEdge(1);
441 Double_t max = fNum->GetXaxis()->GetBinUpEdge(fNum->GetNbinsX());
442 Double_t period = max - min;
443 Double_t middle_bin = fNum->GetXaxis()->FindBin(val);
444 middle_bin = fNum->GetBinCenter(middle_bin);
445 for (int i = 1; i <= fNum->GetNbinsX(); i++) {
446 Double_t x = fNum->GetBinCenter(i);
447 Double_t op_x = 2.0 * val - x;
448 if (op_x > max) op_x = op_x - period;
449 if (op_x <= min) op_x = op_x + period;
450 Int_t j = fNum->FindBin(op_x);
451 Double_t N1 = num_c->GetBinContent(i);
452 Double_t NE1 = num_c->GetBinError(i);
453 Double_t D1 = den_c->GetBinContent(i);
454 Double_t DE1 = den_c->GetBinError(i);
455 Double_t N2 = num_c->GetBinContent(j);
456 Double_t NE2 = num_c->GetBinError(j);
457 Double_t D2 = den_c->GetBinContent(j);
458 Double_t DE2 = den_c->GetBinError(j);
459 Double_t Ncontent = N1 + N2;
460 Double_t Nerror = TMath::Sqrt(NE1 * NE1 + NE2 * NE2);
461 Double_t Dcontent = D1 + D2;
462 Double_t Derror = TMath::Sqrt(DE1 * DE1 + DE2 * DE2);
463 fNum->SetBinContent(i, Ncontent);
464 fNum->SetBinError(i, Nerror);
465 fDen->SetBinContent(i, Dcontent);
466 fDen->SetBinError(i, Derror);
467 }
468 delete num_c;
469 delete den_c;
470 }
471
472 DividedHisto1D::~DividedHisto1D() {
473 if (fNum && fNumIsCloned) { delete fNum; }
474 if (fDen && fDenIsCloned) { delete fDen; }
475 if (fNormMin) delete[] fNormMin;
476 if (fNormMax) delete[] fNormMax;
477 if (fLabels) delete fLabels;
478 }
479
480 void DividedHisto2D::Normalize(TH1* h) const {
481 if (fScale == 0) {
482 for (int i = 0; i < fDim; i++) {
483 Double_t scale = 1.0;
484 if (fNormMin[i] == fNormMax[i]) {
485 Cout::PrintInfo(Form("Same norm min and norm max on %i axis", i), Hal::EInfo::kLowWarning);
486 } else {
487 Int_t binMinX = fNum->GetXaxis()->FindBin(fNormMin[0]);
488 Int_t binMaxX = fNum->GetXaxis()->FindBin(fNormMax[0]);
489 Int_t binMinY = fNum->GetYaxis()->FindBin(fNormMin[1]);
490 Int_t binMaxY = fNum->GetYaxis()->FindBin(fNormMax[1]);
491 Double_t num = ((TH2*) fNum)->Integral(binMinX, binMaxX, binMinY, binMaxY);
492 binMinX = fDen->GetXaxis()->FindBin(fNormMin[0]);
493 binMaxX = fDen->GetXaxis()->FindBin(fNormMax[0]);
494 binMinY = fDen->GetYaxis()->FindBin(fNormMin[1]);
495 binMaxY = fDen->GetYaxis()->FindBin(fNormMax[1]);
496 Double_t den = ((TH2*) fDen)->Integral(binMinX, binMaxX, binMinY, binMaxY);
497 if (den == 0) {
498 Cout::PrintInfo("Numerator integral is zero, set scale to one", Hal::EInfo::kLowWarning);
499 } else {
500 scale = den / num;
501 }
502 }
503 h->Scale(scale);
504 }
505 } else {
506 h->Scale(fScale);
507 }
508 }
509
511
512 TString DividedHisto2D::HTMLExtract(Int_t counter, TString dir) const {
513 TString path = Form("%s/divided_%i", dir.Data(), counter);
514 gSystem->MakeDirectory(path);
515 Bool_t batch = gROOT->IsBatch();
516 gROOT->SetBatch(kTRUE);
517 TCanvas* c1 = new TCanvas("canvas", "canvas", 0, 0, 1000, 1500);
518 c1->Divide(2, 2);
519 c1->cd(1);
520 TH2* cf = (TH2*) GetHist(kTRUE);
521 cf->Draw("colz");
522 c1->cd(3);
523 TH2* num = (TH2*) GetNum()->Clone();
524 TH2* den = (TH2*) GetDen()->Clone();
525 den->Draw("colz");
526 c1->cd(4);
527 num->Draw("colz");
528 c1->cd(2);
529 c1->SaveAs(Form("%s/divided_%i/divided.root", dir.Data(), counter));
530 // create text file
531 gROOT->SetBatch(batch);
532 TString page = CommonExtract(counter, dir);
533 return page;
534 }
535
536 void DividedHisto2D::SetHistoName(TH1* h) const { h->GetZaxis()->SetName(fAxisName); }
537
540 return *this;
541 }
542
543 void DividedHisto2D::Folding1D(Double_t val, Option_t* axis) {
544 TAxis *X, *Y;
545 TString opt = axis;
546 Int_t order = 0;
547 if (Hal::Std::FindParam(opt, "x")) {
548 X = fNum->GetXaxis();
549 Y = fNum->GetYaxis();
550 } else {
551 order = 1;
552 Y = fNum->GetXaxis();
553 X = fNum->GetYaxis();
554 }
555 TH2* num_c = (TH2*) fNum->Clone();
556 TH2* den_c = (TH2*) fDen->Clone();
557 Double_t min = X->GetBinLowEdge(1);
558 Double_t max = X->GetBinUpEdge(X->GetNbins());
559 Double_t period = max - min;
560 Double_t middle_bin = X->FindBin(val);
561 middle_bin = X->GetBinCenter(middle_bin);
562 for (int i = 1; i <= X->GetNbins(); i++) {
563 Double_t x = X->GetBinCenter(i);
564 Double_t op_x = 2.0 * val - x;
565 if (op_x > max) op_x = op_x - period;
566 if (op_x <= min) op_x = op_x + period;
567 Int_t j = X->FindBin(op_x);
568 for (int k = 1; k <= Y->GetNbins(); k++) {
569 Int_t a, b, c, d;
570 if (order == 1) { // around y-axis
571 a = k;
572 b = i;
573 c = k;
574 d = j;
575 } else {
576 a = i;
577 b = k;
578 c = j;
579 d = k;
580 }
581 Double_t N1 = num_c->GetBinContent(a, b);
582 Double_t NE1 = num_c->GetBinError(a, b);
583 Double_t D1 = den_c->GetBinContent(a, b);
584 Double_t DE1 = den_c->GetBinError(a, b);
585 Double_t N2 = num_c->GetBinContent(c, d);
586 Double_t NE2 = num_c->GetBinError(c, d);
587 Double_t D2 = den_c->GetBinContent(c, d);
588 Double_t DE2 = den_c->GetBinError(c, d);
589 Double_t Ncontent = N1 + N2;
590 Double_t Nerror = TMath::Sqrt(NE1 * NE1 + NE2 * NE2);
591 Double_t Dcontent = D1 + D2;
592 Double_t Derror = TMath::Sqrt(DE1 * DE1 + DE2 * DE2);
593 fNum->SetBinContent(a, b, Ncontent);
594 fNum->SetBinError(a, b, Nerror);
595 fDen->SetBinContent(a, b, Dcontent);
596 fDen->SetBinError(a, b, Derror);
597 }
598 }
599 delete num_c;
600 delete den_c;
601 }
602
603 void DividedHisto2D::Folding2D(Double_t valX, Double_t valY, Option_t* opt) {
604 TString option = opt;
605 if (option.Length() == 1) {
606 Fold1D(valX, "x");
607 Fold1D(valY, "y");
608 return;
609 }
610 TH2* num_c = (TH2*) fNum->Clone();
611 TH2* den_c = (TH2*) fDen->Clone();
612 Double_t minX = fNum->GetXaxis()->GetBinLowEdge(1);
613 Double_t maxX = fNum->GetXaxis()->GetBinUpEdge(fNum->GetNbinsX());
614 Double_t minY = fNum->GetYaxis()->GetBinLowEdge(1);
615 Double_t maxY = fNum->GetYaxis()->GetBinUpEdge(fNum->GetNbinsY());
616 Double_t periodX = maxX - minX;
617 Double_t periodY = maxY - minY;
618 for (int i = 1; i <= fNum->GetNbinsX(); i++) {
619 for (int j = 1; j <= fNum->GetNbinsY(); j++) {
620 Double_t x = fNum->GetXaxis()->GetBinCenter(i);
621 Double_t y = fNum->GetYaxis()->GetBinCenter(j);
622 Double_t op_x = 2.0 * valX - x;
623 Double_t op_y = 2.0 * valY - y;
624 if (op_x > maxX) op_x = op_x - periodX;
625 if (op_x <= minX) op_x = op_x + periodX;
626 if (op_y > maxY) op_y = op_y - periodY;
627 if (op_y <= minY) op_y = op_y + periodY;
628 Int_t k = fNum->GetXaxis()->FindBin(op_x);
629 Int_t l = fNum->GetYaxis()->FindBin(op_y);
630 Double_t N1 = num_c->GetBinContent(i, j);
631 Double_t NE1 = num_c->GetBinError(i, j);
632 Double_t D1 = den_c->GetBinContent(i, j);
633 Double_t DE1 = den_c->GetBinError(i, j);
634 Double_t N2 = num_c->GetBinContent(k, l);
635 Double_t NE2 = num_c->GetBinError(k, l);
636 Double_t D2 = den_c->GetBinContent(k, l);
637 Double_t DE2 = den_c->GetBinError(k, l);
638 Double_t Ncontent = N1 + N2;
639 Double_t Nerror = TMath::Sqrt(NE1 * NE1 + NE2 * NE2);
640 Double_t Dcontent = D1 + D2;
641 Double_t Derror = TMath::Sqrt(DE1 * DE1 + DE2 * DE2);
642 fNum->SetBinContent(i, j, Ncontent);
643 fNum->SetBinError(i, j, Nerror);
644 fDen->SetBinContent(i, j, Dcontent);
645 fDen->SetBinError(i, j, Derror);
646 }
647 }
648 delete num_c;
649 delete den_c;
650 }
651
652 void DividedHisto2D::Fold2D(Double_t valX, Double_t valY, Option_t* opt) {
653 TString option = opt;
654 if (fNum == NULL || fDen == NULL) {
655 Cout::PrintInfo("lack of numerator and/or denominator folding cannot be done", Hal::EInfo::kLowWarning);
656 return;
657 }
658 TAxis* ax = NULL;
659 TAxis* ay = NULL;
660 if (option.EqualTo("yx")) {
661 Double_t temp = valX;
662 valX = valY;
663 valY = temp;
664 option = "xy";
665 ax = fNum->GetXaxis();
666 ay = fNum->GetYaxis();
667 }
668 if (option.EqualTo("zy")) {
669 Double_t temp = valX;
670 valX = valY;
671 valY = temp;
672 option = "yz";
673 ax = fNum->GetYaxis();
674 ay = fNum->GetZaxis();
675 }
676 if (option.EqualTo("zx")) {
677 Double_t temp = valX;
678 valX = valY;
679 valY = temp;
680 option = "xz";
681 ax = fNum->GetXaxis();
682 ay = fNum->GetZaxis();
683 }
684 if (option.EqualTo("x")) {
685 ax = fNum->GetYaxis();
686 ay = fNum->GetZaxis();
687 }
688 if (option.EqualTo("y")) {
689 ax = fNum->GetXaxis();
690 ay = fNum->GetZaxis();
691 }
692 if (option.EqualTo("z")) {
693 ax = fNum->GetXaxis();
694 ay = fNum->GetYaxis();
695 }
696 if (option.EqualTo("xy")) {
697 ax = fNum->GetXaxis();
698 ay = fNum->GetYaxis();
699 }
700 if (option.EqualTo("yz")) {
701 ax = fNum->GetYaxis();
702 ay = fNum->GetZaxis();
703 }
704 if (option.EqualTo("xz")) {
705 ax = fNum->GetXaxis();
706 ay = fNum->GetZaxis();
707 }
708 if (ax == NULL || ay == NULL) {
709 Cout::PrintInfo("Wrong option in 2D folding", Hal::EInfo::kLowWarning);
710 return;
711 }
712 Int_t binX = ax->FindBin(valX);
713 Int_t binY = ay->FindBin(valY);
714 Double_t centX = ax->GetBinCenter(binX);
715 Double_t lowX = ax->GetBinLowEdge(binX);
716 Double_t centY = ay->GetBinCenter(binY);
717 Double_t lowY = ay->GetBinLowEdge(binY);
718 if (TMath::Abs(centX - valX) > 1E-10 && TMath::Abs(lowX - valX) > 1E-10) {
719 Cout::PrintInfo("value of folding on X axis is not equal a bin center of bin edge, "
720 "folding cannot be done",
721 Hal::EInfo::kLowWarning);
722 return;
723 }
724 if (TMath::Abs(centY - valY) > 1E-10 && TMath::Abs(lowY - valY) > 1E-10) {
725 Cout::PrintInfo("value of folding on Y axis is not equal a bin center of bin edge, "
726 "folding cannot be done",
727 Hal::EInfo::kLowWarning);
728 return;
729 }
730 Folding2D(valX, valY, option);
731 }
732
733 void DividedHisto2D::Browse(TBrowser* b) {
734 gPad->Clear();
735 TVirtualPad* c1 = gPad;
736 if (!this->InheritsFrom("Hal::DividedHisto3D")) { // not 3d
737 Draw("all+colz");
738 } else {
739 Draw();
740 }
741 gPad = c1;
742 b->Add(fNum);
743 b->Add(fDen);
744 }
745
746 DividedHisto2D::~DividedHisto2D() {}
747
748 void DividedHisto3D::Normalize(TH1* h) const {
749 if (fScale == 0) {
750 for (int i = 0; i < fDim; i++) {
751 Double_t scale = 1.0;
752 if (fNormMin[i] == fNormMax[i]) {
753 Cout::PrintInfo(Form("Same norm min and norm max on %i axis", i), Hal::EInfo::kLowWarning);
754 scale = 1.0;
755 } else {
756 Int_t binMinX = fNum->GetXaxis()->FindBin(fNormMin[0]);
757 Int_t binMaxX = fNum->GetXaxis()->FindBin(fNormMax[0]);
758 Int_t binMinY = fNum->GetYaxis()->FindBin(fNormMin[1]);
759 Int_t binMaxY = fNum->GetYaxis()->FindBin(fNormMax[1]);
760 Int_t binMinZ = fNum->GetZaxis()->FindBin(fNormMin[2]);
761 Int_t binMaxZ = fNum->GetZaxis()->FindBin(fNormMax[2]);
762
763 Double_t num = ((TH3*) fNum)->Integral(binMinX, binMaxX, binMinY, binMaxY, binMinZ, binMaxZ);
764
765 binMinX = fDen->GetXaxis()->FindBin(fNormMin[0]);
766 binMaxX = fDen->GetXaxis()->FindBin(fNormMax[0]);
767 binMinY = fDen->GetYaxis()->FindBin(fNormMin[1]);
768 binMaxY = fDen->GetYaxis()->FindBin(fNormMax[1]);
769 binMinZ = fDen->GetZaxis()->FindBin(fNormMin[2]);
770 binMaxZ = fDen->GetZaxis()->FindBin(fNormMax[2]);
771
772 Double_t den = ((TH3*) fDen)->Integral(binMinX, binMaxX, binMinY, binMaxY, binMinZ, binMaxZ);
773
774 if (den == 0) {
775 Cout::PrintInfo("Numerator integral is zero, set scale to one", Hal::EInfo::kLowWarning);
776 scale = 1.0;
777 } else {
778 scale = den / num;
779 }
780 }
781 h->Scale(scale);
782 }
783 } else {
784 h->Scale(fScale);
785 }
786 }
787
789
791 const HistogramAxisConf& axX,
792 const HistogramAxisConf& axY,
793 const HistogramAxisConf& axZ,
794 Char_t type) :
795 DividedHisto3D(name,
796 axX.GetNBins(),
797 axX.GetMin(),
798 axX.GetMax(),
799 axY.GetNBins(),
800 axY.GetMin(),
801 axY.GetMax(),
802 axZ.GetNBins(),
803 axZ.GetMin(),
804 axZ.GetMax(),
805 type) {
806 TString nameX = axX.GetName();
807 if (nameX.Length()) {
808 GetNum()->GetXaxis()->SetTitle(nameX);
809 GetDen()->GetXaxis()->SetTitle(nameX);
810 }
811 TString nameY = axY.GetName();
812 if (nameY.Length()) {
813 GetNum()->GetYaxis()->SetTitle(nameY);
814 GetDen()->GetYaxis()->SetTitle(nameY);
815 }
816 TString nameZ = axZ.GetName();
817 if (nameZ.Length()) {
818 GetNum()->GetZaxis()->SetTitle(nameZ);
819 GetDen()->GetZaxis()->SetTitle(nameZ);
820 }
821 }
822
823 TString DividedHisto3D::HTMLExtract(Int_t counter, TString dir) const {
824 TString path = Form("%s/divided_%i", dir.Data(), counter);
825 gSystem->MakeDirectory(path);
826 TString page = CommonExtract(counter, dir);
827 Bool_t batch = gROOT->IsBatch();
828 gROOT->SetBatch(kTRUE);
829 TString name = "Divided 1D";
830 TCanvas* c1 = new TCanvas(name, name, 0, 0, 800, 600);
831 c1->Divide(2, 3);
832 c1->cd(1);
833 TH3* num = (TH3*) GetNum();
834 num->Draw(Draw_3D_option);
835 c1->cd(2);
836 // c1->GetPad(2)->SetLogy(1);
837 num->Draw(Draw_3D_option);
838 c1->cd(3);
839 TH3* den = (TH3*) GetDen();
840 den->Draw(Draw_3D_option);
841 c1->cd(4);
842 // c1->GetPad(4)->SetLogy(1);
843 den->Draw(Draw_3D_option);
844 c1->cd(5);
845 TH3* fun = (TH3*) GetHist(kTRUE);
846 fun->Draw(Draw_3D_option);
847 c1->cd(6);
848 // c1->GetPad(6)->SetLogy(1);
849 fun->Draw(Draw_3D_option);
850 c1->SaveAs(Form("%s/divided_%i/divided.png", dir.Data(), counter));
851 gROOT->SetBatch(batch);
852 return page;
853 }
854
855 TH2D* DividedHisto3D::GetProjection2D(Double_t min, Double_t max, Option_t* opt) const {
856 TString option = opt;
857 TH3* histo = NULL;
858 if (Hal::Std::FindParam(option, "num", kTRUE)) {
859 histo = (TH3*) GetNum()->Clone();
860 } else if (Hal::Std::FindParam(option, "den", kTRUE)) {
861 histo = (TH3*) GetDen()->Clone();
862 } else {
863 histo = (TH3*) GetHist(kTRUE);
864 }
865 TH2D* proj = Hal::Std::GetProjection2D(histo, min, max, opt);
866 delete histo;
867 return proj;
868 }
869
872 return *this;
873 }
874
875 void DividedHisto3D::Folding1D(Double_t val, Option_t* axis) {
876 TAxis *X, *Y, *Z;
877 TString opt = axis;
878 Int_t order = 0;
879 if (Hal::Std::FindParam(opt, "x")) {
880 X = fNum->GetXaxis();
881 Y = fNum->GetYaxis();
882 Z = fNum->GetZaxis();
883 } else if (Hal::Std::FindParam(opt, "y")) {
884 order = 1;
885 Y = fNum->GetXaxis();
886 X = fNum->GetYaxis();
887 Z = fNum->GetZaxis();
888 } else if (Hal::Std::FindParam(opt, "z")) {
889 order = 2;
890 X = fNum->GetZaxis();
891 Y = fNum->GetYaxis();
892 Z = fNum->GetXaxis();
893 } else {
894 X = Y = Z = NULL;
895 Cout::PrintInfo("Wrong option in HalDividedHisto3D Folding", Hal::EInfo::kLowWarning);
896 }
897 Double_t middle_bin = X->FindBin(val);
898 middle_bin = X->GetBinCenter(middle_bin);
899 TH3* num_c = (TH3*) fNum->Clone();
900 TH3* den_c = (TH3*) fDen->Clone();
901 Double_t min = X->GetBinLowEdge(1);
902 Double_t max = X->GetBinUpEdge(fNum->GetNbinsX());
903 Double_t period = max - min;
904 for (int i = 1; i <= X->GetNbins(); i++) {
905 Double_t x = X->GetBinCenter(i);
906 Double_t op_x = 2.0 * val - x;
907 if (op_x > max) op_x = op_x - period;
908 if (op_x <= min) op_x = op_x + period;
909 Int_t j = X->FindBin(op_x);
910 for (int k = 1; k <= Y->GetNbins(); k++) {
911 for (int l = 1; l <= Z->GetNbins(); l++) {
912 Int_t a, b, c, d, e, f;
913 if (order == 0) {
914 a = i;
915 b = k;
916 c = l;
917 d = j;
918 e = k;
919 f = l;
920 } else if (order == 1) {
921 a = k;
922 b = i;
923 c = l;
924 d = k;
925 e = j;
926 f = l;
927 } else {
928 a = k;
929 b = l;
930 c = i;
931 d = k;
932 e = l;
933 f = j;
934 }
935 Double_t N1 = num_c->GetBinContent(a, b, c);
936 Double_t NE1 = num_c->GetBinError(a, b, c);
937 Double_t D1 = den_c->GetBinContent(a, b, c);
938 Double_t DE1 = den_c->GetBinError(a, b, c);
939 Double_t N2 = num_c->GetBinContent(d, e, f);
940 Double_t NE2 = num_c->GetBinError(d, e, f);
941 Double_t D2 = den_c->GetBinContent(d, e, f);
942 Double_t DE2 = den_c->GetBinError(d, e, f);
943 Double_t Ncontent = N1 + N2;
944 Double_t Nerror = TMath::Sqrt(NE1 * NE1 + NE2 * NE2);
945 Double_t Dcontent = D1 + D2;
946 Double_t Derror = TMath::Sqrt(DE1 * DE1 + DE2 * DE2);
947 fNum->SetBinContent(a, b, c, Ncontent);
948 fNum->SetBinError(a, b, c, Nerror);
949 fDen->SetBinContent(a, b, c, Dcontent);
950 fDen->SetBinError(a, b, c, Derror);
951 }
952 }
953 }
954 delete num_c;
955 delete den_c;
956 }
957
958 void DividedHisto3D::Folding2D(Double_t valX, Double_t valY, Option_t* opt) {
959 TAxis *X, *Y, *Z;
960 Int_t order = 0;
961 X = fNum->GetXaxis();
962 Y = fNum->GetYaxis();
963 Z = fNum->GetZaxis();
964 TString option = opt;
965 if (option.EqualTo("xy")) {
966 order = 1;
967 } else if (option.EqualTo("xz")) {
968 order = 2;
969 } else if (option.EqualTo("yz")) {
970 order = 3;
971 } else if (option.EqualTo("x")) {
972 Fold1D(valX, "y");
973 Fold1D(valY, "z");
974 return;
975 } else if (option.EqualTo("y")) {
976 Fold1D(valX, "x");
977 Fold1D(valY, "z");
978 return;
979 } else if (option.EqualTo("z")) {
980 Fold1D(valX, "x");
981 Fold1D(valY, "y");
982 return;
983 } else {
984 return;
985 }
986 switch (opt[0]) {
987 case 'x': X = fNum->GetXaxis(); break;
988 case 'y': X = fNum->GetYaxis(); break;
989 case 'z': X = fNum->GetZaxis(); break;
990 }
991 switch (opt[1]) {
992 case 'y': Y = fNum->GetYaxis(); break;
993 case 'z': Y = fNum->GetZaxis(); break;
994 }
995 TH3* num_c = (TH3*) fNum->Clone();
996 TH3* den_c = (TH3*) fDen->Clone();
997 Double_t minX = X->GetBinLowEdge(1);
998 Double_t maxX = X->GetBinUpEdge(X->GetNbins());
999 Double_t periodX = maxX - minX;
1000 Double_t minY = Y->GetBinLowEdge(1);
1001 Double_t maxY = Y->GetBinUpEdge(Y->GetNbins());
1002 Double_t periodY = maxY - minY;
1003 for (int i = 1; i <= X->GetNbins(); i++) {
1004 Double_t x = X->GetBinCenter(i);
1005 Double_t op_x = 2.0 * valX - x;
1006 if (op_x > maxX) op_x = op_x - periodX;
1007 if (op_x <= minX) op_x = op_x + periodX;
1008 Int_t u = X->FindBin(op_x);
1009 for (int j = 1; j <= Y->GetNbins(); j++) {
1010 Double_t y = Y->GetBinCenter(j);
1011 Double_t op_y = 2.0 * valY - y;
1012 if (op_y > maxY) op_y = op_y - periodY;
1013 if (op_y <= minY) op_y = op_y + periodY;
1014 Int_t v = Y->FindBin(op_y);
1015 for (int k = 1; k <= Z->GetNbins(); k++) {
1016 Int_t a, b, c, d, e, f;
1017 switch (order) {
1018 case 1: { // xy
1019 a = i;
1020 b = j;
1021 c = k;
1022
1023 d = u;
1024 e = v;
1025 f = k;
1026 } break;
1027 case 2: { // xz
1028 a = i;
1029 b = k;
1030 c = j;
1031
1032 d = u;
1033 e = k;
1034 f = v;
1035 } break;
1036 case 3: { // yz
1037 a = k;
1038 b = i;
1039 c = j;
1040
1041 d = k;
1042 e = u;
1043 f = v;
1044 } break;
1045 }
1046 Double_t N1 = num_c->GetBinContent(a, b, c);
1047 Double_t NE1 = num_c->GetBinError(a, b, c);
1048 Double_t D1 = den_c->GetBinContent(a, b, c);
1049 Double_t DE1 = den_c->GetBinError(a, b, c);
1050 Double_t N2 = num_c->GetBinContent(d, e, f);
1051 Double_t NE2 = num_c->GetBinError(d, e, f);
1052 Double_t D2 = den_c->GetBinContent(d, e, f);
1053 Double_t DE2 = den_c->GetBinError(d, e, f);
1054
1055 Double_t Ncontent = N1 + N2;
1056 Double_t Nerror = TMath::Sqrt(NE1 * NE1 + NE2 * NE2);
1057 Double_t Dcontent = D1 + D2;
1058 Double_t Derror = TMath::Sqrt(DE1 * DE1 + DE2 * DE2);
1059 fNum->SetBinContent(a, b, c, Ncontent);
1060 fNum->SetBinError(a, b, c, Nerror);
1061 fDen->SetBinContent(a, b, c, Dcontent);
1062 fDen->SetBinError(a, b, c, Derror);
1063 }
1064 }
1065 }
1066 delete num_c;
1067 delete den_c;
1068 }
1069
1070 void DividedHisto3D::Fold3D(Double_t valX, Double_t valY, Double_t valZ, Option_t* opt) {
1071 TString option = opt;
1072 if (option.EqualTo("xyz")) {
1073 TH3* num_c = (TH3*) fNum->Clone();
1074 TH3* den_c = (TH3*) fDen->Clone();
1075 TAxis* X = fNum->GetXaxis();
1076 TAxis* Y = fNum->GetYaxis();
1077 TAxis* Z = fNum->GetZaxis();
1078 Double_t minX = X->GetBinLowEdge(1);
1079 Double_t maxX = X->GetBinUpEdge(X->GetNbins());
1080 Double_t periodX = maxX - minX;
1081 Double_t minY = Y->GetBinLowEdge(1);
1082 Double_t maxY = Y->GetBinUpEdge(Y->GetNbins());
1083 Double_t periodY = maxY - minY;
1084 Double_t minZ = Z->GetBinLowEdge(1);
1085 Double_t maxZ = Z->GetBinUpEdge(Z->GetNbins());
1086 Double_t periodZ = maxZ - minZ;
1087 for (int i = 1; i <= X->GetNbins(); i++) {
1088 Double_t x = X->GetBinCenter(i);
1089 Double_t op_x = 2.0 * valX - x;
1090 if (op_x > maxX) op_x = op_x - periodX;
1091 if (op_x <= minX) op_x = op_x + periodX;
1092 Int_t u = X->FindBin(op_x);
1093 for (int j = 1; j <= Y->GetNbins(); j++) {
1094 Double_t y = Y->GetBinCenter(j);
1095 Double_t op_y = 2.0 * valY - y;
1096 if (op_y > maxY) op_y = op_y - periodY;
1097 if (op_y <= minY) op_y = op_y + periodY;
1098 Int_t v = Y->FindBin(op_y);
1099 for (int k = 1; k <= Z->GetNbins(); k++) {
1100 Double_t z = Z->GetBinCenter(k);
1101 Double_t op_z = 2.0 * valZ - z;
1102 if (op_z > maxZ) op_z = op_z - periodZ;
1103 if (op_z <= minZ) op_z = op_z + periodZ;
1104 Int_t w = Z->FindBin(op_z);
1105
1106 Double_t N1 = num_c->GetBinContent(i, j, k);
1107 Double_t NE1 = num_c->GetBinError(i, j, k);
1108 Double_t D1 = den_c->GetBinContent(i, j, k);
1109 Double_t DE1 = den_c->GetBinError(i, j, k);
1110 Double_t N2 = num_c->GetBinContent(u, v, w);
1111 Double_t NE2 = num_c->GetBinError(u, v, w);
1112 Double_t D2 = den_c->GetBinContent(u, v, w);
1113 Double_t DE2 = den_c->GetBinError(u, v, w);
1114
1115 Double_t Ncontent = N1 + N2;
1116 Double_t Nerror = TMath::Sqrt(NE1 * NE1 + NE2 * NE2);
1117 Double_t Dcontent = D1 + D2;
1118 Double_t Derror = TMath::Sqrt(DE1 * DE1 + DE2 * DE2);
1119 fNum->SetBinContent(i, j, k, Ncontent);
1120 fNum->SetBinError(i, j, k, Nerror);
1121 fDen->SetBinContent(i, j, k, Dcontent);
1122 fDen->SetBinError(i, j, k, Derror);
1123 }
1124 }
1125 }
1126 delete num_c;
1127 delete den_c;
1128 } else {
1129 Fold1D(valX, "x");
1130 Fold1D(valY, "y");
1131 Fold1D(valZ, "z");
1132 }
1133 }
1134
1135 DividedHisto3D::~DividedHisto3D() {}
1136
1137 void DividedHisto3D::SetHistoName(TH1* /*h*/) const {}
1138
1139 void DividedHisto1D::SetComment(TString comment) { fComment = comment; }
1140
1141 void DividedHisto1D::Rebin(Int_t ngroup, Option_t* opt) {
1142 TString option = opt;
1143 Bool_t num_rebin = Hal::Std::FindParam(option, "num");
1144 Bool_t den_rebin = Hal::Std::FindParam(option, "den");
1145 if ((!num_rebin) && (!den_rebin)) {
1146 num_rebin = kTRUE;
1147 den_rebin = kTRUE;
1148 }
1149 if (num_rebin && fNum != NULL) { fNum->Rebin(ngroup); }
1150 if (den_rebin && fDen != NULL) { fDen->Rebin(ngroup); }
1151 }
1152
1153 void DividedHisto2D::Rebin(Int_t ngroup, Option_t* opt) {
1154 TString option = opt;
1155 Bool_t num_rebin = Hal::Std::FindParam(option, "num");
1156 Bool_t den_rebin = Hal::Std::FindParam(option, "den");
1157 Bool_t rebin_x = kFALSE;
1158 Bool_t rebin_y = kFALSE;
1159 if (Hal::Std::FindParam(option, "y")) {
1160 rebin_y = kTRUE;
1161 } else {
1162 rebin_x = kTRUE;
1163 }
1164 if ((!num_rebin) && (!den_rebin)) {
1165 num_rebin = kTRUE;
1166 den_rebin = kTRUE;
1167 }
1168 if (num_rebin && fNum != NULL) {
1169 if (rebin_x) ((TH2*) fNum)->RebinX(ngroup);
1170 if (rebin_y) ((TH2*) fNum)->RebinY(ngroup);
1171 }
1172 if (den_rebin && fDen != NULL) {
1173 if (rebin_x) ((TH2*) fDen)->RebinX(ngroup);
1174 if (rebin_y) ((TH2*) fDen)->RebinY(ngroup);
1175 }
1176 }
1177
1178 void DividedHisto3D::Rebin(Int_t ngroup, Option_t* opt) {
1179 TString option = opt;
1180 Bool_t num_rebin = Hal::Std::FindParam(option, "num");
1181 Bool_t den_rebin = Hal::Std::FindParam(option, "den");
1182 Bool_t rebin_x = kFALSE;
1183 Bool_t rebin_y = kFALSE;
1184 Bool_t rebin_z = kFALSE;
1185 if (option.Contains("x") || option.Contains("X")) { rebin_x = kTRUE; }
1186 if (option.Contains("y") || option.Contains("Y")) { rebin_y = kTRUE; }
1187 if (option.Contains("z") || option.Contains("Z")) { rebin_z = kTRUE; }
1188 if (rebin_y == kFALSE && rebin_z == kFALSE) { rebin_x = kTRUE; }
1189 if ((!num_rebin) && (!den_rebin)) {
1190 num_rebin = kTRUE;
1191 den_rebin = kTRUE;
1192 }
1193 if (num_rebin && fNum != NULL) {
1194 if (rebin_x) ((TH3*) fNum)->RebinX(ngroup);
1195 if (rebin_y) ((TH3*) fNum)->RebinY(ngroup);
1196 if (rebin_z) ((TH3*) fNum)->RebinZ(ngroup);
1197 }
1198 if (den_rebin && fDen != NULL) {
1199 if (rebin_x) ((TH3*) fDen)->RebinX(ngroup);
1200 if (rebin_y) ((TH3*) fDen)->RebinY(ngroup);
1201 if (rebin_z) ((TH3*) fDen)->RebinZ(ngroup);
1202 }
1203 }
1204
1205 void DividedHisto1D::AddLabel(TString label) { fLabels->AddLast(new TObjString(label)); }
1206
1208 if (this != &other) {
1209 if (fNum) delete fNum;
1210 if (fDen) delete fDen;
1211 if (other.fNum != NULL) fNum = (TH1*) other.fNum->Clone();
1212 if (other.fDen != NULL) fDen = (TH1*) other.fDen->Clone();
1213 for (int i = 0; i < fDim; i++) {
1214 fNormMin[i] = other.fNormMin[i];
1215 fNormMax[i] = other.fNormMax[i];
1216 }
1217 fScale = other.fScale;
1218 fAxisName = other.fAxisName;
1219 fComment = other.fComment;
1220 if (fLabels) {
1221 delete fLabels;
1222 fLabels = NULL;
1223 }
1224 if (other.fLabels) fLabels = (TList*) other.fLabels->Clone();
1225 }
1226 return *this;
1227 }
1228
1229 void DividedHisto1D::Draw(Option_t* opt) {
1230 if (gPad == NULL) new TCanvas();
1231 TString option = opt;
1232 if (Hal::Std::FindParam(option, "num", kTRUE)) {
1233 fNum->Draw(option);
1234 } else if (Hal::Std::FindParam(option, "den", kTRUE)) {
1235 fDen->Draw(option);
1236 } else if (Hal::Std::FindParam(option, "all", kTRUE)) {
1237 TVirtualPad* c1 = gPad;
1238 if (gPad->GetListOfPrimitives()->GetEntries() < 4) gPad->Divide(2, 2);
1239 c1->cd(1);
1240 fNum->Draw(option);
1241 c1->cd(2);
1242 fDen->Draw(option);
1243 c1->cd(3);
1244 TH1* h = this->GetHist(kTRUE);
1245 h->SetTitle(fComment);
1246 h->Draw(option);
1247 c1->cd(4);
1248 TH1* num_copy = (TH1*) fNum->Clone("temp_num");
1249 TH1* den_copy = (TH1*) fDen->Clone("temp_den");
1250 num_copy->GetYaxis()->SetTitle(Form("%s (scaled)", fNum->GetYaxis()->GetTitle()));
1251 num_copy->SetLineColor(kRed + 1);
1252 den_copy->SetLineColor(kGreen + 1);
1253 Double_t num_scale = 1.0 / num_copy->GetEntries();
1254 Double_t den_scale = 1.0 / den_copy->GetEntries();
1255 num_copy->Scale(num_scale);
1256 den_copy->Scale(den_scale);
1257 num_copy->SetStats(kFALSE);
1258 TLegend* leg = new TLegend(0.4, 0.8, 0.9, 0.9);
1259 leg->AddEntry(num_copy, Form("Numerator x %E", num_scale), "LPM");
1260 leg->AddEntry(den_copy, Form("Denominator x %E", den_scale), "LPM");
1261 num_copy->Draw(option);
1262 den_copy->Draw("SAME");
1263 leg->SetFillStyle(3002);
1264 leg->Draw("SAME");
1265 gPad = c1;
1266 } else {
1267 GetHist(kTRUE)->Draw(opt);
1268 }
1269 }
1270
1272 if (fLabels == NULL)
1273 return 0;
1274 else
1275 return fLabels->GetEntries();
1276 }
1277
1278 TString DividedHisto1D::GetPic() const {
1279 return Form("<img border=\"0\" src=\"divided.png\" width=\"996\" height=\"1472\">");
1280 }
1281
1282 DividedHisto1D::DividedHisto1D(TString name, Int_t nbins, Double_t min, Double_t max, Char_t type) : DividedHisto1D(name, 1) {
1283 std::vector<TString> out;
1284 out = Hal::Std::ExplodeString(name, ';');
1285 TString x_axis, y_axis;
1286 if (out.size() > 1) {
1287 name = out[0];
1288 x_axis = out[1];
1289 }
1290 if (out.size() > 2) y_axis = out[2];
1291 switch (type) {
1292 case 'D':
1293 fNum = new TH1D(name + "_Num", name + "_Num", nbins, min, max);
1294 fDen = new TH1D(name + "_Den", name + "_Den", nbins, min, max);
1295 break;
1296 case 'I':
1297 fNum = new TH1I(name + "_Num", name + "_Num", nbins, min, max);
1298 fDen = new TH1I(name + "_Den", name + "_Den", nbins, min, max);
1299 break;
1300 case 'F':
1301 fNum = new TH1F(name + "_Num", name + "_Num", nbins, min, max);
1302 fDen = new TH1F(name + "_Den", name + "_Den", nbins, min, max);
1303 break;
1304 default:
1305 fNum = new TH1D(name + "_Num", name + "_Num", nbins, min, max);
1306 fDen = new TH1D(name + "_Den", name + "_Den", nbins, min, max);
1307 break;
1308 }
1309 if (out.size() > 1) {
1310 fNum->GetXaxis()->SetTitle(x_axis);
1311 fDen->GetXaxis()->SetTitle(x_axis);
1312 }
1313 if (out.size() > 2) {
1314 fNum->GetYaxis()->SetTitle(y_axis);
1315 fDen->GetYaxis()->SetTitle(y_axis);
1316 }
1317 fNumIsCloned = fDenIsCloned = kTRUE;
1318 }
1319
1320 DividedHisto1D::DividedHisto1D(TString name, const HistogramAxisConf& conf, Char_t type) :
1321 DividedHisto1D(name, conf.GetNBins(), conf.GetMin(), conf.GetMax(), type) {
1322 TString nameX = conf.GetName();
1323 if (nameX.Length()) {
1324 GetNum()->GetXaxis()->SetTitle(nameX);
1325 GetDen()->GetXaxis()->SetTitle(nameX);
1326 }
1327 }
1328
1330 Int_t binsX,
1331 Double_t minX,
1332 Double_t maxX,
1333 Int_t binsY,
1334 Double_t minY,
1335 Double_t maxY,
1336 Char_t type) :
1337 DividedHisto2D(name, 2) {
1338 std::vector<TString> out = Hal::Std::ExplodeString(name, ';');
1339 TString x_axis, y_axis, z_axis;
1340 if (out.size() > 1) {
1341 name = out[0];
1342 x_axis = out[1];
1343 }
1344 if (out.size() > 2) y_axis = out[2];
1345 if (out.size() > 3) z_axis = out[3];
1346 switch (type) {
1347 case 'D':
1348 fNum = new TH2D(name + "_Num", name + "_Num", binsX, minX, maxX, binsY, minY, maxY);
1349 fDen = new TH2D(name + "_Den", name + "_Den", binsX, minX, maxX, binsY, minY, maxY);
1350 break;
1351 case 'I':
1352 fNum = new TH2I(name + "_Num", name + "_Num", binsX, minX, maxX, binsY, minY, maxY);
1353 fDen = new TH2I(name + "_Den", name + "_Den", binsX, minX, maxX, binsY, minY, maxY);
1354 break;
1355 case 'F':
1356 fNum = new TH2F(name + "_Num", name + "_Num", binsX, minX, maxX, binsY, minY, maxY);
1357 fDen = new TH2F(name + "_Den", name + "_Den", binsX, minX, maxX, binsY, minY, maxY);
1358 break;
1359 default:
1360 fNum = new TH2D(name + "_Num", name + "_Num", binsX, minX, maxX, binsY, minY, maxY);
1361 fDen = new TH2D(name + "_Den", name + "_Den", binsX, minX, maxX, binsY, minY, maxY);
1362 break;
1363 }
1364 if (out.size() > 1) {
1365 fNum->GetXaxis()->SetTitle(x_axis);
1366 fDen->GetXaxis()->SetTitle(x_axis);
1367 }
1368 if (out.size() > 2) {
1369 fNum->GetYaxis()->SetTitle(y_axis);
1370 fDen->GetYaxis()->SetTitle(y_axis);
1371 }
1372 if (out.size() > 3) {
1373 fNum->GetZaxis()->SetTitle(z_axis);
1374 fDen->GetZaxis()->SetTitle(z_axis);
1375 }
1376 SetOwner(kTRUE);
1377 }
1378 DividedHisto2D::DividedHisto2D(TString name, const HistogramAxisConf& axX, const HistogramAxisConf& axY, Char_t type) :
1379 DividedHisto2D(name, axX.GetNBins(), axX.GetMin(), axX.GetMax(), axY.GetNBins(), axY.GetMin(), axY.GetMax(), type) {
1380 TString nameX = axX.GetName();
1381 if (nameX.Length()) {
1382 GetNum()->GetXaxis()->SetTitle(nameX);
1383 GetDen()->GetXaxis()->SetTitle(nameX);
1384 }
1385 TString nameY = axY.GetName();
1386 if (nameY.Length()) {
1387 GetNum()->GetYaxis()->SetTitle(nameY);
1388 GetDen()->GetYaxis()->SetTitle(nameY);
1389 }
1390 }
1391
1393 Int_t binsX,
1394 Double_t minX,
1395 Double_t maxX,
1396 Int_t binsY,
1397 Double_t minY,
1398 Double_t maxY,
1399 Int_t binsZ,
1400 Double_t minZ,
1401 Double_t maxZ,
1402 Char_t type) :
1403 DividedHisto3D(name) {
1404 std::vector<TString> out = Hal::Std::ExplodeString(name, ';');
1405 TString x_axis, y_axis, z_axis;
1406 if (out.size() > 1) {
1407 name = out[0];
1408 x_axis = out[1];
1409 }
1410 if (out.size() > 2) y_axis = out[2];
1411 if (out.size() > 3) z_axis = out[3];
1412 switch (type) {
1413 case 'D':
1414 fNum = new TH3D(name + "_Num", name + "_Num", binsX, minX, maxX, binsY, minY, maxY, binsZ, minZ, maxZ);
1415 fDen = new TH3D(name + "_Den", name + "_Den", binsX, minX, maxX, binsY, minY, maxY, binsZ, minZ, maxZ);
1416 break;
1417 case 'I':
1418 fNum = new TH3I(name + "_Num", name + "_Num", binsX, minX, maxX, binsY, minY, maxY, binsZ, minZ, maxZ);
1419 fDen = new TH3I(name + "_Den", name + "_Den", binsX, minX, maxX, binsY, minY, maxY, binsZ, minZ, maxZ);
1420 break;
1421 case 'F':
1422 fNum = new TH3F(name + "_Num", name + "_Num", binsX, minX, maxX, binsY, minY, maxY, binsZ, minZ, maxZ);
1423 fDen = new TH3F(name + "_Den", name + "_Den", binsX, minX, maxX, binsY, minY, maxY, binsZ, minZ, maxZ);
1424 break;
1425 default:
1426 fNum = new TH3D(name + "_Num", name + "_Num", binsX, minX, maxX, binsY, minY, maxY, binsZ, minZ, maxZ);
1427 fDen = new TH3D(name + "_Den", name + "_Den", binsX, minX, maxX, binsY, minY, maxY, binsZ, minZ, maxZ);
1428 break;
1429 }
1430 if (out.size() > 1) {
1431 fNum->GetXaxis()->SetTitle(x_axis);
1432 fDen->GetXaxis()->SetTitle(x_axis);
1433 }
1434 if (out.size() > 2) {
1435 fNum->GetYaxis()->SetTitle(y_axis);
1436 fDen->GetYaxis()->SetTitle(y_axis);
1437 }
1438 if (out.size() > 3) {
1439 fNum->GetZaxis()->SetTitle(z_axis);
1440 fDen->GetZaxis()->SetTitle(z_axis);
1441 }
1442 SetOwner(kTRUE);
1443 }
1444
1446
1448
1449 Double_t DividedHisto1D::CalculateNorm(Double_t min, Double_t max) const {
1450 Int_t binA = fNum->GetXaxis()->FindBin(min);
1451 Int_t binB = fNum->GetXaxis()->FindBin(max);
1452 Double_t N = 0;
1453 Double_t D = 0;
1454 for (int i = binA; i <= binB; i++) {
1455 N += fNum->GetBinContent(i);
1456 D += fDen->GetBinContent(i);
1457 }
1458 return D / N;
1459 }
1460
1461 Double_t DividedHisto2D::CalculateNorm(Double_t minX, Double_t maxX, Double_t minY, Double_t maxY) const {
1462 Int_t binXm = fNum->GetXaxis()->FindBin(minX);
1463 Int_t binXM = fNum->GetXaxis()->FindBin(maxX);
1464 Int_t binYm = fNum->GetYaxis()->FindBin(minY);
1465 Int_t binYM = fNum->GetYaxis()->FindBin(maxY);
1466 Double_t N = 0;
1467 Double_t D = 0;
1468 for (int i = binXm; i <= binXM; i++) {
1469 for (int j = binYm; j <= binYM; j++) {
1470 N += fNum->GetBinContent(i, j);
1471 D += fDen->GetBinContent(i, j);
1472 }
1473 }
1474 return D / N;
1475 }
1476
1477 Double_t
1478 DividedHisto3D::CalculateNorm(Double_t minX, Double_t maxX, Double_t minY, Double_t maxY, Double_t minZ, Double_t maxZ) const {
1479 Int_t binXm = fNum->GetXaxis()->FindBin(minX);
1480 Int_t binXM = fNum->GetXaxis()->FindBin(maxX);
1481 Int_t binYm = fNum->GetYaxis()->FindBin(minY);
1482 Int_t binYM = fNum->GetYaxis()->FindBin(maxY);
1483 Int_t binZm = fNum->GetZaxis()->FindBin(minZ);
1484 Int_t binZM = fNum->GetZaxis()->FindBin(maxZ);
1485 Double_t N = 0;
1486 Double_t D = 0;
1487 for (int i = binXm; i <= binXM; i++) {
1488 for (int j = binYm; j <= binYM; j++) {
1489 for (int k = binZm; k <= binZM; k++) {
1490 N += fNum->GetBinContent(i, j, k);
1491 D += fDen->GetBinContent(i, j, k);
1492 }
1493 }
1494 }
1495 return D / N;
1496 }
1497
1498 TH1D* DividedHisto2D::Projection2DTo1D(Double_t min, Double_t max, Option_t* opt) const {
1499 TString option = opt;
1500 TH2* histo = NULL;
1501 if (Hal::Std::FindParam(option, "num", kTRUE)) {
1502 histo = (TH2*) GetNum()->Clone();
1503 } else if (Hal::Std::FindParam(option, "den", kTRUE)) {
1504 histo = (TH2*) GetDen()->Clone();
1505 } else {
1506 histo = (TH2*) GetHist(kTRUE);
1507 }
1508 TH1D* projection = Hal::Std::GetProjection1D(histo, min, max, option);
1509 delete histo;
1510 return projection;
1511 }
1512
1513 TH1D* DividedHisto3D::Projection3DTo1D(Double_t min1, Double_t max1, Double_t min2, Double_t max2, Option_t* opt) const {
1514 TString option = opt;
1515 TH3* histo = NULL;
1516 if (Hal::Std::FindParam(option, "num", kTRUE)) {
1517 histo = (TH3*) GetNum()->Clone();
1518 } else if (Hal::Std::FindParam(option, "den", kTRUE)) {
1519 histo = (TH3*) GetDen()->Clone();
1520 } else if (Hal::Std::FindParam(option, "separate", kTRUE)) {
1521 TH1D* num1d = Hal::Std::GetProjection1D((TH3*) fNum, min1, max1, min2, max2, option);
1522 TH1D* den1d = Hal::Std::GetProjection1D((TH3*) fDen, min1, max1, min2, max2, option);
1523 num1d->Divide(den1d);
1524 delete den1d;
1525 return num1d;
1526 } else {
1527 histo = (TH3*) GetHist(kTRUE);
1528 }
1529 TH1D* projection = Hal::Std::GetProjection1D(histo, min1, max1, min2, max2, option);
1530 delete histo;
1531 return projection;
1532 }
1533
1534 void DividedHisto1D::Print(Option_t* /*opt*/) const {
1535 Cout::Text(this->ClassName(), "R", kWhite);
1536 if (GetNum() != nullptr) {
1537 const TAxis* x = GetNum()->GetXaxis();
1538 TString bins = Form("Nbins: %i", x->GetNbins());
1539 TString min = Form("Xmin:%4.3f", x->GetBinLowEdge(1));
1540 TString max = Form("Xmax:%4.3f", x->GetBinUpEdge(x->GetNbins()));
1541 Cout::Text(bins, "R", kWhite);
1542 Cout::Text(min, "R", kWhite);
1543 Cout::Text(max, "R", kWhite);
1544 if (fDim > 1) {
1545 const TAxis* y = GetNum()->GetYaxis();
1546 bins = Form("Nbins: %i", y->GetNbins());
1547 min = Form("Ymin:%4.3f", y->GetBinLowEdge(1));
1548 max = Form("Ymax:%4.3f", y->GetBinUpEdge(y->GetNbins()));
1549 Cout::Text(bins, "R", kWhite);
1550 Cout::Text(min, "R", kWhite);
1551 Cout::Text(max, "R", kWhite);
1552 }
1553 if (fDim > 2) {
1554 const TAxis* z = GetNum()->GetZaxis();
1555 bins = Form("Nbins: %i", z->GetNbins());
1556 min = Form("Zmin:%4.3f", z->GetBinLowEdge(1));
1557 max = Form("Zmax:%4.3f", z->GetBinUpEdge(z->GetNbins()));
1558 Cout::Text(bins, "R", kWhite);
1559 Cout::Text(min, "R", kWhite);
1560 Cout::Text(max, "R", kWhite);
1561 }
1562 }
1563 }
1564
1565 void DividedHisto1D::AddScaled(const DividedHisto1D& other, Double_t scale) {
1566 fNum->Add(other.fNum, scale);
1567 fDen->Add(other.fDen, scale);
1568 TString className = ClassName();
1569 TString otherClassName = other.ClassName();
1570 if (className != otherClassName) {
1572 Form("%s %i: probably wrong AddScaled %s + %s", __FILE__, __LINE__, className.Data(), other.ClassName()), EInfo::kError);
1573 }
1574 }
1575
1577 h.Apply(*fNum);
1578 h.Apply(*fDen);
1579 }
1580
1581 void DividedHisto1D::SetDirectory(TDirectory* dir) {
1582 if (fNum && fNumIsCloned) fNum->SetDirectory(dir);
1583 if (fDen && fDenIsCloned) fDen->SetDirectory(dir);
1584 }
1585
1586} // namespace Hal
static void Text(TString text, TString option="L", Color_t color=-1)
Definition Cout.cxx:92
static void PrintInfo(TString text, Hal::EInfo status)
Definition Cout.cxx:370
void SetScale(Double_t scale)
void SetAxisName(TString name)
virtual TString CommonExtract(Int_t counter, TString dir) const
virtual void AddNum(TH1 *num, Option_t *opt="")
Double_t CalculateNorm(Double_t min, Double_t max) const
Double_t GetNormMax(Int_t no=0) const
TString GetLabel(Int_t i) const
virtual void ApplyStyle(const HistoStyle &h)
void SetComment(TString comment)
Double_t GetNormMin(Int_t no=0) const
virtual TString HTMLExtract(Int_t counter=0, TString dir=" ") const
virtual void SetDirectory(TDirectory *dir=nullptr)
Int_t GetNDim() const
virtual void Normalize(TH1 *h) const
void AddLabel(TString label)
TString GetExtractType() const
virtual void SetHistoName(TH1 *h) const
virtual void Rebin(Int_t ngroup, Option_t *opt)
Int_t GetLabelsNo() const
void Fold1D(Double_t val, Option_t *opt="x")
void SetNorm(Double_t min, Double_t max, Int_t axis=0)
void SetOwner(Bool_t own)
Double_t GetScale() const
virtual void Draw(Option_t *opt="all")
TH1 * GetHist(Bool_t normalized=kTRUE) const
virtual void AddNumDen(TH1 *num, TH1 *den, Option_t *opt="")
virtual void Add(const Object *h)
virtual TString GetPic() const
virtual void AddScaled(const DividedHisto1D &other, Double_t scale=1)
virtual void AddDen(TH1 *den, Option_t *opt="")
virtual void Browse(TBrowser *b)
virtual void Folding1D(Double_t val, Option_t *axis)
void SetBinomial(Bool_t binomial)
virtual Long64_t Merge(TCollection *collection)
DividedHisto1D & operator=(const DividedHisto1D &other)
virtual void Folding1D(Double_t val, Option_t *axis)
virtual void Rebin(Int_t ngroup, Option_t *opt)
void Fold2D(Double_t valX, Double_t valY, Option_t *opt="xy")
Double_t CalculateNorm(Double_t minX, Double_t maxX, Double_t minY, Double_t maxY) const
virtual void Normalize(TH1 *h) const
virtual void Browse(TBrowser *b)
virtual TString HTMLExtract(Int_t counter=0, TString dir=" ") const
TH1D * Projection2DTo1D(Double_t min, Double_t max, Option_t *opt) const
virtual void SetHistoName(TH1 *h) const
virtual void Folding2D(Double_t valX, Double_t valY, Option_t *opt)
DividedHisto2D & operator=(const DividedHisto2D &other)
virtual void Folding1D(Double_t val, Option_t *axis)
virtual void SetHistoName(TH1 *h) const
virtual void Fold3D(Double_t valX, Double_t valY, Double_t valZ, Option_t *opt)
virtual TH1D * Projection3DTo1D(Double_t min1, Double_t max1, Double_t min2, Double_t max2, Option_t *opt) const
virtual void Folding2D(Double_t valX, Double_t valY, Option_t *opt)
virtual TString HTMLExtract(Int_t counter=0, TString dir=" ") const
virtual TH2D * GetProjection2D(Double_t min, Double_t max, Option_t *opt) const
Double_t CalculateNorm(Double_t minX, Double_t maxX, Double_t minY, Double_t maxY, Double_t minZ, Double_t maxZ) const
virtual void Rebin(Int_t ngroup, Option_t *opt)
DividedHisto3D & operator=(const DividedHisto3D &other)
void Normalize(TH1 *h) const
static TString GetUrl(TString adress, TString text)
Definition HtmlCore.cxx:164
static TString GetJsDiv(TString root_file, TString object_name, TString draw_opt="colz", TString draw_div_name="drawing")
Definition HtmlCore.cxx:226
virtual void AddContent(const HtmlObject &obj)
Definition HtmlTable.cxx:29
virtual void AddContent(const HtmlObject &obj)
Definition HtmlTable.cxx:17