Fork me on GitHub

source: git/external/Hector/H_TransportMatrices.cc@ cc5bcb7

ImprovedOutputFile Timing dual_readout llp
Last change on this file since cc5bcb7 was 3c40083, checked in by pavel <pavel@…>, 11 years ago

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