GFDstarHistsAnalysisMc

// Author: Gero Flucke <mailto:flucke@mail.desy.de>
//____________________________________
// GFDstarHistsAnalysisMc
//
//   Author:      Gero Flucke
//   Date:        June 2nd, 2002
//   last update: $Date: 2005/11/16 00:39:02 $  
//   by:          $Author: flucke $
//

#include <iostream>
#include <ctype.h>

#include <TFile.h>
#include <TStyle.h>
#include <TH1.h>
#include <TH2.h>
#include <TH3.h>
#include <TF1.h>
#include <TGraphAsymmErrors.h>
#include <TLine.h>
#include <TCanvas.h>
#include <TString.h>
#include <TObjString.h>
#include <TAxis.h>
#include <TObjArray.h>

#include "/afs/desy.de/user/f/flucke/h1/root/DstarDmFitter/DstarDmFitter.h"

#include "H1Skeleton/H1RunDependent.h"

#include "GFAnalysis/GFDstarHistsAnalysisMc.h"
#include "GFAnalysis/GFAxisVariables.h"
#include "GFUtils/GFHistManager.h"
#include "GFUtils/GFHistArray.h"
#include "GFUtils/GFDstarRunPeriods.h"
#include "GFUtils/GFHistManip.h"
#include "GFUtils/GFMath.h"

using namespace std;

ClassImp(GFDstarHistsAnalysisMc)

const Double_t GFDstarHistsAnalysisMc::kHighestDmMc = .1699999; // outdated
const Double_t GFDstarHistsAnalysisMc::kBranchRatio = .0257;  // PDG 2002
const Double_t GFDstarHistsAnalysisMc::kBranchRatioMc = .025; // Jeannine's thesis
//const Double_t GFDstarHistsAnalysisMc::kXSecCorrFactor = 23.5/30.;// c->D* correction Jeannine
const Double_t GFDstarHistsAnalysisMc::kXSecCorrFactor = 23.5/29.7;// c->D* corr. Lluis


// // now in 'nb'!
// const Double_t loss = 0.979;// due exclusion of YECL-less runs log book: 7.4.03
// // const Double_t GFDstarHistsAnalysisMc::kLumiMcRes = 2613.84 * 1000. * loss; // A+B
// const Double_t GFDstarHistsAnalysisMc::kLumiMcRes = 1306.93 * 1000. * loss;  // A

// const Double_t GFDstarHistsAnalysisMc::kLumiMcDir = 4. * 174.97 *1000. * loss; // A-D

// const Double_t GFDstarHistsAnalysisMc::kLumiMcExc = 4. * 27.92*1000. * loss // old files
// + 0.26172 * 712.07 * 1000.; // new files, cf. logbook 11.3.2004

// // const Double_t GFDstarHistsAnalysisMc::kLumiMcCasc = 55.62*1000. * loss; // A, B, or C
// const Double_t GFDstarHistsAnalysisMc::kLumiMcCasc = 3.* 55.62*1000. * loss; // A-C

//   } else if(fileName.Contains("DirBeauty")) {
//     fLumiMc = 13362.61 * 1000. * loss;
// //     fPeriod = "DirBeauty";
//   } else if(fileName.Contains("ResBeauty")) {
//     fLumiMc = 73859.6 * 1000. * loss;
// //     fPeriod = "ResBeauty";
//   } else if(fileName.Contains("BeautyEx")) {
//     fLumiMc = (878.3+878.25)* 1000. * loss;
// //     fPeriod = "BeautyEx";


//____________________________________
 GFDstarHistsAnalysisMc::GFDstarHistsAnalysisMc
  this->InitMC(period);
}

//____________________________________
 GFDstarHistsAnalysisMc::GFDstarHistsAnalysisMc(const char* mcFileIn)
  : GFDstarHistsAnalysis(mcFileIn)
{
  TString fileString(mcFileIn);
  if(fileString.Contains("/")){ // remove al directory information
    fileString.Remove(0, fileString.Last('/')+1);
  }
  if(kDefDir && fileString.Contains(kDefName)){
    fileString.Remove(0, fileString.Index(kDefName) + strlen(kDefName));
  }
  fileString.ReplaceAll(".root", NULL);
  
  this->InitMC(fileString.Data()); // hopefully now just 'DirCharmB' or so...
}


//____________________________________
 void GFDstarHistsAnalysisMc::InitMC(const char* period)
{
//     fMultiDmFitter->SetFitOrNotFlag(1);
//  4th digit:    Mean
//  3rd :         Sigma
//  2nd :         BackgrExp
//  1st :         BackgrSqr
//  so
//     mode = 1011 means: fix Mean, Sigma and BackgrSqr,
//                        but fit BackgrExp
//                       (and fit BackgrNorm, NDstar)
  fMultiDmFitter->SetFitModeAll(1000);
  fMultiDmFitter->SetFitModeBins(1011);
  fMultiDmFitter->SetStartUsqr(33.);
  TString fitOption = fMultiDmFitter->GetDstarDmFitter()->GetFitOption();
  if(!fitOption.Contains("LL")) {
    fitOption.ReplaceAll("L", NULL);
    fMultiDmFitter->GetDstarDmFitter()->SetFitOption(fitOption);
  }


  fPeriod = (period && strlen(period) ? period : "anyMC");


}

//____________________________________
 GFDstarHistsAnalysisMc::~GFDstarHistsAnalysisMc()
{
  // nothing to be done

  // fFileOut taken care of in base class
}

//____________________________________
 void GFDstarHistsAnalysisMc::Draw(TH1* (GFDstarHistsAnalysisMc::*funcPtr) (const char*))
{
  const Bool_t wasBatch = this->SetBatch();
  
  GFHistArray hists;
  for(Int_t i = 0; i < fVariables->GetEntriesFast(); ++i){
    hists.Add((this->*funcPtr)(fVariables->At(i)->GetName()));
  }
  this->SetBatch(wasBatch);
  fHistManager->Clear();
  fHistManager->AddHists(&hists);
  fHistManager->Draw();
}

//____________________________________
 void GFDstarHistsAnalysisMc::DrawRecEffAll()
{
  this->Draw(&GFDstarHistsAnalysisMc::CreateRecEff);
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateTrigEff(const char* var, Int_t TE)
{
  TH1* histTrig = this->GetHistTrig(var, TE);
  TH1* histRec = this->GetHistNonTrig(var);

  return this->CreateTrigEff(histRec, histTrig, TE, var, "D*:");
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateTrigEff(const char *var1, const char *var2,
					   Int_t firstBin2, Int_t lastBin2, Int_t TE)
{
  TString nameAdd("HistTrig");
  if(TE >= 0) nameAdd += TE;

  Bool_t flip = kFALSE;
  TH2 *hTrig2D = static_cast<TH2*>(this->GetHist2Diff(flip, var1, var2, nameAdd, "trig"));
  nameAdd = "HistNonTrig";
  TH2 *hRec2D  = static_cast<TH2*>(this->GetHist2Diff(flip, var1, var2, nameAdd, "trig"));

  if(!hRec2D || !hTrig2D) return NULL;

  TH1 *hTrig =
    GFHistManip::CreateProjection(hTrig2D, Form("%s%s", hTrig2D->GetName(), (flip?"Y":"X")),
				  firstBin2, lastBin2, !flip);
  this->AddOutOfRange(hTrig2D, firstBin2, lastBin2, flip, hTrig);
  TH1 *hRec =
    GFHistManip::CreateProjection(hRec2D, Form("%s%s",hRec2D->GetName(), (flip?"Y":"X")),
				  firstBin2, lastBin2, !flip);
  this->AddOutOfRange(hRec2D, firstBin2, lastBin2, flip, hRec);
  TAxis *axis = (flip ? hTrig2D->GetXaxis() : hTrig2D->GetYaxis());

  if(!GFHistManip::BinsAreSubRange(axis, firstBin2, lastBin2)){
    // could be extended to tolerate under/overflow?
    this->Error("CreateTrigEff", "problem with bin range %d to %d", firstBin2, lastBin2);
    delete hRec2D;
    delete hTrig2D;
    delete hRec;
    delete hTrig;
    return NULL;
  }

  TString titPrefix(Form("D*, %.2f #leq %s < %.2f:", axis->GetBinLowEdge(firstBin2), 
			 GFAxisVariables::GetAxisLabel(var2), axis->GetBinUpEdge(lastBin2)));
  delete hRec2D;
  delete hTrig2D;
  return this->CreateTrigEff(hRec, hTrig, TE, var1, titPrefix);
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateTrigEffDsJet(const char* variable, 
						const char* region, Int_t TE)
{
  // var: Pt/Eta+Ds/Jet or Deta, Dphi, Dpt
  // region: 'Back', 'Forw' or ''
  // TE: 19, 31 or -1 (=default: 19&&31)
  TString histName("dstar1Jet");
  (histName += variable) += region;
  TString trigName = histName + "Trig";
  if(TE >= 0) trigName += TE;
  TH1* hTrig = this->GetHist(trigName, "trigJet");
  TH1* hNonTrig = this->GetHist(histName + "NonTrig", "trigJet");

  return this->CreateTrigEff(hNonTrig, hTrig, TE, variable, Form("D* + 1 jet %s:", region));
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateTrigEffDsDiJet(const char* var, Int_t TE)
{
  // var: Pt/Eta + DsJet/OthJet/Ds, DetaDijet, DphiDijet, PtDijetDijet, CosThDijet, xGamDijet,
  //      DeltaR, DeltaRNoDjet
  // TE: 19, 31 or -1 (=default: 19&&31)
  TString histName("dstarDiJet");
  histName += var;
  TString trigName = histName + "Trig";
  if(TE >= 0) trigName += TE;
  TH1* hTrig = this->GetHist(trigName, "trigDiJet");
  TH1* hNonTrig = this->GetHist(histName + "NonTrig", "trigDiJet");

  return this->CreateTrigEff(hNonTrig, hTrig, TE, var, "D* + dijet");
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateTrigEff(TH1* histRec, TH1* histTrig, Int_t TE,
					   const char* var, const char* titlePrefix)
{
  if(!histTrig || !histRec) {
    this->Error("CreateTrigEff", "Missed one hist for %s/%s, TE %d", var, titlePrefix, TE);
    return NULL;
  }
  
  fHistManager->Clear();

  GFHistManip::MergeOverToUnderFlow(histRec);
  GFHistManip::MergeOverToUnderFlow(histTrig);
  
  TString effName("trigEff"); if(TE >= 0) effName+= TE;
  effName+=var; effName += titlePrefix;
  this->MakeStringName(effName+=fTrigger);

  TGraphAsymmErrors * grPtr = NULL;
  TH1* effHist = GFHistManip::CreateRatioBinomErrors(histTrig, histRec, NULL,NULL,NULL, &grPtr);
  if(effHist){
    effHist->SetName(effName);
  } else {
    this->Warning("CreateTrigEff", "Cannot do correct errors from binomial!");
    effHist = static_cast<TH1*>(histTrig->Clone(effName));
    effHist->Reset();
    effHist->Divide(histTrig, histRec, 1., 1., "B");
  }

  histTrig->SetName(effName+"NumTrig");
  fHistManager->AddHist(histTrig);
  histRec->SetName(effName+"NumRec");
  fHistManager->AddHist(histRec);

  TString title(Form("%s #epsilon_{trig}, TE ", titlePrefix));
  if(TE >= 0) title += TE; else title+="19+31";
  effHist->SetTitle(title);
  effHist->SetYTitle("#epsilon_{trig}");

  fHistManager->AddHist(effHist,1);

  fHistManager->SetHistsXTitle(GFAxisVariables::GetAxisLabel(var));
  if(grPtr) {
    fHistManager->AddObject(grPtr, 1, 0, "P");
    effHist->SetOption("HIST");
    effHist->TAttLine::Copy(*grPtr);
    effHist->TAttMarker::Copy(*grPtr);
  } else {
    effHist->SetOption("E");
  }
  this->TmpGraph(grPtr);
  fHistManager->Draw();

  GFHistManip::PrintTable(effHist, Form("%s from %s in bins of %s", title.Data(), 
					fPeriod.Data(), var));
  return effHist;
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateDetAcceptance(const char* variable)
{
  //pt,eta,phi,wGammaP,...
  TH1* hGen = this->GetHistGenWeight(variable);
  TH1* hDetAcc = this->GetHist(Form("%sHistGenWeightDetAcc", variable),"gen");
  
  if(!hGen || ! hDetAcc) return NULL;
  fHistManager->Clear();
  TString name("detAcc");
  this->MakeStringName(name+=fTrigger);
  TGraphAsymmErrors* graphPtr = NULL;
  TH1* hResult = GFHistManip::CreateRatioBinomErrors(hDetAcc, hGen, NULL, NULL, NULL,&graphPtr);
  if(hResult){
    hResult->SetName(name);
  } else {
    this->Warning("CreateDetAcceptance", "Cannot do correct errors from binomial!");
    hResult = GFHistManip::CreateAnalog(hDetAcc, name, kFALSE, kFALSE);
    hResult->Divide(hDetAcc, hGen, 1., 1., "B");
  }

  hResult->SetTitle("detector acceptance");
  hResult->SetYTitle("A_{det}");
  fHistManager->Clear();

  fHistManager->AddHist(hGen);
  fHistManager->AddHist(hDetAcc);
  fHistManager->AddHist(hResult, 1);
  fHistManager->SetHistsXTitle(GFAxisVariables::GetAxisLabel(variable));
  if(graphPtr) {
    for(Int_t i = 1; i <= hResult->GetNbinsX(); ++i){
      hResult->SetBinError(i, 1.e-10);
    }
    fHistManager->AddObject(graphPtr, 1, 0, "P");
    hResult->TAttLine::Copy(*graphPtr);
    hResult->TAttMarker::Copy(*graphPtr);
  }
  this->TmpGraph(graphPtr);

  fHistManager->Draw();

  return hResult;
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateHadCorr(const char* var, const char* forwBack)
{
  // var is PtDs, PtJet,EtaJet, EtaDs, Deta, Dphi etc.
  TH1* histGen = this->GetHist(Form("dstar1Jet%s%sGen", var, forwBack), "partonJet");
  TH1* histParton = this->GetHist(Form("dstar1Jet%s%sParton", var, forwBack), "partonJet");
  TString titlePrefix("D*+jet ");
  titlePrefix+=forwBack;
  if(!histGen || !histParton) {
    this->Error("CreateHadCorr", "Missed one hist for %s/%s", var, titlePrefix.Data());
    return NULL;
  }

  return this->CreateHadCorr(histGen, histParton, var, titlePrefix.Data());
}


//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateHadCorrDiJet(const char* var)
{

  TH1* histGen = this->GetHist(Form("dstarDiJet%sGen", var), "partonDiJet");
  TH1* histParton = this->GetHist(Form("dstarDiJet%sParton", var), "partonDiJet");
  const char *titlePrefix = "D* dijet ";
  if(!histGen || !histParton) {
    this->Error("CreateHadCorrDiJet", "Missed one hist for %s/%s", var, titlePrefix);
    return NULL;
  }

  return this->CreateHadCorr(histGen, histParton, var, titlePrefix);
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateHadCorr(TH1 *histGen, TH1 *histParton, 
					   const char *var, const char *titlePrefix)
{
  
  if (!histGen || !histParton) return NULL;

  fHistManager->Clear();
  
  TString hName("hadCorr");
  hName+=var; hName += titlePrefix;
  this->MakeStringName(hName+=fTrigger);

  this->Info("CreateHadCorr", "Do not do correct errors from binomial! (?)");
  TH1* corrHist = static_cast<TH1*>(histGen->Clone(hName));
  corrHist->Reset();
  corrHist->Divide(histGen, histParton, 1., 1., "B");

  histGen->SetName(hName+"NumGen");
  fHistManager->AddHist(histGen);
  histParton->SetName(hName+"NumParton");
  fHistManager->AddHist(histParton);

  TString title(Form("%s had. correction ", titlePrefix));
  corrHist->SetTitle(title);
  corrHist->SetYTitle("C_{had}");

  fHistManager->AddHist(corrHist,1);

  fHistManager->SetHistsXTitle(GFAxisVariables::GetAxisLabel(var));
  corrHist->SetOption("HIST");

  fHistManager->Draw();

  GFHistManip::PrintTable(corrHist, Form("%s from %s in bins of %s", title.Data(), 
					 fPeriod.Data(), var));
  return corrHist;
}


//____________________________________
 TArrayD GFDstarHistsAnalysisMc::TotalHadCorr(const char *forwBack)
{
  // for D* jet
  TH1* histGen = this->GetHist(Form("dstar1JetEtaDs%sGen", forwBack), "partonJet");
  TH1* histParton = this->GetHist(Form("dstar1JetEtaDs%sParton", forwBack), "partonJet");
  if(!histGen || !histParton) {
    this->Error("TotalHadCorr", "Missed one hist");
    return TArrayD(0);
  }

  const Double_t nGen = histGen->Integral();
  const Double_t nParton = histParton->Integral();
  Stat_t s[11];
  histParton->GetStats(s); // now s[1] is sum of squares of weights, s[0] is sum of weights
  const Double_t nEffParton = s[0]*s[0]/s[1];
  const Double_t nFillParton = histParton->GetEntries();

  delete histGen;
  delete histParton;
  
  TArrayD hadCorr = GFMath::Efficiency(nGen, nParton, nEffParton, nFillParton);

  TString text(Form("D*+o. jet %s had. corrections for S%d from %s", 
		    forwBack, fTrigger, fPeriod.Data()));
  this->PrintPlusMinus(hadCorr, text);

  return hadCorr;
}


//____________________________________
 TArrayD GFDstarHistsAnalysisMc::TotalHadCorrDiJet()
{
  // for D* dijet
  TH1* histGen = this->GetHist("dstarDiJetEtaDJetGen", "partonDiJet");
  TH1* histParton = this->GetHist("dstarDiJetEtaDsParton", "partonDiJet");
  if(!histGen || !histParton) {
    this->Error("TotalHadCorrDiJet", "Missed one hist");
    return TArrayD(0);
  }

  const Double_t nGen = histGen->Integral();
  const Double_t nParton = histParton->Integral();
  Stat_t s[11];
  histParton->GetStats(s); // now s[1] is sum of squares of weights, s[0] is sum of weights
  const Double_t nEffParton = s[0]*s[0]/s[1];
  const Double_t nFillParton = histParton->GetEntries();

  delete histGen;
  delete histParton;
  
  TArrayD hadCorr = GFMath::Efficiency(nGen, nParton, nEffParton, nFillParton);

  TString text(Form("D* dijet had. corrections for S%d from %s", fTrigger, fPeriod.Data()));
  this->PrintPlusMinus(hadCorr, text);

  return hadCorr;
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateRecEffDs1Jet(const char* var, 
						const char* forwBack, UInt_t dmRefFlag)
{
  // rec eff for D* + 1 jet
  // var: Pt/Eta+Ds/Jet (e.g. PtDs), Dpt, Deta, Dphi
  // forwBack: 'Forw', 'Back', '' (not NULL) 
  // if(dmRefFlag != 0): First fit a reference hist with all params free, then fit with fFitMode
  // >=2 [default]: D* inclusive
  // 1          : D* + 1 jet inclusive

  TH2* hRec2D =static_cast<TH2*>(this->GetHist(Form("dstar1Jet%s%s",var,forwBack),"dstarJets"));
  TH1* histGen = this->GetHist(Form("dstar1Jet%s%sGenWeight", var, forwBack), "genJet");

  TH1* hRecRef = this->GetHistDs1JetRef(forwBack, dmRefFlag);

  if(hRec2D && histGen && hRecRef){
    return this->CreateRecEff(var, hRec2D, hRecRef, histGen, 
			      TString(Form("D*%s + 1 jet:", forwBack)));
  } else {
    this->Error("CreateRecEffDs1Jet", "Did not find one of the hists!");
    return NULL;
  }
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateRecEffDsDiJet(const char* var, 
						UInt_t dmRefFlag)
{

  // rec eff for D* + dijet
  // var: Pt/Eta + DsJet/OthJet/Ds, Deta, Dphi, PtDijet, CosTh, xGam,
  //      DeltaR, DeltaRNoDjet
  // First fit a reference hist with all params free, then fit with fFitMode
  // >=2 [default]: D* inclusive
  // 0,1          : D* + dijet

  TH2* hRec2D =(TH2*) this->GetHist(TString(Form("dstarDiJet%s", var)),"dsDiJets");
  TH1* histGen = this->GetHist(TString(Form("dstarDiJet%sGenWeight", var)),"genDiJet");

  TH1* hRecRef = NULL;
  switch(dmRefFlag){
  case 0:
  case 1:
    hRecRef = this->GetHist("dstarDiJet" ,"dsDiJets");
    break;
  default: 
    this->Warning("CreateRecEffDsDiJet", "Don't know dmRefFlag %d, use 2",dmRefFlag);
  case 2:
    hRecRef = this->GetHistPtCut();
  }

  if(hRec2D && histGen && hRecRef){
    return this->CreateRecEff(var, hRec2D, hRecRef, histGen, "D* + dijet:");
  } else {
    this->Error("CreateRecEffDsDiJet", "Did not find one of the hists!");
    return NULL;
  }
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateRecEff(const char* variable)
{
  // rec eff for D* hists: pt, eta, phi, wGammaP, ...
  TH2* histRec2D = this->GetHist2D(variable);
  TH1* histGen = this->GetHistGenWeight(variable);
  TH1* histRef = this->GetHistPtCut();
  if(histRec2D && histGen){
    return this->CreateRecEff(variable, histRec2D, histRef, histGen, "D*:");
  } else {
    this->Error("CreateRecEff", "Did not find both hists!");
    return NULL;
  }
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateRecEff(const char* var1, const char* var2,
					  Int_t firstBin2, Int_t lastBin2)
{
  // 

  // get rec hist:
  Bool_t flip = kFALSE;
  TH3 *hRec3D = static_cast<TH3*>(this->GetHist2Diff(flip, var1, var2, "Hist2D", NULL));
  if(!hRec3D) return NULL;
  TAxis *axis = (flip ? hRec3D->GetYaxis() : hRec3D->GetZaxis());
  if(!GFHistManip::BinsAreSubRange(axis, firstBin2, lastBin2)){
    // could be extended to tolerate under/overflow?
    this->Error("CreateRecEff", "problem with bin range %d to %d", firstBin2, lastBin2);
    delete hRec3D;
    return NULL;
  }
  axis->SetRange(firstBin2, lastBin2);
  TH2 *hRec2D = static_cast<TH2*>(hRec3D->Project3D(flip?"zxe":"yxe"));
  this->AddOutOfRange(hRec3D, firstBin2, lastBin2, flip, hRec2D);

  hRec2D->SetTitle(Form("%.2f #leq %s < %.2f", axis->GetBinLowEdge(firstBin2), 
			GFAxisVariables::GetAxisLabel(var2), axis->GetBinUpEdge(lastBin2)));
  delete hRec3D;

  // get gen hist:
  TH2 *hGen2D = static_cast<TH2*>(this->GetHist2Diff(flip, var1, var2, "HistGenWeight", "gen"));
  if(!hGen2D) return NULL;
  TH1 *hGen =
    GFHistManip::CreateProjection(hGen2D, Form("%s%s",hGen2D->GetName(), (flip?"Y":"X")),
				  firstBin2, lastBin2, !flip);
  this->AddOutOfRange(hGen2D, firstBin2, lastBin2, flip, hGen);

  delete hGen2D;
  // ref hist for fit:
  TH1* hRef = this->GetHistPtCut();
  
  TString titPrefix(Form("D*, %s:", hRec2D->GetTitle()));
  return this->CreateRecEff(var1, hRec2D, hRef, hGen, titPrefix);
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateRecEff(const char* var, TH2* hRec2D, TH1 *hRecRef,
					  TH1* hGen, const char* titlePrefix)
{
  // hGen:   hist with generator level distribution
  // hRec2D: hist with reconstruction level distribution on y-axis and delta(m) on x-axis
  // var:    name of variable

  fHistManager->Clear();

  TString name("recEff");
  name+=var; name += titlePrefix;
  this->MakeStringName(name+=fTrigger);

  // N(gen)-hist:
//   this->DeleteObject(name+"NumGen");
//   TH1* histGen =  static_cast<TH1*>(histGen->Clone(name+"NumGen"));
  GFHistManip::MergeOverToUnderFlow(hGen);
  hGen->SetName(name+"NumGen");
  fHistManager->AddHist(hGen);

  // N(rec)-hist:
  TH1* recHist = fMultiDmFitter->Fit(hRec2D, hRecRef);//NULL, NULL);
//   this->DeleteObject(name+"NumRec");
  recHist->SetName(name+"NumRec");
  fHistManager->AddHists(fMultiDmFitter->GetHists1D(), 1);
  delete fMultiDmFitter->GetHistSignalBack(); // little (!) clean up...

  fHistManager->AddHist(recHist);

  TH1* effHist = GFHistManip::CreateRatioBinomErrors(recHist, hGen);
  if(effHist){
    effHist->SetName(name);
  } else {
    this->Warning("CreateRecEff", "Cannot do correct errors from binomial!");
    effHist = static_cast<TH1*>(recHist->Clone(name));
    effHist->Reset();
    effHist->Divide(recHist, hGen, 1., 1., "B");
  }

  effHist->SetTitle(Form("%s #epsilon'_{rec}", titlePrefix));
  effHist->SetYTitle("#epsilon'_{rec}");
  effHist->SetOption("e");
  fHistManager->AddHist(effHist,2);

  fHistManager->SetHistsXTitle(GFAxisVariables::GetAxisLabel(var));
  fHistManager->Draw();

  GFHistManip::PrintTable(effHist, Form("%s rec eff from %s in bins of %s", 
					titlePrefix, fPeriod.Data(), var));
  return effHist;
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateAccept(const char* var)
{
  TH1* histGenAcc = this->GetHistGenWeight(var); // contains also acceptance as weight
  TH1* histGen    = this->GetHistGenNoAccWeight(var); // all weight besides acceptance

  return this->CreateAccept(histGenAcc, histGen, var, "D*:");
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateAccept(const char* var1, const char* var2,
					  Int_t firstBin2, Int_t lastBin2)
{
  // pt, eta, phi, wGammaP, etc. for D*

  Bool_t flip = kFALSE;
  TH2 *hGenAcc2D = static_cast<TH2*>(this->GetHist2Diff(flip, var1,var2,"HistGenWeight","gen"));
  TH2 *hGen2D =static_cast<TH2*>(this->GetHist2Diff(flip,var1,var2,"HistGenNoAccWeight","gen"));

  if(!hGen2D || !hGenAcc2D) return NULL;

  TH1 *hGen =
    GFHistManip::CreateProjection(hGen2D, Form("%s%s", hGen2D->GetName(), (flip?"Y":"X")),
				  firstBin2, lastBin2, !flip);
  this->AddOutOfRange(hGen2D, firstBin2, lastBin2, flip, hGen);
  TH1 *hGenAcc =
    GFHistManip::CreateProjection(hGenAcc2D, Form("%s%s",hGenAcc2D->GetName(), (flip?"Y":"X")),
				  firstBin2, lastBin2, !flip);
  this->AddOutOfRange(hGenAcc2D, firstBin2, lastBin2, flip, hGenAcc);
  TAxis *axis = (flip ? hGen2D->GetXaxis() : hGen2D->GetYaxis());

  if(!GFHistManip::BinsAreSubRange(axis, firstBin2, lastBin2)){
    // could be extended to tolerate under/overflow?
    this->Error("CreateAccept", "problem with bin range %d to %d", firstBin2, lastBin2);
    delete hGen2D;
    delete hGenAcc2D;
    delete hGen;
    delete hGenAcc;
    return NULL;
  }

  TString titPrefix(Form("D*, %.2f #leq %s < %.2f:", axis->GetBinLowEdge(firstBin2), 
			 GFAxisVariables::GetAxisLabel(var2), axis->GetBinUpEdge(lastBin2)));
  delete hGen2D;
  delete hGenAcc2D;
  
  return this->CreateAccept(hGenAcc, hGen, var1, titPrefix);
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateAcceptDsJet(const char* var, const char* regio)
{
  // var: Pt/Eta+Ds/Jet or Deta, Dphi, Dpt
  // region: 'Back', 'Forw' or ''
  TString histName("dstar1Jet");
  (histName += var) += regio;

  TH1* histGenAcc = this->GetHist(histName + "GenWeight", "genJet");// also acceptance as weight
  TH1* histGen    = this->GetHist(histName + "GenNoAccWeight", "genJet");// weight. w/o accept.

  return this->CreateAccept(histGenAcc, histGen, var, Form("D*%s + 1 jet:",regio));
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateAcceptDsDiJet(const char* var)
{
  // var: Pt/Eta + DsJet/OthJet/Ds, Deta, Dphi, PtDijet, CosTh, xGam,
  //      DeltaR, DeltaRNoDjet
  TString histName("dstarDiJet");
  histName += var;

  TH1* histGenAcc = this->GetHist(histName + "GenWeight", "genDiJet");// also acceptance as weight
  TH1* histGen    = this->GetHist(histName + "GenNoAccWeight", "genDiJet");// weight. w/o accept.

  return this->CreateAccept(histGenAcc, histGen, var, "D* + dijet:");
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateAccept(TH1* histGenAcc, TH1* histGen, 
					  const char* variable, const char* titlePrefix)
{
  // histGenAcc: generator level, weighted with tagger acceptance
  // histgen: generator level

  if(!histGenAcc || !histGen) {
    this->Error("CreateAccept", "Missed one hist for %s/%s", variable, titlePrefix);
    return NULL;
  }

  GFHistManip::MergeOverToUnderFlow(histGenAcc);
  GFHistManip::MergeOverToUnderFlow(histGen);

  fHistManager->Clear();
  TString name("acceptMc");
  name += variable; name += titlePrefix;
  this->MakeStringName(name+=fTrigger);
//   this->DeleteObject(name);

  TGraphAsymmErrors * grPtr = NULL;
  TH1* histAccept = GFHistManip::CreateRatioBinomErrors(histGenAcc, histGen, NULL, NULL, NULL,
							&grPtr);
  if(histAccept){
    histAccept->SetName(name);
  } else {
    this->Warning("CreateAccept", "Cannot do correct errors from binomial!");
    histAccept = static_cast<TH1*>(histGen->Clone(name));
    histAccept->Reset();
    histAccept->Divide(histGenAcc, histGen, 1., 1., "B");
  }
  TString title(Form("%s e-tag acceptance (MC)", titlePrefix));
  histAccept->SetTitle(title);
  histAccept->SetYTitle(Form("A_{%s}", this->GetSampleName()));
    
  fHistManager->AddHist(histGenAcc);
  fHistManager->AddHist(histGen);
  fHistManager->AddHist(histAccept,1);


  if(grPtr) {
    fHistManager->AddObject(grPtr, 1, 0, "P");
    histAccept->SetOption("HIST");
    histAccept->TAttLine::Copy(*grPtr);
    histAccept->TAttMarker::Copy(*grPtr);
  } else {
    histAccept->SetOption("E");
  }
  this->TmpGraph(grPtr);

    
  fHistManager->SetHistsXTitle(GFAxisVariables::GetAxisLabel(variable));
//   fHistManager->AddHistsOption("e");
  fHistManager->Draw();

  GFHistManip::PrintTable(histAccept, Form("%s from %s in bins of %s", title.Data(), 
					   fPeriod.Data(), variable)); 
  return histAccept;
}

// before 'per lumi'
// //____________________________________
// TH1* GFDstarHistsAnalysisMc::CreateAcceptOld(const char* variable)//, Bool_t averagePeriods)
// {
//   // outdated!
//   fHistNumDstar = this->GetHistAccept(variable);
//   fHist1D = this->GetHistTrig(variable);

//   fHistManager->Clear();
//   TH1* histAccept = this->CombinedAccept(variable);
//   if(!histAccept){
//     cout << "Using NON lumi averaged acceptance from MC" << endl;
//     TString name("acceptMc");
//     TString axisLabel(GFAxisVariables::GetAxisLabel(variable));
//     this->MakeStringName((name+=axisLabel)+=fTrigger);
    
// //     this->DeleteObject(name);
//     TArrayD bins(fHist1D->GetNbinsX()+1);
//     for(Int_t i = 0; i < bins.GetSize(); ++i){
//       bins[i] = fHist1D->GetBinLowEdge(i+1); //FIXME: correct?
//     }
//     histAccept = new TH1F(name, "e-tagger acceptance (MC)", 
// 			       bins.GetSize()-1, bins.GetArray());
    
//     // if there are other values than'1.', 
//     // ==>  look in html-docu of TH1::Divide(...) in 3.02_07 and 3.03_? 
//     histAccept->Divide(fHistNumDstar, fHist1D, 1.,1.,"B");
//     //   histAccept->SetTitle("e-tagger acceptance (MC)");
    
//     fHistManager->AddHist(fHistNumDstar);
//     fHistManager->AddHist(fHist1D);
//     fHistManager->AddHist(histAccept,1);
    
//     fHistManager->SetHistsXTitle(axisLabel);
//     fHistManager->AddHistsOption("e");
//     fHistManager->Draw();
// //     fHistManager->WriteHistos(fFileOut);
//   }
//   this->SetHistPointerNull();
//   TString title("tagger acceptance from ");
//   (title += fPeriod) += " MC in bins of ";
//   GFHistManip::PrintTable(histAccept, title += variable); 
//   return histAccept;
// }

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CombinedAccept(const char* variable)
{
  if(fTrigger == 83) {
    this->Warning("CombinedAccept","No combining foreseen for ET33!");
    return NULL;
  } else {
    this->Error("CombinedAccept", "Not implemented anymore after 'per lumi'!");
    return NULL;
  }

}


//____________________________________
 TH1* GFDstarHistsAnalysisMc::XSectionTotHist(Bool_t draw, Bool_t acceptanceWeighted)
{
  if(draw) fHistManager->Clear();

  TArrayD xSectAndError = this->XSectionTot(acceptanceWeighted);

  TString name("sigmaTotMc");
  TString title("#sigma_{vis} ep #rightarrow D* + X [nb] MC "); // pb

  TH1* hTemp = this->GetHist2D("wGammaP");
  TH1* histX = new TH1F((name+=fTrigger)+=this->GetPeriod(), 
			title+=this->GetSampleName(), 
			1, hTemp->GetYaxis()->GetXmin(), hTemp->GetYaxis()->GetXmax());

  histX->SetBinContent(1, xSectAndError[0]);
  histX->SetBinError(1, xSectAndError[1]);
  histX->SetXTitle("W_{#gamma p} [GeV]");

  if(draw) {
    fHistManager->AddHist(histX);
    fHistManager->Draw();
  }
  return histX;
}

//____________________________________
 TArrayD GFDstarHistsAnalysisMc::XSectionTot(Bool_t acceptanceWeighted)
{
  
  const Double_t numDstar = acceptanceWeighted ? this->GetHistGenAccept("eta")->Integral()
                                                : this->GetHistGen("eta")->Integral();

  TArrayD xSectAndError(2);
  xSectAndError[0] = kBranchRatioMc ? numDstar * kXSecCorrFactor/kBranchRatioMc : 0.;
  //WRONG:
//xSectAndError[1] = numDstar ? xSectAndError[0] / TMath::Sqrt(numDstar) : 0;

  this->PrintPlusMinus(xSectAndError, "Total X-Section of MC " + fPeriod);
  return xSectAndError;
}

//____________________________________
 TArrayD GFDstarHistsAnalysisMc::XSectionTotDsJet(const char* forwBack)
{
  // total cross section for D* + 1 jet
  // forwback: Forw, Back, ''
  const Double_t nDstar = this->GetHist(Form("dstar1Jet%sGen",forwBack), "genJet")->Integral();

  TArrayD xSectAndError(2);
  xSectAndError[0] = kBranchRatioMc ? nDstar * kXSecCorrFactor/kBranchRatioMc : 0.;
  //WRONG:
//   xSectAndError[1] = nDstar ? xSectAndError[0] / TMath::Sqrt(nDstar) : 0.;

  this->PrintPlusMinus(xSectAndError, Form("D* + 1 jet %s: X-Section of %s (S%d)", 
					   forwBack, fPeriod.Data(), fTrigger));
  return xSectAndError;
}

//____________________________________
 TArrayD GFDstarHistsAnalysisMc::XSectionTotDsDiJet()
{
  // total cross section for D* dijet
  const Double_t nDstar = this->GetHist("dstarDiJetGen", "genDiJet")->Integral();

  TArrayD xSectAndError(2);
  xSectAndError[0] = kBranchRatioMc ? nDstar * kXSecCorrFactor/kBranchRatioMc : 0.;
  //WRONG:
//   xSectAndError[1] = nDstar ? xSectAndError[0] / TMath::Sqrt(nDstar) : 0.;

  this->PrintPlusMinus(xSectAndError, Form("D* dijet: X-Section of %s (S%d)", 
					   fPeriod.Data(), fTrigger));
  return xSectAndError;
}


//____________________________________
 TH1* GFDstarHistsAnalysisMc::XSectionPt(Bool_t acceptanceWeighted)
{
//   fHist1D = acceptanceWeighted ? this->GetHistGenAccept("pt") : this->GetHistGen("pt");
  TH1* hist = this->XSection("pt", acceptanceWeighted);
  fHistManager->SetLogY();
  return hist;
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::XSectionEta(Bool_t acceptanceWeighted)
{
//   fHist1D = acceptanceWeighted ? this->GetHistGenAccept("eta") :this->GetHistGen("eta");
  return this->XSection("eta", acceptanceWeighted);
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::XSectionPhi(Bool_t acceptanceWeighted)
{
//   fHist1D = acceptanceWeighted ? this->GetHistGenAccept("phi") : this->GetHistGen("phi");
  return this->XSection("phi", acceptanceWeighted);
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::XSectionWgammaP(Bool_t acceptanceWeighted)
{
  return this->XSection("wGammaP", acceptanceWeighted);
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::XSection(const char* var, Bool_t acceptanceWeighted)
{
  TH1* h  = acceptanceWeighted ? this->GetHistGenAccept(var) 
                               : this->GetHistGen(var);
  return this->XSection(h, var, (acceptanceWeighted ? "D* (acc. weight)" : "D*"));
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::XSection(TH1* h, const char* var,const char* titPrefix)
{
  if(!h){
    this->Error("XSection", "No histogram");
    return NULL;
  }
  fHistManager->Clear();
  const char* varName = GFAxisVariables::GetAxisLabel(var);

  TString histName("mcXsection");
  (histName += var)+= fTrigger;
  this->MakeStringName(histName);

  TH1* histX = static_cast<TH1*>(h->Clone(histName));
  GFHistManip::CorrectForBinWidth(histX);
  histX->Scale(kXSecCorrFactor/kBranchRatioMc);

  TString title(Form("d#sigma_{vis}(ep#rightarrow %s + X) / d", titPrefix));
  (title += varName)+=" [nb]"; // pb
  title.ReplaceAll("[GeV]", 0);
  histX->SetTitle(title);
  histX->SetXTitle(varName);

  fHistManager->AddHist(histX);
  fHistManager->Draw();

  //  this->SetHistPointerNull();
//   title  = "X-section from  ";
//   (title += fPeriod) += " MC in bins of ";
  GFHistManip::PrintTable(histX, Form("%s from %s (S%d)",title.Data(),fPeriod.Data(),fTrigger));

  return histX;
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::XSection(const char* var1, const char* var2, 
				      Int_t firstBin2, Int_t lastBin2)
{
  Bool_t flip = kFALSE;
  TH2* h2D = static_cast<TH2*>(this->GetHist2Diff(flip, var1, var2, "HistGen", "gen"));
  if(!h2D) return NULL;

  TAxis* axis = (flip ? h2D->GetXaxis() : h2D->GetYaxis());
  if(!GFHistManip::BinsAreSubRange(axis, firstBin2, lastBin2)){
    // could be extended to tolerate under/overflow?
    this->Error("XSection", "problem with bin range %d to %d", firstBin2, lastBin2);
    delete h2D;
    return NULL;
  }
  TH1 *h = GFHistManip::CreateProjection(h2D, Form("%s%s", h2D->GetName(), (flip?"Y":"X")),
					 firstBin2, lastBin2, !flip);

  TH1* hXSec = this->XSection(h, var1, "D*");
  hXSec->SetTitle(Form("%s (%.2f #leq %s < %.2f)", hXSec->GetTitle(), 
		       axis->GetBinLowEdge(firstBin2), 
		       GFAxisVariables::GetAxisLabel(var2), axis->GetBinUpEdge(lastBin2)));
  delete h2D;
  return hXSec;
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::XSectionDsJet(const char* var, const char* forwBack)
{
  TH1* h  = this->GetHist(Form("dstar1Jet%s%sGen", var, forwBack), "genJet");
  return this->XSection(h, var, Form("D* + 1 jet %s", forwBack));
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::XSectionDsDiJet(const char* var)
{
  TH1* h  = this->GetHist(Form("dstarDiJet%sGen", var), "genDiJet");
  return this->XSection(h, var, "D* dijet");
}


//____________________________________
 TH1* GFDstarHistsAnalysisMc::XSectionRec(const char* variable)
{
  //pt,eta,phi,wGammaP,y,nJet,eta/ptJet1/2,xGamma
  TH1* hNDstar = this->DrawDm(variable);
  TH1* hRecEff = this->CreateRecEff(variable);
  //  TH1* hTrigEff = this->CreateTrigEff(variable);//FIXME: put in as DrawDm now requires TEs! 
  TH1* hAccept = this->CreateAccept(variable);
  TH1* hXSecGen = this->XSection(variable);

  return this->XSectionRecInt(hNDstar, hRecEff, //hTrigEff, 
			      hAccept, hXSecGen, 
			      this->GetHistGenNoAccWeight(variable), 
			      this->GetHistGen(variable));
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::XSectionDsJetRec(const char* var, const char* foBack)
{
  TH1* hNDstar = this->DrawDmDs1Jet(var, foBack);
  TH1* hRecEff = this->CreateRecEffDs1Jet(var, foBack);
  //  TH1* hTrigEff = this->CreateTrigEff(variable);//FIXME: put in as DrawDm now requires TEs!
  TH1* hAccept = this->CreateAcceptDsJet(var, foBack);
  TH1* hXSecGen = this->XSectionDsJet(var, foBack);

  TString histName(Form("dstar1Jet%s%s", var, foBack));
  return this->XSectionRecInt(hNDstar, hRecEff, //hTrigEff, 
			      hAccept, hXSecGen, 
			      this->GetHist(histName + "GenNoAccWeight", "genJet"),
			      this->GetHist(histName + "Gen", "genJet"));
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::XSectionRecInt(TH1* hNDstar, TH1* hRecEff,//TH1* hTrigEff,
					    TH1* hAccept, TH1* hXSecGen, 
					    TH1* hGenNoAccW, TH1* hGen)
{
  // hNDstar already corrected for binwidth...!

  //  GFHistManip::CorrectForBinWidth(hNDstar);

  this->Error("XSectionRecInt", "trig eff not included!");

  TH1* hXSec = static_cast<TH1*>(hNDstar->Clone("xsec"));
  hXSec->Reset();
  hXSec->Divide(hNDstar, hRecEff, kXSecCorrFactor, kBranchRatioMc); // after 'per lumi'
//   hXSec->Divide(hTrigEff);//FIXME: put in as DrawDm now requires TEs!
  hXSec->Divide(hAccept);
  hGenNoAccW->Divide(hGen); // this is now the average inverse of the prescale-weight
  hXSec->Divide(hGenNoAccW); // correct for that: It's not in hNDstar!

  for(Int_t i = 1; i <= hNDstar->GetNbinsX(); ++i){
    const Double_t cont = hNDstar->GetBinContent(i);
    if(cont){
      hXSec->SetBinError(i, hNDstar->GetBinError(i)/cont*hXSec->GetBinContent(i));
    } else {
      this->Warning("XSectionRecInt", "bin %d is empty, no relative error", i);
    }
  }

  fHistManager->Clear();
  fHistManager->AddHist(hNDstar);
  fHistManager->AddHist(hRecEff);
  fHistManager->AddHist(hAccept);
  fHistManager->AddHist(hGenNoAccW);
  hGenNoAccW->SetTitle("weight not in N_{gen} (i.e. mean inverse prescale, pt(c),...)");

  fHistManager->AddHist(hXSec, 1, "rec", "lp");
  fHistManager->AddHistSame(hXSecGen, 1, 0, "gen", "lp");
  fHistManager->AddLegend(1, 0, "Cross Section");

  fHistManager->Draw();

  return hXSec;
}


//____________________________________
 void GFDstarHistsAnalysisMc::DrawEflow(const char* region)
{
  // Forw, ForwBar, Barrel, Spacal, "" (sum) or NULL: all 3+1 regions;
  
  TObjArray regions;
  if(region){
    regions.Add(new TObjString(region));
  } else {
    regions.Add(new TObjString("Barrel"));
    regions.Add(new TObjString("ForwBar"));
    regions.Add(new TObjString("Forw"));
    regions.Add(new TObjString("Spacal"));
    regions.Add(new TObjString(""));
  }
  const char directory[] = "Eflow";
  
  fHistManager->Clear();

  for(Int_t i = 0; i < regions.GetEntriesFast(); ++i){
    TString name("ErecEgen");
    name += regions.At(i)->GetName();
    TH2* hist2D = this->GetHist2D(name, directory);
    hist2D->SetXTitle("E_{rec} [GeV]");
    hist2D->SetYTitle("E_{gen} [GeV]");
    hist2D->SetOption("BOX");
    TH1* histRatio = this->GetHist(name + "Hist", directory);
//     histRatio->Fit("gaus", "0","", histRatio->GetMean() - .35, histRatio->GetMean() + .25);
//     histRatio->GetFunction("gaus")->ResetBit(1 << 9);
    histRatio->SetXTitle("E_{rec}/E_{gen}");

    TH1* histRelDiff = this->GetHist(name + "relHist", directory);
//     histRelDiff->Fit("gaus","0","", histRelDiff->GetMean() - 0.3, histRelDiff->GetMean() + 0.2);
//     histRelDiff->GetFunction("gaus")->ResetBit(1 << 9);
    TH1* histDiff = this->GetHist(name + "diffHist", directory);

    fHistManager->AddHist(hist2D, i/2);
    fHistManager->AddHist(histRatio, i/2);
    fHistManager->AddHist(histDiff, i/2);
    fHistManager->AddHist(histRelDiff, i/2);
  }
  fHistManager->SetNumHistsXY(4,2);

  for(Int_t i = 0; i < fHistManager->GetNumLayers(); ++i){
    const Int_t nHistsLayer = fHistManager->GetNumHistsOf(i);
    //    for(Int_t j = 1; j < (region ? 4 : 7); j += 3){
    for(Int_t j = 1; j < nHistsLayer; j += 4){
      TH1* hist = fHistManager->GetHistsOf(i, j-1)->At(0);
      if(hist) {
	TLine *line = new TLine(hist->GetXaxis()->GetXmin(), hist->GetYaxis()->GetXmin(),
				hist->GetXaxis()->GetXmax(), hist->GetYaxis()->GetXmax());
	line->SetLineWidth(3);
	line->SetLineColor(kRed);
	fHistManager->AddObject(line, i, j-1);
      }
      hist = fHistManager->GetHistsOf(i, j)->At(0);
      if(hist){
	TLine *line = new TLine(1., 0., 1., hist->GetMaximum()*1.05);
	line->SetLineWidth(3);
	line->SetLineColor(kRed);
	fHistManager->AddObject(line, i, j);

      }
      hist = fHistManager->GetHistsOf(i, j+1)->At(0);
      if(hist){
	TLine *line = new TLine(0., 0., 0., hist->GetMaximum()*1.05);
	line->SetLineWidth(3);
	line->SetLineColor(kRed);
	fHistManager->AddObject(line, i, j+1);
      }
    }
  }
  fHistManager->Draw();
  regions.Delete();
}

//____________________________________
 void GFDstarHistsAnalysisMc::DrawJetQuark(const char* jet)
{
  // jet1, jet2
  fHistManager->Clear();
  TObjArray vars;
  vars.Add(new TObjString("E"));
  vars.Add(new TObjString("Et"));
  vars.Add(new TObjString("phi"));
  vars.Add(new TObjString("eta"));
  TObjArray endings;
  endings.Add(new TObjString("cVs"));
  endings.Add(new TObjString("cDiff"));
  endings.Add(new TObjString("cOver"));

  for(Int_t iVar = 0; iVar < vars.GetEntriesFast(); ++iVar){
    for(Int_t iEnd = 0; iEnd < endings.GetEntriesFast(); ++iEnd){
      if(iVar < 2 || iEnd < 2){
	TH1* hist = this->GetHist(Form("%s%s%s%s", vars[iVar]->GetName(), "jet1",
				       vars[iVar]->GetName(), endings[iEnd]->GetName()), 
				  "jetquark");
	fHistManager->AddHist(hist, iVar);
	hist = this->GetHist(Form("%s%s%s%s", vars[iVar]->GetName(), "jet2",
				  vars[iVar]->GetName(), endings[iEnd]->GetName()), "jetquark");
	fHistManager->AddHist(hist, iVar);
      }
    }
  }
  fHistManager->AddHist(this->GetHist("distjet1Charm", "jetquark"), vars.GetEntriesFast());
  fHistManager->AddHist(this->GetHist("distjet2Charm", "jetquark"), vars.GetEntriesFast());
  fHistManager->AddHist(this->GetHist("sumDistsJetCharm", "jetquark"), vars.GetEntriesFast());
  fHistManager->AddHist(this->GetHist("distCloserJetCharm", "jetquark"), vars.GetEntriesFast());
  fHistManager->AddHist(this->GetHist("distFarerJetCharm", "jetquark"), vars.GetEntriesFast());

  fHistManager->AddHistsOption("BOX");
  fHistManager->Draw();

  vars.Delete();
  endings.Delete();
}



//____________________________________
 TH1* GFDstarHistsAnalysisMc::DrawDsJetQuark(const char *var, const char *dsOrJet, Int_t mode)
{

  const char *add = NULL;
  TString titles;
  switch(mode){
  case 0:
    add = "Vs";
    titles = Form(";%s(%s);%s(c/#bar{c})", var, dsOrJet, var);
    break;
  case 1:
    add = "Diff";
    titles = Form(";%s(%s)-%s(c/#bar{c})", var, dsOrJet, var);
    break;
  case 2:
    add = "Over";
    titles = Form(";%s(%s)/%s(c/#bar{c})", var, dsOrJet, var);
    break;
  default:
    this->Error("DrawDsJetQuark", "mode must be 0, 1 or 2");
    return NULL;
  }
  TH1 *h = this->GetHist(Form("%s%s%sQu%s", var, dsOrJet, var, add), "dsjetquark");
  if(!h) return NULL;
  if(h->GetDimension() == 1) h->SetMinimum(0.);
  
  titles = h->GetTitle() + titles;
  h->SetTitle(titles);

  fHistManager->Clear();
  fHistManager->AddHist(h, 0);

  fHistManager->Draw();

  return h;
}

//////////////////////////////////////////////////////
// not maintained anymore
//////////////////////////////////////////////////////
//
// //____________________________________
// TH1* GFDstarHistsAnalysisMc::CreateJetPurity(const char* variable)
// {
//   // n, phi, eta, pt
//   TString var(variable);
//   TH1* hAllRec = this->GetHist(var + "JetRec", "jetPurStab");
//   TH1* hGenRec = this->GetHist(var + "JetGenRec", "jetPurStab");

//   TH1* hPurity = this->CreatePurity(hAllRec, hGenRec, var+="^{jet}");

//   return hPurity;
// }

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateDsJetPurity(const char* variable)
{
  //Dphi Deta xGam PtJet PtDs EtaJet

  TString var(variable);
  TH1* hAllRec = this->GetHist(var + "Rec", "dsJetPurStab");
  TH1* hGenRec = this->GetHist(var + "RecGen", "dsJetPurStab");
  
  TH1* hPur = this->CreatePurity(hAllRec, hGenRec, GFAxisVariables::GetAxisLabel(variable), 
				 kFALSE);

  return hPur;
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateDsJetPurityHadCorr(const char* variable)
{
  //Dphi Deta xGam PtJet EtaJet

  TString var(variable);
  TH1* hAllRec = this->GetHist(var + "Gen", "partonJet");
  TH1* hGenRec = this->GetHist(var + "GenPs", "partonJet");
  
  TH1* hPur = this->CreatePurity(hAllRec, hGenRec, GFAxisVariables::GetAxisLabel(variable), 
				 kTRUE);

  return hPur;
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateDsDiJetPurity(const char* variable)
{
  //Dphi/Deta/xGam/Pt/Pt2 + Dijet Pt/Eta+Ds/DJetOthJet
  TString var(variable);
  TH1* hAllRec = this->GetHist(var + "Rec", "dsDiJetPurStab");
  TH1* hGenRec = this->GetHist(var + "RecGen", "dsDiJetPurStab");
  
  TH1* hPur = this->CreatePurity(hAllRec, hGenRec, GFAxisVariables::GetAxisLabel(variable), 
				 kFALSE);

  return hPur;
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreatePurity(TH1* hAllRec,TH1* hGenRec, const char* var,
					  Bool_t flagHadCorr)
{
  // purity as hGenRec/hAllRec, some 
  if(!hAllRec || !hGenRec) {
    return NULL;
  }
  if(!hAllRec->GetSumw2N()) hAllRec->Sumw2();
  if(!hGenRec->GetSumw2N()) hGenRec->Sumw2();
  
  TString newName(var);
  this->MakeStringName(newName += "Purity");
  TGraphAsymmErrors *grPtr = NULL;
  TH1* hPurity = GFHistManip::CreateRatioBinomErrors(hGenRec, hAllRec, NULL, NULL, NULL,&grPtr);
  if(!hPurity){
    this->Error("CreatePurity", "could not do binomial errors");
    hPurity = static_cast<TH1*>(hGenRec->Clone(newName));
    hPurity->Divide(hAllRec);
  } else {
    hPurity->SetName(newName.Data());
  }

  if(flagHadCorr){
    hPurity->SetTitle("purity: N^{gen}_{part}/N_{gen}");
  } else {
    hPurity->SetTitle("purity: N^{rec}_{gen}/N_{rec}");
  }

  fHistManager->Clear();
  fHistManager->AddHist(hAllRec);
  fHistManager->AddHist(hGenRec);
  fHistManager->AddHist(hPurity, 1);
  if(grPtr) fHistManager->AddObject(grPtr, 1, 0,"P");
  this->TmpGraph(grPtr);
  fHistManager->SetHistsXTitle(var);

  fHistManager->Draw();
  GFHistManip::PrintTable(hPurity, Form("%s: purity in bins of %s", fPeriod.Data(), var));

  return hPurity;
}

//////////////////////////////////////////////////////
// not maintained anymore
//////////////////////////////////////////////////////
//
// //____________________________________
// TH1* GFDstarHistsAnalysisMc::CreateJetStability(const char* variable)
// {
//   // n, phi, eta, pt
//
//   //   const Int_t oldTrig = this->GetTrigger();
//   //  this->SetTrigger(-1);
//
//   TString var(variable);
//   TH1* hGenRec = this->GetHist(var + "JetGenRec", "jetPurStab");
//   TH1* hGenNotRec = this->GetHist(var + "JetGenNonRec", "jetPurStab");
//
//   TH1* hStability = this->CreateStability(hGenRec, hGenNotRec, var+="^{jet}");
//
// //   this->SetTrigger(oldTrig);
//   return hStability;
// }

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateDsJetStability(const char* variable)
{
  //Dphi Deta xGam PtJet PtDs EtaJet

  TString var(variable);
  TH1* hGenRec = this->GetHist(var + "RecGen", "dsJetPurStab");
  TH1* hGenNotRec = this->GetHist(var + "NRecGen", "dsJetPurStab");
  
  TH1* hStability = this->CreateStability(hGenRec, hGenNotRec, 
					  GFAxisVariables::GetAxisLabel(variable), kFALSE);

  return hStability;
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateDsJetStabilityHadCorr(const char* variable)
{
  //xGam Deta Dphi EtaJet PtJet PtJet

  TString var(variable);
  TH1* hGenRec = this->GetHist(var + "GenPs", "partonJet");
  TH1* hGenNotRec = this->GetHist(var + "NGenPs", "partonJet");
  
  TH1* hStability = this->CreateStability(hGenRec, hGenNotRec, 
					  GFAxisVariables::GetAxisLabel(variable), kTRUE);

  return hStability;
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateDsDiJetStability(const char* variable)
{
  //Dphi/Deta/xGam/Pt/Pt2 + Dijet Pt/Eta+Ds/DJetOthJet

  TString var(variable);
  TH1* hGenRec = this->GetHist(var + "RecGen", "dsDiJetPurStab");
  TH1* hGenNotRec = this->GetHist(var + "NRecGen", "dsDiJetPurStab");
  
  TH1* hStability = this->CreateStability(hGenRec, hGenNotRec, 
					  GFAxisVariables::GetAxisLabel(variable), kFALSE);

  return hStability;
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::CreateStability(TH1* hGenRec, TH1* hGenNotRec, 
					     const char* variable, Bool_t flagHadCorr)
{
  if(!hGenNotRec || !hGenRec) {
    return NULL;
  }
  if(!hGenRec->GetSumw2N()) hGenRec->Sumw2();
  if(!hGenNotRec->GetSumw2N()) hGenNotRec->Sumw2();
  
  TString newName(variable);
  this->MakeStringName(newName += "Stability");
  TH1* hGen = GFHistManip::CreateAnalog(hGenRec,newName+"Tmp",kTRUE);
  hGen->Add(hGenNotRec, hGenRec);
  TGraphAsymmErrors *grPtr = NULL;
  TH1* hStability = GFHistManip::CreateRatioBinomErrors(hGenRec, hGen, NULL, NULL, NULL,&grPtr);
  if(!hStability){
    this->Error("CreateStability", "could not do binomial errors");
    hStability = static_cast<TH1*>(hGenRec->Clone(newName));
      hStability->Divide(hGenRec, hGen);//, 1., 1., "B");
  } else {
    hStability->SetName(newName);
  }
  delete hGen;
  if(flagHadCorr){
    hStability->SetTitle("stability: N^{gen}_{part}/(N^{gen}_{part}+N^{not gen}_{part})");
  } else {
    hStability->SetTitle("stability: N^{rec}_{gen}/(N^{rec}_{gen}+N^{not rec}_{gen})");
  }

  fHistManager->Clear();
  fHistManager->AddHist(hGenRec);
  fHistManager->AddHist(hGenNotRec);
  fHistManager->AddHist(hStability, 1);
  if(grPtr) fHistManager->AddObject(grPtr, 1, 0,"P");
  this->TmpGraph(grPtr);

  fHistManager->SetHistsXTitle(variable);

  fHistManager->Draw();
  GFHistManip::PrintTable(hStability, Form("%s: stability in bins of %s", fPeriod.Data(),
					   variable));

  return hStability;
}

// //____________________________________
// TH1* GFDstarHistsAnalysisMc::DrawJetGenVsRec(const char* variable, Int_t etaBin)
// {
//   // n, phi, eta, pt; for pt may choose etaBin
//   TString var(variable);
//   TString name(variable);
//   if(etaBin) name += Form("Eta%d", etaBin);
//   TH1* hGenVsRec = this->GetHist(name.Data(), "jetPurStab");
//   hGenVsRec->SetXTitle(var+"_{gen}");
//   hGenVsRec->SetYTitle(var+"_{rec}");

//   fHistManager->Clear();
//   fHistManager->AddHist(hGenVsRec);
//   fHistManager->Draw();

//   return hGenVsRec;
// }

//____________________________________
// TH1* GFDstarHistsAnalysisMc::DrawJetRes(const char* variable, const char* rel, 
// 					Int_t etaBin)
// {
//   // n, phi, eta, pt, possible combined with "Rel" (default: rel = "")
//   // for pt may choose etaBin 
//   TString var(variable);
//   var += "JetRes";
//   var += rel;
//   if(etaBin) var += Form("Eta%d", etaBin);

//   TH1* h = this->GetHist(var, "jetPurStab");

//   fHistManager->Clear();
//   fHistManager->AddHist(h);
//   fHistManager->Draw();

//   return h;
// }

//____________________________________
 TH1* GFDstarHistsAnalysisMc::DrawDsJetRes(const char* variable, const char* rel)
{
  // variable without 'Res', 'ResRel', if Rel needed: rel = "Rel"
  // EtaJetRes EtaJetPt<n>Res PtJetEta<n>ResRel PtJetResRel 
  // PtDsJetResRel PtDsResRel 
  TString var(variable);
  var += "Res";
  if(rel) var += rel;

  TH1* h = this->GetHist(var, "dsJetPurStab");
  fHistManager->Clear();
  fHistManager->AddHist(h);
  fHistManager->Draw();

  return h;
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::DrawDsJetGenVsRec(const char* var, Bool_t hadGen)
{
  // return last drawn hist

  TObjArray vars;
  if(var){
    vars.Add(new TObjString(var));
  } else {
    vars.Add(new TObjString("Dphi"));
    vars.Add(new TObjString("Deta"));
    vars.Add(new TObjString("xGam"));
    vars.Add(new TObjString("PtJet"));
    vars.Add(new TObjString("PtDs"));
    vars.Add(new TObjString("EtaJet"));
    vars.Add(new TObjString("PhiJet"));
    vars.Add(new TObjString("DaughJet"));
    vars.Add(new TObjString("ChFrac"));
    vars.Add(new TObjString("PtDsJet"));
    vars.Add(new TObjString("MDsJet"));
    vars.Add(new TObjString("CosTh"));
  }

  fHistManager->Clear();
  TString directory(hadGen ? "partonJet" : "dsJetPurStab");
  TString xLabel(hadGen ? "Ps" : "Gen");
  TString yLabel(hadGen ? "Gen" : "Rec");

  TH1* h = NULL;
  for(Int_t i = 0; i < vars.GetEntries(); ++i){
    h = this->GetHist(Form("%s%sVs%s", vars[i]->GetName(), yLabel.Data(), xLabel.Data()), 
		      directory.Data());
    if(h) {
      TString axisTitle = GFAxisVariables::GetAxisLabel(vars[i]->GetName());
      TString y = yLabel; y.ToLower();
      TString x = xLabel; x.ToLower();
      h->SetXTitle(axisTitle+Form("^{%s}", x.Data()));
      h->SetYTitle(axisTitle+Form("^{%s}", y.Data()));
      fHistManager->AddHist(h);
    }
  }

  fHistManager->AddHistsOption("COLZ");

  fHistManager->Draw();

  vars.Delete();
  return h;
}


//____________________________________
 TH1* GFDstarHistsAnalysisMc::DrawDsDiJetGenVsRec(const char* var)
{
  // return last drawn hist

  TObjArray vars;
  if(var){
    vars.Add(new TObjString(var));
  } else {
    vars.Add(new TObjString("PtDJet"));
    vars.Add(new TObjString("PtOthJet"));
    vars.Add(new TObjString("EtaDJet"));
    vars.Add(new TObjString("EtaOthJet"));
    vars.Add(new TObjString("DetaDijet"));
    vars.Add(new TObjString("DphiDijet"));
    vars.Add(new TObjString("PtDijet"));
    vars.Add(new TObjString("Pt2Dijet"));
    vars.Add(new TObjString("CosThDijet"));
    vars.Add(new TObjString("xGamDijet"));
    //   vars.Add(new TObjString("DeltaRDijet"));
  }

  fHistManager->Clear();
  TString directory("dsDiJetPurStab");
  TString xLabel("Gen");
  TString yLabel("Rec");

  TH1* h = NULL;
  for(Int_t i = 0; i < vars.GetEntries(); ++i){
    h = this->GetHist(Form("%s%sVs%s", vars[i]->GetName(), yLabel.Data(), xLabel.Data()), 
		      directory.Data());
    if(h) {
      TString axisTitle = GFAxisVariables::GetAxisLabel(vars[i]->GetName());
      TString y = yLabel; y.ToLower();
      TString x = xLabel; x.ToLower();
      h->SetXTitle(axisTitle+Form("^{%s}", x.Data()));
      h->SetYTitle(axisTitle+Form("^{%s}", y.Data()));
      fHistManager->AddHist(h);
    }
  }

  fHistManager->AddHistsOption("COLZ");

  fHistManager->Draw();

  vars.Delete();
  return h;
}

//____________________________________
 TH2* GFDstarHistsAnalysisMc::DrawEffCorrWeights(const char* var)
{
  // upper case Pt, Eta...
  TH2 *h = static_cast<TH2*>(this->GetHist(Form("weightHist%sEffCor", var), "weights"));
  const Int_t oldStat = gStyle->GetOptStat();
  if(h){
    gStyle->SetOptStat(110100);
    fHistManager->Clear();
    fHistManager->AddHist(h);
    h->SetOption("BOX");
    const char *label = GFAxisVariables::GetAxisLabel(var);
    h->SetTitle(Form("%s;%s;efficiency correction", h->GetTitle(), label));
    GFHistArray *hists1D = GFHistManip::CreateHistArray1D
      (h, label, /*name*/NULL, /*noOverUnderFlow*/kTRUE, /*projX*/kFALSE);//,mergeNBins = 1);
    for(Int_t i = 0; i < hists1D->GetEntriesFast(); ++i){
      hists1D->At(i)->SetTitle(Form("%s;efficiency correction;arbitrary units", 
				    hists1D->At(i)->GetTitle()));
      hists1D->At(i)->SetMinimum(0.);
      hists1D->At(i)->SetOption("HIST");
      GFHistManip::Scale(hists1D->At(i), 100000.);
    }
    fHistManager->AddHists(hists1D, 1);
    delete hists1D;

    fHistManager->Draw();
    fHistManager->Update();// to make the style option work...
    if(gPad->GetMother()) gPad->GetMother()->cd(); // to get it work even in last pad...
  }
  gStyle->SetOptStat(oldStat);
  return h;
}

//____________________________________
 TArrayD GFDstarHistsAnalysisMc::TotalRecEff()
{
  //  returning effRec and its error
  // cf. log book 13.9.2002
  const Double_t numGen = this->GetHistGenWeight("eta")->Integral();
  TString name(Form("tmpRecEff%s", this->GetPeriod()));
  GFHistManip::MakeUniqueName(name);
  TH1* hRec = this->GetHistPtCut();//2D("eta")->ProjectionX(name,-1,-1,"E");
  
  TArrayD result = this->TotalRecEff(hRec, numGen, "D*:");
  delete hRec;
  return result;
}

//____________________________________
 TArrayD GFDstarHistsAnalysisMc::TotalRecEffDsJet(const char* region, 
						 UInt_t dmRefFlag)
{

  TString name("dstar1Jet");
  name += region;
  
  const Double_t numGen = this->GetHist(name + "GenWeight", "genJet")->Integral();
//   TString nameTmp(Form("tmpRecEffDsJet%s%s", region, this->GetPeriod()));
//   GFHistManip::MakeUniqueName(nameTmp);
  TH1* hRec = this->GetHistDs1Jet(region);
  TH1* hRecRef = this->GetHistDs1JetRef(region, dmRefFlag);

  TArrayD result = this->TotalRecEff(hRec, numGen, "D* + 1 jet:", hRecRef);
  delete hRec;
  return result;
}
//____________________________________
 TArrayD GFDstarHistsAnalysisMc::TotalRecEffDsDiJet(UInt_t dmRefFlag)
{
// if dmRefFlag == 2 (default): reference fit from inclusive D*

  TString name("dstarDiJet");
  
  const Double_t numGen = this->GetHist(name + "GenWeight", "genDiJet")->Integral();
//   TString nameTmp(Form("tmpRecEffDsDiJet%s", this->GetPeriod()));
//   GFHistManip::MakeUniqueName(nameTmp);
  TH1 *hRec = this->GetHist(name, "dsDiJets");

  TH1* hRecRef = NULL;
  if (dmRefFlag == 2) hRecRef = this->GetHistPtCut();

  TArrayD result = this->TotalRecEff(hRec, numGen, "D* dijet:", hRecRef);
  delete hRec;
  return result;
}

//____________________________________
 TArrayD GFDstarHistsAnalysisMc::TotalRecEff(TH1 *hRec, Double_t nGen, 
					    const char* prefix, TH1* hRecRef)
{
  //  returning effRec and its error
  // cf. log book 13.9.2002
  
  TArrayD effErr(2);
//   effErr[0] = nGen ? nRec/nGen : 0.;
//   effErr[1] = nGen ? TMath::Sqrt(effErr[0] * (1.- effErr[0])/nGen) : 0.;

  TArrayD nDstar;
  if(hRecRef){
    fMultiDmFitter->FitAll(hRecRef);
    nDstar = fMultiDmFitter->FitSingle(hRec, NULL);
  } else {
    nDstar = fMultiDmFitter->FitAll(hRec);
  }
  effErr[0] = (nGen ? nDstar[0]/nGen : 0.);
  effErr[1] = (nGen ? nDstar[1]/nGen : 0.);
 
  this->PrintPlusMinus(effErr, 
		       Form("%s rec eff for S%d from %s", prefix, fTrigger, fPeriod.Data()));
  return effErr;
}

//____________________________________
 TArrayD GFDstarHistsAnalysisMc::TotalTrigEff(Int_t TE)
{
  // D*: 19, 31, -1 (= 19 && 31)
  TH1* hTrig = this->GetHistTrig("eta", TE);
  TH1* hNonTrig = this->GetHistNonTrig("eta");

  const Double_t numTrig = hTrig->Integral();
  const Double_t numRec = hNonTrig->Integral();
  Stat_t s[11];
  hNonTrig->GetStats(s); // now s[1] is sum of squares of weights, s[0] is sum of weights
  const Double_t nEffRec = s[0]*s[0]/s[1];
  const Double_t nFillRec = hNonTrig->GetEntries();

  delete hTrig;
  delete hNonTrig;

  return this->TotalTrigEff(numTrig, numRec, nEffRec, nFillRec, TE, "D*:");
}

//____________________________________
 TArrayD GFDstarHistsAnalysisMc::TotalTrigEffDsJet(const char* region, Int_t TE)
{
  // D* + jet: returning effTrig and its error
  // Back, Forw, ''; 19, 31, default: both

  TString histName("dstar1Jet");
  histName += region;
  TString trigName = histName + "Trig";
  if(TE >= 0) trigName += TE;

  TH1* hTrig = this->GetHist(trigName, "trigJet");
  TH1* hNonTrig = this->GetHist(histName + "NonTrig", "trigJet");

  const Double_t numTrig = hTrig->Integral();
  const Double_t numRec = hNonTrig->Integral();
  Stat_t s[11];
  hNonTrig->GetStats(s); // now s[1] is sum of squares of weights, s[0] is sum of weights
  const Double_t nEffRec = s[0]*s[0]/s[1];
  const Double_t nFillRec = hNonTrig->GetEntries();

  delete hTrig;
  delete hNonTrig;

  return this->TotalTrigEff(numTrig, numRec, nEffRec,nFillRec,TE,Form("D* + 1 jet %s:",region));
}

//____________________________________
 TArrayD GFDstarHistsAnalysisMc::TotalTrigEffDsDiJet(Int_t TE)
{
  // D* dijet: returning effTrig and its error
  // Back, Forw, ''; 19, 31, default: both

  TString histName("dstarDiJet");
  TString trigName = histName + "Trig";
  if(TE >= 0) trigName += TE;

  TH1* hTrig = this->GetHist(trigName, "trigDiJet");
  TH1* hNonTrig = this->GetHist(histName + "NonTrig", "trigDiJet");

  const Double_t numTrig = hTrig->Integral();
  const Double_t numRec = hNonTrig->Integral();
  Stat_t s[11];
  hNonTrig->GetStats(s); // now s[1] is sum of squares of weights, s[0] is sum of weights
  const Double_t nEffRec = s[0]*s[0]/s[1];
  const Double_t nFillRec = hNonTrig->GetEntries();

  delete hTrig;
  delete hNonTrig;

  return this->TotalTrigEff(numTrig, numRec, nEffRec, nFillRec, TE, "D* dijet:");
}

//____________________________________
 TArrayD GFDstarHistsAnalysisMc::TotalTrigEff(Double_t nTrig, Double_t nRec, 
					     Double_t nEffRec, Double_t nFillRec, Int_t TE,
					     const char* prefix)
{
  // returning effTrig and its error: 
  // nTrig: trig & rec, nRec: all rec, nEffRec: effective number entries for nRec
  // nFillRec: hist fill actions for nRec
  // if(nEffRec/nRec >= 0.7): do asymm. binomial error GFMath::BinomialError
  // then result[2/3] are the asymmetric errors
  // result[1]: symmetric easy binomial error
  
  TArrayD effErr = GFMath::Efficiency(nTrig, nRec, nEffRec, nFillRec);

  TString text(Form("%s trig eff for S%d from %s, TE ", prefix, fTrigger, fPeriod.Data()));
  if(TE >= 0) text += TE; else text += "19+31";
  this->PrintPlusMinus(effErr, text);

  return effErr;
}

//____________________________________
 TArrayD GFDstarHistsAnalysisMc::TotalAccept()
{
  //  returning acceptance (and - not yet implemented - its error)

  const Double_t sumAccept = this->GetHistGenWeight("eta")->Integral(); // cf. comment in 
//   Double_t sumAccept = this->GetHistGenAccept("eta")->Integral();     // CreateAccept
  const Double_t numGen = this->GetHistGenNoAccWeight("eta")->Integral();

  return this->TotalAccept(sumAccept, numGen, "D*:");
}

//____________________________________
 TArrayD GFDstarHistsAnalysisMc::TotalAccept(Double_t sumAccept, Double_t numGen, 
					    const char* prefix)
{
  // sumAccept: sum of acceptances of generated events, numGen: num of generated events
  // error not implemented!
  TArrayD acceptErr(2);
  acceptErr[0] = numGen ? sumAccept/numGen : 0.;
  acceptErr[1] = 0;
  this->Warning("TotalAccept", "Error not yet implemented (set 0.)!");

  this->PrintPlusMinus(acceptErr, Form("%s tag accept for S%d from %s", prefix, fTrigger, 
				       fPeriod.Data()));
  return acceptErr;
}

//____________________________________
 TArrayD GFDstarHistsAnalysisMc::TotalAcceptDsJet(const char* region)
{
  //  D* + 1 jet: returning acceptance (and - not yet implemented - its error)
  // region: Forw, Back, ''

  TString histName("dstar1Jet");
  histName += region;

  // which hists: cf. CreateAcceptDsJet(...)
  const Double_t sumAccept = this->GetHist(histName + "GenWeight", "genJet")->Integral();
  const Double_t numGen = this->GetHist(histName + "GenNoAccWeight", "genJet")->Integral();

  return this->TotalAccept(sumAccept, numGen, Form("D* + 1 jet %s:", region));
}


//____________________________________
 TArrayD GFDstarHistsAnalysisMc::TotalAcceptDsDiJet()
{
  //  D* dijet: returning acceptance (and - not yet implemented - its error)

  TString histName("dstarDiJet");

  // which hists: cf. CreateAcceptDsDiJet(...)
  const Double_t sumAccept = this->GetHist(histName + "GenWeight", "genDiJet")->Integral();
  const Double_t numGen = this->GetHist(histName + "GenNoAccWeight", "genDiJet")->Integral();

  return this->TotalAccept(sumAccept, numGen, "D* dijet:");
}


//____________________________________
//____________________________________
//____________________________________
 TH1* GFDstarHistsAnalysisMc::GetHistGen(const char* variable)
{
  TString var(variable);
  return this->GetHist(var+="HistGen", "gen");
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::GetHistGenWeight(const char* variable)
{
  TString var(variable);
  return this->GetHist(var+="HistGenWeight", "gen");
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::GetHistGenNoAccWeight(const char* variable)
{
  TString var(variable);
  return this->GetHist(var+="HistGenNoAccWeight", "gen");
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::GetHistGenAccept(const char* variable)
{
  TString var(variable);
  return this->GetHist(var+="HistGenAcc", "gen");
}


//____________________________________
 TH1* GFDstarHistsAnalysisMc::GetHistTrig(const char* variable, Int_t TE)
{
  TString var(variable);
  var+="HistTrig";
  if(TE >= 0) var += TE;
  return this->GetHist(var, "trig");
}


//____________________________________
 TH1* GFDstarHistsAnalysisMc::GetHistNonTrig(const char* variable)
{
  TString var(variable);
  var+="HistNonTrig";
  return this->GetHist(var, "trig");
}

//____________________________________
 TH1* GFDstarHistsAnalysisMc::GetHistAccept(const char* variable)
{
  TString var(variable);

  return this->GetHist(var+="HistAcc");
}
ROOT page - Class index - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.