#pragma once

#include "rec/util.hh"
#include "rec/JppShowerPdf.hh"

const uint test_shift( const uint precision = 10 ){

  uint error_count = 0;
  
  for( double E = 1e3; E <= 1e8; E *= 2 ){
    
    for( double D = 0.5; D <= 1000; D *= 2 ){
      
      for( double cd = -0.99; cd < 1.0; cd += 0.01 ){

	double test_cd = get_vertex_cd( get_shower_max_D( D, cd, E ), get_shower_max_cd( D, cd, E ), E );
	double test_D  = get_vertex_D ( get_shower_max_D( D, cd, E ), get_shower_max_cd( D, cd, E ), E );
				
	if( not are_emerged_equal( test_cd, cd, precision ) ){
	  cout << std::setprecision( precision + 5 ) << "cosine angles do not match, original cd: " << cd << ", new cd: " << test_cd << endl;
	  cout << "diff: " << test_cd - cd << endl;
	  cout << "E: " << E << ", D: " << D << ", cd: " << cd << endl;
	  error_count++;
	}
	
	if( not are_emerged_equal( test_D, D, precision ) ){
	  cout << std::setprecision( precision + 5 ) << "distances do not match, original D: " << D << ", new D: " << test_D << endl;
	  cout << "diff: " << test_D - D << endl;
	  cout<< "E: " << E << ", D: " << D << ", cd: " << cd << endl;
	  error_count++;
	}
      }
    }
  }
  return error_count;
};


vector<double> time_pdfs( JppShowerPdf& pdf, const double E = 1e6 )
{
  JppShowerNpe npef( pdf );
  
  vector<double> call_times  = {0.0, 0.0, 0.0};
  vector<double> call_sigmas = {0.0, 0.0, 0.0};
  
  uint call_count = 0;
  
  for( double D = 0.5; D < 1000; D *= 2 ){
    
    for( double cd = -0.99; cd < 1.0; cd += 0.01 ){
      
      for( double theta = 0.0; theta < pi; theta += 0.1 ){
	
	for( double phi = 0.0; phi < pi; phi += 0.1 ){

	  //cout << "D: " << D << ", cd: " << cd << ", theta: " << theta << ", phi: " << phi << endl; 
	  double timer0 = std::clock();
	  pdf.npe_tot<0>( 0, E, D, cd, theta, phi );
	  double timer1 = std::clock();
	  double timer2 = std::clock();
	  npef.npe_tot<0>( 0, E, D, cd, theta, phi );
	  double timer3 = std::clock();

	  double time_pdf  = 1e9 * ( timer1 - timer0 ) / (double) CLOCKS_PER_SEC;
	  double time_npef = 1e9 * ( timer3 - timer2 ) / (double) CLOCKS_PER_SEC;
	  
	  call_times[0] += time_pdf;
	  call_times[2] += time_npef;
	  
	  call_sigmas[0] += pow( time_pdf, 2);
	  call_sigmas[2] += pow( time_npef, 2);
	  
	  call_count++;
	}
      }
    }
  }

  for( uint i = 0; i != 3; i++ ){
    call_times[i]  /= call_count;
    call_sigmas[i] /= call_count;
    call_sigmas[i] -= pow( call_times[i], 2);
    call_sigmas[i] = sqrt( call_sigmas[i] );
  }
  
  cout << "Average call times for" << endl;
  cout << "pdf: " << call_times[0] << " +/- " << call_sigmas[0] << " ns" << endl;
  cout << "npe: " << call_times[2] << " +/- " << call_sigmas[2] << " ns" << endl;
  
  return call_times;
 
}