source: PSPA/Interface_Web/branches/insertionsElements/src/particleBeam.cc @ 432

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

ajout d'un element 'beam' (transport)

File size: 14.6 KB
Line 
1#include "particleBeam.h"
2#include "mathematicalConstants.h"
3#include "PhysicalConstants.h"
4//#include <string>
5#include <stdio.h>
6#include <algorithm>
7
8#include "environmentVariables.h"
9
10
11
12particleBeam::particleBeam()  {
13  // rij_transportMoments_.resize(6);
14  // unsigned dim=0;
15  // unsigned k;
16  // for ( k=0; k < 6; k++){
17  //   rij_transportMoments_.at(k).resize(++dim);
18  // }
19  P0Transport_ = 0.0;
20  particleRepresentationOk_ = false;
21  momentRepresentationOk_ = false;
22}
23
24void particleBeam::clear() {
25  unsigned k,j;
26  goodPartic_.clear();
27  rij_.raz();
28  // for ( k=0; k < 6; k++){
29  //   for ( j=0; j <= k; j++) {
30  //     ( rij_transportMoments_.at(k) ).at(j) = 0.0;
31  //   }
32  // }
33  P0Transport_ = 0.0;
34  particleRepresentationOk_ = false;
35  momentRepresentationOk_ = false;
36}
37
38int particleBeam::getNbParticles() const {
39  return goodPartic_.size();
40}
41
42const transportMoments&   particleBeam::getTransportMoments() const  { 
43    return rij_;
44}
45double particleBeam::getP0Transport() const { 
46  return P0Transport_;
47}
48bool particleBeam::particleRepresentationOk() const {
49  return particleRepresentationOk_;
50}
51bool particleBeam::momentRepresentationOk() const {
52  return momentRepresentationOk_;
53}
54
55void  particleBeam::addParticle( bareParticle p)
56{
57  goodPartic_.push_back(p);
58}
59
60const vector<bareParticle>& particleBeam::getParticleVector() const
61{
62  return goodPartic_;
63}
64vector<bareParticle>& particleBeam::getParticleVector() 
65{
66  return goodPartic_;
67}
68
69  // void particleBeam::setFromBeamParameters(double x, double xp, double y, double yp, double dl, double del, double P0)
70  // {
71  //   clear();
72  //   rij_ = transportMoments(x, xp, y, yp, dl, del);
73  //   P0Transport_ = P0;
74  //   momentRepresentationOk_ = true;
75  // }
76
77
78bool  particleBeam::setFromTransport(string elLab, const nomdElements elem)
79{
80  string elementLabel = elLab;
81  // transformer le label en majuscules ; je ne suis pas sur que ca
82  // marchera a tous les coups (glm)
83  std::transform(elementLabel.begin(), elementLabel.end(), elementLabel.begin(), (int (*)(int))std::toupper);
84
85  cout << " particleBeam::setFromTransport on cherche element " << elementLabel << endl;
86  string buf;
87  ifstream infile;
88  string nameIn = WORKINGAREA + "transport.output";
89  infile.open(nameIn.c_str(), ios::in);
90  if (!infile) {
91    cerr << " particleBeam::setFromTransport : error re-opening transport output stream " << nameIn << endl;
92    return false;
93  }
94  //  else cout << " particleBeam::setFromTransport() : ouverture du fichier " << nameIn << endl;
95
96  string::size_type nn = string::npos;
97  while ( getline(infile, buf) ) {
98    //      cout << " buf= "  << buf << endl;
99    nn = buf.find(elementLabel);
100    //     cout << " string::npos= " << string::npos << " nn= " << nn << endl;
101    if ( nn != string::npos ) {
102      //          cout << " particleBeam::setFromTransport : element " << elementLabel << " trouve " << endl;
103      break;
104    }
105  }
106
107  if ( nn == string::npos ) {
108    cerr << " particleBeam::setFromTransport : element " << elementLabel << " non trouve dans le fichier  " << nameIn << endl;
109    return false;
110  }
111  if (elem.getElementType() == bend ) {
112    getline(infile, buf);
113    getline(infile, buf);
114  }
115  readTransportMoments(infile);
116  //  impressionDesMoments();
117  infile.close();
118  momentRepresentationOk_ = true;
119  return true;
120}
121
122
123bool particleBeam::setFromParmela(unsigned numeroElement, double referencefrequency) {
124  unsigned  k;
125  FILE* filefais;
126  string nomfilefais = WORKINGAREA + "parmdesz";
127  cout << " nom fichier desz : " << nomfilefais << endl;
128  filefais = fopen(nomfilefais.c_str(), "r");
129
130  if ( filefais == (FILE*)0 ) {
131    cerr << " particleBeam::setFromParmela() erreur a l'ouverture du fichier" << nomfilefais  << endl;; 
132    return false;
133  }
134  else cout << " particleBeam::setFromParmela() : ouverture du fichier " << nomfilefais << endl; 
135
136  parmelaParticle partic;
137  std::vector<parmelaParticle> faisceau;
138
139  cout << " particleBeam::setFromParmela : numeroElement = " << numeroElement << endl;
140  numeroElement--;
141
142
143
144  int testNombrePartRef =0;
145  double phaseRef;
146
147  while( partic.readFromParmelaFile(filefais) > 0 ) {
148    if ( partic.ne == numeroElement )
149      {
150        faisceau.push_back(partic);
151
152        if ( partic.np == 1 ) {
153          // en principe on est sur la particule de reference
154          if ( fabs(partic.xx) > EPSILON || fabs(partic.yy) > EPSILON || fabs(partic.xxp) > EPSILON  || fabs(partic.yyp) > EPSILON) {
155            printf(" ATTENTION part. reference douteuse  \n");
156            partic.imprim();
157          }
158          phaseRef = partic.phi;
159          TRIDVECTOR  posRef(partic.xx,partic.yy,0.0);
160          TRIDVECTOR betagammaRef(partic.xxp*partic.begamz, partic.yyp*partic.begamz, partic.begamz);
161          referenceParticle_ = bareParticle(posRef, betagammaRef);
162          testNombrePartRef++;
163          if ( testNombrePartRef != 1 ) {
164            cerr << " TROP DE PART. DE REF : " << testNombrePartRef << " !! " << endl;
165            return false;
166          }
167        }
168      }
169  }
170
171  if ( faisceau.size() == 0) 
172    {
173      cerr << " particleBeam::setFromParmela echec lecture  element " << numeroElement+1 << endl;
174      return false;
175    }
176 
177  // facteur  c/ 360. pour calculer (c dphi) / (360.freq)
178  // avec freq en Mhz et dphi en degres et résultat en cm:
179  double FACTEUR =  83.3333;  // ameliorer la precision
180
181
182
183  // pour l'instant on choisit un centroid nul;
184  centroid_ = vector<double>(6,0.0);
185
186  goodPartic_.clear();
187  goodPartic_.resize(faisceau.size(), bareParticle());
188  double x,xp,y,yp;
189  double betagammaz;
190  double betaz;
191  double deltaz;
192  double dephas;
193  double g;
194  TRIDVECTOR  pos;
195  TRIDVECTOR betagamma;
196  // contrairement a ce qu'indique la notice PARMELA, dans parmdesz, les xp et yp
197  // sont donnes en radians
198  for ( k=0; k < faisceau.size(); k++) {
199    x=faisceau.at(k).xx;
200    xp=faisceau.at(k).xxp;
201    y=faisceau.at(k).yy;
202    yp=faisceau.at(k).yyp;
203
204    // dephasage par rapport a la reference 
205    dephas = faisceau.at(k).phi - phaseRef; // degrés
206    g = faisceau.at(k).wz/ERESTMeV;
207    betagammaz = faisceau.at(k).begamz;
208    betaz = betagammaz/(g+1.0);
209    deltaz = FACTEUR * betaz * dephas / referencefrequency;
210    x += xp * deltaz;
211    y += yp * deltaz;
212    pos.setComponents(x,y,deltaz);
213    betagamma.setComponents(xp*betagammaz, yp*betagammaz, betagammaz);
214    goodPartic_.at(k) = bareParticle(pos,betagamma);
215  }
216  particleRepresentationOk_ = true;
217  return true;
218}
219
220
221void particleBeam::getVariance(double& varx, double& vary, double& varz) const {
222  unsigned int k;
223  double x,y,z;
224  double xav = 0.;
225  double yav = 0.;
226  double zav = 0.;
227  double xavsq = 0.;
228  double yavsq = 0.;
229  double zavsq = 0.;
230
231  TRIDVECTOR pos;
232
233
234  for ( k = 0 ; k < goodPartic_.size(); k++) {
235    pos = goodPartic_.at(k).getPosition();
236    pos.getComponents(x,y,z);
237    //      partic_[k].getXYZ(x,y,z);
238    xav += x;
239    xavsq += x*x;
240    yav += y;
241    yavsq += y*y;
242    zav += z;
243    zavsq += z*z;
244  }
245
246  double aginv = double (goodPartic_.size());
247  aginv = 1.0/aginv;
248
249  varx =  aginv * ( xavsq - xav*xav*aginv );
250  vary =  aginv * ( yavsq - yav*yav*aginv );
251  varz =  aginv * ( zavsq - zav*zav*aginv );
252}
253
254
255void particleBeam::printAllXYZ() const {
256  cout << " dump du faisceau : " << endl;
257  cout <<  goodPartic_.size() << " particules " << endl;
258  unsigned int k;
259  for ( k = 0 ; k < goodPartic_.size(); k++)
260    {
261      double xx,yy,zz;
262      goodPartic_.at(k).getPosition().getComponents(xx,yy,zz);
263      cout << " part. numero " << k << "  x= " << xx << " y= " << yy  << " z= " << zz << endl;
264    }
265}
266
267
268
269void particleBeam::Zrange(double& zmin, double& zmax) const {
270  double z;
271  zmin = GRAND;
272  zmax = -zmin;
273
274  unsigned int k;
275  for ( k = 0 ; k < goodPartic_.size(); k++)
276    {
277      z = goodPartic_.at(k).getZ();
278      if ( z < zmin ) zmin = z;
279      else if ( z > zmax) zmax = z;         
280    }
281}
282
283
284
285string particleBeam::FileOutputFlow() const {
286  ostringstream sortie;
287  unsigned int k;
288  for ( k = 0 ; k < goodPartic_.size(); k++)
289    {
290      sortie << goodPartic_.at(k).FileOutputFlow() << endl;
291    }
292  sortie << endl;
293  return sortie.str();
294}
295
296bool particleBeam::FileInput( ifstream& ifs) {
297  bool test = true;
298  string dum1, dum2;
299  double dummy;
300  if ( !( ifs >> dum1 >> dum2 >> dummy) ) return false;
301 
302  bareParticle pp;
303  while ( pp.FileInput(ifs) )
304    {
305      addParticle( pp);
306    }
307  return test;
308}
309
310void particleBeam::buildMomentRepresentation() {
311
312  unsigned k,j,m;
313  double auxj, auxm;
314  if ( !particleRepresentationOk_)
315    {
316      cerr << " particleBeam::buildMomentRepresentation() vecteur de particules invalide" << endl;
317      return;
318    }
319
320  cout << " buildMomentRepresentation " << endl;
321  //  printAllXYZ();
322
323  double gref = referenceParticle_.getGamma() - 1.0;
324  double P_reference_MeV_sur_c = sqrt( gref*(gref+2) );
325
326  // initialisation des moments
327  razDesMoments();
328
329  // accumulation
330  TRIDVECTOR pos;
331  TRIDVECTOR begam;
332  double gamma;
333  double begamz;
334  double g;
335  double PMeVsc;
336  double del;
337  vector<double> part(6);
338
339    vector< vector<double> >& matrice = rij_.getMatrix();
340
341
342  for (k=0; k < goodPartic_.size(); k++) {
343    gamma = goodPartic_.at(k).getGamma();
344    pos = goodPartic_.at(k).getPosition();
345    begam= goodPartic_.at(k).getBetaGamma();
346    begamz = begam.getComponent(2);
347    g = gamma -1.0;
348    PMeVsc = sqrt( g*(g+2) );
349    del = 100.0 * ( PMeVsc -  P_reference_MeV_sur_c ) / P_reference_MeV_sur_c ; // en %
350
351    part[0] = pos.getComponent(0);
352    part[1] = begam.getComponent(0)/begamz;
353    part[2] = pos.getComponent(1);
354    part[3] = begam.getComponent(1)/begamz;
355    part[4] = pos.getComponent(2);
356    part[5] = del;
357
358    for ( j = 0; j < 6; j++) {
359      auxj = part.at(j) - centroid_.at(j);
360      for (m=0; m <= j; m++) 
361        {
362          auxm = part.at(m) - centroid_.at(m);
363
364          ( matrice.at(j) ).at(m) += auxj*auxm;
365          //      ( rij_transportMoments_.at(j) ).at(m) += auxj*auxm;
366
367
368          //          cout << " j= " << j << " m= " << m << " rjm= " << ( rij_transportMoments_.at(j) ).at(m) << endl;
369        }
370    }
371  }
372
373
374  // moyenne
375  double facmoy = 1.0/double( goodPartic_.size() );
376  for ( j = 0; j < 6; j++) {
377        ( matrice.at(j) ).at(j) = sqrt(( matrice.at(j) ).at(j) * facmoy );
378    //    ( rij_transportMoments_.at(j) ).at(j) = sqrt(( rij_transportMoments_.at(j) ).at(j) * facmoy );
379  }
380
381  for ( j = 0; j < 6; j++) {
382    auxj =  ( matrice.at(j) ).at(j);
383    //    auxj =  ( rij_transportMoments_.at(j) ).at(j);
384    for (m=0; m < j; m++) {
385      auxm = ( matrice.at(m) ).at(m);
386      (  matrice.at(j) ).at(m) *= facmoy/(auxj * auxm);
387      //      auxm = ( rij_transportMoments_.at(m) ).at(m);
388      // ( rij_transportMoments_.at(j) ).at(m) *= facmoy/(auxj * auxm);
389    }
390  }
391   
392  // les longueurs sont en cm
393  // les angles en radians, on passe en mrad;
394
395  double uniteAngle = 1.0e+3;
396  ( matrice.at(1) ).at(1)  *= uniteAngle;
397  ( matrice.at(3) ).at(3)  *= uniteAngle;
398
399  // ( rij_transportMoments_.at(1) ).at(1)  *= uniteAngle;
400  // ( rij_transportMoments_.at(3) ).at(3)  *= uniteAngle;
401   
402  P0Transport_ = 1.0e-3*ERESTMeV*P_reference_MeV_sur_c;
403
404  //  cout << " buildmomentrepresentation impression des moments " << endl;
405  //  impressionDesMoments();
406
407  momentRepresentationOk_ = true;
408}
409
410void particleBeam::impressionDesMoments() const {
411  rij_.impression();
412  // unsigned j,m;
413  // cout << " impression des moments " << endl;
414  // for ( j = 0; j < 6; j++)
415  //   {
416  //     for (m=0; m <= j; m++)
417  //    {
418  //      cout  << ( rij_transportMoments_.at(j) ).at(m) << " ";
419  //    }
420  //     cout << endl;
421  //   }
422}
423
424void particleBeam::razDesMoments() {
425  rij_.raz();
426  // // initialisation des moments
427  // unsigned j,m;
428  // for ( j = 0; j < rij_transportMoments_.size(); j++)
429  //   {
430  //     for (m=0; m <= j; m++)
431  //    {
432  //      ( rij_transportMoments_.at(j) ).at(m) = 0.0; // element r_jm
433  //    }
434  //   }
435}
436
437
438void particleBeam::readTransportMoments(ifstream& inp) {
439rij_.readFromTransportOutput(inp);
440  // string  bidString;
441  // double bidon;
442
443  // // initialisation des moments
444  // razDesMoments();
445
446  // inp >> bidon >>  bidString >>  bidon >>  ( rij_transportMoments_.at(0) ).at(0) >> bidString;
447  // inp >> bidon >> ( rij_transportMoments_.at(1) ).at(1) >>  bidString >> ( rij_transportMoments_.at(1) ).at(0);
448  // inp >> bidon >> ( rij_transportMoments_.at(2) ).at(2) >>  bidString >> ( rij_transportMoments_.at(2) ).at(0)  >> ( rij_transportMoments_.at(2) ).at(1);
449  // inp >> bidon >> ( rij_transportMoments_.at(3) ).at(3) >>  bidString >> ( rij_transportMoments_.at(3) ).at(0)  >> ( rij_transportMoments_.at(3) ).at(1) >> ( rij_transportMoments_.at(3) ).at(2);
450
451  // inp >> bidon >> ( rij_transportMoments_.at(4) ).at(4) >>  bidString >> ( rij_transportMoments_.at(4) ).at(0)  >> ( rij_transportMoments_.at(4) ).at(1) >> ( rij_transportMoments_.at(4) ).at(2) >> ( rij_transportMoments_.at(4) ).at(3);
452
453  // inp >> bidon >> ( rij_transportMoments_.at(5) ).at(5) >>  bidString >> ( rij_transportMoments_.at(5) ).at(0)  >> ( rij_transportMoments_.at(5) ).at(1) >> ( rij_transportMoments_.at(5) ).at(2) >> ( rij_transportMoments_.at(5) ).at(3) >> ( rij_transportMoments_.at(5) ).at(4);
454
455}
456
457double particleBeam::getXmaxRms() {
458  if ( !momentRepresentationOk_ ) buildMomentRepresentation();
459  return ( rij_.getMatrix().at(0) ).at(0);
460  // return ( rij_transportMoments_.at(0) ).at(0);
461}
462
463void particleBeam::donneesDessinEllipseXxp(vector<double>& xcor, vector<double>& ycor) {
464  int k;
465  double x,y;
466
467  if ( !momentRepresentationOk_ ) return;
468
469  xcor.clear();
470  ycor.clear();
471
472  double xm = ( rij_.getMatrix().at(0) ).at(0);
473  double ym = ( rij_.getMatrix().at(1) ).at(1);
474  double r  = ( rij_.getMatrix().at(1) ).at(0);
475
476  // double xm = ( rij_transportMoments_.at(0) ).at(0);
477  // double ym = ( rij_transportMoments_.at(1) ).at(1);
478  // double r  = ( rij_transportMoments_.at(1) ).at(0);
479
480  int nbintv = 50;
481  if ( xm == 0.0 ) return;
482  double pas = 2.0 * xm / nbintv;
483
484  //  cout << " r= " << r << endl;
485  double rac = (1 - r*r);
486  if ( rac > 0.0 ) 
487    {
488      rac = sqrt(1 - r*r);
489      double alpha = -r / rac;
490      double beta = xm / ( ym * rac);
491      //  double gamma = ym / ( xm * rac );
492      double epsil = xm * ym * rac;
493      double fac1 = -1.0 / ( beta * beta);
494      double fac2 = epsil/beta;
495      double fac3 = -alpha/beta;
496      double aux;
497      for ( k=0; k < nbintv; k++)
498        {
499          x = -xm + k*pas;
500          aux = fac1 * x * x + fac2;
501          //     cout << " aux2= " << aux << endl;
502          if ( aux <= 0.0 )
503            {
504              aux = 0.0;
505            }
506          else aux = sqrt(aux);
507     
508          //        y = fac3*x;
509          y = fac3*x + aux;
510          xcor.push_back(x);
511          ycor.push_back(y);
512        }
513
514      for ( k=0; k <= nbintv; k++)
515        {
516          x = xm - k*pas;
517          aux =  fac1 * x * x + fac2;
518          if ( aux <= 0.0 ) 
519            {
520              aux = 0.0;
521            }
522          else aux = sqrt(aux);
523          //   y = fac3*x;
524          y = fac3*x - aux;
525          xcor.push_back(x);
526          ycor.push_back(y);
527        }
528    }
529  else
530    // cas degenere
531    {
532      double fac = -ym/xm;
533      for ( k=0; k < nbintv; k++)
534        {
535          x = -xm + k*pas;
536          y = fac*x;
537          xcor.push_back(x);
538          ycor.push_back(y);
539        }
540       
541    }
542}
Note: See TracBrowser for help on using the repository browser.