#include <globals.hh>
#include <Randomize.hh>

#include <RAT/MetaSphereSeed.hh>
#include <RAT/Method.hh>
#include <RAT/OptimiserFactory.hh>
#include <RAT/Processor.hh>
#include <RAT/DS/FitResult.hh>
#include <RAT/DB.hh>
using namespace RAT::DS;
using namespace RAT::Optimisers;

#include <CLHEP/Units/PhysicalConstants.h>
#include <string>
using namespace CLHEP;
using namespace std;

void
MetaSphereSeed::Initialise( const std::string& param )
{
  fRadius = 6060.4;
  MetaNSeed::Initialise( param );
  fNSeeds = 12;
}

void
MetaSphereSeed::SetD( const std::string& param,
                      double value )
{
  if( param == string( "radius" ) )
    fRadius = value;
  else
    throw Processor::ParamUnknown( param );
}

FitResult
MetaSphereSeed::NewSeed( const FitResult& startSeed )
{
  FitResult newSeed = startSeed;
  FitVertex vertex = newSeed.GetVertex(0);
  if( vertex.ContainsPosition() )
	{
	  const double cosTheta = -1.0 + 2.0 * G4UniformRand();
	  const double phi = twopi * G4UniformRand();
	  const double r = pow( G4UniformRand(), 1.0 / 3.0 ) * fRadius;
	  vertex.SetPosition( TVector3( r * sin( acos( cosTheta ) ) * cos( phi ),
									 r * sin( acos( cosTheta ) ) * sin( phi ),
									 r * cosTheta ) );
	  vertex.SetPositionErrors( ROOT::Math::XYZVectorF(fRadius, fRadius, fRadius) );
	}
  if( vertex.ContainsTime() )
	{
	  vertex.SetTime( 240.0 + G4UniformRand() * 20.0 );
	  vertex.SetTimeErrors( 20.0 );
	}
  if( vertex.ContainsDirection() )
	{
	  const double cosTheta = -1.0 + 2.0 * G4UniformRand();
	  const double phi = twopi * G4UniformRand();
	  vertex.SetDirection( TVector3( sin( acos( cosTheta ) ) * cos( phi ),
									 sin( acos( cosTheta ) ) * sin( phi ),
									 cosTheta ) );
	  vertex.SetDirectionErrors( TVector3( 1.0, 1.0, 1.0 ) );
	}
  if( vertex.ContainsEnergy() )
	{
	  // Random numbers
	  vertex.SetEnergy( 6.0 * G4UniformRand() );
	  vertex.SetEnergyErrors( 12.0 );
	}
  newSeed.SetVertex( 0, vertex );
  return newSeed;
}