Fork me on GitHub

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

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

first commit
File size: 11.1 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::printProperties() const {
125 cout << typestring;
126 cout << name;
127 cout <<"\t at s = " << fs;
128 cout <<"\t length = "<< element_length;
129 if(fk!=0) cout <<"\t strength = " << fk;
130 if(element_aperture->getType()!=NONE) {
131 cout <<"\t aperture type = " << element_aperture->getTypeString();
132 element_aperture->printProperties();
133 }
134
135 cout<<endl;
136 if(element_length<0) { if(VERBOSE) cout<<"<H_OpticalElement> ERROR : Interpenetration of elements !"<<endl; }
137 if(element_length==0) { if(VERBOSE) cout<<"<H_OpticalElement> WARNING : 0-length "<< typestring << " !" << endl; }
138
139 return;
140}
141
142void H_OpticalElement::showMatrix() const {
143 printMatrix(element_mat);
144 return;
145}
146
147void H_OpticalElement::drawAperture() const {
148 element_aperture->draw();
149 return;
150}
151
152TMatrix H_OpticalElement::getMatrix() {
153 setMatrix(0,MP,1);
154 return element_mat;
155}
156
157TMatrix H_OpticalElement::getMatrix(const float eloss, const float p_mass, const float p_charge) {
158 setMatrix(eloss,p_mass,p_charge);
159 return element_mat;
160}
161
162void H_OpticalElement::draw(const float meight, const float height) const{
163 /// @param meight is the minimal extend of the graph
164 /// @param height is the maximal extend of the graph
165 float x1 = getS();
166 float x2 = getS() + getLength();
167 float y1 = meight;
168 float y2 = height;
169 TPaveLabel* cur_box = new TPaveLabel(x1,y1,x2,y2,"");
170 cur_box->SetBorderSize(1);
171 cur_box->SetFillStyle(1001);
172 cur_box->SetFillColor((int)getType());
173 cur_box->Draw();
174}
175
176TVectorD H_OpticalElement::getHitPosition(TVectorD init_pos, double energy_loss, double mp, double qp) {
177 if(!element_length) {
178 // cout<<"O-length element ("<<getName()<<"), should not appear here !"<<endl;
179 return init_pos;
180 }
181 // some declarations
182 bool inside = false;
183 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};
184 TMatrixD mat_init(1,MDIM,vec);
185 TMatrixD mat_min(1,MDIM,vec);
186 TMatrixD mat_max(1,MDIM,vec);
187 TMatrixD mat_stop(1,MDIM,vec);
188 H_OpticalElement* temp_el = clone();
189 // initialasing boundaries
190 double min_pos = 0;
191 double max_pos = element_length/2.;
192 double max_old = element_length/2.;
193 // number of iterations
194 // (idea : fix precision instead of number of iterations + add security)
195 // (idea : interpolate between max and min and give the error)
196 const int N = 10;
197 // starting search loop
198 for(int i = 0; i < N; i++) {
199 // fixing position to be investigated
200 temp_el->setLength(max_pos);
201 // initialising the vector at the initial vector + possible shift/tilt
202 mat_max[0][0] = mat_init[0][0] - temp_el->getX();
203 mat_max[0][1] = mat_init[0][1] - tan(temp_el->getTX());
204 mat_max[0][2] = mat_init[0][2] - temp_el->getY();
205 mat_max[0][3] = mat_init[0][3] - tan(temp_el->getTY());
206 // propagating
207 mat_max *= temp_el->getMatrix(energy_loss,mp,qp);
208 // compensating for the previously stated shifts/tilts
209 mat_max[0][0] = mat_max[0][0] + temp_el->getX();
210 mat_max[0][1] = mat_max[0][1] + tan(temp_el->getTX());
211 mat_max[0][2] = mat_max[0][2] + temp_el->getY();
212 mat_max[0][3] = mat_max[0][3] + tan(temp_el->getTY());
213 // fixing new boundaries
214 if(temp_el->isInside(mat_max.GetMatrixArray()[0]*URAD,mat_max.GetMatrixArray()[2]*URAD)) {
215 max_old = max_pos;
216 max_pos = max_pos + (max_pos - min_pos)/2.;
217 min_pos = max_old;
218 inside = true;
219 } else {
220 max_pos = min_pos + (max_pos - min_pos)/2.;
221 inside = false;
222 }
223 // end of loop
224 }
225 // if it passes at the last iteration, choosing the other range²
226 if(inside) min_pos = max_old;
227 // here the interpolation method : now we are sure that the intercept is between min_pos and max_pos
228 // getting vector at min_pos (for first boundary) :
229 bool precision_estimate = false;
230 if(precision_estimate) {
231 temp_el->setLength(min_pos);
232 mat_min[0][0] = mat_init[0][0] - temp_el->getX();
233 mat_min[0][1] = mat_init[0][1] - tan(temp_el->getTX());
234 mat_min[0][2] = mat_init[0][2] - temp_el->getY();
235 mat_min[0][3] = mat_init[0][3] - tan(temp_el->getTY());
236 mat_min *= temp_el->getMatrix(energy_loss,mp,qp);
237 mat_min[0][0] = mat_min[0][0] + temp_el->getX();
238 mat_min[0][1] = mat_min[0][1] + tan(temp_el->getTX());
239 mat_min[0][2] = mat_min[0][2] + temp_el->getY();
240 mat_min[0][3] = mat_min[0][3] + tan(temp_el->getTY());
241 mat_min[0][4] = min_pos + init_pos[4];
242 // getting vector at max_pos (for second boundary) :
243 temp_el->setLength(max_pos);
244 mat_max[0][0] = mat_init[0][0] - temp_el->getX();
245 mat_max[0][1] = mat_init[0][1] - tan(temp_el->getTX());
246 mat_max[0][2] = mat_init[0][2] - temp_el->getY();
247 mat_max[0][3] = mat_init[0][3] - tan(temp_el->getTY());
248 mat_max *= temp_el->getMatrix(energy_loss,mp,qp);
249 mat_max[0][0] = mat_max[0][0] + temp_el->getX();
250 mat_max[0][1] = mat_max[0][1] + tan(temp_el->getTX());
251 mat_max[0][2] = mat_max[0][2] + temp_el->getY();
252 mat_max[0][3] = mat_max[0][3] + tan(temp_el->getTY());
253 mat_max[0][4] = max_pos + init_pos[4];
254 }
255 // getting vector in the middle (for estimate) :
256 temp_el->setLength((max_pos+min_pos)/2.);
257 mat_stop[0][0] = mat_init[0][0] - temp_el->getX();
258 mat_stop[0][1] = mat_init[0][1] - tan(temp_el->getTX());
259 mat_stop[0][2] = mat_init[0][2] - temp_el->getY();
260 mat_stop[0][3] = mat_init[0][3] - tan(temp_el->getTY());
261 mat_stop *= temp_el->getMatrix(energy_loss,mp,qp);
262 mat_stop[0][0] = mat_stop[0][0] + temp_el->getX();
263 mat_stop[0][1] = mat_stop[0][1] + tan(temp_el->getTX());
264 mat_stop[0][2] = mat_stop[0][2] + temp_el->getY();
265 mat_stop[0][3] = mat_stop[0][3] + tan(temp_el->getTY());
266 mat_stop[0][4] = (max_pos+min_pos)/2. + init_pos[4];
267
268 double xys[LENGTH_VEC];
269 xys[INDEX_X]= mat_stop[0][0]*URAD;
270 xys[INDEX_TX]= atan(mat_stop[0][1])*URAD;
271 xys[INDEX_Y]= mat_stop[0][2]*URAD;
272 xys[INDEX_TY]= atan(mat_stop[0][3])*URAD;
273 xys[INDEX_S]= mat_stop[0][4] ;
274 TVectorD temp_vec(LENGTH_VEC,xys);
275
276 if(precision_estimate) {
277 cout<<"--- Results and precision estimates ---"<<endl;
278 cout<<"\t Stopping element : "<<getName()<<endl;
279 cout<<"\t hit point s : "<<mat_stop[0][4]<<" m +- "<<(mat_max[0][4]-mat_stop[0][4])*1000.<<" mm"<<endl;
280 cout<<"\t hit point x : "<<mat_stop[0][0]*URAD;
281 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));
282 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;
283 cout<<"\t hit point y : "<<mat_stop[0][2]*URAD;
284 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));
285 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;
286 cout<<"\t hit point tx : "<<mat_stop[0][1]*URAD;
287 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));
288 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;
289 cout<<"\t hit point ty : "<<mat_stop[0][3]*URAD;
290 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));
291 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;
292
293 }
294
295 delete temp_el;
296 return temp_vec;
297}
Note: See TracBrowser for help on using the repository browser.