00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043 #ifndef G4GHEKinematicsVector_h
00044 #define G4GHEKinematicsVector_h 1
00045
00046 #include "G4ios.hh"
00047
00048 #include <CLHEP/Units/PhysicalConstants.h>
00049
00050 class G4GHEKinematicsVector
00051 {
00052 public:
00053 inline
00054 G4GHEKinematicsVector()
00055 {
00056 momentum.setX( 0.0 );
00057 momentum.setY( 0.0 );
00058 momentum.setZ( 0.0 );
00059 energy = 0.0;
00060 kineticEnergy = 0.0;
00061 mass = 0.0;
00062 charge = 0.0;
00063 timeOfFlight = 0.0;
00064 side = 0;
00065 flag = false;
00066 code = 0;
00067 particleDef = NULL;
00068 }
00069
00070 ~G4GHEKinematicsVector() {}
00071
00072 inline
00073 G4GHEKinematicsVector( const G4GHEKinematicsVector & p )
00074 {
00075 momentum.setX( p.momentum.x() );
00076 momentum.setY( p.momentum.y() );
00077 momentum.setZ( p.momentum.z() );
00078 energy = p.energy;
00079 kineticEnergy = p.kineticEnergy;
00080 mass = p.mass;
00081 charge = p.charge;
00082 timeOfFlight = p.timeOfFlight;
00083 side = p.side;
00084 flag = p.flag;
00085 code = p.code;
00086 particleDef = p.particleDef;
00087 }
00088
00089 inline
00090 G4GHEKinematicsVector & operator = ( const G4GHEKinematicsVector & p )
00091 {
00092 if (this != &p)
00093 {
00094 momentum.setX( p.momentum.x() );
00095 momentum.setY( p.momentum.y() );
00096 momentum.setZ( p.momentum.z() );
00097 energy = p.energy;
00098 kineticEnergy = p.kineticEnergy;
00099 mass = p.mass;
00100 charge = p.charge;
00101 timeOfFlight = p.timeOfFlight;
00102 side = p.side;
00103 flag = p.flag;
00104 code = p.code;
00105 particleDef = p.particleDef;
00106 }
00107 return *this;
00108 }
00109
00110 inline
00111 void SetMomentum( G4ParticleMomentum mom ) { momentum = mom; return; };
00112
00113 inline
00114 void SetMomentumAndUpdate( G4ParticleMomentum mom )
00115 {
00116 momentum = mom;
00117 energy = std::sqrt(mass*mass + momentum.mag2());
00118 kineticEnergy = std::max(0.,energy - mass);
00119 return;
00120 }
00121
00122 inline const
00123 G4ParticleMomentum GetMomentum() const { return momentum; }
00124
00125 inline
00126 void SetMomentum( G4double x, G4double y, G4double z)
00127 {
00128 momentum.setX( x );
00129 momentum.setY( y );
00130 momentum.setZ( z );
00131 return;
00132 }
00133
00134 inline
00135 void SetMomentumAndUpdate( G4double x, G4double y, G4double z )
00136 {
00137 momentum.setX( x );
00138 momentum.setY( y );
00139 momentum.setZ( z );
00140 energy = std::sqrt(mass*mass + momentum.mag2());
00141 kineticEnergy = std::max(0.,energy-mass);
00142 return;
00143 }
00144
00145 inline
00146 void SetMomentum( G4double x, G4double y )
00147 {
00148 momentum.setX( x );
00149 momentum.setY( y );
00150 return;
00151 }
00152
00153 inline
00154 void SetMomentumAndUpdate( G4double x, G4double y )
00155 {
00156 momentum.setX( x );
00157 momentum.setY( y );
00158 energy = std::sqrt(mass*mass + momentum.mag2());
00159 kineticEnergy = std::max(0.,energy-mass);
00160 return;
00161 }
00162
00163 inline
00164 void SetMomentum( G4double z )
00165 {
00166 momentum.setZ( z );
00167 return;
00168 }
00169
00170 inline
00171 void SetMomentumAndUpdate( G4double z )
00172 {
00173 momentum.setZ( z );
00174 energy = std::sqrt(mass*mass + momentum.mag2());
00175 kineticEnergy = std::max(0.,energy-mass);
00176 return;
00177 }
00178
00179 inline
00180 void SetEnergy( G4double e ) { energy = e; return; }
00181
00182 inline
00183 void SetEnergyAndUpdate( G4double e )
00184 {
00185 if (e <= mass)
00186 {
00187 energy = mass;
00188 kineticEnergy = 0.;
00189 momentum.setX( 0.);
00190 momentum.setY( 0.);
00191 momentum.setZ( 0.);
00192 }
00193 else
00194 {
00195 energy = e;
00196 kineticEnergy = energy - mass;
00197 G4double momold = momentum.mag();
00198 G4double momnew = std::sqrt(energy*energy - mass*mass);
00199 if (momold == 0.)
00200 {
00201 G4double cost = 1.0- 2.0*G4UniformRand();
00202 G4double sint = std::sqrt(1. - cost*cost);
00203 G4double phi = CLHEP::twopi* G4UniformRand();
00204 momentum.setX( momnew * sint * std::cos(phi));
00205 momentum.setY( momnew * sint * std::sin(phi));
00206 momentum.setZ( momnew * cost);
00207 }
00208 else
00209 {
00210 momnew /= momold;
00211 momentum.setX(momentum.x()*momnew);
00212 momentum.setY(momentum.y()*momnew);
00213 momentum.setZ(momentum.z()*momnew);
00214 }
00215 }
00216 return;
00217 }
00218
00219 inline
00220 void SetKineticEnergy( G4double ekin ) { kineticEnergy = ekin; return; }
00221
00222 inline
00223 void SetKineticEnergyAndUpdate(G4double ekin)
00224 {
00225 if (ekin <= 0.)
00226 {
00227 energy = mass;
00228 kineticEnergy = 0.;
00229 momentum.setX( 0.);
00230 momentum.setY( 0.);
00231 momentum.setZ( 0.);
00232 }
00233 else
00234 {
00235 energy = ekin + mass;
00236 kineticEnergy = ekin;
00237 G4double momold = momentum.mag();
00238 G4double momnew = std::sqrt(energy*energy - mass*mass);
00239 if (momold == 0.)
00240 {
00241 G4double cost = 1.0-2.0*G4UniformRand();
00242 G4double sint = std::sqrt(1. - cost*cost);
00243 G4double phi = CLHEP::twopi* G4UniformRand();
00244 momentum.setX( momnew * sint * std::cos(phi));
00245 momentum.setY( momnew * sint * std::sin(phi));
00246 momentum.setZ( momnew * cost);
00247 }
00248 else
00249 {
00250 momnew /= momold;
00251 momentum.setX(momentum.x()*momnew);
00252 momentum.setY(momentum.y()*momnew);
00253 momentum.setZ(momentum.z()*momnew);
00254 }
00255 }
00256 return;
00257 }
00258
00259 inline
00260 G4double GetEnergy() {return energy;}
00261
00262 inline
00263 G4double GetKineticEnergy() {return kineticEnergy;}
00264
00265 inline
00266 void SetMass( G4double mas ) { mass = mas; return; }
00267
00268 inline
00269 void SetMassAndUpdate( G4double mas )
00270 {
00271 kineticEnergy = std::max(0., energy - mas);
00272 mass = mas;
00273 energy = kineticEnergy + mass;
00274 G4double momnew = std::sqrt(std::max(0., energy*energy - mass*mass));
00275 if ( momnew == 0.0)
00276 {
00277 momentum.setX( 0.0 );
00278 momentum.setY( 0.0 );
00279 momentum.setZ( 0.0 );
00280 }
00281 else
00282 {
00283 G4double momold = momentum.mag();
00284 if (momold == 0.)
00285 {
00286 G4double cost = 1.-2.*G4UniformRand();
00287 G4double sint = std::sqrt(1.-cost*cost);
00288 G4double phi = CLHEP::twopi*G4UniformRand();
00289 momentum.setX( momnew*sint*std::cos(phi));
00290 momentum.setY( momnew*sint*std::sin(phi));
00291 momentum.setZ( momnew*cost);
00292 }
00293 else
00294 {
00295 momnew /= momold;
00296 momentum.setX( momentum.x()*momnew );
00297 momentum.setY( momentum.y()*momnew );
00298 momentum.setZ( momentum.z()*momnew );
00299 }
00300 }
00301 return;
00302 }
00303
00304 inline
00305 G4double GetMass() { return mass; }
00306
00307 inline
00308 void SetCharge( G4double c ) { charge = c; return; }
00309
00310 inline
00311 G4double GetCharge() {return charge; }
00312
00313 inline
00314 void SetTOF( G4double t ) { timeOfFlight = t; return; }
00315
00316 inline
00317 G4double GetTOF() { return timeOfFlight; }
00318
00319 inline
00320 void SetSide( G4int sid ) { side = sid; return; }
00321
00322 inline
00323 G4int GetSide() { return side; }
00324
00325 inline
00326 void setFlag( G4bool f ) { flag = f; return; }
00327
00328 inline
00329 G4bool getFlag() { return flag; }
00330
00331 inline
00332 void SetCode( G4int c ) { code = c; return; }
00333
00334 inline
00335 void SetParticleDef( G4ParticleDefinition * c ) { particleDef = c; return; }
00336
00337 inline
00338 G4int GetCode() { return code; }
00339
00340 inline
00341 G4ParticleDefinition * GetParticleDef() { return particleDef; }
00342
00343 inline
00344 void SetZero()
00345 {
00346 momentum.setX( 0.0 );
00347 momentum.setY( 0.0 );
00348 momentum.setZ( 0.0 );
00349 energy = 0.0;
00350 kineticEnergy = 0.0;
00351 mass = 0.0;
00352 charge = 0.0;
00353 timeOfFlight = 0.0;
00354 side = 0;
00355 flag = false;
00356 code = 0;
00357 particleDef = NULL;
00358 }
00359
00360 inline
00361 void Add( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2 )
00362 {
00363 momentum = p1.momentum + p2.momentum;
00364 energy = p1.energy + p2.energy;
00365 G4double b = energy*energy - momentum.mag2();
00366 if( b < 0 )
00367 mass = -1. * std::sqrt( -b );
00368 else
00369 mass = std::sqrt( b );
00370 kineticEnergy = std::max(0.,energy - mass);
00371 charge = p1.charge + p2.charge;
00372 code = p1.code + p2.code;
00373 particleDef = p1.particleDef;
00374 }
00375
00376 inline
00377 void Sub( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2 )
00378 {
00379 momentum = p1.momentum - p2.momentum;
00380 energy = p1.energy - p2.energy;
00381 G4double b = energy*energy - momentum.mag2();
00382 if( b < 0 )
00383 mass = -1. * std::sqrt( -b );
00384 else
00385 mass = std::sqrt( b );
00386 kineticEnergy = std::max(0.,energy - mass);
00387 charge = p1.charge - p2.charge;
00388 code = p1.code - p2.code;
00389 particleDef = p1.particleDef;
00390 }
00391
00392 inline
00393 void Lor( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2 )
00394 {
00395 G4double a;
00396 a = ( p1.momentum.dot(p2.momentum)/(p2.energy+p2.mass) - p1.energy ) / p2.mass;
00397 momentum.setX( p1.momentum.x()+a*p2.momentum.x() );
00398 momentum.setY( p1.momentum.y()+a*p2.momentum.y() );
00399 momentum.setZ( p1.momentum.z()+a*p2.momentum.z() );
00400 energy = std::sqrt( sqr(p1.mass) + momentum.mag2() );
00401 mass = p1.mass;
00402 kineticEnergy = std::max(0.,energy - mass);
00403 timeOfFlight = p1.timeOfFlight;
00404 side = p1.side;
00405 flag = p1.flag;
00406 code = p1.code;
00407 particleDef = p1.particleDef;
00408 }
00409
00410 inline
00411 G4double CosAng( const G4GHEKinematicsVector & p )
00412 {
00413 G4double a = std::sqrt( momentum.mag2() * p.momentum.mag2() );
00414 if( a != 0.0 )
00415 {
00416 a = (momentum.x()*p.momentum.x() +
00417 momentum.y()*p.momentum.y() +
00418 momentum.z()*p.momentum.z()) / a;
00419 if( std::fabs(a) > 1.0 ) a<0.0 ? a=-1.0 : a=1.0;
00420 }
00421 return a;
00422 }
00423 inline
00424 G4double Ang(const G4GHEKinematicsVector & p )
00425 {
00426 G4double a = std::sqrt( momentum.mag2() * p.momentum.mag2() );
00427 if( a != 0.0 )
00428 {
00429 a = (momentum.x()*p.momentum.x() +
00430 momentum.y()*p.momentum.y() +
00431 momentum.z()*p.momentum.z()) / a;
00432 if( std::fabs(a) > 1.0 ) a<0.0 ? a=-1.0 : a=1.0;
00433 }
00434 return std::acos(a);
00435 }
00436
00437 inline
00438 G4double Dot4( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
00439 {
00440 return ( p1.energy * p2.energy
00441 - p1.momentum.x() * p2.momentum.x()
00442 - p1.momentum.y() * p2.momentum.y()
00443 - p1.momentum.z() * p2.momentum.z() );
00444 }
00445
00446 inline
00447 G4double Impu( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
00448 {
00449 return ( - sqr( p1.energy - p2.energy)
00450 + sqr(p1.momentum.x() - p2.momentum.x())
00451 + sqr(p1.momentum.y() - p2.momentum.y())
00452 + sqr(p1.momentum.z() - p2.momentum.z()) );
00453 }
00454
00455 inline
00456 void Add3( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
00457 {
00458 momentum.setX( p1.momentum.x() + p2.momentum.x());
00459 momentum.setY( p1.momentum.y() + p2.momentum.y());
00460 momentum.setZ( p1.momentum.z() + p2.momentum.z());
00461 return;
00462 }
00463
00464 inline
00465 void Sub3( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
00466 {
00467 momentum.setX( p1.momentum.x() - p2.momentum.x());
00468 momentum.setY( p1.momentum.y() - p2.momentum.y());
00469 momentum.setZ( p1.momentum.z() - p2.momentum.z());
00470 return;
00471 }
00472
00473 inline
00474 void Cross( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
00475 {
00476 G4double px, py, pz;
00477 px = p1.momentum.y() * p2.momentum.z() - p1.momentum.z() * p2.momentum.y();
00478 py = p1.momentum.z() * p2.momentum.x() - p1.momentum.x() * p2.momentum.z();
00479 pz = p1.momentum.x() * p2.momentum.y() - p1.momentum.y() * p2.momentum.x();
00480 momentum.setX( px );
00481 momentum.setY( py );
00482 momentum.setZ( pz );
00483 return;
00484 }
00485
00486 inline
00487 G4double Dot( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
00488 {
00489 return ( p1.momentum.x() * p2.momentum.x()
00490 + p1.momentum.y() * p2.momentum.y()
00491 + p1.momentum.z() * p2.momentum.z() );
00492 }
00493
00494 inline
00495 void Smul( const G4GHEKinematicsVector & p, G4double h)
00496 {
00497 momentum.setX( h * p.momentum.x());
00498 momentum.setY( h * p.momentum.y());
00499 momentum.setZ( h * p.momentum.z());
00500 return;
00501 }
00502
00503 inline
00504 void SmulAndUpdate( const G4GHEKinematicsVector & p, G4double h)
00505 {
00506 momentum.setX( h * p.momentum.x());
00507 momentum.setY( h * p.momentum.y());
00508 momentum.setZ( h * p.momentum.z());
00509 mass = p.mass;
00510 energy = std::sqrt(momentum.mag2() + mass*mass);
00511 kineticEnergy = energy - mass;
00512 charge = p.charge;
00513 timeOfFlight = p.timeOfFlight;
00514 side = p.side;
00515 flag = p.flag;
00516 code = p.code;
00517 particleDef = p.particleDef;
00518 return;
00519 }
00520
00521 inline
00522 void Norz( const G4GHEKinematicsVector & p )
00523 {
00524 G4double a = p.momentum.mag2();
00525 if (a > 0.0) a = 1./std::sqrt(a);
00526 momentum.setX( a * p.momentum.x() );
00527 momentum.setY( a * p.momentum.y() );
00528 momentum.setZ( a * p.momentum.z() );
00529 mass = p.mass;
00530 energy = std::sqrt(momentum.mag2() + mass*mass);
00531 kineticEnergy = energy - mass;
00532 charge = p.charge;
00533 timeOfFlight = p.timeOfFlight;
00534 side = p.side;
00535 flag = p.flag;
00536 code = p.code;
00537 particleDef = p.particleDef;
00538 return;
00539 }
00540
00541 inline
00542 G4double Length()
00543 {
00544 return momentum.mag() ;
00545 }
00546
00547 inline
00548 void Exch( G4GHEKinematicsVector & p1)
00549 {
00550 G4GHEKinematicsVector mx = *this;
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573 *this = p1;
00574 p1 = mx;
00575 return;
00576 }
00577
00578 inline
00579 void Defs1( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
00580 {
00581 G4double pt2 = sqr(p1.momentum.x()) + sqr(p1.momentum.y());
00582 if (pt2 > 0.0)
00583 {
00584 G4double ph, px, py, pz;
00585 G4double cost = p2.momentum.z()/p2.momentum.mag();
00586 G4double sint = 0.5 * ( std::sqrt(std::fabs((1.-cost)*(1.+cost)))
00587 + std::sqrt(pt2)/p2.momentum.mag());
00588 (p2.momentum.y() < 0.) ? ph = 1.5*CLHEP::pi : ph = CLHEP::halfpi;
00589 if( p2.momentum.x() != 0.0)
00590 ph = std::atan2(p2.momentum.y(),p2.momentum.x());
00591 px = cost*std::cos(ph)*p1.momentum.x() - std::sin(ph)*p1.momentum.y()
00592 + sint*std::cos(ph)*p1.momentum.z();
00593 py = cost*std::sin(ph)*p1.momentum.x() + std::cos(ph)*p1.momentum.y()
00594 + sint*std::sin(ph)*p1.momentum.z();
00595 pz = - sint *p1.momentum.x()
00596 + cost *p1.momentum.z();
00597 momentum.setX( px );
00598 momentum.setY( py );
00599 momentum.setZ( pz );
00600 }
00601 else
00602 {
00603 momentum = p1.momentum;
00604 }
00605 }
00606
00607 inline
00608 void Defs( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2,
00609 G4GHEKinematicsVector & my, G4GHEKinematicsVector & mz )
00610 {
00611 my = p1;
00612 mz = p2;
00613 momentum.setX( my.momentum.y()*mz.momentum.z()
00614 - my.momentum.z()*mz.momentum.y());
00615 momentum.setY( my.momentum.z()*mz.momentum.x()
00616 - my.momentum.x()*mz.momentum.z());
00617 momentum.setZ( my.momentum.x()*mz.momentum.y()
00618 - my.momentum.y()*mz.momentum.x());
00619 my.momentum.setX( mz.momentum.y()*momentum.z()
00620 - mz.momentum.z()*momentum.y());
00621 my.momentum.setY( mz.momentum.z()*momentum.x()
00622 - mz.momentum.x()*momentum.z());
00623 my.momentum.setZ( mz.momentum.x()*momentum.y()
00624 - mz.momentum.y()*momentum.x());
00625 G4double pp;
00626 pp = momentum.mag();
00627 if (pp > 0.)
00628 {
00629 pp = 1./pp;
00630 momentum.setX( momentum.x()*pp );
00631 momentum.setY( momentum.y()*pp );
00632 momentum.setZ( momentum.z()*pp );
00633 }
00634 pp = my.momentum.mag();
00635 if (pp > 0.)
00636 {
00637 pp = 1./pp;
00638 my.momentum.setX( my.momentum.x()*pp );
00639 my.momentum.setY( my.momentum.y()*pp );
00640 my.momentum.setZ( my.momentum.z()*pp );
00641 }
00642 pp = mz.momentum.mag();
00643 if (pp > 0.)
00644 {
00645 pp = 1./pp;
00646 mz.momentum.setX( mz.momentum.x()*pp );
00647 mz.momentum.setY( mz.momentum.y()*pp );
00648 mz.momentum.setZ( mz.momentum.z()*pp );
00649 }
00650 return;
00651 }
00652
00653 inline
00654 void Trac( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & mx,
00655 const G4GHEKinematicsVector & my, const G4GHEKinematicsVector & mz)
00656 {
00657 double px, py, pz;
00658 px = mx.momentum.x()*p1.momentum.x()
00659 + mx.momentum.y()*p1.momentum.y()
00660 + mx.momentum.z()*p1.momentum.z();
00661 py = my.momentum.x()*p1.momentum.x()
00662 + my.momentum.y()*p1.momentum.y()
00663 + my.momentum.z()*p1.momentum.z();
00664 pz = mz.momentum.x()*p1.momentum.x()
00665 + mz.momentum.y()*p1.momentum.y()
00666 + mz.momentum.z()*p1.momentum.z();
00667 momentum.setX( px );
00668 momentum.setY( py );
00669 momentum.setZ( pz );
00670 return;
00671 }
00672
00673 inline
00674 void Print( G4int L)
00675 {
00676 G4cout << "G4GHEKinematicsVector: "
00677 << L << " " << momentum.x() << " " << momentum.y() << " " << momentum.z() << " "
00678 << energy << " " << kineticEnergy << " " << mass << " " << charge << " "
00679 << timeOfFlight << " " << side << " " << flag << " " << code << particleDef << G4endl;
00680 return;
00681 }
00682
00683 G4ParticleMomentum momentum;
00684 G4double energy;
00685 G4double kineticEnergy;
00686 G4double mass;
00687 G4double charge;
00688 G4double timeOfFlight;
00689 G4int side;
00690 G4bool flag;
00691 G4int code;
00692 G4ParticleDefinition * particleDef;
00693 };
00694
00695 #endif