source: trunk/source/processes/electromagnetic/muons/src/G4EnergyLossForExtrapolator.cc@ 1036

Last change on this file since 1036 was 1007, checked in by garnier, 17 years ago

update to geant4.9.2
File size: 18.1 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// $Id: G4EnergyLossForExtrapolator.cc,v 1.18 2008/11/13 14:14:07 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-02 $
28//
29//---------------------------------------------------------------------------
30//
31// ClassName: G4EnergyLossForExtrapolator
32//
33// Description: This class provide calculation of energy loss, fluctuation,
34// and msc angle
35//
36// Author: 09.12.04 V.Ivanchenko
37//
38// Modification:
39// 08-04-05 Rename Propogator -> Extrapolator (V.Ivanchenko)
40// 16-03-06 Add muon tables and fix bug in units (V.Ivanchenko)
41// 21-03-06 Add verbosity defined in the constructor and Initialisation
42// start only when first public method is called (V.Ivanchenko)
43// 03-05-06 Remove unused pointer G4Material* from number of methods (VI)
44// 12-05-06 SEt linLossLimit=0.001 (VI)
45//
46//----------------------------------------------------------------------------
47//
48
49//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
50
51#include "G4EnergyLossForExtrapolator.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 "G4ParticleTable.hh"
62#include "G4LossTableBuilder.hh"
63#include "G4MollerBhabhaModel.hh"
64#include "G4BetheBlochModel.hh"
65#include "G4eBremsstrahlungModel.hh"
66#include "G4MuPairProductionModel.hh"
67#include "G4MuBremsstrahlungModel.hh"
68#include "G4ProductionCuts.hh"
69#include "G4LossTableManager.hh"
70#include "G4WentzelVIModel.hh"
71
72//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73
74G4EnergyLossForExtrapolator::G4EnergyLossForExtrapolator(G4int verb)
75 :maxEnergyTransfer(DBL_MAX),verbose(verb),isInitialised(false)
76{}
77
78//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
79
80G4EnergyLossForExtrapolator:: ~G4EnergyLossForExtrapolator()
81{
82 for(G4int i=0; i<nmat; i++) {delete couples[i];}
83 delete dedxElectron;
84 delete dedxPositron;
85 delete dedxProton;
86 delete rangeElectron;
87 delete rangePositron;
88 delete rangeProton;
89 delete invRangeElectron;
90 delete invRangePositron;
91 delete invRangeProton;
92 delete mscElectron;
93 delete cuts;
94}
95
96//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
97
98G4double G4EnergyLossForExtrapolator::EnergyAfterStep(G4double kinEnergy,
99 G4double stepLength,
100 const G4Material* mat,
101 const G4ParticleDefinition* part)
102{
103 if(!isInitialised) Initialisation();
104 G4double kinEnergyFinal = kinEnergy;
105 if(SetupKinematics(part, mat, kinEnergy)) {
106 G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
107 G4double r = ComputeRange(kinEnergy,part);
108 if(r <= step) {
109 kinEnergyFinal = 0.0;
110 } else if(step < linLossLimit*r) {
111 kinEnergyFinal -= step*ComputeDEDX(kinEnergy,part);
112 } else {
113 G4double r1 = r - step;
114 kinEnergyFinal = ComputeEnergy(r1,part);
115 }
116 }
117 return kinEnergyFinal;
118}
119
120//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
121
122G4double G4EnergyLossForExtrapolator::EnergyBeforeStep(G4double kinEnergy,
123 G4double stepLength,
124 const G4Material* mat,
125 const G4ParticleDefinition* part)
126{
127 if(!isInitialised) Initialisation();
128 G4double kinEnergyFinal = kinEnergy;
129
130 if(SetupKinematics(part, mat, kinEnergy)) {
131 G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
132 G4double r = ComputeRange(kinEnergy,part);
133
134 if(step < linLossLimit*r) {
135 kinEnergyFinal += step*ComputeDEDX(kinEnergy,part);
136 } else {
137 G4double r1 = r + step;
138 kinEnergyFinal = ComputeEnergy(r1,part);
139 }
140 }
141 return kinEnergyFinal;
142}
143
144//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
145
146G4double G4EnergyLossForExtrapolator::TrueStepLength(G4double kinEnergy,
147 G4double stepLength,
148 const G4Material* mat,
149 const G4ParticleDefinition* part)
150{
151 G4double res = stepLength;
152 if(!isInitialised) Initialisation();
153 if(SetupKinematics(part, mat, kinEnergy)) {
154 if(part == electron || part == positron) {
155 G4double x = stepLength*ComputeValue(kinEnergy, mscElectron);
156 if(x < 0.2) res *= (1.0 + 0.5*x + x*x/3.0);
157 else if(x < 0.9999) res = -std::log(1.0 - x)*stepLength/x;
158 else res = ComputeRange(kinEnergy,part);
159
160 } else {
161 res = ComputeTrueStep(mat,part,kinEnergy,stepLength);
162 }
163 }
164 return res;
165}
166
167//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
168
169G4bool G4EnergyLossForExtrapolator::SetupKinematics(const G4ParticleDefinition* part,
170 const G4Material* mat,
171 G4double kinEnergy)
172{
173 if(!part || !mat || kinEnergy < keV) return false;
174 if(!isInitialised) Initialisation();
175 G4bool flag = false;
176 if(part != currentParticle) {
177 flag = true;
178 currentParticle = part;
179 mass = part->GetPDGMass();
180 G4double q = part->GetPDGCharge()/eplus;
181 charge2 = q*q;
182 }
183 if(mat != currentMaterial) {
184 G4int i = mat->GetIndex();
185 if(i >= nmat) {
186 G4cout << "### G4EnergyLossForExtrapolator WARNING:index i= "
187 << i << " is out of table - NO extrapolation" << G4endl;
188 } else {
189 flag = true;
190 currentMaterial = mat;
191 electronDensity = mat->GetElectronDensity();
192 radLength = mat->GetRadlen();
193 index = i;
194 }
195 }
196 if(flag || kinEnergy != kineticEnergy) {
197 kineticEnergy = kinEnergy;
198 G4double tau = kinEnergy/mass;
199
200 gam = tau + 1.0;
201 bg2 = tau * (tau + 2.0);
202 beta2 = bg2/(gam*gam);
203 tmax = kinEnergy;
204 if(part == electron) tmax *= 0.5;
205 else if(part != positron) {
206 G4double r = electron_mass_c2/mass;
207 tmax = 2.0*bg2*electron_mass_c2/(1.0 + 2.0*gam*r + r*r);
208 }
209 if(tmax > maxEnergyTransfer) tmax = maxEnergyTransfer;
210 }
211 return true;
212}
213
214//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
215
216void G4EnergyLossForExtrapolator::Initialisation()
217{
218 isInitialised = true;
219 if(verbose>1)
220 G4cout << "### G4EnergyLossForExtrapolator::Initialisation" << G4endl;
221 currentParticle = 0;
222 currentMaterial = 0;
223 kineticEnergy = 0.0;
224 electron = G4Electron::Electron();
225 positron = G4Positron::Positron();
226 proton = G4Proton::Proton();
227 muonPlus = G4MuonPlus::MuonPlus();
228 muonMinus= G4MuonMinus::MuonMinus();
229
230 currentParticleName = "";
231
232 linLossLimit = 0.001;
233 emin = 1.*MeV;
234 emax = 10.*TeV;
235 nbins = 70;
236
237 nmat = G4Material::GetNumberOfMaterials();
238 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
239 cuts = new G4ProductionCuts();
240
241 const G4MaterialCutsCouple* couple;
242 for(G4int i=0; i<nmat; i++) {
243 couple = new G4MaterialCutsCouple((*mtable)[i],cuts);
244 couples.push_back(couple);
245 }
246
247 dedxElectron = PrepareTable();
248 dedxPositron = PrepareTable();
249 dedxMuon = PrepareTable();
250 dedxProton = PrepareTable();
251 rangeElectron = PrepareTable();
252 rangePositron = PrepareTable();
253 rangeMuon = PrepareTable();
254 rangeProton = PrepareTable();
255 invRangeElectron = PrepareTable();
256 invRangePositron = PrepareTable();
257 invRangeMuon = PrepareTable();
258 invRangeProton = PrepareTable();
259 mscElectron = PrepareTable();
260
261 G4LossTableBuilder builder;
262
263 if(verbose>1)
264 G4cout << "### G4EnergyLossForExtrapolator Builds electron tables" << G4endl;
265
266 ComputeElectronDEDX(electron, dedxElectron);
267 builder.BuildRangeTable(dedxElectron,rangeElectron);
268 builder.BuildInverseRangeTable(rangeElectron, invRangeElectron);
269
270 if(verbose>1)
271 G4cout << "### G4EnergyLossForExtrapolator Builds positron tables" << G4endl;
272
273 ComputeElectronDEDX(positron, dedxPositron);
274 builder.BuildRangeTable(dedxPositron, rangePositron);
275 builder.BuildInverseRangeTable(rangePositron, invRangePositron);
276
277 if(verbose>1)
278 G4cout << "### G4EnergyLossForExtrapolator Builds muon tables" << G4endl;
279
280 ComputeMuonDEDX(muonPlus, dedxMuon);
281 builder.BuildRangeTable(dedxMuon, rangeMuon);
282 builder.BuildInverseRangeTable(rangeMuon, invRangeMuon);
283
284 if(verbose>1)
285 G4cout << "### G4EnergyLossForExtrapolator Builds proton tables" << G4endl;
286
287 ComputeProtonDEDX(proton, dedxProton);
288 builder.BuildRangeTable(dedxProton, rangeProton);
289 builder.BuildInverseRangeTable(rangeProton, invRangeProton);
290
291 ComputeTrasportXS(electron, mscElectron);
292}
293
294//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
295
296G4PhysicsTable* G4EnergyLossForExtrapolator::PrepareTable()
297{
298 G4PhysicsTable* table = new G4PhysicsTable();
299
300 for(G4int i=0; i<nmat; i++) {
301
302 G4PhysicsVector* v = new G4PhysicsLogVector(emin, emax, nbins);
303 v->SetSpline(G4LossTableManager::Instance()->SplineFlag());
304 table->push_back(v);
305 }
306 return table;
307}
308
309//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
310
311const G4ParticleDefinition* G4EnergyLossForExtrapolator::FindParticle(const G4String& name)
312{
313 const G4ParticleDefinition* p = 0;
314 if(name != currentParticleName) {
315 p = G4ParticleTable::GetParticleTable()->FindParticle(name);
316 if(!p) {
317 G4cout << "### G4EnergyLossForExtrapolator WARNING:FindParticle fails to find "
318 << name << G4endl;
319 }
320 } else {
321 p = currentParticle;
322 }
323 return p;
324}
325
326//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
327
328G4double G4EnergyLossForExtrapolator::ComputeDEDX(G4double kinEnergy,
329 const G4ParticleDefinition* part)
330{
331 G4double x = 0.0;
332 if(part == electron) x = ComputeValue(kinEnergy, dedxElectron);
333 else if(part == positron) x = ComputeValue(kinEnergy, dedxPositron);
334 else if(part == muonPlus || part == muonMinus) {
335 x = ComputeValue(kinEnergy, dedxMuon);
336 } else {
337 G4double e = kinEnergy*proton_mass_c2/mass;
338 x = ComputeValue(e, dedxProton)*charge2;
339 }
340 return x;
341}
342
343//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
344
345G4double G4EnergyLossForExtrapolator::ComputeRange(G4double kinEnergy,
346 const G4ParticleDefinition* part)
347{
348 G4double x = 0.0;
349 if(part == electron) x = ComputeValue(kinEnergy, rangeElectron);
350 else if(part == positron) x = ComputeValue(kinEnergy, rangePositron);
351 else if(part == muonPlus || part == muonMinus)
352 x = ComputeValue(kinEnergy, rangeMuon);
353 else {
354 G4double massratio = proton_mass_c2/mass;
355 G4double e = kinEnergy*massratio;
356 x = ComputeValue(e, rangeProton)/(charge2*massratio);
357 }
358 return x;
359}
360
361//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
362
363G4double G4EnergyLossForExtrapolator::ComputeEnergy(G4double range,
364 const G4ParticleDefinition* part)
365{
366 G4double x = 0.0;
367 if(part == electron) x = ComputeValue(range, invRangeElectron);
368 else if(part == positron) x = ComputeValue(range, invRangePositron);
369 else if(part == muonPlus || part == muonMinus)
370 x = ComputeValue(range, invRangeMuon);
371 else {
372 G4double massratio = proton_mass_c2/mass;
373 G4double r = range*massratio*charge2;
374 x = ComputeValue(r, invRangeProton)/massratio;
375 }
376 return x;
377}
378
379//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
380
381void G4EnergyLossForExtrapolator::ComputeElectronDEDX(const G4ParticleDefinition* part,
382 G4PhysicsTable* table)
383{
384 G4DataVector v;
385 G4MollerBhabhaModel* ioni = new G4MollerBhabhaModel();
386 G4eBremsstrahlungModel* brem = new G4eBremsstrahlungModel();
387 ioni->Initialise(part, v);
388 brem->Initialise(part, v);
389
390 mass = electron_mass_c2;
391 charge2 = 1.0;
392 currentParticle = part;
393
394 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
395 if(0<verbose) {
396 G4cout << "G4EnergyLossForExtrapolator::ComputeElectronDEDX for "
397 << part->GetParticleName()
398 << G4endl;
399 }
400 for(G4int i=0; i<nmat; i++) {
401
402 const G4Material* mat = (*mtable)[i];
403 if(1<verbose)
404 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
405 const G4MaterialCutsCouple* couple = couples[i];
406 G4PhysicsVector* aVector = (*table)[i];
407
408 for(G4int j=0; j<nbins; j++) {
409
410 G4double e = aVector->GetLowEdgeEnergy(j);
411 G4double dedx = ioni->ComputeDEDX(couple,part,e,e) + brem->ComputeDEDX(couple,part,e,e);
412 if(1<verbose) {
413 G4cout << "j= " << j
414 << " e(MeV)= " << e/MeV
415 << " dedx(Mev/cm)= " << dedx*cm/MeV
416 << " dedx(Mev.cm2/g)= " << dedx/((MeV*mat->GetDensity())/(g/cm2)) << G4endl;
417 }
418 aVector->PutValue(j,dedx);
419 }
420 }
421 delete ioni;
422 delete brem;
423}
424
425//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
426
427void G4EnergyLossForExtrapolator::ComputeMuonDEDX(const G4ParticleDefinition* part,
428 G4PhysicsTable* table)
429{
430 G4DataVector v;
431 G4BetheBlochModel* ioni = new G4BetheBlochModel();
432 G4MuPairProductionModel* pair = new G4MuPairProductionModel();
433 G4MuBremsstrahlungModel* brem = new G4MuBremsstrahlungModel();
434 ioni->Initialise(part, v);
435 pair->Initialise(part, v);
436 brem->Initialise(part, v);
437
438 mass = part->GetPDGMass();
439 charge2 = 1.0;
440 currentParticle = part;
441
442 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
443
444 if(0<verbose) {
445 G4cout << "G4EnergyLossForExtrapolator::ComputeMuonDEDX for " << part->GetParticleName()
446 << G4endl;
447 }
448
449 for(G4int i=0; i<nmat; i++) {
450
451 const G4Material* mat = (*mtable)[i];
452 if(1<verbose)
453 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
454 const G4MaterialCutsCouple* couple = couples[i];
455 G4PhysicsVector* aVector = (*table)[i];
456 for(G4int j=0; j<nbins; j++) {
457
458 G4double e = aVector->GetLowEdgeEnergy(j);
459 G4double dedx = ioni->ComputeDEDX(couple,part,e,e) +
460 pair->ComputeDEDX(couple,part,e,e) +
461 brem->ComputeDEDX(couple,part,e,e);
462 aVector->PutValue(j,dedx);
463 if(1<verbose) {
464 G4cout << "j= " << j
465 << " e(MeV)= " << e/MeV
466 << " dedx(Mev/cm)= " << dedx*cm/MeV
467 << " dedx(Mev/(g/cm2)= " << dedx/((MeV*mat->GetDensity())/(g/cm2)) << G4endl;
468 }
469 }
470 }
471 delete ioni;
472}
473
474//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
475
476void G4EnergyLossForExtrapolator::ComputeProtonDEDX(const G4ParticleDefinition* part,
477 G4PhysicsTable* table)
478{
479 G4DataVector v;
480 G4BetheBlochModel* ioni = new G4BetheBlochModel();
481 ioni->Initialise(part, v);
482
483 mass = part->GetPDGMass();
484 charge2 = 1.0;
485 currentParticle = part;
486
487 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
488
489 if(0<verbose) {
490 G4cout << "G4EnergyLossForExtrapolator::ComputeProtonDEDX for " << part->GetParticleName()
491 << G4endl;
492 }
493
494 for(G4int i=0; i<nmat; i++) {
495
496 const G4Material* mat = (*mtable)[i];
497 if(1<verbose)
498 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
499 const G4MaterialCutsCouple* couple = couples[i];
500 G4PhysicsVector* aVector = (*table)[i];
501 for(G4int j=0; j<nbins; j++) {
502
503 G4double e = aVector->GetLowEdgeEnergy(j);
504 G4double dedx = ioni->ComputeDEDX(couple,part,e,e);
505 aVector->PutValue(j,dedx);
506 if(1<verbose) {
507 G4cout << "j= " << j
508 << " e(MeV)= " << e/MeV
509 << " dedx(Mev/cm)= " << dedx*cm/MeV
510 << " dedx(Mev.cm2/g)= " << dedx/((mat->GetDensity())/(g/cm2)) << G4endl;
511 }
512 }
513 }
514 delete ioni;
515}
516
517//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
518
519void G4EnergyLossForExtrapolator::ComputeTrasportXS(const G4ParticleDefinition* part,
520 G4PhysicsTable* table)
521{
522 G4DataVector v;
523 G4WentzelVIModel* msc = new G4WentzelVIModel();
524 msc->SetPolarAngleLimit(CLHEP::pi);
525 msc->Initialise(part, v);
526
527 mass = part->GetPDGMass();
528 charge2 = 1.0;
529 currentParticle = part;
530
531 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
532
533 if(0<verbose) {
534 G4cout << "G4EnergyLossForExtrapolator::ComputeProtonDEDX for " << part->GetParticleName()
535 << G4endl;
536 }
537
538 for(G4int i=0; i<nmat; i++) {
539
540 const G4Material* mat = (*mtable)[i];
541 if(1<verbose)
542 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
543 G4PhysicsVector* aVector = (*table)[i];
544 for(G4int j=0; j<nbins; j++) {
545
546 G4double e = aVector->GetLowEdgeEnergy(j);
547 G4double xs = msc->CrossSectionPerVolume(mat,part,e);
548 aVector->PutValue(j,xs);
549 if(1<verbose) {
550 G4cout << "j= " << j << " e(MeV)= " << e/MeV
551 << " xs(1/mm)= " << xs*mm << G4endl;
552 }
553 }
554 }
555 delete msc;
556}
557
558//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
559
Note: See TracBrowser for help on using the repository browser.