source: trunk/source/processes/hadronic/util/include/G4GHEKinematicsVector.hh@ 819

Last change on this file since 819 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 20.5 KB
RevLine 
[819]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: felawka@alph04.triumf.ca
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
46#include "G4ios.hh"
47
48#include <CLHEP/Units/PhysicalConstants.h>
49
50class G4GHEKinematicsVector
51 {
52 public:
53 inline
54 G4GHEKinematicsVector()
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
70 ~G4GHEKinematicsVector() {}
71
72 inline
73 G4GHEKinematicsVector( const G4GHEKinematicsVector & p )
74 {
75 momentum.setX( p.momentum.x() );
76 momentum.setY( p.momentum.y() );
77 momentum.setZ( p.momentum.z() );
78 energy = p.energy;
79 kineticEnergy = p.kineticEnergy;
80 mass = p.mass;
81 charge = p.charge;
82 timeOfFlight = p.timeOfFlight;
83 side = p.side;
84 flag = p.flag;
85 code = p.code;
86 particleDef = p.particleDef;
87 }
88
89 inline
90 G4GHEKinematicsVector & operator = ( const G4GHEKinematicsVector & p )
91 {
92 momentum.setX( p.momentum.x() );
93 momentum.setY( p.momentum.y() );
94 momentum.setZ( p.momentum.z() );
95 energy = p.energy;
96 kineticEnergy = p.kineticEnergy;
97 mass = p.mass;
98 charge = p.charge;
99 timeOfFlight = p.timeOfFlight;
100 side = p.side;
101 flag = p.flag;
102 code = p.code;
103 particleDef = p.particleDef;
104 return *this;
105 }
106
107 inline
108 void SetMomentum( G4ParticleMomentum mom ) { momentum = mom; return; };
109
110 inline
111 void SetMomentumAndUpdate( G4ParticleMomentum mom )
112 {
113 momentum = mom;
114 energy = std::sqrt(mass*mass + momentum.mag2());
115 kineticEnergy = std::max(0.,energy - mass);
116 return;
117 }
118
119 inline const
120 G4ParticleMomentum GetMomentum() const { return momentum; }
121
122 inline
123 void SetMomentum( G4double x, G4double y, G4double z)
124 {
125 momentum.setX( x );
126 momentum.setY( y );
127 momentum.setZ( z );
128 return;
129 }
130
131 inline
132 void SetMomentumAndUpdate( G4double x, G4double y, G4double z )
133 {
134 momentum.setX( x );
135 momentum.setY( y );
136 momentum.setZ( z );
137 energy = std::sqrt(mass*mass + momentum.mag2());
138 kineticEnergy = std::max(0.,energy-mass);
139 return;
140 }
141
142 inline
143 void SetMomentum( G4double x, G4double y )
144 {
145 momentum.setX( x );
146 momentum.setY( y );
147 return;
148 }
149
150 inline
151 void SetMomentumAndUpdate( G4double x, G4double y )
152 {
153 momentum.setX( x );
154 momentum.setY( y );
155 energy = std::sqrt(mass*mass + momentum.mag2());
156 kineticEnergy = std::max(0.,energy-mass);
157 return;
158 }
159
160 inline
161 void SetMomentum( G4double z )
162 {
163 momentum.setZ( z );
164 return;
165 }
166
167 inline
168 void SetMomentumAndUpdate( G4double z )
169 {
170 momentum.setZ( z );
171 energy = std::sqrt(mass*mass + momentum.mag2());
172 kineticEnergy = std::max(0.,energy-mass);
173 return;
174 }
175
176 inline
177 void SetEnergy( G4double e ) { energy = e; return; }
178
179 inline
180 void SetEnergyAndUpdate( G4double e )
181 {
182 if (e <= mass)
183 {
184 energy = mass;
185 kineticEnergy = 0.;
186 momentum.setX( 0.);
187 momentum.setY( 0.);
188 momentum.setZ( 0.);
189 }
190 else
191 {
192 energy = e;
193 kineticEnergy = energy - mass;
194 G4double momold = momentum.mag();
195 G4double momnew = std::sqrt(energy*energy - mass*mass);
196 if (momold == 0.)
197 {
198 G4double cost = 1.0- 2.0*G4UniformRand();
199 G4double sint = std::sqrt(1. - cost*cost);
200 G4double phi = twopi* G4UniformRand();
201 momentum.setX( momnew * sint * std::cos(phi));
202 momentum.setY( momnew * sint * std::sin(phi));
203 momentum.setZ( momnew * cost);
204 }
205 else
206 {
207 momnew /= momold;
208 momentum.setX(momentum.x()*momnew);
209 momentum.setY(momentum.y()*momnew);
210 momentum.setZ(momentum.z()*momnew);
211 }
212 }
213 return;
214 }
215
216 inline
217 void SetKineticEnergy( G4double ekin ) { kineticEnergy = ekin; return; }
218
219 inline
220 void SetKineticEnergyAndUpdate(G4double ekin)
221 {
222 if (ekin <= 0.)
223 {
224 energy = mass;
225 kineticEnergy = 0.;
226 momentum.setX( 0.);
227 momentum.setY( 0.);
228 momentum.setZ( 0.);
229 }
230 else
231 {
232 energy = ekin + mass;
233 kineticEnergy = ekin;
234 G4double momold = momentum.mag();
235 G4double momnew = std::sqrt(energy*energy - mass*mass);
236 if (momold == 0.)
237 {
238 G4double cost = 1.0-2.0*G4UniformRand();
239 G4double sint = std::sqrt(1. - cost*cost);
240 G4double phi = twopi* G4UniformRand();
241 momentum.setX( momnew * sint * std::cos(phi));
242 momentum.setY( momnew * sint * std::sin(phi));
243 momentum.setZ( momnew * cost);
244 }
245 else
246 {
247 momnew /= momold;
248 momentum.setX(momentum.x()*momnew);
249 momentum.setY(momentum.y()*momnew);
250 momentum.setZ(momentum.z()*momnew);
251 }
252 }
253 return;
254 }
255
256 inline
257 G4double GetEnergy() {return energy;}
258
259 inline
260 G4double GetKineticEnergy() {return kineticEnergy;}
261
262 inline
263 void SetMass( G4double m ) { mass = m; return; }
264
265 inline
266 void SetMassAndUpdate( G4double m )
267 {
268 kineticEnergy = std::max(0., energy - m);
269 mass = m;
270 energy = kineticEnergy + mass;
271 G4double momnew = std::sqrt(std::max(0., energy*energy - mass*mass));
272 if ( momnew == 0.0)
273 {
274 momentum.setX( 0.0 );
275 momentum.setY( 0.0 );
276 momentum.setZ( 0.0 );
277 }
278 else
279 {
280 G4double momold = momentum.mag();
281 if (momold == 0.)
282 {
283 G4double cost = 1.-2.*G4UniformRand();
284 G4double sint = std::sqrt(1.-cost*cost);
285 G4double phi = twopi*G4UniformRand();
286 momentum.setX( momnew*sint*std::cos(phi));
287 momentum.setY( momnew*sint*std::sin(phi));
288 momentum.setZ( momnew*cost);
289 }
290 else
291 {
292 momnew /= momold;
293 momentum.setX( momentum.x()*momnew );
294 momentum.setY( momentum.y()*momnew );
295 momentum.setZ( momentum.z()*momnew );
296 }
297 }
298 return;
299 }
300
301 inline
302 G4double GetMass() { return mass; }
303
304 inline
305 void SetCharge( G4double c ) { charge = c; return; }
306
307 inline
308 G4double GetCharge() {return charge; }
309
310 inline
311 void SetTOF( G4double t ) { timeOfFlight = t; return; }
312
313 inline
314 G4double GetTOF() { return timeOfFlight; }
315
316 inline
317 void SetSide( G4int s ) { side = s; return; }
318
319 inline
320 G4int GetSide() { return side; }
321
322 inline
323 void setFlag( G4bool f ) { flag = f; return; }
324
325 inline
326 G4bool getFlag() { return flag; }
327
328 inline
329 void SetCode( G4int c ) { code = c; return; }
330
331 inline
332 void SetParticleDef( G4ParticleDefinition * c ) { particleDef = c; return; }
333
334 inline
335 G4int GetCode() { return code; }
336
337 inline
338 G4ParticleDefinition * GetParticleDef() { return particleDef; }
339
340 inline
341 void SetZero()
342 {
343 momentum.setX( 0.0 );
344 momentum.setY( 0.0 );
345 momentum.setZ( 0.0 );
346 energy = 0.0;
347 kineticEnergy = 0.0;
348 mass = 0.0;
349 charge = 0.0;
350 timeOfFlight = 0.0;
351 side = 0;
352 flag = false;
353 code = 0;
354 particleDef = NULL;
355 }
356
357 inline
358 void Add( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2 )
359 {
360 momentum = p1.momentum + p2.momentum;
361 energy = p1.energy + p2.energy;
362 G4double b = energy*energy - momentum.mag2();
363 if( b < 0 )
364 mass = -1. * std::sqrt( -b );
365 else
366 mass = std::sqrt( b );
367 kineticEnergy = std::max(0.,energy - mass);
368 charge = p1.charge + p2.charge;
369 code = p1.code + p2.code;
370 particleDef = p1.particleDef;
371 }
372
373 inline
374 void Sub( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2 )
375 {
376 momentum = p1.momentum - p2.momentum;
377 energy = p1.energy - p2.energy;
378 G4double b = energy*energy - momentum.mag2();
379 if( b < 0 )
380 mass = -1. * std::sqrt( -b );
381 else
382 mass = std::sqrt( b );
383 kineticEnergy = std::max(0.,energy - mass);
384 charge = p1.charge - p2.charge;
385 code = p1.code - p2.code;
386 particleDef = p1.particleDef;
387 }
388
389 inline
390 void Lor( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2 )
391 {
392 G4double a;
393 a = ( p1.momentum.dot(p2.momentum)/(p2.energy+p2.mass) - p1.energy ) / p2.mass;
394 momentum.setX( p1.momentum.x()+a*p2.momentum.x() );
395 momentum.setY( p1.momentum.y()+a*p2.momentum.y() );
396 momentum.setZ( p1.momentum.z()+a*p2.momentum.z() );
397 energy = std::sqrt( sqr(p1.mass) + momentum.mag2() );
398 mass = p1.mass;
399 kineticEnergy = std::max(0.,energy - mass);
400 timeOfFlight = p1.timeOfFlight;
401 side = p1.side;
402 flag = p1.flag;
403 code = p1.code;
404 particleDef = p1.particleDef;
405 }
406
407 inline
408 G4double CosAng( const G4GHEKinematicsVector & p )
409 {
410 G4double a = std::sqrt( momentum.mag2() * p.momentum.mag2() );
411 if( a != 0.0 )
412 {
413 a = (momentum.x()*p.momentum.x() +
414 momentum.y()*p.momentum.y() +
415 momentum.z()*p.momentum.z()) / a;
416 if( std::fabs(a) > 1.0 ) a<0.0 ? a=-1.0 : a=1.0;
417 }
418 return a;
419 }
420 inline
421 G4double Ang(const G4GHEKinematicsVector & p )
422 {
423 G4double a = std::sqrt( momentum.mag2() * p.momentum.mag2() );
424 if( a != 0.0 )
425 {
426 a = (momentum.x()*p.momentum.x() +
427 momentum.y()*p.momentum.y() +
428 momentum.z()*p.momentum.z()) / a;
429 if( std::fabs(a) > 1.0 ) a<0.0 ? a=-1.0 : a=1.0;
430 }
431 return std::acos(a);
432 }
433
434 inline
435 G4double Dot4( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
436 {
437 return ( p1.energy * p2.energy
438 - p1.momentum.x() * p2.momentum.x()
439 - p1.momentum.y() * p2.momentum.y()
440 - p1.momentum.z() * p2.momentum.z() );
441 }
442
443 inline
444 G4double Impu( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
445 {
446 return ( - sqr( p1.energy - p2.energy)
447 + sqr(p1.momentum.x() - p2.momentum.x())
448 + sqr(p1.momentum.y() - p2.momentum.y())
449 + sqr(p1.momentum.z() - p2.momentum.z()) );
450 }
451
452 inline
453 void Add3( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
454 {
455 momentum.setX( p1.momentum.x() + p2.momentum.x());
456 momentum.setY( p1.momentum.y() + p2.momentum.y());
457 momentum.setZ( p1.momentum.z() + p2.momentum.z());
458 return;
459 }
460
461 inline
462 void Sub3( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
463 {
464 momentum.setX( p1.momentum.x() - p2.momentum.x());
465 momentum.setY( p1.momentum.y() - p2.momentum.y());
466 momentum.setZ( p1.momentum.z() - p2.momentum.z());
467 return;
468 }
469
470 inline
471 void Cross( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
472 {
473 G4double px, py, pz;
474 px = p1.momentum.y() * p2.momentum.z() - p1.momentum.z() * p2.momentum.y();
475 py = p1.momentum.z() * p2.momentum.x() - p1.momentum.x() * p2.momentum.z();
476 pz = p1.momentum.x() * p2.momentum.y() - p1.momentum.y() * p2.momentum.x();
477 momentum.setX( px );
478 momentum.setY( py );
479 momentum.setZ( pz );
480 return;
481 }
482
483 inline
484 G4double Dot( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
485 {
486 return ( p1.momentum.x() * p2.momentum.x()
487 + p1.momentum.y() * p2.momentum.y()
488 + p1.momentum.z() * p2.momentum.z() );
489 }
490
491 inline
492 void Smul( const G4GHEKinematicsVector & p, G4double h)
493 {
494 momentum.setX( h * p.momentum.x());
495 momentum.setY( h * p.momentum.y());
496 momentum.setZ( h * p.momentum.z());
497 return;
498 }
499
500 inline
501 void SmulAndUpdate( const G4GHEKinematicsVector & p, G4double h)
502 {
503 momentum.setX( h * p.momentum.x());
504 momentum.setY( h * p.momentum.y());
505 momentum.setZ( h * p.momentum.z());
506 mass = p.mass;
507 energy = std::sqrt(momentum.mag2() + mass*mass);
508 kineticEnergy = energy - mass;
509 charge = p.charge;
510 timeOfFlight = p.timeOfFlight;
511 side = p.side;
512 flag = p.flag;
513 code = p.code;
514 particleDef = p.particleDef;
515 return;
516 }
517
518 inline
519 void Norz( const G4GHEKinematicsVector & p )
520 {
521 G4double a = p.momentum.mag2();
522 if (a > 0.0) a = 1./std::sqrt(a);
523 momentum.setX( a * p.momentum.x() );
524 momentum.setY( a * p.momentum.y() );
525 momentum.setZ( a * p.momentum.z() );
526 mass = p.mass;
527 energy = std::sqrt(momentum.mag2() + mass*mass);
528 kineticEnergy = energy - mass;
529 charge = p.charge;
530 timeOfFlight = p.timeOfFlight;
531 side = p.side;
532 flag = p.flag;
533 code = p.code;
534 particleDef = p.particleDef;
535 return;
536 }
537
538 inline
539 G4double Length()
540 {
541 return momentum.mag() ;
542 }
543
544 inline
545 void Exch( G4GHEKinematicsVector & p1)
546 {
547 G4GHEKinematicsVector mx = *this;
548// mx.momentum.SetX( momentum.x());
549// mx.momentum.SetY( momentum.y());
550// mx.momentum.SetZ( momentum.z());
551// mx.energy = energy;
552// mx.kineticEnergy = kineticEnergy;
553// mx.mass = mass;
554// mx.charge = charge;
555// mx.timeOfFlight = timeOfFlight;
556// mx.side = side;
557// mx.flag = flag;
558// mx.code = code;
559// momentum.setX( p1.momentum.x());
560// momentum.setY( p1.momentum.y());
561// momentum.setZ( p1.momentum.z());
562// energy = p1.energy;
563// kineticEnergy = p1.kineticEnergy;
564// mass = p1.mass;
565// charge = p1.charge;
566// timeOfFlight = p1.timeOfFlight;
567// side = p1.side
568// flag = p1.flag;
569// code = p1.code;
570 *this = p1;
571 p1 = mx;
572 return;
573 }
574
575 inline
576 void Defs1( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
577 {
578 G4double pt2 = sqr(p1.momentum.x()) + sqr(p1.momentum.y());
579 if (pt2 > 0.0)
580 {
581 G4double ph, px, py, pz;
582 G4double cost = p2.momentum.z()/p2.momentum.mag();
583 G4double sint = 0.5 * ( std::sqrt(std::fabs((1.-cost)*(1.+cost)))
584 + std::sqrt(pt2)/p2.momentum.mag());
585 (p2.momentum.y() < 0.) ? ph = 1.5*pi : ph = halfpi;
586 if( p2.momentum.x() != 0.0)
587 ph = std::atan2(p2.momentum.y(),p2.momentum.x());
588 px = cost*std::cos(ph)*p1.momentum.x() - std::sin(ph)*p1.momentum.y()
589 + sint*std::cos(ph)*p1.momentum.z();
590 py = cost*std::sin(ph)*p1.momentum.x() + std::cos(ph)*p1.momentum.y()
591 + sint*std::sin(ph)*p1.momentum.z();
592 pz = - sint *p1.momentum.x()
593 + cost *p1.momentum.z();
594 momentum.setX( px );
595 momentum.setY( py );
596 momentum.setZ( pz );
597 }
598 else
599 {
600 momentum = p1.momentum;
601 }
602 }
603
604 inline
605 void Defs( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2,
606 G4GHEKinematicsVector & my, G4GHEKinematicsVector & mz )
607 {
608 my = p1;
609 mz = p2;
610 momentum.setX( my.momentum.y()*mz.momentum.z()
611 - my.momentum.z()*mz.momentum.y());
612 momentum.setY( my.momentum.z()*mz.momentum.x()
613 - my.momentum.x()*mz.momentum.z());
614 momentum.setZ( my.momentum.x()*mz.momentum.y()
615 - my.momentum.y()*mz.momentum.x());
616 my.momentum.setX( mz.momentum.y()*momentum.z()
617 - mz.momentum.z()*momentum.y());
618 my.momentum.setY( mz.momentum.z()*momentum.x()
619 - mz.momentum.x()*momentum.z());
620 my.momentum.setZ( mz.momentum.x()*momentum.y()
621 - mz.momentum.y()*momentum.x());
622 G4double pp;
623 pp = momentum.mag();
624 if (pp > 0.)
625 {
626 pp = 1./pp;
627 momentum.setX( momentum.x()*pp );
628 momentum.setY( momentum.y()*pp );
629 momentum.setZ( momentum.z()*pp );
630 }
631 pp = my.momentum.mag();
632 if (pp > 0.)
633 {
634 pp = 1./pp;
635 my.momentum.setX( my.momentum.x()*pp );
636 my.momentum.setY( my.momentum.y()*pp );
637 my.momentum.setZ( my.momentum.z()*pp );
638 }
639 pp = mz.momentum.mag();
640 if (pp > 0.)
641 {
642 pp = 1./pp;
643 mz.momentum.setX( mz.momentum.x()*pp );
644 mz.momentum.setY( mz.momentum.y()*pp );
645 mz.momentum.setZ( mz.momentum.z()*pp );
646 }
647 return;
648 }
649
650 inline
651 void Trac( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & mx,
652 const G4GHEKinematicsVector & my, const G4GHEKinematicsVector & mz)
653 {
654 double px, py, pz;
655 px = mx.momentum.x()*p1.momentum.x()
656 + mx.momentum.y()*p1.momentum.y()
657 + mx.momentum.z()*p1.momentum.z();
658 py = my.momentum.x()*p1.momentum.x()
659 + my.momentum.y()*p1.momentum.y()
660 + my.momentum.z()*p1.momentum.z();
661 pz = mz.momentum.x()*p1.momentum.x()
662 + mz.momentum.y()*p1.momentum.y()
663 + mz.momentum.z()*p1.momentum.z();
664 momentum.setX( px );
665 momentum.setY( py );
666 momentum.setZ( pz );
667 return;
668 }
669
670 inline
671 void Print( G4int L)
672 {
673 G4cout << "G4GHEKinematicsVector: "
674 << L << " " << momentum.x() << " " << momentum.y() << " " << momentum.z() << " "
675 << energy << " " << kineticEnergy << " " << mass << " " << charge << " "
676 << timeOfFlight << " " << side << " " << flag << " " << code << particleDef << G4endl;
677 return;
678 }
679
680 G4ParticleMomentum momentum;
681 G4double energy;
682 G4double kineticEnergy;
683 G4double mass;
684 G4double charge;
685 G4double timeOfFlight;
686 G4int side;
687 G4bool flag;
688 G4int code;
689 G4ParticleDefinition * particleDef;
690};
691
692#endif
Note: See TracBrowser for help on using the repository browser.