//////////////////////////////////////////////////////////////////////////
//                                                                      //
// GFDstarZmvfns                                                   //
//
//   Author:      Gero Flucke
//   Date:        December 20th, 2004
//   last update: $Date: 2005/11/07 18:45:26 $ (UTC)
//   by:          $Author: flucke $
//                                                                      //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

#include "GFDstarZmvfns.h"
#include "GFAnalysis/GFAxisVariables.h"
#include "GFUtils/GFHistArray.h"
#include "GFUtils/GFHistManip.h"
#include "GFUtils/GFHistManager.h"

#include <TH1.h>
#include <TObjString.h>
#include <TLegend.h>
#include <TObjArray.h>

#include <iostream>
using std::cout;
using std::endl;

ClassImp(GFDstarZmvfns)

//_________________________________________________________
 GFDstarZmvfns::GFDstarZmvfns(const char *fileAdd, const char *filePath, 
			     const char *fileStartDir, const char *fileStartRes) : 
  GFDstarVirtualTheory(filePath),fFileStartDirNlo(fileStartDir), fFileStartResNlo(fileStartRes),
  fFileAdd(fileAdd), fName("ZMVFNS"), fNbFactor(1./1000.)
{
  if(fileAdd && strlen(fileAdd)){
//     if(strcmp(fileAdd, "_cteq5_grv") != 0) fName += fileAdd;
    if(strcmp(fileAdd, "_mrst04_grv") != 0) fName += fileAdd;
  } else {
    fName += "mrst03_afg04";
  }
}

//_________________________________________________________
 GFDstarZmvfns::~GFDstarZmvfns()
{
}

//_________________________________________________________
 void GFDstarZmvfns::DrawHists(const char *var, Int_t flag, Bool_t rel)
{
  // flag 0: final
  //      1: scale variation only (i.e. final!)
  TObjArray *vars = this->CreateVarArray(var);

  const Bool_t isBatch = fHistManager->SetBatch();
  fHistManager->Delete();
  for(Int_t i = 0; i < vars->GetEntriesFast(); ++i){
    GFHistArray *hists1 = NULL, *hists2 = NULL;
    TString vari;
    switch(flag){
    case 0:
    case 1:
      hists1 = this->CreateHistsFinal(vars->At(i)->GetName(), 0);
      if(!rel) hists2 = this->CreateHistsFinal(vars->At(i)->GetName(), -1);
      vari = "final";
      break;
    }
    if(hists1){
      hists1->First()->SetTitle(Form("%s, %s variation",hists1->At(0)->GetTitle(),vari.Data()));
      const Bool_t log = !rel && this->IsLog(vars->At(i)->GetName());
      const Int_t nHist = this->DrawTheHists(hists1, hists2, log, kTRUE);
      const Int_t layer = fHistManager->GetNumLayers()-1;
      TLegend *leg = fHistManager->AddLegend(layer, nHist, this->GetName(), kFALSE);
      leg->SetTextAlign(22);
      leg->AddEntry(hists1->At(1), "dir. + res.", "l");
      leg->AddEntry(hists1->At(0), Form("%s up", vari.Data()), "l");
      leg->AddEntry(hists1->At(2), Form("%s down", vari.Data()), "l");
      if(rel){
	hists1->At(0)->Divide(hists1->At(1));
	hists1->At(2)->Divide(hists1->At(1));
	hists1->At(1)->Divide(hists1->At(1));
	hists1->At(0)->SetYTitle("#sigma/#sigma(dir. + res.)");
	hists1->At(0)->GetYaxis()->CenterTitle();
      } else if(hists2){
	leg->AddEntry(hists2->At(0), "res. only (up)", "l");
	leg->AddEntry(hists2->At(1), "res. only", "l");
	leg->AddEntry(hists2->At(2), "res. only (down)", "l");
	delete hists2;
      }
      delete hists1;
    }
  }
  fHistManager->SetBatch(isBatch);
  fHistManager->Draw();
//   fHistManager->Update();// for the TExec added for log scale
  delete vars;
}

//_________________________________________________________
 GFHistArray* GFDstarZmvfns::CreateHistsFinal(const char *var, Int_t dirResFl)
{
  // < 0 res, > 0 dir
  // take max deviation from variation of one of mass and mu
  GFHistArray *histsScale = this->CreateHistsScale(var, dirResFl);
  if(!histsScale) return NULL;

  GFHistArray histsAllVar;
  histsAllVar.AddAll(histsScale);
  TH1 *hUpper = GFHistManip::CreateHistMinMax(&histsAllVar,  1);
  TH1 *hLower = GFHistManip::CreateHistMinMax(&histsAllVar, -1);

  GFHistArray *result = NULL;
  if(hUpper && hLower){
    hUpper->SetDirectory(NULL);
    hLower->SetDirectory(NULL);
    hUpper->SetLineStyle(4);
    hLower->SetLineStyle(3);
    TH1 *hCenter = histsScale->At(1);
    histsScale->AddAt(NULL, 1); // otherwise hCenter will be deleted at clean up below!
    result = new GFHistArray(3);
    result->AddAt(hUpper, 0);
    result->AddAt(hCenter, 1);
    result->AddAt(hLower, 2);
  } else {
    delete hUpper;
    delete hLower;
  }

  histsScale->Delete();
  delete histsScale;

  return result;
}

//_________________________________________________________
 GFHistArray* GFDstarZmvfns::CreateHistsScale(const char *var,Int_t dirResFlag)
{
  if (fFileAdd.Length() == 0 || fFileAdd == "_cteq5_grv") {
    return this->CreateHistsUpDown(var, "_mu2m125", "_mu1m4", dirResFlag);
  } else if  (fFileAdd == "_mrst04_grv") {
//     return this->CreateHistsUpDown(var, "_mu2m125", "_mu1m4", dirResFlag);
    return this->CreateHistsUpDown(var, "_mu2m1mf15", "_mu1m4", dirResFlag);
  } else {
    this->Warning("CreateHistsScale", "no variations: fFileAdd = %s, dirResFlag = %d", 
 		  fFileAdd.Data(), dirResFlag);
    return this->CreateHistsUpDown(var, "", "", dirResFlag);
  }
}

//_________________________________________________________
 TH1* GFDstarZmvfns::CreateHist(const char *var, const char *note)
{
  // get final dir+res for var with file-notifyer note (per bin width)
  TH1 *hDir = this->CreateHistDir(var, note);
  TH1 *hRes = this->CreateHistRes(var, note);

  if(!hDir || !hRes){
    delete hDir;
    delete hRes;
    return NULL;
  }
  const char *dummy = NULL;
  const char *hName = this->GetHistName(var, dummy);
  //                                                                no sumw2, but fill!
  TH1* hResult = GFHistManip::CreateAnalog(hDir, Form("%sSum", hName), kFALSE, kTRUE);
  hResult->SetDirectory(NULL);
  hResult->Add(hRes);

  delete hDir;
  delete hRes;

  return hResult;
}

//_________________________________________________________
 TH1* GFDstarZmvfns::CreateHistDir(const char *var, const char *note)
{
  // note ignored so far
  const char *fileExt = NULL;
  const char *hName = this->GetHistName(var, fileExt);

  TH1 *hDir = this->GetHist(hName, Form("%s%s%s%s.root", fFileStartDirNlo.Data(),
					fileExt, fFileAdd.Data(), note));

  if(hDir){
    hDir->SetName(Form("%sDir", hName));
    hDir->SetTitle(GFAxisVariables::GetAxisLabel(var));
    hDir->SetXTitle(GFAxisVariables::GetAxisLabel(var));

    this->Scale(hDir, var);
  
    if(strcmp(var, "Dphi") == 0 || strcmp(var, "DphiNlo") == 0) {
//     GFHistManip::RadToDegXaxis(hDir, kTRUE);
      this->SetUpperHistEgde(hDir, 180.);
    } else if(strcmp(var, "xGam") == 0) {
      this->SetUpperHistEgde(hDir, 1.0);
    } else {
// get rid of error bars                               sumw2,  fill
      TH1* hNew = GFHistManip::CreateAnalog(hDir, hDir->GetName(), kFALSE, kTRUE);
      hNew->SetDirectory(NULL);
      delete hDir;
      hDir = hNew;
    }
  }
  return hDir;
}


//_________________________________________________________
 TH1* GFDstarZmvfns::CreateHistRes(const char *var, const char *note)
{
  // note ignored so far
  const char *fileExt = NULL;
  const char *hName = this->GetHistName(var, fileExt);

  TH1 *hRes = this->GetHist(hName, Form("%s%s%s%s.root", fFileStartResNlo.Data(),
					fileExt, fFileAdd.Data(), note));

  if(hRes){
    hRes->SetName(Form("%sRes", hName));
    hRes->SetTitle(GFAxisVariables::GetAxisLabel(var));
    hRes->SetXTitle(GFAxisVariables::GetAxisLabel(var));

    this->Scale(hRes, var);
  
    if(strcmp(var, "Dphi") == 0 || strcmp(var, "DphiNlo") == 0) {
//     GFHistManip::RadToDegXaxis(hRes, kTRUE);
      this->SetUpperHistEgde(hRes, 180.);
    } else if(strcmp(var, "xGam") == 0) {
      this->SetUpperHistEgde(hRes, 1.0);
    } else {
// get rid of error bars                              sumw2,  fill
      TH1* hNew = GFHistManip::CreateAnalog(hRes, hRes->GetName(), kFALSE, kTRUE);
      hNew->SetDirectory(NULL);
      delete hRes;
      hRes = hNew;
    }
  }
  return hRes;
}

//_________________________________________________________
 GFHistArray* GFDstarZmvfns::CreateHistsFinalNorm(const char *var, Int_t dirResFlag)
{
  // < 0 res, > 0 dir 
  TObjArray upDowns;
  upDowns.SetOwner();
  upDowns.Add(new TObjString(""));
// following ifs as in TotalCrossSec
  if (fFileAdd.Length() == 0 || fFileAdd == "_cteq5_grv") {
    upDowns.Add(new TObjString("_mu2m125"));
    upDowns.Add(new TObjString("_mu1m4"));
  } else if (fFileAdd == "_mrst04_grv") {
    upDowns.Add(new TObjString("_mu2m1mf15"));
    upDowns.Add(new TObjString("_mu1m4"));
  } else {
    this->Warning("GFDstarZmvfns", "no variations, fFileAdd = %s, dirResFlag = %d",
		  fFileAdd.Data(), dirResFlag);
    upDowns.Add(new TObjString(""));
    upDowns.Add(new TObjString(""));
  }

  return this->CreateHistsFinalNormInt(var, upDowns, dirResFlag);
}

//_________________________________________________________
 TArrayD GFDstarZmvfns::TotalCrossSec(ECrossSecType dstarJet, Int_t dirResFlag)
{
  // returning array of length 3: total cross section, upper and lower bound
  // 
  // dirResFlag < 0: resolved only
  //            > 0: direct only
  //            = 0: sum

  TObjArray upDowns;
  upDowns.SetOwner();
  upDowns.Add(new TObjString(""));
// following ifs as in CreateHistsFinalNorm
  if (fFileAdd.Length() == 0 || fFileAdd == "_cteq5_grv") {
    upDowns.Add(new TObjString("_mu2m125"));
    upDowns.Add(new TObjString("_mu1m4"));
  } else if (fFileAdd == "_mrst04_grv") {
    upDowns.Add(new TObjString("_mu2m1mf15"));
    upDowns.Add(new TObjString("_mu1m4"));
  } else {
    this->Warning("TotalCrossSec", "no variations, fFileAdd = %s, dirResFlag = %d",
		  fFileAdd.Data(), dirResFlag);
    upDowns.Add(new TObjString(""));
    upDowns.Add(new TObjString(""));
  }

  return this->TotalCrossSecInt(dstarJet, dirResFlag, upDowns);
}

//_________________________________________________________
 const char* GFDstarZmvfns::GetHistName(const char *variable, 
					    const char *& fileExt) const
{
  TString var(variable);

  const char *result = NULL;
  // single diff.
  if (var == "eta")          result = "h20";
  else if (var == "pt")     result = "h21";
  else if (var == "wGammaP")result = "h22";
  else if (var == "zDs")    result = "h23";
  else if (var == "etapt11")result = "h24";
  else if (var == "etapt22")result = "h25";
  else if (var == "etapt33")result = "h26";
  if (result) {
    fileExt = "_incl";
    return result;
  }

  // double diff.
  if (var == "EtaDs")        result =  "h20";
  else if (var == "EtaJet") result =  "h21";
  else if (var == "PtDs")   result =  "h22";
  else if (var == "PtJet")  result =  "h23";
  else if (var == "Dphi" || var == "DphiNlo")   result =  "h24";
  else if (var == "PtDsJet")result =  "h25";
  else if (var == "Deta")   result =  "h26";
  else if (var == "MDsJet") result =  "h27";
  else if (var == "xGam")   result =  "h28";
  if (result) {
    fileExt = "_jeta";
    return result;
  }
  
  this->Error("GetHistName", "variable %s not known", variable);
  fileExt = NULL;
  return NULL;
}


//_________________________________________________________
 Bool_t GFDstarZmvfns::Scale(TH1 *hist, const char *var) const
{
  // scale hist to fit cross section in nb

  if(!hist) return kFALSE;

  GFHistManip::Scale(hist, fNbFactor);
  hist->SetYTitle(this->YAxisTitleBinWidth("#sigma [nb]", GFAxisVariables::GetAxisLabel(var)));

  return kTRUE;
}

//_________________________________________________________
 const char* GFDstarZmvfns::GetName() const
{
  return fName.Data();
//   if(fFileAdd.Length() == 0){
//     return "ZMVFNS";//"NLO (ZMVFNS)";
//   } else {
//     static TString buf("ZMVFNS" + fFileAdd);
//     return buf;
//   }
}

//_________________________________________________________
 void GFDstarZmvfns::SetFileAdd(const char *add)
{
  fFileAdd = add; 
  fName = "ZMVFNS";

  if(add && strlen(add)){
//     if(strcmp(add, "_cteq5_grv") != 0) fName += add;
    if(strcmp(add, "_mrst04_grv") != 0) fName += add;
  } else {
    fName += "mrst03_afg04";
  }
}

//_________________________________________________________
 void GFDstarZmvfns::SetUpperHistEgde(TH1 *&hist, Axis_t upper, Bool_t sw2) const
{
  if(!hist) return;
  if(hist->GetDimension() != 1) {
    this->Error("SetUpperHistEdge", "%d-D not yet foreseen, bit %s is!", 
		hist->GetDimension(), hist->GetName());
    return;
  }
  Axis_t *bins = new Axis_t[hist->GetNbinsX() + 1];
  for(Int_t i = 0; i < hist->GetNbinsX(); ++i){
    bins[i] = hist->GetBinLowEdge(i+1);
  }
  if(upper > bins[hist->GetNbinsX()-1]) {
    bins[hist->GetNbinsX()] = upper;
    
    TDirectory *oldHistsDir = hist->GetDirectory();
    hist->SetDirectory(NULL);
    TString tit(Form("%s;%s;%s", hist->GetTitle(), hist->GetXaxis()->GetTitle(),
		     hist->GetYaxis()->GetTitle()));
    TH1 *hnew = new TH1F(hist->GetName(), tit, hist->GetNbinsX(), bins);
    for(Int_t i = 1; i <= hist->GetNbinsX(); ++i){
      hnew->SetBinContent(i, hist->GetBinContent(i));
      if(sw2 && hist->GetSumw2N()) hnew->SetBinError(i, hist->GetBinError(i));
    }
    hnew->SetDirectory(oldHistsDir);
    delete hist;
    hist = hnew;
  } else {
    this->Error("SetUpperHistEdge", "%f is too low!", upper); 
  }
  delete bins;
}