Geant4  10.01.p02
G4BohrFluctuations.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4BohrFluctuations.cc 72048 2013-07-04 12:39:58Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4BohrFluctuations
34 //
35 // Author: Vladimir Ivanchenko
36 //
37 // Creation date: 02.04.2003
38 //
39 // Modifications:
40 //
41 // 23-05-03 Add control on parthalogical cases (V.Ivanchenko)
42 // 16-10-03 Changed interface to Initialisation (V.Ivanchenko)
43 //
44 // Class Description: Sampling of Gaussion fluctuations
45 //
46 // -------------------------------------------------------------------
47 //
48 
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51 
52 #include "G4BohrFluctuations.hh"
53 #include "G4PhysicalConstants.hh"
54 #include "G4SystemOfUnits.hh"
55 #include "Randomize.hh"
56 #include "G4Poisson.hh"
57 #include "G4ParticleDefinition.hh"
58 #include "G4MaterialCutsCouple.hh"
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
61 
62 using namespace std;
63 
66  particle(0),
67  minNumberInteractionsBohr(2.0),
68  minFraction(0.2),
69  xmin(0.2),
70  minLoss(0.001*eV)
71 {
72  particleMass = proton_mass_c2;
73  chargeSquare = 1.0;
74  kineticEnergy = 0.0;
75  beta2 = 0.0;
76 }
77 
78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
79 
81 {}
82 
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84 
86 {
87  particle = part;
88  particleMass = part->GetPDGMass();
89  G4double q = part->GetPDGCharge()/eplus;
90  chargeSquare = q*q;
91 }
92 
93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
94 
95 G4double
97  const G4DynamicParticle* dp,
98  G4double tmax,
99  G4double length,
100  G4double meanLoss)
101 {
102  if(meanLoss <= minLoss) { return meanLoss; }
103  const G4Material* material = couple->GetMaterial();
104  G4double siga = Dispersion(material,dp,tmax,length);
105  G4double loss = meanLoss;
106 
107  G4double navr = meanLoss*meanLoss/siga;
108  //G4cout << "### meanLoss= " << meanLoss << " navr= " << navr << G4endl;
109  if (navr >= minNumberInteractionsBohr) {
110 
111  // Increase fluctuations for big fractional energy loss
112  if ( meanLoss > minFraction*kineticEnergy ) {
113  G4double gam = (kineticEnergy - meanLoss)/particleMass + 1.0;
114  G4double b2 = 1.0 - 1.0/(gam*gam);
115  if(b2 < xmin*beta2) b2 = xmin*beta2;
116  G4double x = b2/beta2;
117  G4double x3 = 1.0/(x*x*x);
118  siga *= 0.25*(1.0 + x)*(x3 + (1.0/b2 - 0.5)/(1.0/beta2 - 0.5) );
119  }
120  siga = sqrt(siga);
121  G4double twomeanLoss = meanLoss + meanLoss;
122  //G4cout << "siga= " << siga << " 2edp= " << twomeanLoss <<G4endl;
123 
124  if(twomeanLoss < siga) {
125  G4double x;
126  do {
127  loss = twomeanLoss*G4UniformRand();
128  x = (loss - meanLoss)/siga;
129  } while (1.0 - 0.5*x*x < G4UniformRand());
130  } else {
131  do {
132  loss = G4RandGauss::shoot(meanLoss,siga);
133  } while (0.0 > loss || loss > twomeanLoss);
134  }
135 
136  // Poisson fluctuations
137  } else {
138  G4double n = (G4double)(G4Poisson(navr));
139  loss = meanLoss*n/navr;
140  }
141  //G4cout << "loss= " << loss << G4endl;
142 
143  return loss;
144 }
145 
146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
147 
149  const G4DynamicParticle* dp,
150  G4double tmax,
151  G4double length)
152 {
153  if(!particle) { InitialiseMe(dp->GetDefinition()); }
154 
155  G4double electronDensity = material->GetElectronDensity();
158  beta2 = kineticEnergy*(kineticEnergy + 2.0*particleMass)/(etot*etot);
159  G4double siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
160  * electronDensity * chargeSquare;
161 
162  return siga;
163 }
164 
165 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
166 
167 
ThreeVector shoot(const G4int Ap, const G4int Af)
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
G4double GetKineticEnergy() const
G4double Dispersion(const G4Material *, const G4DynamicParticle *, G4double, G4double)
const G4ParticleDefinition * particle
G4ParticleDefinition * GetDefinition() const
G4double minNumberInteractionsBohr
#define G4UniformRand()
Definition: Randomize.hh:93
G4double GetElectronDensity() const
Definition: G4Material.hh:217
static const G4double b2
const G4int n
G4BohrFluctuations(const G4String &nam="BohrFluc")
static const double eV
Definition: G4SIunits.hh:194
G4double GetPDGMass() const
void InitialiseMe(const G4ParticleDefinition *)
double G4double
Definition: G4Types.hh:76
static const double eplus
Definition: G4SIunits.hh:178
G4double GetPDGCharge() const
const G4Material * GetMaterial() const
G4double SampleFluctuations(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double, G4double, G4double)