PairAnalysisHistos.cxx

///////////////////////////////////////////////////////////////////////////
//
// Generic Histogram container with support for groups and filling of groups
// by passing a vector of data
//
// Authors:
//   Julian Book   <Julian.Book@cern.ch>
/*

  TOADD: reserved words, MC signals, hist classes, MC weighting, plotting!

  Histograms such as THxF,D, TProfile,2D,3D and n-dimensonal objects such as
  (THn, THnSparse) can be added via the various template functions:
  - AddHistogram
  - AddProfile + different error option (Mean,Rms,.. see: TProfile::BuildOptions(..))
  - AddSparse
  In addition histograms can be filled with weights.

  Different kind of binnings can be passed (linear, arbitrary, logarithmic)
  either via the PairAnalysisHelper functionalities:
  - PairAnalysisHelper::MakeLinBinning(nbins,low,high)
  - PairAnalysisHelper::MakeLogBinning(nbins,low,high)
  - PairAnalysisHelper::MakeArbitraryBinning("1.,4.,10.,..")

  The 'name' and 'title;titleX;titleY;...' of the objects and all axis labels,
  units and names and object names are set in a consistent way and
  it is ensured that they are unique.

  Variables are added via the enumerator of PairAnalysisVarManager. Any
  computation of variables (max=10) are implemented according to TFormula
  using the TMath functions.

  The arguments and options in the function DrawSame and DrawTaskSame provide various and 
  powerful possibilities to plot and post-process histograms. 
  
  Examples:

  // 1 dimensional
  AddHistogram("Event", PairAnalysisHelper::MakeLinBinning(200,-1.,1.),  PairAnalysisVarManager::kXvPrim);

  AddHistogram("Event", PairAnalysisHelper::MakeLinBinning(200,0.0,10.0), "VtxChi/VtxNDF");          // using formulas

  AddHistogram("Track", PairAnalysisHelper::MakeLogBinning(400,0,4.),  PairAnalysisVarManager::kPt); // logarithmic

  AddHistogram("Hit.TRD",PairAnalysisHelper::MakeLinBinning(105,-1.e+2,5.),  "ElossdEdx-ElossMC");   // use MC truth

  // 2 dimensional
  AddHistogram("Track", PairAnalysisHelper::MakeLinBinning(100,-0.5,3.), PairAnalysisVarManager::kY,
                        PairAnalysisHelper::MakeLinBinning(125,0,5.), PairAnalysisVarManager::kPt);

  AddProfile(  "Track", PairAnalysisHelper::MakeLogBinning(50,0.05,10.), PairAnalysisVarManager::kP,
                        PairAnalysisVarManager::kTRDHits, "I"); // Y-average

  ...

*/
//                                                                       //
///////////////////////////////////////////////////////////////////////////

#include <typeinfo>

#include <TAxis.h>
#include <TCanvas.h>
#include <TCollection.h>
#include <TError.h>
#include <TFile.h>
#include <TFormula.h>
#include <TGaxis.h>
#include <TH1.h>
#include <TH1D.h>
#include <TH1F.h>
#include <TH2.h>
#include <TH3.h>
#include <THStack.h>
#include <THashList.h>
#include <THn.h>
#include <THnBase.h>
#include <THnSparse.h>
#include <TKey.h>
#include <TLatex.h>
#include <TLegend.h>
#include <TLegendEntry.h>
#include <TMath.h>
#include <TObjArray.h>
#include <TObjString.h>
#include <TProfile.h>
#include <TProfile2D.h>
#include <TProfile3D.h>
#include <TROOT.h>
#include <TString.h>
#include <TVectorD.h>
#include <TVirtualPS.h>

#include "PairAnalysis.h"
#include "PairAnalysisHelper.h"
#include "PairAnalysisHistos.h"
#include "PairAnalysisMetaData.h"
#include "PairAnalysisSignalMC.h"
#include "PairAnalysisStyler.h"
#include "PairAnalysisVarManager.h"

ClassImp(PairAnalysisHistos)


  PairAnalysisHistos::PairAnalysisHistos()
  :  //   TCollection(),
  TNamed("PairAnalysisHistos", "PairAnalysis Histogram Container")
  , fHistoList()
  , fMetaData(0x0)
  , fList(0x0)
  , fUsedVars(new TBits(PairAnalysisVarManager::kNMaxValues))
  , fReservedWords(new TString("Hit;Track;Pair"))
  , fPrecision(kFloat) {
  //
  // Default constructor
  //
  fHistoList.SetOwner(kTRUE);
  fHistoList.SetName("PairAnalysis_Histos");
  PairAnalysisVarManager::InitFormulas();
  PairAnalysisStyler::LoadStyle();
}

//_____________________________________________________________________________
PairAnalysisHistos::PairAnalysisHistos(const char* name, const char* title)
  :  //   TCollection(),
  TNamed(name, title)
  , fHistoList()
  , fMetaData(0x0)
  , fList(0x0)
  , fUsedVars(new TBits(PairAnalysisVarManager::kNMaxValues))
  , fReservedWords(new TString("Hit;Track;Pair"))
  , fPrecision(kFloat) {
  //
  // TNamed constructor
  //
  fHistoList.SetOwner(kTRUE);
  fHistoList.SetName(name);
  PairAnalysisVarManager::InitFormulas();
  PairAnalysisStyler::LoadStyle();
}

//_____________________________________________________________________________
PairAnalysisHistos::~PairAnalysisHistos() {
  //
  // Destructor
  //
  fHistoList.Clear();
  if (fUsedVars) delete fUsedVars;
  if (fList) fList->Clear();
  if (fMetaData) delete fMetaData;
  delete fReservedWords;
}

//_____________________________________________________________________________
void PairAnalysisHistos::UserHistogram(const char* histClass,
                                       Int_t ndim,
                                       TObjArray* limits,
                                       UInt_t* vars,
                                       UInt_t valTypeW) {
  //
  // Histogram creation n>3 dimension only with non-linear binning
  //

  Bool_t isOk = kTRUE;
  isOk &= (ndim < 21 && ndim > 3);
  if (!isOk) {
    Warning(
      "UserHistogram",
      "Array sizes should be between 3 and 20. Not adding Histogram to '%s'.",
      histClass);
    return;
  }
  isOk &= (ndim == limits->GetEntriesFast());
  if (!isOk) return;

  // set automatic histo name
  TString name;
  for (Int_t iv = 0; iv < ndim; iv++)
    name += Form("%s_", PairAnalysisVarManager::GetValueName(vars[iv]));
  name.Resize(name.Length() - 1);

  isOk &= IsHistogramOk(histClass, name);

  THnD* hist;
  Int_t bins[ndim];
  if (isOk) {
    // get number of bins
    for (Int_t idim = 0; idim < ndim; idim++) {
      TVectorD* vec = (TVectorD*) limits->At(idim);
      bins[idim]    = vec->GetNrows() - 1;
    }

    hist = new THnD(name.Data(), "", ndim, bins, 0x0, 0x0);

    // set binning
    for (Int_t idim = 0; idim < ndim; idim++) {
      TVectorD* vec = (TVectorD*) limits->At(idim);
      hist->SetBinEdges(idim, vec->GetMatrixArray());
    }

    // store variales in axes
    StoreVariables(hist, vars);
    hist->SetUniqueID(valTypeW);  // store weighting variable

    // store which variables are used
    for (Int_t i = 0; i < ndim; i++)
      fUsedVars->SetBitNumber(vars[i], kTRUE);
    fUsedVars->SetBitNumber(valTypeW, kTRUE);

    Bool_t isReserved = fReservedWords->Contains(histClass);
    if (isReserved)
      UserHistogramReservedWords(histClass, hist);
    else
      UserHistogram(histClass, hist);
  }
}

//_____________________________________________________________________________
void PairAnalysisHistos::AddSparse(const char* histClass,
                                   Int_t ndim,
                                   TObjArray* limits,
                                   UInt_t* vars,
                                   UInt_t valTypeW) {
  //
  // THnSparse creation with non-linear binning
  //

  Bool_t isOk = kTRUE;
  isOk &= (ndim == limits->GetEntriesFast());
  if (!isOk) return;

  // set automatic histo name
  TString name;
  for (Int_t iv = 0; iv < ndim; iv++)
    name += Form("%s_", PairAnalysisVarManager::GetValueName(vars[iv]));
  name.Resize(name.Length() - 1);

  isOk &= IsHistogramOk(histClass, name);

  //  THnSparseD *hist;
  PairAnalysisHn* hist;
  Int_t bins[ndim];
  if (isOk) {
    // get number of bins
    for (Int_t idim = 0; idim < ndim; idim++) {
      TVectorD* vec = (TVectorD*) limits->At(idim);
      bins[idim]    = vec->GetNrows() - 1;
    }

    //    hist=new THnSparseD(name.Data(),"", ndim, bins, 0x0, 0x0);
    hist = new PairAnalysisHn(name.Data(), "", ndim, bins, 0x0, 0x0);

    // set binning
    for (Int_t idim = 0; idim < ndim; idim++) {
      TVectorD* vec = (TVectorD*) limits->At(idim);
      hist->SetBinEdges(idim, vec->GetMatrixArray());
    }

    // store variales in axes
    StoreVariables(hist, vars);
    hist->SetUniqueID(valTypeW);  // store weighting variable

    // store which variables are used
    for (Int_t i = 0; i < ndim; i++)
      fUsedVars->SetBitNumber(vars[i], kTRUE);
    fUsedVars->SetBitNumber(valTypeW, kTRUE);

    Bool_t isReserved = fReservedWords->Contains(histClass);
    if (isReserved)
      UserHistogramReservedWords(histClass, hist);
    else
      UserHistogram(histClass, hist);
  }
}

//_____________________________________________________________________________
void PairAnalysisHistos::AddSparse(const char* histClass,
                                   Int_t ndim,
                                   TObjArray* limits,
                                   TFormula** vars,
                                   UInt_t valTypeW) {
  //
  // THnSparse creation with non-linear binning
  //

  Bool_t isOk = kTRUE;
  isOk &= (ndim == limits->GetEntriesFast());
  if (!isOk) return;

  // set automatic histo name
  TString name;
  for (Int_t iv = 0; iv < ndim; iv++)
    name += Form("%s_", (PairAnalysisHelper::GetFormulaName(vars[iv])).Data());
  name.Resize(name.Length() - 1);
  name.ReplaceAll("f(", "");
  name.ReplaceAll(")", "");
  name.Prepend("f(");
  name.Append(")");

  isOk &= IsHistogramOk(histClass, name);

  //  THnSparseD *hist;
  PairAnalysisHn* hist;
  Int_t bins[ndim];
  if (isOk) {
    // get number of bins
    for (Int_t idim = 0; idim < ndim; idim++) {
      TVectorD* vec = (TVectorD*) limits->At(idim);
      bins[idim]    = vec->GetNrows() - 1;
    }

    //    hist=new THnSparseD(name.Data(),"", ndim, bins, 0x0, 0x0);
    hist = new PairAnalysisHn(name.Data(), "", ndim, bins, 0x0, 0x0);

    // set binning
    for (Int_t idim = 0; idim < ndim; idim++) {
      TVectorD* vec = (TVectorD*) limits->At(idim);
      hist->SetBinEdges(idim, vec->GetMatrixArray());
    }

    // add formulas to list of functions
    for (Int_t idim = 0; idim < ndim; idim++) {
      vars[idim]->SetName(Form("axis%dFormula", idim));
      hist->GetListOfFunctions()->Add(vars[idim]);
    }

    // store variales in axes
    //    StoreVariables(hist, vars);
    hist->SetUniqueID(valTypeW);  // store weighting variable

    // adapt the name and title of the histogram in case they are empty
    AdaptNameTitle(hist, histClass);

    // store which variables are used
    //for(Int_t i=0; i<ndim; i++)   fUsedVars->SetBitNumber(vars[i],kTRUE);
    fUsedVars->SetBitNumber(valTypeW, kTRUE);

    Bool_t isReserved = fReservedWords->Contains(histClass);
    if (isReserved)
      UserHistogramReservedWords(histClass, hist);
    else
      UserHistogram(histClass, hist);
  }
}

//_____________________________________________________________________________
TString PairAnalysisHistos::UserHistogram(const char* histClass,
                                          TObject* hist) {
  //
  // Add any type of user histogram
  //

  //special case for the calss Pair. where histograms will be created for all pair classes
  Bool_t isReserved = fReservedWords->Contains(histClass);
  if (isReserved) {
    UserHistogramReservedWords(histClass, hist);
    return hist->GetName();
  }

  // if (!IsHistogramOk(histClass,hist->GetName())) return hist->GetName();
  // THashList *classTable=(THashList*)fHistoList.FindObject(histClass);
  //  //  hist->SetDirectory(0);

  // store variables axis
  UInt_t valType[20] = {0};
  // extract variables from axis
  FillVarArray(hist, valType);
  // TODO: implement THn/PairAnalysisHn
  // change to mc truth variables if available
  TString hclass = histClass;
  if (hclass.Contains("MCtruth") && !(hist->IsA() == PairAnalysisHn::Class())) {
    for (Int_t i = 0; i < 2; i++) {
      /// protection for changes of variable enum (in the postprocessing) --> use axistname indentification
      //      UInt_t valTypeFromTitle = 0;
      if (!i)
        valType[i] = PairAnalysisVarManager::GetValueType(
          ((TH1*) hist)->GetXaxis()->GetName());
      else
        valType[i] = PairAnalysisVarManager::GetValueType(
          ((TH1*) hist)->GetYaxis()->GetName());
      //Printf("SWITCH TO MC: before: %d %s, (via axis name %d) ---->",valType[i],PairAnalysisVarManager::GetValueName(valType[i]),valTypeFromTitle);
      valType[i] = PairAnalysisVarManager::GetValueTypeMC(valType[i]);
      // if theres no corresponding MCtruth variable, skip adding this histogram
      //      if(valType[i] < PairAnalysisVarManager::kNMaxValues && valType[i]>0) return hist->GetName();
      if (valType[i] < PairAnalysisVarManager::kPairMax && valType[i] > 0)
        return hist->GetName();
      // request filling of mc variable
      fUsedVars->SetBitNumber(valType[i], kTRUE);
      //      Printf("after: %d %s",valType[i],PairAnalysisVarManager::GetValueName(valType[i]));
    }
    StoreVariables(hist, valType);
    hist->SetUniqueID(valType[19]);  // store weighting variable
    // check for formulas
    if (hist->InheritsFrom(TH1::Class())) {
      TIter next(((TH1*) hist)->GetListOfFunctions());
      TFormula* f = 0;
      while ((f = dynamic_cast<TFormula*>(next()))) {
        for (Int_t i = 0; i < f->GetNpar(); i++) {
          Int_t parMC =
            PairAnalysisVarManager::GetValueTypeMC(f->GetParameter(i));
          // if theres none corresponding MCtruth variable, skip adding this histogram
          if (parMC < PairAnalysisVarManager::kNMaxValues)
            return hist->GetName();
          f->SetParameter(i, parMC);
          f->SetParName(i, PairAnalysisVarManager::GetValueName(parMC));
          fUsedVars->SetBitNumber(parMC, kTRUE);
        }
      }
      // change histogram key according to mctruth information
      AdaptNameTitle((TH1*) hist, histClass);
    }
  } else {
    StoreVariables(hist, valType);
    hist->SetUniqueID(valType[19]);  // store weighting variable
  }

  // add histogram to class
  if (!IsHistogramOk(histClass, hist->GetName())) return hist->GetName();
  THashList* classTable = (THashList*) fHistoList.FindObject(histClass);
  //  hist->SetDirectory(0);
  if (classTable) classTable->Add(hist);
  return hist->GetName();
}

//_____________________________________________________________________________
TH1* PairAnalysisHistos::GetTHist(const char* histClass,
                                  const char* name,
                                  const char* title,
                                  const TVectorD* const binsX,
                                  const TVectorD* const binsY,
                                  const TVectorD* const binsZ) {
  //
  // retrieve n-dimensional Hist depending on arguments
  //
  Bool_t isOk = kTRUE;
  isOk &= IsHistogramOk(histClass, name);
  isOk &= (binsX != 0x0);
  if (!isOk) return 0x0;
  switch (fPrecision) {
    case kFloat:
      if (!binsY)
        return (new TH1F(
          name, title, binsX->GetNrows() - 1, binsX->GetMatrixArray()));
      else if (!binsZ)
        return (new TH2F(name,
                         title,
                         binsX->GetNrows() - 1,
                         binsX->GetMatrixArray(),
                         binsY->GetNrows() - 1,
                         binsY->GetMatrixArray()));
      else
        return (new TH3F(name,
                         title,
                         binsX->GetNrows() - 1,
                         binsX->GetMatrixArray(),
                         binsY->GetNrows() - 1,
                         binsY->GetMatrixArray(),
                         binsZ->GetNrows() - 1,
                         binsZ->GetMatrixArray()));
      break;
    case kDouble:
      if (!binsY)
        return (new TH1D(
          name, title, binsX->GetNrows() - 1, binsX->GetMatrixArray()));
      else if (!binsZ)
        return (new TH2D(name,
                         title,
                         binsX->GetNrows() - 1,
                         binsX->GetMatrixArray(),
                         binsY->GetNrows() - 1,
                         binsY->GetMatrixArray()));
      else
        return (new TH3D(name,
                         title,
                         binsX->GetNrows() - 1,
                         binsX->GetMatrixArray(),
                         binsY->GetNrows() - 1,
                         binsY->GetMatrixArray(),
                         binsZ->GetNrows() - 1,
                         binsZ->GetMatrixArray()));
      break;
    default: return 0x0; break;
  }
}

//_____________________________________________________________________________
TH1* PairAnalysisHistos::GetTProf(const char* histClass,
                                  const char* name,
                                  const char* title,
                                  const TVectorD* const binsX,
                                  const TVectorD* const binsY,
                                  const TVectorD* const binsZ,
                                  TString option) {
  //
  // retrieve n-dimensional profile histogram with error options depending on arguments
  //
  Bool_t isOk = kTRUE;
  isOk &= IsHistogramOk(histClass, name);
  isOk &= (binsX != 0x0);
  if (!isOk) return 0x0;
  // parse error option
  TString opt   = "";
  Double_t pmin = 0., pmax = 0.;
  if (!option.IsNull()) {
    TObjArray* arr = option.Tokenize(";");
    arr->SetOwner();
    opt = ((TObjString*) arr->At(0))->GetString();
    if (arr->GetEntriesFast() > 1)
      pmin = (((TObjString*) arr->At(1))->GetString()).Atof();
    if (arr->GetEntriesFast() > 2)
      pmax = (((TObjString*) arr->At(2))->GetString()).Atof();
    delete arr;
  }
  // build profile with error options and return it
  TH1* prof = 0x0;
  if (!binsY)
    return (new TProfile(name,
                         title,
                         binsX->GetNrows() - 1,
                         binsX->GetMatrixArray(),
                         pmin,
                         pmax,
                         opt.Data()));
  else if (!binsZ) {
    prof = new TProfile2D(name,
                          title,
                          binsX->GetNrows() - 1,
                          binsX->GetMatrixArray(),
                          binsY->GetNrows() - 1,
                          binsY->GetMatrixArray());
    ((TProfile2D*) prof)->BuildOptions(pmin, pmax, opt.Data());
    return prof;
  } else {
    prof = new TProfile3D(name,
                          title,
                          binsX->GetNrows() - 1,
                          binsX->GetMatrixArray(),
                          binsY->GetNrows() - 1,
                          binsY->GetMatrixArray(),
                          binsZ->GetNrows() - 1,
                          binsZ->GetMatrixArray());
    ((TProfile3D*) prof)->BuildOptions(pmin, pmax, opt.Data());
    return prof;
  }
}

//_____________________________________________________________________________
TFormula* PairAnalysisHistos::GetFormula(const char* name,
                                         const char* formula) {
  //
  // build a TFormula object
  //
  TFormula* form = new TFormula(name, formula);
  // compile function
  if (form->Compile()) return 0x0;
  //set parameter/variable identifier
  for (Int_t i = 0; i < form->GetNpar(); i++) {
    form->SetParName(
      i, PairAnalysisVarManager::GetValueName(form->GetParameter(i)));
    fUsedVars->SetBitNumber((Int_t) form->GetParameter(i), kTRUE);
  }
  return form;
}

//_____________________________________________________________________________
void PairAnalysisHistos::AddClass(const char* histClass) {
  //
  // Add a class of histograms
  // Several classes can be added by separating them by a ';' e.g. 'class1;class2;class3'
  //
  TString hists(histClass);
  TObjArray* arr = hists.Tokenize(";");
  TIter next(arr);
  TObject* o = 0;
  while ((o = next())) {
    if (fHistoList.FindObject(o->GetName())) {
      Warning(
        "AddClass", "Cannot create class '%s' it already exists.", histClass);
      continue;
    }
    if (fReservedWords->Contains(o->GetName())) {
      Error("AddClass", "Pair is a reserved word, please use another name");
      continue;
    }
    THashList* table = new THashList;
    table->SetOwner(kTRUE);
    table->SetName(o->GetName());
    fHistoList.Add(table);
  }
  delete arr;
}

//_____________________________________________________________________________
void PairAnalysisHistos::FillClass(TString histClass, const Double_t* values) {
  //
  // Fill class 'histClass' (by name)
  //
  THashList* classTable = (THashList*) fHistoList.FindObject(histClass.Data());
  if (!classTable) {
    //    Warning("FillClass","Cannot fill class '%s' its not defined.",histClass.Data());
    return;
  }

  TIter nextHist(classTable);
  TObject* obj = 0;
  while ((obj = (TObject*) nextHist()))
    FillValues(obj, values);

  return;
}

//_____________________________________________________________________________
void PairAnalysisHistos::UserHistogramReservedWords(const char* histClass,
                                                    const TObject* hist) {
  //
  // Creation of histogram for all pair or track types
  //
  TString title(hist->GetTitle());

  TIter nextClass(&fHistoList);
  THashList* l = 0;
  while ((l = static_cast<THashList*>(nextClass()))) {
    TString name(l->GetName());
    if (name.Contains(histClass)) {
      TObject* h = hist->Clone();
      // Tobject has no function SetDirectory, didn't we need this???
      //      h->SetDirectory(0);
      if (h->InheritsFrom(TH1::Class()))
        ((TH1*) h)->SetTitle(Form("%s %s", title.Data(), l->GetName()));
      else
        ((THnBase*) h)->SetTitle(Form("%s %s", title.Data(), l->GetName()));

      UserHistogram(l->GetName(), h);
    }
  }
  delete hist;

  return;
}

//_____________________________________________________________________________
void PairAnalysisHistos::DumpToFile(const char* file) {
  //
  // Dump the histogram list to a newly created root file
  //
  TFile f(file, "recreate");
  fHistoList.Write(fHistoList.GetName(), TObject::kSingleKey);
  f.Close();
}

//_____________________________________________________________________________
TObject* PairAnalysisHistos::GetHist(const char* histClass,
                                     const char* name) const {
  //
  // return object 'name' in 'histClass'
  //
  THashList* classTable = (THashList*) fHistoList.FindObject(histClass);
  if (!classTable) return 0x0;
  return classTable->FindObject(name);
}

//_____________________________________________________________________________
TH1* PairAnalysisHistos::GetHistogram(const char* histClass,
                                      const char* name) const {
  //
  // return histogram 'name' in 'histClass'
  //
  return ((TH1*) GetHist(histClass, name));
}

//_____________________________________________________________________________
TObject* PairAnalysisHistos::GetHist(const char* cutClass,
                                     const char* histClass,
                                     const char* name) const {
  //
  // return object from list of list of histograms
  // this function is thought for retrieving histograms if a list of PairAnalysisHistos is set
  //

  if (!fList) return 0x0;
  THashList* h = dynamic_cast<THashList*>(fList->FindObject(cutClass));
  if (!h) return 0x0;
  THashList* classTable = dynamic_cast<THashList*>(h->FindObject(histClass));
  if (!classTable) return 0x0;
  return classTable->FindObject(name);
}

//_____________________________________________________________________________
TH1* PairAnalysisHistos::GetHistogram(const char* cutClass,
                                      const char* histClass,
                                      const char* name) const {
  //
  // return histogram from list of list of histograms
  // this function is thought for retrieving histograms if a list of PairAnalysisHistos is set
  //
  return ((TH1*) GetHist(cutClass, histClass, name));
}

//_____________________________________________________________________________
void PairAnalysisHistos::Draw(const Option_t* option) {
  //
  // Draw histograms
  //

  TString drawStr(option);
  TObjArray* arr = drawStr.Tokenize(";");
  arr->SetOwner();
  TIter nextOpt(arr);

  TString drawClasses;
  TObjString* ostr = 0x0;

  TString currentOpt;
  TString testOpt;
  while ((ostr = (TObjString*) nextOpt())) {
    currentOpt = ostr->GetString();
    currentOpt.Remove(TString::kBoth, '\t');
    currentOpt.Remove(TString::kBoth, ' ');

    testOpt = "classes=";
    if (currentOpt.Contains(testOpt.Data())) {
      drawClasses = currentOpt(testOpt.Length(), currentOpt.Length());
    }
  }

  delete arr;
  drawStr.ToLower();
  //optionsfList
  //   Bool_t same=drawOpt.Contains("same"); //FIXME not yet implemented

  TCanvas* c = 0x0;
  if (gVirtualPS) {
    if (!gPad) {
      Error("Draw",
            "When writing to a file you have to create a canvas before opening "
            "the file!!!");
      return;
    }
    c = gPad->GetCanvas();
    c->cd();
    //     c=new TCanvas;
  }

  TIter nextClass(&fHistoList);
  THashList* classTable = 0;
  //   Bool_t first=kTRUE;
  while ((classTable = (THashList*) nextClass())) {
    //test classes option
    if (!drawClasses.IsNull() && !drawClasses.Contains(classTable->GetName()))
      continue;
    //optimised division
    Int_t nPads = classTable->GetEntries();
    Int_t nCols = (Int_t) TMath::Ceil(TMath::Sqrt(nPads));
    Int_t nRows = (Int_t) TMath::Ceil((Double_t) nPads / (Double_t) nCols);

    //create canvas
    if (!gVirtualPS) {
      TString canvasName;
      canvasName.Form("c%s_%s", GetName(), classTable->GetName());
      c = (TCanvas*) gROOT->FindObject(canvasName.Data());
      if (!c)
        c = new TCanvas(canvasName.Data(),
                        Form("%s: %s", GetName(), classTable->GetName()));
      c->Clear();
    } else {
      //       if (first){
      //         first=kFALSE;
      //         if (nPads>1) gVirtualPS->NewPage();
      //       } else {
      if (nPads > 1) c->Clear();
      //       }
    }
    if (nCols > 1 || nRows > 1) c->Divide(nCols, nRows);

    //loop over histograms and draw them
    TIter nextHist(classTable);
    Int_t iPad = 0;
    TH1* h     = 0;
    while ((h = (TH1*) nextHist())) {
      TString drawOpt;
      if ((h->InheritsFrom(TH2::Class()))) drawOpt = "colz";
      if (nCols > 1 || nRows > 1) c->cd(++iPad);
      if (TMath::Abs(h->GetXaxis()->GetBinWidth(1)
                     - h->GetXaxis()->GetBinWidth(2))
          > 1e-10)
        gPad->SetLogx();
      if (TMath::Abs(h->GetYaxis()->GetBinWidth(1)
                     - h->GetYaxis()->GetBinWidth(2))
          > 1e-10)
        gPad->SetLogy();
      if (TMath::Abs(h->GetZaxis()->GetBinWidth(1)
                     - h->GetZaxis()->GetBinWidth(2))
          > 1e-10)
        gPad->SetLogz();
      TString histOpt = h->GetOption();
      histOpt.ToLower();
      if (histOpt.Contains("logx")) gPad->SetLogx();
      if (histOpt.Contains("logy")) gPad->SetLogy();
      if (histOpt.Contains("logz")) gPad->SetLogz();
      histOpt.ReplaceAll("logx", "");
      histOpt.ReplaceAll("logy", "");
      histOpt.ReplaceAll("logz", "");
      h->Draw(drawOpt.Data());
    }
    if (gVirtualPS) { c->Update(); }
  }
  //   if (gVirtualPS) delete c;
}

//_____________________________________________________________________________
void PairAnalysisHistos::Print(const Option_t* option) const {
  //
  // Print classes and histograms
  //
  TString optString(option);

  if (optString.IsNull()) PrintStructure();
}

//_____________________________________________________________________________
void PairAnalysisHistos::PrintStructure() const {
  //
  // Print classes and histograms in the class to stdout
  //
  if (!fList) {
    TIter nextClass(&fHistoList);
    THashList* classTable = 0;
    while ((classTable = (THashList*) nextClass())) {
      TIter nextHist(classTable);
      TObject* o = 0;
      Printf("+ %s\n", classTable->GetName());
      while ((o = nextHist()))
        Printf("| ->%s\n", o->GetName());
    }
  } else {
    TIter nextCutClass(fList);
    THashList* cutClass = 0x0;
    while ((cutClass = (THashList*) nextCutClass())) {
      TString cla = cutClass->GetName();
      if (cla.Contains("QAcuts")) continue;
      Printf("+ %s\n", cutClass->GetName());
      TIter nextClass(cutClass);
      THashList* classTable = 0;
      while ((classTable = (THashList*) nextClass())) {
        TIter nextHist(classTable);
        TObject* o = 0;
        Printf("|  + %s\n", classTable->GetName());
        while ((o = nextHist()))
          Printf("|  | ->%s\n", o->GetName());
      }
    }
  }
}

//_____________________________________________________________________________
void PairAnalysisHistos::SetHistogramList(THashList& list,
                                          Bool_t setOwner /*=kTRUE*/) {
  //
  // set histogram classes and histograms to this instance. It will take onwnership!
  //
  ResetHistogramList();
  TString name(GetName());
  //  if (name == "PairAnalysisHistos")
  SetName(list.GetName());
  TIter next(&list);
  TObject* o;
  while ((o = next())) {
    fHistoList.Add(o);
  }
  if (setOwner) {
    list.SetOwner(kFALSE);
    fHistoList.SetOwner(kTRUE);
    fHistoList.SetName(list.GetName());
  } else {
    fHistoList.SetOwner(kFALSE);
    fHistoList.SetName(list.GetName());
  }
}

//_____________________________________________________________________________
Bool_t PairAnalysisHistos::SetCutClass(const char* cutClass) {
  //
  // Assign histogram list according to cutClass
  //

  if (!fList) return kFALSE;
  ResetHistogramList();
  THashList* h = dynamic_cast<THashList*>(fList->FindObject(cutClass));
  if (!h) {
    Warning("SetCutClass", "cutClass '%s' not found", cutClass);
    return kFALSE;
  }
  SetHistogramList(*h, kFALSE);
  return kTRUE;
}
//_____________________________________________________________________________
Bool_t PairAnalysisHistos::IsHistogramOk(const char* histClass,
                                         const char* name) {
  //
  // check whether the histogram class exists and the histogram itself does not exist yet
  //
  Bool_t isReserved = fReservedWords->Contains(histClass);
  if (!fHistoList.FindObject(histClass) && !isReserved) {
    Warning("IsHistogramOk",
            "Cannot create histogram. Class '%s' not defined. Please create it "
            "using AddClass before.",
            histClass);
    return kFALSE;
  }
  if (GetHist(histClass, name)) {
    Warning("IsHistogramOk",
            "Cannot create histogram '%s' in class '%s': It already exists!",
            name,
            histClass);
    return kFALSE;
  }
  return kTRUE;
}

// //_____________________________________________________________________________
// TIterator* PairAnalysisHistos::MakeIterator(Bool_t dir) const
// {
//   //
//   //
//   //
//   return new TListIter(&fHistoList, dir);
// }

//_____________________________________________________________________________
void PairAnalysisHistos::ReadFromFile(const char* file,
                                      const char* task,
                                      const char* config) {
  //
  // Read histos from file
  //
  TFile f(file);
  TIter nextKey(f.GetListOfKeys());
  TKey* key = 0;
  while ((key = (TKey*) nextKey())) {
    TString name = key->GetName();
    // check for meta data
    if (name.Contains(Form("PairAnalysisMetaData_%s", task))) {
      fMetaData = new PairAnalysisMetaData();
      fMetaData->SetMetaData(*dynamic_cast<TList*>(f.Get(key->GetName())),
                             kFALSE);
    }
    // check for histos
    if (!name.Contains(Form("PairAnalysisHistos_%s", task))) continue;
    if (!strlen(task) && !name.Contains(task)) continue;
    TObject* o  = f.Get(key->GetName());
    TList* list = dynamic_cast<TList*>(o);
    if (!list)
      continue;
    else
      fList = list;
    THashList* listCfg = dynamic_cast<THashList*>(list->FindObject(config));
    if (!listCfg) continue;
    SetHistogramList(*listCfg);
    break;
  }
  f.Close();
  // load style
  PairAnalysisStyler::LoadStyle();
}

//_____________________________________________________________________________
TObjArray* PairAnalysisHistos::DrawTaskSame(TString histName,
                                            TString opt,
                                            TString histClassDenom,
                                            TString taskDenom) {
  ///
  /// Draw histograms of different tasks into the same canvas
  ///
  /// additional plotting options to DrawSame:
  /// "selcfg":     option to (de-)select certain task/config specified in "histClassDenom",
  ///               delimiters ':;,' are recognized, exclusion done by prepending '!'
  /// "cutstep":    histograms of the cut steps are compared for a single task (see AnalysisFilter::AddHistos)
  ///
  /// if argument "taskDenom" is provided histogram ratios to this task/config are calculated


  TObjArray* selections = histClassDenom.Tokenize(":;,");

  opt.ToLower();
  TString optString(opt);
  Bool_t optGoff    = optString.Contains("goff");
  Bool_t optEff     = optString.Contains("eff");
  Bool_t optCutStep = optString.Contains("cutstep");
  Bool_t optSelCfg  = optString.Contains("selcfg");
  opt.ReplaceAll("selcfg", "");
  fHistoList.SetOwner(kFALSE);

  /// output array
  TObjArray* arr = NULL;
  if (optGoff) Info("DrawTaskSame", "graphics option off, collect an array");
  {
    //    arr=(TObjArray*)gROOT->FindObject(Form("%s",histName.Data()));
    //arr=(TObjArray*)gROOT->FindObject(GetName());
    //    if(arr) { Warning("DrawTaskSame","Clear existing array %s",GetName()); arr->ls(); arr->Clear(); }
    //    else