#/*********************************************************************
 comparison between lumi and beam angles                  
**********************************************************************/

void lumi_beam_comp(TString dirName) {

  TString filePath = dirName;
  TNtuple *lnt = makeGpNtuple(filePath+ "lumi.ee.out");
  TNtuple *b1nt = makeGpNtuple(filePath+ "beam1.dat");
  TNtuple *b2nt = makeGpNtuple(filePath+ "beam2.dat");

  TH1D *lnt_vx = new TH1D("lnt_vx","lnt_vx",100,-0.001,0.001);
  TH1D *lnt_vy = new TH1D("lnt_vy","lnt_vy",100,-0.001,0.001);
  TH1D *b1nt_vx = new TH1D("b1nt_vx","b1nt_vx",100,-0.001,0.001);
  TH1D *b1nt_vy = new TH1D("b1nt_vy","b1nt_vy",100,-0.001,0.001);

  lnt_vx->Sumw2();
  lnt_vy->Sumw2();
  b1nt_vx->Sumw2();
  b1nt_vy->Sumw2();

  lnt->Draw("vx1 >> lnt_vx","","");
  lnt->Draw("vy1 >> lnt_vy","","");
  b1nt->Draw("vx >> b1nt_vx");
  b1nt->Draw("vy >> b1nt_vy","","");

  TCanvas *can = new TCanvas("c1","c1",600,400);
  lnt_vx->Scale(1./lnt_vx->Integral());
  lnt_vy->Scale(1./lnt_vy->Integral());
  b1nt_vx->Scale(1./b1nt_vx->Integral());
  b1nt_vy->Scale(1./b1nt_vy->Integral());

  lnt_vx->SetMarkerStyle(20);
  lnt_vy->SetMarkerStyle(24);
  lnt_vx->SetMaximum(TMath::Max(lnt_vy->GetMaximum(),b1nt_vy->GetMaximum())*1.2);
  lnt_vx->SetMinimum(0.00001);
  lnt_vx->Fit("gaus");

   gPad->SetLogy();
  lnt_vx->Draw("e1");
  lnt_vy->Draw("e1same");
  b1nt_vx->Draw("histsame");
  b1nt_vy->Draw("histsame");

  lnt_vx->GetXaxis()->SetTitle("#phi_{x,y}[radians]");
  lnt_vx->GetYaxis()->SetTitle("events");

  can->Print("../tmp-results/lumi_beam_comp.eps");
}

/*********************************************************************
 Convert angles to polar and plot
**********************************************************************/

void lumi_polar(TString fName) {

  TNtuple *lnt = makeGpNtuple(fName.Data());

  TH1D *theta1 = new TH1D("theta1","theta1",100,0,0.0005);
  TH1D *theta2 = new TH1D("theta2","theta2",100,0,0.0005);
  TH1D *phi1   = new TH1D("phi1","phi1",100,-TMath::Pi(),TMath::Pi());
  TH1D *phi2   = new TH1D("phi2","phi2",100,-TMath::Pi(),TMath::Pi());

  lnt->Draw("asin(sqrt(pow(sin(vx1),2)+pow(sin(vy1),2))) >> theta1","","");
  lnt->Draw("atan2(sin(vy1),sin(vx1)) >> phi1","","");
  lnt->Draw("asin(sqrt(pow(sin(vx2),2)+pow(sin(vy2),2))) >> theta2","","");
  lnt->Draw("atan2(sin(vy2),sin(vx2)) >> phi2","","");

  TCanvas *can = new TCanvas("c1","c1",600,800);
  can->Divide(1,2);

  theta1->SetMarkerStyle(20);
  phi1->SetMarkerStyle(20);

  can->cd(1);
  theta1->GetXaxis()->SetTitle("#theta_{1,2}[radians]");
  theta1->GetYaxis()->SetTitle("events");
  theta1->Draw("e1");
  theta2->Draw("histsame");

  can->cd(2);
  phi1->GetXaxis()->SetTitle("#phi_{1,2}[radians]");
  phi1->GetYaxis()->SetTitle("events");
  phi1->Draw("e1");
  phi2->Draw("histsame");

  can->Print("../tmp-results/lumi_polar.eps");
}

/*********************************************************************
 Convert angles to polar and compute acolinearity and acoplanarity
**********************************************************************/

void lumi_aco(TString fName) {



  TNtuple *lnt = makeGpNtuple(fName);

  TH1D *theta = new TH1D("theta","theta",100,-1,1);
  //TH1D *theta = new TH1D("theta","theta",100,-0.0005,0.0005);
  TH1D *phi   = new TH1D("phi","phi",100,-TMath::Pi(),TMath::Pi());

  lnt->Draw("asin(sqrt(pow(sin(vx1),2)+pow(sin(vy1),2)))-asin(sqrt(pow(sin(vx2),2)+pow(sin(vy2),2))) >> theta","","");
  lnt->Draw("atan2(sin(vy1),sin(vx1))-atan2(sin(vy2),sin(vx2)) >> phi","","");

  TCanvas *can = new TCanvas("c1","c1",600,800);
  can->Divide(1,2);

  can->cd(1);
  theta->SetMarkerStyle(20);
  theta->GetXaxis()->SetTitle("#Delta#theta[radians]");
  theta->GetYaxis()->SetTitle("events");
  theta->Draw("e1");

  can->cd(2);
  phi->SetMarkerStyle(20);
  phi->GetXaxis()->SetTitle("#Delta#phi[radians]");
  phi->GetYaxis()->SetTitle("events");
  phi->Draw("e1");

  can->Print("../tmp-results/lumi_aco.eps");
}

/*********************************************************************
 angle correlation
**********************************************************************/

void lumi_angle_correl(TString fName) {
  TNtuple *lnt = makeGpNtuple(fName.Data());
  lnt->Draw("2*sqrt(e1*e2) : (atan2(sin(vy1),sin(vx1))-atan2(sin(vy2),sin(vx2))) ","","prof");
  //  lnt->Draw("2*sqrt(e1*e2) : (atan2(sin(vy1),sin(vx1))-atan2(sin(vy2),sin(vx2))) ","","");
}


/*********************************************************************
 missing pt
**********************************************************************/

void lumi_missing_pt(TString fName) {
  TNtuple *lnt = makeGpNtuple(fName.Data());
  Float_t e1,e2,vx1,vy1,vx2,vy2;
  lnt->SetBranchAddress("e1",&e1);
  lnt->SetBranchAddress("e2",&e2);
  lnt->SetBranchAddress("vx1",&vx1);
  lnt->SetBranchAddress("vy1",&vy1);
  lnt->SetBranchAddress("vx2",&vx2);
  lnt->SetBranchAddress("vy2",&vy2);

  TH1D *ptmiss = new TH1D("ptmiss","ptmiss",100,0,1.0);

  for(Int_t i=0;i<lnt->GetEntries();i++) {
    lnt->GetEntry(i);
    Float_t theta1 = asin(sqrt(pow(sin(vx1),2)+pow(sin(vy1),2)));
    Float_t theta2 = asin(sqrt(pow(sin(vx2),2)+pow(sin(vy2),2)));
    Float_t phi1   = atan2(sin(vy1),sin(vx1));
    Float_t phi2   = atan2(sin(vy2),sin(vx2));

    TLorentzVector p1 = TLorentzVector(e1*sin(theta1)*cos(phi1),
                   e1*sin(theta1)*sin(phi1),
                   e1*cos(theta1),e1);

    TLorentzVector p2 = TLorentzVector(e2*sin(theta2)*cos(phi2),
                   e2*sin(theta2)*sin(phi2),
                   e2*cos(theta2),e2);

    ptmiss->Fill((p1+p2).Pt());
  }

  TCanvas *can = new TCanvas("c1","c1",600,400);
  ptmiss->Draw();

}