HARPO  5.1.1
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MakeNiceHisto.cxx
Go to the documentation of this file.
1 #include <iostream>
2 
3 #include "MakeNiceHisto.h"
4 
5 #include "TStyle.h"
6 #ifndef ROOT_TFile
7 #include "TFile.h"
8 #endif
9 #ifndef ROOT_TROOT
10 #include "TROOT.h"
11 #endif
12 #ifndef ROOT_TMath
13 #include "TMath.h"
14 #endif
15 #ifndef ROOT_TH1
16 #include "TH1.h"
17 #endif
18 #ifndef ROOT_TH2
19 #include "TH2.h"
20 #endif
21 #ifndef ROOT_TVirtualPad
22 #include "TVirtualPad.h"
23 #endif
24 #ifndef ROOT_TCanvas
25 #include "TCanvas.h"
26 #endif
27 #ifndef ROOT_TF1
28 #include "TF1.h"
29 #endif
30 #ifndef ROOT_TObject
31 #include "TObject.h"
32 #endif
33 // #ifndef ROOT_TMinuit
34 // #include "TMinuit.h"
35 // #endif
36 #ifndef ROOT_TGraphErrors
37 #include "TGraphErrors.h"
38 #endif
39 #ifndef ROOT_TError
40 #include "TError.h"
41 #endif
42 
43 #ifndef __IOSTREAM__
44 #include <iostream>
45 #endif
46 
47 #include "RVersion.h"
48 
49 using std::cout;
50 using std::endl;
51 
52 
53 // TH1F* HistInvert(TH1* hist);
54 // const Double_t *MakeLogBinning(Int_t nbins, Double_t xmin, Double_t xmax);
55 // void MakeNiceHisto(TH1* hist, TVirtualPad* c1 = 0, const char* opt = "");
56 // void MakeNice1dHisto(TH1* hist, TVirtualPad* c1 = 0, const char* opt = "");
57 // void MakeNice2dHisto(TH2* hist, TVirtualPad* c1 = 0, const char* opt = "");
58 // TCanvas* FindCanvas(const char* canvasName,
59 // Int_t xwidth = 700, Int_t ywidth = 500);
60 // void IntegralWithError(TF1* func, Double_t &integral, Double_t &error,
61 // Double_t xlow, Double_t xhigh,
62 // Double_t stepsize = 0.000001);
63 // Double_t IntegralDerivative(TF1* func, Double_t xlow, Double_t xhigh,
64 // Double_t stepsize, Int_t parameter,
65 // Double_t epsilon);
66 // void FunctionWithError(TF1* func, Double_t x, Double_t &value,
67 // Double_t &error);
68 // Double_t FunctionDerivative(TF1* func, Double_t x,
69 // Int_t parameter, Double_t epsilon);
70 // void SetHistMaximum(TH1* hist, Float_t factor);
71 // TH1F* GetHistogramFromGraph(TGraphErrors* graph, const char* histName);
72 // TGraph* GetGraphFromHistogram(TH1* hist);
73 
74 //_______________________________________________________________________
75 TH1F* HistInvert(TH1* hist)
76 {
77  TH1F* histNew = new TH1F(Form("%s_inv", hist->GetName()),
78  Form("%s (inverted)", hist->GetTitle()),
79  hist->GetNbinsX(),
80  hist->GetXaxis()->GetXmin(),
81  hist->GetXaxis()->GetXmax());
82 
83  const Int_t nBins = hist->GetNbinsX();
84 
85  for(Int_t i = 1; i <= nBins; i++) {
86 
87  if(hist->GetBinContent(i) == 0)
88  continue;
89 
90  histNew->SetBinContent(i, 1.0/hist->GetBinContent(i));
91  histNew->SetBinError(i, hist->GetBinError(i)/hist->GetBinContent(i)/hist->GetBinContent(i));
92  }
93 
94  return histNew;
95 }
96 
97 
98 //_______________________________________________________________________
99 TH1F* HistLog(const char* name, const char* title, Int_t nBins, Double_t xMin, Double_t xMax)
100 {
101  if(xMin==0) xMin = xMax*1e-6;
102  const Int_t nTot = nBins+1;
103  Double_t fac = TMath::Power(xMax/xMin,1./(nTot));
104  Double_t xbins[nTot];
105  //cout << "nBins = " << nBins << ", nTot = " << nTot << endl;
106  for (Int_t i=0;i<nTot;i++){
107  xbins[i] = xMin;
108  //cout << "xbins[" << i << "] = " << xMin << endl;
109  xMin *= fac;
110  // cout << "xbins[" << i << "] = " << xbins[i] << endl;
111  }
112 
113  TH1F* h = new TH1F(name,title,nBins,xbins);
114 
115  return h;
116 }
117 
118 
119 //_______________________________________________________________________
120 TH2F* HistLogLog(const char* name, const char* title, Int_t nBinsX, Double_t xMinX, Double_t xMaxX, Int_t nBinsY, Double_t xMinY, Double_t xMaxY)
121 {
122  if(xMinX==0) xMinX = xMaxX*1e-6;
123  Double_t facX = TMath::Power(xMaxX/xMinX,1./(nBinsX+1));
124  Double_t xbinsX[nBinsX+1];
125  for (Int_t i=0;i<=nBinsX;i++){
126  xbinsX[i] = xMinX;
127  xMinX *= facX;
128  // cout << "xbinsX[" << i << "] = " << xbinsX[i] << endl;
129  }
130 
131  if(xMinY==0) xMinY = xMaxY*1e-6;
132  Double_t facY = TMath::Power(xMaxY/xMinY,1./(nBinsY+1));
133  Double_t xbinsY[nBinsY+1];
134  for (Int_t i=0;i<=nBinsY;i++){
135  xbinsY[i] = xMinY;
136  xMinY *= facY;
137  // cout << "xbinsY[" << i << "] = " << xbinsY[i] << endl;
138  }
139 
140  TH2F* h = new TH2F(name,title,nBinsX,xbinsX,nBinsY,xbinsY);
141 
142  return h;
143 }
144 
145 
146 //_______________________________________________________________________
147 TH2F* HistLinLog(const char* name, const char* title, Int_t nBinsX, Double_t xMinX, Double_t xMaxX, Int_t nBinsY, Double_t xMinY, Double_t xMaxY)
148 {
149  if(xMinY==0) xMinY = xMaxY*1e-6;
150  Double_t fac = TMath::Power(xMaxY/xMinY,1./(nBinsY+1));
151  Double_t xbinsY[nBinsY+1];
152  for (Int_t i=0;i<=nBinsY;i++){
153  xbinsY[i] = xMinY;
154  xMinY *= fac;
155  // cout << "xbinsY[" << i << "] = " << xbinsY[i] << endl;
156  }
157 
158  TH2F* h = new TH2F(name,title,nBinsX,xMinX,xMaxX,nBinsY,xbinsY);
159 
160  return h;
161 }
162 
163 //_______________________________________________________________________
164 TH2F* HistLogLin(const char* name, const char* title, Int_t nBinsX, Double_t xMinX, Double_t xMaxX, Int_t nBinsY, Double_t xMinY, Double_t xMaxY)
165 {
166  if(xMinX==0) xMinX = xMaxX*1e-6;
167  Double_t fac = TMath::Power(xMaxX/xMinX,1./(nBinsX+1));
168  Double_t xbinsX[nBinsX+1];
169  for (Int_t i=0;i<=nBinsX;i++){
170  xbinsX[i] = xMinX;
171  xMinX *= fac;
172  // cout << "xbinsX[" << i << "] = " << xbinsX[i] << endl;
173  }
174 
175  TH2F* h = new TH2F(name,title,nBinsX,xbinsX,nBinsY,xMinY,xMaxY);
176 
177  return h;
178 }
179 
180 
181 
182 //_______________________________________________________________________
183 void SetHistMaximum(TH1* hist, Float_t factor)
184 {
185  hist->SetMaximum(hist->GetBinContent(hist->GetMaximumBin())*factor);
186 }
187 
188 //_______________________________________________________________________
189 void SetTitles(TH1* hist, const char* title, const char* xtitle, const char* ytitle)
190 {
191  hist->SetTitle(title);
192  hist->GetXaxis()->SetTitle(xtitle);
193  hist->GetYaxis()->SetTitle(ytitle);
194 }
195 
196 
197 //_______________________________________________________________________
198 TH1F* MakeNiceGraph(TGraph* g, TVirtualPad* c1, const char* opt)
199 {
200 
201  Double_t *x = g->GetX();
202  Double_t *y = g->GetY();
203  Int_t n = g->GetN();
204  if(n<=1) return 0;
205  Double_t xmin = TMath::MinElement(n,x);
206  Double_t xmax = TMath::MaxElement(n,x);
207  Double_t ymin = TMath::MinElement(n,y);
208  Double_t ymax = TMath::MaxElement(n,y);
209 
210  TH1F* h = new TH1F(Form("%s_hist",g->GetName()),"",100,xmin,xmax);
211  h->SetMinimum(ymin);
212  h->SetMaximum(ymax);
213  MakeNiceHisto(h,c1,"",0);
214  g->Draw(Form("%s same",opt));
215 
216  return h;
217 
218 }
219 
220 //_______________________________________________________________________
221 TVirtualPad* MakeNiceHisto(TH1* hist, TVirtualPad* c1, const char* opt, Bool_t copy)
222 {
223  if(c1==0)
224  c1 = new TCanvas();
225 
226  Bool_t is2dHist = hist->InheritsFrom("TH2");
227 
228  SetHistStyle(hist);
229 
230  if(is2dHist)
231  MakeNice2dHisto((TH2*)hist, c1, opt, copy);
232  else
233  MakeNice1dHisto(hist, c1, opt, copy);
234 
235  return c1;
236 }
237 
238 
239 
240 void SetHistStyle(TH1* hist)
241 {
242 
243  gStyle->SetTitleH(0.08);
244  gStyle->SetTitleW(0.78);
245  gStyle->SetTitleX(0.16);
246  gStyle->SetTitleY(0.99);
247  gStyle->SetTitleFont(132," ");
248  gStyle->SetTitleSize(0.09," ");
249 
250  hist->GetXaxis()->SetTitleSize(0.07);
251  hist->GetYaxis()->SetTitleSize(0.07);
252  hist->GetZaxis()->SetTitleSize(0.07);
253  // hist->GetXaxis()->SetTitleOffset(0.9);
254  hist->GetXaxis()->SetTitleOffset(1.1);
255  // hist->GetYaxis()->SetTitleOffset(0.9);
256  hist->GetYaxis()->SetTitleOffset(1.2);
257  hist->GetZaxis()->SetTitleOffset(0.9);
258  hist->GetXaxis()->SetTitleFont(132);
259  hist->GetYaxis()->SetTitleFont(132);
260  hist->GetZaxis()->SetTitleFont(132);
261  hist->GetXaxis()->SetLabelFont(132);
262  hist->GetYaxis()->SetLabelFont(132);
263  hist->GetZaxis()->SetLabelFont(132);
264  hist->GetXaxis()->SetLabelSize(0.06);
265  hist->GetYaxis()->SetLabelSize(0.06);
266  hist->GetZaxis()->SetLabelSize(0.06);
267 
268  hist->GetXaxis()->SetLabelOffset(0.01);
269  hist->GetYaxis()->SetLabelOffset(0.015);
270  hist->GetZaxis()->SetLabelOffset(0.015);
271 
272 
273 }
274 
275 
276 //_______________________________________________________________________
277 void MakeNice1dHisto(TH1* hist, TVirtualPad* c1, const char* opt, Bool_t copy)
278 {
279  Bool_t hasTitle = kTRUE;
280  const char* title = hist->GetTitle();
281  if(title && title[0] == '\0')
282  hasTitle = kFALSE;
283 
284  if(c1) {
285 
286  // c1->SetLeftMargin(0.15);
287  c1->SetLeftMargin(0.18);
288  c1->SetRightMargin(0.03);
289  // c1->SetBottomMargin(0.15);
290  c1->SetBottomMargin(0.17);
291  if(hasTitle)
292  c1->SetTopMargin(0.09);
293  else
294  c1->SetTopMargin(0.03);
295  c1->cd();
296  if(copy)
297  hist->DrawCopy(opt);
298  else
299  hist->Draw(opt);
300  }
301 }
302 
303 //_______________________________________________________________________
304 void MakeNice2dHisto(TH2* hist, TVirtualPad* c1, const char* option, Bool_t copy)
305 {
306 
307  TString opt(option);
308  opt.ToUpper();
309 
310  Bool_t colz = kFALSE;
311  if (opt.Contains("COLZ"))
312  colz = kTRUE;
313 
314  Bool_t hasTitle = kTRUE;
315  const char* title = hist->GetTitle();
316  if(title && title[0] == '\0')
317  hasTitle = kFALSE;
318 
319  if(colz)
320  gStyle->SetTitleW(0.68);
321 
322 
323  if(c1) {
324 
325  // c1->SetLeftMargin(0.15);
326  c1->SetLeftMargin(0.18);
327  if(!colz)
328  c1->SetRightMargin(0.03);
329  else{
330  if(strcmp(hist->GetZaxis()->GetTitle(),"")){
331  c1->SetRightMargin(0.17);
332  }
333  else
334  c1->SetRightMargin(0.14);
335  }
336  //c1->SetBottomMargin(0.15);
337  c1->SetBottomMargin(0.17);
338  if(hasTitle)
339  c1->SetTopMargin(0.09);
340  else
341  c1->SetTopMargin(0.03);
342  c1->cd();
343  if(copy)
344  hist->DrawCopy(option);
345  else
346  hist->Draw(option);
347  }
348 }
349 
350 
351 
352 //_______________________________________________________________________
353 void MakeNiceHistoPad(TH1* hist, TVirtualPad* c1, Double_t scaleX, Double_t scaleY, const char* opt, Bool_t copy)
354 {
355  if(c1==0)
356  return;
357  // c1 = new TPad();
358 
359  Bool_t is2dHist = hist->InheritsFrom("TH2");
360 
361  SetHistStyle(hist);
362 
363  if(is2dHist)
364  MakeNice2dHistoPad((TH2*)hist, c1, scaleX, scaleY, opt, copy);
365  else
366  MakeNice1dHistoPad(hist, c1, scaleX, scaleY, opt, copy);
367 }
368 
369 
370 //_______________________________________________________________________
371 void MakeNiceHistoPad(TH1* hist, TVirtualPad* c, Int_t i, Int_t j, Double_t scaleX, Double_t scaleY, const char* opt, Bool_t copy)
372 {
373  if(c==0)
374  return;
375  // c1 = new TPad();
376  Char_t name[64];
377  sprintf(name,"%s_%d_%d",c->GetName(),i,j);
378  cout << name << endl;
379  TPad* c1 = ((TPad*) gROOT->FindObject(name));
380  //c->cd(0);
381  //c1->Draw();
382  Int_t w = c->GetWw();
383  Double_t wPad = c1->GetWNDC();
384  Int_t h = c->GetWh();
385  Double_t hPad = c1->GetHNDC();
386 
387  cout << w << " " << wPad << " " << h << " " << hPad << endl;
388 
389  // scaleY = 1;
390  // scaleX = w*1./h;
391 
392  // cout << scaleX << " " << scaleY << endl;
393 
394  Bool_t is2dHist = hist->InheritsFrom("TH2");
395 
396  SetHistStyle(hist);
397  // hist->GetXaxis()->SetTitleOffset(0.5);
398  // hist->GetYaxis()->SetTitleOffset(0.5);
399 
400  if(is2dHist)
401  MakeNice2dHistoPad((TH2*)hist, c1, scaleX, scaleY, opt, copy);
402  else
403  MakeNice1dHistoPad(hist, c1, scaleX, scaleY, opt, copy);
404 }
405 
406 
407 
408 
409 //_______________________________________________________________________
410 void MakeNice1dHistoPad(TH1* hist, TVirtualPad* c1, Double_t scaleX, Double_t scaleY, const char* opt, Bool_t copy)
411 {
412  if(!c1) return;
413 
414  cout << "MakeNice1dHistoPad " << c1->GetName() << endl;
415 
416  // Bool_t hasTitle = kTRUE;
417  // const char* title = hist->GetTitle();
418  // if(title && title[0] == '\0')
419  // hasTitle = kFALSE;
420 
421  Double_t ratioX = c1->GetWNDC()*scaleX;
422  Double_t ratioY = c1->GetHNDC()*scaleY;
423 
424  // std::cout << ratioX << " " << ratioY << std::endl;
425 
426  hist->GetXaxis()->SetLabelSize(hist->GetXaxis()->GetLabelSize()/ratioY);
427  hist->GetYaxis()->SetLabelSize(hist->GetYaxis()->GetLabelSize()/ratioY);
428  hist->GetZaxis()->SetLabelSize(hist->GetZaxis()->GetLabelSize()/ratioY);
429 
430  hist->GetXaxis()->SetTitleSize(hist->GetXaxis()->GetTitleSize()/ratioY);
431  hist->GetYaxis()->SetTitleSize(hist->GetYaxis()->GetTitleSize()/ratioY);
432  hist->GetZaxis()->SetTitleSize(hist->GetZaxis()->GetTitleSize()/ratioY);
433 
434  hist->GetXaxis()->SetLabelOffset(hist->GetXaxis()->GetLabelOffset()*ratioX);
435  hist->GetYaxis()->SetLabelOffset(hist->GetYaxis()->GetLabelOffset()*ratioY);
436  hist->GetZaxis()->SetLabelOffset(hist->GetZaxis()->GetLabelOffset()*ratioY);
437 
438  hist->GetXaxis()->SetTitleOffset(hist->GetXaxis()->GetTitleOffset()*ratioX);
439  hist->GetYaxis()->SetTitleOffset(hist->GetYaxis()->GetTitleOffset()*ratioY);
440  hist->GetZaxis()->SetTitleOffset(hist->GetZaxis()->GetTitleOffset()*ratioY);
441 
442 
443  c1->cd();
444  if(copy)
445  hist->DrawCopy(opt);
446  else
447  hist->Draw(opt);
448 }
449 
450 //_______________________________________________________________________
451 void MakeNice2dHistoPad(TH2* hist, TVirtualPad* c1, Double_t scaleX, Double_t scaleY, const char* option, Bool_t copy)
452 {
453  if(!c1) return;
454 
455  TString opt(option);
456  opt.ToUpper();
457 
458  // Bool_t colz = kFALSE;
459  // if (opt.Contains("COLZ"))
460  // colz = kTRUE;
461  // if(colz)
462  // gStyle->SetTitleW(0.68);
463 
464 
465  // Bool_t hasTitle = kTRUE;
466  // const char* title = hist->GetTitle();
467  // if(title && title[0] == '\0')
468  // hasTitle = kFALSE;
469 
470  // Double_t ratioX = c1->GetWNDC()/c1->GetCanvas()->GetWNDC();
471  // Double_t ratioY = c1->GetHNDC()/c1->GetCanvas()->GetHNDC();
472  // Double_t ratioX = (c1->GetWNDC() - c1->GetLeftMargin() - c1->GetRightMargin())/c1->GetWNDC();
473  // Double_t ratioY = (c1->GetHNDC() - c1->GetBottomMargin() - c1->GetTopMargin())/c1->GetHNDC();
474  Double_t ratioX = c1->GetWNDC()*scaleX;
475  Double_t ratioY = c1->GetHNDC()*scaleY;
476  // std::cout << c1->GetHNDC() << " - " << c1->GetBottomMargin() << " - " << c1->GetTopMargin() << std::endl;
477  // std::cout << "MakeNice2dHistoPad: " << ratioX << " " << ratioY << std::endl;
478 
479  hist->GetXaxis()->SetLabelSize(hist->GetXaxis()->GetLabelSize()/ratioX);
480  hist->GetYaxis()->SetLabelSize(hist->GetYaxis()->GetLabelSize()/ratioY);
481  hist->GetZaxis()->SetLabelSize(hist->GetZaxis()->GetLabelSize()/ratioY);
482 
483  hist->GetXaxis()->SetTitleSize(hist->GetXaxis()->GetTitleSize()/ratioX);
484  hist->GetYaxis()->SetTitleSize(hist->GetYaxis()->GetTitleSize()/ratioY);
485  hist->GetZaxis()->SetTitleSize(hist->GetZaxis()->GetTitleSize()/ratioY);
486 
487  hist->GetXaxis()->SetLabelOffset(hist->GetXaxis()->GetLabelOffset()*ratioX);
488  hist->GetYaxis()->SetLabelOffset(hist->GetYaxis()->GetLabelOffset()*ratioY);
489  hist->GetZaxis()->SetLabelOffset(hist->GetZaxis()->GetLabelOffset()*ratioY);
490 
491  hist->GetXaxis()->SetTitleOffset(hist->GetXaxis()->GetTitleOffset()*ratioX);
492  hist->GetYaxis()->SetTitleOffset(hist->GetYaxis()->GetTitleOffset()*ratioY);
493  hist->GetZaxis()->SetTitleOffset(hist->GetZaxis()->GetTitleOffset()*ratioY);
494 
495 
496  // if(c1) {
497 
498  // c1->SetLeftMargin(0.12);
499  // if(!colz)
500  // c1->SetRightMargin(0.03);
501  // else{
502  // if(strcmp(hist->GetZaxis()->GetTitle(),"")){
503  // c1->SetRightMargin(0.17);
504  // }
505  // else
506  // c1->SetRightMargin(0.14);
507  // }
508  // c1->SetBottomMargin(0.15);
509  // if(hasTitle)
510  // c1->SetTopMargin(0.09);
511  // else
512  // c1->SetTopMargin(0.03);
513  // c1->cd();
514  if(copy)
515  hist->DrawCopy(option);
516  else
517  hist->Draw(option);
518  // }
519 }
520 
521 
522 //______________________________________________________________________
523 const Double_t *MakeLogBinning(Int_t nbins, Double_t xmin, Double_t xmax)
524 {
525 
526  if(xmin<=0){
527  cout << "Log scale cannot start from 0 or below" << endl;
528  return 0;
529  }
530 
531  Double_t logxmin = TMath::Log(xmin);
532  Double_t logxmax = TMath::Log(xmax);
533  Double_t binwidth = (logxmax-logxmin)/nbins;
534  const Int_t Nbins = nbins+1;
535  Double_t *xbins;
536  xbins = (Double_t *)malloc(Nbins*sizeof(*xbins));
537  xbins[0] = xmin;
538  for (Int_t i = 1; i<=nbins; i++){
539  xbins[i] = TMath::Exp(logxmin+i*binwidth);
540  cout << xbins[i] << endl;
541  }
542 
543  return xbins;
544 
545 }
546 
547 
548 
549 
550 //______________________________________________________________________
551 TCanvas* FindCanvas(const char* canvasName,
552  Int_t xwidth, Int_t ywidth)
553 {
554 
555  TCanvas *c1 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject(canvasName);
556  if(!c1)
557  c1 = new TCanvas(canvasName, canvasName, xwidth, ywidth);
558  else
559  c1->SetWindowSize(xwidth, ywidth);
560 
561  return c1;
562 }
563 
564 #ifdef ROOT_TMinuit
565 //______________________________________________________________________
566 void IntegralWithError(TF1* func, Double_t &integral, Double_t &error,
567  Double_t xlow, Double_t xhigh,
568  Double_t stepsize)
569 {
570  // set range for derivative method
571 
572  // maybe check that parameter is not fixed????
573  func->SetRange(xlow, xhigh);
574  const Double_t* dummy = 0;
575 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,0,0)
576  integral = func->Integral(xlow, xhigh, stepsize);
577 #else
578  integral = func->Integral(xlow, xhigh, dummy, stepsize);
579 #endif
580 
581  const Int_t npar = func->GetNpar();
582 
583  Double_t covariance[npar][npar];
584  gMinuit->mnemat(&covariance[0][0], func->GetNpar());
585 
586  Double_t partialDiff[npar];
587 
588  for(Int_t i = 0; i < npar; i++) {
589 
590  Double_t par_err = TMath::Sqrt(covariance[i][i]);
591 // if(par_err < 0.0) {
592 
593 // func->Error("IntegralWithError", "Parameter error < 0!!!");
594 // return;
595 // }
596 
597  partialDiff[i] = IntegralDerivative(func, xlow, xhigh, stepsize,
598  i, 0.1*par_err);
599  }
600 
601  Double_t sigma2 = 0;
602 
603  for(Int_t i = 0; i < npar; i++) {
604  for(Int_t j = 0; j < npar; j++) {
605 
606  cout << "sigma2(" << i << ", " << j << ") = "
607  << partialDiff[i]*partialDiff[j]*covariance[i][j]
608  << endl;
609 
610  sigma2+=partialDiff[i]*partialDiff[j]*covariance[i][j];
611  }
612  }
613 
614  error = TMath::Sqrt(sigma2);
615 }
616 
617 //______________________________________________________________________
618 void FunctionWithError(TF1* func, Double_t x, Double_t &value, Double_t &error)
619 {
620  value = func->Eval(x);
621  cout << "Value : " << value << endl;
622 
623  Int_t npar = func->GetNpar();
624 
625  Double_t covariance[npar][npar];
626  gMinuit->mnemat(&covariance[0][0], func->GetNpar());
627 
628  Double_t partialDiff[npar];
629 
630  for(Int_t i = 0; i < npar; i++) {
631 
632  Double_t par_err = TMath::Sqrt(covariance[i][i]);
633  cout << "par_err : " << par_err << endl;
634  partialDiff[i] = FunctionDerivative(func, x, i, 0.1*par_err);
635  }
636 
637  Double_t sigma2 = 0;
638 
639  for(Int_t i = 0; i < npar; i++) {
640  for(Int_t j = 0; j < npar; j++) {
641 
642  cout << "sigma2(" << i << ", " << j << ") = "
643  << partialDiff[i]*partialDiff[j]*covariance[i][j]
644  << endl;
645 
646  sigma2+=partialDiff[i]*partialDiff[j]*covariance[i][j];
647  }
648  }
649 
650  error = TMath::Sqrt(sigma2);
651 }
652 #endif
653 
654 //______________________________________________________________________
655 Double_t IntegralDerivative(TF1* func, Double_t xlow, Double_t xhigh,
656  Double_t stepsize, Int_t parameter,
657  Double_t epsilon)
658 {
659  //*Return derivative of integral of function with respect to parameter*
660  //
661  // The derivative is computed by computing the value of the integral
662  // using all the parameters as they are except parameter where one
663  // uses (parameter + epsilon) for Int1 and (parameter - epsiolon) for
664  // Int2 and then the derivative is
665  // d_Int/d_parameter = (Int1-Int2)/(2*epsilon)
666  //
667  // Should perhaps check if parameter gets out of range....
668 
669  Double_t parameterSave = func->GetParameter(parameter);
670 
671  func->SetParameter(parameter, parameterSave+epsilon);
672  const Double_t* dummy = 0;
673 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,0,0)
674  Double_t int1 = func->Integral(xlow, xhigh, stepsize);
675 #else
676  Double_t int1 = func->Integral(xlow, xhigh, dummy, stepsize);
677 #endif
678 
679  func->SetParameter(parameter, parameterSave-epsilon);
680 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,0,0)
681  Double_t int2 = func->Integral(xlow, xhigh, stepsize);
682 #else
683  Double_t int2 = func->Integral(xlow, xhigh, dummy, stepsize);
684 #endif
685 
686  func->SetParameter(parameter, parameterSave);
687 
688  return (int1-int2)/(2*epsilon);
689 }
690 
691 //______________________________________________________________________
692 Double_t FunctionDerivative(TF1* func, Double_t x,
693  Int_t parameter, Double_t epsilon)
694 {
695  //*Return derivative of function with respect to parameter*
696  //
697  // The derivative is computed by computing the value
698  // using all the parameters as they are except parameter where one
699  // uses (parameter + epsilon) for func1 and (parameter - epsiolon) for
700  // func2 and then the derivative is
701  // d_func/d_parameter = (func1-func2)/(2*epsilon)
702  //
703  // Should perhaps check if parameter gets out of range....
704 
705  Double_t parameterSave = func->GetParameter(parameter);
706 
707  cout << "parameterSaved : " << parameterSave << endl;
708 
709  func->SetParameter(parameter, parameterSave+epsilon);
710  Double_t func1 = func->Eval(x);
711 
712  func->SetParameter(parameter, parameterSave-epsilon);
713  Double_t func2 = func->Eval(x);
714 
715  func->SetParameter(parameter, parameterSave);
716  func->Eval(x); // has to do this to reset the
717  // parameters for x1fit and x2fit
718 
719  return (func1-func2)/(2*epsilon);
720 }
721 
722 //_______________________________________________________________
723 TH1F* GetHistogramFromGraph(TGraphErrors* graph, const char* histName)
724 {
725  Int_t n = graph->GetN();
726 
727  R__ASSERT(n>1);
728 
729  Float_t xmin = graph->GetX()[0];
730  Float_t xmax = graph->GetX()[n-1];
731  Float_t dx = (xmax-xmin)/Float_t(n-1);
732 
733  R__ASSERT(dx>0);
734 
735  TH1F* hist = new TH1F(histName, 0, n, xmin, xmax);
736 
737  for(Int_t i = 0; i < n; i++) {
738 
739  Float_t xl = xmin+i*dx-0.5*dx;
740  Float_t xu = xmin+i*dx+0.5*dx;
741 
742  Float_t x = graph->GetX()[i];
743  Float_t y = graph->GetY()[i];
744  Float_t ey = graph->GetEY()[i];
745 
746  if(x<xl || x>xu)
747  cout << "ARGGHHHHHH!!!!: " << xl << ", " << x << ", " << xu << endl;
748 
749  hist->SetBinContent(i+1, y);
750  hist->SetBinError(i+1, ey);
751  }
752 
753  cout << xmin << ", " << xmax << endl;
754 
755  return hist;
756 }
757 
758 //_______________________________________________________________
759 TGraph* GetGraphFromHistogram(TH1* hist)
760 {
761  const Int_t n = hist->GetNbinsX();;
762 
763  R__ASSERT(n>0);
764 
765  TGraph* graph = new TGraph(n);
766 
767  for(Int_t i = 0; i < n; i++) {
768 
769  graph->SetPoint(i, hist->GetXaxis()->GetBinCenter(i+1),
770  hist->GetBinContent(i));
771  }
772 
773  return graph;
774 }
775 
776 
777 
778 
779 
780 //_______________________________________________________________
781 TClonesArray* CanvasPartition(TVirtualPad *C,const Int_t Nx,const Int_t Ny,
782  Float_t lMargin, Float_t rMargin,
783  Float_t bMargin, Float_t tMargin)
784 {
785  if (!C) return 0;
786 
787  TClonesArray* arr = new TClonesArray("TPad",Nx*Ny);
788 
789  // Setup Pad layout:
790  Float_t vSpacing = 0.0;
791  Float_t vStep = (1.- bMargin - tMargin - (Ny-1) * vSpacing) / Ny;
792 
793  Float_t hSpacing = 0.0;
794  Float_t hStep = (1.- lMargin - rMargin - (Nx-1) * hSpacing) / Nx;
795 
796  //Float_t vposd,vposu,vmard,vmaru,vfactor;
797  //Float_t hposl,hposr,hmarl,hmarr,hfactor;
798  for (Int_t i=0;i<Nx;i++) {
799 
800  Float_t hposl,hposr,hmarl,hmarr,hfactor;
801  if (i==0)
802  hposl = 0.0;
803  else
804  hposl = lMargin + i*(hStep + hSpacing);
805 
806  hposr = hposl + hStep + hSpacing;
807  if(i == 0)
808  hposr += lMargin;
809  if(i == Nx - 1)
810  hposr += rMargin;
811 
812  hfactor = hposr-hposl;
813  hmarl = 0.0;
814  hmarr = 0.0;
815  if (i==0)
816  hmarl = lMargin / hfactor;
817  if (i == Nx-1)
818  hmarr = rMargin / hfactor;
819 
820  for (Int_t j=0;j<Ny;j++) {
821 
822  Float_t vposd,vposu,vmard,vmaru,vfactor;
823  if (j==0)
824  vposd = 0.0;
825  else
826  vposd = bMargin + j*(vStep + vSpacing);
827 
828  vposu = vposd + vStep + vSpacing;
829  if(j == 0)
830  vposu += bMargin;
831  if(j == Ny - 1)
832  vposu += tMargin;
833 
834  vfactor = vposu-vposd;
835  vmard = 0.0;
836  vmaru = 0.0;
837  if (j==0)
838  vmard = bMargin / vfactor;
839  if (j == Ny-1)
840  vmaru = tMargin / vfactor;
841 
842 
843  C->cd(0);
844 
845  char padname[16];
846  sprintf(padname,"%s_%i_%i",C->GetName(),i,j);
847  TPad *pad = (TPad*) gROOT->FindObject(padname);
848  if (pad) delete pad;
849  std::cout << "NAME = " << padname << " " << hposl << " " << vposd << " " << hposr << " " << vposu << std::endl;
850  // pad = new TPad(name,"",hposl,vposd,hposr,vposu);
851  pad = new((*arr)[i + Nx*j]) TPad(padname,"",hposl,vposd,hposr,vposu);
852  pad->SetFillStyle(4000);
853  pad->SetLeftMargin(hmarl);
854  pad->SetRightMargin(hmarr);
855  pad->SetBottomMargin(vmard);
856  pad->SetTopMargin(vmaru);
857 
858  pad->SetFrameBorderMode(0);
859  pad->SetBorderMode(0);
860  pad->SetBorderSize(0);
861 
862  pad->Draw();
863  }
864  }
865 
866 
867 
868 
869 
870 
871 
872  return arr;
873 }
874 
875 
876 
877 
878 //_______________________________________________________________
879 void ForceScale(TH1* h, Double_t scale)
880 {
881  Bool_t is2dHist = h->InheritsFrom("TH2");
882  Int_t nbinsX = h->GetXaxis()->GetNbins();
883  Int_t nbinsY = h->GetYaxis()->GetNbins();
884  for(Int_t i = 1; i<=nbinsX; i++){
885  if(is2dHist){
886  for(Int_t j = 1; j<=nbinsY; j++){
887  h->SetBinContent(i,j,h->GetBinContent(i,j)*scale);
888  h->SetBinError(i,j,h->GetBinError(i,j)*scale);
889  }
890  }else{
891  Double_t val = h->GetBinContent(i);
892  Double_t err = h->GetBinError(i);
893  //cout << i << ": " << val << " +- " << err << " * " << scale << " = " << val*scale << " +- " << err*scale << endl;
894  h->SetBinContent(i,val*scale);
895  h->SetBinError(i,err*scale);
896  }
897  }
898 }
899 
900 
901 
902 
903 //_______________________________________________________________
904 TH1F* ForceRatio(TH1F* h1, TH1F* h2)
905 {
906 
907  Int_t nbins1 = h1->GetXaxis()->GetNbins();
908  Int_t nbins2 = h2->GetXaxis()->GetNbins();
909  if(nbins1 != nbins2) return h1;
910  // Double_t xmin = h1->GetXaxis()->GetXmin();
911  // Double_t xmax = h1->GetXaxis()->GetXmax();
912 
913  // TH1F* hout = new TH1F(Form("%s_over_%s",h1->GetName(),h2->GetName()),h1->GetTitle(),nbins1,xmin,xmax);
914  TH1F* hout = (TH1F*)h1->Clone(Form("%s_over_%s",h1->GetName(),h2->GetName()));
915  hout->Reset();
916 
917  for(Int_t i = 1; i<=nbins1; i++){
918  Double_t val1 = h1->GetBinContent(i);
919  Double_t err1 = h1->GetBinError(i);
920  Double_t val2 = h2->GetBinContent(i);
921  Double_t err2 = h2->GetBinError(i);
922  if(val2==0) continue;
923  if(val1==0) continue;
924  // cout << i << ": " << val1 << " / " << val2 << " " << val1/val2 << endl;
925  hout->SetBinContent(i,val1/val2);
926  hout->SetBinError(i,val1/val2*(err1/val1 + err2/val2));
927  // hout->SetBinError(i,TMath::Sqrt(err1*err1 + err2*val1/val2*err2*val1/val2)/val2);
928  }
929 
930  return hout;
931 }
932 
933 
934 //_______________________________________________________________
935 void NormaliseBins(TH1* h)
936 {
937  Bool_t is2dHist = h->InheritsFrom("TH2");
938  Int_t nbinsX = h->GetXaxis()->GetNbins();
939  Int_t nbinsY = h->GetYaxis()->GetNbins();
940  for(Int_t i = 1; i<=nbinsX; i++){
941  Double_t widthX = h->GetXaxis()->GetBinWidth(i);
942  if(is2dHist){
943  for(Int_t j = 1; j<=nbinsY; j++){
944  Double_t widthY = h->GetYaxis()->GetBinWidth(i);
945  h->SetBinContent(i,j,h->GetBinContent(i,j)/widthX/widthY);
946  h->SetBinError(i,j,h->GetBinError(i,j)/widthX/widthY);
947  }
948  }else{
949  Double_t val = h->GetBinContent(i);
950  Double_t err = h->GetBinError(i);
951  //cout << i << ": " << val << " +- " << err << " / " << widthX << " = " << val/widthX << " +- " << err/widthX << endl;
952  h->SetBinContent(i,val/widthX);
953  h->SetBinError(i,err/widthX);
954  }
955  }
956 }
957 
958 
959 //_______________________________________________________________
960 TGraphErrors* DivideGraph(TGraph* g1, TGraph* g2)
961 {
962  TGraphErrors* gout = new TGraphErrors();
963  TGraphErrors* g2ex = ExtrapolateGraph(g2,g1);
964  if(!g2ex) return 0;
965  Int_t n1 = g1->GetN();
966  Double_t *x1 = g1->GetX();
967  Double_t *y1 = g1->GetY();
968  Double_t *ex1 = g1->GetEX();
969  Double_t *ey1 = g1->GetEY();
970  Int_t n2 = g2ex->GetN();
971  Double_t *x2 = g2ex->GetX();
972  Double_t *y2 = g2ex->GetY();
973  //Double_t *ex2 = g2ex->GetEX();
974  Double_t *ey2 = g2ex->GetEY();
975  Int_t i = 0;
976  for(Int_t i1 = 0; i1<n1; i1++){
977  Int_t i2 = 0;
978  while(i2<n2 && x1[i1]!=x2[i2]) i2++;
979  if(y2[i2] == 0) continue;
980  gout->SetPoint(i,x1[i],y1[i1]/y2[i2]);
981  if(g1->InheritsFrom("TGraphErrors") && ey1[i1] != 0 && ey2[i2] != 0)
982  gout->SetPointError(i,ex1[i1],y1[i1]/y2[i2]*(y1[i1]/ey1[i1]+y2[i2]/ey2[i2]));
983  i++;
984  }
985 
986  return gout;
987 
988 }
989 
990 //_______________________________________________________________
991 TGraphErrors* ScaleGraph(TGraph* g1, Double_t scale)
992 {
993  TGraphErrors* gout = new TGraphErrors();
994  if(!g1) return 0;
995  Int_t n1 = g1->GetN();
996  Double_t *x1 = g1->GetX();
997  Double_t *y1 = g1->GetY();
998  Double_t *ex1 = g1->GetEX();
999  Double_t *ey1 = g1->GetEY();
1000  for(Int_t i1 = 0; i1<n1; i1++){
1001  gout->SetPoint(i1,x1[i1],y1[i1]*scale);
1002  if(g1->InheritsFrom("TGraphErrors"))
1003  gout->SetPointError(i1,ex1[i1],ey1[i1]*scale);
1004  }
1005 
1006  return gout;
1007 
1008 }
1009 
1010 //_______________________________________________________________
1011 TGraphErrors* AddGraph(TGraph* g1, TGraph* g2, Double_t val)
1012 {
1013  TGraphErrors* gout = new TGraphErrors();
1014  TGraphErrors* g2ex = ExtrapolateGraph(g2,g1);
1015  if(!g2ex) return 0;
1016  Int_t n1 = g1->GetN();
1017  Double_t *x1 = g1->GetX();
1018  Double_t *y1 = g1->GetY();
1019  Double_t *ex1 = g1->GetEX();
1020  Double_t *ey1 = g1->GetEY();
1021  Int_t n2 = g2ex->GetN();
1022  Double_t *x2 = g2ex->GetX();
1023  Double_t *y2 = g2ex->GetY();
1024  // Double_t *ex2 = g2ex->GetEX();
1025  Double_t *ey2 = g2ex->GetEY();
1026  Int_t i = 0;
1027  for(Int_t i1 = 0; i1<n1; i1++){
1028  Int_t i2 = 0;
1029  while(i2<n2 && x1[i1]!=x2[i2]) i2++;
1030  if(i2 == n2) continue;
1031  gout->SetPoint(i,x1[i1],y1[i1]+val*y2[i2]);
1032  if(g1->InheritsFrom("TGraphErrors"))
1033  gout->SetPointError(i,ex1[i1],ey1[i1]+TMath::Abs(val)*ey2[i2]);
1034  i++;
1035  }
1036 
1037  return gout;
1038 
1039 }
1040 
1041 //_______________________________________________________________
1042 TGraphErrors* DistGraph(TGraph* g1, TGraph* g2)
1043 {
1044  TGraphErrors* gout = new TGraphErrors();
1045  // TGraphErrors* g2ex = ExtrapolateGraph(g2,g1);
1046  Int_t n1 = g1->GetN();
1047  Double_t *x1 = g1->GetX();
1048  Double_t *y1 = g1->GetY();
1049  //Double_t *ex1 = g1->GetEX();
1050  //Double_t *ey1 = g1->GetEY();
1051  Int_t n2 = g2->GetN();
1052  Double_t *x2 = g2->GetX();
1053  Double_t *y2 = g2->GetY();
1054  //Double_t *ex2 = g2->GetEX();
1055  //Double_t *ey2 = g2->GetEY();
1056  Int_t i = 0;
1057  for(Int_t i1 = 0; i1<n1; i1++){
1058  Double_t dist = 10000;
1059  for(Int_t i2 = 0; i2<n2; i2++){
1060  Double_t d = (x1[i1]-x2[i2])*(x1[i1]-x2[i2])+(y1[i1]-y2[i2])*(y1[i1]-y2[i2]);
1061  if(d<dist) dist = d;
1062  }
1063  gout->SetPoint(i,x1[i1],dist);
1064  i++;
1065  }
1066 
1067  return gout;
1068 
1069 }
1070 
1071 
1072 
1073 
1074 //_______________________________________________________________
1075 TGraphErrors* ExtrapolateGraph(TGraph* g1, TGraph* g2)
1076 {
1077 
1078 
1079  Int_t n1 = g1->GetN();
1080  Double_t *x1 = g1->GetX();
1081  Double_t *y1 = g1->GetY();
1082  Double_t *ex1 = g1->GetEX();
1083  Double_t *ey1 = g1->GetEY();
1084  Int_t n2 = g2->GetN();
1085  Double_t *x2 = g2->GetX();
1086 
1087  if(n1 <= 0) return 0;
1088  if(n2 <= 0) return 0;
1089  TGraphErrors* gout = new TGraphErrors();
1090  Int_t N = 0;
1091 
1092  Double_t xi1 = x1[0];
1093  Double_t xf1 = x1[n1-1];
1094  for(Int_t i2 = 0; i2<n2; i2++){
1095  Double_t X = x2[i2];
1096  if(X<xi1) continue;
1097  if(X>xf1) continue;
1098  // Double_t xm = xi1, xp = xf1, ym, yp, exm, exp, eym, exp;
1099  Int_t ip = -1, im = -1;
1100  for(Int_t i1 = 0; i1<n1; i1++){
1101  if(x1[i1]<=X) im = i1;
1102  if(x1[i1]>=X) ip = i1;
1103  }
1104  if(im < 0) continue;
1105  if(ip < 0) continue;
1106  Double_t Y = 0, EX = 0, EY = 0;
1107  if(x1[im] == x1[ip]){
1108  Y = y1[im];
1109  if(g1->InheritsFrom("TGraphErrors")){
1110  EX = ex1[im];
1111  EY = ey1[im];
1112  }
1113  }else{
1114  Double_t fac = (X-x1[im])/(x1[ip]-x1[im]);
1115  Y = y1[im] + fac*(y1[ip]-y1[im]);
1116  if(g1->InheritsFrom("TGraphErrors")){
1117  EX = ex1[im] + fac*(ex1[ip]-ex1[im]);
1118  EY = ey1[im] + fac*(ey1[ip]-ey1[im]);
1119  }else{
1120  EX = 0;
1121  EY = 0;
1122  }
1123  }
1124  gout->SetPoint(N,X,Y);
1125  gout->SetPointError(N,EX,EY);
1126  N++;
1127  }
1128 
1129 
1130 
1131 
1132  return gout;
1133 
1134 }
TH1F * MakeNiceGraph(TGraph *g, TVirtualPad *c1, const char *opt)
TH2F * HistLogLog(const char *name, const char *title, Int_t nBinsX, Double_t xMinX, Double_t xMaxX, Int_t nBinsY, Double_t xMinY, Double_t xMaxY)
TH2F * HistLogLin(const char *name, const char *title, Int_t nBinsX, Double_t xMinX, Double_t xMaxX, Int_t nBinsY, Double_t xMinY, Double_t xMaxY)
void MakeNice1dHistoPad(TH1 *hist, TVirtualPad *c1, Double_t scaleX, Double_t scaleY, const char *opt, Bool_t copy)
TGraphErrors * DivideGraph(TGraph *g1, TGraph *g2)
TCanvas * FindCanvas(const char *canvasName, Int_t xwidth, Int_t ywidth)
void MakeNiceHistoPad(TH1 *hist, TVirtualPad *c1, Double_t scaleX, Double_t scaleY, const char *opt, Bool_t copy)
TGraphErrors * ScaleGraph(TGraph *g1, Double_t scale)
TH1F * HistLog(const char *name, const char *title, Int_t nBins, Double_t xMin, Double_t xMax)
TClonesArray * CanvasPartition(TVirtualPad *C, const Int_t Nx, const Int_t Ny, Float_t lMargin, Float_t rMargin, Float_t bMargin, Float_t tMargin)
void SetTitles(TH1 *hist, const char *title, const char *xtitle, const char *ytitle)
void SetHistMaximum(TH1 *hist, Float_t factor)
TH1F * ForceRatio(TH1F *h1, TH1F *h2)
void MakeNice1dHisto(TH1 *hist, TVirtualPad *c1, const char *opt, Bool_t copy)
void NormaliseBins(TH1 *h)
const Double_t * MakeLogBinning(Int_t nbins, Double_t xmin, Double_t xmax)
Double_t IntegralDerivative(TF1 *func, Double_t xlow, Double_t xhigh, Double_t stepsize, Int_t parameter, Double_t epsilon)
TGraphErrors * ExtrapolateGraph(TGraph *g1, TGraph *g2)
TVirtualPad * MakeNiceHisto(TH1 *hist, TVirtualPad *c1, const char *opt, Bool_t copy)
void ForceScale(TH1 *h, Double_t scale)
void MakeNice2dHistoPad(TH2 *hist, TVirtualPad *c1, Double_t scaleX, Double_t scaleY, const char *option, Bool_t copy)
void SetHistStyle(TH1 *hist)
TH1F * HistInvert(TH1 *hist)
void MakeNice2dHisto(TH2 *hist, TVirtualPad *c1, const char *option, Bool_t copy)
Double_t FunctionDerivative(TF1 *func, Double_t x, Int_t parameter, Double_t epsilon)
TGraph * GetGraphFromHistogram(TH1 *hist)
TH1F * GetHistogramFromGraph(TGraphErrors *graph, const char *histName)
TH2F * HistLinLog(const char *name, const char *title, Int_t nBinsX, Double_t xMinX, Double_t xMaxX, Int_t nBinsY, Double_t xMinY, Double_t xMaxY)
TGraphErrors * DistGraph(TGraph *g1, TGraph *g2)
TGraphErrors * AddGraph(TGraph *g1, TGraph *g2, Double_t val)