Commit e79828a5 authored by cortes_c's avatar cortes_c
Browse files

bug fixed in Distribution.cpp

parent 5b3a63e3
...@@ -408,6 +408,8 @@ inline typename ClosedOrbitFinder<Value_type, Size_type, Stepper>::container_typ ...@@ -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" * The momentum in \beta * \gamma is obtained by dividing by "a"
*/ */
value_type factor = 1.0 / acon_m(wo_m); value_type factor = 1.0 / acon_m(wo_m);
std::for_each(pr.begin(), pr.end(), [factor](value_type p) { return p * factor; }); std::for_each(pr.begin(), pr.end(), [factor](value_type p) { return p * factor; });
return pr; return pr;
...@@ -463,8 +465,8 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc ...@@ -463,8 +465,8 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
pr_m.resize(N_m+1); pr_m.resize(N_m+1);
// store acon and bcon locally // store acon and bcon locally
value_type acon = acon_m(wo_m); // [acon] = m 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 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 // helper constants
value_type p2; // p^2 = p*p value_type p2; // p^2 = p*p
......
...@@ -1300,18 +1300,18 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub ...@@ -1300,18 +1300,18 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
const double fmLowE = CyclotronElement->getFMLowE(); const double fmLowE = CyclotronElement->getFMLowE();
const double fmHighE = CyclotronElement->getFMHighE(); const double fmHighE = CyclotronElement->getFMHighE();
double lE,hE; if ( fmLowE < 0 || fmHighE < 0 ) {
lE = fmLowE; throw OpalException("Distribution::CreateMatchedGaussDistribution()",
hE = fmHighE; "Missing attributes 'FMLOWE' and/or 'FMHIHGE' in "
"'CYCLOTRON' definition.");
if ((lE<0) || (hE<0)) {
lE = E_m*1E-6;
hE = E_m*1E-6;
} }
int Nint = 1440; int Nint = 1000;
double scaleFactor = 1.0; double scaleFactor = 1.0;
bool writeMap = true; 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; typedef SigmaGenerator<double, unsigned int> sGenerator_t;
sGenerator_t *siggen = new sGenerator_t(I_m, sGenerator_t *siggen = new sGenerator_t(I_m,
...@@ -1322,8 +1322,8 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub ...@@ -1322,8 +1322,8 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
E_m*1E-6, E_m*1E-6,
CyclotronElement->getCyclHarm(), CyclotronElement->getCyclHarm(),
massIneV*1E-6, massIneV*1E-6,
lE, fmLowE,
hE, fmHighE,
(int)Attributes::getReal(itsAttr[Attrib::Distribution::MAGSYM]), (int)Attributes::getReal(itsAttr[Attrib::Distribution::MAGSYM]),
Nint, Nint,
Attributes::getString(itsAttr[Attrib::Distribution::FMAPFN]), Attributes::getString(itsAttr[Attrib::Distribution::FMAPFN]),
...@@ -1366,11 +1366,31 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub ...@@ -1366,11 +1366,31 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
*/ */
//if (Options::cloTuneOnly) if (Options::cloTuneOnly)
// throw OpalException("Do only CLO and tune calculation","... "); throw OpalException("Do only CLO and tune calculation","... ");
auto sigma = siggen->getSigma(); auto sigma = siggen->getSigma();
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[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_p0| = |p - p_0|/beta*gamma --> l'^2 = (delta*beta*gamma)^2 - (x'^2 + z'^2)
double lprime = sigma(5, 5)*(beta*gamma)*(beta*gamma) - (sigmaP_m[0]*sigmaP_m[0] + sigmaP_m[1]*sigmaP_m[1]);
sigmaP_m[1] = std::sqrt(std::abs(lprime));
/*
// change units from mm to m // change units from mm to m
for (unsigned int i = 0; i < 3; ++ i) { for (unsigned int i = 0; i < 3; ++ i) {
for (unsigned int j = 0; j < 6; ++ j) { for (unsigned int j = 0; j < 6; ++ j) {
...@@ -1393,7 +1413,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub ...@@ -1393,7 +1413,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
sigmaP_m[i] = converteVToBetaGamma(sigmaP_m[i], massIneV); sigmaP_m[i] = converteVToBetaGamma(sigmaP_m[i], massIneV);
} }
} }
*/
correlationMatrix_m(1, 0) = sigma(0, 1) / (sqrt(sigma(0, 0) * sigma(1, 1))); 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))); correlationMatrix_m(3, 2) = sigma(2, 3) / (sqrt(sigma(2, 2) * sigma(3, 3)));
correlationMatrix_m(5, 4) = sigma(4, 5) / (sqrt(sigma(4, 4) * sigma(5, 5))); correlationMatrix_m(5, 4) = sigma(4, 5) / (sqrt(sigma(4, 4) * sigma(5, 5)));
...@@ -2477,11 +2497,11 @@ void Distribution::generateGaussZ(size_t numberOfParticles) { ...@@ -2477,11 +2497,11 @@ void Distribution::generateGaussZ(size_t numberOfParticles) {
*gmsg << " \\\\" << endl; *gmsg << " \\\\" << endl;
} }
#endif #endif
int errcode = gsl_linalg_cholesky_decomp(corMat);
//Sets the GSL error handler off, exception will be handled internaly with a renormalization method //Sets the GSL error handler off, exception will be handled internaly with a renormalization method
gsl_set_error_handler_off(); gsl_set_error_handler_off();
double rn = 1e-12;
int errcode = gsl_linalg_cholesky_decomp(corMat);
double rn = 1e-9;
while (errcode == GSL_EDOM) { while (errcode == GSL_EDOM) {
...@@ -2501,12 +2521,23 @@ void Distribution::generateGaussZ(size_t numberOfParticles) { ...@@ -2501,12 +2521,23 @@ void Distribution::generateGaussZ(size_t numberOfParticles) {
} }
//Just to be sure that the renormalization worked //Just to be sure that the renormalization worked
errcode = gsl_linalg_cholesky_decomp(corMat); errcode = gsl_linalg_cholesky_decomp(corMat);
if(errcode != GSL_EDOM) *gmsg << "WARNING: Correlation matrix had to be renormalized"<<endl; if(errcode != GSL_EDOM) *gmsg << "* WARNING: Correlation matrix had to be renormalized"<<endl;
else rn *= 10; else rn *= 10;
} }
//Sets again the standard GSL error handler on //Sets again the standard GSL error handler on
gsl_set_error_handler(NULL); gsl_set_error_handler(NULL);
if (errcode == GSL_EDOM) {
throw OpalException("Distribution::GenerateGaussZ",
"gsl_linalg_cholesky_decomp failed");
}
// so make the upper part zero.
for (int i = 0; i < 5; ++ i) {
for (int j = i+1; j < 6 ; ++ j) {
gsl_matrix_set (corMat, i, j, 0.0);
}
}
#define DISTDBG2 #define DISTDBG2
#ifdef DISTDBG2 #ifdef DISTDBG2
*gmsg << "* m after gsl_linalg_cholesky_decomp" << endl; *gmsg << "* m after gsl_linalg_cholesky_decomp" << endl;
......
...@@ -316,6 +316,11 @@ class SigmaGenerator ...@@ -316,6 +316,11 @@ class SigmaGenerator
* Mcycs is a vector of all cyclotron maps * Mcycs is a vector of all cyclotron maps
*/ */
void updateSigma(const std::vector<matrix_type>&, const std::vector<matrix_type>&); void updateSigma(const std::vector<matrix_type>&, const std::vector<matrix_type>&);
/// Writes the one turn matrix
/*
* @param One turn Matrix and iteration number
*/
void writeOneTurn(const matrix_type&, const value_type);
/// Writes the Matched Distributions /// Writes the Matched Distributions
/* /*
...@@ -498,7 +503,7 @@ template<typename Value_type, typename Size_type> ...@@ -498,7 +503,7 @@ template<typename Value_type, typename Size_type>
tunes = cof.getTunes(); tunes = cof.getTunes();
ravg = cof.getAverageRadius(); // average radius ravg = cof.getAverageRadius(); // average radius
container_type peo = cof.getMomentum(angle); container_type peo = cof.getMomentum(angle); //radial Momentum
// write to terminal // write to terminal
*gmsg << "* ----------------------------" << endl *gmsg << "* ----------------------------" << endl
...@@ -506,7 +511,7 @@ template<typename Value_type, typename Size_type> ...@@ -506,7 +511,7 @@ template<typename Value_type, typename Size_type>
<< "*" << endl << "*" << endl
<< "* average radius: " << ravg << " [m]" << endl << "* average radius: " << ravg << " [m]" << endl
<< "* initial radius: " << r_turn[0] << " [m]" << endl << "* initial radius: " << r_turn[0] << " [m]" << endl
<< "* initial momentum: " << peo[0] << " [Beta Gamma]" << endl << "* initial momentum: " << peo[0] << " [m] (Gordon Units)" << endl
<< "* frequency error: " << cof.getFrequencyError() << endl << "* frequency error: " << cof.getFrequencyError() << endl
<< "* horizontal tune: " << tunes.first << endl << "* horizontal tune: " << tunes.first << endl
<< "* vertical tune: " << tunes.second << endl << "* vertical tune: " << tunes.second << endl
...@@ -733,7 +738,7 @@ void SigmaGenerator<Value_type, Size_type>::initialize(value_type nuz, value_typ ...@@ -733,7 +738,7 @@ void SigmaGenerator<Value_type, Size_type>::initialize(value_type nuz, value_typ
// Remark: We multiply with 10^{6} (= mega) to convert emittances back. // Remark: We multiply with 10^{6} (= mega) to convert emittances back.
// 1 m^{2} = 10^{6} mm^{2} // 1 m^{2} = 10^{6} mm^{2}
matrix_type sigma = boost::numeric::ublas::zero_matrix<value_type>(6); matrix_type sigma = boost::numeric::ublas::zero_matrix<value_type>(6);
// formula (30), [sigma(0,0)] = m^2 rad * 10^{6} = mm^{2} rad = mm mrad // formula (30), [sigma(0,0)] = m^2 * 10^{-6} = mm^{2}
sigma(0,0) = invAB * (B * ex / Omega + A * ez / omega) * mega; sigma(0,0) = invAB * (B * ex / Omega + A * ez / omega) * mega;
// [sigma(0,5)] = [sigma(5,0)] = m rad * 10^{6} = mm mrad // 1000: m --> mm and 1000 to promille // [sigma(0,5)] = [sigma(5,0)] = m rad * 10^{6} = mm mrad // 1000: m --> mm and 1000 to promille
sigma(0,5) = sigma(5,0) = invAB * (ex / Omega + ez / omega) * mega; sigma(0,5) = sigma(5,0) = invAB * (ex / Omega + ez / omega) * mega;
...@@ -811,6 +816,8 @@ template<typename Value_type, typename Size_type> ...@@ -811,6 +816,8 @@ template<typename Value_type, typename Size_type>
mapgen.combine(Mscs,Mcycs); mapgen.combine(Mscs,Mcycs);
matrix_type Mturn = mapgen.getMap(); matrix_type Mturn = mapgen.getMap();
writeOneTurn(Mturn, niterations_m);
while (error_m > accuracy && !stop) { while (error_m > accuracy && !stop) {
/* /*
* Version 2.0 decouple was substituded by setEigenVectors * Version 2.0 decouple was substituded by setEigenVectors
...@@ -847,6 +854,8 @@ template<typename Value_type, typename Size_type> ...@@ -847,6 +854,8 @@ template<typename Value_type, typename Size_type>
mapgen.combine(Mscs,Mcycs); mapgen.combine(Mscs,Mcycs);
Mturn = mapgen.getMap(); Mturn = mapgen.getMap();
writeOneTurn(Mturn, niterations_m);
// check if number of iterations has maxit exceeded or if the eigenvector matrix could be built // check if number of iterations has maxit exceeded or if the eigenvector matrix could be built
stop = (niterations_m++ > maxit) || !isBuild; stop = (niterations_m++ > maxit) || !isBuild;
...@@ -1069,8 +1078,28 @@ typename SigmaGenerator<Value_type, Size_type>::value_type SigmaGenerator<Value_ ...@@ -1069,8 +1078,28 @@ typename SigmaGenerator<Value_type, Size_type>::value_type SigmaGenerator<Value_
return std::sqrt(std::inner_product(diff.data().begin(),diff.data().end(),diff.data().begin(),0.0)); return std::sqrt(std::inner_product(diff.data().begin(),diff.data().end(),diff.data().begin(),0.0));
} }
template<typename Value_type, typename Size_type> void SigmaGenerator<Value_type, Size_type>::writeOneTurn(const matrix_type& M, const value_type iteration){
std::stringstream ss;
ss << iteration;
std::string it = ss.str();
std::string energy = float2string(E_m);
std::ofstream oneTurn("data/oneTurn"+energy+"MeV.txt", std::ios::app);
oneTurn<<"Iteration: "+it<<std::endl;
for(size_type i = 0; i < 6; i++){
for(size_type j = 0; j < 6; j++){
oneTurn<<std::left<<std::setprecision(5)
<<std::setw(12)<<M(i,j);
}
oneTurn<<std::endl;
}
oneTurn<<std::endl;
oneTurn.close();
}
template<typename Value_type, typename Size_type> void SigmaGenerator<Value_type, Size_type>::writeMatched(const matrix_type& match, const size_type permutation,const std::vector<matrix_type>& Mcycs, const std::vector<matrix_type>& Mscs, const matrix_type& Mcyc){ template<typename Value_type, typename Size_type> void SigmaGenerator<Value_type, Size_type>::writeMatched(const matrix_type& match, const size_type permutation,const std::vector<matrix_type>& Mcycs, const std::vector<matrix_type>& Mscs, const matrix_type& Mcyc){
std::string energy = float2string(E_m); std::string energy = float2string(E_m);
std::ofstream writeSigmas("data/Matches"+energy+"MeV.txt", std::ios::app); std::ofstream writeSigmas("data/Matches"+energy+"MeV.txt", std::ios::app);
writeSigmas<<std::left<<std::setprecision(5) writeSigmas<<std::left<<std::setprecision(5)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment