Fork me on GitHub

source: svn/trunk/Utilities/Hector/src/H_BeamLine.cc@ 425

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

new Hector version
File size: 9.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_BeamLine.cc
20/// \brief Reads external files and retrieves features of the real beam optical elements
21
22// c++ #includes
23#include <iostream>
24#include <fstream>
25#include <sstream>
26
27// local #includes
28#include "H_BeamLine.h"
29#include "H_BeamLineParser.h"
30#include "H_RectEllipticAperture.h"
31#include "H_RectangularAperture.h"
32#include "H_EllipticAperture.h"
33#include "H_CircularAperture.h"
34#include "H_RectangularDipole.h"
35#include "H_SectorDipole.h"
36#include "H_HorizontalQuadrupole.h"
37#include "H_VerticalQuadrupole.h"
38#include "H_HorizontalKicker.h"
39#include "H_VerticalKicker.h"
40#include "H_RectangularCollimator.h"
41#include "H_Marker.h"
42
43using namespace std;
44
45// Caution : mad conventions for vertically(horizontally) focusing quadrupoles
46// are opposite to Wille's. See fill() for more details.
47
48H_BeamLine::H_BeamLine(): H_AbstractBeamLine(),
49 direction(1), ips(0), ipx(0), ipy(0), iptx(0), ipty(0) {
50}
51
52H_BeamLine::H_BeamLine(const int si, const float length) : H_AbstractBeamLine(length),
53 direction((si >= abs(si)) ? 1 : -1), ips(0), ipx(0), ipy(0), iptx(0), ipty(0) {
54}
55
56H_BeamLine::H_BeamLine(const H_BeamLine& beam) : H_AbstractBeamLine(beam),
57 direction(beam.direction), ips(beam.ips), ipx(beam.ipx), ipy(beam.ipy), iptx(beam.iptx), ipty(beam.ipty) {
58}
59
60H_BeamLine& H_BeamLine::operator=(const H_BeamLine& beam) {
61 if(this==&beam) return *this;
62 H_AbstractBeamLine::operator=(beam); // call the mother's operator=
63 direction = beam.direction;
64 ips = beam.ips;
65 ipx = beam.ipx;
66 ipy = beam.ipy;
67 iptx = beam.iptx;
68 ipty = beam.ipty;
69 return *this;
70}
71
72void H_BeamLine::findIP(const string& filename, const string& ipname) {
73 // searches for the IP position in the extended table.
74 ifstream tabfile(filename.c_str());
75 if (! tabfile.is_open()) cout << "<H_BeamLine> ERROR: I Can't open \"" << filename << "\"" << endl;
76 bool found = false;
77 int N_col=0;
78 string headers[40]; // table of all column headers
79 int header_type[40]; // table of all column types
80
81 string temp_string;
82 H_BeamLineParser e;
83 istringstream curstring;
84
85 while ( getline(tabfile,temp_string)) {
86 curstring.clear(); // needed when using istringstream::str(string) several times !
87 curstring.str(temp_string);
88
89 // gets the features of each element
90 if (found) {
91 for (int col=0; col < N_col; col++) e.setProperties(curstring,header_type[col]);
92 if(strstr(e.name.c_str(),ipname.c_str())) {
93 ips = e.s; //e.printProperties();
94 ipx = e.x;
95 ipy = e.y;
96 iptx = e.px*URAD;
97 ipty = e.py*URAD;
98 }
99 } else if (strstr(temp_string.c_str(),"K0L")) { //searches for the beginning of the element list.
100 found = true;
101
102 // reads the title of each column
103 while(curstring.good()) { curstring >> headers[N_col]; if(headers[N_col]!="*") N_col++;}
104 if(VERBOSE) cout << N_col << " columns identified" << endl;
105
106 // identifies each column
107 for (int col=0; col< N_col; col++) header_type[col] = column_identification(headers[col]);
108
109 getline(tabfile,temp_string); // skip one line
110 } // if(... "K0L")
111 } // while (!eof)
112
113 if (!found) cout << "<H_BeamLine> ERROR ! IP not found." << endl;
114 tabfile.close();
115}
116
117double* H_BeamLine::getIPProperties() {
118 double* temp = new double[4];
119 temp[0] = ipx;
120 temp[1] = ipy;
121 temp[2] = iptx;
122 temp[3] = ipty;
123 return temp;
124}
125
126
127void H_BeamLine::fill(const string& filename, const int dir, const string& ipname) {
128 string headers[40]; // table of all column headers
129 int header_type[40]; // table of all column types
130 findIP(filename,ipname);
131 ifstream tabfile(filename.c_str());
132 if (! tabfile.is_open()) cout << "<H_BeamLine> ERROR: I Can't open \"" << filename << "\"" << endl;
133 if(VERBOSE) {
134 cout<<"Using file : "<< filename <<" in the "<<((direction>0)?"positive":"negative")<<" direction."<<endl;
135 cout<<"IP was found at position : "<<ips<<endl<<endl;
136 }
137 bool found = false; // found = true when the beginning of the element list is reached
138 int type, N_col=0;
139 float previous_betax =0, previous_betay=0; // needed when direction ==1
140 float previous_dx =0, previous_dy=0; // needed when direction ==1
141 float previous_x =0, previous_y=0; // needed when direction ==1
142
143 string temp_string;
144 H_BeamLineParser e;
145 istringstream curstring;
146 H_OpticalElement *el = 0;
147 H_Aperture * ap = 0;
148
149
150 while (getline(tabfile,temp_string)) {
151 curstring.clear(); // needed when using several tims istringstream::str(string)
152 curstring.str(temp_string);
153
154 // gets the features of each element
155 if (found) {
156 // reads each column
157 for (int col=0; col < N_col; col++) {
158 e.setProperties(curstring,header_type[col]);
159 }
160 //e.printProperties();
161
162 type =0; //init
163 if(e.k0l!=0) { if(strstr(e.name.c_str(),"MB.")) type = SDIPOLE; else type = RDIPOLE;} //all SBEND seems to be called MB.xxxxx.xx
164 else if(e.k1l!=0) type = (e.k1l>0)?HQUADRUPOLE:VQUADRUPOLE;
165 else if(e.hkick!=0) type = HKICKER;
166 else if(e.vkick!=0) type = VKICKER;
167 else if(strstr(e.name.c_str(),"DRIFT")) type=DRIFT;
168 else if((strstr(e.name.c_str(),"DFB"))||(strstr(e.name.c_str(),"TAS"))||(strstr(e.name.c_str(),"TAN"))||(strstr(e.name.c_str(),"TCL"))) type = RCOLLIMATOR;
169 else if(strstr(e.name.c_str(),ipname.c_str())) { type=IP; }
170 else type=MARKER;
171
172 e.s = (e.s > ips) ? (e.s -ips - e.l)*(direction) : (e.s-ips)*(direction);
173
174 el=0; //init
175 switch(type) {
176 case RDIPOLE: { el = new H_RectangularDipole(e.name.c_str(),e.s,dir*e.k0l/e.l,e.l); } break;
177 case SDIPOLE: { el = new H_SectorDipole(e.name.c_str(),e.s,dir*e.k0l/e.l,e.l); } break;
178 case VQUADRUPOLE: { el = new H_VerticalQuadrupole(e.name.c_str(),e.s,-(e.k1l/e.l),e.l); } break;
179 case HQUADRUPOLE: { el = new H_HorizontalQuadrupole(e.name.c_str(),e.s,-(e.k1l/e.l),e.l); } break;
180 case VKICKER: { el = new H_VerticalKicker(e.name.c_str(),e.s,e.vkick,e.l); } break;
181 case HKICKER: { el = new H_HorizontalKicker(e.name.c_str(),e.s,e.hkick,e.l); } break;
182 case RCOLLIMATOR: { el = new H_RectangularCollimator(e.name.c_str(),e.s,e.l); } break;
183 case DRIFT: {previous_betax = e.betx; previous_betay = e.bety; previous_dx = e.dx; previous_dy = e.dy; previous_x = e.x; previous_y = e.y;} break;
184 case IP: { el = new H_Marker(e.name.c_str(),e.s); } break;
185 default: break;
186 }
187
188 if(el!=0) {
189 ap = 0; //init
190 if(e.aper_1 !=0 || e.aper_2 !=0 || e.aper_3 !=0 || e.aper_4 != 0) {
191 e.aper_1 *= URAD;
192 e.aper_2 *= URAD;
193 e.aper_3 *= URAD;
194 e.aper_4 *= URAD; // in [m] in the tables !
195
196 if(strstr(e.apertype.c_str(),"RECTELLIPSE"))
197 ap = new H_RectEllipticAperture(e.aper_1,e.aper_2,e.aper_3,e.aper_4,0,0);
198 else if(strstr(e.apertype.c_str(),"CIRCLE"))
199 ap = new H_CircularAperture(e.aper_1,0,0);
200 else if(strstr(e.apertype.c_str(),"RECTANGLE"))
201 ap = new H_RectangularAperture(e.aper_1,e.aper_2,0,0);
202 else if(strstr(e.apertype.c_str(),"ELLIPSE"))
203 ap = new H_EllipticAperture(e.aper_1,e.aper_2,0,0);
204 }
205
206 if (direction<0) {
207 el->setBetaX(e.betx); el->setBetaY(e.bety);
208 el->setDX(e.dx); el->setDY(e.dy);
209 el->setRelX(e.x); el->setRelY(e.y);
210 } else {
211 el->setBetaX(previous_betax); el->setBetaY(previous_betay);
212 el->setDX(previous_dx); el->setDY(previous_dy);
213 el->setRelX(previous_x); el->setRelY(previous_y);
214 }
215 if(ap!=0) {
216 el->setAperture(ap);
217 delete ap;
218 // if memory is allocated, i.e. if ap!=0,
219 // it wont be deallocated in H_OpticalElement::~H_OpticalElement
220 // ap should then be deleted here
221 } // memory leak-free!
222
223 /// Parses all the elements from the beginning of the file,
224 /// but only keeps the ones from the IP till the desired length
225 if(e.s>=0 && e.s<beam_length) add(el);
226 else { delete el;}
227 // NB : if "el" is added to the beamline, it will be
228 // deleted by H_AbstractBeamLine::~H_AbstractBeamLine
229 // Otherwise, it should be deleted explicitly
230
231 } // if(el!=0)
232 } // if (found)
233 else if(strstr(temp_string.c_str(),"K0L")) { // if (!found)
234 //searches for the beginning of the element list.
235 found = true;
236
237 // reads the title of each column
238 while(curstring.good()) { curstring >> headers[N_col]; if(headers[N_col]!="*") N_col++;}
239 if(VERBOSE) cout << N_col << " columns identified" << endl;
240
241 // identifies each column
242 for (int col=0; col< N_col; col++) {
243 header_type[col] = column_identification(headers[col]);
244 }
245 getline(tabfile,temp_string);
246 } // if temp_string <-> "K0L"
247
248 } // while (!eof)
249 tabfile.close();
250}
251
Note: See TracBrowser for help on using the repository browser.