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 #ifndef G4KineticTrack_h
00036 #define G4KineticTrack_h 1
00037
00038 #include <CLHEP/Units/PhysicalConstants.h>
00039
00040 #include "globals.hh"
00041 #include "G4ios.hh"
00042
00043
00044 #include "Randomize.hh"
00045 #include "G4ThreeVector.hh"
00046 #include "G4LorentzVector.hh"
00047 #include "G4VKineticNucleon.hh"
00048 #include "G4Nucleon.hh"
00049 #include "G4ParticleDefinition.hh"
00050 #include "G4VDecayChannel.hh"
00051
00052
00053
00054 class G4KineticTrackVector;
00055
00056
00057
00058
00059
00060 class G4KineticTrack : public G4VKineticNucleon
00061 {
00062 public:
00063
00064 G4KineticTrack();
00065
00066 G4KineticTrack(const G4KineticTrack& right);
00067
00068 G4KineticTrack(G4ParticleDefinition* aDefinition,
00069 G4double aFormationTime,
00070 G4ThreeVector aPosition,
00071 G4LorentzVector& a4Momentum);
00072 G4KineticTrack(G4Nucleon * nucleon,
00073 G4ThreeVector aPosition,
00074 G4LorentzVector& a4Momentum);
00075
00076 ~G4KineticTrack();
00077
00078 G4KineticTrack& operator=(const G4KineticTrack& right);
00079
00080 G4int operator==(const G4KineticTrack& right) const;
00081
00082 G4int operator!=(const G4KineticTrack& right) const;
00083
00084
00085
00086
00087 G4ParticleDefinition* GetDefinition() const;
00088 void SetDefinition(G4ParticleDefinition* aDefinition);
00089
00090 G4double GetFormationTime() const;
00091 void SetFormationTime(G4double aFormationTime);
00092
00093 const G4ThreeVector& GetPosition() const;
00094 void SetPosition(const G4ThreeVector aPosition);
00095
00096 const G4LorentzVector& Get4Momentum() const;
00097 void Set4Momentum(const G4LorentzVector& a4Momentum);
00098 void Update4Momentum(G4double aEnergy);
00099 void Update4Momentum(const G4ThreeVector & aMomentum);
00100 void SetTrackingMomentum(const G4LorentzVector& a4Momentum);
00101 void UpdateTrackingMomentum(G4double aEnergy);
00102 void UpdateTrackingMomentum(const G4ThreeVector & aMomentum);
00103
00104 const G4LorentzVector& GetTrackingMomentum() const;
00105
00106 G4double SampleResidualLifetime();
00107
00108 void Hit();
00109 void SetNucleon(G4Nucleon * aN) {theNucleon = aN;}
00110
00111 G4bool IsParticipant() const;
00112
00113 G4KineticTrackVector* Decay();
00114
00115
00116 G4double* GetActualWidth() const;
00117
00118 G4double GetActualMass() const;
00119 G4int GetnChannels() const;
00120
00121
00122 enum CascadeState {undefined, outside, going_in, inside,
00123 going_out, gone_out, captured, miss_nucleus };
00124
00125 CascadeState SetState(const CascadeState new_state);
00126 CascadeState GetState() const;
00127 void SetProjectilePotential(const G4double aPotential);
00128 G4double GetProjectilePotential() const;
00129
00130
00131 private:
00132
00133
00134 void SetnChannels(const G4int aChannel);
00135
00136 void SetActualWidth(G4double* anActualWidth);
00137
00138 G4double EvaluateTotalActualWidth();
00139
00140 G4double EvaluateCMMomentum (const G4double mass,
00141 const G4double* m_ij) const;
00142
00143 G4double IntegrateCMMomentum(const G4double lowerLimit) const;
00144
00145 G4double IntegrateCMMomentum(const G4double lowerLimit ,const G4double polemass) const;
00146
00147 G4double IntegrateCMMomentum2() const;
00148
00149 public:
00150
00151 G4double BrWig(const G4double Gamma,
00152 const G4double rmass,
00153 const G4double mass) const;
00154
00155 private:
00156 G4double IntegrandFunction1 (G4double xmass) const;
00157 G4double IntegrandFunction2 (G4double xmass) const;
00158 G4double IntegrandFunction3 (G4double xmass) const;
00159 G4double IntegrandFunction4 (G4double xmass) const;
00160 public:
00161
00162
00163
00164
00165 private:
00166
00167 G4ParticleDefinition* theDefinition;
00168
00169 G4double theFormationTime;
00170
00171 G4ThreeVector thePosition;
00172
00173 G4LorentzVector the4Momentum;
00174 G4LorentzVector theFermi3Momentum;
00175 G4LorentzVector theTotal4Momentum;
00176
00177 G4Nucleon * theNucleon;
00178
00179 G4int nChannels;
00180
00181 G4double theActualMass;
00182
00183 G4double* theActualWidth;
00184
00185
00186
00187 G4double* theDaughterMass;
00188 G4double* theDaughterWidth;
00189
00190 CascadeState theStateToNucleus;
00191
00192 G4double theProjectilePotential;
00193 };
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213 inline G4ParticleDefinition* G4KineticTrack::GetDefinition() const
00214 {
00215 return theDefinition;
00216 }
00217
00218 inline void G4KineticTrack::SetDefinition(G4ParticleDefinition* aDefinition)
00219 {
00220 theDefinition = aDefinition;
00221 }
00222
00223
00224
00225 inline G4double G4KineticTrack::GetFormationTime() const
00226 {
00227 return theFormationTime;
00228 }
00229
00230 inline void G4KineticTrack::SetFormationTime(G4double aFormationTime)
00231 {
00232 theFormationTime = aFormationTime;
00233 }
00234
00235
00236
00237 inline const G4ThreeVector& G4KineticTrack::GetPosition() const
00238 {
00239 return thePosition;
00240 }
00241
00242 inline void G4KineticTrack::SetPosition(const G4ThreeVector aPosition)
00243 {
00244 thePosition = aPosition;
00245 }
00246
00247
00248 inline const G4LorentzVector& G4KineticTrack::Get4Momentum() const
00249 {
00250 return theTotal4Momentum;
00251 }
00252
00253 inline const G4LorentzVector& G4KineticTrack::GetTrackingMomentum() const
00254 {
00255 return the4Momentum;
00256 }
00257
00258 inline void G4KineticTrack::Set4Momentum(const G4LorentzVector& a4Momentum)
00259 {
00260
00261
00262 theTotal4Momentum=a4Momentum;
00263 the4Momentum = theTotal4Momentum;
00264 theFermi3Momentum=G4LorentzVector(0);
00265 }
00266
00267 inline void G4KineticTrack::Update4Momentum(G4double aEnergy)
00268 {
00269
00270
00271 G4double newP(0);
00272 G4double mass2=theTotal4Momentum.mag2();
00273 if ( sqr(aEnergy) > mass2 )
00274 {
00275 newP = std::sqrt(sqr(aEnergy) - mass2 );
00276 } else
00277 {
00278 aEnergy=std::sqrt(mass2);
00279 }
00280 Set4Momentum(G4LorentzVector(newP*the4Momentum.vect().unit(), aEnergy));
00281 }
00282
00283 inline void G4KineticTrack::Update4Momentum(const G4ThreeVector & aMomentum)
00284 {
00285
00286
00287 G4double newE=std::sqrt(theTotal4Momentum.mag2() + aMomentum.mag2());
00288 Set4Momentum(G4LorentzVector(aMomentum, newE));
00289 }
00290
00291 inline void G4KineticTrack::SetTrackingMomentum(const G4LorentzVector& aMomentum)
00292 {
00293
00294
00295 the4Momentum = aMomentum;
00296 theTotal4Momentum=the4Momentum+theFermi3Momentum;
00297
00298 G4double mass2 = aMomentum.mag2();
00299 G4double p2=theTotal4Momentum.vect().mag2();
00300 theTotal4Momentum.setE(std::sqrt(mass2+p2));
00301 }
00302
00303 inline void G4KineticTrack::UpdateTrackingMomentum(G4double aEnergy)
00304 {
00305
00306
00307 G4double newP(0);
00308 G4double mass2=theTotal4Momentum.mag2();
00309 if ( sqr(aEnergy) > mass2 )
00310 {
00311 newP = std::sqrt(sqr(aEnergy) - mass2 );
00312 } else
00313 {
00314 aEnergy=std::sqrt(mass2);
00315 }
00316 SetTrackingMomentum(G4LorentzVector(newP*the4Momentum.vect().unit(), aEnergy));
00317 }
00318
00319 inline void G4KineticTrack::UpdateTrackingMomentum(const G4ThreeVector & aMomentum)
00320 {
00321
00322
00323 G4double newE=std::sqrt(theTotal4Momentum.mag2() + aMomentum.mag2());
00324 SetTrackingMomentum(G4LorentzVector(aMomentum, newE));
00325 }
00326
00327
00328
00329
00330 inline G4double G4KineticTrack::GetActualMass() const
00331 {
00332 return std::sqrt(std::abs(the4Momentum.mag2()));
00333 }
00334
00335
00336
00337 inline G4int G4KineticTrack::GetnChannels() const
00338 {
00339 return nChannels;
00340 }
00341
00342 inline void G4KineticTrack::SetnChannels(const G4int numberOfChannels)
00343 {
00344 nChannels = numberOfChannels;
00345 }
00346
00347
00348
00349 inline G4double* G4KineticTrack::GetActualWidth() const
00350 {
00351 return theActualWidth;
00352 }
00353
00354 inline void G4KineticTrack::SetActualWidth(G4double* anActualWidth)
00355 {
00356 theActualWidth = anActualWidth;
00357 }
00358
00359
00360
00361 inline G4double G4KineticTrack::EvaluateTotalActualWidth()
00362 {
00363 G4int index;
00364 G4double theTotalActualWidth = 0.0;
00365 for (index = nChannels - 1; index >= 0; index--)
00366 {
00367 theTotalActualWidth += theActualWidth[index];
00368 }
00369 return theTotalActualWidth;
00370 }
00371
00372
00373
00374 inline G4double G4KineticTrack::SampleResidualLifetime()
00375 {
00376 G4double theTotalActualWidth = this->EvaluateTotalActualWidth();
00377 G4double tau = CLHEP::hbar_Planck * (-1.0 / theTotalActualWidth);
00378 G4double theResidualLifetime = tau * std::log(G4UniformRand());
00379 return theResidualLifetime*the4Momentum.gamma();
00380 }
00381
00382
00383
00384 inline G4double G4KineticTrack::EvaluateCMMomentum(const G4double mass,
00385 const G4double* m_ij) const
00386 {
00387 G4double theCMMomentum;
00388 if((m_ij[0]+m_ij[1])<mass)
00389 theCMMomentum = 1 / (2 * mass) *
00390 std::sqrt (((mass * mass) - (m_ij[0] + m_ij[1]) * (m_ij[0] + m_ij[1])) *
00391 ((mass * mass) - (m_ij[0] - m_ij[1]) * (m_ij[0] - m_ij[1])));
00392 else
00393 theCMMomentum=0.;
00394
00395 return theCMMomentum;
00396 }
00397
00398 inline G4double G4KineticTrack::BrWig(const G4double Gamma, const G4double rmass, const G4double mass) const
00399 {
00400 G4double Norm = CLHEP::twopi;
00401 return (Gamma/((mass-rmass)*(mass-rmass)+Gamma*Gamma/4.))/Norm;
00402 }
00403
00404 inline
00405 void G4KineticTrack::Hit()
00406 {
00407 if(theNucleon)
00408 {
00409 theNucleon->Hit(1);
00410 }
00411 }
00412
00413 inline
00414 G4bool G4KineticTrack::IsParticipant() const
00415 {
00416 if(!theNucleon) return true;
00417 return theNucleon->AreYouHit();
00418 }
00419
00420 inline
00421 G4KineticTrack::CascadeState G4KineticTrack::GetState() const
00422 {
00423 return theStateToNucleus;
00424 }
00425
00426 inline
00427 G4KineticTrack::CascadeState G4KineticTrack::SetState(const CascadeState new_state)
00428 {
00429 CascadeState old_state=theStateToNucleus;
00430 theStateToNucleus=new_state;
00431 return old_state;
00432 }
00433
00434 inline
00435 void G4KineticTrack::SetProjectilePotential(G4double aPotential)
00436 {
00437 theProjectilePotential = aPotential;
00438 }
00439 inline
00440 G4double G4KineticTrack::GetProjectilePotential() const
00441 {
00442 return theProjectilePotential;
00443 }
00444
00445 #endif
00446
00447
00448