#include<iostream>
#include<fstream>
#include<cmath>

using namespace std;
//par0=norm, par1=mean, x = x[0]

double poissonf(Double_t*x,Double_t*par) //fitting function
{
double formula;
formula= (par[0]*TMath::Power(par[1], x[0])*TMath::Exp(-par[1]))/(TMath::Factorial(x[0]));

return (formula);
}




void ChargedKStar8()
{


  
TH1D *hist = new TH1D("hist","hist",90,0,12);
TH1D *histi = new TH1D("histi","histi",900,0.2,10.0);
TH1D *histo = new TH1D("histo","histo",200,0,2);
TH1D *histo2 = new TH1D("histo2","histo2",200,0,2);
TH1D *histcos = new TH1D("histcos","histcos",20,-1,1);
 
 float energyK4,massK4,ptK4,etaK4,phiK4,pzK4,energyK3,massK3,ptK3,etaK3,phiK3,pzK3;

TLorentzVector KE;
TLorentzVector *Rho;
TLorentzVector *Kaon;
TLorentzVector *PI1;
TLorentzVector *PI2;
TVector3 v1;
TVector3 v2;

TLorentzVector rho;
TLorentzVector kaon;
TLorentzVector pi1;
TLorentzVector pi2;
TLorentzVector pi;
TLorentzVector rk;

TGenPhaseSpace event;
Double_t weight;
TGenPhaseSpace event2;
Double_t weight2;
Double_t rho_px;



TF1 *fE=new TF1("fE",poissonf,0,10,2);
fE->FixParameter(0,1);
fE->FixParameter(1,1.864);

fE->FixParameter(0,1.864);
fE->FixParameter(1,0.005);
fE->FixParameter(2,1);

TF1 *f=new TF1("f",poissonf,0,10,2);
f->FixParameter(0,4);
f->FixParameter(1,4);


 TF1 *fPT1=new TF1("fPT1","exp(-x*x)",0,5);

TH1D* histPt=new TH1D("histPt","myhist",1000,0,10);
TH1D* histEta=new TH1D("histEta","myhistnew",1000,-1,1);

 double massKE,ptKE,etaKE,phiKE,momentumKE,energyKE,cosdipangle,cosopen,massrho,cosfir,pxKE,pyKE,pzKE,pxrho,pyrho,pzrho,pxkaon,pykaon,pzkaon,InvMass,cosfiro;
TRandom *eventGenerator = new TRandom();

for (int i=0;i<20;i++)
 
{


massKE=fE->GetRandom();



etaKE=gRandom->Uniform(-1.0,1.0);


pxKE=eventGenerator->Gaus(0.0,2.0);
pyKE=eventGenerator->Gaus(0.0,2.0);
pzKE=eventGenerator->Gaus(0.0,2.0);
pxrho=eventGenerator->Gaus(0.0,2.0);

pyrho=eventGenerator->Gaus(0.0,2.0);
pzrho=eventGenerator->Gaus(0.0,2.0);


pyKE=eventGenerator->Gaus(0.0,2.0);
pzKE=eventGenerator->Gaus(0.0,2.0);
 momentumKE=sqrt(pxKE*pxKE+pyKE*pyKE+pzKE*pzKE);
 energyKE=sqrt((momentumKE*momentumKE) + (massKE*massKE));



 KE.SetPxPyPzE(pxKE,pyKE,pzKE,energyKE);  





Double_t masses4[2] = {0.139,0.139};//masses of decaying particle 
event.SetDecay(KE,2,masses4);
weight= event.Generate();
Rho = event.GetDecay(0);
Kaon  = event.GetDecay(1);
rho=*Rho;
kaon=*Kaon;

rk=rho+kaon;


TVector3 v1(rho.Px(),rho.Py(),rho.Pz());
TVector3 v2(kaon.Px(),kaon.Py(),kaon.Pz());

 
//hist->Fill(sqrt(KE*KE));
//hist->Fill(massKE);
histi->Fill(massrho);
 if (rho.Pt()<0.2)
 histo->Fill(kaon.Pt());
 if (kaon.Pt()<0.2)
 histo2->Fill(rho.Pt());

 cout<<KE.Pt()<<"--"<<rho.Px()<<"--"<<kaon.Px()<<"--"<<rk.Px()<<":Loop:"<<i<<"--"<<"Diff"<<KE.Pt()-rho.Pt()<<endl;


histPt->Fill(KE.Pt());
histEta->Fill(etaKE);
//hist->Fill(rk.M());
//histi->Fill(rho.px());
 rho_px= KE.Px();
hist->Fill(KE.Px());
hist->Draw();






 if (i%10000==0)
   cout<<"Event loop:"<<i<<endl;


 }



 TH1D *h_dummy=new TH1D("h_dummy","h_dummy",1,0.0,2.0);
 h_dummy->SetMinimum(0);
 h_dummy->SetMaximum(2000);
 
 
    
    
 TH1D *h_pt_point5=Plothisto(histo,kBlue, 24); //v2default
 TH1D *h_pt_1=Plothisto(histo2,kRed, 24); //v2default

 TCanvas *c11 = DrawCanvas("c11");
 //c11->Divide(2,1);
 c11->cd(1);
 TPad *pad = MyPad(0.01,0.01,0.99, 0.99, 0.15,0.02,0.02,0.13, "pad");
 pad->cd();
  
  


// h_dummy->Draw("X0");
// h_pt_point5->Draw("same");
// h_pt_1->Draw("same");
 
// c11->cd(2);
 //hist->Draw();
 //histEta->Draw();

//histi->Draw();
 

}




TCanvas * DrawCanvas(TString opt="c")
{
  TCanvas *c1 = new TCanvas(opt.Data(),opt.Data(),10,10,600,600);
  c1->cd(1);
  c1->SetLeftMargin(0.15);
  c1->SetRightMargin(0.05);
  c1->SetTopMargin(0.05);
  c1->SetBottomMargin(0.15);
  return c1;
}
TPad *MyPad(Double_t x1 = 0.01,Double_t y1 = 0.01,Double_t x2 = 0.49, Double_t y2= 0.49, Float_t lm  =0.2,Float_t rm = 0.05,Float_t tm =0.05,Float_t bm = 0.2,TString name = "pad")
{
  TPad *c1_1 = new TPad(name.Data(), name.Data(),x1,y1,x2,y2);
  c1_1->Draw();
  c1_1->cd();
  c1_1->Range(0,0,1,1);
  c1_1->SetBorderSize(2);
  c1_1->SetBorderMode(0);
  c1_1->SetFillColor(10);
  c1_1->SetFrameFillColor(10);
  c1_1->SetFrameLineWidth(2);
  c1_1->SetLeftMargin(lm);
  c1_1->SetRightMargin(rm);
  c1_1->SetTopMargin(tm);
  c1_1->SetBottomMargin(bm);
  c1_1->SetTicks(1,1);
  return c1_1;
}
TGraphErrors * PlotGraph(Int_t NdataPoint, Int_t MarkerColor = 1, Int_t MarkerStyle = 20, Double_t *X, Double_t *ErX,Double_t *Y, Double_t *ErY)
{
  TGraphErrors *gr = new TGraphErrors(NdataPoint, X ,Y , ErX, ErY);
  gr->SetTitle("");
  gr->SetMarkerStyle(MarkerStyle);
  gr->SetMarkerColor(MarkerColor);
  gr->SetMarkerSize(1.2);
  gr->SetLineWidth(2);
  gr->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
  gr->GetXaxis()->CenterTitle(false);
  gr->GetXaxis()->SetNdivisions(506);
  gr->GetYaxis()->SetNdivisions(505);
  gr->GetXaxis()->SetLabelOffset(0.015);
  gr->GetXaxis()->SetLabelFont(42);
  gr->GetXaxis()->SetTitleFont(42);
  gr->GetXaxis()->SetLabelSize(0.05);
  gr->GetXaxis()->SetTitleSize(0.05);
  gr->GetXaxis()->SetTickLength(0.04);
  gr->GetXaxis()->SetTitleOffset(1.25);
  gr->GetYaxis()->SetTitleOffset(1.45);
  gr->GetYaxis()->CenterTitle(false);
  gr->GetYaxis()->SetDecimals(false);
  //gr->GetYaxis()->SetNdivisions(310);
  gr->GetYaxis()->SetLabelOffset(0.015);
  gr->GetYaxis()->SetLabelFont(42);
  gr->GetYaxis()->SetLabelSize(0.05);
  gr->GetYaxis()->SetTickLength(0.04);
  gr->GetYaxis()->SetTitleSize(0.05);
  gr->GetYaxis()->SetTitle("#it{#rho}_{00}"); 
  gr->GetYaxis()->SetTitleFont(42);
  return gr;
}


TH1D * Plothisto(TH1D *h1, Int_t MarkerColor = 1, Int_t MarkerStyle = 20)
{
  
  h1->SetTitle("");
  h1->SetMarkerStyle(MarkerStyle);
  h1->SetMarkerColor(MarkerColor);
  h1->SetLineColor(MarkerColor);
  
  h1->SetMarkerSize(1.2);
  h1->SetLineWidth(2);
  h1->GetXaxis()->SetTitle("#it{p}_{T}");
  h1->GetXaxis()->CenterTitle(false);
  h1->GetXaxis()->SetNdivisions(506);
  h1->GetYaxis()->SetNdivisions(505);
  h1->GetXaxis()->SetLabelOffset(0.015);
  h1->GetXaxis()->SetLabelFont(42);
  h1->GetXaxis()->SetTitleFont(42);
  h1->GetXaxis()->SetLabelSize(0.05);
  h1->GetXaxis()->SetTitleSize(0.05);
  h1->GetXaxis()->SetTickLength(0.04);
  h1->GetXaxis()->SetTitleOffset(1.15);
  h1->GetYaxis()->SetTitleOffset(1.15);
  h1->GetYaxis()->CenterTitle(true);
  h1->GetYaxis()->SetDecimals(false);
  //h1->GetYaxis()->SetNdivisions(310);
  h1->GetYaxis()->SetLabelOffset(0.015);
  h1->GetYaxis()->SetLabelFont(42);
  h1->GetYaxis()->SetLabelSize(0.05);
  h1->GetYaxis()->SetTickLength(0.04);
  h1->GetYaxis()->SetTitleSize(0.06);
  h1->GetYaxis()->SetTitle("Counts");
  //h1->GetYaxis()->SetTitle("Y"); 
  h1->GetYaxis()->SetTitleFont(42);
  return h1;
}






//==================================
TLegend *DrawLegend(Double_t x1,Double_t y1,Double_t x2,Double_t y2)
{
  TLegend *legend = new TLegend(x1,y1,x2,y2);
  legend->SetTextFont(42);
  legend->SetTextSize(0.03);
  legend->SetLineColor(0);
  legend->SetShadowColor(0);
  return legend;
}
//=====================================
TLatex *DrawText(Double_t x = 0, Double_t y = 0,Int_t tColor = 2,TString name)
{
  TLatex* tex = new TLatex(x,y,name.Data());
  tex->SetTextSize(0.04);
  tex->SetTextColor(tColor);
  tex->SetTextFont(42);
  //tex->Draw();
  return tex;
}
//==================================
void MyStyle()
{
  gStyle->SetOptTitle(0);
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(0);
  gStyle->SetStatColor(0);
  gStyle->SetCanvasColor(0);
  gStyle->SetPadColor(0);
  gStyle->SetPadBorderMode(0);
  gStyle->SetCanvasBorderMode(0);
  gStyle->SetFrameBorderMode(0);
  gStyle->SetFillColor(0);
  gStyle->SetLineColor(1);
  gStyle->SetTitleTextColor(1);
  gStyle->SetTitleColor(1);
  gStyle->SetLineWidth(2);
  gStyle->SetErrorX(0.2);

}

TGraphErrors *DrawFrame(TGraphErrors *h, Int_t MCol,Int_t MSty, Bool_t mrk = 0)
{
  h->GetXaxis()->SetTitle("<#it{N}_{part}>");
  if(mrk){
    h->SetMarkerColor(MCol);
    h->SetMarkerStyle(MSty);
    h->SetMarkerSize(1.4);
  }
  h->SetLineColor(MCol);
  h->SetLineWidth(2);
  h->GetXaxis()->CenterTitle(false);
  h->GetXaxis()->SetNdivisions(505);
  h->GetYaxis()->SetNdivisions(505);
  h->GetXaxis()->SetLabelOffset(0.017);
  h->GetXaxis()->SetLabelFont(42);
  h->GetXaxis()->SetTitleFont(42);
  h->GetXaxis()->SetLabelSize(0.04);
  h->GetXaxis()->SetTitleSize(0.05);
  h->GetXaxis()->SetTickLength(0.06);
  h->GetXaxis()->SetTitleOffset(1.0);
  h->GetYaxis()->SetTitleOffset(1.15);
  h->GetYaxis()->CenterTitle(false);
  h->GetYaxis()->SetDecimals(false);
  h->GetYaxis()->SetLabelOffset(0.015);
  h->GetYaxis()->SetLabelFont(42);
  h->GetYaxis()->SetLabelSize(0.04);
  h->GetYaxis()->SetTickLength(0.04);
  h->GetYaxis()->SetTitleSize(0.05);
  h->GetYaxis()->SetTitle("#it{#rho_{00}}"); 
  h->GetYaxis()->SetTitleFont(42);
  //h->Draw("AP");
  return h;
}