Plot/src/PlotTask.cc

// $Id: PlotTask.cc,v 1.7 2011/05/27 14:21:11 bendavid 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;

const TH1D *PlotTask::sPuWeights = 0;

//--------------------------------------------------------------------------------------------------
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     (65.),
  fYLegend     (94.),
  fNBins       (100),
  fPuTarget    (0)
{
  // 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;
    }
    
    //set pileup weights
    if (fPuTarget) {
      if (sPuWeights) {
        delete sPuWeights;
      }
      
      TH1D *pusource = (TH1D*)dirTmp->Get("hNPU")->Clone();
      pusource->Scale(1.0/pusource->GetSumOfWeights());
      
      sPuWeights = new TH1D( (*fPuTarget) / (*pusource) );
      
    }
      
    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;
}

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