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 #define INCLXX_IN_GEANT4_MODE 1
00034
00035 #include "globals.hh"
00036
00037
00038
00039
00040
00041
00042
00043
00044 #ifndef PARTICLE_HH_
00045 #define PARTICLE_HH_
00046
00047 #include "G4INCLThreeVector.hh"
00048 #include "G4INCLParticleTable.hh"
00049 #include "G4INCLParticleType.hh"
00050 #include "G4INCLParticleSpecies.hh"
00051 #include "G4INCLLogger.hh"
00052 #include <list>
00053 #include <sstream>
00054 #include <string>
00055 #include <algorithm>
00056
00057 namespace G4INCL {
00058
00059 class Particle;
00060
00061 typedef std::list<G4INCL::Particle*> ParticleList;
00062 typedef std::list<G4INCL::Particle*>::const_iterator ParticleIter;
00063
00064 class Particle {
00065 public:
00066 Particle();
00067 Particle(ParticleType t, G4double energy, ThreeVector momentum, ThreeVector position);
00068 Particle(ParticleType t, ThreeVector momentum, ThreeVector position);
00069 virtual ~Particle() {}
00070
00075 Particle(const Particle &rhs) :
00076 theZ(rhs.theZ),
00077 theA(rhs.theA),
00078 theParticipantType(rhs.theParticipantType),
00079 theType(rhs.theType),
00080 theEnergy(rhs.theEnergy),
00081 theFrozenEnergy(rhs.theFrozenEnergy),
00082 theMomentum(rhs.theMomentum),
00083 theFrozenMomentum(rhs.theFrozenMomentum),
00084 thePosition(rhs.thePosition),
00085 nCollisions(rhs.nCollisions),
00086 nDecays(rhs.nDecays),
00087 thePotentialEnergy(rhs.thePotentialEnergy),
00088 theHelicity(rhs.theHelicity),
00089 emissionTime(rhs.emissionTime),
00090 outOfWell(rhs.outOfWell),
00091 theMass(rhs.theMass)
00092 {
00093 if(rhs.thePropagationEnergy == &(rhs.theFrozenEnergy))
00094 thePropagationEnergy = &theFrozenEnergy;
00095 else
00096 thePropagationEnergy = &theEnergy;
00097 if(rhs.thePropagationMomentum == &(rhs.theFrozenMomentum))
00098 thePropagationMomentum = &theFrozenMomentum;
00099 else
00100 thePropagationMomentum = &theMomentum;
00101
00102 ID = nextID++;
00103 }
00104
00105 protected:
00107 void swap(Particle &rhs) {
00108 std::swap(theZ, rhs.theZ);
00109 std::swap(theA, rhs.theA);
00110 std::swap(theParticipantType, rhs.theParticipantType);
00111 std::swap(theType, rhs.theType);
00112 if(rhs.thePropagationEnergy == &(rhs.theFrozenEnergy))
00113 thePropagationEnergy = &theFrozenEnergy;
00114 else
00115 thePropagationEnergy = &theEnergy;
00116 std::swap(theEnergy, rhs.theEnergy);
00117 std::swap(theFrozenEnergy, rhs.theFrozenEnergy);
00118 if(rhs.thePropagationMomentum == &(rhs.theFrozenMomentum))
00119 thePropagationMomentum = &theFrozenMomentum;
00120 else
00121 thePropagationMomentum = &theMomentum;
00122 std::swap(theMomentum, rhs.theMomentum);
00123 std::swap(theFrozenMomentum, rhs.theFrozenMomentum);
00124 std::swap(thePosition, rhs.thePosition);
00125 std::swap(nCollisions, rhs.nCollisions);
00126 std::swap(nDecays, rhs.nDecays);
00127 std::swap(thePotentialEnergy, rhs.thePotentialEnergy);
00128
00129
00130 std::swap(theHelicity, rhs.theHelicity);
00131 std::swap(emissionTime, rhs.emissionTime);
00132 std::swap(outOfWell, rhs.outOfWell);
00133
00134 std::swap(theMass, rhs.theMass);
00135 }
00136
00137 public:
00138
00143 Particle &operator=(const Particle &rhs) {
00144 Particle temporaryParticle(rhs);
00145 swap(temporaryParticle);
00146 return *this;
00147 }
00148
00153 G4INCL::ParticleType getType() const {
00154 return theType;
00155 };
00156
00158 virtual G4INCL::ParticleSpecies getSpecies() const {
00159 return ParticleSpecies(theType);
00160 };
00161
00162 void setType(ParticleType t) {
00163 theType = t;
00164 switch(theType)
00165 {
00166 case DeltaPlusPlus:
00167 theA = 1;
00168 theZ = 2;
00169 break;
00170 case Proton:
00171 case DeltaPlus:
00172 theA = 1;
00173 theZ = 1;
00174 break;
00175 case Neutron:
00176 case DeltaZero:
00177 theA = 1;
00178 theZ = 0;
00179 break;
00180 case DeltaMinus:
00181 theA = 1;
00182 theZ = -1;
00183 break;
00184 case PiPlus:
00185 theA = 0;
00186 theZ = 1;
00187 break;
00188 case PiZero:
00189 theA = 0;
00190 theZ = 0;
00191 break;
00192 case PiMinus:
00193 theA = 0;
00194 theZ = -1;
00195 break;
00196 case Composite:
00197
00198 theA = 0;
00199 theZ = 0;
00200 break;
00201 case UnknownParticle:
00202 theA = 0;
00203 theZ = 0;
00204 ERROR("Trying to set particle type to Unknown!" << std::endl);
00205 break;
00206 }
00207
00208 if( !isResonance() && t!=Composite )
00209 setINCLMass();
00210 }
00211
00215 G4bool isNucleon() const {
00216 if(theType == G4INCL::Proton || theType == G4INCL::Neutron)
00217 return true;
00218 else
00219 return false;
00220 };
00221
00222 ParticipantType getParticipantType() const {
00223 return theParticipantType;
00224 }
00225
00226 void setParticipantType(ParticipantType const p) {
00227 theParticipantType = p;
00228 }
00229
00230 G4bool isParticipant() const {
00231 return (theParticipantType==Participant);
00232 }
00233
00234 G4bool isTargetSpectator() const {
00235 return (theParticipantType==TargetSpectator);
00236 }
00237
00238 G4bool isProjectileSpectator() const {
00239 return (theParticipantType==ProjectileSpectator);
00240 }
00241
00242 virtual void makeParticipant() {
00243 theParticipantType = Participant;
00244 }
00245
00246 virtual void makeTargetSpectator() {
00247 theParticipantType = TargetSpectator;
00248 }
00249
00250 virtual void makeProjectileSpectator() {
00251 theParticipantType = ProjectileSpectator;
00252 }
00253
00255 G4bool isPion() const { return (theType == PiPlus || theType == PiZero || theType == PiMinus); }
00256
00258 inline G4bool isResonance() const { return isDelta(); }
00259
00261 inline G4bool isDelta() const {
00262 return (theType==DeltaPlusPlus || theType==DeltaPlus ||
00263 theType==DeltaZero || theType==DeltaMinus);
00264 }
00265
00267 G4int getA() const { return theA; }
00268
00270 G4int getZ() const { return theZ; }
00271
00272 G4double getBeta() const {
00273 const G4double P = theMomentum.mag();
00274 return P/theEnergy;
00275 }
00276
00283 ThreeVector boostVector() const {
00284 return theMomentum / theEnergy;
00285 }
00286
00293 void boost(const ThreeVector &aBoostVector) {
00294 const G4double beta2 = aBoostVector.mag2();
00295 const G4double gamma = 1.0 / std::sqrt(1.0 - beta2);
00296 const G4double bp = theMomentum.dot(aBoostVector);
00297 const G4double alpha = (gamma*gamma)/(1.0 + gamma);
00298
00299 theMomentum = theMomentum + aBoostVector * (alpha * bp - gamma * theEnergy);
00300 theEnergy = gamma * (theEnergy - bp);
00301 }
00302
00311 void lorentzContract(const ThreeVector &aBoostVector, const ThreeVector &refPos) {
00312 const G4double beta2 = aBoostVector.mag2();
00313 const G4double gamma = 1.0 / std::sqrt(1.0 - beta2);
00314 const ThreeVector theRelativePosition = thePosition - refPos;
00315 const ThreeVector transversePosition = theRelativePosition - aBoostVector * (theRelativePosition.dot(aBoostVector) / aBoostVector.mag2());
00316 const ThreeVector longitudinalPosition = theRelativePosition - transversePosition;
00317
00318 thePosition = refPos + transversePosition + longitudinalPosition / gamma;
00319 }
00320
00322 inline G4double getMass() const { return theMass; }
00323
00325 inline G4double getINCLMass() const {
00326 switch(theType) {
00327 case Proton:
00328 case Neutron:
00329 case PiPlus:
00330 case PiMinus:
00331 case PiZero:
00332 return ParticleTable::getINCLMass(theType);
00333 break;
00334
00335 case DeltaPlusPlus:
00336 case DeltaPlus:
00337 case DeltaZero:
00338 case DeltaMinus:
00339 return theMass;
00340 break;
00341
00342 case Composite:
00343 return ParticleTable::getINCLMass(theA,theZ);
00344 break;
00345
00346 default:
00347 ERROR("Particle::getINCLMass: Unknown particle type." << std::endl);
00348 return 0.0;
00349 break;
00350 }
00351 }
00352
00354 inline virtual G4double getTableMass() const {
00355 switch(theType) {
00356 case Proton:
00357 case Neutron:
00358 case PiPlus:
00359 case PiMinus:
00360 case PiZero:
00361 return ParticleTable::getTableParticleMass(theType);
00362 break;
00363
00364 case DeltaPlusPlus:
00365 case DeltaPlus:
00366 case DeltaZero:
00367 case DeltaMinus:
00368 return theMass;
00369 break;
00370
00371 case Composite:
00372 return ParticleTable::getTableMass(theA,theZ);
00373 break;
00374
00375 default:
00376 ERROR("Particle::getTableMass: Unknown particle type." << std::endl);
00377 return 0.0;
00378 break;
00379 }
00380 }
00381
00383 inline G4double getRealMass() const {
00384 switch(theType) {
00385 case Proton:
00386 case Neutron:
00387 case PiPlus:
00388 case PiMinus:
00389 case PiZero:
00390 return ParticleTable::getRealMass(theType);
00391 break;
00392
00393 case DeltaPlusPlus:
00394 case DeltaPlus:
00395 case DeltaZero:
00396 case DeltaMinus:
00397 return theMass;
00398 break;
00399
00400 case Composite:
00401 return ParticleTable::getRealMass(theA,theZ);
00402 break;
00403
00404 default:
00405 ERROR("Particle::getRealMass: Unknown particle type." << std::endl);
00406 return 0.0;
00407 break;
00408 }
00409 }
00410
00412 void setRealMass() { setMass(getRealMass()); }
00413
00415 void setTableMass() { setMass(getTableMass()); }
00416
00418 void setINCLMass() { setMass(getINCLMass()); }
00419
00431 G4double getEmissionQValueCorrection(const G4int AParent, const G4int ZParent) const {
00432 const G4int ADaughter = AParent - theA;
00433 const G4int ZDaughter = ZParent - theZ;
00434
00435
00436 G4double theQValue;
00437 if(isCluster())
00438 theQValue = -ParticleTable::getTableQValue(theA, theZ, ADaughter, ZDaughter);
00439 else {
00440 const G4double massTableParent = ParticleTable::getTableMass(AParent,ZParent);
00441 const G4double massTableDaughter = ParticleTable::getTableMass(ADaughter,ZDaughter);
00442 const G4double massTableParticle = getTableMass();
00443 theQValue = massTableParent - massTableDaughter - massTableParticle;
00444 }
00445
00446 const G4double massINCLParent = ParticleTable::getINCLMass(AParent,ZParent);
00447 const G4double massINCLDaughter = ParticleTable::getINCLMass(ADaughter,ZDaughter);
00448 const G4double massINCLParticle = getINCLMass();
00449
00450
00451 return theQValue - (massINCLParent-massINCLDaughter-massINCLParticle);
00452 }
00453
00469 G4double getTransferQValueCorrection(const G4int AFrom, const G4int ZFrom, const G4int ATo, const G4int ZTo) const {
00470 const G4int AFromDaughter = AFrom - theA;
00471 const G4int ZFromDaughter = ZFrom - theZ;
00472 const G4int AToDaughter = ATo + theA;
00473 const G4int ZToDaughter = ZTo + theZ;
00474 const G4double theQValue = ParticleTable::getTableQValue(AToDaughter,ZToDaughter,AFromDaughter,ZFromDaughter,AFrom,ZFrom);
00475
00476 const G4double massINCLTo = ParticleTable::getINCLMass(ATo,ZTo);
00477 const G4double massINCLToDaughter = ParticleTable::getINCLMass(AToDaughter,ZToDaughter);
00478
00479
00480
00481
00482
00483 const G4double massINCLParticle = getTableMass();
00484
00485
00486 return theQValue - (massINCLToDaughter-massINCLTo-massINCLParticle);
00487 }
00488
00494 G4double getInvariantMass() const {
00495 const G4double mass = std::pow(theEnergy, 2) - theMomentum.dot(theMomentum);
00496 if(mass < 0.0) {
00497 ERROR("E*E - p*p is negative." << std::endl);
00498 return 0.0;
00499 } else {
00500 return std::sqrt(mass);
00501 }
00502 };
00503
00505 inline G4double getKineticEnergy() const { return theEnergy - theMass; }
00506
00508 inline G4double getPotentialEnergy() const { return thePotentialEnergy; }
00509
00511 inline void setPotentialEnergy(G4double v) { thePotentialEnergy = v; }
00512
00516 G4double getEnergy() const
00517 {
00518 return theEnergy;
00519 };
00520
00524 void setMass(G4double mass)
00525 {
00526 this->theMass = mass;
00527 }
00528
00532 void setEnergy(G4double energy)
00533 {
00534 this->theEnergy = energy;
00535 };
00536
00540 const G4INCL::ThreeVector &getMomentum() const
00541 {
00542 return theMomentum;
00543 };
00544
00546 virtual G4INCL::ThreeVector getAngularMomentum() const
00547 {
00548 return thePosition.vector(theMomentum);
00549 };
00550
00554 virtual void setMomentum(const G4INCL::ThreeVector &momentum)
00555 {
00556 this->theMomentum = momentum;
00557 };
00558
00562 const G4INCL::ThreeVector &getPosition() const
00563 {
00564 return thePosition;
00565 };
00566
00567 virtual void setPosition(const G4INCL::ThreeVector &position)
00568 {
00569 this->thePosition = position;
00570 };
00571
00572 G4double getHelicity() { return theHelicity; };
00573 void setHelicity(G4double h) { theHelicity = h; };
00574
00575 void propagate(G4double step) {
00576 thePosition += ((*thePropagationMomentum)*(step/(*thePropagationEnergy)));
00577 };
00578
00580 G4int getNumberOfCollisions() const { return nCollisions; }
00581
00583 void setNumberOfCollisions(G4int n) { nCollisions = n; }
00584
00586 void incrementNumberOfCollisions() { nCollisions++; }
00587
00589 G4int getNumberOfDecays() const { return nDecays; }
00590
00592 void setNumberOfDecays(G4int n) { nDecays = n; }
00593
00595 void incrementNumberOfDecays() { nDecays++; }
00596
00605 void setOutOfWell() { outOfWell = true; }
00606
00608 G4bool isOutOfWell() const { return outOfWell; }
00609
00610 void setEmissionTime(G4double t) { emissionTime = t; }
00611 G4double getEmissionTime() { return emissionTime; };
00612
00614 ThreeVector getTransversePosition() const {
00615 return thePosition - getLongitudinalPosition();
00616 }
00617
00619 ThreeVector getLongitudinalPosition() const {
00620 return *thePropagationMomentum * (thePosition.dot(*thePropagationMomentum)/thePropagationMomentum->mag2());
00621 }
00622
00624 const ThreeVector &adjustMomentumFromEnergy();
00625
00627 G4double adjustEnergyFromMomentum();
00628
00630 G4bool isInList(ParticleList const &l) const {
00631 return (std::find(l.begin(), l.end(), this)!=l.end());
00632 }
00633
00634 G4bool isCluster() const {
00635 return (theType == Composite);
00636 }
00637
00639 void setFrozenMomentum(const ThreeVector &momentum) { theFrozenMomentum = momentum; }
00640
00642 void setFrozenEnergy(const G4double energy) { theFrozenEnergy = energy; }
00643
00645 ThreeVector getFrozenMomentum() const { return theFrozenMomentum; }
00646
00648 G4double getFrozenEnergy() const { return theFrozenEnergy; }
00649
00651 ThreeVector getPropagationVelocity() const { return (*thePropagationMomentum)/(*thePropagationEnergy); }
00652
00659 void freezePropagation() {
00660 thePropagationMomentum = &theFrozenMomentum;
00661 thePropagationEnergy = &theFrozenEnergy;
00662 }
00663
00670 void thawPropagation() {
00671 thePropagationMomentum = &theMomentum;
00672 thePropagationEnergy = &theEnergy;
00673 }
00674
00680 virtual void rotate(const G4double angle, const ThreeVector &axis) {
00681 thePosition.rotate(angle, axis);
00682 theMomentum.rotate(angle, axis);
00683 theFrozenMomentum.rotate(angle, axis);
00684 }
00685
00686 std::string print() const {
00687 std::stringstream ss;
00688 ss << "Particle (ID = " << ID << ") type = ";
00689 ss << ParticleTable::getName(theType);
00690 ss << std::endl
00691 << " energy = " << theEnergy << std::endl
00692 << " momentum = "
00693 << theMomentum.print()
00694 << std::endl
00695 << " position = "
00696 << thePosition.print()
00697 << std::endl;
00698 return ss.str();
00699 };
00700
00701 std::string dump() const {
00702 std::stringstream ss;
00703 ss << "(particle " << ID << " ";
00704 ss << ParticleTable::getName(theType);
00705 ss << std::endl
00706 << thePosition.dump()
00707 << std::endl
00708 << theMomentum.dump()
00709 << std::endl
00710 << theEnergy << ")" << std::endl;
00711 return ss.str();
00712 };
00713
00714 long getID() const { return ID; };
00715
00719 ParticleList const *getParticles() const {
00720 WARN("Particle::getParticles() method was called on a Particle object" << std::endl);
00721 return 0;
00722 }
00723
00724 protected:
00725 G4int theZ, theA;
00726 ParticipantType theParticipantType;
00727 G4INCL::ParticleType theType;
00728 G4double theEnergy;
00729 G4double *thePropagationEnergy;
00730 G4double theFrozenEnergy;
00731 G4INCL::ThreeVector theMomentum;
00732 G4INCL::ThreeVector *thePropagationMomentum;
00733 G4INCL::ThreeVector theFrozenMomentum;
00734 G4INCL::ThreeVector thePosition;
00735 G4int nCollisions;
00736 G4int nDecays;
00737 G4double thePotentialEnergy;
00738 long ID;
00739
00740 private:
00741 G4double theHelicity;
00742 G4double emissionTime;
00743 G4bool outOfWell;
00744
00745 G4double theMass;
00746 static long nextID;
00747
00748 };
00749 }
00750
00751 #endif