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

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

complements graphiques, legendes et unification unites
File size: 19.1 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  goodPartic_.clear();
22  rij_.raz();
23  P0Transport_ = 0.0;
24  particleRepresentationOk_ = false;
25  momentRepresentationOk_ = false;
26}
27
28int particleBeam::getNbParticles() const {
29  return goodPartic_.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  goodPartic_.clear();
102  goodPartic_ = particles;
103  cout << " particleBeam::setWithParticles taille vect. part. ENREGISTRE " << goodPartic_.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  goodPartic_.push_back(p);
117}
118
119const vector<bareParticle>& particleBeam::getParticleVector() const
120{
121  return goodPartic_;
122}
123
124vector<bareParticle>& particleBeam::getParticleVector() 
125{
126  return goodPartic_;
127}
128
129void 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 <<  goodPartic_.size() << " particules " << endl;
166  unsigned int k;
167  for ( k = 0 ; k < goodPartic_.size(); k++)
168    {
169      double xx,yy,zz;
170      goodPartic_.at(k).getPosition().getComponents(xx,yy,zz);
171      double betgamx, betgamy, betgamz;
172      goodPartic_.at(k).getBetaGamma().getComponents(betgamx, betgamy, betgamz);
173      cout << " part. numero " << k << "  x= " << xx << " y= " << yy  << " z= " << zz << "  betgamx= " << betgamx << " betgamy= " << betgamy  << " betgamz= " << betgamz << endl;
174    }
175}
176
177
178
179void particleBeam::Zrange(double& zmin, double& zmax) const {
180  double z;
181  zmin = GRAND;
182  zmax = -zmin;
183
184  unsigned int k;
185  for ( k = 0 ; k < goodPartic_.size(); k++)
186    {
187      z = goodPartic_.at(k).getZ();
188      if ( z < zmin ) zmin = z;
189      else if ( z > zmax) zmax = z;         
190    }
191}
192
193
194
195string particleBeam::FileOutputFlow() const {
196  ostringstream sortie;
197  unsigned int k;
198  for ( k = 0 ; k < goodPartic_.size(); k++)
199    {
200      sortie << goodPartic_.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
222  unsigned k,j,m;
223  double auxj, auxm;
224  if ( !particleRepresentationOk_)
225    {
226      cerr << " particleBeam::buildMomentRepresentation() vecteur de particules invalide" << endl;
227      return;
228    }
229
230  cout << " buildMomentRepresentation " << endl;
231  //  printAllXYZ();
232
233  double gref = referenceParticle_.getGamma() - 1.0;
234  double P_reference_MeV_sur_c = sqrt( gref*(gref+2) );
235
236  cout << " gref = " << gref << " P_reference_MeV_sur_c = " << P_reference_MeV_sur_c << endl;
237
238
239  // initialisation des moments
240  razDesMoments();
241
242  // accumulation
243  TRIDVECTOR pos;
244  TRIDVECTOR begam;
245  double gamma;
246  double begamz;
247  double g;
248  double PMeVsc;
249  double del;
250  vector<double> part(6, 0.0);
251
252    vector< vector<double> >& matrice = rij_.getMatrix();
253
254
255  for (k=0; k < goodPartic_.size(); k++) {
256    gamma = goodPartic_.at(k).getGamma();
257    pos = goodPartic_.at(k).getPosition();
258    begam= goodPartic_.at(k).getBetaGamma();
259    begamz = begam.getComponent(2);
260    g = gamma -1.0;
261    PMeVsc = sqrt( g*(g+2) );
262    del = 100.0 * ( PMeVsc -  P_reference_MeV_sur_c ) / P_reference_MeV_sur_c ; // en %
263
264    part[0] = pos.getComponent(0);
265    part[1] = begam.getComponent(0)/begamz;
266    part[2] = pos.getComponent(1);
267    part[3] = begam.getComponent(1)/begamz;
268    part[4] = pos.getComponent(2);
269    part[5] = del;
270
271    for ( j = 0; j < 6; j++) {
272      auxj = part.at(j) - centroid_.at(j);
273      for (m=0; m <= j; m++) 
274        {
275          auxm = part.at(m) - centroid_.at(m);
276
277          ( matrice.at(j) ).at(m) += auxj*auxm;
278          //      ( rij_transportMoments_.at(j) ).at(m) += auxj*auxm;
279
280
281          //          cout << " j= " << j << " m= " << m << " rjm= " << ( rij_transportMoments_.at(j) ).at(m) << endl;
282        }
283    }
284  }
285
286
287  // moyenne
288  double facmoy = 1.0/double( goodPartic_.size() );
289  for ( j = 0; j < 6; j++) {
290        ( matrice.at(j) ).at(j) = sqrt(( matrice.at(j) ).at(j) * facmoy );
291  }
292
293  for ( j = 0; j < 6; j++) {
294    auxj =  ( matrice.at(j) ).at(j);
295    for (m=0; m < j; m++) {
296      auxm = ( matrice.at(m) ).at(m);
297      (  matrice.at(j) ).at(m) *= facmoy/(auxj * auxm);
298    }
299  }
300   
301  ////////////////// si C21 = 1 , transport plante ! a voir //////////
302cout << " valeur initiale de  C21: " << ( matrice.at(1) ).at(0) << endl;
303  if ( ( matrice.at(1) ).at(0) >0.999999  ) {
304    ( matrice.at(1) ).at(0) = 0.999999;
305    cout << " j'ai fait la correction C21: " << ( matrice.at(1) ).at(0) << endl;
306  }
307 
308
309  // les longueurs sont en cm
310  // les angles en radians, on passe en mrad;
311
312  double uniteAngle = 1.0e+3;
313  ( matrice.at(1) ).at(1)  *= uniteAngle;
314  ( matrice.at(3) ).at(3)  *= uniteAngle;
315  P0Transport_ = 1.0e-3*EREST_MeV*P_reference_MeV_sur_c;
316
317  //  cout << " buildmomentrepresentation impression des moments " << endl;
318  //  impressionDesMoments();
319
320  momentRepresentationOk_ = true;
321}
322
323void particleBeam::impressionDesMoments() const {
324  rij_.impression();
325}
326
327void particleBeam::razDesMoments() {
328  rij_.raz();
329}
330
331
332// void particleBeam::readTransportMoments(ifstream& inp) {
333// rij_.readFromTransportOutput(inp);
334// }
335
336// void particleBeam::readTransportMoments(stringstream& inp) {
337// rij_.readFromTransportOutput(inp);
338// }
339
340double particleBeam::getXmaxRms() {
341  if ( !momentRepresentationOk_ ) buildMomentRepresentation();
342  return ( rij_.getMatrix().at(0) ).at(0);
343  // return ( rij_transportMoments_.at(0) ).at(0);
344}
345
346
347void particleBeam::particlesPhaseSpaceData(vector<double>& xcor, vector<double>& ycor, vector<string>& legende, string namex, string namey) {
348
349  unsigned indexAbs, indexOrd;
350  indexAbs = pspaCoorIndexFromString(namex);
351  indexOrd = pspaCoorIndexFromString(namey);
352 
353  particlesPhaseSpaceComponent(xcor, indexAbs);
354  particlesPhaseSpaceComponent(ycor, indexOrd);
355  legende.clear();
356    legende.push_back( "phase space " + namex + "," + namey);
357    legende.push_back( "particle number : " +   mixedTools::intToString(getNbParticles()));
358}
359
360void particleBeam::particlesPhaseSpaceComponent(vector<double>& coord, unsigned index) 
361{
362  if ( !particleRepresentationOk_ ) return;
363 
364  coord.clear();
365  coord.resize(goodPartic_.size(), 0.0 );
366  cout << " particleBeam::particlesPhaseSpaceComponent index = " << index << endl;
367
368  if ( index <= 2 ) {
369    for (unsigned i = 0; i < goodPartic_.size(); ++i) {
370      coord.at(i) =  10.*goodPartic_.at(i).getPosition().getComponent(index);  // en mm
371    }
372    return;
373  }
374 
375  if ( index >  2 && index < 5 ) {
376    for (unsigned i = 0; i < goodPartic_.size(); ++i) {
377      double begamz = goodPartic_.at(i).getBetaGamma().getComponent(2);
378      if ( begamz != 0.0) {
379        coord.at(i) =  1000.*goodPartic_.at(i).getBetaGamma().getComponent(index - 3)/begamz; // milliradians
380      } else {
381        coord.at(i) = 0.0;
382      }
383    }
384    return;
385  }
386
387  if ( index == 5 ) {
388    double gamma0 = referenceParticle_.getGamma();
389    if ( gamma0 == 0.0 ) return;
390    for (unsigned i = 0; i < goodPartic_.size(); ++i) {
391      coord.at(i) =  100.*(goodPartic_.at(i).getGamma() - gamma0)/gamma0;  // en %
392    }
393    return;
394  }
395}
396
397unsigned particleBeam::indexFromPspaToTransport(unsigned index) const {
398  cout << " indexFromPspaToTransport entree : " << index << endl;
399  switch ( index ) {
400  case 0 : return  0; // x
401  case 1 : return  2;  // y
402  case 2 : return  4; // z -> l
403  case 3 : return  1;  // xp
404  case 4 : return  3;  // yp
405  case 5 : return  5; // de/E
406  default : { 
407    cout << " particleBeam::indexFromPspaToTransport : coordinate index ERROR inital index :  "<< index << endl;
408    return 99;
409  }
410  }
411}
412
413unsigned particleBeam::pspaCoorIndexFromString(string nameCoor) const {
414  if ( nameCoor == "x" ) {
415    return 0;
416  } else if ( nameCoor == "y" ) {
417    return 1;
418  } else if ( nameCoor == "dz" ) {
419    return 2;
420  } else if ( nameCoor == "xp" ) {
421    return 3;
422  } else if ( nameCoor == "yp" ) {
423    return 4;
424  } else if ( nameCoor == "dE/E" ) {
425    return 5;
426  } else {
427    return 99;
428  }
429}
430
431unsigned particleBeam::transportCoorIndexFromString(string nameCoor) const {
432  if ( nameCoor == "x" ) {
433    return 0;
434  } else if ( nameCoor == "y" ) {
435    return 2;
436  } else if ( nameCoor == "dz" ) {
437    return 4;
438  } else if ( nameCoor == "xp" ) {
439    return 1;
440  } else if ( nameCoor == "yp" ) {
441    return 3;
442  } else if ( nameCoor == "dE/E" ) {
443    return 5;
444  } else {
445    return 99;
446  }
447}
448
449//void particleBeam::donneesDessinEllipse(vector<double>& xcor, vector<double>& ycor, vector<string>& legende, unsigned indexAbs, unsigned indexOrd) {
450
451void particleBeam::donneesDessinEllipse(vector<double>& xcor, vector<double>& ycor, vector<string>& legende, string namex, string namey) {
452  int k;
453  double x,y;
454  if ( !momentRepresentationOk_ ) buildMomentRepresentation();
455
456  //  if ( !momentRepresentationOk_ ) return;
457
458  unsigned indexAbs, indexOrd;
459
460
461  // index TRANSPORT
462  indexAbs = transportCoorIndexFromString(namex);
463  indexOrd = transportCoorIndexFromString(namey);
464  cout << " namex= " << namex << " indexAbs= " << indexAbs << " namey= " << namey << " indexOrd= " << indexOrd << endl; 
465    if ( indexAbs > 5 || indexOrd > 5 ) return;
466
467  xcor.clear();
468  ycor.clear();
469
470
471
472  double dimensionalFactorX, dimensionalFactorY; // to mm, if necessary
473  dimensionalFactorX = dimensionalFactorFromTransportToGraphics(indexAbs);
474  dimensionalFactorY = dimensionalFactorFromTransportToGraphics(indexOrd);
475
476
477
478  legende.clear();
479  //  if ( indexAbs == 0 && indexOrd == 1 ) {
480    string namx = transportVariableName(indexAbs);
481    string namy = transportVariableName(indexOrd);
482    double em = getUnnormalizedTranspPhaseSpaceArea(indexAbs,indexOrd) ;
483    string  emitt = mixedTools::doubleToString(getUnnormalizedTranspPhaseSpaceArea(indexAbs,indexOrd));
484    string xmax = namx + "max= ";
485    string valXmax = mixedTools::doubleToString(dimensionalFactorX*getSigmaTransportij(indexAbs,indexAbs)); 
486    string ymax = namy + "max= ";
487    string valYmax = mixedTools::doubleToString(dimensionalFactorY*getSigmaTransportij(indexOrd,indexOrd)); // mm
488    string correl = " correlation ";
489    string valCorrel = mixedTools::doubleToString(getSigmaTransportij(1,0));
490    string xunit = graphicTransportUnitName(indexAbs);
491    string yunit = graphicTransportUnitName(indexOrd);
492    legende.push_back( "emittance" + namx + "," + namy + ": " + emitt + " pi." + xunit + "." + yunit);
493    legende.push_back( xmax + valXmax + xunit);
494    legende.push_back( ymax + valYmax + yunit);
495  // } else {
496  //   legende.push_back(" text of legend not yet programmed ");
497  // }
498
499  cout << " index x" << indexAbs << " index y " << indexOrd << endl;
500  double xm = dimensionalFactorX*( rij_.getMatrix().at(indexAbs) ).at(indexAbs);
501  double ym = dimensionalFactorY*( rij_.getMatrix().at(indexOrd) ).at(indexOrd);
502  double r;
503  if ( indexOrd > indexAbs ) {
504    r  = ( rij_.getMatrix().at(indexOrd) ).at(indexAbs);
505  } else {
506    r  = ( rij_.getMatrix().at(indexAbs) ).at(indexOrd);
507  }
508
509  cout <<  " racs11= " << xm << " racs22= " << ym << " r12= " << r << endl;
510
511
512  int nbintv = 50;
513  if ( xm == 0.0 ) return;
514  double pas = 2.0 * xm / nbintv;
515
516  //  cout << " r= " << r << endl;
517  double rac = (1 - r*r);
518  if ( rac > 0.0 ) 
519    {
520      cout << " cas rac > " << endl;
521      rac = sqrt(1 - r*r);
522      double alpha = -r / rac;
523      double beta = xm / ( ym * rac);
524      //  double gamma = ym / ( xm * rac );
525      double epsil = xm * ym * rac;
526      double fac1 = -1.0 / ( beta * beta);
527      double fac2 = epsil/beta;
528      double fac3 = -alpha/beta;
529      double aux;
530      for ( k=0; k < nbintv; k++)
531        {
532          x = -xm + k*pas;
533          aux = fac1 * x * x + fac2;
534          //     cout << " aux2= " << aux << endl;
535          if ( aux <= 0.0 )
536            {
537              aux = 0.0;
538            }
539          else aux = sqrt(aux);
540     
541          //        y = fac3*x;
542          y = fac3*x + aux;
543          xcor.push_back(x);
544          ycor.push_back(y);
545        }
546
547      for ( k=0; k <= nbintv; k++)
548        {
549          x = xm - k*pas;
550          aux =  fac1 * x * x + fac2;
551          if ( aux <= 0.0 ) 
552            {
553              aux = 0.0;
554            }
555          else aux = sqrt(aux);
556          //   y = fac3*x;
557          y = fac3*x - aux;
558          xcor.push_back(x);
559          ycor.push_back(y);
560        }
561    }
562  else
563    // cas degenere
564    {
565      cout << " cas degenere " << endl;
566      double fac = ym/xm;
567      for ( k=0; k < nbintv; k++)
568        {
569          x = -xm + k*pas;
570          y = fac*x;
571          xcor.push_back(x);
572          ycor.push_back(y);
573        }
574       
575    }
576}
577
578void particleBeam::histogramme(unsigned int iabs,vector<double>& xcor,vector<int>& hist,double out[3])
579{
580  // out[0]= nbre de particules, out[1]= moyenne, out[2]= ecart-type
581  vector<double> vshf;
582  particlesPhaseSpaceComponent(vshf,iabs); 
583  histogramInitialize(iabs,vshf,out);
584  histogramPacking(out[2],vshf,xcor,hist);
585
586  if(iabs == 5) {
587    out[1] *= EREST_MeV; // moyenne en MeV
588    out[2] *= EREST_keV; // ecart-type en KeV
589  }
590}
591
592void particleBeam::histogramInitialize(unsigned int index,vector<double>& vshf,double out[3])
593{
594  double vmin= GRAND;
595  double vmax= -vmin;
596  double vmoy= 0.0;
597  double ecatyp= 0.0;
598 
599  for (unsigned int k = 0; k < goodPartic_.size(); k++) {
600    if (vshf[k] < vmin) vmin = vshf[k];
601    else if (vshf[k] > vmax) vmax = vshf[k];
602    vmoy += vshf[k];
603    ecatyp += vshf[k]*vshf[k];
604  }
605
606  double sum= (float)goodPartic_.size();
607  out[0]= sum; 
608  vmoy /= sum;
609  out[1]= vmoy;
610  ecatyp /= sum;
611  ecatyp = sqrt(abs(ecatyp-vmoy*vmoy));
612  out[2]= ecatyp;
613
614  if(index == 0) {
615    cout << "position x -moyenne " << vmoy << " cm " << "-mini " << vmin << " cm " << "-maxi " << vmax << " cm " << "ecart type " << ecatyp << " cm" << endl;
616  } 
617  if(index == 1) {
618    cout << "position y -moyenne " << vmoy << " cm " << "-mini " << vmin << " cm " << "-maxi " << vmax << " cm " << "ecart type " << ecatyp << " cm" << endl;
619  } 
620  if(index == 2) {
621    cout << "position z -moyenne " << vmoy << " cm " << "-mini " << vmin << " cm " << "-maxi " << vmax << " cm " << "ecart type " << ecatyp << " cm" << endl;
622  }
623  if(index == 3) {
624    cout << "divergence xp -moyenne " << vmoy << " mrad " << "-mini " << vmin << " mrad " << "-maxi " << vmax << " mrad " << "ecart type " << ecatyp << " mrad" << endl;
625  } 
626  if(index == 4) {
627    cout << "divergence yp -moyenne " << vmoy << " mrad " << "-mini " << vmin << " mrad " << "-maxi " << vmax << " mrad " << "ecart type " << ecatyp << " mrad" << endl;
628  } 
629  if(index == 5) {
630    double gmin = vmin-1.0;
631    double gmax = vmax-1.0;
632    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;
633  }
634
635  for (unsigned int k = 0; k < goodPartic_.size(); k++) {
636    vshf[k] -= vmoy;
637  }
638}
639
640void particleBeam::histogramPacking(double ecatyp,vector<double>vshf,vector<double>&xcor,vector<int>& hist)
641{
642  // demi fenetre hfene= max(3.*ecatyp-vmoy,vmoy-3.*ecatyp);
643  double hfene= 3.*ecatyp;
644  // et pas de l'histogramme
645  double hpas = hfene/25.;
646 
647  cout << "demi fenetre " << hfene << ", hpas= " << hpas << endl;
648
649  double vmin = -hfene;
650  double dfen = 2.*hfene;
651  int ihist = dfen/hpas;
652  double phist = ihist*hpas;
653  double dpas = hpas-(dfen-phist);
654  if(dpas <= hpas*1.e-03) {
655    ihist++;
656    phist= ihist*hpas;
657  }
658  double vmax= vmin+hpas*ihist;
659 
660  cout << "xAxisNumberOfBins= " << ihist <<", xAxisMinimum= " << vmin << ", xAxisMaximum= " << vmax << ", NParticules= " << vshf.size() << ", phist " << phist << endl;
661 
662  if(!xcor.empty()) xcor.clear();
663  xcor.resize(ihist+1); 
664  for (int i = 0; i <= ihist; ++i) {
665    xcor[i] = vmin+i*hpas;
666  }
667
668  /////////////////////////////////////
669
670  if(!hist.empty()) hist.clear();
671  hist.resize(ihist,0); 
672  for (unsigned int i = 0; i < vshf.size(); ++i) {
673    double var= vshf[i]-vmin;
674    if(var < 0 ) {
675      cout<<"lesser that the minimum "<<var<<", ("<< i<<")"<< endl;
676      hist.at(0)++;
677    } else if(var >= phist) {
678      cout<<"greater that the maximum "<<var<<", ("<< i<<")"<< endl;
679      hist.at(ihist-1)++;
680    } else {
681      int kk= (int)floor(var/hpas);
682      hist.at(kk)++;
683    }
684  }
685}
686
Note: See TracBrowser for help on using the repository browser.