Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4GHEKinematicsVector.hh
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27//
28// ------------------------------------------------------------
29// GEANT 4 class header file --- Copyright CERN 1998
30// CERN Geneva Switzerland
31//
32// History: first implementation, based on object model of
33// 2nd December 1995, G.Cosmo
34// ------------ G4GHEKinematicsVector utility class ------
35// by Larry Felawka (TRIUMF), March 1997
36// E-mail: [email protected]
37// ************************************************************
38//-----------------------------------------------------------------------------
39
40// Store, Retrieve and manipulate particle data.
41// Based on "G4GHEVector" class of H. Fesefeldt.
42
43#ifndef G4GHEKinematicsVector_h
44#define G4GHEKinematicsVector_h 1
45
47#include "G4ios.hh"
48#include "G4ParticleMomentum.hh"
49
51 {
52 public:
53 inline
55 {
56 momentum.setX( 0.0 );
57 momentum.setY( 0.0 );
58 momentum.setZ( 0.0 );
59 energy = 0.0;
60 kineticEnergy = 0.0;
61 mass = 0.0;
62 charge = 0.0;
63 timeOfFlight = 0.0;
64 side = 0;
65 flag = false;
66 code = 0;
67 particleDef = NULL;
68 }
69
71
72 inline
74 {
75 momentum.setX( p.momentum.x() );
76 momentum.setY( p.momentum.y() );
77 momentum.setZ( p.momentum.z() );
78 energy = p.energy;
80 mass = p.mass;
81 charge = p.charge;
83 side = p.side;
84 flag = p.flag;
85 code = p.code;
87 }
88
89 inline
91 {
92 if (this != &p)
93 {
94 momentum.setX( p.momentum.x() );
95 momentum.setY( p.momentum.y() );
96 momentum.setZ( p.momentum.z() );
97 energy = p.energy;
99 mass = p.mass;
100 charge = p.charge;
102 side = p.side;
103 flag = p.flag;
104 code = p.code;
106 }
107 return *this;
108 }
109
110 inline
111 void SetMomentum( G4ParticleMomentum mom ) { momentum = mom; return; };
112
113 inline
115 {
116 momentum = mom;
117 energy = std::sqrt(mass*mass + momentum.mag2());
118 kineticEnergy = std::max(0.,energy - mass);
119 return;
120 }
121
122 inline const
124
125 inline
127 {
128 momentum.setX( x );
129 momentum.setY( y );
130 momentum.setZ( z );
131 return;
132 }
133
134 inline
136 {
137 momentum.setX( x );
138 momentum.setY( y );
139 momentum.setZ( z );
140 energy = std::sqrt(mass*mass + momentum.mag2());
141 kineticEnergy = std::max(0.,energy-mass);
142 return;
143 }
144
145 inline
147 {
148 momentum.setX( x );
149 momentum.setY( y );
150 return;
151 }
152
153 inline
155 {
156 momentum.setX( x );
157 momentum.setY( y );
158 energy = std::sqrt(mass*mass + momentum.mag2());
159 kineticEnergy = std::max(0.,energy-mass);
160 return;
161 }
162
163 inline
165 {
166 momentum.setZ( z );
167 return;
168 }
169
170 inline
172 {
173 momentum.setZ( z );
174 energy = std::sqrt(mass*mass + momentum.mag2());
175 kineticEnergy = std::max(0.,energy-mass);
176 return;
177 }
178
179 inline
180 void SetEnergy( G4double e ) { energy = e; return; }
181
182 inline
184 {
185 if (e <= mass)
186 {
187 energy = mass;
188 kineticEnergy = 0.;
189 momentum.setX( 0.);
190 momentum.setY( 0.);
191 momentum.setZ( 0.);
192 }
193 else
194 {
195 energy = e;
197 G4double momold = momentum.mag();
198 G4double momnew = std::sqrt(energy*energy - mass*mass);
199 if (momold == 0.)
200 {
201 G4double cost = 1.0- 2.0*G4UniformRand();
202 G4double sint = std::sqrt(1. - cost*cost);
203 G4double phi = CLHEP::twopi* G4UniformRand();
204 momentum.setX( momnew * sint * std::cos(phi));
205 momentum.setY( momnew * sint * std::sin(phi));
206 momentum.setZ( momnew * cost);
207 }
208 else
209 {
210 momnew /= momold;
211 momentum.setX(momentum.x()*momnew);
212 momentum.setY(momentum.y()*momnew);
213 momentum.setZ(momentum.z()*momnew);
214 }
215 }
216 return;
217 }
218
219 inline
220 void SetKineticEnergy( G4double ekin ) { kineticEnergy = ekin; return; }
221
222 inline
224 {
225 if (ekin <= 0.)
226 {
227 energy = mass;
228 kineticEnergy = 0.;
229 momentum.setX( 0.);
230 momentum.setY( 0.);
231 momentum.setZ( 0.);
232 }
233 else
234 {
235 energy = ekin + mass;
236 kineticEnergy = ekin;
237 G4double momold = momentum.mag();
238 G4double momnew = std::sqrt(energy*energy - mass*mass);
239 if (momold == 0.)
240 {
241 G4double cost = 1.0-2.0*G4UniformRand();
242 G4double sint = std::sqrt(1. - cost*cost);
243 G4double phi = CLHEP::twopi* G4UniformRand();
244 momentum.setX( momnew * sint * std::cos(phi));
245 momentum.setY( momnew * sint * std::sin(phi));
246 momentum.setZ( momnew * cost);
247 }
248 else
249 {
250 momnew /= momold;
251 momentum.setX(momentum.x()*momnew);
252 momentum.setY(momentum.y()*momnew);
253 momentum.setZ(momentum.z()*momnew);
254 }
255 }
256 return;
257 }
258
259 inline
261
262 inline
264
265 inline
266 void SetMass( G4double mas ) { mass = mas; return; }
267
268 inline
270 {
271 kineticEnergy = std::max(0., energy - mas);
272 mass = mas;
274 G4double momnew = std::sqrt(std::max(0., energy*energy - mass*mass));
275 if ( momnew == 0.0)
276 {
277 momentum.setX( 0.0 );
278 momentum.setY( 0.0 );
279 momentum.setZ( 0.0 );
280 }
281 else
282 {
283 G4double momold = momentum.mag();
284 if (momold == 0.)
285 {
286 G4double cost = 1.-2.*G4UniformRand();
287 G4double sint = std::sqrt(1.-cost*cost);
288 G4double phi = CLHEP::twopi*G4UniformRand();
289 momentum.setX( momnew*sint*std::cos(phi));
290 momentum.setY( momnew*sint*std::sin(phi));
291 momentum.setZ( momnew*cost);
292 }
293 else
294 {
295 momnew /= momold;
296 momentum.setX( momentum.x()*momnew );
297 momentum.setY( momentum.y()*momnew );
298 momentum.setZ( momentum.z()*momnew );
299 }
300 }
301 return;
302 }
303
304 inline
305 G4double GetMass() { return mass; }
306
307 inline
308 void SetCharge( G4double c ) { charge = c; return; }
309
310 inline
312
313 inline
314 void SetTOF( G4double t ) { timeOfFlight = t; return; }
315
316 inline
318
319 inline
320 void SetSide( G4int sid ) { side = sid; return; }
321
322 inline
323 G4int GetSide() { return side; }
324
325 inline
326 void setFlag( G4bool f ) { flag = f; return; }
327
328 inline
329 G4bool getFlag() { return flag; }
330
331 inline
332 void SetCode( G4int c ) { code = c; return; }
333
334 inline
336
337 inline
338 G4int GetCode() { return code; }
339
340 inline
342
343 inline
344 void SetZero()
345 {
346 momentum.setX( 0.0 );
347 momentum.setY( 0.0 );
348 momentum.setZ( 0.0 );
349 energy = 0.0;
350 kineticEnergy = 0.0;
351 mass = 0.0;
352 charge = 0.0;
353 timeOfFlight = 0.0;
354 side = 0;
355 flag = false;
356 code = 0;
357 particleDef = NULL;
358 }
359
360 inline
361 void Add( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2 )
362 {
363 momentum = p1.momentum + p2.momentum;
364 energy = p1.energy + p2.energy;
366 if( b < 0 )
367 mass = -1. * std::sqrt( -b );
368 else
369 mass = std::sqrt( b );
370 kineticEnergy = std::max(0.,energy - mass);
371 charge = p1.charge + p2.charge;
372 code = p1.code + p2.code;
374 }
375
376 inline
377 void Sub( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2 )
378 {
379 momentum = p1.momentum - p2.momentum;
380 energy = p1.energy - p2.energy;
382 if( b < 0 )
383 mass = -1. * std::sqrt( -b );
384 else
385 mass = std::sqrt( b );
386 kineticEnergy = std::max(0.,energy - mass);
387 charge = p1.charge - p2.charge;
388 code = p1.code - p2.code;
390 }
391
392 inline
393 void Lor( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2 )
394 {
395 G4double a;
396 a = ( p1.momentum.dot(p2.momentum)/(p2.energy+p2.mass) - p1.energy ) / p2.mass;
397 momentum.setX( p1.momentum.x()+a*p2.momentum.x() );
398 momentum.setY( p1.momentum.y()+a*p2.momentum.y() );
399 momentum.setZ( p1.momentum.z()+a*p2.momentum.z() );
400 energy = std::sqrt( sqr(p1.mass) + momentum.mag2() );
401 mass = p1.mass;
402 kineticEnergy = std::max(0.,energy - mass);
404 side = p1.side;
405 flag = p1.flag;
406 code = p1.code;
408 }
409
410 inline
412 {
413 G4double a = std::sqrt( momentum.mag2() * p.momentum.mag2() );
414 if( a != 0.0 )
415 {
416 a = (momentum.x()*p.momentum.x() +
417 momentum.y()*p.momentum.y() +
418 momentum.z()*p.momentum.z()) / a;
419 if( std::abs(a) > 1.0 ) a<0.0 ? a=-1.0 : a=1.0;
420 }
421 return a;
422 }
423 inline
425 {
426 G4double a = std::sqrt( momentum.mag2() * p.momentum.mag2() );
427 if( a != 0.0 )
428 {
429 a = (momentum.x()*p.momentum.x() +
430 momentum.y()*p.momentum.y() +
431 momentum.z()*p.momentum.z()) / a;
432 if( std::abs(a) > 1.0 ) a<0.0 ? a=-1.0 : a=1.0;
433 }
434 return std::acos(a);
435 }
436
437 inline
439 {
440 return ( p1.energy * p2.energy
441 - p1.momentum.x() * p2.momentum.x()
442 - p1.momentum.y() * p2.momentum.y()
443 - p1.momentum.z() * p2.momentum.z() );
444 }
445
446 inline
448 {
449 return ( - sqr( p1.energy - p2.energy)
450 + sqr(p1.momentum.x() - p2.momentum.x())
451 + sqr(p1.momentum.y() - p2.momentum.y())
452 + sqr(p1.momentum.z() - p2.momentum.z()) );
453 }
454
455 inline
457 {
458 momentum.setX( p1.momentum.x() + p2.momentum.x());
459 momentum.setY( p1.momentum.y() + p2.momentum.y());
460 momentum.setZ( p1.momentum.z() + p2.momentum.z());
461 return;
462 }
463
464 inline
466 {
467 momentum.setX( p1.momentum.x() - p2.momentum.x());
468 momentum.setY( p1.momentum.y() - p2.momentum.y());
469 momentum.setZ( p1.momentum.z() - p2.momentum.z());
470 return;
471 }
472
473 inline
475 {
476 G4double px, py, pz;
477 px = p1.momentum.y() * p2.momentum.z() - p1.momentum.z() * p2.momentum.y();
478 py = p1.momentum.z() * p2.momentum.x() - p1.momentum.x() * p2.momentum.z();
479 pz = p1.momentum.x() * p2.momentum.y() - p1.momentum.y() * p2.momentum.x();
480 momentum.setX( px );
481 momentum.setY( py );
482 momentum.setZ( pz );
483 return;
484 }
485
486 inline
488 {
489 return ( p1.momentum.x() * p2.momentum.x()
490 + p1.momentum.y() * p2.momentum.y()
491 + p1.momentum.z() * p2.momentum.z() );
492 }
493
494 inline
496 {
497 momentum.setX( h * p.momentum.x());
498 momentum.setY( h * p.momentum.y());
499 momentum.setZ( h * p.momentum.z());
500 return;
501 }
502
503 inline
505 {
506 momentum.setX( h * p.momentum.x());
507 momentum.setY( h * p.momentum.y());
508 momentum.setZ( h * p.momentum.z());
509 mass = p.mass;
510 energy = std::sqrt(momentum.mag2() + mass*mass);
512 charge = p.charge;
514 side = p.side;
515 flag = p.flag;
516 code = p.code;
518 return;
519 }
520
521 inline
522 void Norz( const G4GHEKinematicsVector & p )
523 {
524 G4double a = p.momentum.mag2();
525 if (a > 0.0) a = 1./std::sqrt(a);
526 momentum.setX( a * p.momentum.x() );
527 momentum.setY( a * p.momentum.y() );
528 momentum.setZ( a * p.momentum.z() );
529 mass = p.mass;
530 energy = std::sqrt(momentum.mag2() + mass*mass);
532 charge = p.charge;
534 side = p.side;
535 flag = p.flag;
536 code = p.code;
538 return;
539 }
540
541 inline
543 {
544 return momentum.mag() ;
545 }
546
547 inline
549 {
550 G4GHEKinematicsVector mx = *this;
551// mx.momentum.SetX( momentum.x());
552// mx.momentum.SetY( momentum.y());
553// mx.momentum.SetZ( momentum.z());
554// mx.energy = energy;
555// mx.kineticEnergy = kineticEnergy;
556// mx.mass = mass;
557// mx.charge = charge;
558// mx.timeOfFlight = timeOfFlight;
559// mx.side = side;
560// mx.flag = flag;
561// mx.code = code;
562// momentum.setX( p1.momentum.x());
563// momentum.setY( p1.momentum.y());
564// momentum.setZ( p1.momentum.z());
565// energy = p1.energy;
566// kineticEnergy = p1.kineticEnergy;
567// mass = p1.mass;
568// charge = p1.charge;
569// timeOfFlight = p1.timeOfFlight;
570// side = p1.side
571// flag = p1.flag;
572// code = p1.code;
573 *this = p1;
574 p1 = mx;
575 return;
576 }
577
578 inline
580 {
581 G4double pt2 = sqr(p1.momentum.x()) + sqr(p1.momentum.y());
582 if (pt2 > 0.0)
583 {
584 G4double ph, px, py, pz;
585 G4double cost = p2.momentum.z()/p2.momentum.mag();
586 G4double sint = 0.5 * ( std::sqrt(std::abs((1.-cost)*(1.+cost)))
587 + std::sqrt(pt2)/p2.momentum.mag());
588 (p2.momentum.y() < 0.) ? ph = 1.5*CLHEP::pi : ph = CLHEP::halfpi;
589 if( p2.momentum.x() != 0.0)
590 ph = std::atan2(p2.momentum.y(),p2.momentum.x());
591 px = cost*std::cos(ph)*p1.momentum.x() - std::sin(ph)*p1.momentum.y()
592 + sint*std::cos(ph)*p1.momentum.z();
593 py = cost*std::sin(ph)*p1.momentum.x() + std::cos(ph)*p1.momentum.y()
594 + sint*std::sin(ph)*p1.momentum.z();
595 pz = - sint *p1.momentum.x()
596 + cost *p1.momentum.z();
597 momentum.setX( px );
598 momentum.setY( py );
599 momentum.setZ( pz );
600 }
601 else
602 {
603 momentum = p1.momentum;
604 }
605 }
606
607 inline
610 {
611 my = p1;
612 mz = p2;
613 momentum.setX( my.momentum.y()*mz.momentum.z()
614 - my.momentum.z()*mz.momentum.y());
615 momentum.setY( my.momentum.z()*mz.momentum.x()
616 - my.momentum.x()*mz.momentum.z());
617 momentum.setZ( my.momentum.x()*mz.momentum.y()
618 - my.momentum.y()*mz.momentum.x());
619 my.momentum.setX( mz.momentum.y()*momentum.z()
620 - mz.momentum.z()*momentum.y());
621 my.momentum.setY( mz.momentum.z()*momentum.x()
622 - mz.momentum.x()*momentum.z());
623 my.momentum.setZ( mz.momentum.x()*momentum.y()
624 - mz.momentum.y()*momentum.x());
625 G4double pp;
626 pp = momentum.mag();
627 if (pp > 0.)
628 {
629 pp = 1./pp;
630 momentum.setX( momentum.x()*pp );
631 momentum.setY( momentum.y()*pp );
632 momentum.setZ( momentum.z()*pp );
633 }
634 pp = my.momentum.mag();
635 if (pp > 0.)
636 {
637 pp = 1./pp;
638 my.momentum.setX( my.momentum.x()*pp );
639 my.momentum.setY( my.momentum.y()*pp );
640 my.momentum.setZ( my.momentum.z()*pp );
641 }
642 pp = mz.momentum.mag();
643 if (pp > 0.)
644 {
645 pp = 1./pp;
646 mz.momentum.setX( mz.momentum.x()*pp );
647 mz.momentum.setY( mz.momentum.y()*pp );
648 mz.momentum.setZ( mz.momentum.z()*pp );
649 }
650 return;
651 }
652
653 inline
655 const G4GHEKinematicsVector & my, const G4GHEKinematicsVector & mz)
656 {
657 double px, py, pz;
658 px = mx.momentum.x()*p1.momentum.x()
659 + mx.momentum.y()*p1.momentum.y()
660 + mx.momentum.z()*p1.momentum.z();
661 py = my.momentum.x()*p1.momentum.x()
662 + my.momentum.y()*p1.momentum.y()
663 + my.momentum.z()*p1.momentum.z();
664 pz = mz.momentum.x()*p1.momentum.x()
665 + mz.momentum.y()*p1.momentum.y()
666 + mz.momentum.z()*p1.momentum.z();
667 momentum.setX( px );
668 momentum.setY( py );
669 momentum.setZ( pz );
670 return;
671 }
672
673 inline
674 void Print( G4int LLL)
675 {
676 G4cout << "G4GHEKinematicsVector: "
677 << LLL << " " << momentum.x() << " " << momentum.y() << " " << momentum.z() << " "
678 << energy << " " << kineticEnergy << " " << mass << " " << charge << " "
679 << timeOfFlight << " " << side << " " << flag << " " << code << particleDef << G4endl;
680 return;
681 }
682
693};
694
695#endif
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:52
double z() const
double x() const
void setY(double)
double mag2() const
double y() const
double dot(const Hep3Vector &) const
void setZ(double)
double mag() const
void setX(double)
void SetMomentumAndUpdate(G4double z)
void SetMomentumAndUpdate(G4double x, G4double y)
void SetEnergyAndUpdate(G4double e)
void SetMomentumAndUpdate(G4double x, G4double y, G4double z)
void Defs(const G4GHEKinematicsVector &p1, const G4GHEKinematicsVector &p2, G4GHEKinematicsVector &my, G4GHEKinematicsVector &mz)
G4GHEKinematicsVector(const G4GHEKinematicsVector &p)
G4ParticleDefinition * GetParticleDef()
const G4ParticleMomentum GetMomentum() const
void Trac(const G4GHEKinematicsVector &p1, const G4GHEKinematicsVector &mx, const G4GHEKinematicsVector &my, const G4GHEKinematicsVector &mz)
void SmulAndUpdate(const G4GHEKinematicsVector &p, G4double h)
void SetMassAndUpdate(G4double mas)
void Sub3(const G4GHEKinematicsVector &p1, const G4GHEKinematicsVector &p2)
void Exch(G4GHEKinematicsVector &p1)
void SetParticleDef(G4ParticleDefinition *c)
void SetMomentum(G4double x, G4double y)
G4double Ang(const G4GHEKinematicsVector &p)
void SetMomentum(G4ParticleMomentum mom)
void SetMomentumAndUpdate(G4ParticleMomentum mom)
void Sub(const G4GHEKinematicsVector &p1, const G4GHEKinematicsVector &p2)
G4ParticleDefinition * particleDef
void Lor(const G4GHEKinematicsVector &p1, const G4GHEKinematicsVector &p2)
void SetMomentum(G4double x, G4double y, G4double z)
void Cross(const G4GHEKinematicsVector &p1, const G4GHEKinematicsVector &p2)
G4double Dot(const G4GHEKinematicsVector &p1, const G4GHEKinematicsVector &p2)
G4double Impu(const G4GHEKinematicsVector &p1, const G4GHEKinematicsVector &p2)
void Norz(const G4GHEKinematicsVector &p)
G4double CosAng(const G4GHEKinematicsVector &p)
G4GHEKinematicsVector & operator=(const G4GHEKinematicsVector &p)
void Add3(const G4GHEKinematicsVector &p1, const G4GHEKinematicsVector &p2)
void SetKineticEnergyAndUpdate(G4double ekin)
void Defs1(const G4GHEKinematicsVector &p1, const G4GHEKinematicsVector &p2)
G4double Dot4(const G4GHEKinematicsVector &p1, const G4GHEKinematicsVector &p2)
void SetKineticEnergy(G4double ekin)
void Smul(const G4GHEKinematicsVector &p, G4double h)
void Add(const G4GHEKinematicsVector &p1, const G4GHEKinematicsVector &p2)
Definition: inftrees.h:24
T sqr(const T &x)
Definition: templates.hh:128