#include "Trk.hh"
#include "Vec.hh"
#include "Hit.hh"

#include "../utilJpp.hh"
#include "TH1D.h"
#include "TCanvas.h"
#include "TFile.h"

TFile f("htheta.root");
TH1D htheta( "hitparams - params theta", "hitparams - params theta", 1000, -0.1, 3.2);

int main( ){
  
  long double E = 1e5;
  long double D = 10;
  long double cd = 0.1;
  long double phi = 0.1;
  long double dt = 1.0;

  Trk trk;
  trk.E = E;
  trk.pos = Vec( 0, 0, 0 );
  trk.dir = Vec( 0, 0, 1 );
  
  Hit hit;
  hit.pos = trk.pos + Vec().set_angles( acos(cd), 0 ).normalize() * D;
  hit.t = D / 0.22;

  for( uint thetaBin = 1; thetaBin != htheta.GetNbinsX(); thetaBin++ ){

    hit.dir.set_angles( htheta.GetBinCenter( thetaBin ), phi );

    //Method 1
    const long double u = hit.dir.dot( trk.dir );
    long double theta_1 = acos(u);

    //Method 2
    JPosition3D Jpos = getPosition( trk );
    JDirection3D Jdir = getDirection( trk );

    const JRotation3D Rot( Jdir );    
    Jpos.rotate(Rot);
    JAxis3D axis( JPosition3D(  hit.pos.x, hit.pos.y, hit.pos.z ),
		  JDirection3D( hit.dir.x, hit.dir.y, hit.dir.z ) );
    axis.transform( Rot, Jpos);
    long double theta_2 = axis.getTheta();

    long double diff = abs( theta_1 - theta_2 );
    htheta.SetBinContent( thetaBin, diff );      
  }

  TCanvas c1( "c1", "c1", 800, 800 );
  htheta.GetYaxis()->SetTitle( "|HitParams theta - Params theta|" );
  htheta.GetXaxis()->SetTitle( "theta" );
  htheta.Draw();
  c1.SetLogy();
  c1.SaveAs( "compare_theta.pdf" );
  f.Write();

}