#define DEBUG 0
#define SAME 0
// SAME 0 needs to be set, if 2 diff. backgrounds are to be used, e.g. OS vs LSg
void plot_sLSa_mOS(void)
{

    gROOT->SetBatch();

    // CHANGE HERE ---------------------------------------
    // file title
    TString sTitle="pol1_Prec"; // must match the next lines!!!

    TString sHist1="pol1_norm";
    TString sHist2="pol1_norm";
    
    /*
      pol0:      N--     : e-e-
      pol1:      OS      : e+e-
      pol2:      N++     : e+e+
      back_0:    LSg     : 2*sqrt(N++ * N--)
      back_1:    LSa     : 0.5*(N-- + N++)
    */

    // ----------------------------------------------------
    TString sCut="cut7_"; // cuts, c.f. pairs.cc
    TString sCutE="cut10_"; // eff offset 16
    //
    // old files
    //    TString sFile1="hist_exp_aug04_101006_1.root";//black
    //    TString sFile2="hist_mixing_aug04_full_1.root";//red

    //new files with correct normalization for aug04
    TString sFile1="aug04_pairs_EXP_FULL_BADhv_010407.root";//black
    TString sFile2="aug04_eventMixing_FULL_badHV.root";//red
    // ---------------------------------------------------

    // NOTHING TO CHANGE
    TString cTitle1="EXP_AUG04_olay_"+sTitle;
    TString cTitle2="EXP_AUG04_ratio_"+sTitle;

    TString cTitle3="EXP_AUG04_olay_"+sTitle+"_eff";
    TString cTitle4="EXP_AUG04_ratio_"+sTitle+"_eff";
    TString cTitle5="EXP_AUG04_all_"+sTitle+"_eff";
    TString cTitle6="EXP_AUG04_ratio_all_"+sTitle+"_eff";

    // NOTHING TO CHANGE
    TString sMass ="hmass_";
    TString sOpang="hoangle_";
    TString sPt   ="hpt_";
    TString sRap  ="hrap_";

    ///////////// BUILD HISTOGRAM NAMES
    // file 1
    TString sHMass1  = sMass  + sCut + sHist1;
    TString sHOpang1 = sOpang + sCut + sHist1;
    TString sHPt1    = sPt    + sCut + sHist1;
    TString sHRap1   = sRap   + sCut + sHist1;

    TString sHMassE1  = sMass  + sCutE + sHist1;
    TString sHOpangE1 = sOpang + sCutE + sHist1;
    TString sHPtE1    = sPt    + sCutE + sHist1;
    TString sHRapE1   = sRap   + sCutE + sHist1;

    // file2
#if SAME
    TString sHMass2  = sHMass1;
    TString sHOpang2 = sHOpang1;
    TString sHPt2    = sHPt1;
    TString sHRap2   = sHRap1;

    TString sHMassE2  = sHMassE1;
    TString sHOpangE2 = sHOpangE1;
    TString sHPtE2    = sHPtE1;
    TString sHRapE2   = sHRapE1;
#else
    TString sHMass2  = sMass  + sCut + sHist2;
    TString sHOpang2 = sOpang + sCut + sHist2;
    TString sHPt2    = sPt    + sCut + sHist2;
    TString sHRap2   = sRap   + sCut + sHist2;

    TString sHMassE2  = sMass  + sCutE + sHist2;
    TString sHOpangE2 = sOpang + sCutE + sHist2;
    TString sHPtE2    = sPt    + sCutE + sHist2;
    TString sHRapE2   = sRap   + sCutE + sHist2;
#endif


TFile file1(sFile1.Data());

 TH1F* h1m = (TH1F*) file1.Get(sHMass1.Data());
 TH1F* h1o = (TH1F*) file1.Get(sHOpang1.Data());
 TH1F* h1p = (TH1F*) file1.Get(sHPt1.Data());
 TH1F* h1r = (TH1F*) file1.Get(sHRap1.Data());

 h1o->Rebin(4);
 h1r->Rebin(2);
 
 TH1F* h1me = (TH1F*) file1.Get(sHMassE1.Data());
 TH1F* h1oe = (TH1F*) file1.Get(sHOpangE1.Data());
 TH1F* h1pe = (TH1F*) file1.Get(sHPtE1.Data());
 TH1F* h1re = (TH1F*) file1.Get(sHRapE1.Data());

// h1oe->Rebin(4);
// h1re->Rebin(2);

TFile file2(sFile2.Data());

 TH1F* h2m = (TH1F*) file2.Get(sHMass2.Data());
 TH1F* h2o = (TH1F*) file2.Get(sHOpang2.Data());
 TH1F* h2p = (TH1F*) file2.Get(sHPt2.Data());
 TH1F* h2r = (TH1F*) file2.Get(sHRap2.Data());

 h2o->Rebin(4);
 h2r->Rebin(2);
 
 TH1F* h2me = (TH1F*) file2.Get(sHMassE2.Data());
 TH1F* h2oe = (TH1F*) file2.Get(sHOpangE2.Data());
 TH1F* h2pe = (TH1F*) file2.Get(sHPtE2.Data());
 TH1F* h2re = (TH1F*) file2.Get(sHRapE2.Data());

// h2oe->Rebin(4);
// h2re->Rebin(2);



// MAKE CANVAS

TCanvas c1(cTitle1.Data(),cTitle1.Data(),800,800);
TCanvas c2(cTitle2.Data(),cTitle2.Data(),800,800);
TCanvas c3(cTitle3.Data(),cTitle3.Data(),800,800);
TCanvas c4(cTitle4.Data(),cTitle4.Data(),800,800);
TCanvas c5(cTitle5.Data(),cTitle5.Data(),800,800);
TCanvas c6(cTitle6.Data(),cTitle6.Data(),800,800);

c1.Divide(2,2);
c2.Divide(2,2);
c3.Divide(2,2);
c4.Divide(2,2);
c5.Divide(2,2);
c6.Divide(2,2);
c1.SetFillColor(10);
c2.SetFillColor(10);
c3.SetFillColor(10);
c4.SetFillColor(10);
c5.SetFillColor(10);
c6.SetFillColor(10);
// NO BIN WIDTH MARKER

gStyle->SetErrorX(0);
gStyle->SetOptStat(0);
gStyle->SetCanvasColor(10);
gStyle->SetPadColor(10);

///////////////////////////////////////////////////////

// INVARIANT MASS

 setHistoProps(h1me,21,1,1,2.e-7,1.e-12);
 setHistoProps(h1m ,21,1,1,5.e-7,1.e-12);


 h1m->GetXaxis()->SetTitle("M [MeV/c^{2}]");
 h1m->GetYaxis()->SetTitle("dN/dM [1 / MeV/c^{2}]");
 h1m->GetYaxis()->SetTitleOffset(1.4);
 h1m->SetTitle("invariant mass");

 h1me->GetXaxis()->SetTitle("M [MeV/c^{2}]");
 h1me->GetYaxis()->SetTitle("dN/dM [1 / MeV/c^{2}]");
 h1me->GetYaxis()->SetTitleOffset(1.4);
 h1me->SetTitle("invariant mass -- eff. corr.");


 Int_t mbins[4];
 Int_t mbinsE[4];

// Normalization range
 checkBinRanges(h1m,h2m,150.,400.,mbins);
 checkBinRanges(h1me,h2me,150.,400.,mbinsE);

// Integrate
 Stat_t fIntegral = h1m->Integral(mbins[0],
				   mbins[1],
				   "width");
 Stat_t fIntegral2 = h2m->Integral(mbins[2],
				     mbins[3],
				     "width");

 Stat_t fIntegralE = h1me->Integral(mbinsE[0],
				   mbinsE[1],
				   "width");
 Stat_t fIntegral2E = h2me->Integral(mbinsE[2],
				     mbinsE[3],
				     "width");

// Normalize mixed event CB to same integral
h2m  -> Scale(fIntegral/fIntegral2);
h2me -> Scale(fIntegralE/fIntegral2E);

// colors and markers for the overlayed, scaled histo
setHistoProps(h2m,22,2,2);
setHistoProps(h2me,22,2,2);

// draw on respective pad on resp. canvas
c1.cd(1);
 gPad->SetLogy(1);
 h1m->Draw();
 h2m->Draw("same");
 setLegend("Comb. Background","same event",h1m,"mixed event",h2m);
c3.cd(1);
 gPad->SetLogy(1);
 h1me->Draw();
 h2me->Draw("same");
 setLegend("Comb. Background","same event",h1me,"mixed event",h2me);

 // RATIO ////////////////////////////////////////////

 TH1F* h1mc  = (TH1F*) h1m->Clone();
 TH1F* h1mec = (TH1F*) h1me->Clone();

 h1mc->Divide(h2m);
 h1mec->Divide(h2me);

 setHistoProps(h1mc  ,21,1,1,3.,0.3);
 setHistoProps(h1mec ,21,1,1,3.,0.3);

 h1mc->GetXaxis()->SetTitle("M [MeV/c^{2}]");
 h1mc->GetYaxis()->SetTitle("ratio s / m");
 h1mc->GetYaxis()->SetTitleOffset(1.4);
 h1mc->SetTitle("inv. mass");

 h1mec->GetXaxis()->SetTitle("M [MeV/c^{2}]");
 h1mec->GetYaxis()->SetTitle("ratio s / m");
 h1mec->GetYaxis()->SetTitleOffset(1.4);
 h1mec->SetTitle("inv. mass - eff. corr.");

c2.cd(1);
 gPad->SetLogy(0);
 h1mc->Draw();
c4.cd(1);
 gPad->SetLogy(0);
 h1mec->Draw();

///////////////////////////////////////////////////////
// OPENING ANGLE

 setHistoProps(h1oe,21,1,1,1.e-8,0.001);
 setHistoProps(h1o ,21,1,1,1.e-7,0.001);


 h1o->GetXaxis()->SetTitle("#alpha [#circ]");
 h1o->GetYaxis()->SetTitle("dN/d#alpha [1 / 4#circ]");
 h1o->GetYaxis()->SetTitleOffset(1.4);
 h1o->SetTitle("opening angle");

 h1oe->GetXaxis()->SetTitle("op.angle #alpha [#circ]");
 h1oe->GetYaxis()->SetTitle("dN/d#alpha [1 / 4#circ]");
 h1oe->GetYaxis()->SetTitleOffset(1.4);
 h1oe->SetTitle("opening angle -- eff. corr.");


 Int_t mbinso[4];
 Int_t mbinsoE[4];

// Normalization range
 checkBinRanges(h1o,h2o,20.,120.,mbinso);
 checkBinRanges(h1oe,h2oe,20.,120.,mbinsoE);

// Integrate
  fIntegral = h1o->Integral(mbinso[0],
				   mbinso[1],
				   "width");
  fIntegral2 = h2o->Integral(mbinso[2],
				     mbinso[3],
				     "width");

  fIntegralE = h1oe->Integral(mbinsoE[0],
				   mbinsoE[1],
				   "width");
  fIntegral2E = h2oe->Integral(mbinsoE[2],
				     mbinsoE[3],
				     "width");

// Normalize mixed event CB to same integral
h2o  -> Scale(fIntegral/fIntegral2);
h2oe -> Scale(fIntegralE/fIntegral2E);

// colors and markers for the overlayed, scaled histo
setHistoProps(h2o,22,2,2);
setHistoProps(h2oe,22,2,2);

// draw on respective pad on resp. canvas
c1.cd(2);
 gPad->SetLogy(0);
 h1o->Draw();
 h2o->Draw("same");
 setLegend("Comb. Background","same event",h1o,"mixed event",h2o);

c3.cd(2);
 gPad->SetLogy(0);
 gPad->SetFillColor(0);
 h1oe->Draw();
 h2oe->Draw("same");
 setLegend("Comb. Background","same event",h1oe,"mixed event",h2oe);


 // RATIO ////////////////////////////////////////////

 TH1F* h1oc  = (TH1F*) h1o->Clone();
 TH1F* h1oec = (TH1F*) h1oe->Clone();

 h1oc->Divide(h2o);
 h1oec->Divide(h2oe);

 setHistoProps(h1oc  ,21,1,1,3.,0.3);
 setHistoProps(h1oec ,21,1,1,3.,0.3);

 h1oc->GetXaxis()->SetTitle("op.angle #alpha [#circ]");
 h1oc->GetYaxis()->SetTitle("ratio s / m");
 h1oc->GetYaxis()->SetTitleOffset(1.4);
 h1oc->SetTitle("opening angle");

 h1oec->GetXaxis()->SetTitle("op.angle #alpha [#circ]");
 h1oec->GetYaxis()->SetTitle("ratio s / m");
 h1oec->GetYaxis()->SetTitleOffset(1.4);
 h1oec->SetTitle("opening angle - eff. corr.");

c2.cd(2);
 gPad->SetLogy(0);
 h1oc->Draw();
c4.cd(2);
 gPad->SetLogy(0);
 h1oec->Draw();



///////////////////////////////////////////////////////
//
// TRANSVERSE MOMENTUM

 setHistoProps(h1pe,21,1,1,1.e-8,1.e-12);
 setHistoProps(h1p ,21,1,1,1.e-8,1.e-12);


 h1p->GetXaxis()->SetTitle("P_{#perp}  [MeV/c]");
 h1p->GetYaxis()->SetTitle("dN/dP_{#perp}  [1/ 4 MeV/c]");
 h1p->GetYaxis()->SetTitleOffset(1.4);
 h1p->SetTitle("trans. pair mom.");

 h1pe->GetXaxis()->SetTitle("P_{#perp}  [MeV/c]");
 h1pe->GetYaxis()->SetTitle("dN/dP_{#perp}  [1/ 4 MeV/c]");
 h1pe->GetYaxis()->SetTitleOffset(1.4);
 h1pe->SetTitle("trans. pair mom. -- eff. corr.");

 Int_t mbinsp[4];
 Int_t mbinspE[4];

// Normalization range
 checkBinRanges(h1p,h2p,40.,600.,mbinsp);
 checkBinRanges(h1pe,h2pe,40.,600.,mbinspE);

// Integrate
  fIntegral = h1p->Integral(mbinsp[0],
				   mbinsp[1],
				   "width");
  fIntegral2 = h2p->Integral(mbinsp[2],
				     mbinsp[3],
				     "width");

  fIntegralE = h1pe->Integral(mbinspE[0],
				   mbinspE[1],
				   "width");
  fIntegral2E = h2pe->Integral(mbinspE[2],
				     mbinspE[3],
				     "width");

// Normalize mixed event CB to same integral
h2p  -> Scale(fIntegral/fIntegral2);
h2pe -> Scale(fIntegralE/fIntegral2E);

// colors and markers for the overlayed, scaled histo
setHistoProps(h2p,22,2,2);
setHistoProps(h2pe,22,2,2);

// draw on respective pad on resp. canvas
c1.cd(3);
 gPad->SetLogy(1);
 h1p->Draw();
 h2p->Draw("same");
 setLegend("Comb. Background","same event",h1p,"mixed event",h2p);

c3.cd(3);
 gPad->SetLogy(1);
 h1pe->Draw();
 h2pe->Draw("same");
 setLegend("Comb. Background","same event",h1pe,"mixed event",h2pe);


 // RATIO ////////////////////////////////////////////

 TH1F* h1pc  = (TH1F*) h1p->Clone();
 TH1F* h1pec = (TH1F*) h1pe->Clone();

 h1pc->Divide(h2p);
 h1pec->Divide(h2pe);

 setHistoProps(h1pc  ,21,1,1,3.,0.3);
 setHistoProps(h1pec ,21,1,1,3.,0.3);

 h1pc->GetXaxis()->SetTitle("P_{#perp}  [MeV/c]");
 h1pc->GetYaxis()->SetTitle("ratio s / m");
 h1pc->GetYaxis()->SetTitleOffset(1.4);
 h1pc->SetTitle("transv. mom.");

 h1pec->GetXaxis()->SetTitle("P_{#perp}  [MeV/c]");
 h1pec->GetYaxis()->SetTitle("ratio s / m");
 h1pec->GetYaxis()->SetTitleOffset(1.4);
 h1pec->SetTitle("transv. mom. - eff. corr.");

c2.cd(3);
 gPad->SetLogy(0);
 h1pc->Draw();
c4.cd(3);
 gPad->SetLogy(0);
 h1pec->Draw();


///////////////////////////////////////////////////////
// RAPIDITY

 setHistoProps(h1re, 21, 1, 1, 8.e-7, 0.0001);
 setHistoProps(h1r , 21, 1, 1, 3.e-6, 0.0001);

 h1r->GetXaxis()->SetTitle("Y");
 h1r->GetXaxis()->SetRangeUser(-0.1,2.1);
 h1r->GetYaxis()->SetTitle("dN/dY");
 h1r->GetYaxis()->SetTitleOffset(1.4);
 h1r->SetTitle("pair rapidity");

 h1re->GetXaxis()->SetTitle("Y");
 h1re->GetXaxis()->SetRangeUser(-0.1,2.1);
 h1re->GetYaxis()->SetTitle("dN/dY");
 h1re->GetYaxis()->SetTitleOffset(1.4);
 h1re->SetTitle("pair rapidity -- eff. corr.");

 Int_t mbinsr[4];
 Int_t mbinsrE[4];

// Normalization range
 checkBinRanges(h1r,h2r,0.1,1.9.,mbinsr);
 checkBinRanges(h1re,h2re,0.1,1.9,mbinsrE);

// Integrate
  fIntegral = h1r->Integral(mbinsr[0],
				   mbinsr[1],
				   "width");
  fIntegral2 = h2r->Integral(mbinsr[2],
				     mbinsr[3],
				     "width");

  fIntegralE = h1re->Integral(mbinsrE[0],
				   mbinsrE[1],
				   "width");
  fIntegral2E = h2re->Integral(mbinsrE[2],
				     mbinsrE[3],
				     "width");

// Normalize mixed event CB to same integral
h2r  -> Scale(fIntegral/fIntegral2);
h2re -> Scale(fIntegralE/fIntegral2E);

// colors and markers for the overlayed, scaled histo
setHistoProps(h2r,22,2,2);
setHistoProps(h2re,22,2,2);

// draw on respective pad on resp. canvas
c1.cd(4);
 gPad->SetLogy(0);
 h1r->Draw();
 h2r->Draw("same");
 setLegend("Comb. Background","same event",h1r,"mixed event",h2r);
c3.cd(4);
 gPad->SetLogy(0);
 h1re->Draw();
 h2re->Draw("same");
 setLegend("Comb. Background","same event",h1re,"mixed event",h2re);





 // RATIO ////////////////////////////////////////////

 TH1F* h1rc  = (TH1F*) h1r->Clone();
 TH1F* h1rec = (TH1F*) h1re->Clone();

 h1rc->Divide(h2r);
 h1rec->Divide(h2re);

 setHistoProps(h1rc  ,21,1,1,3.,0.3);
 setHistoProps(h1rec ,21,1,1,3.,0.3);

 
 h1rc->GetXaxis()->SetTitle("Y");
 h1rc->GetXaxis()->SetRangeUser(-0.1,2.1);
 h1rc->GetYaxis()->SetTitle("ratio s / m");
 h1rc->GetYaxis()->SetTitleOffset(1.4);
 h1rc->SetTitle("rapidity");

 h1rec->GetXaxis()->SetTitle("Y");
 h1rec->GetXaxis()->SetRangeUser(-0.1,2.1);
 h1rec->GetYaxis()->SetTitle("ratio s / m");
 h1rec->GetYaxis()->SetTitleOffset(1.4);
 h1rec->SetTitle("rapidity - eff. corr.");

c2.cd(4);
 gPad->SetLogy(0);
 gPad->SetFillColor(10);
 h1rc->Draw();
c4.cd(4);
 gPad->SetLogy(0);
 h1rec->Draw();



 c5.cd(1);
 gPad->SetFillColor(10);
 gPad->SetLogy();
 h1m->DrawCopy();
 h2m->DrawCopy("same");
 setLegend("e+e-","same event",h1m,"mixed event",h2m);
 c5.cd(2);
 h1oe->SetTitle("M_{ee} > 0.15 GeV/c^{2}");
 h1oe->DrawCopy();
 h2oe->DrawCopy("same");
 c5.cd(3);
 gPad->SetLogy();
 h1pe->SetTitle("M_{ee} > 0.15 GeV/c^{2}");
 h1pe->DrawCopy();
 h2pe->DrawCopy("same");
 c5.cd(4);
// gPad->SetLogy();
 h1re->SetTitle("M_{ee} > 0.15 GeV/c^{2}");
 h1re->DrawCopy();
 h2re->DrawCopy("same");


 c6.cd(1);
 gPad->SetFillColor(10);
 h1mc->DrawCopy();
 c6.cd(2);
 h1oec->SetTitle("M_{ee} > 0.15 GeV/c^{2}");
 h1oec->DrawCopy();
 c6.cd(3);
 h1pec->SetTitle("M_{ee} > 0.15 GeV/c^{2}");
 h1pec->DrawCopy();
 c6.cd(4);
// gPad->SetLogy();
 h1rec->SetTitle("M_{ee} > 0.15 GeV/c^{2}");
 h1rec->DrawCopy();


 // --------------------------------------------------
 // ------- OUTPUT -----------------------------------

 TString cTitle1c = cTitle1;
 TString cTitle2c = cTitle2;
 TString cTitle3c = cTitle3;
 TString cTitle4c = cTitle4;
 TString cTitle5c = cTitle5;
 TString cTitle6c = cTitle6;

 cTitle1.Append(".ps");
 cTitle2.Append(".ps");
 cTitle3.Append(".ps");
 cTitle4.Append(".ps");
 cTitle5.Append(".ps");
 cTitle6.Append(".ps");

 cTitle1c.Append(".gif");
 cTitle2c.Append(".gif");
 cTitle3c.Append(".gif");
 cTitle4c.Append(".gif");
 cTitle5c.Append(".gif");
 cTitle6c.Append(".gif");

c1.Print(cTitle1.Data());
c2.Print(cTitle2.Data());
c3.Print(cTitle3.Data());
c4.Print(cTitle4.Data());
c5.Print(cTitle5.Data());
c6.Print(cTitle6.Data());

c1.Print(cTitle1c.Data());
c2.Print(cTitle2c.Data());
c3.Print(cTitle3c.Data());
c4.Print(cTitle4c.Data());
c5.Print(cTitle5c.Data());
c6.Print(cTitle6c.Data());

TFile *out = new TFile("output.root","recreate");
c1.Write();
c2.Write();
c3.Write();
c4.Write();
c5.Write();
c6.Write();
out->Close();

}

void checkHistoDims(TH1F* h1, TH1F* h2)
{
    if (h1->GetNbinsX()!=h2->GetNbinsX()) 
    {
	cout << "ERROR: histogram dimensions different" << endl;
	cout << "h1: " << h1->GetNbinsX() << "   "
	     << "h2: " << h2->GetNbinsX() << endl;
	exit(1);
    }
#if DEBUG    
    else cout << "Histo dimensions equal!" << endl;
#endif
}

void checkRange(TH1F* h,Int_t nL, Int_t nU)
{
    // Lower bin edge included, upper bin edge not included !
    if (nL<1) 
    {
	cout << "ERROR: "
	     <<h->GetTitle()
	  <<"  : underflow bin "
	     <<nL<<" included !" <<endl;
	exit(1);
    }
    if (nU>h->GetNbinsX()) 
    {
	cout << "ERROR: "
	     <<h->GetTitle()
	  <<"  : overflow bin "
	     <<nU<<" included !" <<endl;
	exit(1);
    }
}


void setHistoProps(TH1F* h,Int_t nMarkerStyle=21, Int_t nMarkerColor=1, 
		   Int_t nLineColor=1,
		   Float_t fMax=1.e-6,Float_t fMin=1.e-10)
{
    gStyle->SetCanvasColor(10);
    gStyle->SetPadColor(10);
    h->SetMarkerStyle(nMarkerStyle);
    h->SetMarkerSize(0.8);
    h->SetMaximum(fMax);
    h->SetMinimum(fMin);
    h->SetMarkerColor(nMarkerColor);
    h->SetLineColor(nLineColor);

    h->GetXaxis()->SetLabelSize(0.03);
    h->GetYaxis()->SetLabelSize(0.03);
    h->GetXaxis()->SetLabelColor(1);
    h->GetYaxis()->SetLabelColor(1);
    h->GetXaxis()->SetLabelOffset(0.01);
    h->GetYaxis()->SetLabelOffset(0.01);
    h->GetXaxis()->SetTitleColor(1);
    h->GetYaxis()->SetTitleColor(1);
    h->SetFillColor(10);
    gPad->SetTopMargin(0.2);
    gPad->SetFillColor(10);
    gPad->SetBottomMargin(0.15);
    gPad->SetLeftMargin(0.15);


}
 void checkBinRanges(TH1F* h1,TH1F* h2, Float_t fLeft, Float_t fRight,
		     Int_t* nBins)

     {
#if DEBUG 	 
	 printHisto(h1);
	 printHisto(h2);
#endif

	 nBins[0] = h1->FindBin(fLeft);
	 nBins[1] = h1->FindBin(fRight);
	 nBins[2] = h2->FindBin(fLeft);
	 nBins[3] = h2->FindBin(fRight);
#if DEBUG
	 cout << "left :\t "
	      << nBins[0] << " - " << nBins[1] << endl;
	 cout << "right :\t "
	      << nBins[2] << " - " << nBins[3] <<endl;
#endif
	 checkHistoDims(h1,h2);
	 checkRange(h1,nBins[0],nBins[1]);
	 if (nBins[0]!=nBins[2] ||
	     nBins[1]!=nBins[3]   )
	 {
	     cout << "ERROR: different bin limits" << endl;
	     exit(1);
	 }
     }

 void printHisto(TH1F* h)
 {
     // print contents of a histogram
     cout<<"Title of Histogram: "<<h->GetTitle()<<endl;
     Int_t nbXbins = h->GetNbinsX();
     cout<<"Number of bins in histo: "<<nbXbins<<endl;
     for (Int_t j=0;j<nbXbins+2;j++) // loop over bins in X dir
     {
 	Int_t i = h->GetBin(j); // return bin, redundant here
 	Float_t cc = h->GetBinCenter(i);
 	Float_t errn= h->GetBinError(i);
 	Float_t width=h->GetBinWidth(i);
 	Float_t con=h->GetBinContent(i);
 	//	if (cc > 1000.) break; 
 	Float_t relerr=-1.;
 	if (con>0.) relerr = errn/con*100.;
  	cout<<"Index: "<<j<<" Bin: "<<i<<" Center: "<<cc<<" Content:  "
 	    <<con<<" Width: "<<width<<" Error: "<<errn<<" rel. Err.: "
 	    <<relerr<<"%"<<" (sqrt content: "<<TMath::Sqrt(con)<<")"<<endl;
     }
 }


void setLegend( char* header, char* cH1, TH1F* h1, char* cH2, TH1F* h2)
{
	 TLegend *legmass = new TLegend(0.5,0.7,1.0,1.0);
	 legmass->SetFillColor(0);
	 legmass->SetTextSize(0.046);
	 legmass->SetTextFont(42);
	 legmass->SetHeader(header);
	 if (h1) legmass->AddEntry(h1,cH1);
	 if (h2) legmass->AddEntry(h2,cH2);
	 legmass->Draw("same");
}