Fork me on GitHub

source: svn/trunk/external/Hector/H_AbstractBeamLine.cc@ 1364

Last change on this file since 1364 was 1360, checked in by Pavel Demin, 11 years ago

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