#include "GFDstarDrawXSec.h"
#include "GFDstarFmnr.h"
#include "GFDstarZmvfns.h"
#include "GFDstarCascade.h"
#include "GFDstarPythia.h"
#include "GFAnalysis/GFDstarHistsAnalysisMc.h"
#include "GFUtils/GFHistManager.h"
#include "GFUtils/GFHistArray.h"
#include "GFUtils/GFHistManip.h"

#include <TFile.h>
#include <TH1.h>
#include <TStyle.h>
#include <TBox.h>
#include <TExec.h>
#include <TGaxis.h>
#include <TList.h>
#include <TLine.h>
#include <TLegend.h>
#include <TLegendEntry.h>
#include <TLatex.h>
#include <TPaveText.h>

#include <iostream>
#include <string.h>

using namespace std;

ClassImp(GFDstarDrawXSec)

//___________________________________________________________
 GFDstarDrawXSec::GFDstarDrawXSec
//   fAnalysisMc->SetBatch(kTRUE);
  fAnalysisMc = fPythia->GetHistsAnalysisMc();

  fHistManager->DrawDiffStyle(kFALSE);

  fFileOutput = NULL;
}

//___________________________________________________________
 GFDstarDrawXSec::~GFDstarDrawXSec()
{
  this->Clear();

  delete fFileData;
  delete fFmnr;
  delete fZmvfns;
  delete fCascade;
  delete fPythia;
  //  delete fAnalysisMc;
  delete fHistManager;

  delete fFileOutput;
}

//___________________________________________________________
 TFile* GFDstarDrawXSec::SetFileOutput(const char *name)
{
// opens file with given name and prepares to store theory curves...
  delete fFileOutput;
  fFileOutput = NULL;
  if (name && strlen(name)) {
    fFileOutput = TFile::Open(name, "UPDATE");
  }

  return fFileOutput;
}

//___________________________________________________________
 void GFDstarDrawXSec::Clear(Option_t* opt)
{
  //  fHistManager->Delete(); problem: hists are in the list twice!
  fHistManager->Clear();
}

//___________________________________________________________
 TH1 *GFDstarDrawXSec::DrawNloDstar(const char *variables, GFDstarVirtualTheory *nlo, 
				   GFDstarVirtualTheory *nlo2, Bool_t norm)
{
// one ore several variables, e.g. variables = "1D,xGam"
// return pointer to first drawn hist of last variable
  if(!nlo) return NULL;
  this->Clear();

  const TObjArray *vars = nlo->CreateVarArray(variables);
  TH1 *hResult = NULL;
  Bool_t anyLog = kFALSE;
  for (Int_t i = 0; i < vars->GetEntriesFast(); ++i) {
    const char *var = vars->At(i)->GetName();
    GFHistArray *nloHists = (norm ? nlo->CreateHistsFinalNorm(var) : 
			     nlo->CreateHistsFinal(var));
    TH1 *hData = this->GetHistData(var, (norm ? "n" : ""));
    if (norm) this->Warning("DrawNloDstar", "no sys. errror yet for normalised!");
    TH1 *hDataAllErr = this->GetHistData(var, (norm ? "n" : "e")); //FIXME: 'ne', not 'n'
    if(!nloHists || !nloHists->At(1) || !hData || !hDataAllErr){
      this->Error("DrawNloDstar", "unexpected input for %s", var);
      if (nloHists) nloHists->Delete();
      delete nloHists; delete hData; delete hDataAllErr;
      continue;
    }

    Float_t min, max;
    this->GetMinMax(var, min, max, norm);
    hDataAllErr->SetMinimum(min);
    hDataAllErr->SetMaximum(max);
    
    nloHists->At(1)->SetName(var);
    const Int_t layer = this->AddHistsNlo(hDataAllErr, nloHists, nlo->GetName());
    delete nloHists;

    if(nlo2){
      GFHistArray *nloHists2 = (norm ? nlo2->CreateHistsFinalNorm(var) : 
				nlo2->CreateHistsFinal(var));
      nloHists2->At(1)->SetName(var);
      this->AddHistsNlo2(nloHists2, nlo2->GetName());
      delete nloHists2;
    }

    GFDstarVirtualTheory::ECrossSecType type = nlo->IsDsJetOrDs(var);
    TString header(nlo->GetHeader(type));
    if (type == GFDstarVirtualTheory::kInclDoubleDiff) {
      TString title(hData->GetTitle());
      const char *label = hData->GetXaxis()->GetTitle(); // needed for 'norm' only
      const Ssiz_t pos = title.Index(norm ? Form("d%s", label) : "[nb]");
      if(pos != kNPOS){
	title.Remove(0, pos + (norm ? strlen(label) + 4 : 6));
	title.Remove(title.Length()-1, 1);
      } else {
	title = "double differential";
      }
      (header += ": ") += title;
    }
  
    const Bool_t isLog = nlo->IsLog(var);
    if (isLog) anyLog = kTRUE;
    TLegend *legend = this->AddHistsData(hData, hDataAllErr, layer, isLog, header);
    this->SetLegendPos(legend, var);
    hResult = hDataAllErr;
  }

  fHistManager->Draw();
  if (anyLog) fHistManager->Update();
  
  return hResult;
}

//___________________________________________________________
 TH1 *GFDstarDrawXSec::DrawNloDstarJet(const char *variables, const char *forwBack, GFDstarVirtualTheory *nlo, GFDstarVirtualTheory *nlo2, Bool_t norm)
{
// one ore several variables, e.g. variables = "1D,xGam"
// return pointer to first drawn hist of last variable

  if(!nlo) return NULL;
  this->Clear();

  const TObjArray *vars = nlo->CreateVarArray(variables);
  TH1 *hResult = NULL;
  Bool_t anyLog = kFALSE;
  for (Int_t i = 0; i < vars->GetEntriesFast(); ++i) {
    const char *var = vars->At(i)->GetName();
    TString mergeVar(Form("%s%s", var, forwBack));
    GFHistArray *nloHists = (norm ? nlo->CreateHistsFinalNorm(mergeVar) 
			     : nlo->CreateHistsFinal(mergeVar));
    TH1 *hHadCorr = NULL;
    if(nlo->ShouldApplyHadCor() || (nlo2 && nlo2->ShouldApplyHadCor())){
      if (nlo->IsDsJetOrDs(var) == GFDstarVirtualTheory::kDiJet) {
	hHadCorr = fAnalysisMc->CreateHadCorrDiJet(var);
	if (norm) GFHistManip::Scale(hHadCorr, 1./fAnalysisMc->TotalHadCorrDiJet()[0]);
      } else {
	hHadCorr = fAnalysisMc->CreateHadCorr(var, forwBack);
	if (norm) GFHistManip::Scale(hHadCorr, 1./fAnalysisMc->TotalHadCorr(forwBack)[0]);
      }
    }
    TH1 *hData = this->GetHistData(var, (norm ? "n" : ""));
    if (norm) this->Warning("DrawNloDstarJet", "no sys. errror yet for normalised!");
    TH1 *hDataAllErr = this->GetHistData(var, (norm ? "n" : "e")); //FIXME: 'ne', not 'n'
    if(!nloHists || !nloHists->At(1) || !hData || !hDataAllErr){
      this->Error("DrawNloDstarJet", "unexpected input for %s", var);
      if (nloHists) nloHists->Delete();
      delete nloHists; delete hData; delete hDataAllErr; delete hHadCorr;
      continue;
    }
    if(hHadCorr){
      hHadCorr->SetDirectory(NULL);//otherwise hist name re-use
      // get rid of errors of hHadCorr:                      no sumw2 but fill
      TH1 *hTmp = GFHistManip::CreateAnalog(hHadCorr, hHadCorr->GetName(), kFALSE, kTRUE);
      delete hHadCorr;
      hHadCorr = hTmp; hTmp = NULL;
      hHadCorr->SetDirectory(NULL);
    }

    TH1 *hNloPlain = NULL;
    TString nloTitle(nlo->GetName());
    if(nlo->ShouldApplyHadCor()) {
      hNloPlain = GFHistManip::CreateAnalog(nloHists->At(1), 
					    Form("%sPlain", nloHists->At(1)->GetName()),
					    kFALSE, kTRUE); // no sumw2 but fill
      hNloPlain->SetDirectory(NULL);
      for(Int_t i = 0; i < nloHists->GetEntriesFast(); ++i){
	if(nloHists->At(i)){
	  nloHists->At(i)->Multiply(hHadCorr);
	  nloHists->At(i)->SetName(Form("%sHadCorr", nloHists->At(i)->GetName()));
	}
      }
      nloTitle = Form("%s #otimes Had", nlo->GetName());
    }

    Float_t min, max;
    this->GetMinMax(var, min, max, norm);
    hDataAllErr->SetMinimum(min);
    hDataAllErr->SetMaximum(max);
    nloHists->At(1)->SetName(var);
    const Int_t layer = this->AddHistsNlo(hDataAllErr, nloHists, nloTitle, hNloPlain,
					  nlo->GetName());
    delete nloHists;
    
    if(nlo2){
      GFHistArray *nloHists2 = (norm ? nlo2->CreateHistsFinalNorm(mergeVar)
				: nlo2->CreateHistsFinal(mergeVar));
      TH1 *hNloPlain2 = NULL;
      TString nloTitle2(nlo2->GetName());
      if(nlo2->ShouldApplyHadCor()) {
	TH1 *h = nloHists2->At(1);
	if(!nlo->ShouldApplyHadCor()){// without had corr. only once needed  /// no sumw2 but fill
	  hNloPlain2 = GFHistManip::CreateAnalog(h, Form("%sPlain", h->GetName()),kFALSE,kTRUE);
	  hNloPlain2->SetDirectory(NULL);
	}
	for(Int_t i = 0; i < nloHists2->GetEntriesFast(); ++i){
	  if(nloHists2->At(i)){
	    nloHists2->At(i)->Multiply(hHadCorr);
	    nloHists2->At(i)->SetName(Form("%sHadCorr", nloHists2->At(i)->GetName()));
	  }
	}
	nloTitle2 = Form("%s #otimes Had", nlo2->GetName());
      }
      nloHists2->At(1)->SetName(var);
      this->AddHistsNlo2(nloHists2, nloTitle2.Data(), hNloPlain2, nlo2->GetName());
      delete nloHists2;
    }
    delete hHadCorr;
    
    GFDstarVirtualTheory::ECrossSecType type = nlo->IsDsJetOrDs(var);
    TString header(nlo->GetHeader(type));
    if (type != GFDstarVirtualTheory::kDstarJet && type != GFDstarVirtualTheory::kDiJet) {
      header = "unknown or non-jet variable";
    }

    const Bool_t isLog = nlo->IsLog(var);
    if (isLog) anyLog = isLog;
    TLegend *legend = this->AddHistsData(hData, hDataAllErr, layer, isLog, header);
    this->SetLegendPos(legend, var);
    hResult = hDataAllErr;
  }

  fHistManager->Draw();
  if (anyLog) fHistManager->Update();

  return hResult;
}

//___________________________________________________________
 TH1 *GFDstarDrawXSec::DrawFmnrDstar(const char *var)
{
  return this->DrawNloDstar(var, fFmnr);
}

//___________________________________________________________
 TH1 *GFDstarDrawXSec::DrawFmnrDstarJet(const char *var, const char *forwBack)
{
  return this->DrawNloDstarJet(var, forwBack, fFmnr);
}

//___________________________________________________________
 TH1 *GFDstarDrawXSec::DrawMasslessDstar(const char *var)
{
  return this->DrawNloDstar(var, fZmvfns);
}

//___________________________________________________________
 TH1 *GFDstarDrawXSec::DrawMasslessDstarJet(const char *var, const char *forwBack)
{
  return this->DrawNloDstarJet(var, forwBack, fZmvfns);
}

//___________________________________________________________
 TH1 *GFDstarDrawXSec::DrawCascade(const char *var)
{
  // with and without jet!
  return this->DrawNloDstar(var, fCascade);
}

//___________________________________________________________
 TH1 *GFDstarDrawXSec::DrawCascadePythia(const char *var, Bool_t direct)
{
  // with and without jet!
  const Bool_t oldBatch = fHistManager->SetBatch(kTRUE);
  TH1 *h = this->DrawNloDstar(var, fCascade, fPythia);
  TH1 *hPy2 = (direct ? fPythia->CreateHistDir(var) : fPythia->CreateHistRes(var));
  this->AddExtraHist(hPy2, Form("%s (%s)", fPythia->GetName(), (direct ? "dir." : "res.")));

  if (fFileOutput) {
    TDirectory *oldDir = gDirectory;
    fFileOutput->cd();
    TString name(Form("%s%s%s",var,fPythia->GetName(), (direct?"Dir":"Res")));
    GFHistManip::MakeUniqueName(name);
    hPy2->Write(name);
    oldDir->cd();
  }

  fHistManager->SetBatch(oldBatch);
  fHistManager->Draw();
  fHistManager->Update(); // needed?

  return h;
}

//___________________________________________________________
 TH1 *GFDstarDrawXSec::DrawPythia(const char *var, Bool_t direct)
{
  // with and without jet!
  const Bool_t oldBatch = fHistManager->SetBatch(kTRUE);
  TH1 *h = this->DrawNloDstar(var, fPythia);
  TH1 *hPy2 = (direct ? fPythia->CreateHistDir(var) : fPythia->CreateHistRes(var));
  this->AddExtraHist(hPy2, Form("%s (%s)", fPythia->GetName(), (direct ? "dir." : "res.")));

  if (fFileOutput) {
    TDirectory *oldDir = gDirectory;
    fFileOutput->cd();
    TString name(Form("%s%s%s",var, fPythia->GetName(), (direct?"Dir":"Res")));
    GFHistManip::MakeUniqueName(name);
    hPy2->Write(name);
    oldDir->cd();
  }

  fHistManager->SetBatch(oldBatch);
  fHistManager->Draw();
  fHistManager->Update(); // needed?

  return h;
}

//___________________________________________________________
 TH1 *GFDstarDrawXSec::DrawFmnrPythia(const char *var, Bool_t direct)
{
  // with and without jet!
  const Bool_t oldBatch = fHistManager->SetBatch(kTRUE);
  TH1 *h = NULL;
  switch(fPythia->IsDsJetOrDs(var)){
  case GFDstarVirtualTheory::kDstarJet:
  case GFDstarVirtualTheory::kDiJet:
    h = this->DrawNloDstarJet(var, "", fFmnr, fPythia);
    break;
  case GFDstarVirtualTheory::kInclusive:
    h = this->DrawNloDstar(var, fFmnr, fPythia);
    break;
  case GFDstarVirtualTheory::kInclDoubleDiff:
  case GFDstarVirtualTheory::kUnknown:
    this->Error("DrawFmnrPythia", "IsDsJetOrDs returns %d", fPythia->IsDsJetOrDs(var));
  }

  fHistManager->SetBatch(oldBatch);
  if(h){
    TString title(h->GetTitle());
    Ssiz_t pos = title.Index("]");
    if(pos != kNPOS) {
      title.Replace(0, pos+3, "");
      pos = title.Last(')');
      if(pos != kNPOS) title.Replace(pos, 1, "");
    }
    h->SetTitle(title);
    TH1 *hPy2 = (direct ? fPythia->CreateHistDir(var) : fPythia->CreateHistRes(var));
    this->AddExtraHist(hPy2, Form("%s (%s)", fPythia->GetName(), (direct ? "dir." : "res.")));
    
    fHistManager->Draw();
    fHistManager->Update(); // needed?
  }

  return h;
}

//___________________________________________________________
 TH1 *GFDstarDrawXSec::DrawNloDirRes(const char *var, GFDstarVirtualTheory *nlo, 
				    Bool_t direct)
{
  // draw just the NLO and the direct or resolved part of it
  // with and without jet!
  if(!nlo) return NULL;
  const Bool_t oldBatch = fHistManager->SetBatch(kTRUE);
  TH1 *h = NULL;
  switch(nlo->IsDsJetOrDs(var)){
  case GFDstarVirtualTheory::kDstarJet:
  case GFDstarVirtualTheory::kDiJet:
    h = this->DrawNloDstarJet(var, "", nlo);
    break;
  case GFDstarVirtualTheory::kInclusive:
  case GFDstarVirtualTheory::kInclDoubleDiff:
    h = this->DrawNloDstar(var, nlo);
    break;
  case GFDstarVirtualTheory::kUnknown:
    this->Error("DrawNloDirRes", "%s is of unknown type", var);
  }

  fHistManager->SetBatch(oldBatch);
  if(h){
    TString title(h->GetTitle());
    Ssiz_t pos = title.Index("]");
    if(pos != kNPOS) {
      title.Replace(0, pos+3, "");
      pos = title.Last(')');
      if(pos != kNPOS) title.Replace(pos, 1, "");
    }
    h->SetTitle(title);
    TH1 *hNlo2 = (direct ? nlo->CreateHistDir(var) : nlo->CreateHistRes(var));
    this->AddExtraHist(hNlo2, Form("%s (%s)", nlo->GetName(), (direct ? "dir." : "res.")));
    
    fHistManager->Draw();
    fHistManager->Update(); // needed?
  }

  return h;
}

//___________________________________________________________
 TH1 *GFDstarDrawXSec::DrawFmnrCascDstar(const char *var)
{
  return this->DrawNloDstar(var, fFmnr, fCascade);
}

//___________________________________________________________
 TH1 *GFDstarDrawXSec::DrawFmnrCascDstarJet(const char *var, const char *forwBack)
{
  return this->DrawNloDstarJet(var, forwBack, fFmnr, fCascade);
}

//___________________________________________________________
 TH1 *GFDstarDrawXSec::DrawFmnrMasslessDstar(const char *var)
{
  return this->DrawNloDstar(var, fFmnr, fZmvfns);
}

//___________________________________________________________
 TH1 *GFDstarDrawXSec::DrawFmnrMasslessDstarJet(const char *var, const char *forwBack)
{
  return this->DrawNloDstarJet(var, forwBack, fFmnr, fZmvfns);
}

//___________________________________________________________
 TArrayD GFDstarDrawXSec::TotalDstarJet(GFDstarVirtualTheory *nlo) const
{
  if(!nlo) return TArrayD(4);
  TArrayD result(nlo->TotalCrossSec(GFDstarVirtualTheory::kDstarJet));//D*+jet
  if(nlo->ShouldApplyHadCor()){
    const Double_t hadCorr = fAnalysisMc->TotalHadCorr("").At(0);
    for(Int_t i = 0; i < result.GetSize(); ++i){
      result[i] *= hadCorr;
    }

    cout << "---------------------------------------------------------------n"
	 << "Total X-sec D*+o. jet (had. corr.) = " 
	 << result[0] << " + " << result[1]-result[0] 
	 << " - " << result[0]-result[2] << "n" 
	 << "---------------------------------------------------------------n"
	 << endl;
  }

  return result;
}

//___________________________________________________________
 TArrayD GFDstarDrawXSec::TotalDstarDiJet(GFDstarVirtualTheory *nlo) const
{
  if(!nlo) return TArrayD(4);
  TArrayD result(nlo->TotalCrossSec(GFDstarVirtualTheory::kDiJet));//D*+2 jet
  if(nlo->ShouldApplyHadCor()){
    const Double_t hadCorr = fAnalysisMc->TotalHadCorrDiJet().At(0);
    for(Int_t i = 0; i < result.GetSize(); ++i){
      result[i] *= hadCorr;
    }

    cout << "---------------------------------------------------------------n"
	 << "Total X-sec D* dijet (had. corr.) = " 
	 << result[0] << " + " << result[1]-result[0] 
	 << " - " << result[0]-result[2] << "n" 
	 << "---------------------------------------------------------------n"
	 << endl;
  }

  return result;
}

//___________________________________________________________
 TArrayD GFDstarDrawXSec::TotalFmnrDstarJet()
{
  return this->TotalDstarJet(fFmnr);

//   if(!fFmnr) return TArrayD(4);
//   TArrayD result(fFmnr->TotalCrossSec(kTRUE));// true: D*+jet
//   const Double_t hadCorr = fAnalysisMc->TotalHadCorr("").At(0);
//   for(Int_t i = 0; i < result.GetSize(); ++i){
//     result[i] *= hadCorr;
//   }
  
//   cout << "---------------------------------------------------------------n"
//        << "Total X-sec D*+jet (had. corr.) = " 
//        << result[0] << " + " << result[1]-result[0] 
//        << " - " << result[0]-result[2] << "n" 
//        << "---------------------------------------------------------------n"
//        << endl;

//   return result;
}

//___________________________________________________________
 void GFDstarDrawXSec::DrawNloDataRatio(const char *vars, GFDstarVirtualTheory *nlo, 
				       GFDstarVirtualTheory *nlo2, Bool_t hadCorr,
				       Int_t dirResFlag)
{
  if(!nlo) return;
  this->Clear();

  const Int_t layer = fHistManager->GetNumLayers();
  TObjArray *arr = nlo->CreateVarArray(vars);

  for(Int_t nHist = 0; nHist < arr->GetEntriesFast(); ++nHist){ 
    const char *var = arr->At(nHist)->GetName();
    GFHistArray *nloHists = nlo->CreateHistsFinal(var, dirResFlag);
    TH1 *hDataAllErr = static_cast<TH1*>(fFileData->Get(Form("%sAllErr", var)));
    if(!nloHists || !nloHists->At(1) || !hDataAllErr){
      this->Error("DrawNloDataRatio", "unexpected input for %s", var);
      continue;
    }
    hDataAllErr->SetDirectory(NULL);
    if(strcmp(var, "Dphi") == 0) GFHistManip::MergeBins(hDataAllErr, 5, 6, "widthe");

    TH1 *hHadCorr = NULL;
    const Bool_t hadCorr1 = hadCorr && nlo->ShouldApplyHadCor() 
      && nlo->IsDsJetOrDs(var) == GFDstarVirtualTheory::kDstarJet;
    const Bool_t hadCorr2 = hadCorr && nlo2 && nlo2->ShouldApplyHadCor()
      && nlo2->IsDsJetOrDs(var) == GFDstarVirtualTheory::kDstarJet;;
    if (hadCorr1 || hadCorr2) {
      hHadCorr = fAnalysisMc->CreateHadCorr(var, "");
      TH1 *hTmp = GFHistManip::CreateAnalog(hHadCorr, Form("%sHadCor",var), kFALSE, kTRUE);
      delete hHadCorr;
      hHadCorr = hTmp;
    }

    if(nloHists->At(2)) nloHists->At(2)->Divide(nloHists->At(1));
    if(nloHists->At(0)) nloHists->At(0)->Divide(nloHists->At(1));
    if(hadCorr1) nloHists->At(1)->Multiply(hHadCorr);
    
    Double_t *relErr = new Double_t[hDataAllErr->GetNbinsX()];
    for(Int_t i = 1; i <= hDataAllErr->GetNbinsX(); ++i){
      relErr[i-1] = hDataAllErr->GetBinError(i)/hDataAllErr->GetBinContent(i);
    }
    hDataAllErr->Divide(nloHists->At(1));
    for(Int_t i = 1; i <= hDataAllErr->GetNbinsX(); ++i){
      hDataAllErr->SetBinError(i, relErr[i-1]*hDataAllErr->GetBinContent(i));
    }
    delete [] relErr; 
    relErr = NULL;
    
    GFHistArray *nlo2Hists = NULL;
    if (nlo2) {
      nlo2Hists = nlo2->CreateHistsFinal(var, dirResFlag);
      if(!nlo2Hists || !nlo2Hists->At(1)){
	this->Error("DrawNloDataRatio", "unexpected input for %s", var);
	delete nlo2Hists;
	nlo2Hists = NULL;
      } else {
	for (Int_t i = 0; i < nlo2Hists->GetEntriesFast(); ++i) {
	  if(nlo2Hists->At(i)) {
	    nlo2Hists->At(i)->Divide(nloHists->At(1));
	    if(hadCorr2) nlo2Hists->At(i)->Multiply(hHadCorr);
	  }
	}
      }
    }
    
    nloHists->At(1)->Divide(nloHists->At(1));
    
    TString nloName(Form("%s%s", nlo->GetName(), (hadCorr1?" #otimes Had":"")));
    hDataAllErr->SetMarkerStyle(8);// 8 means dot
    fHistManager->AddHist(hDataAllErr, layer, "data", "lp");
    hDataAllErr->SetTitle("data / theory");
    hDataAllErr->SetYTitle(Form("data / theory (%s)", nloName.Data()));
			   
    nloHists->At(1)->SetLineColor(kRed);
    if(nloHists->At(0) &&  nloHists->At(2)){
      for(Int_t i = 1; i <= nloHists->At(1)->GetNbinsX(); ++i){
	const Float_t x1 = nloHists->At(1)->GetBinLowEdge(i);
	const Float_t x2 = x1 + nloHists->At(1)->GetBinWidth(i);
	const Float_t y1 = nloHists->At(2)->GetBinContent(i);
	const Float_t y2 = nloHists->At(0)->GetBinContent(i);
	TBox *box = new TBox(x1, y1, x2, y2);
	box->SetFillStyle(1001);
	box->SetFillColor(kYellow);
	nloHists->At(1)->TAttLine::Copy(*box);
	fHistManager->AddObject(box, layer, nHist);
	if(i == 1) fHistManager->GetLegendOf(layer, nHist)->AddEntry(box, nloName, "lf");
      }
      fHistManager->AddObject(nloHists->At(1), layer, nHist, "SAME");
    } else {
      fHistManager->AddHistSame(nloHists->At(1), layer, nHist, nloName, "lf");
    }
    
    if (nlo2Hists) {
      TString nlo2Name(Form("%s%s", nlo2->GetName(), (hadCorr2?" #otimes Had":"")));
      if (nlo2Hists->At(0) && nlo2Hists->At(2)) {
	for(Int_t i = 1; i <= nlo2Hists->At(1)->GetNbinsX(); ++i){
	  const Float_t x1 = nlo2Hists->At(1)->GetBinLowEdge(i);
	  const Float_t x2 = x1 + nlo2Hists->At(1)->GetBinWidth(i);
	  const Float_t y1 = nlo2Hists->At(2)->GetBinContent(i);
	  const Float_t y2 = nlo2Hists->At(0)->GetBinContent(i);
	  TBox *box = new TBox(x1, y1, x2, y2);
#if ROOT_VERSION_CODE >= ROOT_VERSION(4,1,4)
	  box->SetFillStyle(3445);
#else
	  box->SetFillStyle(3004); // problem with gv!
#endif
	  box->SetFillColor(kBlue);
	  fHistManager->AddObject(box, layer, nHist);
	  if(i == 1) {
	    fHistManager->GetLegendOf(layer, nHist)->AddEntry(box, nlo2Name, "lf");
	    nlo2Hists->At(1)->SetLineColor(box->GetFillColor());
	  }
	  nlo2Hists->At(1)->TAttLine::Copy(*box);
	}
	fHistManager->AddObject(nlo2Hists->At(1), layer, nHist, "SAME");
      } else {
	fHistManager->AddObject(nlo2Hists->At(1), layer, nHist, "SAME");
	fHistManager->GetLegendOf(layer, nHist)->AddEntry(nlo2Hists->At(1), nlo2Name, "l");
      }
      delete nlo2Hists;
    }
    fHistManager->AddObject(hDataAllErr, layer, nHist, "SAME");// draw again on top
    
    delete nloHists;
    delete hHadCorr;
  }
  
  fHistManager->Draw();
  delete arr;
}



//___________________________________________________________
 Int_t GFDstarDrawXSec::AddHistsNlo(TH1 *hData, GFHistArray *nloHists, 
				   const char *nloLabel, TH1 *hPlain,const char *plainLabel)
{
  if (fFileOutput) {
    TDirectory *oldDir = gDirectory;
    fFileOutput->cd();
    TString name(nloHists->At(1)->GetName());
    name += nloLabel;
    GFHistManip::MakeUniqueName(name);
    nloHists->At(1)->Write(name);
    if (nloHists->At(0)) nloHists->At(0)->Write(name+"Upper");
    if (nloHists->At(2)) nloHists->At(2)->Write(name+"Lower");
    oldDir->cd();
  }

//   const Int_t layer = fHistManager->GetNumLayers();
  const Int_t layer = TMath::Max(fHistManager->GetNumLayers() - 1, 0);

//   fHistManager->AddHist(hData, layer, "Data", "lp");//for being the first: legend, min/max,...
  fHistManager->AddHist(hData, layer, "H1", "lp");//for being the first: legend, min/max,...
  const Int_t nHist = fHistManager->GetNumHistsOf(layer)-1;
  hData->SetOption(""); // overwrite old default option from file (it was E1?)

  const Color_t white = 10;
  const Color_t colErrBand = (fColourCode == 1 ? kYellow : 7);//3=lightgreen, 8=darkgreen
  const Color_t colLine = (fColourCode == 1 ? kRed : kBlue); // 7 = cyan
  Bool_t errorBand = kFALSE;
  if(nloHists->At(0) && nloHists->At(2)){
    errorBand = kTRUE;
    nloHists->At(0)->SetFillStyle(1001);
    nloHists->At(0)->SetFillColor(colErrBand);
    nloHists->At(0)->SetLineColor(colErrBand);
    nloHists->At(0)->SetLineWidth(1);
    nloHists->At(0)->SetLineStyle(1);
   
    nloHists->At(2)->SetFillStyle(1001);
    nloHists->At(2)->SetFillColor(white);
    nloHists->At(0)->TAttLine::Copy(*nloHists->At(2));

    nloHists->At(1)->SetFillColor(colErrBand);
  }
  nloHists->At(1)->SetLineColor(colLine);

  if(nloHists->At(0)) fHistManager->AddHistSame(nloHists->At(0), layer, nHist);
  fHistManager->AddHistSame(nloHists->At(1), layer, nHist, nloLabel, (errorBand ? "lf" : "l"));
  if(nloHists->At(2)) fHistManager->AddHistSame(nloHists->At(2), layer, nHist);
  if(hPlain) {
    hPlain->SetLineStyle(2);
    hPlain->SetLineColor(nloHists->At(1)->GetLineColor());
    fHistManager->AddHistSame(hPlain, layer, nHist, plainLabel, "l");
  }
  TH1 *hTmp = GFHistManip::CreateAnalog(nloHists->At(1), 
					Form("%s2", nloHists->At(1)->GetName()),
					kFALSE, kTRUE); // no sumw2 but fill
  hTmp->SetDirectory(NULL);
  nloHists->At(1)->TAttLine::Copy(*hTmp);
  fHistManager->AddHistSame(hTmp, layer, nHist);// redraw on top of 


  return layer;
}

//___________________________________________________________
 Int_t GFDstarDrawXSec::AddHistsNlo2(GFHistArray *nloHists, const char *nloLabel,
				    TH1 *hPlain,const char *plainLabel)
{
  // second calculation drawn via TBoxes
  if (fFileOutput) {
    TDirectory *oldDir = gDirectory;
    fFileOutput->cd();
    TString name(nloHists->At(1)->GetName());
    name += nloLabel;
    GFHistManip::MakeUniqueName(name);
    nloHists->At(1)->Write(name);
    if (nloHists->At(0)) nloHists->At(0)->Write(name+"Upper");
    if (nloHists->At(2)) nloHists->At(2)->Write(name+"Lower");
    oldDir->cd();
  }

  const Int_t layer = fHistManager->GetNumLayers() -1;// go on with same layer!
  if(layer < 0) return layer;
  const Int_t nHist = fHistManager->GetNumHistsOf(layer)-1;
  const Color_t color = (fColourCode == 1 ? kBlue : kRed);
  nloHists->At(1)->SetLineColor(color);
  nloHists->At(1)->SetLineStyle(2);

  if(nloHists->At(0) && nloHists->At(2)){
#if ROOT_VERSION_CODE >= ROOT_VERSION(4,1,4)
    const Style_t fillStyle = 3445;
#else
    const Style_t fillStyle = 3004; // problem with gv!
#endif
    
    nloHists->At(0)->SetLineColor(color);
    nloHists->At(0)->SetLineWidth(1);
    nloHists->At(0)->SetLineStyle(1);
    
    nloHists->At(0)->TAttLine::Copy(*nloHists->At(2));

    fHistManager->AddHistsSame(nloHists, layer, nHist);

    for(Int_t i = 1; i <= nloHists->At(1)->GetNbinsX(); ++i){
      const Float_t x1 = nloHists->At(1)->GetBinLowEdge(i);
      const Float_t x2 = x1 + nloHists->At(1)->GetBinWidth(i);
      const Float_t y1 = nloHists->At(2)->GetBinContent(i);
      const Float_t y2 = nloHists->At(0)->GetBinContent(i);
      TBox *box = new TBox(x1, y1, x2, y2);// log scale has been fixed in ROOT cvs 22.12.2004
      box->SetFillStyle(fillStyle);
      box->SetFillColor(color);
      box->SetLineColor(color);
      box->SetLineStyle(nloHists->At(1)->GetLineStyle());//1);
      box->SetLineWidth(nloHists->At(1)->GetLineWidth());
      fHistManager->AddObject(box, layer, nHist);
      if(i == 1) fHistManager->GetLegendOf(layer, nHist)->AddEntry(box, nloLabel, "lf");
    }
  } else {
    fHistManager->AddHistSame(nloHists->At(1), layer, nHist, nloLabel, "l");
  }

  if(hPlain) {
    hPlain->SetLineStyle(2);
    hPlain->SetLineColor(nloHists->At(1)->GetLineColor());
    fHistManager->AddHistSame(hPlain, layer, nHist, plainLabel, "l");
  }

  return layer;
}

//___________________________________________________________
 Int_t GFDstarDrawXSec::AddExtraHist(TH1 *h, const char *label)
{
  // add h in current layer with legend label <label>
  if(!h) return -1;
  const Int_t layer = fHistManager->GetNumLayers() -1;// go on with same layer!
  if(layer < 0) return layer;
  const Int_t nHist = fHistManager->GetNumHistsOf(layer)-1;

  GFHistArray *prevHists = fHistManager->GetHistsOf(layer, nHist);
  if(!prevHists || !prevHists->Last()) return -2;
  TH1 *hDirPlusRes = prevHists->Last();

  hDirPlusRes->TAttLine::Copy(*h);
  hDirPlusRes->TAttFill::Copy(*h);
  hDirPlusRes->TAttMarker::Copy(*h);
  h->SetLineStyle(3);
  
  TLegend *leg = fHistManager->GetLegendOf(layer, nHist);
  if(leg){
    const Int_t size = leg->GetListOfPrimitives()->GetSize();
    leg->SetY1(leg->GetY1() - (leg->GetY2()-leg->GetY1())/size);
  }

  fHistManager->AddHistSame(h, layer, nHist, label, "l");

  return layer;
}

//___________________________________________________________
 TLegend* GFDstarDrawXSec::AddHistsData(TH1 *hData, TH1 *hDataAllErr, Int_t layer,
				       Bool_t log, const char *header)
{
  if(!hData || !hDataAllErr) return NULL;

  const Int_t nHist = fHistManager->GetNumHistsOf(layer)-1;
  if(log){
    TExec *exe = new TExec("log", "TVirtualPad::Pad()->SetLogy(kTRUE);");
    fHistManager->AddObject(exe, layer, nHist);
#if ROOT_VERSION_CODE >= ROOT_VERSION(4,1,4)
  } else {
    hDataAllErr->SetMinimum(0.); // why is that needed?
#endif
  }

  // draw the hists via AddObject not via AddHistsSame to paint them on top of the TBoxes
  TExec *exe1 = new TExec("enderrorsize", "gStyle->SetErrorX(0.)");
  fHistManager->AddObject(exe1, layer, nHist);

  hDataAllErr->SetMarkerStyle(8);
//   hDataAllErr->SetMarkerSize(1.6); //  to be set outside via gStyle
  hDataAllErr->SetMarkerColor(hData->GetLineColor());
  hDataAllErr->SetTitle(hData->GetTitle());
  hData->TAttLine::Copy(*hDataAllErr);
  fHistManager->AddObject(hData, layer, nHist, "SAMEE1");

  TExec *exe2 = new TExec("enderrorsizeBack", "gStyle->SetErrorX()");
  fHistManager->AddObject(exe2, layer, nHist);

  fHistManager->AddObject(hDataAllErr, layer, nHist, "SAME");

  // now we need to redraw axes, do via TGaxis
  // since gPad->RedrawAxis() results in ugly screen and .(e)ps (superimposed labels)
  const Float_t xMin = hDataAllErr->GetXaxis()->GetXmin();
  const Float_t xMax = hDataAllErr->GetXaxis()->GetXmax();
  const Float_t yMin = (log ? hDataAllErr->GetMinimum() : TMath::Min(0., hDataAllErr->GetMinimum()));
  const Float_t yMax = hDataAllErr->GetMaximum();

  TGaxis *xAxis = new TGaxis(xMin, yMin, xMax, yMin, xMin, xMax, 510, "U");
  TGaxis *yAxis = new TGaxis(xMin, yMin, xMin, yMax, yMin, yMax, 510, (log ? "UG" : "U"));
  TLine  *line  = new TLine (xMax, yMin, xMax, yMax);//could use hist-draw option "][" instead ?

  fHistManager->AddObject(xAxis, layer, nHist);
  fHistManager->AddObject(yAxis, layer, nHist);
  fHistManager->AddObject(line,  layer, nHist);

  Double_t yText = gStyle->GetPadBottomMargin()+0.04;
  if(TString(hDataAllErr->GetName()).Contains("Dphi")) {
    yText += 0.16;
  }
//   TLatex *h1 = new TLatex(gStyle->GetPadLeftMargin()+0.04, yText, "H1");
//   h1->SetTextFont(62);
//   h1->SetNDC();
//   //  h1->SetTextSize(0.11);
//   h1->SetTextSize(0.10);
//   h1->SetTextColor(12);
//   fHistManager->AddObject(h1, layer, nHist);

//   TLatex *prel = new TLatex(gStyle->GetPadLeftMargin()+0.15, h1->GetY(), "(prel.)"); 
// //   prel->SetTextFont(132);
//   prel->SetTextFont(h1->GetTextFont());
//   prel->SetNDC();
//   prel->SetTextSize(0.07);
//   prel->SetTextColor(h1->GetTextColor());
//   fHistManager->AddObject(prel, layer, nHist);

  if(header){
//     TText *title = new TLatex(0.34, 0.925, header);
    TText *title = new TLatex(0.5, 0.95, header);
    title->SetTextFont(12);
    title->SetTextAlign(22);
    title->SetNDC();
    title->SetTextSize(0.065);
    fHistManager->AddObject(title, layer, nHist);
  }
  TLegend *legend = fHistManager->GetLegendOf(layer, nHist);
  if(legend) legend->SetBorderSize(1);

  return legend; 
}

/////////////////////////////////////////////////////////////////
 Bool_t GFDstarDrawXSec::GetMinMax(const char *var, Float_t &min, Float_t &max, 
				  Bool_t norm) const
{
  // need to give minimum for log hists due to usage of TGaxis...
  min = max = -1111.; // resets
  TString v(var);
  if(v == "eta"){
    max = (norm ? 0.9 : 9.);
  } else if(v == "pt"){
    max =  (norm ? 1.0 : 13.5);
    min = (norm ? 0.001 : 0.004);
  } else if(v == "wGammaP"){
    max =  (norm ? 0.023 : 0.23);
  } else if(v == "zDs"){
    max =  (norm ? 3.5 : 22.); 
  } else if(v == "PtDs"){
    min = (norm ? 0.002 : 0.005);
    max = (norm ? 0.6 : 2.5);
  } else if(v == "PtJet"){
    min = (norm ? 0.002 : 0.005);
    max = (norm ? 0.6 : 2.5);
  } else if(v == "EtaDs"){
    max =  (norm ? 0.75 : 2.5); //  3.0;
  } else if(v == "EtaJet"){
    max =  (norm ? 0.75 : 2.5); //  3.0;
  } else if(v == "Deta"){
    max =  (norm ? 0.4 :1.25);
  } else if(v == "Dphi" || v == "DphiNlo"){
    max =  (norm ? 0.04 : 0.11);//0.115;
  } else if(v == "PtDsJet"){
    max =  1.15;
    if (norm) this->Warning("GetMinMax", "ignore request for normed %s!", var);
  } else if(v == "MDsJet"){
    max =  0.43;
    if (norm) this->Warning("GetMinMax", "ignore request for normed %s!", var);
  } else if(v == "xGam"){
    max =  (norm ? 2.5 : 7.5);//9.4;
  } else if(v.BeginsWith("pteta")){
    max =  2.5;
    min = 0.025;
    if (norm) this->Warning("GetMinMax", "ignore request for normed beginning with %s!", var);
  } else if(v == "etapt11"){
    max = (norm ? 0.6 : 6.0); //= 3.5; high for ZMVFNS
  } else if(v == "etapt22"){
    max = (norm ? 0.3 : 2.2); //= 1.7; high for ZMVFNS
  } else if(v == "etapt33"){
    max = (norm ? 0.11 : 0.55);
//   } else if(v.BeginsWith("etapt")){
//     max = 3.3;
  } else if(v == "xGamDijet"){
    max =  (norm ? 4. : 8.7);//9.4;
  } else if(v == "PtDJet" || v == "PtOthJet"){
    min = 0.004;
    max = 1.5;
    if (norm) this->Warning("GetMinMax", "ignore request for normed %s!", var);
  } else if(v == "EtaDJet" || v == "EtaOthJet"){
    max = 1.9;
    if (norm) this->Warning("GetMinMax", "ignore request for normed %s!", var);
  } else if(v == "DetaDijet"){
    max = 1.1;
    if (norm) this->Warning("GetMinMax", "ignore request for normed %s!", var);
  } else if(v == "PtDijet"){
    max = 0.95;
    if (norm) this->Warning("GetMinMax", "ignore request for normed %s!", var);
  } else {
    return kFALSE;
  }

  return kTRUE;
}

//____________________________________________________________
 void GFDstarDrawXSec::SetLegendPos(TLegend *leg, const char *var) const
{
  // set coordinates
  if(!leg) return;

  const Float_t legWidth = 0.366;
  const TString v(var);
  if (v.BeginsWith("xGam") || v.BeginsWith("Dphi")) {
    leg->SetX1(gStyle->GetPadLeftMargin() + 0.04);
    leg->SetX2(leg->GetX1() + legWidth);
  } else {    
    leg->SetX2(1.-gStyle->GetPadRightMargin());
    leg->SetX1(leg->GetX2() - legWidth);
  }
  leg->SetY2(1.-gStyle->GetPadTopMargin());

  const Int_t numLegEnt = leg->GetListOfPrimitives()->GetSize();
  leg->SetY1(leg->GetY2() - TMath::Max(1,numLegEnt)*0.06);

  leg->SetTextSize(0.041);
  
//   if(fColourCode == 1) {  // FIXME...
//     const Int_t layer = fHistManager->GetNumLayers() -1;// go on with same layer!
//     if(layer < 0) {
//       this->Error("SetLegendPos", "problem with layer = %d", layer);
//     } else {
//       const Int_t nHist = fHistManager->GetNumHistsOf(layer)-1;
//       TPaveText *pdf = new TPaveText(leg->GetX1()+0.01, leg->GetY1()-0.01, 
// 				     leg->GetX2()-0.01, leg->GetY1()-0.11,"NDC");
//       pdf->SetBorderSize(0);
// //       pdf->SetFillColor(kWhite);
//       pdf->SetFillStyle(0);
//       pdf->SetTextFont(132);
//       pdf->AddText("calculations with");
//       pdf->AddText("CTEQ5M");
//       fHistManager->AddObject(pdf, layer, nHist);
//     }
//   }
}

//____________________________________________________________
 Bool_t GFDstarDrawXSec::Divide(const GFHistArray *hsToDivide, const TH1 *hDiv) const
{
  // divide all hists in hsToDivide (may have empty slots) by hDiv, 
  // relative errors of hsToDivide are not changed
  // false if a problem (given pointers NULL, no hists to divide or a hist has 
  // a bin with zero content but non-zero error)

  if (!hsToDivide || hsToDivide->IsEmpty() || !hDiv) return kFALSE;

  Bool_t result = kTRUE;
  for (Int_t iHist = 0; iHist < hsToDivide->GetEntriesFast(); ++iHist) {
    TH1 *h = hsToDivide->UncheckedAt(iHist);
    if (!h) continue;

    if (h->GetSumw2N() == 0) h->Sumw2();
    TArrayD relError(h->GetSumw2N());
    for (Int_t iBin = 0; iBin < relError.GetSize(); ++iBin){ // fine even for TH2/3 (?)
      relError[iBin] = h->GetBinError(iBin);
      const Stat_t binCont = h->GetBinContent(iBin);
      if (binCont) relError[iBin] /= binCont;
      else if (relError[iBin] != 0.) result = kFALSE; // 0./0. might be OK
    }
    h->Divide(hDiv);
    for (Int_t iBin = 0; iBin < relError.GetSize(); ++iBin) {// fine even for TH2/3 (?)
      h->SetBinContent(iBin, relError[iBin] * h->GetBinContent(iBin));
    }
  }

  return result;
}

//____________________________________________________________
 TH1* GFDstarDrawXSec::GetHistData(const char *var, Option_t *option)
{
// hist with name 'var' from data file
// if option contains "e": add "AllErr" to name
// if option contains "n": add "Norm" (after eventual adding of "AllErr")
  if (!var || !fFileData) return NULL;

  const TString opt(option);
  TString varStr(var);
  if (opt.Contains("e", TString::kIgnoreCase)) varStr += "AllErr";
  if (opt.Contains("n", TString::kIgnoreCase)) varStr += "Norm";


  TH1 *result = NULL;
#if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,8)
  fFileData->GetObject(varStr.Data(), result);
#else
  result = static_cast<TH1*>(fFileData->Get(varStr.Data()));
#endif
  return result;
}