ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/FGolf/Tools/histtools.cc
Revision: 1.3
Committed: Sat Jul 14 23:02:47 2012 UTC (12 years, 9 months ago) by fgolf
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +70 -0 lines
Log Message:
add function to combine hists, part of preliminary method to make stack plots

File Contents

# Content
1 #include "TCanvas.h"
2 #include "TFile.h"
3 #include "TH1.h"
4 #include "TH2.h"
5 #include "THStack.h"
6 #include "TKey.h"
7 #include "TLegend.h"
8 #include "TRegexp.h"
9 #include "histtools.h"
10 #include <TIterator.h>
11 #include "TClass.h"
12
13 #include <iostream>
14 #include <math.h>
15
16 using namespace std;
17
18 namespace hist {
19
20 void add(const char* outHistName, const char* patORpfx);
21
22 //Add all histograms whose names match one of ten possible regular expression
23 //patterns or begin with one of ten possible given prefixes. Feel free to
24 //mix and match regular expression patterns and prefixes. If the final hist-
25 //ogram named outHistName does not exist it is created.
26
27 void add(const char* outHistName, const char* patORpfx0, const char* patORpfx1, const char* patORpfx2, const char* patORpfx3, const char* patORpfx4, const char* patORpfx5, const char* patORpfx6, const char* patORpfx7, const char* patORpfx8, const char* patORpfx9)
28 {
29 add(outHistName, patORpfx0);
30 add(outHistName, patORpfx1);
31 if (patORpfx2) add(outHistName, patORpfx2);
32 if (patORpfx3) add(outHistName, patORpfx3);
33 if (patORpfx4) add(outHistName, patORpfx4);
34 if (patORpfx5) add(outHistName, patORpfx5);
35 if (patORpfx6) add(outHistName, patORpfx6);
36 if (patORpfx7) add(outHistName, patORpfx7);
37 if (patORpfx8) add(outHistName, patORpfx8);
38 if (patORpfx9) add(outHistName, patORpfx9);
39 }
40
41 //Add all histograms whose names match the given regular expression pattern
42 //or begin with the given prefix. If the final histogram named outHistName
43 //does not exist it is created.
44
45 void add(const char* outHistName, const char* patORpfx) {
46 TRegexp reg(patORpfx, kFALSE);
47
48 TList* list = gDirectory->GetList() ;
49 TIterator* iter = list->MakeIterator();
50
51 TObject* obj = 0;
52 TObject* hist = 0;
53 Bool_t makeOutHist = false;
54
55 hist = gDirectory->Get(outHistName);
56 //If out hist does not exist, remember to create it
57 if (! hist) makeOutHist = true;
58
59 while ( (obj = iter->Next()) ) {
60 if (! obj->InheritsFrom(TH1::Class())) continue;
61
62 TString name = obj->GetName();
63 //Don't add out hist
64 if (name == TString(outHistName)) continue;
65
66 if (TString(patORpfx).MaybeRegexp()) {
67 if (TString(obj->GetName()).Index(reg) < 0 ) continue;
68 } else if (! name.BeginsWith(patORpfx)) continue;
69
70 if (makeOutHist) {
71 hist = obj->Clone(outHistName);
72
73 if (hist->InheritsFrom(TH2::Class()))
74 ((TH2*)hist)->Reset();
75 else
76 ((TH1*)hist)->Reset();
77
78 ((TH1*)hist)->SetTitle(outHistName);
79 ((TH1*)hist)->Sumw2();
80 makeOutHist = false;
81 }
82
83 ((TH1*)hist)->Add((TH1*)obj);
84 }
85 }
86
87 //For all histograms whose names match the given regular expression pattern
88 //or begin with the given prefix, set the fill, line and marker colors to the
89 //given value.
90
91 void color(const char* patORpfx, Color_t color) {
92 TRegexp reg(patORpfx, kFALSE);
93
94 TList* list = gDirectory->GetList() ;
95 if (!list) {
96 cout << "Failed to set color for " << patORpfx << endl;
97 return;
98 }
99 TIterator* iter = list->MakeIterator();
100
101 TObject* obj = 0;
102
103 while ( (obj = iter->Next()) ) {
104 if (! obj->InheritsFrom(TH1::Class())) continue;
105
106 TString name = obj->GetName();
107
108 if (TString(patORpfx).MaybeRegexp()) {
109 if (TString(obj->GetName()).Index(reg) < 0 ) continue;
110 } else if (! name.BeginsWith(patORpfx)) continue;
111
112 ((TH1*)obj)->SetFillColor(color);
113 ((TH1*)obj)->SetLineColor(color);
114 ((TH1*)obj)->SetMarkerColor(color);
115 }
116 }
117
118 //Return a pointer to a TLegend with an entry for each histogram drawn on a
119 //given TCanvas. Display either the line, point or fill values. Optionally
120 //apply colors to all histograms. By default, entry labels are the names of
121 //their respective histograms. Optionally, if histogram names are of the
122 //form XX_YY_ZZ_WW, entry labels can be XX (token=0), YY (token=1), etc.
123
124 TLegend* legend(TCanvas* canvas, Option_t* option, Bool_t addColor, Int_t token,
125 Float_t xmin, Float_t ymin, Float_t xmax, Float_t ymax) {
126 if(! canvas) return 0;
127
128 TLegend* leg = new TLegend(xmin, ymin, xmax, ymax);
129 TList* list = canvas->GetListOfPrimitives();
130 TIterator* iter = list->MakeIterator();
131
132 TObject* obj = 0;
133
134 //Hist color iterator
135 Int_t colorIt = 1;
136
137 while ( (obj = iter->Next()) ) {
138 if (! obj->InheritsFrom(TH1::Class())) continue;
139
140 if (addColor) {
141 hist::color(obj->GetName(), colorIt);
142 ++colorIt;
143 }
144
145 if (token == -1)
146 leg->AddEntry(obj, obj->GetName(), option);
147 else {
148 TString name(obj->GetName());
149 TObjArray* a = name.Tokenize("_");
150 if (a->GetEntries() <= token)
151 leg->AddEntry(obj, obj->GetName(), option);
152 else
153 leg->AddEntry(obj, a->At(token)->GetName(), option);
154 }
155 }
156
157 return leg;
158 }
159
160 //Normalize to one all histograms whose names match the given regular exp-
161 //ression pattern or begin with the given prefix.
162
163 void normalize(const char* patORpfx) {
164 TRegexp reg(patORpfx, kFALSE);
165
166 TList* list = gDirectory->GetList() ;
167 TIterator* iter = list->MakeIterator();
168
169 TObject* obj = 0;
170
171 while ( (obj = iter->Next()) ) {
172 if (! obj->InheritsFrom(TH1::Class())) continue;
173
174 TString name = obj->GetName();
175
176 if (TString(patORpfx).MaybeRegexp()) {
177 if (TString(obj->GetName()).Index(reg) < 0 ) continue;
178 } else if (! name.BeginsWith(patORpfx)) continue;
179
180 Double_t integral = 0;
181
182 if (obj->InheritsFrom(TH2::Class()))
183 integral = ((TH2*)obj)->Integral();
184 else
185 integral = ((TH1*)obj)->Integral();
186
187 if (integral) {
188 ((TH1*)obj)->Sumw2();
189 ((TH1*)obj)->Scale(1./integral);
190 }
191 }
192 }
193
194 //Scale by the given value all histograms whose names match the given regular
195 //expression pattern or begin with the given prefix.
196
197 void scale(const char* patORpfx, Double_t scale) {
198 TRegexp reg(patORpfx, kFALSE);
199
200 TList* list = gDirectory->GetList() ;
201 TIterator* iter = list->MakeIterator();
202
203 TObject* obj = 0;
204
205 while ( (obj = iter->Next()) ) {
206 if (! obj->InheritsFrom(TH1::Class())) continue;
207
208 TString name = obj->GetName();
209
210 if (TString(patORpfx).MaybeRegexp()) {
211 if (TString(obj->GetName()).Index(reg) < 0 ) continue;
212 } else if (! name.BeginsWith(patORpfx)) continue;
213
214 if (!(((TH1*)obj)->GetSumw2N()))
215 ((TH1*)obj)->Sumw2();
216 ((TH1*)obj)->Scale(scale);
217 }
218 }
219
220 //Don't you hate it when you draw multiple histograms on the same canvas only
221 //to find that the bottom histogram's range does not encompass those of the
222 //histograms drawn on top? This method determines the maximum and minimum y
223 //range of all the histograms drawn on a given TCanvas and appropriately re-
224 //sizes the bottom histogram.
225
226 void setrangey(TCanvas* canvas) {
227 if(! canvas) return;
228
229 TList* list = canvas->GetListOfPrimitives();
230 TIterator* iter = list->MakeIterator();
231
232 TObject* obj = 0;
233 TObject* top = 0;
234
235 //Extremes
236 Double_t maxy = -999999;
237 Double_t miny = 999999;
238
239 while ( (obj = iter->Next()) ) {
240 if (! obj->InheritsFrom(TH1::Class())) continue;
241
242 if (! top) top = obj;
243
244 if (((TH1*)obj)->GetMaximum() > maxy) maxy = ((TH1*)obj)->GetMaximum();
245 if (((TH1*)obj)->GetMinimum() < miny) miny = ((TH1*)obj)->GetMinimum();
246 }
247
248 ((TH1*)top)->SetMaximum(maxy*1.9);
249 //Protect against log scale
250 if (canvas->GetLogy() && ! miny)
251 ((TH1*)top)->SetMinimum(1E-4);
252 else
253 ((TH1*)top)->SetMinimum(miny*0.7);
254 }
255
256 //Create a stacked histogram consisting of all histograms whose names match
257 //the given regular expression pattern or begin with the given prefix. If
258 //the THStack named stackHistName does not exist it is created. Optionally
259 //apply colors to all histograms. Set drawOption to "nostack" if you do not
260 //want to stack, to "hist" to display histograms without errors, to "histe"
261 //to display histograms with errors, etc.
262
263 void stack(const char* stackHistName, const char* patORpfx, Bool_t addColor, Option_t* drawOption) {
264 TRegexp reg(patORpfx, kFALSE);
265
266 TList* list = gDirectory->GetList() ;
267 TIterator* iter = list->MakeIterator();
268
269 TObject* obj = 0;
270 TObject* stack = 0;
271 Bool_t makeStackHist = false;
272
273 stack = gDirectory->Get(stackHistName);
274 //If stack hist does not exist, remember to create it
275 if (! stack) makeStackHist = true;
276
277 //Hist color iterator
278 Int_t colorIt = 1;
279
280 while ( (obj = iter->Next()) ) {
281 if (! obj->InheritsFrom(TH1::Class())) continue;
282
283 TString name = obj->GetName();
284
285 if (TString(patORpfx).MaybeRegexp()) {
286 if (TString(obj->GetName()).Index(reg) < 0 ) continue;
287 } else if (! name.BeginsWith(patORpfx)) continue;
288
289 if (makeStackHist) {
290 stack = new THStack(stackHistName, stackHistName);
291 makeStackHist = false;
292 }
293
294 if (addColor) {
295 hist::color(obj->GetName(), colorIt);
296 ++colorIt;
297 }
298
299 ((THStack*)stack)->Add((TH1*)obj, drawOption);
300 }
301
302 // Currently breaks .ls
303 //gDirectory->Append(stack);
304 }
305
306 //Set the x-axis title of all histograms whose names match the given regular
307 //expression pattern or begin with the given prefix.
308
309 void xaxis(const char* patORpfx, const char* title) {
310 TRegexp reg(patORpfx, kFALSE);
311
312 TList* list = gDirectory->GetList() ;
313 TIterator* iter = list->MakeIterator();
314
315 TObject* obj = 0;
316
317 while ( (obj = iter->Next()) ) {
318 if (! (obj->InheritsFrom(TH1::Class()) || obj->InheritsFrom(THStack::Class()))) continue;
319
320 TString name = obj->GetName();
321
322 if (TString(patORpfx).MaybeRegexp()) {
323 if (TString(obj->GetName()).Index(reg) < 0 ) continue;
324 } else if (! name.BeginsWith(patORpfx)) continue;
325
326 if (obj->InheritsFrom(TH1::Class()))
327 ((TH1*)obj)->GetXaxis()->SetTitle(title);
328 if (obj->InheritsFrom(THStack::Class())) {
329 ((THStack*)obj)->Draw();
330 ((THStack*)obj)->GetXaxis()->SetTitle(title);
331 }
332 }
333 }
334
335 //Set the y-axis title of all histograms whose names match the given regular
336 //expression pattern or begin with the given prefix.
337
338 void yaxis(const char* patORpfx, const char* title) {
339 TRegexp reg(patORpfx, kFALSE);
340
341 TList* list = gDirectory->GetList() ;
342 TIterator* iter = list->MakeIterator();
343
344 TObject* obj = 0;
345
346 while ( (obj = iter->Next()) ) {
347 if (! (obj->InheritsFrom(TH1::Class()) || obj->InheritsFrom(THStack::Class()))) continue;
348
349 TString name = obj->GetName();
350
351 if (TString(patORpfx).MaybeRegexp()) {
352 if (TString(obj->GetName()).Index(reg) < 0 ) continue;
353 } else if (! name.BeginsWith(patORpfx)) continue;
354
355 if (obj->InheritsFrom(TH1::Class()))
356 ((TH1*)obj)->GetYaxis()->SetTitle(title);
357 if (obj->InheritsFrom(THStack::Class())) {
358 ((THStack*)obj)->Draw();
359 ((THStack*)obj)->GetYaxis()->SetTitle(title);
360 }
361 }
362 }
363
364 // Input: 4 histogram
365 // Output: one histogram which is the BG subtracted efficiency:
366 // h1 : TOTAL NUMBER OF EVENTS, SIGNAL REGION
367 // h2 : NUMBER OF EVENTS THAT PASS, SIGNAL REGION
368 // h3 : TOTAL NUMBER OF EVENTS, SIDE BAND
369 // h4 : NUMBER OF EVENTS THAT PASS, SIDE BAND
370
371
372 TH1F* eff_bg(TH1F* h1, TH1F* h2, TH1F* h3, TH1F* h4, const char* name){
373
374 // first, verify that all histograms have same binning
375 // nx is the number of visible bins
376 // nxtot = nx+2 includes underflow and overflow
377 Int_t nx = h1->GetNbinsX();
378 Int_t nxtot = nx + 2;
379 if (h2->GetNbinsX() != nx) {
380 cout << "Histograms must have same number of bins" << endl;
381 return 0;
382 }
383 if (h3->GetNbinsX() != nx) {
384 cout << "Histograms must have same number of bins" << endl;
385 return 0;
386 }
387 if (h4->GetNbinsX() != nx) {
388 cout << "Histograms must have same number of bins" << endl;
389 return 0;
390 }
391
392 // get the new histogram
393 TH1F* temp = (TH1F*) h1->Clone(name);
394 temp->SetTitle(name);
395 temp->Reset();
396 temp->Sumw2();
397
398 // Loop over bins, calculate efficiency and error, put it in histogram
399 for (Int_t i=0; i<nxtot; i++) {
400 Double_t x1 = h1->GetBinContent(i);
401 Double_t x2 = h2->GetBinContent(i);
402 Double_t x3 = h3->GetBinContent(i);
403 Double_t x4 = h4->GetBinContent(i);
404 Double_t denom = x1 - x3;
405 Double_t eff;
406 if (denom == 0.) {
407 eff = 0;
408 } else {
409 eff = (x2-x4)/denom;
410 }
411 Double_t failSig = x1 - x2;
412 Double_t failBg = x3 - x4;
413 Double_t blah = (1-eff)*(1-eff)*(x2+x4) + eff*eff*(failSig+failBg);
414 if (blah <= 0.) blah=0.0;
415 Double_t err;
416 if (denom == 0) {
417 err = 0.;
418 } else {
419 err = sqrt(blah)/denom;
420 }
421 temp->SetBinContent(i,eff);
422 temp->SetBinError(i,err);
423 }
424
425 // Done
426 return temp;
427 }
428
429 void deleteHistos() {
430 // Delete all existing histograms in memory
431 TObject* obj;
432 TList* list = gDirectory->GetList() ;
433 TIterator* iter = list->MakeIterator();
434 while ( (obj=iter->Next()) ) {
435 if (obj->IsA()->InheritsFrom(TH1::Class()) ||
436 obj->IsA()->InheritsFrom(TH2::Class()) ) {delete obj;}
437 }
438 }
439
440 void histio()
441 {
442 }
443
444 void saveHist(const char* filename, const char* pat)
445 {
446 TList* list = gDirectory->GetList() ;
447 TIterator* iter = list->MakeIterator();
448
449 TRegexp re(pat,kTRUE) ;
450
451 TFile outf(filename,"RECREATE") ;
452 while(TObject *obj=iter->Next()) {
453 if (TString(obj->GetName()).Index(re)>=0) {
454 obj->Write() ;
455 cout << "." ;
456 cout.flush() ;
457 }
458 }
459 cout << endl ;
460 outf.Close() ;
461
462 delete iter ;
463 }
464
465
466 void loadHist(const char* filename, const char* pfx, const char* pat, Bool_t verbose, Bool_t doAdd)
467 {
468 TFile inf(filename) ;
469 //inf.ReadAll() ;
470 TList* list = inf.GetListOfKeys() ;
471 TIterator* iter = list->MakeIterator();
472
473 TRegexp re(pat,kTRUE) ;
474 if (verbose)
475 std::cout << "pat = " << pat << std::endl;
476
477 gDirectory->cd("Rint:") ;
478
479 TObject* obj ;
480 TKey* key ;
481 if (verbose) {
482 std::cout << "doAdd = " << (doAdd?"T":"F") << std::endl;
483 std::cout << "loadHist: reading.";
484 }
485 while( (key=(TKey*)iter->Next()) ) {
486
487 Int_t ridx = TString(key->GetName()).Index(re) ;
488 if (ridx==-1) {
489 continue ;
490 }
491
492 obj = inf.Get(key->GetName()) ;
493 TObject* clone ;
494 if (pfx) {
495
496 // Find existing TH1-derived objects
497 TObject* oldObj = 0 ;
498 if (doAdd){
499 oldObj = gDirectory->Get(Form("%s_%s",pfx,obj->GetName())) ;
500 if (oldObj && !oldObj->IsA()->InheritsFrom(TH1::Class())) {
501 oldObj = 0 ;
502 }
503 }
504 if (oldObj) {
505 clone = oldObj ;
506 ((TH1*)clone)->Add((TH1*)obj) ;
507 } else {
508 clone = obj->Clone(Form("%s_%s",pfx,obj->GetName())) ;
509 }
510
511
512 } else {
513
514 // Find existing TH1-derived objects
515 TObject* oldObj = 0 ;
516 if (doAdd){
517 oldObj = gDirectory->Get(key->GetName()) ;
518 if (oldObj && !oldObj->IsA()->InheritsFrom(TH1::Class())) {
519 oldObj = 0 ;
520 }
521 }
522
523 if (oldObj) {
524 clone = oldObj ;
525 ((TH1*)clone)->Add((TH1*)obj) ;
526 } else {
527 clone = obj->Clone() ;
528 }
529 }
530 if (!gDirectory->GetList()->FindObject(clone)) {
531 gDirectory->Append(clone) ;
532 }
533 if (verbose) {
534 std::cout << ".";
535 std::cout.flush();
536 }
537 }
538 if (verbose)
539 std::cout << std::endl;
540 inf.Close() ;
541 delete iter ;
542 }
543
544 //takes in a vector of prefixes, and a output prefix. Combines all histograms
545 //whose prefixes matches those in the input vector of TStrings into a new histogram
546 //with the the outPrefix
547 void combineHists(const std::vector<TString> v_prfxsToCombine, const TString outPrefix) {
548
549
550 if(v_prfxsToCombine.size() == 0) {
551 cout << "Input vector size is 0" << endl;
552 }
553 //histogram of common suffixes
554 vector<TString> v_suffixes;
555 //get the list of histograms that match the first entry in the prefix
556 TString temp = v_prfxsToCombine.at(0);
557 TList *rootlist = gDirectory->GetList();
558
559 for(int i = 0; i < rootlist->GetSize(); i++) {
560
561 TObject *obj = rootlist->At(i);
562 TString name = obj->GetName();
563 if(!obj->InheritsFrom(TH1::Class()))
564 continue;
565 if(!name.BeginsWith(temp+"_"))
566 continue;
567
568 name.ReplaceAll(temp+"_", "_");
569
570
571 if(obj->InheritsFrom(TH2::Class())) {
572 TH2F *h = dynamic_cast<TH2F*>(obj->Clone());
573 h->SetName((outPrefix+name).Data());
574 h->SetTitle((outPrefix+name).Data());
575 for(unsigned int j = 1; j < v_prfxsToCombine.size(); j++) {
576 TH2F *htemp = dynamic_cast<TH2F*>(gDirectory->Get((v_prfxsToCombine.at(j) + name).Data()));
577 h->Add(htemp);
578 }
579 } else if(obj->InheritsFrom(TH1::Class())) {
580 TH1F *h = dynamic_cast<TH1F*>(obj->Clone());
581 h->SetName((outPrefix+name).Data());
582 h->SetTitle((outPrefix+name).Data());
583 for(unsigned int j = 1; j < v_prfxsToCombine.size(); j++) {
584 TH1F *htemp = dynamic_cast<TH1F*>(gDirectory->Get((v_prfxsToCombine.at(j) + name).Data()));
585 //cout << "TH1F: " << v_prfxsToCombine.at(j) + name << endl;
586 if(htemp == NULL)
587 cout << "htemp is NULL" << endl;
588 h->Add(htemp);
589 }
590 }
591 }//rootlist loop
592
593
594 //now delete the histos that match the prfxs to combine
595 for(unsigned int i = 0; i < v_prfxsToCombine.size(); i++ ) {
596
597 // Delete all existing histograms in memory
598 TObject* obj;
599 TList* list = gDirectory->GetList() ;
600 TIterator* iter = list->MakeIterator();
601 while ((obj=iter->Next())) {
602 TString name = obj->GetName();
603 if(name.BeginsWith(outPrefix+"_"))
604 continue;
605 if(TString(obj->GetName()).BeginsWith(v_prfxsToCombine.at(i).Data())) {
606 if (obj->IsA()->InheritsFrom(TH1::Class()) ||
607 obj->IsA()->InheritsFrom(TH2::Class()) )
608 delete obj;
609 }
610 }
611 }//loop over prefixes
612
613 }//fnc declaration
614 }