ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/FGolf/Tools/histtools.C
Revision: 1.1
Committed: Sun Aug 29 18:39:12 2010 UTC (14 years, 8 months ago) by fgolf
Content type: text/plain
Branch: MAIN
CVS Tags: synchMay2011v1, ss20May2011, ss2010_10_06, fkwsynch3pbv2, fkw_3pb_samesign_synch
Log Message:
some more useful tools for a variety of studies

File Contents

# User Rev Content
1 fgolf 1.1 #ifndef HISTTOOLS_H
2     #define HISTTOOLS_H
3    
4     #include "TList.h"
5     #include "TObjArray.h"
6     #include "TH1.h"
7     #include "TH2.h"
8     #include "TLegend.h"
9     #include "THStack.h"
10     #include "TCanvas.h"
11     #include "TFile.h"
12     #include "TRegexp.h"
13     #include "TKey.h"
14     #include <iostream>
15     #include <vector>
16    
17     #include "histtools.h"
18     namespace hist {
19    
20     //Add all histograms whose names match the given regular expression pattern
21     //or begin with the given prefix. If the final histogram named outHistName
22     //does not exist it is created.
23    
24     void add(const char* outHistName, const char* patORpfx) {
25     TRegexp reg(patORpfx, kFALSE);
26    
27     TList* list = gDirectory->GetList() ;
28     TIterator* iter = list->MakeIterator();
29    
30     TObject* obj = 0;
31     TObject* hist = 0;
32     Bool_t makeOutHist = false;
33    
34     hist = gDirectory->Get(outHistName);
35     //If out hist does not exist, remember to create it
36     if (! hist) makeOutHist = true;
37    
38     while ((obj = iter->Next())) {
39     if (! obj->InheritsFrom(TH1::Class())) continue;
40    
41     TString name = obj->GetName();
42     //Don't add out hist
43     if (name == TString(outHistName)) continue;
44    
45     if (TString(patORpfx).MaybeRegexp()) {
46     if (TString(obj->GetName()).Index(reg) < 0 ) continue;
47     } else if (! name.BeginsWith(patORpfx)) continue;
48    
49     if (makeOutHist) {
50     hist = obj->Clone(outHistName);
51    
52     if (hist->InheritsFrom(TH2::Class()))
53     ((TH2*)hist)->Reset();
54     else
55     ((TH1*)hist)->Reset();
56    
57     ((TH1*)hist)->SetTitle(outHistName);
58     ((TH1*)hist)->Sumw2();
59     makeOutHist = false;
60     }
61    
62     ((TH1*)hist)->Add((TH1*)obj);
63     }
64     }
65    
66     //Add all histograms whose names match one of ten possible regular expression
67     //patterns or begin with one of ten possible given prefixes. Feel free to
68     //mix and match regular expression patterns and prefixes. If the final hist-
69     //ogram named outHistName does not exist it is created.
70    
71     void add(const char* outHistName, const char* patORpfx0, const char* patORpfx1, const char* patORpfx2,
72     const char* patORpfx3, const char* patORpfx4, const char* patORpfx5, const char* patORpfx6,
73     const char* patORpfx7, const char* patORpfx8, const char* patORpfx9)
74     {
75     add(outHistName, patORpfx0);
76     add(outHistName, patORpfx1);
77     if (patORpfx2) add(outHistName, patORpfx2);
78     if (patORpfx3) add(outHistName, patORpfx3);
79     if (patORpfx4) add(outHistName, patORpfx4);
80     if (patORpfx5) add(outHistName, patORpfx5);
81     if (patORpfx6) add(outHistName, patORpfx6);
82     if (patORpfx7) add(outHistName, patORpfx7);
83     if (patORpfx8) add(outHistName, patORpfx8);
84     if (patORpfx9) add(outHistName, patORpfx9);
85     }
86    
87    
88     //For all histograms whose names match the given regular expression pattern
89     //or begin with the given prefix, set the fill, line and marker colors to the
90     //given value.
91    
92     void color(const char* patORpfx, Color_t color) {
93    
94     TRegexp reg(patORpfx, kFALSE);
95    
96     TList* list = gDirectory->GetList() ;
97    
98     TIterator* iter = list->MakeIterator();
99    
100     TObject* obj = 0;
101    
102     while ((obj = iter->Next())) {
103     if (! obj->InheritsFrom(TH1::Class())) continue;
104    
105     TString name = obj->GetName();
106    
107     if (TString(patORpfx).MaybeRegexp()) {
108     if (TString(obj->GetName()).Index(reg) < 0 ) continue;
109     } else if (! name.BeginsWith(patORpfx)) continue;
110    
111     ((TH1*)obj)->SetFillColor(color);
112     ((TH1*)obj)->SetLineColor(color);
113     ((TH1*)obj)->SetMarkerColor(color);
114     }
115     }
116    
117    
118     void pattern(const char* patORpfx, Color_t pattern) {
119     TRegexp reg(patORpfx, kFALSE);
120    
121     TList* list = gDirectory->GetList() ;
122     TIterator* iter = list->MakeIterator();
123    
124     TObject* obj = 0;
125    
126     while ((obj = iter->Next())) {
127     if (! obj->InheritsFrom(TH1::Class())) continue;
128    
129     TString name = obj->GetName();
130    
131     if (TString(patORpfx).MaybeRegexp()) {
132     if (TString(obj->GetName()).Index(reg) < 0 ) continue;
133     } else if (! name.BeginsWith(patORpfx)) continue;
134    
135     ((TH1*)obj)->SetFillStyle(pattern);
136     }
137     }
138    
139    
140    
141    
142     //Return a pointer to a TLegend with an entry for each histogram drawn on a
143     //given TCanvas. Display either the line, point or fill values. Optionally
144     //apply colors to all histograms. By default, entry labels are the names of
145     //their respective histograms. Optionally, if histogram names are of the
146     //form XX_YY_ZZ_WW, entry labels can be XX (token=0), YY (token=1), etc.
147    
148     TLegend* legend(TCanvas* canvas, Option_t* option, Bool_t addColor, Int_t token,
149     Float_t xmin, Float_t ymin, Float_t xmax, Float_t ymax) {
150     if(! canvas) return 0;
151    
152     TLegend* leg = new TLegend(xmin, ymin, xmax, ymax);
153     TList* list = canvas->GetListOfPrimitives();
154     // TIterator* iter = list->MakeIterator();
155    
156     TObject* obj = 0;
157    
158     //Hist color iterator
159     Int_t colorIt = 1;
160    
161     //while (obj = iter->Next()) {
162     for(int i = list->GetSize() - 1 ; i >=0 ; i--) {
163     obj = list->At(i);
164     if (! obj->InheritsFrom(TH1::Class())) continue;
165    
166     if (addColor) {
167     hist::color(obj->GetName(), colorIt);
168     ++colorIt;
169     }
170    
171    
172     if (token == -1)
173     leg->AddEntry(obj, obj->GetName(), option);
174     else {
175     TString name(obj->GetName());
176     TObjArray* a = name.Tokenize("_");
177     if (a->GetEntries() <= token)
178     leg->AddEntry(obj, obj->GetName(), option);
179     else
180     leg->AddEntry(obj, a->At(token)->GetName(), option);
181     }
182     }
183    
184     return leg;
185     }
186    
187     //Return a pointer to a TLegend with an entry for each histogram added to a
188     //given THStack. Display either the line, point or fill values. Optionally
189     //apply colors to all histograms. By default, entry labels are the names of
190     //their respective histograms. Optionally, if histogram names are of the
191     //form XX_YY_ZZ_WW, entry labels can be XX (token=0), YY (token=1), etc.
192    
193     TLegend* legend(THStack* stack, Option_t* option, Bool_t addColor, Int_t token,
194     Float_t xmin, Float_t ymin, Float_t xmax, Float_t ymax,
195     TH1F *hdata) {
196    
197     if(! stack) return 0;
198    
199     TLegend* leg = new TLegend(xmin, ymin, xmax, ymax);
200     leg->SetFillColor(kWhite);
201     TList* list = stack->GetHists();
202     // TIterator* iter = list->MakeIterator();
203    
204     TObject* obj = 0;
205    
206    
207     if(hdata != NULL)
208     leg->AddEntry(hdata, "Data", "P");
209    
210    
211     //Hist color iterator
212     Int_t colorIt = 1;
213    
214     //while (obj = iter->Next()) {
215     for(int i = list->GetSize() - 1 ; i >=0 ; i--) {
216     obj = list->At(i);
217     if (! obj->InheritsFrom(TH1::Class())) continue;
218     TString name(obj->GetName());
219    
220     if(name.Contains("QCD"))
221     continue;
222    
223     if (addColor) {
224     hist::color(obj->GetName(), colorIt);
225     ++colorIt;
226     }
227    
228     if (token == -1)
229     leg->AddEntry(obj, obj->GetName(), option);
230     else {
231     TString name(obj->GetName());
232     TObjArray* a = name.Tokenize("_");
233     TString entry(a->GetEntries() <= token ? obj->GetName() : a->At(token)->GetName());
234    
235     if (entry == "t"){
236     entry = "Single top";
237     }else if (entry == "ttbar"){
238     entry = "#font[12]{t#bar{t}}";
239     } else if (entry == "ttdil"){
240     entry = "#font[12]{t#bar{t}} signal";
241     } else if (entry == "ttotr"){
242     entry = "#font[12]{t#bar{t}} other";
243     } else if (entry == "VV"){
244     //entry = "#font[12]{ZZ}, #font[12]{WZ}, #font[12]{ZZ}";
245     entry = "#font[12]{VV}";
246     } else if (entry == "DYtautau"){
247     //entry = "DY#rightarrow #tau#tau";
248     entry = "Z/#gamma*#rightarrow#tau^{+}#tau^{-}";
249     } else if (entry == "DYeemm"){
250     //entry = "DY#rightarrow #font[12]{ee}, #mu#mu";
251     entry = "Z/#gamma*#rightarrowl^{+}l^{-}";
252     } else if (entry == "wjets"){
253     //entry = "#font[12]{W}+jets";
254     entry = "#font[12]{W}#rightarrowl#nu";
255     }else {
256     entry.ReplaceAll("tautau", "#tau#tau");
257     entry.ReplaceAll("mm", "#mu#mu");
258     }
259     leg->AddEntry(obj, entry, option);
260     }
261     }
262    
263     return leg;
264     }
265    
266     //Normalize to one all histograms whose names match the given regular exp-
267     //ression pattern or begin with the given prefix.
268    
269     void normalize(const char* patORpfx) {
270     TRegexp reg(patORpfx, kFALSE);
271    
272     TList* list = gDirectory->GetList() ;
273     TIterator* iter = list->MakeIterator();
274    
275     TObject* obj = 0;
276    
277     while ((obj = iter->Next())) {
278     if (! obj->InheritsFrom(TH1::Class())) continue;
279    
280     TString name = obj->GetName();
281    
282     if (TString(patORpfx).MaybeRegexp()) {
283     if (TString(obj->GetName()).Index(reg) < 0 ) continue;
284     } else if (! name.BeginsWith(patORpfx)) continue;
285    
286     Double_t integral = 0;
287    
288     if (obj->InheritsFrom(TH2::Class()))
289     integral = ((TH2*)obj)->Integral();
290     else
291     integral = ((TH1*)obj)->Integral();
292    
293     if (integral) {
294     ((TH1*)obj)->Sumw2();
295     ((TH1*)obj)->Scale(1./integral);
296     }
297     }
298     }
299    
300     //Scale by the given value all histograms whose names match the given regular
301     //expression pattern or begin with the given prefix.
302    
303     void scale(const char* patORpfx, Double_t scale) {
304     TRegexp reg(patORpfx, kFALSE);
305    
306     TList* list = gDirectory->GetList() ;
307     TIterator* iter = list->MakeIterator();
308    
309     TObject* obj = 0;
310    
311     while ((obj = iter->Next())) {
312     if (! obj->InheritsFrom(TH1::Class())) continue;
313    
314     TString name = obj->GetName();
315    
316     if (TString(patORpfx).MaybeRegexp()) {
317     if (TString(obj->GetName()).Index(reg) < 0 ) continue;
318     } else if (! name.BeginsWith(patORpfx)) continue;
319    
320    
321     if ((((TH1*)obj)->GetNbinsX()+2) != ((TH1*)obj)->GetSumw2N())
322     ((TH1*)obj)->Sumw2();
323    
324     ((TH1*)obj)->Scale(scale);
325     }
326     }
327    
328     //Don't you hate it when you draw multiple histograms on the same canvas only
329     //to find that the bottom histogram's range does not encompass those of the
330     //histograms drawn on top? This method determines the maximum and minimum y
331     //range of all the histograms drawn on a given TCanvas and appropriately re-
332     //sizes the bottom histogram.
333    
334     void setrangey(TCanvas* canvas) {
335     if(! canvas) return ;
336    
337     TList* list = canvas->GetListOfPrimitives();
338     TIterator* iter = list->MakeIterator();
339    
340     TObject* obj = 0;
341     TObject* top = 0;
342    
343     //Extremes
344     Double_t maxy = -999999;
345     Double_t miny = 999999;
346    
347     while ((obj = iter->Next())) {
348     if (! obj->InheritsFrom(TH1::Class())) continue;
349    
350     if (! top) top = obj;
351    
352     if (((TH1*)obj)->GetMaximum() > maxy) maxy = ((TH1*)obj)->GetMaximum();
353     if (((TH1*)obj)->GetMinimum() < miny) miny = ((TH1*)obj)->GetMinimum();
354     }
355    
356     ((TH1*)top)->SetMaximum(maxy*1.3);
357     //Protect against log scale
358     if (canvas->GetLogy() && ! miny)
359     ((TH1*)top)->SetMinimum(1E-3);
360     else
361     ((TH1*)top)->SetMinimum(miny*0.7);
362     }
363    
364     //Create a stacked histogram consisting of all histograms whose names match
365     //the given regular expression pattern or begin with the given prefix. If
366     //the THStack named stackHistName does not exist it is created. Optionally
367     //apply colors to all histograms. Set drawOption to "nostack" if you do not
368     //want to stack, to "hist" to display histograms without errors, to "histe"
369     //to display histograms with errors, etc.
370    
371     void stack(const char* stackHistName, const char* patORpfx, Bool_t addColor, Option_t* drawOption, Int_t orderScheme,
372     const char* bsmName, bool doRefPats, std::string prefixToSkip){//, float yMinimum) {
373     TRegexp reg(patORpfx, kTRUE);
374    
375     TList* list = gDirectory->GetList() ;
376    
377     TObject* obj = 0;
378     TObject* stack = 0;
379     Bool_t makeStackHist = false;
380    
381     stack = gDirectory->Get(stackHistName);
382     //If stack hist does not exist, remember to create it
383     if (! stack) makeStackHist = true;
384    
385     //Hist color iterator
386     Int_t colorIt = 1;
387    
388     std::vector<TString*> pats(0);
389     if (orderScheme == 0){
390     pats.push_back(new TString(""));
391     //pats.push_back(new TString(""));
392     //pats.push_back(new TString("ttdil_"));
393     //pats.push_back(new TString("ttotr_"));
394     //pats.push_back(new TString("DYeemm_"));
395     //pats.push_back(new TString("DYtautau_"));
396     //pats.push_back(new TString("wjets_"));
397     //pats.push_back(new TString("VV_"));
398     //pats.push_back(new TString("tw_"));
399     //pats.push_back(new TString("data_"));
400     } else if (orderScheme == 1) {
401     //this requires some care: only the histograms beginning with what's in the pattern will be processed
402     pats.push_back(new TString("ttdil_"));
403     pats.push_back(new TString("ttotr_"));
404     pats.push_back(new TString("t_"));
405     pats.push_back(new TString("VV_"));
406     pats.push_back(new TString("DYtautau_"));
407     pats.push_back(new TString("DYeemm_"));
408     pats.push_back(new TString("wjets_"));
409     pats.push_back(new TString("QCD_"));
410     // pats.push_back(new TString("Vgamma_"));
411    
412     if (doRefPats){
413     pats.push_back(new TString("ref_ttdil_"));
414     pats.push_back(new TString("ref_ttotr_"));
415     pats.push_back(new TString("ref_t_"));
416     pats.push_back(new TString("ref_VV_"));
417     pats.push_back(new TString("ref_DYtautau_"));
418     pats.push_back(new TString("ref_DYeemm_"));
419     pats.push_back(new TString("ref_wjets_"));
420     pats.push_back(new TString("ref_QCD_"));
421     // pats.push_back(new TString("ref_Vgamma_"));
422     }
423     } else if (orderScheme == 2) {
424     //this requires some care: only the histograms beginning with what's in the pattern will be processed
425     pats.push_back(new TString("ttdil_"));
426     pats.push_back(new TString("ttotr_"));
427     pats.push_back(new TString("t_"));
428     pats.push_back(new TString("VV_"));
429     pats.push_back(new TString("DYtautau_"));
430     pats.push_back(new TString("DYeemm_"));
431     pats.push_back(new TString("wjets_"));
432     pats.push_back(new TString("QCD_"));
433     pats.push_back(new TString("Vgamma_"));
434     if (bsmName!=0) pats.push_back(new TString(Form("%s_",bsmName)));
435    
436     if (doRefPats){
437     pats.push_back(new TString("ref_ttdil_"));
438     pats.push_back(new TString("ref_ttotr_"));
439     pats.push_back(new TString("ref_t_"));
440     pats.push_back(new TString("ref_VV_"));
441     pats.push_back(new TString("ref_DYtautau_"));
442     pats.push_back(new TString("ref_DYeemm_"));
443     pats.push_back(new TString("ref_wjets_"));
444     pats.push_back(new TString("ref_QCD_"));
445     pats.push_back(new TString("ref_Vgamma_"));
446     if(bsmName!=0) pats.push_back(new TString(Form("ref_%s_",bsmName)));
447     }
448     } else {
449     pats.push_back(new TString(""));
450     }
451    
452     for (unsigned int iPat = 0; iPat < pats.size(); ++iPat) {
453     TIterator* iter = list->MakeIterator();
454     while ((obj = iter->Next())) {
455     if (! obj->InheritsFrom(TH1::Class())) continue;
456     //((TH1F*)obj)->SetMinimum(yMinimum);
457     TString name = obj->GetName();
458     if(prefixToSkip != "" && name.BeginsWith(prefixToSkip.c_str())) continue;
459     //if(name.Contains("data"))
460     if (!( name.BeginsWith(pats[iPat]->Data()))) continue;
461    
462     if (TString(patORpfx).MaybeRegexp()) {
463     if (TString(obj->GetName()).Index(reg) < 0 ) continue;
464     } else if (! name.BeginsWith(patORpfx)) continue;
465    
466     if (makeStackHist) {
467     stack = new THStack(stackHistName, stackHistName);
468     makeStackHist = false;
469     }
470    
471     if (addColor) {
472     hist::color(obj->GetName(), colorIt);
473     ++colorIt;
474     }
475     ((THStack*)stack)->Add((TH1*)obj, drawOption);
476    
477     }
478     // delete pats[iPat];//cleanup
479     }
480    
481    
482     // Currently breaks .ls
483     //gDirectory->Append(stack);
484     }
485    
486     //Set the x-axis title of all histograms whose names match the given regular
487     //expression pattern or begin with the given prefix.
488    
489     void xaxis(const char* patORpfx, const char* title) {
490     TRegexp reg(patORpfx, kFALSE);
491    
492     TList* list = gDirectory->GetList() ;
493     TIterator* iter = list->MakeIterator();
494    
495     TObject* obj = 0;
496    
497     while ((obj = iter->Next())) {
498     if (! (obj->InheritsFrom(TH1::Class()) || obj->InheritsFrom(THStack::Class()))) continue;
499    
500     TString name = obj->GetName();
501    
502     if (TString(patORpfx).MaybeRegexp()) {
503     if (TString(obj->GetName()).Index(reg) < 0 ) continue;
504     } else if (! name.BeginsWith(patORpfx)) continue;
505    
506     if (obj->InheritsFrom(TH1::Class()))
507     ((TH1*)obj)->GetXaxis()->SetTitle(title);
508     if (obj->InheritsFrom(THStack::Class())) {
509     ((THStack*)obj)->Draw();
510     ((THStack*)obj)->GetXaxis()->SetTitle(title);
511     }
512     }
513     }
514    
515     //Set the y-axis title of all histograms whose names match the given regular
516     //expression pattern or begin with the given prefix.
517    
518     void yaxis(const char* patORpfx, const char* title) {
519     TRegexp reg(patORpfx, kFALSE);
520    
521     TList* list = gDirectory->GetList() ;
522     TIterator* iter = list->MakeIterator();
523    
524     TObject* obj = 0;
525    
526     while ((obj = iter->Next())) {
527     if (! (obj->InheritsFrom(TH1::Class()) || obj->InheritsFrom(THStack::Class()))) continue;
528    
529     TString name = obj->GetName();
530    
531     if (TString(patORpfx).MaybeRegexp()) {
532     if (TString(obj->GetName()).Index(reg) < 0 ) continue;
533     } else if (! name.BeginsWith(patORpfx)) continue;
534    
535     if (obj->InheritsFrom(TH1::Class()))
536     ((TH1*)obj)->GetYaxis()->SetTitle(title);
537     if (obj->InheritsFrom(THStack::Class())) {
538     ((THStack*)obj)->Draw();
539     ((THStack*)obj)->GetYaxis()->SetTitle(title);
540     }
541     }
542     }
543    
544     // Input: 2 histogram
545     // Output: one histogram which is the efficiency:
546     // h1 : TOTAL NUMBER OF EVENTS
547     // h2 : NUMBER OF EVENTS THAT PASS
548    
549     #include "TH1.h"
550    
551     // Method by pointer
552     TH1F* eff(TH1F* h1, TH1F* h2, const char* name){
553    
554     // first, verify that all histograms have same binning
555     // nx is the number of visible bins
556     // nxtot = nx+2 includes underflow and overflow
557     Int_t nx = h1->GetNbinsX();
558     if (h2->GetNbinsX() != nx) {
559     cout << "Histograms must have same number of bins" << endl;
560     return 0;
561     }
562    
563     // get the new histogram
564     TH1F* temp = (TH1F*) h1->Clone(name);
565     temp->SetTitle(name);
566     temp->Reset();
567     temp->Sumw2();
568    
569     // Do the calculation
570     temp->Divide(h2,h1,1.,1.,"B");
571    
572     // Done
573     return temp;
574     }
575    
576    
577     // Method by name
578     TH1F* eff(const char* name1, const char* name2, const char* name){
579    
580     // Get a list of object and their iterator
581     TList* list = gDirectory->GetList() ;
582     TIterator* iter = list->MakeIterator();
583    
584     // Loop over objects, set the pointers
585     TObject* obj;
586     TH1F* h1=0;
587     TH1F* h2=0;
588     TString str1 = Form("%s",name1);
589     TString str2 = Form("%s",name2);
590     while((obj=iter->Next())) {
591     TString objName = obj->GetName();
592     if (objName == str1) h1 = (TH1F*) obj;
593     if (objName == str2) h2 = (TH1F*) obj;
594     }
595    
596     // quit if not found
597     if (h1 == 0) {
598     cout << "Histogram " << name1 << " not found" << endl;
599     return 0;
600     }
601     if (h2 == 0) {
602     cout << "Histogram " << name2 << " not found" << endl;
603     return 0;
604     }
605    
606     // Call the method by pointer
607     TH1F* temp = eff(h1, h2, name);
608     return temp;
609     }
610     // Input: 4 histogram
611     // Output: one histogram which is the BG subtracted efficiency:
612     // h1 : TOTAL NUMBER OF EVENTS, SIGNAL REGION
613     // h2 : NUMBER OF EVENTS THAT PASS, SIGNAL REGION
614     // h3 : TOTAL NUMBER OF EVENTS, SIDE BAND
615     // h4 : NUMBER OF EVENTS THAT PASS, SIDE BAND
616    
617     #include "TH1.h"
618    
619    
620     TH1F* eff_bg(TH1F* h1, TH1F* h2, TH1F* h3, TH1F* h4, const char* name){
621    
622     // first, verify that all histograms have same binning
623     // nx is the number of visible bins
624     // nxtot = nx+2 includes underflow and overflow
625     Int_t nx = h1->GetNbinsX();
626     Int_t nxtot = nx + 2;
627     if (h2->GetNbinsX() != nx) {
628     cout << "Histograms must have same number of bins" << endl;
629     return 0;
630     }
631     if (h3->GetNbinsX() != nx) {
632     cout << "Histograms must have same number of bins" << endl;
633     return 0;
634     }
635     if (h3->GetNbinsX() != nx) {
636     cout << "Histograms must have same number of bins" << endl;
637     return 0;
638     }
639    
640     // get the new histogram
641     TH1F* temp = (TH1F*) h1->Clone(name);
642     temp->SetTitle(name);
643     temp->Reset();
644     temp->Sumw2();
645    
646     // Loop over bins, calculate efficiency and error, put it in histogram
647     for (Int_t i=0; i<nxtot; i++) {
648     Double_t x1 = h1->GetBinContent(i);
649     Double_t x2 = h2->GetBinContent(i);
650     Double_t x3 = h3->GetBinContent(i);
651     Double_t x4 = h4->GetBinContent(i);
652     Double_t denom = x1 - x3;
653     Double_t eff;
654     if (denom == 0.) {
655     eff = 0;
656     } else {
657     eff = (x2-x4)/denom;
658     }
659     Double_t failSig = x1 - x2;
660     Double_t failBg = x3 - x4;
661     Double_t blah = (1-eff)*(1-eff)*(x2+x4) + eff*eff*(failSig+failBg);
662     if (blah <= 0.) blah=0.0;
663     Double_t err;
664     if (denom == 0) {
665     err = 0.;
666     } else {
667     err = sqrt(blah)/denom;
668     }
669     temp->SetBinContent(i,eff);
670     temp->SetBinError(i,err);
671     }
672    
673     // Done
674     return temp;
675     }
676    
677     #include <TList.h>
678     #include <TIterator.h>
679    
680     void deleteHistos() {
681     // Delete all existing histograms in memory
682     TObject* obj;
683     TList* list = gDirectory->GetList() ;
684     TIterator* iter = list->MakeIterator();
685     while ((obj=iter->Next())) {
686     if (obj->IsA()->InheritsFrom(TH1::Class()) ||
687     obj->IsA()->InheritsFrom(TH2::Class()) ) {delete obj;}
688     }
689     }
690     void histio()
691     {
692     }
693    
694     void saveHist(const char* filename, const char* pat)
695     {
696     TList* list = gDirectory->GetList() ;
697     TIterator* iter = list->MakeIterator();
698    
699     TRegexp re(pat,kTRUE) ;
700    
701     TFile outf(filename,"RECREATE") ;
702     while(TObject* obj=iter->Next()) {
703     if (TString(obj->GetName()).Index(re)>=0) {
704     obj->Write() ;
705     cout << "." ;
706     cout.flush() ;
707     }
708     }
709     cout << endl ;
710     outf.Close() ;
711    
712     delete iter ;
713     }
714    
715    
716     void loadHist(const char* filename, const char* pfx, const char* pat, Bool_t doAdd)
717     {
718     TFile inf(filename) ;
719     //inf.ReadAll() ;
720     TList* list = inf.GetListOfKeys() ;
721     TIterator* iter = list->MakeIterator();
722    
723     TRegexp re(pat,kTRUE) ;
724     if (pat!=0) cout << "pat = " << pat << endl ;
725     else cout<<"no pattern: read all"<<endl;
726    
727     gDirectory->cd("Rint:") ;
728    
729     TObject* obj ;
730     TKey* key ;
731     cout << "doAdd = " << (doAdd?"T":"F") << endl ;
732     cout << "loadHist: reading." ;
733     while((key=(TKey*)iter->Next())) {
734    
735     // Int_t ridx = TString(key->GetName()).Index(re) ;
736     // cout<<"reading "<<key->GetName()<<std::endl;
737     if (pat!=0){
738     // cout<<"Check patt "<<pat<<endl;
739     if (TString(key->GetName()).Index(re)==-1) {
740     continue ;
741     }
742     }
743    
744     obj = inf.Get(key->GetName()) ;
745     TObject* clone ;
746     if (pfx) {
747    
748     // Find existing TH1-derived objects
749     TObject* oldObj = 0 ;
750     if (doAdd){
751     oldObj = gDirectory->Get(Form("%s_%s",pfx,obj->GetName())) ;
752     if (oldObj && !oldObj->IsA()->InheritsFrom(TH1::Class())) {
753     oldObj = 0 ;
754     }
755     }
756     if (oldObj) {
757     clone = oldObj ;
758     ((TH1*)clone)->Add((TH1*)obj) ;
759     } else {
760     clone = obj->Clone(Form("%s_%s",pfx,obj->GetName())) ;
761     }
762    
763    
764     } else {
765    
766     // Find existing TH1-derived objects
767     TObject* oldObj = 0 ;
768     if (doAdd){
769     oldObj = gDirectory->Get(key->GetName()) ;
770     if (oldObj && !oldObj->IsA()->InheritsFrom(TH1::Class())) {
771     oldObj = 0 ;
772     }
773     }
774    
775     if (oldObj) {
776     clone = oldObj ;
777     ((TH1*)clone)->Add((TH1*)obj) ;
778     } else {
779     clone = obj->Clone() ;
780     }
781     }
782     if (!gDirectory->GetList()->FindObject(clone)) {
783     gDirectory->Append(clone) ;
784     }
785     //cout << "." ;
786     cout.flush() ;
787     }
788     cout << endl;
789     inf.Close() ;
790     delete iter ;
791     }
792    
793    
794     //takes in a vector of prefixes, and a output prefix. Combines all histograms
795     //whose prefixes matches those in the input vector of TStrings into a new histogram
796     //with the the outPrefix
797     void combineHists(const std::vector<TString> v_prfxsToCombine, const TString outPrefix) {
798    
799    
800     if(v_prfxsToCombine.size() == 0) {
801     cout << "Input vector size is 0" << endl;
802     }
803     //histogram of common suffixes
804     vector<TString> v_suffixes;
805     //get the list of histograms that match the first entry in the prefix
806     TString temp = v_prfxsToCombine.at(0);
807     TList *rootlist = gDirectory->GetList();
808    
809     for(int i = 0; i < rootlist->GetSize(); i++) {
810    
811     TObject *obj = rootlist->At(i);
812     TString name = obj->GetName();
813     if(!obj->InheritsFrom(TH1::Class()))
814     continue;
815     if(!name.BeginsWith(temp+"_"))
816     continue;
817    
818     name.ReplaceAll(temp+"_", "_");
819    
820    
821     if(obj->InheritsFrom(TH2::Class())) {
822     TH2F *h = dynamic_cast<TH2F*>(obj->Clone());
823     h->SetName((outPrefix+name).Data());
824     h->SetTitle((outPrefix+name).Data());
825     for(unsigned int j = 1; j < v_prfxsToCombine.size(); j++) {
826     TH2F *htemp = dynamic_cast<TH2F*>(gDirectory->Get((v_prfxsToCombine.at(j) + name).Data()));
827     h->Add(htemp);
828     }
829     } else if(obj->InheritsFrom(TH1::Class())) {
830     TH1F *h = dynamic_cast<TH1F*>(obj->Clone());
831     h->SetName((outPrefix+name).Data());
832     h->SetTitle((outPrefix+name).Data());
833     for(unsigned int j = 1; j < v_prfxsToCombine.size(); j++) {
834     TH1F *htemp = dynamic_cast<TH1F*>(gDirectory->Get((v_prfxsToCombine.at(j) + name).Data()));
835     //cout << "TH1F: " << v_prfxsToCombine.at(j) + name << endl;
836     if(htemp == NULL)
837     cout << "htemp is NULL" << endl;
838     h->Add(htemp);
839     }
840     }
841     }//rootlist loop
842    
843    
844     //now delete the histos that match the prfxs to combine
845     for(unsigned int i = 0; i < v_prfxsToCombine.size(); i++ ) {
846    
847     // Delete all existing histograms in memory
848     TObject* obj;
849     TList* list = gDirectory->GetList() ;
850     TIterator* iter = list->MakeIterator();
851     while ((obj=iter->Next())) {
852     TString name = obj->GetName();
853     if(name.BeginsWith(outPrefix+"_"))
854     continue;
855     if(TString(obj->GetName()).BeginsWith(v_prfxsToCombine.at(i).Data())) {
856     if (obj->IsA()->InheritsFrom(TH1::Class()) ||
857     obj->IsA()->InheritsFrom(TH2::Class()) )
858     delete obj;
859     }
860     }
861     }//loop over prefixes
862    
863     }//fnc declaration
864    
865     }
866     #endif