Fork me on GitHub

source: svn/trunk/Utilities/Hector/src/H_OpticalElement.cc@ 249

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

typos
File size: 11.7 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_OpticalElement.cc
13/// \brief Class aiming at describing any beam optical element.
14
15// c++ #includes
16#include <iostream>
17#include <string>
18
19// ROOT #includes
20#include "TPaveLabel.h"
21
22// local #includes
23#include "H_Parameters.h"
24#include "H_TransportMatrices.h"
25#include "H_Aperture.h"
26#include "H_OpticalElement.h"
27using namespace std;
28
29/// called by the constructors
30void H_OpticalElement::init(const string& nameE, const int typeE, const double s, const double k, const double l) {
31 // this is called by the constructors
32 // must be in public section !
33 // xpos and ypos are vectors with n point. They define the aperture shape of the optical element.
34 name = nameE;
35 fs = s;
36 fk = k;
37 xpos = 0;
38 ypos = 0;
39 txpos = 0;
40 typos = 0;
41 element_length = l;
42 type = typeE;
43 element_mat.ResizeTo(MDIM,MDIM);
44 element_mat = driftmat(l);
45
46 // do NOT use setAperture for the initialisation ! there are protections there
47
48 if(element_length<0) { if(VERBOSE) cout<<"<H_OpticalElement> ERROR : Interpenetration of elements !"<<endl; }
49 if(element_length==0) { if(VERBOSE) cout<<"<H_OpticalElement> WARNING : 0-length element ! (" << name << ") " << " at " << fs << endl; }
50 betax =0;
51 betay =0;
52}
53
54H_OpticalElement::H_OpticalElement() : element_aperture(new H_Aperture()) {
55 init("",DRIFT,0.,0.,0.1);
56}
57
58H_OpticalElement::H_OpticalElement(const string& nameE, const int typeE, const double s, const double k, const double l) : element_aperture(new H_Aperture()) {
59 init(nameE,typeE,s,k,l);
60}
61
62H_OpticalElement::H_OpticalElement(const string& nameE, const int typeE, const double s, const double k, const double l, H_Aperture* the_app) : element_aperture(the_app->clone()) {
63 init(nameE,typeE,s,k,l);
64}
65
66H_OpticalElement::H_OpticalElement(const int typeE, const double s, const double k, const double l, H_Aperture* the_app) : element_aperture(the_app->clone()) {
67 init("",typeE,s,k,l);
68}
69
70H_OpticalElement::H_OpticalElement(const int typeE, const double s, const double k, const double l) : element_aperture(new H_Aperture()) {
71 init("",typeE,s,k,l);
72}
73
74H_OpticalElement::H_OpticalElement(const H_OpticalElement& el) {
75 fs = el.fs;
76 element_length = el.element_length;
77 fk = el.fk;
78 xpos = el.xpos;
79 ypos = el.ypos;
80 txpos = el.txpos;
81 typos = el.typos;
82 betax = el.betax;
83 betay = el.betay;
84 type = el.type;
85 name = el.name;
86 typestring = el.typestring;
87 element_mat.ResizeTo(MDIM,MDIM);
88 element_mat = el.element_mat;
89 element_aperture = el.element_aperture->clone();
90}
91
92H_OpticalElement& H_OpticalElement::operator=(const H_OpticalElement& el) {
93 if(this==&el) return *this;
94 fs = el.fs;
95 element_length = el.element_length;
96 fk = el.fk;
97 xpos = el.xpos;
98 ypos = el.ypos;
99 txpos = el.txpos;
100 typos = el.typos;
101 betax = el.betax;
102 betay = el.betay;
103 type = el.type;
104 name = el.name;
105 typestring = el.typestring;
106 element_mat.ResizeTo(MDIM,MDIM);
107 element_mat = el.element_mat;
108 element_aperture = el.element_aperture->clone();
109 return *this;
110}
111
112void H_OpticalElement::setAperture(const H_Aperture* ap) {
113 // do NOT use setAperture in your constructor, as element_aperture is not initialized
114 // this function do not take into account ap if ap=0
115 // do nothing if element_mat = ap
116 if (!ap) {cout << "<H_OpticalElement> Trying to set an empty pointer for the aperture ! Nothing done.\n"; return;}
117 if (element_aperture != ap) {
118 delete element_aperture;
119 element_aperture = ap->clone();
120 }
121 return;
122}
123/*
124void H_OpticalElement::setAperture(H_Aperture* ap) {
125 // do NOT use setAperture in your constructor, as element_aperture is not initialized
126 // this function do not take into account ap if ap=0
127 // do nothing if element_mat = ap
128 if (!ap) {cout << "<H_OpticalElement> Trying to set an empty pointer for the aperture ! Nothing done.\n"; return;}
129 if (element_aperture != ap) {
130 delete element_aperture;
131 element_aperture = ap; //->clone();
132 }
133 return;
134}
135*/
136
137
138
139void H_OpticalElement::printProperties() const {
140 cout << typestring;
141 cout << name;
142 cout <<"\t at s = " << fs;
143 cout <<"\t length = "<< element_length;
144 if(fk!=0) cout <<"\t strength = " << fk;
145 if(element_aperture->getType()!=NONE) {
146 cout <<"\t aperture type = " << element_aperture->getTypeString();
147 element_aperture->printProperties();
148 }
149
150 cout<<endl;
151 if(element_length<0) { if(VERBOSE) cout<<"<H_OpticalElement> ERROR : Interpenetration of elements !"<<endl; }
152 if(element_length==0) { if(VERBOSE) cout<<"<H_OpticalElement> WARNING : 0-length "<< typestring << " !" << endl; }
153
154 return;
155}
156
157void H_OpticalElement::showMatrix() const {
158 printMatrix(element_mat);
159 return;
160}
161
162void H_OpticalElement::drawAperture() const {
163 element_aperture->draw();
164 return;
165}
166
167TMatrix H_OpticalElement::getMatrix() {
168 setMatrix(0,MP,1);
169 return element_mat;
170}
171
172TMatrix H_OpticalElement::getMatrix(const float eloss, const float p_mass, const float p_charge) {
173 setMatrix(eloss,p_mass,p_charge);
174 return element_mat;
175}
176
177void H_OpticalElement::draw(const float meight, const float height) const{
178 /// @param meight is the minimal extend of the graph
179 /// @param height is the maximal extend of the graph
180 float x1 = getS();
181 float x2 = getS() + getLength();
182 float y1 = meight;
183 float y2 = height;
184 TPaveLabel* cur_box = new TPaveLabel(x1,y1,x2,y2,"");
185 cur_box->SetBorderSize(1);
186 cur_box->SetFillStyle(1001);
187 cur_box->SetFillColor((int)getType());
188 cur_box->Draw();
189}
190
191TVectorD H_OpticalElement::getHitPosition(TVectorD init_pos, double energy_loss, double mp, double qp) {
192 if(!element_length) {
193 // cout<<"O-length element ("<<getName()<<"), should not appear here !"<<endl;
194 return init_pos;
195 }
196 // some declarations
197 bool inside = false;
198 double vec[MDIM] = {init_pos[INDEX_X]/URAD,tan(init_pos[INDEX_TX]/URAD),init_pos[INDEX_Y]/URAD,tan(init_pos[INDEX_TY]/URAD),-energy_loss,1};
199 TMatrixD mat_init(1,MDIM,vec);
200 TMatrixD mat_min(1,MDIM,vec);
201 TMatrixD mat_max(1,MDIM,vec);
202 TMatrixD mat_stop(1,MDIM,vec);
203 H_OpticalElement* temp_el = clone();
204 // initialasing boundaries
205 double min_pos = 0;
206 double max_pos = element_length/2.;
207 double max_old = element_length/2.;
208 // number of iterations
209 // (idea : fix precision instead of number of iterations + add security)
210 // (idea : interpolate between max and min and give the error)
211 const int N = 10;
212 // starting search loop
213 for(int i = 0; i < N; i++) {
214 // fixing position to be investigated
215 temp_el->setLength(max_pos);
216 // initialising the vector at the initial vector + possible shift/tilt
217 mat_max[0][0] = mat_init[0][0] - temp_el->getX();
218 mat_max[0][1] = mat_init[0][1] - tan(temp_el->getTX());
219 mat_max[0][2] = mat_init[0][2] - temp_el->getY();
220 mat_max[0][3] = mat_init[0][3] - tan(temp_el->getTY());
221 // propagating
222 mat_max *= temp_el->getMatrix(energy_loss,mp,qp);
223 // compensating for the previously stated shifts/tilts
224 mat_max[0][0] = mat_max[0][0] + temp_el->getX();
225 mat_max[0][1] = mat_max[0][1] + tan(temp_el->getTX());
226 mat_max[0][2] = mat_max[0][2] + temp_el->getY();
227 mat_max[0][3] = mat_max[0][3] + tan(temp_el->getTY());
228 // fixing new boundaries
229 if(temp_el->isInside(mat_max.GetMatrixArray()[0]*URAD,mat_max.GetMatrixArray()[2]*URAD)) {
230 max_old = max_pos;
231 max_pos = max_pos + (max_pos - min_pos)/2.;
232 min_pos = max_old;
233 inside = true;
234 } else {
235 max_pos = min_pos + (max_pos - min_pos)/2.;
236 inside = false;
237 }
238 // end of loop
239 }
240 // if it passes at the last iteration, choosing the other range²
241 if(inside) min_pos = max_old;
242 // here the interpolation method : now we are sure that the intercept is between min_pos and max_pos
243 // getting vector at min_pos (for first boundary) :
244 bool precision_estimate = false;
245 if(precision_estimate) {
246 temp_el->setLength(min_pos);
247 mat_min[0][0] = mat_init[0][0] - temp_el->getX();
248 mat_min[0][1] = mat_init[0][1] - tan(temp_el->getTX());
249 mat_min[0][2] = mat_init[0][2] - temp_el->getY();
250 mat_min[0][3] = mat_init[0][3] - tan(temp_el->getTY());
251 mat_min *= temp_el->getMatrix(energy_loss,mp,qp);
252 mat_min[0][0] = mat_min[0][0] + temp_el->getX();
253 mat_min[0][1] = mat_min[0][1] + tan(temp_el->getTX());
254 mat_min[0][2] = mat_min[0][2] + temp_el->getY();
255 mat_min[0][3] = mat_min[0][3] + tan(temp_el->getTY());
256 mat_min[0][4] = min_pos + init_pos[4];
257 // getting vector at max_pos (for second boundary) :
258 temp_el->setLength(max_pos);
259 mat_max[0][0] = mat_init[0][0] - temp_el->getX();
260 mat_max[0][1] = mat_init[0][1] - tan(temp_el->getTX());
261 mat_max[0][2] = mat_init[0][2] - temp_el->getY();
262 mat_max[0][3] = mat_init[0][3] - tan(temp_el->getTY());
263 mat_max *= temp_el->getMatrix(energy_loss,mp,qp);
264 mat_max[0][0] = mat_max[0][0] + temp_el->getX();
265 mat_max[0][1] = mat_max[0][1] + tan(temp_el->getTX());
266 mat_max[0][2] = mat_max[0][2] + temp_el->getY();
267 mat_max[0][3] = mat_max[0][3] + tan(temp_el->getTY());
268 mat_max[0][4] = max_pos + init_pos[4];
269 }
270 // getting vector in the middle (for estimate) :
271 temp_el->setLength((max_pos+min_pos)/2.);
272 mat_stop[0][0] = mat_init[0][0] - temp_el->getX();
273 mat_stop[0][1] = mat_init[0][1] - tan(temp_el->getTX());
274 mat_stop[0][2] = mat_init[0][2] - temp_el->getY();
275 mat_stop[0][3] = mat_init[0][3] - tan(temp_el->getTY());
276 mat_stop *= temp_el->getMatrix(energy_loss,mp,qp);
277 mat_stop[0][0] = mat_stop[0][0] + temp_el->getX();
278 mat_stop[0][1] = mat_stop[0][1] + tan(temp_el->getTX());
279 mat_stop[0][2] = mat_stop[0][2] + temp_el->getY();
280 mat_stop[0][3] = mat_stop[0][3] + tan(temp_el->getTY());
281 mat_stop[0][4] = (max_pos+min_pos)/2. + init_pos[4];
282
283 double xys[LENGTH_VEC];
284 xys[INDEX_X]= mat_stop[0][0]*URAD;
285 xys[INDEX_TX]= atan(mat_stop[0][1])*URAD;
286 xys[INDEX_Y]= mat_stop[0][2]*URAD;
287 xys[INDEX_TY]= atan(mat_stop[0][3])*URAD;
288 xys[INDEX_S]= mat_stop[0][4] ;
289 TVectorD temp_vec(LENGTH_VEC,xys);
290
291 if(precision_estimate) {
292 cout<<"--- Results and precision estimates ---"<<endl;
293 cout<<"\t Stopping element : "<<getName()<<endl;
294 cout<<"\t hit point s : "<<mat_stop[0][4]<<" m +- "<<(mat_max[0][4]-mat_stop[0][4])*1000.<<" mm"<<endl;
295 cout<<"\t hit point x : "<<mat_stop[0][0]*URAD;
296 cout<<" + "<<fabs(((mat_min[0][0]<mat_max[0][0])?(mat_max[0][0]-mat_stop[0][0])*URAD:(mat_min[0][0]-mat_stop[0][0])*URAD));
297 cout<<" - "<<fabs(((mat_min[0][0]<mat_max[0][0])?(mat_stop[0][0]-mat_min[0][0])*URAD:(mat_stop[0][0]-mat_max[0][0])*URAD))<<" µm"<<endl;
298 cout<<"\t hit point y : "<<mat_stop[0][2]*URAD;
299 cout<<" + "<<fabs(((mat_min[0][0]<mat_max[0][2])?(mat_max[0][2]-mat_stop[0][2])*URAD:(mat_min[0][2]-mat_stop[0][2])*URAD));
300 cout<<" - "<<fabs(((mat_min[0][0]<mat_max[0][2])?(mat_stop[0][2]-mat_min[0][2])*URAD:(mat_stop[0][2]-mat_max[0][2])*URAD))<<" µm"<<endl;
301 cout<<"\t hit point tx : "<<mat_stop[0][1]*URAD;
302 cout<<" + "<<fabs(((mat_min[0][0]<mat_max[0][1])?(mat_max[0][1]-mat_stop[0][1])*URAD:(mat_min[0][1]-mat_stop[0][1])*URAD));
303 cout<<" - "<<fabs(((mat_min[0][0]<mat_max[0][1])?(mat_stop[0][1]-mat_min[0][1])*URAD:(mat_stop[0][1]-mat_max[0][1])*URAD))<<" µrad"<<endl;
304 cout<<"\t hit point ty : "<<mat_stop[0][3]*URAD;
305 cout<<" + "<<fabs(((mat_min[0][0]<mat_max[0][3])?(mat_max[0][3]-mat_stop[0][3])*URAD:(mat_min[0][3]-mat_stop[0][3])*URAD));
306 cout<<" - "<<fabs(((mat_min[0][0]<mat_max[0][3])?(mat_stop[0][3]-mat_min[0][3])*URAD:(mat_stop[0][3]-mat_max[0][3])*URAD))<<" µrad"<<endl;
307
308 }
309
310 delete temp_el;
311 return temp_vec;
312}
Note: See TracBrowser for help on using the repository browser.