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

Timestamp:

Mar 1, 2013, 3:38:12 PM (11 years ago)

Author:

lemeur

Message:

complement dessin phase space + rationalisation

File:

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

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

@@ -333 +333 @@
 }
 
-void particleBeam::donneesDessinEllipseXxp(vector<double>& xcor, vector<double>& ycor) {
+// void particleBeam::donneesDessinEllipseXxp(vector<double>& xcor, vector<double>& ycor) {
+//   int k;
+//   double x,y;
+
+//   if ( !momentRepresentationOk_ ) return;
+
+//   xcor.clear();
+//   ycor.clear();
+
+//   double xm = ( rij_.getMatrix().at(0) ).at(0);
+//   double ym = ( rij_.getMatrix().at(1) ).at(1);
+//   double r  = ( rij_.getMatrix().at(1) ).at(0);
+
+//   cout <<  " racs11= " << xm << " racs22= " << ym << " r12= " << r << endl;
+
+
+//   int nbintv = 50;
+//   if ( xm == 0.0 ) return;
+//   double pas = 2.0 * xm / nbintv;
+
+//   //  cout << " r= " << r << endl;
+//   double rac = (1 - r*r);
+//   if ( rac > 0.0 ) 
+//     {
+//       cout << " cas rac > " << endl;
+//       rac = sqrt(1 - r*r);
+//       double alpha = -r / rac;
+//       double beta = xm / ( ym * rac);
+//       //  double gamma = ym / ( xm * rac );
+//       double epsil = xm * ym * rac;
+//       double fac1 = -1.0 / ( beta * beta);
+//       double fac2 = epsil/beta;
+//       double fac3 = -alpha/beta;
+//       double aux;
+//       for ( k=0; k < nbintv; k++)
+//      {
+//        x = -xm + k*pas;
+//        aux = fac1 * x * x + fac2;
+//        //     cout << " aux2= " << aux << endl;
+//        if ( aux <= 0.0 )
+//          {
+//            aux = 0.0;
+//          }
+//        else aux = sqrt(aux);
+      
+//        //        y = fac3*x;
+//        y = fac3*x + aux;
+//        xcor.push_back(x);
+//        ycor.push_back(y);
+//      }
+
+//       for ( k=0; k <= nbintv; k++)
+//      {
+//        x = xm - k*pas;
+//        aux =  fac1 * x * x + fac2;
+//        if ( aux <= 0.0 ) 
+//          {
+//            aux = 0.0;
+//          }
+//        else aux = sqrt(aux);
+//        //   y = fac3*x;
+//        y = fac3*x - aux;
+//        xcor.push_back(x);
+//        ycor.push_back(y);
+//      }
+//     }
+//   else
+//     // cas degenere
+//     {
+//       cout << " cas degenere " << endl;
+//       double fac = ym/xm;
+//       for ( k=0; k < nbintv; k++)
+//      {
+//        x = -xm + k*pas;
+//        y = fac*x;
+//        xcor.push_back(x);
+//        ycor.push_back(y);
+//      }
+        
+//     }
+// }
+
+void particleBeam::particlesPhaseSpaceData(vector<double>& xcor, vector<double>& ycor, unsigned indexAbs, unsigned indexOrd) {
+  particlesPhaseSpaceComponent(xcor, indexAbs);
+  particlesPhaseSpaceComponent(ycor, indexOrd);
+}
+
+void particleBeam::particlesPhaseSpaceComponent(vector<double>& coord, unsigned index) {
+  if ( !particleRepresentationOk_ ) return;
+  coord.clear();
+  coord.resize(goodPartic_.size(), 0.0 );
+  cout << " particleBeam::particlesPhaseSpaceComponent index = " << index << endl;
+  if ( index <= 2 ) {
+    for (unsigned i = 0; i < goodPartic_.size(); ++i) {
+      coord.at(i) =  goodPartic_.at(i).getPosition().getComponent(index);
+    }
+    return;
+  }
+
+  if ( index >  2 && index < 5 ) {
+    for (unsigned i = 0; i < goodPartic_.size(); ++i) {
+      double begamz = goodPartic_.at(i).getBetaGamma().getComponent(2);
+      if ( begamz != 0.0) {
+        coord.at(i) =  1000.*goodPartic_.at(i).getBetaGamma().getComponent(index - 3)/begamz; // mimmiradians
+      } else {
+        coord.at(i) = 0.0;
+      }
+    }
+    return;
+  }
+
+  if ( index == 5 ) {
+    double gamma0 = referenceParticle_.getGamma();
+    cout << " gamma0 = " << gamma0 << endl;
+    if ( gamma0 == 0.0 ) return;
+    for (unsigned i = 0; i < goodPartic_.size(); ++i) {
+      coord.at(i) =  100.*(goodPartic_.at(i).getGamma() - gamma0)/gamma0;  // en %
+      cout << " gamma0 = " << gamma0 << " gamma = " << goodPartic_.at(i).getGamma() << endl;
+    }
+    return;
+  }
+}
+
+void particleBeam::donneesDessinEllipse(vector<double>& xcor, vector<double>& ycor, unsigned indexAbs, unsigned indexOrd) {
   int k;
   double x,y;
@@ -339 +462 @@
   if ( !momentRepresentationOk_ ) return;
 
+  if ( indexAbs > 5 || indexOrd > 5 ) return;
+
   xcor.clear();
   ycor.clear();
-
-  double xm = ( rij_.getMatrix().at(0) ).at(0);
-  double ym = ( rij_.getMatrix().at(1) ).at(1);
-  double r  = ( rij_.getMatrix().at(1) ).at(0);
+  // les index sont dans l'ordre x,y,z,xp,yp, de/E
+  // on traduit en TRANSPORT
+  if ( indexAbs == 1 ) indexAbs = 2;  // y
+  if ( indexAbs == 2 ) indexAbs = 4;  // z -> l
+  if ( indexAbs == 3 ) indexAbs = 1;  // xp
+  if ( indexAbs == 4 ) indexAbs = 3;  // yp
+
+  if ( indexOrd == 1 ) indexOrd = 2;  // y
+  if ( indexOrd == 2 ) indexOrd = 4;  // z -> l
+  if ( indexOrd == 3 ) indexOrd = 1;  // xp
+  if ( indexOrd == 4 ) indexOrd = 3;  // yp
+
+  cout << " index x" << indexAbs << " index y " << indexOrd << endl;
+
+  double xm = ( rij_.getMatrix().at(indexAbs) ).at(indexAbs);
+  double ym = ( rij_.getMatrix().at(indexOrd) ).at(indexOrd);
+  double r;
+  if ( indexOrd > indexAbs ) {
+    r  = ( rij_.getMatrix().at(indexOrd) ).at(indexAbs);
+  } else {
+    r  = ( rij_.getMatrix().at(indexAbs) ).at(indexOrd);
+  }
 
   cout <<  " racs11= " << xm << " racs22= " << ym << " r12= " << r << endl;
@@ -415 +558 @@
 }
 
-
-
-void particleBeam::donneesDessinEllipse(vector<double>& xcor, vector<double>& ycor, unsigned indexAbs, unsigned indexOrd) {
-  int k;
-  double x,y;
-
-  if ( !momentRepresentationOk_ ) return;
-
-  if ( indexAbs > 5 || indexOrd > 5 ) return;
-
-  xcor.clear();
-  ycor.clear();
-  // les index sont dans l'ordre x,y,z,xp,yp, de/E
-  // on traduit en TRANSPORT
-  if ( indexAbs == 1 ) indexAbs = 2;  // y
-  if ( indexAbs == 2 ) indexAbs = 4;  // z -> l
-  if ( indexAbs == 3 ) indexAbs = 1;  // xp
-  if ( indexAbs == 4 ) indexAbs = 3;  // yp
-
-  if ( indexOrd == 1 ) indexOrd = 2;  // y
-  if ( indexOrd == 2 ) indexOrd = 4;  // z -> l
-  if ( indexOrd == 3 ) indexOrd = 1;  // xp
-  if ( indexOrd == 4 ) indexOrd = 3;  // yp
-
-  cout << " index x" << indexAbs << " index y " << indexOrd << endl;
-
-  double xm = ( rij_.getMatrix().at(indexAbs) ).at(indexAbs);
-  double ym = ( rij_.getMatrix().at(indexOrd) ).at(indexOrd);
-  double r;
-  if ( indexOrd > indexAbs ) {
-    r  = ( rij_.getMatrix().at(indexOrd) ).at(indexAbs);
-  } else {
-    r  = ( rij_.getMatrix().at(indexAbs) ).at(indexOrd);
-  }
-
-  cout <<  " racs11= " << xm << " racs22= " << ym << " r12= " << r << endl;
-
-
-  int nbintv = 50;
-  if ( xm == 0.0 ) return;
-  double pas = 2.0 * xm / nbintv;
-
-  //  cout << " r= " << r << endl;
-  double rac = (1 - r*r);
-  if ( rac > 0.0 ) 
-    {
-      cout << " cas rac > " << endl;
-      rac = sqrt(1 - r*r);
-      double alpha = -r / rac;
-      double beta = xm / ( ym * rac);
-      //  double gamma = ym / ( xm * rac );
-      double epsil = xm * ym * rac;
-      double fac1 = -1.0 / ( beta * beta);
-      double fac2 = epsil/beta;
-      double fac3 = -alpha/beta;
-      double aux;
-      for ( k=0; k < nbintv; k++)
-        {
-          x = -xm + k*pas;
-          aux = fac1 * x * x + fac2;
-          //     cout << " aux2= " << aux << endl;
-          if ( aux <= 0.0 )
-            {
-              aux = 0.0;
-            }
-          else aux = sqrt(aux);
-      
-          //        y = fac3*x;
-          y = fac3*x + aux;
-          xcor.push_back(x);
-          ycor.push_back(y);
-        }
-
-      for ( k=0; k <= nbintv; k++)
-        {
-          x = xm - k*pas;
-          aux =  fac1 * x * x + fac2;
-          if ( aux <= 0.0 ) 
-            {
-              aux = 0.0;
-            }
-          else aux = sqrt(aux);
-          //   y = fac3*x;
-          y = fac3*x - aux;
-          xcor.push_back(x);
-          ycor.push_back(y);
-        }
-    }
-  else
-    // cas degenere
-    {
-      cout << " cas degenere " << endl;
-      double fac = ym/xm;
-      for ( k=0; k < nbintv; k++)
-        {
-          x = -xm + k*pas;
-          y = fac*x;
-          xcor.push_back(x);
-          ycor.push_back(y);
-        }
-        
-    }
-}
-
 void particleBeam::histogramme(vector<double>&xcor,vector<int>& hist,int& cnts,double out[3])
 {