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