/*****************************************
Compare gp output: 
phi, theta, vx, vy distributions in 2-D 
density histograms.
Arguments (can be concatenated): 

1.)opt = "vs"  plot  vx,vy
    =  "lumi"   plot  theta, phi
    =  "comp" plot theta1 vs theta2 etc

2.)rs_nom = give the center of mass energy

 ****************************************/

void lumi_angles_comp(TString fName, TString opt="", Double_t rs_nom=350.0) {

  const Double_t pi=TMath::Pi();

  TNtuple *lnt = makeGpNtuple(fName.Data());
  cout << "finished loading file" << endl;


  Float_t e1,e2,vx1,vy1,vx2,vy2;
  Float_t ievent,x,y,z;
  lnt->SetBranchAddress("ievent",&ievent);
  lnt->SetBranchAddress("e1",&e1);
  lnt->SetBranchAddress("e2",&e2);
  lnt->SetBranchAddress("x",&x);
  lnt->SetBranchAddress("y",&y);
  lnt->SetBranchAddress("z",&z);
  lnt->SetBranchAddress("vx1",&vx1);
  lnt->SetBranchAddress("vy1",&vy1);
  lnt->SetBranchAddress("vx2",&vx2);
  lnt->SetBranchAddress("vy2",&vy2);


  TH1D *hvx1 = new TH1D("hvx1","vx1 histogram",101,-0.0005,0.0005);
  TH1D *hvx2 = new TH1D("hvx2","vx2 histogram",101,-0.0005,0.0005);
  TH1D *hvy1 = new TH1D("hvy1","vy1 histogram",101,-0.0002,0.0002);
  TH1D *hvy2 = new TH1D("hvy2","vy2 histogram",101,-0.0002,0.0002);

  //these are the mirror images histograms
  TH1D *hvx1s = new TH1D("hvx1s","vx1 histogram",101,-0.0005,0.0005);
  TH1D *hvx2s = new TH1D("hvx2s","vx2 histogram",101,-0.0005,0.0005);
  TH1D *hvy1s = new TH1D("hvy1s","vy1 histogram",101,-0.0002,0.0002);
  TH1D *hvy2s = new TH1D("hvy2s","vy2 histogram",101,-0.0002,0.0002);

  TH1D *hth1 = new TH1D("hth1","theta1 histogram",100,0,0.0005);
  TH1D *hth2 = new TH1D("hth2","theta2 histogram",100,0,0.0005);
  TH1D *hph1 = new TH1D("hph1","phi1 histogram",100,-pi,pi);
  TH1D *hph2 = new TH1D("hph2","phi2 histogram",100,-pi,pi);


  TH2D *th1th2 = new TH2D("th1th2","theta1 vs theta2",100,-0.00001,0.0005,100,-0.00001,0.0005);
  TH2D *ph1ph2 = new TH2D("ph1ph2","phi1 vs phi2",100,-pi,pi,100,-pi,pi);
  TH2D *vx1vx2 = new TH2D("vx1vx2","vx1 vs vx2",100,-0.0005,0.0005,100,-0.0005,0.0005);
  TH2D *vx1vy1 = new TH2D("vx1vy1","vx1 vs vy1",100,-0.0005,0.0005,100,-0.0005,0.0005);

  for(Int_t i=0;i<lnt->GetEntries();i++) {
    lnt->GetEntry(i);
    Double_t theta1 = asin(sqrt(pow(sin(vx1),2)+pow(sin(vy1),2)));
    Double_t theta2 = asin(sqrt(pow(sin(vx2),2)+pow(sin(vy2),2)));
    Double_t phi1   = atan2(sin(vy1),sin(vx1));
    Double_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);

    hvx1->Fill(vx1);
    hvx2->Fill(vx2);
    hvy1->Fill(vy1);
    hvy2->Fill(vy2);

    hvx1s->Fill(-vx1);
    hvx2s->Fill(-vx2);
    hvy1s->Fill(-vy1);
    hvy2s->Fill(-vy2);

    hth1->Fill(theta1);
    hth2->Fill(theta2);
    hph1->Fill(phi1);
    hph2->Fill(phi2);

    th1th2->Fill(theta1,theta2);
    ph1ph2->Fill(phi1,phi2);
    vx1vx2->Fill(vx1,vx2);
    vx1vy1->Fill(vx1,vy1);

  }
  cout << "finished filling histograms" << endl;

  //draw the stuff, with root-style if()
  if(opt.Contains("surf")){
    TCanvas *can2 = new TCanvas("can2","can2",600,600);
    gStyle->SetOptTitle(1);
    gStyle->SetOptStat(0);
    can2->Divide(2,2);
    can2->cd(1);
    vx1vx2->Draw();
    can2->cd(2);
    vx1vy1->Draw();
    can2->cd(3);
    th1th2->Draw();
    can2->cd(4);
    ph1ph2->Draw();
  }
  if(opt.Contains("vs")){
    TCanvas *can1 = new TCanvas("can1","can1",600,600);
    gStyle->SetOptTitle(1);
    gStyle->SetOptStat(0);
    can1->Divide(2,2);

    can1->cd(1);
    hvx1->Draw();
    hvx1s->SetLineColor(3);
    hvx1s->SetLineStyle(3);
    hvx1s->Draw("same");

    can1->cd(2);
    hvx2->Draw();
    hvx2s->SetLineColor(3);
    hvx2s->SetLineStyle(3);
    hvx2s->Draw("same");

    can1->cd(3);
    hvy1->Draw();
    hvy1s->SetLineColor(3);
    hvy1s->SetLineStyle(3);
    hvy1s->Draw("same");

    can1->cd(4);
    hvy2->Draw();
    hvy2s->SetLineColor(3);
    hvy2s->SetLineStyle(3);
    hvy2s->Draw("same");


  }
  if(opt.Contains("angles")){
    TCanvas *can3 = new TCanvas("can3","can3",600,600);
    gStyle->SetOptTitle(1);
    gStyle->SetOptStat(0);
    can3->Divide(2,2);
    can3->cd(1);
    hth1->Draw();
    hth2->SetLineColor(3);
    hth2->SetLineStyle(3);
    hth2->Draw("same");
    can3->cd(2);
    hth2->Draw();
    can3->cd(3);
    hph1->Draw();
    hph2->SetLineColor(3);
    hph2->SetLineStyle(3);
    hph2->Draw("same");
    can3->cd(4);
    hph2->Draw();
  }

}