Fork me on GitHub

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

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

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