// Author: Gero Flucke <mailto:flucke@mail.desy.de>
//____________________________________
// GFDstarSysErrors
//   Author:      Gero Flucke
//   Date:        April 5th, 2003
//   last update: $Date: 2006/01/11 11:30:24 $
//   by:          $Author: flucke $
//

#include <iostream>

#include <TROOT.h>
// #include <TFile.h>
#include <TH1.h>
// #include <TH2.h>
#include <TString.h>
#include <TMath.h>
// #include <TAxis.h>
// #include <TObjArray.h>

#include "GFAnalysis/GFDstarSysErrors.h"
#include "GFAnalysis/GFDstarHistsAnalysisData.h"
#include "GFAnalysis/GFDstarHistsAnalysisMc.h"
#include "GFAnalysis/GFDstarXSections.h"
#include "GFAnalysis/GFAxisVariables.h"

#include "GFUtils/GFHistManager.h"
#include "GFUtils/GFHistManip.h"
#include "GFUtils/GFHistArray.h"
#include "GFUtils/GFMath.h"

using namespace std;

ClassImp(GFDstarSysErrors)

 GFDstarSysErrors::GFDstarSysErrors(const char* prefix)
{
  this->Init(prefix);
}

//____________________________________
 GFDstarSysErrors::~GFDstarSysErrors()
{
  delete fHistManager;

  delete fXSections;

  delete fAnalysisDataUp;
  delete fAnalysisDataDown;

  delete fAnalysisMcLeft;
  delete fAnalysisMcRight;
  fHist2 = NULL; // no delete...
  delete fOutFile;
}

//____________________________________
 void GFDstarSysErrors::Init(const char* prefix)
{
  fHistManager  = new GFHistManager;

  fXSections = new GFDstarXSections(prefix);
  fXSections->SetBatch();
//   fXSections = new GFDstarXSections("/data/usr/flucke/dstar/data/oo/histFiles/"
// 				    "weight3DDstarAnalysisHistsAll.root",
// 				    "/data/usr/flucke/dstar/data/oo/histFiles/"
// 				    "weight3DDstarAnalysisHistsDirCharm.root",
// 				    "/data/usr/flucke/dstar/data/oo/histFiles/"
// 				    "weight3DDstarAnalysisHistsResCharm.root",
// 				    NULL);
//   fXSections->GetDstarHistsAnalysisData()->SetBatch(kFALSE);
//   fXSections->GetDstarHistsAnalysisMc()->SetBatch(kFALSE);

  fAnalysisDataUp = NULL;//new GFDstarHistsAnalysisData(prefix, "AllCaloUp");
//   fAnalysisDataUp->SetBatch();
  fAnalysisDataDown = NULL;// = new GFDstarHistsAnalysisData(prefix, "AllCaloDown");
//   fAnalysisDataDown->SetBatch();


  fAnalysisMcLeft = NULL;//new GFDstarHistsAnalysisMcMore("left");
//   fAnalysisMcLeft->SetBatch();
  fAnalysisMcRight = NULL; //new GFDstarHistsAnalysisMcMore("right");
//   fAnalysisMcRight->SetBatch();

  fHist2 = NULL;
  fOutFile = NULL;

  fSysBranch = 0.025; // BRanching dstar->Kpipi
  fSysLumi   = 0.015; 
  fSysL4   = 0.02;//0.04;//0.035;
  fSysTrigEff = 0.04;
  fSysTagTotal83=0.058;
  fSysReflect = 0.015;

  fSysParticleID = 0.02;

//   fSysTrackEff = 0.08;//0.11; cf. mail to Sebastian etc. 29th of June 2005
  fSysTrackEff = 0.06;//cf. 
  fSysFit = 0.03;//5;

  fSysDsJetTrig = 0.04;
//   fSysDsJetScale = 0.02;// cf. log book 16/17.3.04
//   fSysDsJetScale = 0.015;// cf. revised log book entry and thesis
  fSysDsJetScale = 0.02;// tested by scaling jet-pt according to its fraction
                        // is pt, xgamma, dphi dependent: 1-3 % 
}

//____________________________________
 void GFDstarSysErrors::SetTrigger(Int_t trigger)
{
  fXSections->SetTrigger(trigger);
  if(fAnalysisMcRight) fAnalysisMcRight->SetTrigger(trigger);
  if(fAnalysisMcLeft) fAnalysisMcLeft->SetTrigger(trigger);
}

//____________________________________
 void GFDstarSysErrors::SetBatch(Bool_t yesNo)
{
  fHistManager->SetBatch(yesNo);
  fXSections->SetBatch(yesNo);
  if(fAnalysisMcRight) fAnalysisMcRight->SetBatch(yesNo);
  if(fAnalysisMcLeft) fAnalysisMcLeft->SetBatch(yesNo);
}

//____________________________________
 void GFDstarSysErrors::SetOutFile(const char *fileName)
{
  if(fOutFile) delete fOutFile;
  TDirectory *oldDir = gDirectory;
  fOutFile = TFile::Open(fileName, "UPDATE");
  if(oldDir) oldDir->cd();
  else gROOT->cd();
}

//____________________________________
 Double_t GFDstarSysErrors::GetSysErrorTotal()
{
  Double_t total2 = 0.;

  const Double_t tri = this->GetSysErrorTrigEff();
  total2 += (tri*tri);

  const Double_t l4 = this->GetSysErrorL4();
  total2 += (l4*l4);

  const Double_t ad = this->GetSysErrorAcceptData();
  total2 += (ad*ad);

  const Double_t tr = this->GetSysErrorTrackEff();
  total2 += (tr*tr);

  const Double_t id = this->GetSysErrorParticleID();
  total2 += (id*id);

  total2 += (fSysBranch*fSysBranch); //BRanching dstar->Kpipi, correlated!
  total2 += (fSysReflect*fSysReflect); // correlated!

  const Double_t fit = this->GetSysErrorFit();
  total2 += (fit*fit);

  const Double_t mod = this->GetSysErrorModel();
  total2 += (mod*mod);

  total2 += (fSysLumi*fSysLumi); //lumi, correlated

  return TMath::Sqrt(total2);
}

//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorTotal(const char* var, Bool_t norm)
{
  GFHistArray uncertHists;

  const TH1* hModel = this->GetSysErrorModel(var, norm);
  uncertHists.Add(hModel->Clone());
  
  if (!norm) {
    uncertHists.Add(this->GetSysErrorTrigEff(var, hModel));
    uncertHists.Add(this->GetSysErrorL4(var, hModel));
    uncertHists.Add(this->GetSysErrorTrackEff(var, hModel));
    uncertHists.Add(this->GetSysErrorParticleID(var, hModel));
    uncertHists.Add(this->GetSysErrorBranching(var, hModel));  // not binwise?
    uncertHists.Add(this->GetSysErrorReflections(var, hModel));// not binwise?
    uncertHists.Add(this->GetSysErrorFit(var, hModel)); // norm???
    uncertHists.Add(this->GetSysErrorLumi(var, hModel));       // not binwise?
    uncertHists.Add(this->GetSysErrorAcceptData(var, hModel)); // fixme for norm W!
  } 
  if (!norm || TString(var) == "wGammaP") {
    uncertHists.Add(this->GetSysErrorAcceptData(var, hModel)); // fixme for norm W!
  }

  return this->CombineErrorHists(&uncertHists, var);
}

//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorFit(const char* var, const TH1* hBins)
{
  TH1* h = NULL;
  if(hBins){
    h = GFHistManip::CreateAnalog(hBins, Form("sysErrFit%s%d", var, fXSections->GetTrigger()), 
				  kFALSE);
  } else {
    this->Info("GetSysErrorFit", "no hBins given, take GetHist(%s,...)", var);
    h = this->GetHist(var, "sysErrFit");
    if(!h) return NULL;
  }
  h->SetTitle("Systematic uncertainty of N(D*)");
  Int_t nBins = h->GetNbinsX();
//   if(fXSections->GetTrigger() == 83 && TString(var) == "eta"){
//     h->SetBinContent(nBins, 0.1); // last eta bins are bad...
//     h->SetBinContent(nBins-1, 0.1);
//     nBins -= 2;
//   }
  for(Int_t bin = 1; bin <= nBins; ++bin){
    h->SetBinContent(bin, this->GetSysErrorFit());
  }
  this->PrintTable(h, h->GetName());
  return h;
}

//____________________________________
 Double_t GFDstarSysErrors::GetSysErrorFit()
{
  return fSysFit;
}

//____________________________________
 TH1* GFDstarSysErrors::CalcSysErrorFit(const char* var, Int_t fitOrNotFlag)
{
  // normal fit vs fitting WC background and subtracting its integral from signal region

  GFDstarHistsAnalysisData *data = fXSections->GetDstarHistsAnalysisData();
  TH1* histNorm = data->DrawDm(var);
  const Int_t oldFlag = data->GetMultiDmFitter()->GetFitOrNotFlag();
  TH1* histCheck = NULL;
  switch(fitOrNotFlag){
  case -1:
    histCheck = data->DrawDmFitWc(var);
    break;
  case 0: case 1: case 2: case 3: case 4: case 5:
    data->GetMultiDmFitter()->SetFitOrNotFlag(fitOrNotFlag);
    histCheck = data->DrawDm(var);
    break;
  default:
    this->Warning("CalcSysErrorFit", "Don't know fitOrNotFlag = %d, use 5", fitOrNotFlag);
    data->GetMultiDmFitter()->SetFitOrNotFlag(5);
    histCheck = data->DrawDm(var);
  }

  data->GetMultiDmFitter()->SetFitOrNotFlag(oldFlag);

//   histMc->SetTitle("D*: A_{33}*#epsilon_{rec} * #epsilon_{L1}");

  TH1* hSysErr = this->GetRelDeviation(histNorm, histCheck, var, "normal fit", 
				       Form("fit flag %d",fitOrNotFlag), "R");
  return hSysErr;
}

//____________________________________
 Double_t GFDstarSysErrors::CalcSysErrorFit(Int_t fitOrNotFlag)
{
  TArrayD normArr = fXSections->GetDstarHistsAnalysisData()->DrawDm();
  const Int_t oldFlag = fXSections->GetDstarHistsAnalysisData()->GetMultiDmFitter()
    ->GetFitOrNotFlag();
  TArrayD checkArr;
  switch(fitOrNotFlag){
  case -1:
    checkArr = fXSections->GetDstarHistsAnalysisData()->DrawDmFitWc();
    break;
  case 0: case 1: case 2: case 3: case 4: case 5:
    fXSections->GetDstarHistsAnalysisData()->GetMultiDmFitter()->SetFitOrNotFlag(fitOrNotFlag);
    checkArr = fXSections->GetDstarHistsAnalysisData()->DrawDm();
    break;
  default:
    this->Warning("CalcSysErrorFit", "Don't know fitOrNotFlag = %d, use 5", fitOrNotFlag);
    fXSections->GetDstarHistsAnalysisData()->GetMultiDmFitter()->SetFitOrNotFlag(5);
    checkArr = fXSections->GetDstarHistsAnalysisData()->DrawDm();
  }

  fXSections->GetDstarHistsAnalysisData()->GetMultiDmFitter()->SetFitOrNotFlag(oldFlag);

//   histMc->SetTitle("D*: A_{33}*#epsilon_{rec} * #epsilon_{L1}");

  const Double_t res = this->GetRelDeviation(normArr[0], checkArr[0]);
  char tmp[] = "------------------------------------------------n";
  cout << tmp << "Rel. deviation (flag " << fitOrNotFlag << ") of N(D*) is " << res 
       << ", rel. stat. error is " << normArr[1]/normArr[0] << "n"
       << tmp << endl;

  return res;
}


//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorRecEff(const char* var)
{
  // comparing different MC's for eff rec

  TH1* histMc = fXSections->GetDstarHistsAnalysisMc()->CreateRecEff(var);
//   TH1* histMcRef = fAnalysisMcRef->CreateRecEff(var);
  TH1* histMcRef = fXSections->GetDstarHistsAnalysisOtherMc()->CreateRecEff(var);
  return this->GetRelDeviation(histMc, histMcRef, var);
}

//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorTrackEff(const char* var, const TH1* hBins)
{
  // 11% for all bins
  TH1* h = NULL;
  if(hBins){
    h = GFHistManip::CreateAnalog(hBins,Form("sysErrTrackEff%s%d",var,fXSections->GetTrigger()),
				  kFALSE);
  } else {
    this->Info("GetSysErrorTrackEff", "no hBins given, take GetHist(%s,...)", var);
    h = this->GetHist(var, "sysErrTrackEff");
    if(!h) return NULL;
  }
  h->SetTitle("Syst. uncert. of track eff");
  
  Double_t total = this->GetSysErrorTrackEff();
  for(Int_t bin = 1; bin <= h->GetNbinsX(); ++bin){
    h->SetBinContent(bin, total);
  }
  this->PrintTable(h, h->GetName());
  return h;
}

//____________________________________
 Double_t GFDstarSysErrors::GetSysErrorTrackEff()
{
  return fSysTrackEff;
}


// //____________________________________
// TH1* GFDstarSysErrors::GetSysErrorAcceptMc(const char* var)
// {
//   // comparing different MC's for accept

//   TH1* histMc = fXSections->GetDstarHistsAnalysisMc()->CreateAccept(var);
//   TH1* histMcRef = fXSections->GetDstarHistsAnalysisOtherMc()->CreateAccept(var);
//   return this->GetRelDeviation(histMc, histMcRef, var);
// }

//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorAcceptData(const char* var, const TH1* hBins)
{
  // hardcoded values: ET33: bins for wGammaP, rest end ET44: fix values 5.8%/10% (4.8%)
  TH1* h = NULL;
  if(hBins){
    h = GFHistManip::CreateAnalog(hBins, Form("sysErrAcceptData%s%d", var, 
					      fXSections->GetTrigger()), kFALSE);
  } else {
    this->Info("GetSysErrorAcceptData", "no hBins given, take GetHist(%s,...)", var);
    h = this->GetHist(var, "sysErrAcceptData");
    if(!h) return NULL;
  }
  h->SetTitle("Syst. uncert. of A_{33/44}");

  if(fXSections->GetTrigger() == 83 && !strcmp(var, "wGammaP")){
    // cf. log book begin of april 2003
    h->SetBinContent(1, 0.119);
    h->SetBinContent(2, 0.038);
    h->SetBinContent(3, 0.029);
    h->SetBinContent(4, 0.044);
  } else {
    Double_t overall = this->GetSysErrorAcceptData();
    for(Int_t bin = 1; bin <= h->GetNbinsX(); ++bin){
      h->SetBinContent(bin, overall);
    }
  }
  this->PrintTable(h, h->GetName());
  return h;
}

//____________________________________
 Double_t GFDstarSysErrors::GetSysErrorAcceptData()
{
  if(fXSections->GetTrigger() == 83){
    return fSysTagTotal83;
  } else if(fXSections->GetTrigger() == 84){
    // 0.048 // from Vladimirs old mail
    return 0.10; // but leave out curve variation...(cf. log book 9.4.03)
  } else {
    return 111.;
  }
}


//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorAcceptCurve(const char* var)
{
  // difference from max deviation of left/right moved acceptacne curve (c.f. logbook 8.4.03)
  if(fXSections->GetTrigger() != 84){
    return NULL;
  }
  fHistManager->Clear();
  
  TH1* hleft = fAnalysisMcLeft->CreateAccept(var);
  TH1* hright= fAnalysisMcRight->CreateAccept(var);
  TH1* hnorm = fXSections->GetDstarHistsAnalysisMc()->CreateAccept(var);

  fHistManager->AddHist(hright, 0, "left", "l");
  fHistManager->AddHistSame(hleft, 0, 0, "right", "l");
  fHistManager->AddHistSame(hnorm, 0, 0, "norm", "l");
  //man->AddLegend(0,0, "efficiency");
  fHistManager->SetHistsOption("e");

  TH1* hRatioLeft = (TH1*) hleft->Clone("ratioLeft");
  hRatioLeft->SetOption("HIST");
  hRatioLeft->Add(hnorm, -1.);
  hRatioLeft->Divide(hnorm);
  //  hRatioLeft->SetTitle("(left - normal) / normal");
  fHistManager->AddHist(hRatioLeft,0, "left");
  fHistManager->AddLegend(0,0, "ET44: A_{44}");
  
  TH1* hRatioRight = (TH1*) hright->Clone("ratioRight");
  hRatioRight->SetOption("HIST");
  hRatioRight->Add(hnorm, -1.);
  hRatioRight->Divide(hnorm);
  //hRatioRight->SetTitle("(right - normal) / normal");
  fHistManager->AddHistSame(hRatioRight,0, 1, "right");
  fHistManager->AddLegend(0,1, "ET44: A_{44} rel. var.");


  const Double_t fudgeFactor = 0.5;
  TH1* result = this->GetHist(var, "sysErrAcceptCurve", "syst. uncert.: extreme A_{44} curves");
  for(Int_t bin = 1; bin <= result->GetNbinsX(); ++bin){
    Double_t absDevRight = TMath::Abs(hRatioRight->GetBinContent(bin));
    Double_t absDevLeft = TMath::Abs(hRatioLeft->GetBinContent(bin));
    result->SetBinContent(bin, fudgeFactor * TMath::Max(absDevLeft, absDevRight));
  }
  this->PrintTable(result, TString(result->GetName())+=var);
  fHistManager->Draw();
  return result;
}


//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorTrigEff(const char* var, const TH1* hBins)
{
  TH1* h = NULL;
  if(hBins){
    h = GFHistManip::CreateAnalog(hBins, Form("sysErrTrigEff%s%d", var, 
					      fXSections->GetTrigger()), kFALSE);
  } else {
    this->Info("GetSysErrorTrigEff", "no hBins given, take GetHist(%s,...)", var);
    h = this->GetHist(var, "sysErrTrigEff");
    if(!h) return NULL;
  }
  h->SetTitle("Syst. uncert. of #epsilon_{trig}");

  const Double_t overall = this->GetSysErrorTrigEff();
  for(Int_t bin = 1; bin <= h->GetNbinsX(); ++bin){
    h->SetBinContent(bin, overall);
  }
  
  return h;
//   TH1* histData = fXSections->GetDstarHistsAnalysisData()->CreateTrigEff(var);
//   TH1* histMc = fXSections->GetDstarHistsAnalysisMc()->CreateTrigEff(var);
//   return this->GetRelDeviation(histMc, histData, var, 
// 			       fXSections->GetDstarHistsAnalysisMc()->GetPeriod(), 
// 			       fXSections->GetDstarHistsAnalysisData()->GetPeriod());
  
//   // ET33: statistic of trig eff from data
//   // ET44: comparing different MC's for eff trig

//   if(fAnalysisMcRef->GetTrigger() == 83){
//     return this->GetSysErrorTrigEffData(var);
//   } else if(fAnalysisMcRef->GetTrigger() == 84){
//     return this->GetSysErrorTrigEff84(var);
//   } else {
//     return NULL;
//   }
}

//____________________________________
 Double_t GFDstarSysErrors::GetSysErrorTrigEff()
{
  return fSysTrigEff;
}

//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorL4(const char* var, const TH1* hBins)
{
  TH1* h = NULL;
  if(hBins){
    h = GFHistManip::CreateAnalog(hBins, Form("sysErrL4%s%d", var, fXSections->GetTrigger()), 
				  kFALSE);
  } else {
    this->Info("GetSysErrorL4", "no hBins given, take GetHist(%s,...)", var);
    h = this->GetHist(var, "sysErrL4");
    if(!h) return NULL;
  }
  h->SetTitle("Syst. uncert. of #epsilon_{L4}");

  const Double_t overall = this->GetSysErrorL4();
  for(Int_t bin = 1; bin <= h->GetNbinsX(); ++bin){
    h->SetBinContent(bin, overall);
  }
  
  return h;
}

//____________________________________
 Double_t GFDstarSysErrors::GetSysErrorL4()
{
  return fSysL4;
}

//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorParticleID(const char* var, const TH1* hBins)
{
  TH1* h = NULL;
  if(hBins){
    h = GFHistManip::CreateAnalog(hBins, Form("sysErrParticleID%s%d", var, 
					      fXSections->GetTrigger()), kFALSE);
  } else {
    this->Info("GetSysErrorParticleID", "no hBins given, take GetHist(%s,...)", var);
    h = this->GetHist(var, "sysErrParticleID");
    if(!h) return NULL;
  }
  h->SetTitle("Syst. uncert. of particle ID");

  const Double_t overall = this->GetSysErrorParticleID();
  for(Int_t bin = 1; bin <= h->GetNbinsX(); ++bin){
    h->SetBinContent(bin, overall);
  }
  
  return h;
}

//____________________________________
 Double_t GFDstarSysErrors::GetSysErrorParticleID()
{
  return fSysParticleID;
}

// //____________________________________
// TH1* GFDstarSysErrors::GetSysErrorTrigEffData(const char* var)
// {
//   // statistic of trig eff from data (to be for ET33)
  
//   TH1* histData = fXSections->GetDstarHistsAnalysisData()->CreateTrigEff(var);
//   TH1* result = this->GetHist(var, "relTrigErrData", "Syst. trigger uncert. (Data)");
//   for(Int_t bin = 1; bin <= result->GetNbinsX(); ++bin){
//     Double_t content = histData->GetBinContent(bin);
//     Double_t error = histData->GetBinError(bin);
//     result->SetBinContent(bin, content ? error/content: 10000.);
//   }
//   this->PrintTable(result,TString("relErrorTrigEffData ")+=var);
//   return result;
// }

//____________________________________
// Double_t GFDstarSysErrors::GetSysErrorTrigEffData()
// {
//   // statistic of trig eff from data (to be for ET33)
//   TArrayD effErr = fXSections->GetDstarHistsAnalysisData()->TotalTrigEff();
//   return effErr[1] / (effErr[0] ? effErr[0] : 0.0000001);
// }


// //____________________________________
// TH1* GFDstarSysErrors::GetSysErrorTrigEffMc(const char* var)
// {
//   // comparing different MC's for trig rec

//   TH1* histMc = fXSections->GetDstarHistsAnalysisMc()->CreateTrigEff(var);
//   TH1* histMcRef = fXSections->GetDstarHistsAnalysisOtherMc()->CreateTrigEff(var);
//   return this->GetRelDeviation(histMc, histMcRef, var);
// }


//____________________________________
// TH1* GFDstarSysErrors::GetSysErrorTrigEff84(const char* var)
// {
//   // mixin from MC (TE19) and data (TE31)
//   const Int_t oldTrig = fXSections->GetTrigger();
//   if(oldTrig != 84)  this->SetTrigger(84);
//   TH1* hMcEff19   = fXSections->GetDstarHistsAnalysisMc()->CreateTrigEff(var, 19);
//   TH1* hDataEff31 = fXSections->GetDstarHistsAnalysisData()->CreateTrigEff(var, 31);
//   Double_t dEffMean19 = fXSections->GetDstarHistsAnalysisData()->TotalTrigEff(19).At(0);
//   Double_t dEffMean31 = fXSections->GetDstarHistsAnalysisData()->TotalTrigEff(31).At(0);

//   TH1* result = this->GetHist(var, "sysUncertMix", 
// 			      "mixed data/mc sys. uncert, #epsilon_{trig}");
//   for(Int_t bin = 1; bin <= result->GetNbinsX(); ++bin){
//     Double_t mcEff19 = hMcEff19->GetBinContent(bin);
//     Double_t relSysErr19 = (mcEff19 - dEffMean19) / (mcEff19 ? mcEff19 :.00001);// abs! but '^2'
    
//     Double_t dEff31    = hDataEff31->GetBinContent(bin);
//     Double_t dEffErr31 =  hDataEff31->GetBinError(bin);
//     if(!dEffErr31) {
//       this->Warning("GetSysErrorTrigEff84", "error 0. take diff to mean (bin %d, %s)!", 
// 		    bin, var);
//       dEffErr31 = TMath::Abs(dEff31 - dEffMean31);
//     }
//     Double_t relSysErr31 = dEffErr31/ (dEff31 ? dEff31 : .000001);

//     result->SetBinContent(bin, TMath::Sqrt(relSysErr19*relSysErr19 + relSysErr31*relSysErr31));
//   }
//   if(oldTrig != 84) this->SetTrigger(oldTrig);
//   this->PrintTable(result, TString("sys Err effTrig84 ")+=var);

//   return result;

  //  return this->GetRelDeviation(histMc, histMcRef, var);
// }


//____________________________________
// Double_t GFDstarSysErrors::GetSysErrorTrigEff84()
// {
//   // mixin from MC (TE19) and data (TE31)
//   const Int_t oldTrig = fXSections->GetTrigger();
//   if(oldTrig != 84)  this->SetTrigger(84);

//   TArrayD mcEff19 = fXSections->GetDstarHistsAnalysisMc()->TotalTrigEff(19);
//   TArrayD dataEff19 = fXSections->GetDstarHistsAnalysisData()->TotalTrigEff(19);
  
//   Double_t err19 = this->GetRelDeviation(mcEff19[0], dataEff19[0]);

//   TArrayD dataEff31 = fXSections->GetDstarHistsAnalysisData()->TotalTrigEff(31);
//   Double_t err31 = dataEff31[1] / (dataEff31[0] ? dataEff31[0] : 0.0000001); 
  
//   if(oldTrig != 84) this->SetTrigger(oldTrig);
//   return TMath::Sqrt(err19*err19 + err31*err31);
// }

//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorBranching(const char* var, const TH1* hBins)
{
  TH1* h = NULL;
  if(hBins){
    h = GFHistManip::CreateAnalog(hBins, Form("sysErrBranching%s%d", var, 
					      fXSections->GetTrigger()), kFALSE);
  } else {
    this->Info("GetSysErrorBranching", "no hBins given, take GetHist(%s,...)", var);
    h = this->GetHist(var, "sysErrBranching");
    if(!h) return NULL;
  }
  h->SetTitle("Syst. uncert. due to BR(D*#toK#pi#pi_{s})");

  const Double_t overall = fSysBranch;
  for(Int_t bin = 1; bin <= h->GetNbinsX(); ++bin){
    h->SetBinContent(bin, overall);
  }
  
  return h;
}

//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorReflections(const char* var, const TH1* hBins)
{
  TH1* h = NULL;
  if(hBins){
    h = GFHistManip::CreateAnalog(hBins, Form("sysErrReflections%s%d", var, 
					      fXSections->GetTrigger()), kFALSE);
  } else {
    this->Info("GetSysErrorReflections", "no hBins given, take GetHist(%s,...)", var);
    h = this->GetHist(var, "sysErrReflections");
  }
  h->SetTitle("Syst. uncert. of reflection correction");

  const Double_t overall = fSysReflect;
  for(Int_t bin = 1; bin <= h->GetNbinsX(); ++bin){
    h->SetBinContent(bin, overall);
  }
  
  return h;
}

//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorLumi(const char* var, const TH1* hBins)
{
  TH1* h = NULL;
  if(hBins){
    h = GFHistManip::CreateAnalog(hBins, Form("sysErrLumi%s%d", var, 
					      fXSections->GetTrigger()), kFALSE);
  } else {
    this->Info("GetSysErrorLumi", "no hBins given, take GetHist(%s,...)", var);
    h = this->GetHist(var, "sysErrLumi");
    if(!h) return NULL;
  }
  h->SetTitle("Syst. uncert. of luminosity");
  const Double_t overall = fSysLumi;
  for(Int_t bin = 1; bin <= h->GetNbinsX(); ++bin){
    h->SetBinContent(bin, overall);
  }
  
  return h;
}

//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorModel(const char* var, Bool_t norm)
{
  // GetSysErrorRecEff, GetSysErrorAccept, GetSysErrorTrigEff

  // acceptance
  GFDstarHistsAnalysisMc *mc = fXSections->GetDstarHistsAnalysisMc();
  GFDstarHistsAnalysisMc *mcRef = fXSections->GetDstarHistsAnalysisOtherMc();
  TH1* histMc = mc->CreateAccept(var);
  TH1* histMcRef = mcRef->CreateAccept(var);

  if (!histMc || !histMcRef) return NULL;

  // reconstruction
  histMc->Multiply(mc->CreateRecEff(var));
  histMcRef->Multiply(mcRef->CreateRecEff(var));
  // trigger L1
  histMc->Multiply(mc->CreateTrigEff(var));
  histMcRef->Multiply(mcRef->CreateTrigEff(var));

  histMc->SetName(Form("accRecTrig%d%s", fXSections->GetTrigger(), 
		       mc->GetPeriod()));
  histMc->SetName(Form("accRecTrig%d%s", fXSections->GetTrigger(), 
		       mcRef->GetPeriod()));
  histMc->SetTitle("D*: A_{33} * #epsilon_{rec} * #epsilon_{L1}");
  if (norm) {
    histMc->SetTitle(Form("%s (norm.)", histMc->GetTitle()));
    const Float_t tot    = mc   ->TotalRecEff()[0] * mc   ->TotalTrigEff()[0] 
      * mc   ->TotalAccept()[0];
    const Float_t totRef = mcRef->TotalRecEff()[0] * mcRef->TotalTrigEff()[0] 
      * mcRef->TotalAccept()[0];
    if (tot && totRef) {
      histMc   ->Scale(1./tot   );
      histMcRef->Scale(1./totRef);
    } else {
      this->Error("GetSysErrorModel", "problems with totals for normalisation");
      delete histMc;
      delete histMcRef;
      return NULL;
    }

  }

  const TString tit(Form("D*, bins of %s", var));
  return this->GetSysErrorModel(histMc, histMcRef, tit, var, norm);
}

//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorModel(TH1 *hMc, TH1 *hMcRef, const char *tit, 
					const char *var, Bool_t norm)
{
  // half of deviation is taken, but minimum is GetSysErrorModel() for total sample!
  if(!hMc || !hMcRef) return NULL;
  const Double_t min = (norm ? this->GetSysErrorModel() : 0.01);
  TH1* hDev = this->GetRelDeviation(hMc, hMcRef, var);
  TH1* hSysErr = 
    GFHistManip::CreateAnalog(hDev, Form("sysErrModel%s%d", var,
					 fXSections->GetTrigger()), kFALSE, kTRUE);
  hSysErr->SetTitle(Form("syst. uncertainty of model%s", (norm ? " (norm.)" : "")));
  hSysErr->Scale(0.5);
  for (Int_t i = 0; i <= hSysErr->GetNbinsX(); ++i) {// underflow contains all else
    if(min > hSysErr->GetBinContent(i)){          // overflow is empty
      hSysErr->SetBinContent(i, min);
    }
  }

  this->PrintTable(hSysErr, Form("systematic uncertainty of model for %s%s", tit,
				 (norm ? " (norm.)" : "")));
  fHistManager->AddHist(hSysErr, 2);
  fHistManager->Draw();

  return hSysErr;
}

//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorModel(const char* var1, const char* var2, 
					Int_t firstBin2, Int_t lastBin2, Bool_t norm)
{ 
  //pt,eta,wGammaP

  GFDstarHistsAnalysisMc *mc = fXSections->GetDstarHistsAnalysisMc();
  GFDstarHistsAnalysisMc *mcRef = fXSections->GetDstarHistsAnalysisOtherMc();
  // acceptance
  TH1* histMc = mc->CreateAccept(var1, var2, firstBin2, lastBin2);
  TH1* histMcRef = mcRef->CreateAccept(var1, var2, firstBin2,lastBin2);
  // reconstruction
  histMc->Multiply(mc->CreateRecEff(var1, var2, firstBin2, lastBin2));
  histMcRef->Multiply(mcRef->CreateRecEff(var1, var2, firstBin2, lastBin2));
  // trigger L1
  histMc->Multiply(mc->CreateTrigEff(var1, var2, firstBin2, lastBin2));
  histMcRef->Multiply(mcRef->CreateTrigEff(var1, var2, firstBin2, lastBin2));

  histMc->SetName(Form("accRecTrig%d%s", fXSections->GetTrigger(),
		       mc->GetPeriod()));
  histMc->SetName(Form("accRecTrig%d%s", fXSections->GetTrigger(),
		       mcRef->GetPeriod()));
  histMc->SetTitle("D* double diff: A_{33}*#epsilon_{rec} * #epsilon_{L1}");

  if (norm) {
    histMc->SetTitle(Form("%s (norm.)", histMc->GetTitle()));
    const Float_t tot    = mc   ->TotalRecEff()[0] * mc   ->TotalTrigEff()[0] 
      * mc   ->TotalAccept()[0];
    const Float_t totRef = mcRef->TotalRecEff()[0] * mcRef->TotalTrigEff()[0] 
      * mcRef->TotalAccept()[0];
    if (tot && totRef) {
      histMc   ->Scale(1./tot   );
      histMcRef->Scale(1./totRef);
    } else {
      this->Error("GetSysErrorModel", "problems with totals for normalisation");
      delete histMc;
      delete histMcRef;
      return NULL;
    }

  }

  const TString tit(Form("D*, bins of %s %s", var1, var2));
  return this->GetSysErrorModel(histMc, histMcRef, tit, var1, norm);
}


//____________________________________
 Double_t GFDstarSysErrors::GetSysErrorModel()
{
  TArrayD mcRec = fXSections->GetDstarHistsAnalysisMc()->TotalRecEff();
  TArrayD mcRecRef = fXSections->GetDstarHistsAnalysisOtherMc()->TotalRecEff();

  TArrayD mcAcc = fXSections->GetDstarHistsAnalysisMc()->TotalAccept();
  TArrayD mcAccRef = fXSections->GetDstarHistsAnalysisOtherMc()->TotalAccept();

  TArrayD mcTrig = fXSections->GetDstarHistsAnalysisMc()->TotalTrigEff();
  TArrayD mcTrigRef = fXSections->GetDstarHistsAnalysisOtherMc()->TotalTrigEff();

  const Double_t deviation = this->GetRelDeviation(mcAcc[0]*mcRec[0]*mcTrig[0], 
						   mcAccRef[0]*mcRecRef[0]*mcTrigRef[0]);

  return deviation/2.;
}

// //____________________________________
// TH1* GFDstarSysErrors::GetSysErrorTrigDsJet(const TH1* histBins, 
// 					    const char* var,const char* forwBack)
// {
//   return this->CreateConstHist(histBins, fSysDsJetTrig, Form("trig: %s%s", var, forwBack));
// }

// //____________________________________
// TH1* GFDstarSysErrors::GetSysErrorFitDsJet(const TH1* histBins, 
// 					   const char* var,const char* forwBack)
// {
//   return this->CreateConstHist(histBins, fSysFit, Form("fit: %s%s", var, forwBack));
// }


//____________________________________
 TH1* GFDstarSysErrors::CalcSysErrorFitDsJet(const char* var,const char* forwBack, 
					    Int_t fitOrNotFlag)
{
  // normal fit vs fitting WC background and subtracting its integral from signal region
  TH1* histNorm = fXSections->GetDstarHistsAnalysisData()->DrawDmDs1Jet(var, forwBack);
  const Int_t oldFlag = fXSections->GetDstarHistsAnalysisData()->GetMultiDmFitter()
    ->GetFitOrNotFlag();
  TH1* histCheck = NULL;
  switch(fitOrNotFlag){
  case -1:
    this->Error("CalcSysErrorFitDsJet", "There is no fit with WC, Neo!");
//     histCheck = fXSections->GetDstarHistsAnalysisData()->DrawDmFitWc(var);
    break;
  case 0: case 1: case 2: case 3: case 4: case 5:
    fXSections->GetDstarHistsAnalysisData()->GetMultiDmFitter()->SetFitOrNotFlag(fitOrNotFlag);
    histCheck = fXSections->GetDstarHistsAnalysisData()->DrawDmDs1Jet(var, forwBack);
    break;
  default:
    this->Warning("CalcSysErrorFit", "Don't know fitOrNotFlag = %d, use 5", fitOrNotFlag);
    fXSections->GetDstarHistsAnalysisData()->GetMultiDmFitter()->SetFitOrNotFlag(5);
    histCheck = fXSections->GetDstarHistsAnalysisData()->DrawDmDs1Jet(var, forwBack);
  }

  fXSections->GetDstarHistsAnalysisData()->GetMultiDmFitter()->SetFitOrNotFlag(oldFlag);

  TH1* hSysErr = this->GetRelDeviation(histNorm, histCheck, var, "normal fit", 
				       Form("fit flag %d",fitOrNotFlag), "R");
  return hSysErr;
}

//____________________________________
 TH1* GFDstarSysErrors::CalcSysErrorFitDsDiJet(const char *var, Int_t fitOrNotFlag)
{
  // normal fit vs fitting WC background and subtracting its integral from signal region
  TH1* histNorm = fXSections->GetDstarHistsAnalysisData()->DrawDmDsDiJet(var);
  const Int_t oldFlag = fXSections->GetDstarHistsAnalysisData()->GetMultiDmFitter()
    ->GetFitOrNotFlag();
  TH1* histCheck = NULL;
  switch(fitOrNotFlag){
  case -1:
    this->Error("CalcSysErrorFitDsJet", "There is no fit with WC, Neo!");
//     histCheck = fXSections->GetDstarHistsAnalysisData()->DrawDmFitWc(var);
    break;
  case 0: case 1: case 2: case 3: case 4: case 5:
    fXSections->GetDstarHistsAnalysisData()->GetMultiDmFitter()->SetFitOrNotFlag(fitOrNotFlag);
    histCheck = fXSections->GetDstarHistsAnalysisData()->DrawDmDsDiJet(var);
    break;
  default:
    this->Warning("CalcSysErrorFitDsDiJet","Don't know fitOrNotFlag = %d, use 5",fitOrNotFlag);
    fXSections->GetDstarHistsAnalysisData()->GetMultiDmFitter()->SetFitOrNotFlag(5);
    histCheck = fXSections->GetDstarHistsAnalysisData()->DrawDmDsDiJet(var);
  }

  fXSections->GetDstarHistsAnalysisData()->GetMultiDmFitter()->SetFitOrNotFlag(oldFlag);

  TH1* hSysErr = this->GetRelDeviation(histNorm, histCheck, var, "normal fit", 
				       Form("fit flag %d",fitOrNotFlag), "R");
  return hSysErr;
}

//____________________________________
 Double_t GFDstarSysErrors::CalcSysErrorFitDsJet(const char *forwBack, 
						Int_t fitOrNotFlag)
{
  TArrayD normArr = fXSections->GetDstarHistsAnalysisData()->DrawDmDsJet(forwBack);
  const Int_t oldFlag = fXSections->GetDstarHistsAnalysisData()->GetMultiDmFitter()
    ->GetFitOrNotFlag();
  TArrayD checkArr;
  switch(fitOrNotFlag){
  case -1:
    this->Error("CalcSysErrorFitDsJet", "There is no fit with WC, Neo!");
//     checkArr = fXSections->GetDstarHistsAnalysisData()->DrawDmFitWc();
    break;
  case 0: case 1: case 2: case 3: case 4: case 5:
    fXSections->GetDstarHistsAnalysisData()->GetMultiDmFitter()->SetFitOrNotFlag(fitOrNotFlag);
    checkArr = fXSections->GetDstarHistsAnalysisData()->DrawDmDsJet(forwBack);
    break;
  default:
    this->Warning("CalcSysErrorFit", "Don't know fitOrNotFlag = %d, use 5", fitOrNotFlag);
    fXSections->GetDstarHistsAnalysisData()->GetMultiDmFitter()->SetFitOrNotFlag(5);
    checkArr = fXSections->GetDstarHistsAnalysisData()->DrawDmDsJet(forwBack);
  }

  fXSections->GetDstarHistsAnalysisData()->GetMultiDmFitter()->SetFitOrNotFlag(oldFlag);

  const Double_t res = this->GetRelDeviation(normArr[0], checkArr[0]);
  char tmp[] = "------------------------------------------------n";
  cout << tmp << "Rel. deviation (flag " << fitOrNotFlag << ") of N(D*+" << forwBack 
       << "jet) is " << res 
       << ", rel. stat. error is " << normArr[1]/normArr[0] << "n"
       << tmp << endl;

  return res;
}

//____________________________________
 TH1* GFDstarSysErrors::CalcSysErrorScaleUpDsJet(const char* var,
					      const char* forwBack)
{
  TH1* histNorm = fXSections->GetDstarHistsAnalysisData()->DrawDmDs1Jet(var, forwBack);
  TH1* histUp = fAnalysisDataUp->DrawDmDs1Jet(var, forwBack);

  return this->GetRelDeviation(histUp, histNorm, var, "calo scale up", "normal scale", "R");
}

//____________________________________
 TH1* GFDstarSysErrors::CalcSysErrorScaleDownDsJet(const char* var,
						  const char* forwBack)
{
  TH1* histNorm = fXSections->GetDstarHistsAnalysisData()->DrawDmDs1Jet(var, forwBack);
  TH1* histDown = fAnalysisDataDown->DrawDmDs1Jet(var, forwBack);

  return this->GetRelDeviation(histDown, histNorm, var, "calo scale down", "normal scale", "R");
}

//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorScaleDsJet(const TH1* histBins, 
					     const char* var,const char* forwBack)
{
  return this->CreateConstHist(histBins, fSysDsJetScale, Form("scale: %s%s", var, forwBack));
}

//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorScaleDsDiJet(const TH1* histBins, const char* var)
{
  return this->CreateConstHist(histBins, fSysDsJetScale, Form("scale: %s", var));
}

//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorTotalDsJet(const char* var,const char* forwBack, 
					     Bool_t norm)
{
  GFHistArray uncertHists;

  TString combVar(Form("%s%s", var, forwBack));

  const TH1* hModel = this->GetSysErrorModelDsJet(var, forwBack);
  uncertHists.Add(hModel->Clone());
  
  if (!norm) {
    uncertHists.Add(this->GetSysErrorTrigEff(combVar, hModel));
    uncertHists.Add(this->GetSysErrorL4(combVar, hModel));
    uncertHists.Add(this->GetSysErrorAcceptData(combVar, hModel));
    uncertHists.Add(this->GetSysErrorTrackEff(combVar, hModel));
    uncertHists.Add(this->GetSysErrorParticleID(combVar, hModel));
    uncertHists.Add(this->GetSysErrorBranching(combVar, hModel));  // not binwise?
    uncertHists.Add(this->GetSysErrorReflections(combVar, hModel));// not binwise?
    uncertHists.Add(this->GetSysErrorFit(combVar, hModel));
  //  uncertHists.Add(this->CalcSysErrorFitDsJet(var,forwBack));// or this?
    uncertHists.Add(this->GetSysErrorLumi(combVar, hModel));       // not binwise?
  // energy scale:
    uncertHists.Add(this->GetSysErrorScaleDsJet(hModel,var,forwBack));
  } else {
    this->Warning("GetSysErrorTotalDsJet", "still ignore E-scale");
  }

  return this->CombineErrorHists(&uncertHists, combVar);
}

//____________________________________
 Double_t GFDstarSysErrors::GetSysErrorTotalDsJet(const char* forwBack)
{
  Double_t total2 = 0.;

  const Double_t tri = this->GetSysErrorTrigEff();
  total2 += (tri*tri);

  const Double_t l4 = this->GetSysErrorL4();
  total2 += (l4*l4);

  const Double_t ad = this->GetSysErrorAcceptData();
  total2 += (ad*ad);

  const Double_t tr = this->GetSysErrorTrackEff();
  total2 += (tr*tr);

  const Double_t id = this->GetSysErrorParticleID();
  total2 += (id*id);

  total2 += (fSysBranch*fSysBranch); //BRanching dstar->Kpipi, correlated!
  total2 += (fSysReflect*fSysReflect); // correlated!

  const Double_t fit = this->GetSysErrorFit();
  total2 += (fit*fit);

  const Double_t mod = this->GetSysErrorModelDsJet(forwBack);
  total2 += (mod*mod);

  total2 += (fSysLumi*fSysLumi); //lumi, correlated

  total2 += (fSysDsJetScale*fSysDsJetScale);

  return TMath::Sqrt(total2);
}

//____________________________________
 Double_t GFDstarSysErrors::GetSysErrorTotalDsDiJet()
{
  Double_t total2 = 0.;

  const Double_t tri = this->GetSysErrorTrigEff();
  total2 += (tri*tri);

  const Double_t l4 = this->GetSysErrorL4();
  total2 += (l4*l4);

  const Double_t ad = this->GetSysErrorAcceptData();
  total2 += (ad*ad);

  const Double_t tr = this->GetSysErrorTrackEff();
  total2 += (tr*tr);

  const Double_t id = this->GetSysErrorParticleID();
  total2 += (id*id);

  total2 += (fSysBranch*fSysBranch); //BRanching dstar->Kpipi, correlated!
  total2 += (fSysReflect*fSysReflect); // correlated!

  const Double_t fit = this->GetSysErrorFit();
  total2 += (fit*fit);

  const Double_t mod = this->GetSysErrorModelDsDiJet();
  total2 += (mod*mod);

  total2 += (fSysLumi*fSysLumi); //lumi, correlated

  total2 += (fSysDsJetScale*fSysDsJetScale);

  return TMath::Sqrt(total2);
}


//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorTotal(const char* var1, const char* var2,
					Int_t firstBin2, Int_t lastBin2, Bool_t norm)
{
  GFHistArray uncertHists;

  TString combVar(Form("%s%s%d%d", var1, var2, firstBin2, lastBin2));

  const TH1* hModel = this->GetSysErrorModel(var1, var2, firstBin2, lastBin2);
  uncertHists.Add(hModel->Clone());
  
  if (!norm) {
    uncertHists.Add(this->GetSysErrorTrigEff(combVar, hModel));
    uncertHists.Add(this->GetSysErrorL4(combVar, hModel));
    uncertHists.Add(this->GetSysErrorAcceptData(combVar, hModel));
    uncertHists.Add(this->GetSysErrorTrackEff(combVar, hModel));
    uncertHists.Add(this->GetSysErrorParticleID(combVar, hModel));
    uncertHists.Add(this->GetSysErrorBranching(combVar, hModel));  // not binwise?
    uncertHists.Add(this->GetSysErrorReflections(combVar, hModel));// not binwise?
    uncertHists.Add(this->GetSysErrorFit(combVar, hModel));
  //  uncertHists.Add(this->CalcSysErrorFitDsJet(var,forwBack));// or this?
    uncertHists.Add(this->GetSysErrorLumi(combVar, hModel));       // not binwise?


    if(strcmp(var1, "wGammaP") == 0 || strcmp(var2, "wGammaP") == 0){
      this->Error("GetSysErrorTotal",
		  "apply systematic uncertainty of total W-range!");
    }
  }

  return this->CombineErrorHists(&uncertHists, combVar);
}

//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorTotalDsDiJet(const char* var, Bool_t norm)
{
  GFHistArray uncertHists;

  const TH1* hModel = this->GetSysErrorModelDsDiJet(var);
  uncertHists.Add(hModel->Clone());
  
  if (!norm) {
    uncertHists.Add(this->GetSysErrorTrigEff(var, hModel));
    uncertHists.Add(this->GetSysErrorL4(var, hModel));
    uncertHists.Add(this->GetSysErrorAcceptData(var, hModel));//FIXME
    uncertHists.Add(this->GetSysErrorTrackEff(var, hModel));
    uncertHists.Add(this->GetSysErrorParticleID(var, hModel));
    uncertHists.Add(this->GetSysErrorBranching(var, hModel));  // not binwise?
    uncertHists.Add(this->GetSysErrorReflections(var, hModel));// not binwise?
    uncertHists.Add(this->GetSysErrorFit(var, hModel));//FIXME
    //  uncertHists.Add(this->CalcSysErrorFitDsDiJet(var));
    uncertHists.Add(this->GetSysErrorLumi(var, hModel));       // not binwise?
    // energy scale:
    uncertHists.Add(this->GetSysErrorScaleDsDiJet(hModel,var));//FIXME
  } else {
    this->Warning("GetSysErrorTotalDsDiJet", "still ignore E-scale");
  }

  return this->CombineErrorHists(&uncertHists, var);
}

//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorModelDsJet(const char* var,const char* forwBack, 
					     Bool_t norm)
{
  // GetSysErrorRecEff, GetSysErrorAccept, GetSysErrorTrigEff

  GFDstarHistsAnalysisMc *mc = fXSections->GetDstarHistsAnalysisMc();
  GFDstarHistsAnalysisMc *mcRef = fXSections->GetDstarHistsAnalysisOtherMc();
  // accept
  TH1* histMc = mc->CreateAcceptDsJet(var, forwBack);
  TH1* histMcRef = mcRef->CreateAcceptDsJet(var, forwBack);
  // rec eff
  histMc->Multiply(mc->CreateRecEffDs1Jet(var, forwBack));
  histMcRef->Multiply(mcRef->CreateRecEffDs1Jet(var, forwBack));
  // trig eff
  histMc->Multiply(mc->CreateTrigEffDsJet(var, forwBack));
  histMcRef->Multiply(mcRef->CreateTrigEffDsJet(var, forwBack));

  histMc->SetTitle(Form("D* %s + 1 jet: #epsilon_{rec}' * #epsilon_{L1} * A_{33}", 
			forwBack));
  if (norm) {
    histMc->SetTitle(Form("%s (norm.)", histMc->GetTitle()));
    const Float_t tot    = mc   ->TotalRecEffDsJet(forwBack)[0] *
      mc   ->TotalTrigEffDsJet(forwBack)[0] * mc   ->TotalAcceptDsJet(forwBack)[0];
    const Float_t totRef = mcRef->TotalRecEffDsJet(forwBack)[0] *
      mcRef->TotalTrigEffDsJet(forwBack)[0] * mcRef->TotalAcceptDsJet(forwBack)[0];
    if (tot && totRef) {
      histMc   ->Scale(1./tot   );
      histMcRef->Scale(1./totRef);
    } else {
      this->Error("GetSysErrorModelDsJet", "problems with totals for normalisation");
      delete histMc;
      delete histMcRef;
      return NULL;
    }

  }


  const TString tit(Form("D*+jet%s %s", forwBack, var));
  return this->GetSysErrorModel(histMc, histMcRef, tit, var, norm);
}

//____________________________________
 Double_t GFDstarSysErrors::GetSysErrorModelDsJet(const char* forwBack)
{
  // from accept, rec and trig eff

  // accept
  Double_t mc = fXSections->GetDstarHistsAnalysisMc()->TotalAcceptDsJet(forwBack).At(0);
  Double_t mcRef = fXSections->GetDstarHistsAnalysisOtherMc()->TotalAcceptDsJet(forwBack).At(0);
  // rec eff
  mc *= fXSections->GetDstarHistsAnalysisMc()->TotalRecEffDsJet(forwBack).At(0);
  mcRef *= fXSections->GetDstarHistsAnalysisOtherMc()->TotalRecEffDsJet(forwBack).At(0);
  // trige eff
  mc *= fXSections->GetDstarHistsAnalysisMc()->TotalTrigEffDsJet(forwBack).At(0);
  mcRef *= fXSections->GetDstarHistsAnalysisOtherMc()->TotalTrigEffDsJet(forwBack).At(0);

  const Double_t deviation = this->GetRelDeviation(mc, mcRef);
  // FIXME: take minimum compared with inclusive D*?
  return deviation/2.;
}

//____________________________________
 TH1* GFDstarSysErrors::GetSysErrorModelDsDiJet(const char* var, Bool_t norm)
{
  // GetSysErrorRecEff, GetSysErrorAccept, GetSysErrorTrigEff

  GFDstarHistsAnalysisMc *mc    = fXSections->GetDstarHistsAnalysisMc();
  GFDstarHistsAnalysisMc *mcRef = fXSections->GetDstarHistsAnalysisOtherMc();
  // accept
  TH1* histMc = mc->CreateAcceptDsDiJet(var);
  TH1* histMcRef = mcRef->CreateAcceptDsDiJet(var);
  // rec eff
  histMc->Multiply(mc->CreateRecEffDsDiJet(var));
  histMcRef->Multiply(mcRef->CreateRecEffDsDiJet(var));
  // trig eff
  histMc->Multiply(mc->CreateTrigEffDsDiJet(var));
  histMcRef->Multiply(mcRef->CreateTrigEffDsDiJet(var));

  histMc->SetTitle(Form("D* dijet: #epsilon_{rec}' * #epsilon_{L1} * A_{33}"));

  if (norm) {
    histMc->SetTitle(Form("%s (norm.)", histMc->GetTitle()));
    const Float_t tot    = mc   ->TotalRecEffDsDiJet()[0] *
      mc   ->TotalTrigEffDsDiJet()[0] * mc   ->TotalAcceptDsDiJet()[0];
    const Float_t totRef = mcRef->TotalRecEffDsDiJet()[0] *
      mcRef->TotalTrigEffDsDiJet()[0] * mcRef->TotalAcceptDsDiJet()[0];
    if (tot && totRef) {
      histMc   ->Scale(1./tot   );
      histMcRef->Scale(1./totRef);
    } else {
      this->Error("GetSysErrorModelDsDiJet","problems: tot. for normalisation?");
      delete histMc;
      delete histMcRef;
      return NULL;
    }
  }

  const TString tit(Form("D* dijet %s", var));
  return this->GetSysErrorModel(histMc, histMcRef, tit, var, norm);
}

//____________________________________
 Double_t GFDstarSysErrors::GetSysErrorModelDsDiJet()
{
  // from accept, rec and trig eff

  // accept
  Double_t mc = fXSections->GetDstarHistsAnalysisMc()->TotalAcceptDsDiJet().At(0);
  Double_t mcRef = fXSections->GetDstarHistsAnalysisOtherMc()->TotalAcceptDsDiJet().At(0);
  // rec eff
  mc *= fXSections->GetDstarHistsAnalysisMc()->TotalRecEffDsDiJet().At(0);
  mcRef *= fXSections->GetDstarHistsAnalysisOtherMc()->TotalRecEffDsDiJet().At(0);
  // trige eff
  mc *= fXSections->GetDstarHistsAnalysisMc()->TotalTrigEffDsDiJet().At(0);
  mcRef *= fXSections->GetDstarHistsAnalysisOtherMc()->TotalTrigEffDsDiJet().At(0);

  const Double_t deviation = this->GetRelDeviation(mc, mcRef);
  // FIXME: take minimum compared with inclusive D*?
  return deviation/2.;
}


// //____________________________________
// TH1* GFDstarSysErrors::GetSysErrorRecDsJet(const char* var,const char* forwBack)
// {
//   TH1* histMc = fXSections->GetDstarHistsAnalysisMc()->CreateRecEffDs1Jet(var, forwBack);
//   TH1* histMcRef = fXSections->GetDstarHistsAnalysisOtherMc()->CreateRecEffDs1Jet(var,forwBack);

//   return this->GetRelDeviation(histMc, histMcRef, Form("%s%s", var, forwBack));
// }

// //____________________________________
// TH1* GFDstarSysErrors::GetSysErrorTrigDsJetMc(const char* var,const char*forwBack)
// {
//   TH1* histMc = fXSections->GetDstarHistsAnalysisMc()->CreateTrigEffDsJet(var, forwBack);
//   TH1* histMcRef = fXSections->GetDstarHistsAnalysisOtherMc()->CreateTrigEffDsJet(var,forwBack);

//   return this->GetRelDeviation(histMc, histMcRef, Form("%s%s", var, forwBack));
// }

// //____________________________________
// TH1* GFDstarSysErrors::GetSysErrorAccDsJet(const char* var,const char* forwBack)
// {
//   TH1* histMc = fXSections->GetDstarHistsAnalysisMc()->CreateAcceptDsJet(var, forwBack);
//   TH1* histMcRef = fXSections->GetDstarHistsAnalysisOtherMc()->CreateAcceptDsJet(var, forwBack);

//   return this->GetRelDeviation(histMc, histMcRef, Form("%s%s", var, forwBack));
// }


//____________________________________
 TH1* GFDstarSysErrors::XSection(const char* var, Bool_t norm)
{
  TH1 *hCross = (norm ? fXSections->XSectionNorm(var) : fXSections->XSection(var));
  if (!hCross) return NULL;

  GFHistArray mcHists; //  FIXME: not OK for norm = kTRUE
  mcHists.Add(fXSections->GetHistManager()->GetHistsOf(1, 0)->At(1));
  mcHists.Add(fXSections->GetHistManager()->GetHistsOf(1, 0)->At(2));
  
  TH1* sysErrHist = this->GetSysErrorTotal(var, norm);
  TH1* h = this->XSection(hCross, sysErrHist, var, &mcHists, NULL, norm);

  return h;
}

// //____________________________________
// TH1* GFDstarSysErrors::XSectionPt()
// {
//   TH1* h = this->XSection("pt");
//   fHistManager->SetLogY(1);
//   return h;
// }

// //____________________________________
// TH1* GFDstarSysErrors::XSectionEta()
// {
//   return this->XSection("eta");
// }

// //____________________________________
// TH1* GFDstarSysErrors::XSectionWgammaP()
// {
//   return this->XSection("wGammaP");
// }

//____________________________________
 TH1* GFDstarSysErrors::XSection(const char* var1, const char* var2, Int_t firstBin2, 
				Int_t lastBin2, Bool_t norm)
{
  TH1* crossHist = (norm ? fXSections->XSectionNorm(var1, var2, firstBin2, lastBin2)
		    : fXSections->XSection(var1, var2, firstBin2, lastBin2));
  if(!crossHist) return NULL;

  GFHistArray mcHists; //  FIXME: not OK for norm = kTRUE
  mcHists.Add(fXSections->GetHistManager()->GetHistsOf(1, 0)->At(1));
  mcHists.Add(fXSections->GetHistManager()->GetHistsOf(1, 0)->At(2));
  
  TH1* sysErrHist = this->GetSysErrorTotal(var1, var2, firstBin2, lastBin2, norm);
  TString varBins(Form("%s%d%d", var2, firstBin2,lastBin2));
  TH1* h = this->XSection(crossHist, sysErrHist, var1, &mcHists, varBins, norm);

  return h;
}

//____________________________________
 TH1* GFDstarSysErrors::XSection(TH1* crossHist, TH1* sysErrHist, 
				const char* var, GFHistArray *mcHists, const char* var2,
				Bool_t norm)
{
  if(!crossHist || !sysErrHist) return NULL;
  fHistManager->Clear();
  fHistManager->AddHist(sysErrHist, 0);
  fHistManager->AddHist(crossHist, 1, (mcHists ? "Data" : NULL), "LP");
  crossHist->SetOption("E1");

  TH1* hCrossAllErr = GFHistManip::CreateAnalog
    (crossHist, Form("xSecAllErr%s%d", var, fXSections->GetTrigger()), kFALSE, kTRUE);
  hCrossAllErr->SetTitle("#sigma, all uncertainties");
  TH1* hStatErr = GFHistManip::CreateAnalog
    (crossHist, Form("xSecStatErr%s%d",var, fXSections->GetTrigger()), kFALSE);
  hStatErr->SetTitle("rel. stat. uncertainties");

  for (Int_t bin = 0; bin <= hCrossAllErr->GetNbinsX()+1; ++bin){
    const Double_t xSec = crossHist->GetBinContent(bin);
    //    hCrossAllErr->SetBinContent(bin, xSec); // done in CreateAnalog
    const Double_t statErrAbs = crossHist->GetBinError(bin);
    const Double_t statErrRel = xSec ? statErrAbs/xSec : 0.;
    const Double_t sysErrRel = sysErrHist->GetBinContent(bin);
    hCrossAllErr->SetBinError(bin, TMath::Sqrt(statErrAbs*statErrAbs 
					       + sysErrRel*sysErrRel*xSec*xSec));
    hStatErr->SetBinContent(bin, statErrRel);
  }
  fHistManager->AddHist(hStatErr,0);
  fHistManager->AddHistSame(hCrossAllErr, 1, 0);
  if (!norm && mcHists) { // GFDstarXSection does not (yet?) foresee MC for XSectionNorm...
    TString label(fXSections->GetDstarHistsAnalysisMc()->GetPeriod());
    if (!norm && !GFMath::IsEqualFloat(fXSections->GetMcFactors()->At(0), 1)) {
      label = Form("%.2f #times %s", fXSections->GetMcFactors()->At(0), 
		   label.Data());
    }
    fHistManager->AddHistSame(mcHists->At(0), 1, 0, label);
    mcHists->At(0)->SetOption("HIST");
    if(fXSections->GetDstarHistsAnalysisOtherMc()) {
      label = fXSections->GetDstarHistsAnalysisOtherMc()->GetPeriod();
      if (!norm && !GFMath::IsEqualFloat(fXSections->GetMcFactors()->At(1), 1.)) {
	label = Form("%.2f #times %s", fXSections->GetMcFactors()->At(1), 
		     label.Data());
      }
      fHistManager->AddHistSame(mcHists->At(1), 1, 0, label); 
      mcHists->At(1)->SetOption("HIST");
    }
  }

  if(fOutFile){
    GFHistArray hists;
    hists.Add(crossHist);
    hists.Add(hCrossAllErr);
//     if(mcHists){
//       hists.AddAll(mcHists);
//     }
    TString label(var);
    if (var2) label += var2;
    if (norm) label += "Norm";
    this->WriteToFile(hists, label);
  }  
//   const Bool_t old = fHistManager->DrawDiffStyle(kFALSE);
  fHistManager->Draw();
//   fHistManager->DrawDiffStyle(old);

  sysErrHist->TAttLine::Copy(*hCrossAllErr);
  hCrossAllErr->SetMarkerColor(sysErrHist->GetMarkerColor());
  fHist2 = hCrossAllErr;
  return crossHist;
}

//____________________________________
 Bool_t GFDstarSysErrors::WriteToFile(const GFHistArray &hists, const char *label)
{
  if(!fOutFile) return kFALSE;

  TDirectory *oldDir = gDirectory;
  fOutFile->cd();
  if(hists[0]) hists[0]->Write(label);
  if(hists[1]) hists[1]->Write(Form("%sAllErr", label));
  if(hists[2]) {
    hists[2]->SetTitle(Form("%s #times %.2f", hists[2]->GetTitle(), 
			    fXSections->GetMcFactors()->At(0)));
    hists[2]->Write(Form("%s%s", label, fXSections->GetDstarHistsAnalysisMc()->GetPeriod()));
  }
  if(hists[3]) {
    hists[3]->SetTitle(Form("%s #times %.2f", hists[3]->GetTitle(), 
			    fXSections->GetMcFactors()->At(1)));
    hists[3]->Write(Form("%s%s",label,fXSections->GetDstarHistsAnalysisOtherMc()->GetPeriod()));
  }

  if(oldDir) oldDir->cd();
  else gROOT->cd();

  return kTRUE;
}

//____________________________________
 TArrayD GFDstarSysErrors::XSectionTot()
{
  TArrayD xSecErr = fXSections->XSectionTot();
  xSecErr.Set(3);
  xSecErr[2] = this->GetSysErrorTotal() * xSecErr[0];
  const char line[] = "------------------------------------------------------------------------------n";
  cout << line << "| systematic error = " << xSecErr[2] << "n" << line << endl;

  return xSecErr;
}

//____________________________________
 TH1* GFDstarSysErrors::XSectionDsJet(const char* var, const char* forwBack, Bool_t norm)
{
  TH1* crossHist = (norm ? fXSections->XSectionDsJetNorm(var, forwBack) : 
		    fXSections->XSectionDsJet(var, forwBack));
  if(!crossHist) return NULL;

  TH1* sysErrHist = this->GetSysErrorTotalDsJet(var, forwBack, norm);

  GFHistArray mcHists; // FIXME: not OK for norm = kTRUE
  mcHists.Add(fXSections->GetHistManager()->GetHistsOf(1, 0)->At(1));
  mcHists.Add(fXSections->GetHistManager()->GetHistsOf(1, 0)->At(2));

  return this->XSection(crossHist, sysErrHist, Form("%s%s", var, forwBack), &mcHists, 
			NULL, norm);
}

//____________________________________
 TArrayD GFDstarSysErrors::XSectionDsJetTot(const char *forwBack)
{
  TArrayD xSecErr = fXSections->XSectionDsJetTot(forwBack);
  xSecErr.Set(3);
  xSecErr[2] = this->GetSysErrorTotalDsJet(forwBack) * xSecErr[0];
  const char line[] = "------------------------------------------------------------------------------n";
  cout << line << "| systematic error = " << xSecErr[2] << "n" << line << endl;

  return xSecErr;
}

//____________________________________
 TH1* GFDstarSysErrors::XSectionDsDiJet(const char* var, Bool_t norm)
{

  TH1 *crossHist = (norm ? fXSections->XSectionDsDiJetNorm(var) :
		    fXSections->XSectionDsDiJet(var));
  if (!crossHist) return NULL;

  TH1 *sysErrHist = this->GetSysErrorTotalDsDiJet(var, norm);

  GFHistArray mcHists; // FIXME: not OK for norm = kTRUE
  mcHists.Add(fXSections->GetHistManager()->GetHistsOf(1, 0)->At(1));
  mcHists.Add(fXSections->GetHistManager()->GetHistsOf(1, 0)->At(2));

  return this->XSection(crossHist, sysErrHist, var, &mcHists, NULL, norm);
}

//____________________________________
 TArrayD GFDstarSysErrors::XSectionDsDiJetTot()
{
  TArrayD xSecErr = fXSections->XSectionDsDiJetTot();
  xSecErr.Set(3);
  xSecErr[2] = this->GetSysErrorTotalDsDiJet() * xSecErr[0];
  const char line[] = "------------------------------------------------------------------------------n";
  cout << line << "| systematic error = " << xSecErr[2] << "n" << line << endl;

  return xSecErr;
}


//____________________________________
//____________________________________
//____________________________________
//____________________________________
 TH1* GFDstarSysErrors::GetRelDeviation(TH1* histMc, TH1* histMcRef, const char* var,
				       const char* legendEntry1, const char* legendEntry2,
				       Option_t* opt)
{
  // draw both hists on top of each other, with legend and legend entries
  // if they are not given, they are taken as MC->GetPeriod()
  // also drawing and returning a hist with the relative deviations
  // if opt contains r or R, the relative error of histMc is plotted on top

  if(!histMc || !histMcRef) return NULL;
  // this could be more flexible: up to several histMcRef, taking max...
  // Then to be used in GetSysErrorAcceptCurve

  const Bool_t relOpt = TString(opt).Contains("R", TString::kIgnoreCase);

  fHistManager->Clear();
  // clone hist without errors, but fill it!
  TH1* histRelDiff = GFHistManip::CreateAnalog(histMc, Form("relDev%s%d", var, fXSections->
							    GetTrigger()), kFALSE, kTRUE);

  if(!legendEntry1) legendEntry1 = fXSections->GetDstarHistsAnalysisMc()->GetPeriod();
  if(!legendEntry2) legendEntry2 = fXSections->GetDstarHistsAnalysisOtherMc()->GetPeriod();
  TString title = legendEntry1;// fXSections->GetDstarHistsAnalysisMc()->GetPeriod();
  (title += " vs ") += legendEntry2;//fXSections->GetDstarHistsAnalysisOtherMc()->GetPeriod();
  histRelDiff->SetTitle(title += ": rel. deviation");

  histRelDiff->Add(histMcRef, -1.);
  histRelDiff->Divide(histMc);
  for(Int_t bin = 0; bin <= histRelDiff->GetNbinsX()+1; ++bin){
    if(histRelDiff->GetBinContent(bin) < 0.){
      histRelDiff->SetBinContent(bin, -histRelDiff->GetBinContent(bin));
    }
  }

  fHistManager->AddHist(histMc, 0, legendEntry1);
  fHistManager->AddHistSame(histMcRef, 0, 0, legendEntry2);
  fHistManager->AddHist(histRelDiff, 0, (relOpt ? "rel. dev." : NULL));
  fHistManager->AddLegend(0, 0, fXSections->GetDstarHistsAnalysisOtherMc()->GetSampleName());
  histRelDiff->SetOption("HIST");
  if(relOpt){
    TH1* hRelStat = GFHistManip::CreateRelErrorHist(histMc);
    fHistManager->AddHistSame(hRelStat, 0, 1, "rel. stat. error");
    this->PrintTable(hRelStat, Form("rel. stat (%s)", GFAxisVariables::GetAxisLabel(var)));
  }
  fHistManager->Draw();

  this->PrintTable(histRelDiff, (title+=" in bins of ") += GFAxisVariables::GetAxisLabel(var));
  return histRelDiff;
}

//____________________________________
Double_t GFDstarSysErrors::GetRelDeviation(Double_t ref, Double_t other)//mcAccRef[0]*mcRecRef[0]);
{
  // releative deviation of the 2nd
  return TMath::Abs((ref-other)/(ref ? ref : 0.0000001));
}

//____________________________________
 TH1* GFDstarSysErrors::CombineErrorHists(const GFHistArray* hists, const char* var)
{
  if(!hists || !hists->GetEntriesFast()) return NULL;
  fHistManager->Clear();
  TH1* h = //this->GetHist(var, "sysErrTot", "Total systematic uncertainty");
    GFHistManip::CreateAnalog(hists->First(),Form("sysErrTot%s%d",var,fXSections->GetTrigger()),
			      kFALSE);
  h->SetTitle("Total systematic uncertainty");

  // add all in squares
  Int_t nBins = h->GetNbinsX();
  for(Int_t i = 0; i < hists->GetEntriesFast(); ++i){
    for(Int_t bin = 1; bin <= nBins; ++bin){
      Double_t errorNew = (*hists)[i]->GetBinContent(bin);
      Double_t error2Sum = h->GetBinContent(bin);
      h->SetBinContent(bin, error2Sum + errorNew*errorNew);
    }
    fHistManager->AddHist((*hists)[i]);
  }

// take square-root
  for(Int_t bin = 1; bin <= nBins; ++bin){
    Double_t err = h->GetBinContent(bin);
    // impossible...    if(err < 0.){ this->Error()...
    h->SetBinContent(bin, TMath::Sqrt(err));
  }
  fHistManager->AddHist(h,1);
  fHistManager->Draw();
  this->PrintTable(h, Form("%s in bins of %s", h->GetName(), var));

  return h;
}


//____________________________________
 TH1* GFDstarSysErrors::GetHist(const char* var, const char* name, const char* title)
{
  // returns hist in bin of 'var', setting name and title if given
  // (adding 'var' and trigger (33/44) info)
  TString v(var);
  TH1* h = fXSections->GetDstarHistsAnalysisMc()->GetHist(v+"HistGenAcc","gen");
  if(!h) h = fXSections->GetDstarHistsAnalysisMc()->GetHist("dstar1Jet" +v+"GenAcc", "genJet");
  if(!h) {
    this->Error("GetHist", "No hist found for '%s'", var);
    return NULL;
  }
  if(name){
    TString n(name);
    h->SetName(n+=v+=fXSections->GetTrigger());
  }
  if(title) h->SetTitle(title+(" "+v));
  h->Reset();
  h->SetMaximum(-1111.);// reset to 'unset'
  h->SetMinimum(-1111.);// reset to 'unset'
  //  h->SetDrawOption("HIST");

  return h;
}

//____________________________________
 TH1* GFDstarSysErrors::CreateConstHist(const TH1* h, Double_t val, const char* tit)
{
  // create a hist with bins as h, fill all bins with val and use title
  if(!h || h->GetDimension() != 1) {
    this->Error("CreateConstHist", "No or wrong dimension (%d) of hist", h?h->GetDimension():0);
    return NULL;
  }
  TString title(Form("Syst. uncert. %s", tit)); 
  TString name(tit);
  fXSections->GetDstarHistsAnalysisData()->MakeStringName(name);
  TH1* result = GFHistManip::CreateAnalog(h, name, kFALSE);
  result->SetTitle(title);

  const Int_t nBins = result->GetNbinsX();// * result->GetNbinsY() * result->GetNbinsZ();
  for(Int_t i = 1; i <= nBins; ++i) result->SetBinContent(i, val);

  return result;
}

//____________________________________
 void GFDstarSysErrors::PrintTable(const TH1 *h, const char* title) const
{
//   fXSections->GetDstarHistsAnalysisData()->PrintTable(h, title);
  GFHistManip::PrintTable(h, title);
}