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

Last change on this file since 371 was 371, checked in by touze, 12 years ago

nvx histos

File size: 19.2 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, unsigned indexAbs, unsigned indexOrd) {
348  particlesPhaseSpaceComponent(xcor, indexAbs);
349  particlesPhaseSpaceComponent(ycor, indexOrd);
350}
351
352void particleBeam::particlesPhaseSpaceComponent(vector<double>& coord, unsigned index) 
353{
354  if ( !particleRepresentationOk_ ) return;
355 
356  coord.clear();
357  coord.resize(goodPartic_.size(), 0.0 );
358  cout << " particleBeam::particlesPhaseSpaceComponent index = " << index << endl;
359  if ( index <= 2 ) {
360    for (unsigned i = 0; i < goodPartic_.size(); ++i) {
361      coord.at(i) =  goodPartic_.at(i).getPosition().getComponent(index);
362    }
363    return;
364  }
365
366  if ( index >  2 && index < 5 ) {
367    for (unsigned i = 0; i < goodPartic_.size(); ++i) {
368      double begamz = goodPartic_.at(i).getBetaGamma().getComponent(2);
369      if ( begamz != 0.0) {
370        coord.at(i) =  1000.*goodPartic_.at(i).getBetaGamma().getComponent(index - 3)/begamz; // mimmiradians
371      } else {
372        coord.at(i) = 0.0;
373      }
374    }
375    return;
376  }
377
378  if ( index == 5 ) {
379    double gamma0 = referenceParticle_.getGamma();
380    //    cout << " gamma0 = " << gamma0 << endl;
381    if ( gamma0 == 0.0 ) return;
382    for (unsigned i = 0; i < goodPartic_.size(); ++i) {
383      coord.at(i) =  100.*(goodPartic_.at(i).getGamma() - gamma0)/gamma0;  // en %
384      //      cout << " gamma0 = " << gamma0 << " gamma = " << goodPartic_.at(i).getGamma() << endl;
385    }
386    return;
387  }
388}
389
390unsigned particleBeam::indexFromPspaToTransport(unsigned index) const {
391  cout << " indexFromPspaToTransport entree : " << index << endl;
392  switch ( index ) {
393  case 0 : return  0; // x
394  case 1 : return  2;  // y
395  case 2 : return  4; // z -> l
396  case 3 : return  1;  // xp
397  case 4 : return  3;  // yp
398  case 5 : return  5; // de/E
399  default : { 
400    cout << " particleBeam::indexFromPspaToTransport : coordinate index ERROR inital index :  "<< index << endl;
401    return 99;
402  }
403  }
404}
405
406
407void particleBeam::donneesDessinEllipse(vector<double>& xcor, vector<double>& ycor, vector<string>& legende, unsigned indexAbs, unsigned indexOrd) {
408  int k;
409  double x,y;
410  cout << " particleBeam::donneesDessinEllipse index recus x" << indexAbs << " index y " << indexOrd << endl;
411
412  if ( !momentRepresentationOk_ ) return;
413
414  if ( indexAbs > 5 || indexOrd > 5 ) return;
415
416  xcor.clear();
417  ycor.clear();
418  // les index sont dans l'ordre x,y,z,xp,yp, de/E
419  // on traduit en TRANSPORT
420
421  indexAbs = indexFromPspaToTransport(indexAbs);
422  indexOrd = indexFromPspaToTransport(indexOrd);
423
424
425  double dimensionalFactorX, dimensionalFactorY; // to mm, if necessary
426  dimensionalFactorX = dimensionalFactorFromTransportToGraphics(indexAbs);
427  dimensionalFactorY = dimensionalFactorFromTransportToGraphics(indexOrd);
428
429
430  // a completer
431  legende.clear();
432  //  if ( indexAbs == 0 && indexOrd == 1 ) {
433    string namx = transportVariableName(indexAbs);
434    string namy = transportVariableName(indexOrd);
435    double em = getUnnormalizedTranspPhaseSpaceArea(indexAbs,indexOrd) ;
436    string  emitt = mixedTools::doubleToString(getUnnormalizedTranspPhaseSpaceArea(indexAbs,indexOrd));
437    string xmax = namx + "max= ";
438    string valXmax = mixedTools::doubleToString(dimensionalFactorX*getSigmaTransportij(indexAbs,indexAbs)); 
439    string ymax = namy + "max= ";
440    string valYmax = mixedTools::doubleToString(dimensionalFactorY*getSigmaTransportij(indexOrd,indexOrd)); // mm
441    string correl = " correlation ";
442    string valCorrel = mixedTools::doubleToString(getSigmaTransportij(1,0));
443    string xunit = graphicTransportUnitName(indexAbs);
444    string yunit = graphicTransportUnitName(indexOrd);
445    legende.push_back( "emittance" + namx + "," + namy + ": " + emitt + " pi." + xunit + "." + yunit);
446    legende.push_back( xmax + valXmax + xunit);
447    legende.push_back( ymax + valYmax + yunit);
448  // } else {
449  //   legende.push_back(" text of legend not yet programmed ");
450  // }
451
452  cout << " index x" << indexAbs << " index y " << indexOrd << endl;
453  double xm = dimensionalFactorX*( rij_.getMatrix().at(indexAbs) ).at(indexAbs);
454  double ym = dimensionalFactorY*( rij_.getMatrix().at(indexOrd) ).at(indexOrd);
455  double r;
456  if ( indexOrd > indexAbs ) {
457    r  = ( rij_.getMatrix().at(indexOrd) ).at(indexAbs);
458  } else {
459    r  = ( rij_.getMatrix().at(indexAbs) ).at(indexOrd);
460  }
461
462  cout <<  " racs11= " << xm << " racs22= " << ym << " r12= " << r << endl;
463
464
465  int nbintv = 50;
466  if ( xm == 0.0 ) return;
467  double pas = 2.0 * xm / nbintv;
468
469  //  cout << " r= " << r << endl;
470  double rac = (1 - r*r);
471  if ( rac > 0.0 ) 
472    {
473      cout << " cas rac > " << endl;
474      rac = sqrt(1 - r*r);
475      double alpha = -r / rac;
476      double beta = xm / ( ym * rac);
477      //  double gamma = ym / ( xm * rac );
478      double epsil = xm * ym * rac;
479      double fac1 = -1.0 / ( beta * beta);
480      double fac2 = epsil/beta;
481      double fac3 = -alpha/beta;
482      double aux;
483      for ( k=0; k < nbintv; k++)
484        {
485          x = -xm + k*pas;
486          aux = fac1 * x * x + fac2;
487          //     cout << " aux2= " << aux << endl;
488          if ( aux <= 0.0 )
489            {
490              aux = 0.0;
491            }
492          else aux = sqrt(aux);
493     
494          //        y = fac3*x;
495          y = fac3*x + aux;
496          xcor.push_back(x);
497          ycor.push_back(y);
498        }
499
500      for ( k=0; k <= nbintv; k++)
501        {
502          x = xm - k*pas;
503          aux =  fac1 * x * x + fac2;
504          if ( aux <= 0.0 ) 
505            {
506              aux = 0.0;
507            }
508          else aux = sqrt(aux);
509          //   y = fac3*x;
510          y = fac3*x - aux;
511          xcor.push_back(x);
512          ycor.push_back(y);
513        }
514    }
515  else
516    // cas degenere
517    {
518      cout << " cas degenere " << endl;
519      double fac = ym/xm;
520      for ( k=0; k < nbintv; k++)
521        {
522          x = -xm + k*pas;
523          y = fac*x;
524          xcor.push_back(x);
525          ycor.push_back(y);
526        }
527       
528    }
529}
530
531void particleBeam::histogramme(unsigned int iabs,vector<double>&xcor,vector<int>& hist,double out[3])
532{
533  vector<double> vshf;
534  histogramComponent(vshf,iabs);
535  histogramInitialize(iabs,vshf,out);
536  histogramPacking(out[2],vshf,xcor,hist);
537
538  // out[0]= nbre de particules, out[1]= moyenne, out[2]= ecart-type
539  if(iabs == 0 || iabs == 1 || iabs == 2) {
540    out[2] *= 1000.0; // ecart-type en mm
541  }
542
543  if(iabs == 5) {
544    for (unsigned int i = 0; i <= xcor.size(); ++i) {
545      xcor[i] = 100.*xcor[i]/out[1]; // en %
546    }
547    out[1] *= EREST_MeV; // moyenne en MeV
548    out[2] *= EREST_keV; // ecart-type en KeV
549  }
550}
551
552void particleBeam::histogramComponent(vector<double>& coord,unsigned int index) 
553{ 
554  coord.clear();
555  coord.resize(goodPartic_.size(),0.0);
556 
557  if(index == 0 || index == 1 || index == 2) 
558    {
559      // collecte d'une coordonnee x(index= 0), y(index= 1) ou z(index= 2) des particules
560      for (unsigned int k = 0; k < goodPartic_.size(); k++) {
561        coord.at(k) = goodPartic_.at(k).getPosition().getComponent(index);
562      }
563    }
564
565  if(index == 3 || index == 4)
566    {
567      // collecte de begamx (index= 3), begamy (index= 4)
568      for (unsigned int k = 0; k < goodPartic_.size(); k++) {
569        double begamz = goodPartic_.at(k).getBetaGamma().getComponent(2);
570        if (begamz == 0.0) continue;
571        double begam = goodPartic_.at(k).getBetaGamma().getComponent(index-3);
572        coord.at(k) = 1000.0*begam/begamz;  //pi.mm.mradians
573      }
574    }
575 
576  if(index == 5)
577    {
578      // collecte de l'energie acquise par les particules
579      for (unsigned int k = 0; k < goodPartic_.size(); k++) {
580        double gamma = goodPartic_.at(k).getGamma();
581        coord.at(k) = gamma-1.0;
582      }
583    }
584}
585
586void particleBeam::histogramInitialize(unsigned int index,vector<double>& vshf,double out[3])
587{
588  double vmin= GRAND;
589  double vmax= -vmin;
590  double vmoy= 0.0;
591  double ecatyp= 0.0;
592 
593  for (unsigned int k = 0; k < goodPartic_.size(); k++) {
594    if (vshf[k] < vmin) vmin = vshf[k];
595    else if (vshf[k] > vmax) vmax = vshf[k];
596    vmoy += vshf[k];
597    ecatyp += vshf[k]*vshf[k];
598  }
599
600  double sum= (float)goodPartic_.size();
601  out[0]= sum; 
602  vmoy /= sum;
603  out[1]= vmoy;
604  ecatyp /= sum;
605  ecatyp = sqrt(abs(ecatyp-vmoy*vmoy));
606  out[2]= ecatyp;
607
608  if(index == 0) {
609    cout << "position x -moyenne " << vmoy << " m " << "-mini " << vmin << " m " << "-maxi " << vmax << " m " << "ecart type " << ecatyp*1000.0 << " mm" << endl;
610  } 
611  if(index == 1) {
612    cout << "position y -moyenne " << vmoy << " m " << "-mini " << vmin << " m " << "-maxi " << vmax << " m " << "ecart type " << ecatyp*1000.0 << " mm" << endl;
613  } 
614  if(index == 2) {
615    cout << "position z -moyenne " << vmoy << " m " << "-mini " << vmin << " m " << "-maxi " << vmax << " m " << "ecart type " << ecatyp*1000.0 << " mm" << endl;
616  }
617  if(index == 3) {
618    cout << "divergence xp -moyenne " << vmoy << " pi.mm.mrad " << "-mini " << vmin << " pi.mm.mrad " << "-maxi " << vmax << " pi.mm.mrad " << "ecart type " << ecatyp << " pi.mm.mrad" << endl;
619  } 
620  if(index == 4) {
621    cout << "divergence yp -moyenne " << vmoy << " pi.mm.mrad " << "-mini " << vmin << " pi.mm.mrad " << "-maxi " << vmax << " pi.mm.mrad " << "ecart type " << ecatyp << " pi.mm.mrad" << endl;
622  } 
623  if(index == 5) {
624    double gmin = vmin-1.0;
625    double gmax = vmax-1.0;
626    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;
627  }
628
629  for (unsigned int k = 0; k < goodPartic_.size(); k++) {
630    vshf[k] -= vmoy;
631  }
632}
633
634void particleBeam::histogramPacking(double ecatyp,vector<double>vshf,vector<double>&xcor,vector<int>& hist)
635{
636  // demi fenetre hfene= max(3.*ecatyp-vmoy,vmoy-3.*ecatyp);
637  double hfene= 3.*ecatyp;
638  // et pas de l'histogramme
639  double hpas = hfene/25.;
640 
641  cout << "demi fenetre " << hfene << ", hpas= " << hpas << endl;
642
643  double vmin = -hfene;
644  double dfen = 2.*hfene;
645  int ihist = dfen/hpas;
646  double phist = ihist*hpas;
647  double dpas = hpas-(dfen-phist);
648  if(dpas <= hpas*1.e-03) {
649    ihist++;
650    phist= ihist*hpas;
651  }
652  double vmax= vmin+hpas*ihist;
653 
654  cout << "xAxisNumberOfBins= " << ihist <<", xAxisMinimum= " << vmin << ", xAxisMaximum= " << vmax << ", NParticules= " << vshf.size() << ", phist " << phist << endl;
655 
656  if(!xcor.empty()) xcor.clear();
657  xcor.resize(ihist+1); 
658  for (int i = 0; i <= ihist; ++i) {
659    xcor[i] = vmin+i*hpas;
660  }
661
662  /////////////////////////////////////
663
664  if(!hist.empty()) hist.clear();
665  hist.resize(ihist,0); 
666  for (unsigned int i = 0; i < vshf.size(); ++i) {
667    double var= vshf[i]-vmin;
668    if(var < 0 ) {
669      cout<<"lesser that the minimum "<<var<<", ("<< i<<")"<< endl;
670      hist[0]++;
671    }
672
673    if(var >= phist) {
674      cout<<"greater that the maximum "<<var<<", ("<< i<<")"<< endl;
675      hist[ihist-1]++;
676    }
677
678    int kk= (int)floor(var/hpas);
679    hist[kk]++;
680  }
681}
682
Note: See TracBrowser for help on using the repository browser.