StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
histutil.h
1 //#ifndef HISTUTIL_H
2 //#define HISTUTIL_H
3 
4 #include <cstring>
5 
6 // 3d to 2d
7 TH2D* HistSlice(TH3* h3,const char* basename, const char* basetitle,
8  const char* slicename, Axis_t min, Axis_t max,
9  const char* slicetype,
10  Option_t* opt="");
11 
12 
13 // 3d to 2d but keep varying width histos
14 TH2D* HistSlice(TH3* h3,const char* basename, const char* basetitle,
15  const char* slicename, const char* slicetype);
16 
17 // 2d to 1d
18 TH1D* HistSlice(TH2* h2,const char* basename, const char* basetitle,
19  const char* slicename,Axis_t min, Axis_t max,
20  const char* projaxis,
21  Option_t* opt="");
22 
23 // 3d to 1d
24 // if project to x then slicename1="y",slicename2="z"
25 // if project to y then slicename2="z",slicename2="x"
26 // if project to z then slicename1="x",slicename2="y"
27 
28 TH1D* HistSlice(TH3* h3,const char* basename, const char* basetitle,
29  const char* slicename1,Axis_t min1, Axis_t max1,
30  const char* slicename2,Axis_t min2, Axis_t max2,
31  const char* projecttype,
32  Option_t* opt="");
33 
34 TProfile* Profile(TH3* h3,const char* basename, const char* basetitle,
35  const char* slicename, Axis_t min, Axis_t max,
36  const char* slicetype,
37  const char* projecttype,
38  Option_t* opt="");
39 
40 TProfile* Profile(TH2* h2,const char* basename,const char* basetitle,
41  const char* slicename, Axis_t min, Axis_t max,
42  const char* profaxis,
43  Option_t* opt="");
44 
45 TH1D* Rms(TH3* h3,const char* basename, const char* basetitle,
46  const char* slicename, Axis_t min, Axis_t max,
47  const char* slicetype,
48  const char* projecttype);
49 TH1D* Rms(TH2* h2,const char* basename,
50  const char* projecttype);
51 // rebin an axis of a 3D histogram (with varying bins);
52 
53 void Rebin(TH3* h3,int axisType,int nBin,Axis_t* ary);
54 
55 Int_t FindMinBin(const TAxis* axis,Stat_t val);
56 Int_t FindMaxBin(const TAxis* axis,Stat_t val);
57 void SetRange(const TAxis* axis,Stat_t min,Stat_t max);
58 Stat_t FindMax(TH1* ha,TH1* hb);
59 Stat_t FindMin(TH1* ha,TH1* hb);
60 Stat_t FindMax(TH1* ha,TH1* hb,TH1* hc);
61 void Divide(TCanvas* c1,int nx,int ny,const char* title=0,const char* file=0);
62 void Print(TCanvas* c1, const char* outDir, TString basename);
63 
64 void SetMinMax(TH1* h);
65 void SetMinMax(TH1* h,float min, float max);
66 void GetMeanRms(TH1* ha,float& cmean,float& crms);
67 void GetMeanRmsYield(TH1* ha,float& cmean,float& crms);
68 void PrintMeanRms(TH1* ha,float x, float y,float size=0.04);
69 void PrintMeanRmsYield(TH1* ha,float x, float y,float size=0.04);
70 void ReplaceTitle(TH1* h,char* current, char* replace);
71 void DrawLine(TH1* h,float y=0);
72 
73 void showHistogramValues(TH1D* h1, char* more=0);
74 void showTGraphValues(TGraphAsymmErrors* graph,char* more=0);
75 void dump(TH1* h1);
76 void dump(TGraph* g);
77 void dump(TGraphAsymmErrors* g);
78 void scale2D(TH2*);
79 //-----------------------------------------------
80 
81 void DrawLine(TH1* h,float y)
82 {
83  TLine* line=new TLine; TAxis* axis=h->GetXaxis();
84  line->DrawLine(axis->GetBinLowEdge(axis->GetFirst()),y,
85  axis->GetBinUpEdge(axis->GetLast()),y);
86 
87 }
88 
89 void ReplaceTitle(TH1* h,char* current, char* replace)
90 {
91  TString s=h->GetTitle();
92  s.ReplaceAll(current,replace);
93  h->SetTitle(s.Data());
94 }
95 
96 void SetMinMax(TH1* h,float min, float max)
97 {
98  h->SetMinimum(min);
99  h->SetMaximum(max);
100 }
101 
102 void SetMinMax(TH1* h)
103 {
104  h->SetMinimum(h->GetMinimum()*0.8);
105  h->SetMaximum(h->GetMaximum()*1.2);
106 }
107 
108 void GetMeanRms(TH1* ha,float& mean,float& rms)
109 {
110  TAxis* axis=ha->GetXaxis();
111  mean=0;rms=0;
112  int n=0;
113 
114  for(int i=axis->GetFirst(); i<=axis->GetLast(); i++){
115  if(ha->GetBinContent(i)>0){
116  float val = (ha->GetBinCenter(i)*ha->GetBinContent(i));
117  mean += val;
118  rms += val*val;
119 
120  }
121  }
122  n=ha->Integral(axis->GetFirst(),axis->GetLast());
123  //n=ha->Integral(1,axis->GetNbins());
124  if(n>0){
125  mean /= n; rms /= n; rms=sqrt(rms - mean*mean);
126  }
127 }
128 
129 void GetMeanRmsYield(TH1* ha,float& mean,float& rms)
130 {
131  TAxis* axis=ha->GetXaxis();
132  mean=0;rms=0;
133  int n=0;
134 
135  // cout << ha->GetName() << endl;
136 
137  for(int i=axis->GetFirst(); i<=axis->GetLast(); i++){
138  if(ha->GetBinContent(i)>0){
139  float val = ha->GetBinContent(i);
140  mean += val;
141  rms += val*val;
142  // cout << "bin=" << i << ",val=" << val << endl;
143 
144  }
145  }
146  n=axis->GetLast()-axis->GetFirst()+1;
147  // cout << "n=" << n << endl;
148  //n=ha->Integral(1,axis->GetNbins());
149  if(n>0){
150  mean /= n; rms /= n; rms=sqrt(rms - mean*mean);
151  }
152  // cout << mean << ", " << rms << endl;
153 }
154 
155 void PrintMeanRms(TH1* ha,float x,float y,float size)
156 {
157  TText* text=new TText; float mean(0),rms(0);
158  text->SetTextSize(size);
159  char buf[100];
160  mean = ha->GetMean(); rms = ha->GetRMS();
161  // if(compute)GetMeanRms(ha,mean,rms);
162  sprintf(buf,"mean: %.3f",mean);
163  text->DrawTextNDC(x,y,buf);
164  sprintf(buf,"rms: %.3f",rms);
165  text->DrawTextNDC(x,y-0.05,buf);
166 }
167 
168 void PrintMeanRmsYield(TH1* ha,float x,float y,float size)
169 {
170  TText* text=new TText; float mean(0),rms(0);
171  text->SetTextSize(size);
172  char buf[100];
173  //mean = ha->GetMean(); rms = ha->GetRMS();
174  GetMeanRmsYield(ha,mean,rms);
175  sprintf(buf,"mean: %.3f",mean);
176  text->DrawTextNDC(x,y,buf);
177  sprintf(buf,"rms: %.3f",rms);
178  text->DrawTextNDC(x,y-0.05,buf);
179 }
180 
181 // divide a canvas and print a title
182 
183 void Divide(TCanvas *c1,int nx,int ny,const char* title,const char* file)
184 {
185  c1->Clear();
186  TString cTitle = (title) ? title : c1->GetTitle();
187  if(title) c1->SetTitle(title);
188 
189  TPaveLabel *mP = new TPaveLabel(.02,.95,.5,.98,title);
190  mP->SetFillColor(0);
191  mP->SetBorderSize(0);
192  mP->Draw();
193 
194  if(file){
195  TText tFile;
196  tFile.SetTextSize(0.02);
197  tFile.DrawTextNDC(0.1,0.0,file);
198  }
199 
200  char buf[10];
201 
202  Double_t ymax = .94;
203  Double_t dy = ymax/ny;
204  Double_t dx = 1./nx;
205  Double_t ymargin = 0.02;
206  double xmargin = 0.02;
207  int n=0;
208  for (Int_t iy=0;iy<ny;iy++) {
209  double yhigh = ymax - iy*dy - ymargin;
210  double ylow = yhigh - dy + 3*ymargin;
211  if (ylow < 0) ylow = 0;
212  if (ylow > yhigh) continue;
213  for (Int_t ix=0;ix<nx;ix++) {
214  n++;
215  double xlow = ix*dx + xmargin;
216  double xhigh = xlow +dx -2*xmargin;
217  if (xlow > xhigh) continue;
218 
219  TPad* p = new TPad(buf,buf,xlow,ylow,xhigh,yhigh);
220  p->SetNumber(n);
221  p->Draw();
222  }
223  }
224 
225 }
226 
227 // given a maximum value, find it's bin
228 Int_t
229 FindMinBin(const TAxis* axis,Stat_t val)
230 {
231  Int_t bin = axis->FindBin(val);
232  if(bin<0) return 1;
233  return (axis->GetBinCenter(bin)<val) ? ++bin : bin;
234 
235 }
236 
237 // given a minimum value, find it's bin
238 Int_t FindMaxBin(const TAxis* axis,Stat_t val)
239 {
240  Int_t bin = axis->FindBin(val);
241  if(bin>axis->GetNbins()) return axis->GetNbins();
242  return (axis->GetBinCenter(bin)>val) ? --bin : bin;
243 }
244 
245 void SetRange(const TAxis* axis, Stat_t min, Stat_t max)
246 {
247  axis->SetRange(FindMinBin(axis,min),FindMaxBin(axis,max));
248 }
249 Stat_t
250 FindMax(TH1* ha, TH1* hb)
251 {
252  Stat_t maxa = ha->GetMaximum();
253  Stat_t maxb = hb->GetMaximum();
254  return (maxa > maxb) ? maxa : maxb;
255 
256 }
257 
258 Stat_t
259 FindMin(TH1* ha, TH1* hb)
260 {
261  Stat_t mina = ha->GetMinimum();
262  Stat_t minb = hb->GetMinimum();
263  return (mina > minb) ? mina : minb;
264 
265 }
266 
267 Stat_t
268 FindMax(TH1* ha,TH1* hb,TH1* hc)
269 {
270  Stat_t max1 = FindMax(ha,hb);
271  Stat_t max2 = FindMax(hb,hc);
272  return (max1>max2) ? max1 : max2;
273 }
274 
275 
276 // 3d to 2d
277 
278 TH2D* HistSlice(TH3* h3,
279  const char* basename,const char* basetitle,
280  const char* slicename, Axis_t min, Axis_t max,
281  const char* slicetype,
282  Option_t* opt)
283 {
284  char name[200], title[200],option[20],buf[100];
285  Int_t minBin, maxBin;
286  Axis_t mmin,mmax;
287  TH2D* h2;
288 
289  if(!h3){
290  cout << "Null pointer for h3 : " << basename << endl;
291  return 0;
292  }
293 
294  // set the range
295  // x slice
296 
297  TAxis* axis = 0; TAxis* yAxis=0,*xAxis=0;
298  if(strcmp(slicetype,"yz")==0 || strcmp(slicetype,"zy")==0 ){
299  axis = h3->GetXaxis();
300  if(strcmp(slicetype,"yz")==0){ yAxis=h3->GetYaxis(); xAxis=h3->GetZaxis();}
301  else { yAxis=h3->GetZaxis(); xAxis=h3->GetYaxis(); }
302  }
303  // y slice
304  else if (strcmp(slicetype,"xz")==0 || strcmp(slicetype,"zx")==0){
305  axis = h3->GetYaxis();
306  if(strcmp(slicetype,"xz")==0){ yAxis=h3->GetXaxis(); xAxis=h3->GetZaxis();}
307  else { yAxis=h3->GetZaxis(); xAxis=h3->GetXaxis(); }
308  }
309  // z slice
310  else if (strcmp(slicetype,"xy")==0 || strcmp(slicetype,"yx")==0){
311  axis = h3->GetZaxis();
312  if(strcmp(slicetype,"xy")==0){ yAxis=h3->GetXaxis(); xAxis=h3->GetYaxis();}
313  else { yAxis=h3->GetYaxis(); xAxis=h3->GetXaxis(); }
314  }
315  else {
316  cerr << "Wrong value for slice axis for " << basename << endl;
317  return 0;
318  }
319 
320  if(min>=max){
321  minBin = 1; maxBin = axis->GetNbins();
322  }
323  else{
324  minBin = FindMinBin(axis,min);
325  maxBin = FindMaxBin(axis,max);
326  }
327  sprintf(option,"%s%s",slicetype,opt);
328 
329  if(!slicename){
330  strcpy(buf,axis->GetTitle()); slicename=buf;
331  }
332 
333 
334  axis->SetRange(minBin,maxBin);
335 
336  h2 = (TH2D*) h3->Project3D(option);
337  h2->SetYTitle(yAxis->GetTitle());
338  h2->SetXTitle(xAxis->GetTitle());
339 
340 
341  mmin = axis->GetBinLowEdge(minBin);
342  mmax = axis->GetBinUpEdge(maxBin);
343 
344  // set the name and title
345  sprintf(name,"%s(%.2f<%s<%.2f)",basename,mmin,slicename,mmax);
346  sprintf(title,"%s (%.2f<%s<%.2f)",basetitle,mmin,slicename,mmax);
347 
348  h2->SetName(name); h2->SetTitle(title);
349 
350  axis->SetRange(0,9999999);
351 
352  return h2;
353 }
354 
355 // 3d to 2d var bins
356 
357 TH2D* HistSlice(TH3* h3,const char* basename, const char* basetitle,
358  const char* slicename, const char* slicetype)
359 {
360  TString opt = slicetype;
361 
362  int projX=-1,projY=-1,projZ=-1;
363 
364  if(opt.Contains("xy")){ projX=1; projY=0; projZ=2; }
365  else if(opt.Contains("yx")){ projX=0; projY=1; projZ=2;}
366  else if(opt.Contains("zy")){ projX=1; projY=2; projZ=0; }
367  else if(opt.Contains("yz")){ projX=2; projY=1; projZ=0; }
368  else if(opt.Contains("xz")){ projX=2; projY=0; projZ=1; }
369  else if(opt.Contains("zx")){ projX=0; projY=2; projZ=1; }
370 
371  TAxis* xAxis=0, *yAxis=0, *zAxis=0;
372  switch(projX){
373  case 0: xAxis = h3->GetXaxis(); break;
374  case 1: xAxis = h3->GetYaxis(); break;
375  case 2: xAxis = h3->GetZaxis(); break;
376  }
377  switch(projY){
378  case 0: yAxis = h3->GetXaxis(); break;
379  case 1: yAxis = h3->GetYaxis(); break;
380  case 2: yAxis = h3->GetZaxis(); break;
381  }
382  switch(projZ){
383  case 0: zAxis = h3->GetXaxis(); break;
384  case 1: zAxis = h3->GetYaxis(); break;
385  case 2: zAxis = h3->GetZaxis(); break;
386  }
387 
388  // create the 2d histogram
389  int ixmin = xAxis->GetFirst();
390  int ixmax = xAxis->GetLast();
391  int iymin = yAxis->GetFirst();
392  int iymax = yAxis->GetLast();
393  int izmin = zAxis->GetFirst();
394  int izmax = zAxis->GetLast();
395  int nx = ixmax-ixmin+1;
396  int ny = iymax-iymin+1;
397  int nz = izmax-izmin+1;
398 
399  char* name = new char[strlen(basename)+200];
400  char* title = new char[strlen(basename)+200];
401 
402  sprintf(name,"%sSbin%dto%d",basename,izmin,izmax);
403  sprintf(title,"%s(%.2f<%s<%.2f)",basetitle,
404  zAxis->GetBinLowEdge(izmin),slicename,zAxis->GetBinUpEdge(izmax));
405 
406  /*
407  cout << "xmin=" << xAxis->GetBinLowEdge(ixmin)
408  << ",xmax=" << xAxis->GetBinUpEdge(ixmax)
409  << endl
410  << "ymin=" << yAxis->GetBinLowEdge(iymin)
411  << ",ymax=" << yAxis->GetBinUpEdge(iymax)
412  << endl
413  << "zmin=" << zAxis->GetBinLowEdge(izmin)
414  << ",zmax=" << zAxis->GetBinUpEdge(izmax)
415  << endl;
416  */
417 
418  double* xAry=0; xAry=xAxis->GetXbins()->GetArray();
419  double* yAry=0; yAry=yAxis->GetXbins()->GetArray();
420 
421  // normal binning
422  TH2D* h2=0;
423  if(!xAry){
424  h2=(TH2D*) h3->Project3D(slicetype);
425  h2->SetName(name); h2->SetTitle(title);
426  return h2;
427  }
428 
429  // varying binning
430  h2=new TH2D(name,title,
431  nx,&xAry[ixmin-1],ny,&yAry[iymin-1]);
432 
433  // fill the new histogram
434 
435  int iBinXreal=0, iBinYreal=0, iBinZreal=0;
436  for(int ix=ixmin; ix<=ixmax; ix++){
437  switch(projX){
438  case 0: iBinXreal = ix; break;
439  case 1: iBinYreal = ix; break;
440  case 2: iBinZreal = ix; break;
441  }
442  for(int iy=iymin; iy<=iymax; iy++){
443  switch(projY){
444  case 0: iBinXreal = iy; break;
445  case 1: iBinYreal = iy; break;
446  case 2: iBinZreal = iy; break;
447  }
448  double content=0;
449  for(int iz=izmin; iz<=izmax; iz++){
450  switch(projZ){
451  case 0: iBinXreal = iz; break;
452  case 1: iBinYreal = iz; break;
453  case 2: iBinZreal = iz; break;
454  }
455  int bin = h3->GetBin(iBinXreal,iBinYreal,iBinZreal);
456  h2->Fill(xAxis->GetBinCenter(ix),
457  yAxis->GetBinCenter(iy),h3->GetBinContent(bin));
458  }
459  }
460  }
461 
462  return h2;
463 }
464 
465 
466 // 2d to 1d
467 
468 TH1D* HistSlice(TH2* h2,const char* basename, const char* basetitle,
469  const char* slicename,Axis_t min, Axis_t max,
470  const char* projaxis,
471  Option_t* opt)
472 {
473  char name[200], title[200];
474  Axis_t mmin(0),mmax(0);
475  Int_t minBin(0), maxBin(0);
476  TH1D* h1;
477  TAxis* axis = 0;
478  char buf[100];
479 
480  if(strcmp(projaxis,"x")==0){
481  axis = h2->GetYaxis();
482  }
483  else if(strcmp(projaxis,"y")==0){
484  axis = h2->GetXaxis();
485  }
486  else{
487  cout << "The argument must be either 'x' or 'y' " << endl;
488  return 0;
489  }
490 
491  if(min>=max){
492  minBin = 1; maxBin = axis->GetNbins();
493 
494  }
495  else{
496  minBin = FindMinBin(axis,min);
497  maxBin = FindMaxBin(axis,max);
498  }
499 
500  if(strcmp(projaxis,"x")==0){
501  h1= h2->ProjectionX("dummy",minBin,maxBin,opt);
502  if(!slicename){
503  strcpy(buf,axis->GetTitle()); slicename=buf;
504  }
505  h1->SetXTitle(h2->GetXaxis()->GetTitle());
506  }
507  else if(strcmp(projaxis,"y")==0){
508  h1= h2->ProjectionY("dummy",minBin,maxBin,opt);
509  if(!slicename){
510  strcpy(buf,axis->GetTitle()); slicename=buf;
511  }
512  h1->SetXTitle(h2->GetYaxis()->GetTitle());
513  }
514 
515  mmin = axis->GetBinLowEdge(minBin);
516  mmax = axis->GetBinUpEdge(maxBin);
517 
518 
519  //
520  // change the name and title
521  //
522  sprintf(name,"%s(%.2f<%s<%.2f)",basename,mmin,slicename,mmax);
523  sprintf(title,"%s (%.2f<%s<%.2f)",basetitle,mmin,slicename,mmax);
524  h1->SetTitle(title);
525  h1->SetName(name);
526 
527  return h1;
528 }
529 
530 TH1D* HistSlice(TH3* h3,const char* basename, const char* basetitle,
531  const char* slicename1,Axis_t min1, Axis_t max1,
532  const char* slicename2,Axis_t min2, Axis_t max2,
533  const char* projecttype,
534  Option_t* opt)
535 {
536  char slicetype[10];
537  char buf1[100],buf2[100];
538  TAxis* axis1, *axis2;
539  if(strcmp(projecttype,"x")==0){ // slice out y
540  strcpy(slicetype,"zx");
541  axis1=h3->GetYaxis(); axis2=h3->GetZaxis();
542  }
543  else if(strcmp(projecttype,"y")==0){ // slice out z
544  strcpy(slicetype,"xy");
545  axis1=h3->GetZaxis(); axis2=h3->GetXaxis();
546  }
547  else if(strcmp(projecttype,"z")==0){ // slice out x
548  strcpy(slicetype,"yz");
549  axis1=h3->GetXaxis(); axis2=h3->GetYaxis();
550  }
551  else{
552  cout << "error: wrong project type " << basename << endl;
553  return;
554  }
555  if(!slicename1){
556  strcpy(buf1,axis1->GetTitle()); slicename1=buf1;
557  }
558  if(!slicename2){
559  strcpy(buf2,axis2->GetTitle()); slicename2=buf2;
560  }
561 
562  TH2D* h2 =
563  (TH2D*) HistSlice(h3,basename,basetitle,
564  slicename1,min1,max1,
565  slicetype,opt);
566 
567  return HistSlice(h2,h2->GetName(),h2->GetTitle(),
568  slicename2,min2,max2,
569  "x",opt);
570 
571 }
572 
573 TProfile* Profile(TH3* h3,const char* basename, const char* basetitle,
574  const char* slicename, Axis_t min, Axis_t max,
575  const char* slicetype,
576  const char* projecttype,
577  Option_t* opt)
578 {
579  // first project to 2d
580  TH2D* h2 = HistSlice(h3,basename,basetitle,slicename,min,max,
581  slicetype);
582  char name[200];
583  TProfile* p;
584  // profile
585  if(strstr(projecttype,"x")){
586  sprintf(name,"%s_profx",h2->GetName());
587  p=h2->ProfileX(name,0,9999999,opt);
588  }
589  else{
590  sprintf(name,"%s_profy",h2->GetName());
591  p=h2->ProfileY(name,0,9999999,opt);
592  }
593  p->SetTitle(h2->GetTitle());
594  return p;
595 }
596 
597 TProfile* Profile(TH2* h2,const char* basename,const char* basetitle,
598  const char* slicename, Axis_t min, Axis_t max,
599  const char* profaxis,
600  Option_t* opt)
601 {
602  TProfile* p;
603  char name[200], title[200];
604  Axis_t mmin(0),mmax(0);
605  Int_t minBin(0), maxBin(0);
606  TAxis* axis = 0;
607  char buf[100];
608 
609  if(strcmp(profaxis,"x")==0){
610  axis = h2->GetYaxis();
611  }
612  else if(strcmp(profaxis,"y")==0){
613  axis = h2->GetXaxis();
614  }
615  else{
616  cout << "The argument must be either 'x' or 'y' " << endl;
617  return 0;
618  }
619 
620  if(min>=max){
621  minBin = 1; maxBin = axis->GetNbins();
622 
623  }
624  else{
625  minBin = FindMinBin(axis,min);
626  maxBin = FindMaxBin(axis,max);
627  }
628 
629  if(strcmp(profaxis,"x")==0){
630  p= h2->ProfileX("dummy",minBin,maxBin,opt);
631  if(!slicename){
632  strcpy(buf,axis->GetTitle()); slicename=buf;
633  }
634  p->SetXTitle(h2->GetXaxis()->GetTitle());
635  }
636  else if(strcmp(profaxis,"y")==0){
637  p= h2->ProfileY("dummy",minBin,maxBin,opt);
638  if(!slicename){
639  strcpy(buf,axis->GetTitle()); slicename=buf;
640  }
641  p->SetXTitle(h2->GetYaxis()->GetTitle());
642  }
643 
644  mmin = axis->GetBinLowEdge(minBin);
645  mmax = axis->GetBinUpEdge(maxBin);
646 
647 
648  //
649  // change the name and title
650  //
651  sprintf(name,"%s(%.2f<%s<%.2f)",basename,mmin,slicename,mmax);
652  sprintf(title,"%s (%.2f<%s<%.2f)",basetitle,mmin,slicename,mmax);
653  p->SetTitle(title);
654  p->SetName(name);
655 
656  return p;
657 
658 }
659 
660 TProfile* Rms(TH2* h2,const char* basename,
661  const char* projecttype)
662 {
663  char name[100], xtitle[100];
664  TProfile* p;
665  if(strstr(projecttype,"x")){
666  p=h2->ProfileX("dummy",1,h2->GetNbinsY(),"s");
667  strcpy(xtitle,h2->GetXaxis()->GetTitle());
668  }
669  else if(strstr(projecttype,"y")){
670  p=h2->ProfileY("dummy",1,h2->GetNbinsX(),"s");
671  strcpy(xtitle,h2->GetYaxis()->GetTitle());
672  }
673  else{
674  cout << "wrong type " << projecttype << endl; exit(1);
675  }
676  sprintf(name,"%s%s_%s",h2->GetName(),basename,projecttype);
677  p->SetName(name);
678  sprintf(name,"%s_%s_rms",h2->GetName(),basename,projecttype);
679  Int_t nBin= p->GetNbinsX();
680  TH1* h1 = new TH1D(name,name,nBin,
681  p->GetXaxis()->GetBinLowEdge(1),
682  p->GetXaxis()->GetBinUpEdge(nBin));
683  // now loop over the tprofile and fill the histogram
684  for(int i=1;i<=nBin;i++){
685  h1->SetBinContent(i,p->GetBinError(i));
686  // set the error as 1% of the content
687  h1->SetBinError(i,p->GetBinContent(i)*.001);
688  }
689  h1->SetXTitle(xtitle);
690  return h1;
691 }
692 
693 TH1D* Rms(TH3* h3,const char* basename, const char* basetitle,
694  const char* slicename, Axis_t min, Axis_t max,
695  const char* slicetype,
696  const char* projecttype)
697 {
698  char name[200]; sprintf(name,"%s_prms",basename);
699  TProfile* p = Profile(h3,name,basetitle,slicename,
700  min,max,slicetype,projecttype,"s");
701 
702  // make a 1 d histogram
703 
704  sprintf(name,"%s_h1",p->GetName());
705  Int_t nBin= p->GetNbinsX();
706  TH1D* h1 = new TH1D(name,p->GetTitle(),nBin,
707  p->GetXaxis()->GetBinLowEdge(1),
708  p->GetXaxis()->GetBinUpEdge(nBin));
709 
710  // now loop over the tprofile and fill the histogram
711  for(int i=1;i<=nBin;i++){
712  h1->SetBinContent(i,p->GetBinError(i));
713  // set the error as 1% of the content
714  h1->SetBinError(i,p->GetBinContent(i)*.001);
715  }
716  h1->SetXTitle(p->GetXaxis()->GetTitle());
717  return h1;
718 
719 }
720 
721 void
722 Printps(TCanvas* c1, const char* outDir, TString basename)
723 {
724 
725  if(!basename.EndsWith(".ps")) basename += ".ps";
726  TString s = outDir; s += "/"; s += basename;
727  c1->Print(s.Data());
728 
729 }
730 
731 void
732 Printgif(TCanvas* c1, const char* outDir, TString basename)
733 {
734 
735  if(!basename.EndsWith(".gif")) basename += ".gif";
736  TString s = outDir; s += "/"; s += basename;
737  c1->Print(s.Data());
738 
739 }
740 
741 // errors are not reset...
742 void
743 Rebin(TH3* h3,int axisType,int nBin, Axis_t* ary)
744 {
745  if(!h3) { cout << "null pointer?" << endl; return; }
746 
747 
748  // clone the histogram
749  TH3* hClone = (TH3*)h3->Clone();
750  TString name = h3->GetName(); name += "Clone";
751 
752  hClone->SetName(name.Data());
753 
754  TAxis *rebinAxis=0, *aAxis=0, *bAxis=0;
755  TAxis *xAxis=hClone->GetXaxis(),
756  *yAxis=hClone->GetYaxis(), *zAxis=hClone->GetZaxis();
757  TAxis *newAxis,*newAaxis,*newBaxis;
758 
759  switch(axisType){
760  case 0:
761  rebinAxis=xAxis; aAxis=yAxis; bAxis=zAxis;
762  newAxis=h3->GetXaxis(); newAaxis=h3->GetYaxis(); newBaxis=h3->GetZaxis();
763  break;
764  case 1:
765  rebinAxis=yAxis; aAxis=xAxis; bAxis=zAxis;
766  newAxis=h3->GetYaxis(); newAaxis=h3->GetXaxis(); newBaxis=h3->GetZaxis();
767  break;
768  case 2:
769  rebinAxis=zAxis; aAxis=xAxis; bAxis=yAxis;
770  newAxis=h3->GetZaxis(); newAxis=h3->GetXaxis(); newBaxis=h3->GetYaxis();
771  break;
772  default: cout << "Wrong axis type: " << axisType << endl; exit(-1);
773  }
774 
775 
776  // set new axis of the original histogram
777  cout << "old axis " << rebinAxis->GetNbins() << endl;
778  newAxis->Set(nBin,ary);
779  cout << "new axis " << newAxis->GetNbins() << endl;
780 
781  int aNbin=aAxis->GetNbins(), bNbin=bAxis->GetNbins();
782  double* aAry = new double[aAxis->GetNbins()+1];
783  double* bAry = new double[bAxis->GetNbins()+1];
784 
785  for(int i=0; i<aNbin; i++){
786  int bin=i+1;
787  aAry[i]=aAxis->GetBinLowEdge(bin);
788  }
789  aAry[aNbin] = aAxis->GetBinUpEdge(aNbin);
790  newAaxis->Set(aNbin,aAry);
791 
792  for(int i=0; i<bNbin; i++){
793  int bin=i+1;
794  bAry[i]=bAxis->GetBinLowEdge(bin);
795  }
796  bAry[bNbin] = bAxis->GetBinUpEdge(bNbin);
797  newBaxis->Set(bNbin,bAry);
798 
799  delete aAry; delete bAry;
800 
801  // now fill it
802 
803  double newval=0;
804  int iBinXreal=0, iBinYreal=0, iBinZreal=0, bin=0;
805 
806  for(int aBin=1; aBin<=aAxis->GetNbins(); aBin++){
807  for(int bBin=1; bBin<=bAxis->GetNbins(); bBin++){
808  for(int newBin=1; newBin<=nBin; newBin++){
809  newval=0;
810 
811  for(int rebinBin=1; rebinBin<=rebinAxis->GetNbins(); rebinBin++){
812  if(newAxis->FindBin(rebinAxis->GetBinCenter(rebinBin))==newBin){
813  switch(axisType){
814  case 0: iBinXreal=rebinBin; iBinYreal=aBin; iBinZreal=bBin; break;
815  case 1: iBinXreal=aBin; iBinYreal=rebinBin; iBinZreal=bBin; break;
816  case 2: iBinXreal=aBin; iBinYreal=bBin; iBinZreal=rebinBin; break;
817  }
818  bin = hClone->GetBin(iBinXreal, iBinYreal, iBinZreal);
819  newval += hClone->GetBinContent(bin);
820  }
821  }
822  // set new bin content
823  switch(axisType){
824  case 0: iBinXreal=newBin; iBinYreal=aBin; iBinZreal=bBin; break;
825  case 1: iBinXreal=aBin; iBinYreal=newBin; iBinZreal=bBin; break;
826  case 2: iBinXreal=aBin; iBinYreal=bBin; iBinZreal=newBin; break;
827  }
828  bin = h3->GetBin(iBinXreal,iBinYreal,iBinZreal);
829  h3->SetBinContent(bin,newval);
830 
831  } // newBin
832  }
833  }
834  // reality check
835  /*
836  TH1* h1Clone=hClone->Project3D("x");
837  TH1* h1New =h3->Project3D("x");
838 
839  for(int i=1; i<=h1New->GetNbinsX(); i++){
840  cout << "new bin " << i
841  << ",low = "<< h1New->GetXaxis()->GetBinLowEdge(i)
842  << ",high = " << h1New->GetXaxis()->GetBinUpEdge(i)
843  << ",val = " << h1New->GetBinContent(i) << endl;
844  double sum=0;
845  for(int j=1; j<=h1Clone->GetNbinsX(); j++){
846  if(h1New->GetXaxis()->FindBin(h1Clone->GetXaxis()->GetBinCenter(j))==i){
847  cout << "\told bin " << j
848  << ",low = " << h1Clone->GetXaxis()->GetBinLowEdge(j)
849  << ",high = " << h1Clone->GetXaxis()->GetBinUpEdge(j)
850  << ",val = " << h1Clone->GetBinContent(j) << endl;
851  sum+= h1Clone->GetBinContent(j);
852  }
853  }
854  cout << "\tsum = " << sum << endl;
855  }
856  */
857 
858  // get rid of the clone
859  delete hClone; delete aAry; delete bAry;
860 
861 }
862 //#endif
863 
864 //____________________
865 
866 void
867 dump(TH1* h1)
868 {
869  cout << h1->GetName() << endl;
870  Int_t nBin = h1->GetNbinsX();
871 
872  for(Int_t i=1; i<=nBin; i++){
873  float errorFrac = (h1->GetBinContent(i)>0)?
874  h1->GetBinError(i)/h1->GetBinContent(i) : 0;
875 
876  printf("\t%d : %.4f -- %.4f \t%.4f \t%.4f \t%.4f\n",
877  i,
878  h1->GetXaxis()->GetBinLowEdge(i),
879  h1->GetXaxis()->GetBinUpEdge(i),
880  h1->GetBinContent(i),
881  h1->GetBinError(i),
882  errorFrac);
883  }
884 }
885 
886 void
887 dump(TGraph* g)
888 {
889  cout << g->GetName() << endl;
890  double *x = g->GetX();
891  double *y = g->GetY();
892 
893  for(int i=0; i<g->GetN(); i++){
894  printf("\t%d : %.4f \t %.4f\n",i,x[i],y[i]);
895  }
896 }
897 
898 void
899 dump(TGraphAsymmErrors* graph)
900 {
901  double* xValues = graph->GetX();
902  double* yValues = graph->GetY();
903  double* errXLow = graph->GetEXlow();
904  double* errXHigh = graph->GetEXhigh();
905 
906  double* errYLow = graph->GetEYlow();
907  double* errYHigh = graph->GetEYhigh();
908 
909  for(int i=0; i<graph->GetN(); i++){
910  printf("\t%d : %.4f < %.4f <%.4f \t %.3e %.3e\n",
911  i,xValues[i]-errXLow[i],xValues[i],xValues[i]+errXHigh[i],
912  yValues[i],errYLow[i]);
913  }
914 
915 
916 }
917 
918 
919 void
920 showHistogramValues(TH1D* h1, char* more)
921 {
922  *ofVal << "#################################################" << endl;
923  *ofVal << h1->GetName();
924  if(more) *ofVal << ":::::"<<more << endl;
925  else *ofVal << endl;
926 
927  Int_t nBin = h1->GetNbinsX();
928 
929  Double_t errorFrac;
930 
931  // no under and over flows
932 
933  for(Int_t i=1; i<=nBin; i++){
934  errorFrac = (h1->GetBinContent(i)>0)?
935  h1->GetBinError(i)/h1->GetBinContent(i) : 0;
936 
937  *ofVal << "\t" << h1->GetXaxis()->GetBinLowEdge(i) << " < pt < "
938  << h1->GetXaxis()->GetBinUpEdge(i) << " "
939  << "\t" << h1->GetBinContent(i) << " \t" << h1->GetBinError(i)
940  << "\t" << errorFrac << endl;
941  }
942 
943 }
944 
945 //____________________
946 
947 void
948 showTGraphValues(TGraphAsymmErrors* graph,char* more)
949 {
950  *ofVal << "#################################################" << endl;
951  *ofVal << graph->GetName() << ":::::" << more << endl;
952 
953  const Int_t nBin = graph->GetN();
954  Double_t* xValues = graph->GetX();
955  Double_t* yValues = graph->GetY();
956  Double_t* errYHigh = graph->GetEYhigh();
957  Double_t* errXLow = graph->GetEXlow();
958  Double_t* errXHigh = graph->GetEXhigh();
959 
960  Double_t lowEdge, upEdge, errorFrac;
961 
962  for(Int_t i=0; i<nBin; i++){
963  lowEdge = xValues[i] - errXLow[i];
964  upEdge = xValues[i] + errXHigh[i];
965 
966  errorFrac = errYHigh[i]/yValues[i];
967 
968  *ofVal << "\t" << lowEdge << "<" << xValues[i] << "<" << upEdge
969  << "\t" << yValues[i] << "\t" << errYHigh[i] << "\t" << errorFrac
970  << endl;
971  }
972 }
973 // for each x slice, divide the y bins by the integral of the x slice
974 void scale2D(TH2* h2)
975 {
976  TAxis *xAxis=h2->GetXaxis();
977  TAxis *yAxis=h2->GetYaxis();
978 
979  for(int ix=1; ix<=xAxis->GetNbins(); ix++){
980  // first get the sum of this slice
981  float sum=0;
982  for(int iy=1; iy<=yAxis->GetNbins(); iy++){
983  int bin=h2->GetBin(ix,iy);
984  sum += h2->GetBinContent(bin);
985  }
986  // now scale
987  for(int iy=1; iy<=yAxis->GetNbins(); iy++){
988  int bin=h2->GetBin(ix,iy);
989  if(sum)
990  h2->SetBinContent(bin,h2->GetBinContent(bin)/sum);
991  }
992  }
993 }
Generated by   doxygen 1.8.5