KaliVeda  1.12/06
Heavy-Ion Analysis Toolkit
KVHistoManipulator.cpp
Go to the documentation of this file.
1 //Created by KVClassFactory on Thu Oct 18 11:48:18 2007
2 //Author: Eric Bonnet
3 
4 #include "KVHistoManipulator.h"
5 #include "Riostream.h"
6 
7 #include "TProfile.h"
8 #include "TH1.h"
9 #include "TH2.h"
10 
11 #include "TF1.h"
12 #include "TCutG.h"
13 #include "TList.h"
14 #include "TString.h"
15 #include "TClass.h"
16 #include "TRandom3.h"
17 #include "TGraph.h"
18 #include "TROOT.h"
19 #include "TMethodCall.h"
20 #include "TMath.h"
21 
22 #include "KVNumberList.h"
23 #include "KVList.h"
24 
25 #include "TSpline.h"
26 #include "TStyle.h"
27 #include "TCanvas.h"
28 #include "TMultiGraph.h"
29 
30 #include <TObjString.h>
31 
32 using namespace std;
33 
35 
37 
38 
41 
43 {
44  //Default constructor
45  init();
46  gHistoManipulator = this;
47  fVDCanvas = nullptr;
48 }
49 
50 
51 
52 
54 
56 {
57  if (gHistoManipulator == this)
58  gHistoManipulator = nullptr;
59 }
60 
61 
62 //###############################################################################################################"
63 //-------------------------------------------------
64 
77 
79 {
80 //-------------------------------------------------
81  // IMPORTANT l'histo est modifie
82  //
83  // Efface les bins dont la statistique est hors de l 'intervalle ]stat_min,stat_max[
84  // l'erreur et le contenu sont mis a zero
85  // le nombre d entree (GetEntries) reste inchange
86  // pour les TH1 cela correspond a GetBinContent(xx)
87  // pour les TH2 cela correspond a GetBinContent(xx,yy) --> cellules
88  // pour les TProfile cela correspond a GetBinEntries(xx)
89  // la fonction renvoie le nbre de bins ou cellules mis a zero
90  // si stat_min ou stat_max sont egales à -1 la borne associée n'est pas prise en compte
91  if (!hh) {
92  cout << "pointeur histogramme nul" << endl;
93  return -1;
94  }
95  Int_t nb_raz = 0;
96  Bool_t raz = kFALSE;
97  if (hh->InheritsFrom("TH2")) {
98 
99  for (Int_t nx = 1; nx <= hh->GetNbinsX(); nx += 1) {
100  for (Int_t ny = 1; ny <= hh->GetNbinsY(); ny += 1) {
101  raz = kFALSE;
102  if (stat_min != -1 && hh->GetBinContent(nx, ny) < stat_min) raz = kTRUE;
103  if (stat_max != -1 && hh->GetBinContent(nx, ny) > stat_max) raz = kTRUE;
104  if (raz) {
105  hh->SetBinContent(nx, ny, 0);
106  hh->SetBinError(nx, ny, 0);
107  nb_raz += 1;
108  }
109  }
110  }
111  }
112  else if (hh->InheritsFrom("TProfile")) {
113  TProfile* prof = (TProfile*)hh;
114  for (Int_t nx = 1; nx <= prof->GetNbinsX(); nx += 1) {
115  raz = kFALSE;
116  if (stat_min != -1 && prof->GetBinEntries(nx) < stat_min) raz = kTRUE;
117  if (stat_max != -1 && prof->GetBinEntries(nx) > stat_max) raz = kTRUE;
118  if (raz) {
119  prof->SetBinContent(nx, 0);
120  prof->SetBinEntries(nx, 0);
121  prof->SetBinError(nx, 0);
122  nb_raz += 1;
123  }
124  }
125  }
126  else if (hh->InheritsFrom("TH1")) {
127  for (Int_t nx = 1; nx <= hh->GetNbinsX(); nx += 1) {
128  raz = kFALSE;
129  if (stat_min != -1 && hh->GetBinContent(nx) < stat_min) raz = kTRUE;
130  if (stat_max != -1 && hh->GetBinContent(nx) > stat_max) raz = kTRUE;
131  if (raz) {
132  hh->SetBinContent(nx, 0);
133  hh->SetBinError(nx, 0);
134  nb_raz += 1;
135  }
136  }
137  }
138 
139  return nb_raz;
140 }
141 
142 
143 //###############################################################################################################"
144 //-------------------------------------------------
145 
155 
157 {
158 //-------------------------------------------------
159  // IMPORTANT l'histo est modifie
160  //
161  // Applique une selection sur un histo 2D a partir d'un TCutG;
162  // Si mode=In seules les cellules comprises dans le TCutG sont gardes
163  // mode=Out --> Inverse
164  // la fonction renvoie le nbre de cellules mis a zero
165  // Attention le test ne se fait que sur les valeurs centrales des cellules (GetBinCenter())
166 
167  if (!hh) {
168  cout << "pointeur histogramme nul" << endl;
169  return -1;
170  }
171  Int_t nb_raz = 0;
172  Bool_t raz = kFALSE;
173 
174  for (Int_t nx = 1; nx <= hh->GetNbinsX(); nx += 1) {
175  for (Int_t ny = 1; ny <= hh->GetNbinsY(); ny += 1) {
176  raz = kFALSE;
177  Double_t valx = hh->GetXaxis()->GetBinCenter(nx);
178  Double_t valy = hh->GetYaxis()->GetBinCenter(ny);
179  if (mode == "in" && !cut->IsInside(valx, valy)) raz = kTRUE;
180  if (mode == "out" && cut->IsInside(valx, valy)) raz = kTRUE;
181  if (raz) {
182  hh->SetBinContent(nx, ny, 0);
183  hh->SetBinError(nx, ny, 0);
184  nb_raz += 1;
185  }
186  }
187  }
188 
189  return nb_raz;
190 }
191 
192 
193 //###############################################################################################################
194 //-------------------------------------------------
195 
221 
222 TH1* KVHistoManipulator::ScaleHisto(TH1* hh, TF1* fx, TF1* fy, Int_t nx, Int_t ny, Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax, Option_t* norm)
223 {
224 //-------------------------------------------------
225  // Applique les transformations correspondantes aux fonctions (TF1 *) donnees en argument
226  // et donne en sortie l histogramme transforme celui ci a pour nom "nom_de_histo_input"_scaled
227  //
228  // IMPORTANT l'histo d'entree n'est pas modifie
229  //
230  // pour TH1 et TProfile n'est pris en compte que TF1 *fx;
231  // les parametres de la fonction doivent etre initialises avant.
232  // Si la fonction est un pointeur NULL, aucune transformation n est appliquee et l axe reste tel quel.
233  // L'intervalle de validite des fonctions TF1::SetRange est determine a partir des limites de l histo apres transformations
234  // nx xmin et xmax sont les parametres de binnage des histo
235  // si nx=-1, ceux ci sont recalcules automatiquement nx = nx_initial, xmin = fx(xmin), xmax = fx(xmax)
236  // si nx!=-1, xmin et xmax doivent etre aussi specifies
237  // il en est de meme pour l'axe y
238  //
239  // OPTION: norm
240  // norm = "" (default) : no adjustment is made for the change in bin width due to the transformation
241  // norm = "width" : bin contents are adjusted for width change, so that the integral of the histogram
242  // contents taking into account the bin width (i.e. TH1::Integral("width")) is the same.
243  //
244  // NOTE for 1D histos:
245  // If the binning of the scaled histogram is imposed by the user (nx!=-1), then in order
246  // to achieve a continuous scaled distribution we have to randomize X within each bin.
247  // This will only work if the bin contents of 'hh' are integer numbers, i.e. unnormalised raw data histogram.
248 
249  Bool_t width = !strcmp(norm, "width");
250  TH1* gg = nullptr;
251  Double_t abs;
252  Bool_t fixed_bins = (nx != -1);
253  if (fx) {
254  if (nx == -1) {
255  nx = hh->GetNbinsX();
256  abs = hh->GetXaxis()->GetBinLowEdge(1);
257  xmin = fx->Eval(abs);
258  abs = hh->GetXaxis()->GetBinUpEdge(hh->GetNbinsX());
259  xmax = fx->Eval(abs);
260  }
261  fx->SetRange(hh->GetXaxis()->GetBinLowEdge(1), hh->GetXaxis()->GetBinUpEdge(nx));
262  fx->SetNpx(nx + 1);
263  }
264  else {
265  nx = hh->GetNbinsX();
266  xmin = hh->GetXaxis()->GetBinLowEdge(1);
267  xmax = hh->GetXaxis()->GetBinUpEdge(hh->GetNbinsX());
268  }
269 
270  if (hh->InheritsFrom("TH2")) {
271  if (fy) {
272  if (ny == -1) {
273  ny = hh->GetNbinsY();
274  abs = hh->GetYaxis()->GetBinLowEdge(1);
275  ymin = fy->Eval(abs);
276  abs = hh->GetYaxis()->GetBinUpEdge(hh->GetNbinsY());
277  ymax = fy->Eval(abs);
278  }
279  fy->SetRange(hh->GetYaxis()->GetBinLowEdge(1), hh->GetYaxis()->GetBinUpEdge(ny));
280  fy->SetNpx(ny + 1);
281  }
282  else {
283  ny = hh->GetNbinsY();
284  ymin = hh->GetYaxis()->GetBinLowEdge(1);
285  ymax = hh->GetYaxis()->GetBinUpEdge(hh->GetNbinsY());
286  }
287  }
288 
289  TClass* clas = TClass::GetClass(hh->ClassName());
290  gg = (TH1*) clas->New();
291  if (!gg) return nullptr;
292  TString hname;
293  hname.Form("%s_scaled", hh->GetName());
294  gg->SetNameTitle(hname.Data(), hh->GetTitle());
295 
296  if (hh->InheritsFrom("TH2")) gg->SetBins(nx, xmin, xmax, ny, ymin, ymax);
297  else gg->SetBins(nx, xmin, xmax);
298 
299  // if option norm="width', take account of change of X bin width between original and scaled histograms
300  Double_t Xbin_width_corr = 1.0, Ybin_width_corr = 1.0;
301  if (width) {
302  Double_t orig_Xbin_width = (hh->GetXaxis()->GetXmax() - hh->GetXaxis()->GetXmin()) / hh->GetNbinsX();
303  Double_t new_Xbin_width = (gg->GetXaxis()->GetXmax() - gg->GetXaxis()->GetXmin()) / gg->GetNbinsX();
304  Xbin_width_corr = orig_Xbin_width / new_Xbin_width;
305  if (hh->InheritsFrom("TH2")) {
306  // Take account of change of X bin width between original and scaled histograms
307  Double_t orig_Ybin_width = (hh->GetYaxis()->GetXmax() - hh->GetYaxis()->GetXmin()) / hh->GetNbinsY();
308  Double_t new_Ybin_width = (gg->GetYaxis()->GetXmax() - gg->GetYaxis()->GetXmin()) / gg->GetNbinsY();
309  Ybin_width_corr = orig_Ybin_width / new_Ybin_width;
310  }
311  }
312 
313  for (Int_t xx = 1; xx <= hh->GetNbinsX(); xx += 1) {
314  Double_t bmin = hh->GetXaxis()->GetBinLowEdge(xx);
315  Double_t bmax = hh->GetXaxis()->GetBinUpEdge(xx);
316  abs = gRandom->Uniform(bmin, bmax);
317  if (abs == bmax) abs = bmin;
318  Double_t resx = abs;
319  if (fx) resx = fx->Eval(abs);
320  if (hh->InheritsFrom("TH2")) {
321  for (Int_t yy = 1; yy <= hh->GetNbinsY(); yy += 1) {
322  if (hh->GetBinContent(xx, yy) > 0) {
323  bmin = hh->GetYaxis()->GetBinLowEdge(yy);
324  bmax = hh->GetYaxis()->GetBinUpEdge(yy);
325  abs = gRandom->Uniform(bmin, bmax);
326  if (abs == bmax) abs = bmin;
327  Double_t resy = abs;
328  if (fy) resy = fy->Eval(abs);
329  gg->SetBinContent(gg->GetXaxis()->FindBin(resx),
330  gg->GetYaxis()->FindBin(resy),
331  hh->GetBinContent(xx, yy)*Xbin_width_corr * Ybin_width_corr);
332  //gg->SetBinError(gg->GetXaxis()->FindBin(resx),gg->GetYaxis()->FindBin(resy),hh->GetBinError(xx,yy));
333  }
334  }
335  }
336  else {
337  // 1-D histogram
338  if (fixed_bins) {
339  // histogram binning imposed by user. we need to randomise inside bins of original histogram
340  // otherwise scaled histogram will be discontinuously filled.
341  Int_t nmax = (Int_t)hh->GetBinContent(xx);
342  for (int i = 0; i < nmax; i++) {
343  abs = gRandom->Uniform(bmin, bmax);
344  Double_t resx = fx->Eval(abs);
345  gg->Fill(resx, Xbin_width_corr);
346  //cout << resx << " " << Xbin_width_corr << endl;
347  }
348  }
349  else {
350  // range and binning calculated from function.
351  Double_t resy = hh->GetBinContent(xx);
352  if (fy) resy = fy->Eval(resy);
353  gg->SetBinContent(gg->GetXaxis()->FindBin(resx), resy * Xbin_width_corr);
354  }
355  }
356  }
357  return gg;
358 }
359 
360 
361 //###############################################################################################################
362 //-------------------------------------------------
363 
373 
375 {
376 //-------------------------------------------------
377  // Applique les transformations correspondantes aux fonctions (TF1 *) donnees en argument
378  // et donne en sortie le graph transforme celui ci a pour nom "nom_de_histo_input"_scaled
379  //
380  // IMPORTANT le graph d'entree n'est pas modifie
381  //
382  // les parametres de la fonction doivent etre initialises avant.
383  // Si la fonction est un pointeur NULL, aucune transformation n est appliquee et l axe reste tel quel.
384 
385  TGraph* gg = nullptr;
386  TClass* clas = TClass::GetClass(hh->ClassName());
387  gg = (TGraph*) clas->New();
388  if (!gg) return nullptr;
389  TString hname;
390  hname.Form("%s_scaled", hh->GetName());
391  gg->SetNameTitle(hname.Data(), hh->GetTitle());
392 
393  Int_t np = hh->GetN();
394  for (Int_t nn = 0; nn < np; nn += 1) {
395  Double_t xx1, yy1;
396  hh->GetPoint(nn, xx1, yy1);
397  Double_t xx2 = xx1;
398  if (fx) xx2 = fx->Eval(xx1);
399  Double_t yy2 = yy1;
400  if (fy) yy2 = fy->Eval(yy1);
401  gg->SetPoint(nn, xx2, yy2);
402  if (gg->InheritsFrom("TGraphErrors")) {
403  // transformation of errors
404  // if f = f(x) the error on f, e_f, for a given error on x, e_x, is
405  // e_f = abs0(df/dx) * e_x
406  Double_t e_x = ((TGraphErrors*)hh)->GetErrorX(nn);
407  Double_t e_y = ((TGraphErrors*)hh)->GetErrorY(nn);
408  if (fx) e_x = TMath::Abs(fx->Derivative(xx1)) * e_x;
409  if (fy) e_y = TMath::Abs(fy->Derivative(yy1)) * e_y;
410  ((TGraphErrors*)gg)->SetPointError(nn, e_x, e_y);
411  }
412  }
413 
414  return gg;
415 
416 }
417 
418 //###############################################################################################################
419 //-------------------------------------------------
420 
425 
427 {
428 //-------------------------------------------------
429  //Exemple d utilisation de la methode KVHistoManipulator::ScaleHisto avec ici
430  //L'obtention des distributions centrees reduites
431 
432  if (!hh) {
433  cout << "pointeur histogramme nul" << endl;
434  return hh;
435  }
436 
437  TString expression;
438  expression.Form("(x-%lf)/%lf", hh->GetMean(1), hh->GetRMS(1));
439  TF1 fx("fx", expression.Data());
440  unique_ptr<TF1> fy;
441  if (hh->InheritsFrom("TH2")) {
442  expression.Form("(x-%lf)/%lf", hh->GetMean(2), hh->GetRMS(2));
443  fy.reset(new TF1("fy", expression.Data()));
444  }
445 
446  TH1* gg = ScaleHisto(hh, &fx, fy.get(), nx, ny, xmin, xmax, ymin, ymax);
447  TString hname;
448  hname.Form("%s_centred", hh->GetName());
449  gg->SetName(hname.Data());
450 
451  return gg;
452 }
453 
454 //###############################################################################################################
455 //-------------------------------------------------
456 
461 
463 {
464 //-------------------------------------------------
465  //Exemple d utilisation de la methode KVHistoManipulator::ScaleHisto avec ici
466  //L'obtention des distributions centrees reduites
467 
468  if (!hh) {
469  cout << "pointeur histogramme nul" << endl;
470  return hh;
471  }
472 
473  TString expression;
474  expression.Form("(x-%lf)/%lf", hh->GetMean(1), hh->GetRMS(1));
475  TF1 fx("fx", expression.Data());
476  TH2* gg = (TH2*)ScaleHisto(hh, &fx, nullptr, nx, -1, xmin, xmax, -1., -1.);
477  TString hname;
478  hname.Form("%s_centred_X", hh->GetName());
479  gg->SetName(hname.Data());
480 
481  return gg;
482 }
483 
484 //###############################################################################################################
485 //-------------------------------------------------
486 
491 
493 {
494 //-------------------------------------------------
495  //Exemple d utilisation de la methode KVHistoManipulator::ScaleHisto avec ici
496  //L'obtention des distributions centrees reduites
497 
498  if (!hh) {
499  cout << "pointeur histogramme nul" << endl;
500  return hh;
501  }
502 
503  TString expression;
504  expression.Form("(x-%lf)/%lf", hh->GetMean(1), hh->GetRMS(1));
505  TF1 fy("fy", expression.Data());
506  TH2* gg = (TH2*)ScaleHisto(hh, nullptr, &fy, -1, ny, -1., -1., ymin, ymax);
507  TString hname;
508  hname.Form("%s_centred_Y", hh->GetName());
509  gg->SetName(hname.Data());
510 
511  return gg;
512 }
513 
514 
515 //###############################################################################################################
516 //-------------------------------------------------
517 
532 
534 {
535 //-------------------------------------------------
536  // Renormalisation de l histogramme 2D sous la contrainte d'une ditribution plate suivant X ou Y
537  // (signal de bimodalite par exemple)
538  // et donne en sortie l histogramme transforme celui ci a pour nom "nom_de_histo_input"_normalised
539  //
540  // IMPORTANT l'histo d'entree n'est pas modifie
541  //
542  // bmin, bmax : intervalle de bins ou l'on effectue la renormalisation (par defaut -1,-1 --> toute la plage)
543  // les contenus hors de l intervalle [bmin,bmax] sont remis a zero
544  // Si on veut une distribution plate en X ils indiquent l'intervalle sur cet axe x
545  // valref : valeur en stat de la distribution plate (par defaut 1)
546  //
547  // ATTENTION la propagation des erreurs n est pas (pour l instant) implementee
548 
549  if (!hh) {
550  cout << "pointeur histogramme nul" << endl;
551  return hh;
552  }
553 
554  if (bmin == -1) bmin = 1;
555  TString hname;
556  hname.Form("%s_normalised", hh->GetName());
557  TH2* clone = 0;
558  if ((clone = (TH2*)gDirectory->Get(hname.Data()))) delete clone;
559  clone = (TH2*)hh->Clone(hname.Data());
560  clone->Reset();
561 
562  if (axis == "X") {
563  if (bmax == -1) {
564  bmax = hh->GetNbinsX();
565  }
566  for (Int_t nx = bmin; nx <= bmax; nx += 1) {
567  Double_t integ = 0;
568  for (Int_t ny = 1; ny <= hh->GetNbinsY(); ny += 1) integ += hh->GetBinContent(nx, ny);
569  if (integ > 0) {
570  for (Int_t ny = 1; ny <= hh->GetNbinsY(); ny += 1) {
571  clone->SetBinContent(nx, ny, hh->GetBinContent(nx, ny)*valref / integ);
572  if (hh->GetBinContent(nx, ny) > 0) {
573  Double_t erreur = clone->GetBinContent(nx, ny) * TMath::Sqrt(1. / hh->GetBinContent(nx, ny) + 1. / integ);
574  clone->SetBinError(nx, ny, erreur);
575  }
576  }
577  }
578  }
579  }
580  else if (axis == "Y") {
581  if (bmax == -1) {
582  bmax = hh->GetNbinsY();
583  }
584  for (Int_t ny = bmin; ny <= bmax; ny += 1) {
585  Double_t integ = 0;
586  for (Int_t nx = 1; nx <= hh->GetNbinsX(); nx += 1) integ += hh->GetBinContent(nx, ny);
587  if (integ > 0) {
588  for (Int_t nx = 1; nx <= hh->GetNbinsX(); nx += 1) {
589  clone->SetBinContent(nx, ny, hh->GetBinContent(nx, ny)*valref / integ);
590  Double_t erreur = clone->GetBinContent(nx, ny) * TMath::Sqrt(1. / hh->GetBinContent(nx, ny) + 1. / integ);
591  clone->SetBinError(nx, ny, erreur);
592  }
593  }
594  }
595  }
596  else {
597  cout << "l option TString axis doit etre X ou Y" << endl;
598  }
599  return clone;
600 }
601 
602 //###############################################################################################################
603 //-------------------------------------------------
604 
610 
612 {
613 //-------------------------------------------------
614  // Les bornes valmin, valmax definissent l'intervalle ou l on effectue la renormalisation
615  // On en derive les valeurs de bins
616  // pour la methode KVHistoManipulator::RenormaliseHisto(TH1 *hh,TString axis,Double_t valref,Int_t bmin,Int_t bmax)
617 
618  if (!hh) {
619  cout << "pointeur histogramme nul" << endl;
620  return hh;
621  }
622  if (axis == "X") return RenormaliseHisto(hh, hh->GetXaxis()->FindBin(valmin), hh->GetXaxis()->FindBin(valmax), axis, valref);
623  else if (axis == "Y") return RenormaliseHisto(hh, hh->GetYaxis()->FindBin(valmin), hh->GetYaxis()->FindBin(valmax), axis, valref);
624  else {
625  cout << "l option TString axis doit etre X ou Y" << endl;
626  return nullptr;
627  }
628 
629 }
630 
631 
632 //###############################################################################################################
633 //-------------------------------------------------
634 
652 
654 {
655 //-------------------------------------------------
656  // Cumule le contenu de l histo hh entre bin bmin et bmax et retourne l histo correspondant
657  // Si direction="C" (default):
658  // SetBinContent(n) = GetBinContent(bmin)+GetBinContent(bmin+1)+ ... + GetBinContent(n)
659  // Si direction="D":
660  // SetBinContent(n) = GetBinContent(bmax)+GetBinContent(bmax-1)+ ... + GetBinContent(n)
661  //
662  // Donne en sortie l histogramme transforme celui ci a pour nom "nom_de_histo_input"_cumulated
663  //
664  // IMPORTANT l'histo d'entree n'est pas modifie
665  //
666  // si bmin=-1, bmin=1
667  // si bmax=-1, bmax=GetNbinsX()
668  //
669  // Avec norm = "surf" (default) l'integral de l histo cumule est egale a 1
670  // Avec norm = "max" le contenu de l'histogram est renormalise de facon a ce que le maximum soit 1
671 
672  if (!hh) {
673  cout << "pointeur histogramme nul" << endl;
674  return hh;
675  }
676  direction.ToUpper();
677  if (direction != "C" && direction != "D") {
678  cout << "l option TString direction doit etre C ou D" << endl;
679  return nullptr;
680  }
681  if (hh->GetDimension() == 1) {
682  if (bmin < 1) bmin = 1;
683  if (bmax < 1) bmax = hh->GetNbinsX();
684  TString hname;
685  hname.Form("%s_cumulated", hh->GetName());
686  TH1* clone = (TH1*)hh->Clone(hname.Data());
687  clone->Reset();
688  Double_t sum = 0;
689  Double_t big_sum = 0;
690  if (direction == "C") {
691  /* "standard" cumulative histogram, from bin 'bmin' to bin 'bmax' */
692  for (Int_t nx = 1; nx <= hh->GetNbinsX(); ++nx) {
693  if (nx < bmin) clone->SetBinContent(nx, 0);
694  else if (nx > bmax) clone->SetBinContent(nx, sum);
695  else {
696  sum += hh->GetBinContent(nx);
697  clone->SetBinContent(nx, sum);
698  }
699  big_sum += sum;
700  }
701  }
702  else {
703  /* "reverse" cumulative histogram, from bin 'bmax' to bin 'bmin' */
704  for (Int_t nx = hh->GetNbinsX(); nx >= 1; --nx) {
705  if (nx > bmax) clone->SetBinContent(nx, 0);
706  else if (nx < bmin) clone->SetBinContent(nx, sum);
707  else {
708  sum += hh->GetBinContent(nx);
709  clone->SetBinContent(nx, sum);
710  }
711  big_sum += sum;
712  }
713  }
714  // normalisation
715  if (!strcmp(norm, "surf")) {
716  clone->Scale(1. / big_sum);
717  }
718  else if (!strcmp(norm, "max")) {
719  clone->Scale(1. / sum);
720  }
721  return clone;
722  }
723  else {
724  cout << "cette methode n est pas prevue pour les TH2, TH3" << endl;
725  return nullptr;
726  }
727 
728 }
729 
730 
731 //###############################################################################################################
732 //-------------------------------------------------
733 
746 
748 {
749 //-------------------------------------------------
750  // retourne la derivee d ordre 0, 1 ou 2 d'un histogramme
751  // 0 -> correspond a un lissage (smooth) de la distribution
752  // 1 et 2 correspondent aux derivees premieres et deuxiemes
753  //
754  // 0 -> derivee zero Yi = 1/35*(-3yi-2 + 12yi-1 +17yi +12yi1 -3yi2)
755  // 1 -> derivee premiere Y'i = 1/12h*(yi-2 - 8yi-1 +8yi1 -1yi2)
756  // 2 -> derivee seconde Y''i = 1/7h/h*(2yi-2 -1 yi-1 -2yi -1yi1 +2yi2)
757  // les derivees pour les bins 1,2 et n-1,n ne sont pas calculees
758  //
759  // IMPORTANT l'histo d'entree n'est pas modifie
760 
761  if (!hh) {
762  cout << "pointeur histogramme nul" << endl;
763  return hh;
764  }
765  if (!(0 <= order && order <= 2)) {
766  cout << "ordre " << order << "n est pas implemente" << endl;
767  return nullptr;
768  }
769  if (hh->GetDimension() == 1) {
770 
771  TString hname;
772  hname.Form("%s_derivated_%d", hh->GetName(), order);
773  TH1* clone = 0;
774 
775  if (hh->InheritsFrom("TProfile")) clone = new TH1F(hname.Data(), hh->GetTitle(), hh->GetNbinsX(), hh->GetBinLowEdge(1), hh->GetBinLowEdge(hh->GetNbinsX() + 1));
776  else clone = (TH1*)hh->Clone(hname.Data());
777 
778  clone->Reset();
779 
780  for (Int_t nx = 3; nx <= hh->GetNbinsX() - 2; nx += 1) {
781  Double_t dev = 0;
782  if (order == 0) {
783  dev = 1. / 35.*(
784  -3 * hh->GetBinContent(nx - 2)
785  + 12 * hh->GetBinContent(nx - 1)
786  + 17 * hh->GetBinContent(nx)
787  + 12 * hh->GetBinContent(nx + 1)
788  - 3 * hh->GetBinContent(nx + 2)
789  );
790  }
791  else if (order == 1) {
792  Double_t h = hh->GetBinWidth(1);
793  dev = 1 / 12. / h * (
794  1 * hh->GetBinContent(nx - 2)
795  - 8 * hh->GetBinContent(nx - 1)
796  + 0 * hh->GetBinContent(nx)
797  + 8 * hh->GetBinContent(nx + 1)
798  - 1 * hh->GetBinContent(nx + 2)
799  );
800  }
801  else {
802  Double_t h2 = pow(hh->GetBinWidth(1), 2.);
803  dev = 1 / 7. / h2 * (
804  2 * hh->GetBinContent(nx - 2)
805  - 1 * hh->GetBinContent(nx - 1)
806  - 2 * hh->GetBinContent(nx)
807  - 1 * hh->GetBinContent(nx + 1)
808  + 2 * hh->GetBinContent(nx + 2)
809  );
810  }
811  clone->SetBinContent(nx, dev);
812  }
813  return clone;
814  }
815  else {
816  cout << "cette methode n est pas prevue pour les TH2, TH3" << endl;
817  return nullptr;
818  }
819 
820 }
821 
822 
823 //###############################################################################################################
824 //-------------------------------------------------
825 
842 
844 {
845 //-------------------------------------------------
846  // Renvoie un TGraph
847  // A partir d'un 2D histo, permet de tracer l'evolution d'un moment en fonction d'un autre moment
848  // ou d'un moment d'un axe en fonction de l'aure axe.
849  // Les TString momentx et momenty doivent etre des "Get like" methodes connues de TH1
850  // comme GetMean, GetRMS ou GetKurtosis :
851  // - Si momenty="" -> on obtient l'evolution du moment momentx en fonction de
852  // la variable associée a l'autre axe.
853  // - Si momenty!="" -> on obtient l'evolution du moment momenty en fonction du momentx
854  //
855  // Si TString axis = X la variable en X est le parametre d echantillonage et vice-versa
856  //
857  // Exemple :
858  // GetMomentEvolution(histo,"GetMean","GetSkewness","X") -> Evolution de la skewness en fonction de la moyenne de l'observable
859  // en Y par tranche de l'observable X
860 
861  if (axis != "X" && axis != "Y") {
862  cout << "GetMomentEvolution(TH2*,TString ,TString ,TString) Mauvaise syntaxe pour TString axis (X ou Y) " << endl;
863  return nullptr;
864  }
865  TMethodCall cmx;
866  cmx.InitWithPrototype(TClass::GetClass("TH1D"), Form("%s", momentx.Data()), "int");
867  if (!cmx.IsValid()) {
868  cout << "GetMomentEvolution(TH2*,TString ,TString ,TString) TString momentx n'est pas une methode valide " << momentx.Data() << endl;
869  return nullptr;
870  }
871  unique_ptr<TMethodCall> Ecmx(new TMethodCall());
872  Ecmx->InitWithPrototype(TClass::GetClass("TH1D"), Form("%sError", momentx.Data()), "int");
873 
874  unique_ptr<TMethodCall> cmy, Ecmy;
875  if (momenty != "") {
876  cmy.reset(new TMethodCall());
877  cmy->InitWithPrototype(TClass::GetClass("TH1D"), Form("%s", momenty.Data()), "int");
878  if (!cmy->IsValid()) {
879  cout << "GetMomentEvolution(TH2*,TString ,TString ,TString) TString momenty n'est pas une methode valide " << momenty.Data() << endl;
880  return nullptr;
881  }
882  Ecmy.reset(new TMethodCall());
883  Ecmy->InitWithPrototype(TClass::GetClass("TH1D"), Form("%sError", momenty.Data()), "int");
884  }
885 
886  TString fmt_histo;
887  fmt_histo.Form("GetMomentEvolution_%s", hh->GetName());
888  KVNumberList lbins = "";
889  Int_t nmax = 0;
890  if (axis == "Y") nmax = hh->GetNbinsX();
891  else nmax = hh->GetNbinsY();
892 
893  Int_t npts = 0;
894  for (Int_t nn = 1; nn <= nmax; nn += 1) {
895  Double_t stat = 0;
896  if (axis == "Y") stat = ((TH1D*)hh->ProjectionY(fmt_histo.Data(), nn, nn))->Integral();
897  else stat = ((TH1D*)hh->ProjectionX(fmt_histo.Data(), nn, nn))->Integral();
898  if (stat > stat_min) {
899  npts += 1;
900  lbins.Add(nn);
901  }
902  gDirectory->Delete(fmt_histo.Data());
903  }
904  if (npts > 0) {
905  TGraphErrors* gr = new TGraphErrors(npts);
906  TH1D* hp;
907  npts = 0;
908  Double_t valx, valy, Evaly = 0, Evalx = 0;
909  lbins.Begin();
910  while (!lbins.End()) {
911  Int_t bin = lbins.Next();
912  if (axis == "Y") hp = (TH1D*)hh->ProjectionY(fmt_histo.Data(), bin, bin);
913  else hp = (TH1D*)hh->ProjectionX(fmt_histo.Data(), bin, bin);
914  cmx.Execute(hp, "1", valx);
915  if (Ecmx->IsValid()) Ecmx->Execute(hp, "1", Evalx);
916  if (momenty != "") {
917  cmy->Execute(hp, "1", valy);
918  if (Ecmy->IsValid()) Ecmy->Execute(hp, "1", Evaly);
919  gr->SetPoint(npts, valx, valy);
920  if (Evalx != 0 && Evaly != 0) gr->SetPointError(npts, Evalx, Evaly);
921  npts += 1;
922  }
923  else {
924  if (axis == "Y") {
925  valy = hh->GetXaxis()->GetBinCenter(bin);
926  Evaly = hh->GetXaxis()->GetBinWidth(bin) / 2.;
927  }
928  else {
929  valy = hh->GetYaxis()->GetBinCenter(bin);
930  Evaly = hh->GetYaxis()->GetBinWidth(bin) / 2.;
931  }
932  gr->SetPoint(npts, valy, valx);
933  if (Evalx != 0) gr->SetPointError(npts, Evaly, Evalx);
934  npts += 1;
935  }
936  gDirectory->Delete(fmt_histo.Data());
937  }
938  return gr;
939  }
940  else {
941  cout << "GetMomentEvolution(TH2*,TString ,TString ,TString) Aucun point dans le TGraph" << endl;
942  return nullptr;
943  }
944 
945 }
946 
947 
948 
949 //-------------------------------------------------
950 
957 
959 {
960 //-------------------------------------------------
961  // A partir de deux graphs ayant le meme nombre de points et le meme axe X,
962  // cette methode produit un graph correspondant a la correlation entre les deux axes Y
963  // Les inputs peuvent etre aussi bien des TGraph que des TGraphErrors dans ce dernier cas
964  // les barres d erreurs sont prises en compte
965 
966  Double_t* yy = gry->GetY();
967  Double_t* xx = grx->GetY();
968  Double_t* ey = 0;
969  Double_t* ex = 0;
970 
971  if (grx->GetN() != gry->GetN()) {
972  printf("ERREUR : KVHistoManipulator::LinkGraphs : les deux graphs n ont pas le meme nbre de points\n");
973  return 0;
974  }
975  Int_t npoints = grx->GetN();
976 
977  TGraph* corre = 0;
978  if (grx->InheritsFrom("TGraphErrors") || gry->InheritsFrom("TGraphErrors")) {
979  if (grx->InheritsFrom("TGraphErrors")) ex = grx->GetEX();
980  if (gry->InheritsFrom("TGraphErrors")) ey = gry->GetEY();
981  corre = new TGraphErrors(npoints, xx, yy, ex, ey);
982  }
983  else corre = new TGraph(npoints, xx, yy);
984 
985  TString name;
986  name.Form("corre_%s_VS_%s", gry->GetName(), grx->GetName());
987  corre->SetNameTitle(name.Data(), grx->GetTitle());
988 
989  return corre;
990 
991 }
992 
993 
994 
995 //###############################################################################################################
996 //-------------------------------------------------
997 
1002 
1004 {
1005 //-------------------------------------------------
1006  // Calcul le coefficient de correlation entre les variables X et Y
1007  // d'un bidim... Equivalent a TH2F::GetCorrelationFactor() de ROOT
1008  Double_t sumxy = 0;
1009  Double_t sumx = 0, sumx2 = 0;
1010  Double_t sumy = 0, sumy2 = 0;
1011 
1012  Double_t compt = 0;
1013  for (Int_t nx = 1; nx <= hh->GetNbinsX(); nx += 1) {
1014  for (Int_t ny = 1; ny <= hh->GetNbinsY(); ny += 1) {
1015  compt += hh->GetBinContent(nx, ny);
1016  sumxy += hh->GetBinContent(nx, ny) * hh->GetXaxis()->GetBinCenter(nx) * hh->GetYaxis()->GetBinCenter(ny);
1017  sumx += hh->GetBinContent(nx, ny) * hh->GetXaxis()->GetBinCenter(nx);
1018  sumy += hh->GetBinContent(nx, ny) * hh->GetYaxis()->GetBinCenter(ny);
1019  sumx2 += hh->GetBinContent(nx, ny) * pow(hh->GetXaxis()->GetBinCenter(nx), 2.);
1020  sumy2 += hh->GetBinContent(nx, ny) * pow(hh->GetYaxis()->GetBinCenter(ny), 2.);
1021  }
1022  }
1023 
1024  Double_t meanxy = sumxy / compt;
1025  Double_t meanx = sumx / compt;
1026  Double_t meany = sumy / compt;
1027  Double_t meanx2 = sumx2 / compt;
1028  Double_t sigmax = sqrt(meanx2 - pow(meanx, 2.));
1029  Double_t meany2 = sumy2 / compt;
1030  Double_t sigmay = sqrt(meany2 - pow(meany, 2.));
1031 
1032  Double_t rho = (meanxy - meanx * meany) / (sigmax * sigmay);
1033  return rho;
1034 
1035 }
1036 
1037 //###############################################################################################################"
1038 //-------------------------------------------------
1039 
1045 
1047 {
1048 //-------------------------------------------------
1049  // Retourne une liste contenant les projections par tranche de l'axe (TString axis="X" ou "Y")
1050  // remplissant une condition leur integral qui doit etre superieur à MinIntegral (=-1 par defaut)
1051  // si axis="X", les projections sur l'axe Y de l'histogramme est fait pour chaque bin de l'axe X
1052 
1053  if (!hh) {
1054  cout << "pointeur histogramme nul" << endl;
1055  return NULL;
1056  }
1057  TH1D* h1d = 0;
1058  TString proj_name;
1059  if (hh->InheritsFrom("TH2")) {
1060  KVList* list = new KVList();
1061  cout << "TH2" << endl;
1062  if (axis == "X") {
1063  for (Int_t nx = 1; nx <= hh->GetNbinsX(); nx += 1) {
1064  Double_t integ = 0;
1065  for (Int_t ny = 1; ny <= hh->GetNbinsY(); ny += 1)
1066  integ += hh->GetBinContent(nx, ny);
1067 
1068  if (integ > MinIntegral) {
1069  proj_name.Form("%s_bX_%d", hh->GetName(), nx);
1070  h1d = hh->ProjectionY(proj_name.Data(), nx, nx);
1071  h1d->SetTitle(Form("%lf", hh->GetXaxis()->GetBinCenter(nx)));
1072  list->Add(h1d);
1073  }
1074 
1075  }
1076  }
1077  else if (axis == "Y") {
1078  for (Int_t ny = 1; ny <= hh->GetNbinsY(); ny += 1) {
1079  Double_t integ = 0;
1080  for (Int_t nx = 1; nx <= hh->GetNbinsX(); nx += 1)
1081  integ += hh->GetBinContent(nx, ny);
1082 
1083  if (integ > MinIntegral) {
1084  proj_name.Form("%s_bY_%d", hh->GetName(), ny);
1085  h1d = hh->ProjectionX(proj_name.Data(), ny, ny);
1086  h1d->SetTitle(Form("%lf", hh->GetYaxis()->GetBinCenter(ny)));
1087  list->Add(h1d);
1088  }
1089  }
1090  }
1091  else {
1092  cout << "l option TString axis doit etre X ou Y" << endl;
1093  }
1094 
1095  return list;
1096  }
1097  else {
1098  cout << "cette methode n est prevue que pour les TH2 and sons" << endl;
1099  return NULL;
1100  }
1101 
1102 }
1103 
1104 
1105 //-------------------------------------------------
1106 
1115 
1117 {
1118 //-------------------------------------------------
1119  // Donne les valeurs permettant de decouper l'integral
1120  // d'un histograme en tranche de 10% de la stat total
1121  // retourne une KVNumberList indiquant les delimitation des tranches
1122  //
1123  // L'utilisateur doit effacer cette liste apres utilisation
1124  //
1125 
1126  if (!hh) {
1127  cout << "pointeur histogramme nul" << endl;
1128  return NULL;
1129  }
1130 
1131  TString proj_name;
1132  Double_t integral = 0;
1133  for (Int_t nx = 1; nx <= hh->GetXaxis()->GetNbins(); nx += 1) {
1134  integral += hh->GetBinContent(nx);
1135  }
1136 
1137  KVNumberList* nl = new KVNumberList();
1138 
1139  Double_t tranche = 0;
1140  printf("integral du spectre %lf -> tranche de %lf\n", integral, integral / ntranches);
1141 
1142  Int_t nt = 0;
1143  for (Int_t nx = hh->GetXaxis()->GetNbins(); nx >= 1; nx -= 1) {
1144  tranche += hh->GetBinContent(nx);
1145  //printf("nx=%d, nt=%d, tranche[nt]=%lf -> %lf\n",nx,nt,tranche,integral/ntranches);
1146  if (tranche >= integral / ntranches) {
1147 
1148  nt += 1;
1149  //printf("\tnouvelle tranche %d %d\n",nx,nt);
1150  nl->Add(nx);
1151  tranche = 0;
1152 
1153  }
1154  }
1155 
1156  //Verif
1157  /*
1158  Int_t sup=1;
1159  Int_t inf = 1;
1160  nl->Begin();
1161  while (!nl->End()){
1162  sup = nl->Next();
1163  printf("%d %d - %lf -> %lf\n",inf,sup,hh->Integral(inf,sup),hh->Integral(inf,sup)/integral*100);
1164  inf = sup+1;
1165  }
1166  sup = hh->GetXaxis()->GetNbins();
1167  printf("%d %d - %lf -> %lf\n",inf,sup,hh->Integral(inf,sup),hh->Integral(inf,sup)/integral*100);
1168  */
1169  return nl;
1170 }
1171 
1172 
1173 //-------------------------------------------------
1174 
1181 
1183 {
1184 //-------------------------------------------------
1185  // Cree un histo 2D en intervertissant les axes
1186  //
1187  // L'utilisateur doit effacer cet histo apres utilisation
1188  //
1189 
1190  if (!hh) {
1191  cout << "pointeur histogramme nul" << endl;
1192  return NULL;
1193  }
1194  if (!hh->InheritsFrom("TH2")) {
1195  Error("PermuteAxis", "methode definie uniquement pour les classes TH2 et filles");
1196  }
1197  Int_t nx = hh->GetNbinsX();
1198  Int_t ny = hh->GetNbinsY();
1199 
1200  TH2F* gg = new TH2F(
1201  Form("%s_perm", hh->GetName()),
1202  hh->GetTitle(),
1203  ny,
1204  hh->GetYaxis()->GetBinLowEdge(1),
1205  hh->GetYaxis()->GetBinLowEdge(ny + 1),
1206  nx,
1207  hh->GetXaxis()->GetBinLowEdge(1),
1208  hh->GetXaxis()->GetBinLowEdge(nx + 1)
1209  );
1210 
1211  for (Int_t xx = 1; xx <= nx; xx += 1) {
1212  for (Int_t yy = 1; yy <= ny; yy += 1) {
1213 
1214  gg->SetBinContent(yy, xx, hh->GetBinContent(xx, yy));
1215 
1216  }
1217  }
1218  return gg;
1219 }
1220 
1221 
1222 //-------------------------------------------------
1223 
1230 
1232 {
1233 //-------------------------------------------------
1234  // Cree un TGraph en intervertissant les axes
1235  //
1236  // L'utilisateur doit effacer ce graph apres utilisation
1237  //
1238 
1239  if (!gr) {
1240  cout << "pointeur graph nul" << endl;
1241  return NULL;
1242  }
1243  if (!gr->InheritsFrom("TGraph")) {
1244  Error("PermuteAxis", "methode definie uniquement pour les classes TGraph et filles");
1245  }
1246 
1247  TGraphErrors* gr2 = new TGraphErrors();
1248  for (Int_t nn = 0; nn < gr->GetN(); nn += 1) {
1249  Double_t px, py;
1250  gr->GetPoint(nn, px, py);
1251  gr2->SetPoint(nn, py, px);
1252  if (gr->InheritsFrom("TGraphErrors")) {
1253  gr2->SetPointError(nn, ((TGraphErrors*)gr)->GetErrorY(nn), ((TGraphErrors*)gr)->GetErrorX(nn));
1254  }
1255  }
1256 
1257  return gr2;
1258 }
1259 
1260 
1261 //-------------------------------------------------
1262 
1269 
1271 {
1272 //-------------------------------------------------
1273  // Cree un graph à partir d un histo
1274  //
1275  // L'utilisateur doit effacer ce TGraph apres utilisation
1276  //
1277 
1278  if (!pf) {
1279  cout << "pointeur histogramme nul" << endl;
1280  return NULL;
1281  }
1282  if (!pf->InheritsFrom("TProfile")) {
1283  //Error("MakeGraphFrom","methode definie uniquement pour les classes TProfile et filles");
1284  return NULL;
1285  }
1286  Int_t nx = pf->GetNbinsX();
1287 
1288  TGraphErrors* gr = new TGraphErrors();
1289  for (Int_t xx = 1; xx <= nx; xx += 1) {
1290  if (pf->GetBinEntries(xx) > 0) {
1291  gr->SetPoint(gr->GetN(), pf->GetBinCenter(xx), pf->GetBinContent(xx));
1292  if (Error)
1293  gr->SetPointError(gr->GetN() - 1, pf->GetBinWidth(xx) / 2, pf->GetBinError(xx));
1294  }
1295  }
1296 
1297  return gr;
1298 }
1299 
1300 
1301 
1313 
1315 {
1316  // Extract mean and standard deviation (if generated with mode "profs") or error on mean (with "prof")
1317  // from TProfile and put their evolution in 2 TGraphs.
1318  // \returns Pointer to graph containing mean values
1319  // \param sigma Pointer to graph containing standard deviation
1320  //
1321  // Example of use:
1322  //~~~~{.cpp}
1323  // TGraph* sig_graf;
1324  // TGraph* mean_graf = HM.ExtractMeanAndSigmaFromProfile(pf, sig_graf);
1325  //~~~~
1326 
1327  Int_t nx = pf->GetNbinsX();
1328 
1329  TGraph* gr = new TGraph;
1330  sigma = new TGraph;
1331  int i = 0;
1332  for (Int_t xx = 1; xx <= nx; xx += 1) {
1333  if (pf->GetBinEntries(xx) > 0) {
1334  gr->SetPoint(i, pf->GetBinCenter(xx), pf->GetBinContent(xx));
1335  sigma->SetPoint(i, pf->GetBinCenter(xx), pf->GetBinError(xx));
1336  ++i;
1337  }
1338  }
1339 
1340  return gr;
1341 }
1342 
1343 
1344 // //-------------------------------------------------
1345 // TGraphErrors* KVHistoManipulator::MakeGraphFrom(TH1* pf, Bool_t Error)
1346 // {
1347 // //-------------------------------------------------
1348 // // Cree un graph à partir d un histo
1349 // //
1350 // // L'utilisateur doit effacer ce TGraph apres utilisation
1351 // //
1352 //
1353 // if (!pf) {
1354 // cout << "pointeur histogramme nul" << endl;
1355 // return NULL;
1356 // }
1357 // if (!pf->InheritsFrom("TH1")) {
1358 // //Error("MakeGraphFrom","methode definie uniquement pour les classes TProfile et filles");
1359 // }
1360 // else if (pf->InheritsFrom("TProfile"))
1361 // {
1362 // return
1363 // }
1364 // Int_t nx = pf->GetNbinsX();
1365 //
1366 // TGraphErrors* gr = new TGraphErrors();
1367 // for (Int_t xx = 1; xx <= nx; xx += 1) {
1368 // if (pf->GetBinEntries(xx) > 0) {
1369 // gr->SetPoint(gr->GetN(), pf->GetBinCenter(xx), pf->GetBinContent(xx));
1370 // if (Error)
1371 // gr->SetPointError(gr->GetN() - 1, pf->GetBinWidth(xx) / 2, pf->GetBinError(xx));
1372 // }
1373 // }
1374 //
1375 // return gr;
1376 // }
1377 
1378 //###############################################################################################################"
1379 //-------------------------------------------------
1380 
1383 
1384 void KVHistoManipulator::DefinePattern(TH1* ob, TString titleX, TString titleY, TString labelX, TString labelY)
1385 {
1386 //-------------------------------------------------
1387  DefinePattern(ob->GetXaxis(), titleX, labelX);
1388  DefinePattern(ob->GetYaxis(), titleY, labelY);
1389 
1390 }
1391 
1392 //###############################################################################################################"
1393 //-------------------------------------------------
1394 
1397 
1398 void KVHistoManipulator::DefinePattern(TGraph* ob, TString titleX, TString titleY, TString labelX, TString labelY)
1399 {
1400 //-------------------------------------------------
1401  DefinePattern(ob->GetXaxis(), titleX, labelX);
1402  DefinePattern(ob->GetYaxis(), titleY, labelY);
1403 
1404 }
1405 
1406 //###############################################################################################################"
1407 //-------------------------------------------------
1408 
1411 
1412 void KVHistoManipulator::DefinePattern(TF1* ob, TString titleX, TString titleY, TString labelX, TString labelY)
1413 {
1414 //-------------------------------------------------
1415  DefinePattern(ob->GetXaxis(), titleX, labelX);
1416  DefinePattern(ob->GetYaxis(), titleY, labelY);
1417 
1418 }
1419 
1420 //###############################################################################################################"
1421 //-------------------------------------------------
1422 
1425 
1427 {
1428 //-------------------------------------------------
1429 
1430  TObjArray* tok = NULL;
1431 
1432  if (!title.IsNull()) {
1433  tok = title.Tokenize(" ");
1434  if (tok->GetEntries() == 3) {
1435  ax->SetTitleFont(((TObjString*)tok->At(0))->GetString().Atoi());
1436  ax->SetTitleSize(((TObjString*)tok->At(1))->GetString().Atof());
1437  ax->SetTitleOffset(((TObjString*)tok->At(2))->GetString().Atof());
1438  }
1439  }
1440 
1441  if (!label.IsNull()) {
1442  tok = label.Tokenize(" ");
1443  if (tok->GetEntries() == 3) {
1444  ax->SetLabelFont(((TObjString*)tok->At(0))->GetString().Atoi());
1445  ax->SetLabelSize(((TObjString*)tok->At(1))->GetString().Atof());
1446  ax->SetLabelOffset(((TObjString*)tok->At(2))->GetString().Atof());
1447  }
1448  }
1449 
1450  if (tok) delete tok;
1451 
1452 }
1453 
1454 
1455 //###############################################################################################################"
1456 //-------------------------------------------------
1457 
1460 
1462 {
1463 //-------------------------------------------------
1464 
1465  TObjArray* tok = NULL;
1466  if (ob->IsA()->InheritsFrom("TAttLine")) {
1467  if (!line.IsNull()) {
1468  tok = line.Tokenize(" ");
1469  if (tok->GetEntries() == 3) {
1470  ob->SetLineColor(((TObjString*)tok->At(0))->GetString().Atoi());
1471  ob->SetLineStyle(((TObjString*)tok->At(1))->GetString().Atoi());
1472  ob->SetLineWidth(((TObjString*)tok->At(2))->GetString().Atoi());
1473  }
1474  }
1475  }
1476  if (tok) delete tok;
1477 
1478 }
1479 
1480 
1481 //###############################################################################################################"
1482 //-------------------------------------------------
1483 
1486 
1488 {
1489 //-------------------------------------------------
1490 
1491  TObjArray* tok = NULL;
1492  if (ob->IsA()->InheritsFrom("TAttMarker")) {
1493  if (!marker.IsNull()) {
1494  tok = marker.Tokenize(" ");
1495  if (tok->GetEntries() == 3) {
1496  ob->SetMarkerColor(((TObjString*)tok->At(0))->GetString().Atoi());
1497  ob->SetMarkerStyle(((TObjString*)tok->At(1))->GetString().Atoi());
1498  ob->SetMarkerSize(((TObjString*)tok->At(2))->GetString().Atof());
1499  }
1500  }
1501  }
1502  if (tok) delete tok;
1503 
1504 }
1505 
1506 
1507 //###############################################################################################################"
1508 //-------------------------------------------------
1509 
1512 
1514 {
1515 //-------------------------------------------------
1516 
1517  DefineLineStyle((TAttLine*)ob, line);
1518  DefineMarkerStyle((TAttMarker*)ob, marker);
1519 
1520 }
1521 
1522 
1523 //###############################################################################################################"
1524 //-------------------------------------------------
1525 
1528 
1530 {
1531 //-------------------------------------------------
1532 
1533  ob->GetXaxis()->SetTitle(xtit);
1534  ob->GetYaxis()->SetTitle(ytit);
1535 
1536 }
1537 
1538 //###############################################################################################################"
1539 //-------------------------------------------------
1540 
1543 
1545 {
1546 //-------------------------------------------------
1547 
1548  ob->GetXaxis()->SetTitle(xtit);
1549  ob->GetYaxis()->SetTitle(ytit);
1550 
1551 }
1552 
1553 //###############################################################################################################"
1554 //-------------------------------------------------
1555 
1558 
1560 {
1561 //-------------------------------------------------
1562 
1563  ob->GetXaxis()->SetTitle(xtit);
1564  ob->GetYaxis()->SetTitle(ytit);
1565 
1566 }
1567 
1568 //###############################################################################################################"
1569 //-------------------------------------------------
1570 
1580 
1582 {
1583  // Return value of abscissa X for which the interpolated value
1584  // of the histogram contents is equal to the given value, val.
1585  // We use the false position method (which should always converge...)
1586  // eps is required precision, i.e. convergence condition is that no further change
1587  // in result greater than eps is found.
1588  // nmax = maximum number of iterations
1589  // A solution is searched for X between the limits Xmin and Xmax of the X axis of the histo
1590  // unless arguments (xmin,xmax) are given to bracket the search
1591 
1592  TSpline5* spline = new TSpline5(ob);
1593  Double_t Xmax;
1594  Double_t Xmin;
1595  if (xmin == xmax) {
1596  Xmax = ob->GetXaxis()->GetXmax();
1597  Xmin = ob->GetXaxis()->GetXmin();
1598  }
1599  else {
1600  Xmax = xmax;
1601  Xmin = xmin;
1602  }
1603 
1604  Int_t n, side = 0;
1605  Double_t r = 0, fr, fs, s, ft, t;
1606  s = Xmin;
1607  t = Xmax;
1608  fs = spline->Eval(s) - val;
1609  ft = spline->Eval(t) - val;
1610 
1611  for (n = 1; n <= nmax; n++) {
1612  r = (fs * t - ft * s) / (fs - ft);
1613  if (fabs(t - s) < eps * fabs(t + s)) break;
1614  fr = spline->Eval(r) - val;
1615 
1616  if (fr * ft > 0) {
1617  t = r;
1618  ft = fr;
1619  if (side == -1) fs /= 2;
1620  side = -1;
1621  }
1622  else if (fs * fr > 0) {
1623  s = r;
1624  fs = fr;
1625  if (side == +1) ft /= 2;
1626  side = +1;
1627  }
1628  else break;
1629  }
1630  delete spline;
1631  return r;
1632 }
1633 
1634 //###############################################################################################################"
1635 //-------------------------------------------------
1636 
1674 
1675 TF1* KVHistoManipulator::RescaleX(TH1* hist1, TH1* hist2, Int_t degree, Double_t* params,
1676  Int_t npoints, const Char_t* direction, Double_t xmin, Double_t xmax, Double_t qmin, Double_t qmax, Double_t eps)
1677 {
1678  // Find the polynomial transformation of the X-axis of 1-D histogram hist1
1679  // so that the distribution ressembles that of histogram hist2.
1680  // Returns a TF1 which can be used to rescale hist1 (hist1 and hist2 are not modified).
1681  // User's responsibility to delete the TF1.
1682  // After fitting, the array params is filled with:
1683  // params[0] = a0
1684  // params[1] = a1
1685  // ...
1686  // params[degree] = an
1687  // params[degree+1] = Chisquare/NDF
1688  // Therefore params MUST BE at least of dimension (degree+2)
1689  //
1690  // npoints : by default (npoints=-1), we use npoints=degree+2 values of comparison
1691  // of the two cumulative distributions (see below).
1692  // more comparison points can be used by setting npoints>=degree+2.
1693  // direction : by default ("C") we use the cumulative histogram summed from low x to high x.
1694  // if direction="D", we use the cumulative histogram summed from high x to low x.
1695  // xmin, xmax : range of values of abscissa used to build cumulative histograms. default
1696  // (xmin=xmax=-1) include all values.
1697  // qmin, qmax : minimum & maximum values of cumulative histograms used for the
1698  // comparison (see below). by default qmin=0.05, qmax=0.95.
1699  //
1700  // METHOD
1701  // ======
1702  // 'hist1' contains the distribution P1 of variable X1
1703  // 'hist2' contains the distribution P2 of variable X2
1704  // Supposing that we can write X2=f_n(X1), with f_n(x)=a0 + a1*x + a2*x^2 + ... + an*x^n,
1705  // what are the parameters of the polynomial for which P1(f_n(X1))=P2(X2) ?
1706  //
1707  // Consider the cumulative distributions, C1(X1) and C2(X2).
1708  // For npoints=2 we compare the 2 X1 & X2 values for which C1=C2=qmin, qmax
1709  // For npoints=3 we compare the 3 X1 & X2 values for which C1=C2=qmin, qmin+(qmax-qmin)/2, qmax
1710  // For npoints=4 we compare the 4 X1 & X2 values for which C1=C2=qmin, qmin+(qmax-qmin)/3, qmin+2*(qmax-qmin)/3, qmax
1711  // etc. etc.
1712  // In each case we fit the npoints couples (X1,X2) with f_n
1713  //
1714 
1715  int i;
1716  // calculate comparison points
1717  npoints = TMath::Max(npoints, degree + 2);
1718  TString func_name = Form("pol%d", degree);
1719  TF1* fonc = new TF1("f", func_name.Data());
1720  fonc->SetName(Form("RescaleX-%s", func_name.Data()));
1721  RescaleX(hist1, hist2, fonc, npoints, direction, xmin, xmax, qmin, qmax, eps);
1722  for (i = 0; i < degree + 1; i++) {
1723  params[i] = fonc->GetParameter(i);
1724  }
1725  Double_t chisquare = fonc->GetChisquare();
1726  if (fonc->GetNDF() > 0.0) chisquare /= fonc->GetNDF();
1727  params[degree + 1] = chisquare;
1728 
1729  return fonc;
1730 }
1731 
1732 //###############################################################################################################"
1733 //-------------------------------------------------
1734 
1800 
1802  Option_t* opt, Int_t npoints, const Char_t* direction, Double_t xmin, Double_t xmax, Double_t qmin, Double_t qmax,
1803  Double_t eps)
1804 {
1805  // Uses RescaleX(TH1* hist1, TH1* hist2, Int_t degree, Double_t* params, Double_t eps)
1806  // to find a polynomial transformation of 'hist1' abscissa so that its
1807  // distribution resembles that of 'hist2', then generates a new rescaled version of 'hist1'.
1808  //
1809  // degree = degree of polynomial to use
1810  // params = array of dimension [degree+2], after rescaling it will contain
1811  // the values of fitted parameters of polynomial, plus the Chi2/NDF of the fit:
1812  // params[0] = a0
1813  // params[1] = a1
1814  // ...
1815  // params[degree] = an
1816  // params[degree+1] = Chisquare/NDF = Chisquare (NDF is always equal to 1)
1817  //
1818  // npoints : by default (npoints=-1), we use npoints=degree+2 values of comparison
1819  // of the two cumulative distributions (see method RescaleX).
1820  // more comparison points can be used by setting npoints>=degree+2.
1821  //
1822  // direction : by default ("C") we use the cumulative histogram summed from low x to high x.
1823  // if direction="D", we use the cumulative histogram summed from high x to low x.
1824  // xmin, xmax : range of values of abscissa used to build cumulative histograms. default
1825  // (xmin=xmax=-1) include all values.
1826  // qmin, qmax : minimum & maximum values of cumulative histograms used for the
1827  // comparison (see method RescaleX). by default qmin=0.05, qmax=0.95.
1828  //
1829  // eps = relative precision used to find comparison points of histos (default = 1.e-07)
1830  //
1831  // OPTIONS
1832  // =======
1833  // opt = "norm" : rescaled histogram normalised to have same integral as hist2
1834  // opt = "bins" : rescaled histogram will have same number of bins & same limits as hist2
1835  // opt = "normbins" |
1836  // opt = "binsnorm" |--> rescaled histogram can be superposed and added to hist2
1837  //
1838  // EXAMPLE OF USE
1839  // =======================
1840  // In the following example, we fill two histograms with different numbers of random
1841  // values drawn from two Gaussian distributions with different centroids and widths.
1842  // We also add to each histogram a 'pedestal' peak which is unrelated to the 'physical'
1843  // distributions, and significant enough so that it does not permit a correct scaling of
1844  // the physical part of the distribution.
1845  //
1846  // We can overcome the problem of the 'pedestal' by
1847  // using the 'inverse cumulated distribution' and excluding the channels
1848  // where the pedestals are present. The correct scaling of the physical
1849  // distributions is recovered, as shown below:
1850  /*
1851  BEGIN_MACRO(source)
1852  {
1853  TH1F* h10 = new TH1F("h10","gaussian",4096,-.5,4095.5);
1854  TF1 g10("g10","gaus(0)",0,4100);
1855  g10.SetParameters(1,1095,233);
1856  h10->FillRandom("g10",56130);
1857  h10->SetBinContent(85,13425);
1858  TH1F* h20 = new TH1F("h20","gaussian",4096,-.5,4095.5);
1859  g10.SetParameters(1,1673,487);
1860  h20->FillRandom("g10",21370);
1861  h20->SetBinContent(78,17900);
1862  KVHistoManipulator HM2;
1863  HM2.SetVisDebug(); // turn on visual debugging -> create canvas showing rescaling procedure
1864  Double_t params0[10];
1865  // make rescaling using inverse cumulative distribution, limit to x=[100,4095]
1866  TH1* sc30 =HM2.MakeHistoRescaleX(h10,h20,1,params0,"binsnorm",5,"D",100,4095);
1867  return gROOT->GetListOfCanvases()->FindObject("VDCanvas");
1868  }
1869  END_MACRO
1870  */
1871 
1872  TF1* scalefunc = RescaleX(hist1, hist2, degree, params, npoints, direction, xmin, xmax, qmin, qmax, eps);
1873  TString options(opt);
1874  options.ToUpper();
1875  Bool_t norm = options.Contains("NORM");
1876  Bool_t bins = options.Contains("BINS");
1877  Int_t nx = (bins ? hist2->GetNbinsX() : -1);
1878  Double_t nxmin = (bins ? hist2->GetXaxis()->GetXmin() : -1);
1879  Double_t nxmax = (bins ? hist2->GetXaxis()->GetXmax() : -1);
1880  TH1* scalehisto = ScaleHisto(hist1, scalefunc, 0, nx, -1, nxmin, nxmax, -1.0, -1.0, "width");
1881  if (norm) scalehisto->Scale(hist2->Integral("width") / scalehisto->Integral("width"));
1882  if (kVisDebug) {
1883  fVDCanvas->cd(4);
1884  scalehisto->DrawCopy()->SetLineColor(kRed);
1885  hist2->DrawCopy("same")->SetLineColor(kBlack);
1886  gPad->SetLogy(kTRUE);
1887  gPad->Modified();
1888  gPad->Update();
1889  }
1890  delete scalefunc;
1891  return scalehisto;
1892 }
1893 
1894 //###############################################################################################################"
1895 //-------------------------------------------------
1896 
1921 
1922 void KVHistoManipulator::RescaleX(TH1* hist1, TH1* hist2, TF1* scale_func, Int_t npoints,
1923  const Char_t* direction, Double_t xmin, Double_t xmax, Double_t qmin, Double_t qmax, Double_t eps)
1924 {
1925  // Find the transformation of the X-axis of 1-D histogram hist1
1926  // so that the distribution ressembles that of histogram hist2.
1927  // The user provides a function f(x) (TF1* scale_func) which is supposed to
1928  // transform the abscissa X1 of 'hist1' in such a way that P1(f(X1)) = P2(X2).
1929  // We fit 'npoints' comparison points (see below), npoints>=2.
1930  //
1931  // direction : by default ("C") we use the cumulative histogram summed from low x to high x.
1932  // if direction="D", we use the cumulative histogram summed from high x to low x.
1933  // xmin, xmax : range of values of abscissa used to build cumulative histograms. default
1934  // (xmin=xmax=-1) include all values.
1935  // qmin, qmax : minimum & maximum values of cumulative histograms used for the
1936  // comparison (see below). by default qmin=0.05, qmax=0.95.
1937  // METHOD
1938  // ======
1939  // 'hist1' contains the distribution P1 of variable X1
1940  // 'hist2' contains the distribution P2 of variable X2
1941  // Supposing that we can write X2=f(X1), we compare the abscissa of different
1942  // points of the two cumulative distributions, C1(X1) and C2(X2).
1943  // For npoints=2 we compare the 2 X1 & X2 values for which C1=C2=qmin, qmax
1944  // For npoints=3 we compare the 3 X1 & X2 values for which C1=C2=qmin, qmin+(qmax-qmin)/2, qmax
1945  // For npoints=4 we compare the 4 X1 & X2 values for which C1=C2=qmin, qmin+(qmax-qmin)/3, qmin+2*(qmax-qmin)/3, qmax
1946  // etc. etc.
1947  // In each case we fit the 'npoints' couples (X1,X2) with the TF1 pointed to by 'scale_func'
1948 
1949  if (kVisDebug) {
1950  if (!fVDCanvas) fVDCanvas = new TCanvas("VDCanvas", "KVHistoManipulator::RescaleX");
1951  gStyle->SetOptStat("");
1952  fVDCanvas->Clear();
1953  fVDCanvas->Divide(2, 2);
1954  fVDCanvas->cd(1);
1955  hist1->DrawCopy()->SetLineColor(kBlue);
1956  hist2->DrawCopy("same")->SetLineColor(kBlack);
1957  gPad->SetLogy(kTRUE);
1958  gPad->Modified();
1959  gPad->Update();
1960  }
1961  int i;
1962  npoints = TMath::Max(2, npoints);
1963  Info("RescaleX", "Calculating transformation of histo %s using reference histo %s, %d points of comparison",
1964  hist1->GetName(), hist2->GetName(), npoints);
1965  TH1* cum1 = 0;
1966  TH1* cum2 = 0;
1967  if (xmin > -1 && xmax > -1) {
1968  cum1 = CumulatedHisto(hist1, xmin, xmax, direction, "max");
1969  cum2 = CumulatedHisto(hist2, xmin, xmax, direction, "max");
1970  }
1971  else {
1972  cum1 = CumulatedHisto(hist1, direction, -1, -1, "max");
1973  cum2 = CumulatedHisto(hist2, direction, -1, -1, "max");
1974  }
1975  if (kVisDebug) {
1976  fVDCanvas->cd(2);
1977  cum1->DrawCopy()->SetLineColor(kBlue);
1978  cum2->DrawCopy("same")->SetLineColor(kBlack);
1979  gPad->Modified();
1980  gPad->Update();
1981  }
1982  // calculate comparison points
1983  Double_t* quantiles = new Double_t[npoints];
1984  Double_t delta_q = (qmax - qmin) / (1.0 * (npoints - 1));
1985  for (i = 0; i < npoints; i++) quantiles[i] = qmin + i * delta_q;
1986  // get X1 and X2 values corresponding to quantiles
1987  Double_t* X1 = new Double_t[npoints];
1988  Double_t* X2 = new Double_t[npoints];
1989  for (i = 0; i < npoints; i++) {
1990  X1[i] = GetX(cum1, quantiles[i], eps);
1991  X2[i] = GetX(cum2, quantiles[i], eps);
1992  }
1993  for (i = 0; i < npoints; i++) {
1994  printf("COMPARISON: i=%d quantile=%f X1=%f X2=%f\n",
1995  i, quantiles[i], X1[i], X2[i]);
1996  }
1997  // fill TGraph with points to fit
1998  TGraph* fitgraph = new TGraph(npoints, X1, X2);
1999  TString fitoptions = "0N";
2000  if (kVisDebug) fitoptions = "";
2001  if (fitgraph->Fit(scale_func, fitoptions.Data()) != 0) {
2002  Error("RescaleX", "Fitting with function %s failed to converge",
2003  scale_func->GetName());
2004  }
2005  if (kVisDebug) {
2006  fVDCanvas->cd(3);
2007  fitgraph->SetMarkerStyle(20);
2008  gStyle->SetOptStat(1011);
2009  fitgraph->DrawClone("ap");
2010  gPad->Modified();
2011  gPad->Update();
2012  }
2013  delete cum1;
2014  delete cum2;
2015  delete fitgraph;
2016  delete [] quantiles;
2017  delete [] X1;
2018  delete [] X2;
2019 }
2020 
2021 //###############################################################################################################"
2022 //-------------------------------------------------
2023 
2047 
2048 TH1* KVHistoManipulator::MakeHistoRescaleX(TH1* hist1, TH1* hist2, TF1* scale_func, Int_t npoints,
2049  Option_t* opt, const Char_t* direction, Double_t xmin, Double_t xmax, Double_t qmin, Double_t qmax, Double_t eps)
2050 {
2051  // Uses RescaleX(TH1* hist1, TH1* hist2, TF1* scale_func, Int_t npoints, Double_t eps)
2052  // to transform 'hist1' abscissa using TF1 'scale_func' so that its
2053  // distribution resembles that of 'hist2', then generates a new rescaled version of 'hist1'.
2054  //
2055  // npoints = number of points of comparison between the two histograms.
2056  // Make sure this is sufficient for the TF1 used in the transformation.
2057  // i.e. for a polynomial of degree 1 (a+bx), 2 points are enough,
2058  // 3 will give a meaningful Chi^2 value.
2059  // direction : by default ("C") we use the cumulative histogram summed from low x to high x.
2060  // if direction="D", we use the cumulative histogram summed from high x to low x.
2061  // xmin, xmax : range of values of abscissa used to build cumulative histograms. default
2062  // (xmin=xmax=-1) include all values.
2063  // qmin, qmax : minimum & maximum values of cumulative histograms used for the
2064  // comparison (see method RescaleX). by default qmin=0.05, qmax=0.95.
2065  // eps = relative precision used to find comparison points of histos (default = 1.e-07)
2066  //
2067  // OPTIONS
2068  // =======
2069  // opt = "norm" : rescaled histogram normalised to have same integral as hist2
2070  // opt = "bins" : rescaled histogram will have same number of bins & same limits as hist2
2071  // opt = "normbins" |
2072  // opt = "binsnorm" |--> rescaled histogram can be superposed and added to hist2
2073 
2074  RescaleX(hist1, hist2, scale_func, npoints, direction, xmin, xmax, qmin, qmax, eps);
2075  TString options(opt);
2076  options.ToUpper();
2077  Bool_t norm = options.Contains("NORM");
2078  Bool_t bins = options.Contains("BINS");
2079  Int_t nx = (bins ? hist2->GetNbinsX() : -1);
2080  Double_t nxmin = (bins ? hist2->GetXaxis()->GetXmin() : -1);
2081  Double_t nxmax = (bins ? hist2->GetXaxis()->GetXmax() : -1);
2082  TH1* scalehisto = ScaleHisto(hist1, scale_func, 0, nx, -1, nxmin, nxmax, -1.0, -1.0, "width");
2083  if (norm) scalehisto->Scale(hist2->Integral("width") / scalehisto->Integral("width"));
2084  if (kVisDebug) {
2085  fVDCanvas->cd(4);
2086  scalehisto->DrawCopy()->SetLineColor(kRed);
2087  hist2->DrawCopy("same")->SetLineColor(kBlack);
2088  gPad->SetLogy(kTRUE);
2089  gPad->Modified();
2090  gPad->Update();
2091  }
2092  return scalehisto;
2093 }
2094 
2095 //-------------------------------------------------
2096 
2100 
2101 TH1* KVHistoManipulator::CumulatedHisto(TH1* hh, Double_t xmin, Double_t xmax, const TString& direction, Option_t* norm)
2102 {
2103  // Cumule le contenu de l histo hh entre xmin et xmax et retourne l histo correspondant
2104  // Voir CumulatedHisto(TH1* ,TString ,Int_t ,Int_t , Option_t*).
2105  Int_t bmin = hh->FindBin(xmin);
2106  Int_t bmax = hh->FindBin(xmax);
2107  if (bmax > hh->GetNbinsX()) bmax = hh->GetNbinsX();
2108  return CumulatedHisto(hh, direction, bmin, bmax, norm);
2109 }
2110 
2111 
2112 //-------------------------------------------------
2113 
2121 
2123 {
2124  //Camcul du chi2 entre un histogramme et une fonction donnée
2125  //Warning : ne prend en compte que les bins avec une stat>0
2126  //de l histogramme
2127  // norm = kTRUE (default), normalise la valeur du Chi2 au nombre de bin pris en compte
2128  // err = kTRUE (default), prend en compte les erreurs du contenu des bins dans le calcul
2129  // (si celle ci est >0)
2130  Double_t chi2 = 0;
2131  Int_t nbre = 0;
2132  if (h1->InheritsFrom("TH2")) {
2133  Double_t xx[2];
2134  for (Int_t nx = 1; nx < h1->GetNbinsX(); nx += 1)
2135  for (Int_t ny = 1; ny <= h1->GetNbinsY(); ny += 1) {
2136  Double_t hval = h1->GetBinContent(nx, ny);
2137  if (hval > 0) {
2138  if (err) {
2139  Double_t herr = h1->GetBinError(nx, ny);
2140  if (herr > 0) {
2141  nbre += 1;
2142  xx[0] = h1->GetXaxis()->GetBinCenter(nx);
2143  xx[1] = h1->GetYaxis()->GetBinCenter(ny);
2144  Double_t fval = f1->EvalPar(xx, para);
2145  chi2 += TMath::Power((hval - fval) / herr, 2.);
2146  }
2147  }
2148  else {
2149  nbre += 1;
2150  xx[0] = h1->GetXaxis()->GetBinCenter(nx);
2151  xx[1] = h1->GetYaxis()->GetBinCenter(ny);
2152  Double_t fval = f1->EvalPar(xx, para);
2153  chi2 += TMath::Power((hval - fval), 2.);
2154  }
2155  }
2156  }
2157  }
2158  else {
2159  Double_t xx[1];
2160  for (Int_t nx = 1; nx < h1->GetNbinsX(); nx += 1) {
2161  Double_t hval = h1->GetBinContent(nx);
2162  if (hval > 0) {
2163  if (err) {
2164  Double_t herr = h1->GetBinError(nx);
2165  if (herr > 0) {
2166  nbre += 1;
2167  xx[0] = h1->GetXaxis()->GetBinCenter(nx);
2168  Double_t fval = f1->EvalPar(xx, para);
2169  chi2 += TMath::Power((hval - fval) / herr, 2.);
2170  }
2171  }
2172  else {
2173  nbre += 1;
2174  xx[0] = h1->GetXaxis()->GetBinCenter(nx);
2175  Double_t fval = f1->EvalPar(xx, para);
2176  chi2 += TMath::Power((hval - fval), 2.);
2177  }
2178  }
2179  }
2180  }
2181 
2182  if (nbre == 0) {
2183  printf("Warning, KVHistoManipulator::GetChisquare :\n\taucune cellule price en compte dans le calcul du Chi2 ...\n");
2184  return -1;
2185  }
2186  return (norm ? chi2 / nbre : chi2);
2187 
2188 
2189 }
2190 
2191 //-------------------------------------------------
2192 
2200 
2202 {
2203  //Calcul du chi2 entre un histogramme et une fonction donnée
2204  //Warning : ne prend en compte que les bins avec une stat>0
2205  //de l histogramme
2206  // norm = kTRUE (default), normalise la valeur du Chi2 au nombre de bin pris en compte
2207  // err = kTRUE (default), prend en compte les erreurs du contenu des bins dans le calcul
2208  // (si celle ci est >0)
2209  Double_t chi2 = 0;
2210  Int_t nbre = 0;
2211  if (h1->InheritsFrom("TH2")) {
2212  Double_t xx[2];
2213  for (Int_t nx = 1; nx < h1->GetNbinsX(); nx += 1)
2214  for (Int_t ny = 1; ny <= h1->GetNbinsY(); ny += 1) {
2215  Double_t hval = h1->GetBinContent(nx, ny);
2216  if (hval > 0) {
2217  nbre += 1;
2218  xx[0] = h1->GetXaxis()->GetBinCenter(nx);
2219  xx[1] = h1->GetYaxis()->GetBinCenter(ny);
2220  Double_t fval = f1->EvalPar(xx, para);
2221  Double_t logfval = TMath::Log(fval);
2222  chi2 += fval - 1 * hval * logfval;
2223  }
2224  }
2225  }
2226  else {
2227  Double_t xx[1];
2228  for (Int_t nx = 1; nx < h1->GetNbinsX(); nx += 1) {
2229  Double_t hval = h1->GetBinContent(nx);
2230  if (hval > 0) {
2231  nbre += 1;
2232  xx[0] = h1->GetXaxis()->GetBinCenter(nx);
2233  Double_t fval = f1->EvalPar(xx, para);
2234  Double_t logfval = TMath::Log(fval);
2235  chi2 += fval - 1 * hval * logfval;
2236  }
2237  }
2238  }
2239 
2240  if (nbre == 0) {
2241  printf("Warning, KVHistoManipulator::GetChisquare :\n\taucune cellule price en compte dans le calcul du Chi2 ...\n");
2242  return -1;
2243  }
2244  return (norm ? chi2 / nbre : chi2);
2245 
2246 }
2247 
2248 
2249 
2256 
2258 {
2259  // Create and fill a TGraph containing, for each point in G1 and G2,
2260  // the value of (y1/y2).
2261  // G1 & G2 should have the same number of points, with the same x-coordinates
2262  // i.e. x1 = x2 for all points
2263  // if any y2 value = 0, we set the corresponding point's y=0
2264 
2265  Int_t npoints = G1->GetN();
2266  if (G2->GetN() != npoints) {
2267  Error("DivideGraphs", "Graphs must have same number of points");
2268  return nullptr;
2269  }
2270  // make copy of G1
2271  AUTO_NEW_CTOR(TGraph, Gdiv)(*G1);
2272  Gdiv->SetName(Form("%s_divided_by_%s", G1->GetName(), G2->GetName()));
2273  Gdiv->SetTitle(Form("%s divided by %s", G1->GetTitle(), G2->GetTitle()));
2274  Double_t* X = G1->GetX();
2275  Double_t* Y1 = G1->GetY();
2276  Double_t* Y2 = G2->GetY();
2277  for (int i = 0; i < npoints; i++) {
2278  if (Y2[i] != 0) Gdiv->SetPoint(i, X[i], Y1[i] / Y2[i]);
2279  else Gdiv->SetPoint(i, X[i], 0.);
2280  }
2281  return Gdiv;
2282 }
2283 
2284 
2285 
2287 
2289 {
2290  Int_t npoints = g0->GetN();
2291  if (g1->GetN() != npoints) {
2292  Error("ComputeNewGraphFrom", "Graphs must have same number of points %d != %d", npoints, g1->GetN());
2293  return nullptr;
2294  }
2295 
2296  TF1 f1("func_ComputeNewGraphFrom", formula, 0, 1);
2297  if (f1.IsZombie() || f1.GetNpar() != 2) {
2298  Error("ComputeNewGraphFrom", "formula %s for the operation is not valid or has not 2 parameters", formula.Data());
2299  return nullptr;
2300  }
2301 
2302  AUTO_NEW_CTOR(TGraph, gfinal)();
2303  gfinal->SetName(Form("from_%s_%s", g0->GetName(), g1->GetName()));
2304  Double_t* x0 = g0->GetX();
2305  Double_t* y0 = g0->GetY();
2306 
2307  Double_t* x1 = g1->GetX();
2308  Double_t* y1 = g1->GetY();
2309 
2310  for (Int_t ii = 0; ii < npoints; ii++) {
2311  f1.SetParameters(y0[ii], y1[ii]);
2312  if (x1[ii] != x0[ii])
2313  Warning("ComputeNewGraphFrom", "X values are different for the same point %d : %lf %lf", ii, x0[ii], x1[ii]);
2314  Double_t result = f1.Eval(x0[ii]);
2315  gfinal->SetPoint(ii, x0[ii], result);
2316  }
2317  return gfinal;
2318 }
2319 
2320 
2321 
2327 
2328 std::vector<Double_t> KVHistoManipulator::GetLimits(TGraph* G1)
2329 {
2330  /*
2331  xmin -> limits[0];
2332  ymin -> limits[1];
2333  xmax -> limits[2];
2334  ymax -> limits[3];
2335  */
2336  std::vector<Double_t> limits(4);
2337  Double_t xx, yy;
2338  for (Int_t ii = 0; ii < G1->GetN(); ii += 1) {
2339  G1->GetPoint(ii, xx, yy);
2340  if (ii == 0) {
2341  limits[0] = limits[2] = xx;
2342  limits[1] = limits[3] = yy;
2343  }
2344  else {
2345  if (xx < limits[0]) limits[0] = xx;
2346  if (yy < limits[1]) limits[1] = yy;
2347  if (xx > limits[2]) limits[2] = xx;
2348  if (yy > limits[3]) limits[3] = yy;
2349  }
2350  }
2351 
2352  return limits;
2353 }
2354 
2355 
2356 
2362 
2363 std::vector<Double_t> KVHistoManipulator::GetLimits(TMultiGraph* mgr)
2364 {
2365 
2366  /*
2367  xmin -> limits[0];
2368  ymin -> limits[1];
2369  xmax -> limits[2];
2370  ymax -> limits[3];
2371  */
2372  std::vector<Double_t> limits, temp;
2373 
2374  TList* lg = mgr->GetListOfGraphs();
2375  TGraph* gr = nullptr;
2376  for (Int_t ii = 0; ii < lg->GetEntries(); ii += 1) {
2377  gr = dynamic_cast<TGraph*>(lg->At(ii));
2378  if (ii == 0) {
2379  limits = GetLimits(gr);
2380  }
2381  else {
2382  temp = GetLimits(gr);
2383  if (temp[0] < limits[0]) limits[0] = temp[0];
2384  if (temp[1] < limits[1]) limits[1] = temp[1];
2385  if (temp[2] > limits[2]) limits[2] = temp[2];
2386  if (temp[3] > limits[3]) limits[3] = temp[3];
2387  }
2388  }
2389 
2390  return limits;
2391 
2392 }
2393 
2394 
2395 
2401 
2402 std::vector<Double_t> KVHistoManipulator::GetLimits(TProfile* G1)
2403 {
2404  /*
2405  xmin -> limits[0];
2406  ymin -> limits[1];
2407  xmax -> limits[2];
2408  ymax -> limits[3];
2409  */
2410  std::vector<Double_t> limits(4);
2411  Double_t xx, yy;
2412  Bool_t first = kTRUE;
2413  for (Int_t ii = 1; ii <= G1->GetNbinsX(); ii += 1) {
2414  Double_t stat = G1->GetBinEntries(ii);
2415  if (stat > 0) {
2416  xx = G1->GetBinCenter(ii);
2417  yy = G1->GetBinContent(ii);
2418  if (first) {
2419  first = kFALSE;
2420  limits[0] = limits[2] = xx;
2421  limits[1] = limits[3] = yy;
2422  }
2423  else {
2424  if (xx < limits[0]) limits[0] = xx;
2425  if (yy < limits[1]) limits[1] = yy;
2426  if (xx > limits[2]) limits[2] = xx;
2427  if (yy > limits[3]) limits[3] = yy;
2428  }
2429  }
2430  }
2431 
2432  return limits;
2433 
2434 }
2435 
2436 
2437 
2443 
2445 {
2446 
2447  /*
2448  xmin -> limits[0];
2449  ymin -> limits[1];
2450  xmax -> limits[2];
2451  ymax -> limits[3];
2452  */
2453  std::vector<Double_t> temp, limits;
2454 
2455  TProfile* gr = nullptr;
2456  for (Int_t ii = 0; ii < mgr->GetEntries(); ii += 1) {
2457  gr = dynamic_cast<TProfile*>(mgr->At(ii));
2458  if (ii == 0) {
2459  limits = GetLimits(gr);
2460  }
2461  else {
2462  temp = GetLimits(gr);
2463  if (temp[0] < limits[0]) limits[0] = temp[0];
2464  if (temp[1] < limits[1]) limits[1] = temp[1];
2465  if (temp[2] > limits[2]) limits[2] = temp[2];
2466  if (temp[3] > limits[3]) limits[3] = temp[3];
2467  }
2468  }
2469 
2470  return limits;
2471 
2472 }
2473 
2474 
2475 
2479 
2481 {
2482 
2483  //Getthe limits of the histogram in the current pad
2484  //and apply them to the others histogram drawn on the others pads
2485 
2486  if (!gPad) return;
2487  TObject* obj = nullptr;
2488  TVirtualPad* tmp = gPad;
2489  TIter nextp(gPad->GetListOfPrimitives());
2490  TH1* h1 = nullptr;
2491  while ((obj = nextp())) {
2492  if (obj->InheritsFrom("TH1")) {
2493  h1 = dynamic_cast<TH1*>(obj);
2494  }
2495  }
2496  if (h1) {
2497  Int_t x1 = h1->GetXaxis()->GetFirst();
2498  Int_t x2 = h1->GetXaxis()->GetLast();
2499  Double_t y1, y2;
2500  Double_t z1, z2;
2501 
2502  if (h1->GetDimension() == 1) { //TH1 ou TProfile
2503  y1 = h1->GetMinimum();
2504  y2 = h1->GetMaximum();
2505  }
2506  else {
2507  y1 = h1->GetYaxis()->GetFirst();
2508  y2 = h1->GetYaxis()->GetLast();
2509  }
2510 
2511  if (h1->GetDimension() == 2) { //TH2 ou TProfile2D
2512  z1 = h1->GetMinimum();
2513  z2 = h1->GetMaximum();
2514  }
2515  else {
2516  z1 = h1->GetZaxis()->GetFirst();
2517  z2 = h1->GetZaxis()->GetLast();
2518  }
2519 
2520  std::cout << x1 << " " << x2 << " - " << y1 << " " << y2 << " - " << z1 << " " << z2 << std::endl;
2521 
2522  Int_t nc = 1;
2523  TCanvas* cc = gPad->GetCanvas();
2524  TVirtualPad* pad = nullptr;
2525  while ((pad = cc->GetPad(nc))) {
2526  if (tmp != pad) {
2527  pad->cd();
2528  if (AlsoLog) {
2529  gPad->SetLogx(tmp->GetLogx());
2530  gPad->SetLogy(tmp->GetLogy());
2531  gPad->SetLogz(tmp->GetLogz());
2532  }
2533  TIter nextq(gPad->GetListOfPrimitives());
2534  TH1* h1 = nullptr;
2535  while ((obj = nextq())) {
2536  if (obj->InheritsFrom("TH1")) {
2537  h1 = dynamic_cast<TH1*>(obj);
2538 
2539  h1->GetXaxis()->SetRange(x1, x2);
2540  if (h1->GetDimension() == 1) {
2541  h1->SetMinimum(y1);
2542  h1->SetMaximum(y2);
2543  }
2544  else {
2545  h1->GetYaxis()->SetRange(y1, y2);
2546  }
2547  if (h1->GetDimension() == 2) {
2548  h1->SetMinimum(z1);
2549  h1->SetMaximum(z2);
2550  }
2551  else {
2552  h1->GetZaxis()->SetRange(y1, y2);
2553  }
2554  }
2555  }
2556  gPad->Update();
2557  }
2558  nc += 1;
2559  }
2560  cc->Paint();
2561  cc->Update();
2562  }
2563  gPad = tmp;
2564 }
2565 
2566 
int Int_t
KVHistoManipulator * gHistoManipulator
ClassImp(KVPartitionList) void KVPartitionList
Initialisation.
ROOT::R::TRInterface & r
char Char_t
const Bool_t kFALSE
bool Bool_t
double Double_t
const Bool_t kTRUE
const char Option_t
kRed
kBlack
kBlue
#define gDirectory
include TDocParser_001 C image html pict1_TDocParser_001 png width
float xmin
float ymin
float xmax
float ymax
double pow(double, double)
R__EXTERN TRandom * gRandom
char * Form(const char *fmt,...)
R__EXTERN TStyle * gStyle
#define gPad
Toolkit for various operations on histograms & graphs not provided by ROOT.
TF1 * RescaleX(TH1 *hist1, TH1 *hist2, Int_t degree, Double_t *params, Int_t npoints=-1, const Char_t *direction="C", Double_t xmin=-1, Double_t xmax=-1, Double_t qmin=0.05, Double_t qmax=0.95, Double_t eps=1.e-07)
void DefineMarkerStyle(TAttMarker *ob, TString marker)
TH1 * MakeHistoRescaleX(TH1 *hist1, TH1 *hist2, Int_t degree, Double_t *params, Option_t *opt="", Int_t npoints=-1, const Char_t *direction="C", Double_t xmin=-1, Double_t xmax=-1, Double_t qmin=0.05, Double_t qmax=0.95, Double_t eps=1.e-07)
TH2 * CentreeReduiteX(TH2 *hh, Int_t nx=-1, Double_t xmin=-1., Double_t xmax=-1.)
std::vector< Double_t > GetLimits(TGraph *G1)
TGraph * ComputeNewGraphFrom(TGraph *g0, TGraph *g1, const TString &formula)
TH2 * CentreeReduiteY(TH2 *hh, Int_t ny=-1, Double_t ymin=-1., Double_t ymax=-1.)
void DefineStyle(TObject *ob, TString line, TString marker)
Double_t GetX(TH1 *ob, Double_t val, Double_t eps=1.e-07, Int_t nmax=50, Double_t xmin=-1.0, Double_t xmax=-1.0)
Double_t GetChisquare(TH1 *h1, TF1 *f1, Bool_t norm=kTRUE, Bool_t err=kTRUE, Double_t *para=nullptr)
KVNumberList * Saucisson(TH1 *hh, Int_t ntranches=10)
TGraph * ScaleGraph(TGraph *hh, TF1 *fx, TF1 *fy)
TGraph * LinkGraphs(TGraph *grx, TGraph *gry)
void DefinePattern(TH1 *ob, TString titleX="42 0.08 0.8", TString titleY="42 0.07 1.2", TString labelX="42 0.05 0.005", TString labelY="42 0.05 0.006")
virtual ~KVHistoManipulator(void)
TH2 * RenormaliseHisto(TH2 *hh, Int_t bmin=-1, Int_t bmax=-1, TString axis="X", Double_t valref=1)
Int_t CutStatBin(TH1 *hh, Int_t stat_min=-1, Int_t stat_max=-1)
TGraphErrors * GetMomentEvolution(TH2 *hh, TString momentx, TString momenty, TString axis="Y", Double_t stat_min=0)
void DefineTitle(TH1 *ob, TString xtit, TString ytit)
TGraph * ExtractMeanAndSigmaFromProfile(TProfile *pf, TGraph *&sigma)
void ApplyCurrentLimitsToAllCanvas(Bool_t AlsoLog=kFALSE)
TH1 * CentreeReduite(TH1 *hh, Int_t nx=-1, Int_t ny=-1, Double_t xmin=-1., Double_t xmax=-1., Double_t ymin=-1., Double_t ymax=-1.)
TH1 * ScaleHisto(TH1 *hh, TF1 *fx, TF1 *fy=NULL, Int_t nx=-1, Int_t ny=-1, Double_t xmin=-1., Double_t xmax=-1., Double_t ymin=-1., Double_t ymax=-1., Option_t *norm="")
TGraph * DivideGraphs(TGraph *G1, TGraph *G2)
Double_t GetLikelihood(TH1 *h1, TF1 *f1, Bool_t norm=kTRUE, Double_t *para=nullptr)
KVList * Give_ProjectionList(TH2 *hh, Double_t MinIntegral=-1, TString axis="X")
TGraphErrors * MakeGraphFrom(TProfile *pf, Bool_t Error=kTRUE)
TH1 * CumulatedHisto(TH1 *hh, TString direction="C", Int_t bmin=-1, Int_t bmax=-1, Option_t *norm="surf")
TH1 * GetDerivative(TH1 *hh, Int_t order)
void DefineLineStyle(TAttLine *ob, TString line)
Double_t GetCorrelationFactor(TH2 *hh)
Int_t Apply_TCutG(TH2 *hh, TCutG *cut, TString mode="in")
Extended TList class which owns its objects by default.
Definition: KVList.h:27
Strings used to represent a set of ranges of values.
Definition: KVNumberList.h:83
Bool_t End(void) const
Definition: KVNumberList.h:196
void Begin(void) const
void Add(Int_t)
Add value 'n' to the list.
Int_t Next(void) const
virtual void Add(TObject *obj)
virtual void SetTitleOffset(Float_t offset=1)
virtual void SetLabelSize(Float_t size=0.04)
virtual void SetTitleFont(Style_t font=62)
virtual void SetLabelOffset(Float_t offset=0.005)
virtual void SetLabelFont(Style_t font=62)
virtual void SetTitleSize(Float_t size=0.04)
virtual void SetLineStyle(Style_t lstyle)
virtual void SetLineWidth(Width_t lwidth)
virtual void SetLineColor(Color_t lcolor)
virtual void SetMarkerColor(Color_t mcolor=1)
virtual void SetMarkerStyle(Style_t mstyle=1)
virtual void SetMarkerSize(Size_t msize=1)
virtual Int_t FindBin(const char *label)
virtual Double_t GetBinCenter(Int_t bin) const
Double_t GetXmax() const
virtual Double_t GetBinLowEdge(Int_t bin) const
Int_t GetLast() const
Double_t GetXmin() const
Int_t GetNbins() const
virtual void SetRange(Int_t first=0, Int_t last=0)
virtual Double_t GetBinWidth(Int_t bin) const
virtual Double_t GetBinUpEdge(Int_t bin) const
Int_t GetFirst() const
void * New(ENewType defConstructor=kClassNew, Bool_t quiet=kFALSE) const
static TClass * GetClass(Bool_t load=kTRUE, Bool_t silent=kFALSE)
virtual Int_t GetEntries() const
virtual Int_t GetNDF() const
virtual Double_t GetParameter(const TString &name) const
TAxis * GetYaxis() const
Double_t GetChisquare() const
virtual void SetRange(Double_t xmin, Double_t xmax)
virtual Int_t GetNpar() const
virtual Double_t Derivative(Double_t x, Double_t *params=0, Double_t epsilon=0.001) const
virtual void SetNpx(Int_t npx=100)
virtual Double_t EvalPar(const Double_t *x, const Double_t *params=0)
virtual void SetParameters(const Double_t *params)
virtual Double_t Eval(Double_t x, Double_t y=0, Double_t z=0, Double_t t=0) const
TAxis * GetXaxis() const
virtual void SetPointError(Double_t ex, Double_t ey)
virtual Int_t IsInside(Double_t x, Double_t y) const
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
virtual Double_t * GetEX() const
virtual Double_t * GetEY() const
virtual void SetNameTitle(const char *name="", const char *title="")
Int_t GetN() const
virtual TFitResultPtr Fit(const char *formula, Option_t *option="", Option_t *goption="", Axis_t xmin=0, Axis_t xmax=0)
Double_t * GetX() const
TAxis * GetXaxis() const
TAxis * GetYaxis() const
Double_t * GetY() const
virtual Int_t GetPoint(Int_t i, Double_t &x, Double_t &y) const
virtual void SetTitle(const char *title)
virtual Double_t GetBinCenter(Int_t bin) const
virtual Int_t GetNbinsY() const
virtual Double_t GetBinError(Int_t bin) const
TAxis * GetXaxis()
virtual Double_t GetMean(Int_t axis=1) const
virtual Int_t GetDimension() const
virtual void Reset(Option_t *option="")
TAxis * GetYaxis()
TObject * Clone(const char *newname=0) const
virtual Double_t GetMaximum(Double_t maxval=FLT_MAX) const
virtual Int_t GetNbinsX() const
virtual void SetMaximum(Double_t maximum=-1111)
virtual void SetBinError(Int_t bin, Double_t error)
virtual Int_t Fill(const char *name, Double_t w)
virtual void SetMinimum(Double_t minimum=-1111)
Double_t GetRMS(Int_t axis=1) const
virtual void SetBinContent(Int_t bin, Double_t content)
virtual Double_t GetBinLowEdge(Int_t bin) const
virtual TH1 * DrawCopy(Option_t *option="", const char *name_postfix="_copy") const
virtual Double_t Integral(Int_t binx1, Int_t binx2, Option_t *option="") const
virtual void SetName(const char *name)
virtual Double_t GetBinContent(Int_t bin) const
virtual void SetBins(Int_t nx, const Double_t *xBins)
virtual void SetNameTitle(const char *name, const char *title)
virtual Double_t GetBinWidth(Int_t bin) const
virtual void Scale(Double_t c1=1, Option_t *option="")
virtual Int_t FindBin(Double_t x, Double_t y=0, Double_t z=0)
TAxis * GetZaxis()
virtual Double_t GetMinimum(Double_t minval=-FLT_MAX) const
TH1D * ProjectionY(const char *name="_py", Int_t firstxbin=0, Int_t lastxbin=-1, Option_t *option="") const
virtual void Reset(Option_t *option="")
TH1D * ProjectionX(const char *name="_px", Int_t firstybin=0, Int_t lastybin=-1, Option_t *option="") const
virtual Double_t GetBinContent(Int_t bin) const
virtual void SetBinContent(Int_t bin, Double_t content)
virtual TObject * At(Int_t idx) const
void InitWithPrototype(const char *function, const char *proto, ROOT::EFunctionMatchMode mode=ROOT::kConversionMatch)
Bool_t IsValid() const
void Execute()
TList * GetListOfGraphs() const
virtual const char * GetName() const
virtual void SetTitle(const char *title="")
virtual const char * GetTitle() const
virtual void SetName(const char *name)
Int_t GetEntries() const
TObject * At(Int_t idx) const
virtual const char * ClassName() const
virtual TObject * DrawClone(Option_t *option="") const
R__ALWAYS_INLINE Bool_t IsZombie() const
virtual Bool_t InheritsFrom(const char *classname) const
virtual void SetBinEntries(Int_t bin, Double_t w)
virtual Double_t GetBinEntries(Int_t bin) const
virtual Double_t GetBinError(Int_t bin) const
virtual Double_t GetBinContent(Int_t bin) const
virtual Double_t Uniform(Double_t x1, Double_t x2)
virtual TObject * At(Int_t idx) const=0
void ToUpper()
TObjArray * Tokenize(const TString &delim) const
const char * Data() const
Bool_t IsNull() const
void Form(const char *fmt,...)
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
void SetOptStat(Int_t stat=1)
virtual Int_t GetLogz() const=0
virtual TVirtualPad * cd(Int_t subpadnumber=0)=0
virtual Int_t GetLogy() const=0
virtual Int_t GetLogx() const=0
TLine * line
VecExpr< UnaryOp< Sqrt< T >, SVector< T, D >, T >, T, D > sqrt(const SVector< T, D > &rhs)
VecExpr< UnaryOp< Fabs< T >, SVector< T, D >, T >, T, D > fabs(const SVector< T, D > &rhs)
const Double_t sigma
const Int_t n
Double_t ey[n]
TGraphErrors * gr
Double_t ex[n]
TH1F * h1
TF1 * f1
TH1 * h
const long double s
Definition: KVUnits.h:94
const long double cc
volumes
Definition: KVUnits.h:83
void Info(const char *location, const char *va_(fmt),...)
void Error(const char *location, const char *va_(fmt),...)
void Warning(const char *location, const char *va_(fmt),...)
Double_t Power(Double_t x, Double_t y)
Double_t Log(Double_t x)
Double_t Sqrt(Double_t x)
Double_t Abs(Double_t d)
Double_t Max(Double_t a, Double_t b)
void spline(double x[], double y[], int n, double yp1, double ypn, double *y2)