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

Last change on this file since 846 was 819, checked in by garnier, 16 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.