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

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

correction softwaregenerator
File size: 22.4 KB
Line 
1
2#include "particleBeam.h"
3#include "mathematicalConstants.h"
4#include "PhysicalConstants.h"
5#include "mathematicalTools.h"
6#include "mixedTools.h"
7
8#include <stdio.h>
9#include <algorithm>
10#include <sstream>
11
12using namespace std;
13
14particleBeam::particleBeam()  {
15  P0Transport_ = 0.0;
16  particleRepresentationOk_ = false;
17  momentRepresentationOk_ = false;
18}
19
20void particleBeam::clear() {
21  relativePartic_.clear();
22  rij_.raz();
23  P0Transport_ = 0.0;
24  particleRepresentationOk_ = false;
25  momentRepresentationOk_ = false;
26}
27
28int particleBeam::getNbParticles() const {
29  return relativePartic_.size();
30}
31
32const beam2Moments& particleBeam::getTransportMoments() const  { 
33  return rij_;
34}
35
36double particleBeam::getSigmaTransportij(unsigned indexI, unsigned indexJ)  {
37  if (  indexI > 5 ||  indexJ > 5 ) {
38    cerr << " particleBeam::getSigmaTransportij() indices out of range  " << endl;
39    return 0.0;
40  }
41  if ( !momentRepresentationOk_ ) {
42    cerr << " particleBeam::getSigmaTransportij() beam is not in moment representation " << endl;
43    return 0.0;
44  }
45 if ( indexI >= indexJ ) {
46   return ( rij_.getMatrix().at(indexI) ).at(indexJ);
47 } else {
48   return ( rij_.getMatrix().at(indexJ) ).at(indexI);
49 }
50 
51}
52
53double particleBeam::getUnnormalizedEmittanceX() {
54  double r = getSigmaTransportij(1,0);
55  double rac = (1 - r*r);
56  if ( rac <= 0.0 ) {
57    return 0.0;
58  }
59  rac = sqrt(1 - r*r);
60  return  dimensionalFactorFromTransportToGraphics(0)*getSigmaTransportij(0,0) * getSigmaTransportij(1,1) * rac; // en pi.mm.mrad
61}
62
63double particleBeam::getUnnormalizedTranspPhaseSpaceArea(unsigned TranspIndexAbs, unsigned TranspIndexOrd) {
64  if (  TranspIndexAbs == TranspIndexOrd ) return 0.0;
65  double r = getSigmaTransportij(TranspIndexAbs,TranspIndexOrd);
66  double rac = (1 - r*r);
67  if ( rac <= 0.0 ) {
68    return 0.0;
69  }
70  rac = sqrt(1 - r*r);
71  return  dimensionalFactorFromTransportToGraphics(TranspIndexAbs)*getSigmaTransportij(TranspIndexAbs,TranspIndexAbs) * 
72             dimensionalFactorFromTransportToGraphics(TranspIndexOrd)*getSigmaTransportij(TranspIndexOrd,TranspIndexOrd) * rac; // en pi.mm.mrad
73}
74
75
76
77double particleBeam::getP0Transport() const { 
78  return P0Transport_;
79}
80
81double particleBeam::referenceKineticEnergyMeV() const {
82  if ( particleRepresentationOk_ ) {
83    return (referenceParticle_.getGamma() -1.) * EREST_MeV;
84  } else {
85    double P0Norm = 1000.0 * P0Transport_ / EREST_MeV;
86    double gamma = sqrt(1.0 +  P0Norm * P0Norm);
87    return (gamma - 1.0) * EREST_MeV;
88  }
89}
90
91void particleBeam::set2Moments(beam2Moments& moments) {
92  rij_ = moments;
93  momentRepresentationOk_ = true;
94}
95
96void particleBeam::setWithParticles(vector<double>& centroid, bareParticle& referencePart, vector<bareParticle>& particles) {
97  cout << " particleBeam::setWithParticles taille vect. part. " << particles.size() << endl;
98  centroid_.clear();
99  centroid_ = centroid;
100  referenceParticle_ = referencePart;
101  relativePartic_.clear();
102  relativePartic_ = particles;
103  cout << " particleBeam::setWithParticles taille vect. part. ENREGISTRE " << relativePartic_.size() << endl;
104  particleRepresentationOk_ = true;
105}
106
107bool particleBeam::particleRepresentationOk() const {
108  return particleRepresentationOk_;
109}
110bool particleBeam::momentRepresentationOk() const {
111  return momentRepresentationOk_;
112}
113
114void  particleBeam::addParticle( bareParticle p)
115{
116  relativePartic_.push_back(p);
117}
118
119const vector<bareParticle>& particleBeam::getParticleVector() const
120{
121  return relativePartic_;
122}
123
124vector<bareParticle>& particleBeam::getParticleVector() 
125{
126  return relativePartic_;
127}
128
129// void particleBeam::getVariance(double& varx, double& vary, double& varz) const {
130//   unsigned int k;
131//   double x,y,z;
132//   double xav = 0.;
133//   double yav = 0.;
134//   double zav = 0.;
135//   double xavsq = 0.;
136//   double yavsq = 0.;
137//   double zavsq = 0.;
138
139//   TRIDVECTOR pos;
140
141
142//   for ( k = 0 ; k < goodPartic_.size(); k++) {
143//     pos = goodPartic_.at(k).getPosition();
144//     pos.getComponents(x,y,z);
145//     //      partic_[k].getXYZ(x,y,z);
146//     xav += x;
147//     xavsq += x*x;
148//     yav += y;
149//     yavsq += y*y;
150//     zav += z;
151//     zavsq += z*z;
152//   }
153
154//   double aginv = double (goodPartic_.size());
155//   aginv = 1.0/aginv;
156
157//   varx =  aginv * ( xavsq - xav*xav*aginv );
158//   vary =  aginv * ( yavsq - yav*yav*aginv );
159//   varz =  aginv * ( zavsq - zav*zav*aginv );
160// }
161
162
163void particleBeam::printAllXYZ() const {
164  cout << " dump du faisceau : " << endl;
165  cout <<  relativePartic_.size() << " particules " << endl;
166  unsigned int k;
167  for ( k = 0 ; k < relativePartic_.size(); k++)
168    {
169      double xx,yy,zz;
170      relativePartic_.at(k).getPosition().getComponents(xx,yy,zz);
171      double betgamx, betgamy, betgamz;
172      relativePartic_.at(k).getBetaGamma().getComponents(betgamx, betgamy, betgamz);
173      cout << " part. numero " << k << "  x= " << xx << " y= " << yy  << " dphas (c.dt, cm) = " << zz << "  betgamx= " << betgamx << " betgamy= " << betgamy  << " betgamz= " << betgamz << endl;
174    }
175}
176
177
178// extension en phase (cm)
179void particleBeam::ZrangeCdt(double& cdtmin, double& cdtmax) const {
180  double z;
181  cdtmin = GRAND;
182  cdtmax = -cdtmin;
183
184  unsigned int k;
185  for ( k = 0 ; k < relativePartic_.size(); k++)
186    {
187      z = relativePartic_.at(k).getZ(); // ce z est un dephasage, c.dt, en cm
188      if ( z < cdtmin ) cdtmin = z;
189      else if ( z > cdtmax) cdtmax = z;     
190    }
191}
192
193
194
195string particleBeam::fileOutputFlow() const {
196  ostringstream sortie;
197  unsigned int k;
198  for ( k = 0 ; k < relativePartic_.size(); k++)
199    {
200      sortie << relativePartic_.at(k).FileOutputFlow() << endl;
201    }
202  sortie << endl;
203  return sortie.str();
204}
205
206bool particleBeam::FileInput( ifstream& ifs) {
207  bool test = true;
208  string dum1, dum2;
209  double dummy;
210  if ( !( ifs >> dum1 >> dum2 >> dummy) ) return false;
211 
212  bareParticle pp;
213  while ( pp.FileInput(ifs) )
214    {
215      addParticle( pp);
216    }
217  return test;
218}
219
220void particleBeam::buildMomentRepresentation() {
221  // le faisceau cense etre represente en un z donne (reference) sera
222  // ici deploye spatialement (faisceau en un temps donne)
223  unsigned k,j,m;
224  double auxj, auxm;
225  if ( !particleRepresentationOk_)
226    {
227      cerr << " particleBeam::buildMomentRepresentation() vecteur de particules invalide" << endl;
228      return;
229    }
230
231  cout << " buildMomentRepresentation " << endl;
232  //  printAllXYZ();
233
234  double gref = referenceParticle_.getGamma() - 1.0;
235  double P_reference_MeV_sur_c = sqrt( gref*(gref+2) );
236
237  cout << " gref = " << gref << " P_reference_MeV_sur_c = " << P_reference_MeV_sur_c << endl;
238
239
240  // initialisation des moments
241  razDesMoments();
242
243  // accumulation
244  TRIDVECTOR pos;
245  TRIDVECTOR begam;
246  double gamma;
247  double begamz;
248  double g;
249  double PMeVsc;
250  double del;
251  //  double xp,yp,dz,cdt;
252  vector<double> part(6, 0.0);
253
254    vector< vector<double> >& matrice = rij_.getMatrix();
255
256    TRIDVECTOR positionDeployee;
257
258  for (k=0; k < relativePartic_.size(); k++) {
259
260    positionDeployee = coordonneesDeployees(k);
261    gamma = relativePartic_.at(k).getGamma();
262
263    //    pos = relativePartic_.at(k).getPosition();
264    //    cdt = pos.getComponent(2);
265    begam= relativePartic_.at(k).getBetaGamma();
266    begamz = begam.getComponent(2);
267    g = gamma -1.0;
268    PMeVsc = sqrt( g*(g+2) );
269    del = 100.0 * ( PMeVsc -  P_reference_MeV_sur_c ) / P_reference_MeV_sur_c ; // en %
270
271    //    dz = begamz * cdt / gamma;
272    //    xp = begam.getComponent(0)/begamz;
273    //    yp = begam.getComponent(1)/begamz;
274
275    part[0] = positionDeployee.getComponent(0);
276    part[1] = begam.getComponent(0)/begamz;
277    part[2] = positionDeployee.getComponent(1);
278    part[3] = begam.getComponent(1)/begamz;
279    part[4] = positionDeployee.getComponent(2);
280    part[5] = del;
281
282    for ( j = 0; j < 6; j++) {
283      auxj = part.at(j) - centroid_.at(j);
284      for (m=0; m <= j; m++) 
285        {
286          auxm = part.at(m) - centroid_.at(m);
287
288          ( matrice.at(j) ).at(m) += auxj*auxm;
289          //      ( rij_transportMoments_.at(j) ).at(m) += auxj*auxm;
290
291
292          //          cout << " j= " << j << " m= " << m << " rjm= " << ( rij_transportMoments_.at(j) ).at(m) << endl;
293        }
294    }
295  }
296
297
298  // moyenne
299  double facmoy = 1.0/double( relativePartic_.size() );
300  for ( j = 0; j < 6; j++) {
301        ( matrice.at(j) ).at(j) = sqrt(( matrice.at(j) ).at(j) * facmoy );
302  }
303
304  for ( j = 0; j < 6; j++) {
305    auxj =  ( matrice.at(j) ).at(j);
306    for (m=0; m < j; m++) {
307      auxm = ( matrice.at(m) ).at(m);
308      (  matrice.at(j) ).at(m) *= facmoy/(auxj * auxm);
309    }
310  }
311   
312  ////////////////// si C21 = 1 , transport plante ! a voir //////////
313cout << " valeur initiale de  C21: " << ( matrice.at(1) ).at(0) << endl;
314  if ( ( matrice.at(1) ).at(0) >0.999999  ) {
315    ( matrice.at(1) ).at(0) = 0.999999;
316    cout << " j'ai fait la correction C21: " << ( matrice.at(1) ).at(0) << endl;
317  }
318 
319
320  // les longueurs sont en cm
321  // les angles en radians, on passe en mrad;
322
323  double uniteAngle = 1.0e+3;
324  ( matrice.at(1) ).at(1)  *= uniteAngle;
325  ( matrice.at(3) ).at(3)  *= uniteAngle;
326  P0Transport_ = 1.0e-3*EREST_MeV*P_reference_MeV_sur_c;
327
328  //  cout << " buildmomentrepresentation impression des moments " << endl;
329  //  impressionDesMoments();
330
331  momentRepresentationOk_ = true;
332}
333
334void particleBeam::impressionDesMoments() const {
335  rij_.impression();
336}
337
338void particleBeam::razDesMoments() {
339  rij_.raz();
340}
341
342
343// void particleBeam::readTransportMoments(ifstream& inp) {
344// rij_.readFromTransportOutput(inp);
345// }
346
347// void particleBeam::readTransportMoments(stringstream& inp) {
348// rij_.readFromTransportOutput(inp);
349// }
350
351double particleBeam::getXmaxRms() {
352  if ( !momentRepresentationOk_ ) buildMomentRepresentation();
353  return ( rij_.getMatrix().at(0) ).at(0);
354  // return ( rij_transportMoments_.at(0) ).at(0);
355}
356
357
358void particleBeam::particlesPhaseSpaceData(vector<double>& xcor, vector<double>& ycor, vector<string>& legende, string namex, string namey) {
359
360  unsigned indexAbs, indexOrd;
361  indexAbs = pspaCoorIndexFromString(namex);
362  indexOrd = pspaCoorIndexFromString(namey);
363 
364  particlesPhaseSpaceComponent(xcor, indexAbs);
365  particlesPhaseSpaceComponent(ycor, indexOrd);
366  legende.clear();
367    legende.push_back( "phase space " + namex + "," + namey);
368    legende.push_back( "particle number : " +   mixedTools::intToString(getNbParticles()));
369}
370
371// renvoie (dans le vecteur coord) la liste des coordonnées d'index 'index' :
372// index = 0 , 1, 2 -> x,y,z (en coordonnees deployees)
373// index = 3,4 -> x' = betax/betaz , y' = betay/betaz
374// index = 5 -> Ec : energie cinetique
375void particleBeam::particlesPhaseSpaceComponent(vector<double>& coord, unsigned index) 
376{
377
378  if ( !particleRepresentationOk_ ) return;
379 
380  coord.clear();
381  coord.resize(relativePartic_.size(), 0.0 );
382  //  cout << " particleBeam::particlesPhaseSpaceComponent index = " << index << endl;
383
384  if ( index <= 2 ) {
385
386    for (unsigned i = 0; i < relativePartic_.size(); ++i) {
387
388      coord.at(i) =  10.*coordonneesDeployees(i).getComponent(index);
389      //      coord.at(i) =  10.*relativePartic_.at(i).getPosition().getComponent(index);  // en mm
390    }
391    return;
392  }
393 
394  if ( index >  2 && index < 5 ) {
395    //    cout << " particleBeam::particlesPhaseSpaceComponent traitement vitesses " << endl;
396    for (unsigned i = 0; i < relativePartic_.size(); ++i) {
397      double begamz = relativePartic_.at(i).getBetaGamma().getComponent(2);
398      if ( begamz != 0.0) {
399        coord.at(i) =  1000.*relativePartic_.at(i).getBetaGamma().getComponent(index - 3)/begamz; // milliradians
400      } else {
401        coord.at(i) = 0.0;
402      }
403    }
404    //    cout << " particleBeam::particlesPhaseSpaceComponent traitement vitesses TERMINE " << endl;
405    return;
406  }
407
408  if ( index == 5 ) {
409    double gamma0 = referenceParticle_.getGamma();
410    if ( gamma0 == 0.0 ) return;
411    for (unsigned i = 0; i < relativePartic_.size(); ++i) {
412      coord.at(i) =  100.*(relativePartic_.at(i).getGamma() - gamma0)/gamma0;  // en %
413    }
414    return;
415  }
416}
417
418unsigned particleBeam::indexFromPspaToTransport(unsigned index) const {
419  cout << " indexFromPspaToTransport entree : " << index << endl;
420  switch ( index ) {
421  case 0 : return  0; // x
422  case 1 : return  2;  // y
423  case 2 : return  4; // z -> l
424  case 3 : return  1;  // xp
425  case 4 : return  3;  // yp
426  case 5 : return  5; // de/E
427  default : { 
428    cout << " particleBeam::indexFromPspaToTransport : coordinate index ERROR inital index :  "<< index << endl;
429    return 99;
430  }
431  }
432}
433
434unsigned particleBeam::pspaCoorIndexFromString(string nameCoor) const {
435  if ( nameCoor == "x" ) {
436    return 0;
437  } else if ( nameCoor == "y" ) {
438    return 1;
439  } else if ( nameCoor == "dz" ) {
440    return 2;
441  } else if ( nameCoor == "xp" ) {
442    return 3;
443  } else if ( nameCoor == "yp" ) {
444    return 4;
445  } else if ( nameCoor == "dE/E" ) {
446    return 5;
447  } else {
448    return 99;
449  }
450}
451
452unsigned particleBeam::transportCoorIndexFromString(string nameCoor) const {
453  if ( nameCoor == "x" ) {
454    return 0;
455  } else if ( nameCoor == "y" ) {
456    return 2;
457  } else if ( nameCoor == "dz" ) {
458    return 4;
459  } else if ( nameCoor == "xp" ) {
460    return 1;
461  } else if ( nameCoor == "yp" ) {
462    return 3;
463  } else if ( nameCoor == "dE/E" ) {
464    return 5;
465  } else {
466    return 99;
467  }
468}
469
470//void particleBeam::donneesDessinEllipse(vector<double>& xcor, vector<double>& ycor, vector<string>& legende, unsigned indexAbs, unsigned indexOrd) {
471
472void particleBeam::donneesDessinEllipse(vector<double>& xcor, vector<double>& ycor, vector<string>& legende, string namex, string namey) {
473  int k;
474  double x,y;
475  if ( !momentRepresentationOk_ ) buildMomentRepresentation();
476
477  //  if ( !momentRepresentationOk_ ) return;
478
479  unsigned indexAbs, indexOrd;
480
481
482  // index TRANSPORT
483  indexAbs = transportCoorIndexFromString(namex);
484  indexOrd = transportCoorIndexFromString(namey);
485  cout << " namex= " << namex << " indexAbs= " << indexAbs << " namey= " << namey << " indexOrd= " << indexOrd << endl; 
486    if ( indexAbs > 5 || indexOrd > 5 ) return;
487
488  xcor.clear();
489  ycor.clear();
490
491
492
493  double dimensionalFactorX, dimensionalFactorY; // to mm, if necessary
494  dimensionalFactorX = dimensionalFactorFromTransportToGraphics(indexAbs);
495  dimensionalFactorY = dimensionalFactorFromTransportToGraphics(indexOrd);
496
497
498
499  legende.clear();
500  //  if ( indexAbs == 0 && indexOrd == 1 ) {
501    string namx = transportVariableName(indexAbs);
502    string namy = transportVariableName(indexOrd);
503    double em = getUnnormalizedTranspPhaseSpaceArea(indexAbs,indexOrd) ;
504    string  emitt = mixedTools::doubleToString(getUnnormalizedTranspPhaseSpaceArea(indexAbs,indexOrd));
505    string xmax = namx + "max= ";
506    string valXmax = mixedTools::doubleToString(dimensionalFactorX*getSigmaTransportij(indexAbs,indexAbs)); 
507    string ymax = namy + "max= ";
508    string valYmax = mixedTools::doubleToString(dimensionalFactorY*getSigmaTransportij(indexOrd,indexOrd)); // mm
509    string correl = " correlation ";
510    string valCorrel = mixedTools::doubleToString(getSigmaTransportij(1,0));
511    string xunit = graphicTransportUnitName(indexAbs);
512    string yunit = graphicTransportUnitName(indexOrd);
513    legende.push_back( "emittance" + namx + "," + namy + ": " + emitt + " pi." + xunit + "." + yunit);
514    legende.push_back( xmax + valXmax + xunit);
515    legende.push_back( ymax + valYmax + yunit);
516  // } else {
517  //   legende.push_back(" text of legend not yet programmed ");
518  // }
519
520  cout << " index x" << indexAbs << " index y " << indexOrd << endl;
521  double xm = dimensionalFactorX*( rij_.getMatrix().at(indexAbs) ).at(indexAbs);
522  double ym = dimensionalFactorY*( rij_.getMatrix().at(indexOrd) ).at(indexOrd);
523  double r;
524  if ( indexOrd > indexAbs ) {
525    r  = ( rij_.getMatrix().at(indexOrd) ).at(indexAbs);
526  } else {
527    r  = ( rij_.getMatrix().at(indexAbs) ).at(indexOrd);
528  }
529
530  cout <<  " racs11= " << xm << " racs22= " << ym << " r12= " << r << endl;
531
532
533  int nbintv = 50;
534  if ( xm == 0.0 ) return;
535  double pas = 2.0 * xm / nbintv;
536
537  //  cout << " r= " << r << endl;
538  double rac = (1 - r*r);
539  if ( rac > 0.0 ) 
540    {
541      cout << " cas rac > " << endl;
542      rac = sqrt(1 - r*r);
543      double alpha = -r / rac;
544      double beta = xm / ( ym * rac);
545      //  double gamma = ym / ( xm * rac );
546      double epsil = xm * ym * rac;
547      double fac1 = -1.0 / ( beta * beta);
548      double fac2 = epsil/beta;
549      double fac3 = -alpha/beta;
550      double aux;
551      for ( k=0; k < nbintv; k++)
552        {
553          x = -xm + k*pas;
554          aux = fac1 * x * x + fac2;
555          //     cout << " aux2= " << aux << endl;
556          if ( aux <= 0.0 )
557            {
558              aux = 0.0;
559            }
560          else aux = sqrt(aux);
561     
562          //        y = fac3*x;
563          y = fac3*x + aux;
564          xcor.push_back(x);
565          ycor.push_back(y);
566        }
567
568      for ( k=0; k <= nbintv; k++)
569        {
570          x = xm - k*pas;
571          aux =  fac1 * x * x + fac2;
572          if ( aux <= 0.0 ) 
573            {
574              aux = 0.0;
575            }
576          else aux = sqrt(aux);
577          //   y = fac3*x;
578          y = fac3*x - aux;
579          xcor.push_back(x);
580          ycor.push_back(y);
581        }
582    }
583  else
584    // cas degenere
585    {
586      cout << " cas degenere " << endl;
587      double fac = ym/xm;
588      for ( k=0; k < nbintv; k++)
589        {
590          x = -xm + k*pas;
591          y = fac*x;
592          xcor.push_back(x);
593          ycor.push_back(y);
594        }
595       
596    }
597}
598
599void particleBeam::histogramme(unsigned int iabs,vector<double>& xcor,vector<int>& hist, vector<string>& legende)
600{
601  double out[3];
602  // out[0]= nbre de particules, out[1]= moyenne, out[2]= ecart-type
603  vector<double> vshf;
604  particlesPhaseSpaceComponent(vshf,iabs); 
605  histogramInitialize(iabs,vshf,out);
606  histogramPacking(out[2],vshf,xcor,hist);
607
608  if(iabs == 5) {
609    out[1] *= EREST_MeV; // moyenne en MeV
610    out[2] *= EREST_keV; // ecart-type en KeV
611  }
612  // legendes
613  string unites[2];
614  if(iabs == 0 || iabs == 1 || iabs == 2) {
615    unites[0]= unites[1]= " mm";
616  }
617  if(iabs == 3 || iabs == 4) {
618    unites[0]= unites[1]= " mrad";
619  }
620  if(iabs == 5) {
621    unites[0]= " MeV"; unites[1]= " KeV";
622  }
623
624  legende.clear();
625  legende.push_back(" entries : "+ mixedTools::intToString((int)out[0]) );
626  legende.push_back(" mean : "+ mixedTools::doubleToString(out[1])+unites[0]);
627  legende.push_back(" sigma rms : "+ mixedTools::doubleToString(out[2])+unites[1]);
628}
629
630void particleBeam::histogramInitialize(unsigned int index,vector<double>& vshf,double out[3])
631{
632  double vmin= GRAND;
633  double vmax= -vmin;
634  double vmoy= 0.0;
635  double ecatyp= 0.0;
636 
637  for (unsigned int k = 0; k < relativePartic_.size(); k++) {
638    if (vshf[k] < vmin) vmin = vshf[k];
639    else if (vshf[k] > vmax) vmax = vshf[k];
640    vmoy += vshf[k];
641    ecatyp += vshf[k]*vshf[k];
642  }
643
644  double sum= (float)relativePartic_.size();
645  out[0]= sum; 
646  vmoy /= sum;
647  out[1]= vmoy;
648  ecatyp /= sum;
649  ecatyp = sqrt(abs(ecatyp-vmoy*vmoy));
650  out[2]= ecatyp;
651
652  if(index == 0) {
653    cout << "position x -moyenne " << vmoy << " cm " << "-mini " << vmin << " cm " << "-maxi " << vmax << " cm " << "ecart type " << ecatyp << " cm" << endl;
654  } 
655  if(index == 1) {
656    cout << "position y -moyenne " << vmoy << " cm " << "-mini " << vmin << " cm " << "-maxi " << vmax << " cm " << "ecart type " << ecatyp << " cm" << endl;
657  } 
658  if(index == 2) {
659    cout << "position z -moyenne " << vmoy << " cm " << "-mini " << vmin << " cm " << "-maxi " << vmax << " cm " << "ecart type " << ecatyp << " cm" << endl;
660  }
661  if(index == 3) {
662    cout << "divergence xp -moyenne " << vmoy << " mrad " << "-mini " << vmin << " mrad " << "-maxi " << vmax << " mrad " << "ecart type " << ecatyp << " mrad" << endl;
663  } 
664  if(index == 4) {
665    cout << "divergence yp -moyenne " << vmoy << " mrad " << "-mini " << vmin << " mrad " << "-maxi " << vmax << " mrad " << "ecart type " << ecatyp << " mrad" << endl;
666  } 
667  if(index == 5) {
668    double gmin = vmin-1.0;
669    double gmax = vmax-1.0;
670    cout << "energie cinetique -moyenne " << vmoy*EREST_MeV << " Mev " << "-mini " << gmin*EREST_MeV << " Mev " << "-maxi " << gmax*EREST_MeV << " Mev " << "ecart type " << ecatyp*EREST_MeV << " Kev" << endl;
671  }
672
673  for (unsigned int k = 0; k < relativePartic_.size(); k++) {
674    vshf[k] -= vmoy;
675  }
676}
677
678void particleBeam::histogramPacking(double ecatyp,vector<double>vshf,vector<double>&xcor,vector<int>& hist)
679{
680  // demi fenetre hfene= max(5.*ecatyp-vmoy,vmoy-5.*ecatyp);
681  double hfene= 5.*ecatyp;
682  // et pas de l'histogramme
683  double hpas = hfene/25.;
684 
685  cout << "demi fenetre " << hfene << ", hpas= " << hpas << endl;
686
687  double vmin = -hfene;
688  double dfen = 2.*hfene;
689  int ihist = dfen/hpas;
690  double phist = ihist*hpas;
691  double dpas = hpas-(dfen-phist);
692  if(dpas <= hpas*1.e-03) {
693    ihist++;
694    phist= ihist*hpas;
695  }
696  double vmax= vmin+hpas*ihist;
697 
698  cout << "xAxisNumberOfBins= " << ihist <<", xAxisMinimum= " << vmin << ", xAxisMaximum= " << vmax << ", NParticules= " << vshf.size() << ", phist " << phist << endl;
699 
700  if(!xcor.empty()) xcor.clear();
701  xcor.resize(ihist+1); 
702  for (int i = 0; i <= ihist; ++i) {
703    xcor[i] = vmin+i*hpas;
704  }
705
706  /////////////////////////////////////
707
708  if(!hist.empty()) hist.clear();
709  hist.resize(ihist,0); 
710  for (unsigned int i = 0; i < vshf.size(); ++i) {
711    double var= vshf[i]-vmin;
712    if(var < 0 ) {
713      cout<<"lesser that the minimum "<<var<<", ("<< i<<")"<< endl;
714      hist.at(0)++;
715    } else if(var >= phist) {
716      cout<<"greater that the maximum "<<var<<", ("<< i<<")"<< endl;
717      hist.at(ihist-1)++;
718    } else {
719      int kk= (int)floor(var/hpas);
720      hist.at(kk)++;
721    }
722  }
723}
724
725 // coordonnees d'une particule dans le faisceau deploye ( passage
726 // d'une representation z=cte a une representation t=cte)
727TRIDVECTOR particleBeam::coordonneesDeployees(unsigned particleIndex, double cdtShift) {
728    double gamma = relativePartic_.at(particleIndex).getGamma();
729    //    TRIDVECTOR pos = relativePartic_.at(particleIndex).getPosition();
730    // cout << " ------------------ particleBeam::coordonneesDeployees ----- " << endl;
731    //   cout << " particleBeam::coordonneesDeployees " << relativePartic_.at(particleIndex).FileOutputFlow() << endl;
732    //   cout << " par reference x= " << relativePartic_.at(particleIndex).getReferenceToPosition().getComponent(0);
733    // cout << " ----------------------------------------------------- " << endl;
734
735
736
737    double xx = relativePartic_.at(particleIndex).getReferenceToPosition().getComponent(0);
738    double yy = relativePartic_.at(particleIndex).getReferenceToPosition().getComponent(1);
739    double cdt = relativePartic_.at(particleIndex).getReferenceToPosition().getComponent(2);
740    cdt += cdtShift;
741    //    TRIDVECTOR begam= goodPartic_.at(particleIndex).getBetaGamma();
742    double begamz = relativePartic_.at(particleIndex).getReferenceToBetaGamma().getComponent(2);
743    double fac = cdt / gamma;
744    //    double dz = begamz * fac;
745    xx += relativePartic_.at(particleIndex).getReferenceToBetaGamma().getComponent(0) * fac;
746    yy += relativePartic_.at(particleIndex).getReferenceToBetaGamma().getComponent(1) * fac;
747    return TRIDVECTOR(xx, yy, begamz * fac);
748 }
Note: See TracBrowser for help on using the repository browser.