Changeset 166 in PSPA


Ignore:
Timestamp:
Dec 10, 2012, 7:06:24 PM (12 years ago)
Author:
lemeur
Message:

reecriture de particleBeam

Location:
Interface_Web/trunk/pspaWT
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • Interface_Web/trunk/pspaWT/include/bareParticle.h

    r56 r166  
    3434~bareParticle() {;}
    3535
    36  inline bareParticle& operator = (const bareParticle& bp)
    37    {
    38      position_ = bp.position_;
    39      betagamma_ = bp.betagamma_;
    40      gamma_ = bp.gamma_;
    41      return *this;
    42    }
    4336
    44  bareParticle(const TRIDVECTOR&  pos , const TRIDVECTOR& bg)
    45    {
    46      position_ = pos;
    47      betagamma_ = bg;
    48      gamma_ = sqrt( 1.0 + betagamma_.norm2() );
    49    }
     37 bareParticle(const TRIDVECTOR&  pos , const TRIDVECTOR& bg);
    5038
    51  inline void resetDynamics(const bareParticle& bp)
    52   {
    53     position_ = bp.position_;
    54     betagamma_ = bp.betagamma_;
    55     gamma_ = bp.gamma_;
     39 bareParticle(bareParticle& bp);
    5640
    57   }
     41 bareParticle(const bareParticle& bp);
     42 
     43 void resetDynamics(const bareParticle& bp);
     44
     45 bareParticle& operator = (const bareParticle& bp);
    5846
    5947
     48 const TRIDVECTOR& getReferenceToPosition() const;
    6049
    61 inline  const TRIDVECTOR& getReferenceToPosition() const
    62   {
    63     return position_;
    64   }
     50 TRIDVECTOR getPosition() const;
    6551
    66 inline  TRIDVECTOR getPosition() const
    67   {
    68     return position_;
    69   }
     52 TRIDVECTOR& getReferenceToPosition();
    7053
    71  inline  TRIDVECTOR& getReferenceToPosition()
    72   {
    73     return position_;
    74   }
    75 
    76 inline  double getZ() const
    77   {
    78     return position_.getComponent(2);
    79   }
     54 double getZ() const;
    8055
    8156
    82 inline void setZ(double z)
    83   {
    84     position_.setComponent(2, z);
    85   }
     57 void setZ(double z);
    8658
    87  inline  void incrementZ( double dz)
    88   {
    89     position_.incrementComponent(2, dz);
    90   }
     59 void incrementZ( double dz);
    9160
    92 inline   void setX(double x)
    93   {
    94     position_.setComponent(0, x);
    95   }
     61 void setX(double x);
    9662
    9763
    98 inline   double getRadius() const
    99   {
    100     double auxx = position_.getComponent(0);
    101     double auxy = position_.getComponent(1);
    102 
    103     return sqrt(auxx * auxx + auxy * auxy);
    104   }
     64 double getRadius() const;
    10565
    10666
    107 inline  TRIDVECTOR getBetaGamma() const
    108   {
    109     return betagamma_;
    110   }
     67 TRIDVECTOR getBetaGamma() const;
    11168
    11269
    113 inline TRIDVECTOR& getReferenceToBetaGamma()
    114   {
    115     return betagamma_;
    116   }
     70 TRIDVECTOR& getReferenceToBetaGamma();
    11771 
    118 inline void setBetaGamma(const TRIDVECTOR& btg)
    119   {
    120     betagamma_ = btg;
    121      gamma_ = sqrt( 1.0 + betagamma_.norm2() );
    122   }
     72 void setBetaGamma(const TRIDVECTOR& btg);
    12373
    124 inline double getBetaz() const
    125   {
    126     return betagamma_.getComponent(2)/gamma_;
    127   }
     74 double getBetaz() const;
    12875
    12976
    130 inline  double getGamma() const
    131   {
    132     return gamma_;
    133   }
     77 double getGamma() const;
    13478
    13579
    … …  
    13781
    13882
    139  inline string FileOutputFlow() const
    140   {
    141     ostringstream sortie;
    142     sortie << position_.output_flow()  << betagamma_.output_flow() << " " <<  gamma_;
    143     return sortie.str();
    144   }
     83 string FileOutputFlow() const;
    14584
    146 virtual inline bool FileInput( ifstream& ifs)
    147  {
    148     bool test = false;
    149     if ( position_.input_flow(ifs) && betagamma_.input_flow(ifs) )
    150       {
    151         if (  ifs >> gamma_ )
    152           {
    153             test = true;
    154           }
    155       }
    156     return test;
    157  }
     85 virtual bool FileInput( ifstream& ifs);
    15886
    15987

  • Interface_Web/trunk/pspaWT/include/mathematicalTools.h

    r52 r166  
    9191 inline const double& operator() (int index) const { return vec_[index]; }
    9292
    93 inline void operator = (const TRIDVECTOR& triv)
    94 {
    95   setComponents( triv.vec_[0],triv.vec_[1], triv.vec_[2]);
     93inline TRIDVECTOR& operator= (const TRIDVECTOR& triv)
     94{
     95  setComponents( triv.vec_[0],triv.vec_[1], triv.vec_[2]);
     96  return *this;
    9697}
    9798

  • Interface_Web/trunk/pspaWT/include/particleBeam.h

    r153 r166  
    1111#include "bareParticle.h"
    1212#include "mathematicalTools.h"
    13 #include "mathematicalConstants.h"
    14 #include "PhysicalConstants.h"
    1513#include "nomdElements.h"
    1614
    1715using namespace std;
    1816
    19 struct particle
     17typedef struct
    2018{
    2119  float xx, xxp, begamx,yy,yyp,begamy,z, begamz,phi,wz;
    … …  
    3129    printf( " %e %e %e %e %e %e %e %e %e %e %d %d %d, %d %e %e %e %e \n", xx, xxp, begamx,yy,yyp,begamy,z, begamz,phi,wz,ne,np,ngood,npart,phi0, ksi1,ksi2,ksi3);
    3230  }
    33 };
     31} parmelaParticle;
    3432
    3533
    … …  
    4543  vector<bareParticle> goodPartic_;
    4644
    47       vector< vector<double> > rij_transportMoments_;
    48       vector<double> centroid_;
    49       double P0Transport_;
     45  vector< vector<double> > rij_transportMoments_;
     46  vector<double> centroid_;
     47  double P0Transport_;
    5048
    5149
    52       void readTransportMoments(ifstream& inp);
     50  void readTransportMoments(ifstream& inp);
    5351
    54       void impressionDesMoments() const;
     52  void impressionDesMoments() const;
     53  void razDesMoments();
    5554
    5655
    5756 public:
    5857
    59  particleBeam() 
    60    {
    61      rij_transportMoments_.resize(6);
    62      unsigned dim=0;
    63      unsigned k;
    64      for ( k=0; k < 6; k++)
    65        {
    66          rij_transportMoments_.at(k).resize(++dim);
    67        }
    68      P0Transport_ = 0.0;
    69      particleRepresentationOk_ = false;
    70      momentRepresentationOk_ = false;
    71 
    72    }
     58  particleBeam();
    7359
    7460  ~particleBeam() {;}
    7561
    7662  bool setFromParmela(unsigned numeroElement,double referencefrequency);
    77       void buildMomentRepresentation();
     63  void buildMomentRepresentation();
    7864
    79       //       bool setFromTransport(ifstream& inp, unsigned nblignes);
    80           bool  setFromTransport(string elementLabel, const nomdElements elem);
     65  bool  setFromTransport(string elementLabel, const nomdElements elem);
    8166
    82   inline void clear()
    83   {
    84     goodPartic_.clear();
    85      P0Transport_ = 0.0;
    86      particleRepresentationOk_ = false;
    87      momentRepresentationOk_ = false;
    88   }
     67  void clear();
    8968
    90 inline  int getNbParticles() const
    91   {
    92     return goodPartic_.size();
    93   }
     69  int getNbParticles() const;
    9470
    95  inline const vector< vector<double> >&  getTransportMoments() const  { return rij_transportMoments_;}
    96  inline double getP0Transport() const { return P0Transport_;}
    97   inline bool particleRepresentationOk() const {return particleRepresentationOk_;}
    98  inline bool momentRepresentationOk() const {return momentRepresentationOk_;}
     71  const vector< vector<double> >&  getTransportMoments() const;
     72  double getP0Transport() const;
     73  bool particleRepresentationOk() const;
     74  bool momentRepresentationOk() const;
    9975
    10076
    101 inline void  addParticle( bareParticle p)
    102     {
    103       goodPartic_.push_back(p);
    104     }
     77  void  addParticle( bareParticle p);
    10578
    10679
    10780
    108 inline const vector<bareParticle>& getParticleVector() const
    109  {
    110    return goodPartic_;
    111  }
     81  const vector<bareParticle>& getParticleVector() const;
    11282
    113 inline vector<bareParticle>& getParticleVector()
    114  {
    115    return goodPartic_;
    116  }
     83  vector<bareParticle>& getParticleVector();
    11784
    118  double getXmaxRms();
     85  double getXmaxRms();
    11986
    120  void getVariance(double& varx, double& vary, double& varz) const;
     87  void getVariance(double& varx, double& vary, double& varz) const;
    12188 
    12289
    123  void printAllXYZ() const;
     90  void printAllXYZ() const;
    12491
    12592
    126  void Zrange(double& zmin, double& zmax) const;
     93  void Zrange(double& zmin, double& zmax) const;
    12794
    128  void donneesDessinEllipseXxp(vector<double>& xcor, vector<double>& ycor);
     95  void donneesDessinEllipseXxp(vector<double>& xcor, vector<double>& ycor);
    12996
    13097
    131  virtual string FileOutputFlow() const;
     98  virtual string FileOutputFlow() const;
    13299
    133100
    134  virtual bool FileInput(ifstream& ifs);
     101  virtual bool FileInput(ifstream& ifs);
    135102
    136103

  • Interface_Web/trunk/pspaWT/src/GWt_pspaApplication.cc

    r161 r166  
    349349
    350350  if (!areDataCoherent()) {
    351     GWt_dialog warningDialog("PSPA : Vérification des sections", " données incohérentes !", GWt_dialog::Error,true,true);
     351    GWt_dialog warningDialog("PSPA : Vérification des sections", " donnees incoherentes !", GWt_dialog::Error,true,true);
    352352    warningDialog.exec();
    353353  }
    … …  
    476476    {
    477477      if ( beam->particleRepresentationOk() ) faireDessinParmela(toto_, beam);
    478       else addConsoleMessage("the beam state does not allow providing a drawing");
    479     }
    480   else addConsoleMessage("type of  drawing not programmed");
     478      else {
     479        addConsoleMessage("the beam state does not allow providing a drawing");
     480        GWt_dialog warningBeamState(" graphical analysis", "the beam state does not allow providing a drawing with macroparticles !", GWt_dialog::Error, false,true);
     481          warningBeamState.exec();
     482      }
     483    }
     484  else {
     485    addConsoleMessage("type of  drawing not programmed");
     486        GWt_dialog warningTypeDrawing(" graphical analysis", "type of  drawing not programmed !", GWt_dialog::Error, false,true);
     487           warningTypeDrawing.exec();
     488 }
    481489  //////////////////////////////////////////
    482490}
    … …  
    492500      faireDessinEnveloppe(toto_, "x");
    493501    }
    494   else addConsoleMessage("type of enveloppe drawing not programmed");
     502  else {
     503    addConsoleMessage("type of enveloppe drawing not programmed");
     504        GWt_dialog warningTypeEnveloppe(" graphical analysis", "type of enveloppe drawing not programmed !", GWt_dialog::Error, false,true);
     505           warningTypeEnveloppe.exec();
     506  }
    495507  //////////////////////////////////////////
    496508}
    … …  
    591603    model->setData(i, 0, xcor.at(i));
    592604    model->setData(i, 1, ycor.at(i));
    593     cout << " PspaApplication::scatterPlot1D el= " << i+1 << " x= " << xcor.at(i) << " y= " << ycor.at(i) << endl;
     605    //    cout << " PspaApplication::scatterPlot1D el= " << i+1 << " x= " << xcor.at(i) << " y= " << ycor.at(i) << endl;
    594606  }
    595607

  • Interface_Web/trunk/pspaWT/src/bareParticle.cc

    r54 r166  
    11#include "bareParticle.h"
     2
     3  bareParticle::bareParticle(const TRIDVECTOR&  pos , const TRIDVECTOR& bg)
     4   {
     5     position_ = pos;
     6     betagamma_ = bg;
     7     gamma_ = sqrt( 1.0 + betagamma_.norm2() );
     8     cout << " bareParticle constructeur " << endl;
     9   }
     10
     11 bareParticle::bareParticle(bareParticle& bp)
     12   {
     13    position_ = bp.position_;
     14    betagamma_ = bp.betagamma_;
     15    gamma_ = bp.gamma_;
     16     cout << " bareParticle copie simple " << endl;
     17   }
     18
     19 bareParticle::bareParticle(const bareParticle& bp)
     20   {
     21    position_ = bp.position_;
     22    betagamma_ = bp.betagamma_;
     23    gamma_ = bp.gamma_;
     24     cout << " bareParticle copie const " << endl;
     25   }
     26
     27 void bareParticle::resetDynamics(const bareParticle& bp)
     28  {
     29    position_ = bp.position_;
     30    betagamma_ = bp.betagamma_;
     31    gamma_ = bp.gamma_;
     32  }
     33
     34 bareParticle& bareParticle::operator = (const bareParticle& bp)
     35   {
     36     position_ = bp.position_;
     37     betagamma_ = bp.betagamma_;
     38     gamma_ = bp.gamma_;
     39     cout << " bareParticle operateur = " << endl;
     40     return *this;
     41   }
     42
     43const TRIDVECTOR&  bareParticle::getReferenceToPosition() const
     44  {
     45    return position_;
     46  }
     47
     48TRIDVECTOR  bareParticle::getPosition() const
     49  {
     50    return position_;
     51  }
     52
     53 TRIDVECTOR&  bareParticle::getReferenceToPosition()
     54  {
     55    return position_;
     56  }
     57
     58double  bareParticle::getZ() const
     59  {
     60    return position_.getComponent(2);
     61  }
     62
     63
     64void  bareParticle::setZ(double z)
     65  {
     66    position_.setComponent(2, z);
     67  }
     68
     69void  bareParticle::incrementZ( double dz)
     70  {
     71    position_.incrementComponent(2, dz);
     72  }
     73
     74void  bareParticle::setX(double x)
     75  {
     76    position_.setComponent(0, x);
     77  }
     78
     79
     80double  bareParticle::getRadius() const
     81  {
     82    double auxx = position_.getComponent(0);
     83    double auxy = position_.getComponent(1);
     84
     85    return sqrt(auxx * auxx + auxy * auxy);
     86  }
     87
     88
     89TRIDVECTOR  bareParticle::getBetaGamma() const
     90  {
     91    return betagamma_;
     92  }
     93
     94
     95TRIDVECTOR&  bareParticle::getReferenceToBetaGamma()
     96  {
     97    return betagamma_;
     98  }
     99 
     100void  bareParticle::setBetaGamma(const TRIDVECTOR& btg)
     101  {
     102    betagamma_ = btg;
     103     gamma_ = sqrt( 1.0 + betagamma_.norm2() );
     104  }
     105
     106double  bareParticle::getBetaz() const
     107  {
     108    return betagamma_.getComponent(2)/gamma_;
     109  }
     110
     111
     112double  bareParticle::getGamma() const
     113  {
     114    return gamma_;
     115  }
     116
     117string bareParticle::FileOutputFlow() const
     118  {
     119    ostringstream sortie;
     120    sortie << position_.output_flow()  << betagamma_.output_flow() << " " <<  gamma_;
     121    return sortie.str();
     122  }
     123
     124bool bareParticle::FileInput( ifstream& ifs)
     125 {
     126    bool test = false;
     127    if ( position_.input_flow(ifs) && betagamma_.input_flow(ifs) )
     128      {
     129        if (  ifs >> gamma_ )
     130          {
     131            test = true;
     132          }
     133      }
     134    return test;
     135 }
     136
    2137
    3138
    … …  
    13148    betas.print();
    14149    cout << " gamma = " << gamma_  << endl;
    15 
    16150  }

  • Interface_Web/trunk/pspaWT/src/dataManager.cc

    r155 r166  
    362362  removeFile("transport.input");
    363363  removeFile("transport.output");
    364   //  currentBeam_.clear();
     364  diagnosticBeam_.clear();
    365365  currentBeam_ = NULL;
    366366  clearSectionToExecute();

  • Interface_Web/trunk/pspaWT/src/nomdElements.cc

    r160 r166  
    4141{
    4242  switch(eType) {
    43   case RFgun : return "RFgun";
     43  case RFgun : return "rfgun";
    4444  case drift :  return "drift";
    4545  case cell :  return "cell";

  • Interface_Web/trunk/pspaWT/src/particleBeam.cc

    r153 r166  
    22#include "mathematicalConstants.h"
    33#include "PhysicalConstants.h"
    4 //#include "nomdElements.h"
    54//#include <string>
    65#include <stdio.h>
    … …  
    98#include "environmentVariables.h"
    109
     10
     11
     12particleBeam::particleBeam()  {
     13  rij_transportMoments_.resize(6);
     14  unsigned dim=0;
     15  unsigned k;
     16  for ( k=0; k < 6; k++){
     17    rij_transportMoments_.at(k).resize(++dim);
     18  }
     19  P0Transport_ = 0.0;
     20  particleRepresentationOk_ = false;
     21  momentRepresentationOk_ = false;
     22}
     23
     24void particleBeam::clear() {
     25  unsigned k,j;
     26  goodPartic_.clear();
     27  for ( k=0; k < 6; k++){
     28    for ( j=0; j <= k; j++) {
     29      ( rij_transportMoments_.at(k) ).at(j) = 0.0;
     30    }
     31  }
     32  P0Transport_ = 0.0;
     33  particleRepresentationOk_ = false;
     34  momentRepresentationOk_ = false;
     35}
     36
     37int particleBeam::getNbParticles() const {
     38  return goodPartic_.size();
     39}
     40
     41const vector< vector<double> >&   particleBeam::getTransportMoments() const  {
     42  return rij_transportMoments_;
     43}
     44double particleBeam::getP0Transport() const {
     45  return P0Transport_;
     46}
     47bool particleBeam::particleRepresentationOk() const {
     48  return particleRepresentationOk_;
     49}
     50bool particleBeam::momentRepresentationOk() const {
     51  return momentRepresentationOk_;
     52}
     53
     54void  particleBeam::addParticle( bareParticle p)
     55{
     56  goodPartic_.push_back(p);
     57}
     58
     59const vector<bareParticle>& particleBeam::getParticleVector() const
     60{
     61  return goodPartic_;
     62}
     63vector<bareParticle>& particleBeam::getParticleVector()
     64{
     65  return goodPartic_;
     66}
     67
     68
    1169bool  particleBeam::setFromTransport(string elLab, const nomdElements elem)
    1270{
    13    string elementLabel = elLab;
     71  string elementLabel = elLab;
    1472  // transformer le label en majuscules ; je ne suis pas sur que ca
    1573  // marchera a tous les coups (glm)
    16 
    1774  std::transform(elementLabel.begin(), elementLabel.end(), elementLabel.begin(), (int (*)(int))std::toupper);
     75
    1876  cout << " particleBeam::setFromTransport on cherche element " << elementLabel << endl;
    1977  string buf;
    … …  
    2179  string nameIn = WORKINGAREA + "transport.output";
    2280  infile.open(nameIn.c_str(), ios::in);
    23   if (!infile)
    24     {
    25       cerr << " particleBeam::setFromTransport : error re-opening transport output stream " << nameIn << endl;
    26       return false;
    27     }
     81  if (!infile) {
     82    cerr << " particleBeam::setFromTransport : error re-opening transport output stream " << nameIn << endl;
     83    return false;
     84  }
    2885  //  else cout << " particleBeam::setFromTransport() : ouverture du fichier " << nameIn << endl;
    2986
    3087  string::size_type nn = string::npos;
    31   while ( getline(infile, buf) )
    32     {
    33       //      cout << " buf= "  << buf << endl;
    34       nn = buf.find(elementLabel);
    35       //     cout << " string::npos= " << string::npos << " nn= " << nn << endl;
    36       if ( nn != string::npos )
    37         {
    38           //      cout << " particleBeam::setFromTransport : element " << elementLabel << " trouve " << endl;
    39           break;
    40         }
    41     }
    42   if ( nn == string::npos )
    43     {
    44       cerr << " particleBeam::setFromTransport : element " << elementLabel << " non trouve dans le fichier  " << nameIn << endl;
    45       return false;
    46     }
    47   if (elem.getElementType() == bend )
    48     {
    49       getline(infile, buf);
    50       getline(infile, buf);
    51     }
     88  while ( getline(infile, buf) ) {
     89    //      cout << " buf= "  << buf << endl;
     90    nn = buf.find(elementLabel);
     91    //     cout << " string::npos= " << string::npos << " nn= " << nn << endl;
     92    if ( nn != string::npos ) {
     93      //          cout << " particleBeam::setFromTransport : element " << elementLabel << " trouve " << endl;
     94      break;
     95    }
     96  }
     97
     98  if ( nn == string::npos ) {
     99    cerr << " particleBeam::setFromTransport : element " << elementLabel << " non trouve dans le fichier  " << nameIn << endl;
     100    return false;
     101  }
     102  if (elem.getElementType() == bend ) {
     103    getline(infile, buf);
     104    getline(infile, buf);
     105  }
    52106  readTransportMoments(infile);
    53107  //  impressionDesMoments();
    54    infile.close();
    55    momentRepresentationOk_ = true;
     108  infile.close();
     109  momentRepresentationOk_ = true;
    56110  return true;
    57111}
    58112
    59113
    60 bool particleBeam::setFromParmela(unsigned numeroElement, double referencefrequency)
    61 {
    62   unsigned int k;
     114bool particleBeam::setFromParmela(unsigned numeroElement, double referencefrequency) {
     115  unsigned  k;
    63116  FILE* filefais;
    64117  string nomfilefais = WORKINGAREA + "parmdesz";
    … …  
    66119  filefais = fopen(nomfilefais.c_str(), "r");
    67120
    68   if ( filefais == (FILE*)0 )
    69     {
    70       cerr << " particleBeam::setFromParmela() erreur a l'ouverture du fichier" << nomfilefais  << endl;;
    71       return false;
    72     }
     121  if ( filefais == (FILE*)0 ) {
     122    cerr << " particleBeam::setFromParmela() erreur a l'ouverture du fichier" << nomfilefais  << endl;;
     123    return false;
     124  }
    73125  else cout << " particleBeam::setFromParmela() : ouverture du fichier " << nomfilefais << endl;
    74126
    75   struct particle partic;
    76   struct particle* partRefPtr=NULL;
    77 
    78   std::vector<struct particle> faisceau;
    79 
    80   //  int numElem=-1;
    81   partRefPtr=NULL;
     127  parmelaParticle partic;
     128  std::vector<parmelaParticle> faisceau;
     129
    82130  cout << " particleBeam::setFromParmela : numeroElement = " << numeroElement << endl;
    83131  numeroElement--;
    84   while( partic.readFromParmelaFile(filefais) > 0 )
    85 
    86     {
    87       // // on ne va conserver que les resultats a la sortie du dernier element
    88       // if ( partic.ne != numElem )
    89       //        {
    90       //          faisceau.clear();
    91       //          partRefPtr=NULL;
    92       //          numElem = partic.ne;
    93       //        }
    94       //      cout << " partic.ne = " << partic.ne << endl;
    95       if ( partic.ne == numeroElement )
    96         {
    97           //      cout << " trouve particule " << endl;
    98           faisceau.push_back(partic);
    99           if ( partic.np == 1 )
    100             {
    101               if ( fabs(partic.xx) > EPSILON || fabs(partic.yy) > EPSILON || fabs(partic.xxp) > EPSILON  || fabs(partic.yyp) > EPSILON)
    102                 {
    103                   printf(" ATTENTION part. reference douteuse  \n");
    104                   partic.imprim();
    105                 }
    106               partRefPtr=&faisceau.back();
    107           //      cout << " indice de la part. ref. " << faisceau.size()-1 << endl;
    108           //          printf("part. reference \n");
    109           //          partRefPtr->imprim();
    110             }
     132
     133
     134
     135  int testNombrePartRef =0;
     136  double phaseRef;
     137
     138  while( partic.readFromParmelaFile(filefais) > 0 ) {
     139    if ( partic.ne == numeroElement )
     140      {
     141        faisceau.push_back(partic);
     142
     143        if ( partic.np == 1 ) {
     144          // en principe on est sur la particule de reference
     145          if ( fabs(partic.xx) > EPSILON || fabs(partic.yy) > EPSILON || fabs(partic.xxp) > EPSILON  || fabs(partic.yyp) > EPSILON) {
     146            printf(" ATTENTION part. reference douteuse  \n");
     147            partic.imprim();
     148          }
     149          phaseRef = partic.phi;
     150          TRIDVECTOR  posRef(partic.xx,partic.yy,0.0);
     151          TRIDVECTOR betagammaRef(partic.xxp*partic.begamz, partic.yyp*partic.begamz, partic.begamz);
     152          referenceParticle_ = bareParticle(posRef, betagammaRef);
     153          testNombrePartRef++;
     154          if ( testNombrePartRef != 1 ) {
     155            cerr << " TROP DE PART. DE REF : " << testNombrePartRef << " !! " << endl;
     156            return false;
     157          }
    111158        }
    112     }
    113 
    114   if ( faisceau.size() == 0) return false;
    115 
    116   //  printf("dernier element %d \n", numElem);
     159      }
     160  }
     161
     162  if ( faisceau.size() == 0)
     163    {
     164      cerr << " particleBeam::setFromParmela echec lecture  element " << numeroElement+1 << endl;
     165      return false;
     166    }
     167 
    117168  // facteur  c/ 360. pour calculer (c dphi) / (360.freq)
    118169  // avec freq en Mhz et dphi en degres et résultat en cm:
    119170  double FACTEUR =  83.3333;  // ameliorer la precision
     171
     172
     173
     174  // pour l'instant on choisit un centroid nul;
     175  centroid_ = vector<double>(6,0.0);
     176
     177  goodPartic_.clear();
     178  goodPartic_.resize(faisceau.size(), bareParticle());
    120179  double x,xp,y,yp;
    121180  double betagammaz;
     181  double betaz;
     182  double deltaz;
     183  double dephas;
     184  double g;
     185  TRIDVECTOR  pos;
     186  TRIDVECTOR betagamma;
    122187  // contrairement a ce qu'indique la notice PARMELA, dans parmdesz, les xp et yp
    123188  // sont donnes en radians
    124   // particule de reference
    125  
    126   x=partRefPtr->xx;
    127   xp=partRefPtr->xxp;
    128   y=partRefPtr->yy;
    129   yp=partRefPtr->yyp;
    130   betagammaz = partRefPtr->begamz;
    131   TRIDVECTOR  posRef(x,y,0.0);
    132   TRIDVECTOR betagammaRef(xp*betagammaz, yp*betagammaz, betagammaz);
    133 
    134   referenceParticle_ = bareParticle(posRef, betagammaRef);
    135 
    136 
    137   // pour l'instant on choisit un centroid nul;
    138   centroid_ = vector<double>(6,0.0);
    139 
    140   goodPartic_.clear();
    141   for ( k=0; k < faisceau.size(); k++)
    142     //  for (int k=0; k < 10; k++)
    143     {
    144       x=faisceau.at(k).xx;
    145       xp=faisceau.at(k).xxp;
    146       y=faisceau.at(k).yy;
    147       yp=faisceau.at(k).yyp;
    148       // dephasage par rapport a la reference 
    149       double dephas = faisceau.at(k).phi - partRefPtr->phi; // degrés
    150       double g = faisceau.at(k).wz/ERESTMeV;
    151       betagammaz = faisceau.at(k).begamz;
    152       double betaz = betagammaz/(g+1.0);
    153       double deltaz = FACTEUR * betaz * dephas / referencefrequency;
    154       x += xp * deltaz;
    155       y += yp * deltaz;
    156       TRIDVECTOR  pos(x,y,deltaz);
    157       TRIDVECTOR betagamma(xp*betagammaz, yp*betagammaz, betagammaz);
    158       bareParticle bp(pos,betagamma);
    159       goodPartic_.push_back(bp);
    160     }
     189  for ( k=0; k < faisceau.size(); k++) {
     190    x=faisceau.at(k).xx;
     191    xp=faisceau.at(k).xxp;
     192    y=faisceau.at(k).yy;
     193    yp=faisceau.at(k).yyp;
     194
     195    // dephasage par rapport a la reference 
     196    dephas = faisceau.at(k).phi - phaseRef; // degrés
     197    g = faisceau.at(k).wz/ERESTMeV;
     198    betagammaz = faisceau.at(k).begamz;
     199    betaz = betagammaz/(g+1.0);
     200    deltaz = FACTEUR * betaz * dephas / referencefrequency;
     201    x += xp * deltaz;
     202    y += yp * deltaz;
     203    pos.setComponents(x,y,deltaz);
     204    betagamma.setComponents(xp*betagammaz, yp*betagammaz, betagammaz);
     205    goodPartic_.at(k) = bareParticle(pos,betagamma);
     206  }
    161207  particleRepresentationOk_ = true;
    162   //  buildMomentRepresentation();
    163208  return true;
    164209}
    165210
    166211
    167 void particleBeam::getVariance(double& varx, double& vary, double& varz) const
    168 {
    169 
     212void particleBeam::getVariance(double& varx, double& vary, double& varz) const {
     213  unsigned int k;
    170214  double x,y,z;
    171215  double xav = 0.;
    … …  
    176220  double zavsq = 0.;
    177221
    178 
    179222  TRIDVECTOR pos;
    180223
    181   unsigned int k;
    182 
    183   for ( k = 0 ; k < goodPartic_.size(); k++)
    184     {
    185       pos = goodPartic_.at(k).getPosition();
    186       pos.getComponents(x,y,z);
    187       //      partic_[k].getXYZ(x,y,z);
    188       xav += x;
    189       xavsq += x*x;
    190       yav += y;
    191       yavsq += y*y;
    192       zav += z;
    193       zavsq += z*z;
    194     }
     224
     225  for ( k = 0 ; k < goodPartic_.size(); k++) {
     226    pos = goodPartic_.at(k).getPosition();
     227    pos.getComponents(x,y,z);
     228    //      partic_[k].getXYZ(x,y,z);
     229    xav += x;
     230    xavsq += x*x;
     231    yav += y;
     232    yavsq += y*y;
     233    zav += z;
     234    zavsq += z*z;
     235  }
    195236
    196237  double aginv = double (goodPartic_.size());
    … …  
    200241  vary =  aginv * ( yavsq - yav*yav*aginv );
    201242  varz =  aginv * ( zavsq - zav*zav*aginv );
    202 
    203 }
    204 
    205 
    206 void particleBeam::printAllXYZ() const
    207 {
     243}
     244
     245
     246void particleBeam::printAllXYZ() const {
    208247  cout << " dump du faisceau : " << endl;
    209 
    210248  cout <<  goodPartic_.size() << " particules " << endl;
    211249  unsigned int k;
    … …  
    220258
    221259
    222 void particleBeam::Zrange(double& zmin, double& zmax) const
    223 {
     260void particleBeam::Zrange(double& zmin, double& zmax) const {
    224261  double z;
    225262  zmin = GRAND;
    … …  
    237274
    238275
    239 string particleBeam::FileOutputFlow() const
    240 {
     276string particleBeam::FileOutputFlow() const {
    241277  ostringstream sortie;
    242278  unsigned int k;
    … …  
    249285}
    250286
    251 bool particleBeam::FileInput( ifstream& ifs)
    252 {
     287bool particleBeam::FileInput( ifstream& ifs) {
    253288  bool test = true;
    254289  string dum1, dum2;
    … …  
    264299}
    265300
    266 void particleBeam::buildMomentRepresentation()
    267 {
     301void particleBeam::buildMomentRepresentation() {
    268302
    269303  unsigned k,j,m;
    270304  double auxj, auxm;
    271 
    272305  if ( !particleRepresentationOk_)
    273306    {
    … …  
    283316
    284317  // initialisation des moments
     318  razDesMoments();
     319
     320  // accumulation
     321  TRIDVECTOR pos;
     322  TRIDVECTOR begam;
     323  double gamma;
     324  double begamz;
     325  double g;
     326  double PMeVsc;
     327  double del;
     328  vector<double> part(6);
     329  for (k=0; k < goodPartic_.size(); k++) {
     330    gamma = goodPartic_.at(k).getGamma();
     331    pos = goodPartic_.at(k).getPosition();
     332    begam= goodPartic_.at(k).getBetaGamma();
     333    begamz = begam.getComponent(2);
     334    g = gamma -1.0;
     335    PMeVsc = sqrt( g*(g+2) );
     336    del = 100.0 * ( PMeVsc -  P_reference_MeV_sur_c ) / P_reference_MeV_sur_c ; // en %
     337
     338    part[0] = pos.getComponent(0);
     339    part[1] = begam.getComponent(0)/begamz;
     340    part[2] = pos.getComponent(1);
     341    part[3] = begam.getComponent(1)/begamz;
     342    part[4] = pos.getComponent(2);
     343    part[5] = del;
     344
     345    for ( j = 0; j < 6; j++) {
     346      auxj = part.at(j) - centroid_.at(j);
     347      for (m=0; m <= j; m++)
     348        {
     349          auxm = part.at(m) - centroid_.at(m);
     350          ( rij_transportMoments_.at(j) ).at(m) += auxj*auxm;
     351          //          cout << " j= " << j << " m= " << m << " rjm= " << ( rij_transportMoments_.at(j) ).at(m) << endl;
     352        }
     353    }
     354  }
     355
     356
     357  // moyenne
     358  double facmoy = 1.0/double( goodPartic_.size() );
     359  for ( j = 0; j < 6; j++) {
     360    ( rij_transportMoments_.at(j) ).at(j) = sqrt(( rij_transportMoments_.at(j) ).at(j) * facmoy );
     361  }
     362
     363  for ( j = 0; j < 6; j++) {
     364    auxj =  ( rij_transportMoments_.at(j) ).at(j);
     365    for (m=0; m < j; m++) {
     366      auxm = ( rij_transportMoments_.at(m) ).at(m);
     367      ( rij_transportMoments_.at(j) ).at(m) *= facmoy/(auxj * auxm);
     368    }
     369  }
     370   
     371  // les longueurs sont en cm
     372  // les angles en radians, on passe en mrad;
     373
     374  double uniteAngle = 1.0e+3;
     375  ( rij_transportMoments_.at(1) ).at(1)  *= uniteAngle;
     376
     377  ( rij_transportMoments_.at(3) ).at(3)  *= uniteAngle;
     378   
     379  P0Transport_ = 1.0e-3*ERESTMeV*P_reference_MeV_sur_c;
     380
     381  //  cout << " buildmomentrepresentation impression des moments " << endl;
     382  //  impressionDesMoments();
     383
     384  momentRepresentationOk_ = true;
     385}
     386
     387void particleBeam::impressionDesMoments() const {
     388  unsigned j,m;
     389  cout << " impression des moments " << endl;
    285390  for ( j = 0; j < 6; j++)
     391    {
     392      for (m=0; m <= j; m++)
     393        {
     394          cout  << ( rij_transportMoments_.at(j) ).at(m) << " ";
     395        }
     396      cout << endl;
     397    }
     398}
     399
     400void particleBeam::razDesMoments() {
     401  // initialisation des moments
     402  unsigned j,m;
     403  for ( j = 0; j < rij_transportMoments_.size(); j++)
    286404    {
    287405      for (m=0; m <= j; m++)
    … …  
    290408        }
    291409    }
    292 
    293   // accumulation
    294   for (k=0; k < goodPartic_.size(); k++)
    295     {
    296       bareParticle bp = goodPartic_.at(k);
    297       double gamma = bp.getGamma();
    298       TRIDVECTOR pos = bp.getPosition();
    299       TRIDVECTOR begam= bp.getBetaGamma();
    300       double begamz = begam.getComponent(2);
    301       double g = gamma -1.0;
    302       double PMeVsc = sqrt( g*(g+2) );
    303       double del = 100.0 * ( PMeVsc -  P_reference_MeV_sur_c ) / P_reference_MeV_sur_c ; // en %
    304 
    305       vector<double> part(6);
    306       part[0] = pos.getComponent(0);
    307       part[1] = begam.getComponent(0)/begamz;
    308       part[2] = pos.getComponent(1);
    309       part[3] = begam.getComponent(1)/begamz;
    310 
    311       part[4] = pos.getComponent(2);
    312 
    313       part[5] = del;
    314 
    315       //      cout << " buildMomentRepresentation part. " << k << " x= " << part[0] << " xp= " << part[1] << " y= " << part[2] << " yp= " << part[3] << endl;
    316 
    317       for ( j = 0; j < 6; j++)
    318         {
    319           auxj = part.at(j) - centroid_.at(j);
    320           for (m=0; m <= j; m++)
    321             {
    322               auxm = part.at(m) - centroid_.at(m);
    323               ( rij_transportMoments_.at(j) ).at(m) += auxj*auxm;
    324               //              cout << " j= " << j << " m= " << m << " rjm= " << ( rij_transportMoments_.at(j) ).at(m) << endl;
    325             }
    326         }
    327     }
    328 
    329 
    330   // moyenne
    331   double facmoy = 1.0/double( goodPartic_.size() );
    332   for ( j = 0; j < 6; j++)
    333     {
    334       ( rij_transportMoments_.at(j) ).at(j) = sqrt(( rij_transportMoments_.at(j) ).at(j) * facmoy );
    335     }
    336 
    337   for ( j = 0; j < 6; j++)
    338     {
    339       auxj =  ( rij_transportMoments_.at(j) ).at(j);
    340       for (m=0; m < j; m++)
    341         {
    342           auxm = ( rij_transportMoments_.at(m) ).at(m);
    343           ( rij_transportMoments_.at(j) ).at(m) *= facmoy/(auxj * auxm);
    344         }
    345     }
    346    
    347   // les longueurs sont en cm
    348   // les angles en radians, on passe en mrad;
    349 
    350   double uniteAngle = 1.0e+3;
    351   ( rij_transportMoments_.at(1) ).at(1)  *= uniteAngle;
    352 
    353   ( rij_transportMoments_.at(3) ).at(3)  *= uniteAngle;
    354    
    355   P0Transport_ = 1.0e-3*ERESTMeV*P_reference_MeV_sur_c;
    356 
    357   // cout << " impression des moments " << endl;
    358   // for ( j = 0; j < 6; j++)
    359   //   {
    360   //     for (m=0; m <= j; m++)
    361   //    {
    362   //          cout  << ( rij_transportMoments_.at(j) ).at(m) << " ";
    363   //    }
    364   //     cout << endl;
    365   //   }
    366 
    367   momentRepresentationOk_ = true;
    368 }
    369 
    370 void particleBeam::impressionDesMoments() const
    371 {
    372   unsigned j,m;
    373   cout << " impression des moments " << endl;
    374   for ( j = 0; j < 6; j++)
    375     {
    376       for (m=0; m <= j; m++)
    377         {
    378               cout  << ( rij_transportMoments_.at(j) ).at(m) << " ";
    379         }
    380       cout << endl;
    381     }
    382 }
    383 
    384 
    385 
    386 void particleBeam::readTransportMoments(ifstream& inp)
    387 {
    388         unsigned j,m;
    389       string  bidString;
    390        double bidon;
     410}
     411
     412
     413void particleBeam::readTransportMoments(ifstream& inp) {
     414  string  bidString;
     415  double bidon;
    391416
    392417  // initialisation des moments
    393   for ( j = 0; j < 6; j++)
    394     {
    395       for (m=0; m <= j; m++)
    396         {
    397           ( rij_transportMoments_.at(j) ).at(m) = 0.0; // element r_jm
    398         }
    399     }
    400 
    401         inp >> bidon >>  bidString >>  bidon >>  ( rij_transportMoments_.at(0) ).at(0) >> bidString;
    402         inp >> bidon >> ( rij_transportMoments_.at(1) ).at(1) >>  bidString >> ( rij_transportMoments_.at(1) ).at(0);
    403         inp >> bidon >> ( rij_transportMoments_.at(2) ).at(2) >>  bidString >> ( rij_transportMoments_.at(2) ).at(0)  >> ( rij_transportMoments_.at(2) ).at(1);
    404         inp >> bidon >> ( rij_transportMoments_.at(3) ).at(3) >>  bidString >> ( rij_transportMoments_.at(3) ).at(0)  >> ( rij_transportMoments_.at(3) ).at(1) >> ( rij_transportMoments_.at(3) ).at(2);
    405 
    406         inp >> bidon >> ( rij_transportMoments_.at(4) ).at(4) >>  bidString >> ( rij_transportMoments_.at(4) ).at(0)  >> ( rij_transportMoments_.at(4) ).at(1) >> ( rij_transportMoments_.at(4) ).at(2) >> ( rij_transportMoments_.at(4) ).at(3);
    407 
    408         inp >> bidon >> ( rij_transportMoments_.at(5) ).at(5) >>  bidString >> ( rij_transportMoments_.at(5) ).at(0)  >> ( rij_transportMoments_.at(5) ).at(1) >> ( rij_transportMoments_.at(5) ).at(2) >> ( rij_transportMoments_.at(5) ).at(3) >> ( rij_transportMoments_.at(5) ).at(4);
    409 
    410 }
    411 
    412 double particleBeam::getXmaxRms()
    413 {
     418  razDesMoments();
     419
     420  inp >> bidon >>  bidString >>  bidon >>  ( rij_transportMoments_.at(0) ).at(0) >> bidString;
     421  inp >> bidon >> ( rij_transportMoments_.at(1) ).at(1) >>  bidString >> ( rij_transportMoments_.at(1) ).at(0);
     422  inp >> bidon >> ( rij_transportMoments_.at(2) ).at(2) >>  bidString >> ( rij_transportMoments_.at(2) ).at(0)  >> ( rij_transportMoments_.at(2) ).at(1);
     423  inp >> bidon >> ( rij_transportMoments_.at(3) ).at(3) >>  bidString >> ( rij_transportMoments_.at(3) ).at(0)  >> ( rij_transportMoments_.at(3) ).at(1) >> ( rij_transportMoments_.at(3) ).at(2);
     424
     425  inp >> bidon >> ( rij_transportMoments_.at(4) ).at(4) >>  bidString >> ( rij_transportMoments_.at(4) ).at(0)  >> ( rij_transportMoments_.at(4) ).at(1) >> ( rij_transportMoments_.at(4) ).at(2) >> ( rij_transportMoments_.at(4) ).at(3);
     426
     427  inp >> bidon >> ( rij_transportMoments_.at(5) ).at(5) >>  bidString >> ( rij_transportMoments_.at(5) ).at(0)  >> ( rij_transportMoments_.at(5) ).at(1) >> ( rij_transportMoments_.at(5) ).at(2) >> ( rij_transportMoments_.at(5) ).at(3) >> ( rij_transportMoments_.at(5) ).at(4);
     428
     429}
     430
     431double particleBeam::getXmaxRms() {
    414432  if ( !momentRepresentationOk_ ) buildMomentRepresentation();
    415433  return ( rij_transportMoments_.at(0) ).at(0);
    416434}
    417435
    418 void particleBeam::donneesDessinEllipseXxp(vector<double>& xcor, vector<double>& ycor)
    419 {
    420       int k;
    421       double x,y;
     436void particleBeam::donneesDessinEllipseXxp(vector<double>& xcor, vector<double>& ycor) {
     437  int k;
     438  double x,y;
    422439
    423440  if ( !momentRepresentationOk_ ) return;

  • Interface_Web/trunk/pspaWT/workingArea/pspa.save

    r149 r166  
    222998.65 1 100000 10
    330
    4 beam01
     4rfgun01
    55500  3.6  0.02
    661e-06 1e-06
Note: See TracChangeset for help on using the changeset viewer.