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

Last change on this file since 1199 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 5.6 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.8 2009/09/29 11:33:22 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
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
71//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
72
73G4BohrFluctuations::~G4BohrFluctuations()
74{}
75
76//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77
78void G4BohrFluctuations::InitialiseMe(const G4ParticleDefinition* part)
79{
80 particle = part;
81 particleMass = part->GetPDGMass();
82 G4double q = part->GetPDGCharge()/eplus;
83 chargeSquare = q*q;
84}
85
86//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
87
88G4double G4BohrFluctuations::SampleFluctuations(const G4Material* material,
89 const G4DynamicParticle* dp,
90 G4double& tmax,
91 G4double& length,
92 G4double& meanLoss)
93{
94 if(meanLoss <= minLoss) return meanLoss;
95 G4double siga = Dispersion(material,dp,tmax,length);
96 G4double loss = meanLoss;
97
98 G4double navr = meanLoss*meanLoss/siga;
99 //G4cout << "### meanLoss= " << meanLoss << " navr= " << navr << G4endl;
100 if (navr >= minNumberInteractionsBohr) {
101
102 // Increase fluctuations for big fractional energy loss
103 if ( meanLoss > minFraction*kineticEnergy ) {
104 G4double gam = (kineticEnergy - meanLoss)/particleMass + 1.0;
105 G4double b2 = 1.0 - 1.0/(gam*gam);
106 if(b2 < xmin*beta2) b2 = xmin*beta2;
107 G4double x = b2/beta2;
108 G4double x3 = 1.0/(x*x*x);
109 siga *= 0.25*(1.0 + x)*(x3 + (1.0/b2 - 0.5)/(1.0/beta2 - 0.5) );
110 }
111 siga = sqrt(siga);
112 G4double twomeanLoss = meanLoss + meanLoss;
113 //G4cout << "siga= " << siga << " 2edp= " << twomeanLoss <<G4endl;
114
115 if(twomeanLoss < siga) {
116 G4double x;
117 do {
118 loss = twomeanLoss*G4UniformRand();
119 x = (loss - meanLoss)/siga;
120 } while (1.0 - 0.5*x*x < G4UniformRand());
121 } else {
122 do {
123 loss = G4RandGauss::shoot(meanLoss,siga);
124 } while (0.0 > loss || loss > twomeanLoss);
125 }
126
127 // Poisson fluctuations
128 } else {
129 G4double n = (G4double)(G4Poisson(navr));
130 loss = meanLoss*n/navr;
131 }
132 //G4cout << "loss= " << loss << G4endl;
133
134 return loss;
135}
136
137//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
138
139G4double G4BohrFluctuations::Dispersion(const G4Material* material,
140 const G4DynamicParticle* dp,
141 G4double& tmax,
142 G4double& length)
143{
144 if(!particle) InitialiseMe(dp->GetDefinition());
145
146 G4double electronDensity = material->GetElectronDensity();
147 kineticEnergy = dp->GetKineticEnergy();
148 G4double etot = kineticEnergy + particleMass;
149 beta2 = kineticEnergy*(kineticEnergy + 2.0*particleMass)/(etot*etot);
150 G4double siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
151 * electronDensity * chargeSquare;
152
153 return siga;
154}
155
156//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
157
158
Note: See TracBrowser for help on using the repository browser.