Plot/src/PlotTask.cc

// $Id: PlotTask.cc,v 1.2 2011/01/25 11:30:30 paus Exp $

#include <vector>
#include <TROOT.h>
#include <TSystem.h>
#include <TMath.h>
#include <TPad.h>
#include <TFile.h>
#include <TH1.h>
#include <TH1D.h>
#include <TBox.h>
#include <TMarker.h>
#include <TLatex.h>
#include "MitAna/DataUtil/interface/Debug.h"
#include "MitPlots/Style/interface/MitStyle.h"
#include "MitPlots/Plot/interface/PlotTask.h"

ClassImp(mithep::PlotTask)

using namespace std;
using namespace mithep;

//--------------------------------------------------------------------------------------------------
PlotTask::PlotTask(const TaskSamples *taskSamples, const double lumi) :
  fTask        (taskSamples),
  fHistStyles  (0),
  fTargetLumi  (lumi),
  fIdxHistMax  (-1),
  fNRebin      (1),
  fEmptyHist   (0),
  fDataHist    (0),
  fHistMinimum (0),
  fHistMaximum (0),
  fHistXMinimum(0),
  fHistXMaximum(0),
  fAxisTitleX  (""),
  fAxisTitleY  ("Number of Events"),
  fXLegend     (70.),
  fYLegend     (98.)
{
  // Constructor
}

//--------------------------------------------------------------------------------------------------
PlotTask::~PlotTask()
{
  // Destructor

  if (fEmptyHist)
    delete fEmptyHist;
  if (fDataHist)
    delete fDataHist;
}

//--------------------------------------------------------------------------------------------------
void PlotTask::SetAxisTitles(const char* xtit, const char* ytit)
{
  fAxisTitleX = TString(xtit);
  fAxisTitleY = TString(ytit);

  return;
}

//--------------------------------------------------------------------------------------------------
void PlotTask::PlotContributions(const char* dir, const char* hist)
{
  // Show present list of defined samples

  // scale the histograms
  ScaleHistograms(dir,hist);
  FindHistMaximum();

  // say what we are doing
  printf("\n ==== Plotting Contributions -- %s ====\n\n",fTask->Name()->Data());
  printf("  index of highest histogram %d (0-%d)\n\n",fIdxHistMax,fHists.size()-1);

  // check x-axis title
  if (fAxisTitleX == TString(""))
    fAxisTitleX = TString(hist);

  // initialize the largest histogram properly and draw it
  TH1D* hmax = 0;
  if (fIdxHistMax == fHists.size())
    hmax = fDataHist;
  else
    hmax = fHists[fHists.size()-1];
  MitStyle::InitHist(hmax,fAxisTitleX.Data(),fAxisTitleY.Data(),kBlack);
  hmax->GetXaxis()->SetNdivisions(505);
  if (fHistMinimum != 0)
    hmax->SetMinimum(fHistMinimum);
  if (fIdxHistMax == fHists.size())
    hmax->Draw("E");
  else
    hmax->Draw("hist");

  // loop through samples and draw all histograms
  int iCol  = (EColor) kBlack;
  int iFill = 3001;
  for (UInt_t i=0; i<*fTask->NSamples(); i++) {
    //const Sample *s = fTask->GetSample(i);
    MitStyle::InitHist(fHists[i],hist,"Number of Events",(EColor) iCol);

    fHists[i]->SetLineColor(iCol);
    //fHists[i]->SetFillColor(iCol);
    //fHists[i]->SetFillStyle(iFill);
    fHists[i]->Draw("same;Hist");

    iFill++; iCol++;
  }

  // overlay a frame to put the text and boxes on
  OverlayFrame();

  return;
}

//--------------------------------------------------------------------------------------------------
void PlotTask::PlotStack(const char* dir, const char* hist)
{
  // Show present list of defined samples

  // scale the histograms
  ScaleHistograms(dir,hist);
  FindHistMaximum();

  // ensure the histogram styles are ready
  if (! fHistStyles)
    fHistStyles = new HistStyles();

  // check basics
  if (fStackedHists.size() < 1) {
    printf(" WARNING -- no histograms with name: %s. EXIT!\n",hist);
    return;
  }

  // determine the width of one bin
  double binW = fEmptyHist->GetBinWidth(1);
  char binWStr[7];
  sprintf(binWStr," /%7.3f",binW);

  // say what we are doing
  printf("\n ==== Plotting Stack -- %s ====\n",fTask->Name()->Data());
  printf("\n   bin width: %7.3f  in  xMin: %7.3f  <--> xMax: %7.3f\n\n",
         binW,fHistXMinimum,fHistXMaximum);

  // check x-axis title
  if (fAxisTitleX == TString(""))
    fAxisTitleX = TString(hist);

  // create the requested frame and initialize the empty histogram properly and draw it
  TH1D *hFrame = 0;
  if (fHistXMinimum != 0. || fHistXMaximum != 0.) {
    hFrame = new TH1D("Frame",fEmptyHist->GetTitle(),1,fHistXMinimum,fHistXMaximum);
  }
  else
    hFrame = fEmptyHist;

  MitStyle::InitHist(hFrame,fAxisTitleX.Data(),(fAxisTitleY+TString(binWStr)).Data(),kBlack);
  hFrame->GetXaxis()->SetNdivisions(505);
  if (fHistMinimum != 0)
    hFrame->SetMinimum(fHistMinimum);
  hFrame->SetMaximum(fHistMaximum*1.1);
  hFrame->Draw("");

  // loop through samples and draw all histograms
  fHistStyles->ResetStyle();
  for (UInt_t i=fStackedHists.size(); i>0; i--) {
    const Sample *s = fTask->GetSample(i);
    MitStyle::InitHist(fStackedHists[i-1],fAxisTitleX.Data(),fAxisTitleY.Data());
    if (i == fStackedHists.size()) {
      fHistStyles->ApplyCurrentStyle(fStackedHists[i-1]);
      fStackedHists[i-1]->Draw("same;Hist");
      fHistStyles->NextStyle();
    }
    else {
      if (*s->Legend() != TString(" ")) {
        // white out the histogram first
        fStackedHists[i-1]->SetFillColor(10);
        fStackedHists[i-1]->SetFillStyle(1001);
        fStackedHists[i-1]->DrawCopy("same;Hist");
        // overlay the histogram with proper hatching
        fHistStyles->ApplyCurrentStyle(fStackedHists[i-1]);
        fStackedHists[i-1]->Draw("same;Hist");
        fHistStyles->NextStyle();
      }
    }
  }

  if (fDataHist) {
    fHistStyles->ApplyDataStyle(fDataHist);
    fDataHist->Draw("same;E");
  }

  // overlay a frame to put the text and boxes on
  OverlayFrame();

  return;
}

//--------------------------------------------------------------------------------------------------
void PlotTask::ScaleHistograms(const char* dir, const char* hist)
{
  // Scale the histograms according to the cross section and the desired lumi and store them
  // for later use

  if (fHists.size() > 0) {
    printf(" WARNING - scaled histograms already exist. EXIT.");
    return;
  }

  fIdxHistMax = -1;

  // some useful definitions
  TString *dirFwk  = new TString("AnaFwkMod");
  TString *allEvts = new TString("hDAllEvents");
  TString  slash  ("/");

  // say what we are doing
  printf("\n ==== Extracting and Scaling Contributions -- %s ====\n\n",fTask->Name()->Data());
  printf("  target luminosity %f 1/fb\n\n",fTargetLumi/1000.);

  //------------------------------------------------------------------------------------------------
  // go through the Monte Carlo samples
  //------------------------------------------------------------------------------------------------
  // loop through samples and determine maximum
  printf("\n Monte Carlo \n");
  double nTotRaw = 0.0, nTot = 0.0, nTot2 = 0.0;
  for (UInt_t i=0; i<*fTask->NSamples(); i++) {
    const Sample *s = fTask->GetSample(i);
    // open file belonging to this sample
    TFile *fif = new TFile((*fTask->Dir()+slash+*s->File()).Data());
    // make sure the file exists
    if (fif->IsOpen() == kFALSE) {
      printf(" WARNING -- sample  %s  does not have a histogram file. Continue without!\n",
             s->Name()->Data());
      continue;
    }
    // read and determine general properties of this sample
    TDirectory *dirTmp = (TDirectory*) gROOT->FindObject(dirFwk->Data());
    if (dirTmp)
      fif->cd(dirFwk->Data());
    TH1D *hAllEvts = (TH1D*) gROOT->FindObject(allEvts->Data());
    if (! hAllEvts) {
      printf(" WARNING -- sample  %s  does not have a framework file. Next sample!\n",
             s->Name()->Data());
      continue;
    }
    double nEvts   = hAllEvts->GetEntries();
    double lumi    = nEvts / *s->Xsec();
    double factor,scale = *s->Scale();
    if      (lumi < 0)
      factor = 1.0;
    else if (lumi > 0)
      factor = fTargetLumi/lumi;
    else
      factor = 0;

    TH1D *h    = (TH1D*) fif->FindObjectAny(hist);
    if (! h) {
      printf(" WARNING -- sample  %s  does not have requested histogram. Next sample!\n",
             s->Name()->Data());
      continue;
    }
    if (! fEmptyHist) {
      fEmptyHist = new TH1D(factor * scale * (*h));
      fEmptyHist->Rebin(fNRebin);
      fEmptyHist->Reset();
    }

    double nEvtsSelRaw = h->GetEntries();
    double nEvtsSel    = nEvtsSelRaw * factor * scale;
    double nEvtsSelErr = TMath::Sqrt(nEvtsSelRaw) * factor * scale;

    printf(" -> %-40s - %14.0f %12.3f +- %8.3f %16.7f: %16.4f (x %f x %f - %p)\n",
           s->Name()->Data(),nEvts,nEvtsSel,nEvtsSelErr,*s->Xsec(),lumi,factor,scale,(void*)h);

    nTotRaw += nEvts;
    nTot    += nEvtsSel;
    nTot2   += nEvtsSelErr*nEvtsSelErr;

    // scale it
    TH1D *hTmp = new TH1D(factor * scale * (*h));
    hTmp->Rebin(fNRebin);
    // put it into our collection
    fHists.push_back(hTmp);
    // make a new histogram from the last added one and add the recent
    if (i == 0)
      fStackedHists.push_back(hTmp);
    else {
      TH1D *hStackedTmp = new TH1D(*fStackedHists[fStackedHists.size()-1]);
      // remember the stored histogram is already rebinned
      hStackedTmp->Add(hTmp);
      fStackedHists.push_back(hStackedTmp);
    }
  }
  // Monte Carlo summary
  printf(" %-40s    - %14.0f %12.3f +- %8.3f %16.7f: %16.4f (x %f x %f)\n",
         "== Monte Carlo Total ==",nTotRaw,nTot,TMath::Sqrt(nTot2),0.0,0.0,1.0,1.0);


  //------------------------------------------------------------------------------------------------
  // go through the data samples
  //------------------------------------------------------------------------------------------------
  // loop through data samples and add them up, straight away
  if (*fTask->NDataSamples() > 0)
    printf("\n Data \n");
  nTotRaw = 0.0;
  nTot    = 0.0;
  nTot2   = 0.0;
  for (UInt_t i=0; i<*fTask->NDataSamples(); i++) {
    const Sample *s = fTask->GetDataSample(i);
    // open file belonging to the data sample
    TFile *fif = new TFile((*fTask->Dir()+slash+*s->File()).Data());
    // make sure the file exists
    if (fif->IsOpen() == kFALSE) {
      printf(" WARNING -- sample  %s  does not have a histogram file. Continue without!\n",
             s->Name()->Data());
      continue;
    }
    else {
      // read and determine general properties of this sample
      TDirectory *dirTmp = (TDirectory*) gROOT->FindObject(dirFwk->Data());
      if (dirTmp)
        fif->cd(dirFwk->Data());
      TH1D *hAllEvts = (TH1D*) gROOT->FindObject(allEvts->Data());
      if (! hAllEvts) {
        printf(" WARNING -- sample  %s  does not have a framework file. Next sample!\n",
               s->Name()->Data());
        continue;
      }
      else {
        double nEvts = hAllEvts->GetEntries();

        TH1D *h = (TH1D*) fif->FindObjectAny(hist);
        if (! h) {
          printf(" WARNING -- sample  %s  does not have requested histogram. Next sample!\n",
                 s->Name()->Data());
          continue;
        }
        else
          // rebin it
          h->Rebin(fNRebin);

        double nEvtsSelRaw = h->GetEntries();
        double nEvtsSel    = nEvtsSelRaw;
        double nEvtsSelErr = TMath::Sqrt(nEvtsSelRaw);

        printf(" -> %-40s - %14.0f %12.3f +- %8.3f %16.7f: %16.4f (x %f)\n",
               s->Name()->Data(),nEvts,nEvtsSel,nEvtsSelErr,*s->Xsec(),fTargetLumi,1.0);

        nTotRaw += nEvts;
        nTot    += nEvtsSel;
        nTot2   += nEvtsSelErr*nEvtsSelErr;

        // construct the complete data histogram
        if (! fDataHist)
          fDataHist = new TH1D(1.0 * (*h));
        else
          fDataHist->Add(h);
      }
    }
  }
  // Data summary
  printf(" %-40s    - %14.0f %12.3f +- %8.3f %16.7f: %16.4f (x %f)\n\n",
         "== Data Total =========",nTotRaw,nTot,TMath::Sqrt(nTot2),0.0,0.0,1.0);

  return;
}

//--------------------------------------------------------------------------------------------------
void PlotTask::FindHistMaximum()
{
  // Find maximum of all histograms

  // first check whether value was overwritten by hand
  if (fHistMaximum>0.)
    return;
  
  // search through the single histograms
  for (UInt_t i=0; i<fHists.size(); i++) {
    if (fHists[i]->GetMaximum() > fHistMaximum) {
      fHistMaximum = fHists[i]->GetMaximum();
      fIdxHistMax  = i;
    }
  }

  // search through the stacked histograms
  fHistMaximum = 0.;
  for (UInt_t i=0; i<fStackedHists.size(); i++) {
    if (fStackedHists[i]->GetMaximum() > fHistMaximum)
      fHistMaximum = fStackedHists[i]->GetMaximum();
  }
  
  // check the data histogram if present
  if (fDataHist && fDataHist->GetMaximum() > fHistMaximum) {
    fHistMaximum = fDataHist->GetMaximum();
    fIdxHistMax  = fHists.size();
  }

  printf(" Histogram maximum is set to be: %f\n",fHistMaximum);

  return;
}

//--------------------------------------------------------------------------------------------------
void PlotTask::OverlayEmptyHist() const
{
  // Overlay an empty histogram onto the picture
  if (fEmptyHist)
    fEmptyHist->Draw("same");

  return;
}

//--------------------------------------------------------------------------------------------------
void PlotTask::OverlayFrame() const
{
  // Overlay a linear frame from user coordinates (0-100,0-100)

  // create new transparent pad for the text
  TPad *transPad = new TPad("transPad","Transparent Pad",0,0,1,1);
  transPad->SetFillStyle(4000);
  transPad->Draw();
  transPad->cd();
  // find out the right normalization to define the new range (histogram: 0,0 -> 100,100)
  double xtot = 1./(1. - transPad->GetLeftMargin() - transPad->GetRightMargin());
  double ytot = 1./(1. - transPad->GetBottomMargin() - transPad->GetTopMargin());
  transPad->Range((xtot*transPad->GetLeftMargin()  *-100      ),
                  (ytot*transPad->GetBottomMargin()*-100      ),
                  (xtot*transPad->GetRightMargin() * 100 + 100),
                  (ytot*transPad->GetTopMargin()   * 100 + 100));
  MDB(kGeneral,1)
    printf(" Range: %f %f %f %f\n",
           (xtot*transPad->GetLeftMargin()  *-100      ),
           (ytot*transPad->GetBottomMargin()*-100      ),
           (xtot*transPad->GetRightMargin() * 100 + 100),
           (ytot*transPad->GetTopMargin()   * 100 + 100) );

  // define the basic text to be used
  TLatex* text = new TLatex();
  text->SetTextFont (42);
  text->SetTextAlign(12);
  text->SetTextSize (0.03);

  // define the basic box to be used
  TBox *box = new TBox();
  box->SetLineWidth(2);
  // base coordinates for the text/boxes
  double xCorner = fXLegend, yCorner = fYLegend, yDelLine = 4.;
  double xIndent = 1.5*yDelLine;
  // count samples with non empty legend
  int nLegends = 0;
  for (UInt_t i=*fTask->NSamples(); i>0; i--) {
    // attach to the specific sample
    const Sample *s = fTask->GetSample(i-1);
    if (*s->Legend() != TString(" "))
      nLegends++;
  }
  int iDat = 0;
  if (fDataHist) {
    iDat = 1;
    nLegends++;
  }

  // set start values
  fHistStyles->ResetStyle();
  int iLeg  = 0;
  // loop through the sampels
  for (UInt_t i=*fTask->NSamples(); i>0; i--) {
    // attach to the specific sample
    const Sample *s = fTask->GetSample(i-1);
    // calculate corners for the text
    double xText = xCorner+xIndent, yText = yCorner-float((iLeg+iDat)+0.5)*yDelLine;
    // say what goes where
    MDB(kGeneral,1)
      printf(" Adding box at: (x1,y1,x2,y2) = (%6.2f,%6.2f,%6.2f,%6.2f)\n",
             xCorner+0.15*xIndent,yText-0.3*yDelLine,
             xCorner+0.85*xIndent,yText+0.4*yDelLine);
    // plot the box
    if (*s->Legend() != TString(" ")) {
      const HistStyle *hStyle = fHistStyles->CurrentStyle();
      box->SetFillStyle(0);
      box->SetFillColor(hStyle->Color());
      box->SetLineColor(hStyle->Color());
      box->DrawBox(xCorner+0.15*xIndent,yText-0.3*yDelLine,
                   xCorner+0.85*xIndent,yText+0.4*yDelLine);
      box->SetFillStyle(hStyle->FillStyle());
      box->DrawBox(xCorner+0.15*xIndent,yText-0.3*yDelLine,
                   xCorner+0.85*xIndent,yText+0.4*yDelLine);

      TString pText = *s->Legend();
      MDB(kGeneral,1)
        printf(" Adding text \"%s\" at: (x,y) = (%6.2f,%6.2f)\n",pText.Data(),xText,yText);
      // set the proper text color
      text->SetTextColor(hStyle->Color());
      // plot the text
      text->SetTextAlign(12);
      text->DrawLatex(xText,yText,pText.Data());

      fHistStyles->NextStyle();
      // keep track of how many were printed
      iLeg++;
    }
  }

  // attach to the data sample
  if (fDataHist) {
    const Sample *s = fTask->GetDataSample(0);
    // calculate corners for the text
    double xText = xCorner+xIndent, yText = yCorner-float(1-0.7)*yDelLine;
    // say what goes where
    MDB(kGeneral,1)
      printf(" Adding marker at: (x,y) = (%6.2f,%6.2f)\n",
             xCorner+0.5*xIndent,yText);
    // plot the marker
    if (*s->Legend() != TString(" ")) {
      const HistStyle *hStyle = fHistStyles->GetDataStyle();

      // draw marker
      TMarker *m = new TMarker(xCorner+0.5*xIndent,yText,20);
      m->SetMarkerColor(hStyle->Color      ());
      m->SetMarkerSize (hStyle->MarkerSize ());
      m->SetMarkerStyle(hStyle->MarkerStyle());
      m->Draw();
      // draw text
      char l[1024];
      sprintf(l,"%4.2f",fTargetLumi);
      //TString pText = *s->Legend() + TString(" (#int L = ") + TString(l) + TString("/pb)");
      TString pText = *s->Legend() + TString(" (") + TString(l) + TString("/pb)");
      MDB(kGeneral,1)
        printf(" Adding text \"%s\" at: (x,y) = (%6.2f,%6.2f)\n",pText.Data(),xText,yText);
      // set the proper text color
      text->SetTextColor(hStyle->Color());
      // plot the text
      text->SetTextAlign(12);
      text->DrawLatex(xText,yText,pText.Data());
    }
  }

  // Make sure the histogram frame is nice
  box->SetFillStyle(0);
  box->SetLineColor(kBlack);
  box->DrawBox(0,0,100,100);

  delete text;
  delete box;

  return;
}
Revision:	1.3
Committed:	Thu Jan 27 14:10:01 2011 UTC (14 years, 3 months ago) by paus
Content type:	text/plain
Branch:	MAIN
Changes since 1.2:	+23 -6 lines
Log Message:	Tuning.
#	User	Rev	Content
1	paus	1.3	// $Id: PlotTask.cc,v 1.2 2011/01/25 11:30:30 paus Exp $
2	paus	1.2
3			#include <vector>
4			#include <TROOT.h>
5			#include <TSystem.h>
6			#include <TMath.h>
7			#include <TPad.h>
8			#include <TFile.h>
9			#include <TH1.h>
10			#include <TH1D.h>
11			#include <TBox.h>
12			#include <TMarker.h>
13			#include <TLatex.h>
14			#include "MitAna/DataUtil/interface/Debug.h"
15			#include "MitPlots/Style/interface/MitStyle.h"
16			#include "MitPlots/Plot/interface/PlotTask.h"
17
18			ClassImp(mithep::PlotTask)
19
20			using namespace std;
21			using namespace mithep;
22
23			//--------------------------------------------------------------------------------------------------
24			PlotTask::PlotTask(const TaskSamples *taskSamples, const double lumi) :
25			fTask (taskSamples),
26			fHistStyles (0),
27			fTargetLumi (lumi),
28			fIdxHistMax (-1),
29			fNRebin (1),
30			fEmptyHist (0),
31			fDataHist (0),
32			fHistMinimum (0),
33			fHistMaximum (0),
34			fHistXMinimum(0),
35			fHistXMaximum(0),
36			fAxisTitleX (""),
37			fAxisTitleY ("Number of Events"),
38			fXLegend (70.),
39			fYLegend (98.)
40			{
41			// Constructor
42			}
43
44			//--------------------------------------------------------------------------------------------------
45			PlotTask::~PlotTask()
46			{
47			// Destructor
48
49			if (fEmptyHist)
50			delete fEmptyHist;
51			if (fDataHist)
52			delete fDataHist;
53			}
54
55			//--------------------------------------------------------------------------------------------------
56			void PlotTask::SetAxisTitles(const char* xtit, const char* ytit)
57			{
58			fAxisTitleX = TString(xtit);
59			fAxisTitleY = TString(ytit);
60
61			return;
62			}
63
64			//--------------------------------------------------------------------------------------------------
65			void PlotTask::PlotContributions(const char* dir, const char* hist)
66			{
67			// Show present list of defined samples
68
69			// scale the histograms
70			ScaleHistograms(dir,hist);
71			FindHistMaximum();
72
73			// say what we are doing
74			printf("\n ==== Plotting Contributions -- %s ====\n\n",fTask->Name()->Data());
75			printf(" index of highest histogram %d (0-%d)\n\n",fIdxHistMax,fHists.size()-1);
76
77			// check x-axis title
78			if (fAxisTitleX == TString(""))
79			fAxisTitleX = TString(hist);
80
81			// initialize the largest histogram properly and draw it
82			TH1D* hmax = 0;
83			if (fIdxHistMax == fHists.size())
84			hmax = fDataHist;
85			else
86			hmax = fHists[fHists.size()-1];
87			MitStyle::InitHist(hmax,fAxisTitleX.Data(),fAxisTitleY.Data(),kBlack);
88			hmax->GetXaxis()->SetNdivisions(505);
89			if (fHistMinimum != 0)
90			hmax->SetMinimum(fHistMinimum);
91			if (fIdxHistMax == fHists.size())
92			hmax->Draw("E");
93			else
94			hmax->Draw("hist");
95
96			// loop through samples and draw all histograms
97			int iCol = (EColor) kBlack;
98			int iFill = 3001;
99			for (UInt_t i=0; i<*fTask->NSamples(); i++) {
100			//const Sample *s = fTask->GetSample(i);
101			MitStyle::InitHist(fHists[i],hist,"Number of Events",(EColor) iCol);
102
103			fHists[i]->SetLineColor(iCol);
104			//fHists[i]->SetFillColor(iCol);
105			//fHists[i]->SetFillStyle(iFill);
106			fHists[i]->Draw("same;Hist");
107
108			iFill++; iCol++;
109			}
110
111			// overlay a frame to put the text and boxes on
112			OverlayFrame();
113
114			return;
115			}
116
117			//--------------------------------------------------------------------------------------------------
118			void PlotTask::PlotStack(const char* dir, const char* hist)
119			{
120			// Show present list of defined samples
121
122			// scale the histograms
123			ScaleHistograms(dir,hist);
124			FindHistMaximum();
125
126			// ensure the histogram styles are ready
127			if (! fHistStyles)
128			fHistStyles = new HistStyles();
129
130			// check basics
131			if (fStackedHists.size() < 1) {
132			printf(" WARNING -- no histograms with name: %s. EXIT!\n",hist);
133			return;
134			}
135
136			// determine the width of one bin
137			double binW = fEmptyHist->GetBinWidth(1);
138			char binWStr[7];
139			sprintf(binWStr," /%7.3f",binW);
140
141			// say what we are doing
142			printf("\n ==== Plotting Stack -- %s ====\n",fTask->Name()->Data());
143			printf("\n bin width: %7.3f in xMin: %7.3f <--> xMax: %7.3f\n\n",
144			binW,fHistXMinimum,fHistXMaximum);
145
146			// check x-axis title
147			if (fAxisTitleX == TString(""))
148			fAxisTitleX = TString(hist);
149
150			// create the requested frame and initialize the empty histogram properly and draw it
151			TH1D *hFrame = 0;
152			if (fHistXMinimum != 0. \|\| fHistXMaximum != 0.) {
153			hFrame = new TH1D("Frame",fEmptyHist->GetTitle(),1,fHistXMinimum,fHistXMaximum);
154			}
155			else
156			hFrame = fEmptyHist;
157
158			MitStyle::InitHist(hFrame,fAxisTitleX.Data(),(fAxisTitleY+TString(binWStr)).Data(),kBlack);
159			hFrame->GetXaxis()->SetNdivisions(505);
160			if (fHistMinimum != 0)
161			hFrame->SetMinimum(fHistMinimum);
162			hFrame->SetMaximum(fHistMaximum*1.1);
163			hFrame->Draw("");
164
165			// loop through samples and draw all histograms
166			fHistStyles->ResetStyle();
167			for (UInt_t i=fStackedHists.size(); i>0; i--) {
168			const Sample *s = fTask->GetSample(i);
169			MitStyle::InitHist(fStackedHists[i-1],fAxisTitleX.Data(),fAxisTitleY.Data());
170			if (i == fStackedHists.size()) {
171			fHistStyles->ApplyCurrentStyle(fStackedHists[i-1]);
172			fStackedHists[i-1]->Draw("same;Hist");
173			fHistStyles->NextStyle();
174			}
175			else {
176			if (*s->Legend() != TString(" ")) {
177			// white out the histogram first
178			fStackedHists[i-1]->SetFillColor(10);
179			fStackedHists[i-1]->SetFillStyle(1001);
180			fStackedHists[i-1]->DrawCopy("same;Hist");
181			// overlay the histogram with proper hatching
182			fHistStyles->ApplyCurrentStyle(fStackedHists[i-1]);
183			fStackedHists[i-1]->Draw("same;Hist");
184			fHistStyles->NextStyle();
185			}
186			}
187			}
188
189			if (fDataHist) {
190			fHistStyles->ApplyDataStyle(fDataHist);
191			fDataHist->Draw("same;E");
192			}
193
194			// overlay a frame to put the text and boxes on
195			OverlayFrame();
196
197			return;
198			}
199
200			//--------------------------------------------------------------------------------------------------
201			void PlotTask::ScaleHistograms(const char* dir, const char* hist)
202			{
203			// Scale the histograms according to the cross section and the desired lumi and store them
204			// for later use
205
206			if (fHists.size() > 0) {
207			printf(" WARNING - scaled histograms already exist. EXIT.");
208			return;
209			}
210
211			fIdxHistMax = -1;
212
213			// some useful definitions
214			TString *dirFwk = new TString("AnaFwkMod");
215			TString *allEvts = new TString("hDAllEvents");
216			TString slash ("/");
217
218			// say what we are doing
219			printf("\n ==== Extracting and Scaling Contributions -- %s ====\n\n",fTask->Name()->Data());
220			printf(" target luminosity %f 1/fb\n\n",fTargetLumi/1000.);
221
222			//------------------------------------------------------------------------------------------------
223			// go through the Monte Carlo samples
224			//------------------------------------------------------------------------------------------------
225			// loop through samples and determine maximum
226	paus	1.3	printf("\n Monte Carlo \n");
227			double nTotRaw = 0.0, nTot = 0.0, nTot2 = 0.0;
228	paus	1.2	for (UInt_t i=0; i<*fTask->NSamples(); i++) {
229			const Sample *s = fTask->GetSample(i);
230			// open file belonging to this sample
231			TFile fif = new TFile((fTask->Dir()+slash+*s->File()).Data());
232			// make sure the file exists
233			if (fif->IsOpen() == kFALSE) {
234			printf(" WARNING -- sample %s does not have a histogram file. Continue without!\n",
235			s->Name()->Data());
236			continue;
237			}
238			// read and determine general properties of this sample
239			TDirectory dirTmp = (TDirectory) gROOT->FindObject(dirFwk->Data());
240			if (dirTmp)
241			fif->cd(dirFwk->Data());
242			TH1D hAllEvts = (TH1D) gROOT->FindObject(allEvts->Data());
243			if (! hAllEvts) {
244			printf(" WARNING -- sample %s does not have a framework file. Next sample!\n",
245			s->Name()->Data());
246			continue;
247			}
248			double nEvts = hAllEvts->GetEntries();
249			double lumi = nEvts / *s->Xsec();
250			double factor,scale = *s->Scale();
251			if (lumi < 0)
252			factor = 1.0;
253			else if (lumi > 0)
254			factor = fTargetLumi/lumi;
255			else
256			factor = 0;
257
258			TH1D h = (TH1D) fif->FindObjectAny(hist);
259			if (! h) {
260			printf(" WARNING -- sample %s does not have requested histogram. Next sample!\n",
261			s->Name()->Data());
262			continue;
263			}
264			if (! fEmptyHist) {
265			fEmptyHist = new TH1D(factor * scale * (*h));
266			fEmptyHist->Rebin(fNRebin);
267			fEmptyHist->Reset();
268			}
269
270			double nEvtsSelRaw = h->GetEntries();
271			double nEvtsSel = nEvtsSelRaw * factor * scale;
272			double nEvtsSelErr = TMath::Sqrt(nEvtsSelRaw) * factor * scale;
273
274	paus	1.3	printf(" -> %-40s - %14.0f %12.3f +- %8.3f %16.7f: %16.4f (x %f x %f - %p)\n",
275	paus	1.2	s->Name()->Data(),nEvts,nEvtsSel,nEvtsSelErr,s->Xsec(),lumi,factor,scale,(void)h);
276
277	paus	1.3	nTotRaw += nEvts;
278			nTot += nEvtsSel;
279			nTot2 += nEvtsSelErr*nEvtsSelErr;
280
281	paus	1.2	// scale it
282			TH1D hTmp = new TH1D(factor scale * (*h));
283			hTmp->Rebin(fNRebin);
284			// put it into our collection
285			fHists.push_back(hTmp);
286			// make a new histogram from the last added one and add the recent
287			if (i == 0)
288			fStackedHists.push_back(hTmp);
289			else {
290			TH1D hStackedTmp = new TH1D(fStackedHists[fStackedHists.size()-1]);
291			// remember the stored histogram is already rebinned
292			hStackedTmp->Add(hTmp);
293			fStackedHists.push_back(hStackedTmp);
294			}
295			}
296	paus	1.3	// Monte Carlo summary
297			printf(" %-40s - %14.0f %12.3f +- %8.3f %16.7f: %16.4f (x %f x %f)\n",
298			"== Monte Carlo Total ==",nTotRaw,nTot,TMath::Sqrt(nTot2),0.0,0.0,1.0,1.0);
299
300	paus	1.2
301			//------------------------------------------------------------------------------------------------
302			// go through the data samples
303			//------------------------------------------------------------------------------------------------
304			// loop through data samples and add them up, straight away
305			if (*fTask->NDataSamples() > 0)
306	paus	1.3	printf("\n Data \n");
307			nTotRaw = 0.0;
308			nTot = 0.0;
309			nTot2 = 0.0;
310	paus	1.2	for (UInt_t i=0; i<*fTask->NDataSamples(); i++) {
311			const Sample *s = fTask->GetDataSample(i);
312			// open file belonging to the data sample
313			TFile fif = new TFile((fTask->Dir()+slash+*s->File()).Data());
314			// make sure the file exists
315			if (fif->IsOpen() == kFALSE) {
316			printf(" WARNING -- sample %s does not have a histogram file. Continue without!\n",
317			s->Name()->Data());
318			continue;
319			}
320			else {
321			// read and determine general properties of this sample
322			TDirectory dirTmp = (TDirectory) gROOT->FindObject(dirFwk->Data());
323			if (dirTmp)
324			fif->cd(dirFwk->Data());
325			TH1D hAllEvts = (TH1D) gROOT->FindObject(allEvts->Data());
326			if (! hAllEvts) {
327			printf(" WARNING -- sample %s does not have a framework file. Next sample!\n",
328			s->Name()->Data());
329			continue;
330			}
331			else {
332			double nEvts = hAllEvts->GetEntries();
333
334			TH1D h = (TH1D) fif->FindObjectAny(hist);
335			if (! h) {
336			printf(" WARNING -- sample %s does not have requested histogram. Next sample!\n",
337			s->Name()->Data());
338			continue;
339			}
340			else
341			// rebin it
342			h->Rebin(fNRebin);
343
344			double nEvtsSelRaw = h->GetEntries();
345			double nEvtsSel = nEvtsSelRaw;
346			double nEvtsSelErr = TMath::Sqrt(nEvtsSelRaw);
347
348	paus	1.3	printf(" -> %-40s - %14.0f %12.3f +- %8.3f %16.7f: %16.4f (x %f)\n",
349	paus	1.2	s->Name()->Data(),nEvts,nEvtsSel,nEvtsSelErr,*s->Xsec(),fTargetLumi,1.0);
350
351	paus	1.3	nTotRaw += nEvts;
352			nTot += nEvtsSel;
353			nTot2 += nEvtsSelErr*nEvtsSelErr;
354	paus	1.2
355			// construct the complete data histogram
356			if (! fDataHist)
357			fDataHist = new TH1D(1.0 * (*h));
358			else
359			fDataHist->Add(h);
360			}
361			}
362			}
363	paus	1.3	// Data summary
364			printf(" %-40s - %14.0f %12.3f +- %8.3f %16.7f: %16.4f (x %f)\n\n",
365			"== Data Total =========",nTotRaw,nTot,TMath::Sqrt(nTot2),0.0,0.0,1.0);
366	paus	1.2
367			return;
368			}
369
370			//--------------------------------------------------------------------------------------------------
371			void PlotTask::FindHistMaximum()
372			{
373			// Find maximum of all histograms
374
375			// first check whether value was overwritten by hand
376			if (fHistMaximum>0.)
377			return;
378
379			// search through the single histograms
380			for (UInt_t i=0; i<fHists.size(); i++) {
381			if (fHists[i]->GetMaximum() > fHistMaximum) {
382			fHistMaximum = fHists[i]->GetMaximum();
383			fIdxHistMax = i;
384			}
385			}
386
387			// search through the stacked histograms
388			fHistMaximum = 0.;
389			for (UInt_t i=0; i<fStackedHists.size(); i++) {
390			if (fStackedHists[i]->GetMaximum() > fHistMaximum)
391			fHistMaximum = fStackedHists[i]->GetMaximum();
392			}
393
394			// check the data histogram if present
395			if (fDataHist && fDataHist->GetMaximum() > fHistMaximum) {
396			fHistMaximum = fDataHist->GetMaximum();
397			fIdxHistMax = fHists.size();
398			}
399
400			printf(" Histogram maximum is set to be: %f\n",fHistMaximum);
401
402			return;
403			}
404
405			//--------------------------------------------------------------------------------------------------
406			void PlotTask::OverlayEmptyHist() const
407			{
408			// Overlay an empty histogram onto the picture
409			if (fEmptyHist)
410			fEmptyHist->Draw("same");
411
412			return;
413			}
414
415			//--------------------------------------------------------------------------------------------------
416			void PlotTask::OverlayFrame() const
417			{
418			// Overlay a linear frame from user coordinates (0-100,0-100)
419
420			// create new transparent pad for the text
421			TPad *transPad = new TPad("transPad","Transparent Pad",0,0,1,1);
422			transPad->SetFillStyle(4000);
423			transPad->Draw();
424			transPad->cd();
425			// find out the right normalization to define the new range (histogram: 0,0 -> 100,100)
426			double xtot = 1./(1. - transPad->GetLeftMargin() - transPad->GetRightMargin());
427			double ytot = 1./(1. - transPad->GetBottomMargin() - transPad->GetTopMargin());
428			transPad->Range((xtottransPad->GetLeftMargin() -100 ),
429			(ytottransPad->GetBottomMargin()-100 ),
430			(xtottransPad->GetRightMargin() 100 + 100),
431			(ytottransPad->GetTopMargin() 100 + 100));
432			MDB(kGeneral,1)
433			printf(" Range: %f %f %f %f\n",
434			(xtottransPad->GetLeftMargin() -100 ),
435			(ytottransPad->GetBottomMargin()-100 ),
436			(xtottransPad->GetRightMargin() 100 + 100),
437			(ytottransPad->GetTopMargin() 100 + 100) );
438
439			// define the basic text to be used
440			TLatex* text = new TLatex();
441			text->SetTextFont (42);
442			text->SetTextAlign(12);
443			text->SetTextSize (0.03);
444
445			// define the basic box to be used
446			TBox *box = new TBox();
447			box->SetLineWidth(2);
448			// base coordinates for the text/boxes
449			double xCorner = fXLegend, yCorner = fYLegend, yDelLine = 4.;
450			double xIndent = 1.5*yDelLine;
451			// count samples with non empty legend
452			int nLegends = 0;
453			for (UInt_t i=*fTask->NSamples(); i>0; i--) {
454			// attach to the specific sample
455			const Sample *s = fTask->GetSample(i-1);
456			if (*s->Legend() != TString(" "))
457			nLegends++;
458			}
459			int iDat = 0;
460			if (fDataHist) {
461			iDat = 1;
462			nLegends++;
463			}
464
465			// set start values
466			fHistStyles->ResetStyle();
467			int iLeg = 0;
468			// loop through the sampels
469			for (UInt_t i=*fTask->NSamples(); i>0; i--) {
470			// attach to the specific sample
471			const Sample *s = fTask->GetSample(i-1);
472			// calculate corners for the text
473			double xText = xCorner+xIndent, yText = yCorner-float((iLeg+iDat)+0.5)*yDelLine;
474			// say what goes where
475			MDB(kGeneral,1)
476			printf(" Adding box at: (x1,y1,x2,y2) = (%6.2f,%6.2f,%6.2f,%6.2f)\n",
477			xCorner+0.15xIndent,yText-0.3yDelLine,
478			xCorner+0.85xIndent,yText+0.4yDelLine);
479			// plot the box
480			if (*s->Legend() != TString(" ")) {
481			const HistStyle *hStyle = fHistStyles->CurrentStyle();
482			box->SetFillStyle(0);
483			box->SetFillColor(hStyle->Color());
484			box->SetLineColor(hStyle->Color());
485			box->DrawBox(xCorner+0.15xIndent,yText-0.3yDelLine,
486			xCorner+0.85xIndent,yText+0.4yDelLine);
487			box->SetFillStyle(hStyle->FillStyle());
488			box->DrawBox(xCorner+0.15xIndent,yText-0.3yDelLine,
489			xCorner+0.85xIndent,yText+0.4yDelLine);
490
491			TString pText = *s->Legend();
492			MDB(kGeneral,1)
493			printf(" Adding text \"%s\" at: (x,y) = (%6.2f,%6.2f)\n",pText.Data(),xText,yText);
494			// set the proper text color
495			text->SetTextColor(hStyle->Color());
496			// plot the text
497			text->SetTextAlign(12);
498			text->DrawLatex(xText,yText,pText.Data());
499
500			fHistStyles->NextStyle();
501			// keep track of how many were printed
502			iLeg++;
503			}
504			}
505
506			// attach to the data sample
507			if (fDataHist) {
508			const Sample *s = fTask->GetDataSample(0);
509			// calculate corners for the text
510			double xText = xCorner+xIndent, yText = yCorner-float(1-0.7)*yDelLine;
511			// say what goes where
512			MDB(kGeneral,1)
513			printf(" Adding marker at: (x,y) = (%6.2f,%6.2f)\n",
514			xCorner+0.5*xIndent,yText);
515			// plot the marker
516			if (*s->Legend() != TString(" ")) {
517			const HistStyle *hStyle = fHistStyles->GetDataStyle();
518
519			// draw marker
520			TMarker m = new TMarker(xCorner+0.5xIndent,yText,20);
521			m->SetMarkerColor(hStyle->Color ());
522			m->SetMarkerSize (hStyle->MarkerSize ());
523			m->SetMarkerStyle(hStyle->MarkerStyle());
524			m->Draw();
525			// draw text
526			char l[1024];
527			sprintf(l,"%4.2f",fTargetLumi);
528			//TString pText = *s->Legend() + TString(" (#int L = ") + TString(l) + TString("/pb)");
529			TString pText = *s->Legend() + TString(" (") + TString(l) + TString("/pb)");
530			MDB(kGeneral,1)
531			printf(" Adding text \"%s\" at: (x,y) = (%6.2f,%6.2f)\n",pText.Data(),xText,yText);
532			// set the proper text color
533			text->SetTextColor(hStyle->Color());
534			// plot the text
535			text->SetTextAlign(12);
536			text->DrawLatex(xText,yText,pText.Data());
537			}
538			}
539
540			// Make sure the histogram frame is nice
541			box->SetFillStyle(0);
542			box->SetLineColor(kBlack);
543			box->DrawBox(0,0,100,100);
544
545			delete text;
546			delete box;
547
548			return;
549			}