77#include < cfloat>
88#include < cstdlib>
99
10- /* * Constructor */
11- BunchTuneAnalysis::BunchTuneAnalysis (): CppPyWrapper(NULL )
12- {
13- betax = 0 ;
14- alphax = 0 ;
15- etax = 0 ;
16- etapx = 0 ;
17- betay = 0 ;
18- alphay = 0 ;
10+
11+ BunchTuneAnalysis::BunchTuneAnalysis (): CppPyWrapper(NULL ) {
12+ double matrix[6 ][6 ] = {
13+ {1.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 },
14+ {0.0 , 1.0 , 0.0 , 0.0 , 0.0 , 0.0 },
15+ {0.0 , 0.0 , 1.0 , 0.0 , 0.0 , 0.0 },
16+ {0.0 , 0.0 , 0.0 , 1.0 , 0.0 , 0.0 },
17+ {0.0 , 0.0 , 0.0 , 0.0 , 1.0 , 0.0 },
18+ {0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 1.0 }
19+ };
1920}
2021
21- /* * Destructor */
22- BunchTuneAnalysis::~BunchTuneAnalysis ()
23- {
22+
23+ BunchTuneAnalysis::~BunchTuneAnalysis () {}
24+
25+ void BunchTuneAnalysis::setNormMatrixElement (int i, int j, double value) {
26+ matrix[i][j] = value;
2427}
2528
26- void BunchTuneAnalysis::assignTwiss (double bx, double ax, double dx, double dpx, double by, double ay){
27- betax = bx;
28- alphax = ax;
29- etax = dx;
30- etapx = dpx;
31- betay = by;
32- alphay = ay;
29+ double BunchTuneAnalysis::getNormMatrixElement (int i, int j) {
30+ return matrix[i][j];
31+ }
32+
33+ void BunchTuneAnalysis::assignTwiss (double betax, double alphax, double etax, double etapx, double betay, double alphay) {
34+ // Set V_{-1} = I
35+ for (int i = 0 ; i < 6 ; i++) {
36+ for (int j = 0 ; j < 6 ; j++) {
37+ matrix[i][j] = 0.0 ;
38+ }
39+ }
40+ for (int i = 0 ; i < 6 ; i++) {
41+ matrix[i][i] = 1.0 ;
42+ }
43+
44+ // 2D normalization (x-x')
45+ matrix[0 ][0 ] = 1.0 / sqrt (betax);
46+ matrix[0 ][1 ] = 0.0 ;
47+ matrix[1 ][0 ] = alphax / sqrt (betax);
48+ matrix[1 ][1 ] = sqrt (betax);
49+
50+ // 2D normalization (y-y')
51+ matrix[2 ][2 ] = 1.0 / sqrt (betay);
52+ matrix[2 ][3 ] = 0.0 ;
53+ matrix[3 ][2 ] = alphay / sqrt (betay);
54+ matrix[3 ][3 ] = sqrt (betay);
55+
56+ // Dispersion (x-x')
57+ matrix[0 ][5 ] = -etax / sqrt (betax);
58+ matrix[1 ][5 ] = -etax * (alphax / sqrt (betax)) - etapx * sqrt (betax);
59+
60+ // Dispersion (y-y')
61+ double etay = 0.0 ;
62+ double etapy = 0.0 ;
63+ matrix[2 ][5 ] = -etay / sqrt (betay);
64+ matrix[3 ][5 ] = -etay * (alphay / sqrt (betay)) - etapy * sqrt (betay);
65+
3366}
3467
35- /* * Performs the Tune analysis of the bunch */
3668void BunchTuneAnalysis::analyzeBunch (Bunch* bunch){
3769
38- // initialization
3970 bunch->compress ();
4071 SyncPart* syncPart = bunch->getSyncPart ();
4172 double beta = syncPart->getBeta ();
@@ -53,144 +84,56 @@ void BunchTuneAnalysis::analyzeBunch(Bunch* bunch){
5384 bunch->addParticleAttributes (" ParticlePhaseAttributes" , tunemap);
5485 }
5586
56- if (bunch->hasParticleAttributes (" ParticlePhaseAttributes" )){
87+ if (bunch->hasParticleAttributes (" ParticlePhaseAttributes" )){
5788 for (int i=0 ; i < bunch->getSize (); i++)
5889 {
90+ // Extract phase space coordinates
5991 double x = part_coord_arr[i][0 ];
6092 double xp = part_coord_arr[i][1 ];
6193 double y = part_coord_arr[i][2 ];
6294 double yp = part_coord_arr[i][3 ];
63- double Etot = syncPart->getEnergy () + syncPart->getMass ();
64- double dpp = 1 /(beta*beta)*part_coord_arr[i][5 ]/Etot;
65-
66- double xval = (x - etax * dpp)/sqrt (betax);
67- double xpval = (xp - etapx * dpp) * sqrt (betax) + xval * alphax;
68- double yval = y / sqrt (betay);
69- double ypval = (yp + y * alphay/betay) * sqrt (betay);
70-
71- double angle = atan2 (xpval, xval);
72- if (angle < 0 .) angle += (2.0 *OrbitConst::PI);
73- double xPhase = angle;
74- double xPhaseOld = bunch->getParticleAttributes (" ParticlePhaseAttributes" )->attValue (i,0 );
75- double xTune = (xPhaseOld - xPhase) / (2.0 *OrbitConst::PI);
76- if (xTune < 0 .) xTune += 1 .;
77- bunch->getParticleAttributes (" ParticlePhaseAttributes" )->attValue (i, 0 ) = xPhase;
78- bunch->getParticleAttributes (" ParticlePhaseAttributes" )->attValue (i, 2 ) = xTune;
79-
80- angle = atan2 (ypval, yval);
81- if (angle < 0 .) angle += (2.0 *OrbitConst::PI);
82- double yPhase = angle;
83- double yPhaseOld = bunch->getParticleAttributes (" ParticlePhaseAttributes" )->attValue (i,1 );
84- double yTune = (yPhaseOld - yPhase) / (2.0 *OrbitConst::PI);
85- if (yTune < 0 .) yTune += 1 .;
86- bunch->getParticleAttributes (" ParticlePhaseAttributes" )->attValue (i, 1 ) = yPhase;
87- bunch->getParticleAttributes (" ParticlePhaseAttributes" )->attValue (i, 3 ) = yTune;
88-
89- double xcanonical = x - etax * dpp;
90- double ycanonical = y;
91- double xpfac = xp - etapx * dpp;
92- double ypfac = yp;
93- double pxcanonical = xpfac + xcanonical * (alphax/betax);
94- double pycanonical = ypfac + ycanonical * (alphay/betay);
95- double xAction = xcanonical * xcanonical / betax + pxcanonical * pxcanonical * betax;
96- double yAction = ycanonical * ycanonical / betay + pycanonical * pycanonical * betay;
97-
98- bunch->getParticleAttributes (" ParticlePhaseAttributes" )->attValue (i, 4 ) = xAction;
99- bunch->getParticleAttributes (" ParticlePhaseAttributes" )->attValue (i, 5 ) = yAction;
100- }
101- }
102-
103- }
104-
105-
106-
107- BunchTuneAnalysis4D::BunchTuneAnalysis4D (): CppPyWrapper(NULL )
108- {
109- double matrix[4 ][4 ] = {
110- {1.0 , 0.0 , 0.0 , 0.0 },
111- {0.0 , 1.0 , 0.0 , 0.0 },
112- {0.0 , 0.0 , 1.0 , 0.0 },
113- {0.0 , 0.0 , 0.0 , 1.0 }
114- };
115- }
116-
117- BunchTuneAnalysis4D::~BunchTuneAnalysis4D ()
118- {
119- }
120-
121- void BunchTuneAnalysis4D::setNormMatrixElement (int i, int j, double value) {
122- matrix[i][j] = value;
123- }
124-
125- double BunchTuneAnalysis4D::getNormMatrixElement (int i, int j) {
126- return matrix[i][j];
127- }
128-
129- void BunchTuneAnalysis4D::analyzeBunch (Bunch* bunch){
130-
131- bunch->compress ();
132- SyncPart* syncPart = bunch->getSyncPart ();
133- double beta = syncPart->getBeta ();
134- double ** part_coord_arr = bunch->coordArr ();
135-
136- if (!bunch->hasParticleAttributes (" ParticlePhaseAttributes" )){
137- cerr << " Adding particle phase information attribute.\n " ;
138- std::map<std::string, double > tunemap;
139- tunemap.insert (std::make_pair (" phase1" , 0 ));
140- tunemap.insert (std::make_pair (" phase2" , 0 ));
141- tunemap.insert (std::make_pair (" tune1" , 0 ));
142- tunemap.insert (std::make_pair (" tune2" , 0 ));
143- tunemap.insert (std::make_pair (" action1" , 0 ));
144- tunemap.insert (std::make_pair (" action2" , 0 ));
145- bunch->addParticleAttributes (" ParticlePhaseAttributes" , tunemap);
146- }
147-
148- if (bunch->hasParticleAttributes (" ParticlePhaseAttributes" )) {
149- for (int i=0 ; i < bunch->getSize (); i++) {
150- // Get phase space coordinates
151- double x = part_coord_arr[i][0 ];
152- double xp = part_coord_arr[i][1 ];
153- double y = part_coord_arr[i][2 ];
154- double yp = part_coord_arr[i][3 ];
95+ double z = part_coord_arr[i][4 ];
15596 double Etot = syncPart->getEnergy () + syncPart->getMass ();
15697 double dpp = 1.0 / (beta * beta) * part_coord_arr[i][5 ] / Etot;
15798
158- // Normalize coordinates
159- double xval = matrix[0 ][0 ] * x + matrix[0 ][1 ] * xp + matrix[0 ][2 ] * y + matrix[0 ][3 ] * yp;
160- double xpval = matrix[1 ][0 ] * x + matrix[1 ][1 ] * xp + matrix[1 ][2 ] * y + matrix[1 ][3 ] * yp;
161- double yval = matrix[2 ][0 ] * x + matrix[2 ][1 ] * xp + matrix[2 ][2 ] * y + matrix[2 ][3 ] * yp;
162- double ypval = matrix[3 ][0 ] * x + matrix[3 ][1 ] * xp + matrix[3 ][2 ] * y + matrix[3 ][3 ] * yp;
99+ // Normalize phase space coordinates
100+ double xval = matrix[0 ][0 ] * x + matrix[0 ][1 ] * xp + matrix[0 ][2 ] * y + matrix[0 ][3 ] * yp + matrix[ 0 ][ 4 ] * z + matrix[ 0 ][ 5 ] * dpp ;
101+ double xpval = matrix[1 ][0 ] * x + matrix[1 ][1 ] * xp + matrix[1 ][2 ] * y + matrix[1 ][3 ] * yp + matrix[ 1 ][ 4 ] * z + matrix[ 1 ][ 5 ] * dpp ;
102+ double yval = matrix[2 ][0 ] * x + matrix[2 ][1 ] * xp + matrix[2 ][2 ] * y + matrix[2 ][3 ] * yp + matrix[ 2 ][ 4 ] * z + matrix[ 2 ][ 5 ] * dpp ;
103+ double ypval = matrix[3 ][0 ] * x + matrix[3 ][1 ] * xp + matrix[3 ][2 ] * y + matrix[3 ][3 ] * yp + matrix[ 3 ][ 4 ] * z + matrix[ 3 ][ 5 ] * dpp ;
163104
164- // Compute phases (mode 1)
105+ // Compute phase advance in normalized x-x'
165106 double angle = atan2 (xpval, xval);
166107 if (angle < 0.0 ) {
167108 angle += (2.0 * OrbitConst::PI);
168109 }
169110 double xPhase = angle;
170111 double xPhaseOld = bunch->getParticleAttributes (" ParticlePhaseAttributes" )->attValue (i, 0 );
171- double xTune = (xPhaseOld - xPhase) / (2.0 * OrbitConst::PI);
172- if (xTune < 0.0 ) {
112+ double xTune = (xPhaseOld - xPhase) / (2.0 * OrbitConst::PI);
113+ if (xTune < 0.0 ) {
173114 xTune += 1.0 ;
174115 }
175116 bunch->getParticleAttributes (" ParticlePhaseAttributes" )->attValue (i, 0 ) = xPhase;
176117 bunch->getParticleAttributes (" ParticlePhaseAttributes" )->attValue (i, 2 ) = xTune;
177118
119+ // Compute phase advance in normalized y-y'
178120 angle = atan2 (ypval, yval);
179121 if (angle < 0.0 ) {
180122 angle += (2.0 * OrbitConst::PI);
181123 }
182124 double yPhase = angle;
183125 double yPhaseOld = bunch->getParticleAttributes (" ParticlePhaseAttributes" )->attValue (i, 1 );
184- double yTune = (yPhaseOld - yPhase) / (2.0 * OrbitConst::PI);
126+ double yTune = (yPhaseOld - yPhase) / (2.0 * OrbitConst::PI);
185127 if (yTune < 0.0 ) {
186128 yTune += 1.0 ;
187129 }
188130 bunch->getParticleAttributes (" ParticlePhaseAttributes" )->attValue (i, 1 ) = yPhase;
189131 bunch->getParticleAttributes (" ParticlePhaseAttributes" )->attValue (i, 3 ) = yTune;
190132
133+ // Compute actions
191134 double xAction = xval * xval + xpval * xpval;
192135 double yAction = yval * yval + ypval * ypval;
193-
136+
194137 bunch->getParticleAttributes (" ParticlePhaseAttributes" )->attValue (i, 4 ) = xAction;
195138 bunch->getParticleAttributes (" ParticlePhaseAttributes" )->attValue (i, 5 ) = yAction;
196139 }
0 commit comments