FGolf/Tools/histtools.C

#include "TCanvas.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "THStack.h"
#include "TKey.h"
#include "TLegend.h"
#include "TRegexp.h"

#include <iostream>

using namespace std;

namespace hist {

     void add(const char* outHistName, const char* patORpfx);

   //Add all histograms whose names match one of ten possible regular expression
   //patterns or begin with one of ten possible given prefixes.  Feel free to
   //mix and match regular expression patterns and prefixes.  If the final hist-
   //ogram named outHistName does not exist it is created.

   void add(const char* outHistName, const char* patORpfx0, const char* patORpfx1, const char* patORpfx2 = 0, const char* patORpfx3 = 0, const char* patORpfx4 = 0, const char* patORpfx5 = 0, const char* patORpfx6 = 0, const char* patORpfx7 = 0, const char* patORpfx8 = 0, const char* patORpfx9 = 0)
   {
      add(outHistName, patORpfx0);
      add(outHistName, patORpfx1);
      if (patORpfx2) add(outHistName, patORpfx2);
      if (patORpfx3) add(outHistName, patORpfx3);
      if (patORpfx4) add(outHistName, patORpfx4);
      if (patORpfx5) add(outHistName, patORpfx5);
      if (patORpfx6) add(outHistName, patORpfx6);
      if (patORpfx7) add(outHistName, patORpfx7);
      if (patORpfx8) add(outHistName, patORpfx8);
      if (patORpfx9) add(outHistName, patORpfx9);
   }

   //Add all histograms whose names match the given regular expression pattern
   //or begin with the given prefix.  If the final histogram named outHistName
   //does not exist it is created.

   void add(const char* outHistName, const char* patORpfx) {
      TRegexp reg(patORpfx, kFALSE);

      TList* list = gDirectory->GetList() ;
      TIterator* iter = list->MakeIterator();

      TObject* obj = 0;
      TObject* hist = 0;
      Bool_t makeOutHist = false;

      hist = gDirectory->Get(outHistName);
      //If out hist does not exist, remember to create it
      if (! hist) makeOutHist = true;

      while ( (obj = iter->Next()) ) {
         if (! obj->InheritsFrom(TH1::Class())) continue;

         TString name = obj->GetName();
         //Don't add out hist
         if (name == TString(outHistName)) continue;

         if (TString(patORpfx).MaybeRegexp()) {
            if (TString(obj->GetName()).Index(reg) < 0 ) continue;
         } else if (! name.BeginsWith(patORpfx)) continue;

         if (makeOutHist) {
            hist = obj->Clone(outHistName);

            if (hist->InheritsFrom(TH2::Class()))
               ((TH2*)hist)->Reset();
            else
               ((TH1*)hist)->Reset();

            ((TH1*)hist)->SetTitle(outHistName);
            ((TH1*)hist)->Sumw2();
            makeOutHist = false;
         }

         ((TH1*)hist)->Add((TH1*)obj);
      }
   }

   //For all histograms whose names match the given regular expression pattern
   //or begin with the given prefix, set the fill, line and marker colors to the
   //given value.

   void color(const char* patORpfx, Color_t color) {
      TRegexp reg(patORpfx, kFALSE);

      TList* list = gDirectory->GetList() ;
      if (!list) {
        cout << "Failed to set color for " << patORpfx << endl;
        return;
      }
      TIterator* iter = list->MakeIterator();

      TObject* obj = 0;

      while ( (obj = iter->Next()) ) {
         if (! obj->InheritsFrom(TH1::Class())) continue;

         TString name = obj->GetName();

         if (TString(patORpfx).MaybeRegexp()) {
            if (TString(obj->GetName()).Index(reg) < 0 ) continue;
         } else if (! name.BeginsWith(patORpfx)) continue;

         ((TH1*)obj)->SetFillColor(color);
         ((TH1*)obj)->SetLineColor(color);
         ((TH1*)obj)->SetMarkerColor(color);
      }
   }

   //Return a pointer to a TLegend with an entry for each histogram drawn on a
   //given TCanvas.  Display either the line, point or fill values.  Optionally
   //apply colors to all histograms.  By default, entry labels are the names of
   //their respective histograms.  Optionally, if histogram names are of the
   //form XX_YY_ZZ_WW, entry labels can be XX (token=0), YY (token=1), etc.

     TLegend* legend(TCanvas* canvas, Option_t* option = "lpf", Bool_t addColor = kFALSE, Int_t token = -1,
                     Float_t xmin = 0.75, Float_t ymin = 0.75, Float_t xmax = 0.99, Float_t ymax = 0.99) {
      if(! canvas) return 0;

      TLegend* leg = new TLegend(xmin, ymin, xmax, ymax);
      TList* list = canvas->GetListOfPrimitives();
      TIterator* iter = list->MakeIterator();

      TObject* obj = 0;

      //Hist color iterator
      Int_t colorIt = 1;

      while ( (obj = iter->Next()) ) {
         if (! obj->InheritsFrom(TH1::Class())) continue;

         if (addColor) {
            hist::color(obj->GetName(), colorIt);
            ++colorIt;
         }

         if (token == -1)
            leg->AddEntry(obj, obj->GetName(), option);
         else {
            TString name(obj->GetName());
            TObjArray* a = name.Tokenize("_");
            if (a->GetEntries() <= token)
               leg->AddEntry(obj, obj->GetName(), option);
            else
               leg->AddEntry(obj, a->At(token)->GetName(), option);
         }
      }

      return leg;
   }

   //Return a pointer to a TLegend with an entry for each histogram added to a
   //given THStack.  Display either the line, point or fill values.  Optionally
   //apply colors to all histograms.  By default, entry labels are the names of
   //their respective histograms.  Optionally, if histogram names are of the
   //form XX_YY_ZZ_WW, entry labels can be XX (token=0), YY (token=1), etc.

   TLegend* legend(THStack* stack, Option_t* option = "lpf", Bool_t addColor = kFALSE, Int_t token = -1,
                   Float_t xmin = 0.75, Float_t ymin = 0.75, Float_t xmax = 0.99, Float_t ymax = 0.99) {
      if(! stack) return 0;

      TLegend* leg = new TLegend(xmin, ymin, xmax, ymax);
      TList* list = stack->GetHists();
      TIterator* iter = list->MakeIterator();

      TObject* obj = 0;

      //Hist color iterator
      Int_t colorIt = 1;

      while ( (obj = iter->Next()) ) {
         if (! obj->InheritsFrom(TH1::Class())) continue;

         if (addColor) {
            hist::color(obj->GetName(), colorIt);
            ++colorIt;
         }

         if (token == -1)
            leg->AddEntry(obj, obj->GetName(), option);
         else {
            TString name(obj->GetName());
            TObjArray* a = name.Tokenize("_");
            if (a->GetEntries() <= token)
               leg->AddEntry(obj, obj->GetName(), option);
            else
               leg->AddEntry(obj, a->At(token)->GetName(), option);
         }
      }

      return leg;
   }

   //Normalize to one all histograms whose names match the given regular exp-
   //ression pattern or begin with the given prefix.

   void normalize(const char* patORpfx) {
      TRegexp reg(patORpfx, kFALSE);

      TList* list = gDirectory->GetList() ;
      TIterator* iter = list->MakeIterator();

      TObject* obj = 0;

      while ( (obj = iter->Next()) ) {
         if (! obj->InheritsFrom(TH1::Class())) continue;

         TString name = obj->GetName();

         if (TString(patORpfx).MaybeRegexp()) {
            if (TString(obj->GetName()).Index(reg) < 0 ) continue;
         } else if (! name.BeginsWith(patORpfx)) continue;

         Double_t integral = 0;

         if (obj->InheritsFrom(TH2::Class()))
            integral = ((TH2*)obj)->Integral();
         else
            integral = ((TH1*)obj)->Integral();

         if (integral) {
            ((TH1*)obj)->Sumw2();
            ((TH1*)obj)->Scale(1./integral);
         }
      }
   }

   //Scale by the given value all histograms whose names match the given regular
   //expression pattern or begin with the given prefix.

   void scale(const char* patORpfx, Double_t scale) {
      TRegexp reg(patORpfx, kFALSE);

      TList* list = gDirectory->GetList() ;
      TIterator* iter = list->MakeIterator();

      TObject* obj = 0;

      while ( (obj = iter->Next()) ) {
         if (! obj->InheritsFrom(TH1::Class())) continue;

         TString name = obj->GetName();

         if (TString(patORpfx).MaybeRegexp()) {
            if (TString(obj->GetName()).Index(reg) < 0 ) continue;
         } else if (! name.BeginsWith(patORpfx)) continue;

         ((TH1*)obj)->Sumw2();
         ((TH1*)obj)->Scale(scale);
      }
   }

   //Don't you hate it when you draw multiple histograms on the same canvas only
   //to find that the bottom histogram's range does not encompass those of the
   //histograms drawn on top?  This method determines the maximum and minimum y
   //range of all the histograms drawn on a given TCanvas and appropriately re-
   //sizes the bottom histogram.

   void setrangey(TCanvas* canvas) {
      if(! canvas) return;

      TList* list = canvas->GetListOfPrimitives();
      TIterator* iter = list->MakeIterator();

      TObject* obj = 0;
      TObject* top = 0;

      //Extremes
      Double_t maxy = -999999;
      Double_t miny = 999999;

      while ( (obj = iter->Next()) ) {
         if (! obj->InheritsFrom(TH1::Class())) continue;

         if (! top) top = obj;

         if (((TH1*)obj)->GetMaximum() > maxy) maxy = ((TH1*)obj)->GetMaximum();
         if (((TH1*)obj)->GetMinimum() < miny) miny = ((TH1*)obj)->GetMinimum();
      }

      ((TH1*)top)->SetMaximum(maxy*1.3);
      //Protect against log scale
      if (canvas->GetLogy() && ! miny)
         ((TH1*)top)->SetMinimum(1E-4);
      else
         ((TH1*)top)->SetMinimum(miny*0.7);
   }

   //Create a stacked histogram consisting of all histograms whose names match
   //the given regular expression pattern or begin with the given prefix.  If
   //the THStack named stackHistName does not exist it is created.  Optionally
   //apply colors to all histograms.  Set drawOption to "nostack" if you do not
   //want to stack, to "hist" to display histograms without errors, to "histe"
   //to display histograms with errors, etc.

   void stack(const char* stackHistName, const char* patORpfx, Bool_t addColor = kFALSE, Option_t* drawOption = "") {
      TRegexp reg(patORpfx, kFALSE);

      TList* list = gDirectory->GetList() ;
      TIterator* iter = list->MakeIterator();

      TObject* obj = 0;
      TObject* stack = 0;
      Bool_t makeStackHist = false;

      stack = gDirectory->Get(stackHistName);
      //If stack hist does not exist, remember to create it
      if (! stack) makeStackHist = true;

      //Hist color iterator
      Int_t colorIt = 1;

      while ( (obj = iter->Next()) ) {
         if (! obj->InheritsFrom(TH1::Class())) continue;

         TString name = obj->GetName();

         if (TString(patORpfx).MaybeRegexp()) {
            if (TString(obj->GetName()).Index(reg) < 0 ) continue;
         } else if (! name.BeginsWith(patORpfx)) continue;

         if (makeStackHist) {
            stack = new THStack(stackHistName, stackHistName);
            makeStackHist = false;
         }

         if (addColor) {
            hist::color(obj->GetName(), colorIt);
            ++colorIt;
         }

         ((THStack*)stack)->Add((TH1*)obj, drawOption);
      }

      // Currently breaks .ls
      //gDirectory->Append(stack);
   }

   //Set the x-axis title of all histograms whose names match the given regular
   //expression pattern or begin with the given prefix.

   void xaxis(const char* patORpfx, const char* title) {
      TRegexp reg(patORpfx, kFALSE);

      TList* list = gDirectory->GetList() ;
      TIterator* iter = list->MakeIterator();

      TObject* obj = 0;

      while ( (obj = iter->Next()) ) {
         if (! (obj->InheritsFrom(TH1::Class()) || obj->InheritsFrom(THStack::Class()))) continue;

         TString name = obj->GetName();

         if (TString(patORpfx).MaybeRegexp()) {
            if (TString(obj->GetName()).Index(reg) < 0 ) continue;
         } else if (! name.BeginsWith(patORpfx)) continue;

         if (obj->InheritsFrom(TH1::Class()))
            ((TH1*)obj)->GetXaxis()->SetTitle(title);
         if (obj->InheritsFrom(THStack::Class())) {
            ((THStack*)obj)->Draw();
            ((THStack*)obj)->GetXaxis()->SetTitle(title);
         }
      }
   }

   //Set the y-axis title of all histograms whose names match the given regular
   //expression pattern or begin with the given prefix.

   void yaxis(const char* patORpfx, const char* title) {
      TRegexp reg(patORpfx, kFALSE);

      TList* list = gDirectory->GetList() ;
      TIterator* iter = list->MakeIterator();

      TObject* obj = 0;

      while ( (obj = iter->Next()) ) {
         if (! (obj->InheritsFrom(TH1::Class()) || obj->InheritsFrom(THStack::Class()))) continue;

         TString name = obj->GetName();

         if (TString(patORpfx).MaybeRegexp()) {
            if (TString(obj->GetName()).Index(reg) < 0 ) continue;
         } else if (! name.BeginsWith(patORpfx)) continue;

         if (obj->InheritsFrom(TH1::Class()))
            ((TH1*)obj)->GetYaxis()->SetTitle(title);
         if (obj->InheritsFrom(THStack::Class())) {
            ((THStack*)obj)->Draw();
            ((THStack*)obj)->GetYaxis()->SetTitle(title);
         }
      }
   }

}
// Input:  2 histogram
// Output: one histogram which is the efficiency:
// h1 :  TOTAL NUMBER OF EVENTS
// h2 :  NUMBER OF EVENTS THAT PASS

#include "TH1.h"

// Method by pointer
TH1F* eff(TH1F* h1, TH1F* h2, const char* name="eff"){

  // first, verify that all histograms have same binning
  // nx is the number of visible bins
  // nxtot = nx+2 includes underflow and overflow
  Int_t nx = h1->GetNbinsX();
  if (h2->GetNbinsX() != nx) {
    cout << "Histograms must have same number of bins" << endl;
    return 0;
  }

  // get the new histogram
  TH1F* temp = (TH1F*) h1->Clone(name);
  temp->SetTitle(name);
  temp->Reset();
  temp->Sumw2();

  // Do the calculation
  temp->Divide(h2,h1,1.,1.,"B");

  // Done
  return temp;
}


// Method by name
TH1F* eff(const char* name1, const char* name2, const char* name="eff"){

  // Get a list of object and their iterator
  TList* list = gDirectory->GetList() ;
  TIterator* iter = list->MakeIterator();

  // Loop over objects, set the pointers
  TObject* obj;
  TH1F* h1=0;
  TH1F* h2=0;
  TString str1 = Form("%s",name1);
  TString str2 = Form("%s",name2);
  while( (obj=iter->Next()) ) {
    TString objName = obj->GetName();
    if (objName == str1) h1 = (TH1F*) obj;
    if (objName == str2) h2 = (TH1F*) obj;
  }

  // quit if not found
  if (h1 == 0) {
    cout << "Histogram " << name1 << " not found" << endl;
    return 0;
  }
  if (h2 == 0) {
    cout << "Histogram " << name2 << " not found" << endl;
    return 0;
  }

  // Call the method by pointer
  TH1F* temp = eff(h1, h2, name);
  return temp;
}
// Input:  4 histogram
// Output: one histogram which is the BG subtracted efficiency:
// h1 :  TOTAL NUMBER OF EVENTS, SIGNAL REGION
// h2 :  NUMBER OF EVENTS THAT PASS, SIGNAL REGION
// h3 :  TOTAL NUMBER OF EVENTS, SIDE BAND
// h4 :  NUMBER OF EVENTS THAT PASS, SIDE BAND

#include "TH1.h"


TH1F* eff_bg(TH1F* h1, TH1F* h2, TH1F* h3, TH1F* h4, const char* name="eff"){

  // first, verify that all histograms have same binning
  // nx is the number of visible bins
  // nxtot = nx+2 includes underflow and overflow
  Int_t nx = h1->GetNbinsX();
  Int_t nxtot = nx + 2;
  if (h2->GetNbinsX() != nx) {
    cout << "Histograms must have same number of bins" << endl;
    return 0;
  }
  if (h3->GetNbinsX() != nx) {
    cout << "Histograms must have same number of bins" << endl;
    return 0;
  }
  if (h3->GetNbinsX() != nx) {
    cout << "Histograms must have same number of bins" << endl;
    return 0;
  }

  // get the new histogram
  TH1F* temp = (TH1F*) h1->Clone(name);
  temp->SetTitle(name);
  temp->Reset();
  temp->Sumw2();

  // Loop over bins, calculate efficiency and error, put it in histogram
  for (Int_t i=0; i<nxtot; i++) {
    Double_t x1 = h1->GetBinContent(i);
    Double_t x2 = h2->GetBinContent(i);
    Double_t x3 = h3->GetBinContent(i);
    Double_t x4 = h4->GetBinContent(i);
    Double_t denom = x1 - x3;
    Double_t eff;
    if (denom == 0.) {
      eff = 0;
    } else {
      eff   = (x2-x4)/denom;
    }
    Double_t failSig = x1 - x2;
    Double_t failBg  = x3 - x4;
    Double_t blah    = (1-eff)*(1-eff)*(x2+x4) + eff*eff*(failSig+failBg);
    if (blah <= 0.) blah=0.0;
    Double_t err;
    if (denom == 0) {
      err = 0.;
    } else {
      err = sqrt(blah)/denom;
    }
    temp->SetBinContent(i,eff);
    temp->SetBinError(i,err);
  }

  // Done
  return temp;
}

#include <TList.h>
#include <TIterator.h>

void deleteHistos() {
   // Delete all existing histograms in memory
   TObject* obj;
   TList* list = gDirectory->GetList() ;
   TIterator* iter = list->MakeIterator();
   while ( (obj=iter->Next()) ) {
     if (obj->IsA()->InheritsFrom(TH1::Class()) ||
         obj->IsA()->InheritsFrom(TH2::Class()) ) {delete obj;}
   }
}

void histio()
{
}

void saveHist(const char* filename, const char* pat="*")
{
   TList* list = gDirectory->GetList() ;
   TIterator* iter = list->MakeIterator();

   TRegexp re(pat,kTRUE) ;

   TFile outf(filename,"RECREATE") ;
   while(TObject *obj=iter->Next()) {
      if (TString(obj->GetName()).Index(re)>=0) {
         obj->Write() ;
         cout << "." ;
         cout.flush() ;
      }
   }
   cout << endl ;
   outf.Close() ;

   delete iter ;
}


void loadHist(const char* filename, const char* pfx=0, const char* pat="*", Bool_t doAdd=kFALSE)
{
   TFile inf(filename) ;
   //inf.ReadAll() ;
   TList* list = inf.GetListOfKeys() ;
   TIterator* iter = list->MakeIterator();

   TRegexp re(pat,kTRUE) ;
   cout << "pat = " << pat << endl ;

   gDirectory->cd("Rint:") ;

   TObject* obj ;
   TKey* key ;
   cout << "doAdd = " << (doAdd?"T":"F") << endl ;
   cout << "loadHist: reading." ;
   while( (key=(TKey*)iter->Next()) ) {

      Int_t ridx = TString(key->GetName()).Index(re) ;
      if (ridx==-1) {
         continue ;
      }

      obj = inf.Get(key->GetName()) ;
      TObject* clone ;
      if (pfx) {

         // Find existing TH1-derived objects
         TObject* oldObj = 0 ;
         if (doAdd){
            oldObj = gDirectory->Get(Form("%s_%s",pfx,obj->GetName())) ;
            if (oldObj && !oldObj->IsA()->InheritsFrom(TH1::Class())) {
               oldObj = 0 ;
            }
         }
         if (oldObj) {
            clone = oldObj ;
            ((TH1*)clone)->Add((TH1*)obj) ;
         } else {
            clone = obj->Clone(Form("%s_%s",pfx,obj->GetName())) ;
         }


      } else {

         // Find existing TH1-derived objects
         TObject* oldObj = 0 ;
         if (doAdd){
            oldObj = gDirectory->Get(key->GetName()) ;
            if (oldObj && !oldObj->IsA()->InheritsFrom(TH1::Class())) {
               oldObj = 0 ;
            }
         }

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