Ignore:
Timestamp:
Jun 10, 2013, 4:52:52 PM (11 years ago)
Author:
lemeur
Message:

mise a jour des manipulations de faisceau

File:
1 edited

Legend:

Unmodified
Added
Removed
  • Interface_Web/trunk/pspaWT/sources/controler/src/particleBeam.cc

    r414 r417  
    1919
    2020void particleBeam::clear() {
    21   goodPartic_.clear();
     21  relativePartic_.clear();
    2222  rij_.raz();
    2323  P0Transport_ = 0.0;
     
    2727
    2828int particleBeam::getNbParticles() const {
    29   return goodPartic_.size();
     29  return relativePartic_.size();
    3030}
    3131
     
    9999  centroid_ = centroid;
    100100  referenceParticle_ = referencePart;
    101   goodPartic_.clear();
    102   goodPartic_ = particles;
    103   cout << " particleBeam::setWithParticles taille vect. part. ENREGISTRE " << goodPartic_.size() << endl;
     101  relativePartic_.clear();
     102  relativePartic_ = particles;
     103  cout << " particleBeam::setWithParticles taille vect. part. ENREGISTRE " << relativePartic_.size() << endl;
    104104  particleRepresentationOk_ = true;
    105105}
     
    114114void  particleBeam::addParticle( bareParticle p)
    115115{
    116   goodPartic_.push_back(p);
     116  relativePartic_.push_back(p);
    117117}
    118118
    119119const vector<bareParticle>& particleBeam::getParticleVector() const
    120120{
    121   return goodPartic_;
     121  return relativePartic_;
    122122}
    123123
    124124vector<bareParticle>& particleBeam::getParticleVector()
    125125{
    126   return goodPartic_;
    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 }
     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// }
    161161
    162162
    163163void particleBeam::printAllXYZ() const {
    164164  cout << " dump du faisceau : " << endl;
    165   cout <<  goodPartic_.size() << " particules " << endl;
     165  cout <<  relativePartic_.size() << " particules " << endl;
    166166  unsigned int k;
    167   for ( k = 0 ; k < goodPartic_.size(); k++)
     167  for ( k = 0 ; k < relativePartic_.size(); k++)
    168168    {
    169169      double xx,yy,zz;
    170       goodPartic_.at(k).getPosition().getComponents(xx,yy,zz);
     170      relativePartic_.at(k).getPosition().getComponents(xx,yy,zz);
    171171      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 
    179 void particleBeam::Zrange(double& zmin, double& zmax) const {
     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 {
    180180  double z;
    181   zmin = GRAND;
    182   zmax = -zmin;
     181  cdtmin = GRAND;
     182  cdtmax = -cdtmin;
    183183
    184184  unsigned int k;
    185   for ( k = 0 ; k < goodPartic_.size(); k++)
     185  for ( k = 0 ; k < relativePartic_.size(); k++)
    186186    {
    187       z = goodPartic_.at(k).getZ();
    188       if ( z < zmin ) zmin = z;
    189       else if ( z > zmax) zmax = z;         
    190     }
    191 }
    192 
    193 
    194 
    195 string particleBeam::FileOutputFlow() const {
     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 {
    196196  ostringstream sortie;
    197197  unsigned int k;
    198   for ( k = 0 ; k < goodPartic_.size(); k++)
     198  for ( k = 0 ; k < relativePartic_.size(); k++)
    199199    {
    200       sortie << goodPartic_.at(k).FileOutputFlow() << endl;
     200      sortie << relativePartic_.at(k).FileOutputFlow() << endl;
    201201    }
    202202  sortie << endl;
     
    219219
    220220void particleBeam::buildMomentRepresentation() {
    221 
     221  // le faisceau cense etre represente en un z donne (reference) sera
     222  // ici deploye spatialement (faisceau en un temps donne)
    222223  unsigned k,j,m;
    223224  double auxj, auxm;
     
    248249  double PMeVsc;
    249250  double del;
     251  //  double xp,yp,dz,cdt;
    250252  vector<double> part(6, 0.0);
    251253
    252254    vector< vector<double> >& matrice = rij_.getMatrix();
    253255
    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();
     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();
    259266    begamz = begam.getComponent(2);
    260267    g = gamma -1.0;
     
    262269    del = 100.0 * ( PMeVsc -  P_reference_MeV_sur_c ) / P_reference_MeV_sur_c ; // en %
    263270
    264     part[0] = pos.getComponent(0);
     271    //    dz = begamz * cdt / gamma;
     272    //    xp = begam.getComponent(0)/begamz;
     273    //    yp = begam.getComponent(1)/begamz;
     274
     275    part[0] = positionDeployee.getComponent(0);
    265276    part[1] = begam.getComponent(0)/begamz;
    266     part[2] = pos.getComponent(1);
     277    part[2] = positionDeployee.getComponent(1);
    267278    part[3] = begam.getComponent(1)/begamz;
    268     part[4] = pos.getComponent(2);
     279    part[4] = positionDeployee.getComponent(2);
    269280    part[5] = del;
    270281
     
    286297
    287298  // moyenne
    288   double facmoy = 1.0/double( goodPartic_.size() );
     299  double facmoy = 1.0/double( relativePartic_.size() );
    289300  for ( j = 0; j < 6; j++) {
    290301        ( matrice.at(j) ).at(j) = sqrt(( matrice.at(j) ).at(j) * facmoy );
     
    359370
    360371// renvoie (dans le vecteur coord) la liste des coordonnées d'index 'index' :
    361 // index = 0 , 1, 2 -> x,y,z
     372// index = 0 , 1, 2 -> x,y,z (en coordonnees deployees)
    362373// index = 3,4 -> x' = betax/betaz , y' = betay/betaz
    363374// index = 5 -> Ec : energie cinetique
     
    367378 
    368379  coord.clear();
    369   coord.resize(goodPartic_.size(), 0.0 );
     380  coord.resize(relativePartic_.size(), 0.0 );
    370381  cout << " particleBeam::particlesPhaseSpaceComponent index = " << index << endl;
    371382
    372383  if ( index <= 2 ) {
    373     for (unsigned i = 0; i < goodPartic_.size(); ++i) {
    374       coord.at(i) =  10.*goodPartic_.at(i).getPosition().getComponent(index);  // en mm
     384
     385    for (unsigned i = 0; i < relativePartic_.size(); ++i) {
     386
     387      coord.at(i) =  10.*coordonneesDeployees(i).getComponent(index);
     388      //      coord.at(i) =  10.*relativePartic_.at(i).getPosition().getComponent(index);  // en mm
    375389    }
    376390    return;
     
    378392 
    379393  if ( index >  2 && index < 5 ) {
    380     for (unsigned i = 0; i < goodPartic_.size(); ++i) {
    381       double begamz = goodPartic_.at(i).getBetaGamma().getComponent(2);
     394    for (unsigned i = 0; i < relativePartic_.size(); ++i) {
     395      double begamz = relativePartic_.at(i).getBetaGamma().getComponent(2);
    382396      if ( begamz != 0.0) {
    383         coord.at(i) =  1000.*goodPartic_.at(i).getBetaGamma().getComponent(index - 3)/begamz; // milliradians
     397        coord.at(i) =  1000.*relativePartic_.at(i).getBetaGamma().getComponent(index - 3)/begamz; // milliradians
    384398      } else {
    385399        coord.at(i) = 0.0;
     
    392406    double gamma0 = referenceParticle_.getGamma();
    393407    if ( gamma0 == 0.0 ) return;
    394     for (unsigned i = 0; i < goodPartic_.size(); ++i) {
    395       coord.at(i) =  100.*(goodPartic_.at(i).getGamma() - gamma0)/gamma0;  // en %
     408    for (unsigned i = 0; i < relativePartic_.size(); ++i) {
     409      coord.at(i) =  100.*(relativePartic_.at(i).getGamma() - gamma0)/gamma0;  // en %
    396410    }
    397411    return;
     
    618632  double ecatyp= 0.0;
    619633 
    620   for (unsigned int k = 0; k < goodPartic_.size(); k++) {
     634  for (unsigned int k = 0; k < relativePartic_.size(); k++) {
    621635    if (vshf[k] < vmin) vmin = vshf[k];
    622636    else if (vshf[k] > vmax) vmax = vshf[k];
     
    625639  }
    626640
    627   double sum= (float)goodPartic_.size();
     641  double sum= (float)relativePartic_.size();
    628642  out[0]= sum;
    629643  vmoy /= sum;
     
    654668  }
    655669
    656   for (unsigned int k = 0; k < goodPartic_.size(); k++) {
     670  for (unsigned int k = 0; k < relativePartic_.size(); k++) {
    657671    vshf[k] -= vmoy;
    658672  }
     
    706720}
    707721
     722 // coordonnees d'une particule dans le faisceau deploye ( passage
     723 // d'une representation z=cte a une representation t=cte)
     724TRIDVECTOR particleBeam::coordonneesDeployees(unsigned particleIndex, double cdtShift) {
     725    double gamma = relativePartic_.at(particleIndex).getGamma();
     726    //    TRIDVECTOR pos = relativePartic_.at(particleIndex).getPosition();
     727    double xx = relativePartic_.at(particleIndex).getReferenceToPosition().getComponent(0);
     728    double yy = relativePartic_.at(particleIndex).getReferenceToPosition().getComponent(1);
     729    double cdt = relativePartic_.at(particleIndex).getReferenceToPosition().getComponent(2);
     730    cdt += cdtShift;
     731    //    TRIDVECTOR begam= goodPartic_.at(particleIndex).getBetaGamma();
     732    double begamz = relativePartic_.at(particleIndex).getReferenceToBetaGamma().getComponent(2);
     733    double dz = begamz * cdt / gamma;
     734    double xp = relativePartic_.at(particleIndex).getReferenceToBetaGamma().getComponent(0)/begamz;
     735    double yp = relativePartic_.at(particleIndex).getReferenceToBetaGamma().getComponent(1)/begamz;
     736    xx += xp * dz;
     737    yy += yp * dz;   
     738    return TRIDVECTOR(xx, yy, dz);
     739 }
Note: See TracChangeset for help on using the changeset viewer.