Fork me on GitHub

source: git/external/Hector/H_AbstractBeamLine.cc@ 5af205e

Last change on this file since 5af205e was 3c40083, checked in by pavel <pavel@…>, 11 years ago

switch to a more stable Hector version
  • Property mode set to 100644
File size: 20.5 KB
RevLine 
[3c40083]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
[5b822e5]12
13/// \file H_AbstractBeamLine.cc
14/// \brief Class describing ideal beamline.
15///
16/// Units : angles [rad], distances [m], energies [GeV], c=[1].
17
18// c++ #includes
19#include <iostream>
20#include <cmath>
21
22// ROOT #includes
[3c40083]23//#include "TPaveLabel.h"
24//#include "TLine.h"
25//#include "TGaxis.h"
26//#include "TLegend.h"
27//#include "TF1.h"
28//#include "TROOT.h"
[5b822e5]29
30// local #includes
31#include "H_Parameters.h"
32#include "H_TransportMatrices.h"
33#include "H_Drift.h"
34#include "H_AbstractBeamLine.h"
[3c40083]35#include "H_BeamParticle.h"
[5b822e5]36#include "H_RomanPot.h"
37using namespace std;
38
39void H_AbstractBeamLine::init(const float length) {
[3c40083]40 beam_mat = new TMatrix(MDIM,MDIM);
[5b822e5]41 beam_length = length;
42 H_Drift * drift0 = new H_Drift("Drift0",0.,length);
43 add(drift0);
44 return;
45}
46
[3c40083]47H_AbstractBeamLine::H_AbstractBeamLine(const H_AbstractBeamLine& beamline) {
48 elements = beamline.elements;
49 matrices = beamline.matrices;
50 beam_mat = new TMatrix(*(beamline.beam_mat));
[5b822e5]51 beam_length = beamline.beam_length;
52}
53
54H_AbstractBeamLine& H_AbstractBeamLine::operator=(const H_AbstractBeamLine& beamline) {
55 if(this== &beamline) return *this;
[3c40083]56 elements = beamline.elements;
[5b822e5]57 matrices = beamline.matrices;
[3c40083]58 beam_mat = new TMatrix(*(beamline.beam_mat));
[5b822e5]59 beam_length = beamline.beam_length;
60 return *this;
61}
62
63H_AbstractBeamLine::~H_AbstractBeamLine() {
64 vector<H_OpticalElement*>::iterator element_i;
65 for (element_i = elements.begin(); element_i<elements.end(); element_i++) {
66 delete (*element_i);
67 }
68 elements.clear();
69 matrices.clear();
[3c40083]70 delete beam_mat;
[5b822e5]71}
72
73void H_AbstractBeamLine::add(H_OpticalElement * newElement) {
74 /// @param newElement is added to the beamline
75// H_OpticalElement * el = new H_OpticalElement(*newElement);
76// H_OpticalElement * el = const_cast<H_OpticalElement*> newElement;
77// H_OpticalElement * el = newElement;
78 elements.push_back(newElement);
79 float a = newElement->getS()+newElement->getLength();
80 if (a > beam_length) {
81 beam_length = a;
[3c40083]82 if(VERBOSE) cout<<"WARNING : element ("<< newElement->getName()<<") too far away. The beam length has been extended to "<< beam_length << ". "<<endl;
83 }
84 calcSequence();
85 calcMatrix();
86}
87
88void H_AbstractBeamLine::add(H_OpticalElement & newElement) {
89 /// @param newElement is added to the beamline
90// H_OpticalElement * el = new H_OpticalElement(newElement);
91// elements.push_back(el);
92 elements.push_back(&newElement);
93 float a = newElement.getS()+newElement.getLength();
94 if (a > beam_length) {
95 beam_length = a;
96 if(VERBOSE) cout<<"WARNING : element ("<< newElement.getName()<<") too far away. The beam length has been extended to "<< beam_length << ". "<<endl;
[5b822e5]97 }
98 calcSequence();
99 calcMatrix();
100}
101
[3c40083]102const TMatrix * H_AbstractBeamLine::getBeamMatrix() const {
103 TMatrix * mat = new TMatrix(*beam_mat);
104 return mat;
[5b822e5]105}
106
[3c40083]107const TMatrix * H_AbstractBeamLine::getBeamMatrix(const float eloss,const float p_mass, const float p_charge) {
[5b822e5]108
109 vector<H_OpticalElement*>::iterator element_i;
110 TMatrix calc_mat(MDIM,MDIM);
111
112 // initialization
113 calc_mat.UnitMatrix();
114
115 // multiplies the matrix of each beam's element
116 // and add each product matrix to the list of matrices.
117 for(element_i = elements.begin(); element_i < elements.end(); element_i++) {
118 calc_mat *= (*element_i)->getMatrix(eloss,p_mass,p_charge);
119 }
[3c40083]120 const TMatrix* bmat = new TMatrix(calc_mat);
121 return bmat;
[5b822e5]122}
123
[3c40083]124const TMatrix * H_AbstractBeamLine::getPartialMatrix(const string elname, const float eloss, const float p_mass, const float p_charge) {
[5b822e5]125
126 vector<H_OpticalElement*>::iterator element_i;
127 TMatrix calc_mat(MDIM,MDIM);
128
129 calc_mat.UnitMatrix();
130
131 for(element_i = elements.begin(); element_i < elements.end(); element_i++) {
132 calc_mat *= (*element_i)->getMatrix(eloss,p_mass,p_charge);
133 if(elname==(*element_i)->getName()) {
[3c40083]134 const TMatrix* bmat = new TMatrix(calc_mat);
135 return bmat;
[5b822e5]136 }
137 }
[3c40083]138 cout<<"Element "<<elname<<" desn't exist. Returning full beam matrix"<<endl;
139 const TMatrix* bmat = new TMatrix(calc_mat);
140 return bmat;
[5b822e5]141}
142
[3c40083]143const TMatrix * H_AbstractBeamLine::getPartialMatrix(const unsigned int element_position) const {
[5b822e5]144 //const int N = (element_position<0)?0:(( (element_position)>elements.size()-1)?elements.size()-1:element_position);
145 const int N = (element_position>elements.size()-1)?elements.size()-1:element_position;
[3c40083]146 return &(*(matrices.begin()+N)); // //for optimization of the code :same as return &matrices[N];
[5b822e5]147}
148
[3c40083]149const TMatrix * H_AbstractBeamLine::getPartialMatrix(const H_OpticalElement * element) const{
[5b822e5]150 // returns the transport matrix to transport until the end of the specified element
151 // !!! 2 elements should never have the same name in "elements" !!!
152
153 vector<H_OpticalElement*>::const_iterator element_i;
154 vector<TMatrix>::const_iterator matrix_i;
[3c40083]155 TMatrix * calc_mat = new TMatrix(MDIM,MDIM);
[5b822e5]156
157 // parses the list of optical elements and find the searched one
158 for(element_i = elements.begin(),matrix_i = matrices.begin(); element_i < elements.end(); element_i++, matrix_i++) {
159 if(element->getName() == (*element_i)->getName()) {
160 // element has been found
[3c40083]161 calc_mat = const_cast<TMatrix*>( &(*matrix_i));
[5b822e5]162 }
163 }
[3c40083]164 return (const TMatrix*) calc_mat;
[5b822e5]165}
166
167H_OpticalElement * H_AbstractBeamLine::getElement(const unsigned int element_position) {
168 const unsigned int N = (element_position>elements.size())?elements.size():element_position;
169 return *(elements.begin()+N);//for optimization of the code :same as return &elements[N];
170}
171
[3c40083]172const H_OpticalElement * H_AbstractBeamLine::getElement(const unsigned int element_position) const {
[5b822e5]173 const unsigned int N = (element_position>elements.size())?elements.size():element_position;
174 return *(elements.begin()+N);//for optimization of the code :same as return &elements[N];
175}
176
177
[3c40083]178H_OpticalElement * H_AbstractBeamLine::getElement(const string el_name) {
[5b822e5]179 for(unsigned int i=0; i < elements.size(); i++) {
180 if( (*(elements.begin()+i))->getName() == el_name )
181 return *(elements.begin()+i);
182 } // if found -> return ; else : not found at all !
[3c40083]183 cout<<"Element "<<el_name<<" not found"<<endl;
[5b822e5]184 return *(elements.begin()+1);
185}
186
[3c40083]187const H_OpticalElement * H_AbstractBeamLine::getElement(const string el_name) const {
[5b822e5]188 for(unsigned int i=0; i < elements.size(); i++) {
189 if( (*(elements.begin()+i))->getName() == el_name)
190 return *(elements.begin()+i);
191 } // if found -> return ; else : not found at all !
[3c40083]192 cout<<"Element "<<el_name<<" not found"<<endl;
[5b822e5]193 return *(elements.begin()+1);
194}
195
[3c40083]196void H_AbstractBeamLine::printProperties() const {
197 vector<H_OpticalElement*>::const_iterator element_i;
198 cout << "Pointeurs des elements du faisceau" << endl;
199 for (element_i = elements.begin(); element_i < elements.end(); element_i++) {
200 cout << (int)(element_i-elements.begin()) << "\t" << (*element_i)->getName() << "\t" << (*element_i)->getS() << endl;
201 }
202 return;
[5b822e5]203}
204
205void H_AbstractBeamLine::showElements() const{
206 vector<H_OpticalElement*>::const_iterator element_i;
207 for(element_i = elements.begin(); element_i < elements.end(); element_i++) {
208 (*element_i)->printProperties();
209 }
210 cout << "Beam length = " << beam_length << endl;
211 cout << "Number of elements (including drifts) = " << getNumberOfElements() << endl;
212 return;
213}
214
215void H_AbstractBeamLine::showElements(const int type_el) const{
216 vector<H_OpticalElement*>::const_iterator element_i;
217 for(element_i = elements.begin(); element_i < elements.end(); element_i++) {
218 if ((*element_i)->getType()==type_el)
219 (*element_i)->printProperties();
220 }
221 return;
222}
223
224void H_AbstractBeamLine::showMatrix() const {
225 cout << "Transport matrix for the whole beam : " << endl;
226 cout << "(x,x',...) = (x*,x'*,...) M " <<endl;
227 printMatrix(beam_mat);
228 return;
229}
230
231void H_AbstractBeamLine::showMatrices() const{
232 // prints the list of all transport matrices, from the whole beam.
233
234 vector<TMatrix>::const_iterator matrix_i;
235 vector<H_OpticalElement*>::const_iterator element_i;
236 TMatrix temp(MDIM,MDIM);
237
238 for(matrix_i = matrices.begin(), element_i = elements.begin(); matrix_i < matrices.end(); matrix_i++, element_i++) {
239 temp = *matrix_i;
240 cout << "Matrix for transport until s=" << (*element_i)->getS() + (*element_i)->getLength() << "m (" << (*element_i)->getName() << "). " << endl;
[3c40083]241 printMatrix(&temp);
[5b822e5]242 cout << endl;
243 }
244 return ;
245}
246
247void H_AbstractBeamLine::calcSequence() {
248 // reorders the elements, computes the drifts;
249
250 vector<H_OpticalElement*> temp_elements;
251 vector<H_OpticalElement*>::iterator element_i;
252 // element_i is a pointer to elements[i]
253
254 if(elements.size()==1) { return; }
255
256 // getting rid of drifts before calculating
257 for(element_i = elements.begin(); element_i < elements.end(); element_i++) {
[3c40083]258 if((*element_i)->getType() == DRIFT) {elements.erase(element_i); }
[5b822e5]259 }
260
261 // ordering the elements in position
262 sort(elements.begin(),elements.end(),ordering());
263 // inserting the drifts before the other elements
264 float current_pos = 0;
265 float drift_length=0;
266
267 for(element_i=elements.begin(); element_i < elements.end(); element_i++) {
268 drift_length = (*element_i)->getS() - current_pos;
269 if(drift_length>0) {
270 H_Drift *dr = new H_Drift(current_pos,drift_length);
271 temp_elements.push_back(dr);
272 }
273 temp_elements.push_back(*element_i);
274 current_pos = (*element_i)->getS() + (*element_i)->getLength();
275 }
276
277 //adding the last drift
278 drift_length = beam_length - current_pos;
279 if (drift_length>0) {
280 H_Drift *dr = new H_Drift(current_pos,drift_length);
281 temp_elements.push_back(dr);
282 }
283 elements.clear();
284 for(element_i=temp_elements.begin(); element_i < temp_elements.end(); element_i++) {
285 elements.push_back(*element_i);
286 }
287}
288
289void H_AbstractBeamLine::calcMatrix() {
290 // computes the transport matrix for the beam upto here...
291 vector<H_OpticalElement*>::iterator element_i;
292 TMatrix calc_mat(MDIM,MDIM);
293
294 // initialization
295 matrices.clear();
296 calc_mat.UnitMatrix();
297
298 // multiplies the matrix of each beam's element
299 // and add each product matrix to the list of matrices.
300 for(element_i = elements.begin(); element_i < elements.end(); element_i++) {
301 calc_mat *= (*element_i)->getMatrix();
302 matrices.push_back(calc_mat);
303 }
304
[3c40083]305 *beam_mat = calc_mat;
[5b822e5]306 return;
307}
308
309float qh(float k) {
310 float beta = (log((float)10.0))/0.05;
311 // put (std::log((float)10.0)) instead of log(10) to avoid compilation errors
312 return 0.8*(1-exp(-beta*fabs(k)));
313}
314
315float dh(float k) {
316 float psi = (log((float)10.0))/0.002;
317 // put (std::log((float)10.0)) instead of log(10) to avoid compilation errors
318 return 0.8*(1-exp(-psi*fabs(k)));
319}
320
[3c40083]321void H_AbstractBeamLine::draw() const{
322/* gROOT->SetStyle("Plain");
323 TLegend* leg = new TLegend(0.85,0.50,1,1,"Legend");
[5b822e5]324 leg->SetBorderSize(1);
325 TBox* b1 = new TBox();
326 TBox* b2 = new TBox(0,0,10,10);
327 TBox* b3 = new TBox(0,0,0,0);
328 TBox* b4 = new TBox(0,0,0,0);
329 TBox* b5 = new TBox(0,0,0,0);
330 TBox* b6 = new TBox(0,0,0,0);
331 TBox* b7 = new TBox(0,0,0,0);
332 b1->SetFillColor(RDIPOLE);
333 b2->SetFillColor(SDIPOLE);
334 b3->SetFillColor(VQUADRUPOLE);
335 b4->SetFillColor(HQUADRUPOLE);
336 b5->SetFillColor(HKICKER);
337 b6->SetFillColor(VKICKER);
338 b7->SetFillColor(RCOLLIMATOR);
339 leg->AddEntry(b1,"R-Dipole");
340 leg->AddEntry(b2,"S-Dipole");
341 leg->AddEntry(b3,"V-Quadrupole");
342 leg->AddEntry(b4,"H-Quadrupole");
343 leg->AddEntry(b5,HKICKERNAME);
344 leg->AddEntry(b6,VKICKERNAME);
345 leg->AddEntry(b7,"RCollimator");
346 leg->Draw();
[3c40083]347*/
[5b822e5]348/* TLine* l1 = new TLine(0.05,0.5,0.95,0.5);
349 TLine* l2 = new TLine(0.1,0.1,0.1,0.9);
350 TLine* l3 = new TLine(0.9,0.1,0.9,0.9);
351 TPaveLabel* p1 = new TPaveLabel(0.05,0.5,0.1,0.6,"IP");
352 TPaveLabel* p2 = new TPaveLabel(0.9,0.5,0.95,0.6,"RP");
353 TGaxis* a1 = new TGaxis(0.1,0.1,0.9,0.1,0,beam_length);
354 a1->SetLabelSize(0.08);
355 p1->SetBorderSize(1);
356 p2->SetBorderSize(1);
357 p1->SetFillColor(0);
358 p2->SetFillColor(0);
359 l1->Draw();
360 l2->Draw();
361 l3->Draw();
362 p1->Draw();
363 p2->Draw();
364 float x1,x2,y1,y2;
365 vector<TPaveLabel*> boxes;
366 vector<H_OpticalElement*>::const_iterator element_i;
367 for(element_i = elements.begin(); element_i < elements.end(); element_i++) {
368 x1 = 0.1 + ((*element_i)->getS()/beam_length)*0.8;
369 x2 = x1 + ((*element_i)->getLength()/beam_length)*0.8;
370 if((*element_i)->getType()>5) {
371 y1 = 0.3;
372 y2 = 0.7;
373 }
374 else if((*element_i)->getType()>3) {
375 y1 = 0.5 - qh((*element_i)->getK()*(*element_i)->getLength())/2.;
376 y2 = 0.5 + qh((*element_i)->getK()*(*element_i)->getLength())/2.;
377 } else {
378 y1 = 0.5 - dh((*element_i)->getK()*(*element_i)->getLength())/2.;
379 y2 = 0.5 + dh((*element_i)->getK()*(*element_i)->getLength())/2.;
380 }
381 TPaveLabel* cur_box = new TPaveLabel(x1,y1,x2,y2,"");
382 cur_box->SetFillStyle(1);
383 cur_box->SetFillColor(((int)(*element_i)->getType()));
384 cur_box->SetBorderSize(1);
385 if((*element_i)->getType()!=DRIFT) boxes.push_back(cur_box);
386 }
387 vector<TPaveLabel*>::iterator box_i;
388 for(box_i = boxes.begin(); box_i < boxes.end(); box_i++) {
389 (*box_i)->Draw();
390 }
391 a1->Draw();
392*/
393 return;
394}
395
[3c40083]396void H_AbstractBeamLine::drawX(const float a_min, const float a_max) const{
[5b822e5]397 /// @param a_min defines the size of the drawing
398 /// @param a_max defines the size of the drawing
399 const int N = getNumberOfElements();
400 for(int i=0;i<N;i++) {
401 float height = fabs(a_max);
402 float meight = fabs(a_min);
403 float size = (height>meight)?meight:height;
[3c40083]404 float middle = getElement(i)->getX()*URAD;
[5b822e5]405 if(getElement(i)->getType()!=DRIFT) getElement(i)->draw(middle+size/2.,middle-size/2.);
406 }
407}
408
409void H_AbstractBeamLine::drawY(const float a_min, const float a_max) const{
410 /// @param a_min defines the size of the drawing
411 /// @param a_max defines the size of the drawing
412 const int N = getNumberOfElements();
413 for(int i=0;i<N;i++) {
414 float height = fabs(a_max);
415 float meight = fabs(a_min);
416 float size = (height>meight)?meight:height;
417 float middle = getElement(i)->getY()*URAD;
418 if(getElement(i)->getType()!=DRIFT) getElement(i)->draw(middle+size/2.,middle-size/2.);
419 }
420}
421
[3c40083]422void H_AbstractBeamLine::moveElement(const string name, const float new_s) {
[5b822e5]423 /// @param name identifies the element to move
424 /// @param new_s is where to put it
425 vector<H_OpticalElement*>::iterator element_i;
426 for(element_i = elements.begin(); element_i < elements.end(); element_i++) {
427 if(name==(*element_i)->getName()) { (*element_i)->setS(new_s); }
428 }
429
430 calcSequence();
431 calcMatrix();
432 return;
433}
434
[3c40083]435void H_AbstractBeamLine::alignElement(const string name, const float disp_x, const float disp_y) {
[5b822e5]436 /// @param name identifies the element to move
437 /// @param disp_x identifies the displacement to add in x [\f$ \mu m \f$]
438 /// @param disp_y identifies the displacement to add in y [\f$ \mu m \f$]
439 vector<H_OpticalElement*>::iterator element_i;
440 for(element_i = elements.begin(); element_i < elements.end(); element_i++) {
441 if(name==(*element_i)->getName()) {
442 (*element_i)->setX((*element_i)->getX()+disp_x);
443 (*element_i)->setY((*element_i)->getY()+disp_y);
444 return ;
445 }
446 }
[3c40083]447 cout<<"Element "<<name<<" not found."<<endl;
448 if(VERBOSE) cout<<"Element "<<name<<" not found."<<endl;
[5b822e5]449 return;
450}
451
[3c40083]452void H_AbstractBeamLine::tiltElement(const string name, const float ang_x, const float ang_y) {
[5b822e5]453 /// @param name identifies the element to move
454 /// @param ang_x identifies the angle to add in x
455 /// @param ang_y identifies the angle to add in y
456 vector<H_OpticalElement*>::iterator element_i;
457 for(element_i = elements.begin(); element_i < elements.end(); element_i++) {
458 if(name==(*element_i)->getName()) {
459 (*element_i)->setTX((*element_i)->getTX()+ang_x);
460 (*element_i)->setTY((*element_i)->getTY()+ang_y);
461 return ;
462 }
463 }
[3c40083]464 if(VERBOSE) cout<<"Element "<<name<<" not found."<<endl;
[5b822e5]465 return;
466}
467
468void H_AbstractBeamLine::offsetElements(const float start, const float offset) {
469 /// @param start After this s [m] coordinate, all elements will be offset.
470 /// @param offset In meters
471
472 extern int relative_energy;
473 if(!relative_energy) {
474 vector<H_OpticalElement*>::iterator element_i;
475 for(element_i = elements.begin(); element_i < elements.end(); element_i++) {
476 if((*element_i)->getS() > start ) {
477 (*element_i)->setX(offset);
478 }
479 }
480 }
481}
482
[3c40083]483/*
[5b822e5]484TGraph * H_AbstractBeamLine::getBetaX() const{
485 const int N = elements.size();
486 float * s = new float[N], * b = new float[N], temp;
487 int i=0, n=N;
488
489 vector<H_OpticalElement*>::const_iterator element_i;
490 for(element_i = elements.begin(); element_i < elements.end(); element_i++) {
491 temp=(*element_i)->getBetaX();
492 if (temp !=0) {
493 b[i] = (*element_i)->getBetaX();
494 s[i] = (*element_i)->getS();
495 i++;
496 n=i;
497 }
498 }
499
500 TGraph * betax = new TGraph(n,s,b);
501 betax->SetLineColor(1);
502 betax->SetLineStyle(2);
503 delete [] s;
504 delete [] b;
505 return betax;
506}
507
508TGraph * H_AbstractBeamLine::getBetaY() const{
509 const int N = elements.size();
510 float * s = new float[N], * b = new float[N], temp;
511 int i=0, n=N;
512
513 vector<H_OpticalElement*>::const_iterator element_i;
514 for(element_i = elements.begin(); element_i < elements.end(); element_i++) {
515 temp=(*element_i)->getBetaY();
516 if (temp !=0) {
517 b[i] = (*element_i)->getBetaY();
518 s[i] = (*element_i)->getS();
519 i++;
520 n=i;
521 }
522 }
523
524 TGraph * betay = new TGraph(n,s,b);
525 betay->SetLineColor(2);
526 betay->SetLineStyle(2);
527 delete [] s;
528 delete [] b;
529 return betay;
530}
531
532TGraph * H_AbstractBeamLine::getDX() const{
533 const int N = elements.size();
534 float * s = new float[N], * d = new float[N], temp;
535 int i=0, n=N;
536
537 vector<H_OpticalElement*>::const_iterator element_i;
538 for(element_i = elements.begin(); element_i < elements.end(); element_i++) {
539 temp=(*element_i)->getDX();
540 if (temp !=0) {
541 d[i] = (*element_i)->getDX();
542 s[i] = (*element_i)->getS();
543 i++;
544 n=i;
545 }
546 }
547
548 TGraph * dispx = new TGraph(n,s,d);
549 dispx->SetLineColor(8);
550 dispx->SetLineStyle(2);
551 delete [] s;
552 delete [] d;
553 return dispx;
554}
555
556TGraph * H_AbstractBeamLine::getDY() const{
557 const int N = elements.size();
558 float * s = new float[N], * d = new float[N], temp;
559 int i=0, n=N;
560
561 vector<H_OpticalElement*>::const_iterator element_i;
562 for(element_i = elements.begin(); element_i < elements.end(); element_i++) {
563 temp=(*element_i)->getDY();
564 if (temp !=0) {
565 d[i] = (*element_i)->getDY();
566 s[i] = (*element_i)->getS();
567 i++;
568 n=i;
569 }
570 }
571
572 TGraph * dispy = new TGraph(n,s,d);
573 dispy->SetLineColor(kBlue);
574 dispy->SetLineStyle(2);
575 delete [] s;
576 delete [] d;
577 return dispy;
578}
579
580
581TGraph * H_AbstractBeamLine::getRelX() const{
582 const int N = elements.size();
583 float * s = new float[N], * r = new float[N], temp;
584 int i=0, n=N;
585
586 vector<H_OpticalElement*>::const_iterator element_i;
587 for(element_i = elements.begin(); element_i < elements.end(); element_i++) {
588 temp=(*element_i)->getRelX();
589 if((*element_i)->getType() != DRIFT) {
590 r[i] = (*element_i)->getRelX();
591 s[i] = (*element_i)->getS();
592 i++;
593 n=i;
594 }
595 }
596
597 TGraph * relx = new TGraph(n,s,r);
598 relx->SetLineColor(kBlack);
599 relx->SetMarkerStyle(kOpenSquare);
600 relx->SetMarkerSize(0.6);
601 relx->SetLineStyle(2);
602 delete [] s;
603 delete [] r;
604 return relx;
605}
606
607TGraph * H_AbstractBeamLine::getRelY() const{
608 const int N = elements.size();
609 float * s = new float[N], * r = new float[N], temp;
610 int i=0, n=N;
611
612 vector<H_OpticalElement*>::const_iterator element_i;
613 for(element_i = elements.begin(); element_i < elements.end(); element_i++) {
614 temp=(*element_i)->getRelY();
615 if((*element_i)->getType() != DRIFT) {
616 r[i] = (*element_i)->getRelY();
617 s[i] = (*element_i)->getS();
618 i++;
619 n=i;
620 }
621 }
622
623 TGraph * rely = new TGraph(n,s,r);
624 rely->SetLineColor(kRed);
625 rely->SetMarkerStyle(kOpenSquare);
626 rely->SetMarkerSize(0.6);
627 rely->SetLineStyle(2);
628 delete [] s;
629 delete [] r;
630 return rely;
631}
[3c40083]632*/
Note: See TracBrowser for help on using the repository browser.