Fork me on GitHub

source: git/external/Hector/H_Beam.cc@ 0bbf248

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

switch to a more stable Hector version
  • Property mode set to 100644
File size: 15.9 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_Beam.cc
13/// \brief Describes a set a particles as a beam
14///
15
16// ROOT #includes
17//#include "TGraph.h"
18#include "TRandom.h"
19
20// local #includes
21#include "H_Beam.h"
22using namespace std;
23
24H_Beam::H_Beam() {
25 setPosition(PX,PY,TX+CRANG,TY,PS);
26 setE(BE);
27 setDispersion(SX,SY,STX,STY,SS);
28 setDE(SBE);
29 Nparticles=0;
30}
31
32H_Beam::H_Beam(const H_Beam& be) {
33 beamParticles = be.beamParticles;
34 setPosition(be.fx_ini,be.fy_ini,tx_ini,ty_ini,be.fs_ini);
35 setE(be.fe_ini);
36 setDispersion(be.x_disp,be.y_disp,be.tx_disp,be.ty_disp,be.s_disp);
37 setDE(be.e_disp);
38 Nparticles = be.Nparticles;
39}
40
41H_Beam& H_Beam::operator=(const H_Beam& be) {
42 if(this==&be) return *this;
43 beamParticles = be.beamParticles;
44 setPosition(be.fx_ini,be.fy_ini,tx_ini,ty_ini,be.fs_ini);
45 setE(be.fe_ini);
46 setDispersion(be.x_disp,be.y_disp,be.tx_disp,be.ty_disp,be.s_disp);
47 setDE(be.e_disp);
48 Nparticles = be.Nparticles;
49 return *this;
50}
51
52H_Beam::~H_Beam() {
53 beamParticles.clear();
54 return;
55};
56
57void H_Beam::createBeamParticles(const unsigned int Number_of_particles) {
58 createBeamParticles(Number_of_particles,MP,QP);
59}
60
61void H_Beam::createBeamParticles(const unsigned int Number_of_particles, const double p_mass, const double p_charge) {
62 beamParticles.clear();
63 Nparticles = (Number_of_particles<1) ? 1 : Number_of_particles;
64 for (unsigned int i=0; i<Nparticles; i++) {
65 H_BeamParticle p(p_mass,p_charge);
66 p.setPosition(fx_ini,fy_ini,tx_ini,ty_ini,fs_ini);
67 p.setE(fe_ini);
68 p.smearPos(x_disp,y_disp);
69 p.smearAng(tx_disp,ty_disp);
70 p.smearE(e_disp);
71 p.smearS(s_disp);
72 if (VERBOSE) {if (i==0) cout << " x_ini , tx_ini " << p.getX() << " " << p.getTX() << endl;}
73 beamParticles.push_back(p);
74 }
75}
76
77void H_Beam::createXScanningBeamParticles(const unsigned int Number_of_particles, const float fx_max) {
78 beamParticles.clear();
79 Nparticles = (Number_of_particles<2) ? 2 : Number_of_particles;
80 for (unsigned int i=0; i<Nparticles; i++) {
81 H_BeamParticle p;
82 float fx = fx_ini + i/(float)(Nparticles-1) * (fx_max-fx_ini);
83 p.setPosition(fx,fy_ini,0,0,fs_ini);
84 p.setE(fe_ini);
85 beamParticles.push_back(p);
86 }
87}
88
89void H_Beam::createYScanningBeamParticles(const unsigned int Number_of_particles, const float fy_max) {
90
91 beamParticles.clear();
92 Nparticles = (Number_of_particles<2) ? 2 : Number_of_particles;
93 for (unsigned int i=0; i<Nparticles; i++) {
94 H_BeamParticle p;
95 float fy = fy_ini + i/(float)(Nparticles-1) * (fy_max-fy_ini);
96 p.setPosition(fx_ini,fy,0,0,fs_ini);
97 p.setE(fe_ini);
98 beamParticles.push_back(p);
99 }
100}
101
102void H_Beam::createTXScanningBeamParticles(const unsigned int Number_of_particles, const float tx_max) {
103 beamParticles.clear();
104 Nparticles = (Number_of_particles<2) ? 2 : Number_of_particles;
105 for (unsigned int i=0; i<Nparticles; i++) {
106 H_BeamParticle p;
107 float tx = tx_ini + i/(float)(Nparticles-1) * (tx_max-tx_ini);
108 p.setPosition(fx_ini,fy_ini,tx,ty_ini,fs_ini);
109 p.setE(fe_ini);
110 beamParticles.push_back(p);
111 }
112}
113
114void H_Beam::createTYScanningBeamParticles(const unsigned int Number_of_particles, const float ty_max) {
115 beamParticles.clear();
116 Nparticles = (Number_of_particles<2) ? 2 : Number_of_particles;
117 for (unsigned int i=0; i<Nparticles; i++) {
118 H_BeamParticle p;
119 float ty = ty_ini + i/(float)(Nparticles-1) * (ty_max-ty_ini);
120 p.setPosition(fx_ini,fy_ini,tx_ini,ty,fs_ini);
121 p.setE(fe_ini);
122 beamParticles.push_back(p);
123 }
124}
125
126const H_BeamParticle * H_Beam::getBeamParticle(const unsigned int particle_index) const {
127// const int N = (particle_index<0)?0:(( particle_index>Nparticles)?Nparticles:particle_index);
128 const int N = (particle_index>Nparticles)?Nparticles:particle_index;
129 return &(*(beamParticles.begin()+N));// same as "return &beamParticles[N];" but more efficient
130}
131
132H_BeamParticle * H_Beam::getBeamParticle(const unsigned int particle_index) {
133// const int N = (particle_index<0)?0:(( particle_index>Nparticles)?Nparticles:particle_index);
134 const int N = (particle_index>Nparticles)?Nparticles:particle_index;
135 return &(*(beamParticles.begin()+N));// same as "return &beamParticles[N];" but more efficient
136}
137
138void H_Beam::add(const H_BeamParticle &p) {
139 beamParticles.push_back(p);
140 Nparticles++;
141}
142
143void H_Beam::computePath(const H_AbstractBeamLine * beamline, const bool NonLinear) {
144 vector<H_BeamParticle>::iterator particle_i;
145
146 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
147 particle_i->computePath(beamline,NonLinear);
148 }
149}
150
151void H_Beam::computePath(const H_AbstractBeamLine * beamline) {
152 computePath(beamline,false);
153}
154
155/// Propagates the beam until a given s
156void H_Beam::propagate(const float position) {
157 vector<H_BeamParticle>::iterator particle_i;
158 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
159 particle_i->propagate(position);
160 }
161}
162
163void H_Beam::emitGamma(const double gee, const double gq2) {
164 /// @param gee = \f$ E_{\gamma} \f$ is the photon energy
165 /// @param gq2 = \f$ Q^2 < 0 \f$ is virtuality of photon \f$ Q^{2} = E^{2}-\vec{k}^{2} \f$
166 emitGamma(gee,gq2,0,2*PI);
167}
168
169void H_Beam::emitGamma(const double gee, const double gq2, const double phimin, const double phimax) {
170 /// @param gee = \f$ E_{\gamma} \f$ is the photon energy
171 /// @param gq2 = \f$ Q^2 < 0 \f$ is virtuality of photon \f$ Q^{2} = E^{2}-\vec{k}^{2} \f$
172 /// @param phimin : lower bound for \f$ \phi \f$
173 /// @param phimax : higher bound for \f$ \phi \f$
174 vector<H_BeamParticle>::iterator particle_i;
175 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++)
176 particle_i->emitGamma(gee,gq2,phimin,phimax);
177}
178
179float H_Beam::getBetaX(const float s, float& error_on_betax) {
180 /// @param s is the position [m] to propagate to
181 /// @param error_on_betax : getBetaX(...) returns its error in this variable
182 /// not a const method because does a propagate to s!
183 vector<H_BeamParticle>::iterator particle_i;
184 float EX2=0,dummy, mean=getX(s,dummy), temp;
185
186 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
187 particle_i->propagate(s);
188 temp = particle_i->getX()-mean;
189 EX2 += temp*temp;
190 }
191 EX2 /= (float)Nparticles;
192 float emitx = getEmittanceX();
193 EX2 = (emitx==0)?0:(float) (EX2 /(float) (emitx*URAD))/URAD;
194 error_on_betax = EX2 / (float) sqrt((double)2*Nparticles);
195 return EX2;
196}
197
198float H_Beam::getBetaY(const float s, float& error_on_betay) {
199 /// @param s is the position [m] to propagate to
200 /// @param error_on_betay : getBetaY(...) returns its error in this variable
201 /// not a const method because does a propagate to s!
202 vector<H_BeamParticle>::iterator particle_i;
203 float EY2 =0, dummy, mean=getY(s,dummy), temp;
204
205 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
206 particle_i->propagate(s);
207 temp = particle_i->getY() - mean;
208 EY2 += temp*temp;
209 }
210 EY2 /= (float)Nparticles;
211 float emity = getEmittanceY();
212 EY2 = (emity==0)?0:(float) (EY2 / (float) (emity*URAD))/URAD;
213 error_on_betay = EY2 / (float) sqrt((double)2*Nparticles);
214 return EY2;
215}
216/*
217TGraphErrors * H_Beam::getBetaX(const float length, const unsigned int number_of_points) {
218 /// @param length [m]
219 /// @number_of_points in the graph (typ. 200)
220 const unsigned int N = number_of_points;
221 float * s = new float[N], * b = new float[N], * es = new float[N], * eb = new float[N];
222 for (unsigned int i=0; i<N; i++) {
223 s[i] = (float) fs_ini + i/(float)(N-1) *length;
224 b[i] = getBetaX(s[i],eb[i]);
225 es[i] = 0;
226 }
227 TGraphErrors * betax = new TGraphErrors(N,s,b,es,eb);
228 betax->SetLineColor(kBlack);
229 betax->SetFillColor(kYellow);
230 delete [] s;
231 delete [] b;
232 delete [] es;
233 delete [] eb;
234 return betax;
235}
236
237TGraphErrors * H_Beam::getBetaY(const float length, const unsigned int number_of_points) {
238 /// @param length [m]
239 /// @number_of_points in the graph (typ. 200)
240 const unsigned int N = number_of_points;
241 float * s = new float[N], * b = new float[N], * es = new float[N], *eb = new float[N];
242 for (unsigned int i=0; i<N; i++) {
243 s[i] = (float) fs_ini + i/(float)(N-1) *length;
244 b[i] = getBetaY(s[i],eb[i]);
245 es[i]=0;
246 }
247 TGraphErrors * betay = new TGraphErrors(N,s,b,es,eb);
248 betay->SetLineColor(kRed);
249 betay->SetFillColor(kYellow);
250 delete [] s;
251 delete [] b;
252 delete [] es;
253 delete [] eb;
254 return betay;
255}
256*/
257
258float H_Beam::getX(const float s, float& error_on_posx) {
259 vector<H_BeamParticle>::iterator particle_i;
260 float mean=0;
261
262 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
263 particle_i->propagate(s);
264 mean += particle_i->getX();
265 }
266 mean = mean / (float) Nparticles;
267 error_on_posx = mean / (float) sqrt((double)Nparticles);
268 return mean;
269}
270
271float H_Beam::getY(const float s, float& error_on_posy) {
272 vector<H_BeamParticle>::iterator particle_i;
273 float mean=0;
274
275 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
276 particle_i->propagate(s);
277 mean += particle_i->getY();
278 }
279 mean = mean / (float) Nparticles;
280 error_on_posy = mean / (float) sqrt((double)Nparticles);
281 return mean;
282}
283
284unsigned int H_Beam::getStoppedNumber(const H_AbstractBeamLine * beamline) {
285 int number =0;
286 vector<H_BeamParticle>::iterator particle_i;
287 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
288 if(particle_i->stopped(beamline)) number++;
289 }
290 return number;
291}
292
293void H_Beam::getStoppingElements(const H_AbstractBeamLine * beamline, vector<H_OpticalElement>& list, vector<int>& numb) {
294 vector<H_BeamParticle>::iterator particle_i;
295 vector<H_OpticalElement>::iterator element_i;
296 H_OpticalElement temp_el;
297 vector<int>::iterator n_i;
298 int number =0;
299 bool found;
300
301 list.clear();
302 numb.clear();
303
304 // creates a list of elements where beamParticles have stopped
305 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
306 found = false;
307 if(particle_i->stopped(beamline)) {
308 temp_el = *(particle_i->getStoppingElement());
309 if(list.size()==0) {
310 number=1;
311 list.push_back(temp_el);
312 numb.push_back(number);
313 } else {
314 for (element_i = list.begin(), n_i = numb.begin(); element_i < list.end(); element_i++, n_i++) {
315 string el_i_name = element_i->getName();
316 string temp_el_name = temp_el.getName();
317 if(el_i_name == temp_el_name) {
318 number = *n_i;
319 number++;
320 *n_i = number;
321 found = true;
322 }
323 }
324 if(!found) {
325 number=1;
326 list.push_back(temp_el);
327 numb.push_back(number);
328 }
329 }
330 } // if particle_i->stopped
331 }// for particle_i
332} // H_Beam::getStoppingElements
333
334void H_Beam::printInitialState() const {
335 cout << "Initial parameters of the beam" << endl;
336 cout << "(x,y,s) = (" << fx_ini << "," << fy_ini << "," << fs_ini << ") ";
337 cout << "(theta_x, theta_y) = (" << tx_ini << "," << ty_ini << ") ";
338 cout << "energy = " << fe_ini << endl;
339 cout << endl;
340 cout << "Dispersion on these values : " << endl;
341 cout << "(dx,dy,ds) = (" << x_disp << "," << y_disp << "," << s_disp << ") ";
342 cout << "(dtheta_x, dtheta_y) = (" << tx_disp << "," << ty_disp << ") ";
343 cout << "de = " << e_disp << endl << endl;
344
345 float mean_ini =0;
346 vector<H_BeamParticle>::const_iterator particle_i;
347 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
348 mean_ini += particle_i->getX();
349 }
350 mean_ini /= (float) beamParticles.size();
351 cout << "Mean ini x = " << mean_ini << endl;
352}
353
354void H_Beam::printProperties() const {
355 vector<H_BeamParticle>::const_iterator particle_i;
356 cout << "There are " << Nparticles << " in the beam." << endl;
357 for (particle_i = beamParticles.begin();particle_i < beamParticles.end(); particle_i++) {
358 particle_i->printProperties();
359 }
360}
361
362void H_Beam::printStoppingElements(const vector<H_OpticalElement>& list, const vector<int>& numb) const{
363 /// see also H_Beam::getStoppingElements
364 vector<H_OpticalElement>::const_iterator element_i;
365 vector<int>::const_iterator n_i;
366
367 // prints the list
368 for (element_i=list.begin(), n_i = numb.begin(); element_i < list.end(); element_i++, n_i++) {
369 cout << *n_i << " particules in " << element_i->getName();
370 cout << " (" << element_i->getTypeString() << ") at " << element_i->getS() << "m" << endl;
371 element_i->getAperture()->printProperties();
372 }
373} // H_Beam::printStoppingElements
374/*
375TH2F * H_Beam::drawProfile(const float s) {
376 /// not a const method because does a propagate to s!
377 char title[50];
378 sprintf(title,"Beam profile at %.2f m",s);
379 vector<H_BeamParticle>::iterator particle_i;
380 float xmax, xmin, ymax, ymin;
381 float xx, yy, xborder, yborder;
382
383 particle_i=beamParticles.begin();
384 xmin = particle_i->getX();
385 xmax = particle_i->getX();
386 ymin = particle_i->getY();
387 ymax = particle_i->getY();
388
389 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
390 particle_i->propagate(s);
391 xx = particle_i->getX();
392 yy = particle_i->getY();
393
394 xmax = xx>xmax ? xx : xmax;
395 ymax = yy>ymax ? yy : ymax;
396 xmin = xx<xmin ? xx : xmin;
397 ymin = yy<ymin ? yy : ymin;
398 }
399
400 // in order to avoid some drawing problems, when the beam divergence is null
401 if(xmax == xmin) xmax += 0.1;
402 if(ymax == ymin) ymax += 0.1;
403
404 xborder = (xmax-xmin)*0.2;
405 yborder = (ymax-ymin)*0.2;
406
407 xmax += xborder;
408 xmin -= xborder;
409 ymax += yborder;
410 ymin -= yborder;
411
412 TH2F * profile = new TH2F("profile",title,10000,xmin,xmax,1000,ymin,ymax);
413 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
414 profile->Fill(particle_i->getX(), particle_i->getY());
415 }
416 return profile;
417}*/
418/*
419TMultiGraph * H_Beam::drawBeamX(const int color) const {
420 int mycolor = color;
421 vector<H_BeamParticle>::const_iterator particle_i;
422 TMultiGraph * beam_profile_x = new TMultiGraph("beam_profile_x","");
423
424 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
425 TGraph * ppath_x = particle_i->getPath(0,mycolor);
426 beam_profile_x->Add(ppath_x);
427 }
428 return beam_profile_x;
429}
430
431TMultiGraph * H_Beam::drawBeamY(const int color) const {
432 int mycolor = color;
433 vector<H_BeamParticle>::const_iterator particle_i;
434 TMultiGraph * beam_profile_y = new TMultiGraph("beam_profile_y","");
435
436 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
437 TGraph * ppath_y = particle_i->getPath(1,mycolor);
438 beam_profile_y->Add(ppath_y);
439 }
440 return beam_profile_y;
441}
442*/
Note: See TracBrowser for help on using the repository browser.