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

Last change on this file since 877 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 16.0 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.13 2007/07/28 13:44:25 vnivanch Exp $
27// GEANT4 tag $Name: $
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
70//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71
72G4EnergyLossForExtrapolator::G4EnergyLossForExtrapolator(G4int verb)
73 :maxEnergyTransfer(DBL_MAX),verbose(verb),isInitialised(false)
74{}
75
76//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
77
78G4EnergyLossForExtrapolator:: ~G4EnergyLossForExtrapolator()
79{
80 for(G4int i=0; i<nmat; i++) {delete couples[i];}
81 delete dedxElectron;
82 delete dedxPositron;
83 delete dedxProton;
84 delete rangeElectron;
85 delete rangePositron;
86 delete rangeProton;
87 delete invRangeElectron;
88 delete invRangePositron;
89 delete invRangeProton;
90 delete cuts;
91}
92
93//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
94
95G4double G4EnergyLossForExtrapolator::EnergyAfterStep(G4double kinEnergy,
96 G4double stepLength,
97 const G4Material* mat,
98 const G4ParticleDefinition* part)
99{
100 if(!isInitialised) Initialisation();
101 G4double kinEnergyFinal = kinEnergy;
102 if(mat && part) {
103 G4double step = ComputeTrueStep(mat,part,kinEnergy,stepLength);
104 G4double r = ComputeRange(kinEnergy,part);
105 if(r <= step) {
106 kinEnergyFinal = 0.0;
107 } else if(step < linLossLimit*r) {
108 kinEnergyFinal -= step*ComputeDEDX(kinEnergy,part);
109 } else {
110 G4double r1 = r - step;
111 kinEnergyFinal = ComputeEnergy(r1,part);
112 }
113 }
114 return kinEnergyFinal;
115}
116
117//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
118
119G4double G4EnergyLossForExtrapolator::EnergyBeforeStep(G4double kinEnergy,
120 G4double stepLength,
121 const G4Material* mat,
122 const G4ParticleDefinition* part)
123{
124 if(!isInitialised) Initialisation();
125 G4double kinEnergyFinal = kinEnergy;
126
127 if(mat && part) {
128 G4double step = ComputeTrueStep(mat,part,kinEnergy,stepLength);
129 G4double r = ComputeRange(kinEnergy,part);
130
131 if(step < linLossLimit*r) {
132 kinEnergyFinal += step*ComputeDEDX(kinEnergy,part);
133 } else {
134 G4double r1 = r + step;
135 kinEnergyFinal = ComputeEnergy(r1,part);
136 }
137 }
138 return kinEnergyFinal;
139}
140
141//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
142
143G4double G4EnergyLossForExtrapolator::ComputeTrueStep(const G4Material* mat,
144 const G4ParticleDefinition* part,
145 G4double kinEnergy, G4double stepLength)
146{
147 if(!isInitialised) Initialisation();
148 G4bool flag = false;
149 if(part != currentParticle) {
150 flag = true;
151 currentParticle = part;
152 mass = part->GetPDGMass();
153 G4double q = part->GetPDGCharge()/eplus;
154 charge2 = q*q;
155 }
156 if(mat != currentMaterial) {
157 G4int i = mat->GetIndex();
158 if(i >= nmat) {
159 G4cout << "### G4EnergyLossForExtrapolator WARNING:index i= "
160 << i << " is out of table - NO extrapolation" << G4endl;
161 } else {
162 flag = true;
163 currentMaterial = mat;
164 electronDensity = mat->GetElectronDensity();
165 radLength = mat->GetRadlen();
166 index = i;
167 }
168 }
169 if(flag || kinEnergy != kineticEnergy) {
170 kineticEnergy = kinEnergy;
171 G4double tau = kinEnergy/mass;
172
173 gam = tau + 1.0;
174 bg2 = tau * (tau + 2.0);
175 beta2 = bg2/(gam*gam);
176 tmax = kinEnergy;
177 if(part == electron) tmax *= 0.5;
178 else if(part != positron) {
179 G4double r = electron_mass_c2/mass;
180 tmax = 2.0*bg2*electron_mass_c2/(1.0 + 2.0*gam*r + r*r);
181 }
182 if(tmax > maxEnergyTransfer) tmax = maxEnergyTransfer;
183 }
184 G4double theta = ComputeScatteringAngle(stepLength);
185 return stepLength*std::sqrt(1.0 + 0.625*theta*theta);
186}
187
188//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
189
190void G4EnergyLossForExtrapolator::Initialisation()
191{
192 isInitialised = true;
193 if(verbose>1)
194 G4cout << "### G4EnergyLossForExtrapolator::Initialisation" << G4endl;
195 currentParticle = 0;
196 currentMaterial = 0;
197 kineticEnergy = 0.0;
198 electron = G4Electron::Electron();
199 positron = G4Positron::Positron();
200 proton = G4Proton::Proton();
201 muonPlus = G4MuonPlus::MuonPlus();
202 muonMinus= G4MuonMinus::MuonMinus();
203
204 currentParticleName = "";
205
206 linLossLimit = 0.001;
207 emin = 1.*MeV;
208 emax = 10.*TeV;
209 nbins = 70;
210
211 nmat = G4Material::GetNumberOfMaterials();
212 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
213 cuts = new G4ProductionCuts();
214
215 const G4MaterialCutsCouple* couple;
216 for(G4int i=0; i<nmat; i++) {
217 couple = new G4MaterialCutsCouple((*mtable)[i],cuts);
218 couples.push_back(couple);
219 }
220
221 dedxElectron = PrepareTable();
222 dedxPositron = PrepareTable();
223 dedxMuon = PrepareTable();
224 dedxProton = PrepareTable();
225 rangeElectron = PrepareTable();
226 rangePositron = PrepareTable();
227 rangeMuon = PrepareTable();
228 rangeProton = PrepareTable();
229 invRangeElectron = PrepareTable();
230 invRangePositron = PrepareTable();
231 invRangeMuon = PrepareTable();
232 invRangeProton = PrepareTable();
233
234 G4LossTableBuilder builder;
235
236 if(verbose>1)
237 G4cout << "### G4EnergyLossForExtrapolator Builds electron tables" << G4endl;
238
239 ComputeElectronDEDX(electron, dedxElectron);
240 builder.BuildRangeTable(dedxElectron,rangeElectron);
241 builder.BuildInverseRangeTable(rangeElectron, invRangeElectron);
242
243 if(verbose>1)
244 G4cout << "### G4EnergyLossForExtrapolator Builds positron tables" << G4endl;
245
246 ComputeElectronDEDX(positron, dedxPositron);
247 builder.BuildRangeTable(dedxPositron, rangePositron);
248 builder.BuildInverseRangeTable(rangePositron, invRangePositron);
249
250 if(verbose>1)
251 G4cout << "### G4EnergyLossForExtrapolator Builds muon tables" << G4endl;
252
253 ComputeMuonDEDX(muonPlus, dedxMuon);
254 builder.BuildRangeTable(dedxMuon, rangeMuon);
255 builder.BuildInverseRangeTable(rangeMuon, invRangeMuon);
256
257 if(verbose>1)
258 G4cout << "### G4EnergyLossForExtrapolator Builds proton tables" << G4endl;
259
260 ComputeProtonDEDX(proton, dedxProton);
261 builder.BuildRangeTable(dedxProton, rangeProton);
262 builder.BuildInverseRangeTable(rangeProton, invRangeProton);
263
264}
265
266//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
267
268G4PhysicsTable* G4EnergyLossForExtrapolator::PrepareTable()
269{
270 G4PhysicsTable* table = new G4PhysicsTable();
271
272 for(G4int i=0; i<nmat; i++) {
273
274 G4PhysicsVector* v = new G4PhysicsLogVector(emin, emax, nbins);
275 table->push_back(v);
276 }
277 return table;
278}
279
280//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
281
282const G4ParticleDefinition* G4EnergyLossForExtrapolator::FindParticle(const G4String& name)
283{
284 const G4ParticleDefinition* p = 0;
285 if(name != currentParticleName) {
286 p = G4ParticleTable::GetParticleTable()->FindParticle(name);
287 if(!p) {
288 G4cout << "### G4EnergyLossForExtrapolator WARNING:FindParticle fails to find "
289 << name << G4endl;
290 }
291 } else {
292 p = currentParticle;
293 }
294 return p;
295}
296
297//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
298
299G4double G4EnergyLossForExtrapolator::ComputeDEDX(G4double kinEnergy,
300 const G4ParticleDefinition* part)
301{
302 G4double x = 0.0;
303 if(part == electron) x = ComputeValue(kinEnergy, dedxElectron);
304 else if(part == positron) x = ComputeValue(kinEnergy, dedxPositron);
305 else if(part == muonPlus || part == muonMinus) {
306 x = ComputeValue(kinEnergy, dedxMuon);
307 } else {
308 G4double e = kinEnergy*proton_mass_c2/mass;
309 x = ComputeValue(e, dedxProton)*charge2;
310 }
311 return x;
312}
313
314//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
315
316G4double G4EnergyLossForExtrapolator::ComputeRange(G4double kinEnergy,
317 const G4ParticleDefinition* part)
318{
319 G4double x = 0.0;
320 if(part == electron) x = ComputeValue(kinEnergy, rangeElectron);
321 else if(part == positron) x = ComputeValue(kinEnergy, rangePositron);
322 else if(part == muonPlus || part == muonMinus)
323 x = ComputeValue(kinEnergy, rangeMuon);
324 else {
325 G4double massratio = proton_mass_c2/mass;
326 G4double e = kinEnergy*massratio;
327 x = ComputeValue(e, rangeProton)/(charge2*massratio);
328 }
329 return x;
330}
331
332//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
333
334G4double G4EnergyLossForExtrapolator::ComputeEnergy(G4double range,
335 const G4ParticleDefinition* part)
336{
337 G4double x = 0.0;
338 if(part == electron) x = ComputeValue(range, invRangeElectron);
339 else if(part == positron) x = ComputeValue(range, invRangePositron);
340 else if(part == muonPlus || part == muonMinus)
341 x = ComputeValue(range, invRangeMuon);
342 else {
343 G4double massratio = proton_mass_c2/mass;
344 G4double r = range*massratio*charge2;
345 x = ComputeValue(r, invRangeProton)/massratio;
346 }
347 return x;
348}
349
350//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
351
352void G4EnergyLossForExtrapolator::ComputeElectronDEDX(const G4ParticleDefinition* part,
353 G4PhysicsTable* table)
354{
355 G4DataVector v;
356 G4MollerBhabhaModel* ioni = new G4MollerBhabhaModel();
357 G4eBremsstrahlungModel* brem = new G4eBremsstrahlungModel();
358 ioni->Initialise(part, v);
359 brem->Initialise(part, v);
360
361 mass = electron_mass_c2;
362 charge2 = 1.0;
363 currentParticle = part;
364
365 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
366 if(0<verbose) {
367 G4cout << "G4EnergyLossForExtrapolator::ComputeElectronDEDX for "
368 << part->GetParticleName()
369 << G4endl;
370 }
371 for(G4int i=0; i<nmat; i++) {
372
373 const G4Material* mat = (*mtable)[i];
374 if(1<verbose)
375 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
376 const G4MaterialCutsCouple* couple = couples[i];
377 G4PhysicsVector* aVector = (*table)[i];
378
379 for(G4int j=0; j<nbins; j++) {
380
381 G4double e = aVector->GetLowEdgeEnergy(j);
382 G4double dedx = ioni->ComputeDEDX(couple,part,e,e) + brem->ComputeDEDX(couple,part,e,e);
383 if(1<verbose) {
384 G4cout << "j= " << j
385 << " e(MeV)= " << e/MeV
386 << " dedx(Mev/cm)= " << dedx*cm/MeV
387 << " dedx(Mev.cm2/g)= " << dedx/((MeV*mat->GetDensity())/(g/cm2)) << G4endl;
388 }
389 aVector->PutValue(j,dedx);
390 }
391 }
392 delete ioni;
393 delete brem;
394}
395
396//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
397
398void G4EnergyLossForExtrapolator::ComputeMuonDEDX(const G4ParticleDefinition* part,
399 G4PhysicsTable* table)
400{
401 G4DataVector v;
402 G4BetheBlochModel* ioni = new G4BetheBlochModel();
403 G4MuPairProductionModel* pair = new G4MuPairProductionModel();
404 G4MuBremsstrahlungModel* brem = new G4MuBremsstrahlungModel();
405 ioni->Initialise(part, v);
406 pair->Initialise(part, v);
407 brem->Initialise(part, v);
408
409 mass = part->GetPDGMass();
410 charge2 = 1.0;
411 currentParticle = part;
412
413 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
414
415 if(0<verbose) {
416 G4cout << "G4EnergyLossForExtrapolator::ComputeMuonDEDX for " << part->GetParticleName()
417 << G4endl;
418 }
419
420 for(G4int i=0; i<nmat; i++) {
421
422 const G4Material* mat = (*mtable)[i];
423 if(1<verbose)
424 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
425 const G4MaterialCutsCouple* couple = couples[i];
426 G4PhysicsVector* aVector = (*table)[i];
427 for(G4int j=0; j<nbins; j++) {
428
429 G4double e = aVector->GetLowEdgeEnergy(j);
430 G4double dedx = ioni->ComputeDEDX(couple,part,e,e) +
431 pair->ComputeDEDX(couple,part,e,e) +
432 brem->ComputeDEDX(couple,part,e,e);
433 aVector->PutValue(j,dedx);
434 if(1<verbose) {
435 G4cout << "j= " << j
436 << " e(MeV)= " << e/MeV
437 << " dedx(Mev/cm)= " << dedx*cm/MeV
438 << " dedx(Mev/(g/cm2)= " << dedx/((MeV*mat->GetDensity())/(g/cm2)) << G4endl;
439 }
440 }
441 }
442 delete ioni;
443}
444
445//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
446
447void G4EnergyLossForExtrapolator::ComputeProtonDEDX(const G4ParticleDefinition* part,
448 G4PhysicsTable* table)
449{
450 G4DataVector v;
451 G4BetheBlochModel* ioni = new G4BetheBlochModel();
452 ioni->Initialise(part, v);
453
454 mass = part->GetPDGMass();
455 charge2 = 1.0;
456 currentParticle = part;
457
458 const G4MaterialTable* mtable = G4Material::GetMaterialTable();
459
460 if(0<verbose) {
461 G4cout << "G4EnergyLossForExtrapolator::ComputeProtonDEDX for " << part->GetParticleName()
462 << G4endl;
463 }
464
465 for(G4int i=0; i<nmat; i++) {
466
467 const G4Material* mat = (*mtable)[i];
468 if(1<verbose)
469 G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
470 const G4MaterialCutsCouple* couple = couples[i];
471 G4PhysicsVector* aVector = (*table)[i];
472 for(G4int j=0; j<nbins; j++) {
473
474 G4double e = aVector->GetLowEdgeEnergy(j);
475 G4double dedx = ioni->ComputeDEDX(couple,part,e,e);
476 aVector->PutValue(j,dedx);
477 if(1<verbose) {
478 G4cout << "j= " << j
479 << " e(MeV)= " << e/MeV
480 << " dedx(Mev/cm)= " << dedx*cm/MeV
481 << " dedx(Mev.cm2/g)= " << dedx/((mat->GetDensity())/(g/cm2)) << G4endl;
482 }
483 }
484 }
485 delete ioni;
486}
487
488//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
489
Note: See TracBrowser for help on using the repository browser.