Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4LivermoreBremsstrahlungModel.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/lowenergy/src/G4LivermoreBremsstrahlungModel.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4LivermoreBremsstrahlungModel.cc (Version 9.2.p4)


  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:     G4LivermoreBremsstrahlungMod    
 33 //                                                
 34 // Author:        Vladimir Ivanchenko use inhe    
 35 //                base class implementing ultr    
 36 //                model                           
 37 //                                                
 38 // Creation date: 04.10.2011                      
 39 //                                                
 40 // Modifications:                                 
 41 //                                                
 42 // -------------------------------------------    
 43 //                                                
 44 //....oooOO0OOooo........oooOO0OOooo........oo    
 45 //....oooOO0OOooo........oooOO0OOooo........oo    
 46                                                   
 47 #include "G4LivermoreBremsstrahlungModel.hh"      
 48 #include "G4PhysicalConstants.hh"                 
 49 #include "G4SystemOfUnits.hh"                     
 50 #include "G4Electron.hh"                          
 51 #include "G4Positron.hh"                          
 52 #include "G4Gamma.hh"                             
 53 #include "Randomize.hh"                           
 54 #include "G4AutoLock.hh"                          
 55 #include "G4Material.hh"                          
 56 #include "G4Element.hh"                           
 57 #include "G4ElementVector.hh"                     
 58 #include "G4ProductionCutsTable.hh"               
 59 #include "G4ParticleChangeForLoss.hh"             
 60 #include "G4Generator2BS.hh"                      
 61                                                   
 62 #include "G4Physics2DVector.hh"                   
 63 #include "G4Exp.hh"                               
 64 #include "G4Log.hh"                               
 65                                                   
 66 #include "G4ios.hh"                               
 67 #include <fstream>                                
 68 #include <iomanip>                                
 69                                                   
 70 //....oooOO0OOooo........oooOO0OOooo........oo    
 71                                                   
 72 namespace { G4Mutex LivermoreBremsstrahlungMod    
 73 using namespace std;                              
 74                                                   
 75                                                   
 76 G4Physics2DVector* G4LivermoreBremsstrahlungMo    
 77 G4double G4LivermoreBremsstrahlungModel::ylimi    
 78 G4double G4LivermoreBremsstrahlungModel::expnu    
 79                                                   
 80 static const G4double emaxlog = 4*G4Log(10.);     
 81 static const G4double alpha = CLHEP::twopi*CLH    
 82 static const G4double epeaklimit= 300*CLHEP::M    
 83 static const G4double elowlimit = 20*CLHEP::ke    
 84                                                   
 85 G4LivermoreBremsstrahlungModel::G4LivermoreBre    
 86   const G4ParticleDefinition* p, const G4Strin    
 87   : G4eBremsstrahlungRelModel(p,nam),useBicubi    
 88 {                                                 
 89   SetLowEnergyLimit(10.0*eV);                     
 90   SetAngularDistribution(new G4Generator2BS())    
 91 }                                                 
 92                                                   
 93 //....oooOO0OOooo........oooOO0OOooo........oo    
 94                                                   
 95 G4LivermoreBremsstrahlungModel::~G4LivermoreBr    
 96 {                                                 
 97   if(IsMaster()) {                                
 98     for(size_t i=0; i<101; ++i) {                 
 99       if(dataSB[i]) {                             
100   delete dataSB[i];                               
101   dataSB[i] = nullptr;                            
102       }                                           
103     }                                             
104   }                                               
105 }                                                 
106                                                   
107 //....oooOO0OOooo........oooOO0OOooo........oo    
108                                                   
109 void G4LivermoreBremsstrahlungModel::Initialis    
110             const G4DataVector& cuts)             
111 {                                                 
112   // Access to elements                           
113   if(IsMaster()) {                                
114     // check environment variable                 
115     // Build the complete string identifying t    
116     const char* path = G4FindDataDir("G4LEDATA    
117                                                   
118     const G4ElementTable* theElmTable = G4Elem    
119     size_t numOfElm = G4Element::GetNumberOfEl    
120     if(numOfElm > 0) {                            
121       for(size_t i=0; i<numOfElm; ++i) {          
122   G4int Z = (*theElmTable)[i]->GetZasInt();       
123   if(Z < 1)        { Z = 1; }                     
124   else if(Z > 100) { Z = 100; }                   
125   //G4cout << "Z= " << Z << G4endl;               
126   // Initialisation                               
127   if(!dataSB[Z]) { ReadData(Z, path); }           
128       }                                           
129     }                                             
130   }                                               
131   G4eBremsstrahlungRelModel::Initialise(p, cut    
132 }                                                 
133                                                   
134 //....oooOO0OOooo........oooOO0OOooo........oo    
135                                                   
136 G4String G4LivermoreBremsstrahlungModel::Direc    
137 {                                                 
138   return "/livermore/brem/br";                    
139 }                                                 
140                                                   
141 //....oooOO0OOooo........oooOO0OOooo........oo    
142                                                   
143 void G4LivermoreBremsstrahlungModel::ReadData(    
144 {                                                 
145   if(dataSB[Z]) { return; }                       
146   const char* datadir = path;                     
147                                                   
148   if(nullptr == datadir) {                        
149     datadir = G4FindDataDir("G4LEDATA");          
150     if(!datadir) {                                
151       G4Exception("G4LivermoreBremsstrahlungMo    
152       FatalException,"Environment variable G4L    
153       return;                                     
154     }                                             
155   }                                               
156   std::ostringstream ost;                         
157   ost << datadir << DirectoryPath() << Z;         
158   std::ifstream fin(ost.str().c_str());           
159   if( !fin.is_open()) {                           
160     G4ExceptionDescription ed;                    
161     ed << "Bremsstrahlung data file <" << ost.    
162        << "> is not opened!";                     
163     G4Exception("G4LivermoreBremsstrahlungMode    
164     FatalException,ed,                            
165     "G4LEDATA version should be G4EMLOW8.0 or     
166     return;                                       
167   }                                               
168   //G4cout << "G4LivermoreBremsstrahlungModel     
169   //   << ">" << G4endl;                          
170   G4Physics2DVector* v = new G4Physics2DVector    
171   if(v->Retrieve(fin)) {                          
172     if(useBicubicInterpolation) { v->SetBicubi    
173     dataSB[Z] = v;                                
174     ylimit[Z] = v->Value(0.97, emaxlog, idx, i    
175   } else {                                        
176     G4ExceptionDescription ed;                    
177     ed << "Bremsstrahlung data file <" << ost.    
178        << "> is not retrieved!";                  
179     G4Exception("G4LivermoreBremsstrahlungMode    
180                 FatalException,ed,                
181     "G4LEDATA version should be G4EMLOW8.0 or     
182     delete v;                                     
183   }                                               
184 }                                                 
185                                                   
186 //....oooOO0OOooo........oooOO0OOooo........oo    
187                                                   
188 G4double                                          
189 G4LivermoreBremsstrahlungModel::ComputeDXSecti    
190 {                                                 
191   if(gammaEnergy < 0.0 || fPrimaryKinEnergy <=    
192   G4double x = gammaEnergy/fPrimaryKinEnergy;     
193   G4double y = G4Log(fPrimaryKinEnergy/MeV);      
194   G4int Z = fCurrentIZ;                           
195                                                   
196   //G4cout << "G4LivermoreBremsstrahlungModel:    
197   //   << " x= " << x << " y= " << y << " " <<    
198   if(!dataSB[Z]) { InitialiseForElement(0, Z);    
199                                                   
200   G4double invb2 = fPrimaryTotalEnergy*fPrimar    
201                    *(fPrimaryKinEnergy + 2.*fP    
202   G4double cross = dataSB[Z]->Value(x,y,idx,id    
203                                                   
204   if(!fIsElectron) {                              
205     G4double invbeta1 = sqrt(invb2);              
206     G4double e2 = fPrimaryKinEnergy - gammaEne    
207     if(e2 > 0.0) {                                
208       G4double invbeta2 = (e2 + fPrimaryPartic    
209                           /sqrt(e2*(e2 + 2.*fP    
210       G4double xxx = alpha*fCurrentIZ*(invbeta    
211       if(xxx < expnumlim) { cross = 0.0; }        
212       else { cross *= G4Exp(xxx); }               
213     } else {                                      
214       cross = 0.0;                                
215     }                                             
216   }                                               
217                                                   
218   return cross;                                   
219 }                                                 
220                                                   
221 //....oooOO0OOooo........oooOO0OOooo........oo    
222                                                   
223 void                                              
224 G4LivermoreBremsstrahlungModel::SampleSecondar    
225                                         std::v    
226           const G4MaterialCutsCouple* couple,     
227           const G4DynamicParticle* dp,            
228           G4double cutEnergy,                     
229           G4double maxEnergy)                     
230 {                                                 
231   G4double kineticEnergy = dp->GetKineticEnerg    
232   G4double cut  = std::min(cutEnergy, kineticE    
233   G4double emax = std::min(maxEnergy, kineticE    
234   if(cut >= emax) { return; }                     
235   // sets total energy, kinetic energy and den    
236   SetupForMaterial(fPrimaryParticle, couple->G    
237                                                   
238   const G4Element* elm =                          
239     SelectRandomAtom(couple,fPrimaryParticle,k    
240   fCurrentIZ = elm->GetZasInt();                  
241   G4int Z = fCurrentIZ;                           
242                                                   
243   G4double totMomentum = sqrt(kineticEnergy*(f    
244   /*                                              
245   G4cout << "G4LivermoreBremsstrahlungModel::S    
246    << kineticEnergy/MeV                           
247    << " Z= " << Z << " cut(MeV)= " << cut/MeV     
248    << " emax(MeV)= " << emax/MeV << " corr= "     
249   */                                              
250   G4double xmin = G4Log(cut*cut + fDensityCorr    
251   G4double xmax = G4Log(emax*emax  + fDensityC    
252   G4double y = G4Log(kineticEnergy/MeV);          
253                                                   
254   G4double gammaEnergy, v;                        
255                                                   
256   // majoranta                                    
257   G4double x0 = cut/kineticEnergy;                
258   G4double vmax = dataSB[Z]->Value(x0, y, idx,    
259                                                   
260   // majoranta corrected for e-                   
261   if(fIsElectron && x0 < 0.97 &&                  
262      ((kineticEnergy > epeaklimit) || (kinetic    
263     G4double ylim = std::min(ylimit[Z],1.1*dat    
264     if(ylim > vmax) { vmax = ylim; }              
265   }                                               
266   if(x0 < 0.05) { vmax *= 1.2; }                  
267                                                   
268   do {                                            
269     //++ncount;                                   
270     G4double x = G4Exp(xmin + G4UniformRand()*    
271     if(x < 0.0) { x = 0.0; }                      
272     gammaEnergy = sqrt(x);                        
273     G4double x1 = gammaEnergy/kineticEnergy;      
274     v = dataSB[Z]->Value(x1, y, idx, idy);        
275                                                   
276     // correction for positrons                   
277     if(!fIsElectron) {                            
278       G4double e1 = kineticEnergy - cut;          
279       G4double invbeta1 = (e1 + fPrimaryPartic    
280                           /sqrt(e1*(e1 + 2*fPr    
281       G4double e2 = kineticEnergy - gammaEnerg    
282       G4double invbeta2 = (e2 + fPrimaryPartic    
283                           /sqrt(e2*(e2 + 2*fPr    
284       G4double xxx = twopi*fine_structure_cons    
285                                                   
286       if(xxx < expnumlim) { v = 0.0; }            
287       else { v *= G4Exp(xxx); }                   
288     }                                             
289                                                   
290     if (v > 1.05*vmax && nwarn < 5) {             
291       ++nwarn;                                    
292       G4ExceptionDescription ed;                  
293       ed << "### G4LivermoreBremsstrahlungMode    
294    << v << " > " << vmax << " by " << v/vmax      
295    << " Egamma(MeV)= " << gammaEnergy             
296    << " Ee(MeV)= " << kineticEnergy               
297    << " Z= " << Z << "  " << fPrimaryParticle-    
298                                                   
299       if ( 20 == nwarn ) {                        
300   ed << "\n ### G4LivermoreBremsstrahlungModel    
301       }                                           
302       G4Exception("G4LivermoreBremsstrahlungMo    
303       JustWarning, ed,"");                        
304                                                   
305     }                                             
306   } while (v < vmax*G4UniformRand());             
307                                                   
308   //                                              
309   // angles of the emitted gamma. ( Z - axis a    
310   // use general interface                        
311   //                                              
312                                                   
313   G4ThreeVector gammaDirection =                  
314     GetAngularDistribution()->SampleDirection(    
315                 Z, couple->GetMaterial());        
316                                                   
317   // create G4DynamicParticle object for the G    
318   G4DynamicParticle* gamma =                      
319     new G4DynamicParticle(fGammaParticle,gamma    
320   vdp->push_back(gamma);                          
321                                                   
322   G4ThreeVector direction = (totMomentum*dp->G    
323            - gammaEnergy*gammaDirection).unit(    
324                                                   
325   /*                                              
326   G4cout << "### G4SBModel: v= "                  
327    << " Eg(MeV)= " << gammaEnergy                 
328    << " Ee(MeV)= " << kineticEnergy               
329    << " DirE " << direction << " DirG " << gam    
330    << G4endl;                                     
331   */                                              
332   // energy of primary                            
333   G4double finalE = kineticEnergy - gammaEnerg    
334                                                   
335   // stop tracking and create new secondary in    
336   if(gammaEnergy > SecondaryThreshold()) {        
337     fParticleChange->ProposeTrackStatus(fStopA    
338     fParticleChange->SetProposedKineticEnergy(    
339     G4DynamicParticle* el =                       
340       new G4DynamicParticle(const_cast<G4Parti    
341           direction, finalE);                     
342     vdp->push_back(el);                           
343                                                   
344     // continue tracking                          
345   } else {                                        
346     fParticleChange->SetProposedMomentumDirect    
347     fParticleChange->SetProposedKineticEnergy(    
348   }                                               
349 }                                                 
350                                                   
351 //....oooOO0OOooo........oooOO0OOooo........oo    
352                                                   
353 void G4LivermoreBremsstrahlungModel::Initialis    
354                                      const G4P    
355              G4int Z)                             
356 {                                                 
357   G4AutoLock l(&LivermoreBremsstrahlungModelMu    
358   if(!dataSB[Z]) { ReadData(Z); }                 
359   l.unlock();                                     
360 }                                                 
361                                                   
362 //....oooOO0OOooo........oooOO0OOooo........oo    
363