Fork me on GitHub

source: git/external/Hector/H_BeamParticle.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: 19.8 KB
RevLine 
[3c40083]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*/
[5b822e5]11
12/// \file H_BeamParticle.cc
13/// \brief Class aiming at simulating a particle in the LHC beam
14
15// from IP to RP, with emission of a photon of defined energy and Q.
16// Units : angles [rad], distances [m], energies [GeV], c=[1].
17// !!! no comment statement at the end of a #define line !!!
18
19// c++ #includes
20#include <iostream>
21#include <iomanip>
22#include <vector>
23#include <string>
24#include <cmath>
25
26// ROOT #includes
27#include "H_Parameters.h"
28#include "TRandom.h"
[3c40083]29//#include "TView.h"
30//#include "TPolyLine3D.h"
[5b822e5]31#ifdef _include_pythia_
32#include "TPythia6.h"
[3c40083]33#include "TRandom.h"
[5b822e5]34#endif
35// local #includes
36#include "H_OpticalElement.h"
37#include "H_BeamParticle.h"
38#include "H_Drift.h"
39
40using namespace std;
41
[3c40083]42void H_BeamParticle::init() {
[5b822e5]43 mp = MP;
44 qp = QP;
45 fx = 0;
46 fy = 0;
47 thx = 0;
48 thy = 0;
49 fs = 0;
50 energy = BE;
51 hasstopped = false;
52 hasemitted = false;
53 isphysical = true;
54 addPosition(fx,thx,fy,thy,fs);
55 stop_position = new TVectorD(LENGTH_VEC);
56 for (int i=0; i<LENGTH_VEC; i++) (*stop_position)[i] = -1;
57 stop_element = 0;
58}
59
[3c40083]60H_BeamParticle::H_BeamParticle() {
[5b822e5]61 init();
62}
63
[3c40083]64H_BeamParticle::H_BeamParticle(const H_BeamParticle& p) {
65 mp = p.mp;
66 qp = p.qp;
67 fx = p.fx;
68 fy = p.fy;
69 thx = p.thx;
70 thy = p.thy;
71 fs = p.fs;
72 energy = p.energy;
73 hasstopped = p.hasstopped;
74 hasemitted = p.hasemitted;
75 isphysical = p.isphysical;
76 stop_position = new TVectorD(*(p.stop_position));
[5b822e5]77 if(p.hasstopped) stop_element = new H_OpticalElement(*(p.stop_element));
[3c40083]78 positions = p.positions;
[5b822e5]79}
80
81H_BeamParticle::H_BeamParticle(const double mass, const double charge) {
82 init();
83 mp = mass;
84 qp = (mass==0) ? 0 : charge;
85 /// rejects particles with mass = 0 and charge != 0
86}
87
88H_BeamParticle& H_BeamParticle::operator=(const H_BeamParticle& p) {
89 if(this==&p) return *this;
90 mp = p.mp;
91 qp = p.qp;
92 fx = p.fx;
93 fy = p.fy;
94 thx = p.thx;
95 thy = p.thy;
96 fs = p.fs;
97 energy = p.energy;
98 hasstopped = p.hasstopped;
99 hasemitted = p.hasemitted;
100 isphysical = p.isphysical;
101 stop_position = new TVectorD(*(p.stop_position));
102 if(p.hasstopped) stop_element = new H_OpticalElement(*(p.stop_element));
103 positions = p.positions;
104 return *this;
105}
106
[3c40083]107bool H_BeamParticle::stopped(const H_AbstractBeamLine * beamline) {
[5b822e5]108 vector<TVectorD>::const_iterator position_i;
109 for(position_i = positions.begin(); position_i < positions.end()-1; position_i++) {
110 const unsigned int pos = position_i-positions.begin();
111 if(beamline->getElement(pos)->getAperture()->getType()!=NONE) {
112 bool has_passed_entrance = beamline->getElement(pos)->isInside((*position_i)[INDEX_X],(*position_i)[INDEX_Y]);
[3c40083]113 bool has_passed_exit = beamline->getElement(pos)->isInside((*(position_i+1))[INDEX_X],(*(position_i+1))[INDEX_Y]);
[5b822e5]114 if(!(has_passed_entrance && has_passed_exit)) {
[3c40083]115 // cout << "p =" << (*position_i)[INDEX_X] << "; el=" << beamline->getElement(pos)->getX()<<endl;
116 // beamline->getElement(position_i-positions.begin())->printProperties();
[5b822e5]117 if(VERBOSE) cout<<"particle stopped at "<<(beamline->getElement(pos))->getName();
118 if(VERBOSE) cout<<" (s = "<<(*position_i)[4] << ")" << endl;
[3c40083]119 if(VERBOSE && !has_passed_exit) cout << "Particle stopped inside the element" << endl;
[5b822e5]120 hasstopped=true;
121 stop_element = const_cast<H_OpticalElement*>(beamline->getElement(pos));
122 *stop_position = *position_i;
123 return hasstopped;
124 } // if
[3c40083]125// else cout << "outside aperture " << endl;
[5b822e5]126 } // if
127 } // for
128 return hasstopped;
129}
130
131void H_BeamParticle::addPosition(const double x, const double tx, const double y, const double ty, const double s) {
132 // [x] = [y] = m ; [tx] = [ty] = rad; [s] = m
133 double xys[LENGTH_VEC];
134 xys[INDEX_X]=x;
135 xys[INDEX_TX]=tx;
136 xys[INDEX_Y]=y;
137 xys[INDEX_TY]=ty;
138 xys[INDEX_S]=s;
139
140 TVectorD temp_vec(LENGTH_VEC,xys);
141 positions.push_back(temp_vec);
142}
143
144void H_BeamParticle::smearPos(const double dx,const double dy, TRandom* r) {
145 // the beam is centered on (fx,fy) at IP
146 fx = r->Gaus(fx,dx);
147 fy = r->Gaus(fy,dy);
148 positions.clear();
149 addPosition(fx,thx,fy,thy,fs);
150 return;
151}
152
153void H_BeamParticle::smearAng(const double tx, const double ty, TRandom* r) {
154 // the beam transverse direction is centered on (thx,thy) at IP
155 thx = r->Gaus(thx,tx);
156 thy = r->Gaus(thy,ty);
157 positions.clear();
158 addPosition(fx,thx,fy,thy,fs);
159 return;
160}
161
162void H_BeamParticle::smearE(const double erre, TRandom* r) {
163 energy = r->Gaus(energy,erre);
164 return;
165}
166
167void H_BeamParticle::smearS(const double errs, TRandom* r) {
168 fs= r->Gaus(fs,errs);
169 positions.clear();
170 addPosition(fx,thx,fy,thy,fs);
171 return;
172}
173
174void H_BeamParticle::set4Momentum(const double px, const double py, const double pz, const double ene) {
175 /// @param px, py, pz, ene is \f$ (\vec p , E) [GeV]\f$
176 ///
177 /// Clears the H_BeamParticle::positions vector.
178 positions.clear();
179 if(pz==0) {
180 cout<<" ERROR in H_BeamParticle::set4Momentum : no momentum in the beamline direction !"<<endl;
181 return;
182 }
183 thx = thx + URAD*atan(px/pz);
184 thy = thy + URAD*atan(py/pz);
185 energy = ene;
186 positions.clear();
187 addPosition(fx,thx,fy,thy,fs);
188 return;
189}
190
191void H_BeamParticle::setE(const double ene) {
192 energy = ene;
193 return;
194}
195
196void H_BeamParticle::setPosition(const double x, const double y, const double tx, const double ty, const double s) {
197 /// @param x, y are the transverse positions in \f$ \mu \f$ m
198 /// @param tx, ty are the angles in \f$ \mu \f$ rad
199 /// @param s is the longitudinal coordinate in m
200 // clear positions and sets the initial one.
201 fx=x;
202 fy=y;
203 thx=tx;
204 thy=ty;
205 fs = s;
206 positions.clear();
207 addPosition(fx,thx,fy,thy,s);
208
209 return;
210}
211
212const H_OpticalElement* H_BeamParticle::getStoppingElement() const{
213 if(hasstopped) return stop_element;
214 else { H_OpticalElement * dummy_el = new H_Drift("",0,0); return dummy_el;}
215}
216
[3c40083]217void H_BeamParticle::emitGamma(const double gee, const double gq2) {
218 emitGamma(gee,gq2,0,2*PI);
219 return;
220}
221
[5b822e5]222void H_BeamParticle::emitGamma(const double gee, const double gq2, const double phimin, const double phimax) {
223 /// @param gee = \f$ E_{\gamma} \f$ is the photon energy
224 /// @param gq2 = Q < 0 is virtuality of photon \f$ Q^{2} = E^{2}-\vec{k}^{2} \f$
225 /// @param phimin : lower bound for \f$ phi \f$
226 /// @param phimax : higher bound for \f$ phi \f$
227
228 if(gq2==0) {
229 if(VERBOSE) cout<<"No virtuality : only energy has changed"<<endl;
230 setE(energy-gee);
231 return;
232 }
233
234 double m1 = mp;
235
236 double e1 = energy , e2 = energy - gee; // particle energy : before (1) / after (2)
237
238 double p1 = sqrt(pow(e1,2) - pow(m1,2)), p2 = sqrt(pow(e2,2) - pow(m1,2)); // particle momentum : before (1) / after (2)
239 double q2min = pow(gee,2) - pow(p1+p2,2); // lower bound from kinematics E^2 - (p1 + p2)^2
240 double q2max = -2 * pow(m1*gee/(p1+p2),2) * (1 + (pow(e1,2) + pow(e2,2) -pow(m1,2))/(e1*e2 + p1*p2) );
241 // upper bound from kinematics ; E^2 - (p1-p2)^2; is bad for numerical computations
242
243
244 // if q2min < q2 < q2max is NOT true, there will be mathematical problems (like cos(eta) > 1).
245 // so if q2 > q2max, we force q2 = q2max (-> cos(eta) = 1)
246 // and if q2 < q2min, we force q2 = q2min (-> cos(eta) = 1)
247 // BUT the user knows something was wrong with the value of "H_BeamParticle::isphysical"
248 const double q2 = (gq2 > q2max ) ? q2max : (gq2 < q2min) ? q2min : gq2;
249
250 if( (gq2>q2max) || (gq2<q2min)) {
[3c40083]251 if(VERBOSE) cout<<"Non physical particle ! Q2 (" << q2 << ") and E ("<<gee << ") are not compatible." << endl;
[5b822e5]252 isphysical = false;
253 }
254
[3c40083]255 if(hasemitted) { cout<<"particle has already emitted at least one gamma !"<<endl;}
[5b822e5]256 hasemitted = true;
257 energy = energy - gee;
258 // gkk is k
259 double gkk = sqrt(pow(gee,2)-q2);
260 // eta is the angle between gamma and initial direction of the gamma-emitting particle
261 // ceta = cos(eta) and seta = sin(eta)
262
263 double ceta = sqrt( pow(mp/p1,2) + 1 ) * sqrt( q2/pow(gkk,2) + 1 ) - q2/(2*p1*gkk);
264 double seta = sqrt(1 - ceta*ceta);
265 // theta is the angle between particle and beam
266 double theta = URAD*atan(seta/(BE/gkk - ceta));
267 double phi = phimin + gRandom->Uniform(phimax-phimin);
268 thx = thx + theta*cos(phi);
269 thy = thy - theta*sin(phi);
270
271 // caution : emitting a photon erases all known positions !
272 positions.clear();
273 addPosition(fx,thx,fy,thy,fs);
274 return;
275}
276
277void H_BeamParticle::doInelastic() {
278#ifdef _include_pythia_
279// if(!gROOT->GetClass("TPythia6")) {
280// gROOT->Reset();
281// gSystem->Load("libPythia6");
282// gSystem->Load("libEG");
283// gSystem->Load("libEGPythia6");
284// }
285
286 TPythia6 gen;
287 // select AB -> AX process
288 gen.SetMSEL(0);
289 gen.SetMSUB(93,1);
290 // no showers/decays
291 gen.SetMSTP(111,0);
292 // no printouts
293 gen.SetMSTP(122,0);
294 gen.SetMSTU(12,0);
295 // generator initialization
296 gen.Initialize("CMS","p","p",14000);
297 // event generation
298 gen.GenerateEvent();
299 // list particles
300// gen.Pylist(1);
301 thx = thx + URAD*atan(gen.GetP(5,1)/gen.GetP(5,3));
302 thy = thy + URAD*atan(gen.GetP(5,2)/gen.GetP(5,3));
303 energy = gen.GetP(5,4);
304 positions.clear();
305 addPosition(fx,thx,fy,thy,fs);
306#endif
307 return;
308}
309
[3c40083]310void H_BeamParticle::printProperties() const {
311 cout << " M = " << getM() << "GeV ";
312 cout << " Q = " << getQ() << "e";
313 cout << " fx = " << getX() << "m ";
314 cout << " fy = " << getY() << "m ";
315 cout << " thx = " << getTX() << "rad ";
316 cout << " thy = " << getTY() << "rad ";
317 cout << endl;
318 return;
[5b822e5]319}
320
321/// The phase space vector is (x,x',y,y',E)
322/// [x] = [y] = meters
323/// [x'] = [y'] = 1 with x' = dx/ds = tan (thetaX)
324/// [E] = GeV
325const TMatrixD * H_BeamParticle::getV() const {
326 double vec[MDIM] = {fx/URAD, tan(thx/URAD), fy/URAD, tan(thy/URAD),energy};
327 TMatrixD * mat = new TMatrixD(1,MDIM,vec);
328 return mat;
329}
330
331void H_BeamParticle::printV() const {
332 TMatrixD X(*getV());
333 cout << " x = " << (X.GetMatrixArray())[0] << "m ";
334 cout << " x' = " << (X.GetMatrixArray())[1] << " ";
335 cout << " y = " << (X.GetMatrixArray())[2] << "m ";
336 cout << " y' = " << (X.GetMatrixArray())[3] << " ";
337 cout << endl;
338 return;
339}
340
341void H_BeamParticle::propagate(const double position) {
342 /// @param position is the s coordinate in m to reach
343 /// Does not propagate if position is in the middle of an otics element of the beamline.
344 vector<TVectorD>::const_iterator position_i = positions.begin();
345 double l = 0.;
346 if(position != fs) { // avoid repeating the computation if already done at this position
347 if(position == (*position_i)[INDEX_S]) {
348 fs = position;
349 fx = (*position_i)[INDEX_X];
350 fy = (*position_i)[INDEX_Y];
351 thx= (*position_i)[INDEX_TX];
352 thy= (*position_i)[INDEX_TY];
353 return;
354 } else
355 for(position_i = positions.begin(); position_i < positions.end(); position_i++) {
356 if((*position_i)[INDEX_S]>=position) {
357 if(position_i==positions.begin()) {
[3c40083]358 if(VERBOSE) cout<<"ERROR : non reachable value"<<endl;
[5b822e5]359 return;
360 }
361 l = (*position_i)[INDEX_S] - (*(position_i-1))[INDEX_S];
362 if(l==0) {
[3c40083]363 if(VERBOSE) cout<<"WARNING : no luck in choosing position, no propagation done"<<endl;
[5b822e5]364 return;
365 }
366 fs = position;
367 fx = (*(position_i-1))[INDEX_X] + (position-(*(position_i-1))[INDEX_S])*((*position_i)[INDEX_X] - (*(position_i-1))[INDEX_X])/l;
368 fy = (*(position_i-1))[INDEX_Y] + (position-(*(position_i-1))[INDEX_S])*((*position_i)[INDEX_Y] - (*(position_i-1))[INDEX_Y])/l;
369 thx = (*(position_i-1))[INDEX_TX];
370 thy = (*(position_i-1))[INDEX_TY];
371 return;
372 }
373 }
374 position_i = positions.begin();
375 cout << "Desired position is : " << position << " & positions.begin() is " << (*position_i)[INDEX_S] << endl;
[3c40083]376 cout<<"ERROR : position not reachable"<<endl;
[5b822e5]377 return;
378 }
379}
380
[3c40083]381/// Caution : do not use this method !!!
[5b822e5]382void H_BeamParticle::propagate(const H_AbstractBeamLine * beam, const H_OpticalElement * element) {
383 TMatrixD X(*getV());
[3c40083]384 X *= *(beam->getPartialMatrix(element));
[5b822e5]385 fx = URAD*(X.GetMatrixArray())[0];
386 thx = URAD*atan((X.GetMatrixArray())[1]);
387 fy = URAD*(X.GetMatrixArray())[2];
388 thy = URAD*atan((X.GetMatrixArray())[3]);
389 return;
390}
391
[3c40083]392void H_BeamParticle::propagate(const H_AbstractBeamLine * beam, const string el_name) {
393 propagate(beam->getElement(el_name)->getS());
[5b822e5]394 return;
395}
396
397void H_BeamParticle::propagate(const H_AbstractBeamLine * beam) {
398 TMatrixD X(*getV());
[3c40083]399 X *= *(beam->getBeamMatrix());
[5b822e5]400 fx = URAD*(X.GetMatrixArray())[0];
401 thx = URAD*atan((X.GetMatrixArray())[1]);
402 fy = URAD*(X.GetMatrixArray())[2];
403 thy = URAD*atan((X.GetMatrixArray())[3]);
404 return;
405}
406
407void H_BeamParticle::showPositions() const{
408 vector<TVectorD>::const_iterator position_i;
409 TVectorD temp_vec(LENGTH_VEC);
410
411 for(position_i = positions.begin(); position_i < positions.end(); position_i++) {
412 cout << "Vector (x,y,s) = (" << (*position_i)[INDEX_X] << ", " << (*position_i)[INDEX_Y] << ", " << (*position_i)[INDEX_S] << "). " << endl;
413 }
414 return ;
415}
[3c40083]416/*
[5b822e5]417TGraph * H_BeamParticle::getPath(const int x_or_y, const int color) const{
418 /// @param x_or_y = 0(1) draws the x(y) component;
419
420 const int N = (int) positions.size();
[3c40083]421 int mycolor = color;
422 if(N<2) cout<<"particle positions not calculated : please run computePath"<<endl;
[5b822e5]423 double * s = new double[N], * graph = new double[N];
424
425 int index;
426 if(x_or_y==0) {index = INDEX_X;} else {index = INDEX_Y;}
427
428 vector<TVectorD>::const_iterator position_i;
429 for(position_i = positions.begin(); position_i < positions.end(); position_i++) {
430 graph[(int)(position_i-positions.begin())] = (*position_i)[index];
431 s[(int)(position_i-positions.begin())] = (*position_i)[INDEX_S];
432 }
433
434 TGraph * ppath = new TGraph(N,s,graph);
435 ppath->SetLineColor(mycolor);
436 delete [] s;
437 delete [] graph;
438 return ppath;
439}
440
[3c40083]441TPolyLine3D * H_BeamParticle::getPath3D(const H_AbstractBeamLine * beam, const bool isfirst, const int color, const int side) const{
442 const int N = (int) positions.size();
443 int mycolor = color;
444 if(N<2) cout<<"WARNING : particle positions not calculated. Run computePath"<<endl;
445 double * s = new double[N], * graphx = new double[N], * graphy = new double[N];
446 int direction = (side<0)?-1:1;
[5b822e5]447
448
449 vector<TVectorD>::const_iterator position_i;
450 for(position_i = positions.begin(); position_i < positions.end(); position_i++) {
[3c40083]451 graphx[(int)(position_i-positions.begin())] = (*position_i)[INDEX_X];
452 graphy[(int)(position_i-positions.begin())] = (*position_i)[INDEX_Y];
453 s[(int)(position_i-positions.begin())] = (*position_i)[INDEX_S]*1000*direction;
[5b822e5]454 }
455
[3c40083]456 float coi[3] = {beam->getLength()*(-1000),-10000,-5000};
457 float cof[3] = {beam->getLength()*1000,10000,5000};
458 TView *view = TView::CreateView(11);
459 view->SetRange(coi[0],coi[1],coi[2],cof[0],cof[1],cof[2]);
460 TPolyLine3D* ppath = new TPolyLine3D(N,s,graphx,graphy);
461 ppath->SetLineColor(mycolor);
462
463 if(isfirst) {
464 ppath->Draw();
465 view->ShowAxis();
466 } else {
467 ppath->Draw("same");
468 }
469
470 delete [] s;
471 delete [] graphx;
472 delete [] graphy;
473 return ppath;
474} // getPath3D
475*/
476void H_BeamParticle::computePath(const H_AbstractBeamLine * beam) {
477 computePath(beam,true);
478}
[5b822e5]479
480// should be removed later, to keep only computePath(const H_AbstractBeamLine & , const bool)
481void H_BeamParticle::computePath(const H_AbstractBeamLine * beam, const bool NonLinear) {
482 TMatrixD temp_mat(MDIM,MDIM);
483 double temp_x, temp_y, temp_s, temp_tx, temp_ty;
484
485 temp_x = (positions.front())[INDEX_X];
486 temp_tx = (positions.front())[INDEX_TX];
487 temp_y = (positions.front())[INDEX_Y];
488 temp_ty = (positions.front())[INDEX_TY];
489 temp_s = (positions.front())[INDEX_S];
490
491 double vec[MDIM] = {temp_x/URAD, tan(temp_tx/URAD), temp_y/URAD, tan(temp_ty/URAD),energy,1};
492
493 extern bool relative_energy;
494 if(relative_energy) {
495 vec[4] = energy-BE;
496 } else {
497 vec[4] = energy;
498 }
499
500 TMatrixD mat(1,MDIM,vec);
501
502 const int N =beam->getNumberOfElements();
503 double xys[LENGTH_VEC];
504
505 double energy_loss = NonLinear?BE-energy:0;
506
507 for (int i=0; i<N; i++) {
508 const unsigned pos = i;
[3c40083]509 mat[0][0] = mat[0][0] - beam->getElement(pos)->getX();
510 mat[0][1] = mat[0][1] - tan(beam->getElement(pos)->getTX()/URAD)*URAD;
511 mat[0][2] = mat[0][2] - beam->getElement(pos)->getY();
512 mat[0][3] = mat[0][3] - tan(beam->getElement(pos)->getTY()/URAD)*URAD;
513 mat *= beam->getElement(pos)->getMatrix(energy_loss,mp,qp);
514 mat[0][0] = mat[0][0] + beam->getElement(pos)->getX();
515 mat[0][1] = mat[0][1] + tan(beam->getElement(pos)->getTX()/URAD)*URAD;
516 mat[0][2] = mat[0][2] + beam->getElement(pos)->getY();
517 mat[0][3] = mat[0][3] + tan(beam->getElement(pos)->getTY()/URAD)*URAD;
518 xys[0] = mat.GetMatrixArray()[0]*URAD;
519 xys[1] = atan(mat.GetMatrixArray()[1])*URAD;
520 xys[2] = mat.GetMatrixArray()[2]*URAD;
521 xys[3] = atan(mat.GetMatrixArray()[3])*URAD;
522 xys[4] = beam->getElement(pos)->getS()+beam->getElement(pos)->getLength();
523 addPosition(xys[0],xys[1],xys[2],xys[3],xys[4]);
524 fx = xys[0];
525 fy = xys[2];
526 thx = xys[1];
527 thy = xys[3];
528 }
[5b822e5]529}
530
531void H_BeamParticle::computePath(const H_AbstractBeamLine & beam, const bool NonLinear) {
532 TMatrixD temp_mat(MDIM,MDIM);
533 double temp_x, temp_y, temp_s, temp_tx, temp_ty;
534
535 temp_x = (positions.front())[INDEX_X];
536 temp_tx = (positions.front())[INDEX_TX];
537 temp_y = (positions.front())[INDEX_Y];
538 temp_ty = (positions.front())[INDEX_TY];
539 temp_s = (positions.front())[INDEX_S];
540
541 double vec[MDIM] = {temp_x/URAD, tan(temp_tx/URAD), temp_y/URAD, tan(temp_ty/URAD),energy,1};
542
543 extern bool relative_energy;
544 if(relative_energy) {
545 vec[4] = energy-BE;
546 } else {
547 vec[4] = energy;
548 }
549
550 TMatrixD mat(1,MDIM,vec);
551
552 const int N =beam.getNumberOfElements();
553 double xys[LENGTH_VEC];
554
555 double energy_loss = NonLinear?BE-energy:0;
556
557 for (int i=0; i<N; i++) {
558 const unsigned pos = i;
[3c40083]559 mat[0][0] = mat[0][0] - beam.getElement(pos)->getX();
560 mat[0][1] = mat[0][1] - tan(beam.getElement(pos)->getTX());
561 mat[0][2] = mat[0][2] - beam.getElement(pos)->getY();
562 mat[0][3] = mat[0][3] - tan(beam.getElement(pos)->getTY());
563 mat *= beam.getElement(pos)->getMatrix(energy_loss,mp,qp);
564 mat[0][0] = mat[0][0] + beam.getElement(pos)->getX();
565 mat[0][1] = mat[0][1] + tan(beam.getElement(pos)->getTX());
566 mat[0][2] = mat[0][2] + beam.getElement(pos)->getY();
567 mat[0][3] = mat[0][3] + tan(beam.getElement(pos)->getTY());
568 xys[0] = mat.GetMatrixArray()[0]*URAD;
569 xys[1] = atan(mat.GetMatrixArray()[1])*URAD;
570 xys[2] = mat.GetMatrixArray()[2]*URAD;
571 xys[3] = atan(mat.GetMatrixArray()[3])*URAD;
572 xys[4] = beam.getElement(pos)->getS()+beam.getElement(pos)->getLength();
573 addPosition(xys[0],xys[1],xys[2],xys[3],xys[4]);
574 fx = xys[0];
575 fy = xys[2];
576 thx = xys[1];
577 thy = xys[3];
578 }
[5b822e5]579}
580
581void H_BeamParticle::resetPath() {
582 double temp_x, temp_y, temp_s, temp_tx, temp_ty;
583
584 temp_x = (positions.front())[INDEX_X];
585 temp_tx = (positions.front())[INDEX_TX];
586 temp_y = (positions.front())[INDEX_Y];
587 temp_ty = (positions.front())[INDEX_TY];
588 temp_s = (positions.front())[INDEX_S];
589 positions.clear();
590 addPosition(temp_x,temp_tx,temp_y,temp_ty,temp_s);
591}
592
593const TVectorD * H_BeamParticle::getPosition(const int element_position) const {
594 const int N = (element_position<0)?0:(( ((unsigned int) element_position)>positions.size()-1)?positions.size()-1:element_position);
595 return &(*(positions.begin()+N));// same as "return &positions[N];", but more efficient
596}
Note: See TracBrowser for help on using the repository browser.