Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/muons/src/G4EnergyLossForExtrapolator.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/muons/src/G4EnergyLossForExtrapolator.cc (Version 11.3.0) and /processes/electromagnetic/muons/src/G4EnergyLossForExtrapolator.cc (Version 4.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 // ClassName:    G4EnergyLossForExtrapolator      
 30 //                                                
 31 // Description:  This class provide calculatio    
 32 //               and msc angle                    
 33 //                                                
 34 // Author:       09.12.04 V.Ivanchenko            
 35 //                                                
 36 // Modification:                                  
 37 // 08-04-05 Rename Propogator -> Extrapolator     
 38 // 16-03-06 Add muon tables and fix bug in uni    
 39 // 21-03-06 Add verbosity defined in the const    
 40 //          start only when first public metho    
 41 // 03-05-06 Remove unused pointer G4Material*     
 42 // 12-05-06 SEt linLossLimit=0.001 (VI)           
 43 //                                                
 44 //--------------------------------------------    
 45 //                                                
 46                                                   
 47 //....oooOO0OOooo........oooOO0OOooo........oo    
 48                                                   
 49 #include "G4EnergyLossForExtrapolator.hh"         
 50 #include "G4PhysicalConstants.hh"                 
 51 #include "G4SystemOfUnits.hh"                     
 52 #include "G4ParticleDefinition.hh"                
 53 #include "G4Material.hh"                          
 54 #include "G4MaterialCutsCouple.hh"                
 55 #include "G4Electron.hh"                          
 56 #include "G4Positron.hh"                          
 57 #include "G4Proton.hh"                            
 58 #include "G4MuonPlus.hh"                          
 59 #include "G4MuonMinus.hh"                         
 60 #include "G4ParticleTable.hh"                     
 61                                                   
 62 //....oooOO0OOooo........oooOO0OOooo........oo    
 63                                                   
 64 #ifdef G4MULTITHREADED                            
 65 G4Mutex G4EnergyLossForExtrapolator::extrMutex    
 66 #endif                                            
 67                                                   
 68 G4TablesForExtrapolator* G4EnergyLossForExtrap    
 69                                                   
 70 G4EnergyLossForExtrapolator::G4EnergyLossForEx    
 71   : maxEnergyTransfer(DBL_MAX), verbose(verb)     
 72 {                                                 
 73   emin = 1.*CLHEP::MeV;                           
 74   emax = 100.*CLHEP::TeV;                         
 75 }                                                 
 76                                                   
 77 //....oooOO0OOooo........oooOO0OOooo........oo    
 78                                                   
 79 G4EnergyLossForExtrapolator::~G4EnergyLossForE    
 80 {                                                 
 81   if(isMaster) {                                  
 82     delete tables;                                
 83     tables = nullptr;                             
 84   }                                               
 85 }                                                 
 86                                                   
 87 //....oooOO0OOooo........oooOO0OOooo........oo    
 88                                                   
 89 G4double                                          
 90 G4EnergyLossForExtrapolator::EnergyAfterStep(G    
 91                G4double stepLength,               
 92                const G4Material* mat,             
 93                const G4ParticleDefinition* par    
 94 {                                                 
 95   G4double kinEnergyFinal = kinEnergy;            
 96   if(SetupKinematics(part, mat, kinEnergy)) {     
 97     G4double step = TrueStepLength(kinEnergy,s    
 98     G4double r  = ComputeRange(kinEnergy,part,    
 99     if(r <= step) {                               
100       kinEnergyFinal = 0.0;                       
101     } else if(step < linLossLimit*r) {            
102       kinEnergyFinal -= step*ComputeDEDX(kinEn    
103     } else {                                      
104       G4double r1 = r - step;                     
105       kinEnergyFinal = ComputeEnergy(r1,part,m    
106     }                                             
107   }                                               
108   return kinEnergyFinal;                          
109 }                                                 
110                                                   
111 //....oooOO0OOooo........oooOO0OOooo........oo    
112                                                   
113 G4double                                          
114 G4EnergyLossForExtrapolator::EnergyBeforeStep(    
115                 G4double stepLength,              
116                 const G4Material* mat,            
117                 const G4ParticleDefinition* pa    
118 {                                                 
119   //G4cout << "G4EnergyLossForExtrapolator::En    
120   G4double kinEnergyFinal = kinEnergy;            
121                                                   
122   if(SetupKinematics(part, mat, kinEnergy)) {     
123     G4double step = TrueStepLength(kinEnergy,s    
124     G4double r = ComputeRange(kinEnergy,part,m    
125                                                   
126     if(step < linLossLimit*r) {                   
127       kinEnergyFinal += step*ComputeDEDX(kinEn    
128     } else {                                      
129       G4double r1 = r + step;                     
130       kinEnergyFinal = ComputeEnergy(r1,part,m    
131     }                                             
132   }                                               
133   return kinEnergyFinal;                          
134 }                                                 
135                                                   
136 //....oooOO0OOooo........oooOO0OOooo........oo    
137                                                   
138 G4double                                          
139 G4EnergyLossForExtrapolator::TrueStepLength(G4    
140               G4double stepLength,                
141               const G4Material* mat,              
142               const G4ParticleDefinition* part    
143 {                                                 
144   G4double res = stepLength;                      
145   //G4cout << "## G4EnergyLossForExtrapolator:    
146   //   <<  "  " << part->GetParticleName() <<     
147   if(SetupKinematics(part, mat, kinEnergy)) {     
148     if(part == electron || part == positron) {    
149       const G4double x = stepLength*              
150   ComputeValue(kinEnergy, GetPhysicsTable(fMsc    
151       //G4cout << " x= " << x << G4endl;          
152       if(x < 0.2)         { res *= (1.0 + 0.5*    
153       else if(x < 0.9999) { res = -G4Log(1.0 -    
154       else { res = ComputeRange(kinEnergy, par    
155     } else {                                      
156       res = ComputeTrueStep(mat,part,kinEnergy    
157     }                                             
158   }                                               
159   return res;                                     
160 }                                                 
161                                                   
162 //....oooOO0OOooo........oooOO0OOooo........oo    
163                                                   
164 G4bool                                            
165 G4EnergyLossForExtrapolator::SetupKinematics(c    
166                const G4Material* mat,             
167                G4double kinEnergy)                
168 {                                                 
169   if(mat->GetNumberOfMaterials() != nmat) { In    
170   if(nullptr == part || nullptr == mat || kinE    
171     { return false; }                             
172   if(part != currentParticle) {                   
173     currentParticle = part;                       
174     G4double q = part->GetPDGCharge()/eplus;      
175     charge2 = q*q;                                
176   }                                               
177   if(mat != currentMaterial) {                    
178     size_t i = mat->GetIndex();                   
179     if(i >= nmat) {                               
180       G4cout << "### G4EnergyLossForExtrapolat    
181        << i << " above number of materials " <    
182       return false;                               
183     } else {                                      
184       currentMaterial = mat;                      
185       electronDensity = mat->GetElectronDensit    
186       radLength       = mat->GetRadlen();         
187     }                                             
188   }                                               
189   if(kinEnergy != kineticEnergy) {                
190     kineticEnergy = kinEnergy;                    
191     G4double mass = part->GetPDGMass();           
192     G4double tau  = kinEnergy/mass;               
193                                                   
194     gam   = tau + 1.0;                            
195     bg2   = tau * (tau + 2.0);                    
196     beta2 = bg2/(gam*gam);                        
197     tmax  = kinEnergy;                            
198     if(part == electron) tmax *= 0.5;             
199     else if(part != positron) {                   
200       G4double r = CLHEP::electron_mass_c2/mas    
201       tmax = 2.0*bg2*CLHEP::electron_mass_c2/(    
202     }                                             
203     tmax = std::min(tmax, maxEnergyTransfer);     
204   }                                               
205   return true;                                    
206 }                                                 
207                                                   
208 //....oooOO0OOooo........oooOO0OOooo........oo    
209                                                   
210 const G4ParticleDefinition*                       
211 G4EnergyLossForExtrapolator::FindParticle(cons    
212 {                                                 
213   currentParticle = G4ParticleTable::GetPartic    
214   if(nullptr == currentParticle) {                
215     G4cout << "### G4EnergyLossForExtrapolator    
216      << "FindParticle() fails to find " << nam    
217   }                                               
218   return currentParticle;                         
219 }                                                 
220                                                   
221 //....oooOO0OOooo........oooOO0OOooo........oo    
222                                                   
223 G4double                                          
224 G4EnergyLossForExtrapolator::ComputeDEDX(G4dou    
225            const G4ParticleDefinition* part,      
226                                          const    
227 {                                                 
228   if(mat->GetNumberOfMaterials() != nmat) { In    
229   G4double x = 0.0;                               
230   if(part == electron)  {                         
231     x = ComputeValue(ekin, GetPhysicsTable(fDe    
232   } else if(part == positron) {                   
233     x = ComputeValue(ekin, GetPhysicsTable(fDe    
234   } else if(part == muonPlus || part == muonMi    
235     x = ComputeValue(ekin, GetPhysicsTable(fDe    
236   } else {                                        
237     G4double e = ekin*CLHEP::proton_mass_c2/pa    
238     G4double q = part->GetPDGCharge()/CLHEP::e    
239     x = ComputeValue(e, GetPhysicsTable(fDedxP    
240   }                                               
241   return x;                                       
242 }                                                 
243                                                   
244 //....oooOO0OOooo........oooOO0OOooo........oo    
245                                                   
246 G4double                                          
247 G4EnergyLossForExtrapolator::ComputeRange(G4do    
248             const G4ParticleDefinition* part,     
249             const G4Material* mat)                
250 {                                                 
251   if(mat->GetNumberOfMaterials() != nmat) { In    
252   G4double x = 0.0;                               
253   if(part == electron) {                          
254     x = ComputeValue(ekin, GetPhysicsTable(fRa    
255   } else if(part == positron) {                   
256     x = ComputeValue(ekin, GetPhysicsTable(fRa    
257   } else if(part == muonPlus || part == muonMi    
258     x = ComputeValue(ekin, GetPhysicsTable(fRa    
259   } else {                                        
260     G4double massratio = CLHEP::proton_mass_c2    
261     G4double e = ekin*massratio;                  
262     G4double q = part->GetPDGCharge()/CLHEP::e    
263     x = ComputeValue(e, GetPhysicsTable(fRange    
264       /(q*q*massratio);                           
265   }                                               
266   return x;                                       
267 }                                                 
268                                                   
269 //....oooOO0OOooo........oooOO0OOooo........oo    
270                                                   
271 G4double                                          
272 G4EnergyLossForExtrapolator::ComputeEnergy(G4d    
273              const G4ParticleDefinition* part,    
274              const G4Material* mat)               
275 {                                                 
276   if(mat->GetNumberOfMaterials() != nmat) { In    
277   G4double x = 0.0;                               
278   if(part == electron) {                          
279     x = ComputeValue(range,GetPhysicsTable(fIn    
280   } else if(part == positron) {                   
281     x = ComputeValue(range,GetPhysicsTable(fIn    
282   } else if(part == muonPlus || part == muonMi    
283     x = ComputeValue(range, GetPhysicsTable(fI    
284   } else {                                        
285     G4double massratio = CLHEP::proton_mass_c2    
286     G4double q = part->GetPDGCharge()/CLHEP::e    
287     G4double r = range*massratio*q*q;             
288     x = ComputeValue(r, GetPhysicsTable(fInvRa    
289   }                                               
290   return x;                                       
291 }                                                 
292                                                   
293 //....oooOO0OOooo........oooOO0OOooo........oo    
294                                                   
295 G4double                                          
296 G4EnergyLossForExtrapolator::EnergyDispersion(    
297                 G4double stepLength,              
298                 const G4Material* mat,            
299                 const G4ParticleDefinition* pa    
300 {                                                 
301   G4double sig2 = 0.0;                            
302   if(SetupKinematics(part, mat, kinEnergy)) {     
303     G4double step = ComputeTrueStep(mat,part,k    
304     sig2 = (1.0/beta2 - 0.5)                      
305       *CLHEP::twopi_mc2_rcl2*tmax*step*electro    
306   }                                               
307   return sig2;                                    
308 }                                                 
309                                                   
310 //....oooOO0OOooo........oooOO0OOooo........oo    
311                                                   
312 G4double G4EnergyLossForExtrapolator::AverageS    
313                         G4double kinEnergy,       
314       G4double stepLength,                        
315       const G4Material* mat,                      
316       const G4ParticleDefinition* part)           
317 {                                                 
318   G4double theta = 0.0;                           
319   if(SetupKinematics(part, mat, kinEnergy)) {     
320     G4double t = stepLength/radLength;            
321     G4double y = std::max(0.001, t);              
322     theta = 19.23*CLHEP::MeV*std::sqrt(charge2    
323       /(beta2*gam*part->GetPDGMass());            
324   }                                               
325   return theta;                                   
326 }                                                 
327                                                   
328 //....oooOO0OOooo........oooOO0OOooo........oo    
329                                                   
330 void G4EnergyLossForExtrapolator::Initialisati    
331 {                                                 
332   if(verbose>0) {                                 
333     G4cout << "### G4EnergyLossForExtrapolator    
334      << tables << G4endl;                         
335   }                                               
336   electron = G4Electron::Electron();              
337   positron = G4Positron::Positron();              
338   proton   = G4Proton::Proton();                  
339   muonPlus = G4MuonPlus::MuonPlus();              
340   muonMinus= G4MuonMinus::MuonMinus();            
341                                                   
342   // initialisation for the 1st run               
343   if(nullptr == tables) {                         
344 #ifdef G4MULTITHREADED                            
345     G4MUTEXLOCK(&extrMutex);                      
346     if(nullptr == tables) {                       
347 #endif                                            
348       isMaster = true;                            
349       tables = new G4TablesForExtrapolator(ver    
350       tables->Initialisation();                   
351       nmat = G4Material::GetNumberOfMaterials(    
352       if(verbose > 0) {                           
353         G4cout << "### G4EnergyLossForExtrapol    
354                << nmat << " materials Nbins= "    
355                << nbins << " Emin(MeV)= " << e    
356                << G4endl;                         
357       }                                           
358 #ifdef G4MULTITHREADED                            
359     }                                             
360     G4MUTEXUNLOCK(&extrMutex);                    
361 #endif                                            
362   }                                               
363                                                   
364   // initialisation for the next run              
365   if(isMaster && G4Material::GetNumberOfMateri    
366 #ifdef G4MULTITHREADED                            
367     G4MUTEXLOCK(&extrMutex);                      
368 #endif                                            
369     tables->Initialisation();                     
370 #ifdef G4MULTITHREADED                            
371     G4MUTEXUNLOCK(&extrMutex);                    
372 #endif                                            
373   }                                               
374   nmat = G4Material::GetNumberOfMaterials();      
375 }                                                 
376                                                   
377 //....oooOO0OOooo........oooOO0OOooo........oo    
378