Plot/src/PlotTask.cc

// $Id: PlotTask.cc,v 1.6 2011/05/21 00:23:12 klute 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 <TTree.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.),
  fNBins       (100)
{
  // 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::SetDrawExp(const char* draw, const char* sel)
{
  fDrawExp = TString(draw);
  fSelExp  = TString(sel);

  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,int(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;

    TDirectory *fid = 0;
    TString histname = hist;
    if (TString(hist).Contains("/")) {
      TObjArray *substra = TString(hist).Tokenize("/");
      fid = (TDirectory*)fif->FindObjectAny(((TObjString*)(substra->At(0)))->GetString());
      histname = ((TObjString*)substra->At(1))->GetString();
      delete substra;
    }
    else {
      fid = fif;
    }
        
    
    TH1D *h = dynamic_cast<TH1D*>(fid->FindObjectAny(histname));
    
    //histogram doesn't exist, try to find TTree instead
    if (!h) {
      TTree *htree = dynamic_cast<TTree*>(fid->FindObjectAny(histname));
      if (!htree) {
        printf(" WARNING -- sample  %s  does not have requested histogram. Next sample!\n",
        s->Name()->Data());
        continue;
      }
      
      TString histname("htempmc_");
      histname += i;
      h = new TH1D(histname,histname,fNBins,fHistXMinimum,fHistXMaximum);
      TString drawexp = fDrawExp + TString(">>") + histname;
      htree->Draw(drawexp,fSelExp);
    }
    
    if (! fEmptyHist) {
      fEmptyHist = new TH1D(factor * scale * (*h));
      fEmptyHist->Rebin(fNRebin);
      fEmptyHist->Reset();
    }

    double nEvtsSelRaw = h->GetSumOfWeights();
    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();
        
        TDirectory *fid = 0;
        TString histname = hist;
        if (TString(hist).Contains("/")) {
          TObjArray *substra = TString(hist).Tokenize("/");
          fid = (TDirectory*)fif->FindObjectAny(((TObjString*)(substra->At(0)))->GetString());
          histname = ((TObjString*)substra->At(1))->GetString();
          delete substra;
        }
        else {
          fid = fif;
        }
            
        TH1D *h = dynamic_cast<TH1D*>(fid->FindObjectAny(histname));
        //histogram doesn't exist, try to find TTree instead
        if (!h) {
          TTree *htree = dynamic_cast<TTree*>(fid->FindObjectAny(histname));
          if (!htree) {
            printf(" WARNING -- sample  %s  does not have requested histogram. Next sample!\n",
            s->Name()->Data());
            continue;
          }
          
          TString histname("htempdata_");
          histname += i;
          h = new TH1D(histname,histname,fNBins,fHistXMinimum,fHistXMaximum);
          TString drawexp = fDrawExp + TString(">>") + histname;
          //printf ("Draw(%s, %s);\n",drawexp.Data(),fSelExp.Data());
          htree->Draw(drawexp,fSelExp);
        }
        else
          // rebin it
          h->Rebin(fNRebin);

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