Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4UniversalFluctuation.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/standard/src/G4UniversalFluctuation.cc (Version 11.3.0) and /processes/electromagnetic/standard/src/G4UniversalFluctuation.cc (Version 5.1)


  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 24 // *******************************************    
 25 //                                                
 26 //                                                
 27 // -------------------------------------------    
 28 //                                                
 29 // GEANT4 Class file                              
 30 //                                                
 31 //                                                
 32 // File name:     G4UniversalFluctuation          
 33 //                                                
 34 // Author:        V. Ivanchenko for Laszlo Urb    
 35 //                                                
 36 // Creation date: 03.01.2002                      
 37 //                                                
 38 // Modifications:                                 
 39 //                                                
 40 //                                                
 41                                                   
 42 //....oooOO0OOooo........oooOO0OOooo........oo    
 43 //....oooOO0OOooo........oooOO0OOooo........oo    
 44                                                   
 45 #include "G4UniversalFluctuation.hh"              
 46 #include "G4PhysicalConstants.hh"                 
 47 #include "G4SystemOfUnits.hh"                     
 48 #include "Randomize.hh"                           
 49 #include "G4Poisson.hh"                           
 50 #include "G4Material.hh"                          
 51 #include "G4MaterialCutsCouple.hh"                
 52 #include "G4DynamicParticle.hh"                   
 53 #include "G4ParticleDefinition.hh"                
 54 #include "G4Log.hh"                               
 55                                                   
 56 //....oooOO0OOooo........oooOO0OOooo........oo    
 57                                                   
 58 G4UniversalFluctuation::G4UniversalFluctuation    
 59  :G4VEmFluctuationModel(nam),                     
 60   minLoss(10.*CLHEP::eV)                          
 61 {                                                 
 62   rndmarray = new G4double[sizearray];            
 63 }                                                 
 64                                                   
 65 //....oooOO0OOooo........oooOO0OOooo........oo    
 66                                                   
 67 G4UniversalFluctuation::~G4UniversalFluctuatio    
 68 {                                                 
 69   delete [] rndmarray;                            
 70 }                                                 
 71                                                   
 72 //....oooOO0OOooo........oooOO0OOooo........oo    
 73                                                   
 74 void G4UniversalFluctuation::InitialiseMe(cons    
 75 {                                                 
 76   particle = part;                                
 77   particleMass = part->GetPDGMass();              
 78   const G4double q = part->GetPDGCharge()/CLHE    
 79                                                   
 80   // Derived quantities                           
 81   m_Inv_particleMass = 1.0 / particleMass;        
 82   m_massrate = CLHEP::electron_mass_c2 * m_Inv    
 83   chargeSquare = q*q;                             
 84 }                                                 
 85                                                   
 86 //....oooOO0OOooo........oooOO0OOooo........oo    
 87                                                   
 88 G4double                                          
 89 G4UniversalFluctuation::SampleFluctuations(con    
 90                                            con    
 91                                            con    
 92                                            con    
 93                                            con    
 94                                            con    
 95 {                                                 
 96   // Calculate actual loss from the mean loss.    
 97   // The model used to get the fluctuations is    
 98   // as in Glandz in Geant3 (Cern program libr    
 99   // L. Urban et al. NIM A362, p.416 (1995) an    
100                                                   
101   // shortcut for very small loss or from a st    
102   // (out of validity of the model)               
103   //                                              
104   if (averageLoss < minLoss) { return averageL    
105   meanLoss = averageLoss;                         
106   const G4double tkin  = dp->GetKineticEnergy(    
107   //G4cout<< "Emean= "<< meanLoss<< " tmax= "<    
108                                                   
109   if(dp->GetDefinition() != particle) { Initia    
110                                                   
111   CLHEP::HepRandomEngine* rndmEngineF = G4Rand    
112                                                   
113   const G4double gam   = tkin * m_Inv_particle    
114   const G4double gam2  = gam*gam;                 
115   const G4double beta  = dp->GetBeta();           
116   const G4double beta2 = beta*beta;               
117                                                   
118   G4double loss(0.), siga(0.);                    
119                                                   
120   const G4Material* material = couple->GetMate    
121                                                   
122   // Gaussian regime                              
123   // for heavy particles only and conditions      
124   // for Gauusian fluct. has been changed         
125   //                                              
126   if (particleMass > CLHEP::electron_mass_c2 &    
127       meanLoss >= minNumberInteractionsBohr*tc    
128                                                   
129     siga = std::sqrt((tmax/beta2 - 0.5*tcut)*C    
130                       length*chargeSquare*mate    
131     const G4double sn = meanLoss/siga;            
132                                                   
133     // thick target case                          
134     if (sn >= 2.0) {                              
135                                                   
136       const G4double twomeanLoss = meanLoss +     
137       do {                                        
138   loss = G4RandGauss::shoot(rndmEngineF, meanL    
139   // Loop checking, 03-Aug-2015, Vladimir Ivan    
140       } while  (0.0 > loss || twomeanLoss < lo    
141                                                   
142       // Gamma distribution                       
143     } else {                                      
144                                                   
145       const G4double neff = sn*sn;                
146       loss = meanLoss*G4RandGamma::shoot(rndmE    
147     }                                             
148     //G4cout << "Gauss: " << loss << G4endl;      
149     return loss;                                  
150   }                                               
151                                                   
152   auto ioni = material->GetIonisation();          
153   e0 = ioni->GetEnergy0fluct();                   
154                                                   
155   // very small step or low-density material      
156   if(tcut <= e0) { return meanLoss; }             
157                                                   
158   ipotFluct = ioni->GetMeanExcitationEnergy();    
159   ipotLogFluct = ioni->GetLogMeanExcEnergy();     
160                                                   
161   // width correction for small cuts              
162   const G4double scaling = std::min(1.+0.5*CLH    
163   meanLoss /= scaling;                            
164                                                   
165   w2 = (tcut > ipotFluct) ?                       
166     G4Log(2.*CLHEP::electron_mass_c2*beta2*gam    
167   return SampleGlandz(rndmEngineF, material, t    
168 }                                                 
169                                                   
170 //....oooOO0OOooo........oooOO0OOooo........oo    
171                                                   
172 G4double                                          
173 G4UniversalFluctuation::SampleGlandz(CLHEP::He    
174                                      const G4M    
175                                      const G4d    
176 {                                                 
177   G4double a1(0.0), a3(0.0);                      
178   G4double loss = 0.0;                            
179   G4double e1 = ipotFluct;                        
180                                                   
181   if(tcut > e1) {                                 
182     a1 = meanLoss*(1.-rate)/e1;                   
183     if(a1 < a0) {                                 
184       const G4double fwnow = 0.1+(fw-0.1)*std:    
185       a1 /= fwnow;                                
186       e1 *= fwnow;                                
187     } else {                                      
188       a1 /= fw;                                   
189       e1 *= fw;                                   
190     }                                             
191   }                                               
192                                                   
193   const G4double w1 = tcut/e0;                    
194   a3 = rate*meanLoss*(tcut - e0)/(e0*tcut*G4Lo    
195   if(a1 <= 0.) { a3 /= rate; }                    
196                                                   
197   //'nearly' Gaussian fluctuation if a1>nmaxCo    
198   G4double emean = 0.;                            
199   G4double sig2e = 0.;                            
200                                                   
201   // excitation of type 1                         
202   if(a1 > 0.0) { AddExcitation(rndmEngineF, a1    
203                                                   
204   if(sig2e > 0.0) { SampleGauss(rndmEngineF, e    
205                                                   
206   // ionisation                                   
207   if(a3 > 0.) {                                   
208     emean = 0.;                                   
209     sig2e = 0.;                                   
210     G4double p3 = a3;                             
211     G4double alfa = 1.;                           
212     if(a3 > nmaxCont) {                           
213       alfa = w1*(nmaxCont+a3)/(w1*nmaxCont+a3)    
214       const G4double alfa1  = alfa*G4Log(alfa)    
215       const G4double namean = a3*w1*(alfa-1.)/    
216       emean += namean*e0*alfa1;                   
217       sig2e += e0*e0*namean*(alfa-alfa1*alfa1)    
218       p3 = a3 - namean;                           
219     }                                             
220                                                   
221     const G4double w3 = alfa*e0;                  
222     if(tcut > w3) {                               
223       const G4double w = (tcut-w3)/tcut;          
224       const G4int nnb = (G4int)G4Poisson(p3);     
225       if(nnb > 0) {                               
226         if(nnb > sizearray) {                     
227           sizearray = nnb;                        
228           delete [] rndmarray;                    
229           rndmarray = new G4double[nnb];          
230         }                                         
231         rndmEngineF->flatArray(nnb, rndmarray)    
232         for (G4int k=0; k<nnb; ++k) { loss +=     
233       }                                           
234     }                                             
235     if(sig2e > 0.0) { SampleGauss(rndmEngineF,    
236   }                                               
237   //G4cout << "### loss=" << loss << G4endl;      
238   return loss;                                    
239 }                                                 
240                                                   
241 //....oooOO0OOooo........oooOO0OOooo........oo    
242                                                   
243                                                   
244 G4double G4UniversalFluctuation::Dispersion(      
245                           const G4Material* ma    
246                           const G4DynamicParti    
247                           const G4double tcut,    
248                           const G4double tmax,    
249                           const G4double lengt    
250 {                                                 
251   if(dp->GetDefinition() != particle) { Initia    
252   const G4double beta = dp->GetBeta();            
253   return (tmax/(beta*beta) - 0.5*tcut) * CLHEP    
254     * material->GetElectronDensity() * chargeS    
255 }                                                 
256                                                   
257 //....oooOO0OOooo........oooOO0OOooo........oo    
258                                                   
259 void                                              
260 G4UniversalFluctuation::SetParticleAndCharge(c    
261                                              G    
262 {                                                 
263   if(part != particle) {                          
264     particle = part;                              
265     particleMass = part->GetPDGMass();            
266                                                   
267     // Derived quantities                         
268     m_Inv_particleMass = 1.0 / particleMass;      
269     m_massrate = CLHEP::electron_mass_c2 * m_I    
270   }                                               
271   chargeSquare = q2;                              
272 }                                                 
273                                                   
274 //....oooOO0OOooo........oooOO0OOooo........oo    
275