{
gStyle->SetStatBorderSize(1);
gStyle->SetStatColor(0);
gStyle->SetTitleFillColor(0);
gStyle->SetCanvasColor(0);
gStyle->SetDrawBorder(0);
gStyle->SetPadColor(0);
gStyle->SetPadBorderMode(0);
gStyle->SetFrameBorderMode(0);
gStyle->SetPalette(1,0);
gStyle->SetOptStat(0);
gStyle->SetOptTitle(0);

TFile* ff = new TFile("tree_cu63.root","r");
TTree* T_bl = (TTree*)gROOT->FindObject("T_bl");
TTree* T_dl = (TTree*)gROOT->FindObject("T_dl");
TTree* T_bs = (TTree*)gROOT->FindObject("T_bs");
TTree* T_ds = (TTree*)gROOT->FindObject("T_ds");

TH1F* hbl1 = new TH1F("hbl1","Total momentum of pKpi",500,0,2);
TH1F* hbl2 = new TH1F("hbl2","Total momentum of pKpi",500,0,2);
TH1F* hbs1 = new TH1F("hbs1","Total momentum of pKpi",500,0,2);
TH1F* hbs2 = new TH1F("hbs2","Total momentum of pKpi",500,0,2);
TH1F* hdl = new TH1F("hdl","Total momentum of pKpi",500,0,2);
TH1F* hds = new TH1F("hds","Total momentum of pKpi",500,0,2);

c1 = new TCanvas("c1","",800,600);

T_bl->Draw("sqrt((p1.af.P*cos(p1.af.Ph)*sin(p1.af.Th)+pi.af.P*cos(pi.af.Ph)*sin(pi.af.Th)+k.af.P*cos(k.af.Ph)*sin(k.af.Th))**2+(p1.af.P*sin(p1.af.Th)*sin(p1.af.Ph)+pi.af.P*sin(pi.af.Th)*sin(pi.af.Ph)+k.af.P*sin(k.af.Th)*sin(k.af.Ph))**2+(p1.af.P*cos(p1.af.Th)+pi.af.P*cos(pi.af.Th)+k.af.P*cos(k.af.Th))**2)>>hbl1","w.pf*w.lk*w.ppi*w.nf/50.");
 //*(((p1.af.P*cos(p1.af.Ph)*sin(p1.af.Th)+pi.af.P*cos(pi.af.Ph)*sin(pi.af.Th)+k.af.P*cos(k.af.Ph)*sin(k.af.Th))**2+(p1.af.P*sin(p1.af.Th)*sin(p1.af.Ph)+pi.af.P*sin(pi.af.Th)*sin(pi.af.Ph)+k.af.P*sin(k.af.Th)*sin(k.af.Ph))**2+(p1.af.P*cos(p1.af.Th)+pi.af.P*cos(pi.af.Th)+k.af.P*cos(k.af.Th))**2)<2)");

T_bl->Draw("sqrt((p2.af.P*cos(p2.af.Ph)*sin(p2.af.Th)+pi.af.P*cos(pi.af.Ph)*sin(pi.af.Th)+k.af.P*cos(k.af.Ph)*sin(k.af.Th))**2+(p2.af.P*sin(p2.af.Th)*sin(p2.af.Ph)+pi.af.P*sin(pi.af.Th)*sin(pi.af.Ph)+k.af.P*sin(k.af.Th)*sin(k.af.Ph))**2+(p2.af.P*cos(p2.af.Th)+pi.af.P*cos(pi.af.Th)+k.af.P*cos(k.af.Th))**2)>>hbl2","w.pf*w.lk*w.ppi*w.nf/50.");
 //*(((p2.af.P*cos(p2.af.Ph)*sin(p2.af.Th)+pi.af.P*cos(pi.af.Ph)*sin(pi.af.Th)+k.af.P*cos(k.af.Ph)*sin(k.af.Th))**2+(p2.af.P*sin(p2.af.Th)*sin(p2.af.Ph)+pi.af.P*sin(pi.af.Th)*sin(pi.af.Ph)+k.af.P*sin(k.af.Th)*sin(k.af.Ph))**2+(p2.af.P*cos(p2.af.Th)+pi.af.P*cos(pi.af.Th)+k.af.P*cos(k.af.Th))**2)<2)");

T_bs->Draw("sqrt((p1.af.P*cos(p1.af.Ph)*sin(p1.af.Th)+pi.af.P*cos(pi.af.Ph)*sin(pi.af.Th)+k.af.P*cos(k.af.Ph)*sin(k.af.Th))**2+(p1.af.P*sin(p1.af.Th)*sin(p1.af.Ph)+pi.af.P*sin(pi.af.Th)*sin(pi.af.Ph)+k.af.P*sin(k.af.Th)*sin(k.af.Ph))**2+(p1.af.P*cos(p1.af.Th)+pi.af.P*cos(pi.af.Th)+k.af.P*cos(k.af.Th))**2)>>hbs1","w.pf*w.sk*w.lg*w.ppi*w.nf/50.");
 //*(((p1.af.P*cos(p1.af.Ph)*sin(p1.af.Th)+pi.af.P*cos(pi.af.Ph)*sin(pi.af.Th)+k.af.P*cos(k.af.Ph)*sin(k.af.Th))**2+(p1.af.P*sin(p1.af.Th)*sin(p1.af.Ph)+pi.af.P*sin(pi.af.Th)*sin(pi.af.Ph)+k.af.P*sin(k.af.Th)*sin(k.af.Ph))**2+(p1.af.P*cos(p1.af.Th)+pi.af.P*cos(pi.af.Th)+k.af.P*cos(k.af.Th))**2)<2)");

T_bs->Draw("sqrt((p2.af.P*cos(p2.af.Ph)*sin(p2.af.Th)+pi.af.P*cos(pi.af.Ph)*sin(pi.af.Th)+k.af.P*cos(k.af.Ph)*sin(k.af.Th))**2+(p2.af.P*sin(p2.af.Th)*sin(p2.af.Ph)+pi.af.P*sin(pi.af.Th)*sin(pi.af.Ph)+k.af.P*sin(k.af.Th)*sin(k.af.Ph))**2+(p2.af.P*cos(p2.af.Th)+pi.af.P*cos(pi.af.Th)+k.af.P*cos(k.af.Th))**2)>>hbs2","w.pf*w.sk*w.lg*w.ppi*w.nf/50.");
 //*(((p2.af.P*cos(p2.af.Ph)*sin(p2.af.Th)+pi.af.P*cos(pi.af.Ph)*sin(pi.af.Th)+k.af.P*cos(k.af.Ph)*sin(k.af.Th))**2+(p2.af.P*sin(p2.af.Th)*sin(p2.af.Ph)+pi.af.P*sin(pi.af.Th)*sin(pi.af.Ph)+k.af.P*sin(k.af.Th)*sin(k.af.Ph))**2+(p2.af.P*cos(p2.af.Th)+pi.af.P*cos(pi.af.Th)+k.af.P*cos(k.af.Th))**2)<2)");

T_dl->Draw("sqrt((p1.af.P*cos(p1.af.Ph)*sin(p1.af.Th)+pi.af.P*cos(pi.af.Ph)*sin(pi.af.Th)+k.af.P*cos(k.af.Ph)*sin(k.af.Th))**2+(p1.af.P*sin(p1.af.Th)*sin(p1.af.Ph)+pi.af.P*sin(pi.af.Th)*sin(pi.af.Ph)+k.af.P*sin(k.af.Th)*sin(k.af.Ph))**2+(p1.af.P*cos(p1.af.Th)+pi.af.P*cos(pi.af.Th)+k.af.P*cos(k.af.Th))**2)>>hdl","w.lk*w.ppi/100.");
 //*(((p1.af.P*cos(p1.af.Ph)*sin(p1.af.Th)+pi.af.P*cos(pi.af.Ph)*sin(pi.af.Th)+k.af.P*cos(k.af.Ph)*sin(k.af.Th))**2+(p1.af.P*sin(p1.af.Th)*sin(p1.af.Ph)+pi.af.P*sin(pi.af.Th)*sin(pi.af.Ph)+k.af.P*sin(k.af.Th)*sin(k.af.Ph))**2+(p1.af.P*cos(p1.af.Th)+pi.af.P*cos(pi.af.Th)+k.af.P*cos(k.af.Th))**2)<2)");

T_ds->Draw("sqrt((p1.af.P*cos(p1.af.Ph)*sin(p1.af.Th)+pi.af.P*cos(pi.af.Ph)*sin(pi.af.Th)+k.af.P*cos(k.af.Ph)*sin(k.af.Th))**2+(p1.af.P*sin(p1.af.Th)*sin(p1.af.Ph)+pi.af.P*sin(pi.af.Th)*sin(pi.af.Ph)+k.af.P*sin(k.af.Th)*sin(k.af.Ph))**2+(p1.af.P*cos(p1.af.Th)+pi.af.P*cos(pi.af.Th)+k.af.P*cos(k.af.Th))**2)>>hds","w.sk*w.lg*w.ppi/100.");
 //*(((p1.af.P*cos(p1.af.Ph)*sin(p1.af.Th)+pi.af.P*cos(pi.af.Ph)*sin(pi.af.Th)+k.af.P*cos(k.af.Ph)*sin(k.af.Th))**2+(p1.af.P*sin(p1.af.Th)*sin(p1.af.Ph)+pi.af.P*sin(pi.af.Th)*sin(pi.af.Ph)+k.af.P*sin(k.af.Th)*sin(k.af.Ph))**2+(p1.af.P*cos(p1.af.Th)+pi.af.P*cos(pi.af.Th)+k.af.P*cos(k.af.Th))**2)<2)");

//hbl2->SetTitle("Total momentum of pK^{+}#pi^{-}");
hbl2->SetXTitle("Total momentum: |P_{p}+P_{K^{+}}+P_{#pi^{-}}| (GeV/c)");
hbl2->SetYTitle("Arbitary Scale");
hbl2->GetXaxis()->SetLabelSize(0.035);
hbl2->GetYaxis()->SetLabelSize(0.035);
hbl2->SetLineColor(2);
hbl1->SetLineColor(6);
hbs2->SetLineColor(7);
hbs1->SetLineColor(1);
hdl->SetLineColor(4);
hds->SetLineColor(8);

hbl2->SetLineWidth(1);
hbl1->SetLineWidth(1);
hbs2->SetLineWidth(1);
hbs1->SetLineWidth(1);
hdl->SetLineWidth(1);
hds->SetLineWidth(1);

hbl2->SetMaximum(150);
hbl2->Draw();
hbl1->Draw("same");
hbs1->Draw("same");
hbs2->Draw("same");
hdl->Draw("same");
hds->Draw("same");

TLegend* leg = new TLegend(0.15,0.5,0.40,0.88,"Total momentum of pK^{+}#pi^{-}");
leg->SetFillColor(0);
leg->SetEntrySeparation(0.5);
leg->AddEntry(hbl2,"#phi-N to #Lambda^{o}","l");
leg->AddEntry(hbl1,"#phi-N to #Lambda^{o} w/ recoil p","l");
leg->AddEntry(hbs2,"#phi-N to #Sigma^{o}","l");
leg->AddEntry(hbs1,"#phi-N to #Sigma^{o} w/ recoil p","l");
leg->AddEntry(hdl,"Direct #Lambda^{o} Production","l");
leg->AddEntry(hds,"Direct #Sigma^{o} Production","l");
leg->Draw();
}