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

Last change on this file since 961 was 961, checked in by garnier, 15 years ago

update processes
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-ref-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.