Fork me on GitHub

source: svn/trunk/Utilities/Hector/src/H_Beam.cc@ 108

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

first commit
File size: 19.8 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
77//void H_Beam::particleGun(const unsigned int Number_of_particles, const float E_min=BE, const float E_max=BE, const float fs_min=0, const float fs_max=0, const float fx_min=0, const float fx_max=0, const float fy_min=0, const float fy_max=0, const float tx_min=-PI/2., const float tx_max=PI/2., const float ty_min=-PI/2., const float ty_max=PI/2., const float p_mass=MP, const double p_charge=QP) {
78void H_Beam::particleGun(const unsigned int Number_of_particles, const float E_min, const float E_max, const float fs_min, const float fs_max, const float fx_min, const float fx_max, const float fy_min, const float fy_max, const float tx_min, const float tx_max, const float ty_min, const float ty_max, const float p_mass, const double p_charge, const bool flat) {
79 beamParticles.clear();
80 Nparticles = (Number_of_particles<2) ? 2 : Number_of_particles;
81 float gx,gy,gs,gtx,gty,gE;
82 for (unsigned int i=0; i<Nparticles; i++) {
83 H_BeamParticle p(p_mass,p_charge);
84 if (flat) {
85 gx = gRandom->Uniform(fx_min,fx_max);
86 gy = gRandom->Uniform(fy_min,fy_max);
87 gs = gRandom->Uniform(fs_min,fs_max);
88 gtx = gRandom->Uniform(tx_min,tx_max);
89 gty = gRandom->Uniform(ty_min,ty_max);
90 gE = gRandom->Uniform(E_min,E_max);
91 } else {
92 gx = gRandom->Gaus((fx_min+fx_max)/2,(-fx_min+fx_max)/2);
93 gy = gRandom->Gaus((fy_min+fy_max)/2,(-fy_min+fy_max)/2);
94 gs = gRandom->Gaus((fs_min+fs_max)/2,(-fs_min+fs_max)/2);
95 gtx = gRandom->Gaus((tx_min+tx_max)/2,(-tx_min+tx_max)/2);
96 gty = gRandom->Gaus((ty_min+ty_max)/2,(-ty_min+ty_max)/2);
97 gE = gRandom->Gaus ((E_min+E_max)/2,(-E_min+E_max)/2);
98 }
99 p.setPosition(gx,gy,gtx,gty,gs);
100 p.setE(gE);
101 beamParticles.push_back(p);
102 }
103 return;
104}
105
106
107void H_Beam::createXScanningBeamParticles(const unsigned int Number_of_particles, const float fx_max) {
108 beamParticles.clear();
109 Nparticles = (Number_of_particles<2) ? 2 : Number_of_particles;
110 for (unsigned int i=0; i<Nparticles; i++) {
111 H_BeamParticle p;
112 float fx = fx_ini + i/(float)(Nparticles-1) * (fx_max-fx_ini);
113 p.setPosition(fx,fy_ini,0,0,fs_ini);
114 p.setE(fe_ini);
115 beamParticles.push_back(p);
116 }
117}
118
119void H_Beam::createYScanningBeamParticles(const unsigned int Number_of_particles, const float fy_max) {
120
121 beamParticles.clear();
122 Nparticles = (Number_of_particles<2) ? 2 : Number_of_particles;
123 for (unsigned int i=0; i<Nparticles; i++) {
124 H_BeamParticle p;
125 float fy = fy_ini + i/(float)(Nparticles-1) * (fy_max-fy_ini);
126 p.setPosition(fx_ini,fy,0,0,fs_ini);
127 p.setE(fe_ini);
128 beamParticles.push_back(p);
129 }
130}
131
132void H_Beam::createTXScanningBeamParticles(const unsigned int Number_of_particles, const float tx_max) {
133 beamParticles.clear();
134 Nparticles = (Number_of_particles<2) ? 2 : Number_of_particles;
135 for (unsigned int i=0; i<Nparticles; i++) {
136 H_BeamParticle p;
137 float tx = tx_ini + i/(float)(Nparticles-1) * (tx_max-tx_ini);
138 p.setPosition(fx_ini,fy_ini,tx,ty_ini,fs_ini);
139 p.setE(fe_ini);
140 beamParticles.push_back(p);
141 }
142}
143
144void H_Beam::createTYScanningBeamParticles(const unsigned int Number_of_particles, const float ty_max) {
145 beamParticles.clear();
146 Nparticles = (Number_of_particles<2) ? 2 : Number_of_particles;
147 for (unsigned int i=0; i<Nparticles; i++) {
148 H_BeamParticle p;
149 float ty = ty_ini + i/(float)(Nparticles-1) * (ty_max-ty_ini);
150 p.setPosition(fx_ini,fy_ini,tx_ini,ty,fs_ini);
151 p.setE(fe_ini);
152 beamParticles.push_back(p);
153 }
154}
155
156const H_BeamParticle * H_Beam::getBeamParticle(const unsigned int particle_index) const {
157// const int N = (particle_index<0)?0:(( particle_index>Nparticles)?Nparticles:particle_index);
158 const int N = (particle_index>Nparticles)?Nparticles:particle_index;
159 return &(*(beamParticles.begin()+N));// same as "return &beamParticles[N];" but more efficient
160}
161
162H_BeamParticle * H_Beam::getBeamParticle(const unsigned int particle_index) {
163// const int N = (particle_index<0)?0:(( particle_index>Nparticles)?Nparticles:particle_index);
164 const int N = (particle_index>Nparticles)?Nparticles:particle_index;
165 return &(*(beamParticles.begin()+N));// same as "return &beamParticles[N];" but more efficient
166}
167
168void H_Beam::add(const H_BeamParticle &p) {
169 beamParticles.push_back(p);
170 Nparticles++;
171}
172
173void H_Beam::computePath(const H_AbstractBeamLine * beamline, const bool NonLinear) {
174 vector<H_BeamParticle>::iterator particle_i;
175
176 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
177 particle_i->computePath(beamline,NonLinear);
178 }
179}
180
181void H_Beam::computePath(const H_AbstractBeamLine * beamline) {
182 computePath(beamline,false);
183}
184
185/// Propagates the beam until a given s
186void H_Beam::propagate(const float position) {
187 vector<H_BeamParticle>::iterator particle_i;
188 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
189 particle_i->propagate(position);
190 }
191}
192
193void H_Beam::emitGamma(const double gee, const double gq2) {
194 /// @param gee = \f$ E_{\gamma} \f$ is the photon energy
195 /// @param gq2 = \f$ Q^2 < 0 \f$ is virtuality of photon \f$ Q^{2} = E^{2}-\vec{k}^{2} \f$
196 emitGamma(gee,gq2,0,2*PI);
197}
198
199void H_Beam::emitGamma(const double gee, const double gq2, const double phimin, const double phimax) {
200 /// @param gee = \f$ E_{\gamma} \f$ is the photon energy
201 /// @param gq2 = \f$ Q^2 < 0 \f$ is virtuality of photon \f$ Q^{2} = E^{2}-\vec{k}^{2} \f$
202 /// @param phimin : lower bound for \f$ \phi \f$
203 /// @param phimax : higher bound for \f$ \phi \f$
204 vector<H_BeamParticle>::iterator particle_i;
205 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++)
206 particle_i->emitGamma(gee,gq2,phimin,phimax);
207}
208
209float H_Beam::getBetaX(const float s, float& error_on_betax) {
210 /// @param s is the position [m] to propagate to
211 /// @param error_on_betax : getBetaX(...) returns its error in this variable
212 /// not a const method because does a propagate to s!
213 vector<H_BeamParticle>::iterator particle_i;
214 float EX2=0,dummy, mean=getX(s,dummy), temp;
215
216 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
217 particle_i->propagate(s);
218 temp = particle_i->getX()-mean;
219 EX2 += temp*temp;
220 }
221 EX2 /= (float)Nparticles;
222 float emitx = getEmittanceX();
223 EX2 = (emitx==0)?0:(float) (EX2 /(float) (emitx*URAD))/URAD;
224 error_on_betax = EX2 / (float) sqrt((double)2*Nparticles);
225 return EX2;
226}
227
228float H_Beam::getBetaY(const float s, float& error_on_betay) {
229 /// @param s is the position [m] to propagate to
230 /// @param error_on_betay : getBetaY(...) returns its error in this variable
231 /// not a const method because does a propagate to s!
232 vector<H_BeamParticle>::iterator particle_i;
233 float EY2 =0, dummy, mean=getY(s,dummy), temp;
234
235 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
236 particle_i->propagate(s);
237 temp = particle_i->getY() - mean;
238 EY2 += temp*temp;
239 }
240 EY2 /= (float)Nparticles;
241 float emity = getEmittanceY();
242 EY2 = (emity==0)?0:(float) (EY2 / (float) (emity*URAD))/URAD;
243 error_on_betay = EY2 / (float) sqrt((double)2*Nparticles);
244 return EY2;
245}
246
247TGraphErrors * H_Beam::getBetaX(const float length, const unsigned int number_of_points) {
248 /// @param length [m]
249 /// @number_of_points in the graph (typ. 200)
250 const unsigned int N = number_of_points;
251 float * s = new float[N], * b = new float[N], * es = new float[N], * eb = new float[N];
252 for (unsigned int i=0; i<N; i++) {
253 s[i] = (float) fs_ini + i/(float)(N-1) *length;
254 b[i] = getBetaX(s[i],eb[i]);
255 es[i] = 0;
256 }
257 TGraphErrors * betax = new TGraphErrors(N,s,b,es,eb);
258 betax->SetLineColor(kBlack);
259 betax->SetFillColor(kYellow);
260 delete [] s;
261 delete [] b;
262 delete [] es;
263 delete [] eb;
264 return betax;
265}
266
267TGraphErrors * H_Beam::getBetaY(const float length, const unsigned int number_of_points) {
268 /// @param length [m]
269 /// @number_of_points in the graph (typ. 200)
270 const unsigned int N = number_of_points;
271 float * s = new float[N], * b = new float[N], * es = new float[N], *eb = new float[N];
272 for (unsigned int i=0; i<N; i++) {
273 s[i] = (float) fs_ini + i/(float)(N-1) *length;
274 b[i] = getBetaY(s[i],eb[i]);
275 es[i]=0;
276 }
277 TGraphErrors * betay = new TGraphErrors(N,s,b,es,eb);
278 betay->SetLineColor(kRed);
279 betay->SetFillColor(kYellow);
280 delete [] s;
281 delete [] b;
282 delete [] es;
283 delete [] eb;
284 return betay;
285}
286
287float H_Beam::getX(const float s, float& error_on_posx) {
288 vector<H_BeamParticle>::iterator particle_i;
289 float mean=0;
290
291 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
292 particle_i->propagate(s);
293 mean += particle_i->getX();
294 }
295 mean = mean / (float) Nparticles;
296 error_on_posx = mean / (float) sqrt((double)Nparticles);
297 return mean;
298}
299
300float H_Beam::getY(const float s, float& error_on_posy) {
301 vector<H_BeamParticle>::iterator particle_i;
302 float mean=0;
303
304 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
305 particle_i->propagate(s);
306 mean += particle_i->getY();
307 }
308 mean = mean / (float) Nparticles;
309 error_on_posy = mean / (float) sqrt((double)Nparticles);
310 return mean;
311}
312
313unsigned int H_Beam::getStoppedNumber(const H_AbstractBeamLine * beamline) {
314 int number =0;
315 vector<H_BeamParticle>::iterator particle_i;
316 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
317 if(particle_i->stopped(beamline)) number++;
318 }
319 return number;
320}
321
322vector<TVectorD> H_Beam::getStoppingElements(const H_AbstractBeamLine * beamline, vector<H_OpticalElement>& list, vector<int>& numb) {
323
324 vector<TVectorD> stop_positions;
325 vector<H_BeamParticle>::iterator particle_i;
326 vector<H_OpticalElement>::iterator element_i;
327 H_OpticalElement temp_el;
328 vector<int>::iterator n_i;
329 int number =0;
330 bool found;
331
332 list.clear();
333 numb.clear();
334
335 // creates a list of elements where beamParticles have stopped
336 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
337 found = false;
338 if(particle_i->stopped(beamline)) {
339 temp_el = *(particle_i->getStoppingElement());
340 stop_positions.push_back(*(particle_i->getStopPosition()));
341 if(list.size()==0) {
342 number=1;
343 list.push_back(temp_el);
344 numb.push_back(number);
345
346 } else {
347 for (element_i = list.begin(), n_i = numb.begin(); element_i < list.end(); element_i++, n_i++) {
348 string el_i_name = element_i->getName();
349 string temp_el_name = temp_el.getName();
350 if(el_i_name == temp_el_name) {
351 number = *n_i;
352 number++;
353 *n_i = number;
354 found = true;
355 }
356 }
357 if(!found) {
358 number=1;
359 list.push_back(temp_el);
360 numb.push_back(number);
361 }
362 }
363 } // if particle_i->stopped
364 }// for particle_i
365 return stop_positions;
366} // H_Beam::getStoppingElements
367
368void H_Beam::printInitialState() const {
369 cout << "Initial parameters of the beam" << endl;
370 cout << "(x,y,s) = (" << fx_ini << "," << fy_ini << "," << fs_ini << ") ";
371 cout << "(theta_x, theta_y) = (" << tx_ini << "," << ty_ini << ") ";
372 cout << "energy = " << fe_ini << endl;
373 cout << endl;
374 cout << "Dispersion on these values : " << endl;
375 cout << "(dx,dy,ds) = (" << x_disp << "," << y_disp << "," << s_disp << ") ";
376 cout << "(dtheta_x, dtheta_y) = (" << tx_disp << "," << ty_disp << ") ";
377 cout << "de = " << e_disp << endl << endl;
378
379 float mean_ini =0;
380 vector<H_BeamParticle>::const_iterator particle_i;
381 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
382 mean_ini += particle_i->getX();
383 }
384 mean_ini /= (float) beamParticles.size();
385 cout << "Mean ini x = " << mean_ini << endl;
386}
387
388void H_Beam::printProperties() const {
389 vector<H_BeamParticle>::const_iterator particle_i;
390 cout << "There are " << Nparticles << " in the beam." << endl;
391 for (particle_i = beamParticles.begin();particle_i < beamParticles.end(); particle_i++) {
392 particle_i->printProperties();
393 }
394}
395
396void H_Beam::printStoppingElements(const vector<H_OpticalElement>& list, const vector<int>& numb) const{
397 /// see also H_Beam::getStoppingElements
398 vector<H_OpticalElement>::const_iterator element_i;
399 vector<int>::const_iterator n_i;
400
401 // prints the list
402 for (element_i=list.begin(), n_i = numb.begin(); element_i < list.end(); element_i++, n_i++) {
403 cout << *n_i << " particules in " << element_i->getName();
404 cout << " (" << element_i->getTypeString() << ") at " << element_i->getS() << "m" << endl;
405 element_i->getAperture()->printProperties();
406 }
407} // H_Beam::printStoppingElements
408
409TH2F * H_Beam::drawAngleProfile(const float s) {
410 /// not a const method because does a propagate to s!
411 char title[50];
412 sprintf(title,"Beam profile at %.2f m",s);
413 vector<H_BeamParticle>::iterator particle_i;
414 float xmax, xmin, ymax, ymin;
415 float xx, yy, xborder, yborder;
416
417 particle_i=beamParticles.begin();
418 xmin = particle_i->getTX();
419 xmax = particle_i->getTX();
420 ymin = particle_i->getTY();
421 ymax = particle_i->getTY();
422
423 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
424 particle_i->propagate(s);
425 xx = particle_i->getTX();
426 yy = particle_i->getTY();
427
428 xmax = xx>xmax ? xx : xmax;
429 ymax = yy>ymax ? yy : ymax;
430 xmin = xx<xmin ? xx : xmin;
431 ymin = yy<ymin ? yy : ymin;
432 }
433
434 // in order to avoid some drawing problems, when the beam divergence is null
435 if(!(xmax || xmin)) xmax +=0.1;
436 if(!(ymax || ymin)) xmax +=0.1;
437
438 if(xmax == xmin) xmax *= 1.1;
439 if(ymax == ymin) ymax *= 1.1;
440
441 xborder = (xmax-xmin)*0.2;
442 yborder = (ymax-ymin)*0.2;
443
444 xmax += xborder;
445 xmin -= xborder;
446 ymax += yborder;
447 ymin -= yborder;
448
449 TH2F * profile = new TH2F("profile",title,10000,xmin,xmax,1000,ymin,ymax);
450 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
451 profile->Fill(particle_i->getTX(), particle_i->getTY());
452 }
453 return profile;
454}
455
456
457TH2F * H_Beam::drawProfile(const float s) {
458 /// not a const method because does a propagate to s!
459 char title[50];
460 sprintf(title,"Beam profile at %.2f m",s);
461 vector<H_BeamParticle>::iterator particle_i;
462 float xmax, xmin, ymax, ymin;
463 float xx, yy, xborder, yborder;
464
465 particle_i=beamParticles.begin();
466 xmin = particle_i->getX();
467 xmax = particle_i->getX();
468 ymin = particle_i->getY();
469 ymax = particle_i->getY();
470
471 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
472 particle_i->propagate(s);
473 xx = particle_i->getX();
474 yy = particle_i->getY();
475
476 xmax = xx>xmax ? xx : xmax;
477 ymax = yy>ymax ? yy : ymax;
478 xmin = xx<xmin ? xx : xmin;
479 ymin = yy<ymin ? yy : ymin;
480 }
481
482 // in order to avoid some drawing problems, when the beam divergence is null
483 if(!(xmax || xmin)) xmax +=0.1;
484 if(!(ymax || ymin)) xmax +=0.1;
485
486 if(xmax == xmin) xmax += 0.1;
487 if(ymax == ymin) ymax += 0.1;
488
489 xborder = (xmax-xmin)*0.2;
490 yborder = (ymax-ymin)*0.2;
491
492 xmax += xborder;
493 xmin -= xborder;
494 ymax += yborder;
495 ymin -= yborder;
496
497 TH2F * profile = new TH2F("profile",title,10000,xmin,xmax,1000,ymin,ymax);
498 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
499 profile->Fill(particle_i->getX(), particle_i->getY());
500 }
501 return profile;
502}
503
504TMultiGraph * H_Beam::drawBeamX(const int color) const {
505 int mycolor = color;
506 vector<H_BeamParticle>::const_iterator particle_i;
507 TMultiGraph * beam_profile_x = new TMultiGraph("beam_profile_x","");
508
509 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
510 TGraph * ppath_x = particle_i->getPath(0,mycolor);
511 beam_profile_x->Add(ppath_x);
512 }
513 return beam_profile_x;
514}
515
516TMultiGraph * H_Beam::drawBeamY(const int color) const {
517 int mycolor = color;
518 vector<H_BeamParticle>::const_iterator particle_i;
519 TMultiGraph * beam_profile_y = new TMultiGraph("beam_profile_y","");
520
521 for (particle_i = beamParticles.begin(); particle_i < beamParticles.end(); particle_i++) {
522 TGraph * ppath_y = particle_i->getPath(1,mycolor);
523 beam_profile_y->Add(ppath_y);
524 }
525 return beam_profile_y;
526}
Note: See TracBrowser for help on using the repository browser.