//////////////////////////////////////////////////////////////////////////
//                                                                      //
// GFDstarVirtualTheory                                                      //
//
//   Author:      Gero Flucke
//   Date:        November 3rd, 2004
//   last update: $Date: 2005/12/14 12:57:30 $ (UTC)
//   by:          $Author: flucke $
//                                                                      //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

#include "GFDstarVirtualTheory.h"
#include "GFAnalysis/GFAxisVariables.h"
#include "GFAnalysis/GFDstarHistsAnalysisMc.h"
#include <TArrayD.h>
#include <TH1.h>
#include <TFile.h>
#include <TObjArray.h>
#include <TObjString.h>
#include <TExec.h>
#include <TLegend.h>
#include "GFUtils/GFHistManager.h"
#include "GFUtils/GFHistManip.h"
#include "GFUtils/GFHistArray.h"

#include <iostream>
using namespace std;

ClassImp(GFDstarVirtualTheory)

GFDstarVirtualTheory::GFDstarVirtualTheory(const char *filePath) : 
  fHistManager(new GFHistManager), fFilePath(filePath), fOpenFiles(new TObjArray)
{

}

 GFDstarVirtualTheory::~GFDstarVirtualTheory()
{
  fOpenFiles->Delete(); // close all files
  delete fOpenFiles;
  fOpenFiles = NULL;

  fHistManager->Delete(); // delete hists
  delete fHistManager;
  fHistManager = NULL;
}

//_________________________________________________________
 void GFDstarVirtualTheory::DrawHistsNorm(const char *var, Int_t dirResFlag)
{
  TObjArray *vars = this->CreateVarArray(var);

  const Bool_t isBatchOld = fHistManager->SetBatch();
  Bool_t anyLog = kFALSE;
  fHistManager->Delete();

  for(Int_t i = 0; i < vars->GetEntriesFast(); ++i){
    GFHistArray *hs = this->CreateHistsFinalNorm(vars->At(i)->GetName(), dirResFlag);
    if (hs){
//       hs->First()->SetTitle(Form("%s",hists1->At(0)->GetTitle()));
      const Bool_t log = this->IsLog(vars->At(i)->GetName());
      if (log) anyLog = kTRUE;
      const Int_t nHist = this->DrawTheHists(hs, NULL, log, kTRUE);
      const Int_t layer = fHistManager->GetNumLayers()-1;
      TLegend *leg = fHistManager->AddLegend(layer, nHist, this->GetName(), kFALSE);
      leg->SetTextAlign(22);
      if (hs->At(0)) leg->AddEntry(hs->At(0), "final up", "l");
      leg->AddEntry(hs->At(1), "central", "l");
      if (hs->At(2)) leg->AddEntry(hs->At(2), "final down", "l");
      delete hs;
    }
  }
  fHistManager->SetBatch(isBatchOld);
  fHistManager->Draw();
  if (anyLog) fHistManager->Update();// for the TExec added for log scale
  delete vars;
}

//_________________________________________________________
 void GFDstarVirtualTheory::DrawSimple(const char *var, const char *note)
{
  TObjArray *vars = this->CreateVarArray(var);

  fHistManager->Delete();
  TString header(this->GetName());
  if(note && strlen(note)) header += Form("  (%s)", note);

  for(Int_t i = 0; i < vars->GetEntriesFast(); ++i){
    TH1 *h = this->CreateHist(vars->At(i)->GetName(), note);
    if(!h) continue;
    fHistManager->AddHist(h, 0, "dir. + res.");
    if(this->IsLog(vars->At(i)->GetName())){
      TExec *exe = new TExec("log", "TVirtualPad::Pad()->SetLogy(kTRUE);");
      fHistManager->AddObject(exe, 0, i);
    }
    h = this->CreateHistRes(vars->At(i)->GetName(), note);
    if(h){
      fHistManager->AddHistSame(h, 0, i, "res.");
    }
    fHistManager->AddLegend(0, i, header, kFALSE);
  }
  delete vars;
  fHistManager->DrawDiffStyle(kTRUE);
  fHistManager->Draw();
}

//_________________________________________________________
 TH1* GFDstarVirtualTheory::DrawResRatio(const char *var, const char *note)
{
  // the last ratio hist is returned
  TObjArray *vars = this->CreateVarArray(var);

  fHistManager->Delete();
  TH1 *hLastRatio = NULL;
  for(Int_t i = 0; i < vars->GetEntriesFast(); ++i){
    TH1 *h = this->CreateHist(vars->At(i)->GetName(), note);
    TH1 *hRes = this->CreateHistRes(vars->At(i)->GetName(), note);
    if(h && hRes) {
      hRes->Divide(h);
      hLastRatio = GFHistManip::CreateAnalog(hRes, Form("hResRatio%s", vars->At(i)->GetName()),
					     kFALSE, kTRUE);
      hLastRatio->SetDirectory(NULL);
      hLastRatio->SetYTitle("");
      fHistManager->AddHist(hLastRatio, 0, "resolved fraction");
      fHistManager->AddLegend(0, i, this->GetName(), kFALSE);
    } 
    delete h;
    delete hRes;
   
  }
  delete vars;
  fHistManager->DrawDiffStyle(kTRUE);
  fHistManager->Draw();

  return hLastRatio;
}

//_________________________________________________________
 Int_t GFDstarVirtualTheory::DrawTheHists(GFHistArray *hists, GFHistArray *hists2, 
				      Bool_t log, Bool_t add)
{
  if(!hists || !hists->GetEntriesFast()) return 0;
  const Int_t layer = 0;
  if(!add) fHistManager->Delete();
  const Int_t nHist = (fHistManager->GetNumLayers() ? fHistManager->GetNumHistsOf(layer) : 0);

  fHistManager->AddHist(hists->At(0), layer);
  fHistManager->AddHistSame(hists->At(1), layer, nHist);
  if(hists->At(2)) fHistManager->AddHistSame(hists->At(2), layer, nHist);
  if(hists2){
    fHistManager->AddHistSame(hists2->At(0), layer, nHist);
    TIter next(hists2);
    while(TH1 *h = static_cast<TH1*>(next())) h->SetLineColor(8);
    hists2->At(1)->SetLineStyle(2);//1+hists2->At(1)->GetLineStyle());
    fHistManager->AddHistSame(hists2->At(1), layer, nHist);
    if(hists2->At(2)) fHistManager->AddHistSame(hists2->At(2), layer, nHist);
  }
  if(log && hists->First()){
    TExec *exe = new TExec("log", "TVirtualPad::Pad()->SetLogy(kTRUE);");
    fHistManager->AddObject(exe, layer, nHist);
  }
  fHistManager->DrawDiffStyle(kFALSE);
  fHistManager->Draw();

  return nHist;
}

//_________________________________________________________
 TH1* GFDstarVirtualTheory::CreateHistNorm(const char *var, const char *note, 
					  Int_t dirResFlag)
{
  // normalised cross section with appendix 'note', dir/res or dir+res(default)
  // (normalised to inclusive D* or D*+other jet or D* tagged dijet)

  TH1 *h = this->CreateHist(var, note);
  if (h) {
    ECrossSecType type = this->IsDsJetOrDs(var);
    Double_t tot = this->TotalCrossSec(type, note, dirResFlag);
    if (tot) {
//       if (type == kDstarJet || type == kDiJet) {
// 	if (!mc) this->Warning("CreateHistNorm", "without had. corr. for %s", var);
// 	else {
// 	  TH1 *hHadCorr = NULL;
// 	  if (type == kDstarJet) {
// 	    tot *= mc->TotalHadCorr("").At(0);
// 	    hHadCorr = mc->CreateHadCorr(var, "");
// 	  } else { // dijet!
// 	    tot *= mc->TotalHadCorrDiJet().At(0);
// 	    hHadCorr = mc->CreateHadCorrDiJet(var);
// 	  }
// 	  // get rid of errors of hHadCorr:                   no sumw2 but fill
// 	  TH1 *hCor2 = GFHistManip::CreateAnalog(hHadCorr, "tmpName", kFALSE, kTRUE);
// 	  h->Multiply(hCor2);
// 	  h->SetTitle(Form("%s (had. cor. applied)", h->GetTitle()));
// 	  delete hHadCorr; delete hCor2;
// 	}
//       }
      h->Scale(1./tot);
      h->SetName(Form("%sNorm", h->GetName()));
      TString yTit(Form("1/#sigma_{%s} %s",this->GetHeader(type),h->GetYaxis()->GetTitle()));
      this->RemoveNb(yTit);
      h->SetYTitle(yTit);
    } else {
      delete h;
      h = NULL;
    }
  }

  return h;
}

//_________________________________________________________
 GFHistArray* GFDstarVirtualTheory::CreateHistsFinalNormInt(const char *var, 
						  const TObjArray &upDowns, 
						  Int_t dirResFlag)
{
  // upDowns is array of string of 'notes' (appendices) indicating 
  // the scale variation (given as argument to CreateHistNorm)
  // assuming the first entry is the 'default'
  // if upDowns contains only 1 entry, resulting array has only 1 hist (at 1)

  GFHistArray hists;
  for(Int_t i = 0; i < upDowns.GetEntriesFast(); ++i){
    TH1 *h = this->CreateHistNorm(var, upDowns[i]->GetName());
    if (h) hists.Add(h);
  }
  if (hists.IsEmpty()) return NULL;

  GFHistArray *result = new GFHistArray(3);
  result->AddAt(hists.At(0), 1); // take care not to delete (cf. belowe!)

  if (hists.GetEntriesFast() > 1) {
    TH1 *hUpper = GFHistManip::CreateHistMinMax(&hists,  1);
    hUpper->SetDirectory(NULL);
    hUpper->SetLineStyle(4);
    result->AddAt(hUpper, 0);
    TH1 *hLower = GFHistManip::CreateHistMinMax(&hists, -1);
    hLower->SetDirectory(NULL);
    hLower->SetLineStyle(3);
    result->AddAt(hLower, 2);
  }

  hists.RemoveAt(0); // to prevent deletion of central value put into result array
  hists.Delete();
  return result;
}

//_________________________________________________________
 Bool_t GFDstarVirtualTheory::ShouldApplyHadCor() const
{
  // whether or not hadronisation corrections should be applied or not.
  // default is true
  return kTRUE;
}

//_________________________________________________________
 TH1* GFDstarVirtualTheory::GetHist(const char *hName, const char *fileName)
{
  if (!hName || strlen(hName) == 0) return NULL;

  const TString fullFileName(Form("%s/%s", fFilePath.Data(), fileName));

  TFile *file = NULL;
  TIter nextFile(fOpenFiles);
  while(TObject *obj = nextFile()){
    if(obj->GetName() == fullFileName){
      file = static_cast<TFile*>(obj);
      break;
    }
  }

  if(!file) {
    TDirectory *oldDir = gDirectory;
    TFile *oldFile = gFile;
    file = TFile::Open(fullFileName);
    if(file) fOpenFiles->Add(file);
    gDirectory = oldDir;
    gFile = oldFile;
  }
  if(!file){
    this->Error("GetHist", "No file %s available", fileName);
    return NULL;
  }

  TObject *object = (hName ? file->Get(hName) : NULL); // ROOT crashes on NULL, fixed in 4.01_03
  if(object && object->InheritsFrom(TH1::Class())){
    TH1 *h = static_cast<TH1*>(object);
    h->SetDirectory(NULL);
    return h;
  } else {
    this->Error("GetHist", "%s not in file %s or not a TH1", hName, fullFileName.Data());
    delete object;
    return NULL;
  }
}

//_________________________________________________________
 GFHistArray* GFDstarVirtualTheory::CreateHistsUpDown(const char *var, 
						  const char *flagUp, const char *flagDown,
						  Int_t dirResFlag)
{
  // dirResFlag < 0: resolved only
  //            > 0: direct only
  //            = 0: sum
  TH1 *hCenter = NULL, *hUpper = NULL, *hLower = NULL;
  if(dirResFlag < 0){
    hCenter = this->CreateHistRes(var);
    hUpper  = this->CreateHistRes(var, flagUp);
    hLower  = this->CreateHistRes(var, flagDown);
  } else if(dirResFlag == 0) {
    hCenter = this->CreateHist(var);
    hUpper  = this->CreateHist(var, flagUp);
    hLower  = this->CreateHist(var, flagDown);
  } else if(dirResFlag > 0){
    hCenter = this->CreateHistDir(var);
    hUpper  = this->CreateHistDir(var, flagUp);
    hLower  = this->CreateHistDir(var, flagDown);
  }

  if(!hCenter) { // || !hUpper || !hLower){
    delete hCenter;
    delete hUpper;
    delete hLower;
    return NULL;
  }
  GFHistArray *result = new GFHistArray(3);
  result->AddAt(hUpper, 0);
  result->AddAt(hCenter, 1);
  result->AddAt(hLower, 2);
  if (hUpper) hUpper->SetLineStyle(4);
  if (hLower) hLower->SetLineStyle(3);

  return result;
}

//_________________________________________________________
 Double_t GFDstarVirtualTheory::TotalCrossSec(ECrossSecType dstarJet,
					     const char *note, Int_t dirResFlag)
{
  // total cross section of type 'dstarJet' for appendix 'note'
  // direct (>0), resolved (<0) or dir.+res. (==0,default)

  // FIXME: overwrite in Pythia to avoid print of eta/EtaDs/EtaDJet
  const char *var  = this->GetVarNameForTotal(dstarJet);
  TH1 *h = NULL;
  if(dirResFlag < 0){
    h = this->CreateHistRes(var, note);
  } else if(dirResFlag == 0) {
    h = this->CreateHist(var, note);
  } else if(dirResFlag > 0){
    h = this->CreateHistDir(var, note);
  }

  if (!h) {
    this->Error("TotalCrossSec", "missing hist for %s", var);
    return 0.;
  } else {
    const Double_t result = h->Integral("width");
    delete h;
    return result;
  }
}

//_________________________________________________________
 TArrayD GFDstarVirtualTheory::TotalCrossSecInt(ECrossSecType dstarJet,Int_t dirResFlag,
					       const TObjArray& upDowns)
{
  // returning array of length 3: total cross section, upper and lower bound
  // 
  // dirResFlag < 0: resolved only
  //            > 0: direct only
  //            = 0: sum
  GFHistArray etaHists;
  const char *var  = this->GetVarNameForTotal(dstarJet);

  if(dirResFlag < 0){
    for(Int_t i = 0; i < upDowns.GetEntriesFast(); ++i){
      etaHists.Add(this->CreateHistRes(var, upDowns[i]->GetName()));
    }
  } else if(dirResFlag == 0) {
    for(Int_t i = 0; i < upDowns.GetEntriesFast(); ++i){
      etaHists.Add(this->CreateHist(var, upDowns[i]->GetName()));
    }
  } else if(dirResFlag > 0){
    for(Int_t i = 0; i < upDowns.GetEntriesFast(); ++i){
      etaHists.Add(this->CreateHistDir(var, upDowns[i]->GetName()));
    }
  }// else only if dirResFlag is nan (is that possible for int?)

  TArrayD xSecs(etaHists.GetEntriesFast());
  for(Int_t i = 0; i < xSecs.GetSize(); ++i){
    if(etaHists[i]) xSecs[i] = etaHists[i]->Integral("width");
  }
  etaHists.Delete();

  TArrayI ind(xSecs.GetSize());
  TMath::Sort(xSecs.GetSize(), xSecs.GetArray(), ind.GetArray());

  TArrayD result(3);
  result[0] = xSecs[0];
  result[1] = xSecs[ind[0]];
  result[2] = xSecs[ind[ind.GetSize()-1]];

  cout << "---------------------------------------------------------------n"
       << this->GetName() << ", total X-sec " << this->GetHeader(dstarJet) << " = " 
       << result[0] << " + " << result[1]-result[0] 
       << " - " << result[0]-result[2] << "n" 
       << "---------------------------------------------------------------n"
       << "max from " << upDowns[ind[0]]->GetName() 
       << ", min from " << upDowns[ind[ind.GetSize()-1]]->GetName() << "n"
       << "---------------------------------------------------------------n"
       << endl;

  if(dstarJet != kInclDoubleDiff && dstarJet != kInclusive && this->ShouldApplyHadCor()){
    cout << "This is without hadronisation correction!n"
	 << "---------------------------------------------------------------n"
	 << endl;
  }

  
  return result; 
}

//_________________________________________________________
 TH1* GFDstarVirtualTheory::CreateHistBins(const char *variable) const
{
  TString var(variable);
  TString var2;
  TArrayD xBins;
  if(var == "pt"){
    Double_t help[] = {2.0, 2.5, 3.0, 3.5, 4.25, 5.0, 6.0, 8.5, 12.};
    xBins.Set(sizeof(help)/sizeof(help[0]), help);
  } else if(var == "eta"){
    Double_t help[] = {-1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5};
    xBins.Set(sizeof(help)/sizeof(help[0]), help);
  } else if(var == "wGammaP"){
    Double_t help[] = {172., 192., 212., 232., 256.};
    xBins.Set(sizeof(help)/sizeof(help[0]), help);
  } else if(var == "zDs"){
//     this->Info("CreateHistBins", "bin is 0.46 instead 0.45"); 
//     Double_t help[] = {0., 0.1, 0.2, 0.3, 0.46, 0.6, 0.8};
    Double_t help[] = {0., 0.1, 0.2, 0.3, 0.45, 0.6, 0.8};
    xBins.Set(sizeof(help)/sizeof(help[0]), help);
  } else if(var == "PtJet" || var == "PtDJet" || var == "PtOthJet") {
    Double_t help[] = {3., 5., 7.5, 11., 16.};
    xBins.Set(sizeof(help)/sizeof(help[0]), help);
  } else if(var == "EtaJet" || var == "EtaDJet" || var == "EtaOthJet"){
    Double_t help[] = {-1.5, -1.0, -0.5, 0.1, 0.8, 1.5};
    xBins.Set(sizeof(help)/sizeof(help[0]), help);
  } else if(var == "PtDs"){ 
    Double_t help[] = {2.0, 3.5, 5.0, 7.5, 11.};
    xBins.Set(sizeof(help)/sizeof(help[0]), help);
  } else if(var == "EtaDs"){  
    Double_t help[] = {-1.5, -1., -0.5, 0.1, 0.8, 1.5};
    xBins.Set(sizeof(help)/sizeof(help[0]), help);
  } else if(var == "Deta" || var == "DetaDijet") {   
    Double_t help[] = {-2.8, -1.9, -1., 0., 1., 1.9, 2.8};
    xBins.Set(sizeof(help)/sizeof(help[0]), help);
  } else if(var == "Dphi" || var == "DphiDijet") {   
    Double_t help[] = {0., 86., 114., 138., 154., 170., 180.};
    xBins.Set(sizeof(help)/sizeof(help[0]), help);
  } else if(var == "DphiNlo" || var == "DphiDijetNlo") {   
     Double_t help[] = {0., 86., 114., 138., 154., 180.};// unify last two bins:
    xBins.Set(sizeof(help)/sizeof(help[0]), help);
  } else if(var == "MDsJet" || var == "MDijet"){ 
    Double_t help[] = {5., 8.5, 13., 18., 28.};
    xBins.Set(sizeof(help)/sizeof(help[0]), help);
  } else if(var == "PtDsJet" || var == "PtDijet"){
    Double_t help[] = {0., 1.8, 3.5, 6., 9.};
    xBins.Set(sizeof(help)/sizeof(help[0]), help);
  } else if(var == "xGam" || var == "xGamDijet"){ 
    Double_t help[] = {0.0, 0.3, 0.6, 0.8, 1.};
    xBins.Set(sizeof(help)/sizeof(help[0]), help);
//   } else if(var == "xGam"){
//     Double_t help[] = {0.15, 0.35, 0.625, 0.825, 1.};
//     xBins.Set(sizeof(help)/sizeof(help[0]), help);
//   } else if(var == "xGamDijet"){ 
//     Double_t help[] = {0.2, 0.425, 0.675, 0.85, 1.};
//     xBins.Set(sizeof(help)/sizeof(help[0]), help);
  } else if(var.BeginsWith("pteta")){
    Double_t help[] = {2.0, 3., 4.25, 7.};
    xBins.Set(sizeof(help)/sizeof(help[0]), help);
    var2 = var(2, var.Length()-2);
    var = "pt";
  } else if(var.BeginsWith("etapt")){
    Double_t help[] = {-1.5, -0.9, 0.2, 1.5};
    xBins.Set(sizeof(help)/sizeof(help[0]), help);
    var2 = var(3, var.Length()-3);
    var = "eta";
  }

  if(xBins.GetSize() > 2){
    TString title(Form("%s %s;%s", GFAxisVariables::GetAxisLabel(var), var2.Data(),
		       GFAxisVariables::GetAxisLabel(var)));
    TH1 *h = new TH1F(variable, title, xBins.GetSize()-1, xBins.GetArray());
    h->SetDirectory(NULL);
    return h;
  } else {
    this->Error("CreateHistBins", "Variable %s not found or less than 1 bin", variable);
    return NULL;
  }
}


//_________________________________________________________
 Bool_t GFDstarVirtualTheory::IsLog(const char *variable) const
{
  TString var(variable);

//   return var.Contains("pt", TString::kIgnoreCase) && var != "PtDsJet"
//     && !var.BeginsWith("etapt", TString::kIgnoreCase)
//     && !var.BeginsWith("zDspt", TString::kIgnoreCase)
//     && !var.BeginsWith("wGammaPpt", TString::kIgnoreCase)
  return (var.BeginsWith("pt", TString::kIgnoreCase) && 
	  var != "PtDijet" && var != "PtDsJet");
}

//_________________________________________________________
 TString GFDstarVirtualTheory::YAxisTitleBinWidth(const char *yTitle,
					    const char *xAxisTitle) const
{
  TString tit(Form("d%s/d%s", yTitle, xAxisTitle));
  tit.ReplaceAll(" [nb]", NULL);
  Ssiz_t pos = tit.First('[');
  if(pos != kNPOS){
    tit.Insert(pos+1, "nb/");
  } else {
    tit += " [nb]";
  }

  return tit;
}


//_________________________________________________________
 TObjArray* GFDstarVirtualTheory::CreateVarArray(const char *var) const
{
  // creates array filled with TObjString containing vars (array is owner!)
  // 1D (inclusive D*), 1Ddd (inclusive D* double diff), 2D (D*+jet),
  // dijet (all D* dijet)
  // all (1D+1Ddd+2D+dijet) or a specific variable
  // or any comma separated list thereof

#if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,8)
  TObjArray *vars = TString(var).Tokenize(",");
#else
  TObjArray *vars = new TObjArray;
  TString tmpStr(var);
  if (tmpStr.Index(",") == kNPOS) {
    vars->Add(new TObjString(var));
  } else {
    this->Warning("CreateVarArray", "ignore non-first in %s (old ROOT...)", var);
    tmpStr = tmpStr(0, tmpStr.Index(",") + 1);
    vars->Add(new TObjString(tmpStr));
  }
#endif

  if (vars->Contains("1D") || vars->Contains("all")) {
    while (TObject *o = vars->Remove(vars->FindObject("1D"))) delete o;
    vars->Add(new TObjString("eta"));
    vars->Add(new TObjString("pt"));
    vars->Add(new TObjString("wGammaP"));
    vars->Add(new TObjString("zDs"));
  } 
  if (vars->Contains("1Ddd") || vars->Contains("all")) {
    while (TObject *o = vars->Remove(vars->FindObject("1Ddd"))) delete o;
    vars->Add(new TObjString("etapt11"));
    vars->Add(new TObjString("etapt22"));
    vars->Add(new TObjString("etapt33"));
    //       "pteta",  "etapt", "zDseta", "etazDs", "zDspt", "zDspt"
  } 
  if (vars->Contains("2D") || vars->Contains("all")) {
    while (TObject *o = vars->Remove(vars->FindObject("2D"))) delete o;
    vars->Add(new TObjString("EtaDs"));
    vars->Add(new TObjString("PtDs"));
    vars->Add(new TObjString("EtaJet"));
    vars->Add(new TObjString("PtJet"));
    vars->Add(new TObjString("Deta"));
    vars->Add(new TObjString("Dphi"));
    vars->Add(new TObjString("DphiNlo"));
    vars->Add(new TObjString("PtDsJet"));
    vars->Add(new TObjString("MDsJet"));
    vars->Add(new TObjString("xGam"));
  }
  if (vars->Contains("dijet") || vars->Contains("all")) {
    while (TObject *o = vars->Remove(vars->FindObject("dijet"))) delete o;
    vars->Add(new TObjString("EtaDJet"));
    vars->Add(new TObjString("EtaOthJet"));
    vars->Add(new TObjString("PtDJet"));
    vars->Add(new TObjString("PtOthJet"));
    vars->Add(new TObjString("DetaDijet"));
    vars->Add(new TObjString("DphiDijet"));
    vars->Add(new TObjString("DphiDijetNlo"));
    vars->Add(new TObjString("PtDijet"));
    vars->Add(new TObjString("MDijet"));
    vars->Add(new TObjString("xGamDijet"));
  }
  
  while (TObject *o = vars->Remove(vars->FindObject("all"))) delete o;
  vars->Compress();
  return vars;
}

////////////////////////////////////////////////////////////
GFDstarVirtualTheory::ECrossSecType GFDstarVirtualTheory::IsDsJetOrDs(const char *var) const
{
  // 0 if 1D variable, 1 if D*-jet, 2 if D* double diff, 3 if D* with dijet
  // -1 if not known

  TString thisVar(var);
  ECrossSecType flag = kUnknown;

  TObjArray *arr = this->CreateVarArray("1D");
  if (arr->Contains(var)) {
    flag = kInclusive;
  } else {
    delete arr;
    arr = this->CreateVarArray("2D");
    if (arr->Contains(var)) flag = kDstarJet;
    else {
      delete arr;
      arr = this->CreateVarArray("1Ddd");
      if (arr->Contains(var)) flag = kInclDoubleDiff;
      else {
	delete arr;
	arr = this->CreateVarArray("dijet");
	if (arr->Contains(var)) flag = kDiJet;
      }
    }
  }

  delete arr;
  return flag;
}

//____________________________________________________
 const char* GFDstarVirtualTheory::GetVarNameForTotal(GFDstarVirtualTheory::ECrossSecType type) const //, const char **titlePtr) const
{
  // return name of distribution to be integrated to total integrated
  // cross sections for incl. D*, D*+o. jet, D* dijet

  switch (type) {
  case kDiJet:
    return "EtaDJet";
  case kDstarJet:
    return "EtaDs";
  case kInclDoubleDiff:
  case kInclusive:
    return "eta";
  case kUnknown:
  default:
    this->Warning("GetVarNameForTotal", "unknown type %d", type);
    return "";
  }

}

//____________________________________________________
 const char* GFDstarVirtualTheory::GetHeader(GFDstarVirtualTheory::ECrossSecType type) const
{
  // header according to cross section type: 
  // 'inclusive D*', 'D* + other jet',...

  switch (type) {
  case kDiJet:
    return "D* tagged dijets";
  case kDstarJet:
    return "D* + other jet";
  case kInclDoubleDiff:
  case kInclusive:
    return "inclusive D*";
  case kUnknown:
  default:
    this->Warning("GetHeader", "unknown type %d", type);
    return "D* + ? (unknown)";
  }

}

//____________________________________________________
 Bool_t GFDstarVirtualTheory::RemoveNb(TString &str) const
{
  // remove [nb] or replace [nb/<?>] by [1/<?>]
  // (not very general: await them around the end...)
  // true if one of the patterns is found and str.is changed
  
  if (!str.Contains("[nb")) return kFALSE;
  const Ssiz_t posBracC = str.Last(']');
  if (posBracC == kNPOS) return kFALSE;

  const Ssiz_t posBracO = str.Index("[nb");
  if (posBracO >= posBracC) return kFALSE;

  const Ssiz_t posBy = str.Last('/');
  if (posBy == kNPOS || posBy <= posBracO) str.Remove(posBracO, posBracC - posBracO +1);
  else str.Replace(posBracO+1, posBy-1 - posBracO, "1");

  return kTRUE;
}