Fork me on GitHub

source: svn/trunk/Utilities/Hector/src/H_BeamParticle.cc@ 185

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

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