#include <string>
#include <iostream>
#include <iomanip>

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

#include "JTools/JQuantile.hh"

#include "JPhysics/JRadiation.hh"
#include "JPhysics/JRadiationFunction.hh"
#include "JPhysics/JRadiationSource.hh"
#include "JPhysics/JGeane.hh"
#include "JPhysics/JConstants.hh"

#include "JSirene/JSeaWater.hh"
#include "JLang/JSharedPointer.hh"

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


namespace {

  static const std::string LENGTH = "length";  //!< calculation of interaction length
  static const std::string ENERGY = "energy";  //!< calculation of shower energy
  static const std::string RANGE  = "range";   //!< calculation of muon range
  static const std::string ELOSS  = "eloss";   //!< calculation of b parameter in dE/dx formula

  /**
   * Clone master histogram.
   *
   * \param  h1       pointer to master histogram
   * \param  buffer   name of clone
   * \return          pointer to cloned histogram
   */
  inline TH1* clone(TH1* h1, const char* buffer)
  {
    if (h1 != NULL) 
      return (TH1*) h1->Clone(buffer);
    else
      return NULL;
  }
}


/**
 * \file
 *
 * Example program to histogram radiation cross sections, shower energy, range and b(E)
 * as a function of the muon energy.
 * \author mdejong
 */
int main(int argc, char* argv[])
{
  using namespace std;

  string  outputFile;
  int     numberOfPoints;
  string  option;
  int     debug;

  try {

    JParser<> zap("Example program to histogram radiation cross sections, shower energy, range and b(E).");

    zap['o'] = make_field(outputFile)      = "radiation.root";
    zap['n'] = make_field(numberOfPoints)  = 0;
    zap['O'] = make_field(option)          = LENGTH, ENERGY, RANGE, ELOSS;
    zap['d'] = make_field(debug)           = 0;
    
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  using namespace JPP;


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

  TH1* h0 = NULL;                                                               // master histogram

  if (option == LENGTH)  {
    h0 = new TH1D("total", NULL, 160, 2.0, 10.0);                               // interaction length [m]
  }

  if (option == ENERGY) {
    h0 = new TH1D("total", NULL,  24, 2.0, 10.0);                               // <E> [GeV]
  }

  NOTICE("Setting up radiation tables... " << flush);

  typedef pair<JSharedPointer<JRadiationInterface>, TH1*> tuple;
  typedef vector<tuple>                                   ntuple;

  ntuple radiation;

  const JRadiation hydrogen ( 1.0,  1.0, 40, 0.01, 0.1, 0.1);
  const JRadiation oxygen   ( 8.0, 16.0, 40, 0.01, 0.1, 0.1);
  const JRadiation chlorine (17.0, 35.0, 40, 0.01, 0.1, 0.1);
  const JRadiation sodium   (11.0, 23.0, 40, 0.01, 0.1, 0.1);

  JSharedPointer<JRadiation> Hydrogen (new JRadiationFunction(hydrogen, 300, 0.2, 1.0e11));
  JSharedPointer<JRadiation> Oxygen   (new JRadiationFunction(oxygen,   300, 0.2, 1.0e11));
  JSharedPointer<JRadiation> Chlorine (new JRadiationFunction(chlorine, 300, 0.2, 1.0e11));
  JSharedPointer<JRadiation> Sodium   (new JRadiationFunction(sodium,   300, 0.2, 1.0e11));

  radiation.push_back(tuple(new JRadiationSource(11, Oxygen,   DENSITY_SEA_WATER * JSeaWater::O(),  EErad_t), clone(h0, "[eerad O]" )));
  radiation.push_back(tuple(new JRadiationSource(12, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), EErad_t), clone(h0, "[eerad Cl]")));
  radiation.push_back(tuple(new JRadiationSource(13, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(),  EErad_t), clone(h0, "[eerad H]" )));
  radiation.push_back(tuple(new JRadiationSource(14, Sodium,   DENSITY_SEA_WATER * JSeaWater::Na(), EErad_t), clone(h0, "[eerad Na]" )));
  
  radiation.push_back(tuple(new JRadiationSource(21, Oxygen,   DENSITY_SEA_WATER * JSeaWater::O(),  Brems_t), clone(h0, "[Brems O]" )));
  radiation.push_back(tuple(new JRadiationSource(22, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), Brems_t), clone(h0, "[Brems Cl]")));
  radiation.push_back(tuple(new JRadiationSource(23, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(),  Brems_t), clone(h0, "[Brems H]" )));
  radiation.push_back(tuple(new JRadiationSource(24, Sodium,   DENSITY_SEA_WATER * JSeaWater::Na(), Brems_t), clone(h0, "[Brems Na]" )));

  radiation.push_back(tuple(new JRadiationSource(31, Oxygen,   DENSITY_SEA_WATER * JSeaWater::O(),  GNrad_t), clone(h0, "[gnrad O]" )));
  radiation.push_back(tuple(new JRadiationSource(32, Chlorine, DENSITY_SEA_WATER * JSeaWater::Cl(), GNrad_t), clone(h0, "[gnrad Cl]")));
  radiation.push_back(tuple(new JRadiationSource(33, Hydrogen, DENSITY_SEA_WATER * JSeaWater::H(),  GNrad_t), clone(h0, "[gnrad H]" )));
  radiation.push_back(tuple(new JRadiationSource(34, Sodium,   DENSITY_SEA_WATER * JSeaWater::Na(), GNrad_t), clone(h0, "[gnrad Na]" )));

  NOTICE("OK" << endl);

  
  if (option == LENGTH) {

    for (int i = 1; i <= h0->GetNbinsX(); ++i) {
    
      const double x = h0->GetBinCenter(i);
      const double E = pow(10.0, x); 
    
      STATUS("Energy: " << SCIENTIFIC(8,1) << E << '\r'); DEBUG(endl);

      for (ntuple::iterator p = radiation.begin(); p != radiation.end(); ++p) {

	const double li = p->first->getInverseInteractionLength(E);

	p->second->Fill(x, 1.0/li);
      }
    }
    STATUS(endl);
  }


  if (option == ENERGY) {

    for (int i = 1; i <= h0->GetNbinsX(); ++i) {
    
      const double x = h0->GetBinCenter(i);
      const double E = pow(10.0, x); 

      STATUS("Energy: " << SCIENTIFIC(8,1) << E << '\r'); DEBUG(endl);

      for (ntuple::iterator p = radiation.begin(); p != radiation.end(); ++p) {

	JQuantile Q;

	for (int n = 0; n != numberOfPoints; ++n) {
	  Q.put(p->first->getEnergyOfShower(E));
	}
	
	p->second->SetBinContent(i, Q.getMean());
	//p->second->SetBinError  (i, Q.getSTDev());
      }
    }
    STATUS(endl);
  }


  if (option == RANGE) {
    
    TH1D* Ra = new TH1D("R[analytical]", NULL, 12, 2.0, 8.0);                                            // Range [m]
    TH1D* Rb = new TH1D("R[numerical]",  NULL, 12, 2.0, 8.0);                                            // Range [m]
    
    for (int i = 1; i <= Ra->GetNbinsX(); ++i) {

      const double x  = Ra->GetBinCenter(i);
      const double E0 = pow(10.0, x); 
      const double Z  = gWater(E0);
      
      Ra->SetBinContent(i, Z*1e-3);
    }
    
    for (int j = 1; j <= Rb->GetNbinsX(); ++j) {

      const double x  = Rb->GetBinCenter(j);
      const double E0 = pow(10.0, x); 
    
      STATUS("Energy: " << SCIENTIFIC(8,1) << E0 << '\r'); DEBUG(endl);

      JQuantile Q;

      for (int n = 0; n != numberOfPoints; ++n) {
      
	double E = E0;
	double z = 0.0;
      
	while (E >= 0.5) {
	
	  const int N = radiation.size();
	
	  double li[N];                                       // inverse interaction lengths
	  double ls = 1.0e-5;                                 // minimal total inverse interaction length (100 km)^-1
	
	  for (int i = 0; i != N; ++i) {
	    ls += li[i] = radiation[i].first->getInverseInteractionLength(E);
	  }
	
	  double dz = min(gRandom->Exp(1.0) / ls, gWater(E));

	  E -= gWater.getA() * dz;
	
	  double Es = 0.0;                                            // shower energy [GeV]
	  double y  = gRandom->Uniform(ls);
	
	  for (int i = 0; i != N; ++i) {
	  
	    y -= li[i];
	  
	    if (y < 0.0) {
	      Es = radiation[i].first->getEnergyOfShower(E);          // shower energy [GeV]
	      break;
	    }
	  }
	
	  z += dz;
	  E -= Es;
	}

	Q.put(z);
      }

      Rb->SetBinContent(j, Q.getMean() * 1e-3);
      //Rb->SetBinError  (j, Q.getSTDev()  * 1e-3);
    }
    STATUS(endl);
  }
  
    
  if (option == ELOSS) {

    TH1D* hb = new TH1D("hb", NULL, 12, 2.0, 8.0);                                            // b(E)

    for (int j = 1; j <= hb->GetNbinsX(); ++j) {

      const double x  = hb->GetBinCenter(j);
      const double E  = pow(10.0, x); 

      STATUS("Energy: " << SCIENTIFIC(8,1) << E << '\r'); DEBUG(endl);

      JQuantile Q;

      const int N = radiation.size();
    
      double li[N];                                       // inverse interaction lengths
      double ls = 1.0e-5;                                 // minimal total inverse interaction length (100 km)^-1
    
      for (int i = 0; i != N; ++i) {
	ls += li[i] = radiation[i].first->getInverseInteractionLength(E);
      }
	
      for (int i = 0; i != numberOfPoints; ++i) {
      
	double Es = 0.0;                                            // shower energy [GeV]
	double y  = gRandom->Uniform(ls);
	
	for (int i = 0; i != N; ++i) {
	
	  y -= li[i];
	
	  if (y < 0.0) {
	    Es = radiation[i].first->getEnergyOfShower(E);          // shower energy [GeV]
	    break;
	  }
	}

	Q.put(1e3*ls*Es/E);
      }

      hb->SetBinContent(j, Q.getMean());
      //hb->SetBinError  (j, Q.getSTDev());
    }
    STATUS(endl);
  }

  delete h0;

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