Fork me on GitHub

source: git/external/Hector/H_BeamLine.cc@ cd616b9

ImprovedOutputFile Timing dual_readout llp
Last change on this file since cd616b9 was 3c40083, checked in by pavel <pavel@…>, 11 years ago

switch to a more stable Hector version

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