////////////////////////////////////////////////////////////
//
// Filename: ComparePlots.cpp
// Created:  14-NOV-2007,  Shaohua Fu
//
// Purpose:  Compare histograms from two root files:
//           the first is black, and the second is red
//           If there is only one file, then just plot
//
////////////////////////////////////////////////////////////

#include "ComparePlots.hpp"
#include <iostream>

#include "TStyle.h"
#include "TCanvas.h"
#include "TPostScript.h"
#include "TText.h"

using namespace std;


int main(int argc, char *argv[])
{
  if (argc < 2 || argc > 4) {
    cout<<"Usage: ComparePlots mctest.root (optional)mctest2.root (optional)test.ps"<<endl;
    cout<<"    or ComparePlots mctest.root (optional)test.ps"<<endl;
    return 2;
  }

  TString filename1 = argv[1];
  TString filename2 = "";
  TString psfilename = "test.ps";

  if (argc > 2) filename2 = argv[2];
  if (argc > 3) psfilename = argv[3];
  if (filename2(filename2.Sizeof()-4,3)==".ps")
    {psfilename = filename2; filename2 = "";}
  cout<<"Plots are made in file "<<psfilename<<endl;


  //// set background color to be white ////
  gStyle->SetCanvasColor(10);
  gStyle->SetStatColor(10);
  gStyle->SetTitleColor(10);
  gStyle->SetOptStat(111111);
  gStyle->SetOptFit(2);


  bool _comp = true;
  if (filename2=="") _comp = false;

  HistoFile hist1(filename1);
  HistoFile hist2;
  if (_comp) {
    hist2.Init(filename2); hist2.SetColor(2);
    if (hist2._h1pvx->Integral()!=0) {
      double scale = hist1._h1pvx->Integral()/hist2._h1pvx->Integral();
      hist2.Scale(scale); cout<<"Scale the second set of histograms by "<<scale<<endl;
    }
  }


  TCanvas *cc = new TCanvas("cc","mctest",900,1200);
  TPad *c0 = new TPad("c0","",0.01,0.95,0.99,0.99,10,2);
  c0->Draw();
  TPad *c1 = new TPad("c1","",0.01,0.01,0.99,0.94,10,2);
  c1->Draw();

  // postscript output type
  // type = 111   portrait  ps
  // type = 112   landscape ps
  // type = 113   eps
  TPostScript *ps = new TPostScript(psfilename, 111);
  TText *title = new TText();
  title->SetTextSize(0.6);
  title->SetTextColor(1);
  Float_t x=0.4;
  Float_t y=0.4;

  ps->NewPage();
  c0->cd(); c0->Clear();
  title->DrawTextNDC(x,y,"MC Particles");
  c0->Update();
  c1->Clear(); c1->Divide(2,3);
  c1->cd(1); hist1._h1nmc->Draw();   if (_comp) hist2._h1nmc->Draw("same");
  c1->cd(2); hist1._h1nstab->Draw(); if (_comp) hist2._h1nstab->Draw("same");
  c1->cd(3); hist1._h1mcpt->Draw();  if (_comp) hist2._h1mcpt->Draw("same");
  c1->cd(4); hist1._h1mceta->Draw(); if (_comp) hist2._h1mceta->Draw("same");
  c1->cd(5); hist1._h1mcphi->Draw(); if (_comp) hist2._h1mcphi->Draw("same");
  c1->Update();

  ps->NewPage();
  c0->cd(); c0->Clear();
  title->DrawTextNDC(x,y,"Electrons");
  c0->Update();
  c1->Clear(); c1->Divide(3,3);
  c1->cd(1); hist1._h1nel_all->Draw(); if (_comp) hist2._h1nel_all->Draw("same");
  c1->cd(2); hist1._h1elpt_all->Draw(); if (_comp) hist2._h1elpt_all->Draw("same");
  c1->cd(4); hist1._h1nel->Draw();    if (_comp) hist2._h1nel->Draw("same");
  c1->cd(5); hist1._h1elpt->Draw();   if (_comp) hist2._h1elpt->Draw("same");
  c1->cd(6); hist1._h1eleta->Draw();  if (_comp) hist2._h1eleta->Draw("same");
  c1->cd(7); hist1._h1elphi->Draw();  if (_comp) hist2._h1elphi->Draw("same");
  c1->cd(8); hist1._h1elpid->Draw();  if (_comp) hist2._h1elpid->Draw("same");
  c1->cd(9); hist1._h1elpid2->Draw(); if (_comp) hist2._h1elpid2->Draw("same");
  c1->Update();

  ps->NewPage();
  c0->cd(); c0->Clear();
  title->DrawTextNDC(x,y,"Muons");
  c0->Update();
  c1->Clear(); c1->Divide(3,3);
  c1->cd(1); hist1._h1nmu_all->Draw(); if (_comp) hist2._h1nmu_all->Draw("same");
  c1->cd(2); hist1._h1mupt_all->Draw(); if (_comp) hist2._h1mupt_all->Draw("same");
  c1->cd(4); hist1._h1nmu->Draw();    if (_comp) hist2._h1nmu->Draw("same");
  c1->cd(5); hist1._h1mupt->Draw();   if (_comp) hist2._h1mupt->Draw("same");
  c1->cd(6); hist1._h1mueta->Draw();  if (_comp) hist2._h1mueta->Draw("same");
  c1->cd(7); hist1._h1muphi->Draw();  if (_comp) hist2._h1muphi->Draw("same");
  c1->cd(8); hist1._h1mupid->Draw();  if (_comp) hist2._h1mupid->Draw("same");
  c1->cd(9); hist1._h1mupid2->Draw(); if (_comp) hist2._h1mupid2->Draw("same");
  c1->Update();

  ps->NewPage();
  c0->cd(); c0->Clear();
  title->DrawTextNDC(x,y,"Parton Jets");
  c0->Update();
  c1->Clear(); c1->Divide(2,3);
  c1->cd(1); hist1._h1_jet_n->Draw();   if (_comp) hist2._h1_jet_n->Draw("same");
  c1->cd(3); hist1._h1_jet_pt->Draw();  if (_comp) hist2._h1_jet_pt->Draw("same");
  c1->cd(5)->SetLogy(1); hist1._h1_jet_pt2->Draw(); if (_comp) hist2._h1_jet_pt2->Draw("same");
  c1->cd(2); hist1._h1_jet_n_no_e->Draw(); if (_comp) hist2._h1_jet_n_no_e->Draw("same");
  c1->cd(4); hist1._h1_jet_pt_no_e->Draw(); if (_comp) hist2._h1_jet_pt_no_e->Draw("same");
  c1->cd(6)->SetLogy(1); hist1._h1_jet_pt2_no_e->Draw(); if (_comp) hist2._h1_jet_pt2_no_e->Draw("same");
  c1->Update();

  ps->NewPage();
  c0->cd(); c0->Clear();
  title->DrawTextNDC(x,y,"Parton Jets (continued)");
  c0->Update();
  c1->Clear(); c1->Divide(2,3);
  c1->cd(1); hist1._h1_jet_eta->Draw(); if (_comp) hist2._h1_jet_eta->Draw("same");
  c1->cd(2); hist1._h1_jet_phi->Draw(); if (_comp) hist2._h1_jet_phi->Draw("same");
  c1->cd(3); hist1._h1_jet_mjj->Draw(); if (_comp) hist2._h1_jet_mjj->Draw("same");
  c1->cd(4); hist1._h1_jet_dphi->Draw(); if (_comp) hist2._h1_jet_dphi->Draw("same");
  c1->Update();

  ps->NewPage();
  c0->cd(); c0->Clear();
  title->DrawTextNDC(x,y,"Vertex");
  c0->Update();
  c1->Clear(); c1->Divide(2,3);
  c1->cd(1); hist1._h1pvx->Draw(); if (_comp) hist2._h1pvx->Draw("same");
  c1->cd(2); hist1._h1pvy->Draw(); if (_comp) hist2._h1pvy->Draw("same");
  c1->cd(3); hist1._h1pvz->Draw(); if (_comp) hist2._h1pvz->Draw("same");
  //c1->cd(4); hist1._h1nvtx->Draw(); if (_comp) hist2._h1nvtx->Draw("same");
  c1->Update();

  ps->NewPage();
  c0->cd(); c0->Clear();
  title->DrawTextNDC(x,y,"Event Information");
  c0->Update();
  c1->Clear(); c1->Divide(2,3);
  c1->cd(1); hist1._h1pdfx1->Draw();  if (_comp) hist2._h1pdfx1->Draw("same");
  c1->cd(2); hist1._h1pdfx2->Draw();  if (_comp) hist2._h1pdfx2->Draw("same");
  c1->cd(3); hist1._h1xsec->Draw();   if (_comp) hist2._h1xsec->Draw("same");
  if (hist1._h1_xsec_dist) {
    c1->cd(4); hist1._h1_xsec_dist->Draw();
    if (_comp && hist2._h1_xsec_dist) {c1->cd(5); hist2._h1_xsec_dist->Draw();}
  }
  c1->cd(6); hist1._h1lumi->Draw();   if (_comp) hist2._h1lumi->Draw("same");
  c1->Update();

  ps->NewPage();
  c0->cd(); c0->Clear();
  title->DrawTextNDC(x,y,"Generated W Bosons");
  c0->Update();
  c1->Clear(); c1->Divide(2,3);
  c1->cd(1); hist1._h1_w_n->Draw();   if (_comp) hist2._h1_w_n->Draw("same");
  c1->cd(2); hist1._h1_w_pt->Draw();  if (_comp) hist2._h1_w_pt->Draw("same");
  c1->cd(3)->SetLogy(1); if (hist1._h1_w_pt2->Integral()>0) hist1._h1_w_pt2->Draw();
  if (_comp && hist2._h1_w_pt2->Integral()>0) hist2._h1_w_pt2->Draw("same");
  c1->cd(4); hist1._h1_w_eta->Draw(); if (_comp) hist2._h1_w_eta->Draw("same");
  c1->cd(5); hist1._h1_w_phi->Draw(); if (_comp) hist2._h1_w_phi->Draw("same");
  c1->cd(6); hist1._h1_w_mass->Draw(); if (_comp) hist2._h1_w_mass->Draw("same");
  c1->Update();

  ps->NewPage();
  c0->cd(); c0->Clear();
  title->DrawTextNDC(x,y,"Generated Z Bosons");
  c0->Update();
  c1->Clear(); c1->Divide(2,3);
  c1->cd(1); hist1._h1_z_n->Draw();   if (_comp) hist2._h1_z_n->Draw("same");
  c1->cd(2); hist1._h1_z_pt->Draw();  if (_comp) hist2._h1_z_pt->Draw("same");
  c1->cd(3)->SetLogy(1); if (hist1._h1_z_pt2->Integral()>0) hist1._h1_z_pt2->Draw();
  if (_comp && hist2._h1_z_pt2->Integral()>0) hist2._h1_z_pt2->Draw("same");
  c1->cd(4); hist1._h1_z_eta->Draw(); if (_comp) hist2._h1_z_eta->Draw("same");
  c1->cd(5); hist1._h1_z_phi->Draw(); if (_comp) hist2._h1_z_phi->Draw("same");
  c1->cd(6); hist1._h1_z_mass->Draw(); if (_comp) hist2._h1_z_mass->Draw("same");
  c1->Update();

  ps->NewPage();
  c0->cd(); c0->Clear();
  title->DrawTextNDC(x,y,"Generated t quarks");
  c0->Update();
  c1->Clear(); c1->Divide(2,3);
  c1->cd(1); hist1._h1_t_n->Draw();    if (_comp) hist2._h1_t_n->Draw("same");
  c1->cd(2); hist1._h1_t_pt->Draw();   if (_comp) hist2._h1_t_pt->Draw("same");
  c1->cd(3); hist1._h1_t_eta->Draw();  if (_comp) hist2._h1_t_eta->Draw("same");
  c1->cd(4); hist1._h1_t_phi->Draw();  if (_comp) hist2._h1_t_phi->Draw("same");
  c1->cd(5); hist1._h1_t_mass->Draw(); if (_comp) hist2._h1_t_mass->Draw("same");
  c1->Update();

  ps->NewPage();
  c0->cd(); c0->Clear();
  title->DrawTextNDC(x,y,"Generated b quarks");
  c0->Update();
  c1->Clear(); c1->Divide(3,3);
  c1->cd(1); hist1._h1_b_n->Draw();   if (_comp) hist2._h1_b_n->Draw("same");
  c1->cd(2); hist1._h1_b_eta->Draw(); if (_comp) hist2._h1_b_eta->Draw("same");
  c1->cd(3); hist1._h1_b_phi->Draw(); if (_comp) hist2._h1_b_phi->Draw("same");
  c1->cd(4); hist1._h1_b_pt->Draw();  if (_comp) hist2._h1_b_pt->Draw("same");
  c1->cd(5); hist1._h1_b1_pt->Draw(); if (_comp) hist2._h1_b1_pt->Draw("same");
  c1->cd(6); hist1._h1_b2_pt->Draw(); if (_comp) hist2._h1_b2_pt->Draw("same");
  c1->cd(7); hist1._h1_b3_pt->Draw(); if (_comp) hist2._h1_b3_pt->Draw("same");
  c1->cd(8); hist1._h1_b4_pt->Draw(); if (_comp) hist2._h1_b4_pt->Draw("same");
  c1->cd(9); hist1._h1_b_drmin->Draw(); if (_comp) hist2._h1_b_drmin->Draw("same");
  c1->Update();

  ps->NewPage();
  c0->cd(); c0->Clear();
  title->DrawTextNDC(x,y,"Generated c quarks");
  c0->Update();
  c1->Clear(); c1->Divide(3,3);
  c1->cd(1); hist1._h1_c_n->Draw();   if (_comp) hist2._h1_c_n->Draw("same");
  c1->cd(2); hist1._h1_c_eta->Draw(); if (_comp) hist2._h1_c_eta->Draw("same");
  c1->cd(3); hist1._h1_c_phi->Draw(); if (_comp) hist2._h1_c_phi->Draw("same");
  c1->cd(4); hist1._h1_c_pt->Draw();  if (_comp) hist2._h1_c_pt->Draw("same");
  c1->cd(5); hist1._h1_c1_pt->Draw(); if (_comp) hist2._h1_c1_pt->Draw("same");
  c1->cd(6); hist1._h1_c2_pt->Draw(); if (_comp) hist2._h1_c2_pt->Draw("same");
  c1->cd(7); hist1._h1_c3_pt->Draw(); if (_comp) hist2._h1_c3_pt->Draw("same");
  c1->cd(8); hist1._h1_c4_pt->Draw(); if (_comp) hist2._h1_c4_pt->Draw("same");
  c1->cd(9); hist1._h1_c_drmin->Draw(); if (_comp) hist2._h1_c_drmin->Draw("same");
  c1->Update();

  ps->NewPage();
  c0->cd(); c0->Clear();
  title->DrawTextNDC(x-0.1,y,"Generated light-partons");
  c0->Update();
  c1->Clear(); c1->Divide(2,3);
  c1->cd(1); hist1._h1_lp_n_all->Draw(); if (_comp) hist2._h1_lp_n_all->Draw("same");
  c1->cd(2); hist1._h1_lp_pt_all->Draw(); if (_comp) hist2._h1_lp_pt_all->Draw("same");
  c1->cd(3); hist1._h1_lp_n->Draw();   if (_comp) hist2._h1_lp_n->Draw("same");
  c1->cd(4); hist1._h1_lp_pt->Draw();  if (_comp) hist2._h1_lp_pt->Draw("same");
  c1->cd(5); hist1._h1_lp_eta->Draw(); if (_comp) hist2._h1_lp_eta->Draw("same");
  c1->cd(6); hist1._h1_lp_phi->Draw(); if (_comp) hist2._h1_lp_phi->Draw("same");
  c1->Update();

  ps->NewPage();
  c0->cd(); c0->Clear();
  title->DrawTextNDC(x-0.1,y,"Generated light-partons (continued)");
  c0->Update();
  c1->Clear(); c1->Divide(2,3);
  c1->cd(1); hist1._h1_lp1_pt->Draw(); if (_comp) hist2._h1_lp1_pt->Draw("same");
  c1->cd(2); hist1._h1_lp2_pt->Draw(); if (_comp) hist2._h1_lp2_pt->Draw("same");
  c1->cd(3); hist1._h1_lp3_pt->Draw(); if (_comp) hist2._h1_lp3_pt->Draw("same");
  c1->cd(4); hist1._h1_lp4_pt->Draw(); if (_comp) hist2._h1_lp4_pt->Draw("same");
  c1->cd(5); hist1._h1_lp5_pt->Draw(); if (_comp) hist2._h1_lp5_pt->Draw("same");
  c1->cd(6); hist1._h1_lp_drmin->Draw(); if (_comp) hist2._h1_lp_drmin->Draw("same");
  c1->Update();


  ps->Close();
  return 0;
}


////////////////////////////////////////////////////////////
// Functions for class HistoFile

HistoFile::HistoFile(const TString filename) {
  Init(filename);
}

void HistoFile::Init(const TString filename) {
  _hfile = new TFile(filename, "READ");

  for (int i=0; i<200; ++i) _h1_array[i] = 0;

  int j = 0;

  _h1nmc      = (TH1F*)_hfile->Get("nmc");     _h1_array[j++] = _h1nmc;
  _h1nstab    = (TH1F*)_hfile->Get("nstab");   _h1_array[j++] = _h1nstab;
  _h1mcpt     = (TH1F*)_hfile->Get("mcpt");    _h1_array[j++] = _h1mcpt;
  _h1mceta    = (TH1F*)_hfile->Get("mceta");   _h1_array[j++] = _h1mceta;
  _h1mcphi    = (TH1F*)_hfile->Get("mcphi");   _h1_array[j++] = _h1mcphi;

  _h1nel_all  = (TH1F*)_hfile->Get("nel_all"); _h1_array[j++] = _h1nel_all;
  _h1elpt_all = (TH1F*)_hfile->Get("elpt_all");_h1_array[j++] = _h1elpt_all;
  _h1nel      = (TH1F*)_hfile->Get("nel");     _h1_array[j++] = _h1nel;
  _h1elpt     = (TH1F*)_hfile->Get("elpt");    _h1_array[j++] = _h1elpt;
  _h1eleta    = (TH1F*)_hfile->Get("eleta");   _h1_array[j++] = _h1eleta;
  _h1elphi    = (TH1F*)_hfile->Get("elphi");   _h1_array[j++] = _h1elphi;
  _h1elpid    = (TH1F*)_hfile->Get("elpid");   _h1_array[j++] = _h1elpid;
  _h1elpid2   = (TH1F*)_hfile->Get("elpid2");  _h1_array[j++] = _h1elpid2;

  _h1nmu_all  = (TH1F*)_hfile->Get("nmu_all"); _h1_array[j++] = _h1nmu_all;
  _h1mupt_all = (TH1F*)_hfile->Get("mupt_all");_h1_array[j++] = _h1mupt_all;
  _h1nmu      = (TH1F*)_hfile->Get("nmu");     _h1_array[j++] = _h1nmu;
  _h1mupt     = (TH1F*)_hfile->Get("mupt");    _h1_array[j++] = _h1mupt;
  _h1mueta    = (TH1F*)_hfile->Get("mueta");   _h1_array[j++] = _h1mueta;
  _h1muphi    = (TH1F*)_hfile->Get("muphi");   _h1_array[j++] = _h1muphi;
  _h1mupid    = (TH1F*)_hfile->Get("mupid");   _h1_array[j++] = _h1mupid;
  _h1mupid2   = (TH1F*)_hfile->Get("mupid2");  _h1_array[j++] = _h1mupid2;

  _h1xsec     = (TH1F*)_hfile->Get("xsec");    _h1_array[j++] = _h1xsec;
  _h1pdfx1    = (TH1F*)_hfile->Get("pdfx1");   _h1_array[j++] = _h1pdfx1;
  _h1pdfx2    = (TH1F*)_hfile->Get("pdfx2");   _h1_array[j++] = _h1pdfx2;
  _h1nvtx     = (TH1F*)_hfile->Get("nvtx");    _h1_array[j++] = _h1nvtx;
  _h1pvx      = (TH1F*)_hfile->Get("pvx");     _h1_array[j++] = _h1pvx;
  _h1pvy      = (TH1F*)_hfile->Get("pvy");     _h1_array[j++] = _h1pvy;
  _h1pvz      = (TH1F*)_hfile->Get("pvz");     _h1_array[j++] = _h1pvz;
  _h1lumi     = (TH1F*)_hfile->Get("lumi");    _h1_array[j++] = _h1lumi;

  _h1_w_n     = (TH1F*)_hfile->Get("w_n");     _h1_array[j++] = _h1_w_n;
  _h1_w_pt    = (TH1F*)_hfile->Get("w_pt");    _h1_array[j++] = _h1_w_pt;
  _h1_w_pt2   = (TH1F*)_hfile->Get("w_pt2");   _h1_array[j++] = _h1_w_pt2;
  _h1_w_eta   = (TH1F*)_hfile->Get("w_eta");   _h1_array[j++] = _h1_w_eta;
  _h1_w_phi   = (TH1F*)_hfile->Get("w_phi");   _h1_array[j++] = _h1_w_phi;
  _h1_w_mass  = (TH1F*)_hfile->Get("w_mass");  _h1_array[j++] = _h1_w_mass;
  _h1_z_n     = (TH1F*)_hfile->Get("z_n");     _h1_array[j++] = _h1_z_n;
  _h1_z_pt    = (TH1F*)_hfile->Get("z_pt");    _h1_array[j++] = _h1_z_pt;
  _h1_z_pt2   = (TH1F*)_hfile->Get("z_pt2");   _h1_array[j++] = _h1_z_pt2;
  _h1_z_eta   = (TH1F*)_hfile->Get("z_eta");   _h1_array[j++] = _h1_z_eta;
  _h1_z_phi   = (TH1F*)_hfile->Get("z_phi");   _h1_array[j++] = _h1_z_phi;
  _h1_z_mass  = (TH1F*)_hfile->Get("z_mass");  _h1_array[j++] = _h1_z_mass;

  _h1_t_n     = (TH1F*)_hfile->Get("t_n");     _h1_array[j++] = _h1_t_n;
  _h1_t_pt    = (TH1F*)_hfile->Get("t_pt");    _h1_array[j++] = _h1_t_pt;
  _h1_t_eta   = (TH1F*)_hfile->Get("t_eta");   _h1_array[j++] = _h1_t_eta;
  _h1_t_phi   = (TH1F*)_hfile->Get("t_phi");   _h1_array[j++] = _h1_t_phi;
  _h1_t_mass  = (TH1F*)_hfile->Get("t_mass");  _h1_array[j++] = _h1_t_mass;

  _h1_b_n     = (TH1F*)_hfile->Get("b_n");     _h1_array[j++] = _h1_b_n;
  _h1_b_pt    = (TH1F*)_hfile->Get("b_pt");    _h1_array[j++] = _h1_b_pt;
  _h1_b_eta   = (TH1F*)_hfile->Get("b_eta");   _h1_array[j++] = _h1_b_eta;
  _h1_b_phi   = (TH1F*)_hfile->Get("b_phi");   _h1_array[j++] = _h1_b_phi;
  _h1_b1_pt   = (TH1F*)_hfile->Get("b1_pt");   _h1_array[j++] = _h1_b1_pt;
  _h1_b2_pt   = (TH1F*)_hfile->Get("b2_pt");   _h1_array[j++] = _h1_b2_pt;
  _h1_b3_pt   = (TH1F*)_hfile->Get("b3_pt");   _h1_array[j++] = _h1_b3_pt;
  _h1_b4_pt   = (TH1F*)_hfile->Get("b4_pt");   _h1_array[j++] = _h1_b4_pt;
  _h1_b_drmin = (TH1F*)_hfile->Get("b_drmin"); _h1_array[j++] = _h1_b_drmin;
  _h1_c_n     = (TH1F*)_hfile->Get("c_n");     _h1_array[j++] = _h1_c_n;
  _h1_c_pt    = (TH1F*)_hfile->Get("c_pt");    _h1_array[j++] = _h1_c_pt;
  _h1_c_eta   = (TH1F*)_hfile->Get("c_eta");   _h1_array[j++] = _h1_c_eta;
  _h1_c_phi   = (TH1F*)_hfile->Get("c_phi");   _h1_array[j++] = _h1_c_phi;
  _h1_c1_pt   = (TH1F*)_hfile->Get("c1_pt");   _h1_array[j++] = _h1_c1_pt;
  _h1_c2_pt   = (TH1F*)_hfile->Get("c2_pt");   _h1_array[j++] = _h1_c2_pt;
  _h1_c3_pt   = (TH1F*)_hfile->Get("c3_pt");   _h1_array[j++] = _h1_c3_pt;
  _h1_c4_pt   = (TH1F*)_hfile->Get("c4_pt");   _h1_array[j++] = _h1_c4_pt;
  _h1_c_drmin = (TH1F*)_hfile->Get("c_drmin"); _h1_array[j++] = _h1_c_drmin;

  _h1_lp_n_all = (TH1F*)_hfile->Get("lp_n_all"); _h1_array[j++] = _h1_lp_n_all;
  _h1_lp_pt_all = (TH1F*)_hfile->Get("lp_pt_all"); _h1_array[j++] = _h1_lp_pt_all;
  _h1_lp_n    = (TH1F*)_hfile->Get("lp_n");    _h1_array[j++] = _h1_lp_n;
  _h1_lp_pt   = (TH1F*)_hfile->Get("lp_pt");   _h1_array[j++] = _h1_lp_pt;
  _h1_lp_eta  = (TH1F*)_hfile->Get("lp_eta");  _h1_array[j++] = _h1_lp_eta;
  _h1_lp_phi  = (TH1F*)_hfile->Get("lp_phi");  _h1_array[j++] = _h1_lp_phi;
  _h1_lp1_pt  = (TH1F*)_hfile->Get("lp1_pt");  _h1_array[j++] = _h1_lp1_pt;
  _h1_lp2_pt  = (TH1F*)_hfile->Get("lp2_pt");  _h1_array[j++] = _h1_lp2_pt;
  _h1_lp3_pt  = (TH1F*)_hfile->Get("lp3_pt");  _h1_array[j++] = _h1_lp3_pt;
  _h1_lp4_pt  = (TH1F*)_hfile->Get("lp4_pt");  _h1_array[j++] = _h1_lp4_pt;
  _h1_lp5_pt  = (TH1F*)_hfile->Get("lp5_pt");  _h1_array[j++] = _h1_lp5_pt;
  _h1_lp_drmin = (TH1F*)_hfile->Get("lp_drmin"); _h1_array[j++] = _h1_lp_drmin;

  _h1_jet_n   = (TH1F*)_hfile->Get("jet_n");   _h1_array[j++] = _h1_jet_n;
  _h1_jet_pt  = (TH1F*)_hfile->Get("jet_pt");  _h1_array[j++] = _h1_jet_pt;
  _h1_jet_pt2 = (TH1F*)_hfile->Get("jet_pt2"); _h1_array[j++] = _h1_jet_pt2;
  _h1_jet_n_no_e = (TH1F*)_hfile->Get("jet_n_no_e"); _h1_array[j++] = _h1_jet_n_no_e;
  _h1_jet_pt_no_e = (TH1F*)_hfile->Get("jet_pt_no_e"); _h1_array[j++] = _h1_jet_pt_no_e;
  _h1_jet_pt2_no_e = (TH1F*)_hfile->Get("jet_pt2_no_e"); _h1_array[j++] = _h1_jet_pt2_no_e;
  _h1_jet_eta = (TH1F*)_hfile->Get("jet_eta"); _h1_array[j++] = _h1_jet_eta;
  _h1_jet_phi = (TH1F*)_hfile->Get("jet_phi"); _h1_array[j++] = _h1_jet_phi;
  _h1_jet_mjj = (TH1F*)_hfile->Get("jet_mjj"); _h1_array[j++] = _h1_jet_mjj;
  _h1_jet_dphi = (TH1F*)_hfile->Get("jet_dphi"); _h1_array[j++] = _h1_jet_dphi;

  _h1_xsec_dist = 0;

  // make an extra plot for cross section distribution
  if (_h1xsec) {
    int nxs = 0;
    double xsmin=-1, xsmax=0, xs=0;
    for (int i=0; i<_h1xsec->GetNbinsX(); ++i) {
      double yy = _h1xsec->GetBinContent(i+1);
      if (yy>0) {
	xs += yy;
	++nxs;
	if (yy<xsmin || xsmin<0) xsmin=yy;
	if (yy>xsmax) xsmax=yy;
      }
    }
    if (nxs!=0) {
      xs = xs/(double)nxs;
      double epsilon = fabs(xs-0.5*(xsmin+xsmax));
      if (nxs==1) epsilon = xs;
      double x1 = xsmin-epsilon, x2 = xsmax+epsilon;
      _h1_xsec_dist = new TH1F("xsec_dist","Cross section of last event per run",50,x1,x2);
      _h1_xsec_dist->SetStats(1111);
      for (int i=0; i<_h1xsec->GetNbinsX(); ++i) {
	double yy = _h1xsec->GetBinContent(i+1);
	if (yy>0) _h1_xsec_dist->Fill(yy);
      }
      if (nxs>5) _h1_xsec_dist->Fit("gaus");
      _h1_array[j++] = _h1_xsec_dist;
    }
  }

}


void HistoFile::SetColor(int icolor) {
  for (int i=0; i<200; ++i)
    if (_h1_array[i]) _h1_array[i]->SetLineColor(icolor);
}

void HistoFile::Scale(double scale) {
  for (int i=0; i<200; ++i)
    if (_h1_array[i] && _h1_array[i]!=_h1xsec && _h1_array[i]!=_h1_xsec_dist)
      _h1_array[i]->Scale(scale);
}