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 #include "G4HEVector.hh"
00038 #include "globals.hh"
00039 #include "G4ios.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4SystemOfUnits.hh"
00042 #include "G4ParticleDefinition.hh"
00043
00044 G4HEVector::G4HEVector(const G4HadProjectile * aParticle)
00045 {
00046 G4ThreeVector aMom = 1./GeV*aParticle->Get4Momentum().vect();
00047 px = aMom.x();
00048 py = aMom.y();
00049 pz = aMom.z();
00050 energy = aParticle->GetTotalEnergy()/GeV;
00051 kineticEnergy = aParticle->GetKineticEnergy()/GeV;
00052 mass = aParticle->GetDefinition()->GetPDGMass()/GeV;
00053 charge = aParticle->GetDefinition()->GetPDGCharge()/eplus;
00054 timeOfFlight = 0.0;
00055 side = 0;
00056 flag = false;
00057 code = aParticle->GetDefinition()->GetPDGEncoding();
00058 baryon = aParticle->GetDefinition()->GetBaryonNumber();
00059 particleName = getParticleName(code, baryon);
00060 particleType = aParticle->GetDefinition()->GetParticleType();
00061 strangeness = aParticle->GetDefinition()->GetQuarkContent(3);
00062 }
00063
00064 G4double G4HEVector::Amax(G4double a, G4double b)
00065 {
00066 G4double c = a;
00067 if(b > a) c = b;
00068 return c;
00069 }
00070
00071 G4String G4HEVector::getParticleName(G4int aCode, G4int aBaryon)
00072 {
00073 G4String name;
00074 if(aCode == 211) name = "PionPlus";
00075 else if(aCode == 111) name = "PionZero";
00076 else if(aCode == -211) name = "PionMinus";
00077 else if(aCode == 321) name = "KaonPlus";
00078 else if(aCode == 311) name = "KaonZero";
00079 else if(aCode == -311) name = "AntiKaonZero";
00080 else if(aCode == -321) name = "KaonMinus";
00081 else if(aCode == 310) name = "KaonZeroShort";
00082 else if(aCode == 130) name = "KaonZeroLong";
00083 else if(aCode == 2212) name = "Proton";
00084 else if(aCode == -2212) name = "AntiProton";
00085 else if(aCode == 2112) name = "Neutron";
00086 else if(aCode == -2112) name = "AntiNeutron";
00087 else if(aCode == 3122) name = "Lambda";
00088 else if(aCode == -3122) name = "AntiLambda";
00089 else if(aCode == 3222) name = "SigmaPlus";
00090 else if(aCode == 3212) name = "SigmaZero";
00091 else if(aCode == 3112) name = "SigmaMinus";
00092 else if(aCode == -3222) name = "AntiSigmaPlus";
00093 else if(aCode == -3212) name = "AntiSigmaZero";
00094 else if(aCode == -3112) name = "AntiSigmaMinus";
00095 else if(aCode == 3322) name = "XiZero";
00096 else if(aCode == 3312) name = "XiMinus";
00097 else if(aCode == -3322) name = "AntiXiZero";
00098 else if(aCode == -3312) name = "AntiXiMinus";
00099 else if(aCode == 3334) name = "OmegaMinus";
00100 else if(aCode == -3334) name = "AntiOmegaMinus";
00101 else if(aCode == 0)
00102 {
00103 if(aBaryon==2) name = "Deuteron";
00104 else if(aBaryon==3) name = "Triton";
00105 else if(aBaryon==4) name = "Alpha";
00106 }
00107 else if(aCode == 22) name = "Gamma";
00108 else
00109 {
00110 G4cout << "particle " << aCode << " " <<aBaryon<< " not known in this generator!!" << G4endl;
00111 }
00112 return name;
00113 }
00114
00115
00116 void
00117 G4HEVector::setMomentum(const G4ParticleMomentum mom )
00118 {
00119 px = mom.x();
00120 py = mom.y();
00121 pz = mom.z();
00122 return;
00123 }
00124
00125 void
00126 G4HEVector::setMomentum(const G4ParticleMomentum * mom )
00127 {
00128 px = mom->x();
00129 py = mom->y();
00130 pz = mom->z();
00131 return;
00132 }
00133
00134 void
00135 G4HEVector::setMomentumAndUpdate( const G4ParticleMomentum mom )
00136 {
00137 px = mom.x();
00138 py = mom.y();
00139 pz = mom.z();
00140 energy = std::sqrt(mass*mass + px*px + py*py + pz*pz);
00141 kineticEnergy = Amax(0.,energy - mass);
00142 return;
00143 }
00144
00145 void
00146 G4HEVector::setMomentumAndUpdate( const G4ParticleMomentum * mom )
00147 {
00148 px = mom->x();
00149 py = mom->y();
00150 pz = mom->z();
00151 energy = std::sqrt(mass*mass + px*px + py*py + pz*pz);
00152 kineticEnergy = Amax(0.,energy - mass);
00153 return;
00154 }
00155
00156 const G4ParticleMomentum
00157 G4HEVector::getMomentum() const
00158 {
00159 G4ParticleMomentum mom;
00160 mom.setX(px);
00161 mom.setY(py);
00162 mom.setZ(pz);
00163 return mom;
00164 }
00165
00166 G4double G4HEVector::getTotalMomentum() const
00167 {
00168 return std::sqrt(px*px + py*py + pz*pz);
00169 }
00170
00171 void
00172 G4HEVector::setMomentum(G4double x, G4double y, G4double z)
00173 {
00174 px = x;
00175 py = y;
00176 pz = z;
00177 return;
00178 }
00179
00180 void
00181 G4HEVector::setMomentumAndUpdate(G4double x, G4double y, G4double z)
00182 {
00183 px = x;
00184 py = y;
00185 pz = z;
00186 energy = std::sqrt(mass*mass + px*px + py*py + pz*pz);
00187 kineticEnergy = Amax(0.,energy-mass);
00188 return;
00189 }
00190
00191 void
00192 G4HEVector::setMomentum(G4double x, G4double y)
00193 {
00194 px = x;
00195 py = y;
00196 return;
00197 }
00198
00199 void
00200 G4HEVector::setMomentumAndUpdate( G4double x, G4double y )
00201 {
00202 px = x;
00203 py = y;
00204 energy = std::sqrt(mass*mass + px*px + py*py + pz*pz);
00205 kineticEnergy = Amax(0.,energy-mass);
00206 return;
00207 }
00208
00209 void
00210 G4HEVector::setMomentum( G4double z )
00211 {
00212 pz = z;
00213 return;
00214 }
00215
00216 void
00217 G4HEVector::setMomentumAndUpdate( G4double z )
00218 {
00219 pz = z;
00220 energy = std::sqrt(mass*mass + px*px + py*py + pz*pz);
00221 kineticEnergy = Amax(0.,energy-mass);
00222 return;
00223 }
00224
00225 void
00226 G4HEVector::setEnergy( G4double e )
00227 {
00228 energy = e;
00229 return;
00230 }
00231
00232 void
00233 G4HEVector::setEnergyAndUpdate( G4double e )
00234 {
00235 if (e <= mass)
00236 {
00237 energy = mass;
00238 kineticEnergy = 0.;
00239 px = 0.;
00240 py = 0.;
00241 pz = 0.;
00242 }
00243 else
00244 {
00245 energy = e;
00246 kineticEnergy = energy - mass;
00247 G4double momold = std::sqrt(px*px + py*py + pz*pz);
00248 G4double momnew = std::sqrt(energy*energy - mass*mass);
00249 if (momold == 0.)
00250 {
00251 G4double cost = 1.0- 2.0*G4UniformRand();
00252 G4double sint = std::sqrt(1. - cost*cost);
00253 G4double phi = twopi* G4UniformRand();
00254 px = momnew * sint * std::cos(phi);
00255 py = momnew * sint * std::sin(phi);
00256 pz = momnew * cost;
00257 }
00258 else
00259 {
00260 momnew /= momold;
00261 px *= momnew;
00262 py *= momnew;
00263 pz *= momnew;
00264 }
00265 }
00266 return;
00267 }
00268
00269 void
00270 G4HEVector::setKineticEnergy( G4double ekin )
00271 {
00272 kineticEnergy = ekin;
00273 return;
00274 }
00275
00276 void
00277 G4HEVector::setKineticEnergyAndUpdate(G4double ekin)
00278 {
00279 if (ekin <= 0.)
00280 {
00281 energy = mass;
00282 kineticEnergy = 0.;
00283 px = 0.;
00284 py = 0.;
00285 pz = 0.;
00286 }
00287 else
00288 {
00289 energy = ekin + mass;
00290 kineticEnergy = ekin;
00291 G4double momold = std::sqrt(px*px + py*py + pz*pz);
00292 G4double momnew = std::sqrt(energy*energy - mass*mass);
00293 if (momold == 0.)
00294 {
00295 G4double cost = 1.0-2.0*G4UniformRand();
00296 G4double sint = std::sqrt(1. - cost*cost);
00297 G4double phi = twopi* G4UniformRand();
00298 px = momnew * sint * std::cos(phi);
00299 py = momnew * sint * std::sin(phi);
00300 pz = momnew * cost;
00301 }
00302 else
00303 {
00304 momnew /= momold;
00305 px *= momnew;
00306 py *= momnew;
00307 pz *= momnew;
00308 }
00309 }
00310 return;
00311 }
00312
00313 G4double G4HEVector::getEnergy() const
00314 {
00315 return energy;
00316 }
00317
00318 G4double G4HEVector::getKineticEnergy() const
00319 {
00320 return kineticEnergy;
00321 }
00322
00323
00324 void G4HEVector::setMass(G4double amass)
00325 {
00326 mass = amass;
00327 return;
00328 }
00329
00330
00331 void G4HEVector::setMassAndUpdate(G4double amass)
00332 {
00333 kineticEnergy = Amax(0., energy - mass);
00334 mass = amass;
00335 energy = kineticEnergy + mass;
00336 G4double momnew = std::sqrt(Amax(0., energy*energy - mass*mass));
00337 if (momnew == 0.0) {
00338 px = 0.;
00339 py = 0.;
00340 pz = 0.;
00341 } else {
00342 G4double momold = std::sqrt(px*px + py*py + pz*pz);
00343 if (momold == 0.) {
00344 G4double cost = 1.-2.*G4UniformRand();
00345 G4double sint = std::sqrt(1.-cost*cost);
00346 G4double phi = twopi*G4UniformRand();
00347 px = momnew*sint*std::cos(phi);
00348 py = momnew*sint*std::sin(phi);
00349 pz = momnew*cost;
00350 } else {
00351 momnew /= momold;
00352 px *= momnew ;
00353 py *= momnew ;
00354 pz *= momnew ;
00355 }
00356 }
00357 return;
00358 }
00359
00360
00361 G4double G4HEVector::getMass() const
00362 {
00363 return mass;
00364 }
00365
00366 void
00367 G4HEVector::setCharge( G4double c )
00368 {
00369 charge = c;
00370 return;
00371 }
00372
00373 G4double G4HEVector::getCharge() const
00374 {
00375 return charge;
00376 }
00377
00378 void
00379 G4HEVector::setTOF( G4double t )
00380 {
00381 timeOfFlight = t;
00382 return;
00383 }
00384
00385 G4double
00386 G4HEVector::getTOF()
00387 {
00388 return timeOfFlight;
00389 }
00390
00391
00392 void G4HEVector::setSide(G4int aside)
00393 {
00394 side = aside;
00395 return;
00396 }
00397
00398
00399 G4int
00400 G4HEVector::getSide()
00401 {
00402 return side;
00403 }
00404
00405
00406 void
00407 G4HEVector::setFlag( G4bool f )
00408 {
00409 flag = f;
00410 return;
00411 }
00412
00413 G4bool
00414 G4HEVector::getFlag()
00415 {
00416 return flag;
00417 }
00418
00419 void
00420 G4HEVector::setCode( G4int c )
00421 {
00422 code = c;
00423 return;
00424 }
00425
00426 G4int G4HEVector::getCode() const
00427 {
00428 return code;
00429 }
00430
00431 G4String G4HEVector::getName() const
00432 {
00433 return particleName;
00434 }
00435
00436 G4String G4HEVector::getType() const
00437 {
00438 return particleType;
00439 }
00440
00441 G4int G4HEVector::getBaryonNumber() const
00442 {
00443 return baryon;
00444 }
00445
00446 G4int G4HEVector::getStrangenessNumber() const
00447 {
00448 return strangeness;
00449 }
00450
00451 G4int
00452 G4HEVector::getQuarkContent(G4int flavor)
00453 {
00454 if(flavor > 0 && flavor < 8)
00455 {
00456 G4int check;
00457 check = FillQuarkContent();
00458 if(check != code)
00459 {
00460 return 0;
00461 }
00462 else
00463 {
00464 return theQuarkContent[flavor-1];
00465 }
00466 }
00467 else
00468 {
00469 return 0;
00470 }
00471 }
00472
00473 G4int
00474 G4HEVector::getAntiQuarkContent(G4int flavor)
00475 {
00476 if(flavor > 0 && flavor < 8)
00477 {
00478 G4int check;
00479 check = FillQuarkContent();
00480 if(check != code)
00481 {
00482 return 0;
00483 }
00484 else
00485 {
00486 return theAntiQuarkContent[flavor-1];
00487 }
00488 }
00489 else
00490 {
00491 return 0;
00492 }
00493 }
00494
00495 void
00496 G4HEVector::setZero()
00497 {
00498 px = 0.0;
00499 py = 0.0;
00500 pz = 0.0;
00501 energy = 0.0;
00502 kineticEnergy = 0.0;
00503 mass = 0.0;
00504 charge = 0.0;
00505 timeOfFlight = 0.0;
00506 side = 0;
00507 flag = false;
00508 code = 0;
00509 particleName = "";
00510 particleType = "";
00511 baryon = 0;
00512 strangeness = 0;
00513 }
00514
00515 void
00516 G4HEVector::Add( const G4HEVector & p1, const G4HEVector & p2 )
00517 {
00518 px = p1.px + p2.px;
00519 py = p1.py + p2.py;
00520 pz = p1.pz + p2.pz;
00521 energy = p1.energy + p2.energy;
00522 G4double b = energy*energy - px*px - py*py - pz*pz;
00523 if( b < 0 )
00524 mass = -1. * std::sqrt( -b );
00525 else
00526 mass = std::sqrt( b );
00527 kineticEnergy = Amax(0.,energy - mass);
00528 charge = p1.charge + p2.charge;
00529 code = 0;
00530 particleName = "";
00531 particleType = "";
00532 baryon = 0;
00533 strangeness = 0;
00534 }
00535
00536 void
00537 G4HEVector::Sub( const G4HEVector & p1, const G4HEVector & p2 )
00538 {
00539 px = p1.px - p2.px;
00540 py = p1.py - p2.py;
00541 pz = p1.pz - p2.pz;
00542 energy = p1.energy - p2.energy;
00543 G4double b = energy*energy - px*px - py*py - pz*pz;
00544 if( b < 0 )
00545 mass = -1. * std::sqrt( -b );
00546 else
00547 mass = std::sqrt( b );
00548 kineticEnergy = Amax(0.,energy - mass);
00549 charge = p1.charge - p2.charge;
00550 code = 0;
00551 particleName = "";
00552 particleType = "";
00553 baryon = 0;
00554 strangeness = 0;
00555 }
00556
00557 void
00558 G4HEVector::Lor( const G4HEVector & p1, const G4HEVector & p2 )
00559 {
00560 G4double a;
00561 a = ( Dot(p1,p2)/(p2.energy+p2.mass) - p1.energy ) / p2.mass;
00562 px = p1.px + a*p2.px;
00563 py = p1.py + a*p2.py;
00564 pz = p1.pz + a*p2.pz;
00565 energy = std::sqrt( sqr(p1.mass) + px*px + py*py + pz*pz);
00566 mass = p1.mass;
00567 kineticEnergy = Amax(0.,energy - mass);
00568 timeOfFlight = p1.timeOfFlight;
00569 side = p1.side;
00570 flag = p1.flag;
00571 code = p1.code;
00572 particleName = p1.particleName;
00573 particleType = p1.particleType;
00574 baryon = p1.baryon;
00575 strangeness = p1.strangeness;
00576 }
00577
00578 G4double G4HEVector::CosAng(const G4HEVector& p) const
00579 {
00580 G4double a = std::sqrt( (px*px + py*py + pz*pz)*(p.px*p.px + p.py*p.py + p.pz*p.pz) );
00581 if (a != 0.0) {
00582 a = (px*p.px + py*p.py + pz*p.pz)/a;
00583 if (std::fabs(a) > 1.0) {
00584 if (a < 0.0) a = -1.0;
00585 else a = 1.0;
00586 }
00587 }
00588 return a;
00589 }
00590
00591 G4double
00592 G4HEVector::Ang(const G4HEVector & p )
00593 {
00594 G4double a = std::sqrt( (px*px + py*py + pz*pz)*(p.px*p.px + p.py*p.py + p.pz*p.pz) );
00595 if( a != 0.0 )
00596 {
00597 a = (px*p.px + py*p.py + pz*p.pz)/a;
00598 if( std::fabs(a) > 1.0 )
00599 {
00600 if(a<0.0) a=-1.0;
00601 else a=1.0;
00602 }
00603 }
00604 return std::acos(a);
00605 }
00606
00607 G4double
00608 G4HEVector::Dot4( const G4HEVector & p1, const G4HEVector & p2)
00609 {
00610 return ( p1.energy*p2.energy - p1.px*p2.px - p1.py*p2.py - p1.pz*p2.pz );
00611 }
00612
00613 G4double
00614 G4HEVector::Impu( const G4HEVector & p1, const G4HEVector & p2)
00615 {
00616 return ( - sqr( p1.energy - p2.energy)
00617 + sqr( p1.px - p2.px)
00618 + sqr( p1.py - p2.py)
00619 + sqr( p1.pz - p2.pz) );
00620 }
00621
00622 void
00623 G4HEVector::Add3( const G4HEVector & p1, const G4HEVector & p2)
00624 {
00625 px = p1.px + p2.px;
00626 py = p1.py + p2.py;
00627 pz = p1.pz + p2.pz;
00628 return;
00629 }
00630
00631 void
00632 G4HEVector::Sub3( const G4HEVector & p1, const G4HEVector & p2)
00633 {
00634 px = p1.px - p2.px;
00635 py = p1.py - p2.py;
00636 pz = p1.pz - p2.pz;
00637 return;
00638 }
00639
00640 void
00641 G4HEVector::Cross( const G4HEVector & p1, const G4HEVector & p2)
00642 {
00643 G4double tpx = p1.py * p2.pz - p1.pz * p2.py;
00644 G4double tpy = p1.pz * p2.px - p1.px * p2.pz;
00645 G4double tpz = p1.px * p2.py - p1.py * p2.px;
00646 px=tpx;
00647 py=tpy;
00648 pz=tpz;
00649 return;
00650 }
00651
00652 G4double
00653 G4HEVector::Dot( const G4HEVector & p1, const G4HEVector & p2)
00654 {
00655 return ( p1.px * p2.px + p1.py * p2.py + p1.pz * p2.pz );
00656 }
00657
00658 void
00659 G4HEVector::Smul( const G4HEVector & p, G4double h)
00660 {
00661 px = h * p.px;
00662 py = h * p.py;
00663 pz = h * p.pz;
00664 return;
00665 }
00666
00667 void
00668 G4HEVector::SmulAndUpdate( const G4HEVector & p, G4double h)
00669 {
00670 px = h * p.px;
00671 py = h * p.py;
00672 pz = h * p.pz;
00673 mass = p.mass;
00674 energy = std::sqrt(px*px + py*py + pz*pz + mass*mass);
00675 kineticEnergy = energy - mass;
00676 charge = p.charge;
00677 timeOfFlight = p.timeOfFlight;
00678 side = p.side;
00679 flag = p.flag;
00680 code = p.code;
00681 particleName = p.particleName;
00682 particleType = p.particleType;
00683 baryon = p.baryon;
00684 strangeness = p.strangeness;
00685 return;
00686 }
00687
00688 void
00689 G4HEVector::Norz( const G4HEVector & p )
00690 {
00691 G4double a = p.px*p.px + p.py*p.py + p.pz*p.pz;
00692 if (a > 0.0) a = 1./std::sqrt(a);
00693 px = a * p.px;
00694 py = a * p.py;
00695 pz = a * p.pz;
00696 mass = p.mass;
00697 energy = std::sqrt(px*px + py*py + pz*pz + mass*mass);
00698 kineticEnergy = energy - mass;
00699 charge = p.charge;
00700 timeOfFlight = p.timeOfFlight;
00701 side = p.side;
00702 flag = p.flag;
00703 code = p.code;
00704 particleName = p.particleName;
00705 particleType = p.particleType;
00706 baryon = p.baryon;
00707 strangeness = p.strangeness;
00708 return;
00709 }
00710
00711 G4double G4HEVector::Length() const
00712 {
00713 return std::sqrt(px*px + py*py + pz*pz);
00714 }
00715
00716 void
00717 G4HEVector::Exch( G4HEVector & p1)
00718 {
00719 G4HEVector mx = *this;
00720 *this = p1;
00721 p1 = mx;
00722 return;
00723 }
00724
00725 void
00726 G4HEVector::Defs1( const G4HEVector & p1, const G4HEVector & p2)
00727 {
00728 G4double pt2 = p2.px*p2.px + p2.py*p2.py;
00729 if (pt2 > 0.0)
00730 {
00731 G4double ph, qx, qy, qz;
00732 G4double a = std::sqrt(p2.px*p2.px + p2.py*p2.py + p2.pz*p2.pz);
00733 G4double cost = p2.pz/a;
00734 G4double sint = 0.5 * (std::sqrt(std::fabs((1.-cost)*(1.+cost))) + std::sqrt(pt2)/a);
00735 if(p2.py < 0.) ph = 1.5*pi;
00736 else ph = halfpi;
00737 if( p2.px != 0.0)
00738 ph = std::atan2(p2.py,p2.px);
00739 qx = cost*std::cos(ph)*p1.px - std::sin(ph)*p1.py
00740 + sint*std::cos(ph)*p1.pz;
00741 qy = cost*std::sin(ph)*p1.px + std::cos(ph)*p1.py
00742 + sint*std::sin(ph)*p1.pz;
00743 qz = - sint *p1.px
00744 + cost *p1.pz;
00745 px = qx;
00746 py = qy;
00747 pz = qz;
00748 }
00749 else
00750 {
00751 px = p1.px;
00752 py = p1.py;
00753 pz = p1.pz;
00754 }
00755 }
00756
00757 void
00758 G4HEVector::Defs( const G4HEVector & p1, const G4HEVector & p2,
00759 G4HEVector & my, G4HEVector & mz )
00760 {
00761 my = p1;
00762 mz = p2;
00763 px = my.py*mz.pz - my.pz*mz.py;
00764 py = my.pz*mz.px - my.px*mz.pz;
00765 pz = my.px*mz.py - my.py*mz.px;
00766 my.px = mz.py*pz - mz.pz*py;
00767 my.py = mz.pz*px - mz.px*pz;
00768 my.pz = mz.px*py - mz.py*px;
00769 G4double pp;
00770 pp = std::sqrt(px*px + py*py + pz*pz);
00771 if (pp > 0.)
00772 {
00773 pp = 1./pp;
00774 px = px*pp ;
00775 py = py*pp ;
00776 pz = pz*pp ;
00777 }
00778 pp = std::sqrt(my.px*my.px + my.py*my.py + my.pz*my.pz);
00779 if (pp > 0.)
00780 {
00781 pp = 1./pp;
00782 my.px = my.px*pp ;
00783 my.py = my.py*pp ;
00784 my.pz = my.pz*pp ;
00785 }
00786 pp = std::sqrt(mz.px*mz.px + mz.py*mz.py + mz.pz*mz.pz);
00787 if (pp > 0.)
00788 {
00789 pp = 1./pp;
00790 mz.px = mz.px*pp ;
00791 mz.py = mz.py*pp ;
00792 mz.pz = mz.pz*pp ;
00793 }
00794 return;
00795 }
00796
00797 void
00798 G4HEVector::Trac( const G4HEVector & p1, const G4HEVector & mx,
00799 const G4HEVector & my, const G4HEVector & mz)
00800 {
00801 G4double qx, qy, qz;
00802 qx = mx.px*p1.px + mx.py*p1.py + mx.pz*p1.pz;
00803 qy = my.px*p1.px + my.py*p1.py + my.pz*p1.pz;
00804 qz = mz.px*p1.px + mz.py*p1.py + mz.pz*p1.pz;
00805 px = qx ;
00806 py = qy ;
00807 pz = qz ;
00808 return;
00809 }
00810
00811 void
00812 G4HEVector::setDefinition(G4String name)
00813 {
00814 if(name == "PionPlus")
00815 {
00816 mass = 0.1395700;
00817 charge = 1.;
00818 code = 211;
00819 particleType = "Meson";
00820 particleName = name;
00821 baryon = 0;
00822 strangeness = 0;
00823 }
00824 else if(name == "PionZero")
00825 {
00826 mass = 0.1349764;
00827 charge = 0.;
00828 code = 111;
00829 particleType = "Meson";
00830 particleName = name;
00831 baryon = 0;
00832 strangeness = 0;
00833 }
00834 else if(name == "PionMinus")
00835 {
00836 mass = 0.1395700;
00837 charge = -1.;
00838 code = -211;
00839 particleType = "Meson";
00840 particleName = name;
00841 baryon = 0;
00842 strangeness = 0;
00843 }
00844 else if(name == "KaonPlus")
00845 {
00846 mass = 0.493677;
00847 charge = 1.;
00848 code = 321;
00849 particleType = "Meson";
00850 particleName = name;
00851 baryon = 0;
00852 strangeness = 1;
00853 }
00854 else if(name == "KaonZero")
00855 {
00856 mass = 0.497672;
00857 charge = 0.;
00858 code = 311;
00859 particleType = "Meson";
00860 particleName = name;
00861 baryon = 0;
00862 strangeness = 1;
00863 }
00864 else if(name == "AntiKaonZero")
00865 {
00866 mass = 0.497672;
00867 charge = 0.;
00868 code = -311;
00869 particleType = "Meson";
00870 particleName = name;
00871 baryon = 0;
00872 strangeness =-1;
00873 }
00874 else if(name == "KaonMinus")
00875 {
00876 mass = 0.493677;
00877 charge = -1.;
00878 code = -321;
00879 particleType = "Meson";
00880 particleName = name;
00881 baryon = 0;
00882 strangeness = -1;
00883 }
00884 else if(name == "KaonZeroShort")
00885 {
00886 mass = 0.497672;
00887 charge = 0.;
00888 code = 310;
00889 particleType = "Meson";
00890 particleName = name;
00891 baryon = 0;
00892 strangeness = 0;
00893 }
00894 else if(name == "KaonZeroLong")
00895 {
00896 mass = 0.497672;
00897 charge = 0.;
00898 code = 130;
00899 particleType = "Meson";
00900 particleName = name;
00901 baryon = 0;
00902 strangeness = 0;
00903 }
00904 else if(name == "Proton")
00905 {
00906 mass = 0.9382723;
00907 charge = 1.;
00908 code = 2212;
00909 particleType = "Baryon";
00910 particleName = name;
00911 baryon = 1;
00912 strangeness = 0;
00913 }
00914 else if(name == "AntiProton")
00915 {
00916 mass = 0.9382723;
00917 charge = -1.;
00918 code = -2212;
00919 particleType = "AntiBaryon";
00920 particleName = name;
00921 baryon = -1;
00922 strangeness = 0;
00923 }
00924 else if(name == "Neutron")
00925 {
00926 mass = 0.93956563;
00927 charge = 0.;
00928 code = 2112;
00929 particleType = "Baryon";
00930 particleName = name;
00931 baryon = 1;
00932 strangeness = 0;
00933 }
00934 else if(name == "AntiNeutron")
00935 {
00936 mass = 0.93956563;
00937 charge = 0.;
00938 code = -2112;
00939 particleType = "AntiBaryon";
00940 particleName = name;
00941 baryon = -1;
00942 strangeness = 0;
00943 }
00944 else if(name == "Lambda")
00945 {
00946 mass = 1.115684;
00947 charge = 0.;
00948 code = 3122;
00949 particleType = "Baryon";
00950 particleName = name;
00951 baryon = 1;
00952 strangeness = -1;
00953 }
00954 else if(name == "AntiLambda")
00955 {
00956 mass = 1.115684;
00957 charge = 0.;
00958 code = -3122;
00959 particleType = "AntiBaryon";
00960 particleName = name;
00961 baryon = -1;
00962 strangeness = 1;
00963 }
00964 else if(name == "SigmaPlus")
00965 {
00966 mass = 1.18937;
00967 charge = 1.;
00968 code = 3222;
00969 particleType = "Baryon";
00970 particleName = name;
00971 baryon = 1;
00972 strangeness = -1;
00973 }
00974 else if(name == "SigmaZero")
00975 {
00976 mass = 1.19255;
00977 charge = 0.;
00978 code = 3212;
00979 particleType = "Baryon";
00980 particleName = name;
00981 baryon = 1;
00982 strangeness = -1;
00983 }
00984 else if(name == "SigmaMinus")
00985 {
00986 mass = 1.19744;
00987 charge = -1.;
00988 code = 3112;
00989 particleType = "Baryon";
00990 particleName = name;
00991 baryon = 1;
00992 strangeness = -1;
00993 }
00994 else if(name == "AntiSigmaPlus")
00995 {
00996 mass = 1.18937;
00997 charge = -1.;
00998 code = -3222;
00999 particleType = "AntiBaryon";
01000 particleName = name;
01001 baryon = -1;
01002 strangeness = 1;
01003 }
01004 else if(name == "AntiSigmaZero")
01005 {
01006 mass = 1.19255;
01007 charge = 0.;
01008 code = -3212;
01009 particleType = "AntiBaryon";
01010 particleName = name;
01011 baryon = -1;
01012 strangeness = 1;
01013 }
01014 else if(name == "AntiSigmaMinus")
01015 {
01016 mass = 1.19744;
01017 charge = 1.;
01018 code = -3112;
01019 particleType = "AntiBaryon";
01020 particleName = name;
01021 baryon = -1;
01022 strangeness = 1;
01023 }
01024 else if(name == "XiZero")
01025 {
01026 mass = 1.3149;
01027 charge = 0.;
01028 code = 3322;
01029 particleType = "Baryon";
01030 particleName = name;
01031 baryon = 1;
01032 strangeness = -2;
01033 }
01034 else if(name == "XiMinus")
01035 {
01036 mass = 1.32132;
01037 charge = -1.;
01038 code = 3312;
01039 particleType = "Baryon";
01040 particleName = name;
01041 baryon = 1;
01042 strangeness = -2;
01043 }
01044 else if(name == "AntiXiZero")
01045 {
01046 mass = 1.3149;
01047 charge = 0.;
01048 code = -3322;
01049 particleType = "AntiBaryon";
01050 particleName = name;
01051 baryon = -1;
01052 strangeness = 2;
01053 }
01054 else if(name == "AntiXiMinus")
01055 {
01056 mass = 1.32132;
01057 charge = 1.;
01058 code = -3312;
01059 particleType = "AntiBaryon";
01060 particleName = name;
01061 baryon = -1;
01062 strangeness = 2;
01063 }
01064 else if(name == "OmegaMinus")
01065 {
01066 mass = 1.67245;
01067 charge = -1.;
01068 code = 3334;
01069 particleType = "Baryon";
01070 particleName = name;
01071 baryon = 1;
01072 strangeness = -3;
01073 }
01074 else if(name == "AntiOmegaMinus")
01075 {
01076 mass = 1.67245;
01077 charge = 1.;
01078 code = -3334;
01079 particleType = "AntiBaryon";
01080 particleName = name;
01081 baryon = -1;
01082 strangeness = 3;
01083 }
01084 else if(name == "Deuteron")
01085 {
01086 mass = 1.875613;
01087 charge = 1.;
01088 code = 0;
01089 particleType = "Nucleus";
01090 particleName = name;
01091 baryon = 2;
01092 strangeness = 0;
01093 }
01094 else if(name == "Triton")
01095 {
01096 mass = 2.80925;
01097 charge = 1.;
01098 code = 0;
01099 particleType = "Nucleus";
01100 particleName = name;
01101 baryon = 3;
01102 strangeness = 0;
01103 }
01104 else if(name == "Alpha")
01105 {
01106 mass = 3.727417;
01107 charge = 2.;
01108 code = 0;
01109 particleType = "Nucleus";
01110 particleName = name;
01111 baryon = 4;
01112 strangeness = 0;
01113 }
01114 else if(name == "Gamma")
01115 {
01116 mass = 0.;
01117 charge = 0.;
01118 code = 22;
01119 particleType = "Boson";
01120 particleName = name;
01121 baryon = 0;
01122 strangeness = 0;
01123 }
01124 else
01125 {
01126 G4cout << "particle " << name << " not known in this generator!!" << G4endl;
01127 return;
01128 }
01129 px = 0.;
01130 py = 0.;
01131 pz = 0.;
01132 kineticEnergy = 0.;
01133 energy = mass;
01134 timeOfFlight = 0.;
01135 side = 0;
01136 flag = false;
01137 return;
01138 }
01139
01140 G4int G4HEVector::FillQuarkContent()
01141 {
01142
01143
01144
01145
01146 G4int tempPDGcode = code;
01147 G4double echarge = 1.;
01148
01149 for (G4int flavor=0; flavor<NumberOfQuarkFlavor; flavor++){
01150 theQuarkContent[flavor] =0;
01151 theAntiQuarkContent[flavor] =0;
01152 }
01153
01154 G4int temp = std::abs(tempPDGcode);
01155 G4int multiplet = temp/10000;
01156 temp -= G4int(multiplet*10000);
01157 G4int quark1 = temp/1000;
01158 temp -= G4int(quark1*1000);
01159 G4int quark2 = temp/100;
01160 temp -= G4int(quark2*100);
01161 G4int quark3 = temp/10;
01162 temp -= G4int(quark3*10);
01163
01164
01165
01166 if (particleType =="quark") {
01167 if (tempPDGcode>0){
01168 if (tempPDGcode<=NumberOfQuarkFlavor){
01169 theQuarkContent[tempPDGcode-1] =1;
01170 } else {
01171
01172 tempPDGcode = 0;
01173 }
01174 } else {
01175 G4int temp0 = -1*tempPDGcode;
01176 if (temp0 <= NumberOfQuarkFlavor) {
01177 theAntiQuarkContent[temp0-1] =1;
01178 } else {
01179
01180 tempPDGcode = 0;
01181 }
01182 }
01183
01184 } else if (particleType == "Meson") {
01185
01186
01187
01188
01189 if (tempPDGcode == 130) {
01190
01191 quark2 = 3;
01192 quark3 = 1;
01193 }
01194
01195 if (quark1 !=0)
01196 {
01197 tempPDGcode = 0;
01198 }
01199 if ((quark2==0)||(quark3==0)){
01200 tempPDGcode = 0;
01201 }
01202 if (quark2<quark3) {
01203 tempPDGcode = 0;
01204 }
01205
01206 if (quark2>=NumberOfQuarkFlavor){
01207 tempPDGcode = 0;
01208 }
01209
01210 if (quark2 & 1) {
01211
01212 if (tempPDGcode >0) {
01213 theQuarkContent[quark3-1] =1;
01214 theAntiQuarkContent[quark2-1] =1;
01215 } else {
01216 theQuarkContent[quark2-1] =1;
01217 theAntiQuarkContent[quark3-1] =1;
01218 }
01219 } else {
01220
01221 if (tempPDGcode >0) {
01222 theQuarkContent[quark2-1] =1;
01223 theAntiQuarkContent[quark3-1] =1;
01224 } else {
01225 theQuarkContent[quark3-1] =1;
01226 theAntiQuarkContent[quark2-1] =1;
01227 }
01228 }
01229
01230 G4double totalCharge = 0.0;
01231 for (G4int flavor= 0; flavor<NumberOfQuarkFlavor-1; flavor+=2){
01232 totalCharge += (-1./3.)*echarge*theQuarkContent[flavor];
01233 totalCharge += 1./3.*echarge*theAntiQuarkContent[flavor];
01234 totalCharge += 2./3.*echarge*theQuarkContent[flavor+1];
01235 totalCharge += (-2./3.)*echarge*theAntiQuarkContent[flavor+1];
01236 }
01237 if (std::abs(totalCharge-charge)>0.1*echarge) {
01238 tempPDGcode = 0;
01239 }
01240 } else if (particleType == "baryon"){
01241
01242 if ((quark1==0)||(quark2==0)||(quark3==0)){
01243 tempPDGcode = 0;
01244 }
01245
01246 if (std::abs(tempPDGcode) == 3122) {
01247
01248 quark2 = 2;
01249 quark3 = 1;
01250
01251 } else if (std::abs(tempPDGcode) == 4122) {
01252
01253 quark2 = 2;
01254 quark3 = 1;
01255
01256 }
01257
01258 if ((quark1<quark2)||(quark2<quark3)||(quark1<quark3)) {
01259 tempPDGcode = 0;
01260 }
01261 if (quark1>=NumberOfQuarkFlavor) {
01262 tempPDGcode = 0;
01263 }
01264 if (tempPDGcode >0) {
01265 theQuarkContent[quark1-1] ++;
01266 theQuarkContent[quark2-1] ++;
01267 theQuarkContent[quark3-1] ++;
01268 } else {
01269 theAntiQuarkContent[quark1-1] ++;
01270 theAntiQuarkContent[quark2-1] ++;
01271 theAntiQuarkContent[quark3-1] ++;
01272 }
01273
01274 G4double totalCharge = 0.0;
01275 for (G4int flavor= 0; flavor<NumberOfQuarkFlavor-1; flavor+=2){
01276 totalCharge += (-1./3.)*echarge*theQuarkContent[flavor];
01277 totalCharge += 1./3.*echarge*theAntiQuarkContent[flavor];
01278 totalCharge += 2./3.*echarge*theQuarkContent[flavor+1];
01279 totalCharge += (-2./3.)*echarge*theAntiQuarkContent[flavor+1];
01280 }
01281 if (std::abs(totalCharge - charge) > 0.1*echarge) tempPDGcode = 0;
01282
01283 } else {
01284 }
01285 return tempPDGcode;
01286 }
01287
01288
01289 void G4HEVector::Print(G4int L) const
01290 {
01291 G4cout << "HEV: "
01292 << L << " " << px << " " << py << " " << pz << " "
01293 << energy << " " << mass << " " << charge << " "
01294 << timeOfFlight << " " << side << " " << flag << " "
01295 << code << " " << baryon << " " << particleName << G4endl;
01296 return;
01297 }