//Function To Calculate Eta
Float_t ETA(Float_t px, Float_t py, Float_t pz )
{
  Float_t r,theta,eta;
  r        = TMath::Sqrt(px*px + py*py);
  theta    = TMath::ATan2(r,pz);
  eta   = -TMath::Log(TMath::Tan(0.5*theta));
  return eta;
}
//Function To Calculate Eta
Float_t THETA(Float_t px, Float_t py, Float_t pz )
{
  Float_t r,theta,eta;
  r        = TMath::Sqrt(px*px + py*py);
  theta    = TMath::ATan2(r,pz);
  return theta;
}

Float_t PHI(Float_t px,Float_t py)
{
  Float_t phi,pybypx, phi1;
  Float_t pi = TMath::Pi();
  if(px==0)
    {
      if(py>0) phi = pi/2.;
      if(py<0) phi = (3.*pi/2.);
    }
  if(px != 0)
    {
      pybypx = py/px;
      if(pybypx<0) pybypx = - pybypx;
      phi1 = TMath::ATan(pybypx);
      if(px < 0 && py > 0) phi = pi - phi1;
      if(px < 0 && py < 0) phi = pi + phi1;
      if(px > 0 && py < 0) phi = 2.*pi - phi1;
      if(px > 0 && py > 0) phi = phi1;
    }
  return phi;
}

void TestCode()
{
  //Reset ROOT and connect tree file
  gROOT->Reset();
  TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("Kinematics.root");
  if (!f) {
    f = new TFile("Kinematics.root");
  }
 
  TTree   *tree  = (TTree*)f->Get("Event0/TreeK");//Tree
  TBranch *branch1 = tree->GetBranch("fPx");//Branch
  TBranch *branch2 = tree->GetBranch("fPy");//Branch
  TBranch *branch3 = tree->GetBranch("fPz");//Branch
  TBranch *branch4 = tree->GetBranch("fE");//Branch
  TBranch *branch5 = tree->GetBranch("fCalcMass");//Branch
  TBranch *branch6 = tree->GetBranch("fDaughter[2]");//Branch
  TBranch *branch7 = tree->GetBranch("fPolarPhi");//Branch
  TBranch *branch8 = tree->GetBranch("fPolarTheta");//Branch
    



  TCanvas *c1 = new TCanvas("c1","c1",200,10,800,600);
  c1->Divide(2,2);
  TCanvas *c2 = new TCanvas("c2","c2",200,10,800,600);
  c2->Divide(2,2);


  TH1F *hPx = new TH1F("Px","Px",100,-2.0,4.0); 
  TH1F *hPy = new TH1F("Py","Py",100,-2.0,1.5); 
  TH1F *hPz = new TH1F("Pz","Pz",100,-60.0,140.0); 
  TH1F *hE = new TH1F("E","E",100,0.0,140.0); 

  TH1F *hEta = new TH1F("Eta","Eta",100,-7.5,7.5); 
  TH1F *hTheta = new TH1F("Theta","Theta",100,0.0,3.15); 
  TH1F *hPhi = new TH1F("Phi","Phi",100,0.0,6.5); 
  
  //Declaration of leaves types
  Double_t        fPx;
  Double_t        fPy;
  Double_t        fPz;
  Double_t        fE;
  Double_t        fCalcMass;
  Double_t        fDaughter[2];
  Double_t        fPolarPhi;
  Double_t        fPolarTheta;
  
  Double_t        Eta;
  Double_t        Phi;
  Double_t        Theta;

  //the following call is required when reading a branch of a superbranch
  tree->SetMakeClass(1);
  
  branch1->SetAddress(&fPx);
  branch2->SetAddress(&fPy);
  branch3->SetAddress(&fPz);
  branch4->SetAddress(&fE);
  branch5->SetAddress(&fCalcMass);
  branch6->SetAddress(fDaughter[0]);
  branch6->SetAddress(fDaughter[1]);
  branch7->SetAddress(fPolarPhi);
  branch8->SetAddress(fPolarTheta);
  
  Long64_t nentries = branch1->GetEntries();
  cout<<"Entries = "<<nentries<<endl;
  
  Long64_t nbytes = 0;
  
  for (Long64_t i=0; i<nentries;i++) 
    {
      branch1->GetEntry(i);
      branch2->GetEntry(i);
      branch3->GetEntry(i);
      branch4->GetEntry(i);
      branch5->GetEntry(i);
      branch6->GetEntry(i);
      branch7->GetEntry(i);
      branch8->GetEntry(i);
      /*
      cout<<"Px["<<i+1<<"] = "<<fPx<<endl;
      cout<<"Py["<<i+1<<"] = "<<fPy<<endl;
      cout<<"Pz["<<i+1<<"] = "<<fPz<<endl;
      cout<<"E["<<i+1<<"] = "<<fE<<endl;
      cout<<"Calculated Mass["<<i+1<<"] = "<<fCalcMass<<endl;
      cout<<"Daughter1["<<i+1<<"] = "<<fDaughter[0]<<endl;
      cout<<"Daughter2["<<i+1<<"] = "<<fDaughter[1]<<endl;
      cout<<"PolarPhi["<<i+1<<"] = "<<fPolarPhi<<endl;
      cout<<"PolarTheta["<<i+1<<"] = "<<fPolarTheta<<endl;
      */
      Eta = ETA(fPx,fPy,fPz);
      Theta = THETA(fPx,fPy,fPz);
      Phi = PHI(fPx,fPy);
      
      hPx->Fill(fPx);
      hPy->Fill(fPy);
      hPz->Fill(fPz);
      hE->Fill(fE);

      hEta->Fill(Eta);
      hPhi->Fill(Phi);
      hTheta->Fill(Theta);
    }
  c1->cd(1);
  hPx->Draw();
  c1->cd(2);
  hPy->Draw();
  c1->cd(3);
  hPz->Draw();
  c1->cd(4);
  hE->Draw();

  c2->cd(1);
  hEta->Draw();
  c2->cd(2);
  hTheta->Draw();
  c2->cd(3);
  hPhi->Draw();
}