#include "TFile.h"
#include "TTree.h"
#include "TMath.h"
#include "TLeaf.h"

int convert_carv_to_polar() {

  TFile * fin = new TFile("mnt_density_data.root", "READ");

  TTree * oldtree = (TTree*)fin->Get("tree_mnt");

  int Ntotal = oldtree->GetEntries();

  double x,y,z;
  double r,theta,phi;
  double rho, drho;

  double det_x = -150;
  double det_y =  350;
  double det_z =    0;

  TFile * fout = new TFile("mnt_density_data_converted.root", "RECREATE");
  TTree * tree = new TTree("tree", "");
  tree->Branch("x", &x, "x/D");
  tree->Branch("y", &y, "y/D");
  tree->Branch("z", &z, "z/D");
  tree->Branch("rho", &rho, "rho/D");
  tree->Branch("drho", &drho, "drho/D");
  tree->Branch("theta", &theta, "theta/D");
  tree->Branch("phi", &phi, "phi/D");
  tree->Branch("r", &r, "r/D");

  
  for (int i=0; i<Ntotal; i++) {
    
    oldtree->GetEntry(i);
    
    x = oldtree->GetLeaf("x")->GetValue();
    y = oldtree->GetLeaf("y")->GetValue();
    z = oldtree->GetLeaf("z")->GetValue();
    rho = oldtree->GetLeaf("rho")->GetValue();
    drho = rho-2.7;
    
    x = x - det_x;
    y = y - det_y;
    z = z - det_z;
    r = sqrt(x*x+y*y+z*z);
    
    theta = TMath::ASin(sqrt(x*x+y*y)/r);
    if (z<0)  theta = - theta + TMath::Pi();
    phi = TMath::ATan(y/x);
    if ( x<0 && y>0) phi =   TMath::Pi() - fabs(phi);
    if ( x<0 && y<0) phi = - TMath::Pi() + fabs(phi);
    if ( x>0 && y<0) phi = - fabs(phi);

    tree->Fill();
    
  }
  tree->Write();
  fout->Close();

}