source: PSPA/Interface_Web/trunk/pspaWT/sources/controler/src/softwareGenerator.cc @ 434

#include "softwareGenerator.h"
#include "abstractElement.h"
#include "generatorParticle.h"
#include "mathematicalConstants.h"
#include "PhysicalConstants.h"
#include "dataManager.h"

softwareGenerator::softwareGenerator() : abstractSoftware()
{
  nameOfSoftware_ = nomDeLogiciel("generator");
}

softwareGenerator::softwareGenerator(string inputFileName, globalParameters* globals, dataManager* dt) : abstractSoftware(inputFileName, globals, dt)
{
  nameOfSoftware_ = nomDeLogiciel("generator");
  registerElement(nomdElements::RFgun,TBoolOk);
}

bool softwareGenerator::createInputFile(particleBeam* beamBefore, unsigned int numeroDeb, unsigned int numeroFin, string workingDir)
{
  if ( !initComputationLimits(numeroDeb,numeroFin) ) return false;

  if ( numeroDeb_ != numeroFin_ ) return false;

  abstractElement* elPtr;
  elPtr = dataManager_->getElementPointerFromNumero(numeroDeb_);

  if ( elPtr->getNomdElement().getElementType() != nomdElements::RFgun ) {
    dataManager_->consoleMessage(" softwareGenerator::createInputFile : the element must be rfgun " );
    cerr << " softwareGenerator::createInputFile : the element must be rfgun" << endl;
    return false;
  }

  ofstream outfile;
  string name = workingDir + inputFileName_;
  outfile.open(name.c_str(), ios::out);
  if (!outfile) {
    dataManager_->consoleMessage(" softwareGenerator::createInputFile : error opening output stream " );
    cerr << " error opening output stream " << name << endl;
    return false;
  }
  outfile << "&INPUT" << endl;
  string fichier = "'" + workingDir + "faisceau.ini" + "'";
  outfile << "   FNAME = " << fichier << endl;
  outfile << " Add=FALSE, N_add=0" << endl;
  outfile << " Species='electrons'" << endl;
  //  outfile << elPtr->generatorOutputFlow() << endl;

  outfile << elementsData(elPtr->parametersToSoftware()) << endl;

  outfile << "/" << endl;
  outfile.close();
  dataManager_->consoleMessage("fichier input termine pour GENERATOR");
  return true;
}

bool softwareGenerator::execute(string workingDir) {
  
  bool ExecuteStatus = true;
  
  ostringstream sortie;
  sortie << " EXECUTION DE GENERATOR DE l'ELEMENT " << numeroDeb_ << " A L'ELEMENT " << numeroFin_ << endl;

  string generatorJob = workingDir + "generator";
  generatorJob += string("   ");
  generatorJob += workingDir + inputFileName_;

  string resultOfRun;
  bool success = launchJob(generatorJob,resultOfRun);
  sortie << resultOfRun << endl;
  if ( !success ) {
    sortie << " launching of generator failed " << endl;
    ExecuteStatus = false;
  } else {
    cout << " execution generator MARCHE " << endl;
    sortie << resultOfRun;
    string nameOut = workingDir + "generator.output";
    ofstream outfile;
    outfile.open(nameOut.c_str(), ios::out);
    if (!outfile) {
      sortie << " error first opening transport output stream " << nameOut << endl;
      ExecuteStatus = false;
    } else {
      // on copie la sortie dans un fichier 'generator.out'
      outfile << resultOfRun << endl;
      outfile.close();
    }
  }
  
  dataManager_->consoleMessage(sortie.str());
  return ExecuteStatus;
}

bool  softwareGenerator::buildBeamAfterElements( string workingDir) {

  bool result = true;
  if ( !ComputationLimitsOk() ) return false;

  if ( numeroDeb_ != numeroFin_ ) {
    dataManager_->consoleMessage(" softwareGenerator::buildBeamAfterElements : only one element (rfgun) must be calculated " );
    return false;
  }
  // on initialise une nouvelle sortie diagnostic 
  particleBeam* newDiag = dataManager_->updateCurrentDiagnostic(true);

  //  beams.push_back(particleBeam());
  vector<double> centroid;
  bareParticle refPart;
  vector<bareParticle> particles;
  vector<bareParticle> particlesPassives; // on ne fait rien de ces particules pour l'instant
  if (beamFromGenerator(string("faisceau.ini"),workingDir, centroid, refPart,particles, particlesPassives )) {
    newDiag->setWithParticles(centroid, refPart,particles); 
  } else {
    dataManager_->consoleMessage(" softwareGenerator::buildBeamAfterElements : error  " );
    result = false;
  }
  return result;
}

bool softwareGenerator::beamFromGenerator(string beamFileName, string workingDir, vector<double>& centroid, bareParticle& refPart,vector<bareParticle>& particles, vector<bareParticle>& passiveParticles ) {
  unsigned  k;
  FILE* filefais;
  string nomfilefais = workingDir + beamFileName;
  cout << " nom fichier sortie generator : " << nomfilefais << endl;
  filefais = fopen(nomfilefais.c_str(), "r");
  if ( filefais == (FILE*)0 ) {
    dataManager_->consoleMessage(" softwareGenerator::beamFromGenerator : error openig the file " + nomfilefais );
    cerr << " beamFromGenerator() erreur a l'ouverture du fichier" << nomfilefais  << endl;; 
    return false;
  }
  else cout << " beamFromGenerator() : ouverture du fichier " << nomfilefais << endl; 

  generatorParticle partic;
  std::vector<generatorParticle> faisceau;
  int nbProbPart =0;
  double timeRef = 0.0;
  double betagammaZRef = 0.0;
  // lecture part. de reference 
  if ( partic.readFromGeneratorFile(filefais) > 0 )
    {
         if ( partic.index != 1 ) {
    dataManager_->consoleMessage(" softwareGenerator::beamFromGenerator : particles are not electrons, we have to reconsider this method beamFromGenerator" );
           cout << " ERROR softwareGenerator::beamFromGenerator : particles are not electrons, we have to reconsider this method " << endl;
        return false;
      }
   
      if ( partic.flag != -1 ) {
        cout << " ATTENTION softwareGenerator::beamFromGenerator : flag different de -1 " << endl;
      }
      if (fabs(partic.xx) > EPSILON || fabs(partic.yy) > EPSILON || fabs(partic.px) > EPSILON  || fabs(partic.py) > EPSILON) {
        printf(" ATTENTION softwareGenerator::beamFromGenerator : part. reference douteuse  \n");
        partic.imprim();
      }
      timeRef = partic.clock;
      TRIDVECTOR  posRef(100.*partic.xx,100.*partic.yy,100.*partic.zz); // en cm
      betagammaZRef = partic.pz/EREST_eV;
      // l'impulsion donnee par generator est en eV/c
      TRIDVECTOR betagammaRef(partic.px/EREST_eV , partic.py/EREST_eV, betagammaZRef);
      refPart = bareParticle(posRef, betagammaRef);

      // seule la part. de reference a un pz absolu (les autres pat. on un pz relatif a la ref)
      // je mets ici a zero, pour homogeneiser
      partic.pz = 0.0;
      faisceau.push_back(partic);
    }

  while( partic.readFromGeneratorFile(filefais) > 0 ) {
    faisceau.push_back(partic);
    if ( partic.flag != -1 ) nbProbPart++;
  }

  if ( faisceau.size() == 0) 
    {
    dataManager_->consoleMessage(" softwareGenerator::beamFromGenerator : error no particle found" );
      cerr << " softwareGenerator::beamFromGenerator echec lecture " << endl;
      return false;
    }
  // pour l'instant on choisit un centroid nul;
  centroid.clear();
  centroid = vector<double>(6,0.0);

  particles.clear();
  passiveParticles.clear();
  //  particles.resize(faisceau.size() - nbProbPart, bareParticle());
  //  passiveParticles.resize(nbProbPart, bareParticle());
  double x,y;
  //  double deltaz;
  double cdeltat;
  TRIDVECTOR  pos;
  TRIDVECTOR betagamma;
  double pxPart;
  double pyPart;
  double pzPartRel;
  double deltaPzPart;

  for ( k=0; k < faisceau.size(); k++) {

    pxPart = faisceau.at(k).px;
    pyPart = faisceau.at(k).py;
    pzPartRel = faisceau.at(k).pz;
    deltaPzPart = faisceau.at(k).pz;
    x=faisceau.at(k).xx;
    y=faisceau.at(k).yy;

    // tout ce qui suit sera a clarifier
    double betaGammax = pxPart/EREST_eV;
    double betaGammay = pyPart/EREST_eV;
    double betaGammaz = betagammaZRef + pzPartRel/EREST_eV;
    betagamma.setComponents(betaGammax, betaGammay, betaGammaz);
    double gamma = sqrt(1.0 + betagamma.norm2());

    // decalage temporel par rapport a la reference  ?
    // on prend le faisceau sous la forme "a un instant donne"

    //   deltat = faisceau.at(k).clock - timeRef; // nanoseondes
    cdeltat = CLIGHT_m_per_ns * (faisceau.at(k).clock - timeRef); // metres
    //    double ds = CLIGHT_m_per_ns * deltat /gamma;
    //    x += betaGammax * ds;  // en metres
    //    y += betaGammay * ds;

    // ici on neglige la difference entre gamma de la part. et gamma de la ref.
    //    deltaz = (pzPartRel/EREST_eV) * CLIGHT_m_per_ns * deltat; // en metres

    //    pos.setComponents(100.*x,100.*y,100.*deltaz);  // en cm
    pos.setComponents(100.*x,100.*y,100.*cdeltat);  // en cm
    if ( faisceau.at(k).flag == -1 ) {
      particles.push_back(bareParticle(pos,betagamma));
      //      cout << " generator enregistre " << particles.back().FileOutputFlow() << endl;
    } else {
      passiveParticles.push_back(bareParticle(pos,betagamma));
    }
  }


  return true;
}


string softwareGenerator::elementsData(const vector< pair<string, vector<string> > >& donnees) const {
  unsigned k;
  cout << " PASSAGE softwareGenerator::elementsData " << endl;
  if ( donnees.at(0).first != "labelsGenericSpecific" ) {
    cout << " softwareGenerator::elementsData ERROR : element badly defined " << endl;
    return string();
  }
  string genericName = donnees.at(0).second.at(0);
  if ( genericName == "rfgun" ) return rfgunData(donnees);
  return string();
}

string softwareGenerator::rfgunData(const vector< pair<string, vector<string> > >& donnees) const {

  cout << " PASSAGE softwareGenerator::rfgunData " << endl;
  string nmacrop = "0";
  string sigma_t = "0.0";
  string sigma_r = "0.0";
  string emit_x = "0.0";
  string emit_y = "0.0";
  string E_cin = "0.0";
  string sigma_E = "0.0";
  string phaseStep = "0.0";
  string totalCharge = "0.0";

  unsigned k;
  for ( k=1; k < donnees.size(); k++) {
    if ( donnees.at(k).first == "nbMacroparticles" ) { 
      nmacrop = donnees.at(k).second.at(0);
    } else if ( donnees.at(k).first == "sigmasTR" ) {
      sigma_t = donnees.at(k).second.at(0);
      sigma_r = donnees.at(k).second.at(1);
    } else if ( donnees.at(k).first == "emittancesXY" ) {
      emit_x = donnees.at(k).second.at(0);
      emit_y = donnees.at(k).second.at(1);
    } else if ( donnees.at(k).first == "kineticE" ) {
      E_cin = donnees.at(k).second.at(0);
      sigma_E = donnees.at(k).second.at(1);
    } else if ( donnees.at(k).first == "phaseStep" ) {
      phaseStep = donnees.at(k).second.at(0);
    } else if ( donnees.at(k).first == "totalCharge" ) {
      totalCharge = donnees.at(k).second.at(0);
    }
  }
  ostringstream sortie;


    sortie << "Ipart=" << nmacrop << endl;

    sortie << "Probe=.True." << endl;
    sortie << "Noise_reduc=.T" << endl;
    sortie << "Cathode=.T." << endl;
    sortie << "Q_total=" << totalCharge << endl;
    sortie << "Ref_zpos=0.0" << endl;
    sortie << "Ref_clock=0.0" << endl;
    sortie << "Ref_Ekin=" << E_cin << endl; // tjs en MeV
    sortie << "Dist_z='g'" << endl;
    sortie << "sig_clock=" << 1.0e-3*atof(sigma_t.c_str()) << endl; // passage de ps en ns
    sortie <<  "Dist_pz='g', sig_Ekin=" << 1000.*atof(sigma_E.c_str()) << ",  emit_z=0.0 ,    cor_Ekin=0.0 " << endl;  // passage en keV
    sortie << "Dist_x='gauss',  sig_x=" <<  10.*atof(sigma_r.c_str()) << endl;    // passage en mm
    sortie << "Dist_px='g', Nemit_x=" << emit_x << ",   cor_px=0.0" << endl;
    sortie << "Dist_y='gauss',  sig_y=" <<  10.*atof(sigma_r.c_str()) << endl;    // passage en mm
    sortie << "Dist_py='g', Nemit_y=" << emit_y << ",   cor_py=0.0" << endl;
 
  return sortie.str();
}