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

Last change on this file since 417 was 417, checked in by lemeur, 11 years ago

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