[3] | 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"
|
---|
| 22 | using namespace std;
|
---|
| 23 |
|
---|
| 24 | H_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 |
|
---|
| 32 | H_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 |
|
---|
| 41 | H_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 |
|
---|
| 52 | H_Beam::~H_Beam() {
|
---|
| 53 | beamParticles.clear();
|
---|
| 54 | return;
|
---|
| 55 | };
|
---|
| 56 |
|
---|
| 57 | void H_Beam::createBeamParticles(const unsigned int Number_of_particles) {
|
---|
| 58 | createBeamParticles(Number_of_particles,MP,QP);
|
---|
| 59 | }
|
---|
| 60 |
|
---|
| 61 | void 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) {
|
---|
| 78 | void 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 |
|
---|
| 107 | void 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 |
|
---|
| 119 | void 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 |
|
---|
| 132 | void 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 |
|
---|
| 144 | void 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 |
|
---|
| 156 | const 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 |
|
---|
| 162 | H_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 |
|
---|
| 168 | void H_Beam::add(const H_BeamParticle &p) {
|
---|
| 169 | beamParticles.push_back(p);
|
---|
| 170 | Nparticles++;
|
---|
| 171 | }
|
---|
| 172 |
|
---|
| 173 | void 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 |
|
---|
| 181 | void H_Beam::computePath(const H_AbstractBeamLine * beamline) {
|
---|
| 182 | computePath(beamline,false);
|
---|
| 183 | }
|
---|
| 184 |
|
---|
| 185 | /// Propagates the beam until a given s
|
---|
| 186 | void 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 |
|
---|
| 193 | void 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 |
|
---|
| 199 | void 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 |
|
---|
| 209 | float 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 |
|
---|
| 228 | float 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 |
|
---|
| 247 | TGraphErrors * 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 |
|
---|
| 267 | TGraphErrors * 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 |
|
---|
| 287 | float 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 |
|
---|
| 300 | float 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 |
|
---|
| 313 | unsigned 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 |
|
---|
| 322 | vector<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 |
|
---|
| 368 | void 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 |
|
---|
| 388 | void 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 |
|
---|
| 396 | void 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 |
|
---|
| 409 | TH2F * 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 |
|
---|
| 457 | TH2F * 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 |
|
---|
| 504 | TMultiGraph * 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 |
|
---|
| 516 | TMultiGraph * 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 | }
|
---|