src/Distribution/ClosedOrbitFinder.h

@@ -304,7 +304,7 @@ ClosedOrbitFinder<Value_type,
     
     // read in magnetic fieldmap
     bField_m.read(type, scaleFactor);
-
+    
     // compute closed orbit
     converged_m = findOrbit(accuracy, maxit);
 
@@ -408,6 +408,8 @@ inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_typ
      * The momentum in \beta * \gamma is obtained by dividing by "a"
      */
     value_type factor =  1.0 / acon_m(wo_m);
+    
+    
     std::for_each(pr.begin(), pr.end(), [factor](value_type p) { return p * factor; });
     
     return pr;
@@ -463,8 +465,8 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
     pr_m.resize(N_m+1);
 
     // store acon and bcon locally
-    value_type acon = acon_m(wo_m);               // [acon] = m
-    value_type invbcon = 1.0 / bcon_m(E0_m, wo_m);        // [bcon] = MeV*s/(C*m^2) = 10^6 T = 10^7 kG (kilo Gauss)
+    value_type acon = acon_m(wo_m);                     // [acon] = m
+    value_type invbcon = 1.0 / bcon_m(E0_m, wo_m);      // [bcon] = MeV*s/(C*m^2) = 10^6 T = 10^7 kG (kilo Gauss)
 
     // helper constants
     value_type p2;                                      // p^2 = p*p
@@ -481,7 +483,7 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
         r_m[idx] = y[0];
         pr_m[idx] = y[1];
 
-        // count number of crossings (excluding starting point --> idx>0)
+	// count number of crossings (excluding starting point --> idx>0)
         nxcross_m += (idx > 0) * (y[4] * xold < 0);
         xold = y[4];
         ++idx;
@@ -493,6 +495,7 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
                                                                                               const double theta)
     {
         pr2 = y[1] * y[1];
+		
         if (p2 < pr2)
             throw OpalException("ClosedOrbitFinder::findOrbit()", "p_{r}^2 > p^{2} (defined in Gordon paper) --> Square root of negative number.");
 
@@ -510,7 +513,7 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
         dydt[0] = y[0] * y[1] * invptheta;
         // Gordon, formula (5b)
         dydt[1] = ptheta - y[0] * bint;
-        // Gordon, formulas (9a) and (9b)
+	// Gordon, formulas (9a) and (9b)
         for (size_type i = 2; i < 5; i += 2) {
             dydt[i] = (y[1] * y[i] + y[0] * p2 * y[i+1] * invptheta * invptheta) * invptheta;
             dydt[i+1] = - y[1] * y[i+1] * invptheta - (bint + y[0] * brint) * y[i];

src/Distribution/Distribution.cpp

@@ -1306,9 +1306,12 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
                             "'CYCLOTRON' definition.");
     }
     
-    int Nint = 1000;
+    int Nint = 1440;
     double scaleFactor = 1.0;
     bool writeMap = true;
+    
+    double gamma = E_m / massIneV + 1.0;
+    double beta = std::sqrt(1.0 - 1.0 / (gamma * gamma));
 
     typedef SigmaGenerator<double, unsigned int> sGenerator_t;
     sGenerator_t *siggen = new sGenerator_t(I_m,
@@ -1334,7 +1337,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
                       CyclotronElement->getPHIinit(),
                       Attributes::getReal(itsAttr[Attrib::Distribution::RGUESS]),
                       Attributes::getString(itsAttr[Attrib::Distribution::FMTYPE]),
-                      false))  {
+                      false, false))  {
 
         std::array<double,3> Emit = siggen->getEmittances();
 
@@ -1359,37 +1362,35 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
           Units of the Sigma Matrix:
 
           spatial: mm
-          moment:  rad
+          moment:  mrad
 
         */
 
         if (Options::cloTuneOnly)
             throw OpalException("Do only CLO and tune calculation","... ");
 
-
         auto sigma = siggen->getSigma();
-        // change units from mm to m
-        for (unsigned int i = 0; i < 3; ++ i) {
-            for (unsigned int j = 0; j < 6; ++ j) {
-                sigma(2 * i, j) *= 1e-3;
-                sigma(j, 2 * i) *= 1e-3;
-            }
-        }
-
+        // convertion from mm mrad to m rad 
+        for (unsigned int i = 0; i < 6; ++ i)
+            for (unsigned int j = 0; j < 6; ++ j) sigma(i, j) *= 1e-6;
+                
         for (unsigned int i = 0; i < 3; ++ i) {
             if ( sigma(2 * i, 2 * i) < 0 || sigma(2 * i + 1, 2 * i + 1) < 0 )
                 throw OpalException("Distribution::CreateMatchedGaussDistribution()",
-                                    "Negative value on the diagonal of the sigma matrix.");
-
-            sigmaR_m[i] = std::sqrt(sigma(2 * i, 2 * i));
-            sigmaP_m[i] = std::sqrt(sigma(2 * i + 1, 2 * i + 1));
-        }
-
-        if (inputMoUnits_m == InputMomentumUnitsT::EV) {
-            for (unsigned int i = 0; i < 3; ++ i) {
-                sigmaP_m[i] = converteVToBetaGamma(sigmaP_m[i], massIneV);
-            }
+                                    "Negative value on the diagonal of the sigma matrix.");   
         }
+        
+        sigmaR_m[0] = std::sqrt(sigma(0, 0));
+        sigmaP_m[0] = std::sqrt(sigma(1, 1))*beta*gamma;
+        sigmaR_m[2] = std::sqrt(sigma(2, 2));
+        sigmaP_m[2] = std::sqrt(sigma(3, 3))*beta*gamma;
+        sigmaR_m[1] = std::sqrt(sigma(4, 4));
+        // delta = |p - p_0|/|p_0| = |p - p_0|/beta*gamma --> l'^2 = (delta*beta*gamma)^2 - (x'^2 + z'^2) 
+        //double lprime2 = sigma(5, 5)*(beta*beta*gamma*gamma) - sigmaP_m[0]*sigmaP_m[0] - sigmaP_m[2]*sigmaP_m[2];
+        //sigmaP_m[1] = std::sqrt(std::abs(lprime2));
+        double pl2 = (std::sqrt(sigma(5,5)) + 1)*(std::sqrt(sigma(5,5)) + 1)*beta*gamma*beta*gamma - sigmaP_m[0]*sigmaP_m[0] - sigmaP_m[2]*sigmaP_m[2];
+        double pl = std::sqrt(pl2);
+        sigmaP_m[1] = gamma*(pl - beta*gamma);
 
         correlationMatrix_m(1, 0) = sigma(0, 1) / (sqrt(sigma(0, 0) * sigma(1, 1)));
         correlationMatrix_m(3, 2) = sigma(2, 3) / (sqrt(sigma(2, 2) * sigma(3, 3)));
@@ -2475,8 +2476,35 @@ void Distribution::generateGaussZ(size_t numberOfParticles) {
     }
 #endif
 
+    //Sets the GSL error handler off, exception will be handled internaly with a renormalization method
+    gsl_set_error_handler_off();
     int errcode = gsl_linalg_cholesky_decomp(corMat);
-
+    double rn = 1e-12;
+
+    while (errcode == GSL_EDOM) {
+        
+        // Resets the correlation matrix
+        for (unsigned int i = 0; i < 6; ++ i) {
+            gsl_matrix_set(corMat, i, i, correlationMatrix_m(i, i));
+            for (unsigned int j = 0; j < i; ++ j) {
+                gsl_matrix_set(corMat, i, j, correlationMatrix_m(i, j));
+                gsl_matrix_set(corMat, j, i, correlationMatrix_m(i, j));
+            }
+        }
+        // Applying a renormalization method corMat = corMat + rn*Unitymatrix
+        // This is the renormalization
+        for(int i = 0; i < 6; i++){
+            double corMati = gsl_matrix_get(corMat, i , i);
+            gsl_matrix_set(corMat, i, i, corMati + rn);
+        }
+        //Just to be sure that the renormalization worked
+        errcode = gsl_linalg_cholesky_decomp(corMat);
+        if(errcode != GSL_EDOM) *gmsg << "* WARNING: Correlation matrix had to be renormalized"<<endl;
+        else rn *= 10;
+    }
+    //Sets again the standard GSL error handler on
+    gsl_set_error_handler(NULL);
+    //Just to be sure
     if (errcode == GSL_EDOM) {
         throw OpalException("Distribution::GenerateGaussZ",
                             "gsl_linalg_cholesky_decomp failed");
@@ -3464,11 +3492,11 @@ void Distribution::setAttributes() {
     // itsAttr[Attrib::Distribution::E2]
     //     = Attributes::makeReal("E2", "If E2<Eb, we compute the tunes from the beams energy Eb to E2 with dE=0.25 MeV ", 0.0);
     itsAttr[Attrib::Distribution::RESIDUUM]
-        = Attributes::makeReal("RESIDUUM", "Residuum for the closed orbit finder and sigma matrix generator ", 1e-8);
+        = Attributes::makeReal("RESIDUUM", "Residuum for the closed orbit finder and sigma matrix generator ", 1e-9);
     itsAttr[Attrib::Distribution::MAXSTEPSCO]
         = Attributes::makeReal("MAXSTEPSCO", "Maximum steps used to find closed orbit ", 100);
     itsAttr[Attrib::Distribution::MAXSTEPSSI]
-        = Attributes::makeReal("MAXSTEPSSI", "Maximum steps used to find matched distribution ",2000);
+        = Attributes::makeReal("MAXSTEPSSI", "Maximum steps used to find matched distribution ",1000);
     itsAttr[Attrib::Distribution::ORDERMAPS]
         = Attributes::makeReal("ORDERMAPS", "Order used in the field expansion ", 7);
     itsAttr[Attrib::Distribution::MAGSYM]
@@ -4735,4 +4763,4 @@ void Distribution::adjustPhaseSpace(double massIneV) {
 // mode:c++
 // c-basic-offset: 4
 // indent-tabs-mode:nil
-// End:
\ No newline at end of file
+// End:

src/Distribution/MapGenerator.h

@@ -147,7 +147,7 @@ MapGenerator<Value_type,
     
     for (size_type i = 0; i < 6; ++i) {
         for (size_type j = 0; j < 6; ++j) {
-            map(i,j) = matrix[i][j];
+	    map(i,j) = matrix[i][j];
         }
     }