Geant4 Cross Reference

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


  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 // File name:     G4PenelopeBremsstrahlungAngu    
 30 //                                                
 31 // Author:        Luciano Pandola                 
 32 //                                                
 33 // Creation date: November 2010                   
 34 //                                                
 35 // History:                                       
 36 // -----------                                    
 37 // 23 Nov 2010  L. Pandola       1st implement    
 38 // 24 May 2011  L. Pandola       Renamed (make    
 39 // 13 Mar 2012  L. Pandola       Made a derive    
 40 //                               and update th    
 41 // 18 Jul 2012  L. Pandola       Migrated to t    
 42 //                               Now returns a    
 43 // 03 Oct 2013  L. Pandola       Migrated to M    
 44 // 17 Oct 2013  L. Pandola       Partially rev    
 45 //                                thread-local    
 46 //                                                
 47 //--------------------------------------------    
 48                                                   
 49 #include "G4PenelopeBremsstrahlungAngular.hh"     
 50                                                   
 51 #include "globals.hh"                             
 52 #include "G4PhysicalConstants.hh"                 
 53 #include "G4SystemOfUnits.hh"                     
 54 #include "G4PhysicsFreeVector.hh"                 
 55 #include "G4PhysicsTable.hh"                      
 56 #include "G4Material.hh"                          
 57 #include "Randomize.hh"                           
 58 #include "G4Exp.hh"                               
 59                                                   
 60 G4PenelopeBremsstrahlungAngular::G4PenelopeBre    
 61   G4VEmAngularDistribution("Penelope"), fEffec    
 62   fLorentzTables1(nullptr),fLorentzTables2(nul    
 63                                                   
 64 {                                                 
 65   fDataRead = false;                              
 66   fVerbosityLevel = 0;                            
 67 }                                                 
 68                                                   
 69 //....oooOO0OOooo........oooOO0OOooo........oo    
 70                                                   
 71 G4PenelopeBremsstrahlungAngular::~G4PenelopeBr    
 72 {                                                 
 73   ClearTables();                                  
 74 }                                                 
 75                                                   
 76 //....oooOO0OOooo........oooOO0OOooo........oo    
 77                                                   
 78 void G4PenelopeBremsstrahlungAngular::Initiali    
 79 {                                                 
 80   ClearTables();                                  
 81 }                                                 
 82                                                   
 83 //....oooOO0OOooo........oooOO0OOooo........oo    
 84                                                   
 85 void G4PenelopeBremsstrahlungAngular::ClearTab    
 86 {                                                 
 87   if (fLorentzTables1)                            
 88     {                                             
 89       for (auto j=fLorentzTables1->begin(); j     
 90         {                                         
 91     G4PhysicsTable* tab = j->second;              
 92           tab->clearAndDestroy();                 
 93           delete tab;                             
 94         }                                         
 95       fLorentzTables1->clear();                   
 96       delete fLorentzTables1;                     
 97       fLorentzTables1 = nullptr;                  
 98     }                                             
 99                                                   
100   if (fLorentzTables2)                            
101     {                                             
102       for (auto j=fLorentzTables2->begin(); j     
103         {                                         
104     G4PhysicsTable* tab = j->second;              
105     tab->clearAndDestroy();                       
106           delete tab;                             
107         }                                         
108       fLorentzTables2->clear();                   
109       delete fLorentzTables2;                     
110       fLorentzTables2 = nullptr;                  
111     }                                             
112   if (fEffectiveZSq)                              
113     {                                             
114       delete fEffectiveZSq;                       
115       fEffectiveZSq = nullptr;                    
116     }                                             
117 }                                                 
118                                                   
119 //....oooOO0OOooo........oooOO0OOooo........oo    
120                                                   
121 void G4PenelopeBremsstrahlungAngular::ReadData    
122 {                                                 
123    //Read information from DataBase file          
124   const char* path = G4FindDataDir("G4LEDATA")    
125   if (!path)                                      
126     {                                             
127       G4String excep =                            
128   "G4PenelopeBremsstrahlungAngular - G4LEDATA     
129       G4Exception("G4PenelopeBremsstrahlungAng    
130       "em0006",FatalException,excep);             
131       return;                                     
132     }                                             
133   G4String pathString(path);                      
134   G4String pathFile = pathString + "/penelope/    
135   std::ifstream file(pathFile);                   
136                                                   
137   if (!file.is_open())                            
138     {                                             
139       G4String excep = "G4PenelopeBremsstrahlu    
140       G4Exception("G4PenelopeBremsstrahlungAng    
141       "em0003",FatalException,excep);             
142       return;                                     
143     }                                             
144   G4int i=0,j=0,k=0; // i=index for Z, j=index    
145                                                   
146   for (k=0;k<fNumberofKPoints;k++)                
147     for (i=0;i<fNumberofZPoints;i++)              
148       for (j=0;j<fNumberofEPoints;j++)            
149   {                                               
150     G4double a1,a2;                               
151     G4int ik1,iz1,ie1;                            
152     G4double zr,er,kr;                            
153     file >> iz1 >> ie1 >> ik1 >> zr >> er >> k    
154     //check the data are correct                  
155     if ((iz1-1 == i) && (ik1-1 == k) && (ie1-1    
156       {                                           
157         fQQ1[i][j][k]=a1;                         
158         fQQ2[i][j][k]=a2;                         
159       }                                           
160     else                                          
161       {                                           
162         G4ExceptionDescription ed;                
163         ed << "Corrupted data file " << pathFi    
164         G4Exception("G4PenelopeBremsstrahlungA    
165       "em0005",FatalException,ed);                
166       }                                           
167   }                                               
168   file.close();                                   
169   fDataRead = true;                               
170 }                                                 
171                                                   
172 //....oooOO0OOooo........oooOO0OOooo........oo    
173                                                   
174 void G4PenelopeBremsstrahlungAngular::PrepareT    
175 {                                                 
176   //Unused at the moment: the G4PenelopeBremss    
177   //builds its own version of the tables.         
178                                                   
179   //Check if data file has already been read      
180   if (!fDataRead)                                 
181     {                                             
182       ReadDataFile();                             
183       if (!fDataRead)                             
184   G4Exception("G4PenelopeBremsstrahlungAngular    
185         "em2001",FatalException,"Unable to bui    
186     }                                             
187                                                   
188   if (!fLorentzTables1)                           
189       fLorentzTables1 = new std::map<G4double,    
190   if (!fLorentzTables2)                           
191     fLorentzTables2 = new std::map<G4double,G4    
192                                                   
193   G4double Zmat = CalculateEffectiveZ(material    
194                                                   
195   const G4int reducedEnergyGrid=21;               
196   //Support arrays.                               
197   G4double betas[fNumberofEPoints]; //betas fo    
198   //tables for interpolation                      
199   G4double Q1[fNumberofEPoints][fNumberofKPoin    
200   G4double Q2[fNumberofEPoints][fNumberofKPoin    
201   //expanded tables for interpolation             
202   G4double Q1E[fNumberofEPoints][reducedEnergy    
203   G4double Q2E[fNumberofEPoints][reducedEnergy    
204   G4double pZ[fNumberofZPoints] = {2.0,8.0,13.    
205                                                   
206   G4int i=0,j=0,k=0; // i=index for Z, j=index    
207   //Interpolation in Z                            
208   for (i=0;i<fNumberofEPoints;i++)                
209     {                                             
210       for (j=0;j<fNumberofKPoints;j++)            
211   {                                               
212     G4PhysicsFreeVector* QQ1vector =              
213       new G4PhysicsFreeVector(fNumberofZPoints    
214     G4PhysicsFreeVector* QQ2vector =              
215       new G4PhysicsFreeVector(fNumberofZPoints    
216                                                   
217     //fill vectors                                
218     for (k=0;k<fNumberofZPoints;k++)              
219       {                                           
220         QQ1vector->PutValues(k,pZ[k],G4Log(fQQ    
221         QQ2vector->PutValues(k,pZ[k],fQQ2[k][i    
222       }                                           
223     //Filled: now calculate derivatives           
224     QQ1vector->FillSecondDerivatives();           
225     QQ2vector->FillSecondDerivatives();           
226                                                   
227     Q1[i][j]= G4Exp(QQ1vector->Value(Zmat));      
228     Q2[i][j]=QQ2vector->Value(Zmat);              
229     delete QQ1vector;                             
230     delete QQ2vector;                             
231   }                                               
232     }                                             
233   G4double pE[fNumberofEPoints] = {1.0e-03*MeV    
234           1.0e-01*MeV,5.0e-01*MeV};               
235   G4double pK[fNumberofKPoints] = {0.0,0.6,0.8    
236   G4double ppK[reducedEnergyGrid];                
237                                                   
238   for(i=0;i<reducedEnergyGrid;i++)                
239     ppK[i]=((G4double) i) * 0.05;                 
240                                                   
241                                                   
242   for(i=0;i<fNumberofEPoints;i++)                 
243     betas[i]=std::sqrt(pE[i]*(pE[i]+2*electron    
244                                                   
245                                                   
246   for (i=0;i<fNumberofEPoints;i++)                
247     {                                             
248       for (j=0;j<fNumberofKPoints;j++)            
249   Q1[i][j]=Q1[i][j]/Zmat;                         
250     }                                             
251                                                   
252   //Expanded table of distribution parameters     
253   for (i=0;i<fNumberofEPoints;i++)                
254     {                                             
255       G4PhysicsFreeVector* Q1vector = new G4Ph    
256       G4PhysicsFreeVector* Q2vector = new G4Ph    
257                                                   
258       for (j=0;j<fNumberofKPoints;j++)            
259   {                                               
260     Q1vector->PutValues(j,pK[j],G4Log(Q1[i][j]    
261     Q2vector->PutValues(j,pK[j],Q2[i][j]);        
262   }                                               
263                                                   
264       for (j=0;j<reducedEnergyGrid;j++)           
265   {                                               
266     Q1E[i][j]=Q1vector->Value(ppK[j]);            
267     Q2E[i][j]=Q2vector->Value(ppK[j]);            
268   }                                               
269       delete Q1vector;                            
270       delete Q2vector;                            
271     }                                             
272   //                                              
273   //TABLES to be stored                           
274   //                                              
275   G4PhysicsTable* theTable1 = new G4PhysicsTab    
276   G4PhysicsTable* theTable2 = new G4PhysicsTab    
277   // the table will contain reducedEnergyGrid     
278   // values of k,                                 
279   // Each of the G4PhysicsFreeVectors has a pr    
280   // y vs. E                                      
281   //                                              
282   //reserve space of the vectors.                 
283   for (j=0;j<reducedEnergyGrid;j++)               
284     {                                             
285       G4PhysicsFreeVector* thevec = new G4Phys    
286       theTable1->push_back(thevec);               
287       G4PhysicsFreeVector* thevec2 = new G4Phy    
288       theTable2->push_back(thevec2);              
289     }                                             
290                                                   
291   for (j=0;j<reducedEnergyGrid;j++)               
292     {                                             
293       G4PhysicsFreeVector* thevec = (G4Physics    
294       G4PhysicsFreeVector* thevec2 = (G4Physic    
295       for (i=0;i<fNumberofEPoints;i++)            
296   {                                               
297     thevec->PutValues(i,betas[i],Q1E[i][j]);      
298     thevec2->PutValues(i,betas[i],Q2E[i][j]);     
299   }                                               
300       //Vectors are filled: calculate derivati    
301       thevec->FillSecondDerivatives();            
302       thevec2->FillSecondDerivatives();           
303     }                                             
304                                                   
305   if (fLorentzTables1 && fLorentzTables2)         
306     {                                             
307       fLorentzTables1->insert(std::make_pair(Z    
308       fLorentzTables2->insert(std::make_pair(Z    
309     }                                             
310   else                                            
311     {                                             
312       G4ExceptionDescription ed;                  
313       ed << "Unable to create tables of Lorent    
314       ed << "<Z>= "  << Zmat << " in G4Penelop    
315       delete theTable1;                           
316       delete theTable2;                           
317       G4Exception("G4PenelopeBremsstrahlungAng    
318       "em2005",FatalException,ed);                
319     }                                             
320                                                   
321   return;                                         
322 }                                                 
323                                                   
324 //....oooOO0OOooo........oooOO0OOooo........oo    
325                                                   
326 G4ThreeVector& G4PenelopeBremsstrahlungAngular    
327                 G4double eGamma,                  
328                 G4int,                            
329                 const G4Material* material)       
330 {                                                 
331   if (!material)                                  
332     {                                             
333       G4Exception("G4PenelopeBremsstrahlungAng    
334       "em2040",FatalException,"The pointer to     
335       return fLocalDirection;                     
336     }                                             
337                                                   
338   //Retrieve the effective Z                      
339   G4double Zmat = 0;                              
340                                                   
341   //The model might be initialized incorrectly    
342   //G4PenelopeBremsstrahungModel: make sure it    
343   if (!fEffectiveZSq)                             
344     {                                             
345       G4Exception("G4PenelopeBremsstrahlungAng    
346       "em2040",JustWarning,"EffectiveZSq table    
347       PrepareTables(material,false);              
348       //return fLocalDirection;                   
349     }                                             
350                                                   
351   //found in the table: return it                 
352   if (fEffectiveZSq->count(material))             
353     Zmat = fEffectiveZSq->find(material)->seco    
354   else //this can happen in unit tests or when    
355     //models other than G4PenelopeBremsstrahun    
356     {                                             
357       G4Exception("G4PenelopeBremsstrahlungAng    
358       "em2040",JustWarning,"Material not found    
359       PrepareTables(material,false);              
360       Zmat = fEffectiveZSq->find(material)->se    
361       //      return fLocalDirection;             
362     }                                             
363                                                   
364   if (fVerbosityLevel > 0)                        
365     {                                             
366       G4cout << "Effective <Z> for material :     
367   " = " << Zmat << G4endl;                        
368     }                                             
369                                                   
370   G4double ePrimary = dp->GetKineticEnergy();     
371                                                   
372   G4double beta = std::sqrt(ePrimary*(ePrimary    
373     (ePrimary+electron_mass_c2);                  
374   G4double cdt = 0;                               
375   G4double sinTheta = 0;                          
376   G4double phi = 0;                               
377                                                   
378   //Use a pure dipole distribution for energy     
379   if (ePrimary > 500*keV)                         
380     {                                             
381       cdt = 2.0*G4UniformRand() - 1.0;            
382       if (G4UniformRand() > 0.75)                 
383   {                                               
384     if (cdt<0)                                    
385       cdt = -1.0*std::pow(-cdt,1./3.);            
386     else                                          
387       cdt = std::pow(cdt,1./3.);                  
388   }                                               
389       cdt = (cdt+beta)/(1.0+beta*cdt);            
390       //Get primary kinematics                    
391       sinTheta = std::sqrt(1. - cdt*cdt);         
392       phi  = twopi * G4UniformRand();             
393       fLocalDirection.set(sinTheta* std::cos(p    
394           sinTheta* std::sin(phi),                
395           cdt);                                   
396       //rotate                                    
397       fLocalDirection.rotateUz(dp->GetMomentum    
398       //return                                    
399       return fLocalDirection;                     
400     }                                             
401                                                   
402   if (!(fLorentzTables1->count(Zmat)) || !(fLo    
403     {                                             
404       G4ExceptionDescription ed;                  
405       ed << "Unable to retrieve Lorentz tables    
406       G4Exception("G4PenelopeBremsstrahlungAng    
407       "em2006",FatalException,ed);                
408     }                                             
409                                                   
410   //retrieve actual tables                        
411   const G4PhysicsTable* theTable1 = fLorentzTa    
412   const G4PhysicsTable* theTable2 = fLorentzTa    
413                                                   
414   G4double RK=20.0*eGamma/ePrimary;               
415   G4int ik=std::min((G4int) RK,19);               
416                                                   
417   G4double P10=0,P11=0,P1=0;                      
418   G4double P20=0,P21=0,P2=0;                      
419                                                   
420   //First coefficient                             
421   const G4PhysicsFreeVector* v1 = (G4PhysicsFr    
422   const G4PhysicsFreeVector* v2 = (G4PhysicsFr    
423   P10 = v1->Value(beta);                          
424   P11 = v2->Value(beta);                          
425   P1=P10+(RK-(G4double) ik)*(P11-P10);            
426                                                   
427   //Second coefficient                            
428   const G4PhysicsFreeVector* v3 = (G4PhysicsFr    
429   const G4PhysicsFreeVector* v4 = (G4PhysicsFr    
430   P20=v3->Value(beta);                            
431   P21=v4->Value(beta);                            
432   P2=P20+(RK-(G4double) ik)*(P21-P20);            
433                                                   
434   //Sampling from the Lorenz-trasformed dipole    
435   P1=std::min(G4Exp(P1)/beta,1.0);                
436   G4double betap = std::min(std::max(beta*(1.0    
437                                                   
438   G4double testf=0;                               
439                                                   
440   if (G4UniformRand() < P1)                       
441     {                                             
442       do{                                         
443   cdt = 2.0*G4UniformRand()-1.0;                  
444   testf=2.0*G4UniformRand()-(1.0+cdt*cdt);        
445       }while(testf>0);                            
446     }                                             
447   else                                            
448     {                                             
449       do{                                         
450   cdt = 2.0*G4UniformRand()-1.0;                  
451   testf=G4UniformRand()-(1.0-cdt*cdt);            
452       }while(testf>0);                            
453     }                                             
454   cdt = (cdt+betap)/(1.0+betap*cdt);              
455                                                   
456   //Get primary kinematics                        
457   sinTheta = std::sqrt(1. - cdt*cdt);             
458   phi  = twopi * G4UniformRand();                 
459   fLocalDirection.set(sinTheta* std::cos(phi),    
460           sinTheta* std::sin(phi),                
461           cdt);                                   
462   //rotate                                        
463   fLocalDirection.rotateUz(dp->GetMomentumDire    
464   //return                                        
465   return fLocalDirection;                         
466                                                   
467 }                                                 
468                                                   
469 //....oooOO0OOooo........oooOO0OOooo........oo    
470                                                   
471 G4double G4PenelopeBremsstrahlungAngular::Calc    
472 {                                                 
473   if (!fEffectiveZSq)                             
474     fEffectiveZSq = new std::map<const G4Mater    
475                                                   
476   //Already exists: return it                     
477   if (fEffectiveZSq->count(material))             
478     return fEffectiveZSq->find(material)->seco    
479                                                   
480   //Helper for the calculation                    
481   std::vector<G4double> *StechiometricFactors     
482   G4int nElements = (G4int)material->GetNumber    
483   const G4ElementVector* elementVector = mater    
484   const G4double* fractionVector = material->G    
485   for (G4int i=0;i<nElements;++i)                 
486     {                                             
487       G4double fraction = fractionVector[i];      
488       G4double atomicWeigth = (*elementVector)    
489       StechiometricFactors->push_back(fraction    
490     }                                             
491   //Find max                                      
492   G4double MaxStechiometricFactor = 0.;           
493   for (G4int i=0;i<nElements;++i)                 
494     {                                             
495       if ((*StechiometricFactors)[i] > MaxStec    
496         MaxStechiometricFactor = (*Stechiometr    
497     }                                             
498   //Normalize                                     
499   for (G4int i=0;i<nElements;++i)                 
500     (*StechiometricFactors)[i] /=  MaxStechiom    
501                                                   
502   G4double sumz2 = 0;                             
503   G4double sums = 0;                              
504   for (G4int i=0;i<nElements;++i)                 
505     {                                             
506       G4double Z = (*elementVector)[i]->GetZ()    
507       sumz2 += (*StechiometricFactors)[i]*Z*Z;    
508       sums  += (*StechiometricFactors)[i];        
509     }                                             
510   delete StechiometricFactors;                    
511                                                   
512   G4double ZBR = std::sqrt(sumz2/sums);           
513   fEffectiveZSq->insert(std::make_pair(materia    
514                                                   
515   return ZBR;                                     
516 }                                                 
517