Changeset 417 in PSPA for Interface_Web/trunk/pspaWT/sources/controler/src/particleBeam.cc

Timestamp:

Jun 10, 2013, 4:52:52 PM (11 years ago)

Author:

lemeur

Message:

mise a jour des manipulations de faisceau

File:

• Interface_Web/trunk/pspaWT/sources/controler/src/particleBeam.cc (modified) (16 diffs)

Interface_Web/trunk/pspaWT/sources/controler/src/particleBeam.cc

@@ -19 +19 @@
 
 void particleBeam::clear() {
-  goodPartic_.clear();
+  relativePartic_.clear();
   rij_.raz();
   P0Transport_ = 0.0;
@@ -27 +27 @@
 
 int particleBeam::getNbParticles() const {
-  return goodPartic_.size();
+  return relativePartic_.size();
 }
 
@@ -99 +99 @@
   centroid_ = centroid;
   referenceParticle_ = referencePart;
-  goodPartic_.clear();
-  goodPartic_ = particles;
-  cout << " particleBeam::setWithParticles taille vect. part. ENREGISTRE " << goodPartic_.size() << endl;
+  relativePartic_.clear();
+  relativePartic_ = particles;
+  cout << " particleBeam::setWithParticles taille vect. part. ENREGISTRE " << relativePartic_.size() << endl;
   particleRepresentationOk_ = true;
 }
@@ -114 +114 @@
 void  particleBeam::addParticle( bareParticle p)
 {
-  goodPartic_.push_back(p);
+  relativePartic_.push_back(p);
 }
 
 const vector<bareParticle>& particleBeam::getParticleVector() const
 {
-  return goodPartic_;
+  return relativePartic_;
 }
 
 vector<bareParticle>& particleBeam::getParticleVector() 
 {
-  return goodPartic_;
-}
-
-void particleBeam::getVariance(double& varx, double& vary, double& varz) const {
-  unsigned int k;
-  double x,y,z;
-  double xav = 0.;
-  double yav = 0.;
-  double zav = 0.;
-  double xavsq = 0.;
-  double yavsq = 0.;
-  double zavsq = 0.;
-
-  TRIDVECTOR pos;
-
-
-  for ( k = 0 ; k < goodPartic_.size(); k++) {
-    pos = goodPartic_.at(k).getPosition();
-    pos.getComponents(x,y,z);
-    //      partic_[k].getXYZ(x,y,z);
-    xav += x;
-    xavsq += x*x;
-    yav += y;
-    yavsq += y*y;
-    zav += z;
-    zavsq += z*z;
-  }
-
-  double aginv = double (goodPartic_.size());
-  aginv = 1.0/aginv;
-
-  varx =  aginv * ( xavsq - xav*xav*aginv );
-  vary =  aginv * ( yavsq - yav*yav*aginv );
-  varz =  aginv * ( zavsq - zav*zav*aginv );
-}
+  return relativePartic_;
+}
+
+// void particleBeam::getVariance(double& varx, double& vary, double& varz) const {
+//   unsigned int k;
+//   double x,y,z;
+//   double xav = 0.;
+//   double yav = 0.;
+//   double zav = 0.;
+//   double xavsq = 0.;
+//   double yavsq = 0.;
+//   double zavsq = 0.;
+
+//   TRIDVECTOR pos;
+
+
+//   for ( k = 0 ; k < goodPartic_.size(); k++) {
+//     pos = goodPartic_.at(k).getPosition();
+//     pos.getComponents(x,y,z);
+//     //      partic_[k].getXYZ(x,y,z);
+//     xav += x;
+//     xavsq += x*x;
+//     yav += y;
+//     yavsq += y*y;
+//     zav += z;
+//     zavsq += z*z;
+//   }
+
+//   double aginv = double (goodPartic_.size());
+//   aginv = 1.0/aginv;
+
+//   varx =  aginv * ( xavsq - xav*xav*aginv );
+//   vary =  aginv * ( yavsq - yav*yav*aginv );
+//   varz =  aginv * ( zavsq - zav*zav*aginv );
+// }
 
 
 void particleBeam::printAllXYZ() const {
   cout << " dump du faisceau : " << endl;
-  cout <<  goodPartic_.size() << " particules " << endl;
+  cout <<  relativePartic_.size() << " particules " << endl;
   unsigned int k;
-  for ( k = 0 ; k < goodPartic_.size(); k++)
+  for ( k = 0 ; k < relativePartic_.size(); k++)
     {
       double xx,yy,zz;
-      goodPartic_.at(k).getPosition().getComponents(xx,yy,zz);
+      relativePartic_.at(k).getPosition().getComponents(xx,yy,zz);
       double betgamx, betgamy, betgamz;
-      goodPartic_.at(k).getBetaGamma().getComponents(betgamx, betgamy, betgamz);
-      cout << " part. numero " << k << "  x= " << xx << " y= " << yy  << " z= " << zz << "  betgamx= " << betgamx << " betgamy= " << betgamy  << " betgamz= " << betgamz << endl;
-    }
-}
-
-
-
-void particleBeam::Zrange(double& zmin, double& zmax) const {
+      relativePartic_.at(k).getBetaGamma().getComponents(betgamx, betgamy, betgamz);
+      cout << " part. numero " << k << "  x= " << xx << " y= " << yy  << " dphas (c.dt, cm) = " << zz << "  betgamx= " << betgamx << " betgamy= " << betgamy  << " betgamz= " << betgamz << endl;
+    }
+}
+
+
+// extension en phase (cm)
+void particleBeam::ZrangeCdt(double& cdtmin, double& cdtmax) const {
   double z;
-  zmin = GRAND;
-  zmax = -zmin;
+  cdtmin = GRAND;
+  cdtmax = -cdtmin;
 
   unsigned int k;
-  for ( k = 0 ; k < goodPartic_.size(); k++)
+  for ( k = 0 ; k < relativePartic_.size(); k++)
     {
-      z = goodPartic_.at(k).getZ();
-      if ( z < zmin ) zmin = z;
-      else if ( z > zmax) zmax = z;         
-    }
-}
-
-
-
-string particleBeam::FileOutputFlow() const {
+      z = relativePartic_.at(k).getZ(); // ce z est un dephasage, c.dt, en cm
+      if ( z < cdtmin ) cdtmin = z;
+      else if ( z > cdtmax) cdtmax = z;     
+    }
+}
+
+
+
+string particleBeam::fileOutputFlow() const {
   ostringstream sortie;
   unsigned int k;
-  for ( k = 0 ; k < goodPartic_.size(); k++)
+  for ( k = 0 ; k < relativePartic_.size(); k++)
     {
-      sortie << goodPartic_.at(k).FileOutputFlow() << endl;
+      sortie << relativePartic_.at(k).FileOutputFlow() << endl;
     }
   sortie << endl;
@@ -219 +219 @@
 
 void particleBeam::buildMomentRepresentation() {
-
+  // le faisceau cense etre represente en un z donne (reference) sera 
+  // ici deploye spatialement (faisceau en un temps donne)
   unsigned k,j,m;
   double auxj, auxm;
@@ -248 +249 @@
   double PMeVsc;
   double del;
+  //  double xp,yp,dz,cdt;
   vector<double> part(6, 0.0);
 
     vector< vector<double> >& matrice = rij_.getMatrix();
 
-
-  for (k=0; k < goodPartic_.size(); k++) {
-    gamma = goodPartic_.at(k).getGamma();
-    pos = goodPartic_.at(k).getPosition();
-    begam= goodPartic_.at(k).getBetaGamma();
+    TRIDVECTOR positionDeployee;
+
+  for (k=0; k < relativePartic_.size(); k++) {
+
+    positionDeployee = coordonneesDeployees(k);
+    gamma = relativePartic_.at(k).getGamma();
+
+    //    pos = relativePartic_.at(k).getPosition();
+    //    cdt = pos.getComponent(2);
+    begam= relativePartic_.at(k).getBetaGamma();
     begamz = begam.getComponent(2);
     g = gamma -1.0;
@@ -262 +269 @@
     del = 100.0 * ( PMeVsc -  P_reference_MeV_sur_c ) / P_reference_MeV_sur_c ; // en %
 
-    part[0] = pos.getComponent(0);
+    //    dz = begamz * cdt / gamma;
+    //    xp = begam.getComponent(0)/begamz;
+    //    yp = begam.getComponent(1)/begamz;
+
+    part[0] = positionDeployee.getComponent(0);
     part[1] = begam.getComponent(0)/begamz;
-    part[2] = pos.getComponent(1);
+    part[2] = positionDeployee.getComponent(1);
     part[3] = begam.getComponent(1)/begamz;
-    part[4] = pos.getComponent(2);
+    part[4] = positionDeployee.getComponent(2);
     part[5] = del;
 
@@ -286 +297 @@
 
   // moyenne
-  double facmoy = 1.0/double( goodPartic_.size() );
+  double facmoy = 1.0/double( relativePartic_.size() );
   for ( j = 0; j < 6; j++) {
         ( matrice.at(j) ).at(j) = sqrt(( matrice.at(j) ).at(j) * facmoy );
@@ -359 +370 @@
 
 // renvoie (dans le vecteur coord) la liste des coordonnées d'index 'index' : 
-// index = 0 , 1, 2 -> x,y,z
+// index = 0 , 1, 2 -> x,y,z (en coordonnees deployees)
 // index = 3,4 -> x' = betax/betaz , y' = betay/betaz
 // index = 5 -> Ec : energie cinetique
@@ -367 +378 @@
   
   coord.clear();
-  coord.resize(goodPartic_.size(), 0.0 );
+  coord.resize(relativePartic_.size(), 0.0 );
   cout << " particleBeam::particlesPhaseSpaceComponent index = " << index << endl;
 
   if ( index <= 2 ) {
-    for (unsigned i = 0; i < goodPartic_.size(); ++i) {
-      coord.at(i) =  10.*goodPartic_.at(i).getPosition().getComponent(index);  // en mm
+
+    for (unsigned i = 0; i < relativePartic_.size(); ++i) {
+
+      coord.at(i) =  10.*coordonneesDeployees(i).getComponent(index);
+      //      coord.at(i) =  10.*relativePartic_.at(i).getPosition().getComponent(index);  // en mm
     }
     return;
@@ -378 +392 @@
   
   if ( index >  2 && index < 5 ) {
-    for (unsigned i = 0; i < goodPartic_.size(); ++i) {
-      double begamz = goodPartic_.at(i).getBetaGamma().getComponent(2);
+    for (unsigned i = 0; i < relativePartic_.size(); ++i) {
+      double begamz = relativePartic_.at(i).getBetaGamma().getComponent(2);
       if ( begamz != 0.0) {
-        coord.at(i) =  1000.*goodPartic_.at(i).getBetaGamma().getComponent(index - 3)/begamz; // milliradians
+        coord.at(i) =  1000.*relativePartic_.at(i).getBetaGamma().getComponent(index - 3)/begamz; // milliradians
       } else {
         coord.at(i) = 0.0;
@@ -392 +406 @@
     double gamma0 = referenceParticle_.getGamma();
     if ( gamma0 == 0.0 ) return;
-    for (unsigned i = 0; i < goodPartic_.size(); ++i) {
-      coord.at(i) =  100.*(goodPartic_.at(i).getGamma() - gamma0)/gamma0;  // en %
+    for (unsigned i = 0; i < relativePartic_.size(); ++i) {
+      coord.at(i) =  100.*(relativePartic_.at(i).getGamma() - gamma0)/gamma0;  // en %
     }
     return;
@@ -618 +632 @@
   double ecatyp= 0.0;
   
-  for (unsigned int k = 0; k < goodPartic_.size(); k++) {
+  for (unsigned int k = 0; k < relativePartic_.size(); k++) {
     if (vshf[k] < vmin) vmin = vshf[k];
     else if (vshf[k] > vmax) vmax = vshf[k];
@@ -625 +639 @@
   }
 
-  double sum= (float)goodPartic_.size();
+  double sum= (float)relativePartic_.size();
   out[0]= sum; 
   vmoy /= sum;
@@ -654 +668 @@
   }
 
-  for (unsigned int k = 0; k < goodPartic_.size(); k++) {
+  for (unsigned int k = 0; k < relativePartic_.size(); k++) {
     vshf[k] -= vmoy;
   }
@@ -706 +720 @@
 }
 
+ // coordonnees d'une particule dans le faisceau deploye ( passage 
+ // d'une representation z=cte a une representation t=cte)
+TRIDVECTOR particleBeam::coordonneesDeployees(unsigned particleIndex, double cdtShift) {
+    double gamma = relativePartic_.at(particleIndex).getGamma();
+    //    TRIDVECTOR pos = relativePartic_.at(particleIndex).getPosition();
+    double xx = relativePartic_.at(particleIndex).getReferenceToPosition().getComponent(0);
+    double yy = relativePartic_.at(particleIndex).getReferenceToPosition().getComponent(1);
+    double cdt = relativePartic_.at(particleIndex).getReferenceToPosition().getComponent(2);
+    cdt += cdtShift;
+    //    TRIDVECTOR begam= goodPartic_.at(particleIndex).getBetaGamma();
+    double begamz = relativePartic_.at(particleIndex).getReferenceToBetaGamma().getComponent(2);
+    double dz = begamz * cdt / gamma;
+    double xp = relativePartic_.at(particleIndex).getReferenceToBetaGamma().getComponent(0)/begamz;
+    double yp = relativePartic_.at(particleIndex).getReferenceToBetaGamma().getComponent(1)/begamz;
+    xx += xp * dz;
+    yy += yp * dz;    
+    return TRIDVECTOR(xx, yy, dz);
+ }