source: PSPA/Interface_Web/trunk/pspaWT/sources/controler/src/softwareMadx.cc @ 495

#include "softwareMadx.h"
#include "dataManager.h"

#include "mathematicalConstants.h"
#include "PhysicalConstants.h"

#include <boost/filesystem.hpp>

softwareMadx::softwareMadx() : abstractSoftware()
{
  nameOfSoftware_ = new nomDeLogiciel("madx");
}

// softwareMadx::softwareMadx(string inputFileName,sectionToExecute* sect, dataManager* data) : abstractSoftware(inputFileName,sect,data)
// {
//   nameOfSoftware_ = new nomDeLogiciel("madx");
// }

softwareMadx::softwareMadx(string inputFileName,computingBlock* cmpb, dataManager* data) : abstractSoftware(inputFileName,cmpb,data)
{
  nameOfSoftware_ = new nomDeLogiciel("madx");
}

bool softwareMadx::createInputFile(particleBeam* beamBefore,string workingDir)
{
  cout << "***********************************" << endl;
  cout << " softwareMadx::createInputFile(...) " << endl << endl;
  dataManager_->consoleMessage(" softwareMadx::createInputFile");

  string name = workingDir + inputFileName_;
  ofstream outfile;
  outfile.open(name.c_str(),ios::out);
  if(!outfile) {
    dataManager_->consoleMessage(" softwareMadx::createInputFile : error opening output stream ");
    cerr << " error opening output stream " << name << endl;
    return false;
  }

  //  sector* sector= getSectionToExecute()->getSector();
  //  cout << " softwareMadx::createInputFile sector " << sector->getName() << endl;
  ///////////////////////////////////////////
  
  unsigned firstIndex= 0;
  string sbeam;
  abstractElement* elPtr;
  //  elPtr = getSectionToExecute()->getElements().front();
  elPtr = getComputingBlock()->getFirstElement();
  nomdElements::typedElement eType = elPtr->getNomdElement().getElementType();

  if(eType == nomdElements::RFgun) {
    // le 1er elt est RFgun
    firstIndex = 1;
    dataManager_->consoleMessage(" softwareMadx::createInputFile : set from RF cavity the quantities to be supplied to the madx BEAM command");
    sbeam = RFgunData(elPtr->parametersToSoftware());
  } else if(eType == nomdElements::beam) {
    // le 1er elt est beam
    firstIndex = 1;
    dataManager_->consoleMessage(" softwareMadx::createInputFile : set from beam the quantities to be supplied to the madx BEAM command");
    sbeam = beamData(elPtr->parametersToSoftware());
  } else {
    // on suppose qu'il y a déjà un faisceau
    if(beamBefore == NULL) { 
      // il n'y a pas de faisceau : erreur
      dataManager_->consoleMessage(" softwareMadx::createInputFile : no input beam, input file not created");
      return false;
    } else { 
      // il y a un faisceau : on le met au format madx
      dataManager_->consoleMessage(" softwareMadx::createInputFile : extract from input beam the quantities to be supplied to the madx BEAM command");
      
      if (beamBefore->particleRepresentationOk()) {
        // le faisceau est représenté à la "parmela"
        cout << "passsage de la représentation macroparticules à la representation moments" << endl;
        beamBefore->buildMomentRepresentation();
      } else if(beamBefore->momentRepresentationOk()) {
        // le faisceau est représenté à la "transport"
        cout << "le faisceau est représenté par moments" << endl;
      } else {
        // représentation ni "macroparticules" ni "moments"
        dataManager_->consoleMessage(" softwareMadx::createInputFile : that's unclear: the beam is represented neither by macroparticles neither by moments");
        return false;
      }
      sbeam = beamData(beamBefore);
    }
  }
  ///////////////////////////////////////////

  // element label //////////////////////////
  outfile << endl; // saut de ligne

  ostringstream os;
  os << "L: line=(";
  //  unsigned nElts= getSectionToExecute()->getElements().size();
  unsigned nElts= getComputingBlock()->getNumberOfElements();
  for(unsigned k = firstIndex; k < nElts; k++)
    {
      //      elPtr = getSectionToExecute()->getElements()[k];
      elPtr = getComputingBlock()->getElement(k);
      cout << " debug:: element [" << k << "] " << elPtr->getNomdElement().getExpandedName() << endl;
      vector<statements> v= elPtr->parametersToSoftware();
      outfile << inputFormat(v);

      if(k < nElts-1) 
        os << elPtr->getLabel() << ",";
      else
        os << elPtr->getLabel() << ");" << endl;
    }

  // lines/sublines ///////////////////////////////
  outfile << endl; // saut de ligne
  
  // relection and repetition ///////////////
  //  os << "all: " << "line=(" << sector->getRepetitionNumber() << "*L);";
  int bid = 1;
os << "all: " << "line=(" << bid << "*L);";
  ///////////////////////////////////////////
  
  outfile << os.str() << endl;
  outfile << endl; // saut de ligne
  ///////////////////////////////////////////

  outfile << sbeam << endl; // beam commnd p46
  outfile << "use, period = all;" << endl << endl; // action commands p45
  
  outfile << "select,flag = twiss,column = name,s,betx,bety;" << endl; //p48
  outfile << "twiss,save,centre,file = "+workingDir+"twiss.out;" << endl; //p51
  outfile << "stop;" << endl;

  outfile.close();
  dataManager_->consoleMessage("input file done for madx");
  return true;
}

string softwareMadx::beamData(const vector<statements>& v) const 
{
  // x  (cm)   = v.at(1).second.at(0)
  // xp (mrad) = v.at(1).second.at(1)
  // y  (cm)   = v.at(1).second.at(2)
  // yp (mrad) = v.at(1).second.at(3)
  // dl  (cm)   = v.at(2).second.at(0)
  // del (mrad) = v.at(2).second.at(1)
  // p0 (GeV/c) = v.at(3).second.at(0)
  
  ostringstream os;
  os << "beam, PARTICLE = ELECTRON";

  double PC = atof(v.at(3).second.at(0).c_str());
  os << ", ENERGY = " << PC; // momentum of the central trajectory [GeV/c]

  double x  = atof(v.at(1).second.at(0).c_str());
  double xp = atof(v.at(1).second.at(1).c_str());
  double y  = atof(v.at(1).second.at(2).c_str());
  double yp = atof(v.at(1).second.at(3).c_str());
  double dl = atof(v.at(2).second.at(0).c_str());
  double del= atof(v.at(2).second.at(1).c_str());
  beam2Moments mts(x,xp,y,yp,dl,del);
  os << emittances(mts) << ";";
  return os.str();  
}

string softwareMadx::beamData(particleBeam* beam) const 
{
  // les commandes BEAM de madx sont :
  // PARTICLE (defaut POSITRON)
  // ENERGY (defaut 1 GeV) ou bien PC GeV/c
  // EX et EY (defaut 1 m.rad)
  // ET (defaut 1 m)

  ostringstream os;
  os << "beam, PARTICLE = ELECTRON";

  double PC = beam->getP0Transport();
  os << ", ENERGY = " << PC; // en GeV/c

  const beam2Moments& mts = beam->getTransportMoments();
  os << emittances(mts) << ";";
  return os.str();
}

string softwareMadx::emittances(const beam2Moments& mts) const 
{
  const vector< vector<double> > rij= mts.getMoments();
  ostringstream os;

  // emittance (x,x')
  double r10= rij.at(1).at(0); // r10 et s10 = s01= r10*r00*r11
  double rac= 0.0;
  if(r10*r10 < 1.0) rac= sqrt(1.0-r10*r10);
  double EX = rij.at(0).at(0)*rij.at(1).at(1)*rac; // r00*r11*sqrt(1-s01*s10)
  // r00 en cm (= 1E-02 m) et r11 en mrad (= 1E-03 rad)
  os << ", EX = " << EX*1.E-05;

  // emittance (y,y')
  double r32= rij.at(3).at(2); // r32 et s32 = s23= r32*r22*r33
  rac= 0.0;
  if(r32*r32 < 1.0) rac= sqrt(1.0-r32*r32);
  double EY = rij.at(2).at(2)*rij.at(3).at(3)*rac; // r22*r33*sqrt(1-s23*s32)
  os << ", EY = " << EY*1.E-05;

  return os.str();
}

string softwareMadx::RFgunData(const vector<statements>& v) const 
{
  //NPART = v.at(1).second.at(0);
  //SIGT (m) = v.at(2).second.at(0) (ps)
  //?? = sigma_r = v.at(2).second.at(1); en cm
  //EX (m.rad) = v.at(3).second.at(0) (pi.mm.mrad)
  //EY (m.rad) = v.at(3).second.at(1) (pi.mm.mrad) 
  //ENERGY (GeV) = v.at(4).second.at(0) (MeV)
  //SIGE (GeV) = v.at(4).second.at(1) (MeV)
  //CHARGE = v.at(6).second.at(0);

  ostringstream os;
  os << "beam, PARTICLE = ELECTRON";

  // masse au repos (en GeV) E0 = 1.E-03*EREST_MeV
  // energie cinetique (en GeV) T = 1.E-03*E_cin
  // l'energie totale  (en Gev) W = E0+T
  double W = 1.E-03*(EREST_MeV + atof(v.at(4).second.at(0).c_str()));
  os << ", ENERGY = " << W; // total energy in [Gev]

  // pi*EX = 1E-06*emit_x en pi.m.rad
  double EX = 1.E-06*atof(v.at(3).second.at(0).c_str())/PI;
  os << ", EX = " << EX; // horizontal emittance in [rad.m]
  double EY = 1.E-06*atof(v.at(3).second.at(1).c_str())/PI;
  os << ", EY = " << EY; // vertical emittance in [rad.m]

  // SIGT = c*sigma_t*1E-12
  double SIGT = 1.E-04*CLIGHT_E8*atof(v.at(2).second.at(0).c_str());
  os << ", SIGT = " << SIGT; // bunch length in [m]

  // SIGE = sigma_Ecin/(p0*c) ou p0 impulsion de la paricule de réf. 
  // os << ", SIGE = " << atof(v.at(4).second.at(1).c_str());

  os << ";";
  return os.str();
}

string softwareMadx::inputFormat(const vector<statements>& v) const 
{
  // v.at(0).first = "labelsGenericSpecific"
  // v.at(0).second = vector<string> de 2 elements (type de l'element,label)
  string keyword= v.at(0).second.at(0);
  string label= v.at(0).second.at(1);
  ostringstream os;
  if(keyword == "drift") {
    double length = atof(v.at(1).second.at(0).c_str());
    os << label << ":" << " drift, l=" << 0.01*length << ";" << endl;

  } else if(keyword == "mpole") {
    double x[9]= {0.0}; // 0 <= n <= 9 in user's manual of MAD program  
    int n= atoi(v.at(1).second.at(0).c_str());
    x[n]= atof(v.at(1).second.at(1).c_str());
    os << label << ":" << " multipole, knl={" << x[0];
    for(int i = 1; i <= n; i++) {
      os << "," << x[i];
    }
    os <<"};" << endl;
    
  } else if(keyword == "qpole") {
    double ln = atof(v.at(1).second.at(0).c_str());
    double k1 = atof(v.at(1).second.at(1).c_str());
    os << label << ":" << " quadrupole, l=" << ln << ", k1= " << k1 << ";" << endl;

  } else if(keyword == "spole") {
    double ln = atof(v.at(1).second.at(0).c_str());
    double k2 = atof(v.at(1).second.at(1).c_str());
    os << label << ":" << " sextupole, l=" << ln << ", k2= " << k2 << ";" << endl;

  } else {
    cout << "softwareMadx::inputFormat ERROR : element type= " << keyword << " not defined" << endl;
  }
  
  return os.str();
}

bool softwareMadx::execute(string workingDir)
{
  cout << "***********************************" << endl;
  cout << " softwareMadx::execute(...) " << endl << endl;

  dataManager_->consoleMessage(" softwareMadx::execute");
  ostringstream sortie;
  bool status= true;
  string mjob = workingDir + "madx64";
  
  if(!boost::filesystem::exists(mjob.c_str())) {  
    sortie << "Error:: " << mjob << " does not exist\n";
    status = false;
  } else {
    mjob += string(" <  ");
    mjob += workingDir + inputFileName_;
    sortie << " run " << mjob << endl;
    string resultOfRun;
    bool success = launchJob(mjob,resultOfRun);
    // xx sortie << resultOfRun << endl;
    if ( !success ) {
      //sortie << " launching of madx failed " << endl;
      status = false;
    } else {
      sortie << " madx finished normally" << endl;
      string nameOut = workingDir + "madx.output";
      ofstream outfile;
      outfile.open(nameOut.c_str(),ios::out);
      if (!outfile) {
        sortie << " error in opening madx output stream " << nameOut << endl;
        status = false;
      } else {
        // on copie la sortie dans un fichier 'madx.out'
        outfile << resultOfRun << endl;
        outfile.close();
      }
    }
  }

  dataManager_->consoleMessage(sortie.str());
  return status;
}

bool softwareMadx::buildBeamAfterElements(string workingDir)
{
  cout << "***********************************" << endl;
  cout << " softwareMadx::buildBeamAfterElements(...) " << endl << endl;

  dataManager_->consoleMessage(" softwareMadx::buildBeamAfterElements: not programmed :O((");
  return false;
}