source: trunk/source/processes/electromagnetic/standard/src/G4BohrFluctuations.cc@ 1346

Last change on this file since 1346 was 1340, checked in by garnier, 15 years ago

update ti head
File size: 5.7 KB
Line 
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,v 1.9 2010/10/25 18:23:36 vnivanch Exp $
27// GEANT4 tag $Name: emstand-V09-03-24 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name: G4BohrFluctuations
35//
36// Author: Vladimir Ivanchenko
37//
38// Creation date: 02.04.2003
39//
40// Modifications:
41//
42// 23-05-03 Add control on parthalogical cases (V.Ivanchenko)
43// 16-10-03 Changed interface to Initialisation (V.Ivanchenko)
44//
45// Class Description: Sampling of Gaussion fluctuations
46//
47// -------------------------------------------------------------------
48//
49
50//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
52
53#include "G4BohrFluctuations.hh"
54#include "Randomize.hh"
55#include "G4Poisson.hh"
56#include "G4ParticleDefinition.hh"
57
58//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
59
60using namespace std;
61
62G4BohrFluctuations::G4BohrFluctuations(const G4String& nam)
63 :G4VEmFluctuationModel(nam),
64 particle(0),
65 minNumberInteractionsBohr(2.0),
66 minFraction(0.2),
67 xmin(0.2),
68 minLoss(0.001*eV)
69{
70 particleMass = proton_mass_c2;
71 chargeSquare = 1.0;
72 kineticEnergy = 0.0;
73 beta2 = 0.0;
74}
75
76//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77
78G4BohrFluctuations::~G4BohrFluctuations()
79{}
80
81//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
82
83void G4BohrFluctuations::InitialiseMe(const G4ParticleDefinition* part)
84{
85 particle = part;
86 particleMass = part->GetPDGMass();
87 G4double q = part->GetPDGCharge()/eplus;
88 chargeSquare = q*q;
89}
90
91//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92
93G4double G4BohrFluctuations::SampleFluctuations(const G4Material* material,
94 const G4DynamicParticle* dp,
95 G4double& tmax,
96 G4double& length,
97 G4double& meanLoss)
98{
99 if(meanLoss <= minLoss) { return meanLoss; }
100 G4double siga = Dispersion(material,dp,tmax,length);
101 G4double loss = meanLoss;
102
103 G4double navr = meanLoss*meanLoss/siga;
104 //G4cout << "### meanLoss= " << meanLoss << " navr= " << navr << G4endl;
105 if (navr >= minNumberInteractionsBohr) {
106
107 // Increase fluctuations for big fractional energy loss
108 if ( meanLoss > minFraction*kineticEnergy ) {
109 G4double gam = (kineticEnergy - meanLoss)/particleMass + 1.0;
110 G4double b2 = 1.0 - 1.0/(gam*gam);
111 if(b2 < xmin*beta2) b2 = xmin*beta2;
112 G4double x = b2/beta2;
113 G4double x3 = 1.0/(x*x*x);
114 siga *= 0.25*(1.0 + x)*(x3 + (1.0/b2 - 0.5)/(1.0/beta2 - 0.5) );
115 }
116 siga = sqrt(siga);
117 G4double twomeanLoss = meanLoss + meanLoss;
118 //G4cout << "siga= " << siga << " 2edp= " << twomeanLoss <<G4endl;
119
120 if(twomeanLoss < siga) {
121 G4double x;
122 do {
123 loss = twomeanLoss*G4UniformRand();
124 x = (loss - meanLoss)/siga;
125 } while (1.0 - 0.5*x*x < G4UniformRand());
126 } else {
127 do {
128 loss = G4RandGauss::shoot(meanLoss,siga);
129 } while (0.0 > loss || loss > twomeanLoss);
130 }
131
132 // Poisson fluctuations
133 } else {
134 G4double n = (G4double)(G4Poisson(navr));
135 loss = meanLoss*n/navr;
136 }
137 //G4cout << "loss= " << loss << G4endl;
138
139 return loss;
140}
141
142//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
143
144G4double G4BohrFluctuations::Dispersion(const G4Material* material,
145 const G4DynamicParticle* dp,
146 G4double& tmax,
147 G4double& length)
148{
149 if(!particle) { InitialiseMe(dp->GetDefinition()); }
150
151 G4double electronDensity = material->GetElectronDensity();
152 kineticEnergy = dp->GetKineticEnergy();
153 G4double etot = kineticEnergy + particleMass;
154 beta2 = kineticEnergy*(kineticEnergy + 2.0*particleMass)/(etot*etot);
155 G4double siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
156 * electronDensity * chargeSquare;
157
158 return siga;
159}
160
161//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
162
163
Note: See TracBrowser for help on using the repository browser.