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 |
|
---|
39 | using namespace std;
|
---|
40 |
|
---|
41 | void 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 |
|
---|
59 | H_BeamParticle::H_BeamParticle() {
|
---|
60 | init();
|
---|
61 | }
|
---|
62 |
|
---|
63 | H_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 |
|
---|
80 | H_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 |
|
---|
87 | H_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 |
|
---|
106 | 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 |
|
---|
148 | void 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 | void 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 |
|
---|
170 | void 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 |
|
---|
179 | void 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 |
|
---|
188 | void 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 |
|
---|
197 | void H_BeamParticle::smearE(const double erre) {
|
---|
198 | energy = gRandom->Gaus(energy,erre);
|
---|
199 | return;
|
---|
200 | }
|
---|
201 |
|
---|
202 | void H_BeamParticle::smearE() {
|
---|
203 | energy = gRandom->Gaus(energy,SBE);
|
---|
204 | return;
|
---|
205 | }
|
---|
206 |
|
---|
207 | void 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 |
|
---|
214 | void H_BeamParticle::smearS() {
|
---|
215 | fs = gRandom->Gaus(fs,SS);
|
---|
216 | positions.clear();
|
---|
217 | addPosition(fx,thx,fy,thy,fs);
|
---|
218 | return;
|
---|
219 | }
|
---|
220 |
|
---|
221 | void 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 |
|
---|
238 | void 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 |
|
---|
245 | void H_BeamParticle::setE(const double ene) {
|
---|
246 | energy = ene;
|
---|
247 | return;
|
---|
248 | }
|
---|
249 |
|
---|
250 | void 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 |
|
---|
266 | const 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 |
|
---|
271 | void H_BeamParticle::emitGamma(const double gee, const double gq2) {
|
---|
272 | emitGamma(gee,gq2,0,2*PI);
|
---|
273 | return;
|
---|
274 | }
|
---|
275 |
|
---|
276 | void 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 |
|
---|
331 | void 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 |
|
---|
364 | void 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
|
---|
379 | const 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 |
|
---|
385 | void 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 |
|
---|
395 | void 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) !!!
|
---|
436 | void 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 |
|
---|
446 | void 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 |
|
---|
451 | void 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 |
|
---|
461 | void 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 |
|
---|
471 | TGraph * 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 |
|
---|
495 | void 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)
|
---|
500 | void 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
|
---|
575 | void 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 |
|
---|
646 | void 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 |
|
---|
658 | const 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 | }
|
---|