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