source: PSPA/Interface_Web/trunk/pspaWT/src/particleBeam.cc @ 236

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

epuration particleBeam

File size: 18.0 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#include <sstream>
8//#include "environmentVariables.h"
9#include <Wt/WApplication>
10
11using namespace std;
12
13particleBeam::particleBeam()  {
14  // rij_transportMoments_.resize(6);
15  // unsigned dim=0;
16  // unsigned k;
17  // for ( k=0; k < 6; k++){
18  //   rij_transportMoments_.at(k).resize(++dim);
19  // }
20  P0Transport_ = 0.0;
21  particleRepresentationOk_ = false;
22  momentRepresentationOk_ = false;
23}
24
25void particleBeam::clear() {
26  goodPartic_.clear();
27  rij_.raz();
28  P0Transport_ = 0.0;
29  particleRepresentationOk_ = false;
30  momentRepresentationOk_ = false;
31}
32
33int particleBeam::getNbParticles() const {
34  return goodPartic_.size();
35}
36
37const beam2Moments&   particleBeam::getTransportMoments() const  { 
38    return rij_;
39}
40
41double particleBeam::getSigmaTransportij(unsigned i, unsigned j)  {
42  if ( i < 1 || i > 6 || j < 1 || j > 6 ) {
43    cerr << " particleBeam::getSigmaTransportij() indices out of range  " << endl;
44    return 0.0;
45  }
46  if ( !momentRepresentationOk_ ) {
47    cerr << " particleBeam::getSigmaTransportij() beam is not in moment representation " << endl;
48    return 0.0;
49  }
50
51  i--;
52  j--;
53  if ( j > i ) {
54    unsigned aux = i;
55    i = j;
56    j = aux;
57  }
58  return ( rij_.getMatrix().at(i) ).at(j);
59}
60
61
62double particleBeam::getUnnormalizedEmittanceX() {
63  double r = getSigmaTransportij(2,1);
64  double rac = (1 - r*r);
65  if ( rac <= 0.0 ) {
66    return 0.0;
67  }
68  rac = sqrt(1 - r*r);
69  return  10.*getSigmaTransportij(1,1) * getSigmaTransportij(2,2) * rac; // en pi.mm.mrad
70}
71
72double particleBeam::getP0Transport() const { 
73  return P0Transport_;
74}
75
76void particleBeam::set2Moments(beam2Moments& moments) {
77  rij_ = moments;
78  momentRepresentationOk_ = true;
79}
80
81
82
83bool particleBeam::particleRepresentationOk() const {
84  return particleRepresentationOk_;
85}
86bool particleBeam::momentRepresentationOk() const {
87  return momentRepresentationOk_;
88}
89
90void  particleBeam::addParticle( bareParticle p)
91{
92  goodPartic_.push_back(p);
93}
94
95const vector<bareParticle>& particleBeam::getParticleVector() const
96{
97  return goodPartic_;
98}
99vector<bareParticle>& particleBeam::getParticleVector() 
100{
101  return goodPartic_;
102}
103
104
105// bool  particleBeam::setFromTransport(string workingDir, string elLab, const nomdElements elem)
106// {
107//   string elementLabel = elLab;
108//   // transformer le label en majuscules ; je ne suis pas sur que ca
109//   // marchera a tous les coups (glm)
110//   std::transform(elementLabel.begin(), elementLabel.end(), elementLabel.begin(), (int (*)(int))std::toupper);
111
112//   cout << " particleBeam::setFromTransport on cherche element " << elementLabel << endl;
113//   string buf;
114//   ifstream infile;
115//   //  string nameIn = WORKINGAREA + "transport.output";
116//   string nameIn = workingDir + "transport.output";
117//   infile.open(nameIn.c_str(), ios::in);
118//   if (!infile) {
119//     cerr << " particleBeam::setFromTransport : error re-opening transport output stream " << nameIn << endl;
120//     return false;
121//   }
122//   //  else cout << " particleBeam::setFromTransport() : ouverture du fichier " << nameIn << endl;
123
124//   string::size_type nn = string::npos;
125//   while ( getline(infile, buf) ) {
126//     //      cout << " buf= "  << buf << endl;
127//     nn = buf.find(elementLabel);
128//     //     cout << " string::npos= " << string::npos << " nn= " << nn << endl;
129//     if ( nn != string::npos ) {
130//       //       cout << " particleBeam::setFromTransport : element " << elementLabel << " trouve " << endl;
131//       break;
132//     }
133//   }
134
135//   if ( nn == string::npos ) {
136//     cerr << " particleBeam::setFromTransport : element " << elementLabel << " non trouve dans le fichier  " << nameIn << endl;
137//     return false;
138//   }
139//   if (elem.getElementType() == bend ) {
140//     getline(infile, buf);
141//     getline(infile, buf);
142//   }
143//   readTransportMoments(infile);
144//   impressionDesMoments();
145//   infile.close();
146//   momentRepresentationOk_ = true;
147//   return true;
148// }
149
150// bool  particleBeam::setFromTransport(string workingDir, string elLab, const nomdElements elem)
151// {
152//   string elementLabel = elLab;
153//   // transformer le label en majuscules ; je ne suis pas sur que ca
154//   // marchera a tous les coups (glm)
155//   std::transform(elementLabel.begin(), elementLabel.end(), elementLabel.begin(), (int (*)(int))std::toupper);
156
157//   cout << " particleBeam::setFromTransport on cherche element " << elementLabel << endl;
158//   ifstream infile;
159//   //  string nameIn = WORKINGAREA + "transport.output";
160//   string nameIn = workingDir + "transport.output";
161//   infile.open(nameIn.c_str(), ios::in);
162//   if (!infile) {
163//     cerr << " particleBeam::setFromTransport : error re-opening transport output stream " << nameIn << endl;
164//     return false;
165//   }
166//   //  else cout << " particleBeam::setFromTransport() : ouverture du fichier " << nameIn << endl;
167
168//   string::size_type nn = string::npos;
169//   string  fichier;
170//   string buf;
171//   unsigned compteur = 0;
172//   while ( getline(infile, buf) ) {
173//       fichier += buf+"\n";
174//     nn = buf.find(elementLabel);
175//     if ( nn != string::npos ) {
176//       compteur++;
177//       //      break;
178//     }
179//   }
180//   infile.close();
181//   cout << " compteur= " << compteur << endl;
182//   if ( compteur == 0 ) {
183//     cerr << " particleBeam::setFromTransport : element " << elementLabel << " non trouve dans le fichier  " << nameIn << endl;
184//     return false;
185//   }
186
187
188//   stringstream fichierStream(fichier);
189//   buf.clear();
190//   unsigned relu = 0;
191//   while ( std::getline(fichierStream, buf) ) {
192//     //    cout << " relecture buffer : " << buf << endl;
193//     nn = buf.find(elementLabel);
194//     if ( nn != string::npos ) {
195//       relu++;
196//       if ( relu == compteur ) {
197//      cout << " TROUVE !" << endl;
198//      break;
199//       }
200//     }
201//   }
202//
203//   readTransportMoments(fichierStream);
204//   //  impressionDesMoments();
205//   momentRepresentationOk_ = true;
206//   return true;
207// }
208
209
210bool particleBeam::setFromParmela(string workingDir,unsigned numeroElement, double referencefrequency) {
211  unsigned  k;
212  FILE* filefais;
213  string nomfilefais = workingDir + "parmdesz";
214  cout << " nom fichier desz : " << nomfilefais << endl;
215  filefais = fopen(nomfilefais.c_str(), "r");
216
217  if ( filefais == (FILE*)0 ) {
218    cerr << " particleBeam::setFromParmela() erreur a l'ouverture du fichier" << nomfilefais  << endl;; 
219    return false;
220  }
221  else cout << " particleBeam::setFromParmela() : ouverture du fichier " << nomfilefais << endl; 
222
223  parmelaParticle partic;
224  std::vector<parmelaParticle> faisceau;
225
226  cout << " particleBeam::setFromParmela : numeroElement = " << numeroElement << endl;
227  unsigned indexElement = numeroElement-1;
228 
229
230
231
232  int testNombrePartRef =0;
233  double phaseRef;
234
235  while( partic.readFromParmelaFile(filefais) > 0 ) {
236    if ( partic.ne == (int)indexElement )
237      {
238        faisceau.push_back(partic);
239
240        if ( partic.np == 1 ) {
241          // en principe on est sur la particule de reference
242          if ( fabs(partic.xx) > EPSILON || fabs(partic.yy) > EPSILON || fabs(partic.xxp) > EPSILON  || fabs(partic.yyp) > EPSILON) {
243            printf(" ATTENTION part. reference douteuse  \n");
244            partic.imprim();
245          }
246          phaseRef = partic.phi;
247          TRIDVECTOR  posRef(partic.xx,partic.yy,0.0);
248          TRIDVECTOR betagammaRef(partic.xxp*partic.begamz, partic.yyp*partic.begamz, partic.begamz);
249          referenceParticle_ = bareParticle(posRef, betagammaRef);
250          testNombrePartRef++;
251          if ( testNombrePartRef != 1 ) {
252            cerr << " TROP DE PART. DE REF : " << testNombrePartRef << " !! " << endl;
253            return false;
254          }
255        }
256      }
257  }
258
259  if ( faisceau.size() == 0) 
260    {
261      cerr << " particleBeam::setFromParmela echec lecture  element " << numeroElement << endl;
262      return false;
263    }
264 
265  // facteur  c/ 360. pour calculer (c dphi) / (360.freq)
266  // avec freq en Mhz et dphi en degres et résultat en cm:
267  double FACTEUR =  83.3333;  // ameliorer la precision
268
269
270
271  // pour l'instant on choisit un centroid nul;
272  centroid_ = vector<double>(6,0.0);
273
274  goodPartic_.clear();
275  goodPartic_.resize(faisceau.size(), bareParticle());
276  double x,xp,y,yp;
277  double betagammaz;
278  double betaz;
279  double deltaz;
280  double dephas;
281  double g;
282  TRIDVECTOR  pos;
283  TRIDVECTOR betagamma;
284  // contrairement a ce qu'indique la notice PARMELA, dans parmdesz, les xp et yp
285  // sont donnes en radians
286  for ( k=0; k < faisceau.size(); k++) {
287    x=faisceau.at(k).xx;
288    xp=faisceau.at(k).xxp;
289    y=faisceau.at(k).yy;
290    yp=faisceau.at(k).yyp;
291
292    // dephasage par rapport a la reference 
293    dephas = faisceau.at(k).phi - phaseRef; // degrés
294    g = faisceau.at(k).wz/ERESTMeV;
295    betagammaz = faisceau.at(k).begamz;
296    betaz = betagammaz/(g+1.0);
297    deltaz = FACTEUR * betaz * dephas / referencefrequency;
298    x += xp * deltaz;
299    y += yp * deltaz;
300    pos.setComponents(x,y,deltaz);
301    betagamma.setComponents(xp*betagammaz, yp*betagammaz, betagammaz);
302    goodPartic_.at(k) = bareParticle(pos,betagamma);
303  }
304  particleRepresentationOk_ = true;
305  return true;
306}
307
308
309void particleBeam::getVariance(double& varx, double& vary, double& varz) const {
310  unsigned int k;
311  double x,y,z;
312  double xav = 0.;
313  double yav = 0.;
314  double zav = 0.;
315  double xavsq = 0.;
316  double yavsq = 0.;
317  double zavsq = 0.;
318
319  TRIDVECTOR pos;
320
321
322  for ( k = 0 ; k < goodPartic_.size(); k++) {
323    pos = goodPartic_.at(k).getPosition();
324    pos.getComponents(x,y,z);
325    //      partic_[k].getXYZ(x,y,z);
326    xav += x;
327    xavsq += x*x;
328    yav += y;
329    yavsq += y*y;
330    zav += z;
331    zavsq += z*z;
332  }
333
334  double aginv = double (goodPartic_.size());
335  aginv = 1.0/aginv;
336
337  varx =  aginv * ( xavsq - xav*xav*aginv );
338  vary =  aginv * ( yavsq - yav*yav*aginv );
339  varz =  aginv * ( zavsq - zav*zav*aginv );
340}
341
342
343void particleBeam::printAllXYZ() const {
344  cout << " dump du faisceau : " << endl;
345  cout <<  goodPartic_.size() << " particules " << endl;
346  unsigned int k;
347  for ( k = 0 ; k < goodPartic_.size(); k++)
348    {
349      double xx,yy,zz;
350      goodPartic_.at(k).getPosition().getComponents(xx,yy,zz);
351      cout << " part. numero " << k << "  x= " << xx << " y= " << yy  << " z= " << zz << endl;
352    }
353}
354
355
356
357void particleBeam::Zrange(double& zmin, double& zmax) const {
358  double z;
359  zmin = GRAND;
360  zmax = -zmin;
361
362  unsigned int k;
363  for ( k = 0 ; k < goodPartic_.size(); k++)
364    {
365      z = goodPartic_.at(k).getZ();
366      if ( z < zmin ) zmin = z;
367      else if ( z > zmax) zmax = z;         
368    }
369}
370
371
372
373string particleBeam::FileOutputFlow() const {
374  ostringstream sortie;
375  unsigned int k;
376  for ( k = 0 ; k < goodPartic_.size(); k++)
377    {
378      sortie << goodPartic_.at(k).FileOutputFlow() << endl;
379    }
380  sortie << endl;
381  return sortie.str();
382}
383
384bool particleBeam::FileInput( ifstream& ifs) {
385  bool test = true;
386  string dum1, dum2;
387  double dummy;
388  if ( !( ifs >> dum1 >> dum2 >> dummy) ) return false;
389 
390  bareParticle pp;
391  while ( pp.FileInput(ifs) )
392    {
393      addParticle( pp);
394    }
395  return test;
396}
397
398void particleBeam::buildMomentRepresentation() {
399
400  unsigned k,j,m;
401  double auxj, auxm;
402  if ( !particleRepresentationOk_)
403    {
404      cerr << " particleBeam::buildMomentRepresentation() vecteur de particules invalide" << endl;
405      return;
406    }
407
408  cout << " buildMomentRepresentation " << endl;
409  //  printAllXYZ();
410
411  double gref = referenceParticle_.getGamma() - 1.0;
412  double P_reference_MeV_sur_c = sqrt( gref*(gref+2) );
413
414  // initialisation des moments
415  razDesMoments();
416
417  // accumulation
418  TRIDVECTOR pos;
419  TRIDVECTOR begam;
420  double gamma;
421  double begamz;
422  double g;
423  double PMeVsc;
424  double del;
425  vector<double> part(6);
426
427    vector< vector<double> >& matrice = rij_.getMatrix();
428
429
430  for (k=0; k < goodPartic_.size(); k++) {
431    gamma = goodPartic_.at(k).getGamma();
432    pos = goodPartic_.at(k).getPosition();
433    begam= goodPartic_.at(k).getBetaGamma();
434    begamz = begam.getComponent(2);
435    g = gamma -1.0;
436    PMeVsc = sqrt( g*(g+2) );
437    del = 100.0 * ( PMeVsc -  P_reference_MeV_sur_c ) / P_reference_MeV_sur_c ; // en %
438
439    part[0] = pos.getComponent(0);
440    part[1] = begam.getComponent(0)/begamz;
441    part[2] = pos.getComponent(1);
442    part[3] = begam.getComponent(1)/begamz;
443    part[4] = pos.getComponent(2);
444    part[5] = del;
445
446    for ( j = 0; j < 6; j++) {
447      auxj = part.at(j) - centroid_.at(j);
448      for (m=0; m <= j; m++) 
449        {
450          auxm = part.at(m) - centroid_.at(m);
451
452          ( matrice.at(j) ).at(m) += auxj*auxm;
453          //      ( rij_transportMoments_.at(j) ).at(m) += auxj*auxm;
454
455
456          //          cout << " j= " << j << " m= " << m << " rjm= " << ( rij_transportMoments_.at(j) ).at(m) << endl;
457        }
458    }
459  }
460
461
462  // moyenne
463  double facmoy = 1.0/double( goodPartic_.size() );
464  for ( j = 0; j < 6; j++) {
465        ( matrice.at(j) ).at(j) = sqrt(( matrice.at(j) ).at(j) * facmoy );
466  }
467
468  for ( j = 0; j < 6; j++) {
469    auxj =  ( matrice.at(j) ).at(j);
470    for (m=0; m < j; m++) {
471      auxm = ( matrice.at(m) ).at(m);
472      (  matrice.at(j) ).at(m) *= facmoy/(auxj * auxm);
473    }
474  }
475   
476  ////////////////// si C21 = 1 , transport plante ! a voir //////////
477cout << " valeur initiale de  C21: " << ( matrice.at(1) ).at(0) << endl;
478  if ( ( matrice.at(1) ).at(0) >0.999999  ) {
479    ( matrice.at(1) ).at(0) = 0.999999;
480    cout << " j'ai fait la correction C21: " << ( matrice.at(1) ).at(0) << endl;
481  }
482 
483
484  // les longueurs sont en cm
485  // les angles en radians, on passe en mrad;
486
487  double uniteAngle = 1.0e+3;
488  ( matrice.at(1) ).at(1)  *= uniteAngle;
489  ( matrice.at(3) ).at(3)  *= uniteAngle;
490  P0Transport_ = 1.0e-3*ERESTMeV*P_reference_MeV_sur_c;
491
492  //  cout << " buildmomentrepresentation impression des moments " << endl;
493  //  impressionDesMoments();
494
495  momentRepresentationOk_ = true;
496}
497
498void particleBeam::impressionDesMoments() const {
499  rij_.impression();
500}
501
502void particleBeam::razDesMoments() {
503  rij_.raz();
504}
505
506
507// void particleBeam::readTransportMoments(ifstream& inp) {
508// rij_.readFromTransportOutput(inp);
509// }
510
511// void particleBeam::readTransportMoments(stringstream& inp) {
512// rij_.readFromTransportOutput(inp);
513// }
514
515double particleBeam::getXmaxRms() {
516  if ( !momentRepresentationOk_ ) buildMomentRepresentation();
517  return ( rij_.getMatrix().at(0) ).at(0);
518  // return ( rij_transportMoments_.at(0) ).at(0);
519}
520
521void particleBeam::donneesDessinEllipseXxp(vector<double>& xcor, vector<double>& ycor) {
522  int k;
523  double x,y;
524
525  if ( !momentRepresentationOk_ ) return;
526
527  xcor.clear();
528  ycor.clear();
529
530  double xm = ( rij_.getMatrix().at(0) ).at(0);
531  double ym = ( rij_.getMatrix().at(1) ).at(1);
532  double r  = ( rij_.getMatrix().at(1) ).at(0);
533
534  cout <<  " racs11= " << xm << " racs22= " << ym << " r12= " << r << endl;
535
536
537  int nbintv = 50;
538  if ( xm == 0.0 ) return;
539  double pas = 2.0 * xm / nbintv;
540
541  //  cout << " r= " << r << endl;
542  double rac = (1 - r*r);
543  if ( rac > 0.0 ) 
544    {
545      cout << " cas rac > " << endl;
546      rac = sqrt(1 - r*r);
547      double alpha = -r / rac;
548      double beta = xm / ( ym * rac);
549      //  double gamma = ym / ( xm * rac );
550      double epsil = xm * ym * rac;
551      double fac1 = -1.0 / ( beta * beta);
552      double fac2 = epsil/beta;
553      double fac3 = -alpha/beta;
554      double aux;
555      for ( k=0; k < nbintv; k++)
556        {
557          x = -xm + k*pas;
558          aux = fac1 * x * x + fac2;
559          //     cout << " aux2= " << aux << endl;
560          if ( aux <= 0.0 )
561            {
562              aux = 0.0;
563            }
564          else aux = sqrt(aux);
565     
566          //        y = fac3*x;
567          y = fac3*x + aux;
568          xcor.push_back(x);
569          ycor.push_back(y);
570        }
571
572      for ( k=0; k <= nbintv; k++)
573        {
574          x = xm - k*pas;
575          aux =  fac1 * x * x + fac2;
576          if ( aux <= 0.0 ) 
577            {
578              aux = 0.0;
579            }
580          else aux = sqrt(aux);
581          //   y = fac3*x;
582          y = fac3*x - aux;
583          xcor.push_back(x);
584          ycor.push_back(y);
585        }
586    }
587  else
588    // cas degenere
589    {
590      cout << " cas degenere " << endl;
591      double fac = ym/xm;
592      for ( k=0; k < nbintv; k++)
593        {
594          x = -xm + k*pas;
595          y = fac*x;
596          xcor.push_back(x);
597          ycor.push_back(y);
598        }
599       
600    }
601}
602
603void particleBeam::histogramme(vector<double>&xcor,vector<int>& hist,int& cnts)
604{
605  double gammin= GRAND;
606  double gammax= -gammin;
607  double Emoy= 0.0;
608  double ecatyp= 0.0;
609
610  for (unsigned int k = 0; k < goodPartic_.size(); k++)
611    {
612      double gamma = goodPartic_.at(k).getGamma();
613      double EMev = (gamma-1.0)*ERESTMeV;
614      if (gamma < gammin) gammin = gamma;
615      else if (gamma > gammax) gammax = gamma;
616      Emoy += EMev;
617      ecatyp += EMev*EMev;
618    }
619
620  double sum= (float)goodPartic_.size();
621  Emoy /= sum;
622  ecatyp /= sum;
623  ecatyp = sqrt(abs(ecatyp-Emoy*Emoy));
624
625  double Emin = (gammin-1.0)*ERESTMeV;
626  double Emax = (gammax-1.0)*ERESTMeV;
627  cout << "energie cinetique -moyenne " << Emoy << " Mev " << "-mini " << Emin << " Mev " << "-maxi " << Emax << " Mev " << "ecart type " << ecatyp*1000.0 << " Kev" << endl;
628
629  vector<double> Eshf;
630  for (unsigned int k = 0; k < goodPartic_.size(); k++)
631    {
632      double gamma = goodPartic_.at(k).getGamma();
633      double EMev = (gamma-1.0)*ERESTMeV;
634      Eshf.push_back(EMev-Emoy);
635    }
636
637  //////////////////////////////////////////////////////////////////////////////
638 
639  // demi fenetre en energie, et pas de l'histogramme
640  //  double hfene= max(3.*ecatyp-Emoy,Emoy-3.*ecatyp);
641  double hfene= 3.*ecatyp;
642  double hpas = hfene/25.;
643
644  cout << "demi fenetre " << hfene << ", hpas= " << hpas << endl;
645
646  double vmin = -hfene;
647  double dfen = 2.*hfene;
648  int ihist = dfen/hpas;
649  double phist = ihist*hpas;
650  double dpas = hpas-(dfen-phist);
651  if(dpas <= hpas*1.e-03) {
652    ihist++;
653    phist= ihist*hpas;
654  }
655  double vmax= vmin+hpas*ihist;
656 
657  cout << "'xAxisNumberOfBins= " << ihist <<", xAxisMinimum= " << vmin << ", xAxisMaximum= " << vmax << ", NParticules= " << Eshf.size() << ", phist " << phist << endl;
658
659  xcor= vector<double>(ihist);
660  for (int i = 0; i < ihist; ++i) {
661
662    // on gradue l'abcisse en pourcents
663    xcor[i]= 100.*( vmin+i*hpas );
664  }
665
666  hist= vector<int>(ihist,0);
667  for (unsigned i = 0; i < Eshf.size(); ++i) {
668    double var= Eshf[i]-vmin;
669    if(var < 0 || var >= phist) cout<<"out of range "<<var<<", ("<< i<<")"<< endl;
670    int k= var/hpas;
671    int kk= (int)floor(var/hpas);
672    if(i%20 == 0) cout<<"v("<<i<<")= " <<var<<" ["<<k<<"-"<<kk<<"], "<<endl; 
673    hist[kk]++;
674  }
675
676  cnts= 0;
677  for (int i = 0; i < ihist; ++i) {
678    if(hist.at(i) > 0) cnts++;
679    cout<<"("<<xcor.at(i)<<","<<hist.at(i)<<") ";
680  }
681  cout<< " ... cnts=  " << cnts << endl;
682}
Note: See TracBrowser for help on using the repository browser.