source: PSPA/Interface_Web/trunk/pspaWT/src/particleBeam.cc @ 224

#include "particleBeam.h"
#include "mathematicalConstants.h"
#include "PhysicalConstants.h"
//#include <string>
#include <stdio.h>
#include <algorithm>

#include "environmentVariables.h"



particleBeam::particleBeam()  {
  // rij_transportMoments_.resize(6);
  // unsigned dim=0;
  // unsigned k;
  // for ( k=0; k < 6; k++){
  //   rij_transportMoments_.at(k).resize(++dim); 
  // }
  P0Transport_ = 0.0;
  particleRepresentationOk_ = false;
  momentRepresentationOk_ = false;
}

void particleBeam::clear() {
  goodPartic_.clear();
  rij_.raz();
  P0Transport_ = 0.0;
  particleRepresentationOk_ = false;
  momentRepresentationOk_ = false;
}

int particleBeam::getNbParticles() const {
  return goodPartic_.size();
}

const transportMoments&   particleBeam::getTransportMoments() const  { 
    return rij_;
}

double particleBeam::getSigmaTransportij(unsigned i, unsigned j)  {
  if ( i < 1 || i > 6 || j < 1 || j > 6 ) {
    cerr << " particleBeam::getSigmaTransportij() indices out of range  " << endl;
    return 0.0;
  }
  if ( !momentRepresentationOk_ ) {
    cerr << " particleBeam::getSigmaTransportij() beam is not in moment representation " << endl;
    return 0.0;
  }

  i--;
  j--;
  if ( j > i ) {
    unsigned aux = i;
    i = j;
    j = aux;
  }
  return ( rij_.getMatrix().at(i) ).at(j);
}


double particleBeam::getUnnormalizedEmittanceX() {
  double r = getSigmaTransportij(2,1);
  double rac = (1 - r*r);
  if ( rac <= 0.0 ) {
    return 0.0;
  }
  rac = sqrt(1 - r*r);
  return  10.*getSigmaTransportij(1,1) * getSigmaTransportij(2,2) * rac; // en pi.mm.mrad
}

double particleBeam::getP0Transport() const { 
  return P0Transport_;
}
bool particleBeam::particleRepresentationOk() const {
  return particleRepresentationOk_;
}
bool particleBeam::momentRepresentationOk() const {
  return momentRepresentationOk_;
}

void  particleBeam::addParticle( bareParticle p)
{
  goodPartic_.push_back(p);
}

const vector<bareParticle>& particleBeam::getParticleVector() const
{
  return goodPartic_;
}
vector<bareParticle>& particleBeam::getParticleVector() 
{
  return goodPartic_;
}


bool  particleBeam::setFromTransport(string elLab, const nomdElements elem)
{
  string elementLabel = elLab;
  // transformer le label en majuscules ; je ne suis pas sur que ca 
  // marchera a tous les coups (glm)
  std::transform(elementLabel.begin(), elementLabel.end(), elementLabel.begin(), (int (*)(int))std::toupper);

  cout << " particleBeam::setFromTransport on cherche element " << elementLabel << endl;
  string buf;
  ifstream infile;
  string nameIn = WORKINGAREA + "transport.output";
  infile.open(nameIn.c_str(), ios::in);
  if (!infile) {
    cerr << " particleBeam::setFromTransport : error re-opening transport output stream " << nameIn << endl;
    return false;
  }
  //  else cout << " particleBeam::setFromTransport() : ouverture du fichier " << nameIn << endl; 

  string::size_type nn = string::npos;
  while ( getline(infile, buf) ) {
    //      cout << " buf= "  << buf << endl;
    nn = buf.find(elementLabel);
    //     cout << " string::npos= " << string::npos << " nn= " << nn << endl;
    if ( nn != string::npos ) {
      //          cout << " particleBeam::setFromTransport : element " << elementLabel << " trouve " << endl;
      break;
    }
  }

  if ( nn == string::npos ) {
    cerr << " particleBeam::setFromTransport : element " << elementLabel << " non trouve dans le fichier  " << nameIn << endl;
    return false;
  }
  if (elem.getElementType() == bend ) {
    getline(infile, buf);
    getline(infile, buf);
  }
  readTransportMoments(infile);
  //  impressionDesMoments();
  infile.close();
  momentRepresentationOk_ = true;
  return true;
}


bool particleBeam::setFromParmela(unsigned numeroElement, double referencefrequency) {
  unsigned  k;
  FILE* filefais;
  string nomfilefais = WORKINGAREA + "parmdesz";
  cout << " nom fichier desz : " << nomfilefais << endl;
  filefais = fopen(nomfilefais.c_str(), "r");

  if ( filefais == (FILE*)0 ) {
    cerr << " particleBeam::setFromParmela() erreur a l'ouverture du fichier" << nomfilefais  << endl;; 
    return false;
  }
  else cout << " particleBeam::setFromParmela() : ouverture du fichier " << nomfilefais << endl; 

  parmelaParticle partic;
  std::vector<parmelaParticle> faisceau;

  cout << " particleBeam::setFromParmela : numeroElement = " << numeroElement << endl;
  unsigned indexElement = numeroElement-1;
 



  int testNombrePartRef =0;
  double phaseRef;

  while( partic.readFromParmelaFile(filefais) > 0 ) {
    if ( partic.ne == indexElement )
      {
        faisceau.push_back(partic);

        if ( partic.np == 1 ) {
          // en principe on est sur la particule de reference
          if ( fabs(partic.xx) > EPSILON || fabs(partic.yy) > EPSILON || fabs(partic.xxp) > EPSILON  || fabs(partic.yyp) > EPSILON) {
            printf(" ATTENTION part. reference douteuse  \n");
            partic.imprim();
          }
          phaseRef = partic.phi;
          TRIDVECTOR  posRef(partic.xx,partic.yy,0.0);
          TRIDVECTOR betagammaRef(partic.xxp*partic.begamz, partic.yyp*partic.begamz, partic.begamz);
          referenceParticle_ = bareParticle(posRef, betagammaRef);
          testNombrePartRef++;
          if ( testNombrePartRef != 1 ) {
            cerr << " TROP DE PART. DE REF : " << testNombrePartRef << " !! " << endl;
            return false;
          }
        }
      }
  }

  if ( faisceau.size() == 0) 
    {
      cerr << " particleBeam::setFromParmela echec lecture  element " << numeroElement << endl;
      return false;
    }
  
  // facteur  c/ 360. pour calculer (c dphi) / (360.freq) 
  // avec freq en Mhz et dphi en degres et résultat en cm: 
  double FACTEUR =  83.3333;  // ameliorer la precision



  // pour l'instant on choisit un centroid nul;
  centroid_ = vector<double>(6,0.0);

  goodPartic_.clear();
  goodPartic_.resize(faisceau.size(), bareParticle());
  double x,xp,y,yp;
  double betagammaz;
  double betaz;
  double deltaz;
  double dephas;
  double g;
  TRIDVECTOR  pos;
  TRIDVECTOR betagamma;
  // contrairement a ce qu'indique la notice PARMELA, dans parmdesz, les xp et yp 
  // sont donnes en radians
  for ( k=0; k < faisceau.size(); k++) {
    x=faisceau.at(k).xx;
    xp=faisceau.at(k).xxp;
    y=faisceau.at(k).yy;
    yp=faisceau.at(k).yyp;

    // dephasage par rapport a la reference  
    dephas = faisceau.at(k).phi - phaseRef; // degrés
    g = faisceau.at(k).wz/ERESTMeV;
    betagammaz = faisceau.at(k).begamz;
    betaz = betagammaz/(g+1.0);
    deltaz = FACTEUR * betaz * dephas / referencefrequency;
    x += xp * deltaz;
    y += yp * deltaz;
    pos.setComponents(x,y,deltaz);
    betagamma.setComponents(xp*betagammaz, yp*betagammaz, betagammaz);
    goodPartic_.at(k) = bareParticle(pos,betagamma);
  }
  particleRepresentationOk_ = true;
  return true;
}


void particleBeam::getVariance(double& varx, double& vary, double& varz) const {
  unsigned int k;
  double x,y,z;
  double xav = 0.;
  double yav = 0.;
  double zav = 0.;
  double xavsq = 0.;
  double yavsq = 0.;
  double zavsq = 0.;

  TRIDVECTOR pos;


  for ( k = 0 ; k < goodPartic_.size(); k++) {
    pos = goodPartic_.at(k).getPosition();
    pos.getComponents(x,y,z);
    //      partic_[k].getXYZ(x,y,z);
    xav += x;
    xavsq += x*x;
    yav += y;
    yavsq += y*y;
    zav += z;
    zavsq += z*z;
  }

  double aginv = double (goodPartic_.size());
  aginv = 1.0/aginv;

  varx =  aginv * ( xavsq - xav*xav*aginv );
  vary =  aginv * ( yavsq - yav*yav*aginv );
  varz =  aginv * ( zavsq - zav*zav*aginv );
}


void particleBeam::printAllXYZ() const {
  cout << " dump du faisceau : " << endl;
  cout <<  goodPartic_.size() << " particules " << endl;
  unsigned int k;
  for ( k = 0 ; k < goodPartic_.size(); k++)
    {
      double xx,yy,zz;
      goodPartic_.at(k).getPosition().getComponents(xx,yy,zz);
      cout << " part. numero " << k << "  x= " << xx << " y= " << yy  << " z= " << zz << endl;
    }
}



void particleBeam::Zrange(double& zmin, double& zmax) const {
  double z;
  zmin = GRAND;
  zmax = -zmin;

  unsigned int k;
  for ( k = 0 ; k < goodPartic_.size(); k++)
    {
      z = goodPartic_.at(k).getZ();
      if ( z < zmin ) zmin = z;
      else if ( z > zmax) zmax = z;         
    }
}



string particleBeam::FileOutputFlow() const {
  ostringstream sortie;
  unsigned int k;
  for ( k = 0 ; k < goodPartic_.size(); k++)
    {
      sortie << goodPartic_.at(k).FileOutputFlow() << endl;
    }
  sortie << endl;
  return sortie.str();
}

bool particleBeam::FileInput( ifstream& ifs) {
  bool test = true;
  string dum1, dum2;
  double dummy;
  if ( !( ifs >> dum1 >> dum2 >> dummy) ) return false;
 
  bareParticle pp;
  while ( pp.FileInput(ifs) )
    {
      addParticle( pp);
    }
  return test;
}

void particleBeam::buildMomentRepresentation() {

  unsigned k,j,m;
  double auxj, auxm;
  if ( !particleRepresentationOk_)
    {
      cerr << " particleBeam::buildMomentRepresentation() vecteur de particules invalide" << endl;
      return;
    }

  cout << " buildMomentRepresentation " << endl;
  //  printAllXYZ();

  double gref = referenceParticle_.getGamma() - 1.0;
  double P_reference_MeV_sur_c = sqrt( gref*(gref+2) );

  // initialisation des moments
  razDesMoments();

  // accumulation
  TRIDVECTOR pos;
  TRIDVECTOR begam;
  double gamma;
  double begamz;
  double g;
  double PMeVsc;
  double del;
  vector<double> part(6);

    vector< vector<double> >& matrice = rij_.getMatrix();


  for (k=0; k < goodPartic_.size(); k++) {
    gamma = goodPartic_.at(k).getGamma();
    pos = goodPartic_.at(k).getPosition();
    begam= goodPartic_.at(k).getBetaGamma();
    begamz = begam.getComponent(2);
    g = gamma -1.0;
    PMeVsc = sqrt( g*(g+2) );
    del = 100.0 * ( PMeVsc -  P_reference_MeV_sur_c ) / P_reference_MeV_sur_c ; // en %

    part[0] = pos.getComponent(0);
    part[1] = begam.getComponent(0)/begamz;
    part[2] = pos.getComponent(1);
    part[3] = begam.getComponent(1)/begamz;
    part[4] = pos.getComponent(2);
    part[5] = del;

    for ( j = 0; j < 6; j++) {
      auxj = part.at(j) - centroid_.at(j);
      for (m=0; m <= j; m++) 
        {
          auxm = part.at(m) - centroid_.at(m);

          ( matrice.at(j) ).at(m) += auxj*auxm;
          //      ( rij_transportMoments_.at(j) ).at(m) += auxj*auxm;


          //          cout << " j= " << j << " m= " << m << " rjm= " << ( rij_transportMoments_.at(j) ).at(m) << endl;
        }
    }
  }


  // moyenne
  double facmoy = 1.0/double( goodPartic_.size() );
  for ( j = 0; j < 6; j++) {
        ( matrice.at(j) ).at(j) = sqrt(( matrice.at(j) ).at(j) * facmoy );
  }

  for ( j = 0; j < 6; j++) {
    auxj =  ( matrice.at(j) ).at(j);
    for (m=0; m < j; m++) {
      auxm = ( matrice.at(m) ).at(m);
      (  matrice.at(j) ).at(m) *= facmoy/(auxj * auxm);
    }
  }
    
  // les longueurs sont en cm
  // les angles en radians, on passe en mrad;

  double uniteAngle = 1.0e+3;
  ( matrice.at(1) ).at(1)  *= uniteAngle;
  ( matrice.at(3) ).at(3)  *= uniteAngle;
  P0Transport_ = 1.0e-3*ERESTMeV*P_reference_MeV_sur_c;

  //  cout << " buildmomentrepresentation impression des moments " << endl;
  //  impressionDesMoments();

  momentRepresentationOk_ = true;
}

void particleBeam::impressionDesMoments() const {
  rij_.impression();
}

void particleBeam::razDesMoments() {
  rij_.raz();
}


void particleBeam::readTransportMoments(ifstream& inp) {
rij_.readFromTransportOutput(inp);
}

double particleBeam::getXmaxRms() {
  if ( !momentRepresentationOk_ ) buildMomentRepresentation();
  return ( rij_.getMatrix().at(0) ).at(0);
  // return ( rij_transportMoments_.at(0) ).at(0);
}

void particleBeam::donneesDessinEllipseXxp(vector<double>& xcor, vector<double>& ycor) {
  int k;
  double x,y;

  if ( !momentRepresentationOk_ ) return;

  xcor.clear();
  ycor.clear();

  double xm = ( rij_.getMatrix().at(0) ).at(0);
  double ym = ( rij_.getMatrix().at(1) ).at(1);
  double r  = ( rij_.getMatrix().at(1) ).at(0);

  cout <<  " racs11= " << xm << " racs22= " << ym << " r12= " << r << endl;


  int nbintv = 50;
  if ( xm == 0.0 ) return;
  double pas = 2.0 * xm / nbintv;

  //  cout << " r= " << r << endl;
  double rac = (1 - r*r);
  if ( rac > 0.0 ) 
    {
      cout << " cas rac > " << endl;
      rac = sqrt(1 - r*r);
      double alpha = -r / rac;
      double beta = xm / ( ym * rac);
      //  double gamma = ym / ( xm * rac );
      double epsil = xm * ym * rac;
      double fac1 = -1.0 / ( beta * beta);
      double fac2 = epsil/beta;
      double fac3 = -alpha/beta;
      double aux;
      for ( k=0; k < nbintv; k++)
        {
          x = -xm + k*pas;
          aux = fac1 * x * x + fac2;
          //     cout << " aux2= " << aux << endl;
          if ( aux <= 0.0 )
            {
              aux = 0.0;
            }
          else aux = sqrt(aux);
      
          //        y = fac3*x;
          y = fac3*x + aux;
          xcor.push_back(x);
          ycor.push_back(y);
        }

      for ( k=0; k <= nbintv; k++)
        {
          x = xm - k*pas;
          aux =  fac1 * x * x + fac2;
          if ( aux <= 0.0 ) 
            {
              aux = 0.0;
            }
          else aux = sqrt(aux);
          //   y = fac3*x;
          y = fac3*x - aux;
          xcor.push_back(x);
          ycor.push_back(y);
        }
    }
  else
    // cas degenere
    {
      cout << " cas degenere " << endl;
      double fac = ym/xm;
      for ( k=0; k < nbintv; k++)
        {
          x = -xm + k*pas;
          y = fac*x;
          xcor.push_back(x);
          ycor.push_back(y);
        }
        
    }
}

void particleBeam::histogramme(vector<double>&xcor,vector<int>& hist,int& cnts)
{
  double gammin= GRAND;
  double gammax= -gammin;
  double Emoy= 0.0;
  double ecatyp= 0.0;

  for (unsigned int k = 0; k < goodPartic_.size(); k++)
    {
      double gamma = goodPartic_.at(k).getGamma();
      double EMev = (gamma-1.0)*ERESTMeV;
      if (gamma < gammin) gammin = gamma;
      else if (gamma > gammax) gammax = gamma;
      Emoy += EMev;
      ecatyp += EMev*EMev;
    }

  double sum= (float)goodPartic_.size();
  Emoy /= sum;
  ecatyp /= sum;
  ecatyp = sqrt(abs(ecatyp-Emoy*Emoy));

  double Emin = (gammin-1.0)*ERESTMeV;
  double Emax = (gammax-1.0)*ERESTMeV;
  cout << "energie cinetique -moyenne " << Emoy << " Mev " << "-mini " << Emin << " Mev " << "-maxi " << Emax << " Mev " << "ecart type " << ecatyp*1000.0 << " Kev" << endl;

  vector<double> Eshf;
  for (unsigned int k = 0; k < goodPartic_.size(); k++)
    {
      double gamma = goodPartic_.at(k).getGamma();
      double EMev = (gamma-1.0)*ERESTMeV;
      Eshf.push_back(EMev-Emoy);
    }

  //////////////////////////////////////////////////////////////////////////////
  
  // demi fenetre en energie, et pas de l'histogramme
  double hfene= max(3.*ecatyp-Emoy,Emoy-3.*ecatyp);
  double hpas = hfene/50.;

  cout << "demi fenetre " << hfene << ", hpas= " << hpas << endl;

  double vmin = -hfene;
  double dfen = 2.*hfene;
  int ihist = dfen/hpas;
  double phist = ihist*hpas;
  double dpas = hpas-(dfen-phist);
  if(dpas <= hpas*1.e-03) {
    ihist++;
    phist= ihist*hpas;
  }
  double vmax= vmin+hpas*ihist;
  
  cout << "'xAxisNumberOfBins= " << ihist <<", xAxisMinimum= " << vmin << ", xAxisMaximum= " << vmax << ", NParticules= " << Eshf.size() << ", phist " << phist << endl;

  xcor= vector<double>(ihist);
  for (unsigned i = 0; i < ihist; ++i) {
    xcor[i]= vmin+i*hpas;
  }

  hist= vector<int>(ihist,0);
  for (unsigned i = 0; i < Eshf.size(); ++i) {
    double var= Eshf[i]-vmin;
    if(var < 0 || var >= phist) cout<<"out of range "<<var<<", ("<< i<<")"<< endl;
    int k= var/hpas;
    int kk= (int)floor(var/hpas);
    if(i%20 == 0) cout<<"v("<<i<<")= " <<var<<" ["<<k<<"-"<<kk<<"], "<<endl; 
    hist[kk]++;
  }

  cnts= 0;
  for (unsigned i = 0; i < ihist; ++i) {
    if(hist.at(i) > 0) cnts++;
    cout<<"("<<xcor.at(i)<<","<<hist.at(i)<<") ";
  }
  cout<< " ... cnts=  " << cnts << endl;
}