source: PSPA/Interface_Web/trunk/pspaWT/sources/controler/src/softwareParmela.cc @ 481

#include "softwareParmela.h"
#include "abstractElement.h"
#include "parmelaParticle.h"
#include "mathematicalConstants.h"
#include "PhysicalConstants.h"
#include "dataManager.h"
#include "mixedTools.h"

softwareParmela::softwareParmela() : abstractSoftware()
{
  nameOfSoftware_ = nomDeLogiciel("parmela");
}

softwareParmela::softwareParmela(string inputFileName,sectionToExecute* sect, dataManager* data) : abstractSoftware(inputFileName, sect, data)
{
  nameOfSoftware_ = nomDeLogiciel("parmela");
  registerElement(nomdElements::RFgun,TBoolOk);
  registerElement(nomdElements::drift,TBoolOk);
  registerElement(nomdElements::cell,TBoolOk);
  registerElement(nomdElements::bend,TBoolOk);
  registerElement(nomdElements::soleno,TBoolOk);
  registerElement(nomdElements::fit,TBoolIgnore);
  registerElement(nomdElements::snapshot,TBoolIgnore);
}



bool softwareParmela::createInputFile(particleBeam* beamBefore, string workingDir)
{
  unsigned int k;

//  setRelativeParmelaElementIndices();
  string name = workingDir + inputFileName_;
  ofstream outfile;
  outfile.open(name.c_str(), ios::out);
  if (!outfile) {
    dataManager_->consoleMessage(" softwareParmela::createInputFile : error opening output stream " );
    cerr << " softwareParmela::createInputFile : error opening output stream " << name << endl;
    return false;
  }
  
  abstractElement* elPtr = getSectionToExecute()->getElements().front();
  double initalKineticEnergy = 0.0;
  bool there_is_rfgun = (elPtr->getNomdElement().getElementType() == nomdElements::RFgun );
  
  if ( !there_is_rfgun ) {
    string nameOfInput = workingDir + "parin.input0";
    if ( !beamToParmela(nameOfInput,beamBefore) ) return false;
    initalKineticEnergy = beamBefore->referenceKineticEnergyMeV();
    // les elements de parmela sont indexes de 1 Ã  max, s'il n'y a pas de rfgun
    //    offsetNumElem_ = numeroDeb_ -1;
  } else {
    elPtr->setPhaseStep( dataManager_->getGlobalParameters()->getIntegrationStep() );
    initalKineticEnergy = elPtr->getInitialKineticEnergy();
    // les elements de parmela sont indexes de 0 Ã  max, s'il y a un rfgun
  }
  
  outfile << "TITLE" << endl;
  outfile << " titre provisoire " << endl;
  outfile << "RUN /n0=1 /ip=999 /freq=" << dataManager_->getGlobalParameters()->getFrequency() << "  /z0=0.0 /W0=" << initalKineticEnergy << "  /itype=1" << endl;
  outfile << "OUTPUT 0" << endl;
  if ( there_is_rfgun ) {
    outfile << elementsData(elPtr->parametersToSoftware());
  } else {
    outfile << "INPUT 0 /NP=" << beamBefore->getNbParticles() << endl;
  }
  
 // retrouver le sector !!
  for ( k =1; k < getSectionToExecute()->getElements().size(); k++)
    {
      outfile << elementsData(getSectionToExecute()->getElements()[k]->parametersToSoftware());
    }
    
  outfile << "ZOUT" << endl;
  outfile << "START /wt=0.0 /dwt=" << dataManager_->getGlobalParameters()->getIntegrationStep() << "  /nsteps=" << dataManager_->getGlobalParameters()->getNbSteps() << "  nsc=" << dataManager_->getGlobalParameters()->getScPeriod() << "  /nout=10" << endl;
  outfile << "END" << endl;
  outfile.close();
  dataManager_->consoleMessage("fichier input termine pour PARMELA");
  return true;
}

bool softwareParmela::execute(string workingDir)
{
  bool ExecuteStatus = true;
  
  ostringstream sortie;
  sortie << " EXECUTION DE PARMELA " <<  endl;
    
  string parmelaJob = workingDir + "parmela";
  parmelaJob += string("   ");
  parmelaJob += workingDir;
  
  string resultOfRun;
  bool success = launchJob(parmelaJob,resultOfRun);


  //  sortie << resultOfRun << endl;
  if ( !success ) {
    sortie << " launching of parmela failed " << endl;
    ExecuteStatus = false;
  } else {
    sortie << " successful launching of parmela " << endl;
    cout << " execution parmela MARCHE " << endl;
    string::size_type nn = (resultOfRun).find("normal");
    if ( nn == string::npos )
      {
        sortie << " abnormal exit of parmela " << endl;
        ExecuteStatus = false;
      } 
  }
  
  dataManager_->consoleMessage(sortie.str());
  return ExecuteStatus;
}

bool softwareParmela::buildBeamAfterElements( string workingDir)
{
  bool result = true;

  unsigned curseur;
  for ( unsigned int k=0; k < getSectionToExecute()->getElements().size() ; k++ ) {
    abstractElement* elem = getSectionToExecute()->getElements()[k];
      if ( elem == NULL ) {
        dataManager_->consoleMessage(" softwareParmela::buildBeamAfterElements : null pointer on element " );
        return  false;        
      }

    curseur = k;

    if (!(doAcceptElement(getSectionToExecute()->getElements()[curseur]->getNomdElement().getElementType()) == TBoolOk)) {
//      if ( relativeParmelaElementIndices_.at(curseur) < 0 ) {

        // si l'element doit etre ignore de parmela, on renvoie sur le diag precedent
        particleBeam* lastDiag = dataManager_->updateCurrentDiagnostic(false);

        // if(elem->getNomdElement().getElementType() == snapshot) {
        //   // si cet element est un snapshot, on organise la sortie correspondante
        //   string* param = elem->getParametersString();
        //   string cliche = workingDir + param[2].c_str();
        //   if( beamToParmela(cliche,lastDiag) ) {
        //       //         cout  <<  "["  <<  k  <<  "] : ecrit sur fichier " << cliche << " le contenu de beam["<<avantDernier<<"]"<<endl;
        //     cout  <<  "["  <<  k  <<  "] : ecrit sur fichier " << cliche << " le contenu de beam[ ]"<<endl;
        //   }
        // }
        // si le numero est reconnu de parmela
      } else {
         
        // on initialise une nouvelle sortie diagnostic 
          particleBeam* newDiag = dataManager_->updateCurrentDiagnostic(true);
          vector<double> centroid;
          bareParticle refPart;
          vector<bareParticle> particles;
          double freq= dataManager_->getGlobalParameters()->getFrequency();
          unsigned numeroParmel;
//        numeroParmel = (unsigned)relativeParmelaElementIndices_.at(curseur);
    numeroParmel = curseur;
          cout << " lecture PARMELA el no absolu " << k << " numero relatif " << numeroParmel << " nom " << elem->getNomdElement().getExpandedName() << endl;
        // lecture sortie parmela
        if(!beamFromParmela(workingDir,numeroParmel,freq,centroid,refPart,particles))
          {
            // si echec, fin
            dataManager_->consoleMessage(" softwareParmela::buildBeamAfterElements : failure in reading parmela result for element " + elem->getLabel() + " for unknown reason " );
            return  false;
          } else {
          // si succes, on complete le diagnostic
          newDiag->setWithParticles(centroid,refPart,particles);
        }
      }
    }
  return result;
}
bool softwareParmela::beamFromParmela(string workingDir,unsigned numeroParmel, double referencefrequency, vector<double>& centroid, bareParticle& refPart,vector<bareParticle>& particles )
{    
  string nomfilefais = workingDir + "parmdesz";
  cout << " nom fichier desz : " << nomfilefais << endl;
  FILE *filefais = fopen(nomfilefais.c_str(), "r");
  
  if ( filefais == (FILE*)0 ) {
    dataManager_->consoleMessage(" beamFromParmela() erreur a l'ouverture du fichier 'parmdesz'");
    cerr << " beamFromParmela() erreur a l'ouverture du fichier" << nomfilefais  << endl;
    return false;
  } else
    cout << " beamFromParmela() : ouverture du fichier " << nomfilefais << endl;
  
  parmelaParticle partic;
  std::vector<parmelaParticle> faisceau;
  int testNombrePartRef =0;
  double phaseRef = 0.0;
  
  while( partic.readFromParmelaFile(filefais) > 0 )
    {
      
      if ( partic.ne == (int)numeroParmel ) {
        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;
          
          // le 'z' est 'absolu' (le long de la trajectoire)
          TRIDVECTOR  posRef(partic.xx,partic.yy,partic.z);
          TRIDVECTOR betagammaRef(partic.xxp*partic.begamz, partic.yyp*partic.begamz, partic.begamz);
          refPart = bareParticle(posRef, betagammaRef);
          testNombrePartRef++;
          if ( testNombrePartRef != 1 ) {
            dataManager_->consoleMessage(" beamFromParmela : nombre de part. de ref different de 1 :");
            cerr << " nombre de part. de ref different de 1 : " << testNombrePartRef << " !! " << endl;
            return false;
          }
        }
        faisceau.push_back(partic);
      }
    } //while
  
  if ( faisceau.size() == 0)
    {
      string stringNum = mixedTools::intToString( (int)numeroParmel );
      dataManager_->consoleMessage("beamFromParmela echec lecture  element numero relatif  parmela : " + stringNum);
      cerr << " beamFromParmela echec lecture  element numero relatif  parmela " << numeroParmel << 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.clear();
  centroid = vector<double>(6,0.0);
  
  particles.clear();
  particles.resize(faisceau.size(), bareParticle());
  //  double x,xp,y,yp;
  double xp, yp;
  double betagammaz;
  //  double betaz;
  double cdt;
  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 (unsigned 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/EREST_MeV;
    betagammaz = faisceau.at(k).begamz;
    //    betaz = betagammaz/(g+1.0);
    //    deltaz = FACTEUR * betaz * dephas / referencefrequency;
    cdt = FACTEUR * dephas / referencefrequency;
    //    x += xp * deltaz;
    //    y += yp * deltaz;
    pos.setComponents(faisceau.at(k).xx,faisceau.at(k).yy,cdt);
    betagamma.setComponents(xp*betagammaz, yp*betagammaz, betagammaz);
    particles.at(k) = bareParticle(pos,betagamma);
  }
  
  return true;
}

// sauvegarde d'un 'particleBeam' sur un fichier parmela, en guise d'INPUT
// pour l'instant de nom standard 'parin.input0'
bool softwareParmela::beamToParmela(string nameOfFile,particleBeam* beam)
{
    if ( !beam->particleRepresentationOk() ) {
        dataManager_->consoleMessage("softwareParmela::beamToParmela : beam not in particles form : not yet programmed");
        cout << " softwareParmela::beamToParmela : beam not in particles form : not yet programmed " << endl;
        return false;
    }
    
    ofstream outfile;
    outfile.open(nameOfFile.c_str(),ios::out);
    if (!outfile) {
        dataManager_->consoleMessage(" softwareParmela::beamToParmela : error opening output stream ");
        cerr << " softwareParmela::beamToParmela : error opening output stream " << nameOfFile << endl;
        return false;
    }
    
    //    const vector<bareParticle>& partic = beam->getParticleVector();
    double weight = 1.0;
    double xx,yy,zz;
    double begamx, begamy, begamz;
    //    TRIDVECTOR pos, begam;
    double zmin = GRAND;
    double zmax = -zmin; 
    double cdtmin, cdtmax;
    beam->ZrangeCdt(cdtmin, cdtmax);
    cout << " softwareParmela::beamToParmela cdtmin = " << cdtmin << " cdtmax = " << cdtmax << endl;
  
    for (unsigned k = 0; k < beam->getNbParticles(); k++) {
        // partic.at(k).getPosition().getComponents(xx,yy,zz);
        // partic.at(k).getBetaGamma().getComponents(begamx,begamy,begamz);
      beam->coordonneesDeployees(k, -cdtmax).getComponents(xx,yy,zz);
      beam->betaGamma(k).getComponents(begamx,begamy,begamz);
      if ( zz > zmax) zmax = zz;
      if ( zz < zmin ) zmin = zz;
        outfile << xx << " " << begamx << " " <<  yy << " " << begamy << " " << zz << " " << begamz  << " " << weight << endl;
    }
    outfile.close();
    cout << " softwareParmela::beamToParmela zmn = " << zmin << " zmax = " << zmax << endl;
    return true;
}


string softwareParmela::elementsData(const vector< pair<string, vector<string> > >& donnees) const {
  unsigned k;
  cout << " PASSAGE softwareParmela::elementsData " << endl;
  if ( donnees.at(0).first != "labelsGenericSpecific" ) {
    cout << " softwareParmela::elementsData ERROR : element badly defined " << endl;
    return string();
  }
  string genericName = donnees.at(0).second.at(0);
  if ( genericName == "rfgun" ) return rfgunData(donnees);
  if ( genericName == "cell" ) return cellData(donnees);
  if ( genericName == "drift" ) return driftData(donnees);
  if ( genericName == "solnd" ) return solenoData(donnees);
  if ( genericName == "bend" ) return bendData(donnees);
  return string();
}

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

  cout << " PASSAGE softwareParmela::rfgunData " << endl;
  string nmacrop = "0";
  string sigma_t = "0.0";
  string sigma_r = "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 == "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;
  // on prend les troncatures tmax et rmax à 3 sigmas 
  sortie << "INPUT 10 /np=" << nmacrop << "  /sigt=" << sigma_t << endl;
  sortie << "  /tmax=" << 3.3*atof(sigma_t.c_str()) << " /sigr=" << sigma_r << endl;
  sortie << " /rmax=" << 3.0*atof(sigma_r.c_str()) << " /W0=" << E_cin << " /dw0=" << sigma_E << endl;
  sortie << " /dwt=" << phaseStep << " /ran=2" << endl;

  // on doit entrer le nb vrai de part. (avec signe moins)
  sortie << "SCHEFF /beami=" << -fabs( atof(totalCharge.c_str()) )/ELECTRONANOCOULOMB << " /nprog=2 /point=-1.7"  << endl;
 
  return sortie.str();
}

string softwareParmela::cellData(const vector< pair<string, vector<string> > >& donnees) const {
  cout << " PASSAGE softwareParmela::cellData " << endl;

  string lenght = "0.0";
  string aperture = "0.0";
  string initialPhase = "0.0";
  string acceleratingField = "0.0";
  string phaseStepMax = "0.0";
  string acceleratingShapeFile = "";
  string focusingMagFile = "";
  string offsetMag = "0.0";
  string scaleFactor = "0.0";

  unsigned k;
  for ( k=1; k < donnees.size(); k++) {
    if ( donnees.at(k).first == "lenghtAperture" ) { 
      lenght = donnees.at(k).second.at(0);
      aperture = donnees.at(k).second.at(1);
    } else if ( donnees.at(k).first == "phaseInitialStepmax" ) {
      initialPhase = donnees.at(k).second.at(0);
      phaseStepMax = donnees.at(k).second.at(1);
    } else if ( donnees.at(k).first == "fieldValueFile" ) {
      acceleratingField = donnees.at(k).second.at(0);
      acceleratingShapeFile = donnees.at(k).second.at(1);
    } else if ( donnees.at(k).first == "MagFocusingFileOffsetScale") {
      focusingMagFile = donnees.at(k).second.at(0);
      offsetMag = donnees.at(k).second.at(1);
      scaleFactor = donnees.at(k).second.at(2);
    }
  }
  ostringstream sortie;
    sortie << "CELL /l=" << lenght << "  /aper=" << aperture << endl;
    sortie << "  /iout=1  /phi0=" << initialPhase << " /E0=" << acceleratingField << endl;
    sortie << " /nc=1 /dwtmax=" << phaseStepMax << " /sym=-1" << endl;
    sortie << "CFIELD 1" << endl;
    sortie << acceleratingShapeFile << endl;
    if ( focusingMagFile != "") {
      sortie << "POISSON /zoff=" << offsetMag << " /rmult=" << scaleFactor << endl;
      sortie << focusingMagFile << endl;
    }
  return sortie.str();

}

string softwareParmela::driftData(const vector< pair<string, vector<string> > >& donnees) const {
  cout << " PASSAGE softwareParmela::driftData " << endl;

  string lenght = "0.0";
  string aperture = "0.0";

  unsigned k;
  for ( k=1; k < donnees.size(); k++) {
    if ( donnees.at(k).first == "lenghtAperture" ) { 
      lenght = donnees.at(k).second.at(0);
      aperture = donnees.at(k).second.at(1);
    }
  }
  ostringstream sortie;
  sortie << "DRIFT /l=" << lenght << "  /aper=" << aperture << "  /iout=1" << endl;
  return sortie.str();
}

string softwareParmela::solenoData(const vector< pair<string, vector<string> > >& donnees) const {
  cout << " PASSAGE softwareParmela::solenoData " << endl;
  string lenght = "0.0";
  string aperture = "0.0";
  string B0 = "0.0";

  unsigned k;
  for ( k=1; k < donnees.size(); k++) {
    if ( donnees.at(k).first == "lenghtAperture" ) { 
      lenght = donnees.at(k).second.at(0);
      aperture = donnees.at(k).second.at(1);
    } else if ( donnees.at(k).first == "field" ) {
      B0 = donnees.at(k).second.at(0);
    }
  }
    ostringstream sortie;
    // on passe l'induction magnetique en Gauss 
    sortie << "SOLENOID /l=" << lenght << "  /aper=" << aperture << "  /iout=1 /h=" << 1000.*atof(B0.c_str())  << endl;
    return sortie.str();

}

string softwareParmela::bendData(const vector< pair<string, vector<string> > >& donnees) const {
  cout << " PASSAGE softwareParmela::bendData " << endl;

  string lenght = "0.0";
  string aperture = "0.0";
  string angleDeg = "0.0";
  string momentum = "0.0";
  string beta1 = "0.0";
  string beta2 = "0.0";
  unsigned k;
  for ( k=1; k < donnees.size(); k++) {
    if ( donnees.at(k).first == "lenghtAperture" ) { 
      lenght = donnees.at(k).second.at(0);
      aperture = donnees.at(k).second.at(1);
    } else if ( donnees.at(k).first == "angleDegre" ) {
      angleDeg = donnees.at(k).second.at(0);
    } else if ( donnees.at(k).first == "momentum" ) {
      momentum = donnees.at(k).second.at(0);
     } else if ( donnees.at(k).first == "rotatedFaces" ) {
      beta1 = donnees.at(k).second.at(0);
      beta2 = donnees.at(k).second.at(1);
    }   
  }

  double ecin = atof(momentum.c_str())/EREST_MeV;
    ecin = ecin*ecin + 1.;
    ecin = EREST_MeV*(sqrt(ecin) - 1.0);

    ostringstream sortie;
    // il faut entrer l'energie cinetique
    sortie << "BEND /l=" << lenght << "  / aper=" << aperture << "  / iout=1 / wr=" << ecin << " /alpha=" << angleDeg << " / beta1=" << beta1;
    sortie << " / beta2=" << beta2  << endl;

    return sortie.str();

}