#include <string>
#include <iostream>
#include <cmath>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TRandom.h"

#include "JPhysics/JNeutrino.hh"
#include "JPhysics/JConstants.hh"

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

namespace {

  const char nc_enabled        = 'A';    //!< Neutral-current interactions as simulated.
  const char nc_no_energy_loss = 'C';    //!< Neutral-current interactions with Bjorken-Y set to zero.
  const char nc_no_scattering  = 'D';    //!< Neutral-current interactions with transverse energy set to zero.
  const char nc_disabled       = 'E';    //!< Neutral-current interactions which have no effect.
  const char nc_absorption     = 'F';    //!< Neutral-current interactions lead to absorption.
}


/**
 * \file
 *
 * Mickey-mouse program to propagate neutrinos through the Earth.
 * \author mdejong
 */
int main(int argc, char* argv[])
{
  using namespace std;
  using namespace JPP;

  typedef long long int  counter_type;

  string        outputFile;
  double        E_GeV;
  counter_type  numberOfEvents;
  UInt_t        seed;
  char          option;
  int           debug;

  try {

    JParser<> zap("Mickey-mouse program to propagate neutrinos through the Earth.");

    zap['o'] = make_field(outputFile,         "output file")                    = "earth.root";
    zap['E'] = make_field(E_GeV,              "neutrino energy at source");
    zap['n'] = make_field(numberOfEvents,     "number of generated events");
    zap['S'] = make_field(seed)              = 0;
    zap['O'] = make_field(option, endl
			  << nc_enabled        << " -> " << "Neutral-current interactions as simulated."                       << endl
			  << nc_no_energy_loss << " -> " << "Neutral-current interactions with Bjorken-Y set to zero."         << endl
			  << nc_no_scattering  << " -> " << "Neutral-current interactions with transverse energy set to zero." << endl
			  << nc_disabled       << " -> " << "Neutral-current interactions which have no effect."               << endl
			  << nc_absorption     << " -> " << "Neutral-current interactions lead to absorption."                 << endl) =
      nc_enabled,
      nc_no_energy_loss,
      nc_no_scattering,
      nc_disabled,
      nc_absorption;

    zap['d'] = make_field(debug)             = 0;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  gRandom->SetSeed(seed);


  const double R_SOURCE_M         =  40000.0;                 // [m]
  const double R_TARGET_M         =    500.0;                 // [m]
  const double EMIN_GEV           =     10.0;                 // [GeV]


  const double Zmin = 0.0;                                    // [m]
  const double Zmax = R_EARTH_KM * 2.0e3;                     // [m]

  enum {
    CC = 0,                                                   // charged-current interaction
    NC,                                                       // neutral-current interaction
    N                                                         // number of interactions
  };

  const JCrossSection* sigma[N] = { NULL };                   // cross-sections

  counter_type number_of_events[2][2] = {{ 0 }};              // [source inside/outside][target inside/outside]

  TFile out(outputFile.c_str(), "recreate");

  TH2D hs("hs", NULL, 
	  101, -R_SOURCE_M, +R_SOURCE_M,
	  101, -R_SOURCE_M, +R_SOURCE_M);

  TH2D ht("ht", NULL, 
	  101, -R_SOURCE_M, +R_SOURCE_M,
	  101, -R_SOURCE_M, +R_SOURCE_M);

  TH1D hn("hn", NULL, 11, -0.5, 10.5);

  TH1D ha("ha", NULL, 100, 0.0, 1.0);
  TH1D hb("hb", NULL, 100, 0.0, 1.0);

  TH2D h2("h2", NULL, 
	  101, -0.01, +0.01,
	  101, -0.01, +0.01);


  for (counter_type count = 0; count != numberOfEvents; ++count) {

    if (count%10000 == 0) {
      STATUS("event " << setw(9) << count << '\r'); DEBUG(endl);
    }

    const bool neutrino = (count%2 == 0);                     // [anti-]neutrino ratio 1:1

    if (neutrino) {
      sigma[CC] = &cc_nu;
      sigma[NC] = &nc_nu;
    } else {
      sigma[CC] = &cc_nubar;
      sigma[NC] = &nc_nubar;
    }

    const double r2  = gRandom->Uniform(0.0, R_SOURCE_M * R_SOURCE_M);
    const double r1  = sqrt(r2);
    const double phi = gRandom->Uniform(-PI, +PI);

    // start values

    double x  = r1 * cos(phi);
    double y  = r1 * sin(phi);
    double z  = Zmin;
    double Tx = 0.0;
    double Ty = 0.0;
    double E  = E_GeV;

    hs.Fill(x,y);

    const int i1 = (sqrt(x*x + y*y) <= R_TARGET_M ? 0 : 1);   // inside -> 0; outside -> 1

    int ni[N] = { 0 };                                        // number of interactions

    while (ni[CC] == 0 && E > EMIN_GEV) {

      double li[N];                                           // inverse interaction length [m^-1]
      double ls = 0.0;

      for (int i = 0; i != N; ++i) {
	ls += li[i] = 1.0e+2 * (*sigma[i])(E) * DENSITY_EARTH * AVOGADRO;
      }

      double dz = gRandom->Exp(1.0) / ls;

      if (dz > Zmax - z) {
	dz = Zmax - z;
      }

      // step

      if (option != nc_disabled && option != nc_no_scattering) {
	x += dz * Tx;
	y += dz * Ty;
      }
      z += dz;

      if (z >= Zmax) {

	break;                                                // ready

      } else {

	double y = gRandom->Uniform(ls);

	if ((y-= li[CC]) < 0.0) {                             // charged-current interaction

	  ++ni[CC];

	  break;                                              // stop
	}

	if ((y-= li[NC]) < 0.0) {                             // neutral-current interaction

	  ++ni[NC];

	  if (option == nc_absorption) {
	    break;
	  }
	  
	  // simplified neutrino scattering

	  double s   = E * (MASS_PROTON + MASS_NEUTRON);
	  double ET  = 0.5 * sqrt(s) * gRandom->Uniform(0.0, 1.0);
	  double phi = gRandom->Uniform(-PI, +PI);

	  double By  = gRandom->Uniform(0.0, 1.0);

	  if (!neutrino) {
	    while (gRandom->Rndm() > (1.0 - By) * (1.0 - By)) {
	      By  = gRandom->Uniform(0.0, 1.0);
	    }
	  }

	  if (option == nc_disabled || option == nc_no_scattering) {
	    ET  = 0.0;
	    phi = 0.0;
	  }

	  if (option == nc_disabled || option == nc_no_energy_loss) {
	    By  = 0.0;
	  }
	  
	  if (neutrino)
	    ha.Fill(By);
	  else
	    hb.Fill(By);

	  Tx = ET * cos(phi) / E;
	  Ty = ET * sin(phi) / E;
	  E  = (1.0 - By) * E;
	}
      }
    }

    if (z >= Zmax) {                                          // reach target
      hn.Fill((double) ni[NC]);
      ht.Fill(x,y);
      h2.Fill(Tx,Ty);
    }

    if (z >= Zmax && sqrt(x*x + y*y) <= R_TARGET_M) {         // hit target

      number_of_events[i1][0] += 1;                           // count events

      if (ni[NC] == 0) {
	number_of_events[i1][1] += 1;                         // count events
      }
    }
  }
  STATUS(endl);

  // summary

  const double n00 = number_of_events[0][0];                  // number of events in target from source area covered by target
  const double n10 = number_of_events[1][0];                  // number of events in target from source area not covered by target
  const double n01 = number_of_events[0][1];                  // number of events in target from source area covered by target, without neutral-current interaction 

  if (option == nc_enabled) {
    cout << " " << SCIENTIFIC(7,1) << E_GeV << flush; 
  }
  {
    cout << " & " << FIXED(7,0) << n00 << flush;
  }
  if (option != nc_disabled && option != nc_no_scattering) {
    cout << " & " << FIXED(7,0) << n10 << flush;
  }
  if (option == nc_enabled) {
    cout << " & " << FIXED(7,0) << n01 << " & " << FIXED(5,2) << (n00 + n10) / n01 << flush;
  }
  if (option == nc_disabled) {
    cout << "  \\\\" << endl;
  }

  if (debug >= status_t) {

    cout << ' ' << setw(8) << right << "inside"  << ' ' << setw(8) << right << "outside" << endl;

    for (int row = 0; row != N; ++row) {
      for (int col = 0; col != N; ++col) {
	cout << ' ' << setw(8) << number_of_events[row][col];
      }
      cout << endl;
    }
  }

  out.Write();
  out.Close();
}