Fork me on GitHub

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

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

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