#include "TCanvas.h"
#include "TF1.h"
#include "TFile.h"
#include "TH1.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TLatex.h"
#include "TLine.h"
#include "TObjArray.h"
#include "TObjString.h"
#include "TPad.h"
#include "TROOT.h"
#include "TString.h"
#include "TTree.h"
#include "TArrow.h"
#include "TLegend.h"
#include "TEventList.h"
#include "TStyle.h"

#include <iostream>
#include <vector>
#include <algorithm>

TString savepath   = "fig/ananote/";
TString savepathqa = savepath+"scans/"; 
TString savepref   = "";

int cnt=0;

TString sigcut = "mode>9";
TString bkgcut1 = "mode==0";
TString bkgcut2 = "mode==2";

double S0e   = 98000;
double B0e   = 9580637576;
double B0e2  = 100000;

double BRe   = 0.05971;    // BR(J/psi -> e+ e-)

double S0m   = 100000;
double B0m   = 8871910789;
double B0m2  = 99000;
double B0maf = 2600283835;

double BRm   = 0.0596;     // BR(J/psi -> mu+ mu-)

double SigXs  = 100;         // [nb]
double BkgXs  = 46e6;       // [nb] inelastic
double BkgXst = 59e6;       // [nb] total
double BkgNR  = 1.2;        // [nb] ppbar -> J/psi pi+ pi-

double BRX = 0.05;         // BR(X -> J/psi pi+ pi-)


typedef std::vector<TString> StrVec;

void PandaSmartLabel(TString align="Rprel", Double_t xpos=-1, Double_t ypos=-1, Double_t scale=-1, TString str="-1", Double_t scale2=-1) 
{
  // H as the first parameter triggers printing comamnd options
  if (align == "H" || gPad==0) {
   cout 
      << endl
      << "  USAGE: PandaSmartLabel(\"align\",xpos,ypos,scale,\"str\",scale2);" << endl
      << "  Prints the official BaBar label on the active ROOT pad" << endl
      << "      align   R or L, optionally with additional 'prel' or 'front': str is aligned to the right (default) or left edge of the \"PANDA\" label" << endl
      << "      xpos    X position of the \"PANDA\" label\'s bottom right (align==R)/left (==L) corner, 0 < xpos < 1, defaults to (L) left margin + 0.02 or (R) 1.0 - right margin - 0.02" << endl 
      << "      ypos    Y position of the \"PANDA\" label\'s bottom right (align==R)/left (==L) corner, 0 < ypos < 1, defaults to 1 - top margin - 0.08" << endl
      << "      scale   relative size of the label, defaults to 0.9" << endl 
      << "      str     LaTeX-style text that goes under the PANDA label." << endl
      << "      scale2  relative size of the second line of text, defaults to 0.45" << endl;
   cout <<endl
	  << "    Examples: " << endl
      << "      PandaSmartLabel(\"H\")         -> Prints out this help text" << endl
      << "      PandaSmartLabel()            -> Label 'PANDA//MC simulation' in top right corner" << endl
      << "      PandaSmartLabel(\"L\")         -> Label 'PANDA//MC simulation' in top left corner" << endl
      << "      PandaSmartLabel(\"-1\",-1,-1,-1,\"preliminary\",0.4) -> Label 'PANDA//preliminary' in top right corner " << endl
      << "      PandaSmartLabel(\"Rprel\")     -> Label 'PANDA' in top right corner  with a large 'PRELIMINARY' across the plot behind histograms (works only with gStyle->SetFrameFillStyle(0))" << endl
      << "      PandaSmartLabel(\"Lprelfront\")-> Label 'PANDA' in top left corner  with a large 'PRELIMINARY' across the plot behind histograms (works only with gStyle->SetFrameFillStyle(0))" << endl << endl;
    return;
  }

  if (gPad==0) return;
	
  // Making -1 a placeholder for function's default value 
  if (align=="") align="R";
  
  bool rightalign = !(align.Contains("L"));
  bool preliminary = align.Contains("prel");
  bool front = align.Contains("front");
  
  double mt = gPad->GetTopMargin();
  double mb = gPad->GetBottomMargin();
  double ml = gPad->GetLeftMargin();
  double mr = gPad->GetRightMargin();
  
  if (ypos==-1) ypos = 1.0 - mt - 0.08;
  
  if (rightalign) // Right align 
  {
	  if (xpos==-1) xpos = 1.0 - mr - 0.02;
  }
  else
  {							// Left align
	  if (xpos==-1) xpos = ml + 0.02;
  }
  
  if (scale == -1) scale = 0.9; 
  if (str == "-1") str = "MC simulation"; 
  if (scale2 == -1) scale2 = 0.43;
  

  // Draw the label 
  TLatex *panda = new TLatex();
  panda->SetNDC(kTRUE); 
   double fcx=1., fcy=1., dx=1.0-xpos, dy = 1.0-ypos;
//   
  double padSizeX = gPad->GetWw()*gPad->GetWNDC(); // Get pad's dimensions
  double padSizeY = gPad->GetWh()*gPad->GetHNDC();
//   if (padSizeX>padSizeY) fcx=padSizeY/padSizeX;
//   if (padSizeX<padSizeY) fcy=padSizeX/padSizeY;
//   
//   cout <<padSizeX<<" "<<padSizeY<<" : "<<fcx<<" "<<fcy<<endl;
//   
  panda->SetNDC(kTRUE); 

  if (preliminary)
  {
	  panda->SetTextFont(42);
	  panda->SetTextColor(17);
	  panda->SetTextSize(0.17*(1.-ml-mr));
	  panda->SetTextAlign(22);
	  double sizeX = padSizeX*(1. - ml - mr);
	  double sizeY = padSizeY*(1. - mt - mb);
	  double posX  = (ml+(1.-mr))*0.5;
	  double posY  = (mb+(1.-mt))*0.42;
	  panda->SetTextAngle(atan(sizeY/sizeX)/3.1415*142.);
	  panda->DrawLatex(posX,posY,"PRELIMINARY");
	  
	  TList *ll=gPad->GetListOfPrimitives();
	  if (!front) ll->AddAt(ll->RemoveAt(ll->GetEntries()-1),0);
  }
  
  panda->SetTextColor(13);
  panda->SetTextAngle(0);
  
  panda->SetTextFont(62); // Bold helvetiva
  if (rightalign) panda->SetTextAlign(31); // Right-Bottom
  else panda->SetTextAlign(11);            // Left bottom
  panda->SetTextSize(0.07*scale);
  panda->DrawLatex(xpos*fcy,ypos*fcx,"#bar{P}ANDA");
  
  panda->SetTextFont(42); // Helvetica
  if (rightalign) panda->SetTextAlign(33); // Right-Bottom
  else panda->SetTextAlign(13);            // Left bottom
  panda->SetTextSize(0.08*scale2);
  panda->DrawLatex(xpos*fcx,ypos*fcy*0.99,str);
 
  delete panda;
  
  return;
}

// --------------------------------------------------------------------

int SplitString(TString s, TString delim, StrVec &toks)
{
	toks.clear();
	TObjArray *tok = s.Tokenize(delim);
	int N = tok->GetEntries();	
	for (int i=0;i<N;++i) 
	{
		TString token = (((TObjString*)tok->At(i))->String()).Strip(TString::kBoth);
		toks.push_back(token);
	}
	return toks.size();
}
void setPandaStyle(TString opt="")
{
	//TStyle *pandaStyle=new TStyle("PANDA","PANDA approved plots style");

	// use plain black on white colors
	gStyle->SetFrameBorderMode(0);
	gStyle->SetCanvasBorderMode(0);
	gStyle->SetPadBorderMode(0);
	gStyle->SetPadColor(0);
	gStyle->SetCanvasColor(0);
	gStyle->SetFrameFillColor(0);
	gStyle->SetFrameFillStyle(0);
	//gStyle->SetFillColor(0);
	
	// set the paper & margin sizes
	//double rightmargin = 0.05;
	//if (opt=="colz") rightmargin=0.12;

	gStyle->SetPaperSize(20,26);
	gStyle->SetPadTopMargin(0.05);
	gStyle->SetPadRightMargin(0.05);
	gStyle->SetPadBottomMargin(0.14);
	gStyle->SetPadLeftMargin(0.13);

	gStyle->SetStatColor(0);
	gStyle->SetOptStat(1110) ; 
	gStyle->SetStatBorderSize(1);
	gStyle->SetStatFont(42);
	gStyle->SetStatFontSize(0.04);
	gStyle->SetStatX(0.95);
	gStyle->SetStatY(0.95);
	gStyle->SetStatW(0.16);
	gStyle->SetStatH(0.18);


	// use sans serif fonts
// 	gStyle->SetTextFont(132);
// 	gStyle->SetTextSize(0.08);
	gStyle->SetLabelFont(42,"x");
	gStyle->SetLabelFont(42,"y");
	gStyle->SetLabelFont(42,"z");
	gStyle->SetLabelSize(0.05,"x");
	gStyle->SetTitleSize(0.06,"x");
	gStyle->SetTitleOffset(1.0,"X");
	gStyle->SetLabelSize(0.05,"y");
	gStyle->SetTitleSize(0.06,"y");
	gStyle->SetTitleOffset(1.,"y");
	gStyle->SetLabelSize(0.05,"z");
	gStyle->SetTitleSize(0.06,"z");
	
	gStyle->SetNdivisions(505,"X");
	
	gStyle->SetStripDecimals(false);

	// use bold lines and markers
	gStyle->SetMarkerStyle(20);
//	gStyle->SetHistLineWidth(1.85);
	//gStyle->SetLineStyleString(2,"[12 12]"); // postscript dashes

	// get rid of X error bars and y error bar caps
	gStyle->SetErrorX(0.001);
	
	gStyle->SetHistMinimumZero(1);
	
	// do not display any of the standard histogram decorations
	gStyle->SetOptTitle(0);
	gStyle->SetOptStat(1110);
	gStyle->SetOptFit(1011);
	gStyle->SetFitFormat(".3g");
		
	gStyle->SetPalette(1);

	// put tick marks on top and RHS of plots
 	gStyle->SetPadTickX(1);
 	gStyle->SetPadTickY(1);
	gStyle->SetTickLength(0.020,"xz");
	gStyle->SetTickLength(0.015,"y");
	gStyle->SetLineScalePS(2);

	//pandaStyle->SetLabelOffset(0.02,"xyz");

	//gROOT->Reset();
	//pandaStyle->SetLineWidth(2);
	//gROOT->SetStyle("PANDA");
	//gROOT->ForceStyle();
}


// --------------------------------------------------------------------
// saves a canvas in various types
void saveCanvas(TCanvas *c, TString name, TString ftypes="gif,C,root,pdf")
{
	//cout <<savepath<<endl;
	//return;
	StrVec types;
	SplitString(ftypes,",",types);
	
	for (unsigned int i=0;i<types.size();++i)
		c->SaveAs(savepath+savepref+name+"."+types[i]);
}

// --------------------------------------------------------------------

void confhist(TH1 *h, int col=1, TString tx="m [GeV/c^{2}]", TString ty="", double offy=1.5, int ndiv=505)
{
	double width = (h->GetXaxis()->GetXmax()-h->GetXaxis()->GetXmin())/h->GetNbinsX()*1000.;
	
	//if (width<0.001) 
	//h->SetNdivisions(500);
	
	if (ty=="") ty=TString::Format("counts / %.1f MeV/c^{2}",width);
	
	h->SetXTitle(tx);
	h->SetYTitle(ty);
	h->GetYaxis()->SetTitleOffset(offy);
	h->SetLineColor(col);
	h->SetStats(0);
	h->SetLineWidth(2);
	h->GetXaxis()->SetNdivisions(ndiv);
}

// --------------------------------------------------------------------
// draws cut mass window
void drawwin(TH1* h, double c, double w)
{
	double max = h->GetMaximum();
	double min = h->GetMinimum();
	
	TLine l;
	l.SetLineColor(6);
	l.SetLineStyle(7);
	l.SetLineWidth(2);
	
	l.DrawLine(c-w, min, c-w, max*0.8);
	l.DrawLine(c+w, min, c+w, max*0.8);
}

// --------------------------------------------------------------------

double writeLatex(double x, double y, TString txt, int col=0, double size=0.05, int font=42)
{
	TLatex lat(x,y,txt);
	lat.SetNDC();
	
	lat.SetTextColor(col);
	lat.SetTextFont(font);
	lat.SetTextSize(size);
	
	lat.DrawLatex(x,y,txt);
	
	return lat.GetYsize();
}

// --------------------------------------------------------------------
// draws signal and background distribution for all individual cuts in a selection cut
void prepHisto(TH1F* hs, TString tit="", TString xtit="", int col=1)
{
	hs->GetXaxis()->SetTitleSize(0.05);
	hs->GetXaxis()->SetLabelSize(0.05);
	hs->GetYaxis()->SetLabelSize(0.05);
	
	hs->SetTitle(tit); 
	hs->SetXTitle(xtit);
	
	hs->SetStats(0);
	
	hs->SetLineColor(col);
	hs->SetLineWidth(2);
	
	hs->SetNdivisions(505);
}

void drawVars(TTree *t, TString cut, TCanvas *c, TString precut="", bool drawcut=true)
{
	if (c==0) return;
	
	// separate individual cuts
	StrVec v;
	SplitString( cut, "&&", v);
		
	int nvar = v.size();
	
	// new divison of TCanvas according to number of variables
	//c->Clear();
	//int xdim = 700, ydim=(nvar/2+nvar%2)*350;//, ydim = 700;
	//c->SetWindowSize(xdim, ydim);
	c->Divide(2,nvar/2+nvar%2);
	
	t->SetEventList(0);
	t->Draw(">>el",precut);
	TEventList *el=(TEventList*)gDirectory->Get("el");
	t->SetEventList(el);
	
	
	for (int i=0;i<nvar;++i)
	{
		c->cd(i+1);
		
		bool gt=true;
		int cutdir = v[i].Contains("<") ? -1 : 1 ;
		int pos = v[i].Index(">");
		if (pos<0) {pos=v[i].Index("<");gt=false;}
		if (pos<0) continue;
		
		TString vi = v[i], var = vi(0,pos), sval =vi(pos+1,1000);
		
		double  cutval = sval.Atof();
		
		double xmin = t->GetMinimum(var)-0.001;
		double xmax = t->GetMaximum(var)+0.001;

		cout <<var<<(gt?" > ":" < ")<<cutval<<"  min="<<xmin<<"  max="<<xmax<<endl;
		
		if (xmin==-0.001 && xmax==0.001)
		{
			xmin=1e8; xmax=-1e8;
			t->Draw(var+Form(">>h%d(100)",cnt++),sigcut);
			TH1F *h=(TH1F*)gROOT->FindObject(Form("h%d",cnt-1));
			if (h && h->GetXaxis()->GetXmin()<xmin) xmin = h->GetXaxis()->GetXmin();
			if (h && h->GetXaxis()->GetXmax()>xmax) xmax = h->GetXaxis()->GetXmax();
			
			t->Draw(var+Form(">>h%d(100)",cnt++),bkgcut1);
			h=(TH1F*)gROOT->FindObject(Form("h%d",cnt-1));
			if (h && h->GetXaxis()->GetXmin()<xmin) xmin = h->GetXaxis()->GetXmin();
			if (h && h->GetXaxis()->GetXmax()>xmax) xmax = h->GetXaxis()->GetXmax();
			
			t->Draw(var+Form(">>h%d(100)",cnt++),bkgcut2);
			h=(TH1F*)gROOT->FindObject(Form("h%d",cnt-1));
			if (h && h->GetXaxis()->GetXmin()<xmin) xmin = h->GetXaxis()->GetXmin();
			if (h && h->GetXaxis()->GetXmax()>xmax) xmax = h->GetXaxis()->GetXmax();
			
			xmin -= (xmax-xmin)*.01;
			xmax += (xmax-xmin)*.01;
			
		}

		t->Draw(var+Form(">>h%d(100,%f,%f)",cnt++,xmin,xmax),sigcut);
		t->Draw(var+Form(">>h%d(100,%f,%f)",cnt++,xmin,xmax),bkgcut1);
		t->Draw(var+Form(">>h%d(100,%f,%f)",cnt++,xmin,xmax),bkgcut2);

		TH1F *hsig=(TH1F*)gDirectory->Get(Form("h%d",cnt-3));
		TH1F *hbkg=(TH1F*)gDirectory->Get(Form("h%d",cnt-2));
		TH1F *hbkg2=(TH1F*)gDirectory->Get(Form("h%d",cnt-1));
		
		prepHisto(hsig, var, var,1);
		prepHisto(hbkg, var, var,2);
		prepHisto(hbkg2, var, var,4);
				
		hsig->Scale(1./hsig->Integral(1,100));
		hbkg->Scale(1./hbkg->Integral(1,100));
		hbkg2->Scale(1./hbkg2->Integral(1,100));
		
		double maxi = 1.05*max(max(hsig->GetMaximum(), hbkg->GetMaximum()),hbkg2->GetMaximum());
		hsig->SetMaximum(maxi);
		hbkg->SetMaximum(maxi);
		hbkg2->SetMaximum(maxi);

		hsig->Draw(); hbkg->Draw("same");hbkg2->Draw("same");
		
		// prepare line and arrow for drawing cut 
		TLine ln;
		ln.SetLineColor(6);	ln.SetLineStyle(7); ln.SetLineWidth(2);
		
		TArrow arr;
		arr.SetAngle(40); arr.SetLineWidth(2); arr.SetLineColor(6); arr.SetFillColor(6);

		if (drawcut) 
		{
			ln.DrawLine(cutval,0,cutval,maxi*0.9);
			arr.DrawArrow(cutval,maxi*0.9,cutval+0.1*(xmax-xmin)*cutdir,maxi*0.9,0.01,"|>");
		}
	}
	
	t->SetEventList(0);
}



// --------------------------------------------------------------------
// draw variables sig/bkg with cuts
//void drawVars(TTree *t, TString cut, TCanvas *c, TString precut="", bool drawcut=true)
//{
	//if (c==0) return;
	
	//// separate individual cuts
	//StrVec v;
	//SplitString( cut, "&&", v);
		
	//int nvar = v.size();
	
	//// new divison of TCanvas accroding to number of variables
	//c->Clear();
	//int xdim = (nvar/2+nvar%2)*350, ydim = 700;
	//c->SetWindowSize(xdim, ydim);
	//c->Divide(nvar/2+nvar%2,2);
	
	//t->SetEventList(0);
	//t->Draw(">>el",precut);
	//TEventList *el=(TEventList*)gDirectory->Get("el");
	//t->SetEventList(el);
	
	//for (int i=0;i<nvar;++i)
	//{
		//c->cd(i+1);
		
		//bool gt=true;
		//int cutdir = v[i].Contains("<") ? -1 : 1 ;
		//int pos = v[i].Index(">");
		//if (pos<0) {pos=v[i].Index("<");gt=false;}
		//if (pos<0) continue;
		
		//TString vi = v[i], var = vi(0,pos), sval =vi(pos+1,1000);
		//double  cutval = sval.Atof();
		
		//double xmin = t->GetMinimum(var)-0.001;
		//double xmax = t->GetMaximum(var)+0.001;

		//cout <<var<<(gt?" > ":" < ")<<cutval<<"  min="<<xmin<<"  max="<<xmax<<endl;
		
		//t->Draw(var+Form(">>h%d(100,%f,%f)",cnt++,xmin,xmax),sigcut);
		//t->Draw(var+Form(">>h%d(100,%f,%f)",cnt++,xmin,xmax),bkgcut);
		
		//TH1F *hsig=(TH1F*)gDirectory->Get(Form("h%d",cnt-2));
		//TH1F *hbkg=(TH1F*)gDirectory->Get(Form("h%d",cnt-1));
		
		//prepHisto(hsig, var, var,1);
		//prepHisto(hbkg, var, var,2);
				
		//hsig->Scale(1./hsig->Integral(1,100));
		//hbkg->Scale(1./hbkg->Integral(1,100));
		
		//double maxi = 1.05*max(hsig->GetMaximum(), hbkg->GetMaximum());
		//hsig->SetMaximum(maxi);
		//hbkg->SetMaximum(maxi);

		//hsig->Draw(); hbkg->Draw("same");
		
		//// prepare line and arrow for drawing cut 
		//TLine ln;
		//ln.SetLineColor(6);	ln.SetLineStyle(7); ln.SetLineWidth(2);
		
		//TArrow arr;
		//arr.SetAngle(40); arr.SetLineWidth(2); arr.SetLineColor(6); arr.SetFillColor(6);

		//if (drawcut) 
		//{
			//ln.DrawLine(cutval,0,cutval,maxi*0.9);
			//arr.DrawArrow(cutval,maxi*0.9,cutval+0.1*(xmax-xmin)*cutdir,maxi*0.9,0.01,"|>");
		//}
	//}
	
	//t->SetEventList(0);
//}

// --------------------------------------------------------------------
// counts different events for a cut
int cntEvents(TTree *t, TString cut)
{
	int bev;
	t->SetEventList(0);
	t->Draw(">>el",cut);
	TEventList *el=(TEventList*)gROOT->FindObject("el");
	
	t->SetBranchStatus("*",0);
	t->SetBranchStatus("ev",1);
	t->SetBranchAddress("ev",&bev);
	
	int lev = -1, cnt=0;
	for (int i=0;i<el->GetN();++i)
	{
		t->GetEvent(el->GetEntry(i));
		if (lev!=bev) cnt++;
		lev = bev;
	}
	
	t->SetBranchStatus("*",1);
	
	return cnt;
}

// --------------------------------------------------------------------
// draws signal and background distribution for all individual cuts in a selection cut
void varsQa(TTree *t, TString var, TString cut, TCanvas *c, TString precut="", int bins=50)
{
	if (c==0) return;
	
	// separate individual cuts
	StrVec v;
	SplitString( cut, "&&", v);
		
	int nvar = v.size();
	
	// new divison of TCanvas accroding to number of variables
	//c->Clear();
	//int xdim = 700, ydim=(nvar/2+nvar%2)*350;//, ydim = 700;
	//c->SetWindowSize(xdim, ydim);
	c->Divide(2,nvar/2+nvar%2);

	//int xdim = (nvar/2+nvar%2)*350, ydim = 700;
	//c->SetWindowSize(xdim, ydim);
	//c->Divide(nvar/2+nvar%2,2);
	
	t->SetEventList(0);
	t->Draw(">>el",precut);
	TEventList *el=(TEventList*)gDirectory->Get("el");
	t->SetEventList(el);
	
	for (int i=0;i<nvar;++i)
	{
		c->cd(i+1);
		double xmin = t->GetMinimum(var)+0.001;
		double xmax = t->GetMaximum(var)-0.001;
		
		t->Draw(var+Form(">>h%d(%d,%f,%f)",cnt++,bins,xmin,xmax),"","e");
		t->Draw(var+Form(">>h%d(%d,%f,%f)",cnt++,bins,xmin,xmax),v[i],"same e");
		TH1F *h    = (TH1F*)gDirectory->Get(Form("h%d",cnt-2));
		TH1F *hcut = (TH1F*)gDirectory->Get(Form("h%d",cnt-1));
		
		prepHisto(h,   v[i],var);
		prepHisto(hcut,v[i],var,2);
	}
}

// --------------------------------------------------------------------

void xplots_final()
{
	savepath   = "fig/ananote/";
	savepathqa = savepath+"scans/"; 
	savepref   = "";

	//gROOT->LoadMacro("setPandaStyle.C");
	setPandaStyle();
	
	gStyle->SetTitleSize(0.05,"x");
	gStyle->SetTitleSize(0.05,"y");
	gStyle->SetPadTopMargin(0.05);
	gStyle->SetPadRightMargin(0.04);
	gStyle->SetPadBottomMargin(0.14);
	gStyle->SetPadLeftMargin(0.16);
	gStyle->SetTitleOffset(1.2,"x");
	//gStyle->SetLineScalePS(1);
	
//gStyle->SetPadBorderMode(0);
//gStyle->SetPadBorderSize(0);
//gStyle->SetPadColor(0);

//gStyle->SetDrawBorder(0);
//gStyle->SetCanvasBorderMode(0);
//gStyle->SetCanvasBorderSize(0);
//gStyle->SetCanvasColor(0);

//gStyle->SetFrameBorderMode(0);
//gStyle->SetFrameBorderSize(0);
//gStyle->SetFrameFillColor(0);
//gStyle->SetFrameFillStyle(0);

	
	TF1 *f1=new TF1("f1","gaus(0)");
	f1->SetLineColor(4);
	f1->SetLineWidth(2);

	// canvas to show J/psi mass for e+ e- case
	TCanvas *c1=new TCanvas("c1","c1",10,10,1500,500);
	c1->Divide(3,1);
	
	// canvas to show J/psi mass for mu+ mu- case
	TCanvas *c2=new TCanvas("c2","c2",200,200,1500,500);
	c2->Divide(3,1);
	
	////canvases for displaying the variables
	//TCanvas *c3=new TCanvas("c3","c3",600,900); // e e
	//TCanvas *c4=new TCanvas("c4","c4",600,600); // mu mu
	
	//// canvases for showing individual cut quality
	//TCanvas *c5=new TCanvas("c5","c5",600,900); // e e
	//TCanvas *c6=new TCanvas("c6","c6",600,600); // mu mu
	
	//canvas for filter quality check
	TCanvas *c7=new TCanvas("c7","c7",500,500);  // mu mu
	
	//// canvas to show 4C fitted pbarp mass both channels
	//TCanvas *c8=new TCanvas("c8","c8",10,10,1000,500);
	//c8->Divide(2,1);
	
	// canvas to show correlation J/psi and rho mass
	TCanvas *c9=new TCanvas("c9","c9",200,200,1500,500);
	c9->Divide(3,1);
	
	TCanvas *c9a=new TCanvas("c9a","c9a",220,2200,1500,500);
	c9a->Divide(3,1);
	
	TCanvas *c9b=new TCanvas("c9b","c9b",240,200,1500,500);
	c9b->Divide(3,1);
	
	TCanvas *c9c=new TCanvas("c9c","c9c",260,200,1500,500);
	c9c->Divide(3,1);
	
	// canvas to show correlation J/psi and rho mass
	//TCanvas *c10=new TCanvas("c10","c10",250,250,1000,500);
	//c10->Divide(2,1);

	
	// ----------------------------------
	
	double win1 = 0.02;
	double fjwin = 8;

	TLatex lat1;
	lat1.SetTextFont(42);
	lat1.SetTextSize(0.05);
	
	TLatex lat3;
	lat3.SetTextFont(42);
	lat3.SetTextSize(0.05);
	lat3.SetTextColor(4);
	
	TLatex lat2;
	lat2.SetTextFont(42);
	lat2.SetTextColor(2);
	lat2.SetTextSize(0.05);

	// *****************
	// electron channel
	// *****************
	
	TFile *fe = new TFile("MJ2e_S98k_NR100k_HD9580M_ntp2.root");
	TTree *te = (TTree*)fe->Get("ntp2");
	
	te->SetAlias("msum","(f4cxd0m+f4cxd1m)");
	
	double ecntr = 3.0946, ewin = 0.025; 	
	double ejmin = ecntr-fjwin*ewin, ejmax=ecntr+fjwin*ewin;
	//double ejmin = 2.9, ejmax=3.3;
	
	TString mcute = "abs(f4cxd0m-3.0946)<0.2";
	TString cute  = "msum>3.77&&xd0d0pide>0.95&&xd0d1pide>0.95&&f4cxm>3.867&&f4cxm<3.874&&xd0oang<2.1&&xpcm<0.4";
			
	TH1F *he1pres = new TH1F("he1pres","#bar{p}p 4C fit (J/#psi#rightarrow e^{+}e^{-}), preselection",100,3.872-win1,3.872+win1);
	TH1F *he1s    = new TH1F("he1s","#bar{p}p 4C fit (J/#psi#rightarrow e^{+}e^{-}), full selection",100,3.872-win1,3.872+win1);

	TH1F *he2pres = new TH1F("he2pres","m(J/#psi#rightarrow e^{+}e^{-}) 4C fit, preselection",100,ejmin,ejmax);	
	TH1F *he2s    = new TH1F("he2s","m(J/#psi#rightarrow e^{+}e^{-}) 4C fit, full selection",100,ejmin,ejmax);	
	
	TH1F *he1preb = new TH1F("he1preb","#bar{p}p 4C fit (J/#psi#rightarrow e^{+}e^{-}), preselection (Bg)",100,3.872-win1,3.872+win1);
	TH1F *he1b    = new TH1F("he1b","#bar{p}p 4C fit (J/#psi#rightarrow e^{+}e^{-}), full selection (Bg)",100,3.872-win1,3.872+win1);
	
	TH1F *he2preb = new TH1F("he2preb","m(J/#psi#rightarrow e^{+}e^{-}) 4C fit, preselection (Bg)",100,ejmin,ejmax);	
	TH1F *he2b    = new TH1F("he2b","m(J/#psi#rightarrow e^{+}e^{-}) 4C fit, full selection (Bg)",100,ejmin,ejmax);	
	
	TH1F *he2preb2 = new TH1F("he2preb2","m(J/#psi#rightarrow e^{+}e^{-}) 4C fit, preselection (Bg)",100,ejmin,ejmax);	
	TH1F *he2b2    = new TH1F("he2b2","m(J/#psi#rightarrow e^{+}e^{-}) 4C fit, full selection (Bg)",100,ejmin,ejmax);	
	
	confhist(he1pres,1,"m_{fit}(e^{+}e^{-}#pi^{+}#pi^{-}) [GeV/c^{2}]");
	confhist(he1s,1,"m_{fit}(e^{+}e^{-}#pi^{+}#pi^{-}) [GeV/c^{2}]");
	confhist(he2pres,1,"m_{fit}(e^{+}e^{-}) [GeV/c^{2}]");
	confhist(he2s,1,"m_{fit}(e^{+}e^{-}) [GeV/c^{2}]");
	confhist(he1preb,2);
	confhist(he1b,2);
	confhist(he2preb,2);
	confhist(he2b,2);
	confhist(he2preb2,4);
	confhist(he2b2,4);
	
	//he1preb->SetLineStyle(2);
	//he2preb->SetLineStyle(2);
	//he2preb2->SetLineStyle(7);
	//he1b->SetLineStyle(2);
	//he2b->SetLineStyle(2);
	//he2b2->SetLineStyle(7);
	
	// create the histograms
	te->Project("he1pres","f4cxm",sigcut);
	te->Project("he1s","f4cxm",mcute+"&&"+cute+"&&"+sigcut);

	te->Project("he2pres","f4cxd0m", sigcut);	
	te->Project("he2s","f4cxd0m",cute+"&&"+sigcut);

	te->Project("he1preb","f4cxm", bkgcut1);
	te->Project("he1b","f4cxm",mcute+"&&"+cute+"&&"+bkgcut1);

	te->Project("he2preb","f4cxd0m", bkgcut1);	
	te->Project("he2b","f4cxd0m",cute+"&&"+bkgcut1);
	
	te->Project("he2preb2","f4cxd0m", bkgcut2);	
	te->Project("he2b2","f4cxd0m",cute+"&&"+bkgcut2);

	// determine signal and background counts
	double Sep  = cntEvents(te,mcute+"&&"+sigcut);//te->GetEntries(mcute+"&&"+sigcut);
	double Bep  = cntEvents(te,mcute+"&&"+bkgcut1);//te->GetEntries(mcute+"&&"+bkgcut);
	
	double Se  = cntEvents(te,mcute+"&&"+cute+"&&"+sigcut);//e->GetEntries(mcute+"&&"+cute+"&&"+sigcut);
	double Be  = he2b->Integral(1,100);
	double Be2 = cntEvents(te,mcute+"&&"+cute+"&&"+bkgcut2);
	double dBe = sqrt(he2b->Integral(1,100))/fjwin;
	
	// display all variables with cuts
	//drawVars(te,"msum>3.77&&xd0d0pide>0.95&&f4cxm>3.867&&f4cxm<3.874&&xd0oang<2.1&&xpcm<0.4", c3, "f4cxm>3.83&&f4cxm<3.92&&msum>3&&msum<4&&xpcm<3");
	//TString savepath   = "fig/";

	//saveCanvas(c3,"J2e_cutvars");

	//c3->SaveAs("fig/J2e_cutvars_2.gif");
	//c3->SaveAs("fig/J2e_cutvars_2.C");

	// QA check for bkg distortion cut wise
	//varsQa(te, "f4cxd0m", "msum>3.77&&xd0d0pide>0.95&&f4cxm>3.867&&f4cxm<3.874&&xd0oang<2.1&&xpcm<0.4", c5, mcute+"&&"+bkgcut1,10);
	
	//saveCanvas(c5,"J2e_qaflat");
	//c5->SaveAs("fig/J2e_qaflat_2.gif");
	//c5->SaveAs("fig/J2e_qaflat_2.C");

	// plot histograms, sig+bg before and after cuts
	c1->cd(2); gPad->SetLogy();
	he2pres->Draw();
	
	f1->SetParameters(1000,ecntr,0.0084);
	he2pres->Fit("f1","","",ecntr-ewin/2., ecntr+ewin/2.);
	double sige = f1->GetParameter(2);
	
	//lat3.DrawLatex((ejmax-ejmin)*0.02+ejmin,he2pres->GetMaximum()*0.95,TString::Format("#sigma_{core} = %.1f MeV/c^{2}",sige*1000));
	//lat1.DrawLatex((ejmax-ejmin)*0.02+ejmin,he2s->GetMaximum()*0.5,TString::Format("S = %.0f",Sep));
	//lat2.DrawLatex((ejmax-ejmin)*0.02+ejmin,he2s->GetMaximum()*0.26,TString::Format("B = %.0f",Bep));
	
	he2preb->Draw("same");		
	//he2preb2->Draw("same");		
	//drawwin(he2pres,ecntr, ewin);
	PandaSmartLabel("R");


	c1->cd(3); gPad->SetLogy();
	he2s->Draw(); he2b->Draw("same");he2b2->Draw("same");			
	//drawwin(he2s,ecntr, ewin);
	
	//lat1.DrawLatex((ejmax-ejmin)*0.02+ejmin,he2s->GetMaximum()*0.95,TString::Format("S = %.0f",Se));
	//lat3.DrawLatex((ejmax-ejmin)*0.02+ejmin,he2s->GetMaximum()*0.5,TString::Format("B_{nr} = %.0f",Be2));
	//lat2.DrawLatex((ejmax-ejmin)*0.02+ejmin,he2s->GetMaximum()*0.25,TString::Format("B_{hd} = %.2f",Be));
	
	TLine lsc; lsc.SetLineColor(2); lsc.SetLineWidth(2); lsc.SetLineStyle(9);
	double hdevtscale = BkgXs*Be*S0e/(SigXs*B0e*BRX*BRe)/100;	
	lsc.DrawLine(ejmin,hdevtscale,ejmax,hdevtscale);
	cout <<"ee bkg level scaled = "<<hdevtscale<<endl;
	
	//c1->cd(3);
	//he1preb->Draw(); he1b->Draw("same");		
	
	//c1->cd(4);
	//he2preb->Draw(); he2b->Draw("same");		
	PandaSmartLabel("R");

		
	// 4C fitted pbarp mass
	c1->cd(1); gPad->SetLogy();
	he1pres->SetNdivisions(505);
	he1pres->Draw();
	he1preb->Draw("same");

	PandaSmartLabel("R");

	saveCanvas(c1,"J2e_selection");
	//c1->SaveAs("fig/J2e_selection_2.gif");
	//c1->SaveAs("fig/J2e_selection_2.C");
		

	// *****************
	// muon channel
	// *****************
	
	TFile *fm = new TFile("MJ2mu_S100k_NR99k_HD8870M_ntp2.root");
	TTree *tm = (TTree*)fm->Get("ntp2");
		
	tm->SetAlias("msum","(f4cxd0m+f4cxd1m)");
	
	
	double mcntr = 3.097, mwin = 0.02; 
	double mjmin = mcntr-fjwin*mwin, mjmax=mcntr+fjwin*mwin;
	//double mjmin = 2.88, mjmax=3.3;
	 
	TString mcutm = "abs(f4cxd0m-3.097)<0.16";
	TString cutm("msum>3.78&&xd0d0pidmu>0.99&&xd0d1pidmu>0.99&&xd0oang<1.4&&essph<0.11");
	
	TH1F *hm1pres = new TH1F("hm1pres","#bar{p}p 4C fit (J/#psi#rightarrow #mu^{+}#mu^{-}), preselection",100,3.872-win1,3.872+win1);
	TH1F *hm1s    = new TH1F("hm1s","#bar{p}p 4C fit (J/#psi#rightarrow #mu^{+}#mu^{-}), full selection",100,3.872-win1,3.872+win1);

	TH1F *hm2pres = new TH1F("hm2pres","m(J/#psi#rightarrow #mu^{+}#mu^{-}) 4C fit, preselection",100,mjmin,mjmax);	
	TH1F *hm2s    = new TH1F("hm2s","m(J/#psi#rightarrow #mu^{+}#mu^{-}) 4C fit, full selection",100,mjmin,mjmax);	
	
	TH1F *hm1preb = new TH1F("hm1preb","#bar{p}p 4C fit (J/#psi#rightarrow #mu^{+}#mu^{-}), preselection (Bg)",100,3.872-win1,3.872+win1);
	TH1F *hm1b    = new TH1F("hm1b","#bar{p}p 4C fit (J/#psi#rightarrow #mu^{+}#mu^{-}), full selection (Bg)",100,3.872-win1,3.872+win1);
	
	TH1F *hm2preb = new TH1F("hm2preb","m(J/#psi#rightarrow #mu^{+}#mu^{-}) 4C fit, preselection (Bg)",100,mjmin,mjmax);
	TH1F *hm2b    = new TH1F("hm2b","m(J/#psi#rightarrow #mu^{+}#mu^{-}) 4C fit, full selection (Bg)",100,mjmin,mjmax);

	TH1F *hm2preb2 = new TH1F("hm2preb2","m(J/#psi#rightarrow #mu^{+}#mu^{-}) 4C fit, preselection (Bg)",100,mjmin,mjmax);
	TH1F *hm2b2    = new TH1F("hm2b2","m(J/#psi#rightarrow #mu^{+}#mu^{-}) 4C fit, full selection (Bg)",100,mjmin,mjmax);

	confhist(hm1pres,1,"m_{fit}(#mu^{+}#mu^{-}#pi^{+}#pi^{-}) [GeV/c^{2}]");
	confhist(hm1s,1,"m_{fit}(#mu^{+}mu^{-}#pi^{+}#pi^{-}) [GeV/c^{2}]");
	confhist(hm2pres,1,"m_{fit}(#mu^{+}#mu^{-}) [GeV/c^{2}]");
	confhist(hm2s,1,"m_{fit}(#mu^{+}#mu^{-}) [GeV/c^{2}]");
	confhist(hm1preb,2);
	confhist(hm1b,2);
	confhist(hm2preb,2);
	confhist(hm2b,2);
	confhist(hm2preb2,4);
	confhist(hm2b2,4);


	// create the histograms
	tm->Project("hm1pres","f4cxm",sigcut);
	tm->Project("hm1s","f4cxm",mcutm+"&&"+cutm+"&&"+sigcut);

	tm->Project("hm2pres","f4cxd0m", sigcut);	
	tm->Project("hm2s","f4cxd0m",cutm+"&&"+sigcut);

	tm->Project("hm1preb","f4cxm", bkgcut1);
	tm->Project("hm1b","f4cxm",mcutm+"&&"+cutm+"&&"+bkgcut1);

	tm->Project("hm2preb","f4cxd0m", bkgcut1);	
	tm->Project("hm2b","f4cxd0m",cutm+"&&"+bkgcut1);
		
	tm->Project("hm2preb2","f4cxd0m", bkgcut2);	
	tm->Project("hm2b2","f4cxd0m",cutm+"&&"+bkgcut2);
	
	// determine signal and background counts
	double Smp  = cntEvents(tm,mcutm+"&&"+sigcut);//tm->GetEntries(mcutm+"&&"+sigcut);
	double Bmp  = cntEvents(tm,mcutm+"&&"+bkgcut1);//tm->GetEntries(mcutm+"&&"+bkgcut);
	
	double Sm  = cntEvents(tm,mcutm+"&&"+cutm+"&&"+sigcut);//tm->GetEntries(mcutm+"&&"+cutm+"&&"+sigcut);
	double Bm  = hm2b->Integral(1,100);
	double Bm2 = cntEvents(tm,mcutm+"&&"+cutm+"&&"+bkgcut2);
	double dBm = sqrt(hm2b->Integral(1,100))/fjwin;

	// display all variables with cuts
	//drawVars(tm,"msum>3.78&&xd0d0pidmu>0.99&&xd0oang<1.4&&essph<0.11", c4, "msum>3&&msum<4");
	//saveCanvas(c4,"J2mu_cutvars");
	//c4->SaveAs("fig/J2mu_cutvars_2.gif");
	//c4->SaveAs("fig/J2mu_cutvars_2.C");
	
	// QA check for bkg distortion cut wise
	//varsQa(tm, "f4cxd0m", "msum>3.78&&xd0d0pidmu>0.99&&xd0oang<1.4&&essph<0.11", c6, mcutm+"&&"+bkgcut1);
	//saveCanvas(c6,"J2mu_qaflat");
	//c6->SaveAs("fig/J2mu_qaflat_2.gif");
	//c6->SaveAs("fig/J2mu_qaflat_2.C");

	// plot histograms, sig+bg before and after cuts
	c2->cd(2); gPad->SetLogy();
	hm2pres->Draw(); 
	f1->SetParameters(1000,ecntr,0.007);
	hm2pres->Fit("f1","","",mcntr-mwin/2., mcntr+mwin/2.);
	double sigm = f1->GetParameter(2);
	
	//lat3.DrawLatex((mjmax-mjmin)*0.02+mjmin,hm2pres->GetMaximum()*0.95,TString::Format("#sigma_{core} = %.1f MeV/c^{2}",sigm*1000));
	
	//lat1.DrawLatex((mjmax-mjmin)*0.02+mjmin,hm2s->GetMaximum()*0.5,TString::Format("S = %.0f",Smp));
	//lat2.DrawLatex((mjmax-mjmin)*0.02+mjmin,hm2s->GetMaximum()*0.26,TString::Format("B = %.0f",Bmp));
	
	hm2preb->Draw("same");		
	//hm2preb2->Draw("same");		
	//drawwin(hm2pres,mcntr, mwin);
	PandaSmartLabel("R");

	c2->cd(3); gPad->SetLogy();
	hm2s->Draw(); hm2b->Draw("same");hm2b2->Draw("same");
	//drawwin(hm2s,mcntr, mwin);
			
	//lat1.DrawLatex((mjmax-mjmin)*0.02+mjmin,hm2s->GetMaximum()*0.95,TString::Format("S = %.0f",Sm));
	//lat3.DrawLatex((mjmax-mjmin)*0.02+mjmin,hm2s->GetMaximum()*0.5,TString::Format("B_{nr} = %.0f",Bm2));
	//lat2.DrawLatex((mjmax-mjmin)*0.02+mjmin,hm2s->GetMaximum()*0.25,TString::Format("B_{hd} = %.1f",Bm));

	double hdevtscalem = BkgXs*Bm*S0m/(SigXs*B0m*BRX*BRm)/100;	
	lsc.DrawLine(mjmin,hdevtscalem,mjmax,hdevtscalem);
	cout <<"mumu bkg level scaled = "<<hdevtscalem<<endl;

	// ***********************
	// 4C fitted pbarp mass
	// ***********************
	PandaSmartLabel("R");

	c2->cd(1); gPad->SetLogy();
	hm1pres->SetNdivisions(505);
	hm1pres->Draw();
	hm1preb->Draw("same");
		
	//saveCanvas(c8,"pbarp_4c_preselection");
	//c8->SaveAs("fig/pbarp_4c_preselection_2.gif");
	//c8->SaveAs("fig/pbarp_4c_preselection_2.C");
	PandaSmartLabel("R");
	

	saveCanvas(c2,"J2mu_selection");
	//c2->SaveAs("fig/J2mu_selection_2.gif");
	//c2->SaveAs("fig/J2mu_selection_2.C");
	
	
	
	
	
	
	// ***********************
	// Correlation J and m_J + m_rho for mu mu channel	
	// ***********************

	tm->SetEventList(0);
	
	int ncbin=40;
	
	TLine lll;
	lll.SetLineColor(2);
	lll.SetLineWidth(2);
	lll.SetLineStyle(7);
	
	//TH2F* hcorrs=new TH2F("hcorrs","Signal: m(#mu^{+}#mu^{-}) vs. m(#pi^{+}#pi^{-})",ncbin,0.27,1.,ncbin,2.8,3.3);
	//TH2F* hcorrb=new TH2F("hcorrb","Haronic bkg: m(#mu^{+}#mu^{-}) vs. m(#pi^{+}#pi^{-})",ncbin,0.27,1.,ncbin,2.8,3.3);
	//TH2F* hcorrb2=new TH2F("hcorrb2","J/#psi NR bkg: m(#mu^{+}#mu^{-}) vs. m(#pi^{+}#pi^{-})",ncbin,0.27,1.,ncbin,2.8,3.3);

	TH2F* hcorrs=new TH2F("hcorrs","Signal: m(#mu#mu) vs. m(#mu#mu)+m(#pi#pi)",          ncbin, 3.5, 3.9, ncbin, mcntr-8*mwin, mcntr+8*mwin);
	TH2F* hcorrb=new TH2F("hcorrb","Haronic bkg: m(#mu#mu) vs. m(#mu#mu)+m(#pi#pi)",     ncbin, 3.5, 3.9, ncbin, mcntr-8*mwin, mcntr+8*mwin);
	TH2F* hcorrb2=new TH2F("hcorrb2","J/#psi NR bkg: m(#mu#mu) vs. m(#mu#mu)+m(#pi#pi)", ncbin, 3.5, 3.9, ncbin, mcntr-8*mwin, mcntr+8*mwin);
	
	confhist(hcorrs,1,"m(#pi^{+}#pi^{-})+m(#mu^{+}#mu^{-}) [GeV/c^{2}]","m(#mu^{+}#mu^{-}) [GeV/c^{2}]", 1.75, 505);
	confhist(hcorrb,1,"m(#pi^{+}#pi^{-})+m(#mu^{+}#mu^{-}) [GeV/c^{2}]","m(#mu^{+}#mu^{-}) [GeV/c^{2}]", 1.75, 505);
	confhist(hcorrb2,1,"m(#pi^{+}#pi^{-})+m(#mu^{+}#mu^{-}) [GeV/c^{2}]","m(#mu^{+}#mu^{-}) [GeV/c^{2}]",1.75, 505);
	
	hcorrs->SetFillColor(1);hcorrs->SetFillStyle(1000);
	hcorrb->SetFillColor(1);hcorrs->SetFillStyle(1000);
	hcorrb2->SetFillColor(1);hcorrs->SetFillStyle(1000);
		
	c9->cd(1);
	tm->Project("hcorrs","f4cxd0m:msum",sigcut);
	hcorrs->Draw("box");
	PandaSmartLabel("R");
	//lll.DrawLine(3.5,mcntr-8*mwin,3.5,mcntr+8*mwin);
	
	c9->cd(2);
	tm->Project("hcorrb","f4cxd0m:msum",bkgcut1);
	hcorrb->Draw("box");
	PandaSmartLabel("R");
	//lll.DrawLine(3.5,mcntr-8*mwin,3.5,mcntr+8*mwin);
	
	c9->cd(3);
	tm->Project("hcorrb2","f4cxd0m:msum",bkgcut2);
	hcorrb2->Draw("box");
	PandaSmartLabel("R");
	//lll.DrawLine(3.5,mcntr-8*mwin,3.5,mcntr+8*mwin);
	
	c9->Update();

	saveCanvas(c9,"J2mu_Jrho_corr");
	
	
	
	// ***********************
	// Correlation J and rho for mu mu channel	
	// ***********************
	double mrmin = 0.25, mrmax = 1.0;
	
	TH2F* hcorrs_2=new TH2F("hcorrs_2","Signal: m(#mu#mu) vs. m(#pi#pi)",          ncbin, mrmin, mrmax, ncbin, mcntr-8*mwin, mcntr+8*mwin);
	TH2F* hcorrb_2=new TH2F("hcorrb_2","Haronic bkg: m(#mu#mu) vs. m(#pi#pi)",     ncbin, mrmin, mrmax, ncbin, mcntr-8*mwin, mcntr+8*mwin);
	TH2F* hcorrb2_2=new TH2F("hcorrb2_2","J/#psi NR bkg: m(#mu#mu) vs. m(#pi#pi)", ncbin, mrmin, mrmax, ncbin, mcntr-8*mwin, mcntr+8*mwin);
	
	confhist(hcorrs_2,1,"m(#pi^{+}#pi^{-}) [GeV/c^{2}]","m(#mu^{+}#mu^{-}) [GeV/c^{2}]", 1.75, 505);
	confhist(hcorrb_2,1,"m(#pi^{+}#pi^{-}) [GeV/c^{2}]","m(#mu^{+}#mu^{-}) [GeV/c^{2}]", 1.75, 505);
	confhist(hcorrb2_2,1,"m(#pi^{+}#pi^{-}) [GeV/c^{2}]","m(#mu^{+}#mu^{-}) [GeV/c^{2}]",1.75, 505);
	
	hcorrs_2->SetFillColor(1);hcorrs->SetFillStyle(1000);
	hcorrb_2->SetFillColor(1);hcorrs->SetFillStyle(1000);
	hcorrb2_2->SetFillColor(1);hcorrs->SetFillStyle(1000);
		
	c9a->cd(1);
	tm->Project("hcorrs_2","f4cxd0m:f4cxd1m",sigcut);
	hcorrs_2->Draw("box");
	PandaSmartLabel("R");
	
	c9a->cd(2);
	tm->Project("hcorrb_2","f4cxd0m:f4cxd1m",bkgcut1);
	hcorrb_2->Draw("box");
	PandaSmartLabel("R");
	
	c9a->cd(3);
	tm->Project("hcorrb2_2","f4cxd0m:f4cxd1m",bkgcut2);
	hcorrb2_2->Draw("box");
	PandaSmartLabel("R");
	
	c9a->Update();

	saveCanvas(c9a,"J2mu_Jrho_corr2");
	
	
	
	
	
	// ***********************
	// Correlation J and m_J + m_rho for e e channel	
	// ***********************

	te->SetEventList(0);
	
	//int ncbin=40;
	
	//TLine lll;
	//lll.SetLineColor(2);
	//lll.SetLineWidth(2);
	//lll.SetLineStyle(7);
	
	//TH2F* hcorrs=new TH2F("hcorrs","Signal: m(#mu^{+}#mu^{-}) vs. m(#pi^{+}#pi^{-})",ncbin,0.27,1.,ncbin,2.8,3.3);
	//TH2F* hcorrb=new TH2F("hcorrb","Haronic bkg: m(#mu^{+}#mu^{-}) vs. m(#pi^{+}#pi^{-})",ncbin,0.27,1.,ncbin,2.8,3.3);
	//TH2F* hcorrb2=new TH2F("hcorrb2","J/#psi NR bkg: m(#mu^{+}#mu^{-}) vs. m(#pi^{+}#pi^{-})",ncbin,0.27,1.,ncbin,2.8,3.3);

	TH2F* hcorres=new TH2F("hcorres","Signal: m(ee) vs. m(ee)+m(#pi#pi)",          ncbin, 3.5, 3.9, ncbin, mcntr-8*mwin, mcntr+8*mwin);
	TH2F* hcorreb=new TH2F("hcorreb","Haronic bkg: m(ee) vs. m(ee)+m(#pi#pi)",     ncbin, 3.5, 3.9, ncbin, mcntr-8*mwin, mcntr+8*mwin);
	TH2F* hcorreb2=new TH2F("hcorreb2","J/#psi NR bkg: m(ee) vs. m(ee)+m(#pi#pi)", ncbin, 3.5, 3.9, ncbin, mcntr-8*mwin, mcntr+8*mwin);
	
	confhist(hcorres,1,"m(#pi^{+}#pi^{-})+m(e^{+}e^{-}) [GeV/c^{2}]","m(e^{+}e^{-}) [GeV/c^{2}]", 1.75, 505);
	confhist(hcorreb,1,"m(#pi^{+}#pi^{-})+m(e^{+}e^{-}) [GeV/c^{2}]","m(e^{+}e^{-}) [GeV/c^{2}]", 1.75, 505);
	confhist(hcorreb2,1,"m(#pi^{+}#pi^{-})+m(e^{+}e^{-}) [GeV/c^{2}]","m(e^{+}e^{-}) [GeV/c^{2}]",1.75, 505);
	
	hcorres->SetFillColor(1);hcorres->SetFillStyle(1000);
	hcorreb->SetFillColor(1);hcorreb->SetFillStyle(1000);
	hcorreb2->SetFillColor(1);hcorreb2->SetFillStyle(1000);
		
	c9b->cd(1);
	te->Project("hcorres","f4cxd0m:msum",sigcut);
	hcorres->Draw("box");
	PandaSmartLabel("R");
	//lll.DrawLine(3.5,mcntr-8*mwin,3.5,mcntr+8*mwin);
	
	c9b->cd(2);
	te->Project("hcorreb","f4cxd0m:msum",bkgcut1);
	hcorreb->Draw("box");
	PandaSmartLabel("R");
	//lll.DrawLine(3.5,mcntr-8*mwin,3.5,mcntr+8*mwin);
	
	c9b->cd(3);
	te->Project("hcorreb2","f4cxd0m:msum",bkgcut2);
	hcorreb2->Draw("box");
	PandaSmartLabel("R");
	//lll.DrawLine(3.5,mcntr-8*mwin,3.5,mcntr+8*mwin);
	
	c9b->Update();

	saveCanvas(c9b,"J2e_Jrho_corr");
	
	// ***********************
	// Correlation J and rho for e e channel	
	// ***********************
	//double mrmin = 0.25, mrmax = 1.0;
	
	TH2F* hcorres_2=new TH2F("hcorres_2","Signal: m(ee) vs. m(#pi#pi)",          ncbin, mrmin, mrmax, ncbin, mcntr-8*mwin, mcntr+8*mwin);
	TH2F* hcorreb_2=new TH2F("hcorreb_2","Haronic bkg: m(ee) vs. m(#pi#pi)",     ncbin, mrmin, mrmax, ncbin, mcntr-8*mwin, mcntr+8*mwin);
	TH2F* hcorreb2_2=new TH2F("hcorreb2_2","J/#psi NR bkg: m(ee) vs. m(#pi#pi)", ncbin, mrmin, mrmax, ncbin, mcntr-8*mwin, mcntr+8*mwin);
	
	confhist(hcorres_2,1,"m(#pi^{+}#pi^{-}) [GeV/c^{2}]","m(e^{+}e^{-}) [GeV/c^{2}]", 1.75, 505);
	confhist(hcorreb_2,1,"m(#pi^{+}#pi^{-}) [GeV/c^{2}]","m(e^{+}e^{-}) [GeV/c^{2}]", 1.75, 505);
	confhist(hcorreb2_2,1,"m(#pi^{+}#pi^{-}) [GeV/c^{2}]","m(e^{+}e^{-}) [GeV/c^{2}]",1.75, 505);
	
	hcorres_2->SetFillColor(1);hcorres_2->SetFillStyle(1000);
	hcorreb_2->SetFillColor(1);hcorreb_2->SetFillStyle(1000);
	hcorreb2_2->SetFillColor(1);hcorreb2_2->SetFillStyle(1000);
		
	c9c->cd(1);
	te->Project("hcorres_2","f4cxd0m:f4cxd1m",sigcut);
	hcorres_2->Draw("box");
	PandaSmartLabel("R");
	
	c9c->cd(2);
	te->Project("hcorreb_2","f4cxd0m:f4cxd1m",bkgcut1);
	hcorreb_2->Draw("box");
	PandaSmartLabel("R");
	
	c9c->cd(3);
	te->Project("hcorreb2_2","f4cxd0m:f4cxd1m",bkgcut2);
	hcorreb2_2->Draw("box");
	PandaSmartLabel("R");
	
	c9c->Update();

	saveCanvas(c9c,"J2e_Jrho_corr2");
	
	//c9->SaveAs("fig/J2mu_Jrho_corr_2.gif");
	//c9->SaveAs("fig/J2mu_Jrho_corr_2.C");
	
	//c10->cd(1);
	//tm->Project("hcorrs","f4cxd0m:msum",sigcut);
	//hcorrs->Draw("col");
	
	//c10->cd(2);
	//tm->Project("hcorrb2","f4cxd0m:msum",bkgcut2);
	//hcorrb2->Draw("col");
	
	//c10->Update();

	//saveCanvas(c10,"J2mu_Jrho_corr_NR");
	//c10->SaveAs("fig/J2mu_Jrho_corr_NR_2.gif");
	//c10->SaveAs("fig/J2mu_Jrho_corr_NR_2.C");
	
	
	
	// ***********************
	// Prefilter check
	// ***********************
	
	
	c7->cd();
	gStyle->SetLineScalePS(3);

	TFile *fmaf = new TFile("MDPMAltFiltJ2mu_Bg2600M_ntp2.root");
	TTree *tmaf = (TTree*)fmaf->Get("ntp2");

	TH1F *hFilt28 = new TH1F("hFilt28","m(J/#psi#rightarrow #mu^{+}#mu^{-}) 4C fit, preselection (Bg)",30,mjmin,mjmax);	
	TH1F *hFilt25 = new TH1F("hFilt25","m(J/#psi#rightarrow #mu^{+}#mu^{-}) 4C fit, preselection (Bg)",30,mjmin,mjmax);	
	confhist(hFilt28,2,"m_{fit}(#mu^{+}#mu^{-}) [GeV/c^{2}]");
	confhist(hFilt25,4,"m_{fit}(#mu^{+}#mu^{-}) [GeV/c^{2}]");
	
	hFilt28->SetMarkerColor(2);hFilt28->SetMarkerStyle(20);	
	hFilt25->SetMarkerColor(4);hFilt28->SetMarkerStyle(21);

	tm->SetEventList(0);
	tm->Project("hFilt28","f4cxd0m", bkgcut1);
	hFilt28->Sumw2();
	hFilt28->Scale(B0maf/B0m);

	tmaf->Project("hFilt25","f4cxd0m", bkgcut1);
	hFilt25->SetMinimum(0);

	TF1 *f2=new TF1("f2","pol2(0)");
	f2->SetParameters(1,1,1);
	f2->SetLineColor(4);
	f2->SetLineWidth(2);
	hFilt25->Draw("e");
	hFilt25->Fit(f2);

	double integ25  = f2->Integral(mjmin,mjmax);	
	double integ25r = f2->Integral(mcntr-mwin,mcntr+mwin);	
	
	f2->SetLineColor(2);
	//f2->SetLineStyle(7);
	hFilt28->Draw("same e");
	hFilt28->Fit(f2);
	//((TF1*)hFilt28->GetListOfFunctions()->At(0))->SetLineColor(4);
	
	double integ28  = f2->Integral(mjmin,mjmax);	
	double integ28r = f2->Integral(mcntr-mwin,mcntr+mwin);	
	
	cout <<"Relative integral difference        : "<<fabs(integ25-integ28)/integ28<<endl;
	cout <<"Relative integral difference in ROI : "<<fabs(integ25r-integ28r)/integ28r<<endl;
	
	//drawwin(hFilt25, mcntr, mwin);
	
	TLegend *leg1 = new TLegend(0.45,0.83,0.96,0.95);
	leg1->SetBorderSize(1);
	leg1->SetTextFont(42);
	leg1->AddEntry(hFilt25,"m_{filter}(#mu#mu) > 2.5 GeV/c^{2}","lp");
	leg1->AddEntry(hFilt28,"m_{filter}(#mu#mu) > 2.8 GeV/c^{2}","lp");
	leg1->Draw();
	
	//double hgt = writeLatex(0.18,0.25,Form("#Delta_{tot} = %.1f %%",fabs(integ25-integ28)/integ28*100),1);
	//writeLatex(0.18,0.25-hgt*1.2,Form("#Delta_{roi} = %.1f %%",fabs(integ25r-integ28r)/integ28r*100),6);
	
	cout <<"Filter:"<<fabs(integ25-integ28)/integ28*100<<endl;
	PandaSmartLabel("L",0.2,0.2);

	saveCanvas(c7,"J2mu_filtercheck");
	//c7->SaveAs("fig/J2mu_filtercheck_2.gif");
	//c7->SaveAs("fig/J2mu_filtercheck_2.C");
	
	// ***********************
	// Compute numbers	
	// ***********************
	
	cout <<he2b->Integral(1,100)<<"  "<<hm2b->Integral(1,100)<<endl;

	
	double effSe    = Se/S0e;
	double effBe    = Be/B0e;
	double effBe2   = Be2/B0e2;
	double errEffBe = dBe/B0e; 
	
	double effSm    = Sm/S0m;
	double effBm    = Bm/B0m;
	double effBm2   = Bm2/B0m2;
	double errEffBm = dBm/B0m; 

	double SNe = (effSe*BRe*BRX*SigXs)/(effBe*BkgXs + effBe2*BkgNR*BRe);
	double SNm = (effSm*BRm*BRX*SigXs)/(effBm*BkgXs + effBm2*BkgNR*BRm);
	double SN  = ((effSe*BRe + effSm*BRm)*BRX*SigXs)/((effBe+effBm)*BkgXs + (effBe2*BRe + effBm2*BRm)*BkgNR);
	
	cout <<"sigma_core_e = "<<sige<<"  sigma_core_mu = "<<sigm<<endl;
	printf("effSe = %.3f  effBe = %.3e  effBeNR = %.3f\neffSm = %.3f  effBm = %.3e  effBeNR = %.3f\n",effSe, effBe, effBe2, effSm, effBm, effBm2);
	printf("S_e = %.0f  B_e = %.2f  B_e2 = %.0f | S_m = %.0f  B_m = %.2f  B_m2 = %.0f\n",Se, Be, Be2, Sm, Bm, Bm2);
	printf("SNe = %.1f:1  SNmu = %.1f:1  SNtot = %.1f:1\n",SNe, SNm, SN);  
	
}