Fork me on GitHub

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

Last change on this file since 751 was 547, checked in by severine ovyn, 15 years ago

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