| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id: G4VEmModel.cc,v 1.20 2008/11/13 23:13:18 schaelic Exp $
>> 27 // GEANT4 tag $Name: geant4-09-02-patch-04 $
>> 28 //
26 // ------------------------------------------- 29 // -------------------------------------------------------------------
27 // 30 //
28 // GEANT4 Class file 31 // GEANT4 Class file
29 // 32 //
30 // 33 //
31 // File name: G4VEmModel 34 // File name: G4VEmModel
32 // 35 //
33 // Author: Vladimir Ivanchenko 36 // Author: Vladimir Ivanchenko
34 // 37 //
35 // Creation date: 25.07.2005 38 // Creation date: 25.07.2005
36 // 39 //
37 // Modifications: 40 // Modifications:
38 // 25.10.2005 Set default highLimit=100.TeV (V 41 // 25.10.2005 Set default highLimit=100.TeV (V.Ivanchenko)
39 // 06.02.2006 add method ComputeMeanFreePath() 42 // 06.02.2006 add method ComputeMeanFreePath() (mma)
40 // 16.02.2009 Move implementations of virtual <<
41 // 43 //
42 // 44 //
43 // Class Description: 45 // Class Description:
44 // 46 //
45 // Abstract interface to energy loss models 47 // Abstract interface to energy loss models
46 48
47 // ------------------------------------------- 49 // -------------------------------------------------------------------
48 // 50 //
49 51
50 #include "G4VEmModel.hh" 52 #include "G4VEmModel.hh"
51 #include "G4ElementData.hh" <<
52 #include "G4LossTableManager.hh" 53 #include "G4LossTableManager.hh"
53 #include "G4LossTableBuilder.hh" <<
54 #include "G4ProductionCutsTable.hh" 54 #include "G4ProductionCutsTable.hh"
55 #include "G4ParticleChangeForLoss.hh" <<
56 #include "G4ParticleChangeForGamma.hh" <<
57 #include "G4EmParameters.hh" <<
58 #include "G4SystemOfUnits.hh" <<
59 #include "G4EmUtility.hh" <<
60 #include "G4Log.hh" <<
61 #include "Randomize.hh" <<
62 #include <iostream> <<
63 55
64 //....oooOO0OOooo........oooOO0OOooo........oo 56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
65 //....oooOO0OOooo........oooOO0OOooo........oo 57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
66 58
67 G4VEmModel::G4VEmModel(const G4String& nam): 59 G4VEmModel::G4VEmModel(const G4String& nam):
68 inveplus(1.0/CLHEP::eplus), << 60 fluc(0), name(nam), lowLimit(0.1*keV), highLimit(100.0*TeV),
69 lowLimit(0.1*CLHEP::keV), << 61 polarAngleLimit(0.0),secondaryThreshold(DBL_MAX),theLPMflag(false),
70 highLimit(100.0*CLHEP::TeV), << 62 pParticleChange(0),nuclearStopping(false),nsec(5)
71 polarAngleLimit(CLHEP::pi), <<
72 name(nam) <<
73 { 63 {
74 xsec.resize(nsec); 64 xsec.resize(nsec);
75 fEmManager = G4LossTableManager::Instance(); << 65 nSelectors = 0;
76 fEmManager->Register(this); << 66 G4LossTableManager::Instance()->Register(this);
77 isMaster = fEmManager->IsMaster(); <<
78 <<
79 G4LossTableBuilder* bld = fEmManager->GetTab <<
80 theDensityFactor = bld->GetDensityFactors(); <<
81 theDensityIdx = bld->GetCoupleIndexes(); <<
82 } 67 }
83 68
84 //....oooOO0OOooo........oooOO0OOooo........oo 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
85 70
86 G4VEmModel::~G4VEmModel() 71 G4VEmModel::~G4VEmModel()
87 { 72 {
88 if(localElmSelectors) { << 73 G4LossTableManager::Instance()->DeRegister(this);
89 for(G4int i=0; i<nSelectors; ++i) { << 74 G4int n = elmSelectors.size();
90 delete (*elmSelectors)[i]; << 75 if(n > 0) {
91 } << 76 for(G4int i=0; i<n; i++) {
92 delete elmSelectors; << 77 delete elmSelectors[i];
93 } <<
94 delete anglModel; <<
95 <<
96 if(localTable && xSectionTable != nullptr) { <<
97 xSectionTable->clearAndDestroy(); <<
98 delete xSectionTable; <<
99 xSectionTable = nullptr; <<
100 } <<
101 fEmManager->DeRegister(this); <<
102 } <<
103 <<
104 //....oooOO0OOooo........oooOO0OOooo........oo <<
105 <<
106 G4ParticleChangeForLoss* G4VEmModel::GetPartic <<
107 { <<
108 G4ParticleChangeForLoss* p = nullptr; <<
109 if (pParticleChange != nullptr) { <<
110 p = static_cast<G4ParticleChangeForLoss*>( <<
111 } else { <<
112 p = new G4ParticleChangeForLoss(); <<
113 pParticleChange = p; <<
114 } <<
115 if(fTripletModel != nullptr) { fTripletModel <<
116 return p; <<
117 } <<
118 <<
119 //....oooOO0OOooo........oooOO0OOooo........oo <<
120 <<
121 G4ParticleChangeForGamma* G4VEmModel::GetParti <<
122 { <<
123 G4ParticleChangeForGamma* p = nullptr; <<
124 if (pParticleChange != nullptr) { <<
125 p = static_cast<G4ParticleChangeForGamma*> <<
126 } else { <<
127 p = new G4ParticleChangeForGamma(); <<
128 pParticleChange = p; <<
129 } <<
130 if(fTripletModel != nullptr) { fTripletModel <<
131 return p; <<
132 } <<
133 <<
134 //....oooOO0OOooo........oooOO0OOooo........oo <<
135 <<
136 void G4VEmModel::InitialiseElementSelectors(co <<
137 co <<
138 { <<
139 if(highLimit <= lowLimit) { return; } <<
140 G4EmUtility::InitialiseElementSelectors(this <<
141 localElmSelectors = true; <<
142 } <<
143 <<
144 //....oooOO0OOooo........oooOO0OOooo........oo <<
145 <<
146 void G4VEmModel::InitialiseLocal(const G4Parti <<
147 {} <<
148 <<
149 //....oooOO0OOooo........oooOO0OOooo........oo <<
150 <<
151 void G4VEmModel::InitialiseForMaterial(const G <<
152 const G <<
153 { <<
154 if(material != nullptr) { <<
155 G4int n = (G4int)material->GetNumberOfElem <<
156 for(G4int i=0; i<n; ++i) { <<
157 G4int Z = material->GetElement(i)->GetZa <<
158 InitialiseForElement(part, Z); <<
159 } 78 }
160 } 79 }
161 } 80 }
162 81
163 //....oooOO0OOooo........oooOO0OOooo........oo 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
164 83
165 void G4VEmModel::InitialiseForElement(const G4 << 84 G4double G4VEmModel::CrossSectionPerVolume(const G4Material* material,
166 {} << 85 const G4ParticleDefinition* p,
167 << 86 G4double ekin,
168 //....oooOO0OOooo........oooOO0OOooo........oo << 87 G4double emin,
169 << 88 G4double emax)
170 G4double G4VEmModel::ComputeDEDXPerVolume(cons <<
171 cons <<
172 G4do <<
173 { 89 {
174 return 0.0; << 90 SetupForMaterial(p, material, ekin);
175 } << 91 G4double cross = 0.0;
176 << 92 const G4ElementVector* theElementVector = material->GetElementVector();
177 //....oooOO0OOooo........oooOO0OOooo........oo << 93 const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
178 << 94 G4int nelm = material->GetNumberOfElements();
179 G4double G4VEmModel::CrossSectionPerVolume(con <<
180 con <<
181 G4d <<
182 G4d <<
183 G4d <<
184 { <<
185 SetupForMaterial(p, mat, ekin); <<
186 const G4double* theAtomNumDensityVector = ma <<
187 G4int nelm = (G4int)mat->GetNumberOfElements <<
188 if(nelm > nsec) { 95 if(nelm > nsec) {
189 xsec.resize(nelm); 96 xsec.resize(nelm);
190 nsec = nelm; 97 nsec = nelm;
191 } 98 }
192 G4double cross = 0.0; << 99 for (G4int i=0; i<nelm; i++) {
193 for (G4int i=0; i<nelm; ++i) { <<
194 cross += theAtomNumDensityVector[i]* 100 cross += theAtomNumDensityVector[i]*
195 ComputeCrossSectionPerAtom(p,mat->GetEle << 101 ComputeCrossSectionPerAtom(p,(*theElementVector)[i],ekin,emin,emax);
196 xsec[i] = cross; 102 xsec[i] = cross;
197 } 103 }
198 return cross; 104 return cross;
199 } 105 }
200 106
201 //....oooOO0OOooo........oooOO0OOooo........oo 107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
202 108
203 G4double G4VEmModel::GetPartialCrossSection(co << 109 G4double G4VEmModel::ComputeMeanFreePath(const G4ParticleDefinition* p,
204 co << 110 G4double ekin,
205 G4 << 111 const G4Material* material,
206 { << 112 G4double emin,
207 return 0.0; << 113 G4double emax)
208 } <<
209 <<
210 //....oooOO0OOooo........oooOO0OOooo........oo <<
211 <<
212 void G4VEmModel::StartTracking(G4Track*) <<
213 {} <<
214 <<
215 //....oooOO0OOooo........oooOO0OOooo........oo <<
216 <<
217 const G4Element* G4VEmModel::SelectRandomAtom( <<
218 <<
219 <<
220 <<
221 <<
222 { <<
223 G4int n = (G4int)mat->GetNumberOfElements(); <<
224 fCurrentElement = mat->GetElement(0); <<
225 if (n > 1) { <<
226 const G4double x = G4UniformRand()* <<
227 G4VEmModel::CrossSectionPerVolume(mat,pd <<
228 for(G4int i=0; i<n; ++i) { <<
229 if (x <= xsec[i]) { <<
230 fCurrentElement = mat->GetElement(i); <<
231 break; <<
232 } <<
233 } <<
234 } <<
235 return fCurrentElement; <<
236 } <<
237 <<
238 //....oooOO0OOooo........oooOO0OOooo........oo <<
239 <<
240 const G4Element* G4VEmModel::GetCurrentElement <<
241 { <<
242 const G4Element* elm = fCurrentElement; <<
243 if(nullptr == elm && nullptr != mat) { <<
244 elm = G4EmUtility::SampleRandomElement(mat <<
245 } <<
246 return elm; <<
247 } <<
248 <<
249 //....oooOO0OOooo........oooOO0OOooo........oo <<
250 <<
251 G4int G4VEmModel::SelectRandomAtomNumber(const <<
252 { <<
253 const G4Element* elm = GetCurrentElement(mat <<
254 return (nullptr == elm) ? 0 : elm->GetZasInt <<
255 } <<
256 <<
257 //....oooOO0OOooo........oooOO0OOooo........oo <<
258 <<
259 const G4Isotope* G4VEmModel::GetCurrentIsotope <<
260 { <<
261 const G4Isotope* iso = nullptr; <<
262 const G4Element* el = elm; <<
263 if(nullptr == el && nullptr != fCurrentCoupl <<
264 el = GetCurrentElement(fCurrentCouple->Get <<
265 } <<
266 if(nullptr != el) { <<
267 iso = G4EmUtility::SampleRandomIsotope(el) <<
268 } <<
269 return iso; <<
270 } <<
271 <<
272 //....oooOO0OOooo........oooOO0OOooo........oo <<
273 <<
274 G4int G4VEmModel::SelectIsotopeNumber(const G4 <<
275 { <<
276 auto iso = GetCurrentIsotope(elm); <<
277 return (nullptr != iso) ? iso->GetN() : 0; <<
278 } <<
279 <<
280 //....oooOO0OOooo........oooOO0OOooo........oo <<
281 <<
282 G4double G4VEmModel::ComputeCrossSectionPerAto <<
283 <<
284 <<
285 { <<
286 return 0.0; <<
287 } <<
288 <<
289 //....oooOO0OOooo........oooOO0OOooo........oo <<
290 <<
291 G4double <<
292 G4VEmModel::ComputeCrossSectionPerShell(const <<
293 G4int, <<
294 G4doub <<
295 { <<
296 return 0.0; <<
297 } <<
298 <<
299 //....oooOO0OOooo........oooOO0OOooo........oo <<
300 <<
301 void G4VEmModel::DefineForRegion(const G4Regio <<
302 {} <<
303 <<
304 //....oooOO0OOooo........oooOO0OOooo........oo <<
305 <<
306 void G4VEmModel::FillNumberOfSecondaries(G4int <<
307 G4int <<
308 { <<
309 numberOfTriplets = 0; <<
310 numberOfRecoil = 0; <<
311 } <<
312 <<
313 //....oooOO0OOooo........oooOO0OOooo........oo <<
314 <<
315 G4double G4VEmModel::ChargeSquareRatio(const G <<
316 { <<
317 return GetChargeSquareRatio(track.GetParticl <<
318 track.GetMateria <<
319 } <<
320 <<
321 //....oooOO0OOooo........oooOO0OOooo........oo <<
322 <<
323 G4double G4VEmModel::GetChargeSquareRatio(cons <<
324 cons <<
325 { <<
326 const G4double q = p->GetPDGCharge()*inveplu <<
327 return q*q; <<
328 } <<
329 <<
330 //....oooOO0OOooo........oooOO0OOooo........oo <<
331 <<
332 G4double G4VEmModel::GetParticleCharge(const G <<
333 const G <<
334 { 114 {
335 return p->GetPDGCharge(); << 115 G4double mfp = DBL_MAX;
>> 116 G4double cross = CrossSectionPerVolume(material,p,ekin,emin,emax);
>> 117 if (cross > DBL_MIN) mfp = 1./cross;
>> 118 return mfp;
336 } 119 }
337 120
338 //....oooOO0OOooo........oooOO0OOooo........oo 121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
339 122
340 void G4VEmModel::CorrectionsAlongStep(const G4 << 123 void G4VEmModel::InitialiseElementSelectors(const G4ParticleDefinition* p,
341 const G4 << 124 const G4DataVector& cuts)
342 const G4 <<
343 {} <<
344 <<
345 //....oooOO0OOooo........oooOO0OOooo........oo <<
346 <<
347 G4double G4VEmModel::Value(const G4MaterialCut <<
348 const G4ParticleDef <<
349 { 125 {
350 SetCurrentCouple(couple); << 126 G4int nbins = G4int(std::log10(highLimit/lowLimit) + 0.5);
351 return pFactor*e*e*CrossSectionPerVolume(pBa << 127 if(nbins < 3) nbins = 3;
352 } << 128 G4bool spline = G4LossTableManager::Instance()->SplineFlag();
353 129
354 //....oooOO0OOooo........oooOO0OOooo........oo << 130 G4ProductionCutsTable* theCoupleTable=
>> 131 G4ProductionCutsTable::GetProductionCutsTable();
>> 132 G4int numOfCouples = theCoupleTable->GetTableSize();
355 133
356 G4double G4VEmModel::MinPrimaryEnergy(const G4 << 134 // prepare vector
357 const G4 << 135 if(numOfCouples > nSelectors) {
358 G4double << 136 elmSelectors.resize(numOfCouples);
359 { << 137 nSelectors = numOfCouples;
360 return 0.0; <<
361 } <<
362 <<
363 //....oooOO0OOooo........oooOO0OOooo........oo <<
364 <<
365 G4double G4VEmModel::MinEnergyCut(const G4Part <<
366 const G4Mate <<
367 { <<
368 return 0.0; <<
369 } <<
370 <<
371 //....oooOO0OOooo........oooOO0OOooo........oo <<
372 <<
373 G4double G4VEmModel::MaxSecondaryEnergy(const <<
374 G4doub <<
375 { <<
376 return kineticEnergy; <<
377 } <<
378 <<
379 //....oooOO0OOooo........oooOO0OOooo........oo <<
380 <<
381 void G4VEmModel::SetupForMaterial(const G4Part <<
382 const G4Mate <<
383 { <<
384 GetChargeSquareRatio(p, mat, ekin); <<
385 } <<
386 <<
387 //....oooOO0OOooo........oooOO0OOooo........oo <<
388 <<
389 void <<
390 G4VEmModel::SetParticleChange(G4VParticleChang <<
391 { <<
392 if(p != nullptr && pParticleChange != p) { p <<
393 if(flucModel != f) { flucModel = f; } <<
394 } <<
395 <<
396 //....oooOO0OOooo........oooOO0OOooo........oo <<
397 <<
398 void G4VEmModel::SetCrossSectionTable(G4Physic <<
399 { <<
400 xSectionTable = p; <<
401 localTable = isLocal; <<
402 } <<
403 <<
404 //....oooOO0OOooo........oooOO0OOooo........oo <<
405 <<
406 void G4VEmModel::SetLPMFlag(G4bool) <<
407 { <<
408 if (G4EmParameters::Instance()->Verbose() > <<
409 G4ExceptionDescription ed; <<
410 ed << "The obsolete method SetLPMFlag(..) <<
411 << " is called. Please, use G4EmParamet <<
412 << " instead"; <<
413 G4Exception("G4VEmModel::SetLPMFlag", "em0 <<
414 } 138 }
415 } <<
416 139
417 //....oooOO0OOooo........oooOO0OOooo........oo << 140 // initialise vector
>> 141 for(G4int i=0; i<numOfCouples; i++) {
>> 142 const G4MaterialCutsCouple* couple =
>> 143 theCoupleTable->GetMaterialCutsCouple(i);
>> 144 const G4Material* material = couple->GetMaterial();
>> 145 G4int idx = couple->GetIndex();
418 146
419 void G4VEmModel::SetMasterThread(G4bool) << 147 // selector already exist check if should be deleted
420 {} << 148 G4bool create = true;
>> 149 if(elmSelectors[i]) {
>> 150 if(material == elmSelectors[i]->GetMaterial()) create = false;
>> 151 else delete elmSelectors[i];
>> 152 }
>> 153 if(create) {
>> 154 elmSelectors[i] = new G4EmElementSelector(this,material,nbins,
>> 155 lowLimit,highLimit,spline);
>> 156 }
>> 157 elmSelectors[i]->Initialise(p, cuts[idx]);
>> 158 //elmSelectors[i]->Dump(p);
>> 159 }
>> 160 }
421 161
422 //....oooOO0OOooo........oooOO0OOooo........oo 162 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
423 163
424 void G4VEmModel::ModelDescription(std::ostrea <<
425 { <<
426 outFile << "The description for this model h <<
427 } <<
428 164
429 //....oooOO0OOooo........oooOO0OOooo........oo <<
430 165