Changeset 356 in PSPA

Timestamp:

Mar 2, 2013, 7:54:38 PM (12 years ago)

Author:

lemeur

Message:

bugs dans les dessins

Location:

Interface_Web/trunk/pspaWT/sources

Files:

• controler/include/PhysicalConstants.h (modified) (1 diff)
• controler/include/particleBeam.h (modified) (2 diffs)
• controler/src/elementBend.cc (modified) (1 diff)
• controler/src/particleBeam.cc (modified) (9 diffs)
• controler/src/softwareGenerator.cc (modified) (3 diffs)
• controler/src/softwareParmela.cc (modified) (1 diff)
• userInterface/src/GWt_pspaApplication.cc (modified) (4 diffs)

Interface_Web/trunk/pspaWT/sources/controler/include/PhysicalConstants.h

@@ -8 +8 @@
 #define ELECTRONANOCOULOMB 1.60217653e-10
 #define ELECTRONRADIUS 2.817940325e-13  // cm!
-#define ERESTMeV 0.510998918
-#define ERESTeV 0.510998918e+06
+#define EREST_MeV 0.510998918
+#define EREST_keV 0.510998918e+03
+#define EREST_eV 0.510998918e+06
 #define CLIGHT_E8 2.99792458
 #define CLIGHT_m_per_ns 0.299792458

Interface_Web/trunk/pspaWT/sources/controler/include/particleBeam.h

@@ -41 +41 @@
   void razDesMoments();
   void particlesPhaseSpaceComponent(vector<double>& coord, unsigned index);
+  unsigned indexFromPspaToTransport(unsigned index) const;
 
  public:
@@ -67 +68 @@
   void Zrange(double& zmin, double& zmax) const;
   //  void donneesDessinEllipseXxp(vector<double>& xcor, vector<double>& ycor);
-  void donneesDessinEllipse(vector<double>& xcor, vector<double>& ycor, unsigned indexAbs, unsigned indexOrd);
+  void donneesDessinEllipse(vector<double>& xcor, vector<double>& ycor, vector<string>& legende, unsigned indexAbs, unsigned indexOrd);
   void particlesPhaseSpaceData(vector<double>& xcor, vector<double>& ycor, unsigned indexAbs, unsigned indexOrd);
   virtual string FileOutputFlow() const;

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

@@ -96 +96 @@
     ostringstream sortie;
     // il faut entrer l'energie cinetique
-    double ecin = momentum_/ERESTMeV;
+    double ecin = momentum_/EREST_MeV;
     ecin = ecin*ecin + 1.;
-    ecin = ERESTMeV*(sqrt(ecin) - 1.0);
+    ecin = EREST_MeV*(sqrt(ecin) - 1.0);
     sortie << "BEND /l=" << lenghtElem_ << "  / aper=" << aperture_ << "  / iout=1 / wr=" << ecin << " /alpha=" << angleDeg_ << " / beta1=" << beta1_;
     sortie << " / beta2=" << beta1_;

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

@@ -4 +4 @@
 #include "PhysicalConstants.h"
 #include "mathematicalTools.h"
+#include "mixedTools.h"
 
 #include <stdio.h>
@@ -70 +71 @@
 double particleBeam::referenceKineticEnergyMeV() const {
   if ( particleRepresentationOk_ ) {
-    return (referenceParticle_.getGamma() -1.) * ERESTMeV;
+    return (referenceParticle_.getGamma() -1.) * EREST_MeV;
   } else {
-    double P0Norm = 1000.0 * P0Transport_ / ERESTMeV;
+    double P0Norm = 1000.0 * P0Transport_ / EREST_MeV;
     double gamma = sqrt(1.0 +  P0Norm * P0Norm);
-    return (gamma - 1.0) * ERESTMeV;
+    return (gamma - 1.0) * EREST_MeV;
   }
 }
@@ -302 +303 @@
   ( matrice.at(1) ).at(1)  *= uniteAngle;
   ( matrice.at(3) ).at(3)  *= uniteAngle;
-  P0Transport_ = 1.0e-3*ERESTMeV*P_reference_MeV_sur_c;
+  P0Transport_ = 1.0e-3*EREST_MeV*P_reference_MeV_sur_c;
 
   //  cout << " buildmomentrepresentation impression des moments " << endl;
@@ -365 +366 @@
   if ( index == 5 ) {
     double gamma0 = referenceParticle_.getGamma();
-    cout << " gamma0 = " << gamma0 << endl;
+    //    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;
+      //      cout << " gamma0 = " << gamma0 << " gamma = " << goodPartic_.at(i).getGamma() << endl;
     }
     return;
   }
 }
-
-void particleBeam::donneesDessinEllipse(vector<double>& xcor, vector<double>& ycor, unsigned indexAbs, unsigned indexOrd) {
+unsigned particleBeam::indexFromPspaToTransport(unsigned index) const {
+  cout << " indexFromPspaToTransport entree : " << index << endl;
+  unsigned indexRetour;
+  switch ( index ) {
+  case 0 : return  0; break; // x
+  case 1 : return  2; break; // y
+  case 2 : return  4; break; // z -> l
+  case 3 : return  1; break; // xp
+  case 4 : return  3; break; // yp
+  case 5 : return  5; break; // de/E
+  default : { 
+    cout << " particleBeam::indexFromPspaToTransport : coordinate index ERROR inital index :  "<< index << endl;
+    return 99;
+  }
+  }
+}
+
+void particleBeam::donneesDessinEllipse(vector<double>& xcor, vector<double>& ycor, vector<string>& legende, unsigned indexAbs, unsigned indexOrd) {
   int k;
   double x,y;
+  cout << " particleBeam::donneesDessinEllipse index recus x" << indexAbs << " index y " << indexOrd << endl;
 
   if ( !momentRepresentationOk_ ) return;
@@ -387 +405 @@
   // 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
+
+  indexAbs = indexFromPspaToTransport(indexAbs);
+  indexOrd = indexFromPspaToTransport(indexOrd);
+
+  // a completer 
+  legende.clear();
+  if ( indexAbs == 0 && indexOrd == 1 ) {
+    string namx = " x ";
+    string namy = " xp ";
+    string  emitt = mixedTools::doubleToString(getUnnormalizedEmittanceX());
+    string xmax = "x max= ";
+    string valXmax = mixedTools::doubleToString(10.0*getSigmaTransportij(1,1)); // mm
+    string ymax = "x' max= ";
+    string valYmax = mixedTools::doubleToString(getSigmaTransportij(2,2)); // mm
+    string correl = " correlation ";
+    string valCorrel = mixedTools::doubleToString(getSigmaTransportij(2,1));
+    string xunit = " mm";
+    string yunit = " mrad";
+    legende.push_back( "emittance" + namx + "," + namy + ": " + emitt + " pi." + xunit + "." + yunit);
+    legende.push_back( xmax + valXmax + xunit);
+    legende.push_back( ymax + valYmax + yunit);
+  } else {
+    legende.push_back(" text of legend not yet programmed ");
+  }
 
   cout << " index x" << indexAbs << " index y " << indexOrd << endl;
@@ -481 +514 @@
   // sortie pour la legende: out[0]= entries, out[1]= mean, out[2]= rms
 
+  // on calcule avec les gammas ; les MeV c'est pour les sorties
+
   double gammin= GRAND;
   double gammax= -gammin;
-  double Emoy= 0.0;
+  //  double Emoy= 0.0;
   double ecatyp= 0.0;
-
+  double gmoy=0.0;
   for (unsigned int k = 0; k < goodPartic_.size(); k++)
     {
       double gamma = goodPartic_.at(k).getGamma();
-      double EMev = (gamma-1.0)*ERESTMeV;
+      //      double EMev = (gamma-1.0)*ERESTMeV;
       if (gamma < gammin) gammin = gamma;
       else if (gamma > gammax) gammax = gamma;
-      Emoy += EMev;
-      ecatyp += EMev*EMev;
+      //      Emoy += EMev;
+      double g = gamma-1.0;
+      gmoy += g;
+      //      ecatyp += EMev*EMev;
+      ecatyp += g*g;
     }
 
   double sum= (float)goodPartic_.size();
   out[0]= sum;
-  Emoy /= sum;
-  out[1]= Emoy; //MeV
+  //  Emoy /= sum;
+  gmoy /= sum;
+  //  out[1]= Emoy; //MeV
+  out[1]= gmoy*EREST_MeV ; //MeV
   ecatyp /= sum;
-  out[2]= 1000.0*sqrt(ecatyp); //KeV
-  ecatyp = sqrt(abs(ecatyp-Emoy*Emoy));
-
-  double Emin = (gammin-1.0)*ERESTMeV;
-  double Emax = (gammax-1.0)*ERESTMeV;
-  cout << "energie cinetique -moyenne " << Emoy << " Mev " << "-mini " << Emin << " Mev " << "-maxi " << Emax << " Mev " << "ecart type " << ecatyp*1000.0 << " Kev" << endl;
+  //  out[2]= 1000.0*sqrt(ecatyp); //KeV
+  //  ecatyp = sqrt(abs(ecatyp-Emoy*Emoy));
+
+  ecatyp = sqrt(abs(ecatyp-gmoy*gmoy));
+  out[2]= ecatyp*EREST_keV; //KeV
+
+  double gmin = gammin-1.0;
+  double gmax = gammax-1.0;
+  cout << "energie cinetique -moyenne " << gmoy*EREST_MeV << " Mev " << "-mini " << gmin*EREST_MeV << " Mev " << "-maxi " << gmax*EREST_MeV << " Mev " << "ecart type " << ecatyp*EREST_keV << " Kev" << endl;
 
   vector<double> Eshf;
@@ -512 +555 @@
     {
       double gamma = goodPartic_.at(k).getGamma();
-      double EMev = (gamma-1.0)*ERESTMeV;
-      Eshf.push_back(EMev-Emoy);
+      // double EMev = (gamma-1.0)*ERESTMeV;
+      // Eshf.push_back(EMev-Emoy);
+      // modif glm 2/03/2013
+      double g = (gamma-1.0);
+      Eshf.push_back( (g-gmoy) );
     }
 
@@ -542 +588 @@
 
     // on gradue l'abscisse en pourcents
-    xcor[i]= 100.*( vmin+i*hpas );
+    //   xcor[i]= 100.*( vmin+i*hpas );
+    xcor[i]= 100.*( vmin+i*hpas )/gmoy;
   }
 
@@ -548 +595 @@
   for (unsigned i = 0; i < Eshf.size(); ++i) {
     double var= Eshf[i]-vmin;
-    if(var < 0 || var >= phist) cout<<"out of range "<<var<<", ("<< i<<")"<< endl;
-    int k= var/hpas;
+    //    if(var < 0 || var >= phist) cout<<"out of range "<<var<<", ("<< i<<")"<< endl;
+    if(var < 0 ) {
+      cout<<"lesser that the minimum "<<var<<", ("<< i<<")"<< endl;
+      hist[0]++;
+    }
+    if(var >= phist) {
+      cout<<"greater that the maximum "<<var<<", ("<< i<<")"<< endl;
+      hist[ihist-1]++;
+    }
+    //    int k= var/hpas;
     int kk= (int)floor(var/hpas);
     //if(i%20 == 0) cout<<"v("<<i<<")= " <<var<<" ["<<k<<"-"<<kk<<"], "<<endl; 

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

@@ -140 +140 @@
       timeRef = partic.clock;
       TRIDVECTOR  posRef(100.*partic.xx,100.*partic.yy,100.*partic.zz); // en cm
-      betagammaZRef = partic.pz/ERESTeV;
+      betagammaZRef = partic.pz/EREST_eV;
       // l'impulsion donnee par generator est en eV/c
-      TRIDVECTOR betagammaRef(partic.px/ERESTeV , partic.py/ERESTeV, betagammaZRef);
+      TRIDVECTOR betagammaRef(partic.px/EREST_eV , partic.py/EREST_eV, betagammaZRef);
       refPart = bareParticle(posRef, betagammaRef);
 
@@ -190 +190 @@
 
     // tout ce qui suit sera a clarifier
-    double betaGammax = pxPart/ERESTeV;
-    double betaGammay = pyPart/ERESTeV;
-    double betaGammaz = betagammaZRef + pzPartRel/ERESTeV;
+    double betaGammax = pxPart/EREST_eV;
+    double betaGammay = pyPart/EREST_eV;
+    double betaGammaz = betagammaZRef + pzPartRel/EREST_eV;
     betagamma.setComponents(betaGammax, betaGammay, betaGammaz);
     double gamma = sqrt(1.0 + betagamma.norm2());
@@ -203 +203 @@
 
     // ici on neglige la difference entre gamma de la part. et gamma de la ref.
-    deltaz = (pzPartRel/ERESTeV) * CLIGHT_m_per_ns * deltat; // en metres
+    deltaz = (pzPartRel/EREST_eV) * CLIGHT_m_per_ns * deltat; // en metres
 
     pos.setComponents(100.*x,100.*y,100.*deltaz);  // en cm

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

@@ -243 +243 @@
     // dephasage par rapport a la reference  
     dephas = faisceau.at(k).phi - phaseRef; // degrés
-    g = faisceau.at(k).wz/ERESTMeV;
+    g = faisceau.at(k).wz/EREST_MeV;
     betagammaz = faisceau.at(k).begamz;
     betaz = betagammaz/(g+1.0);

Interface_Web/trunk/pspaWT/sources/userInterface/src/GWt_pspaApplication.cc

@@ -502 +502 @@
     {
       if ( !beam->momentRepresentationOk() ) beam->buildMomentRepresentation();
+      cout << " PspaApplication::dessinerPhaseSpace cood cliquees " << xabs << " " << yord << endl;
       faireDessinTransport(toto_, beam, xabs, yord, nameAbs, nameOrd );
     }
@@ -582 +583 @@
   ellipseDialog->setClosable(true);
   ellipseDialog->show();
-    
-  new WText(nameOfCase_ + " : emittance " + namex + " , " + namey + " : " + mixedTools::doubleToString(beam->getUnnormalizedEmittanceX()) + " pi.mm.mrad" , ellipseDialog->contents());
-  new WBreak(ellipseDialog->contents());
-  new WText(" xmax = " + mixedTools::doubleToString(beam->getSigmaTransportij(1,1)) + " cm ", ellipseDialog->contents());
-  new WBreak(ellipseDialog->contents());
-  new WText(" x' max = " + mixedTools::doubleToString(beam->getSigmaTransportij(2,2)) + " mrad ", ellipseDialog->contents());
-  new WBreak(ellipseDialog->contents());
-  new WText(" corr. 12 = " + mixedTools::doubleToString(beam->getSigmaTransportij(2,1))  + " sqrt(cm.rad) ", ellipseDialog->contents());
+  new WText(nameOfCase_, ellipseDialog->contents());
+  // new WText(nameOfCase_ + " : emittance " + namex + " , " + namey + " : " + mixedTools::doubleToString(beam->getUnnormalizedEmittanceX()) + " pi.mm.mrad" , ellipseDialog->contents());
+  // new WBreak(ellipseDialog->contents());
+  // new WText(" xmax = " + mixedTools::doubleToString(beam->getSigmaTransportij(1,1)) + " cm ", ellipseDialog->contents());
+  // new WBreak(ellipseDialog->contents());
+  // new WText(" x' max = " + mixedTools::doubleToString(beam->getSigmaTransportij(2,2)) + " mrad ", ellipseDialog->contents());
+  // new WBreak(ellipseDialog->contents());
+  // new WText(" corr. 12 = " + mixedTools::doubleToString(beam->getSigmaTransportij(2,1))  + " sqrt(cm.rad) ", ellipseDialog->contents());
   
   vector<double> xcor;
   vector<double> ycor;
-  beam->donneesDessinEllipse(xcor,ycor,indexAbs, indexOrd);
+  vector<string> legende;
+  beam->donneesDessinEllipse(xcor,ycor,legende, indexAbs, indexOrd);
+  unsigned k;
+  for (int k=0 ; k < legende.size(); k++) {
+  new WBreak(ellipseDialog->contents());
+    new WText(legende.at(k), ellipseDialog->contents());
+  }
   //  beam->donneesDessinEllipseXxp(xcor,ycor);
   //  scatterPlot1D(ellipseDialog->contents(),xcor,ycor);  
@@ -863 +870 @@
   new WText(" mean : "+ mixedTools::doubleToString(out[1])+" MeV",w);
   new WBreak(w);
-  new WText(" rms : "+ mixedTools::doubleToString(out[2])+" KeV",w);
+  new WText(" sigma rms : "+ mixedTools::doubleToString(out[2])+" KeV",w);
   
   WCartesianChart *chart = new WCartesianChart(w);
-  chart->setTitle("Histogram of particle energies");
+  chart->setTitle("Histogram of kinetic energy");
   chart->setModel(model);        // set the model
   chart->setXSeriesColumn(0);    // set the column that holds the categories
@@ -881 +888 @@
   axis.setLabelFormat("%.3f");
   axis.setGridLinesEnabled(true);
-  axis.setTitle(WString("legende x"));
+  axis.setTitle(WString(" dEcin/Ecin (%)"));
 
   chart->axis(Y1Axis).setGridLinesEnabled(true);