| Geant4 Cross Reference |
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 // ClassName: G4TablesForExtrapolator
29 //
30 // Description: This class provide calculatio
31 // and msc angle
32 //
33 // Author: 09.12.04 V.Ivanchenko
34 //
35 // Modification:
36 // 08-04-05 Rename Propogator -> Extrapolator
37 // 16-03-06 Add muon tables and fix bug in uni
38 // 21-03-06 Add verbosity defined in the const
39 // start only when first public metho
40 // 03-05-06 Remove unused pointer G4Material*
41 // 12-05-06 SEt linLossLimit=0.001 (VI)
42 //
43 //--------------------------------------------
44 //
45
46 //....oooOO0OOooo........oooOO0OOooo........oo
47
48 #include "G4TablesForExtrapolator.hh"
49 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "G4LossTableManager.hh"
52 #include "G4PhysicsLogVector.hh"
53 #include "G4ParticleDefinition.hh"
54 #include "G4Material.hh"
55 #include "G4MaterialCutsCouple.hh"
56 #include "G4Electron.hh"
57 #include "G4Positron.hh"
58 #include "G4Proton.hh"
59 #include "G4MuonPlus.hh"
60 #include "G4MuonMinus.hh"
61 #include "G4EmParameters.hh"
62 #include "G4MollerBhabhaModel.hh"
63 #include "G4BetheBlochModel.hh"
64 #include "G4eBremsstrahlungRelModel.hh"
65 #include "G4MuPairProductionModel.hh"
66 #include "G4MuBremsstrahlungModel.hh"
67 #include "G4ProductionCuts.hh"
68 #include "G4LossTableBuilder.hh"
69 #include "G4WentzelVIModel.hh"
70 #include "G4ios.hh"
71
72 //....oooOO0OOooo........oooOO0OOooo........oo
73
74 G4TablesForExtrapolator::G4TablesForExtrapolat
75
76 : emin(e1), emax(e2), verbose(verb), nbins(b
77 {
78 electron = G4Electron::Electron();
79 positron = G4Positron::Positron();
80 proton = G4Proton::Proton();
81 muonPlus = G4MuonPlus::MuonPlus();
82 muonMinus= G4MuonMinus::MuonMinus();
83 }
84
85 //....oooOO0OOooo........oooOO0OOooo........oo
86
87 G4TablesForExtrapolator::~G4TablesForExtrapola
88 {
89 if(nullptr != dedxElectron) {
90 dedxElectron->clearAndDestroy();
91 delete dedxElectron;
92 }
93 if(nullptr != dedxPositron) {
94 dedxPositron->clearAndDestroy();
95 delete dedxPositron;
96 }
97 if(nullptr != dedxProton) {
98 dedxProton->clearAndDestroy();
99 delete dedxProton;
100 }
101 if(nullptr != dedxMuon) {
102 dedxMuon->clearAndDestroy();
103 delete dedxMuon;
104 }
105 if(nullptr != rangeElectron) {
106 rangeElectron->clearAndDestroy();
107 delete rangeElectron;
108 }
109 if(nullptr != rangePositron) {
110 rangePositron->clearAndDestroy();
111 delete rangePositron;
112 }
113 if(nullptr != rangeProton) {
114 rangeProton->clearAndDestroy();
115 delete rangeProton;
116 }
117 if(nullptr != rangeMuon) {
118 rangeMuon->clearAndDestroy();
119 delete rangeMuon;
120 }
121 if(nullptr != invRangeElectron) {
122 invRangeElectron->clearAndDestroy();
123 delete invRangeElectron;
124 }
125 if(nullptr != invRangePositron) {
126 invRangePositron->clearAndDestroy();
127 delete invRangePositron;
128 }
129 if(nullptr != invRangeProton) {
130 invRangeProton->clearAndDestroy();
131 delete invRangeProton;
132 }
133 if(nullptr != invRangeMuon) {
134 invRangeMuon->clearAndDestroy();
135 delete invRangeMuon;
136 }
137 if(nullptr != mscElectron) {
138 mscElectron->clearAndDestroy();
139 delete mscElectron;
140 }
141 delete pcuts;
142 delete builder;
143 }
144
145 //....oooOO0OOooo........oooOO0OOooo........oo
146
147 const G4PhysicsTable*
148 G4TablesForExtrapolator::GetPhysicsTable(ExtTa
149 {
150 const G4PhysicsTable* table = nullptr;
151 switch (type)
152 {
153 case fDedxElectron:
154 table = dedxElectron;
155 break;
156 case fDedxPositron:
157 table = dedxPositron;
158 break;
159 case fDedxProton:
160 table = dedxProton;
161 break;
162 case fDedxMuon:
163 table = dedxMuon;
164 break;
165 case fRangeElectron:
166 table = rangeElectron;
167 break;
168 case fRangePositron:
169 table = rangePositron;
170 break;
171 case fRangeProton:
172 table = rangeProton;
173 break;
174 case fRangeMuon:
175 table = rangeMuon;
176 break;
177 case fInvRangeElectron:
178 table = invRangeElectron;
179 break;
180 case fInvRangePositron:
181 table = invRangePositron;
182 break;
183 case fInvRangeProton:
184 table = invRangeProton;
185 break;
186 case fInvRangeMuon:
187 table = invRangeMuon;
188 break;
189 case fMscElectron:
190 table = mscElectron;
191 }
192 return table;
193 }
194
195 //....oooOO0OOooo........oooOO0OOooo........oo
196
197 void G4TablesForExtrapolator::Initialisation()
198 {
199 if(verbose>1) {
200 G4cout << "### G4TablesForExtrapolator::In
201 }
202 G4int num = (G4int)G4Material::GetNumberOfMa
203 if(nmat == num) { return; }
204 nmat = num;
205 cuts.resize(nmat, DBL_MAX);
206 couples.resize(nmat, nullptr);
207
208 const G4MaterialTable* mtable = G4Material::
209 if(!pcuts) { pcuts = new G4ProductionCuts();
210
211 for(G4int i=0; i<nmat; ++i) {
212 couples[i] = new G4MaterialCutsCouple((*mt
213 }
214
215 dedxElectron = PrepareTable(dedxElectron
216 dedxPositron = PrepareTable(dedxPositron
217 dedxMuon = PrepareTable(dedxMuon);
218 dedxProton = PrepareTable(dedxProton);
219 rangeElectron = PrepareTable(rangeElectro
220 rangePositron = PrepareTable(rangePositro
221 rangeMuon = PrepareTable(rangeMuon);
222 rangeProton = PrepareTable(rangeProton)
223 invRangeElectron = PrepareTable(invRangeElec
224 invRangePositron = PrepareTable(invRangePosi
225 invRangeMuon = PrepareTable(invRangeMuon
226 invRangeProton = PrepareTable(invRangeProt
227 mscElectron = PrepareTable(mscElectron)
228
229 builder = new G4LossTableBuilder(true);
230 builder->SetBaseMaterialActive(false);
231
232 if(verbose>1) {
233 G4cout << "### G4TablesForExtrapolator Bui
234 << G4endl;
235 }
236 ComputeElectronDEDX(electron, dedxElectron);
237 builder->BuildRangeTable(dedxElectron,rangeE
238 builder->BuildInverseRangeTable(rangeElectro
239
240 if(verbose>1) {
241 G4cout << "### G4TablesForExtrapolator Bui
242 << G4endl;
243 }
244 ComputeElectronDEDX(positron, dedxPositron);
245 builder->BuildRangeTable(dedxPositron, range
246 builder->BuildInverseRangeTable(rangePositro
247
248 if(verbose>1) {
249 G4cout << "### G4TablesForExtrapolator Bui
250 }
251 ComputeMuonDEDX(muonPlus, dedxMuon);
252 builder->BuildRangeTable(dedxMuon, rangeMuon
253 builder->BuildInverseRangeTable(rangeMuon, i
254
255 if(verbose>2) {
256 G4cout << "DEDX MUON" << G4endl;
257 G4cout << *dedxMuon << G4endl;
258 G4cout << "RANGE MUON" << G4endl;
259 G4cout << *rangeMuon << G4endl;
260 G4cout << "INVRANGE MUON" << G4endl;
261 G4cout << *invRangeMuon << G4endl;
262 }
263 if(verbose>1) {
264 G4cout << "### G4TablesForExtrapolator Bui
265 << G4endl;
266 }
267 ComputeProtonDEDX(proton, dedxProton);
268 builder->BuildRangeTable(dedxProton, rangePr
269 builder->BuildInverseRangeTable(rangeProton,
270
271 ComputeTrasportXS(electron, mscElectron);
272 }
273
274 //....oooOO0OOooo........oooOO0OOooo........oo
275
276 G4PhysicsTable* G4TablesForExtrapolator::Prepa
277 {
278 G4PhysicsTable* table = ptr;
279 if(nullptr == ptr) { table = new G4PhysicsTa
280 G4int n = (G4int)table->length();
281 for(G4int i=n; i<nmat; ++i) {
282 G4PhysicsVector* v = new G4PhysicsLogVecto
283 table->push_back(v);
284 }
285 return table;
286 }
287
288 //....oooOO0OOooo........oooOO0OOooo........oo
289
290 void G4TablesForExtrapolator::ComputeElectronD
291 const G4Part
292 G4PhysicsTable* table)
293 {
294 G4MollerBhabhaModel* ioni = new G4MollerBhab
295 G4eBremsstrahlungRelModel* brem = new G4eBre
296 ioni->Initialise(part, cuts);
297 brem->Initialise(part, cuts);
298 ioni->SetUseBaseMaterials(false);
299 brem->SetUseBaseMaterials(false);
300
301 mass = electron_mass_c2;
302 charge2 = 1.0;
303 currentParticle = part;
304
305 const G4MaterialTable* mtable = G4Material::
306 if(0<verbose) {
307 G4cout << "G4TablesForExtrapolator::Comput
308 << part->GetParticleName()
309 << G4endl;
310 }
311 for(G4int i=0; i<nmat; ++i) {
312
313 const G4Material* mat = (*mtable)[i];
314 if(1<verbose) {
315 G4cout << "i= " << i << " mat= " << mat
316 }
317 G4PhysicsVector* aVector = (*table)[i];
318
319 for(G4int j=0; j<=nbins; ++j) {
320
321 G4double e = aVector->Energy(j);
322 G4double dedx = ioni->ComputeDEDXPerVol
323 + brem->ComputeDEDXPerVolume(mat,part,e,e);
324 if(1<verbose) {
325 G4cout << "j= " << j
326 << " e(MeV)= " << e/MeV
327 << " dedx(Mev/cm)= " << dedx*c
328 << " dedx(Mev.cm2/g)= "
329 << dedx/((MeV*mat->GetDensity())/(g/cm2))
330 << G4endl;
331 }
332 aVector->PutValue(j,dedx);
333 }
334 if(splineFlag) { aVector->FillSecondDeriva
335 }
336 }
337
338 //....oooOO0OOooo........oooOO0OOooo........oo
339
340 void
341 G4TablesForExtrapolator::ComputeMuonDEDX(const
342 G4PhysicsTable* table)
343 {
344 G4BetheBlochModel* ioni = new G4BetheBlochMo
345 G4MuPairProductionModel* pair = new G4MuPair
346 G4MuBremsstrahlungModel* brem = new G4MuBrem
347 ioni->Initialise(part, cuts);
348 pair->Initialise(part, cuts);
349 brem->Initialise(part, cuts);
350 ioni->SetUseBaseMaterials(false);
351 pair->SetUseBaseMaterials(false);
352 brem->SetUseBaseMaterials(false);
353
354 mass = part->GetPDGMass();
355 charge2 = 1.0;
356 currentParticle = part;
357
358 const G4MaterialTable* mtable = G4Material::
359
360 if(0<verbose) {
361 G4cout << "G4TablesForExtrapolator::Comput
362 << part->GetParticleName()
363 << G4endl;
364 }
365
366 for(G4int i=0; i<nmat; ++i) {
367
368 const G4Material* mat = (*mtable)[i];
369 if(1<verbose) {
370 G4cout << "i= " << i << " mat= " << mat
371 }
372 G4PhysicsVector* aVector = (*table)[i];
373 for(G4int j=0; j<=nbins; ++j) {
374
375 G4double e = aVector->Energy(j);
376 G4double dedx = ioni->ComputeDEDXPerVol
377 pair->ComputeDEDXPerVolume(ma
378 brem->ComputeDEDXPerVolume(ma
379 aVector->PutValue(j,dedx);
380 if(1<verbose) {
381 G4cout << "j= " << j
382 << " e(MeV)= " << e/MeV
383 << " dedx(Mev/cm)= " << dedx*c
384 << " dedx(Mev/(g/cm2)= "
385 << dedx/((MeV*mat->GetDensity(
386 << G4endl;
387 }
388 }
389 if(splineFlag) { aVector->FillSecondDeriva
390 }
391 }
392
393 //....oooOO0OOooo........oooOO0OOooo........oo
394
395 void
396 G4TablesForExtrapolator::ComputeProtonDEDX(con
397 G4PhysicsTable*
398 {
399 G4BetheBlochModel* ioni = new G4BetheBlochMo
400 ioni->Initialise(part, cuts);
401 ioni->SetUseBaseMaterials(false);
402
403 mass = part->GetPDGMass();
404 charge2 = 1.0;
405 currentParticle = part;
406
407 const G4MaterialTable* mtable = G4Material::
408
409 if(0<verbose) {
410 G4cout << "G4TablesForExtrapolator::Comput
411 << part->GetParticleName()
412 << G4endl;
413 }
414
415 for(G4int i=0; i<nmat; ++i) {
416
417 const G4Material* mat = (*mtable)[i];
418 if(1<verbose)
419 G4cout << "i= " << i << " mat= " << mat
420 G4PhysicsVector* aVector = (*table)[i];
421 for(G4int j=0; j<=nbins; ++j) {
422
423 G4double e = aVector->Energy(j);
424 G4double dedx = ioni->ComputeDEDXPerVol
425 aVector->PutValue(j,dedx);
426 if(1<verbose) {
427 G4cout << "j= " << j
428 << " e(MeV)= " << e/MeV
429 << " dedx(Mev/cm)= " << dedx*c
430 << " dedx(Mev.cm2/g)= " << ded
431 << G4endl;
432 }
433 }
434 if(splineFlag) { aVector->FillSecondDeriva
435 }
436 }
437
438 //....oooOO0OOooo........oooOO0OOooo........oo
439
440 void
441 G4TablesForExtrapolator::ComputeTrasportXS(con
442 G4PhysicsTable* table)
443 {
444 G4WentzelVIModel* msc = new G4WentzelVIModel
445 msc->SetPolarAngleLimit(CLHEP::pi);
446 msc->Initialise(part, cuts);
447 msc->SetUseBaseMaterials(false);
448
449 mass = part->GetPDGMass();
450 charge2 = 1.0;
451 currentParticle = part;
452
453 const G4MaterialTable* mtable = G4Material::
454
455 if(0<verbose) {
456 G4cout << "G4TablesForExtrapolator::Comput
457 << part->GetParticleName()
458 << G4endl;
459 }
460
461 for(G4int i=0; i<nmat; ++i) {
462
463 const G4Material* mat = (*mtable)[i];
464 msc->SetCurrentCouple(couples[i]);
465 if(1<verbose)
466 G4cout << "i= " << i << " mat= " << mat
467 G4PhysicsVector* aVector = (*table)[i];
468 for(G4int j=0; j<=nbins; ++j) {
469
470 G4double e = aVector->Energy(j);
471 G4double xs = msc->CrossSectionPerVolum
472 aVector->PutValue(j,xs);
473 if(1<verbose) {
474 G4cout << "j= " << j << " e(MeV)= "
475 << " xs(1/mm)= " << xs*mm << G
476 }
477 }
478 if(splineFlag) { aVector->FillSecondDeriva
479 }
480 }
481
482 //....oooOO0OOooo........oooOO0OOooo........oo
483
484