ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4TablesForExtrapolator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4TablesForExtrapolator.cc
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 //---------------------------------------------------------------------------
27 //
28 // ClassName: G4TablesForExtrapolator
29 //
30 // Description: This class provide calculation of energy loss, fluctuation,
31 // and msc angle
32 //
33 // Author: 09.12.04 V.Ivanchenko
34 //
35 // Modification:
36 // 08-04-05 Rename Propogator -> Extrapolator (V.Ivanchenko)
37 // 16-03-06 Add muon tables and fix bug in units (V.Ivanchenko)
38 // 21-03-06 Add verbosity defined in the constructor and Initialisation
39 // start only when first public method is called (V.Ivanchenko)
40 // 03-05-06 Remove unused pointer G4Material* from number of methods (VI)
41 // 12-05-06 SEt linLossLimit=0.001 (VI)
42 //
43 //----------------------------------------------------------------------------
44 //
45 
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47 
49 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "G4PhysicsLogVector.hh"
52 #include "G4ParticleDefinition.hh"
53 #include "G4Material.hh"
54 #include "G4MaterialCutsCouple.hh"
55 #include "G4Electron.hh"
56 #include "G4Positron.hh"
57 #include "G4Proton.hh"
58 #include "G4MuonPlus.hh"
59 #include "G4MuonMinus.hh"
60 #include "G4EmParameters.hh"
61 #include "G4MollerBhabhaModel.hh"
62 #include "G4BetheBlochModel.hh"
66 #include "G4ProductionCuts.hh"
67 #include "G4LossTableBuilder.hh"
68 #include "G4WentzelVIModel.hh"
69 
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71 
74  : builder(nullptr), verbose(verb), nbins(bins), nmat(0), emin(e1), emax(e2)
75 {
84  pcuts = nullptr;
86 }
87 
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89 
91 {
92  for(G4int i=0; i<nmat; i++) { delete couples[i]; }
93 
107 
108  delete dedxElectron;
109  delete dedxPositron;
110  delete dedxProton;
111  delete dedxMuon;
112  delete rangeElectron;
113  delete rangePositron;
114  delete rangeProton;
115  delete rangeMuon;
116  delete invRangeElectron;
117  delete invRangePositron;
118  delete invRangeProton;
119  delete invRangeMuon;
120  delete mscElectron;
121  delete pcuts;
122  delete builder;
123 }
124 
125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
126 
127 const G4PhysicsTable*
129 {
130  const G4PhysicsTable* table = nullptr;
131  switch (type)
132  {
133  case fDedxElectron:
134  table = dedxElectron;
135  break;
136  case fDedxPositron:
137  table = dedxPositron;
138  break;
139  case fDedxProton:
140  table = dedxProton;
141  break;
142  case fDedxMuon:
143  table = dedxMuon;
144  break;
145  case fRangeElectron:
146  table = rangeElectron;
147  break;
148  case fRangePositron:
149  table = rangePositron;
150  break;
151  case fRangeProton:
152  table = rangeProton;
153  break;
154  case fRangeMuon:
155  table = rangeMuon;
156  break;
157  case fInvRangeElectron:
158  table = invRangeElectron;
159  break;
160  case fInvRangePositron:
161  table = invRangePositron;
162  break;
163  case fInvRangeProton:
164  table = invRangeProton;
165  break;
166  case fInvRangeMuon:
167  table = invRangeMuon;
168  break;
169  case fMscElectron:
170  table = mscElectron;
171  }
172  return table;
173 }
174 
175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
176 
178 {
179  if(verbose>1) {
180  G4cout << "### G4TablesForExtrapolator::Initialisation" << G4endl;
181  }
182  currentParticle = nullptr;
183  mass = charge2 = 0.0;
184 
187  if(!pcuts) { pcuts = new G4ProductionCuts(); }
188 
189  G4int i0 = couples.size();
190  if(0 == i0) {
191  couples.reserve(nmat);
192  cuts.reserve(nmat);
193  }
194  for(G4int i=i0; i<nmat; ++i) {
195  couples.push_back(new G4MaterialCutsCouple((*mtable)[i],pcuts));
196  cuts.push_back(DBL_MAX);
197  }
198 
200 
214 
215  if(!builder) { builder = new G4LossTableBuilder(); }
217 
218  if(verbose>1) {
219  G4cout << "### G4TablesForExtrapolator Builds electron tables"
220  << G4endl;
221  }
225 
226  if(verbose>1) {
227  G4cout << "### G4TablesForExtrapolator Builds positron tables"
228  << G4endl;
229  }
233 
234  if(verbose>1) {
235  G4cout << "### G4TablesForExtrapolator Builds muon tables" << G4endl;
236  }
240  /*
241  G4cout << "DEDX MUON" << G4endl
242  G4cout << *dedxMuon << G4endl;
243  G4cout << "RANGE MUON" << G4endl
244  G4cout << *rangeMuon << G4endl;
245  G4cout << "INVRANGE MUON" << G4endl
246  G4cout << *invRangeMuon << G4endl;
247  */
248  if(verbose>1) {
249  G4cout << "### G4TablesForExtrapolator Builds proton tables"
250  << G4endl;
251  }
255 
257 }
258 
259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
260 
262 {
263  G4PhysicsTable* table;
264  G4int i0 = 0;
265  if(ptr) {
266  table = ptr;
267  i0 = ptr->size();
268  } else {
269  table = new G4PhysicsTable();
270  }
271  for(G4int i=i0; i<nmat; ++i) {
273  v->SetSpline(splineFlag);
274  table->push_back(v);
275  }
276  return table;
277 }
278 
279 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
280 
282  const G4ParticleDefinition* part,
283  G4PhysicsTable* table)
284 {
287  ioni->Initialise(part, cuts);
288  brem->Initialise(part, cuts);
289  ioni->SetUseBaseMaterials(false);
290  brem->SetUseBaseMaterials(false);
291 
293  charge2 = 1.0;
295 
297  if(0<verbose) {
298  G4cout << "G4TablesForExtrapolator::ComputeElectronDEDX for "
299  << part->GetParticleName()
300  << G4endl;
301  }
302  for(G4int i=0; i<nmat; ++i) {
303 
304  const G4Material* mat = (*mtable)[i];
305  if(1<verbose) {
306  G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
307  }
308  const G4MaterialCutsCouple* couple = couples[i];
309  G4PhysicsVector* aVector = (*table)[i];
310 
311  for(G4int j=0; j<=nbins; ++j) {
312 
313  G4double e = aVector->Energy(j);
314  G4double dedx = ioni->ComputeDEDX(couple,part,e,e)
315  + brem->ComputeDEDX(couple,part,e,e);
316  if(1<verbose) {
317  G4cout << "j= " << j
318  << " e(MeV)= " << e/MeV
319  << " dedx(Mev/cm)= " << dedx*cm/MeV
320  << " dedx(Mev.cm2/g)= "
321  << dedx/((MeV*mat->GetDensity())/(g/cm2))
322  << G4endl;
323  }
324  aVector->PutValue(j,dedx);
325  }
326  if(splineFlag) { aVector->FillSecondDerivatives(); }
327  }
328  delete ioni;
329  delete brem;
330 }
331 
332 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
333 
334 void
336  G4PhysicsTable* table)
337 {
338  G4BetheBlochModel* ioni = new G4BetheBlochModel();
341  ioni->Initialise(part, cuts);
342  pair->Initialise(part, cuts);
343  brem->Initialise(part, cuts);
344  ioni->SetUseBaseMaterials(false);
345  pair->SetUseBaseMaterials(false);
346  brem->SetUseBaseMaterials(false);
347 
348  mass = part->GetPDGMass();
349  charge2 = 1.0;
351 
353 
354  if(0<verbose) {
355  G4cout << "G4TablesForExtrapolator::ComputeMuonDEDX for "
356  << part->GetParticleName()
357  << G4endl;
358  }
359 
360  for(G4int i=0; i<nmat; ++i) {
361 
362  const G4Material* mat = (*mtable)[i];
363  if(1<verbose) {
364  G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
365  }
366  const G4MaterialCutsCouple* couple = couples[i];
367  G4PhysicsVector* aVector = (*table)[i];
368  for(G4int j=0; j<=nbins; j++) {
369 
370  G4double e = aVector->Energy(j);
371  G4double dedx = ioni->ComputeDEDX(couple,part,e,e) +
372  pair->ComputeDEDX(couple,part,e,e) +
373  brem->ComputeDEDX(couple,part,e,e);
374  aVector->PutValue(j,dedx);
375  if(1<verbose) {
376  G4cout << "j= " << j
377  << " e(MeV)= " << e/MeV
378  << " dedx(Mev/cm)= " << dedx*cm/MeV
379  << " dedx(Mev/(g/cm2)= "
380  << dedx/((MeV*mat->GetDensity())/(g/cm2))
381  << G4endl;
382  }
383  }
384  if(splineFlag) { aVector->FillSecondDerivatives(); }
385  }
386  delete ioni;
387 }
388 
389 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
390 
391 void
393  G4PhysicsTable* table)
394 {
395  G4BetheBlochModel* ioni = new G4BetheBlochModel();
396  ioni->Initialise(part, cuts);
397  ioni->SetUseBaseMaterials(false);
398 
399  mass = part->GetPDGMass();
400  charge2 = 1.0;
402 
404 
405  if(0<verbose) {
406  G4cout << "G4TablesForExtrapolator::ComputeProtonDEDX for "
407  << part->GetParticleName()
408  << G4endl;
409  }
410 
411  for(G4int i=0; i<nmat; ++i) {
412 
413  const G4Material* mat = (*mtable)[i];
414  if(1<verbose)
415  G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
416  const G4MaterialCutsCouple* couple = couples[i];
417  G4PhysicsVector* aVector = (*table)[i];
418  for(G4int j=0; j<=nbins; j++) {
419 
420  G4double e = aVector->Energy(j);
421  G4double dedx = ioni->ComputeDEDX(couple,part,e,e);
422  aVector->PutValue(j,dedx);
423  if(1<verbose) {
424  G4cout << "j= " << j
425  << " e(MeV)= " << e/MeV
426  << " dedx(Mev/cm)= " << dedx*cm/MeV
427  << " dedx(Mev.cm2/g)= " << dedx/((mat->GetDensity())/(g/cm2))
428  << G4endl;
429  }
430  }
431  if(splineFlag) { aVector->FillSecondDerivatives(); }
432  }
433  delete ioni;
434 }
435 
436 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
437 
438 void
440  G4PhysicsTable* table)
441 {
442  G4WentzelVIModel* msc = new G4WentzelVIModel();
444  msc->Initialise(part, cuts);
445  msc->SetUseBaseMaterials(false);
446 
447  mass = part->GetPDGMass();
448  charge2 = 1.0;
450 
452 
453  if(0<verbose) {
454  G4cout << "G4TablesForExtrapolator::ComputeProtonDEDX for "
455  << part->GetParticleName()
456  << G4endl;
457  }
458 
459  for(G4int i=0; i<nmat; ++i) {
460 
461  const G4Material* mat = (*mtable)[i];
462  msc->SetCurrentCouple(couples[i]);
463  if(1<verbose)
464  G4cout << "i= " << i << " mat= " << mat->GetName() << G4endl;
465  G4PhysicsVector* aVector = (*table)[i];
466  for(G4int j=0; j<=nbins; j++) {
467 
468  G4double e = aVector->Energy(j);
469  G4double xs = msc->CrossSectionPerVolume(mat,part,e);
470  aVector->PutValue(j,xs);
471  if(1<verbose) {
472  G4cout << "j= " << j << " e(MeV)= " << e/MeV
473  << " xs(1/mm)= " << xs*mm << G4endl;
474  }
475  }
476  if(splineFlag) { aVector->FillSecondDerivatives(); }
477  }
478  delete msc;
479 }
480 
481 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
482