#include "CalorimeterSimulation.h"
#include "CalorimeterConstants.h"
#include "CalorimeterGeometry.h"
#include "CellAddress.h"

#include "TRandom.h"

#include <cstdlib>
#include <cmath>

CalorimeterSimulation::CalorimeterSimulation()
{
  m_NbShoot = 5000;
}

void CalorimeterSimulation::CalorimeterData(std::vector<CalCell>& caldata) const
{
  for ( std::map<CellAddress,float>::const_iterator it = m_caldata.begin(); 
	it != m_caldata.end(); ++it )
    caldata.push_back(CalCell(it->first,it->second));
}

void CalorimeterSimulation::SimulateShower(float xstart, float ystart, float energy)
{
  using namespace std;
  // constants
  static const double TWOPI = 2.*acos(-1.0);
  static const float X0 = 0.01;  // radiation length
  static const float RM = 0.05;  // Moliere radius
  static const float ShowerSigma = 0.03;

  // first check that the point in within calorimeter volume
  if ( xstart < CalConst::XYMin || xstart > CalConst::XYMax ) return;
  if ( ystart < CalConst::XYMin || ystart > CalConst::XYMax ) return;

  float zstart = CalConst::ZMin+0.001;


  int Nshoot = 0;
  while ( Nshoot < m_NbShoot )
  {
    // generate in local coordinate system at starting point
    float x1=0.,y1=0.,z1=0.;
    bool accept = false;
    while ( !accept ) 
    {
      static const float gammamax = gamma(-1.);
      z1 = gRandom->Uniform(CalConst::ZMax-CalConst::ZMin);
      float height = gRandom->Uniform(0.,gammamax);
      if ( height < gamma(z1/X0) ) accept = true;
    }

    // generate transverse shape 
    float r1 = gRandom->Gaus(0.,ShowerSigma);
    float phi1 = gRandom->Uniform(0.,TWOPI);
    x1 = r1*cos(phi1);
    y1 = r1*sin(phi1);

    // translate to starting point
    double xyz[3];
    xyz[0] = x1 + xstart;
    xyz[1] = y1 + ystart;
    xyz[2] = z1 + zstart;



    CellAddress ca;
    if ( CalorimeterGeometry::IsInside(xyz,ca) ) 
    {
      // add energy to shower
      std::map<CellAddress,float>::iterator where = m_caldata.find(ca);
      if ( where != m_caldata.end() )
      {
	where->second += energy/static_cast<float>(m_NbShoot);
      }
      else
	m_caldata[ca] =  energy/static_cast<float>(m_NbShoot);

      Nshoot++;
    }
  }
}
  

float CalorimeterSimulation::gamma(float t) const
{
  // see Particle Data Book
  static const float b = 0.5; // for electrons 
  static const float a = 4.;   // for a shower max at 6 (X0 units)
  if ( t < 0. ) return gamma((a-1.f)/b);  // maximum
  return exp((a-1.)*log(b*t))*exp(-b*t);
}

std::ostream& operator<<(std::ostream& os, const CalorimeterSimulation &cs)
{
  os << "Calorimeter Cells: " << std::endl;
  for ( std::map<CellAddress,float>::const_iterator it = cs.m_caldata.begin();
	it != cs.m_caldata.end(); it++ ){
    os << CalCell(it->first,it->second) << std::endl;
  }
  return os;
}