Fork me on GitHub

source: git/external/Hector/H_Beam.cc@ c526439

ImprovedOutputFile Timing dual_readout llp
Last change on this file since c526439 was 5b822e5, checked in by pavel <pavel@…>, 11 years ago

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