#ifndef PARTICLEBEAM_SEEN
#define  PARTICLEBEAM_SEEN
#include <iostream>
#include <sstream>
#include <stdio.h>
#include<cmath>
#include <vector>

#include "bareParticle.h"
#include "nomdElements.h"
#include "beam2Moments.h"

using namespace std;

class particleBeam 
{
  bool particleRepresentationOk_;
  bool momentRepresentationOk_;
  beam2Moments rij_;
  double P0Transport_;

  // la representation du faisceau de particules sera a preciser ulterieurement 
  // Pour l'instant une particule courante est definie par : 
  // x
  // y
  // deltaz = z - z0 (z0 est le z de la particule de reference)
  // gamma.betax
  // gamma.betay
  // gamma.betaz
  vector<bareParticle> goodPartic_;

  // par ailleurs, la particules de reference est definie par x,y,z,gamma.betax, gamma.betay, gamma.betaz
  // elle est dupliquee dans le faisceau prorement dit (goodPartic_) avec z = z-z0 et ne doit, par consequent,
  // pas etre prise en compte dans les calculs de trajectoire.
  bareParticle referenceParticle_;
  

  vector<double> centroid_;

  void impressionDesMoments() const;
  void razDesMoments();
  void particlesPhaseSpaceComponent(vector<double>& coord, unsigned index);
  unsigned indexFromPspaToTransport(unsigned index) const;

  void histogramInitialize(unsigned int index,vector<double>& vshf,double out[3]);
  void histogramPacking(double ecatyp,vector<double>vshf,vector<double>&xcor,vector<int>& hist);

  double getSigmaTransportij(unsigned indexI, unsigned indexJ);
double getUnnormalizedTranspPhaseSpaceArea(unsigned TranspIndexAbs, unsigned TranspIndexOrd);


 unsigned pspaCoorIndexFromString(string nameCoor) const;
 unsigned transportCoorIndexFromString(string nameCoor) const;

inline double dimensionalFactorFromTransportToGraphics(unsigned index) {
  if ( index == 0 || index == 2 || index == 4 ) return 10.;   // mm
  else return  1.0;  // mrad
}




 inline string transportVariableName(unsigned index) {
  switch ( index ) {
  case 0 : return  " x ";
  case 1 : return  " x' ";
  case 2 : return " y ";
  case 3 : return " y' ";
  case 4 : return " l ";
  case 5 : return "dp/p";
  default : { 
    cout << " particleBeam::transportVariableName : coordinate index ERROR index :  "<< index << endl;
    return "";
  }
  }
 }




 inline string graphicTransportUnitName(unsigned index) {
  switch ( index ) {
  case 0 : return  " mm ";
  case 1 : return  " mrad ";
  case 2 : return " mm ";
  case 3 : return " mrad ";
  case 4 : return " mm ";
  case 5 : return " % ";
  default : { 
    cout << " particleBeam::transportVariableName : coordinate index ERROR index :  "<< index << endl;
    return "";
  }
  }
 }

 public:

  particleBeam(); 
  virtual ~particleBeam() {;}

  void buildMomentRepresentation();
  void clear();
  int getNbParticles() const;
  const beam2Moments&  getTransportMoments() const;
  double getUnnormalizedEmittanceX();
  double getP0Transport() const;
  double referenceKineticEnergyMeV() const;
  bool particleRepresentationOk() const;
  bool momentRepresentationOk() const;
  void  addParticle( bareParticle p);
  const vector<bareParticle>& getParticleVector() const;
  vector<bareParticle>& getParticleVector(); 
  double getXmaxRms();
  void getVariance(double& varx, double& vary, double& varz) const;
  void set2Moments(beam2Moments& moments);
  void setWithParticles(vector<double>& centroid, bareParticle& referencePart, vector<bareParticle>& particles);
  void printAllXYZ() const;
  void Zrange(double& zmin, double& zmax) const;
  void donneesDessinEllipse(vector<double>& xcor, vector<double>& ycor, vector<string>& legende, string namex, string namey);
  void particlesPhaseSpaceData(vector<double>& xcor, vector<double>& ycor, vector<string>& legende, string namex, string namey);
  virtual string FileOutputFlow() const;
  virtual bool FileInput(ifstream& ifs);

  void histogramme(unsigned int index,vector<double>&xcor,vector<int>& hist,double out[3]);
};
#endif