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

Last change on this file since 398 was 398, checked in by lemeur, 12 years ago

fichier parmout dans workingarea
File size: 12.0 KB
Line 
1#include "softwareParmela.h"
2#include "abstractElement.h"
3#include "parmelaParticle.h"
4#include "mathematicalConstants.h"
5#include "PhysicalConstants.h"
6#include "dataManager.h"
7#include "mixedTools.h"
8//#include "nomDeLogiciel.h"
9softwareParmela::softwareParmela() : abstractSoftware()
10{
11  nameOfSoftware_ = nomDeLogiciel("parmela");
12}
13
14softwareParmela::softwareParmela(string inputFileName, globalParameters* globals, dataManager* dt) : abstractSoftware(inputFileName, globals, dt)
15{
16  nameOfSoftware_ = nomDeLogiciel("parmela");
17}
18
19void softwareParmela::setRelativeParmelaElementIndices() {
20  relativeParmelaElementIndices_.clear();
21  relativeParmelaElementIndices_.resize(numeroFin_ - numeroDeb_ + 1, -1);
22  cout << " setRelativeParmelaElementIndices() taille a priori : " << relativeParmelaElementIndices_.size() << endl;
23  abstractElement* elPtr = dataManager_->getElementPointerFromNumero(numeroDeb_);
24  bool there_is_rfgun = ( elPtr->getNomdElement().getElementType() == RFgun );
25  unsigned offsetNumElem;
26  // les elements de parmela sont indexes de 1 à max, s'il n'y a pas de rfgun
27  if ( !there_is_rfgun ) {
28    offsetNumElem = numeroDeb_ -1;
29    // les elements de parmela sont indexes de 0 à max, s'il y a un rfgun
30  } else {
31    offsetNumElem = numeroDeb_;
32  }
33
34  // index du premier element de parmela
35  int id= numeroDeb_ - offsetNumElem;
36  unsigned k;
37  unsigned curseur = 0;
38  for ( k=numeroDeb_; k <= numeroFin_ ; k++ ) {
39    abstractElement* elem = dataManager_->getElementPointerFromNumero(k);
40    cout << " liste PARMELA no absolu " << k << " relatif provisoire " << relativeParmelaElementIndices_.at(curseur) << endl;
41    if ( elem->is_accepted_by_software(nameOfSoftware_) == TBoolOk ) {
42      relativeParmelaElementIndices_.at(curseur) = id;
43      cout << " mis a " << id << endl;
44      id++;
45    }
46    curseur++;
47  }
48}
49
50
51bool softwareParmela::createInputFile(particleBeam* beamBefore,unsigned int numeroDeb,unsigned int numeroFin,string workingDir)
52{
53  unsigned int k;
54  if ( !initComputationLimits(numeroDeb,numeroFin) ) return false;
55  setRelativeParmelaElementIndices();
56  string name = workingDir + inputFileName_;
57  ofstream outfile;
58  outfile.open(name.c_str(), ios::out);
59  if (!outfile) {
60    dataManager_->consoleMessage(" softwareParmela::createInputFile : error opening output stream " );
61    cerr << " softwareParmela::createInputFile : error opening output stream " << name << endl;
62    return false;
63  }
64 
65  abstractElement* elPtr;
66  double initalKineticEnergy = 0.0;
67  elPtr = dataManager_->getElementPointerFromNumero(numeroDeb_);
68  bool there_is_rfgun = ( elPtr->getNomdElement().getElementType() == RFgun );
69 
70  if ( !there_is_rfgun ) {
71    string nameOfInput = workingDir + "parin.input0";
72    if ( !beamToParmela(nameOfInput,beamBefore) ) return false;
73    initalKineticEnergy = beamBefore->referenceKineticEnergyMeV();
74    // les elements de parmela sont indexes de 1 à max, s'il n'y a pas de rfgun
75    //    offsetNumElem_ = numeroDeb_ -1;
76  } else {
77    elPtr->setPhaseStep( globParamPtr_->getIntegrationStep() );
78    initalKineticEnergy = elPtr->getInitialKineticEnergy();
79    // les elements de parmela sont indexes de 0 à max, s'il y a un rfgun
80    //    offsetNumElem_ = numeroDeb_;
81  }
82 
83  outfile << "TITLE" << endl;
84  outfile << " titre provisoire " << endl;
85  outfile << "RUN /n0=1 /ip=999 /freq=" << globParamPtr_->getFrequency() << "  /z0=0.0 /W0=" << initalKineticEnergy << "  /itype=1" << endl;
86  outfile << "OUTPUT 0" << endl;
87  unsigned int premier = numeroDeb_ ;
88  if ( there_is_rfgun ) {
89    outfile << dataManager_->getElementPointerFromNumero(numeroDeb_)->parmelaOutputFlow();
90    premier++;
91  } else {
92    outfile << "INPUT 0 /NP=" << beamBefore->getNbParticles() << endl;
93  }
94 
95  for ( k = premier; k <= numeroFin_; k++)
96    {
97      elPtr = dataManager_->getElementPointerFromNumero(k);
98      if (elPtr) {
99        outfile << elPtr->parmelaOutputFlow();
100      }
101    }
102   
103  outfile << "ZOUT" << endl;
104  outfile << "START /wt=0.0 /dwt=" << globParamPtr_->getIntegrationStep() << "  /nsteps=" << globParamPtr_->getNbSteps() << "  nsc=" << globParamPtr_->getScPeriod() << "  /nout=10" << endl;
105  outfile << "END" << endl;
106  outfile.close();
107  dataManager_->consoleMessage("fichier input termine pour PARMELA");
108  return true;
109}
110
111bool softwareParmela::execute(string workingDir)
112{
113  bool ExecuteStatus = true;
114 
115  ostringstream sortie;
116  sortie << " EXECUTION DE PARMELA DE l'ELEMENT " << numeroDeb_ << " A L'ELEMENT " << numeroFin_ << endl;
117   
118  string parmelaJob = workingDir + "parmela";
119  parmelaJob += string("   ");
120  parmelaJob += workingDir;
121 
122  string resultOfRun;
123  bool success = launchJob(parmelaJob,resultOfRun);
124
125
126  //  sortie << resultOfRun << endl;
127  if ( !success ) {
128    sortie << " launching of parmela failed " << endl;
129    ExecuteStatus = false;
130  } else {
131    sortie << " successful launching of parmela " << endl;
132    cout << " execution parmela MARCHE " << endl;
133    string::size_type nn = (resultOfRun).find("normal");
134    if ( nn == string::npos )
135      {
136        sortie << " abnormal exit of parmela " << endl;
137        ExecuteStatus = false;
138      } 
139  }
140 
141  dataManager_->consoleMessage(sortie.str());
142  return ExecuteStatus;
143}
144
145bool softwareParmela::buildBeamAfterElements( string workingDir)
146{
147  bool result = true;
148
149  if ( !ComputationLimitsOk() ) return false;
150  unsigned curseur;
151  for(unsigned k = numeroDeb_; k <= numeroFin_; k++)
152    {
153      abstractElement* elem = dataManager_->getElementPointerFromNumero(k);
154      if ( elem == NULL ) {
155        dataManager_->consoleMessage(" softwareParmela::buildBeamAfterElements : null pointer on element " );
156        return  false;       
157      }
158
159      curseur = k - numeroDeb_;
160
161      if ( relativeParmelaElementIndices_.at(curseur) < 0 ) {
162
163        // si l'element doit etre ignore de parmela, on renvoie sur le diag precedent
164        particleBeam* lastDiag = dataManager_->updateCurrentDiagnostic(false);
165
166        // if(elem->getNomdElement().getElementType() == snapshot) {
167        //   // si cet element est un snapshot, on organise la sortie correspondante
168        //   string* param = elem->getParametersString();
169        //   string cliche = workingDir + param[2].c_str();
170        //   if( beamToParmela(cliche,lastDiag) ) {
171        //       //         cout  <<  "["  <<  k  <<  "] : ecrit sur fichier " << cliche << " le contenu de beam["<<avantDernier<<"]"<<endl;
172        //     cout  <<  "["  <<  k  <<  "] : ecrit sur fichier " << cliche << " le contenu de beam[ ]"<<endl;
173        //   }
174        // }
175        // si le numero est reconnu de parmela
176      } else {
177         
178        // on initialise une nouvelle sortie diagnostic
179          particleBeam* newDiag = dataManager_->updateCurrentDiagnostic(true);
180          vector<double> centroid;
181          bareParticle refPart;
182          vector<bareParticle> particles;
183          double freq= globParamPtr_->getFrequency();
184          unsigned numeroParmel;
185          numeroParmel = (unsigned)relativeParmelaElementIndices_.at(curseur);
186          cout << " lecture PARMELA el no absolu " << k << " numero relatif " << numeroParmel << " nom " << elem->getNomdElement().getElementName() << endl;
187        // lecture sortie parmela
188        if(!beamFromParmela(workingDir,numeroParmel,freq,centroid,refPart,particles))
189          {
190            // si echec, fin
191            dataManager_->consoleMessage(" softwareParmela::buildBeamAfterElements : failure in reading parmela result for element " + elem->getLabel() + " for unknown reason " );
192            return  false;
193          } else {
194          // si succes, on complete le diagnostic
195          newDiag->setWithParticles(centroid,refPart,particles);
196        }
197      }
198    }
199  return result;
200}
201bool softwareParmela::beamFromParmela(string workingDir,unsigned numeroParmel, double referencefrequency, vector<double>& centroid, bareParticle& refPart,vector<bareParticle>& particles )
202{   
203  string nomfilefais = workingDir + "parmdesz";
204  cout << " nom fichier desz : " << nomfilefais << endl;
205  FILE *filefais = fopen(nomfilefais.c_str(), "r");
206 
207  if ( filefais == (FILE*)0 ) {
208    dataManager_->consoleMessage(" beamFromParmela() erreur a l'ouverture du fichier 'parmdesz'");
209    cerr << " beamFromParmela() erreur a l'ouverture du fichier" << nomfilefais  << endl;
210    return false;
211  } else
212    cout << " beamFromParmela() : ouverture du fichier " << nomfilefais << endl;
213 
214  parmelaParticle partic;
215  std::vector<parmelaParticle> faisceau;
216  int testNombrePartRef =0;
217  double phaseRef = 0.0;
218 
219  while( partic.readFromParmelaFile(filefais) > 0 )
220    {
221     
222      if ( partic.ne == (int)numeroParmel ) {
223        if ( partic.np == 1 ) {
224          // en principe on est sur la particule de reference
225          if ( fabs(partic.xx) > EPSILON || fabs(partic.yy) > EPSILON || fabs(partic.xxp) > EPSILON  || fabs(partic.yyp) > EPSILON) {
226            printf(" ATTENTION part. reference douteuse  \n");
227            partic.imprim();
228          }
229          phaseRef = partic.phi;
230         
231          // le 'z' est 'absolu' (le long de la trajectoire)
232          TRIDVECTOR  posRef(partic.xx,partic.yy,partic.z);
233          TRIDVECTOR betagammaRef(partic.xxp*partic.begamz, partic.yyp*partic.begamz, partic.begamz);
234          refPart = bareParticle(posRef, betagammaRef);
235          testNombrePartRef++;
236          if ( testNombrePartRef != 1 ) {
237            dataManager_->consoleMessage(" beamFromParmela : nombre de part. de ref different de 1 :");
238            cerr << " nombre de part. de ref different de 1 : " << testNombrePartRef << " !! " << endl;
239            return false;
240          }
241        }
242        faisceau.push_back(partic);
243      }
244    } //while
245 
246  if ( faisceau.size() == 0)
247    {
248      string stringNum = mixedTools::intToString( (int)numeroParmel );
249      dataManager_->consoleMessage("beamFromParmela echec lecture  element numero relatif  parmela : " + stringNum);
250      cerr << " beamFromParmela echec lecture  element numero relatif  parmela " << numeroParmel << endl;
251      return false;
252    }
253 
254  // facteur  c/ 360. pour calculer (c dphi) / (360.freq)
255  // avec freq en Mhz et dphi en degres et résultat en cm:
256  double FACTEUR =  83.3333;  // ameliorer la precision
257 
258  // pour l'instant on choisit un centroid nul;
259  centroid.clear();
260  centroid = vector<double>(6,0.0);
261 
262  particles.clear();
263  particles.resize(faisceau.size(), bareParticle());
264  double x,xp,y,yp;
265  double betagammaz;
266  double betaz;
267  double deltaz;
268  double dephas;
269  double g;
270  TRIDVECTOR  pos;
271  TRIDVECTOR betagamma;
272  // contrairement a ce qu'indique la notice PARMELA, dans parmdesz, les xp et yp
273  // sont donnes en radians
274  for (unsigned k = 0; k < faisceau.size(); k++) {
275    x= faisceau.at(k).xx;
276    xp=faisceau.at(k).xxp;
277    y= faisceau.at(k).yy;
278    yp=faisceau.at(k).yyp;
279    // dephasage par rapport a la reference
280    dephas = faisceau.at(k).phi - phaseRef; // degrés
281    g = faisceau.at(k).wz/EREST_MeV;
282    betagammaz = faisceau.at(k).begamz;
283    betaz = betagammaz/(g+1.0);
284    deltaz = FACTEUR * betaz * dephas / referencefrequency;
285    x += xp * deltaz;
286    y += yp * deltaz;
287    pos.setComponents(x,y,deltaz);
288    betagamma.setComponents(xp*betagammaz, yp*betagammaz, betagammaz);
289    particles.at(k) = bareParticle(pos,betagamma);
290  }
291 
292  return true;
293}
294
295// sauvegarde d'un 'particleBeam' sur un fichier parmela, en guise d'INPUT
296// pour l'instant de nom standard 'parin.input0'
297bool softwareParmela::beamToParmela(string nameOfFile,particleBeam* beam)
298{
299    if ( !beam->particleRepresentationOk() ) {
300        dataManager_->consoleMessage("softwareParmela::beamToParmela : beam not in particles form : not yet programmed");
301        cout << " softwareParmela::beamToParmela : beam not in particles form : not yet programmed " << endl;
302        return false;
303    }
304   
305    ofstream outfile;
306    outfile.open(nameOfFile.c_str(),ios::out);
307    if (!outfile) {
308        dataManager_->consoleMessage(" softwareParmela::beamToParmela : error opening output stream ");
309        cerr << " softwareParmela::beamToParmela : error opening output stream " << nameOfFile << endl;
310        return false;
311    }
312   
313    const vector<bareParticle>& partic = beam->getParticleVector();
314    double weight = 1.0;
315    double xx,yy,zz;
316    double begamx, begamy, begamz;
317    for (unsigned k = 0; k < partic.size(); k++) {
318        partic.at(k).getPosition().getComponents(xx,yy,zz);
319        partic.at(k).getBetaGamma().getComponents(begamx,begamy,begamz);
320        outfile << xx << " " << begamx << " " <<  yy << " " << begamy << " " << zz << " " << begamz  << " " << weight << endl;
321    }
322    outfile.close();
323    return true;
324}
Note: See TracBrowser for help on using the repository browser.