source: PSPA/Interface_Web/trunk/pspaWT/sources/controler/src/particleBeam.cc @ 354

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

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