include/Geant4/G4KineticTrack.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 //
0028 // -----------------------------------------------------------------------------
0029 //      GEANT 4 class header file
0030 //
0031 //      History: first implementation, A. Feliciello, 20th May 1998
0032 // -----------------------------------------------------------------------------
0033
0034 #ifndef G4KineticTrack_h
0035 #define G4KineticTrack_h 1
0036
0037 #include <CLHEP/Units/PhysicalConstants.h>
0038
0039 #include "globals.hh"
0040 #include "G4ios.hh"
0041
0042
0043 #include "Randomize.hh"
0044 #include "G4ThreeVector.hh"
0045 #include "G4LorentzVector.hh"
0046 #include "G4VKineticNucleon.hh"
0047 #include "G4Nucleon.hh"
0048 #include "G4ParticleDefinition.hh"
0049 #include "G4VDecayChannel.hh"
0050 #include "G4Log.hh"
0051
0052 // #include "G4Allocator.hh"
0053
0054 class G4KineticTrackVector;
0055
0056 class G4KineticTrack : public G4VKineticNucleon
0057 {
0058   public:
0059
0060       G4KineticTrack();
0061
0062       G4KineticTrack(const G4KineticTrack& right);
0063
0064       G4KineticTrack(const G4ParticleDefinition* aDefinition,
0065                      G4double aFormationTime,
0066                      const G4ThreeVector& aPosition,
0067                      const G4LorentzVector& a4Momentum);
0068       G4KineticTrack(G4Nucleon * nucleon,
0069                      const G4ThreeVector& aPosition,
0070                      const G4LorentzVector& a4Momentum);
0071
0072       ~G4KineticTrack();
0073
0074       G4KineticTrack& operator=(const G4KineticTrack& right);
0075
0076       G4bool operator==(const G4KineticTrack& right) const;
0077
0078       G4bool operator!=(const G4KineticTrack& right) const;
0079 /*
0080       inline void *operator new(size_t);
0081       inline void operator delete(void *aTrack);
0082 */
0083       const G4ParticleDefinition* GetDefinition() const;
0084       void SetDefinition(const G4ParticleDefinition* aDefinition);
0085
0086       G4double GetFormationTime() const;
0087       void SetFormationTime(G4double aFormationTime);
0088
0089       const G4ThreeVector& GetPosition() const;
0090       void SetPosition(const G4ThreeVector aPosition);
0091
0092       const G4LorentzVector& Get4Momentum() const;
0093       void Set4Momentum(const G4LorentzVector& a4Momentum);
0094       void Update4Momentum(G4double aEnergy);           // update E and p, not changing mass
0095       void Update4Momentum(const G4ThreeVector & aMomentum);    // idem
0096       void SetTrackingMomentum(const G4LorentzVector& a4Momentum);
0097       void UpdateTrackingMomentum(G4double aEnergy);            // update E and p, not changing mass
0098       void UpdateTrackingMomentum(const G4ThreeVector & aMomentum); // idem
0099
0100       const G4LorentzVector& GetTrackingMomentum() const;
0101
0102       G4double SampleResidualLifetime();
0103
0104       void Hit();
0105       void SetNucleon(G4Nucleon * aN) {theNucleon = aN;}
0106
0107       G4bool IsParticipant() const;
0108
0109       G4KineticTrackVector* Decay();
0110
0111   // LB move to public (before was private) LB
0112       G4double* GetActualWidth() const;
0113
0114       G4double GetActualMass() const;
0115       G4int GetnChannels() const;
0116
0117 //   position relativ to nucleus "state"
0118       enum CascadeState {undefined, outside, going_in, inside,
0119                          going_out, gone_out, captured, miss_nucleus };
0120
0121       CascadeState SetState(const CascadeState new_state);
0122       CascadeState GetState() const;
0123       void SetProjectilePotential(const G4double aPotential);
0124       G4double GetProjectilePotential() const;
0125
0126       void SetCreatorModelID(G4int id);
0127       G4int GetCreatorModelID() const;
0128
0129       const G4ParticleDefinition* GetParentResonanceDef() const;
0130       void SetParentResonanceDef(const G4ParticleDefinition* parent);
0131       G4int GetParentResonanceID() const;
0132       void SetParentResonanceID(const G4int parentID);
0133
0134   private:
0135
0136
0137       void SetnChannels(const G4int aChannel);
0138
0139       void SetActualWidth(G4double* anActualWidth);
0140
0141       G4double EvaluateTotalActualWidth();
0142
0143       G4double EvaluateCMMomentum (const G4double mass,
0144                                    const G4double* m_ij) const;
0145
0146       G4double IntegrateCMMomentum(const G4double lowerLimit) const;
0147
0148       G4double IntegrateCMMomentum(const G4double lowerLimit ,const G4double polemass) const;
0149
0150       G4double IntegrateCMMomentum2() const;
0151
0152   public:
0153
0154       G4double BrWig(const G4double Gamma,
0155                      const G4double rmass,
0156                      const G4double mass) const;
0157
0158 private:
0159       G4double IntegrandFunction1 (G4double xmass) const;
0160       G4double IntegrandFunction2 (G4double xmass) const;
0161       G4double IntegrandFunction3 (G4double xmass) const;
0162       G4double IntegrandFunction4 (G4double xmass) const;
0163 public:
0164   //   friend G4double IntegrandFunction3 (G4double xmass);
0165
0166   //   friend G4double IntegrandFunction4 (G4double xmass);
0167
0168   private:
0169
0170       const G4ParticleDefinition* theDefinition;
0171
0172       G4double theFormationTime;
0173
0174       G4ThreeVector thePosition;
0175
0176       G4LorentzVector the4Momentum;
0177       G4LorentzVector theFermi3Momentum;
0178       G4LorentzVector theTotal4Momentum;
0179
0180       G4Nucleon * theNucleon;
0181
0182       G4int nChannels;
0183
0184       G4double theActualMass;
0185
0186       G4double* theActualWidth;
0187
0188      // Temporary storage for daughter masses and widths
0189       // (needed because Integrand Function cannot take > 1 argument)
0190       G4double* theDaughterMass;
0191       G4double* theDaughterWidth;
0192
0193       CascadeState theStateToNucleus;
0194
0195       G4double theProjectilePotential;
0196
0197       G4int theCreatorModel;
0198
0199       const G4ParticleDefinition* theParentResonanceDef = nullptr;
0200       G4int theParentResonanceID;
0201 };
0202
0203 // extern G4Allocator<G4KineticTrack> theKTAllocator;
0204
0205
0206 // Class G4KineticTrack
0207 /*
0208 inline void * G4KineticTrack::operator new(size_t)
0209 {
0210   void * aT;
0211   aT = (void *) theKTAllocator.MallocSingle();
0212   return aT;
0213 }
0214
0215 inline void G4KineticTrack::operator delete(void * aT)
0216 {
0217   theKTAllocator.FreeSingle((G4KineticTrack *) aT);
0218 }
0219 */
0220
0221 inline const G4ParticleDefinition* G4KineticTrack::GetDefinition() const
0222 {
0223   return theDefinition;
0224 }
0225
0226 inline void G4KineticTrack::SetDefinition(const G4ParticleDefinition* aDefinition)
0227 {
0228   theDefinition = aDefinition;
0229 }
0230
0231 inline G4double G4KineticTrack::GetFormationTime() const
0232 {
0233   return theFormationTime;
0234 }
0235
0236 inline void G4KineticTrack::SetFormationTime(G4double aFormationTime)
0237 {
0238   theFormationTime = aFormationTime;
0239 }
0240
0241 inline const G4ThreeVector& G4KineticTrack::GetPosition() const
0242 {
0243   return thePosition;
0244 }
0245
0246 inline void G4KineticTrack::SetPosition(const G4ThreeVector aPosition)
0247 {
0248   thePosition = aPosition;
0249 }
0250
0251 inline const G4LorentzVector& G4KineticTrack::Get4Momentum() const
0252 {
0253   return theTotal4Momentum;
0254 }
0255
0256 inline const G4LorentzVector& G4KineticTrack::GetTrackingMomentum() const
0257 {
0258    return the4Momentum;
0259 }
0260
0261 inline void G4KineticTrack::Set4Momentum(const G4LorentzVector& a4Momentum)
0262 {
0263 //  set the4Momentum and update theTotal4Momentum
0264
0265   theTotal4Momentum=a4Momentum;
0266   the4Momentum = theTotal4Momentum;
0267   theFermi3Momentum=G4LorentzVector(0);
0268 }
0269
0270 inline void G4KineticTrack::Update4Momentum(G4double aEnergy)
0271 {
0272 // update the4Momentum with aEnergy at constant mass (the4Momentum.mag()
0273 //   updates theTotal4Momentum as well.
0274   G4double newP(0);
0275   G4double mass2=theTotal4Momentum.mag2();
0276   if ( sqr(aEnergy) > mass2 )
0277   {
0278       newP = std::sqrt(sqr(aEnergy) - mass2 );
0279   } else
0280   {
0281       aEnergy=std::sqrt(mass2);
0282   }
0283   Set4Momentum(G4LorentzVector(newP*the4Momentum.vect().unit(), aEnergy));
0284 }
0285
0286 inline void G4KineticTrack::Update4Momentum(const G4ThreeVector & aMomentum)
0287 {
0288 // update the4Momentum with aMomentum at constant mass (the4Momentum.mag()
0289 //   updates theTotal4Momentum as well.
0290   G4double newE=std::sqrt(theTotal4Momentum.mag2() + aMomentum.mag2());
0291   Set4Momentum(G4LorentzVector(aMomentum, newE));
0292 }
0293
0294 inline void G4KineticTrack::SetTrackingMomentum(const G4LorentzVector& aMomentum)
0295 {
0296 //  set the4Momentum and update theTotal4Momentum, keep the mass of aMomentum
0297
0298   the4Momentum = aMomentum;
0299   theTotal4Momentum=the4Momentum+theFermi3Momentum;
0300 //     keep mass of aMomentum for the total momentum
0301   G4double mass2 = aMomentum.mag2();
0302   G4double p2=theTotal4Momentum.vect().mag2();
0303   theTotal4Momentum.setE(std::sqrt(mass2+p2));
0304 }
0305
0306 inline void G4KineticTrack::UpdateTrackingMomentum(G4double aEnergy)
0307 {
0308 // update the4Momentum with aEnergy at constant mass (the4Momentum.mag()
0309 //   updates theTotal4Momentum as well.
0310   G4double newP(0);
0311   G4double mass2=theTotal4Momentum.mag2();
0312   if ( sqr(aEnergy) > mass2 )
0313   {
0314       newP = std::sqrt(sqr(aEnergy) - mass2 );
0315   } else
0316   {
0317       aEnergy=std::sqrt(mass2);
0318   }
0319   SetTrackingMomentum(G4LorentzVector(newP*the4Momentum.vect().unit(), aEnergy));
0320 }
0321
0322 inline void G4KineticTrack::UpdateTrackingMomentum(const G4ThreeVector & aMomentum)
0323 {
0324 // update the4Momentum with aMomentum at constant mass (the4Momentum.mag()
0325 //   updates theTotal4Momentum as well.
0326   G4double newE=std::sqrt(theTotal4Momentum.mag2() + aMomentum.mag2());
0327   SetTrackingMomentum(G4LorentzVector(aMomentum, newE));
0328 }
0329
0330 inline G4double G4KineticTrack::GetActualMass() const
0331 {
0332   return std::sqrt(std::abs(the4Momentum.mag2()));
0333 }
0334
0335 inline G4int G4KineticTrack::GetnChannels() const
0336 {
0337   return nChannels;
0338 }
0339
0340 inline void G4KineticTrack::SetnChannels(const G4int numberOfChannels)
0341 {
0342   nChannels = numberOfChannels;
0343 }
0344
0345 inline G4double* G4KineticTrack::GetActualWidth() const
0346 {
0347   return theActualWidth;
0348 }
0349
0350 inline void G4KineticTrack::SetActualWidth(G4double* anActualWidth)
0351 {
0352   theActualWidth = anActualWidth;
0353 }
0354
0355
0356
0357 inline G4double G4KineticTrack::EvaluateTotalActualWidth()
0358 {
0359  G4int index;
0360  G4double theTotalActualWidth = 0.0;
0361  for (index = nChannels - 1; index >= 0; index--)
0362     {
0363      theTotalActualWidth += theActualWidth[index];
0364     }
0365  return theTotalActualWidth;
0366 }
0367
0368 inline G4double G4KineticTrack::SampleResidualLifetime()
0369 {
0370  G4double theTotalActualWidth = this->EvaluateTotalActualWidth();
0371  G4double tau = CLHEP::hbar_Planck * (-1.0 / theTotalActualWidth);
0372  G4double theResidualLifetime = tau * G4Log(G4UniformRand());
0373  return theResidualLifetime*the4Momentum.gamma();
0374 }
0375
0376 inline G4double G4KineticTrack::EvaluateCMMomentum(const G4double mass,
0377                                                  const G4double* m_ij) const
0378 {
0379   G4double theCMMomentum;
0380   if((m_ij[0]+m_ij[1])<mass)
0381    theCMMomentum = 1 / (2 * mass) *
0382           std::sqrt (((mass * mass) - (m_ij[0] + m_ij[1]) * (m_ij[0] + m_ij[1])) *
0383                 ((mass * mass) - (m_ij[0] - m_ij[1]) * (m_ij[0] - m_ij[1])));
0384   else
0385    theCMMomentum=0.;
0386
0387  return theCMMomentum;
0388 }
0389
0390 inline G4double G4KineticTrack::BrWig(const G4double Gamma, const G4double rmass, const G4double mass) const
0391 {
0392   G4double Norm = CLHEP::twopi;
0393   return (Gamma/((mass-rmass)*(mass-rmass)+Gamma*Gamma/4.))/Norm;
0394 }
0395
0396 inline
0397 void G4KineticTrack::Hit()
0398 {
0399   if(theNucleon)
0400   {
0401     theNucleon->Hit(1);
0402   }
0403 }
0404
0405 inline
0406 G4bool G4KineticTrack::IsParticipant() const
0407 {
0408   if(!theNucleon) return true;
0409   return theNucleon->AreYouHit();
0410 }
0411
0412 inline
0413 G4KineticTrack::CascadeState G4KineticTrack::GetState() const
0414 {
0415     return theStateToNucleus;
0416 }
0417
0418 inline
0419 G4KineticTrack::CascadeState G4KineticTrack::SetState(const CascadeState new_state)
0420 {
0421     CascadeState old_state=theStateToNucleus;
0422     theStateToNucleus=new_state;
0423     return old_state;
0424 }
0425
0426 inline
0427 void G4KineticTrack::SetProjectilePotential(G4double aPotential)
0428 {
0429     theProjectilePotential = aPotential;
0430 }
0431 inline
0432 G4double G4KineticTrack::GetProjectilePotential() const
0433 {
0434     return theProjectilePotential;
0435 }
0436
0437 inline
0438 void G4KineticTrack::SetCreatorModelID(G4int id)
0439 {
0440         theCreatorModel = id;
0441 }
0442 inline
0443 G4int G4KineticTrack::GetCreatorModelID() const
0444 {
0445        return theCreatorModel;
0446 }
0447
0448 inline
0449 const G4ParticleDefinition* G4KineticTrack::GetParentResonanceDef() const
0450 {
0451         return theParentResonanceDef;
0452 }
0453
0454 inline
0455 void G4KineticTrack::SetParentResonanceDef(const G4ParticleDefinition* parentDef)
0456 {
0457         theParentResonanceDef = parentDef;
0458 }
0459
0460 inline
0461 G4int G4KineticTrack::GetParentResonanceID() const
0462 {
0463         return theParentResonanceID;
0464 }
0465
0466 inline
0467 void G4KineticTrack::SetParentResonanceID(const G4int parentID)
0468 {
0469         theParentResonanceID = parentID;
0470 }
0471
0472 #endif