Fork me on GitHub

source: svn/trunk/Utilities/Hector/src/H_TransportMatrices.cc@ 552

Last change on this file since 552 was 281, checked in by Xavier Rouby, 16 years ago

new Hector version
File size: 10.6 KB
Line 
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * *
2 * *
3* --<--<-- A fast simulator --<--<-- *
4* / --<--<-- of particle --<--<-- *
5* ----HECTOR----< *
6* \ -->-->-- transport through -->-->-- *
7* -->-->-- generic beamlines -->-->-- *
8* *
9* JINST 2:P09005 (2007) *
10* X Rouby, J de Favereau, K Piotrzkowski (CP3) *
11* http://www.fynu.ucl.ac.be/hector.html *
12* *
13* Center for Cosmology, Particle Physics and Phenomenology *
14* Universite catholique de Louvain *
15* Louvain-la-Neuve, Belgium *
16 * *
17 * * * * * * * * * * * * * * * * * * * * * * * * * * * */
18
19/// \file H_TransportMatrices.cc
20/// \brief Includes the implementation of every transport matrix.
21
22// c++ #includes
23#include <iostream>
24
25// C #includes
26#include <cmath>
27
28// local #includes
29#include "H_Parameters.h"
30#include "H_TransportMatrices.h"
31using namespace std;
32
33bool relative_energy = 1;
34
35// caution : do not change particle mass, not implemented yet.
36
37extern double omega(const double k, const double l) {
38 // [l] = [m] and [k] = [1/m^2] for quadrupoles
39 // [omega] = [1]
40 return sqrt(fabs(k))*l;
41}
42
43extern double radius(const double k) {
44 // [k] = [1/m^2] for quadrupoles
45 // [k] = [1/m] for dipoles
46 // [radius(k)] = [m]
47 if(k==0 && VERBOSE) cout<<"<H_TransportMatrices> ERROR : Dipole has no effect : results will be corrupted"<<endl;
48 // this is protected by the "if(k==0) -> driftmat" in every matrix below (ex vquatmat)
49 return (k==0) ? 1 : 1/k;
50}
51
52extern void printMatrix(TMatrix TMat) {
53 char temp[20];
54 float * el = new float[MDIM*MDIM];
55 el = (TMat.GetMatrixArray());
56
57 cout << endl << "\t";
58 for(int i=0;i<MDIM*MDIM;i++) {
59 if (el[i]<0)
60 {sprintf(temp,"%.5e",el[i]);}
61 else if (el[i]>0)
62 {sprintf(temp," %.5e",el[i]);}
63 else {sprintf(temp," 0 ");}
64
65 cout << temp << " ";
66 if((i+1)%MDIM == 0) { cout << endl << "\t"; }
67 }
68 cout << endl;
69}
70
71extern TMatrix vquadmat(const float l, const float k, const float eloss = 0., const float p_mass=MP, const float p_charge=QP) {
72 // the length l is in [m]
73 // the strength k is in [1/mￜ] for quadrupoles
74 // eloss in [GeV]
75 // ke is the modified field with respect to the eloss
76 // k = e/p * dB/dx with p = mv (and m = MP)
77 // k -> ke = k * p/ (p - dp) <- chromacity
78 // ke -> ke * p_charge / QP <- if not a proton
79 // ke = 0 if charge = 0, whatever the mass
80
81 const double p0 = sqrt( (BE-MP)*(BE+MP) );
82 const double E = BE - eloss;
83 const double p = sqrt( (E-p_mass)*(E+p_mass) );
84 const float ke = (p_charge==0) ? 0 : k* p0/p *p_charge/QP;
85 if (ke==0) {
86 TMatrix drift(driftmat(l));
87 return drift;
88 }
89 // else... :
90 float om = omega(ke,l);
91 float * mat = new float[MDIM*MDIM];
92 float tmat[MDIM*MDIM] = {cosh(om),sqrt(ke)*sinh(om),0.,0., 0.,0.,
93 (1/sqrt(ke))*sinh(om),cosh(om),0.,0., 0.,0.,
94 0.,0.,cos(om),-sqrt(ke)*sin(om), 0.,0.,
95 0.,0.,(1/sqrt(ke))*sin(om),cos(om), 0.,0.,
96 0., 0., 0., 0., 1., 0.,
97 0., 0., 0., 0., 0., 1.
98 };
99 for (int i=0; i<MDIM*MDIM; i++) {mat[i] = tmat[i];}
100 TMatrix TMat(MDIM,MDIM,mat);
101 delete [] mat;
102 return TMat;
103}
104
105extern TMatrix hquadmat(const float l, const float k, const float eloss = 0., const float p_mass=MP, const float p_charge=QP) {
106 // the length l is in [m]
107 // the strength k is in [1/mￜ] for quadrupoles
108 // ke is the modified field with respect to the eloss
109 // k = e/p * dB/dx with p = mv (and m = MP)
110 // k -> ke = k * p/ (p- dp) <- chromacity
111 // ke -> ke *p_charge/QP <- if not a proton
112 const double p0 = sqrt( (BE-MP)*(BE+MP) );
113 const double E = BE - eloss;
114 const double p = sqrt( (E-p_mass)*(E+p_mass) );
115 const float ke = (p_charge==0) ? 0 : fabs(k* p0/p) *p_charge/QP;
116
117 if (ke==0) {
118 TMatrix drift(driftmat(l));
119 return drift;
120 }
121 float om = omega(ke,l);
122 float * mat = new float[MDIM*MDIM];
123 float tmat[MDIM*MDIM] = {cos(om),-sqrt(ke)*sin(om),0.,0., 0., 0.,
124 (1/sqrt(ke))*sin(om),cos(om),0.,0., 0., 0.,
125 0.,0.,cosh(om),sqrt(ke)*sinh(om), 0., 0.,
126 0.,0.,(1/sqrt(ke))*sinh(om),cosh(om), 0., 0.,
127 0., 0., 0., 0., 1., 0.,
128 0., 0., 0., 0., 0., 1.
129 };
130 for(int i=0;i<MDIM*MDIM;i++) { mat[i] = tmat[i]; }
131 TMatrix TMat(MDIM,MDIM,mat);
132 delete [] mat;
133 return TMat;
134}
135
136extern TMatrix rdipmat(const float l, const float k, const float eloss = 0., const float p_mass=MP, const float p_charge=QP) {
137 // the length l is in [m]
138 // the strength k is in [1/m] for dipoles
139 // ke is the modified field with respect to the eloss
140 // k = e/p * dB/dx with p = mv (and m = MP)
141 // k -> ke = k * p/ (p- dp) <- chromacity
142 // ke -> ke * q_mass/QP <- if not a proton
143
144 const double p0 = sqrt( (BE-MP)*(BE+MP) );
145 const double E = BE - eloss;
146 const double p = sqrt( (E-p_mass)*(E+p_mass) );
147 const float ke = (p_charge==0) ? 0 : k* p0/p *p_charge/QP;
148
149 if (ke==0) {
150 TMatrix drift(driftmat(l));
151 return drift;
152 }
153 float r = radius(ke);
154 float * mat = new float[MDIM*MDIM];
155 float * efmat = new float[MDIM*MDIM];
156 double simp = r*2*sin(l/(2*r))*sin(l/(2*r))/BE;
157 double psy = ke*l/2.;
158 float tefmat[MDIM*MDIM] = {1., tan(psy)*ke, 0., 0., 0., 0.,
159 0., 1., 0., 0., 0., 0.,
160 0., 0., 1., -tan(psy)*ke, 0., 0.,
161 0., 0., 0., 1., 0., 0.,
162 0., 0., 0., 0., 1., 0.,
163 0., 0., 0., 0., 0., 1. };
164
165 float tmat[MDIM*MDIM] = {cos(l/r),(-1/r)*sin(l/r),0.,0., 0., 0.,
166 r*sin(l/r),cos(l/r),0.,0., 0., 0.,
167 0.,0.,1.,0., 0., 0.,
168 0.,0.,l,1., 0., 0.,
169 simp, sin(l/r)/BE, 0., 0., 1., 0.,
170 0., 0., 0., 0., 0., 1. };
171 for(int i=0;i<MDIM*MDIM;i++) {
172 mat[i] = tmat[i];
173 efmat[i] = tefmat[i];
174 }
175 TMatrix TMat(MDIM,MDIM,mat);
176 TMatrix TEfmat(MDIM,MDIM,efmat);
177 if(relative_energy) {
178 TMat *= TEfmat;
179 TEfmat *= TMat;
180 }
181
182// if(VERBOSE) cout<<"\t WARNING : RDipoles not implemented and replaced by SDipoles" << endl;
183 delete [] mat;
184 delete [] efmat;
185 if(relative_energy) {
186 return TEfmat;
187 } else {
188 return TMat;
189 }
190}
191
192extern TMatrix sdipmat(const float l, const float k, const float eloss = 0., const float p_mass=MP, const float p_charge=QP) {
193 // the length l is in [m]
194 // the strength k is in [1/m] for dipoles
195 // ke is the modified field with respect to the eloss
196 // k = e/p * dB/dx with p = mv (and m = MP)
197 // k -> ke = k * p/ (p- dp) <- chromacity
198 // ke -> ke * q_mass/QP <- if not a proton
199
200 const double p0 = sqrt( (BE-MP)*(BE+MP) );
201 const double E = BE - eloss;
202 const double p = sqrt( (E-p_mass)*(E+p_mass) );
203 const float ke = (p_charge==0) ? 0 : k* p0/p *p_charge/QP;
204
205 if (ke==0) {
206 TMatrix drift(driftmat(l));
207 return drift;
208 }
209 extern bool relative_energy;
210 float r = radius(ke);
211 float * mat = new float[MDIM*MDIM];
212
213 float simp = 2*r*sin(l/(2*r))*sin(l/(2*r))/BE;
214 float tmat[MDIM*MDIM] = {cos(l/r),(-1/r)*sin(l/r),0.,0., 0., 0.,
215 r*sin(l/r),cos(l/r),0.,0., 0., 0.,
216 0.,0.,1.,0., 0., 0.,
217 0.,0.,l,1., 0., 0.,
218 simp, sin(l/r)/BE, 0., 0., 1., 0.,
219 0., 0., 0., 0., 0., 1.
220 };
221 if(!relative_energy) {
222 tmat[24] = 0;
223 tmat[25] = 0;
224 }
225
226 for(int i=0;i<MDIM*MDIM;i++) { mat[i] = tmat[i]; }
227 TMatrix TMat(MDIM,MDIM,mat);
228 delete [] mat;
229 return TMat;
230}
231
232extern TMatrix driftmat(const float l) {
233 // the length l is in [m]
234 float * mat = new float[MDIM*MDIM];
235 float tmat[MDIM*MDIM] = {1.,0.,0.,0.,0.,0.,
236 l ,1.,0.,0.,0.,0.,
237 0.,0.,1.,0.,0.,0.,
238 0.,0.,l ,1.,0.,0.,
239 0.,0.,0.,0.,1.,0.,
240 0., 0., 0., 0., 0., 1.
241 };
242 for(int i=0;i<MDIM*MDIM;i++) { mat[i] = tmat[i]; }
243 TMatrix TMat(MDIM,MDIM,mat);
244 delete [] mat;
245 return TMat;
246}
247
248
249extern TMatrix hkickmat(const float l, const float k, const float eloss =0., const float p_mass=MP, const float p_charge=QP) {
250 // the length l is in [m]
251 // the strength k is in [rad]
252 const double p0 = sqrt( (BE-MP)*(BE+MP) );
253 const double E = BE - eloss;
254 const double p = sqrt( (E-p_mass)*(E+p_mass) );
255 const float ke = (p_charge==0) ? 0 : -k* p0/p *p_charge/QP;
256
257 if (ke==0) {
258 TMatrix drift(driftmat(l));
259 return drift;
260 }
261 float * mat = new float[MDIM*MDIM];
262 float tmat[MDIM*MDIM] = {1.,0.,0.,0.,0.,0.,
263 l ,1.,0.,0.,0.,0.,
264 0.,0.,1.,0.,0.,0.,
265 0.,0.,l ,1.,0.,0.,
266 0.,0.,0.,0.,1.,0.,
267 l*tan(ke)/2.,ke, 0., 0., 0., 1.
268 };
269
270 for(int i=0;i<MDIM*MDIM;i++) { mat[i] = tmat[i]; }
271 TMatrix TMat(MDIM,MDIM,mat);
272 delete [] mat;
273 return TMat;
274}
275
276extern TMatrix vkickmat(const float l, const float k, const float eloss=0., const float p_mass=MP, const float p_charge=QP) {
277 // the length l is in [m]
278 // the strength k is in [rad]
279 const double p0 = sqrt( (BE-MP)*(BE+MP) );
280 const double E = BE - eloss;
281 const double p = sqrt( (E-p_mass)*(E+p_mass) );
282 const float ke = (p_charge==0) ? 0 : -k* p0/p *p_charge/QP;
283
284 if (ke==0) {
285 TMatrix drift(driftmat(l));
286 return drift;
287 }
288 float * mat = new float[MDIM*MDIM];
289 float tmat[MDIM*MDIM] = {1.,0.,0.,0.,0.,0.,
290 l ,1.,0.,0.,0.,0.,
291 0.,0.,1.,0.,0.,0.,
292 0.,0.,l ,1.,0.,0.,
293 0.,0.,0.,0.,1.,0.,
294 0.,0.,l*tan(ke)/2.,ke, 0., 1.
295 };
296
297 for(int i=0;i<MDIM*MDIM;i++) { mat[i] = tmat[i]; }
298 TMatrix TMat(MDIM,MDIM,mat);
299 delete [] mat;
300 return TMat;
301}
302
Note: See TracBrowser for help on using the repository browser.